使用 tbcontrol 和 python 将框图中的变量相乘

如何解决使用 tbcontrol 和 python 将框图中的变量相乘

我正在尝试使用 sympy 和 block diagram below 通过 python 实现 tbcontrol。我一直在计算 E。我可以得到 P_H2、P_H2O 和 P_{O2 使用 sim_res = diagram.simulate(ts) 运行模拟后，我无法计算 Vcell 或 m 或 DC/AC 逆变器块和反馈。}

我的问题是如何在 E 中包含 diagram 的计算？

我尝试使用函数 E_fun，但它不起作用。

感谢任何帮助。谢谢

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','Ip2',Ip2_ae)

def Ip3_ae(t,I):
    return I*Kr
Gc3 = blocksim.AlgebraicEquation('Gc3','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','Ip4',Ip4_ae)

#################### AE B0 ####################
def i_ae(t,I):
    return B*sym.ln(C*I)
Gc6 = blocksim.AlgebraicEquation('Gc6','Ip5',i_ae)
#################### TRANSFERS ####################
G1 = blocksim.LTI('G1','i1','q_meth',[k3*t3,t3],[t3,0],0)
G2 = blocksim.LTI('G2',CV,[t1*t2,(t1+t2),1],0)
G3 = blocksim.LTI('G3','i3','p_h2',1/KH2,[tH2,0)
G4 = blocksim.LTI('G4','p_h2o',1/KH2O,[tH2O,0)
G5 = blocksim.LTI('G5','i5','p_o2',1/KO2,[tO2,0)
G7 = blocksim.LTI('G7','i7','Vcell',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()