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.algos.opt.nlopt.settings.nlopt_cobyla_settings import NLOPT_COBYLA_Settings
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.
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 settings:
slsqp_settings = SLSQP_Settings(
max_iter=30,
xtol_rel=1e-7,
xtol_abs=1e-7,
ftol_rel=1e-7,
ftol_abs=1e-7,
ineq_tolerance=1e-4,
)
cobyla_settings = NLOPT_COBYLA_Settings(
max_iter=50,
xtol_rel=1e-7,
xtol_abs=1e-7,
ftol_rel=1e-7,
ftol_abs=1e-7,
ineq_tolerance=1e-4,
)
Build a sub-scenario for Propulsion#
This sub-scenario will minimize SFC.
sc_prop = create_scenario(
propu,
"y_34",
design_space.filter("x_3", copy=True),
name="PropulsionScenario",
formulation_name="DisciplinaryOpt",
)
sc_prop.set_algorithm(slsqp_settings)
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,
"y_24",
design_space.filter("x_2", copy=True),
name="AerodynamicsScenario",
maximize_objective=True,
formulation_name="DisciplinaryOpt",
)
sc_aero.set_algorithm(slsqp_settings)
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,
"y_11",
design_space.filter("x_1", copy=True),
name="StructureScenario",
maximize_objective=True,
formulation_name="DisciplinaryOpt",
)
sc_str.add_constraint("g_1", constraint_type="ineq")
sc_str.set_algorithm(slsqp_settings)
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],
"y_4",
design_space.filter("x_shared", copy=True),
apply_cstr_tosub_scenarios=False,
parallel_scenarios=False,
multithread_scenarios=True,
main_mda_settings={"tolerance": 1e-14, "max_mda_iter": 30},
maximize_objective=True,
sub_scenarios_log_level=WARNING,
formulation_name="BiLevel",
)
system_scenario.add_constraint(["g_1", "g_2", "g_3"], constraint_type="ineq")
WARNING - 08:36:38: Unsupported feature 'minItems' in JSONGrammar 'SobieskiMission_discipline_output' for property 'y_4' in conversion to SimpleGrammar.
WARNING - 08:36:38: Unsupported feature 'maxItems' in JSONGrammar 'SobieskiMission_discipline_output' for property 'y_4' in conversion to SimpleGrammar.
Visualize the XDSM#
Generate the XDSM on the fly:
log_workflow_status=Truewill log the status of the workflow in the console,save_html(defaultTrue) will generate a self-contained HTML file, that can be automatically opened usingshow_html=True.
system_scenario.xdsmize(save_html=False, pdf_build=False)
Execute the main scenario#
system_scenario.execute(cobyla_settings)
INFO - 08:36:38:
INFO - 08:36:38: *** Start MDOScenario execution ***
INFO - 08:36:38: MDOScenario
INFO - 08:36:38: Disciplines: AerodynamicsScenario PropulsionScenario SobieskiMission StructureScenario
INFO - 08:36:38: MDO formulation: BiLevel
INFO - 08:36:38: Optimization problem:
INFO - 08:36:38: minimize -y_4(x_shared)
INFO - 08:36:38: with respect to x_shared
INFO - 08:36:38: subject to constraints:
INFO - 08:36:38: g_1_g_2_g_3(x_shared) <= 0
INFO - 08:36:38: over the design space:
INFO - 08:36:38: +-------------+-------------+-------+-------------+-------+
INFO - 08:36:38: | Name | Lower bound | Value | Upper bound | Type |
INFO - 08:36:38: +-------------+-------------+-------+-------------+-------+
INFO - 08:36:38: | x_shared[0] | 0.01 | 0.05 | 0.09 | float |
INFO - 08:36:38: | x_shared[1] | 30000 | 45000 | 60000 | float |
INFO - 08:36:38: | x_shared[2] | 1.4 | 1.6 | 1.8 | float |
INFO - 08:36:38: | x_shared[3] | 2.5 | 5.5 | 8.5 | float |
INFO - 08:36:38: | x_shared[4] | 40 | 55 | 70 | float |
INFO - 08:36:38: | x_shared[5] | 500 | 1000 | 1500 | float |
INFO - 08:36:38: +-------------+-------------+-------+-------------+-------+
INFO - 08:36:38: Solving optimization problem with algorithm NLOPT_COBYLA:
INFO - 08:36:38: 2%|▏ | 1/50 [00:00<00:09, 5.03 it/sec, obj=-553]
WARNING - 08:36:38: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:38: The solution is not feasible.
INFO - 08:36:38: 4%|▍ | 2/50 [00:00<00:09, 5.01 it/sec, obj=-574]
WARNING - 08:36:38: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:38: The solution is not feasible.
INFO - 08:36:38: 6%|▌ | 3/50 [00:00<00:09, 5.05 it/sec, obj=-813]
WARNING - 08:36:38: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:38: The solution is not feasible.
INFO - 08:36:39: 8%|▊ | 4/50 [00:00<00:09, 4.92 it/sec, obj=-751]
WARNING - 08:36:39: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:39: The solution is not feasible.
INFO - 08:36:39: 10%|█ | 5/50 [00:01<00:10, 4.48 it/sec, obj=-734]
WARNING - 08:36:39: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:39: The solution is not feasible.
INFO - 08:36:39: 12%|█▏ | 6/50 [00:01<00:09, 4.61 it/sec, obj=-977]
WARNING - 08:36:39: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:39: The solution is not feasible.
INFO - 08:36:39: 14%|█▍ | 7/50 [00:01<00:09, 4.46 it/sec, obj=-1.05e+3]
WARNING - 08:36:39: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:39: The solution is not feasible.
INFO - 08:36:40: 16%|█▌ | 8/50 [00:01<00:09, 4.43 it/sec, obj=-1.67e+3]
INFO - 08:36:40: 18%|█▊ | 9/50 [00:02<00:09, 4.42 it/sec, obj=-1.73e+3]
WARNING - 08:36:40: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:40: The solution is not feasible.
INFO - 08:36:40: 20%|██ | 10/50 [00:02<00:09, 4.36 it/sec, obj=-2.59e+3]
INFO - 08:36:40: 22%|██▏ | 11/50 [00:02<00:08, 4.36 it/sec, obj=-2.94e+3]
WARNING - 08:36:40: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:40: The solution is not feasible.
INFO - 08:36:41: 24%|██▍ | 12/50 [00:02<00:08, 4.28 it/sec, obj=-2.39e+3]
INFO - 08:36:41: 26%|██▌ | 13/50 [00:03<00:08, 4.28 it/sec, obj=-2.66e+3]
WARNING - 08:36:41: MDAJacobi has reached its maximum number of iterations but the normed residual 4.00581977874369e-13 is still above the tolerance 1e-14.
INFO - 08:36:41: 28%|██▊ | 14/50 [00:03<00:08, 4.26 it/sec, obj=-2.54e+3]
WARNING - 08:36:41: MDAJacobi has reached its maximum number of iterations but the normed residual 2.6870836615612218e-14 is still above the tolerance 1e-14.
INFO - 08:36:41: 30%|███ | 15/50 [00:03<00:08, 4.19 it/sec, obj=-2.51e+3]
WARNING - 08:36:42: MDAJacobi has reached its maximum number of iterations but the normed residual 2.6870836615612218e-14 is still above the tolerance 1e-14.
INFO - 08:36:42: 32%|███▏ | 16/50 [00:03<00:08, 4.14 it/sec, obj=-2.2e+3]
INFO - 08:36:42: 34%|███▍ | 17/50 [00:04<00:07, 4.16 it/sec, obj=-2.76e+3]
INFO - 08:36:42: 36%|███▌ | 18/50 [00:04<00:07, 4.13 it/sec, obj=-2.47e+3]
WARNING - 08:36:42: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:42: The solution is not feasible.
INFO - 08:36:42: 38%|███▊ | 19/50 [00:04<00:07, 4.13 it/sec, obj=-2.97e+3]
INFO - 08:36:43: 40%|████ | 20/50 [00:04<00:07, 4.16 it/sec, obj=-2.51e+3]
INFO - 08:36:43: 42%|████▏ | 21/50 [00:04<00:06, 4.21 it/sec, obj=-2.85e+3]
INFO - 08:36:43: 44%|████▍ | 22/50 [00:05<00:06, 4.21 it/sec, obj=-2.99e+3]
INFO - 08:36:43: 46%|████▌ | 23/50 [00:05<00:06, 4.24 it/sec, obj=-2.7e+3]
INFO - 08:36:43: 48%|████▊ | 24/50 [00:05<00:06, 4.25 it/sec, obj=-3.1e+3]
INFO - 08:36:44: 50%|█████ | 25/50 [00:05<00:05, 4.31 it/sec, obj=-2.8e+3]
INFO - 08:36:44: 52%|█████▏ | 26/50 [00:06<00:05, 4.28 it/sec, obj=-3.07e+3]
INFO - 08:36:44: 54%|█████▍ | 27/50 [00:06<00:05, 4.32 it/sec, obj=-3.15e+3]
WARNING - 08:36:44: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:44: The solution is not feasible.
INFO - 08:36:44: 56%|█████▌ | 28/50 [00:06<00:05, 4.28 it/sec, obj=-2.8e+3]
WARNING - 08:36:44: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:44: The solution is not feasible.
WARNING - 08:36:45: MDAJacobi has reached its maximum number of iterations but the normed residual 2.6870836615612218e-14 is still above the tolerance 1e-14.
INFO - 08:36:45: 58%|█████▊ | 29/50 [00:06<00:04, 4.27 it/sec, obj=-3.14e+3]
WARNING - 08:36:45: Optimization found no feasible point; the least infeasible point is selected.
WARNING - 08:36:45: The solution is not feasible.
INFO - 08:36:45: 60%|██████ | 30/50 [00:06<00:04, 4.32 it/sec, obj=-3.12e+3]
WARNING - 08:36:45: MDAJacobi has reached its maximum number of iterations but the normed residual 1.7359275648513364e-14 is still above the tolerance 1e-14.
INFO - 08:36:45: 62%|██████▏ | 31/50 [00:07<00:04, 4.34 it/sec, obj=-3.12e+3]
WARNING - 08:36:45: MDAJacobi has reached its maximum number of iterations but the normed residual 2.6870836615612218e-14 is still above the tolerance 1e-14.
INFO - 08:36:45: 64%|██████▍ | 32/50 [00:07<00:04, 4.36 it/sec, obj=-3.12e+3]
INFO - 08:36:45: 66%|██████▌ | 33/50 [00:07<00:03, 4.40 it/sec, obj=-3.22e+3]
INFO - 08:36:45: 68%|██████▊ | 34/50 [00:07<00:03, 4.43 it/sec, obj=-3.12e+3]
INFO - 08:36:46: 70%|███████ | 35/50 [00:07<00:03, 4.45 it/sec, obj=-3.13e+3]
INFO - 08:36:46: 72%|███████▏ | 36/50 [00:08<00:03, 4.49 it/sec, obj=-3.11e+3]
INFO - 08:36:46: 74%|███████▍ | 37/50 [00:08<00:02, 4.52 it/sec, obj=-3.14e+3]
INFO - 08:36:46: 76%|███████▌ | 38/50 [00:08<00:02, 4.55 it/sec, obj=-3.16e+3]
INFO - 08:36:46: 78%|███████▊ | 39/50 [00:08<00:02, 4.58 it/sec, obj=-3.2e+3]
WARNING - 08:36:46: MDAJacobi has reached its maximum number of iterations but the normed residual 1.7359275648513364e-14 is still above the tolerance 1e-14.
INFO - 08:36:46: 80%|████████ | 40/50 [00:08<00:02, 4.59 it/sec, obj=-3.16e+3]
INFO - 08:36:47: 82%|████████▏ | 41/50 [00:08<00:01, 4.62 it/sec, obj=-3.16e+3]
INFO - 08:36:47: 84%|████████▍ | 42/50 [00:09<00:01, 4.65 it/sec, obj=-3.15e+3]
INFO - 08:36:47: 86%|████████▌ | 43/50 [00:09<00:01, 4.68 it/sec, obj=-3.17e+3]
INFO - 08:36:47: 88%|████████▊ | 44/50 [00:09<00:01, 4.70 it/sec, obj=-3.18e+3]
INFO - 08:36:47: 90%|█████████ | 45/50 [00:09<00:01, 4.73 it/sec, obj=-3.19e+3]
WARNING - 08:36:47: MDAJacobi has reached its maximum number of iterations but the normed residual 1.7359275648513364e-14 is still above the tolerance 1e-14.
INFO - 08:36:47: 92%|█████████▏| 46/50 [00:09<00:00, 4.73 it/sec, obj=-3.2e+3]
INFO - 08:36:48: 94%|█████████▍| 47/50 [00:09<00:00, 4.76 it/sec, obj=-3.21e+3]
INFO - 08:36:48: 96%|█████████▌| 48/50 [00:10<00:00, 4.78 it/sec, obj=-3.22e+3]
INFO - 08:36:48: 98%|█████████▊| 49/50 [00:10<00:00, 4.80 it/sec, obj=-3.23e+3]
INFO - 08:36:48: 100%|██████████| 50/50 [00:10<00:00, 4.82 it/sec, obj=-3.24e+3]
INFO - 08:36:48: Optimization result:
INFO - 08:36:48: Optimizer info:
INFO - 08:36:48: Status: None
INFO - 08:36:48: Message: Maximum number of iterations reached. GEMSEO stopped the driver.
INFO - 08:36:48: Number of calls to the objective function by the optimizer: 52
INFO - 08:36:48: Solution:
INFO - 08:36:48: The solution is feasible.
INFO - 08:36:48: Objective: -3240.7083588450046
INFO - 08:36:48: Standardized constraints:
INFO - 08:36:48: g_1_g_2_g_3 = [-3.41793260e-12 -1.26719747e-02 -2.55059715e-02 -3.52857326e-02
INFO - 08:36:48: -4.26719747e-02 -1.80140397e-01 -5.98596026e-02 0.00000000e+00
INFO - 08:36:48: -7.67227050e-01 -2.32772950e-01 2.22044605e-16 -1.83255000e-01]
INFO - 08:36:48: Design space:
INFO - 08:36:48: +-------------+-------------+---------------------+-------------+-------+
INFO - 08:36:48: | Name | Lower bound | Value | Upper bound | Type |
INFO - 08:36:48: +-------------+-------------+---------------------+-------------+-------+
INFO - 08:36:48: | x_shared[0] | 0.01 | 0.05999999999999999 | 0.09 | float |
INFO - 08:36:48: | x_shared[1] | 30000 | 60000 | 60000 | float |
INFO - 08:36:48: | x_shared[2] | 1.4 | 1.4 | 1.8 | float |
INFO - 08:36:48: | x_shared[3] | 2.5 | 3.729811315256278 | 8.5 | float |
INFO - 08:36:48: | x_shared[4] | 40 | 70 | 70 | float |
INFO - 08:36:48: | x_shared[5] | 500 | 1500 | 1500 | float |
INFO - 08:36:48: +-------------+-------------+---------------------+-------------+-------+
INFO - 08:36:48: *** End MDOScenario execution (time: 0:00:10.391182) ***
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)
Plot the system optimization history view#
system_scenario.post_process(post_name="OptHistoryView", save=False, show=True)
<gemseo.post.opt_history_view.OptHistoryView object at 0x7f6dc1162700>
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.optimization_problem)
opt_problem.database = database
execute_post(opt_problem, post_name="OptHistoryView", save=False, show=True)
for disc in [propu, aero, mission, struct]:
print(f"{disc.name}: {disc.execution_statistics.n_calls} calls.")
SobieskiPropulsion: 1642 calls.
SobieskiAerodynamics: 1784 calls.
SobieskiMission: 50 calls.
SobieskiStructure: 1956 calls.
Total running time of the script: (0 minutes 14.708 seconds)