BiLevel-based MDO on the Sobieski SSBJ test case

from __future__ import division, unicode_literals

from copy import deepcopy

from future import standard_library

from gemseo.api import configure_logger, create_discipline, create_scenario
from gemseo.problems.sobieski.core import SobieskiProblem

standard_library.install_aliases()
configure_logger()

Out:

<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().read_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",
)
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",
    deepcopy(design_space).filter("x_2"),
    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",
    deepcopy(design_space).filter("x_1"),
    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).

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,
    apply_cstr_tosub_scenarios=False,
    parallel_scenarios=False,
    multithread_scenarios=True,
    tolerance=1e-14,
    max_mda_iter=30,
    maximize_objective=True,
)
system_scenario.add_constraint(["g_1", "g_2", "g_3"], "ineq")

# system_scenario.xdsmize(open_browser=True)
system_scenario.execute(
    {"max_iter": 50, "algo": "NLOPT_COBYLA", "algo_options": algo_options}
)

Out:

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

Plot the history of the MDA residuals

For the first MDA:

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

Plot the optimization history view

system_scenario.post_process("OptHistoryView", show=True, save=False)
for disc in [propu, aero, mission, struct]:
    print("{}: {} calls.".format(disc.name, disc.n_calls))
  • Evolution of the optimization variables
  • Evolution of the objective value
  • Distance to the optimum
  • Evolution of the inequality constraints

Out:

/home/docs/checkouts/readthedocs.org/user_builds/gemseo/conda/3.0.3/lib/python3.8/site-packages/gemseo/post/opt_history_view.py:312: UserWarning: FixedFormatter should only be used together with FixedLocator
  ax1.set_yticklabels(y_labels)
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/conda/3.0.3/lib/python3.8/site-packages/gemseo/post/opt_history_view.py:626: MatplotlibDeprecationWarning: default base will change from np.e to 10 in 3.4.  To suppress this warning specify the base keyword argument.
  norm=SymLogNorm(linthresh=1.0, vmin=-vmax, vmax=vmax),
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/conda/3.0.3/lib/python3.8/site-packages/gemseo/post/opt_history_view.py:619: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  im1 = ax1.imshow(
SobieskiPropulsion: 1274 calls.
SobieskiAerodynamics: 1514 calls.
SobieskiMission: 50 calls.
SobieskiStructure: 1577 calls.

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

Gallery generated by Sphinx-Gallery