### 1.34 Plot the response of the inverted pendulum problem using state space

Problem: Given the inverted pendulum shown below, use state space using one input (the force on the cart) and 2 outputs (the cart horizontal displacement, and the pendulum angle. Plot each output separately for the same input.

Mathematica

Remove["Global*"];
a = {{0, 1, 0, 0},
{0, 0, ((-m)*g)/M, 0},
{0, 0, 0, 1},
{0, 0, ((M + m)*g)/(M*L), 0}};



$\left ( {\begin {array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & -\frac {g m}{M} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & \frac {g (m+M)}{L M} & 0 \\ \end {array}} \right )$

b = {{0}, {1/M}, {0}, {-(1/(M*L))}};
MatrixForm[b]



$\left ( {\begin {array}{c} 0 \\ \frac {1}{M} \\ 0 \\ -\frac {1}{L M} \\ \end {array}} \right )$

c = {{1, 0, 0, 0}, {0, 0, 1, 0}};
MatrixForm[b]



$\left ( {\begin {array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \end {array}} \right )$

sys=StateSpaceModel[{a,b ,c}]



It is now possible to obtain the response of the system as analytical expression or an interpolatingFunction.

It is much more eļ¬cient to obtain the response as interpolatingFunction. This requires that the time domain be given.

Here is example of obtaining the analytical expression of the response

sys=sys/.{g->9.8,m->1,M->20,L->1};
y=Chop[OutputResponse[sys,UnitStep[t],t]]



$${\begin {array}{l} e^{-6.41561 t} (0.0238095 e^{6.41561 t} t^2 \theta (t)+0.000115693 e^{3.2078 t} \theta (t)-0.000231385 e^{6.41561 t} \theta (t)+0.000115693 e^{9.62341 t} \theta (t))\\,e^{-6.41561 t} (-0.00242954 e^{3.2078 t} \theta (t)+0.00485909 e^{6.41561 t} \theta (t)-0.00242954 e^{9.62341 t} \theta (t)) \end {array}}$$

Plot[y[[1]],{t,0,3},
PlotLabel->"y(t) response",
FrameLabel->{"time (sec)","y (meter)"},
Frame->True,
PlotRange->All,
ImageSize->300,
GridLines->Automatic,
GridLinesStyle->Dashed,
RotateLabel->False]



Plot[ y[[2]],{t,0,3},
PlotLabel->"\[Theta](t) response",
Frame->True,
PlotRange->All,
ImageSize->300,
GridLines->Automatic,
GridLinesStyle->Dashed,
RotateLabel->False]



Matlab

clear all; close all;

m=1; M=20; L=1; g=9.8;

A=[0 1      0           0;
0 0   -m*g/M         0;
0 0      0           1;
0 0  (M+m)*g/(M*L)   0
];
B=[0  1/M  0 -1/(M*L)]';
C=[1 0 0 0;
0 0 1 0];
D=[0];

sys=ss(A,B,C(1,:),0);
step(sys,3);
grid on;
set(gcf,'Position',[10,10,320,320]);
title('y(t) due to step input');



figure;
sys=ss(A,B,C(2,:),0);
step(sys,3);
title('\theta(t) due to step input');
grid on;
set(gcf,'Position',[10,10,320,320]);

`