我如何估计稳态反应器的动力学参数

Question

我想求解一个参数未知的 ODE 系统k1,k2,k3:

dC1/dx=-k1*C1
dC2/dx=k1*C1-k2*C2
dC3/dx=k2*C2-k3*C3

并且我有一组实验结果，在x=0（入口点）和x=1（终点）处的值为C1,C2,C3。

我在 x=0 和 x=1 之间没有任何数据可以使用像 ode45() 或 ode23() 这样的 ODE 函数来解决它，然后使用优化功能。我如何在 MATLAB 中解决这个问题？

Answer 1

这是一个包含 3 个常微分方程组的问题。 C1、C2、C3在x=0和x=1处的值是初始条件。您可以使用 dsolve 解析求解。在我的示例中，我假设初始条件是 C1(0)=1; C2(0)=1; C3(1)=1;。您可以根据自己的数据进行修改。

syms C1(x) C2(x) C3(x) k1 k2 k3

% Define the equations using == and represent differentiation using the diff function.
ode1 = diff(C1) == -k1*C1;
ode2 = diff(C2) == k1*C1-k2*C2;
ode3 = diff(C3) == k2*C2-k3*C3;
odes = [ode1; ode2; ode3]

% Define the initial conditions using ==
cond1 = C1(0) == 1;
cond2 = C2(0) == 1;
cond3 = C3(1) == 1;
conds = [cond1; cond2; cond3];

% dsolve function finds values for the constants that satisfy these conditions.
[C1Sol(x) C2Sol(x) C3Sol(x)] = dsolve(odes, conds)

输出

odes(x) =

           D(C1)(x) == -k1*C1(x)
 D(C2)(x) == k1*C1(x) - k2*C2(x)
 D(C3)(x) == k2*C2(x) - k3*C3(x)


C1Sol(x) =

exp(-k1*x)


C2Sol(x) =

-(exp(-k1*x)*exp(-k2*x)*((k1^2*k2*exp(k2*x))/((k1 - k2)*(k1 - k3)) - (k2^2*exp(k1*x)*(2*k1 - k2))/((k1 - k2)*(k2 - k3)) - (k1*k2*k3*exp(k2*x))/((k1 - k2)*(k1 - k3)) + (k2*k3*exp(k1*x)*(2*k1 - k2))/((k1 - k2)*(k2 - k3))))/k2


C3Sol(x) =

-exp(-k1*x)*exp(-k2*x)*exp(-k3*x)*((k2*exp(k1*x)*exp(k3*x)*(2*k1 - k2))/((k1 - k2)*(k2 - k3)) - (k1*k2*exp(k2*x)*exp(k3*x))/((k1 - k2)*(k1 - k3)) + (exp(k1*x)*exp(k2*x)*exp(k3)*(k1*k2^2 - k1^2*k2 - k1*k3^2 + k1^2*k3 + k2*k3^2 - k2^2*k3 + k1*k2^2*exp(-k1) + k1*k2^2*exp(-k2) - 2*k1^2*k2*exp(-k2) - k2^2*k3*exp(-k2) - k1*k2*k3*exp(-k1) + 2*k1*k2*k3*exp(-k2)))/((k1 - k2)*(k1 - k3)*(k2 - k3)))

Answer 2

您可以尝试使用边界值求解器bvp4c（或bvp5c）。可以满足的边界条件的数量是 ODE 的维度和未知参数数量的总和，因此您希望参数固定为 6 个已知值的 6 个可能条件。也就是说，传递给求解器的函数是

function dCdx = odesys(x,C,k)
    dCdx = zeros_like(C)
    dCdx(1) = -k(1)*C(1)
    dCdx(2) =  k(1)*C(1)-k(2)*C(2)
    dCdx(3) =  k(2)*C(2)-k(3)*C(3)

function bc = boundary(C0, CT, k)
    bc = [ C0(1)-c01, C0(2)-c02, C0(3)-c03, CT(1)-cT1, CT(2)-cT2, CT(3)-cT3 ]

作为初始猜测，您可以使用已知值之间的一些线性插值。

关于具有提供的边界值的修改系统

dC1/dx=-k1*C1/(1+k1*C1) ; 
dC2/dx=k1*C1/(1+k1*C1)-k2*C2 ; 
dC3/dx=k2*C2-k3*C3 ;

在脚本中使用类似的 python 求解器（但将 x1 更改为 100）

c0A, c0B, c0C = 293.3 , 2.1414, 3.6884
c1A, c1B, c1C = 208.09, 33.823, 78.561 
x0, x1 = 0.0, 100.0

def odesys(x,C,k):
    dCdx = zeros_like(C)
    dCdx[0] = -k[0]*C[0]/(1+k[0]*C[0]) ;
    dCdx[1] =  k[0]*C[0]/(1+k[0]*C[0])-k[1]*C[1]
    dCdx[2] =  k[1]*C[1]-k[2]*C[2]
    return dCdx

def bc(C0, C1, k):
    return [ C0[0]-c0A, C0[1]-c0B, C0[2]-c0C, C1[0]-c1A, C1[1]-c1B, C1[2]-c1C ]

x_init = linspace(x0,x1,21)
s = (x_init-x0)/(x1-x0)
C_init = [ c0A+s*(c1A-c0A), c0B+s*(c1B-c0B), c0C+s*(c1C-c0C)]

k = [1.0e-2, 1.0, 1.0]

res = solve_bvp(odesys, bc, x_init, C_init, k)
print res.message, res.p

它将参数 [ k1, k2, k3] 的结果打印为

The algorithm converged to the desired accuracy. 
[ 0.02319266  0.02248122 -0.00678455]

第三个参数是非物理负值暗示模型或某些参数不正确。此外，绘制解决方案：

x_sol=np.linspace(x0,x1,301)    
C_sol=res.sol(x_sol)
for k in range(3): plt.subplot(1,3,k+1); plt.plot(x_sol, C_sol[k]); plt.grid()
plt.show()

似乎显示了正确的行为。使用这些参数与另一种积分方法的正向积分证实了解决方案的正确性

C_int = odeint( lambda C,t: odesys(t,C,res.p), res.sol(0), x_sol)
for k in range(3): plt.subplot(1,3,k+1); plt.plot(x_sol, C_int[:,k]); plt.grid()
plt.show()

我如何估计稳态反应器的动力学参数

How can i estimate kinitic parameter of a steady-state reactor

optimization

matlab

model

curve-fitting

ode