solving state equation by ODE 45 by user defined time step and also checking response at each timestep

Hi I want to calculate response x (fixed time step=0.01 sec) by solving state equation using ODE45: xdot=A.x+B.u --> A,B,x0(initial value of x) & u are known to me, but I will decide (or modify) the input u at (i+1)-th time based on the response I am getting at i-th time step. If I use ODE45 directly then I will get the final responses at each time and I will not be able to modify u accordingly. May be some loop is required to do that but HOW!!! Please help me.Thanks in advance.

Answer 1

[Star Strider]

The MATLAB ODE solvers return the current value of your independent variable (here, time) and the current values of the integrated function. If you have a time-dependent input, you can use the time value to change the input. If you are making small adjustments to ‘u’, this should not be a problem.

Note that the ODE solvers don’t deal well with significant discontinuities, so consider stopping the integration and re-starting it with new values for ‘u’ (using the previous final values of ‘x’ as the new initial conditions) if the discontinuities are significant.

[Kamal Bera]

I think you did not get my point. My point are (1) I want to use ODE 45 with fixed time step, not the arbitrary one as taken by default. (2) During solution I want to track the response at each time step so that I can put some constraint on 'u' according to the response of the previous step. But if I use ode45 syntax directly this is not possible. How to do that.

Mr Strider I am not getting your reply. Can you please write the code for me. Thank you.

[Star Strider]

To use a fixed step, instead of using:

tspan = [t_start  t_end];

use:

tspan = linspace(t_start, t_end, 50);

to create a 50-element fixed-step time vector, for example.

[Kamal Bera]

Thank you for your quick response. Now my question no (1) is solved. But what about question (2). I think some 'for'loop is necessary e.g.

    for i=1:50  % for 50-element fixed-step time vector
        [t(i),y(i)]=ode45(@FN,tspan(i),y0);
        y0=y(i);
    end
Something like this. So that from y(i), I can decide my 'u' for next step.( *Above code is not correct*)
Please help me. I am not able to progress further. I will accept all your answer.

[Star Strider]

That will throw an error. Your tspan has to be a vector, so you have to define it using linspace at each step. A single value won’t work. You code would have to become something like:

tspan = linspace(0, t_end, 50);
y0 = . . .;
  for i=1:50  % for 50-element fixed-step time vector
      [t(:,i),y(:,:,i)]=ode45(@FN,tspan,y0);
      y0=y(end,:,i);
      tspan = linspace(t(end,i), t(end,i)+T, 50);
  end

The ‘T’ increment is a value you would have to define. Keep the ‘tspan’ vector to the same length (I chose 50 arbitrarily, but they all have to be the same lengths to save them in your array. You don’t have to save it in a separate array, since it is saved in your ‘t’ matrix.)

Note: This is UNTESTED CODE.

[Kamal Bera]

Hi Mr Strider, I have run your code for a small problem(ydot=A*y; y has two component y1 & y2). Finally I am getting the output y as an array of size 50x2x50; the required solution is y(:,:,1). But I am not able to understand what 'T' is doing, as my required solution not affected with the choice of 'T'. I have accepted your answer.

Now my actual problem is ydot=A*y+B*u; where u is time dependent. With reference to your answer dated 27 Jan 2016 at 14:47, you have written | If you have a time-dependent input, you can use the time value to change the input.|Can you clarify me (1) How to write the function FN which will represent the differential equation ydot=A*y+B*u. (2)How to use the time value to change the input as suggested by you.

It will be very helpful if you write the code for me. Again thank you very much for your timely response.

[Star Strider]

I appreciate your having Accepted my Answer. However I cannot write the exact code for you, in part because if you are doing control theory, you need to learn how to program control systems.

The ‘T’ value is a constant you define that defines the ‘tspan’ vector, and therefore the end of the next time increment. It ensures that all the time increments are the same length and duration. I should have included it in my code outline, so including it now:

tspan = linspace(0, t_end, 50);                 % Initial ‘tspan’ Vector
T = 50;                                         % Constant Chosen To Define The End Time Of Each ‘tspan’ Value
y0 = . . .;
for i=1:50  % for 50-element fixed-step time vector
    [t(:,i),y(:,:,i)]=ode45(@FN,tspan,y0);
    y0=y(end,:,i);
    tspan = linspace(t(end,i), t(end,i)+T, 50);
end

If you have a linear problem, you don’t have to use the ODE solvers. You can use the matrix exponential expm function, but you would have to use the solved version of your control equation. You control theory textbook should explain that in detail (in involves matrix Laplace transforms).

