Surrogate models - Advanced features

A deeper dive into meta-models created by make_surrogate.

Update of sub-system variables at each surrogate model execution

As discussed in the tutorial on surrogate models, a meta-model, when activated, supersedes the behaviour of a system as originally defined by its compute() method. The meta-model is a black box whose inputs are specified by the Design of Experiment (DOE) on which it was trained. But what about its outputs?

At the very least, the output of a system surrogate should match those of the original system. CoSApp runs the extra mile, by adding all outputs and connected inputs of all sub-systems to the list of meta-model outputs. Moreover, at each surrogate model execution, all sub-system variables are synchronized with the computed meta-model outputs.

This way, the meta-modeled system and its sub-systems are kept in a state consistent with its free inputs.

Consider the systems defined in the tutorial on surrogate models:

from cosapp.base import System


class MultiplyByM(System):
    """System computing y = m * x"""
    def setup(self):
        self.add_inward('m', 1.0)
        self.add_inward('x')
        self.add_outward('y')

    def compute(self):
        self.y = self.m * self.x


class AffineSystem(System):
    """System computing y = a.m * b.m * x + 3"""
    def setup(self):
        self.add_outward('y')
        # sub-systems
        a = self.add_child(MultiplyByM('a'), pulling='x')
        b = self.add_child(MultiplyByM('b'))
        # connections
        self.connect(a, b, {'y': 'x'})

    def compute(self):
        self.y = self.b.y + 3

    @property
    def y_exact(self) -> float:
        """Expected value for `self.y`"""
        return self.b.m * (self.a.m * self.x) + 3

import pandas
import itertools

def Cartesian_DoE(axes: dict) -> pandas.DataFrame:
    """Simple Cartesian grid DoE from 1D samples in all axis directions"""
    return pandas.DataFrame(
        itertools.product(*axes.values()),
        columns=axes.keys(),
    )

import numpy

axes = {
    'x': numpy.linspace(-5, 5, 11),
    'a.m': numpy.linspace(-3, 3, 6),
    'b.m': numpy.linspace(-3, 3, 6),
}

doe = Cartesian_DoE(axes)

head = AffineSystem('head')
meta = head.make_surrogate(doe)

print(
    f"Outputs of {head.name!r}:",
    dict(filter(lambda item: len(item[1]) > 0, head.outputs.items())),  # non-empty output ports
    f"\nOutputs of meta-model:",
    meta.synched_outputs,
    sep='\n',
)

Outputs of 'head':
{'outwards': ExtensiblePort: {'y': 48.0}}

Outputs of meta-model:
['outwards.y', 'a.inwards.x', 'a.outwards.y', 'b.inwards.x', 'b.outwards.y']

As can be seen in previous cell, system head has only one output y, whereas meta has five, the last four of which are sub-system variables. Every time head is executed, its meta-model not only updates head.y, but also head.a and head.b variables, as would be the case with the original system.

In cell below, we show the internal state of system head with or without meta-model bypass.

def run_and_print(head: AffineSystem, activate: bool):
    inputs = ['x']
    outputs = ['y', 'a.x', 'a.y', 'b.x', 'b.y']
    for var in outputs:
        head[var] = -999.999  # set to bogus value
    # Set surrogate status and run model
    head.active_surrogate = activate
    head.run_once()
    # Print info on internal variables
    print('\nActivated meta-model:', head.active_surrogate)
    for var in inputs + outputs:
        print(f"{head.name}.{var}: ", head[var], sep='\t')

head.x = 3.14
head.a.m = 0.25
head.b.m = -1

run_and_print(head, activate=False)
run_and_print(head, activate=True)

Activated meta-model: False
head.x:         3.14
head.y:         2.215
head.a.x:       3.14
head.a.y:       0.785
head.b.x:       0.785
head.b.y:       -0.785

Activated meta-model: True
head.x:         3.14
head.y:         2.214842218714314
head.a.x:       3.1471118198553447
head.a.y:       0.7871162151169441
head.b.x:       0.7871162151169441
head.b.y:       -0.7851577812874428

After meta-model activation, top-level output head.y is estimated with a numerical error.

