IDF-based MDO on the Sobieski SSBJ test case

from __future__ import annotations

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.problem import SobieskiProblem

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

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 sobieski idf 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 IDF formulation. We tell the scenario to minimize -y_4 instead of minimizing y_4 (range), which is the default option.

Instantiate the scenario

design_space = SobieskiProblem().design_space
design_space
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


scenario = create_scenario(
    disciplines,
    "IDF",
    objective_name="y_4",
    design_space=design_space,
    maximize_objective=True,
)

Set the design constraints

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

Visualize the XDSM

Generate the XDSM file 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)


Define the algorithm inputs

We set the maximum number of iterations, the optimizer and the optimizer options

algo_options = {
    "ftol_rel": 1e-10,
    "ineq_tolerance": 1e-3,
    "eq_tolerance": 1e-3,
    "normalize_design_space": True,
}
scn_inputs = {"max_iter": 20, "algo": "SLSQP", "algo_options": algo_options}

Execute the scenario

scenario.execute(scn_inputs)
    INFO - 08:24:31:
    INFO - 08:24:31: *** Start MDOScenario execution ***
    INFO - 08:24:31: MDOScenario
    INFO - 08:24:31:    Disciplines: SobieskiAerodynamics SobieskiMission SobieskiPropulsion SobieskiStructure
    INFO - 08:24:31:    MDO formulation: IDF
    INFO - 08:24:31: Optimization problem:
    INFO - 08:24:31:    minimize -y_4(x_shared, y_14, y_24, y_34)
    INFO - 08:24:31:    with respect to x_1, x_2, x_3, x_shared, y_12, y_14, y_21, y_23, y_24, y_31, y_32, y_34
    INFO - 08:24:31:    subject to constraints:
    INFO - 08:24:31:       g_1(x_shared, x_1, y_31, y_21) <= 0.0
    INFO - 08:24:31:       g_2(x_shared, x_2, y_32, y_12) <= 0.0
    INFO - 08:24:31:       g_3(x_shared, x_3, y_23) <= 0.0
    INFO - 08:24:31:       y_31_y_32_y_34(x_shared, x_3, y_23): y_31(x_shared, x_3, y_23) - y_31 == 0.0
    INFO - 08:24:31:                                            y_32(x_shared, x_3, y_23) - y_32 == 0.0
    INFO - 08:24:31:                                            y_34(x_shared, x_3, y_23) - y_34 == 0.0
    INFO - 08:24:31:       y_21_y_23_y_24(x_shared, x_2, y_32, y_12): y_21(x_shared, x_2, y_32, y_12) - y_21 == 0.0
    INFO - 08:24:31:                                                  y_23(x_shared, x_2, y_32, y_12) - y_23 == 0.0
    INFO - 08:24:31:                                                  y_24(x_shared, x_2, y_32, y_12) - y_24 == 0.0
    INFO - 08:24:31:       y_12_y_14(x_shared, x_1, y_31, y_21): y_12(x_shared, x_1, y_31, y_21) - y_12 == 0.0
    INFO - 08:24:31:                                             y_14(x_shared, x_1, y_31, y_21) - y_14 == 0.0
    INFO - 08:24:31:    over the design space:
    INFO - 08:24:31:    +-------------+-------------+--------------------+-------------+-------+
    INFO - 08:24:31:    | name        | lower_bound |       value        | upper_bound | type  |
    INFO - 08:24:31:    +-------------+-------------+--------------------+-------------+-------+
    INFO - 08:24:31:    | x_shared[0] |     0.01    |        0.05        |     0.09    | float |
    INFO - 08:24:31:    | x_shared[1] |    30000    |       45000        |    60000    | float |
    INFO - 08:24:31:    | x_shared[2] |     1.4     |        1.6         |     1.8     | float |
    INFO - 08:24:31:    | x_shared[3] |     2.5     |        5.5         |     8.5     | float |
    INFO - 08:24:31:    | x_shared[4] |      40     |         55         |      70     | float |
    INFO - 08:24:31:    | x_shared[5] |     500     |        1000        |     1500    | float |
    INFO - 08:24:31:    | x_1[0]      |     0.1     |        0.25        |     0.4     | float |
    INFO - 08:24:31:    | x_1[1]      |     0.75    |         1          |     1.25    | float |
    INFO - 08:24:31:    | x_2         |     0.75    |         1          |     1.25    | float |
    INFO - 08:24:31:    | x_3         |     0.1     |        0.5         |      1      | float |
    INFO - 08:24:31:    | y_14[0]     |    24850    |   50606.9741711    |    77100    | float |
    INFO - 08:24:31:    | y_14[1]     |    -7700    |   7306.20262124    |    45000    | float |
    INFO - 08:24:31:    | y_32        |    0.235    | 0.5027962499999999 |    0.795    | float |
    INFO - 08:24:31:    | y_31        |     2960    |   6354.32430691    |    10185    | float |
    INFO - 08:24:31:    | y_24        |     0.44    |     4.15006276     |    11.13    | float |
    INFO - 08:24:31:    | y_34        |     0.44    |     1.10754577     |     1.98    | float |
    INFO - 08:24:31:    | y_23        |     3365    |   12194.2671934    |    26400    | float |
    INFO - 08:24:31:    | y_21        |    24850    |   50606.9741711    |    77250    | float |
    INFO - 08:24:31:    | y_12[0]     |    24850    |     50606.9742     |    77250    | float |
    INFO - 08:24:31:    | y_12[1]     |     0.45    |        0.95        |     1.5     | float |
    INFO - 08:24:31:    +-------------+-------------+--------------------+-------------+-------+
    INFO - 08:24:31: Solving optimization problem with algorithm SLSQP:
    INFO - 08:24:31: ...   0%|          | 0/20 [00:00<?, ?it]
    INFO - 08:24:31: ...   5%|▌         | 1/20 [00:00<00:00, 252.26 it/sec, obj=-536]
    INFO - 08:24:31: ...  10%|█         | 2/20 [00:00<00:00, 57.90 it/sec, obj=-1.49e+3]
    INFO - 08:24:31: ...  15%|█▌        | 3/20 [00:00<00:00, 60.65 it/sec, obj=-3.83e+3]
    INFO - 08:24:31: ...  20%|██        | 4/20 [00:00<00:00, 62.22 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  25%|██▌       | 5/20 [00:00<00:00, 63.30 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  30%|███       | 6/20 [00:00<00:00, 64.20 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  35%|███▌      | 7/20 [00:00<00:00, 64.58 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  40%|████      | 8/20 [00:00<00:00, 64.98 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  45%|████▌     | 9/20 [00:00<00:00, 69.27 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  50%|█████     | 10/20 [00:00<00:00, 69.00 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  55%|█████▌    | 11/20 [00:00<00:00, 72.50 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  60%|██████    | 12/20 [00:00<00:00, 75.65 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  65%|██████▌   | 13/20 [00:00<00:00, 78.54 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  70%|███████   | 14/20 [00:00<00:00, 81.23 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  75%|███████▌  | 15/20 [00:00<00:00, 83.66 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  80%|████████  | 16/20 [00:00<00:00, 85.99 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  85%|████████▌ | 17/20 [00:00<00:00, 88.17 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  90%|█████████ | 18/20 [00:00<00:00, 90.18 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ...  95%|█████████▌| 19/20 [00:00<00:00, 92.06 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: ... 100%|██████████| 20/20 [00:00<00:00, 93.79 it/sec, obj=-3.96e+3]
    INFO - 08:24:31: Optimization result:
    INFO - 08:24:31:    Optimizer info:
    INFO - 08:24:31:       Status: 8
    INFO - 08:24:31:       Message: Positive directional derivative for linesearch
    INFO - 08:24:31:       Number of calls to the objective function by the optimizer: 21
    INFO - 08:24:31:    Solution:
    INFO - 08:24:31:       The solution is feasible.
    INFO - 08:24:31:       Objective: -3963.909125280606
    INFO - 08:24:31:       Standardized constraints:
    INFO - 08:24:31:          g_1 = [-0.01807265 -0.03335477 -0.04425595 -0.0518399  -0.05733055 -0.13720865
    INFO - 08:24:31:  -0.10279135]
    INFO - 08:24:31:          g_2 = 7.162693463458325e-06
    INFO - 08:24:31:          g_3 = [-7.67188159e-01 -2.32811841e-01 -3.73860557e-05 -1.83255000e-01]
    INFO - 08:24:31:          y_12_y_14 = [ 1.12943503e-05  4.28787041e-07  1.13267743e-05 -1.50426345e-05]
    INFO - 08:24:31:          y_21_y_23_y_24 = [-2.77708306e-16  1.33085464e-05 -2.62473913e-05]
    INFO - 08:24:31:          y_31_y_32_y_34 = [-2.35661889e-06 -2.32437808e-06  1.57269783e-05]
    INFO - 08:24:31:       Design space:
    INFO - 08:24:31:       +-------------+-------------+---------------------+-------------+-------+
    INFO - 08:24:31:       | name        | lower_bound |        value        | upper_bound | type  |
    INFO - 08:24:31:       +-------------+-------------+---------------------+-------------+-------+
    INFO - 08:24:31:       | x_shared[0] |     0.01    | 0.06000179067336589 |     0.09    | float |
    INFO - 08:24:31:       | x_shared[1] |    30000    |        60000        |    60000    | float |
    INFO - 08:24:31:       | x_shared[2] |     1.4     |         1.4         |     1.8     | float |
    INFO - 08:24:31:       | x_shared[3] |     2.5     |         2.5         |     8.5     | float |
    INFO - 08:24:31:       | x_shared[4] |      40     |          70         |      70     | float |
    INFO - 08:24:31:       | x_shared[5] |     500     |         1500        |     1500    | float |
    INFO - 08:24:31:       | x_1[0]      |     0.1     |         0.4         |     0.4     | float |
    INFO - 08:24:31:       | x_1[1]      |     0.75    |         0.75        |     1.25    | float |
    INFO - 08:24:31:       | x_2         |     0.75    |         0.75        |     1.25    | float |
    INFO - 08:24:31:       | x_3         |     0.1     |  0.156238904271528  |      1      | float |
    INFO - 08:24:31:       | y_14[0]     |    24850    |  44749.85202975167  |    77100    | float |
    INFO - 08:24:31:       | y_14[1]     |    -7700    |  19351.86291108692  |    45000    | float |
    INFO - 08:24:31:       | y_32        |    0.235    |  0.7328131425977524 |    0.795    | float |
    INFO - 08:24:31:       | y_31        |     2960    |  9437.362336215188  |    10185    | float |
    INFO - 08:24:31:       | y_24        |     0.44    |   8.05763183450019  |    11.13    | float |
    INFO - 08:24:31:       | y_34        |     0.44    |  0.9239115967364436 |     1.98    | float |
    INFO - 08:24:31:       | y_23        |     3365    |   5553.60943830111  |    26400    | float |
    INFO - 08:24:31:       | y_21        |    24850    |  44749.85202975168  |    77250    | float |
    INFO - 08:24:31:       | y_12[0]     |    24850    |  44749.85202975167  |    77250    | float |
    INFO - 08:24:31:       | y_12[1]     |     0.45    |  0.9027908995968825 |     1.5     | float |
    INFO - 08:24:31:       +-------------+-------------+---------------------+-------------+-------+
    INFO - 08:24:31: *** End MDOScenario execution (time: 0:00:00.253011) ***

{'max_iter': 20, 'algo_options': {'ftol_rel': 1e-10, 'ineq_tolerance': 0.001, 'eq_tolerance': 0.001, 'normalize_design_space': True}, 'algo': 'SLSQP'}

Save the optimization history

We can save the whole optimization problem and its history for further post processing:

scenario.save_optimization_history("idf_history.h5", file_format="hdf5")
INFO - 08:24:31: Export optimization problem to file: idf_history.h5

We can also save only calls to functions and design variables history:

scenario.save_optimization_history("idf_history.xml", file_format="ggobi")
INFO - 08:24:31: Export to ggobi for functions: ['-y_4', 'g_1', 'g_2', 'g_3', 'y_12_y_14', 'y_21_y_23_y_24', 'y_31_y_32_y_34']
INFO - 08:24:31: Export to ggobi file: idf_history.xml

Plot the optimization history view

scenario.post_process("OptHistoryView", save=True, show=True)
<gemseo.post.opt_history_view.OptHistoryView object at 0x7f0cbcc30100>

Plot the quadratic approximation of the objective

scenario.post_process("QuadApprox", function="-y_4", save=False, show=True)
  • Hessian matrix SR1 approximation of -y_4
  • plot sobieski idf example
<gemseo.post.quad_approx.QuadApprox object at 0x7f0cbb5265e0>

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

Gallery generated by Sphinx-Gallery