‘Now my actual problem is ydot=A*y+B*u; where u is time dependent.’ You solve this using the ‘y’ vector your ‘FN’ ODE function takes as input, doing the necessary A*y matrix multiplication (creating a column vector), then adding the column vector created by multiplying B*u. This creates the column vector ‘FN’ returns to ode45 (the usual solve for these problems).

For a time-dependent ‘u(t)’, you have to define it. I do not know what form it takes, so I can only say that if it is a continuous function, you can write it in your ‘FN’ function as an anonymous function. If it is a discontinuous step function, it is best to stop the integration, increment ‘u(t)’, and continue. To include an external input to ‘FN’, you would create it as for example:

FN = @(t,y,u) A*y + B*u;

Prototype code for your problem (that I tested to verify that it works):

A = [0.1 -0.2; 0.3 0.4];
B = [0.4; 0.5];
u = 3;
y0 = [0.1; 0.2];
FN = @(t,y,u) A*y + B*u;
tspan = linspace(0, 10, 20);
T = 10;
for i=1:10  % for 50-element fixed-step time vector
    u = i/100;
    [t(:,i),y(:,:,i)]=ode45(@(t,y) FN(t,y,u),tspan,y0);
    y0=y(end,:,i);
    tspan = linspace(t(end,i), t(end,i)+T, 20);
end
tv = reshape(t, [], 1);
yp = permute(y, [1 3 2]);
yv = reshape(yp, [], size(y,2));

figure(1)
plot(tv, yv)
grid

You will have to make the necessary changes in this code to adapt it to your control system.

[Kamal Bera]

Hi Mr Strider, My problem is also in the similar line except in my case I will assign ' u' as some constant multiplied by the response at each step. I have tried a lot but unfortunately I am not getting the controlled response. At i=1, u is assigned as zero, then for i=2, the code is not taking the calculated u value. As a result I am always getting the uncontrolled response. I think some minor adjustment is needed so that it can take new u value at each time step. I am attaching my matlab code. I am requesting you to kindly go through it and make the proper adjustment, so as to obtain the controlled response. I will be highly obliged if you kindly consider my request.

[Star Strider]

I have a reasonable background in modern control theory, but I have no idea what you mean by ‘controlled response’. That could be anything. You have to design your own control system using linear predictive control, pole placement, or however you want to implement it.

Linear control is in general not difficult to implement, but you need to take the time to understand the problem, understand the techniques, and then design your control system to control your plant to produce the result you want.

If you have the Control System Toolbox, this is considerably easier to do.

[Kamal Bera]

Actually I am getting the same response with ' u' and without * u*, where ' u=-k*y', ' k' is the output feedback gain already calculated. But it is not always ' u=-k*y',as can be seen from if else statement inside the for loop. That means there is some additional restriction on control input. In my code irrespective of the value of ' k', I am getting the same response, which is not possible. That means the effect of ' u' is not incorporated when solving equation: xdot=A*x+B*u. I am doing something wrong inside the _ for loop. Please help me.

• I have already solved the whole control problem without using ODE45 (I have used matrix exponential). But as per my supervisor requirement I have to solve it by using Runge-Kutta method which as per him is numerically more accurate and I have to follow his instruction. Otherwise he will not accept the result.

[Star Strider]

I did some experiments with my code (that I posted here earlier) using various stepwise assignments to ‘u’. The effect of ‘u’ is definitely apparent, but you may have to scale it up to see the effect. In my code, defining ‘u’ as:

u = (i^4*(-1)^i);

definitely demonstrates that the function uses ‘u’ appropriately in its calculations, as it should. I am not certain what the scaling factor should be (perhaps 1000?) but it seems you will have to experiment to see the effect.

As I mentioned, I usually use the Control System Toolbox for these sorts of problems, so it’s not obvious to me why ‘u’ is not having more of an effect than it is in this code. I am satisfied that my code is working correctly.

[Kamal Bera]

definitely your code is working. But you are defining 'u' as u = (i^4*(-1)^i); before calling the function FN. In my case from the first step (i=1), I am getting some response & that is used to get the 'u' which will be used in the next step(i=2). So inside the loop I have tell the code to take 'u' obtained from step 1 as 'u' of the step 2. See line no 43 in my code i.e. u=u(:,i). I think this is where I am making the mistake. Please see the code and help. I know it is easy for you.

[Star Strider]

When I searched through your code, I don’t see that you’ve defined ‘B’ anywhere, and certainly not prior to your ‘FN’ function definition. (I don’t know why that would not have thrown an error.) If ‘B’ is undefined (and MATLAB is case-sensitive, so it would have to be assigned as ‘B’ for it to work), that would explain your not seeing any effect of ‘u’, because it is not being transmitted to the rest of your ‘FN’ function.

Since it’s apparently not throwing an error, you need to find out what ‘B’ is assigned as. That could be the problem.

