BiLevel-based DOE on the Sobieski SSBJ test case

from __future__ import annotations

from copy import deepcopy
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.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.

propu, aero, mission, struct = create_discipline(
    [
        "SobieskiPropulsion",
        "SobieskiAerodynamics",
        "SobieskiMission",
        "SobieskiStructure",
    ]
)

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

Then, we build a sub-scenario for each strongly coupled disciplines, using the following algorithm, maximum number of iterations and algorithm options:

algo_options = {
    "xtol_rel": 1e-7,
    "xtol_abs": 1e-7,
    "ftol_rel": 1e-7,
    "ftol_abs": 1e-7,
    "ineq_tolerance": 1e-4,
}
sub_sc_opts = {"max_iter": 30, "algo": "SLSQP", "algo_options": algo_options}

Build a sub-scenario for Propulsion

This sub-scenario will minimize SFC.

sc_prop = create_scenario(
    propu,
    "DisciplinaryOpt",
    "y_34",
    design_space=deepcopy(design_space).filter("x_3"),
    name="PropulsionScenario",
)

Build a sub-scenario for Aerodynamics

This sub-scenario will minimize L/D.

sc_aero = create_scenario(
    aero,
    "DisciplinaryOpt",
    "y_24",
    deepcopy(design_space).filter("x_2"),
    name="AerodynamicsScenario",
    maximize_objective=True,
)

Build a sub-scenario for Structure

This sub-scenario will maximize log(aircraft total weight / (aircraft total weight - fuel weight)).

sc_str = create_scenario(
    struct,
    "DisciplinaryOpt",
    "y_11",
    deepcopy(design_space).filter("x_1"),
    name="StructureScenario",
    maximize_objective=True,
)

Build a scenario for Mission

This scenario is based on the three previous sub-scenarios and on the Mission and aims to maximize the range (Breguet).

sub_disciplines = [sc_prop, sc_aero, sc_str] + [mission]
design_space = deepcopy(design_space).filter("x_shared")
system_scenario = create_scenario(
    sub_disciplines,
    "BiLevel",
    "y_4",
    design_space,
    parallel_scenarios=False,
    reset_x0_before_opt=True,
    scenario_type="DOE",
)

Note

Setting reset_x0_before_opt=True is mandatory when doing a DOE in parallel. If we want reproducible results, don’t reuse previous xopt.

system_scenario.formulation.mda1.warm_start = False
system_scenario.formulation.mda2.warm_start = False

Note

This is mandatory when doing a DOE in parallel if we want always exactly the same results, don’t warm start mda1 to have exactly the same process whatever the execution order and process dispatch.

for sub_sc in sub_disciplines[0:3]:
    sub_sc.default_inputs = {"max_iter": 20, "algo": "L-BFGS-B"}

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().

system_scenario.execute(
    {
        "n_samples": 30,
        "algo": "lhs",
        "algo_options": {"n_processes": 1 if os_name == "nt" else 4},
    }
)

system_scenario.print_execution_metrics()

 INFO - 14:46:14:
 INFO - 14:46:14: *** Start DOEScenario execution ***
 INFO - 14:46:14: DOEScenario
 INFO - 14:46:14:    Disciplines: AerodynamicsScenario PropulsionScenario SobieskiMission StructureScenario
 INFO - 14:46:14:    MDO formulation: BiLevel
 INFO - 14:46:14: Optimization problem:
 INFO - 14:46:14:    minimize y_4(x_shared)
 INFO - 14:46:14:    with respect to x_shared
 INFO - 14:46:14:    over the design space:
 INFO - 14:46:14:    +-------------+-------------+-------+-------------+-------+
 INFO - 14:46:14:    | name        | lower_bound | value | upper_bound | type  |
 INFO - 14:46:14:    +-------------+-------------+-------+-------------+-------+
 INFO - 14:46:14:    | x_shared[0] |     0.01    |  0.05 |     0.09    | float |
 INFO - 14:46:14:    | x_shared[1] |    30000    | 45000 |    60000    | float |
 INFO - 14:46:14:    | x_shared[2] |     1.4     |  1.6  |     1.8     | float |
 INFO - 14:46:14:    | x_shared[3] |     2.5     |  5.5  |     8.5     | float |
 INFO - 14:46:14:    | x_shared[4] |      40     |   55  |      70     | float |
 INFO - 14:46:14:    | x_shared[5] |     500     |  1000 |     1500    | float |
 INFO - 14:46:14:    +-------------+-------------+-------+-------------+-------+
 INFO - 14:46:14: Solving optimization problem with algorithm lhs:
 INFO - 14:46:14: ...   0%|          | 0/30 [00:00<?, ?it]
 INFO - 14:46:14: Running DOE in parallel on n_processes = 4
