I want to simulate the Point Reactor Kinetic Equations in Simulink. Solving the equations using ODE45 in MATLAB works fine but the Simulink model never achieves steady state. Can someone advise?

The equations as used in a MATLAB function are:

    function ndot=neudens(t,y,bet,B,lam,L)
        ndot(1) = ((react-B)/L)*y(1) + (lam(1)*y(2)) + (lam(2)*y(3)) + (lam(3)*y(4)) + (lam(4)*y(5)) + (lam(5)*y(6)) + (lam(6)*y(7));
        ndot(2) = ((bet(1)/L)*y(1)) - (lam(1)*y(2));
        ndot(3) = ((bet(2)/L)*y(1)) - (lam(2)*y(3));
        ndot(4) = ((bet(3)/L)*y(1)) - (lam(3)*y(4));
        ndot(5) = ((bet(4)/L)*y(1)) - (lam(4)*y(5));
        ndot(6) = ((bet(5)/L)*y(1)) - (lam(5)*y(6));
        ndot(7) = ((bet(6)/L)*y(1)) - (lam(6)*y(7));
    endfunction

bet = [0.000229,0.000832,0.000710,0.000852,0.000171,0.000102]';
B   = sum(bet);
lam = [0.0126,0.0337,0.139,0.325,1.13,2.50]';
L = 0.0005;

nt_0 = 1.0; Ct_0 = (bet(i)/(L*lam(i)))*nt; for i=1..6

I use ODE45 to solve the above equation and for react=0 the system achieves steady state, i.e., no change. But when I try solving the same set of equations in Simulink (model attached), the system never reaches steady state. What am I overlooking or doing wrong for it to behave this way?

Answer 1

[Zack Peters]

Hi Vikram,

While the algorithms for the ODE45 are the same across Simulink and MATLAB, they are implemented in different parts of the source code that have 0 overlap and can at times results in some floating point differences. With that said, You're is showing more than a floating point difference and I was able to implement several changes to get the results the same.

1. Now, to begin with I gave the same options to both. Here I used RelTol, AbsTol, and MaxStep with the odeset() function. This helped me get a sense of where the system was at. I saw that Simulink had a significant less data points and was shifted to the left a bit:

2. Simulink handles errors on a state-by-state basis. The model that was sent had 7 integrators, meaning 7 different states. Simulink will measure each one and if any one of them is out of tolerance then that will force a smaller step. With MATLAB ODE45, it acts slightly differently in that it will consider all of the outputs, y, as one vector and then determine the error based on the norm of the vector. To fix this I made Simulink have 1 single state that contains a vector. I also put everything into an EML block so that its easier to model.

After doing this you will see that the outputs "align" and the left-shift is gone, but the MATLAB ODE45 still has significantly more data points/noise in its signal:

3. If you examine the MATLAB ODE45 output closely you'll notice that what it's really doing is dividing every major step by 4. Its especially obvious when you specify an initial step size, like 1 second, and the first 5 output times are: 0, 0.25, 0.50, 0.75, and 1.0. This tells me that it's refining the output. The default is a refine value of 4, meaning it's going to refine if nothing else is specified. So I added Refine value of '1' to the odeset() function and we can now see that the results align perfectly:

I hope that this helps explain some of the nuances of the ODE functions and comparison between MATLAB and Simulink versions.

~Zack