A Monte Carlo demonstration 
¶
Uncertainty is evaluated assuming Gaussian distributions.
For input data :
value + best has a 15 % probability to be reached
value + worst has a 85 % probabiliy to be reached
At interfaces, value is modified by perturbation when going through the connector. Perturbation value is:
+ best with 15 % probability
+ worst with 85 % probability
Design dispersion¶
[1]:
import logging
logging.getLogger().setLevel(logging.INFO)
We will use the circuit define for the multi-points design demonstration.
[2]:
from cosapp.base import System, Port
from cosapp.drivers import NonLinearSolver, RunSingleCase
from cosapp.recorders import DataFrameRecorder
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')
self.add_outward('deltaV', unit='V') # Not mandatory; could be local to compute method
def compute(self):
self.deltaV = self.V_in.V - self.V_out.V
self.I.I = self.deltaV / 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')
self.add_unknown('V') # Iterative variable
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')
self.add_equation('sum_I_in == sum_I_out', name='V')
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)
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)
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.))
p.add_child(Circuit('circuit'))
p.connect(p.source.I_out, p.circuit.I_in)
p.connect(p.ground.V_out, p.circuit.Vg)
solve = p.add_driver(NonLinearSolver('solve', verbose=True))
p.run_drivers()
Connector source 'circuit.inwards.n1_V' is dimensionless, but target 'R1.V_in.V' has physical unit V.
Connector source 'circuit.inwards.n1_V' is dimensionless, but target 'R2.V_in.V' has physical unit V.
Connector source 'circuit.inwards.n2_V' is dimensionless, but target 'R2.V_out.V' has physical unit V.
Connector source 'circuit.inwards.n2_V' is dimensionless, but target 'R3.V_in.V' has physical unit V.
The following function will carry out the execution of the Monte Carlo simulation.
It has two optional parameters:
draws: number of samples to drawlinear: use approximate linear derivative instead of mathematical resolution
[3]:
from cosapp.utils.distributions import Normal
from cosapp.drivers import MonteCarlo
from cosapp.recorders import DataFrameRecorder
def run_analysis(syst, draws=1000, linear=True):
syst.drivers.clear() # Remove all drivers on the System
solver = syst.add_driver(NonLinearSolver('design', verbose=False)) # Add numerical solver
# Add driver to set boundary conditions on point 1
point1 = solver.add_child(RunSingleCase('pt1'))
# Same as previous for a second point
point2 = solver.add_child(RunSingleCase('pt2'))
point1.set_values({
'source.I': 0.08,
'ground.V': 0,
})
point1.design.add_unknown('circuit.R2.R').add_equation('circuit.n2.V == 8')
point2.set_values({
'source.I': 0.15,
'ground.V': 0,
})
point2.design.add_unknown('circuit.R1.R').add_equation('circuit.n1.V == 50')
syst.run_drivers()
# Montecarlo
syst.drivers.clear()
montecarlo = syst.add_driver(MonteCarlo('mc'))
montecarlo.add_recorder(DataFrameRecorder(includes='circuit.R?.R'))
montecarlo.add_child(solver)
montecarlo.draws = draws
montecarlo.linear = linear
# parameters for uncertainties in the data
R_attr = syst.circuit.R3.inwards.get_details('R')
# Set the distribution around the current value
R_attr.distribution = Normal(best=100, worst=-50)
montecarlo.add_random_variable('circuit.R3.R')
montecarlo.add_response(['circuit.R1.R', 'circuit.R2.R'])
# Computation
syst.run_drivers()
return montecarlo.recorder.export_data()
[4]:
results = run_analysis(p)
[5]:
import plotly.graph_objs as go
from plotly.subplots import make_subplots
# List of fields to plot
fields = [n for n in results.columns if 'circuit' in n]
# Create the figure object
fig = make_subplots(rows=1, cols=3)
fig.layout.title = "Probability"
fig.layout.yaxis.title = 'percent'
for i, field in enumerate(fields):
# Add plot
fig.add_trace(
go.Histogram(
x=results[field],
histnorm='percent',
name=field
),
row = 1,
col = i + 1,
)
fig.get_subplot(1, i + 1).xaxis.title = field
fig