Simple circuit¶
Preliminary note:
This case is taken from OpenMDAO. OpenMDAO is an open-source computing platform for systems analysis and multidisciplinary optimization developed by the NASA. Its philosophy shares some of the goals with CoSApp. So this example is also an opportunity to compare both software.
OpenMDAO is licensed under Apache License.
Case description¶
This is a simple circuit having a source providing the circuit with a constant intensity. The resolution is driving the potential V at node1 and node2 so that the current flowing in equals the one flowing out.
Creating elementary bricks¶
First, we define voltage and current ports. For convenience, we also define an abstract model of directional dipole, to be used as base class for resitors and diodes.
[1]:
import abc
from cosapp.base import System, Port
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
Next, we define the following elements:
Resistors
Diodes
Electric nodes
Source
Ground
[2]:
from __future__ import annotations
from cosapp.base import System
import math
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 Diode(Dipole):
"""Regularized diode model
The current intensity flowing through the diode is calculated based on
$ I = I_s \\exp \\left( \\dfrac{V_{in} - V_{out}}{V_t} - 1 \\right) $
"""
def setup(self):
super().setup()
self.add_inward('Is', 1e-15, desc='Saturation current in Amps')
self.add_inward('Vt', 0.025875, desc='Thermal voltage in Volts')
def compute_I(self):
"""Regularized diode model"""
self.I.I = self.Is * math.exp(self.deltaV / self.Vt - 1.)
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, parent, 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
Building the circuit¶
The circuit is built by reusing Node and Resistor elements customized thanks to keywords arguments that will be processed in the setup method of those classes.
[3]:
class Circuit(System):
def setup(self):
R1 = self.add_child(Resistor('R1', R=100.), pulling={'V_out': 'Vg'})
R2 = self.add_child(Resistor('R2', R=10000.))
D1 = self.add_child(Diode('D1'), 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=[D1],
)
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)
Solving the problem¶
[4]:
from cosapp.drivers import NonLinearSolver
p.add_driver(NonLinearSolver('solver'))
p.run_drivers()
print(
f"{p.circuit.n1.V = }",
f"{p.circuit.n2.V = }",
f"{p.circuit.R1.I = }",
f"{p.circuit.R2.I = }",
f"{p.circuit.D1.I = }",
sep="\n",
)
print(f"\nSanity check: {p.circuit.R1.I.I + p.circuit.D1.I.I = } ({p.source.I} expected)")
p.circuit.n1.V = 9.908302820430912
p.circuit.n2.V = 0.7385848635220567
p.circuit.R1.I = Intensity: {'I': 0.09908302820430911}
p.circuit.R2.I = Intensity: {'I': 0.0009169717956908855}
p.circuit.D1.I = Intensity: {'I': 0.0009169717956908876}
Sanity check: p.circuit.R1.I.I + p.circuit.D1.I.I = 0.1 (0.1 expected)
Show mathematical problem:
[5]:
p.drivers['solver'].problem
[5]:
Unknowns
circuit.n1_V = 9.908302820430912
circuit.n2_V = 0.7385848635220567
Equations
circuit.n1: current balance := 0.0
circuit.n2: current balance := -2.0599841277224584e-18