ERROR - 14:46:15: Failed to execute task indexed 0
ERROR - 14:46:15: math domain error
 INFO - 14:46:15: ...   3%|▎         | 1/30 [00:01<00:01, 28.24 it/sec, obj=485]
 INFO - 14:46:16: ...  10%|█         | 3/30 [00:01<00:01, 22.62 it/sec, obj=621]
 INFO - 14:46:16: ...  17%|█▋        | 5/30 [00:01<00:01, 19.14 it/sec, obj=495]
 INFO - 14:46:16: ...  23%|██▎       | 7/30 [00:01<00:01, 17.53 it/sec, obj=367]
 INFO - 14:46:16: ...  27%|██▋       | 8/30 [00:01<00:01, 16.45 it/sec, obj=549]
 INFO - 14:46:16: ...  30%|███       | 9/30 [00:02<00:01, 14.44 it/sec, obj=918]
 INFO - 14:46:17: ...  40%|████      | 12/30 [00:02<00:01, 12.20 it/sec, obj=354]
 INFO - 14:46:17: ...  50%|█████     | 15/30 [00:02<00:01, 11.54 it/sec, obj=337]
 INFO - 14:46:17: ...  57%|█████▋    | 17/30 [00:02<00:01, 10.42 it/sec, obj=2.27e+3]
 INFO - 14:46:17: ...  63%|██████▎   | 19/30 [00:03<00:01,  9.65 it/sec, obj=829]
 INFO - 14:46:18: ...  67%|██████▋   | 20/30 [00:03<00:01,  8.80 it/sec, obj=832]
 INFO - 14:46:18: ...  73%|███████▎  | 22/30 [00:03<00:00,  8.54 it/sec, obj=640]
 INFO - 14:46:18: ...  80%|████████  | 24/30 [00:03<00:00,  7.62 it/sec, obj=1.19e+3]
 INFO - 14:46:18: ...  93%|█████████▎| 28/30 [00:04<00:00,  7.19 it/sec, obj=952]
 INFO - 14:46:19: ...  97%|█████████▋| 29/30 [00:04<00:00,  7.11 it/sec, obj=647]
 INFO - 14:46:19: Optimization result:
 INFO - 14:46:19:    Optimizer info:
 INFO - 14:46:19:       Status: None
 INFO - 14:46:19:       Message: None
 INFO - 14:46:19:       Number of calls to the objective function by the optimizer: 30
 INFO - 14:46:19:    Solution:
 INFO - 14:46:19:       Objective: 292.91178826079175
 INFO - 14:46:19:       Design space:
 INFO - 14:46:19:       +-------------+-------------+---------------------+-------------+-------+
 INFO - 14:46:19:       | name        | lower_bound |        value        | upper_bound | type  |
 INFO - 14:46:19:       +-------------+-------------+---------------------+-------------+-------+
 INFO - 14:46:19:       | x_shared[0] |     0.01    | 0.04666883227502976 |     0.09    | float |
 INFO - 14:46:19:       | x_shared[1] |    30000    |  37789.27932845149  |    60000    | float |
 INFO - 14:46:19:       | x_shared[2] |     1.4     |  1.667044086506944  |     1.8     | float |
 INFO - 14:46:19:       | x_shared[3] |     2.5     |   8.29772323088249  |     8.5     | float |
 INFO - 14:46:19:       | x_shared[4] |      40     |  42.41730480236713  |      70     | float |
 INFO - 14:46:19:       | x_shared[5] |     500     |  724.9388551459947  |     1500    | float |
 INFO - 14:46:19:       +-------------+-------------+---------------------+-------------+-------+
 INFO - 14:46:19: *** End DOEScenario execution (time: 0:00:04.231293) ***
 INFO - 14:46:19: Scenario Execution Statistics
 INFO - 14:46:19:    Discipline: PropulsionScenario
 INFO - 14:46:19:       Executions number: 30
 INFO - 14:46:19:       Execution time: 0.8828918189999513 s
 INFO - 14:46:19:       Linearizations number: 0
 INFO - 14:46:19:    Discipline: AerodynamicsScenario
 INFO - 14:46:19:       Executions number: 30
 INFO - 14:46:19:       Execution time: 1.0742737650004983 s
 INFO - 14:46:19:       Linearizations number: 0
 INFO - 14:46:19:    Discipline: StructureScenario
 INFO - 14:46:19:       Executions number: 30
 INFO - 14:46:19:       Execution time: 4.105284698999412 s
 INFO - 14:46:19:       Linearizations number: 0
 INFO - 14:46:19:    Discipline: SobieskiMission
 INFO - 14:46:19:       Executions number: 29
 INFO - 14:46:19:       Execution time: 0.0016093849994831544 s
 INFO - 14:46:19:       Linearizations number: 0
 INFO - 14:46:19:    Total number of executions calls: 119
 INFO - 14:46:19:    Total number of linearizations: 0

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 the method export_to_dataframe():

dataset = system_scenario.export_to_dataset("a_name_for_my_dataset")
dataframe = dataset.export_to_dataframe()

Plot the optimization history view

system_scenario.post_process("OptHistoryView", show=False, save=False)

• 
• 
• 

<gemseo.post.opt_history_view.OptHistoryView object at 0x7f3bfb0f3f10>

Plot the scatter matrix

system_scenario.post_process(
    "ScatterPlotMatrix", show=False, save=False, variable_names=["y_4", "x_shared"]
)

<gemseo.post.scatter_mat.ScatterPlotMatrix object at 0x7f3bfb14bbe0>

Plot parallel coordinates

system_scenario.post_process("ParallelCoordinates", show=False, save=False)

• 
• 

<gemseo.post.para_coord.ParallelCoordinates object at 0x7f3c106eedc0>

Plot correlations

system_scenario.post_process("Correlations", show=False, save=False)
# Workaround for HTML rendering, instead of ``show=True``
plt.show()

INFO - 14:46:22: Detected 0 correlations > 0.95