Noticeably, it also appears that internal connections are accounted for by the meta-model, as a.y and b.x are equal, for example. Moreover, the behaviour \(y = m\,x\) is also satisfied (within numerical error) for sub-systems a and b, with head.a.m = 0.25 and head.b.m = -1. Precision is rather poor here, because the meta-model was trained on scarce data.

Post-synchronization can be altered at meta-model creation, with optional argument postsynch, specifying a name pattern (or a list thereof) for variables that must be synchronized. By default, this argument is set to '*', meaning “synchronize everything”. An empty list or None value means “do not synchronize anything apart from top-level outputs”.

In the example below, we specify that only variables of sub-system head.a should be synchronized:

head = AffineSystem('head')
meta = head.make_surrogate(doe, postsynch='a.*')

print(
    f"Outputs of meta-model:",
    meta.synched_outputs,
    sep='\n',
)
head.x = 3.14
head.a.m = 0.25
head.b.m = -1

run_and_print(head, activate=False)
run_and_print(head, activate=True)

Outputs of meta-model:
['y', 'a.y']

Activated meta-model: False
head.x:         3.14
head.y:         2.215
head.a.x:       3.14
head.a.y:       0.785
head.b.x:       0.785
head.b.y:       -0.785

Activated meta-model: True
head.x:         3.14
head.y:         2.2144689816859966
head.a.x:       3.14
head.a.y:       0.7873347384770749
head.b.x:       -999.999
head.b.y:       -999.999

As requested, we observe that ports b.x and b.y have not been modified after meta-model execution.

Limiting post-synchronization may also help reduce training time.

Training a sub-system

Methods

• make_surrogate

• dump_surrogate

• load_surrogate

• active_surrogate

can be called at any system level.

Systems above the meta-modeled module (i.e. parent systems) will be unaffected by the procedure. Systems below (child systems), on the other hand, will automatically be deactivated, together with their drivers, if any.

In the example below, we create a meta-model on sub-system head.a alone:

head = AffineSystem('head')

axes = {
    'm': numpy.linspace(-2, 2, 11),
    'x': numpy.linspace(-10, 10, 21),
}

doe = Cartesian_DoE(axes)
head.a.make_surrogate(doe, activate=False)  # create meta-model, but do not activate it

# Sanity check:
print(
    f"{head.has_surrogate = }",
    f"{head.a.has_surrogate = }",
    "",
    sep="\n",
)

head.x = 5.2
head.a.m = 1.55

head.run_drivers()
print(f"{head.a.active_surrogate = }; \t{head.y = }")

head.a.active_surrogate = True
head.run_drivers()
print(f"{head.a.active_surrogate = }; \t{head.y = }")

head.has_surrogate = False
head.a.has_surrogate = True

head.a.active_surrogate = False;        head.y = 11.06
head.a.active_surrogate = True;         head.y = 11.060051256302355

Impact of training space

What happens when the meta-model is trained on a subset of inputs?

In the example below, we again create a meta-model for sub-system head.a, but do not include input m in the DOE:

head = AffineSystem('head')
head.a.m = 1.55

axes = {
    # 'm': numpy.linspace(-2, 2, 11),  # do not train on `m`
    'x': numpy.linspace(-10, 10, 21),
}

doe = Cartesian_DoE(axes)
head.a.make_surrogate(doe, activate=False)  # create meta-model, but do not activate it

head.x = 5.2

head.run_drivers()
print(f"{head.a.active_surrogate = }; \t{head.y = }")

head.a.active_surrogate = True
head.run_drivers()
print(f"{head.a.active_surrogate = }; \t{head.y = }")

head.a.active_surrogate = False;        head.y = 11.06
head.a.active_surrogate = True;         head.y = 11.060211519142264

Here, the surrogate model was trained assuming head.a.m is constant and equals 1.55. The result after activation is correct, as head.a.m was unchanged.

However, results are incorrect when this parameter is modified:

head.a.m = 0.83

head.run_drivers()
print(f"{head.y = } instead of {head.y_exact}")

head.y = 11.060211519142264 instead of 7.316

