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{gm}{M}& 0\\ 0& 0& 0& 1\\ 0& 0& \frac{g(m+M)}{LM}& 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}{LM}\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 efficient 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.41561t}(0.0238095e^{6.41561t}{t}^2\theta (t)+0.000115693e^{3.2078t}\theta (t)-0.000231385e^{6.41561t}\theta (t)+0.000115693e^{9.62341t}\theta (t))\\ ,e^{-6.41561t}(-0.00242954e^{3.2078t}\theta (t)+0.00485909e^{6.41561t}\theta (t)-0.00242954e^{9.62341t}\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", 
 FrameLabel->{"time (sec)","\[Theta] (radians)"}, 
 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]);