Multi-point design

Up to now, systems have been solved in their up-and-running state, that is all inputs are set, and we want to obtain the state of the system for a given set of boundary conditions.

As CoSApp focuses on system design, we will now show how one can free parameters and compute them to satisfy design equations. In complex systems, parameters are designed under most critical constraints, which usually occur in different operating conditions. As will be illustrated here, CoSApp is able to solve such multi-point design calculations.

Concepts

Multi-point design is based on the interaction between a NonLinearSolver driver, and several RunSingleCase sub-drivers (one per design point).

engine = Turbofan("engine")

solver = engine.add_driver(NonLinearSolver("solver"))

# Add design points:
takeoff = solver.add_child(RunSingleCase("takeoff"))
cruise = solver.add_child(RunSingleCase("cruise"))

Each RunSingleCase driver will contain

• Operating conditions;

• Design unknowns and equations to be satisfied;

• Other nonlinear conditions which may be satisfied at design point only.

Operating conditions are specified by a dictionary of the kind {variable_name: value}, using method set_values:

cruise.set_values({
    "altitude": 3000,
    "rho": "1.21 * exp(-altitude / 8e3)",
    # and so on
})

Additional values may be later appended, with method add_values.

Design and off-design problems

Design parameters refer to variables whose value is independent of operating conditions. Geometric parameters, typically, can be used as design parameters. Design unknowns are unique, and will all be computed jointly by the solver. Therefore, by convention, all unknowns declared at solver level will be regarded as design parameters.

For example, imagine we seek to determine two geometric parameters of our turbofan engine, in order to meet conditions at take-off and cruise. This problem is simply declared by:

solver.add_unknown(["fan.diameter", "core.turbine.inlet.area"])

takeoff.add_equation("thrust == 1.2e5")
cruise.add_equation("Mach == 0.8")

In this example, the fan diameter and turbine inlet area used in system-wide computations will be common to both take-off and cruise design points. The constraint on thrust, though, will be satisfied at take-off only, whereas the Mach number target will be met in cruise conditions.

In practice, specific design parameters are often associated to operating criteria. Thus, for convenience, design unknowns can be associated to design equations at design point level, using attribute design of RunSingleCase drivers. The code below generates a multi-point design problem equivalent to previous snippet:

takeoff.design.add_unknown("core.turbine.inlet.area").add_equation("thrust == 1e4")

cruise.design.add_unknown("fan.diameter").add_equation("Mach == 0.8")

Oppositely, off-design unknowns may be declared several times, in various operating conditions. However, they can (and usually do) converge towards different values at different points. Thus, a local clone of each off-design unknown is created in each design point.

Off-design conditions can be specified at operating point level, using the offdesign problem attribute, instead of design. For the sake of convenience and conciseness, though, offdesign may be omitted. For example,

takeoff.add_unknown("fuel_flowrate").add_equation("core.burner.flow_out.T == 1500")

cruise.add_unknown("fuel_flowrate").add_equation("core.burner.c_nox == 1e-2")

will determine the fuel flowrate at take-off corresponding to a target temperature at combustor outlet; fuel consumption in cruise conditions will be computed to satisfy a criterion on NOx concentration. The value assumed in other points (if any) may also be set independently, using set_values, for example.

Notes

• For a RunSingleCase design point case, case.add_unknown(...) (respectively add_equation, extend) is a shortcut to case.offdesign.add_unknown(...).

• Fixing a parameter using case.set_values({"x": "0.123"}) is mathematically equivalent to the trivial constaint case.add_unknown("x").add_equation("x == 0.123"). Thus, the local off-design problem can be viewed as an extension of set_values for non-trivial, nonlinear constraints.

• In single-point design, solver and case problems (either off-design or design) are interchangeable. However, it is good practice to follow multi-point design rules, for the sake of consistency.

• Equations declared at solver level are satisfied independently on each RunSingleCase subdriver. Thus, in multi-point problems, be aware that a single solver equation will eventually result in several equations in the assembled mathematical problem.

Example

Case description

The simple circuit test case having a constant intensity source will be used here. In addition to the off-design resolution driving the potential V at nodes to balance current fluxes, we seek to determine the value of two resistances from two different operating points.

Design variables: R1.R and R2.R.

Operating point	Boundary conditions	Design equation
Point 1	\(I_{source} = 0.08\) A	n2.V == 8 V
Point 2	\(I_{source} = 0.15\) A	n1.V == 50 V

Building the circuit

The circuit is built as an assembly of elementary models. The default case is then solved to initialize all variables.

from __future__ import annotations
from cosapp.base import Port, System
import abc


