Non-Negative Solution with Integrator in Simulink

4 visualizaciones (últimos 30 días)
Rajmohan el 7 de Mayo de 2024
Comentada: Sam Chak el 9 de Mayo de 2024
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?
0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

Paul el 7 de Mayo de 2024
Hi Rajmohan,
Do the parameters Limit Output along with Upper and Lower Saturation Limit do what's needed? Integrator
3 comentariosMostrar 1 comentario más antiguoOcultar 1 comentario más antiguo
Paul el 8 de Mayo de 2024
"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 el 9 de Mayo de 2024
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')

Iniciar sesión para comentar.

Más respuestas (1)

Sam Chak el 8 de Mayo de 2024
Editada: Sam Chak el 8 de Mayo de 2024
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')
1 comentarioMostrar -1 comentarios más antiguosOcultar -1 comentarios más antiguos
Sam Chak el 8 de Mayo de 2024
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.

Iniciar sesión para comentar.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by