I’m seeking help from someone skilled in MATLAB and in chemical engineering to check whether there are any issues with this kinetic model

Thread Content

function nihe_0316 % Parameter estimation for differential equations % clear all; clc; format long; global t0; % k10,E,Ka correspond to k(1), k(2), k(3) k0 =; % Initial values of parameters, important! The results can be repeatedly substituted for iterative calculation! ! x0 = 0.63713; % Initial value; if not available, use the data from the first row. lb = ; % Lower bound of parameters, to be adjusted appropriately based on the fitting results. ub = ; % Upper bound of parameters. % Experimental data: expdata= ... ; t0=expdata(:,1); yexp = expdata(:,2); % yexp: corresponds to experimental data. % ------------------------------------------------------------------ % Use the lsqnonlin() function to estimate the parameters. k = lsqnonlin(@ObjFunc4LNL, k0, lb, ub, x0, yexp); fprintf('\n\nThe parameter values estimated using lsqnonlin() are:\n') fprintf('\tk10 = %.8f\n', k(1)) fprintf('\tE = %.8f\n', k(2)) fprintf('\tKa = %.8f\n', k(3)) % ---------------------------------------------------- tspan = ; x = ode45(@KineticEqs, tspan, x0, k); x = spline(t, xx, t0); % Interpolation to make the number of theoretical and experimental values equal. % Calculate the correlation coefficient. A = corrcoef(yexp, x(:,1)); fprintf('\nThe correlation coefficient R for Xb is:\n') fprintf('\tR = %.8f\n', A(1,2)) % % Calculate the coefficient of determination. r(1) = 1 – sum((yexp – x(:,1)).^2) / sum((yexp – mean(yexp)).^2); % ******************************** Coefficient of determination. fprintf('\nThe coefficient of determination R^2 for Xb is:\n') fprintf('\tR^2 = %.8f\n', r(1)) % Calculate the root mean square error. rmse(1) = sqrt(sum((yexp – x(:,1)).^2) / length(t0)); fprintf('\nThe root mean square error RMSE for Xb is:\n') fprintf('\tRMSE = %.8f\n', rmse(1)) % % Calculate the sum of squared errors (SSE). ssr(1) = sum((yexp(:,1) – x(:,1)).^2); fprintf('\nThe sum of squared errors for Xb is:\n') fprintf('\tSSE = %.8f\n', ssr(1)) % ---------------------------------------------------- figure(1) plot(t, xx(:,1), ‘b-’, t0, yexp, ‘ro’), legend(‘Theoretical value’, ‘Experimental value’, ‘Location’, ‘best’) xlabel(‘Time’); ylabel(‘Experimental value N’); % ------------------------------------------------------------------ function f = ObjFunc4LNL(k, x0, yexp) tspan = ; x = ode45(@KineticEqs, tspan, x0, k); x = spline(t, xx, t0); f = x – yexp; % The difference between theoretical and experimental values, i.e., the residuals. end % ------------------------------------------------------------------ function dXbdt = KineticEqs(t, Xb, k) % Differential equation; k10, E, Ka correspond to 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 % ------------------------------------------------------------------ end The calculated value of X does not change over time. Theoretically, it should be a differential equation of the form dx/dt = f(x). I’ve seen that for equations of the form dc/dt = f(c), the RK method can be used for calculation and fitting. I’m not sure why it doesn’t work in my case. Please help me figure this out