Tip: Before creating a system surrogate, be sure to lookup variables that may have an impact on outputs:

soi = head.a  # system of interest

print(
    "Input ports:", dict(soi.inputs),
    "",
    f"{list(soi.unknowns) = }",
    f"{list(soi.transients) = }",
    sep="\n",
)

Input ports:
{'inwards': ExtensiblePort: {'m': 0.83, 'x': 5.2}, 'modevars_in': ModeVarPort: {}}

list(soi.unknowns) = []
list(soi.transients) = []

Be aware that this technique does not report sub-system inputs that may influence outputs, though.

Surrogate models and solvers

Good practices when the system of interest has unknowns.

Consider a class ClosedSystem, assembling one sub-system usys (of type UnknownSystem) bearing an unknown, and one sub-system esys (of type EquationSystem) containing an equation.

class UnknownSystem(System):
    def setup(self):
        self.add_inward('m', 1.0, desc="Multiplying factor")
        self.add_inward('x', 1.0)
        self.add_outward('y')
        # Unknowns
        self.add_unknown('m')

    def compute(self):
        self.y = self.x * self.m

class EquationSystem(System):
    def setup(self, p_target=25):
        self.add_inward('p', 1.0)
        self.add_property('p_target', p_target)
        self.add_equation("p == p_target")

class ClosedSystem(System):
    def setup(self):
        usys = self.add_child(UnknownSystem('usys'), pulling='x')
        esys = self.add_child(EquationSystem('esys'))

        self.connect(usys, esys, {'y': 'p'})

    @property
    def solution(self) -> float:
        """Exact solution for self.usys.m"""
        try:
            return self.esys.p_target / self.x

        except ZeroDivisionError:
            return numpy.sign(self.esys.p_target) * numpy.inf

from cosapp.drivers import NonLinearSolver
from cosapp.utils.surrogate_models import (
    FloatKrigingSurrogate,
    LinearNearestNeighbor,
)

Case 1: The system of interest is well-posed

In this part, the system of interest is the top assembly, mathematically closed.

If the system has a NonLinearSolver driver before the generation of the meta-model, CoSApp automatically tracks unknowns, and you don’t have to add them to the input training dataset.

head = ClosedSystem("head")
head.add_driver(NonLinearSolver('solver', tol=1e-9))
# Check that assembly system + driver can be solved
print(f"{head.is_standalone() = }")

doe = Cartesian_DoE({
    'x': numpy.linspace(0.1, 10, 20),
})
meta = head.make_surrogate(doe, model=LinearNearestNeighbor)

# Check that `usys.m` is tracked by meta-model
print(f"{meta.synched_outputs = }")

head.is_standalone() = True
meta.synched_outputs = ['usys.m', 'usys.inwards.x', 'usys.outwards.y', 'esys.inwards.p']

import plotly.graph_objects as go
from plotly.subplots import make_subplots

def make_figure(head, meta=None, x_values=numpy.linspace(3, 10, 21)):
    """
    Utility function to compare the behaviour of a `ClosedSystem`
    instance `head` with and without meta-modeling on head[meta].
    """
    fig = make_subplots(
        rows=1,
        cols=2,
        subplot_titles=("Solution", "Relative difference"),
    )

    system = head if meta is None else head[meta]
    meta_results = numpy.zeros_like(x_values)
    normal_results = numpy.zeros_like(x_values)

    for k, head.x in enumerate(x_values):
        system.active_surrogate = True
        head.run_drivers()
        meta_results[k] = head.usys.m

        system.active_surrogate = False
        head.run_drivers()
        normal_results[k] = head.usys.m

    ############################### BUILDING FIGURE ###############################

    def add_scatter(figure, row, col, y, **options):
        figure.add_trace(
            go.Scatter(x=x_values, y=y, **options),
            row=row, col=col,
        )

    add_scatter(fig, 1, 1, meta_results,
        name='Meta-model on',
    )
    add_scatter(fig, 1, 1, normal_results,
        name='Meta-model off',
        mode='markers',
    )
    add_scatter(fig, 1, 2,
        y = numpy.absolute(meta_results / normal_results - 1),
        name='Relative diff',
        mode='markers',
    )

    ############################### GRAPH LAYOUT ###############################
    sysname = f"{head.name}" if meta is None else f"{head.name}.{system.name}"
    fig.update_layout(
        title=dict(
            text=f"Solution with and without meta-model on {sysname!r}",
            y=0.95,
            x=0.5,
            xanchor='center',
            yanchor='top',
        ),
        legend_title_text='Legend:',
        hovermode='x',
    )
    fig.update_xaxes(title_text="head.x", row=1, col=1)
    fig.update_xaxes(title_text="head.x", row=1, col=2)
    fig.update_yaxes(title_text="head.usys.m", row=1, col=1)
    fig.update_yaxes(title_text="|meta / normal - 1|", row=1, col=2)
    fig.update_traces(showlegend=True)

    return fig

