MDF-based DOE on the Sobieski SSBJ test case

from os import name as os_name

from gemseo.api import configure_logger
from gemseo.api import create_discipline
from gemseo.api import create_scenario
from gemseo.api import generate_n2_plot
from gemseo.problems.sobieski.core.problem import SobieskiProblem
from matplotlib import pyplot as plt

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: {'y_23': array([12562.01206488]), 'x_3': array([0.5]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'c_3': array([4360.])}
SobieskiAerodynamics
Default inputs: {'x_2': array([1.]), 'y_32': array([0.50279625]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'y_12': array([5.06069742e+04, 9.50000000e-01]), 'c_4': array([0.01375])}
SobieskiMission
Default inputs: {'y_14': array([50606.9741711 ,  7306.20262124]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'y_24': array([4.15006276]), 'y_34': array([1.10754577])}
SobieskiStructure
Default inputs: {'y_21': array([50606.9741711]), 'y_31': array([6354.32430691]), 'x_1': array([0.25, 1.  ]), 'x_shared': array([5.0e-02, 4.5e+04, 1.6e+00, 5.5e+00, 5.5e+01, 1.0e+03]), 'c_0': array([2000.]), 'c_1': array([25000.]), 'c_2': array([6.])}

You may also be interested in plotting the couplings of your disciplines. A quick way of getting this information is the API 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 = SobieskiProblem().design_space
print(design_space)
Design space:
+----------+-------------+--------------------+-------------+-------+
| name     | lower_bound |       value        | upper_bound | type  |
+----------+-------------+--------------------+-------------+-------+
| x_shared |     0.01    |        0.05        |     0.09    | float |
| x_shared |    30000    |       45000        |    60000    | float |
| x_shared |     1.4     |        1.6         |     1.8     | float |
| x_shared |     2.5     |        5.5         |     8.5     | float |
| x_shared |      40     |         55         |      70     | float |
| x_shared |     500     |        1000        |     1500    | float |
| x_1      |     0.1     |        0.25        |     0.4     | float |
| x_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     |    24850    |   50606.9741711    |    77100    | float |
| y_14     |    -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     |    24850    |     50606.9742     |    77250    | float |
| y_12     |     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")

Execute the scenario

Use provided analytic derivatives

scenario.set_differentiation_method("user")

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. For Python versions < 3.7 and Numpy < 1.20.0, subprocesses may get hung randomly during execution. It is strongly recommended to update your environment to avoid this problem. 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 - 13:45:29:
    INFO - 13:45:29: *** Start DOEScenario execution ***
    INFO - 13:45:29: DOEScenario
    INFO - 13:45:29:    Disciplines: SobieskiAerodynamics SobieskiMission SobieskiPropulsion SobieskiStructure
    INFO - 13:45:29:    MDO formulation: MDF
    INFO - 13:45:29: Optimization problem:
    INFO - 13:45:29:    minimize -y_4(x_shared, x_1, x_2, x_3)
    INFO - 13:45:29:    with respect to x_1, x_2, x_3, x_shared
    INFO - 13:45:29:    subject to constraints:
    INFO - 13:45:29:       g_1(x_shared, x_1, x_2, x_3) <= 0.0
    INFO - 13:45:29:       g_2(x_shared, x_1, x_2, x_3) <= 0.0
    INFO - 13:45:29:       g_3(x_shared, x_1, x_2, x_3) <= 0.0
    INFO - 13:45:29:    over the design space:
    INFO - 13:45:29:    +----------+-------------+-------+-------------+-------+
    INFO - 13:45:29:    | name     | lower_bound | value | upper_bound | type  |
    INFO - 13:45:29:    +----------+-------------+-------+-------------+-------+
    INFO - 13:45:29:    | x_shared |     0.01    |  0.05 |     0.09    | float |
    INFO - 13:45:29:    | x_shared |    30000    | 45000 |    60000    | float |
    INFO - 13:45:29:    | x_shared |     1.4     |  1.6  |     1.8     | float |
    INFO - 13:45:29:    | x_shared |     2.5     |  5.5  |     8.5     | float |
    INFO - 13:45:29:    | x_shared |      40     |   55  |      70     | float |
    INFO - 13:45:29:    | x_shared |     500     |  1000 |     1500    | float |
    INFO - 13:45:29:    | x_1      |     0.1     |  0.25 |     0.4     | float |
    INFO - 13:45:29:    | x_1      |     0.75    |   1   |     1.25    | float |
    INFO - 13:45:29:    | x_2      |     0.75    |   1   |     1.25    | float |
    INFO - 13:45:29:    | x_3      |     0.1     |  0.5  |      1      | float |
    INFO - 13:45:29:    +----------+-------------+-------+-------------+-------+
    INFO - 13:45:29: Solving optimization problem with algorithm lhs:
    INFO - 13:45:29: ...   0%|          | 0/30 [00:00<?, ?it]
    INFO - 13:45:29: Running DOE in parallel on n_processes = 4
    INFO - 13:45:30: ...   3%|▎         | 1/30 [00:00<00:00, 43.26 it/sec, obj=-285]
    INFO - 13:45:30: ...  17%|█▋        | 5/30 [00:01<00:00, 25.80 it/sec, obj=-481]
    INFO - 13:45:31: ...  30%|███       | 9/30 [00:01<00:01, 18.36 it/sec, obj=-419]
    INFO - 13:45:31: ...  43%|████▎     | 13/30 [00:02<00:01, 14.71 it/sec, obj=-574]
    INFO - 13:45:31: ...  53%|█████▎    | 16/30 [00:02<00:00, 14.01 it/sec, obj=-433]
    INFO - 13:45:32: ...  60%|██████    | 18/30 [00:02<00:01, 11.89 it/sec, obj=-481]
    INFO - 13:45:32: ...  70%|███████   | 21/30 [00:02<00:00, 10.26 it/sec, obj=-247]
    INFO - 13:45:32: ...  77%|███████▋  | 23/30 [00:03<00:00,  9.88 it/sec, obj=-273]
    INFO - 13:45:32: ...  83%|████████▎ | 25/30 [00:03<00:00,  9.08 it/sec, obj=-624]
    INFO - 13:45:33: ...  90%|█████████ | 27/30 [00:03<00:00,  8.67 it/sec, obj=-606]
    INFO - 13:45:33: ...  97%|█████████▋| 29/30 [00:03<00:00,  8.11 it/sec, obj=-485]
    INFO - 13:45:33: ... 100%|██████████| 30/30 [00:03<00:00,  7.98 it/sec, obj=-405]
    INFO - 13:45:33: Optimization result:
    INFO - 13:45:33:    Optimizer info:
    INFO - 13:45:33:       Status: None
    INFO - 13:45:33:       Message: None
    INFO - 13:45:33:       Number of calls to the objective function by the optimizer: 30
    INFO - 13:45:33:    Solution:
    INFO - 13:45:33:       The solution is feasible.
    INFO - 13:45:33:       Objective: -485.49229968733584
    INFO - 13:45:33:       Standardized constraints:
    INFO - 13:45:33:          g_1 = [-0.11350951 -0.10812292 -0.1045109  -0.10204971 -0.10028641 -0.01838903
    INFO - 13:45:33:  -0.22161097]
    INFO - 13:45:33:          g_2 = -0.02400000000000002
    INFO - 13:45:33:          g_3 = [-0.33063134 -0.66936866 -0.73821755 -0.07789536]
    INFO - 13:45:33:       Design space:
    INFO - 13:45:33:       +----------+-------------+---------------------+-------------+-------+
    INFO - 13:45:33:       | name     | lower_bound |        value        | upper_bound | type  |
    INFO - 13:45:33:       +----------+-------------+---------------------+-------------+-------+
    INFO - 13:45:33:       | x_shared |     0.01    | 0.05400000000000001 |     0.09    | float |
    INFO - 13:45:33:       | x_shared |    30000    |        46500        |    60000    | float |
    INFO - 13:45:33:       | x_shared |     1.4     |  1.686666666666667  |     1.8     | float |
    INFO - 13:45:33:       | x_shared |     2.5     |         5.2         |     8.5     | float |
    INFO - 13:45:33:       | x_shared |      40     |         66.5        |      70     | float |
    INFO - 13:45:33:       | x_shared |     500     |  583.3333333333334  |     1500    | float |
    INFO - 13:45:33:       | x_1      |     0.1     |        0.185        |     0.4     | float |
    INFO - 13:45:33:       | x_1      |     0.75    |  0.9416666666666667 |     1.25    | float |
    INFO - 13:45:33:       | x_2      |     0.75    |        0.775        |     1.25    | float |
    INFO - 13:45:33:       | x_3      |     0.1     |        0.115        |      1      | float |
    INFO - 13:45:33:       +----------+-------------+---------------------+-------------+-------+
    INFO - 13:45:33: *** End DOEScenario execution (time: 0:00:03.772875) ***

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

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.export_to_dataset() will allow you to export your results to a Dataset, the basic GEMSEO class to store data. From a dataset, you can even obtain a Pandas dataframe with its method export_to_dataframe():

dataset = scenario.export_to_dataset("a_name_for_my_dataset")
dataframe = dataset.export_to_dataframe(
    variable_names=["-y_4", "x_1", "x_2", "x_3", "x_shared"]
)

Plot the optimization history view

scenario.post_process("OptHistoryView", show=False, save=False)
  • 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 0x7f1a500fa6a0>

Tip

Each post-processing method requires different inputs and offers a variety of customization options. Use the API 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", show=False, save=False, variable_names=["y_4", "x_shared"]
)
plot doe sobieski mdf example
<gemseo.post.scatter_mat.ScatterPlotMatrix object at 0x7f1a68d164c0>

Plot correlations

scenario.post_process("Correlations", show=False, save=False)
# Workaround for HTML rendering, instead of ``show=True``
plt.show()
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 - 13:45:36: Detected 10 correlations > 0.95

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

Gallery generated by Sphinx-Gallery