class Voltage(Port):
    def setup(self):
        self.add_variable("V", unit="V")


class Intensity(Port):
    def setup(self):
        self.add_variable("I", unit="A")


class Dipole(System):
    """Abstract directional dipole model computing
    current from end values of electric potential.
    """
    def setup(self):
        self.add_input(Voltage, "V_in")
        self.add_input(Voltage, "V_out")

        self.add_output(Intensity, "I")
        self.add_outward("deltaV", unit="V")

    def compute(self):
        self.deltaV = self.V_in.V - self.V_out.V
        self.compute_I()

    @abc.abstractmethod
    def compute_I(self) -> None:
        pass


class Resistor(Dipole):

    def setup(self, R=1.0):
        super().setup()
        self.add_inward("R", R, desc="Resistance in Ohms")

    def compute_I(self):
        self.I.I = self.deltaV / self.R


class Node(System):
    """Electric node model with `n_in` incoming and `n_out` outgoing currents.
    """
    def setup(self, n_in=1, n_out=1):
        self.add_property("n_in", max(1, int(n_in)))
        self.add_property("n_out", max(1, int(n_out)))

        incoming = tuple(
            self.add_input(Intensity, f"I_in{i}")
            for i in range(self.n_in)
        )
        outgoing = tuple(
            self.add_input(Intensity, f"I_out{i}")
            for i in range(self.n_out)
        )
        self.add_property("incoming_currents", incoming)
        self.add_property("outgoing_currents", outgoing)

        self.add_inward("V", 1.0, unit="V")
        self.add_outward("sum_I_in", 0., unit="A", desc="Sum of all incoming currents")
        self.add_outward("sum_I_out", 0., unit="A", desc="Sum of all outgoing currents")

        self.add_unknown("V")
        self.add_equation("sum_I_in == sum_I_out", name="current balance")

    def compute(self):
        self.sum_I_in = sum(current.I for current in self.incoming_currents)
        self.sum_I_out = sum(current.I for current in self.outgoing_currents)

    @classmethod
    def make(
        cls,
        name: str,
        parent: System,
        incoming: list[Dipole]=[],
        outgoing: list[Dipole]=[],
        pulling=None,
    ) -> Node:
        """Factory creating new node within `parent`, with
        appropriate connections with incoming and outgoing dipoles.
        """
        node = cls(name, n_in=len(incoming), n_out=len(outgoing))
        parent.add_child(node, pulling=pulling)

        for dipole, current in zip(incoming, node.incoming_currents):
            # print(dipole.name, type(dipole), dipole.I)
            parent.connect(dipole.I, current)
            parent.connect(dipole.V_out, node.inwards, "V")

        for dipole, current in zip(outgoing, node.outgoing_currents):
            parent.connect(dipole.I, current)
            parent.connect(dipole.V_in, node.inwards, "V")

        return node


class Ground(System):

    def setup(self, V=0.0):
        self.add_inward("V", V, unit="V")
        self.add_output(Voltage, "V_out", {"V": V})

    def compute(self):
        self.V_out.V = self.V

class Circuit(System):

    def setup(self, I=1.0):
        ground = self.add_child(Ground("ground", V=0.0))

        R1 = self.add_child(Resistor("R1", R=1.00e3))
        R2 = self.add_child(Resistor("R2", R=0.50e3))
        R3 = self.add_child(Resistor("R3", R=0.25e3))

        # Define nodes
        Node.make("n1",
            parent=self,
            pulling={"I_in0": "source"},
            outgoing=[R1, R2],
        )
        Node.make("n2",
            parent=self,
            incoming=[R2],
            outgoing=[R3],
        )

        self.connect(ground.V_out, R1.V_out)
        self.connect(ground.V_out, R3.V_out)

        self.source.I = I  # initial current source

from cosapp.drivers import NonLinearSolver

circuit = Circuit("circuit")

# Set source current
circuit.source.I = 0.1

# Set resistance values
circuit.R1.R = 1000.0
circuit.R2.R = 500.0
circuit.R3.R = 250.0

# Add solver and solve circuit
solver = circuit.add_driver(NonLinearSolver("solver"))
circuit.run_drivers()

# Print out results
print(
    "Resistances:",
    f"  R1, R2, R3 = ({circuit.R1.R:.6}, {circuit.R2.R:.6}, {circuit.R3.R:.6}) Ohm",
    "Currents:",
    f"  I1, I2, I3 = ({circuit.R1.I.I:.4}, {circuit.R2.I.I:.4}, {circuit.R3.I.I:.4}) A",
    "Node voltages:",
    f"  n1: {circuit.n1.V:.4} V",
    f"  n2: {circuit.n2.V:.4} V",
    sep="\n",
)

