Scalable problem

We want to solve the Aerostructure MDO problem by means of the MDF formulation with a higher dimension for the sweep parameter. For that, we use the ScalableProblem class.

from __future__ import annotations

from gemseo import configure_logger
from gemseo import create_discipline
from gemseo import create_scenario
from gemseo.problems.aerostructure.aerostructure_design_space import (
    AerostructureDesignSpace,
)
from gemseo.problems.scalable.data_driven.problem import ScalableProblem

configure_logger()
<RootLogger root (INFO)>

Define the design problem

In a first step, we define the design problem in terms of objective function (to maximize or minimize), design variables (local and global) and constraints (equality and inequality).

design_variables = ["thick_airfoils", "thick_panels", "sweep"]
objective_function = "range"
eq_constraints = ["c_rf"]
ineq_constraints = ["c_lift"]
maximize_objective = True

Create the disciplinary datasets

Then, we create the disciplinary AbstractFullCache datasets based on a DiagonalDOE.

disciplines = create_discipline(["Aerodynamics", "Structure", "Mission"])
datasets = []
for discipline in disciplines:
    design_space = AerostructureDesignSpace()
    design_space.filter(discipline.get_input_data_names())
    output_names = iter(discipline.get_output_data_names())
    scenario = create_scenario(
        discipline,
        "DisciplinaryOpt",
        next(output_names),
        design_space,
        scenario_type="DOE",
    )
    for output_name in output_names:
        scenario.add_observable(output_name)
    scenario.execute({"algo": "DiagonalDOE", "n_samples": 10})
    datasets.append(scenario.to_dataset(name=discipline.name, opt_naming=False))
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/stable/lib/python3.9/site-packages/gemseo/algos/design_space.py:466: ComplexWarning: Casting complex values to real discards the imaginary part
  self.__current_value[name] = array_value.astype(
    INFO - 10:51:56:
    INFO - 10:51:56: *** Start DOEScenario execution ***
    INFO - 10:51:56: DOEScenario
    INFO - 10:51:56:    Disciplines: Aerodynamics
    INFO - 10:51:56:    MDO formulation: DisciplinaryOpt
    INFO - 10:51:56: Optimization problem:
    INFO - 10:51:56:    minimize drag(thick_airfoils, sweep, displ)
    INFO - 10:51:56:    with respect to displ, sweep, thick_airfoils
    INFO - 10:51:56:    over the design space:
    INFO - 10:51:56:       +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:       | Name           | Lower bound | Value | Upper bound | Type  |
    INFO - 10:51:56:       +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:       | thick_airfoils |      5      |   15  |      25     | float |
    INFO - 10:51:56:       | sweep          |      10     |   25  |      35     | float |
    INFO - 10:51:56:       | displ          |    -1000    |  -700 |     1000    | float |
    INFO - 10:51:56:       +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:56: Solving optimization problem with algorithm DiagonalDOE:
    INFO - 10:51:56:     10%|█         | 1/10 [00:00<00:00, 219.47 it/sec, obj=422]
    INFO - 10:51:56:     20%|██        | 2/10 [00:00<00:00, 351.99 it/sec, obj=336]
    INFO - 10:51:56:     30%|███       | 3/10 [00:00<00:00, 448.51 it/sec, obj=250]
    INFO - 10:51:56:     40%|████      | 4/10 [00:00<00:00, 521.45 it/sec, obj=166]
    INFO - 10:51:56:     50%|█████     | 5/10 [00:00<00:00, 578.14 it/sec, obj=82.3]
    INFO - 10:51:56:     60%|██████    | 6/10 [00:00<00:00, 619.51 it/sec, obj=-.0983]
    INFO - 10:51:56:     70%|███████   | 7/10 [00:00<00:00, 651.29 it/sec, obj=-81.6]
    INFO - 10:51:56:     80%|████████  | 8/10 [00:00<00:00, 671.14 it/sec, obj=-162]
    INFO - 10:51:56:     90%|█████████ | 9/10 [00:00<00:00, 697.11 it/sec, obj=-242]
    INFO - 10:51:56:    100%|██████████| 10/10 [00:00<00:00, 716.18 it/sec, obj=-320]
    INFO - 10:51:56: Optimization result:
    INFO - 10:51:56:    Optimizer info:
    INFO - 10:51:56:       Status: None
    INFO - 10:51:56:       Message: None
    INFO - 10:51:56:       Number of calls to the objective function by the optimizer: 10
    INFO - 10:51:56:    Solution:
    INFO - 10:51:56:       Objective: -319.99905478395067
    INFO - 10:51:56:       Design space:
    INFO - 10:51:56:          +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:          | Name           | Lower bound | Value | Upper bound | Type  |
    INFO - 10:51:56:          +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:          | thick_airfoils |      5      |   25  |      25     | float |
    INFO - 10:51:56:          | sweep          |      10     |   35  |      35     | float |
    INFO - 10:51:56:          | displ          |    -1000    |  1000 |     1000    | float |
    INFO - 10:51:56:          +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:56: *** End DOEScenario execution (time: 0:00:00.027084) ***
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/stable/lib/python3.9/site-packages/gemseo/algos/design_space.py:466: ComplexWarning: Casting complex values to real discards the imaginary part
  self.__current_value[name] = array_value.astype(
    INFO - 10:51:56:
    INFO - 10:51:56: *** Start DOEScenario execution ***
    INFO - 10:51:56: DOEScenario
    INFO - 10:51:56:    Disciplines: Structure
    INFO - 10:51:56:    MDO formulation: DisciplinaryOpt
    INFO - 10:51:56: Optimization problem:
    INFO - 10:51:56:    minimize mass(thick_panels, sweep, forces)
    INFO - 10:51:56:    with respect to forces, sweep, thick_panels
    INFO - 10:51:56:    over the design space:
    INFO - 10:51:56:       +--------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:       | Name         | Lower bound | Value | Upper bound | Type  |
    INFO - 10:51:56:       +--------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:       | thick_panels |      1      |   3   |      20     | float |
    INFO - 10:51:56:       | sweep        |      10     |   25  |      35     | float |
    INFO - 10:51:56:       | forces       |    -1000    |  400  |     1000    | float |
    INFO - 10:51:56:       +--------------+-------------+-------+-------------+-------+
    INFO - 10:51:56: Solving optimization problem with algorithm DiagonalDOE:
    INFO - 10:51:56:     10%|█         | 1/10 [00:00<00:00, 228.81 it/sec, obj=100]
    INFO - 10:51:56:     20%|██        | 2/10 [00:00<00:00, 367.05 it/sec, obj=4.48e+4]
    INFO - 10:51:56:     30%|███       | 3/10 [00:00<00:00, 465.12 it/sec, obj=8.94e+4]
    INFO - 10:51:56:     40%|████      | 4/10 [00:00<00:00, 532.64 it/sec, obj=1.34e+5]
    INFO - 10:51:56:     50%|█████     | 5/10 [00:00<00:00, 587.42 it/sec, obj=1.79e+5]
    INFO - 10:51:56:     60%|██████    | 6/10 [00:00<00:00, 631.23 it/sec, obj=2.23e+5]
    INFO - 10:51:56:     70%|███████   | 7/10 [00:00<00:00, 666.37 it/sec, obj=2.68e+5]
    INFO - 10:51:56:     80%|████████  | 8/10 [00:00<00:00, 691.33 it/sec, obj=3.13e+5]
    INFO - 10:51:56:     90%|█████████ | 9/10 [00:00<00:00, 716.19 it/sec, obj=3.57e+5]
    INFO - 10:51:56:    100%|██████████| 10/10 [00:00<00:00, 721.84 it/sec, obj=4.02e+5]
    INFO - 10:51:56: Optimization result:
    INFO - 10:51:56:    Optimizer info:
    INFO - 10:51:56:       Status: None
    INFO - 10:51:56:       Message: None
    INFO - 10:51:56:       Number of calls to the objective function by the optimizer: 10
    INFO - 10:51:56:    Solution:
    INFO - 10:51:56:       Objective: 100.08573388203513
    INFO - 10:51:56:       Design space:
    INFO - 10:51:56:          +--------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:          | Name         | Lower bound | Value | Upper bound | Type  |
    INFO - 10:51:56:          +--------------+-------------+-------+-------------+-------+
    INFO - 10:51:56:          | thick_panels |      1      |   1   |      20     | float |
    INFO - 10:51:56:          | sweep        |      10     |   10  |      35     | float |
    INFO - 10:51:56:          | forces       |    -1000    | -1000 |     1000    | float |
    INFO - 10:51:56:          +--------------+-------------+-------+-------------+-------+
    INFO - 10:51:56: *** End DOEScenario execution (time: 0:00:00.026830) ***
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/stable/lib/python3.9/site-packages/gemseo/algos/design_space.py:466: ComplexWarning: Casting complex values to real discards the imaginary part
  self.__current_value[name] = array_value.astype(
    INFO - 10:51:56:
    INFO - 10:51:56: *** Start DOEScenario execution ***
    INFO - 10:51:56: DOEScenario
    INFO - 10:51:56:    Disciplines: Mission
    INFO - 10:51:56:    MDO formulation: DisciplinaryOpt
    INFO - 10:51:56: Optimization problem:
    INFO - 10:51:56:    minimize range(drag, lift, mass, reserve_fact)
    INFO - 10:51:56:    with respect to drag, lift, mass, reserve_fact
    INFO - 10:51:56:    over the design space:
    INFO - 10:51:56:       +--------------+-------------+--------+-------------+-------+
    INFO - 10:51:56:       | Name         | Lower bound | Value  | Upper bound | Type  |
    INFO - 10:51:56:       +--------------+-------------+--------+-------------+-------+
    INFO - 10:51:56:       | drag         |     100     |  340   |     1000    | float |
    INFO - 10:51:56:       | lift         |     0.1     |  0.5   |      1      | float |
    INFO - 10:51:56:       | mass         |    100000   | 100000 |    500000   | float |
    INFO - 10:51:56:       | reserve_fact |    -1000    |   0    |     1000    | float |
    INFO - 10:51:56:       +--------------+-------------+--------+-------------+-------+
    INFO - 10:51:56: Solving optimization problem with algorithm DiagonalDOE:
    INFO - 10:51:56:     10%|█         | 1/10 [00:00<00:00, 234.32 it/sec, obj=8e+3+j]
    INFO - 10:51:56:     20%|██        | 2/10 [00:00<00:00, 353.65 it/sec, obj=5.54e+3+j]
    INFO - 10:51:56:     30%|███       | 3/10 [00:00<00:00, 448.88 it/sec, obj=4.24e+3+j]
    INFO - 10:51:56:     40%|████      | 4/10 [00:00<00:00, 519.85 it/sec, obj=3.43e+3+j]
    INFO - 10:51:56:     50%|█████     | 5/10 [00:00<00:00, 569.34 it/sec, obj=2.88e+3+j]
    INFO - 10:51:56:     60%|██████    | 6/10 [00:00<00:00, 612.38 it/sec, obj=2.48e+3+j]
    INFO - 10:51:56:     70%|███████   | 7/10 [00:00<00:00, 647.40 it/sec, obj=2.18e+3+j]
    INFO - 10:51:56:     80%|████████  | 8/10 [00:00<00:00, 676.96 it/sec, obj=1.95e+3+j]
    INFO - 10:51:56:     90%|█████████ | 9/10 [00:00<00:00, 685.34 it/sec, obj=1.76e+3+j]
    INFO - 10:51:56:    100%|██████████| 10/10 [00:00<00:00, 705.97 it/sec, obj=1.6e+3+j]
    INFO - 10:51:56: Optimization result:
    INFO - 10:51:56:    Optimizer info:
    INFO - 10:51:56:       Status: None
    INFO - 10:51:56:       Message: None
    INFO - 10:51:56:       Number of calls to the objective function by the optimizer: 10
    INFO - 10:51:56:    Solution:
    INFO - 10:51:56:       Objective: (1600+0j)
    INFO - 10:51:56:       Design space:
    INFO - 10:51:56:          +--------------+-------------+--------+-------------+-------+
    INFO - 10:51:56:          | Name         | Lower bound | Value  | Upper bound | Type  |
    INFO - 10:51:56:          +--------------+-------------+--------+-------------+-------+
    INFO - 10:51:56:          | drag         |     100     |  1000  |     1000    | float |
    INFO - 10:51:56:          | lift         |     0.1     |   1    |      1      | float |
    INFO - 10:51:56:          | mass         |    100000   | 500000 |    500000   | float |
    INFO - 10:51:56:          | reserve_fact |    -1000    |  1000  |     1000    | float |
    INFO - 10:51:56:          +--------------+-------------+--------+-------------+-------+
    INFO - 10:51:56: *** End DOEScenario execution (time: 0:00:00.028322) ***

Instantiate a scalable problem

In a third stage, we instantiate a ScalableProblem from these disciplinary datasets and from the definition of the MDO problem. We also increase the dimension of the sweep parameter.

problem = ScalableProblem(
    datasets,
    design_variables,
    objective_function,
    eq_constraints,
    ineq_constraints,
    maximize_objective,
    sizes={"sweep": 2},
)
print(problem)
MDO problem
   Disciplines: Aerodynamics, Structure, Mission
   Design variables: thick_airfoils, thick_panels, sweep
   Objective function: range (to maximize)
   Inequality constraints: c_lift
   Equality constraints: c_rf
   Sizes: displ (1), sweep (2), thick_airfoils (1), drag (1), forces (1), lift (1), thick_panels (1), mass (1), reserve_fact (1), c_lift (1), c_rf (1), range (1)

Note

We could also provide options to the ScalableModel objects by means of the constructor of ScalableProblem, e.g. fill_factor in the frame of the ScalableDiagonalModel. In this example, we use the standard ones.

Visualize the N2 chart

We can see the coupling between disciplines through this N2 chart:

problem.plot_n2_chart(save=False, show=True)
plot problem

Create an MDO scenario

Lastly, we create an MDOScenario with the MDF formulation and start the optimization at equilibrium, thus ensuring the feasibility of the first iterate.

scenario = problem.create_scenario("MDF", start_at_equilibrium=True)
INFO - 10:51:57: Build a preliminary MDA to start at equilibrium

Note

We could also provide options for the scalable models to the constructor of ScalableProblem, e.g. fill_factor in the frame of the ScalableDiagonalModel. In this example, we use the standard ones.

Once the scenario is created, we can execute it as any scenario. Here, we use the NLOPT_SLSQP optimization algorithm with no more than 100 iterations.

scenario.execute({"algo": "NLOPT_SLSQP", "max_iter": 100})
    INFO - 10:51:57:
    INFO - 10:51:57: *** Start MDOScenario execution ***
    INFO - 10:51:57: MDOScenario
    INFO - 10:51:57:    Disciplines: sdm_Aerodynamics sdm_Mission sdm_Structure
    INFO - 10:51:57:    MDO formulation: MDF
    INFO - 10:51:57: Optimization problem:
    INFO - 10:51:57:    minimize -range(thick_airfoils, thick_panels, sweep)
    INFO - 10:51:57:    with respect to sweep, thick_airfoils, thick_panels
    INFO - 10:51:57:    subject to constraints:
    INFO - 10:51:57:       c_lift(thick_airfoils, thick_panels, sweep) <= [0.74554856]
    INFO - 10:51:57:       c_rf(thick_airfoils, thick_panels, sweep) == 0.49642016361892943
    INFO - 10:51:57:    over the design space:
    INFO - 10:51:57:       +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:57:       | Name           | Lower bound | Value | Upper bound | Type  |
    INFO - 10:51:57:       +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:57:       | thick_airfoils |      0      |  0.5  |      1      | float |
    INFO - 10:51:57:       | thick_panels   |      0      |  0.5  |      1      | float |
    INFO - 10:51:57:       | sweep[0]       |      0      |  0.5  |      1      | float |
    INFO - 10:51:57:       | sweep[1]       |      0      |  0.5  |      1      | float |
    INFO - 10:51:57:       +----------------+-------------+-------+-------------+-------+
    INFO - 10:51:57: Solving optimization problem with algorithm NLOPT_SLSQP:
    INFO - 10:51:57:      1%|          | 1/100 [00:00<00:02, 33.41 it/sec, obj=-.168]
    INFO - 10:51:57:      2%|▏         | 2/100 [00:00<00:11,  8.24 it/sec, obj=-.172]
    INFO - 10:51:57:      3%|▎         | 3/100 [00:00<00:10,  9.28 it/sec, obj=-.2]
    INFO - 10:51:57:      4%|▍         | 4/100 [00:00<00:09,  9.84 it/sec, obj=-.302]
    INFO - 10:51:57:      5%|▌         | 5/100 [00:00<00:09, 10.36 it/sec, obj=-.302]
    INFO - 10:51:57:      6%|▌         | 6/100 [00:00<00:09, 10.34 it/sec, obj=-.303]
    INFO - 10:51:58:      7%|▋         | 7/100 [00:00<00:08, 10.51 it/sec, obj=-.304]
    INFO - 10:51:58:      8%|▊         | 8/100 [00:00<00:08, 10.80 it/sec, obj=-.309]
    INFO - 10:51:58:      9%|▉         | 9/100 [00:00<00:08, 11.03 it/sec, obj=-.309]
   ERROR - 10:51:58: NLopt run failed: NLopt roundoff-limited, RoundoffLimited
Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/stable/lib/python3.9/site-packages/gemseo/algos/opt/lib_nlopt.py", line 498, in _run
    nlopt_problem.optimize(x_0.real)
  File "/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/stable/lib/python3.9/site-packages/nlopt/nlopt.py", line 335, in optimize
    return _nlopt.opt_optimize(self, *args)
nlopt.RoundoffLimited: NLopt roundoff-limited
    INFO - 10:51:58:     10%|█         | 10/100 [00:00<00:07, 12.20 it/sec, obj=-.309]
    INFO - 10:51:58: Optimization result:
    INFO - 10:51:58:    Optimizer info:
    INFO - 10:51:58:       Status: None
    INFO - 10:51:58:       Message:  GEMSEO Stopped the driver
    INFO - 10:51:58:       Number of calls to the objective function by the optimizer: 11
    INFO - 10:51:58:    Solution:
    INFO - 10:51:58:       The solution is feasible.
    INFO - 10:51:58:       Objective: -0.30937609894434914
    INFO - 10:51:58:       Standardized constraints:
    INFO - 10:51:58:          [c_lift+offset] = [-0.42031165]
    INFO - 10:51:58:          [c_rf-0.49642016361892943] = -2.7755575615628914e-16
    INFO - 10:51:58:       Design space:
    INFO - 10:51:58:          +----------------+-------------+--------------------+-------------+-------+
    INFO - 10:51:58:          | Name           | Lower bound |       Value        | Upper bound | Type  |
    INFO - 10:51:58:          +----------------+-------------+--------------------+-------------+-------+
    INFO - 10:51:58:          | thick_airfoils |      0      | 0.3004770228426673 |      1      | float |
    INFO - 10:51:58:          | thick_panels   |      0      |         1          |      1      | float |
    INFO - 10:51:58:          | sweep[0]       |      0      |         1          |      1      | float |
    INFO - 10:51:58:          | sweep[1]       |      0      |         1          |      1      | float |
    INFO - 10:51:58:          +----------------+-------------+--------------------+-------------+-------+
    INFO - 10:51:58: *** End MDOScenario execution (time: 0:00:00.836733) ***

{'max_iter': 100, 'algo': 'NLOPT_SLSQP'}

We can post-process the results. Here, we use the standard OptHistoryView.

scenario.post_process("OptHistoryView", save=False, show=True)
  • 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 equality constraints
<gemseo.post.opt_history_view.OptHistoryView object at 0x7f1de1da2940>

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

Gallery generated by Sphinx-Gallery