Note

Go to the end to download the full example code

BiLevel-based MDO on the Sobieski SSBJ test case

from __future__ import annotations

from copy import deepcopy
from logging import WARNING

from gemseo import configure_logger
from gemseo import create_discipline
from gemseo import create_scenario
from gemseo import execute_post
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.

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 = SobieskiDesignSpace()

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.filter("x_3", copy=True),
    name="PropulsionScenario",
)
sc_prop.default_inputs = sub_sc_opts
sc_prop.add_constraint("g_3", constraint_type="ineq")

Build a sub-scenario for Aerodynamics

This sub-scenario will minimize L/D.

sc_aero = create_scenario(
    aero,
    "DisciplinaryOpt",
    "y_24",
    design_space.filter("x_2", copy=True),
    name="AerodynamicsScenario",
    maximize_objective=True,
)
sc_aero.default_inputs = sub_sc_opts
sc_aero.add_constraint("g_2", constraint_type="ineq")

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",
    design_space.filter("x_1", copy=True),
    name="StructureScenario",
    maximize_objective=True,
)
sc_str.add_constraint("g_1", constraint_type="ineq")
sc_str.default_inputs = sub_sc_opts

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

system_scenario = create_scenario(
    [sc_prop, sc_aero, sc_str, mission],
    "BiLevel",
    "y_4",
    design_space.filter("x_shared", copy=True),
    apply_cstr_tosub_scenarios=False,
    parallel_scenarios=False,
    multithread_scenarios=True,
    tolerance=1e-14,
    max_mda_iter=30,
    maximize_objective=True,
    sub_scenarios_log_level=WARNING,
)
system_scenario.add_constraint(["g_1", "g_2", "g_3"], constraint_type="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.

system_scenario.xdsmize(save_html=False)


Execute the main scenario

system_scenario.execute({
    "max_iter": 50,
    "algo": "NLOPT_COBYLA",
    "algo_options": algo_options,
})

    INFO - 13:54:07:
    INFO - 13:54:07: *** Start MDOScenario execution ***
    INFO - 13:54:07: MDOScenario
    INFO - 13:54:07:    Disciplines: AerodynamicsScenario PropulsionScenario SobieskiMission StructureScenario
    INFO - 13:54:07:    MDO formulation: BiLevel
    INFO - 13:54:07: Optimization problem:
    INFO - 13:54:07:    minimize -y_4(x_shared)
    INFO - 13:54:07:    with respect to x_shared
    INFO - 13:54:07:    subject to constraints:
    INFO - 13:54:07:       g_1_g_2_g_3(x_shared) <= 0.0
    INFO - 13:54:07:    over the design space:
    INFO - 13:54:07:       +-------------+-------------+-------+-------------+-------+
    INFO - 13:54:07:       | Name        | Lower bound | Value | Upper bound | Type  |
    INFO - 13:54:07:       +-------------+-------------+-------+-------------+-------+
    INFO - 13:54:07:       | x_shared[0] |     0.01    |  0.05 |     0.09    | float |
    INFO - 13:54:07:       | x_shared[1] |    30000    | 45000 |    60000    | float |
    INFO - 13:54:07:       | x_shared[2] |     1.4     |  1.6  |     1.8     | float |
    INFO - 13:54:07:       | x_shared[3] |     2.5     |  5.5  |     8.5     | float |
    INFO - 13:54:07:       | x_shared[4] |      40     |   55  |      70     | float |
    INFO - 13:54:07:       | x_shared[5] |     500     |  1000 |     1500    | float |
    INFO - 13:54:07:       +-------------+-------------+-------+-------------+-------+
    INFO - 13:54:07: Solving optimization problem with algorithm NLOPT_COBYLA:
    INFO - 13:54:08:      2%|▏         | 1/50 [00:00<00:17,  2.80 it/sec, obj=-553]
 WARNING - 13:54:08: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:08:       The solution is not feasible.
    INFO - 13:54:08:      4%|▍         | 2/50 [00:00<00:15,  3.18 it/sec, obj=-574]
 WARNING - 13:54:08: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:08:       The solution is not feasible.
    INFO - 13:54:08:      6%|▌         | 3/50 [00:00<00:13,  3.42 it/sec, obj=-813]
 WARNING - 13:54:08: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:08:       The solution is not feasible.
    INFO - 13:54:08:      8%|▊         | 4/50 [00:01<00:13,  3.52 it/sec, obj=-751]
 WARNING - 13:54:09: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:09:       The solution is not feasible.
    INFO - 13:54:09:     10%|█         | 5/50 [00:01<00:13,  3.33 it/sec, obj=-734]
 WARNING - 13:54:09: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:09:       The solution is not feasible.
    INFO - 13:54:09:     12%|█▏        | 6/50 [00:01<00:12,  3.45 it/sec, obj=-977]
 WARNING - 13:54:09: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:09:       The solution is not feasible.
    INFO - 13:54:09:     14%|█▍        | 7/50 [00:02<00:12,  3.43 it/sec, obj=-1.05e+3]
 WARNING - 13:54:10: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:10:       The solution is not feasible.
    INFO - 13:54:10:     16%|█▌        | 8/50 [00:02<00:12,  3.41 it/sec, obj=-1.67e+3]
    INFO - 13:54:10:     18%|█▊        | 9/50 [00:02<00:12,  3.41 it/sec, obj=-1.73e+3]
 WARNING - 13:54:10: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:10:       The solution is not feasible.
    INFO - 13:54:10:     20%|██        | 10/50 [00:02<00:11,  3.36 it/sec, obj=-2.59e+3]
    INFO - 13:54:11:     22%|██▏       | 11/50 [00:03<00:11,  3.38 it/sec, obj=-2.94e+3]
 WARNING - 13:54:11: MDAJacobi has reached its maximum number of iterations but the normed residual 2.189447656198491e-13 is still above the tolerance 1e-14.
    INFO - 13:54:11:     24%|██▍       | 12/50 [00:03<00:11,  3.27 it/sec, obj=-2.64e+3]
 WARNING - 13:54:11: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:11:       The solution is not feasible.
    INFO - 13:54:11:     26%|██▌       | 13/50 [00:03<00:11,  3.31 it/sec, obj=-2.85e+3]
 WARNING - 13:54:11: MDAJacobi has reached its maximum number of iterations but the normed residual 1.8505396985126146e-13 is still above the tolerance 1e-14.
    INFO - 13:54:12:     28%|██▊       | 14/50 [00:04<00:10,  3.29 it/sec, obj=-2.79e+3]
 WARNING - 13:54:12: MDAJacobi has reached its maximum number of iterations but the normed residual 3.0145345591022575e-14 is still above the tolerance 1e-14.
 WARNING - 13:54:12: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:12:       The solution is not feasible.
    INFO - 13:54:12:     30%|███       | 15/50 [00:04<00:10,  3.21 it/sec, obj=-2.4e+3]
    INFO - 13:54:12:     32%|███▏      | 16/50 [00:04<00:10,  3.26 it/sec, obj=-3.07e+3]
 WARNING - 13:54:12: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:12:       The solution is not feasible.
    INFO - 13:54:12:     34%|███▍      | 17/50 [00:05<00:09,  3.30 it/sec, obj=-3.01e+3]
    INFO - 13:54:13:     36%|███▌      | 18/50 [00:05<00:09,  3.34 it/sec, obj=-3.39e+3]
    INFO - 13:54:13:     38%|███▊      | 19/50 [00:05<00:09,  3.38 it/sec, obj=-3.84e+3]
 WARNING - 13:54:13: MDAJacobi has reached its maximum number of iterations but the normed residual 3.8877123894689336e-14 is still above the tolerance 1e-14.
    INFO - 13:54:13:     40%|████      | 20/50 [00:05<00:08,  3.39 it/sec, obj=-3.58e+3]
    INFO - 13:54:13:     42%|████▏     | 21/50 [00:06<00:08,  3.42 it/sec, obj=-3.66e+3]
    INFO - 13:54:14:     44%|████▍     | 22/50 [00:06<00:08,  3.45 it/sec, obj=-3.77e+3]
 WARNING - 13:54:14: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:14:       The solution is not feasible.
    INFO - 13:54:14:     46%|████▌     | 23/50 [00:06<00:07,  3.48 it/sec, obj=-3.75e+3]
    INFO - 13:54:14:     48%|████▊     | 24/50 [00:06<00:07,  3.49 it/sec, obj=-2.94e+3]
 WARNING - 13:54:14: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:14:       The solution is not feasible.
    INFO - 13:54:15:     50%|█████     | 25/50 [00:07<00:07,  3.50 it/sec, obj=-3.71e+3]
    INFO - 13:54:15:     52%|█████▏    | 26/50 [00:07<00:06,  3.56 it/sec, obj=-3.96e+3]
    INFO - 13:54:15:     54%|█████▍    | 27/50 [00:07<00:06,  3.61 it/sec, obj=-3.87e+3]
 WARNING - 13:54:15: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:15:       The solution is not feasible.
    INFO - 13:54:15:     56%|█████▌    | 28/50 [00:07<00:06,  3.65 it/sec, obj=-3.95e+3]
    INFO - 13:54:15:     58%|█████▊    | 29/50 [00:07<00:05,  3.70 it/sec, obj=-3.93e+3]
 WARNING - 13:54:15: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:15:       The solution is not feasible.
    INFO - 13:54:15:     60%|██████    | 30/50 [00:08<00:05,  3.74 it/sec, obj=-3.97e+3]
    INFO - 13:54:16:     62%|██████▏   | 31/50 [00:08<00:05,  3.77 it/sec, obj=-3.96e+3]
    INFO - 13:54:16:     64%|██████▍   | 32/50 [00:08<00:04,  3.81 it/sec, obj=-3.96e+3]
    INFO - 13:54:16:     66%|██████▌   | 33/50 [00:08<00:04,  3.85 it/sec, obj=-3.95e+3]
    INFO - 13:54:16:     68%|██████▊   | 34/50 [00:08<00:04,  3.91 it/sec, obj=-3.94e+3]
 WARNING - 13:54:16: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:16:       The solution is not feasible.
    INFO - 13:54:16:     70%|███████   | 35/50 [00:08<00:03,  3.95 it/sec, obj=-3.97e+3]
 WARNING - 13:54:16: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:16:       The solution is not feasible.
    INFO - 13:54:16:     72%|███████▏  | 36/50 [00:09<00:03,  3.97 it/sec, obj=-3.98e+3]
    INFO - 13:54:17:     74%|███████▍  | 37/50 [00:09<00:03,  4.00 it/sec, obj=-3.96e+3]
    INFO - 13:54:17:     76%|███████▌  | 38/50 [00:09<00:02,  4.06 it/sec, obj=-3.95e+3]
    INFO - 13:54:17:     78%|███████▊  | 39/50 [00:09<00:02,  4.09 it/sec, obj=-3.95e+3]
    INFO - 13:54:17:     80%|████████  | 40/50 [00:09<00:02,  4.13 it/sec, obj=-3.96e+3]
    INFO - 13:54:17:     82%|████████▏ | 41/50 [00:09<00:02,  4.17 it/sec, obj=-3.96e+3]
    INFO - 13:54:17:     84%|████████▍ | 42/50 [00:09<00:01,  4.21 it/sec, obj=-3.96e+3]
    INFO - 13:54:18:     86%|████████▌ | 43/50 [00:10<00:01,  4.23 it/sec, obj=-3.95e+3]
 WARNING - 13:54:18: MDAJacobi has reached its maximum number of iterations but the normed residual 3.8001101574110354e-14 is still above the tolerance 1e-14.
    INFO - 13:54:18:     88%|████████▊ | 44/50 [00:10<00:01,  4.21 it/sec, obj=-3.96e+3]
    INFO - 13:54:18:     90%|█████████ | 45/50 [00:10<00:01,  4.24 it/sec, obj=-3.96e+3]
    INFO - 13:54:18:     92%|█████████▏| 46/50 [00:10<00:00,  4.27 it/sec, obj=-3.96e+3]
    INFO - 13:54:18:     94%|█████████▍| 47/50 [00:10<00:00,  4.30 it/sec, obj=-3.96e+3]
 WARNING - 13:54:18: Optimization found no feasible point !  The least infeasible point is selected.
 WARNING - 13:54:18:       The solution is not feasible.
    INFO - 13:54:18:     96%|█████████▌| 48/50 [00:11<00:00,  4.34 it/sec, obj=-3.96e+3]
    INFO - 13:54:19:     98%|█████████▊| 49/50 [00:11<00:00,  4.37 it/sec, obj=-3.96e+3]
    INFO - 13:54:19:    100%|██████████| 50/50 [00:11<00:00,  4.41 it/sec, obj=-3.96e+3]
    INFO - 13:54:19: Optimization result:
    INFO - 13:54:19:    Optimizer info:
    INFO - 13:54:19:       Status: None
    INFO - 13:54:19:       Message: Maximum number of iterations reached. GEMSEO Stopped the driver
    INFO - 13:54:19:       Number of calls to the objective function by the optimizer: 52
    INFO - 13:54:19:    Solution:
    INFO - 13:54:19:       The solution is feasible.
    INFO - 13:54:19:       Objective: -3963.3800122701787
    INFO - 13:54:19:       Standardized constraints:
    INFO - 13:54:19:          g_1_g_2_g_3 = [-0.01805093 -0.03333915 -0.04424381 -0.05182998 -0.05732217 -0.13720865
    INFO - 13:54:19:  -0.10279135  0.         -0.76718646 -0.23281354  0.         -0.183255  ]
    INFO - 13:54:19:       Design space:
    INFO - 13:54:19:          +-------------+-------------+---------------------+-------------+-------+
    INFO - 13:54:19:          | Name        | Lower bound |        Value        | Upper bound | Type  |
    INFO - 13:54:19:          +-------------+-------------+---------------------+-------------+-------+
    INFO - 13:54:19:          | x_shared[0] |     0.01    | 0.05999999999999999 |     0.09    | float |
    INFO - 13:54:19:          | x_shared[1] |    30000    |        60000        |    60000    | float |
    INFO - 13:54:19:          | x_shared[2] |     1.4     |         1.4         |     1.8     | float |
    INFO - 13:54:19:          | x_shared[3] |     2.5     |         2.5         |     8.5     | float |
    INFO - 13:54:19:          | x_shared[4] |      40     |          70         |      70     | float |
    INFO - 13:54:19:          | x_shared[5] |     500     |         1500        |     1500    | float |
    INFO - 13:54:19:          +-------------+-------------+---------------------+-------------+-------+
    INFO - 13:54:19: *** End MDOScenario execution (time: 0:00:11.361317) ***

{'max_iter': 50, 'algo_options': {'xtol_rel': 1e-07, 'xtol_abs': 1e-07, 'ftol_rel': 1e-07, 'ftol_abs': 1e-07, 'ineq_tolerance': 0.0001}, 'algo': 'NLOPT_COBYLA'}

Plot the history of the MDA residuals

For the first MDA:

system_scenario.formulation.mda1.plot_residual_history(save=False, show=True)
# For the second MDA:
system_scenario.formulation.mda2.plot_residual_history(save=False, show=True)
  • MDAJacobi: residual plot
  • MDAJacobi: residual plot

Plot the system optimization history view

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

Plot the structure optimization histories of the 2 first iterations

struct_databases = system_scenario.formulation.scenario_adapters[2].databases
for database in struct_databases[:2]:
    opt_problem = deepcopy(sc_str.formulation.opt_problem)
    opt_problem.database = database
    execute_post(opt_problem, "OptHistoryView", save=False, show=True)

for disc in [propu, aero, mission, struct]:
    print(f"{disc.name}: {disc.n_calls} calls.")
  • Evolution of the optimization variables
  • Evolution of the objective value
  • Distance to the optimum
  • Hessian diagonal approximation
  • Evolution of the inequality constraints
  • Evolution of the optimization variables
  • Evolution of the objective value
  • Distance to the optimum
  • Hessian diagonal approximation
  • Evolution of the inequality constraints
SobieskiPropulsion: 1290 calls.
SobieskiAerodynamics: 1343 calls.
SobieskiMission: 50 calls.
SobieskiStructure: 1394 calls.

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

Gallery generated by Sphinx-Gallery