Matlab 仿真——单自由度倒立摆（1）系统建模

1. 受控对象与设计要求

该例的系统包含一个装有一个倒立摆的小车。我们的控制目标是通过给小车作用一个力，使顶部的倒立摆不落下来。下图标示出了各个变量的含义。

对于这个例子，我们有以下参数：
(M) 小车质量 0.5 kg
(m) 倒立摆质量 0.2 kg
(b) 小车的摩擦系数 0.1 N/m/sec
(l) 转动点到倒立摆质心的长度 0.3 m
(I) 倒立摆的转动惯量 0.006 kg.m^2
(F) 小车受到的力
(x) 小车的位置坐标
(theta) 倒立摆与垂线的角度

我们会用不同的方法设计该系统的控制器：PID，根轨迹，频域分析，状态空间。由于PID，根轨迹以及频域分析法最适合SISO系统的分析研究，因此用这三种方法的时候我们不考虑小车的位置，仅仅考虑如何控制好倒立摆的角度。我们会设计一款控制器，使得小车收到一个冲击（1Nsec）时仍能保持倒立摆垂直，设计要求是5s内稳定倒立摆，且倒立摆角度与稳定状态下的角度相比不超过0.05 弧度。倒立摆的初始位置垂直于地面。因此，我们可以总结出我们的设计需求：

1. θ的稳定时间 < 5s
2. |θ-θ0| < 0.05 radians

而对于状态空间设计方法，我们得以能够处理多输出系统。因此，在用状态空间设计法的时候我们尝试同时控制倒立摆的角度以及小车的位置。以下是我们的设计需求：

1. x 与 θ 的稳定时间 < 5s
2. x 的上升时间 < 0.5s
3. |θ-θ0| < 0.05 radians （也就是20°）
4. 对于x和θ来说，稳态误差 < 2%

2. 力分析与系统方程

下图是该系统的受力分析图：

水平方向对小车进行受力分析得到（1）式：
$$M\ddot{x}+b\dot{x}+N=F$$
你也可以在竖直方向对小车进行受力分析，但没什么卵用。接着在水平方向对倒立摆进行受力分析，得到（2）式 (这一步推导有点快，请参考这个文档)：
$$N=m\ddot{x}+ml\ddot{\theta }\cos \theta -ml{\dot{\theta }}^2\sin \theta$$
接着把（2）式代入到（1）式得到该系统的第一个动力学公式：
第一个动力学方程： $$(M+m)\ddot{x}+b\dot{x}+ml\ddot{\theta }\cos \theta -ml{\dot{\theta }}^2\sin \theta =F$$
为了得到第二个动力学公式，我们在垂直于倒立摆的方向上列出力平衡方程：
$$P\sin \theta +N\cos \theta -mg\sin \theta =ml\ddot{\theta }+m\ddot{x}\cos \theta$$
左边的三项不难理解，右边的第一项为引起转动加速度的力，第二项为引起平动加速度的力。为了消去P和N，我们再列出关于倒立摆重心的转动惯量方程：
$$-Pl\sin \theta -Nl\cos \theta =I\ddot{\theta }$$
稍微整合一下上述两个方程得到我们要用的第二个动力学方程：
第二个动力学方程： $$(I+m{l}^2)\ddot{\theta }+mgl\sin \theta =-ml\ddot{x}\cos \theta$$
上述两个动力学方程是非线性的，为了采用线性系统的方法去设计控制器，还需要把上述两个动力学方程线性化一下。因此我们假设倒立摆的平衡位置 $\theta$ = $\pi$，且倒立摆只在平衡位置下小范围内摆动。我们用$\varphi$ 代表倒立摆与平衡位置的角度偏差，那么倒立摆在任意时刻的位置 $\theta$ = $\pi$ + $\varphi$。有了上述线性化假设之后，我们可以得到：

1. $$\cos \theta =\cos (\pi +\varphi )\approx -1$$
2. $$\sin \theta =\sin (\pi +\varphi )\approx -\varphi$$
3. $${\dot{\theta }}^2={\dot{\varphi }}^2\approx 0$$

将上述结果带入到我们的动力学方程组中得到：
第一个动力学方程： $$(M+m)\ddot{x}+b\dot{x}-ml\ddot{\varphi }=u$$
第二个动力学方程： $$(I+m{l}^2)\ddot{\varphi }-mgl\varphi =ml\ddot{x}$$

2.1 转换方程

