how to plot control input from ode45

C
C
2011 9 5
3 个回答
更新时间：2018 12 30
4 次查看（30 天）

0 个投票

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.

请先登录，再回答此问题。

回答（3 个）

Walter Roberson
Walter Roberson 2011-9-5

0 个投票

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()
C
C 2011-9-6

0 个投票

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;

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Walter Roberson
Walter Roberson 2011-9-6
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"

请先登录，再进行评论。

C
C 2011-9-6

0 个投票

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?

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Walter Roberson
Walter Roberson 2011-9-6
Sorry, I have never tried that kind of plotting.
Bilal sadiq
Bilal sadiq 2018-7-17
Is any one there who can solve this problem for me ?
Walter Roberson
Walter Roberson 2018-7-17
Bilal sadiq, please clarify what problem it is that you are trying to solve?
Sunil Ojwani
Sunil Ojwani 2018-12-30
bilal sadiq , have you got your answer ??? i have same problem regarding control input plot from ode45

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

提问:

C
C
2011-9-5

评论：

Sunil Ojwani
Sunil Ojwani
2018-12-30

Community Treasure Hunt

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

Start Hunting!

Translated by