how to plot control input from ode45

I'm trying to solve a nonlinear differential equation using ode45 I have a nonlinear system in the following form x = A*x + B*u and u is a function of the states x. Here is my function file,

function xprime = myode(t, x)

global g hs m S Cd Cl r0 rho0 h u

L = 0.5*rho0*exp(-(x(1)-r0)/hs)*x(3)^2*S*Cl/m;

f1 = x(3)*sin(x(6));
f2 = x(5);
f3 = -x(2) - g*sin(x(6));
f4 = x(3)*cos(x(6));
f5 = -x(5)*x(3)*sin(x(6))/hs + x(2)*(x(2) + g*sin(x(6)))*sin(x(6))/hs +...
    x(2)*cos(x(6))^2*(g-x(3)^2/x(1))/hs -...
    2*x(5)*(x(2) + g*sin(x(6)))/x(3) -2*x(2)*(x(2) +...
    g*sin(x(6)))^2/x(3)^2 - 2*x(2)*x(5)/x(3) +...
    2*x(2)*g*cos(x(6))^2*(g-x(3)^2/x(1))/x(3)^2;
f6 = -(g-x(3)^2/x(1))*cos(x(6))/x(3);

x2r = -.1221*t.^2 + 5.9859*t - 5.7866;
x5r =  -.2442*t + 5.9859;

x2rdot = x5r;
x5rdot =  -.2442;

x2rddot = x5rdot;

w1 =  -x(2)*cos(x(6))*(1/hs + 2*g/x(3)^2)*L;
w2 =  -x(2)*cos(x(6))*(1/hs + 2*g/x(3)^2)*L;

W = [w1;w2];

Delta = [h^2/2 0;0 h];

e1 = x(2) - x2r;
e2 = x(5) - x5r;

e = [e1;e2];

z1 = h*x(5) + (h^2/2)*f5;
z2 = h*f5;

z = [z1;z2];

d1 = h*x2rdot + h^2/2*x2rddot;
d2 = h*x5rdot;

d = [d1;d2];

u = -((Delta*W)'*Delta*W)^(-1)*(Delta*W)'*(e + z - d);

if( abs(u) > 1 )         %limits on control input 
        u = sign(u); 
end

xprime = zeros(6,1);
xprime(1) = f1;
xprime(2) = f2;
xprime(3) = f3;
xprime(4) = f4;
xprime(5) = f5 -x(2)*cos(x(6))*(1/hs + 2*g/x(3)^2)*L*u;
xprime(6) = f6 + L*u/x(3);

When I go to my main function file I can easily plot my states by using the command

plot(t,x(:,1)
plot(t,x(:,2)
.
.
.
plot(t,x(:,n))

and so on for n states. Since u is a function of these states I could plot u, i.e, u = blahblah(x(:,1),x(:,2),...,x(:,3)) and use the command plot(t,u)

However, my control input in this case is really long and complex and rewriting my control input in terms of my numerical answers x(:,n) seems time-consuming. Is there an easier way to plot my control input without having to rewrite the whole thing in my main function file? Any help is appreciated. Thanks.

Answer 1

[Walter Roberson]

I am not clear as to what your control input is for this case ?

The built in ways for easier plotting are through the options structure; see in particular http://www.mathworks.com/help/techdoc/ref/odeset.html#f92-1016858 and read down at least as far as the mention of odeplot()

Answer 2

[C]

Thanks for the reply. Here is my code using the options structure to plot my control input. I use options = odeset('OutputFcn',@getplot). And it seems to be plotting the input but only for a short time interval. I'm integrating my nonlinear system from t= 0 to 50 seconds but my control input only plots from t= 46.5 to 50 seconds. Is there anyway I can get around this? I have also marked where my control input is in case you're wondering. Thanks

function status = getplot(t,x,task)

global g hs m S Cd Cl r0 rho0 h u

switch task
    case ''

L = 0.5*rho0*exp(-(x(1)-r0)/hs)*x(3)^2*S*Cl/m;

f1 = x(3)*sin(x(6));
f2 = x(5);
f3 = -x(2) - g*sin(x(6));
f4 = x(3)*cos(x(6));
f5 = -x(5)*x(3)*sin(x(6))/hs + x(2)*(x(2) + g*sin(x(6)))*sin(x(6))/hs +...
    x(2)*cos(x(6))^2*(g-x(3)^2/x(1))/hs -...
    2*x(5)*(x(2) + g*sin(x(6)))/x(3) -2*x(2)*(x(2) +...
    g*sin(x(6)))^2/x(3)^2 - 2*x(2)*x(5)/x(3) +...
    2*x(2)*g*cos(x(6))^2*(g-x(3)^2/x(1))/x(3)^2;
f6 = -(g-x(3)^2/x(1))*cos(x(6))/x(3);

x2r = -.1221*t.^2 + 5.9859*t - 5.7866;
x5r =  -.2442*t + 5.9859;

x2rdot = x5r;
x5rdot =  -.2442;

x2rddot = x5rdot;

w1 =  -x(2)*cos(x(6))*(1/hs + 2*g/x(3)^2)*L;
w2 =  -x(2)*cos(x(6))*(1/hs + 2*g/x(3)^2)*L;

e1 = x(2) - x2r;
e2 = x(5) - x5r;

z1 = h*x(5) + (h^2/2)*f5;
z2 = h*f5;

d1 = h*x2rdot + h^2/2*x2rddot;
d2 = h*x5rdot;

% ******* control input u with w1,w2,e1,e2, etc. defined above ****
u = (h^4/4*w1^2+h^2*w2^2)^(-1)*(h^2/2*w1*(e1+z1-d1) + h*w2*(e2-z2-d2));  

% if( abs(u) > 1 )         %implementing limits on control input 
%         u = sign(u); 
% end

figure(9)
 plot(t,u)
end

status = 0;

[Walter Roberson]

Each plot() call is drawing a new plot overwriting the previous one. You only see the final one because you never give permission to MATLAB to draw the others. Add a drawnow() call after the plot() call to be able to see the intermediate ones, and have a look at the documentation for "hold on"

Answer 3

[C]

Thanks again for the reply. I have one quick question, however. I've made those minor adjustments you suggested but the plot is not looking smooth at all. The plot looks ragged like a piecewise function, is this supposed to happen?

[Walter Roberson]

Sorry, I have never tried that kind of plotting.

[Bilal sadiq]

Is any one there who can solve this problem for me ?

[Walter Roberson]

Bilal sadiq, please clarify what problem it is that you are trying to solve?

[Sunil Ojwani]

bilal sadiq , have you got your answer ??? i have same problem regarding control input plot from ode45