Как я могу оценить кинетический параметр стационарного реактора - PullRequest
Новая
       9

Как я могу оценить кинетический параметр стационарного реактора

0 голосов
/ 30 октября 2018

Я хочу решить систему ODE с некоторыми неизвестными параметрами k1,k2,k3: 

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

И у меня есть набор экспериментальных результатов со значениями C1,C2,C3 при x=0 (точка входа) и x=1 (конечная точка).

У меня нет данных между x=0 и x=1, чтобы решить их с помощью функции ODE, такой как ode45() или ode23(), а затем использовать функцию оптимизации. Как я могу решить эту проблему в MATLAB?

Ответы [ 2 ]

0 голосов
/ 31 октября 2018

Вы можете попробовать использовать решатель граничных значений 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 ;

enter image description here

с использованием аналогичного решателя 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()

enter image description here

, кажется, показывает правильное поведение. Прямая интеграция с использованием этих параметров с другим методом интеграции подтверждает правильность решения 

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()
0 голосов
/ 31 октября 2018

Это система из 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)))
...