Store observables

Introduction

In this example, we will learn how to store the history of state variables using the add_observable() method. This is useful in situations where we wish to access, post-process, or save the values of discipline outputs that are not design variables, constraints or objective functions.

The Sellar problem

We will consider in this example the Sellar problem:

\[\begin{split}\begin{aligned} \text{minimize the objective function }&\text{obj}=x_{\text{local}}^2 + x_{\text{shared},2} +y_1^2+e^{-y_2} \\ \text{with respect to the design variables }&x_{\text{shared}},\,x_{\text{local}} \\ \text{subject to the general constraints } & c_1 \leq 0\\ & c_2 \leq 0\\ \text{subject to the bound constraints } & -10 \leq x_{\text{shared},1} \leq 10\\ & 0 \leq x_{\text{shared},2} \leq 10\\ & 0 \leq x_{\text{local}} \leq 10. \end{aligned}\end{split}\]

where the coupling variables are

\[\text{Discipline 1: } y_1 = \sqrt{x_{\text{shared},1}^2 + x_{\text{shared},2} + x_{local} - 0.2\,y_2},\]

and

\[\text{Discipline 2: }y_2 = |y_1| + x_{\text{shared},1} + x_{\text{shared},2}.\]

and where the general constraints are

\[ \begin{align}\begin{aligned}c_1 = 3.16 - y_1^2\\c_2 = y_2 - 24.\end{aligned}\end{align} \]

Imports

All the imports needed for the tutorials are performed here.

from __future__ import annotations

from gemseo import configure_logger
from gemseo import create_discipline
from gemseo import create_scenario
from gemseo.algos.design_space import DesignSpace
from numpy import array
from numpy import ones

configure_logger()
<RootLogger root (INFO)>

Create the problem disciplines

In this section, we use the available classes Sellar1, Sellar2 and SellarSystem to define the disciplines of the problem. The create_discipline() API function allows us to carry out this task easily, as well as store the instances in a list to be used later on.

disciplines = create_discipline(["Sellar1", "Sellar2", "SellarSystem"])

Create and execute the scenario

Create the design space

In this section, we define the design space which will be used for the creation of the MDOScenario.

design_space = DesignSpace()
design_space.add_variable("x_local", l_b=0.0, u_b=10.0, value=ones(1))
design_space.add_variable(
    "x_shared", 2, l_b=(-10, 0.0), u_b=(10.0, 10.0), value=array([4.0, 3.0])
)

Create the scenario

In this section, we build the MDO scenario which links the disciplines with the formulation, the design space and the objective function.

scenario = create_scenario(
    disciplines, formulation="MDF", objective_name="obj", design_space=design_space
)

Add the constraints

Then, we have to set the design constraints

scenario.add_constraint("c_1", "ineq")
scenario.add_constraint("c_2", "ineq")

Add the observables

Only the design variables, objective function and constraints are stored by default. In order to be able to recover the data from the state variables, y1 and y2, we have to add them as observables. All we have to do is enter the variable name as a string to the add_observable() method. If more than one output name is provided (as a list of strings), the observable function returns a concatenated array of the output values.

scenario.add_observable("y_1")

It is also possible to add the observable with a custom name, using the option observable_name. Let us store the variable y_2 as y2.

scenario.add_observable("y_2", observable_name="y2")

Execute the scenario

Then, we execute the MDO scenario with the inputs of the MDO scenario as a dictionary. In this example, the gradient-based SLSQP optimizer is selected, with 10 iterations at maximum:

