MDF-based DOE on the Sobieski SSBJ test case

from __future__ import annotations

from os import name as os_name

from gemseo import configure_logger
from gemseo import create_discipline
from gemseo import create_scenario
from gemseo import generate_n2_plot
from gemseo.problems.sobieski.core.design_space import SobieskiDesignSpace

configure_logger()
<RootLogger root (INFO)>

Instantiate the disciplines

First, we instantiate the four disciplines of the use case: SobieskiPropulsion, SobieskiAerodynamics, SobieskiMission and SobieskiStructure.

disciplines = create_discipline([
    "SobieskiPropulsion",
    "SobieskiAerodynamics",
    "SobieskiMission",
    "SobieskiStructure",
])

We can quickly access the most relevant information of any discipline (name, inputs, and outputs) with Python’s print() function. Moreover, we can get the default input values of a discipline with the attribute MDODiscipline.default_inputs

for discipline in disciplines:
    print(discipline)
    print(f"Default inputs: {discipline.default_inputs}")
SobieskiPropulsion
Default inputs: {'c_3': array([4360.]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'x_3': array([0.5]), 'y_23': array([12562.01206488])}
SobieskiAerodynamics
Default inputs: {'c_4': array([0.01375]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'y_32': array([0.50279625]), 'y_12': array([5.06069742e+04, 9.50000000e-01]), 'x_2': array([1.])}
SobieskiMission
Default inputs: {'y_24': array([4.15006276]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'y_34': array([1.10754577]), 'y_14': array([50606.9741711 ,  7306.20262124])}
SobieskiStructure
Default inputs: {'y_21': array([50606.9741711]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'y_31': array([6354.32430691]), 'c_2': array([6.]), 'x_1': array([0.25, 1.  ]), 'c_0': array([2000.]), 'c_1': array([25000.])}

You may also be interested in plotting the couplings of your disciplines. A quick way of getting this information is the high-level function generate_n2_plot(). A much more detailed explanation of coupling visualization is available here.

generate_n2_plot(disciplines, save=False, show=True)
plot doe sobieski mdf example

Build, execute and post-process the scenario

Then, we build the scenario which links the disciplines with the formulation and the optimization algorithm. Here, we use the BiLevel formulation. We tell the scenario to minimize -y_4 instead of minimizing y_4 (range), which is the default option.

We need to define the design space.

design_space = SobieskiDesignSpace()
design_space
Sobieski design space:
Name Lower bound Value Upper bound Type
x_shared[0] 0.01 0.05 0.09 float
x_shared[1] 30000 45000 60000 float
x_shared[2] 1.4 1.6 1.8 float
x_shared[3] 2.5 5.5 8.5 float
x_shared[4] 40 55 70 float
x_shared[5] 500 1000 1500 float
x_1[0] 0.1 0.25 0.4 float
x_1[1] 0.75 1 1.25 float
x_2 0.75 1 1.25 float
x_3 0.1 0.5 1 float
y_14[0] 24850 50606.9741711 77100 float
y_14[1] -7700 7306.20262124 45000 float
y_32 0.235 0.5027962499999999 0.795 float
y_31 2960 6354.32430691 10185 float
y_24 0.44 4.15006276 11.13 float
y_34 0.44 1.10754577 1.98 float
y_23 3365 12194.2671934 26400 float
y_21 24850 50606.9741711 77250 float
y_12[0] 24850 50606.9742 77250 float
y_12[1] 0.45 0.95 1.5 float


Instantiate the scenario

scenario = create_scenario(
    disciplines,
    formulation="MDF",
    objective_name="y_4",
    design_space=design_space,
    maximize_objective=True,
    scenario_type="DOE",
)

Set the design constraints

for constraint in ["g_1", "g_2", "g_3"]:
    scenario.add_constraint(constraint, "ineq")

Visualize the XDSM

Generate the XDSM on the fly:

  • log_workflow_status=True will log the status of the workflow in the console,

  • save_html (default True) will generate a self-contained HTML file, that can be automatically opened using show_html=True.

scenario.xdsmize(save_html=False)


Execute the scenario

Use provided analytic derivatives

scenario.set_differentiation_method()

Multiprocessing

It is possible to run a DOE in parallel using multiprocessing, in order to do this, we specify the number of processes to be used for the computation of the samples.

Warning

The multiprocessing option has some limitations on Windows. Due to problems with sphinx, we disable it in this example. The features MemoryFullCache and HDF5Cache are not available for multiprocessing on Windows. As an alternative, we recommend the method DOEScenario.set_optimization_history_backup().

We define the algorithm options. Here the criterion = center option of pyDOE centers the points within the sampling intervals.

algo_options = {
    "criterion": "center",
    # Evaluate gradient of the MDA
    # with coupled adjoint
    "eval_jac": True,
    # Run in parallel on 1 or 4 processors
    "n_processes": 1 if os_name == "nt" else 4,
}

scenario.execute({"n_samples": 30, "algo": "lhs", "algo_options": algo_options})
    INFO - 10:52:33:
    INFO - 10:52:33: *** Start DOEScenario execution ***
    INFO - 10:52:33: DOEScenario
    INFO - 10:52:33:    Disciplines: SobieskiAerodynamics SobieskiMission SobieskiPropulsion SobieskiStructure
    INFO - 10:52:33:    MDO formulation: MDF
    INFO - 10:52:33: Optimization problem:
    INFO - 10:52:33:    minimize -y_4(x_shared, x_1, x_2, x_3)
    INFO - 10:52:33:    with respect to x_1, x_2, x_3, x_shared
    INFO - 10:52:33:    subject to constraints:
    INFO - 10:52:33:       g_1(x_shared, x_1, x_2, x_3) <= 0.0
    INFO - 10:52:33:       g_2(x_shared, x_1, x_2, x_3) <= 0.0
    INFO - 10:52:33:       g_3(x_shared, x_1, x_2, x_3) <= 0.0
    INFO - 10:52:33:    over the design space:
    INFO - 10:52:33:       +-------------+-------------+-------+-------------+-------+
    INFO - 10:52:33:       | Name        | Lower bound | Value | Upper bound | Type  |
    INFO - 10:52:33:       +-------------+-------------+-------+-------------+-------+
    INFO - 10:52:33:       | x_shared[0] |     0.01    |  0.05 |     0.09    | float |
    INFO - 10:52:33:       | x_shared[1] |    30000    | 45000 |    60000    | float |
    INFO - 10:52:33:       | x_shared[2] |     1.4     |  1.6  |     1.8     | float |
    INFO - 10:52:33:       | x_shared[3] |     2.5     |  5.5  |     8.5     | float |
    INFO - 10:52:33:       | x_shared[4] |      40     |   55  |      70     | float |
    INFO - 10:52:33:       | x_shared[5] |     500     |  1000 |     1500    | float |
    INFO - 10:52:33:       | x_1[0]      |     0.1     |  0.25 |     0.4     | float |
    INFO - 10:52:33:       | x_1[1]      |     0.75    |   1   |     1.25    | float |
    INFO - 10:52:33:       | x_2         |     0.75    |   1   |     1.25    | float |
    INFO - 10:52:33:       | x_3         |     0.1     |  0.5  |      1      | float |
    INFO - 10:52:33:       +-------------+-------------+-------+-------------+-------+
    INFO - 10:52:33: Solving optimization problem with algorithm lhs:
    INFO - 10:52:33: Running DOE in parallel on n_processes = 4
    INFO - 10:52:34:      3%|▎         | 1/30 [00:00<00:28,  1.03 it/sec, obj=-222]
    INFO - 10:52:34:      7%|▋         | 2/30 [00:01<00:14,  2.00 it/sec, obj=-285]
    INFO - 10:52:34:     10%|█         | 3/30 [00:01<00:09,  2.96 it/sec, obj=-567]
    INFO - 10:52:34:     13%|█▎        | 4/30 [00:01<00:06,  3.84 it/sec, obj=-344]
    INFO - 10:52:35:     17%|█▋        | 5/30 [00:01<00:07,  3.13 it/sec, obj=-235]
    INFO - 10:52:35:     20%|██        | 6/30 [00:01<00:06,  3.69 it/sec, obj=-519]
    INFO - 10:52:35:     23%|██▎       | 7/30 [00:01<00:05,  4.18 it/sec, obj=-429]
    INFO - 10:52:35:     27%|██▋       | 8/30 [00:01<00:04,  4.48 it/sec, obj=-481]
    INFO - 10:52:35:     30%|███       | 9/30 [00:02<00:05,  4.17 it/sec, obj=-538]
    INFO - 10:52:36:     33%|███▎      | 10/30 [00:02<00:04,  4.31 it/sec, obj=-306]
    INFO - 10:52:36:     37%|███▋      | 11/30 [00:02<00:04,  4.71 it/sec, obj=-419]
    INFO - 10:52:36:     40%|████      | 12/30 [00:02<00:03,  5.04 it/sec, obj=-556]
    INFO - 10:52:36:     43%|████▎     | 13/30 [00:02<00:03,  4.72 it/sec, obj=-495]
    INFO - 10:52:36:     47%|████▋     | 14/30 [00:02<00:03,  4.79 it/sec, obj=-574]
    INFO - 10:52:36:     50%|█████     | 15/30 [00:02<00:02,  5.06 it/sec, obj=-433]
    INFO - 10:52:36:     53%|█████▎    | 16/30 [00:02<00:02,  5.35 it/sec, obj=-621]
    INFO - 10:52:37:     57%|█████▋    | 17/30 [00:03<00:02,  5.09 it/sec, obj=-481]
    INFO - 10:52:37:     60%|██████    | 18/30 [00:03<00:02,  5.06 it/sec, obj=-1.07e+3]
    INFO - 10:52:37:     63%|██████▎   | 19/30 [00:03<00:02,  5.31 it/sec, obj=-287]
    INFO - 10:52:37:     67%|██████▋   | 20/30 [00:03<00:01,  5.52 it/sec, obj=-602]
    INFO - 10:52:37:     70%|███████   | 21/30 [00:04<00:01,  5.23 it/sec, obj=-247]
    INFO - 10:52:37:     73%|███████▎  | 22/30 [00:04<00:01,  5.32 it/sec, obj=-273]
    INFO - 10:52:37:     77%|███████▋  | 23/30 [00:04<00:01,  5.51 it/sec, obj=-414]
    INFO - 10:52:38:     80%|████████  | 24/30 [00:04<00:01,  5.63 it/sec, obj=-1.18e+3]
    INFO - 10:52:38:     83%|████████▎ | 25/30 [00:04<00:00,  5.43 it/sec, obj=-383]
    INFO - 10:52:38:     87%|████████▋ | 26/30 [00:04<00:00,  5.56 it/sec, obj=-624]
    INFO - 10:52:38:     90%|█████████ | 27/30 [00:04<00:00,  5.62 it/sec, obj=-606]
    INFO - 10:52:38:     93%|█████████▎| 28/30 [00:04<00:00,  5.75 it/sec, obj=-351]
    INFO - 10:52:38:     97%|█████████▋| 29/30 [00:05<00:00,  5.69 it/sec, obj=-405]
    INFO - 10:52:38:    100%|██████████| 30/30 [00:05<00:00,  5.88 it/sec, obj=-485]
    INFO - 10:52:38: Optimization result:
    INFO - 10:52:38:    Optimizer info:
    INFO - 10:52:38:       Status: None
    INFO - 10:52:38:       Message: None
    INFO - 10:52:38:       Number of calls to the objective function by the optimizer: 30
    INFO - 10:52:38:    Solution:
    INFO - 10:52:38:       The solution is feasible.
    INFO - 10:52:38:       Objective: -485.49213288049987
    INFO - 10:52:38:       Standardized constraints:
    INFO - 10:52:38:          g_1 = [-0.11350951 -0.10812292 -0.1045109  -0.10204971 -0.10028641 -0.01838903
    INFO - 10:52:38:  -0.22161097]
    INFO - 10:52:38:          g_2 = -0.02400000000000002
    INFO - 10:52:38:          g_3 = [-0.33063169 -0.66936831 -0.73821755 -0.07789536]
    INFO - 10:52:38:       Design space:
    INFO - 10:52:38:          +-------------+-------------+---------------------+-------------+-------+
    INFO - 10:52:38:          | Name        | Lower bound |        Value        | Upper bound | Type  |
    INFO - 10:52:38:          +-------------+-------------+---------------------+-------------+-------+
    INFO - 10:52:38:          | x_shared[0] |     0.01    | 0.05400000000000001 |     0.09    | float |
    INFO - 10:52:38:          | x_shared[1] |    30000    |        46500        |    60000    | float |
    INFO - 10:52:38:          | x_shared[2] |     1.4     |  1.686666666666667  |     1.8     | float |
    INFO - 10:52:38:          | x_shared[3] |     2.5     |         5.2         |     8.5     | float |
    INFO - 10:52:38:          | x_shared[4] |      40     |         66.5        |      70     | float |
    INFO - 10:52:38:          | x_shared[5] |     500     |  583.3333333333334  |     1500    | float |
    INFO - 10:52:38:          | x_1[0]      |     0.1     |        0.185        |     0.4     | float |
    INFO - 10:52:38:          | x_1[1]      |     0.75    |  0.9416666666666667 |     1.25    | float |
    INFO - 10:52:38:          | x_2         |     0.75    |        0.775        |     1.25    | float |
    INFO - 10:52:38:          | x_3         |     0.1     |        0.115        |      1      | float |
    INFO - 10:52:38:          +-------------+-------------+---------------------+-------------+-------+
    INFO - 10:52:38: *** End DOEScenario execution (time: 0:00:05.181279) ***

{'eval_jac': False, 'algo_options': {'criterion': 'center', 'eval_jac': True, 'n_processes': 4}, 'n_samples': 30, 'algo': 'lhs'}

Warning

On Windows, the progress bar may show duplicated instances during the initialization of each subprocess. In some cases it may also print the conclusion of an iteration ahead of another one that was concluded first. This is a consequence of the pickling process and does not affect the computations of the scenario.

Exporting the problem data.

After the execution of the scenario, you may want to export your data to use it elsewhere. The method Scenario.to_dataset() will allow you to export your results to a Dataset, the basic GEMSEO class to store data.

dataset = scenario.to_dataset("a_name_for_my_dataset")

Plot the optimization history view

scenario.post_process("OptHistoryView", save=False, show=True)
  • Evolution of the optimization variables
  • Evolution of the objective value
  • Distance to the optimum
  • Hessian diagonal approximation
  • Evolution of the inequality constraints
<gemseo.post.opt_history_view.OptHistoryView object at 0x7f1dde513250>

Tip

Each post-processing method requires different inputs and offers a variety of customization options. Use the high-level function get_post_processing_options_schema() to print a table with the attributes for any post-processing algo. Or refer to our dedicated page: Post-processing algorithms.

Plot the scatter matrix

scenario.post_process(
    "ScatterPlotMatrix",
    save=False,
    show=True,
    variable_names=["y_4", "x_shared"],
)
plot doe sobieski mdf example
<gemseo.post.scatter_mat.ScatterPlotMatrix object at 0x7f1ddda643d0>

Plot correlations

scenario.post_process("Correlations", save=False, show=True)
R=0.98987, R=0.97511, R=0.99671, R=0.96235, R=0.99119, R=0.99866, R=0.95205, R=0.98584, R=0.99618, R=0.99936
    INFO - 10:52:42: Detected 10 correlations > 0.95

<gemseo.post.correlations.Correlations object at 0x7f1ddd56d4c0>

Total running time of the script: (0 minutes 9.393 seconds)

Gallery generated by Sphinx-Gallery