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(...)
(respectivelyadd_equation
,extend
) is a shortcut tocase.offdesign.add_unknown(...)
.Fixing a parameter using
case.set_values({'x': '0.123'})
is mathematically equivalent to the trivial constaintcase.add_unknown('x').add_equation('x == 0.123')
. Thus, the local off-design problem can be viewed as an extension ofset_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.
Operating point |
Boundary conditions |
Design variable |
Design equation |
---|---|---|---|
Point 1 |
\(I_{source} = 0.08\) A |
|
|
Point 2 |
\(I_{source} = 0.15\) A |
|
|
Building the circuit¶
The circuit is built as an assembly of elementary models. The default case is then solved to initialize all variables.
[1]:
from __future__ import annotations
from cosapp.base import Port, System
import math
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 Source(System):
def setup(self, I=0.1):
self.add_inward('I', I, unit='A')
self.add_output(Intensity, 'I_out', {'I': I})
def compute(self):
self.I_out.I = self.I
class Ground(System):
def setup(self, V=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):
R1 = self.add_child(Resistor('R1', R=1.00e3), pulling={'V_out': 'Vg'})
R2 = self.add_child(Resistor('R2', R=0.50e3))
R3 = self.add_child(Resistor('R3', R=0.25e3), pulling={'V_out': 'Vg'})
# Define nodes
Node.make('n1',
parent=self,
pulling={'I_in0': 'I_in'},
outgoing=[R1, R2],
)
Node.make('n2',
parent=self,
incoming=[R2],
outgoing=[R3],
)
[2]:
class AssembledModel(System):
def setup(self, I=0.1):
self.add_child(Source('source', I=I))
self.add_child(Ground('ground', V=0.0))
self.add_child(Circuit('circuit'))
self.connect(self.source.I_out, self.circuit.I_in)
self.connect(self.ground.V_out, self.circuit.Vg)
[3]:
from cosapp.drivers import NonLinearSolver
model = AssembledModel('model')
solver = model.add_driver(NonLinearSolver('solver'))
model.run_drivers()
# Print out results
print(
"Resistances:",
f" R1, R2, R3 = ({model.circuit.R1.R:.6}, {model.circuit.R2.R:.6}, {model.circuit.R3.R:.6}) Ohm",
"Currents:",
f" I1, I2, I3 = ({model.circuit.R1.I.I:.4}, {model.circuit.R2.I.I:.4}, {model.circuit.R3.I.I:.4}) A",
"Node voltages:",
f" n1: {model.circuit.n1.V:.4} V",
f" n2: {model.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¶
[4]:
from cosapp.utils import get_state, set_state
initial_state = get_state(model)
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.
[5]:
from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder
# Clear all previously defined drivers and add new solver
model.drivers.clear()
set_state(model, initial_state)
solver = model.add_driver(NonLinearSolver('solver'))
# Define design unknowns
solver.add_unknown(['circuit.R1.R', 'circuit.R2.R'])
# Add driver to set boundary conditions on point 1
point1 = solver.add_child(RunSingleCase('point1'))
point1.add_equation('circuit.n2.V == 8')
point1.set_values({
'source.I': 0.08,
'ground.V': 0,
})
# Same as previous for a second point
point2 = solver.add_child(RunSingleCase('point2'))
point2.add_equation('circuit.n1.V == 50')
point2.set_values({
'source.I': 0.15,
'ground.V': 0,
})
solver.add_recorder(
DataFrameRecorder(
includes = ['*.n?.V', '*R', 'source.I'],
excludes = 'ground.*',
)
)
model.run_drivers()
solver.problem
[5]:
Unknowns [6]
circuit.R1.R = 555.5555555512082
circuit.R2.R = 583.3333333305679
point1[circuit.n1_V] = 26.66666666656489
point1[circuit.n2_V] = 8.0
point2[circuit.n1_V] = 50.0
point2[circuit.n2_V] = 15.000000000017456
Equations [6]
point1[circuit.n2.V == 8] := 0.0
point1[circuit.n1: current balance] := -1.6964207816272392e-13
point1[circuit.n2: current balance] := -2.2773449792623524e-14
point2[circuit.n1.V == 50] := 0.0
point2[circuit.n1: current balance] := -9.587886040662852e-13
point2[circuit.n2: current balance] := 1.846856001463948e-13
Export recorded data¶
[6]:
data = solver.recorder.export_data()
data = data.drop(['Section', 'Status', 'Error code'], axis=1)
data
[6]:
Reference | circuit.R1.R | circuit.R2.R | circuit.R3.R | circuit.n1.V | circuit.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.
[7]:
from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder
# Create and initialize new model
model = AssembledModel('model')
set_state(model, initial_state)
# Add solver
solver = model.add_driver(NonLinearSolver('solver'))
# Define design unknwowns.
# A maximum relative step between iterations is requested,
# to stabilize the resolution.
solver.add_unknown(['circuit.R1.R', 'circuit.R2.R'], max_rel_step=0.5)
# Add driver to set boundary conditions on point 1
point1 = solver.add_child(RunSingleCase('point1'))
point1.add_equation('circuit.n2.V == 8')
point1.set_values({
'source.I': 0.08,
'ground.V': 0,
'circuit.R3.R': 0.5e3,
})
# Same as previous for a second point
point2 = solver.add_child(RunSingleCase('point2'))
point2.add_unknown('circuit.R3.R') # local, off-design unknown
point2.add_equation([
'circuit.n1.V == 50',
'circuit.n2.V == 5',
])
point2.set_values({
'source.I': 0.15,
'ground.V': 0,
})
solver.add_recorder(
DataFrameRecorder(
includes = ['*.n?.V', '*R', 'source.I'],
excludes = 'ground.*',
)
)
model.run_drivers()
solver.problem
[7]:
Unknowns [7]
circuit.R1.R = 438.26050234978294
circuit.R2.R = 1253.042009399127
point1[circuit.n1_V] = 28.048672150386107
point1[circuit.n2_V] = 8.0
point2[circuit.R3.R] = 139.22688993337923
point2[circuit.n1_V] = 50.0
point2[circuit.n2_V] = 5.0
Equations [7]
point1[circuit.n2.V == 8] := 0.0
point1[circuit.n1: current balance] := -5.551115123125783e-17
point1[circuit.n2: current balance] := 5.898059818321144e-17
point2[circuit.n1.V == 50] := 0.0
point2[circuit.n2.V == 5] := 0.0
point2[circuit.n1: current balance] := -2.7755575615628914e-16
point2[circuit.n2: current balance] := 3.68594044175552e-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
.
[8]:
data = solver.recorder.export_data()
data = data.drop(['Section', 'Status', 'Error code'], axis=1)
data
[8]:
Reference | circuit.R1.R | circuit.R2.R | circuit.R3.R | circuit.n1.V | circuit.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('circuit.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 circuit.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:
[9]:
from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder
# Create and initialize new model
model = AssembledModel('model')
set_state(model, initial_state)
# Add solver
solver = model.add_driver(NonLinearSolver('solver'))
# Define design unknwowns.
solver.add_unknown(['circuit.R1.R', 'circuit.R2.R'], max_rel_step=0.5)
# Set a target on node 2 voltage (will pertain to all design points)
solver.add_target('circuit.n2.V')
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,
'circuit.R3.R': 0.5e3,
})
# Define design point #2
point2 = solver.add_child(RunSingleCase('point2'))
point2.add_unknown('circuit.R3.R') # local, off-design unknown
point2.add_target('circuit.n1.V') # local target on node 1 voltage
point2.set_values({
'source.I': 0.15,
'ground.V': 0,
})
Once the problem is defined, we can specify point-dependent target values using method set_init
:
[10]:
point1.set_init({
'circuit.n2.V': 8.0,
})
point2.set_init({
'circuit.n1.V': 50.0,
'circuit.n2.V': 5.0,
})
model.run_drivers()
solver.problem
[10]:
Unknowns [7]
circuit.R1.R = 438.26050234978294
circuit.R2.R = 1253.042009399127
point1[circuit.n1_V] = 28.048672150386107
point1[circuit.n2_V] = 8.0
point2[circuit.R3.R] = 139.22688993337923
point2[circuit.n1_V] = 50.0
point2[circuit.n2_V] = 5.0
Equations [7]
point1[circuit.n1: current balance] := -5.551115123125783e-17
point1[circuit.n2: current balance] := 5.898059818321144e-17
point1[circuit.n2.V == 8.0] := 0.0
point2[circuit.n1: current balance] := -2.7755575615628914e-16
point2[circuit.n2: current balance] := 3.68594044175552e-14
point2[circuit.n1.V == 50.0] := 0.0
point2[circuit.n2.V == 5.0] := 0.0
Target values can now be redefined interactively, by simply changing initial values:
[11]:
point1.set_init({
'circuit.n2.V': 12.0,
})
point2.set_init({
'circuit.n1.V': 62.0,
'circuit.n2.V': 7.5,
})
model.run_drivers()
solver.problem
[11]:
Unknowns [7]
circuit.R1.R = 646.2866761049455
circuit.R2.R = 1008.0022442448729
point1[circuit.n1_V] = 36.19205386187696
point1[circuit.n2_V] = 12.0
point2[circuit.R3.R] = 138.71590517131278
point2[circuit.n1_V] = 62.0
point2[circuit.n2_V] = 7.5
Equations [7]
point1[circuit.n1: current balance] := -1.3877787807814457e-17
point1[circuit.n2: current balance] := 6.938893903907228e-18
point1[circuit.n2.V == 12.0] := 0.0
point2[circuit.n1: current balance] := 0.0
point2[circuit.n2: current balance] := -6.938893903907228e-18
point2[circuit.n1.V == 62.0] := 0.0
point2[circuit.n2.V == 7.5] := 0.0