fig = make_figure(head)
fig.show()

Case 2: The system of interest is ill-posed

Here, we assume that the system of interest has an unbalanced number of unknowns and equations, such that it cannot be solved on its own.

This is the case of sub-system usys of assembly ClosedSystem defined above.

head = ClosedSystem("head")
head.add_driver(NonLinearSolver('solver', tol=1e-9))

doe = Cartesian_DoE({
    'm': numpy.linspace(-10, 10, 21),
    'x': numpy.linspace(-10, 10, 20),
})
head.usys.make_surrogate(doe, activate=False)

print(
    f"{head.is_standalone() = }",   # top assembly is solvable
    f"{head.usys.is_standalone() = }",  # sub-system `usys` is not
    sep="\n",
)

head.is_standalone() = True
head.usys.is_standalone() = False

head.x = 4.

print(f"Exact solution: head.usys.m = {head.solution}")

def solve_and_print(activate: bool):
    head.usys.active_surrogate = activate
    head.run_drivers()
    print(
        f"\nMeta-model activated: {head.usys.active_surrogate}",
        f"{head.usys.m = }",
        f"Error: {head.esys.p - head.esys.p_target = }",
        sep="\n",
    )

solve_and_print(activate=False)
solve_and_print(activate=True)

Exact solution: head.usys.m = 6.25

Meta-model activated: False
head.usys.m = 6.25
Error: head.esys.p - head.esys.p_target = 0.0

Meta-model activated: True
head.usys.m = 6.250589017670577
Error: head.esys.p - head.esys.p_target = -2.27018404075352e-12

fig = make_figure(head, meta='usys')
fig.show()

Depending on model type (FloatKrigingSurrogate, LinearNearestNeighbor…), the training precision for your unknowns will vary. If the DoE space is not sufficiently large, NonLinearSolver may fail.

As a general rule of thumb, it is better to have an idea of the solution, so you can train in a range where NonLinearSolver is not likely to fail.

Common Pitfall

A system unknown is left out of the training set:

head = ClosedSystem("head")

doe = Cartesian_DoE({
    # 'm': numpy.linspace(-10, 10, 21),  # `m` left out on purpose
    'x': numpy.linspace(-10, 10, 21),
})

head.usys.make_surrogate(doe)

print(f"{head.usys.active_surrogate = }")

head.usys.active_surrogate = True

/home/docs/checkouts/readthedocs.org/user_builds/cosapp/conda/latest/lib/python3.11/site-packages/cosapp/systems/systemSurrogate.py:348: UserWarning:

The following unknowns/transients are not part of the training set; future attempts to compute them with a driver may fail: ['m']

Since m is not part of the training set, the surrogate model cannot predict its impact on outputs.

As a consequence, NonLinearSolver will fail:

head.add_driver(NonLinearSolver('solver', tol=1e-9))

head.x = 4.5
head.run_drivers()

/home/docs/checkouts/readthedocs.org/user_builds/cosapp/conda/latest/lib/python3.11/site-packages/cosapp/core/numerics/root.py:502: LinAlgWarning:

Diagonal number 1 is exactly zero. Singular matrix.

The solver failed: Singular 1x1 Jacobian matrix
The 1 following parameter(s) have no influence: ['m']
The 1 following residue(s) are not influenced: ['esys: p == p_target']