ODE suite for solving switched systems
11 views (last 30 days)
Show older comments
I have the following code.
The initial conditions for
X = [0;0;0;0]
sw = 0
n =100
h = 1e-9
for k = 1 : n
if (sw(k) == 0 && X(2,k)> 0.7 ) || (sw(k) == 1 && X(1,k) >0)
sw(k+1) = 1;
X(:,k+1) = X(:,k) +h*(A1*X(:,k) + B1*u);
else
sw(k+1) =0;
X(:,k+1) = X(:,k) +h*(A0*X(:,k) + B0*u);
end
end
I need to execute this code in the ODE suite. How do I take check sw and X after each iteration to see which of the two ODEs is to be solved.
0 Comments
Answers (2)
Walter Roberson
on 28 Aug 2016
You do not show two ODEs and you do not show any looping calling an ODE solver, so we do not know what you are trying to do.
Do not use "if" within a single ode objective function, not unless it is an "if" to trap errors or to determine whether you need to calculate the jacobian. If you have a "if" within an ode objective function then chances are that you have a discontinuity in the objective function and will not be able to calculate the integrals. If you have a situation where the calculation depends upon the range of the input t or x values, then you should probably be using "event functions" to signal the termination of the ode.
If you have multiple ode functions for different ranges then you can do the test outside the objective function. For example,
if t > 3 & t < 5
fun_handle = @first_ode;
else
fun_handle = @second_ode;
end
[t, x] = ode15s(fun_handle, [t, t+5], x0, ....);
5 Comments
Walter Roberson
on 29 Aug 2016
h = 1e-9;
tinit = 0;
tfinal = 1e-3;
timestep_list = tinit : h : tfinal;
n = length(timestep_list);
X_step = zeros(4,1);
sw = zeros(1, n);
X = zeros(4, n);
first_fun = @first_ode_fun;
second_fun = @second_ode_fun;
X(:, 1) = X_step;
for k = 1 : n - 1;
if (sw(k) == 0 && X(2,k) > 0.7 ) || (sw(k) == 1 && X(1,k) > 0)
sw(k+1) = 1;
[t, X_step] = ode23tb(first_fun, timestep_list(k:k+1), X_step);
else
sw(k+1) = 0;
[t, X_step] = ode23tb(second_fun, timestep_list(k:k+1), X_step);
end
X_step = X_step(:,end); %see note
X(:, k) = X_step;
end
Note: when you pass in exactly two times for the tspan then you get one output for each internal time that the ode happens to adaptively evaluate at. You only care about the final result, so that one has to be extracted.
Notice that the timestep list as a whole is not being passed to ode23tb to evaluate at. Although that would give you outputs only at those times (provided there were at least 3 of them), the ode*() routines assume continuity exists over the entire time range, and that is something that your system does not have.
This setup only checks the condition at the end of each of your fixed timesteps distance h apart, even though ode23tb() is going to be evaluating at a number of adaptively-determined times for each of your imposed fixed times. This setup will not change state as soon as the internal state of the ode might be able to detect the conditions are met: this setup only potentially changes state at the beginning of each of your imposed fixed steps. (You could rewrite it without difficulty to change state after each of the fixed steps.)
It would be possible for ode23tb() to determine pretty closely where the internal state would correspond to the switch conditions, at some intermediate time instead of at the fixed times: the mechanism for doing that would be to use an event function (and the loop would be written a bit differently.)
Warning: the h and final time that you indicated, h = 1e-9 and run to 1e-3, imply 10^6 imposed timesteps each with a call to ode23tb. The loop might take some time to finish. You would probably want to emit some output or update a graph along the way.
Kwin
on 11 Oct 2016
One way of doing this is by adding the two solution, but multiply the one that is inactive by zero,
dX = ((sw(k) == 0 && X(2,k)> 0.7 ) || (sw(k) == 1 && X(1,k) >0)) * (A1*X(:,k) + B1*u) + ...
~((sw(k) == 0 && X(2,k)> 0.7 ) || (sw(k) == 1 && X(1,k) >0)) * (A0*X(:,k) + B0*u);
2 Comments
Walter Roberson
on 12 Oct 2016
Your ode function should never have an "if" or conditional computation that depends upon the values of either of the first two inputs, the "t" or the "x": if you have such a thing then you almost certainly have a discontinuity in the function. All of the ode solvers assume continuity not just of values but of derivatives, and will break if you have a change of formula part way through the range.
Walter Roberson
on 12 Oct 2016
Edited: Walter Roberson
on 12 Oct 2016
In general there is a hidden problem with using
result = condition1(x) .* expression1(x) + condition2(x) .* expression2(x)
to try to implement the vectorized version of
if condition1(x(K))
result(K) = expression1(x(K));
elseif condition2(x(K))
result(K) = expression2(x(K));
end
The problem is that with the single-statement form, expression1(x(K)) and expression2(x(K)) are always executed even for the positions where the mask is false. If expression1(x(K)) is NaN where condition1(x(K)) is false, then the NaN times the 0 of the "false" will be NaN, and then it does not matter what the other condition produces, the result of the additional will be NaN. Likewise, if expression1(x(K)) is infinite where condition1(x(K)) is false, then the infinity times the 0 of the "false" will be NaN, and then it does not matter what the other condition produces, the result of the additional will be NaN. Similar constraints hold for condition2 and expression2. Therefore you can only use the fully vectorized expression if the result at the locations that you do not want selected are well-defined and not NaN or infinity.
For example you cannot use
(x ~= 0) .* (1./x) + (x -- 0) .* 1
because at 0, the 1./x will still be calculated, giving an infinity, that would then be multiplied by the false of x ~= 0, but 0 * infinity = NaN, and NaN + 1 is NaN
The logical indexing version does not suffer from that problem:
result = zeros(size(x));
mask = condition1(x);
result(mask) = expression1(x(mask));
mask = condition2(x);
result(mask) = expression2(x(mask));
However, in the context of an ODE, you are still at big risk of the ode* routines interpreting the change in slope as singularity.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!