Помогите, пожалуйста, кто-нибудь из специалистов в области химической инженерии, владеющий MATLAB: проверьте, нет ли проблем с этой динамической моделью

Thread Content

function nihe_0316 % Оценка параметров дифференциального уравнения % clear all; clc; format long; global t0; % k10,E,Ka соответствуют k(1), k(2), k(3) k0 =; % Начальные значения параметров, важно! Результат можно неоднократно подставлять для итеративных вычислений! ! x0 = 0.63713; % Начальное значение; если его нет, берется значение из первой строки данных.
lb = ; % Нижняя граница параметров; значения корректируются в соответствии с результатами аппроксимации.
ub = ; % Верхняя граница параметров.
% Экспериментальные данные: expdata= ... ;
t0=expdata(:,1); yexp = expdata(:,2); % yexp: соответствует экспериментальным данным.

% Используется функция lsqnonlin() для оценки параметров:
k=lsqnonlin(@ObjFunc4LNL, k0, lb, ub, x0, yexp);
fprintf('\n\nЗначения параметров, полученные с помощью функции lsqnonlin():\n')
fprintf('\tk10 = %.8f\n', k(1))
fprintf('\tE = %.8f\n', k(2))
fprintf('\tKa = %.8f\n', k(3))

% tspan = ; = ode45(@KineticEqs, tspan, x0, k);
x=spline(t, xx, t0); % Интерполяция для совпадения количества теоретических и экспериментальных значений.
% Вычисление коэффициента корреляции:
A=corrcoef(yexp, x(:,1));
fprintf('\nКоэффициент корреляции R равен:\n')
fprintf('\tR = %.8f\n', A(1,2))

% Вычисление коэффициента детерминации:
r(1)=1-sum((yexp-x(:,1)).^2)./sum((yexp-mean(yexp)).^2);
% ******************************** Коэффициент детерминации
fprintf('\nКоэффициент детерминации R^2 равен:\n')
fprintf('\tR^2 = %.8f\n', r(1))

% Вычисление среднеквадратичной ошибки:
rmse(1)=sqrt(sum((yexp-x(:,1)).^2)/length(t0));
fprintf('\nСреднеквадратичная ошибка RMSE равна:\n')
fprintf('\tRMSE = %.8f\n', rmse(1))

% Вычисление суммы квадратов остатков:
ssr(1)=sum((yexp(:,1)-x(:,1)).^2);
fprintf('\nСумма квадратов остатков для Xb равна:\n')
fprintf('\tSSE = %.8f\n', ssr(1))

% Рисунок:
figure(1) plot(t, xx(:,1), ‘b-’, t0, yexp, ‘ro’), legend(‘Теоретические значения’, ‘Экспериментальные значения’, ‘Location’, ‘best’)
xlabel(‘Время’); ylabel(‘Экспериментальные значения N’);

% Функция ObjFunc4LNL(k, x0, yexp):
function f = ObjFunc4LNL(k, x0, yexp)
tspan = ; = ode45(@KineticEqs, tspan, x0, k);
x=spline(t, xx, t0);
f = x - yexp; % Разница между теоретическими и экспериментальными значениями, остатки.
end

% Функция KineticEqs(t, Xb, k):
function dXbdt = KineticEqs(t, Xb, k)
% Дифференциальное уравнение; k10, E, Ka соответствуют k(1), k(2), k(3).
T=358.15; m=3; N0=0.1467; Cb0=0.5564; Ca=Cb0*(3-Xb); Cb=Cb0*(1-Xb); Cc=Cb0*Xb; Keq=exp(4873/T-13.813); Kb=5.2112*k(3); Kc=2.2695*k(3);
dXbdt = (m*k(1)*exp(-k(2)/8.314*T)*(Ca*Cb-Cc/Keq))/(N0*(1+k(3)*Ca+Kb*Cb+Kc*Cc)^2);
end

Результаты вычислений показывают, что значение X не меняется со временем. Теоретически это уравнение вида dx/dt=f(x). Я читал, что для уравнения вида dc/dt=f(c) можно использовать метод RK для вычислений и аппроксимации. Не понимаю, почему в моем случае это не работает. Пожалуйста, помогите мне разобраться