Resistances:
  R1, R2, R3 = (1000.0, 500.0, 250.0) Ohm
Currents:
  I1, I2, I3 = (0.04286, 0.05714, 0.05714) A
Node voltages:
  n1: 42.86 V
  n2: 14.29 V

Save state for later

from cosapp.utils import get_state, set_state

initial_state = get_state(circuit)

Defining the design case

After purging all drivers, the design case can be defined. First, a numerical solver is attached to the head system. Then, for each design point, a sub-driver RunSingleCase is added to the solver.

from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder

# Clear all previously defined drivers and add new solver
circuit.drivers.clear()
set_state(circuit, initial_state)

solver = circuit.add_driver(NonLinearSolver("solver"))

# Add sub-driver to set boundary conditions on point 1
point1 = solver.add_child(RunSingleCase("point1"))
point1.set_values({
    "source.I": 0.08,
    "ground.V": 0.0,
})

# Same as previous for a second point
point2 = solver.add_child(RunSingleCase("point2"))
point2.set_values({
    "source.I": 0.15,
    "ground.V": 0.0,
})

# Define multi-point design problem
solver.add_unknown(["R1.R", "R2.R"])
point1.add_equation("n2.V == 8")
point2.add_equation("n1.V == 50")

# Add recorder to collect results
solver.add_recorder(
    DataFrameRecorder(
        includes = ["n?.V", "*R", "source.I"],
        excludes = "ground.*",
    )
)

# Solve circuit and show algebraic problem
circuit.run_drivers()

solver.problem

Unknowns [6]
  R1.R = 555.5555555512082
  R2.R = 583.3333333305677
  point1[n1_V] = 26.666666666564883
  point1[n2_V] = 8.0
  point2[n1_V] = 50.0
  point2[n2_V] = 15.00000000001746
Equations [6]
  point1[n2.V == 8] := 0.0
  point1[n1: current balance] := -1.696282003749161e-13
  point1[n2: current balance] := -2.2773449792623524e-14
  point2[n1.V == 50] := 0.0
  point2[n1: current balance] := -9.588163596419008e-13
  point2[n2: current balance] := 1.846994779342026e-13

Export recorded data

data = solver.recorder.export_data()
data = data.drop(["Section", "Status", "Error code"], axis=1)
data

	Reference	R1.R	R2.R	R3.R	n1.V	n2.V	source.I
0	point1	555.555556	583.333333	250.0	26.666667	8.0	0.08
1	point2	555.555556	583.333333	250.0	50.000000	15.0	0.15

Note that resistances R1 and R2 each assume a unique value throughout design points. Resistance R3 is also constant, since it has not been explicitly redefined.

In the next example, we perform a slight variant of the design problem above, with different values of R3, to simulate a potentiometer, say. Additional constraints are:

• At design point 1, we impose R3 = 500 Ohm.

• At design point 2, we seek the value of R3 yielding a potential of 5 V at node 2.

from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder

# Create and initialize new model
circuit = Circuit("circuit")
set_state(circuit, initial_state)

# Add solver
solver = circuit.add_driver(NonLinearSolver("solver"))

# Define design unknwowns, with a maximum relative step
# between iterations to stabilize the resolution
solver.add_unknown(["R1.R", "R2.R"], max_rel_step=0.5)

# Define design point 1
point1 = solver.add_child(RunSingleCase("point1"))

point1.set_values({
    "source.I": 0.08,
    "ground.V": 0,
    "R3.R": 500.0,  # constant value for point 1
})
# Add constraints for point 1
point1.add_equation("n2.V == 8")

# Same as previous for a second point
point2 = solver.add_child(RunSingleCase("point2"))

point2.set_values({
    "source.I": 0.15,
    "ground.V": 0,
})
# Add constraints for point 2
point2.add_equation(["n1.V == 50", "n2.V == 5"])
point2.add_unknown("R3.R")  # local, off-design unknown

# Add recorder to collect results
solver.add_recorder(
    DataFrameRecorder(
        includes = ["n?.V", "*R", "source.I"],
        excludes = "ground.*",
    )
)

circuit.run_drivers()
solver.problem

Unknowns [7]
  R1.R = 438.2605023497829
  R2.R = 1253.0420093991265
  point1[n1_V] = 28.048672150386103
  point1[n2_V] = 8.0
  point2[R3.R] = 139.2268899333792
  point2[n1_V] = 50.0
  point2[n2_V] = 5.0