[Kamal Bera]

Sorry. I forget to write the 'B' matrix in the code that I send to you. Actually it was there in the code that I run and it does not give any error related to ' B'. I am again sending the corrected code. Please check it. Again extremely sorry fir the mistake.

[Star Strider]

I don’t see what the problem is. Your ‘u’ doesn’t change much but it does definitely change between iterations, and adding my modified plot code at the end:

tv = reshape(t, [], 1);
yp = permute(y, [1 3 2]);
yv = reshape(yp, [], size(y,2));

tv = [0:size(yv,2)-1]*0.001;
figure(1)
plot(tv,yv)
grid

shows the output to be a bounded damped sinusoid. I would save ‘u’ to a vector (perhaps ‘uv’), and then plot it against ‘tv’ as well to see how its changes correspond with the output.

[Kamal Bera]

No sir, there is something wrong. Just check with input 'u' as zero. There is no difference between the plots (with 'u' as written in the code vs 'u' as zero). I have checked. Please check it. There must be some difference (whatever small it may be) between two cases. I am attaching the code(' without using ODE45') which I have already run to solve this problem. Kindly see the response. Another doubt regarding your truncated plot through tv = [0:size(yv,2)-1]*0.001;when I am writing ' 1' instead of '0.001', the range of x axis it shows is from 0 to 1201, it should be from 0 to 120 sec as my final time is *120 sec* only. ' * 1201*' is number of time step. May be I have implemented your idea(original code) wrongly.||Please check.

[Star Strider]

That is not the correct implementation of expm in that instance.

Derivation:

x' = A*x + B*u                          % Original System Differential Equation
s*X = A*X + B*U                         % Matrix Laplace Transform
s*X - A*X = B*U                         % Collect Terms
(s*I - A)*X = B*U                       % I = Identity Matrix
X = (s*I - A)^-1 * B*U                  % ‘^-1’ => Inverse
x = expm(A*t)*B*u                       % Take Inverse Laplace Trasnform To Get The Result
                                Q.E.D.

Continuing on:

y = C*x + D*u
y = C*expm(A*t)*B*u + D*u
y = (C*expm(A*t)*B + D)*u

Re-write your equation using that derivation, and it should work out to be the same as that using ode45.

[Kamal Bera]

Please see the following corrected derivation

x' = A*x + B*u                 % Original System Differential Equation
s*X-x(0) = A*X + B*U           % Matrix Laplace Transform, x(0)=x at time t=0;
s*X - A*X = x(0) + B*U         % Collect Terms
(s*I - A)*X = x(0)+B*U         % I = Identity Matrix
X = (s*I - A)^-1 *x(0) +(s*I - A)^-1*B*U                  
Now if you take inverse laplace transform 
x=exp(A*t)*x(0)+integration of [exp(A*(t-tao)*B*u(tao)*dtao]
The second term can be approximated and for a time marching scheme the following expression can be written.
x(n*deltat)=exp(A*deltat)*x((n-1)*deltat)+B*u((n-1)*deltat)*deltat; n=1,2,3.....                    
I am writing the above expression from the book Modern Control Design with MATLAB and SIMULINK by Ashish Tewari. Please see section 4.4 in chapter 4 & equation no. (4.66). This book can be downloaded online.

[Star Strider]

I’ll look at the book tomorrow. It’s the end of my day here.

In control theory as I learned it, all initial conditions are assumed to be zero, by convention. My derivation is correct according to Chen’s text. It needs no correction.

The point however is that if you use my implementation of the solved differential equation, it should closely approximate the result produced by ode45. It may also with the version you quoted. I’ve not checked it. I will leave that to you.

[Kamal Bera]

Please give a last try with the code ContRespNew.m that I have already sent to you on 30 Jan 2016 at 18:56. Change ' k' (that will eventually change the ' u' also) whatever you want.You will get same response each time. I think this is not correct. If correct, can you explain me!!!

[Star Strider]

If that’s the file with the ode45 call, I already looked at it, ran it, and commented on it. It appears to be correct to me.

The other problem was your (incorrect) use of the expm function that I also commented on in some detail.

It is quite possible to quickly converge on a similar-appearing solution. The only way you can determine if it is working correctly is to closely examine what it is doing over time, preferably by plotting the output of the plant and the control inputs over time.

You would have to save the control signals, perhaps as ‘uv’, by inserting this statement at the end of your loop (after defining the new value of ‘u’ in each iteration):

uv(:,i) = u;

Then in your plot, plot ‘uv’ against ‘t’, as well as the signals you are already plotting. I am certain you will see them change appropriately.

Answer 2

[Nithin Chava]

?(?) = ?+6 ?(? 2+4?+3)

[Star Strider]

What?