ODE suite for solving switched systems

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.

Answer 1

[Walter Roberson]

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, ....);

[Hadhiq Khan]

Let me elaborate a bit. I haven't used the ode function yet. The above code is a part of original code. I start with the initial conditions. Since I know the step size (h), & I know timespan, I know the value of n. So, what I am doing is that I begin solving the first differential equation. The condition on sw and X is to be checked after each iteration. Initially, the 'if' statement is true. Gradually, as the solution progresses, the 'else' statement become true. I wanted to write this as an ode function with the conditions as stated. In essence, depending on the 'sw', either of the two ode need to be solved. Hope you get my point.

[Walter Roberson]

None of the ode solvers use a fixed timestep. If you have a fixed timestep you should probably be using some kind of numeric integration such as trapz()

It is not clear which input your sw is in the ode function. Is it the time? Is it one of the inputs typically designated as "x", the second parameter to the objective function ? Is it a third input that is being calculated outside the ode function and is being passed in through the function handle? http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Your ode function should never have an "if" 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. If you need to compute differently based upon the values of the first two inputs, then you should instead use an event function to mark the crossing as terminal to exit the ode*() call, and the call ode*() again with a different ode function.

[Hadhiq Khan]

That is my point, Walter. I have used a fixed step here. I want to used the ode23tb (since the system is stiff). I am solving this differential equation:

Xdot = A*X + B*u.

My input is u. A and B are the system matrices. The system is such that it can either be 'on' or 'off'. This 'on' or 'off' is taken care of by a variable, say 'sw'. Consider the situation at n+1: If the value of X(2,n) and X(1,n) are above a certain threshold and the sw(n) is 1, at (n+1), the system will be 'on' and the following differential equation is to be solved:

Xdot = Aon*X + Bon*u

Similarly, if sw at (n) is 0 X(2,n) and X(1,n) are below that threshold, the system will be 'off' and X(n+1) will be solved based on:

 Xdot = Aoff*X + Boff*u

There has to be continuity in the variable 'X' from the 'on' to 'off' state and vice-versa. The terminal condition of the 'off' state will be used as initial condition of the 'on' state. The same will happen in the other case. The initial conditions on X and sw are: X=zeros(4,1). sw = 0. Using this, the system has to be solved for 1e-3 seconds. I wanted to solve this system using ode23tb. Since the ode will give me the solution over the entire time span, how do I approach this problem? How do I take care of the 'sw' in the event? Any suggestions?

Hope you get what I mean.

[Walter Roberson]

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.

[Hadhiq Khan]

Thanks for your answer. This is what I did:

[t_sol, x_sol] = ode15s('diode_ode',[t_start,t_end], x0)

and

function xdot = diode_ode(t,x);
global A1 B1 A0 B0 u1 u2 flag
if (flag == 0 && x(2,:) > 0.7) || ( flag == 1 && x(1,:) > 0)

    xdot = A1*x + B1*u1;
else

    xdot = A0*x + B0*u1;
end

Used the variable 'flag' for the condition. This did the job for me.

Answer 2

[Kwin]

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);

[Walter Roberson]

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]

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.