Equations [7]
  point1[n2.V == 8] := 0.0
  point1[n1: current balance] := -6.938893903907228e-17
  point1[n2: current balance] := 6.245004513516506e-17
  point2[n1.V == 50] := 0.0
  point2[n2.V == 5] := 0.0
  point2[n1: current balance] := -3.0531133177191805e-16
  point2[n2: current balance] := 3.6866343311459104e-14

This time, circuit.R3.R assumes different values for the two design points: its value is fixed at 500 Ohm in point1, but is is computed as an unknown in point2.

data = solver.recorder.export_data()
data = data.drop(["Section", "Status", "Error code"], axis=1)
data

	Reference	R1.R	R2.R	R3.R	n1.V	n2.V	source.I
0	point1	438.260502	1253.042009	500.00000	28.048672	8.0	0.08
1	point2	438.260502	1253.042009	139.22689	50.000000	5.0	0.15

Multi-point design problems involving targets

Targets offer a convenient way to declare a dynamically changeable target value on a variable (either input or output), rather than imposing a hard-coded right-hand-side value, as in:

point1.add_equation("n2.V == 8")

Further detail on targets may be found in a dedicated tutorial.

In single-point design problems, target values can be redefined by simply reassigning the targetted variables prior to solver execution (see tutorial). In multi-point design problems, though, a variable may be assigned different target values in different design points. In the last circuit design problem, for example, potential n2.V is required to be 8 V in point1, and 5 V in point2. Setting point-wise targets can be achieved with method set_init.

Let us reformulate the same design problem, using targets:

from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder

# Create and initialize new model
circuit = Circuit("circuit")
set_state(circuit, initial_state)

# Add solver
solver = circuit.add_driver(NonLinearSolver("solver"))

# Add recorder to collect results
solver.add_recorder(
    DataFrameRecorder(
        includes = ["n?.V", "*R", "source.I"],
        excludes = "ground.*",
    )
)

# Define design point #1
point1 = solver.add_child(RunSingleCase("point1"))

point1.set_values({
    "source.I": 0.08,
    "ground.V": 0.0,
    "R3.R": 500.0,  # constant value for point 1
})

# Define design point #2
point2 = solver.add_child(RunSingleCase("point2"))

point2.set_values({
    "source.I": 0.15,
    "ground.V": 0.0,
})

# Define design unknwowns, with a maximum relative step
# between iterations to stabilize the resolution
solver.add_unknown(["R1.R", "R2.R"], max_rel_step=0.5)

# Set a target on node 2 voltage (pertains to all design points)
solver.add_target("n2.V")

# Add constraints specific to point 2
point2.add_unknown("R3.R")  # local, off-design unknown
point2.add_target("n1.V")   # local target on node 1 voltage

Unknowns [1]
  R3.R = 250.0
Equations [1]
  n1.V == 42.85714285714286 (target)

Once the problem is defined, we can specify point-dependent target values using method set_init:

point1.set_init({
    "n2.V": 8.0,
})
point2.set_init({
    "n1.V": 50.0,
    "n2.V": 5.0,
})

circuit.run_drivers()

solver.problem

Unknowns [7]
  R1.R = 438.2605023497829
  R2.R = 1253.0420093991265
  point1[n1_V] = 28.048672150386103
  point1[n2_V] = 8.0
  point2[R3.R] = 139.2268899333792
  point2[n1_V] = 50.0
  point2[n2_V] = 5.0
Equations [7]
  point1[n1: current balance] := -6.938893903907228e-17
  point1[n2: current balance] := 6.245004513516506e-17
  point1[n2.V == 8.0] := 0.0
  point2[n1: current balance] := -3.0531133177191805e-16
  point2[n2: current balance] := 3.6866343311459104e-14
  point2[n1.V == 50.0] := 0.0
  point2[n2.V == 5.0] := 0.0

Target values can now be redefined interactively, by simply changing initial values:

point1.set_init({
    "n2.V": 12.0,
})
point2.set_init({
    "n1.V": 62.0,
    "n2.V": 7.5,
})

circuit.run_drivers()

solver.problem

Unknowns [7]
  R1.R = 646.2866761049455
  R2.R = 1008.0022442448729
  point1[n1_V] = 36.19205386187696
  point1[n2_V] = 12.0
  point2[R3.R] = 138.71590517131278
  point2[n1_V] = 62.0
  point2[n2_V] = 7.5
Equations [7]
  point1[n1: current balance] := -1.3877787807814457e-17
  point1[n2: current balance] := 6.938893903907228e-18
  point1[n2.V == 12.0] := 0.0
  point2[n1: current balance] := 0.0
  point2[n2: current balance] := -6.938893903907228e-18
  point2[n1.V == 62.0] := 0.0
  point2[n2.V == 7.5] := 0.0