No explanation of how to install keras

Odeint is no longer the recommended way to integrate ivps. We should rather use solve_ivp.

This template is an excellent version of what is possible.

The idea of finding parameters by matching coefficients is used often in textbooks, for instance to find the matching PID parameters for a DS or IMC design.

Hi,

I'm trying to use tbcontrol to study the dynamical response of a fuel cell. The model I'm trying to implement can be found in Fig. 5 of this paper. Fig. 5 is:

The parts I haven't figured out how to implement are the two yellow boxes in the image above, that is, Vcell and \delta.

I think Vcell is calculated using the simulation results for P_H2, P_H2O and P_O2. The results of Vcell then are used to update m through V_ac and \delta but \delta is updated based on the value for m.

The following code is my attempt. I'm using the algebraicequation E_fun but it is not working.
How can I pass the 3 pressures as inputs to the AlgebraicEquation E_fun?
I hope this makes sense and you can shed light on how to add the calculation of Vcell into the block diagram.
Thanks for the help.

import sympy as sym
from sympy.abc import s,t,x,y,z
from sympy.integrals import inverse_laplace_transform

import tbcontrol
from tbcontrol import blocksim

import numpy as np
import matplotlib.pyplot as plt

#Parameters
T = 343 #K
F = 96484600 # C/kmol
R = 8314.47 #J/(kmol K)
Eo = 0.6 #V
No = 88
Kr = 0.996e-6 #kmol/(s A)
U = 0.8
KH2 = 4.22e-5 #kmol/(s atm)
KH2O = 7.716e-6 #kmol/(s atm)
KO2 = 2.11e-5 #kmol/(s atm)
tH2 = 3.37 #s
tH2O = 18.418 #s
tO2 = 6.74 #s
t1 = 2 #s
t2 = 2 #s
t3 = 2 #s
CV = 2
B = 0.04777 #1/A
C = 0.0136 #V
Rint = 0.00303 #ohm
X = 0.05 #ohm
k3 = 1/(2*CV)
k5 = k6 = 10
Vr = 1 #pu
Q_meth_ref = 0.000015 #kmol/s
r_H_O = 1.168
Vs = 1 #pu
Km = 100 #pu
tm = 10 #pu
def Ip1_ae(t,I):
    return I*No/(2*F*U)
Gc1 = blocksim.AlgebraicEquation('Gc1', 'I', 'Ip1', Ip1_ae)

def Ip2_ae(t,I):
    return I*2*Kr
Gc2 = blocksim.AlgebraicEquation('Gc2', 'I', 'Ip2', Ip2_ae)

def Ip3_ae(t,I):
    return I*Kr
Gc3 = blocksim.AlgebraicEquation('Gc3', 'I', 'Ip3', Ip3_ae)

def qO2_ae(t,q_h2):
    return q_h2/r_H_O
Gc4 = blocksim.AlgebraicEquation('Gc4', 'q_h2', 'qO2', qO2_ae)

def Ip4_ae(t,I):
    return I*Rint
Gc5 = blocksim.AlgebraicEquation('Gc5', 'I', 'Ip4', Ip4_ae)

def i_ae(t,I):
    return B*sym.ln(C*I)
Gc6 = blocksim.AlgebraicEquation('Gc6', 'I', 'Ip5', i_ae)

#################### TRANSFERS ####################
G1 = blocksim.LTI('G1', 'i1', 'q_meth', [k3*t3, t3], [t3,0], 0)
G2 = blocksim.LTI('G2', 'q_meth', 'q_h2', CV, [t1*t2,(t1+t2),1], 0)
G3 = blocksim.LTI('G3', 'i3', 'p_h2', 1/KH2, [tH2,1], 0)
G4 = blocksim.LTI('G4', 'Ip2', 'p_h2o', 1/KH2O, [tH2O,1], 0)
G5 = blocksim.LTI('G5', 'i5', 'p_o2', 1/KO2, [tO2,1], 0)
G7 = blocksim.LTI('G7', 'i7','Vcell',1,1,0)

def E_fun(t,u):
    return No*( Eo+(R*T/2*F)*sym.log(p_h2*sym.sqrt(p_o2)/p_h2o) )
Gc7 = blocksim.AlgebraicEquation('Gc7','u','E',E_fun)   

diagram = blocksim.Diagram([G1, G2, G3, G4, G5, G7, Gc1, Gc2, Gc3, Gc4, Gc5, Gc6,Gc7],
        sums={  'i1': ('+Ip1', '-q_h2'),
                'i3': ('+q_h2','-Ip2'),
                'i5': ('+qO2','-Ip3'),
                'i7': ('+E','-Ip4','-Ip5')
             },
        inputs={'I': lambda t: 0.5*np.sin(0.005*t)})

ts = np.linspace(0,10,100)
sim_res = diagram.simulate(ts)
plt.plot(sim_res['p_h2'],'.')
plt.plot(sim_res['p_o2'],'.')
plt.plot(sim_res['p_h2o'],'.')
plt.show()

The input I is also a function of m and delta but for the time being I'm considering a sine function

In the current rendering on readthedocs [1], Slycot wasn't installed when the example was run.

Slycot 0.4.0 is available from conda-forge for Linux, Windows, and macOS (all 64-bit only).

[1] https://dynamics-and-control.readthedocs.io/en/latest/1_Dynamics/6_Multivariable_system_representation/Transfer%20function%20matrices.html#Conversion-to-state-space

Just wanted to let you know: the Laplace transform got major updates in SymPy 1.11 and 1.12 (and is still being updated for 1.13).

In Laplace Transforms.ipynb, you have an example More complicated inverses which shows hiow sometimes .apart() might help, but that example does not work anymore since laplace_transform will use .apart() by itself if appropriate.

Hi!

I just came across your great collection of notebooks (almost a textbook in fact). While reading, I noticed a typo in the last link of

https://dynamics-and-control.readthedocs.io/en/latest/2_Control/2_Laplace_domain_analysis_of_control_systems/Stability%20analysis.html

How can we determine stability without calculating the roots? See the next notebook.

This link is broken because it should be point to SymPy Routh Array.html (capital case in Sympy + extension). I guess that for the extension there is a tradeoff between having it work in the html export versus in the actual collection of notebook files.

Hi,

The result of simulate is in the time domain already or still in the Laplace (frequency) domain?

Thanks

Hi all,

First of all when I try to download the whole book I can't get the pdf. I do not want the .ipynb separately nor the html,

I prefer one pdf that contains all chapters in DYnamics and Control here:
https://dynamics-and-control.readthedocs.io/en/latest/index.html

I want to print it into a book so I can read it without internet. Thanks before. It is a very helpful to all of us want to learn.

Travis isn't free anymore, we should use GitHub actions

In Cheatsheet I get the following error on Binder:

On CoCalc I get a similar error after importing the complete repo.

Is there a way to import tbcontrol directly from a location on Github?

This appears to need slycot these days, perhaps add instructions on the correct install procedure or get rid of the dependency