scenario.execute(input_data={"max_iter": 10, "algo": "SLSQP"})
    INFO - 16:25:55:
    INFO - 16:25:55: *** Start MDOScenario execution ***
    INFO - 16:25:55: MDOScenario
    INFO - 16:25:55:    Disciplines: Sellar1 Sellar2 SellarSystem
    INFO - 16:25:55:    MDO formulation: MDF
    INFO - 16:25:55: Optimization problem:
    INFO - 16:25:55:    minimize obj(x_local, x_shared)
    INFO - 16:25:55:    with respect to x_local, x_shared
    INFO - 16:25:55:    subject to constraints:
    INFO - 16:25:55:       c_1(x_local, x_shared) <= 0.0
    INFO - 16:25:55:       c_2(x_local, x_shared) <= 0.0
    INFO - 16:25:55:    over the design space:
    INFO - 16:25:55:    +-------------+-------------+-------+-------------+-------+
    INFO - 16:25:55:    | name        | lower_bound | value | upper_bound | type  |
    INFO - 16:25:55:    +-------------+-------------+-------+-------------+-------+
    INFO - 16:25:55:    | x_local     |      0      |   1   |      10     | float |
    INFO - 16:25:55:    | x_shared[0] |     -10     |   4   |      10     | float |
    INFO - 16:25:55:    | x_shared[1] |      0      |   3   |      10     | float |
    INFO - 16:25:55:    +-------------+-------------+-------+-------------+-------+
    INFO - 16:25:55: Solving optimization problem with algorithm SLSQP:
    INFO - 16:25:55: ...   0%|          | 0/10 [00:00<?, ?it]
    INFO - 16:25:55: ...  10%|█         | 1/10 [00:00<00:00, 30.74 it/sec, obj=21.8+j]
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/5.0.0/lib/python3.9/site-packages/scipy/optimize/_slsqp_py.py:422: ComplexWarning: Casting complex values to real discards the imaginary part
  slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw,
    INFO - 16:25:55: ...  20%|██        | 2/10 [00:00<00:00, 24.23 it/sec, obj=5.39+j]
    INFO - 16:25:56: ...  30%|███       | 3/10 [00:00<00:00, 23.32 it/sec, obj=3.41+j]
    INFO - 16:25:56: ...  40%|████      | 4/10 [00:00<00:00, 22.78 it/sec, obj=3.19+j]
    INFO - 16:25:56: ...  50%|█████     | 5/10 [00:00<00:00, 22.44 it/sec, obj=3.18+j]
    INFO - 16:25:56: ...  60%|██████    | 6/10 [00:00<00:00, 22.26 it/sec, obj=3.18+j]
    INFO - 16:25:56: ...  70%|███████   | 7/10 [00:00<00:00, 22.11 it/sec, obj=3.18+j]
    INFO - 16:25:56: ...  80%|████████  | 8/10 [00:00<00:00, 24.41 it/sec, obj=3.18+j]
    INFO - 16:25:56: Optimization result:
    INFO - 16:25:56:    Optimizer info:
    INFO - 16:25:56:       Status: None
    INFO - 16:25:56:       Message: Successive iterates of the objective function are closer than ftol_rel or ftol_abs. GEMSEO Stopped the driver
    INFO - 16:25:56:       Number of calls to the objective function by the optimizer: 9
    INFO - 16:25:56:    Solution:
    INFO - 16:25:56:       The solution is feasible.
    INFO - 16:25:56:       Objective: (3.1833939452393087+0j)
    INFO - 16:25:56:       Standardized constraints:
    INFO - 16:25:56:          c_1 = (6.486713388653698e-09+0j)
    INFO - 16:25:56:          c_2 = (-20.244722236724627+0j)
    INFO - 16:25:56:       Design space:
    INFO - 16:25:56:       +-------------+-------------+-------------------+-------------+-------+
    INFO - 16:25:56:       | name        | lower_bound |       value       | upper_bound | type  |
    INFO - 16:25:56:       +-------------+-------------+-------------------+-------------+-------+
    INFO - 16:25:56:       | x_local     |      0      |         0         |      10     | float |
    INFO - 16:25:56:       | x_shared[0] |     -10     | 1.977638881636786 |      10     | float |
    INFO - 16:25:56:       | x_shared[1] |      0      |         0         |      10     | float |
    INFO - 16:25:56:       +-------------+-------------+-------------------+-------------+-------+
    INFO - 16:25:56: *** End MDOScenario execution (time: 0:00:00.342063) ***

{'max_iter': 10, 'algo': 'SLSQP'}

Access the observable variables

Retrieve observables from a dataset

In order to create a dataset, we use the corresponding OptimizationProblem:

opt_problem = scenario.formulation.opt_problem

We can easily build an OptimizationDataset from this OptimizationProblem: either by separating the design parameters from the functions (default option):

dataset = opt_problem.to_dataset("sellar_problem")
print(dataset)
GROUP       designs                          ...  functions
VARIABLE    x_local  x_shared                ...        obj         y2       y_1
COMPONENT         0         0             1  ...          0          0         0
1          1.000000  4.000000  3.000000e+00  ...  21.757227  11.213931  4.213931
2          0.000000  2.520049  0.000000e+00  ...   5.390534   4.840098  2.320049
3          0.003071  2.040410  1.594470e-08  ...   3.410649   3.881611  1.841201
4          0.000000  1.978711  5.354587e-12  ...   3.187156   3.757421  1.778711
5          0.000000  1.977639  5.040991e-13  ...   3.183395   3.755278  1.777639
6          0.000000  1.977639  3.942508e-13  ...   3.183394   3.755278  1.777639
7          0.000000  1.977639  0.000000e+00  ...   3.183394   3.755278  1.777639
8          0.000000  1.977639  1.971886e-13  ...   3.183394   3.755278  1.777639

[8 rows x 8 columns]

or by considering all features as default parameters:

