Non-Negative Solution with Integrator in Simulink

When solving ode in Matlab, the ode solver has an option that allows setting the solution non-negative property to be true. For instance -

x = odeset; 
x.NonNegative = ones(numStates,1);

Is there an equivalent setting in Simulink? Something that can be adjusted in either model parameters or the [1/s] block?

Answer 1

[Paul]

Hi Rajmohan,

Do the parameters Limit Output along with Upper and Lower Saturation Limit do what's needed? Integrator

[Rajmohan]

Unfortunately, no.

The code equivalent of that would be -

y = ode15s(@callFun(),...);

yhat = min(y,0);

The saturation modifies the output only, it does not limit the state values and as a result the error can keep propagating.

[Paul]

"The saturation modifies the output only, it does not limit the state value"

I'm quite sure that's incorrect. I ran the equivalent of this model in Simulink, with Limit Output checked and Lower Saturation Limit set to 0 in the Integrator parameters and got the same result.

[t,x]=ode45(@(t,x) -t.*sin(t),[0 30],1,odeset('MaxStep',0.01,'NonNegative',1));
figure
plot(t,x)

[Sam Chak]

Hi @Rajmohan

Now that I've had the chance to test it in Simulink, I can confirm that @Paul's suggestion of setting the Lower Saturation Limit to 0 does indeed enforce the nonnegativity constraint and yield the correct solution. This option is also directly available in the Integrator Limited block. Here's a demo using the liquid level drainage system for reference.

If we don't impose the nonnegativity constraint (either mathematically or technically), the Simulink solver will terminate at 2 seconds, and the ode45 solver will provide an incorrect solution.

ode     = @(t, h) - sqrt(h);
[t, h]  = ode45(ode, [0 10], 1);
plot(t, h), grid on, xlabel('t'), ylabel('h(t)')

Warning: Imaginary parts of complex X and/or Y arguments ignored.

title('Solution without nonnegativity constraint')

Answer 2

[Sam Chak]

Hi @Rajmohan

I'm not familiar with all types of settings, codes, and tools in MATLAB/Simulink. Therefore, I usually resolve numerical issues using my own mathematical tricks, whenever possible.

In this case, you can attempt to address the ODE by incorporating a signum function in the Simulink model. You may find the demonstration below helpful, which is obtained from the "Nonnegative ODE Solution" article. However, I have also devised my own solution without relying on the built-in 'NonNegative' option.

Update: I have added a proper non-negative solution using the popular ReLU function in Case 4.

ode      = @(t, y) - abs(y);

% Analytic solution
tspan    = [0 40];
t        = linspace(tspan(1), tspan(2), 1001);
y        = exp(-t);

%% Case 1: Standard solution with ode45
options1 = odeset('Refine',1);
[t1, y1] = ode45(ode, tspan, 1, options1);

%% Case 2: Solution with nonnegative constraint
options2 = odeset(options1, 'NonNegative', 1);
[t2, y2] = ode45(ode, tspan, 1, options2);

%% Case 3: Solution with signum function
odefix   = @(t, y) - abs(y)*sign(y);
[t3, y3] = ode45(odefix, tspan, 1, options1);

%% Case 4: Solution with ReLU function
odeReLU  = @(t, y) - abs((sqrt(y.^2) + y)/2);   % ReLU(x) = (√(x²) + x)/2
[t4, y4] = ode45(odeReLU, tspan, 1, options1);

%% Plot results
plot(t, y, '-', t1, y1, 'o', t2, y2, '*', t3, y3, 'x', t4, y4, 'p'), grid on, ylim([-5 2])
xlabel('t'), ylabel('y(t)')
legend('Exact solution', 'No constraints', 'Nonnegativity', 'Signum Fix', 'ReLU Fix', 'Location', 'best')

[Sam Chak]

Or try using this Detect Rise Nonnegative block

• The output is true (equal to 1) when the input signal is greater than or equal to zero, and its previous value was less than zero.
• The output is false (equal to 0) when the input signal is less than zero, or if the input signal is nonnegative, its previous value was also nonnegative.