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.algos.opt.scipy_local.settings.slsqp import SLSQP_Settings
from gemseo.problems.mdo.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 Discipline.default_input_data

for discipline in disciplines:
    print(discipline)
    print(f"Default inputs: {discipline.default_input_data}")
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 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 = 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


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

Set the design constraints#

for c_name in ["g_1", "g_2", "g_3"]:
    scenario.add_constraint(c_name, constraint_type="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, pdf_build=False)


Define the algorithm inputs#

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

slsqp_settings = SLSQP_Settings(
    max_iter=20,
    ftol_rel=1e-10,
    ineq_tolerance=1e-3,
    eq_tolerance=1e-3,
    normalize_design_space=True,
)

Execute the scenario#

scenario.execute(slsqp_settings)
INFO - 15:40:13: *** Start MDOScenario execution ***
INFO - 15:40:13: MDOScenario
INFO - 15:40:13:    Disciplines: SobieskiAerodynamics SobieskiMission SobieskiPropulsion SobieskiStructure
INFO - 15:40:13:    MDO formulation: IDF
INFO - 15:40:13: Optimization problem:
INFO - 15:40:13:    minimize -y_4(x_shared, y_14, y_24, y_34)
INFO - 15:40:13:    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 - 15:40:13:    subject to constraints:
INFO - 15:40:13:       g_1(x_shared, x_1, y_31, y_21) <= 0
INFO - 15:40:13:       g_2(x_shared, x_2, y_32, y_12) <= 0
INFO - 15:40:13:       g_3(x_shared, x_3, y_23) <= 0
INFO - 15:40:13:       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 - 15:40:13:                                            y_32(x_shared, x_3, y_23) - y_32 == 0.0
INFO - 15:40:13:                                            y_34(x_shared, x_3, y_23) - y_34 == 0.0
INFO - 15:40:13:       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 - 15:40:13:                                                  y_23(x_shared, x_2, y_32, y_12) - y_23 == 0.0
INFO - 15:40:13:                                                  y_24(x_shared, x_2, y_32, y_12) - y_24 == 0.0
INFO - 15:40:13:       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 - 15:40:13:                                             y_14(x_shared, x_1, y_31, y_21) - y_14 == 0.0
INFO - 15:40:13:    over the design space:
INFO - 15:40:13:       +-------------+-------------+--------------------+-------------+-------+
INFO - 15:40:13:       | Name        | Lower bound |       Value        | Upper bound | Type  |
INFO - 15:40:13:       +-------------+-------------+--------------------+-------------+-------+
INFO - 15:40:13:       | x_shared[0] |     0.01    |        0.05        |     0.09    | float |
INFO - 15:40:13:       | x_shared[1] |    30000    |       45000        |    60000    | float |
INFO - 15:40:13:       | x_shared[2] |     1.4     |        1.6         |     1.8     | float |
INFO - 15:40:13:       | x_shared[3] |     2.5     |        5.5         |     8.5     | float |
INFO - 15:40:13:       | x_shared[4] |      40     |         55         |      70     | float |
INFO - 15:40:13:       | x_shared[5] |     500     |        1000        |     1500    | float |
INFO - 15:40:13:       | x_1[0]      |     0.1     |        0.25        |     0.4     | float |
INFO - 15:40:13:       | x_1[1]      |     0.75    |         1          |     1.25    | float |
INFO - 15:40:13:       | x_2         |     0.75    |         1          |     1.25    | float |
INFO - 15:40:13:       | x_3         |     0.1     |        0.5         |      1      | float |
INFO - 15:40:13:       | y_14[0]     |    24850    |   50606.9741711    |    77100    | float |
INFO - 15:40:13:       | y_14[1]     |    -7700    |   7306.20262124    |    45000    | float |
INFO - 15:40:13:       | y_32        |    0.235    | 0.5027962499999999 |    0.795    | float |
INFO - 15:40:13:       | y_31        |     2960    |   6354.32430691    |    10185    | float |
INFO - 15:40:13:       | y_24        |     0.44    |     4.15006276     |    11.13    | float |
INFO - 15:40:13:       | y_34        |     0.44    |     1.10754577     |     1.98    | float |
INFO - 15:40:13:       | y_23        |     3365    |   12194.2671934    |    26400    | float |
INFO - 15:40:13:       | y_21        |    24850    |   50606.9741711    |    77250    | float |
INFO - 15:40:13:       | y_12[0]     |    24850    |     50606.9742     |    77250    | float |
INFO - 15:40:13:       | y_12[1]     |     0.45    |        0.95        |     1.5     | float |
INFO - 15:40:13:       +-------------+-------------+--------------------+-------------+-------+
INFO - 15:40:13: Solving optimization problem with algorithm SLSQP:
INFO - 15:40:13:      5%|▌         | 1/20 [00:00<00:00, 322.47 it/sec, obj=-536]
INFO - 15:40:13:     10%|█         | 2/20 [00:00<00:00, 81.31 it/sec, obj=-1.49e+3]
INFO - 15:40:13:     15%|█▌        | 3/20 [00:00<00:00, 91.43 it/sec, obj=-3.83e+3]
INFO - 15:40:13:     20%|██        | 4/20 [00:00<00:00, 97.52 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     25%|██▌       | 5/20 [00:00<00:00, 101.97 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     30%|███       | 6/20 [00:00<00:00, 105.42 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     35%|███▌      | 7/20 [00:00<00:00, 107.13 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     40%|████      | 8/20 [00:00<00:00, 108.21 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     45%|████▌     | 9/20 [00:00<00:00, 115.47 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     50%|█████     | 10/20 [00:00<00:00, 115.79 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     55%|█████▌    | 11/20 [00:00<00:00, 122.09 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     60%|██████    | 12/20 [00:00<00:00, 127.96 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     65%|██████▌   | 13/20 [00:00<00:00, 133.25 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     70%|███████   | 14/20 [00:00<00:00, 138.27 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     75%|███████▌  | 15/20 [00:00<00:00, 142.87 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     80%|████████  | 16/20 [00:00<00:00, 147.22 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     85%|████████▌ | 17/20 [00:00<00:00, 151.24 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     90%|█████████ | 18/20 [00:00<00:00, 154.67 it/sec, obj=-3.96e+3]
INFO - 15:40:13:     95%|█████████▌| 19/20 [00:00<00:00, 158.18 it/sec, obj=-3.96e+3]
INFO - 15:40:13:    100%|██████████| 20/20 [00:00<00:00, 161.57 it/sec, obj=-3.96e+3]
INFO - 15:40:13: Optimization result:
INFO - 15:40:13:    Optimizer info:
INFO - 15:40:13:       Status: 8
INFO - 15:40:13:       Message: Positive directional derivative for linesearch
INFO - 15:40:13:       Number of calls to the objective function by the optimizer: 21
INFO - 15:40:13:    Solution:
INFO - 15:40:13:       The solution is feasible.
INFO - 15:40:13:       Objective: -3963.909125280606
INFO - 15:40:13:       Standardized constraints:
INFO - 15:40:13:          g_1 = [-0.01807265 -0.03335477 -0.04425595 -0.0518399  -0.05733055 -0.13720865
INFO - 15:40:13:  -0.10279135]
INFO - 15:40:13:          g_2 = 7.162693463458325e-06
INFO - 15:40:13:          g_3 = [-7.67188159e-01 -2.32811841e-01 -3.73860557e-05 -1.83255000e-01]
INFO - 15:40:13:          y_12_y_14 = [ 1.12943503e-05  4.28787041e-07  1.13267743e-05 -1.50426345e-05]
INFO - 15:40:13:          y_21_y_23_y_24 = [-2.77708306e-16  1.33085464e-05 -2.62473913e-05]
INFO - 15:40:13:          y_31_y_32_y_34 = [-2.35661889e-06 -2.32437808e-06  1.57269783e-05]
INFO - 15:40:13:       Design space:
INFO - 15:40:13:          +-------------+-------------+---------------------+-------------+-------+
INFO - 15:40:13:          | Name        | Lower bound |        Value        | Upper bound | Type  |
INFO - 15:40:13:          +-------------+-------------+---------------------+-------------+-------+
INFO - 15:40:13:          | x_shared[0] |     0.01    | 0.06000179067336589 |     0.09    | float |
INFO - 15:40:13:          | x_shared[1] |    30000    |        60000        |    60000    | float |
INFO - 15:40:13:          | x_shared[2] |     1.4     |         1.4         |     1.8     | float |
INFO - 15:40:13:          | x_shared[3] |     2.5     |         2.5         |     8.5     | float |
INFO - 15:40:13:          | x_shared[4] |      40     |          70         |      70     | float |
INFO - 15:40:13:          | x_shared[5] |     500     |         1500        |     1500    | float |
INFO - 15:40:13:          | x_1[0]      |     0.1     |         0.4         |     0.4     | float |
INFO - 15:40:13:          | x_1[1]      |     0.75    |         0.75        |     1.25    | float |
INFO - 15:40:13:          | x_2         |     0.75    |         0.75        |     1.25    | float |
INFO - 15:40:13:          | x_3         |     0.1     |  0.156238904271528  |      1      | float |
INFO - 15:40:13:          | y_14[0]     |    24850    |  44749.85202975167  |    77100    | float |
INFO - 15:40:13:          | y_14[1]     |    -7700    |  19351.86291108692  |    45000    | float |
INFO - 15:40:13:          | y_32        |    0.235    |  0.7328131425977524 |    0.795    | float |
INFO - 15:40:13:          | y_31        |     2960    |  9437.362336215188  |    10185    | float |
INFO - 15:40:13:          | y_24        |     0.44    |   8.05763183450019  |    11.13    | float |
INFO - 15:40:13:          | y_34        |     0.44    |  0.9239115967364436 |     1.98    | float |
INFO - 15:40:13:          | y_23        |     3365    |   5553.60943830111  |    26400    | float |
INFO - 15:40:13:          | y_21        |    24850    |  44749.85202975168  |    77250    | float |
INFO - 15:40:13:          | y_12[0]     |    24850    |  44749.85202975167  |    77250    | float |
INFO - 15:40:13:          | y_12[1]     |     0.45    |  0.9027908995968825 |     1.5     | float |
INFO - 15:40:13:          +-------------+-------------+---------------------+-------------+-------+
INFO - 15:40:13: *** End MDOScenario execution (time: 0:00:00.138194) ***

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 - 15:40:13: Exporting the optimization problem to the 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")

Plot the optimization history view#

scenario.post_process(post_name="OptHistoryView", save=False, show=True)
  • Evolution of the optimization variables
  • Evolution of the objective value
  • Evolution of the distance to the optimum
  • Evolution of the inequality constraints
  • Evolution of the equality constraints
<gemseo.post.opt_history_view.OptHistoryView object at 0x7f51f905f470>

Plot the quadratic approximation of the objective#

scenario.post_process(post_name="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 0x7f51f72a7da0>

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

Gallery generated by Sphinx-Gallery