! python --versionPython 3.11.14
In this project we apply the established Active Inference (AIF) modeling approach (making use of Bond Graphs, Generalized Filters, and Generalized Coordinates of Motion) to an electric system. The latter is not a complicated system. Also, the generalized coordinates part of the experiment needs more tuning. In this first part, we will make use of Bond Graphs to derive the state space equations that will be used in the following parts.
A few coding conventions:
- Bold variable names usually start with an underscore, i.e. ’_’. That means all structured variables including vectors (i.e. tensors). The focus is on the code equivalents of the mathematics, i.e. not on Python structured variables, e.g. lists.
- Global variable names (i.e. used across cells) are usually indicated by the Python ‘global’ keyword
- External (i.e. true) states and parameters: The ‘^x’ symbol is used for code variables, i.e. ‘v^x’. In other words, superscript ‘x’ is used to imitate superscript ’*’ in code. In documentation the standard ‘\(v^*\)’ is used.
- This notebook follows a standardized framework. This is my own framework, developed over a number of years, and is informed by the CRISP-DM framework, the work of Warren Powell (Princeton), the work of Bert De Vries, the work of Karl Friston as well as that of the larger Active Inference community.
import BondGraphTools as bgt
from BondGraphTools import component_manager as cm
import importlib.metadata as md
import sympy as sp
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from matplotlib.gridspec import GridSpec
from IPython.display import Math, display
print("BondGraphTools:", md.version("BondGraphTools"))
print("sympy:", sp.__version__)BondGraphTools: 0.4.6
sympy: 1.14.01 BUSINESS UNDERSTANDING
The current need is to apply an Active Inference (AIF) modeling approach (making use of Bond Graphs, Generalized Filters, and Generalized Coordinates of Motion) to a simple electric system. In this first part, we will make use of Bond Graphs to derive the state space equations that will be used in the following parts.
Bond graphs provide an energy-based, domain-independent way to derive state-space models that is especially useful for complex, multi-physics engineering systems. They turn the messy bookkeeping of variables and interconnections into a structured graphical procedure for obtaining first-order equations in a systematic way.
Using bond graphs, engineers can:
- Work with a unified representation across mechanical, electrical, thermal and other domains, which simplifies deriving coupled state-space equations for multi-domain systems.
- Exploit causality assignment to automatically identify suitable state variables (typically energy storage variables) and algorithmically generate consistent state-space or differential-algebraic equations.
- Gain physical insight into power flow, storage and dissipation, which helps validate models, anticipate structural properties (e.g. controllability issues), and support modular model reuse and extension.
Generalized filtering provides a variational Bayesian framework for estimating and controlling the hidden states of nonlinear dynamical systems in continuous time. By formulating inference as gradient descent on a free-energy functional, it delivers online state and parameter estimates that can be directly coupled to control laws, without the need for backward passes as in standard smoothing approaches.
When applied to systems whose state equations are derived from bond graphs, generalized filters operate on physically meaningful energy-based state variables, allowing control to respect the underlying power flows and conservation laws encoded in the model. This combination enables a principled closed-loop scheme in which bond-graph state equations generate predictions, generalized filtering updates beliefs about hidden states and disturbances, and control actions are optimized by minimizing variational free energy—effectively unifying estimation and control within a single coherent formalism.
1.1 Bond Graph construction
To derive the state equations we will make use of bondgraph modeling. We consider the displacement of electric charge \(q\) in a capacitor and the magnetic flux in an inductor as state variables. The effort \(e(t)\) will be the voltage across a component and the flow \(f(t)\) will be the current flowing through the component. The following diagrams show a circuit of bondgraph components as well as the associated bondgraph.
from IPython.display import Image, display
display(Image(filename="ElectricSystem-CIRCUIT.png"))
from IPython.display import Image, display
display(Image(filename="ElectricSystem-BG.png"))
The signal \(a_5(t)\) will be for control, \(v_7(t)\) is a disturbance or exogenous force, and \(w_1(t)\) is process noise.
print("Libraries:")
for lib_id, lib_desc in cm.get_library_list():
print(f"- {lib_id}: {lib_desc}")
print("\nElements in each library:")
for lib_id, lib_desc in cm.get_library_list():
print(f"\n[{lib_id}] {lib_desc}")
comps = cm.get_components_list(lib_id)
for cid, desc in comps:
print(f" {cid:10s} {desc}")Libraries:
- elec: Electrical Components
- BioChem: Biochemical Components
- base: Basic Components
Elements in each library:
[elec] Electrical Components
Di Diode
BJT Bipolar Junction Transistor
[BioChem] Biochemical Components
Ce Concentration of Chemical Species
Re Biochemical Reaction
Y Stoichiometry
[base] Basic Components
0 Equal Effort Node
1 Equal Flow Node
R Generalised Linear Resistor
C Generalised Linear Capacitor
I Generalised Linear Inductor
Se Effort Source
Sf Flow Source
TF Linear Transformer
GY Linear Gyrator
SS Source Sensor
PH Port Hamiltonianbg = bgt.new(name="Electric System")Se1_sym, R2_sym, C3_sym, Sf5_sym, I6_sym, Sf7_sym = \
sp.symbols('Se1 R2, C3 Sf5, I6 Sf7') ## names
Se1 = bgt.new(component='Se', name='Se1', value=Se1_sym)
R2 = bgt.new(component='R', name='R2', value=R2_sym)
_1_ = bgt.new(component='1', name='_1_')
C3 = bgt.new(component='C', name='C3', value=C3_sym)
Sf5 = bgt.new(component='Sf', name='Sf5', value=Sf5_sym)
_0_ = bgt.new(component='0', name='_0_')
I6 = bgt.new(component='I', name='I6', value=I6_sym)
Sf7 = bgt.new(component='Sf', name='Sf7', value=Sf7_sym)bgt.add(bg, Se1, R2, _1_, C3, Sf5, _0_, I6, Sf7)bgt.connect(Se1, _1_)
bgt.connect(_1_, R2)
bgt.connect(_1_, C3)
bgt.connect(_1_, _0_)
bgt.connect(Sf5, _0_)
bgt.connect(_0_, I6)
bgt.connect(_0_, Sf7)bgt.draw(bg)
plt.show()
1.2 Equations
print("State vars:", bg.state_vars); print()
for i,rel in enumerate(bg.constitutive_relations):
print(i, rel)State vars: {'x_0': (C: C3, 'q_0'), 'x_1': (I: I6, 'p_0')}
0 Sf5 + Sf7 + dx_0 - x_1/I6
1 -R2*Sf5 - R2*Sf7 - Se1 + dx_1 + R2*x_1/I6 + x_0/C3print("State vars:", bg.state_vars); print()
for i,rel in enumerate(bg.constitutive_relations):
print(i, rel)
print()
## internal names (from BondGraphTools)
x0, x1 = sp.symbols('x_0 x_1')
dx0, dx1 = sp.symbols('dx_0 dx_1')
Se1 = sp.symbols('Se1')
Sf5 = sp.symbols('Sf5')
Sf7 = sp.symbols('Sf7')
## external preferred names
q3, p6 = sp.symbols('q3 p6')
dq3, dp6 = sp.symbols('dq3 dp6')
w1 = sp.symbols('w1(t)')
a5 = sp.symbols('a5(t)')
v7 = sp.symbols('v7(t)')
subs_map = {
x0: q3,
x1: p6,
dx0: dq3,
dx1: dp6,
Se1: w1,
Sf5: a5,
Sf7: v7
}
relabeled = [rel.subs(subs_map) for rel in bg.constitutive_relations]
for i,r in enumerate(relabeled):
print(i, r)State vars: {'x_0': (C: C3, 'q_0'), 'x_1': (I: I6, 'p_0')}
0 Sf5 + Sf7 + dx_0 - x_1/I6
1 -R2*Sf5 - R2*Sf7 - Se1 + dx_1 + R2*x_1/I6 + x_0/C3
0 a5(t) + dq3 + v7(t) - p6/I6
1 -R2*a5(t) - R2*v7(t) + dp6 - w1(t) + R2*p6/I6 + q3/C3def simplify_for_state_space(sol):
sol_terms = sp.expand(sol); display(Math(sp.latex(sol_terms))) ## needed even when next statement is used
terms = sol_terms.as_ordered_terms(); display(Math(sp.latex(terms))) ## list of terms in the sum
simplified_terms = []
for term in terms:
## Get numerator and denominator of this term
num, den = sp.fraction(term)
## Factor numerator and denominator
num_f = sp.factor(num)
den_f = sp.factor(den)
## Cancel common factors within this term
term_simpl = sp.cancel(num_f / den_f)
simplified_terms.append(term_simpl)
## Recombine simplified terms into a single expression
expr_simpl = sp.Add(*simplified_terms); display(Math(sp.latex(expr_simpl)))
return expr_simpl
def simplify_final_terms(ex):
terms2 = ex.as_ordered_terms(); display(Math(sp.latex(terms2)))
simplified_terms2 = []
for term2 in terms2:
term_simpl2 = sp.simplify(term2)
simplified_terms2.append(term_simpl2)
expr_simpl2 = sp.Add(*simplified_terms2); display(Math(sp.latex(expr_simpl2)))
return expr_simpl2## setup symbols for solver
R2,C3,I6 = sp.symbols('R2 C3 I6')
t = sp.symbols('t')
w1 = sp.Function('w1')
a5 = sp.Function('a5')
v7 = sp.Function('v7')## solve for dq3
## 0:
expr = a5(t) + dq3 + v7(t) - p6/I6
eq = sp.Eq(expr, 0)
sol = sp.solve(eq, dq3)[0]
ex = simplify_for_state_space(sol)
## ex = sp.collect(ex, [F(t)]); display(Math(sp.latex(ex)))
ex = simplify_final_terms(ex)
## ex = ex.as_ordered_terms()
## sp.factor(sp.denom(ex[0]))\(\displaystyle - a_{5}{\left(t \right)} - v_{7}{\left(t \right)} + \frac{p_{6}}{I_{6}}\)
\(\displaystyle \left[ - a_{5}{\left(t \right)}, \ - v_{7}{\left(t \right)}, \ \frac{p_{6}}{I_{6}}\right]\)
\(\displaystyle - a_{5}{\left(t \right)} - v_{7}{\left(t \right)} + \frac{p_{6}}{I_{6}}\)
\(\displaystyle \left[ - a_{5}{\left(t \right)}, \ - v_{7}{\left(t \right)}, \ \frac{p_{6}}{I_{6}}\right]\)
\(\displaystyle - a_{5}{\left(t \right)} - v_{7}{\left(t \right)} + \frac{p_{6}}{I_{6}}\)
\[ \dot{q}_3(t) = + \frac{1}{I_6} p_{6}(t) - v_7(t) - a_5(t) \]
## solve for dp6
## 1:
expr = -R2*a5(t) - R2*v7(t) + dp6 - w1(t) + R2*p6/I6 + q3/C3
eq = sp.Eq(expr, 0)
sol = sp.solve(eq, dp6)[0]
ex = simplify_for_state_space(sol)
## ex = sp.collect(ex, [F(t)]); display(Math(sp.latex(ex)))
ex = simplify_final_terms(ex)
## ex = ex.as_ordered_terms()
## sp.factor(sp.denom(ex[0]))\(\displaystyle R_{2} a_{5}{\left(t \right)} + R_{2} v_{7}{\left(t \right)} + w_{1}{\left(t \right)} - \frac{R_{2} p_{6}}{I_{6}} - \frac{q_{3}}{C_{3}}\)
\(\displaystyle \left[ R_{2} a_{5}{\left(t \right)}, \ R_{2} v_{7}{\left(t \right)}, \ w_{1}{\left(t \right)}, \ - \frac{R_{2} p_{6}}{I_{6}}, \ - \frac{q_{3}}{C_{3}}\right]\)
\(\displaystyle R_{2} a_{5}{\left(t \right)} + R_{2} v_{7}{\left(t \right)} + w_{1}{\left(t \right)} - \frac{R_{2} p_{6}}{I_{6}} - \frac{q_{3}}{C_{3}}\)
\(\displaystyle \left[ R_{2} a_{5}{\left(t \right)}, \ R_{2} v_{7}{\left(t \right)}, \ w_{1}{\left(t \right)}, \ - \frac{R_{2} p_{6}}{I_{6}}, \ - \frac{q_{3}}{C_{3}}\right]\)
\(\displaystyle R_{2} a_{5}{\left(t \right)} + R_{2} v_{7}{\left(t \right)} + w_{1}{\left(t \right)} - \frac{R_{2} p_{6}}{I_{6}} - \frac{q_{3}}{C_{3}}\)
\[ \dot{p}_6(t) = - \frac{1}{C_3} q_{3}(t) - \frac{R_2}{I_6}p_6(t) + R_2v_7(t) + R_2 a_5(t) + w_1(t) \]
Collect all the equations:
\[\begin{align} \dot{q}_3(t) &= + \frac{1}{I_6} p_{6}(t) - v_7(t) - a_5(t) \\ \dot{p}_6(t) &= - \frac{1}{C_3} q_{3}(t) - \frac{R_2}{I_6}p_6(t) + R_2v_7(t) + R_2 a_5(t) + w_1(t) \end{align}\]
1.3 Simulation
ALPHA = 0.6
ENV_WIDTH = 4. ## 8.
AGT_WIDTH = 2.
AGT_LINESTYLES = ['-', '--', '-.', ':'] ## for agent
ORANGE_HUE_BY_SATURATION = [
"#FFA500", ## 100% bright orange
"#FFB347", ## 80% light vivid orange
"#FFD580", ## 60% softer pastel orange
"#FFE5B4" ## 40% light faded orange
]
ORANGE_HUE_BY_BRIGHTNESS = [
"#FFA500", ## 100% pure bright orange
"#FFC04D", ## 80% light orange
"#FFD699", ## 60% pastel orange
"#FFE5CC" ## 40% very pale orange
]
ORANGE_HUE = ORANGE_HUE_BY_BRIGHTNESS
styles = {
#### envir
"a": {
"color": "crimson",
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"xˣ[0]": {
"color": "black",
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"xˣ[1]": {
"color": "darkgrey",
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"vˣ[0]": {
"color": "blue",
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"vˣ[1]": {
"color": "dodgerblue",
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"y[0]": {
"color": ORANGE_HUE[0],
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"y[1]": {
"color": ORANGE_HUE[1],
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"y[2]": {
"color": ORANGE_HUE[2],
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"y[3]": {
"color": ORANGE_HUE[3],
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
"w": {
"color": "yellow",
"linestyle": "-",
"linewidth": ENV_WIDTH,
"alpha": ALPHA
},
#### agent
"μ_x[0]": {
"color": "black",
"linestyle": AGT_LINESTYLES[0],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"μ_x[1]": {
"color": "black",
"linestyle": AGT_LINESTYLES[1],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"v[0]": {
"color": "blue",
"linestyle": AGT_LINESTYLES[0],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"v[1]": {
"color": "blue",
"linestyle": AGT_LINESTYLES[1],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"μ_v[0]": {
"color": "dodgerblue",
"linestyle": AGT_LINESTYLES[0],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"μ_v[1]": {
"color": "dodgerblue",
"linestyle": AGT_LINESTYLES[1],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"μ_y[0]": {
"color": "orange",
"linestyle": AGT_LINESTYLES[0],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"μ_y[1]": {
"color": "orange",
"linestyle": AGT_LINESTYLES[1],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"μ_y[2]": {
"color": "orange",
"linestyle": AGT_LINESTYLES[2],
"linewidth": AGT_WIDTH,
"alpha": ALPHA
},
"F": {
"color": "#9400D3", ## vivid deep violet
"linestyle": "-",
"linewidth": AGT_WIDTH,
"alpha": ALPHA
}
}def plot_states_and_sources(
t_span,
t_states,
state_series,
state_labels,
state_styles,
t_sources,
source_series,
source_labels,
source_styles,
figsize=(8, 5),
):
"""
state_series: list of 1D arrays (same shape as t_states)
state_labels: list of strings, LaTeX-ready
state_styles: list of dicts, e.g. [styles['xˣ[0]'], styles['xˣ[1]'], ...]
source_series: list of 1D arrays (same shape as t_sources)
source_labels: list of strings
source_styles: list of dicts, e.g. [styles['vˣ[0]'], ...]
"""
ylabelx = -0.3
ylabely = 0.3
fig = plt.figure(figsize=figsize)
grid_rows = 2
gs = GridSpec(grid_rows, 1, figure=fig, height_ratios=[1, 1], hspace=0.2)
ax = [fig.add_subplot(gs[i]) for i in range(grid_rows)]
## Common axis formatting
for axis in ax:
axis.grid(True)
axis.spines['top'].set_visible(False)
axis.spines['right'].set_visible(False)
## Top: states
i = 0
ax[i].set_title('states', fontweight='regular', fontsize=12)
for y, lab, st in zip(state_series, state_labels, state_styles):
ax[i].plot(t_states, y, **st, label=lab)
ax[i].set_xticklabels([])
ax[i].yaxis.set_label_coords(ylabelx, ylabely)
ax[i].legend(loc='upper left')
## Bottom: exogvs (exogenous forces)
i = 1
ax[i].set_title('sources', fontweight='regular', fontsize=12)
for y, lab, st in zip(source_series, source_labels, source_styles):
ax[i].step(t_sources, y, where="post", **st, label=lab)
ax[i].legend(loc='upper left')
## X ticks from t_span
x_ticks = np.arange(t_span[0], t_span[1] + 1, 1.0)
ax[i].set_xticks(x_ticks)
ax[i].set_xlabel(r'$\mathrm{Timestep}\ t$', fontweight='bold', fontsize=12)
plt.subplots_adjust(hspace=0.1)
plt.show()
return fig, ax## build values model for simulation
bg = bgt.new(name="Electric System")
Se1_val = 1.
R2_val = 2.
C3_val = 5.
Sf5_val = 1.
I6_val = 3.
Sf7_val = 1.
Se1 = bgt.new(component='Se', name='Se1', value=Se1_val)
R2 = bgt.new(component='R', name='R2', value=R2_val)
_1_ = bgt.new(component='1', name='_1_')
C3 = bgt.new(component='C', name='C3', value=C3_val)
Sf5 = bgt.new(component='Sf', name='Sf5', value=Sf5_val)
_0_ = bgt.new(component='0', name='_0_')
I6 = bgt.new(component='I', name='I6', value=I6_val)
Sf7 = bgt.new(component='Sf', name='Sf7', value=Sf7_val)
bgt.add(bg, Se1, R2, _1_, C3, Sf5, _0_, I6, Sf7)
bgt.connect(Se1, _1_)
bgt.connect(_1_, R2)
bgt.connect(_1_, C3)
bgt.connect(_1_, _0_)
bgt.connect(Sf5, _0_)
bgt.connect(_0_, I6)
bgt.connect(_0_, Sf7)
print(bg.state_vars)
t_span = (0.0, 10.0)
t_eval = np.linspace(t_span[0], t_span[1], 1001)
x0 = [1., 2.] ## initial conditions for states (adjust length if needed){'x_0': (C: C3, 'q_0'), 'x_1': (I: I6, 'p_0')}params = {
## Se1_val = 1.
'R2': R2_val,
'C3': C3_val,
## Sf5_val = 1.
'I6': I6_val,
## Sf7_val = 1.
}
def make_rhs(funcs, params, **func_kwargs):
## left sides must be symbols in latex equations
R2 = params['R2']
C3 = params['C3']
I6 = params['I6']
## unpack per-source theta dicts, default to empty
v7_theta = func_kwargs.get('v7_theta', {})
a5_theta = func_kwargs.get('a5_theta', {})
w1_theta = func_kwargs.get('w1_theta', {})
def rhs(t, x):
q3, p6 = x
## exogenous inputs
v7 = funcs['v7'](t, **v7_theta)
a5 = funcs['a5'](t, **a5_theta)
w1 = funcs['w1'](t, **w1_theta)
# t1_num = 1; t1_den = I6
# t2_num = 1; t2_den = 1
# t3_num = 1; t3_den = 1
dq3 = +(1/I6)*p6 - v7 - a5
# t1_num = 1; t1_den = C3
# t2_num = R2; t2_den = I6
# t3_num = R2; t3_den = 1
# t3_num = R2; t3_den = 1
dp6 = -(1/C3)*q3 - (R2/I6)*p6 + R2*v7 + R2*a5 + w1
return [dq3, dp6]
return rhs
def f_const(t, A=0.0, **kwargs):
t_arr = np.asarray(t)
val = A * np.ones_like(t_arr, dtype=float)
return float(val) if np.isscalar(t) else val
def f_cos(t, A=1.0, θ1=0.0):
t_arr = np.asarray(t, dtype=float)
val = A * np.cos(θ1 * t_arr)
return float(val) if np.isscalar(t) else val
def f_pulses(t, **kwargs):
"""
Sum of rectangular pulses defined in physical time.
kwargs:
'amps' : list/array of amplitudes [A1, A2, ...]
't_ons' : list/array of start times [t_on1, t_on2, ...]
't_offs' : list/array of end times [t_off1, t_off2, ...]
"""
amps = np.asarray(kwargs['amps'], dtype=float)
t_ons = np.asarray(kwargs['t_ons'], dtype=float)
t_offs = np.asarray(kwargs['t_offs'], dtype=float)
t_arr = np.asarray(t, dtype=float)
val = np.zeros_like(t_arr, dtype=float)
for A, t_on, t_off in zip(amps, t_ons, t_offs):
val += A * (np.heaviside(t_arr - t_on, 0.0) - np.heaviside(t_arr - t_off, 0.0))
return float(val) if np.isscalar(t) else val
def f_lines(t, **kwargs):
"""
Piecewise-linear signal defined by points (t_i, y_i).
kwargs:
'y_pts' : list/array of values [y0, y1, ..., yN]
't_pts' : list/array of time points [t0, t1, ..., tN]
For t between t_i and t_{i+1}, returns the straight line
interpolating between (t_i, y_i) and (t_{i+1}, y_{i+1}).
For t < t0 returns y0, for t > tN returns yN.
"""
y_pts = np.asarray(kwargs['y_pts'], dtype=float)
t_pts = np.asarray(kwargs['t_pts'], dtype=float)
if t_pts.ndim != 1 or y_pts.ndim != 1 or t_pts.size != y_pts.size:
raise ValueError("t_pts and y_pts must be 1D arrays of the same length")
t_arr = np.asarray(t, dtype=float)
## use numpy's 1D linear interpolation
val = np.interp(t_arr, t_pts, y_pts)
return float(val) if np.isscalar(t) else val
rng = np.random.default_rng(1234)
## --- Option 1: White Gaussian noise with std σ ---
def f_normal(t, **kwargs):
"""
Gaussian white noise with std = A*sigma, mean = 0.
Works for scalar or array t.
Expects kwargs: 'A', 'sigma'.
"""
A = kwargs.get('A', 1.0)
sigma = kwargs.get('sigma', 1.0)
t_arr = np.asarray(t)
noise = rng.normal(loc=0.0, scale=A * sigma, size=t_arr.shape)
return float(noise) if np.isscalar(t) else noise
## --- Option 2: Uniform noise in [-A, A] ---
def f_uniform(t, **kwargs):
"""
Uniform white noise in [-A, A].
Works for scalar or array t.
Expects kwargs: 'A'.
"""
A = kwargs.get('A', 1.0)
t_arr = np.asarray(t)
u = rng.random(size=t_arr.shape)
noise = (2 * A) * (u - 0.5)
return float(noise) if np.isscalar(t) else noisesources = {
'v7': f_const,
'a5': f_const,
'w1': f_normal,
}
v7_const_A = 10.
a5_const_A = 2.
w1_normal_A = 1
w1_normal_sigma = 1.
rhs = make_rhs(
sources,
params,
v7_theta={'A': v7_const_A},
a5_theta={'A': a5_const_A},
w1_theta={'A': w1_normal_A, 'sigma': w1_normal_sigma},
)
sol = solve_ivp(rhs, t_span, x0, t_eval=t_eval)
t = sol.t
q3 = sol.y[0] ## x_0
p6 = sol.y[1] ## x_1
v7_vals = f_const(t_eval, A=v7_const_A)
a5_vals = f_const(t_eval, A=a5_const_A)
w1_vals = f_normal(t_eval, A=w1_normal_A, sigma=w1_normal_sigma)
fig, ax = plot_states_and_sources(
t_span=t_span,
t_states=t,
state_series=[q3, p6],
state_labels=[
r'$q_{3}(t)$', r'$p_{6}(t)$'
],
state_styles=[
styles['xˣ[0]'], styles['xˣ[1]']
],
t_sources=t_eval,
source_series=[v7_vals, a5_vals, w1_vals],
source_labels = [r'$v_7(t)$', r'$a_5(t)$', r'$w_1(t)$'],
source_styles=[styles['vˣ[0]'], styles['a'], styles['w']],
)
sources = {
'v7': f_cos,
'a5': f_lines,
'w1': f_normal,
}
v7_cos_A = 5.
v7_cos_θ1 = 3.
a5_lines_y_pts = [20, -3, -9, 10]
a5_lines_t_pts = [0., 2.5, 5, 10]
w1_normal_A = 1
w1_normal_sigma = 1.
rhs = make_rhs(
sources,
params,
v7_theta={'A': v7_cos_A, 'θ1': v7_cos_θ1},
a5_theta={'y_pts': a5_lines_y_pts, 't_pts': a5_lines_t_pts},
w1_theta={'A': w1_normal_A, 'sigma': w1_normal_sigma},
)
sol = solve_ivp(rhs, t_span, x0, t_eval=t_eval)
t = sol.t
q3 = sol.y[0] ## x_0
p6 = sol.y[1] ## x_1
v7_vals = f_cos(t_eval, A=v7_cos_A, θ1=v7_cos_θ1)
a5_vals = f_lines(t_eval, y_pts=a5_lines_y_pts, t_pts=a5_lines_t_pts)
w1_vals = f_normal(t_eval, A=w1_normal_A, sigma=w1_normal_sigma)
fig, ax = plot_states_and_sources(
t_span=t_span,
t_states=t,
state_series=[q3, p6],
state_labels=[
r'$q_{3}(t)$', r'$p_{6}(t)$'
],
state_styles=[
styles['xˣ[0]'], styles['xˣ[1]']
],
t_sources=t_eval,
source_series=[v7_vals, a5_vals, w1_vals],
source_labels = [r'$v_7(t)$', r'$a_5(t)$', r'$w_1(t)$'],
source_styles=[styles['vˣ[0]'], styles['a'], styles['w']],
)
sources = {
'v7': f_cos,
'a5': f_pulses,
'w1': f_normal,
}
v7_cos_A = 1.
v7_cos_θ1 = 3.
a5_pulses_amps = [7., -1., -3, 10]
a5_pulses_t_ons = [0., 1, 2, 5]
a5_pulses_t_offs = [1, 2, 5, 1e10]
w1_normal_A = 1
w1_normal_sigma = .1
rhs = make_rhs(
sources,
params,
v7_theta={'A': v7_cos_A, 'θ1': v7_cos_θ1},
a5_theta={'amps': a5_pulses_amps, 't_ons': a5_pulses_t_ons, 't_offs': a5_pulses_t_offs},
w1_theta={'A': w1_normal_A, 'sigma': w1_normal_sigma},
)
sol = solve_ivp(rhs, t_span, x0, t_eval=t_eval)
t = sol.t
q3 = sol.y[0] ## x_0
p6 = sol.y[1] ## x_1
v7_vals = f_cos(t_eval, A=v7_cos_A, θ1=v7_cos_θ1)
a5_vals = f_pulses(
t_eval,
amps=a5_pulses_amps,
t_ons=a5_pulses_t_ons,
t_offs=a5_pulses_t_offs,
)
w1_vals = f_normal(t_eval, A=w1_normal_A, sigma=w1_normal_sigma)
fig, ax = plot_states_and_sources(
t_span=t_span,
t_states=t,
state_series=[q3, p6],
state_labels=[
r'$q_{3}(t)$', r'$p_{6}(t)$'
],
state_styles=[
styles['xˣ[0]'], styles['xˣ[1]']
],
t_sources=t_eval,
source_series=[v7_vals, a5_vals, w1_vals],
source_labels = [r'$v_7(t)$', r'$a_5(t)$', r'$w_1(t)$'],
source_styles=[styles['vˣ[0]'], styles['a'], styles['w']],
)