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.

simple-circuit

Operating point

Boundary conditions

Design variable

Design equation

Point 1

\(I_{source} = 0.08\) A

R2.R

n2.V == 8 V

Point 2

\(I_{source} = 0.15\) A

R1.R

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.

[1]:

from cosapp.base import System, Port
from cosapp.drivers import NonLinearSolver, RunSingleCase
import numpy as np


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 Resistor(System):

    def setup(self, R = 1.):
        self.add_input(Voltage, 'V_in')
        self.add_input(Voltage, 'V_out')
        self.add_output(Intensity, 'I')

        self.add_inward('R', R, unit='ohm', desc='Resistance in Ohms')


    def compute(self):
        self.I.I = (self.V_in.V - self.V_out.V) / self.R


class Node(System):

    def setup(self, n_in=1, n_out=1):
        self.add_property('n_in', int(n_in))
        self.add_property('n_out', int(n_out))

        if min(self.n_in, self.n_out) < 1:
            raise ValueError("Node needs at least one incoming and one outgoing current")

        for i in range(self.n_in):
            self.add_input(Intensity, f"I_in{i}")
        for i in range(self.n_out):
            self.add_input(Intensity, f"I_out{i}")

        self.add_inward('V', unit = 'V')

        self.add_outward('sum_I_in', 0., desc='Sum of all incoming currents')
        self.add_outward('sum_I_out', 0., desc='Sum of all outgoing currents')

        # Off-design current balance condition
        self.add_unknown('V').add_equation('sum_I_in == sum_I_out', name='current balance')

    def compute(self):
        self.sum_I_in = sum(self[f"I_in{i}.I"] for i in range(self.n_in))
        self.sum_I_out = sum(self[f"I_out{i}.I"] for i in range(self.n_out))


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.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):
        n1 = self.add_child(Node('n1', n_in=1, n_out=2), pulling={'I_in0': 'I_in'})
        n2 = self.add_child(Node('n2'))

        R1 = self.add_child(Resistor('R1', R=1000.), pulling={'V_out': 'Vg'})
        R2 = self.add_child(Resistor('R2', R=500.))
        R3 = self.add_child(Resistor('R3', R=250.), pulling={'V_out': 'Vg'})

        self.connect(R1.V_in, n1.inwards, 'V')
        self.connect(R2.V_in, n1.inwards, 'V')
        self.connect(R1.I, n1.I_out0)
        self.connect(R2.I, n1.I_out1)

        self.connect(R2.V_out, n2.inwards, 'V')
        self.connect(R3.V_in, n2.inwards, 'V')
        self.connect(R2.I, n2.I_in0)
        self.connect(R3.I, n2.I_out0)


p = System('model')
p.add_child(Source('source', I=0.1))
p.add_child(Ground('ground', V=0.0))
p.add_child(Circuit('circuit'))

p.connect(p.source.I_out, p.circuit.I_in)
p.connect(p.ground.V_out, p.circuit.Vg)

solver = p.add_driver(NonLinearSolver('solver'))
p.run_drivers()

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

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

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.

[2]:

from cosapp.recorders import DataFrameRecorder

# Clear all previously defined drivers and add solver
p.drivers.clear()
solver = p.add_driver(NonLinearSolver('solver', tol=1e-9))

solver.add_unknown(['circuit.R1.R', 'circuit.R2.R'])  # design unknowns

# 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 = '*R3*',
    )
)
p.run_drivers()

# Export recorded dataframe
data = solver.recorder.export_data()
assert data.at[0, 'circuit.n2.V'] == 8
assert data.at[1, 'circuit.n1.V'] == 50

data

[2]:
Section Status Error code Reference circuit.R1.R circuit.R2.R circuit.n1.V circuit.n2.V ground.V source.I
0 0 point1 555.555556 583.333333 26.666667 8.0 0 0.08
1 0 point2 555.555556 583.333333 50.000000 15.0 0 0.15

Note that resistances R1 and R2 each assume a unique value throughout design points.

[3]:

print(solver.problem)

Unknowns
  circuit.R1.R = 555.5555555546839
  circuit.R2.R = 583.3333333340137
  point1[circuit.n1_V] = 26.666666666706195
  point1[circuit.n2_V] = 8.0
  point1[circuit.n1.I_out0.I] = 0.048
  point1[circuit.n1.I_out1.I] = 0.032
  point1[circuit.n2.I_in0.I] = 0.032
  point1[circuit.n2.I_out0.I] = 0.032
  point2[circuit.n1_V] = 50.0
  point2[circuit.n2_V] = 14.999999999998831
  point2[circuit.n1.I_out0.I] = 0.09000000000000467
  point2[circuit.n1.I_out1.I] = 0.05999999999999532
  point2[circuit.n2.I_in0.I] = 0.05999999999999532
  point2[circuit.n2.I_out0.I] = 0.05999999999999532
Equations
  point1[circuit.n2.V == 8] := 0.0
  point1[circuit.n1: current balance] := 0.0
  point1[circuit.n2: current balance] := 0.0
  point1[circuit: n1.I_out0.I == R1.I.I (loop)] := -1.4645923362976987e-13
  point1[circuit: n1.I_out1.I == R2.I.I (loop)] := -3.044092755644101e-14
  point1[circuit: n2.I_in0.I == R2.I.I (loop)] := -3.044092755644101e-14
  point1[circuit: n2.I_out0.I == R3.I.I (loop)] := 0.0
  point2[circuit.n1.V == 50] := 0.0
  point2[circuit.n1: current balance] := 0.0
  point2[circuit.n2: current balance] := 0.0
  point2[circuit: n1.I_out0.I == R1.I.I (loop)] := -1.3652967645327863e-13
  point2[circuit: n1.I_out1.I == R2.I.I (loop)] := 6.330352908534564e-14
  point2[circuit: n2.I_in0.I == R2.I.I (loop)] := 6.330352908534564e-14
  point2[circuit: n2.I_out0.I == R3.I.I (loop)] := -6.938893903907228e-18