对动力学方程组分别做拉氏变换
$$(I+m{l}^2)\Phi (s)s^2-mgl\Phi (s)=mlX(s)s^2$$
$$(M+m)X(s)s^2+bX(s)s-ml\Phi (s)s^2=U(s)$$
要得到输入$U(s)$与输出 $\Phi (s)$ 的关系,我们需要消掉$X(s)$,于是经过一系列的化简代入我们最终得到：
$$P_{pend}(s)=\frac{\Phi (s)}{U(s)}=\frac{\frac{ml}{q}s}{s^3+\frac{b(I+m{l}^2)}{q}{s}^2-\frac{(M+m)mgl}{q}s-\frac{bmgl}{q}}[\frac{rad}{N}]$$
$$P_{cart}(s)=\frac{X(s)}{U(s)}=\frac{\frac{(I+m{l}^2)s^2-gml}{q}}{s^4+\frac{b(I+m{l}^2)}{q}{s}^3-\frac{(M+m)mgl}{q}{s}^2-\frac{bmgl}{q}s}[\frac{m}{N}]$$
其中：
$$q=[(M+m)(I+m{l}^2)-(ml{)}^2]$$

2.2 状态空间

有了线性化之后的动力学方程组，一顿化简得到：
$$\left[\begin{array}{c}\dot{x}\\ \ddot{x}\\ \dot{\varphi }\\ \ddot{\varphi }\end{array}\right]=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& \frac{-(I+m{l}^2)b}{I(M+m)+Mm{l}^2}& \frac{m^{2}g{l}^2}{I(M+m)+Mm{l}^2}& 0\\ 0& 0& 0& 1\\ 0& \frac{-mlb}{I(M+m)+Mm{l}^2}& \frac{mgl(M+m)}{I(M+m)+Mm{l}^2}& 0\end{array}\right]\left[\begin{array}{c}x\\ \dot{x}\\ \varphi \\ \dot{\varphi }\end{array}\right]+\left[\begin{array}{c}0\\ \frac{I+m{l}^2}{I(M+m)+Mm{l}^2}\\ 0\\ \frac{ml}{I(M+m)+Mm{l}^2}\end{array}\right]u$$
$$\mathbf{y}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}x\\ \dot{x}\\ \varphi \\ \dot{\varphi }\end{array}\right]+\left[\begin{array}{c}0\\ 0\end{array}\right]u$$
该系统的输出有两个，一是小车的位置$x$ ，二是倒立摆的角度$\theta$ ,因此$y$有两项

3. Matlab表达

现在我们就可以在Matlab里面表达出我们的系统了

3.1 转换方程

M = 0.5;
m = 0.2;
b = 0.1;
I = 0.006;
g = 9.8;
l = 0.3;
q = (M+m)*(I+m*l^2)-(m*l)^2;
s = tf('s');

P_cart = (((I+m*l^2)/q)*s^2 - (m*g*l/q))/(s^4 + (b*(I + m*l^2))*s^3/q - ((M + m)*m*g*l)*s^2/q - b*m*g*l*s/q);

P_pend = (m*l*s/q)/(s^3 + (b*(I + m*l^2))*s^2/q - ((M + m)*m*g*l)*s/q - b*m*g*l/q);

sys_tf = [P_cart ; P_pend];

inputs = {'u'};
outputs = {'x'; 'phi'};

set(sys_tf,'InputName',inputs)
set(sys_tf,'OutputName',outputs)

sys_tf

输出：

3.2 状态空间

M = .5;
m = 0.2;
b = 0.1;
I = 0.006;
g = 9.8;
l = 0.3;

p = I*(M+m)+M*m*l^2; %denominator for the A and B matrices

A = [0      1              0           0;
     0 -(I+m*l^2)*b/p  (m^2*g*l^2)/p   0;
     0      0              0           1;
     0 -(m*l*b)/p       m*g*l*(M+m)/p  0];
B = [     0;
     (I+m*l^2)/p;
          0;
        m*l/p];
C = [1 0 0 0;
     0 0 1 0];
D = [0;
     0];

states = {'x' 'x_dot' 'phi' 'phi_dot'};
inputs = {'u'};
outputs = {'x'; 'phi'};

sys_ss = ss(A,B,C,D,'statename',states,'inputname',inputs,'outputname',outputs)

输出

我们可以通过把状态空间表达转换为转换方程：

sys_tf = tf(sys_ss)

输出

我们发现上述转换方程里面有系数非常接近于零，这是由于Matlab在四舍五入的时候带来的小误差。我们完全可以手动把这些系数改为0。完了之后我们这里转换得到的转换方程和之前我们直接得到转换方程表征的应该是完全一样的系统

4. 引用

https://ctms.engin.umich.edu/CTMS/index.php?example=InvertedPendulum&section=SystemModeling