Stochastic Modelling in fmdtools
This notebook covers the basics of stochastic modelling in fmdtools. Stochastic models are models which run non-determinsitically, meaning that simulating the model at the same inputs results in different outputs because of randomness or uncertainty in the model behavior. Thus, to determine the distribution of outcomes which may result from a simulation, a stochastic model must be run a number of times. This notebook will cover: - How to construct a stochastic model by adding stochastic
states to function blocks and incorporating them in function behavior - How to simulate single scenarios and distributions of scenarios in propagate methods - How to visualize and analyze the results of stochastic model simulations
[1]:
from fmdtools.sim.sample import ParameterSample, ParameterDomain
import fmdtools.sim.propagate as propagate
from fmdtools.analyze.graph import FunctionArchitectureGraph
import inspect # we will be using inspect to look at
This notebook uses Pump model in pump_stochastic.py, which is an adaptation of the original ex_pump.py model with stochastic states added.
[2]:
from pump_stochastic import Pump
Below is the structure, with the same functions/flows as before.
[3]:
mg = FunctionArchitectureGraph(Pump())
fig, ax = mg.draw()
Model Setup
This model has been augmented stochastic states and behaviors to enable stochastic simulation using the Rand class.
See rand documentation: https://nasa.github.io/fmdtools/docs/fmdtools.define.html#fmdtools.define.container.rand.Rand
This is what a Rand object looks like:
[4]:
from fmdtools.define.container.rand import Rand
r = Rand()
r
[4]:
Rand(rng=Generator(PCG64) at 0x22121CCE880, probs=[], probdens=1.0, seed=42, run_stochastic=False)
Below, ImportEERand container is used in the ImportEE function to add two random states (via ImportEERandState): - effstate, the quality of the input voltage (i.e. large fluctuations from a step change in power) - grid_noise, which is meant to be fluctuations in voltage from the power source (i.e., more ordinary noise)
This is then reflected in the behavior, which is defined in the Rand.set_rand_state method.
[5]:
help(Rand.set_rand_state)
Help on function set_rand_state in module fmdtools.define.container.rand:
set_rand_state(self, statename, methodname, *args)
Update the given random state with a given method and arguments.
(if in run_stochastic mode)
Parameters
----------
statename : str
name of the random state defined
methodname :
str name of the numpy method to call in the rng
*args : args
arguments for the numpy method
In this method,
statenameis the name of the random state, whilemethodnamecorresponds to the name of a random distribution to update the state from. These distributions are pulled from a numpy random number generator (see numpy documentation) and include most one might need including normal, beta, uniform, etc.*argsrefers to the corresponding arguments to call the given method with, which should be pulled from the corresponding documentation. Below, we have points where this is called.
Belo we show how this is used in the model:
[6]:
from pump_stochastic import ImportEERand, ImportEERandState, ImportEE
[7]:
print(inspect.getsource(ImportEERand))
print(inspect.getsource(ImportEERandState))
class ImportEERand(Rand, copy_default=True):
s: ImportEERandState = ImportEERandState()
class ImportEERandState(State):
effstate: float = 1.0
grid_noise: float = 1.0
In this setup, the effstate and grid_noise states are to be set by the rng in ImportEERand. Below shows how this is embodied in behavior:
[8]:
print(inspect.getsource(ImportEE.behavior))
def behavior(self, time):
if self.m.has_fault('no_v'):
self.r.s.effstate = 0.0 # an open circuit means no voltage is exported
elif self.m.has_fault('inf_v'):
self.r.s.effstate = 100.0 # a voltage spike means voltage is much higher
else:
if time > self.t.time:
self.r.set_rand_state('effstate', 'triangular', 0.9, 1, 1.1)
if time > self.t.time:
self.r.set_rand_state('grid_noise', 'normal',
1, 0.05*(2+np.sin(np.pi/2*time)))
self.ee_out.s.voltage = self.r.s.grid_noise*self.r.s.effstate * 500
As shown, set_rand_state is used:
First, to pull
effstatefrom a triangular distribution with minimum 0.9, mode 1.0, and max 1.1Parameters are provided here: https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.triangular.html#numpy.random.Generator.triangular
Second, to pull
grid_noisefrom a normal distribution centered on 1 with a standard deviation that varies over time according to a sine wave.Parameters are provided here: https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.normal.html#numpy.random.Generator.normal
Move Water Example
The Move Water function has a similar setup, except the stochasticity is auto-updated at each timestep. This is done by providing an eff_update class varaible to the RandState.
[9]:
from pump_stochastic import MoveWat, MoveWatRandState, MoveWatRand
[10]:
print(inspect.getsource(MoveWatRandState))
class MoveWatRandState(State):
eff: float = 1.0
eff_update = ('normal', (1.0, 0.2))
In this case, the normal distribution will be drawn from with mean 1.0 and standard deviation 0.2.
The advantage of this is that simple stochastic behaviors do not need to be defined in the behavior method, however, it is less flexible than the previous approach, since one is limited to always drawing from the same distribution at each time-step. The corresponding behavior method (below) thus has no set_rand_state call, since the state eff is automaticall updated outside the behavior definition.
[11]:
print(inspect.getsource(MoveWat.behavior))
def behavior(self, time):
self.s.eff = self.r.s.eff
super().behavior(time)
if time > self.t.time:
self.s.inc(total_flow=self.wat_out.s.flowrate)
Model Considerations
To get a distribution of stochastic behavior, the simulation must be run over a number of different random seeds. This seed parameter ensures that each simulation in the same random thread will produce the same results, and that results can be replicated from stochastic simulation by running the model with the same seed.
By default, the seed for the Rand in each of the functions is set to match the Model Rand, meaning that passing a sead at the Model level sets it in all the containing functions also. This can be set by passing an r parameter to the Model.
[12]:
mdl1 = Pump(r={'seed':23})
mdl1.fxns['move_water'].r
[12]:
MoveWatRand(rng=Generator(PCG64) at 0x22121C95620, probs=[], probdens=1.0, seed=23, run_stochastic=False, s=MoveWatRandState(eff=1.0))
[13]:
mdl2 = Pump(r={'seed':10})
mdl2.fxns['move_water'].r
[13]:
MoveWatRand(rng=Generator(PCG64) at 0x22121C95A80, probs=[], probdens=1.0, seed=10, run_stochastic=False, s=MoveWatRandState(eff=1.0))
Simulation and Analysis
With this model set up, we can now simulate it using the methods in propagate.
First, let’s simulate it in the nominal scenario:
[14]:
mdl = Pump()
endclass, mdlhist = propagate.nominal(mdl)
fig, axs = mdlhist.plot_line('fxns.import_ee.r.s.effstate',
'fxns.import_ee.r.s.grid_noise',
'fxns.import_signal.r.s.sig_noise',
'fxns.move_water.r.s.eff')
As shown, even with all of these random states, it doesn’t appear that the behavior actually changes stochastically over time. This is because, by default, the model is set to run deterministically, meaning the stochastic states take their default values. To run stochastically, we use the option run_stochastic=True.
[15]:
mdl = Pump()
endclass, mdlhist = propagate.nominal(mdl, run_stochastic=True)
[16]:
fig, axs = mdlhist.plot_line('fxns.import_ee.r.s.effstate',
'fxns.import_ee.r.s.grid_noise',
'fxns.import_signal.r.s.sig_noise',
'fxns.move_water.r.s.eff')
As shown, this is more what we would expect from a random . Note that this simulation comes from the default model seed, and will this will always be the same. To get a different seed, we can pass new modelparams as an argument to propagate.nominal
[17]:
mdl = Pump()
endclass, mdlhist = propagate.nominal(mdl, run_stochastic=True,
mdl_kwargs=dict(r={'seed':110}))
fig, axs = mdlhist.plot_line('fxns.import_ee.r.s.effstate',
'fxns.import_ee.r.s.grid_noise',
'fxns.import_signal.r.s.sig_noise',
'fxns.move_water.r.s.eff')
We can further simulate faults as we would before.
[18]:
mdl = Pump()
endclass, mdlhist = propagate.one_fault(mdl, 'export_water','block', time=20,
run_stochastic=True, mdl_kwargs=dict(r={'seed':110}))
fig = mdlhist.plot_line('fxns.move_water.r.s.eff',
'fxns.move_water.s.total_flow',
'flows.wat_2.s.flowrate',
'flows.wat_2.s.pressure',
time_slice=[20], legend_loc=False)
As shown, the stochastic states still simulated over time after the fault.
One thing to watch in stochastic models is that off-nominal behavior may be more difficult to track, since stochastic behavior may vary without necessarily being adverse behavior. This can be a problem because many of the visualizations (e.g., in graph) are based on finding differences between faulty and nominal models to visualize degradation, these methods may be less useful/reliable.
[19]:
fig = mdlhist.plot_line('fxns.import_ee.r.s.effstate', 'fxns.import_ee.r.s.grid_noise',
'fxns.import_signal.r.s.sig_noise', 'fxns.move_water.r.s.eff', legend_loc=False)
In the blockage scenario, there is no effect on the underlying stochastic states, since nothing was set up for this in the behavior. In the no_sig fault, on the other hand, we defined the signal noise to go to zero, as shown:
[20]:
mdl = Pump()
endclass, mdlhist = propagate.one_fault(mdl, 'import_signal', 'no_sig', time=20,
run_stochastic=True, mdl_kwargs=dict(r={'seed':110}))
fig = mdlhist.plot_line('fxns.import_ee.r.s.effstate',
'fxns.import_ee.r.s.grid_noise',
'fxns.import_signal.r.s.sig_noise',
'fxns.move_water.r.s.eff',
legend_loc=False,
time_slice=20)
Because stochastic models are non-deterministic, we are often interested not in the results of a single thread, but of the distribution of outcomes that might occur. To perform this kind of assessment, we can use a ParameterSample to instantiate the model with a number of different seeds.
[21]:
help(ParameterSample.add_variable_replicates)
Help on function add_variable_replicates in module fmdtools.sim.sample:
add_variable_replicates(self, x_combos, replicates=1, seed_comb='shared', name='var', weight=1.0)
Add replicates for the given cominations of variables.
Parameters
----------
x_combos : list
List of combinations of variable values. [x_1, x_2...]. If empty, the
default values of x are used.
replicates : int, optional
Number of replicates to add of the variable. The default is 1.
seed_comb : str, optional
How to combine seeds ('shared' or 'independent'). 'shared' uses the same
seed for each value while 'independent' uses different seeds over all
variable values. The default is 'shared'.
name : str, optional
Name prefix for the set of replicates. The default is 'var'.
weight : float, optional
Total weight for the replicates (i.e., total probability to divide between
samples). The default is 1.0.
Examples
--------
>>> ex_ps = ParameterSample(expd)
>>> ex_ps.add_variable_replicates([[1,1], [2,2]])
>>> ex_ps
ParameterSample of scenarios:
- rep0_var_0
- rep0_var_1
>>> scen0 =ex_ps.scenarios()[0]
>>> scen0.prob
0.5
>>> scen0.p['x']
1
>>> scen0.p['y']
1
>>> scen0.inputparams
{0: 1, 1: 1}
>>> scen1 =ex_ps.scenarios()[1]
>>> scen1.prob
0.5
>>> scen1.p['x']
2
>>> scen1.p['y']
2
>>> scen1.inputparams
{0: 2, 1: 2}
The replicates option can also be used in ParameterSample methods to simultaneously change the input parameters and seeds of the model. Below, we create an sample to simulate the model 100 times.
[22]:
ps = ParameterSample()
ps.add_variable_replicates([], replicates=100)
[23]:
endclasses, mdlhists = propagate.parameter_sample(mdl, ps, run_stochastic=True)
SCENARIOS COMPLETE: 100%|██████████| 100/100 [00:04<00:00, 21.05it/s]
To evaluate this behavior over time, we can then use History.plot_line, to plot the distribution of behaviors over time. History.plot_line has a number of different options for visualization of distributions. For example, below we plot model states as individual lines:
[24]:
fig, ax = mdlhists.plot_line('fxns.move_water.r.s.eff', 'fxns.move_water.s.total_flow',
'flows.wat_2.s.flowrate', 'flows.wat_2.s.pressure',
'fxns.import_ee.r.s.effstate', 'fxns.import_ee.r.s.grid_noise',
'flows.ee_1.s.voltage', 'flows.sig_1.s.power',
color='blue', comp_groups={}, alpha=0.1, legend_loc=False)
Or as percentiles:
[25]:
fig = mdlhist.plot_line('fxns.move_water.r.s.eff', 'fxns.move_water.s.total_flow',
'flows.wat_2.s.flowrate', 'flows.wat_2.s.pressure',
'fxns.import_ee.r.s.effstate', 'fxns.import_ee.r.s.grid_noise',
'flows.ee_1.s.voltage', 'flows.sig_1.s.power',
color='blue', comp_groups={}, aggregation='percentile')
Or as a mean with confidence interval:
[26]:
fig, ax = mdlhists.plot_line('fxns.move_water.r.s.eff', 'fxns.move_water.s.total_flow',
'flows.wat_2.s.flowrate', 'flows.wat_2.s.pressure',
'fxns.import_ee.r.s.effstate', 'fxns.import_ee.r.s.grid_noise',
'flows.ee_1.s.voltage', 'flows.sig_1.s.power',
color='blue', comp_groups={}, aggregation='mean_ci')
We can also compare stochastic output over a set of scenarios using plot.hist with the comp_groups parameter, which places the scenarios in different groups.
This can be done by first creating an ParameterDomain with the variable we wish to compare. The pump still has one real parameter–the fault delay.
[27]:
from pump_stochastic import PumpParam
pd = ParameterDomain(PumpParam)
pd.add_variable("delay")
pd
[27]:
ParameterDomain with:
- variables: {'delay': (0, 100)}
- constants: {}
- parameter_initializer: PumpParam
This can then be sampled with a ParameterSample:
[28]:
ps = ParameterSample(pd)
ps.add_variable_replicates([[1]], replicates=100, name="delay1")
ps.add_variable_replicates([[10]], replicates=100, name="delay10")
ps
[28]:
ParameterSample of scenarios:
- rep0_delay1_0
- rep1_delay1_1
- rep2_delay1_2
- rep3_delay1_3
- rep4_delay1_4
- rep5_delay1_5
- rep6_delay1_6
- rep7_delay1_7
- rep8_delay1_8
- rep9_delay1_9
- ... (200 total)
Since this delay only shows up in blockage fault modes, to compare the behaviors, we will first simulate it in a nested approach with only the blockage fault added.
[29]:
import multiprocessing as mp
[30]:
faultdomains = {"fd": (("fault", "export_water", "block"), {})}
faultsamples = {"fs": (("fault_phases", "fd"), {})}
endclasses, mdlhists, apps=propagate.nested_sample(mdl, ps, run_stochastic=True, faultdomains=faultdomains, faultsamples=faultsamples)
NESTED SCENARIOS COMPLETE: 0%| | 1/200 [00:00<00:44, 4.45it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 2%|▏ | 4/200 [00:00<00:44, 4.38it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 3%|▎ | 6/200 [00:01<00:43, 4.43it/s]
Faults found during the nominal run {'import_ee': ['no_v']}
NESTED SCENARIOS COMPLETE: 4%|▍ | 9/200 [00:02<00:43, 4.35it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 6%|▌ | 11/200 [00:02<00:43, 4.34it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 6%|▋ | 13/200 [00:02<00:42, 4.43it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 8%|▊ | 15/200 [00:03<00:40, 4.53it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 8%|▊ | 17/200 [00:03<00:40, 4.54it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 10%|▉ | 19/200 [00:04<00:40, 4.46it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 10%|█ | 21/200 [00:04<00:41, 4.36it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 12%|█▏ | 23/200 [00:05<00:40, 4.37it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 12%|█▎ | 25/200 [00:05<00:40, 4.35it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 14%|█▎ | 27/200 [00:06<00:40, 4.32it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 14%|█▍ | 29/200 [00:06<00:38, 4.48it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 16%|█▌ | 31/200 [00:07<00:38, 4.42it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 16%|█▋ | 33/200 [00:07<00:37, 4.45it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'import_ee': ['no_v']}
NESTED SCENARIOS COMPLETE: 18%|█▊ | 35/200 [00:07<00:37, 4.37it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 18%|█▊ | 37/200 [00:08<00:36, 4.48it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 20%|█▉ | 39/200 [00:08<00:36, 4.37it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 20%|██ | 41/200 [00:09<00:35, 4.52it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 22%|██▏ | 44/200 [00:09<00:33, 4.62it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 23%|██▎ | 46/200 [00:10<00:33, 4.57it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 24%|██▍ | 48/200 [00:10<00:33, 4.59it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 26%|██▌ | 52/200 [00:11<00:33, 4.43it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 28%|██▊ | 55/200 [00:12<00:33, 4.37it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 28%|██▊ | 56/200 [00:12<00:33, 4.32it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 29%|██▉ | 58/200 [00:13<00:31, 4.47it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 30%|███ | 61/200 [00:13<00:30, 4.58it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 32%|███▏ | 64/200 [00:14<00:29, 4.62it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 33%|███▎ | 66/200 [00:14<00:29, 4.51it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 34%|███▍ | 68/200 [00:15<00:29, 4.54it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 35%|███▌ | 70/200 [00:15<00:28, 4.49it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 36%|███▌ | 72/200 [00:16<00:28, 4.56it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 37%|███▋ | 74/200 [00:16<00:27, 4.55it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 40%|████ | 80/200 [00:17<00:27, 4.44it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 42%|████▏ | 84/200 [00:18<00:26, 4.35it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 43%|████▎ | 86/200 [00:19<00:26, 4.36it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 44%|████▍ | 88/200 [00:19<00:25, 4.40it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 46%|████▌ | 91/200 [00:20<00:24, 4.37it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 46%|████▋ | 93/200 [00:20<00:23, 4.50it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 48%|████▊ | 96/200 [00:21<00:22, 4.57it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 49%|████▉ | 98/200 [00:22<00:22, 4.51it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 50%|████▉ | 99/200 [00:22<00:22, 4.49it/s]
Faults found during the nominal run {'move_water': ['mech_break']}
Faults found during the nominal run {'move_water': ['mech_break']}
NESTED SCENARIOS COMPLETE: 53%|█████▎ | 106/200 [00:23<00:20, 4.50it/s]
Faults found during the nominal run {'import_ee': ['no_v']}
NESTED SCENARIOS COMPLETE: 67%|██████▋ | 134/200 [00:30<00:14, 4.41it/s]
Faults found during the nominal run {'import_ee': ['no_v']}
NESTED SCENARIOS COMPLETE: 84%|████████▎ | 167/200 [00:37<00:07, 4.37it/s]
Faults found during the nominal run {'import_ee': ['no_v']}
Faults found during the nominal run {'import_ee': ['no_v']}
NESTED SCENARIOS COMPLETE: 100%|██████████| 200/200 [00:45<00:00, 4.38it/s]
Faults found during the nominal run {'import_ee': ['no_v']}
Next, we get just the scenarios of interest:
[31]:
comp_mdlhists = mdlhists.get_scens('export_water_block_t27p0')
Finally, we create some comparison groups to group the results by:
[32]:
comp_groups = {'delay_1': ps.get_scens(p={'delay': 1}), 'delay_10': ps.get_scens(p={'delay': 10})}
These are the resulting behaviors:
[33]:
figax = comp_mdlhists.plot_line('fxns.move_water.r.s.eff', 'fxns.move_water.s.total_flow',
'flows.wat_2.s.flowrate', 'flows.wat_2.s.pressure',
'fxns.import_ee.r.s.effstate', 'fxns.import_ee.r.s.grid_noise',
'flows.ee_1.s.voltage', 'flows.sig_1.s.power',
comp_groups=comp_groups, aggregation='percentile', time_slice=27)
There a few interesting aspects of this simulation: - First, some of the delay_1 design simulations have a zero flowrate before the fault is injected. This is because the reduced delay means that nominal behavior can easily cause the fault–all that nees to happen is for the pressure to go over 15 for two time-steps, which is entirely forseeable with all the variability in the simulation. This can also be seen in the “total_flow” graph. - Second, the delay_10 gives the nominal behavior
until the fault time, when pressure increases dramatically. This is also observed in the delay_1 design, but with a much smaller delay (1 instead of 1).
Finally, running multiple simulations relies on the NominalApproach, it can be combined fairly easily with existing methods used in conjunction with NominalApproach. (e.g., tabulate.NestedComparison).
[34]:
import numpy as np
from fmdtools.analyze.tabulate import NestedComparison
comp = NestedComparison(endclasses, ps, ["p.delay"], apps, ["fault"], metrics=['cost'], ci_metrics=['cost'], default_stat=np.mean)
[35]:
comp.as_table()
[35]:
| cost | cost_lb | cost_ub | ||
|---|---|---|---|---|
| 10 | block | 15198.134116 | 14514.478099 | 15881.346779 |
| 1 | block | 13590.756110 | 13119.199763 | 14104.400379 |
[36]:
comp.as_plot("cost", color_factor="fault")
[36]:
(<Figure size 600x400 with 1 Axes>, <Axes: xlabel='p.delay'>)