dataset = opt_problem.to_dataset("sellar_problem", categorize=False)
print(dataset)
GROUP     parameters                          ...
VARIABLE     x_local  x_shared                ...        obj         y2       y_1
COMPONENT          0         0             1  ...          0          0         0
0           1.000000  4.000000  3.000000e+00  ...  21.757227  11.213931  4.213931
1           0.000000  2.520049  0.000000e+00  ...   5.390534   4.840098  2.320049
2           0.003071  2.040410  1.594470e-08  ...   3.410649   3.881611  1.841201
3           0.000000  1.978711  5.354587e-12  ...   3.187156   3.757421  1.778711
4           0.000000  1.977639  5.040991e-13  ...   3.183395   3.755278  1.777639
5           0.000000  1.977639  3.942508e-13  ...   3.183394   3.755278  1.777639
6           0.000000  1.977639  0.000000e+00  ...   3.183394   3.755278  1.777639
7           0.000000  1.977639  1.971886e-13  ...   3.183394   3.755278  1.777639

[8 rows x 8 columns]

or by using an input-output naming rather than an optimization naming:

dataset = opt_problem.to_dataset("sellar_problem", opt_naming=False)
print(dataset)
GROUP        inputs                          ...    outputs
VARIABLE    x_local  x_shared                ...        obj         y2       y_1
COMPONENT         0         0             1  ...          0          0         0
0          1.000000  4.000000  3.000000e+00  ...  21.757227  11.213931  4.213931
1          0.000000  2.520049  0.000000e+00  ...   5.390534   4.840098  2.320049
2          0.003071  2.040410  1.594470e-08  ...   3.410649   3.881611  1.841201
3          0.000000  1.978711  5.354587e-12  ...   3.187156   3.757421  1.778711
4          0.000000  1.977639  5.040991e-13  ...   3.183395   3.755278  1.777639
5          0.000000  1.977639  3.942508e-13  ...   3.183394   3.755278  1.777639
6          0.000000  1.977639  0.000000e+00  ...   3.183394   3.755278  1.777639
7          0.000000  1.977639  1.971886e-13  ...   3.183394   3.755278  1.777639

[8 rows x 8 columns]

Access observables by name

We can get the observable data by variable names:

print(dataset.get_view(variable_names=["y_1", "y2"]))
GROUP       outputs
VARIABLE        y_1         y2
COMPONENT         0          0
0          4.213931  11.213931
1          2.320049   4.840098
2          1.841201   3.881611
3          1.778711   3.757421
4          1.777639   3.755278
5          1.777639   3.755278
6          1.777639   3.755278
7          1.777639   3.755278

Use the observables in a post-processing method

Finally, we can generate plots with the observable variables. Have a look at the Basic History plot and the Scatter Plot Matrix:

scenario.post_process(
    "BasicHistory",
    variable_names=["obj", "y_1", "y2"],
    save=False,
    show=True,
)
scenario.post_process(
    "ScatterPlotMatrix",
    variable_names=["obj", "c_1", "c_2", "y2", "y_1"],
    save=False,
    show=True,
)
  • History plot
  • plot store observables
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/5.0.0/lib/python3.9/site-packages/genson/schema/strategies/base.py:42: UserWarning: Schema incompatible. Keyword 'description' has conflicting values ('The width and height of the figure in inches, e.g. ``(w, h)``.\nIf ``None``, use the :attr:`.OptPostProcessor.DEFAULT_FIG_SIZE`\nof the post-processor.' vs. 'The width and height of the figure in inches, e.g. `(w, h)`.\nIf ``None``, use the :attr:`.OptPostProcessor.DEFAULT_FIG_SIZE`\nof the post-processor.'). Using 'The width and height of the figure in inches, e.g. ``(w, h)``.\nIf ``None``, use the :attr:`.OptPostProcessor.DEFAULT_FIG_SIZE`\nof the post-processor.'
  warn(('Schema incompatible. Keyword {0!r} has conflicting '
/home/docs/checkouts/readthedocs.org/user_builds/gemseo/envs/5.0.0/lib/python3.9/site-packages/genson/schema/strategies/base.py:42: UserWarning: Schema incompatible. Keyword 'description' has conflicting values ('The width and height of the figure in inches, e.g. ``(w, h)``.\nIf ``None``, use the :attr:`.OptPostProcessor.DEFAULT_FIG_SIZE`\nof the post-processor.' vs. 'The width and height of the figure in inches, e.g. `(w, h)`.\nIf ``None``, use the :attr:`.OptPostProcessor.DEFAULT_FIG_SIZE`\nof the post-processor.'). Using 'The width and height of the figure in inches, e.g. ``(w, h)``.\nIf ``None``, use the :attr:`.OptPostProcessor.DEFAULT_FIG_SIZE`\nof the post-processor.'
  warn(('Schema incompatible. Keyword {0!r} has conflicting '

<gemseo.post.scatter_mat.ScatterPlotMatrix object at 0x7fea47b36e80>

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

Gallery generated by Sphinx-Gallery