倒立摆的状态反馈控制仿真

系统模型

一阶倒立摆是控制理论学习中常见的被控对象,本文使用它的近似线性模型,状态方程为

\[ \left[ \begin{array}{c} \dot{x}_1\\ \dot{x}_2\\ \dot{x}_3\\ \dot{x}_4\\ \end{array} \right] =\left[ \begin{matrix} 0& 1& 0& 0\\ 0& \frac{-\left( I+ml^2 \right) b}{I\left( M+m \right) +Mml^2}& \frac{-m^2gl^2}{I\left( M+m \right) +Mml^2}& 0\\ 0& 0& 0& 1\\ 0& \frac{mlb}{I\left( M+m \right) +Mml^2}& \frac{mgl\left( M+m \right)}{I\left( M+m \right) +Mml^2}& 0\\ \end{matrix} \right] \left[ \begin{array}{c} x_1\\ x_2\\ x_3\\ x_4\\ \end{array} \right] +\left[ \begin{array}{c} 0\\ \frac{I+ml^2}{I\left( M+m \right) +Mml^2}\\ 0\\ \frac{-ml}{I\left( M+m \right) +Mml^2}\\ \end{array} \right] u \tag{1} \]

其中\(x_1=x, x_2=\dot{x}, x_3=\theta, x_4=\dot{\theta}\)\(x\)是小车的位移,\(\theta\)是摆杆与垂直线的夹角。当\(\theta=0\)时摆杆是垂直的,所以系统平衡点(系统状态等于零)即为目标位置,这种设置更方便。假设传感器测量位移\(x\)和夹角\(\theta\),而且系统输出\(y\)由这两个量构成,则该系统在MATLAB中的表示代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
M = 0.5; m = 0.2; b = 0.1; I = 0.006; g = 9.8; l = 0.3;
den = I*(M+m)+M*m*l^2;
A = [0, 1, 0, 0;
0, -(I+m*l^2)*b/den, -(m^2*g*l^2)/den, 0;
0, 0, 0, 1;
0, (m*l*b)/den, m*g*l*(M+m)/den, 0];
B = [ 0;
(I+m*l^2)/den;
0;
-m*l/den];
C = [1, 0, 0, 0;
0, 0, 1, 0];
D = [0;
0];
sys_open = ss(A, B, C, D);

开环系统的极点即为矩阵\(A\)的特征值,可以通过运行eig(A)命令得到,可知开环系统是不稳定的。运行如下命令可得系统能控矩阵的秩为4,所以系统是能控的。

1
2
co = ctrb(sys_open);
rank(co);

状态反馈控制

假设系统输入信号为\(r(t)\),采用如下状态反馈控制律

\[ u=r-Kx \tag{2} \]

常用的反馈控制律\(u=-Kx\)则是零输入时的特殊情况。可得到闭环系统的状态方程为

\[ \dot{x}=(A-BK)x+Br \tag{3} \]

由于系统是能控的,所以闭环极点(\(A-BK\)矩阵的特征值)能够配置到任意位置。如下命令求解出状态反馈矩阵\(K\),将闭环系统极点配置到指定位置。

1
2
poles = [-5; -10; -5+10j; -5-10j];  % 根据期望性能指标选取极点
K = place(A, B, poles);

如果输入是定值,即\(r(t)=r_{ss}\)时,根据\((3)\)式可知系统状态的稳态值\(x_{ss}\)满足

\[ 0=(A-BK)x_{ss}+Br_{ss} \tag{4} \]

由此可得

\[ x_{ss} = -(A-BK)^{-1}Br_{ss} \tag{5} \]

如果\(D=0\),则有

\[ y_{ss} = -C(A-BK)^{-1}Br_{ss} \tag{6} \]

在经典控制理论中,系统是单输入单输出的,输入\(r(t)\)通常作为系统的给定值,将输出与输入比较得到误差,根据误差计算控制量使输出跟踪输入。但从上面的分析可看出,在状态空间法中,输入\(r(t)\)和控制量\(u(t)\)的维数相同,大部分情况下无法完整表示状态\(x(t)\)或输出\(y(t)\)的目标值,也就无法通过\((2)\)那种形式的控制律达到跟踪控制的效果。

本文讨论的倒立摆模型中,\(u\)是个标量(施加在小车上的力),而且可计算得到\(-C(A-BK)^{-1}B\)矩阵的值为

1
2
3
4
>> -C*inv(A-B*K)*B
ans =
-0.0071
0.0000

因此,可以结合\((6)\)式选取\(r_{ss}\)决定位移\(x\)的稳态值,但角度\(\theta\)会趋于零,无法控制到任意给定值。下面的代码将闭环系统模型表示出来,通过lsim函数得到系统输出。

1
2
3
4
5
6
7
8
9
10
Ac = A-B*K; Bc = B;
Cc = C; Dc = D;
sys_closed = ss(Ac, Bc, Cc, Dc);
t = 0:0.01:1.5;
% 位移x的稳态值为0.1时的输入信号
rss = 0.1/(-0.0071);
r = rss*ones(size(t));
x0 = [-0.2; 0; 0.1; 0]; % 系统初始状态
y = lsim(sys_closed, r, t, x0);
plot(t, y(:, 1:2), 'LineWidth', 1.5); grid on

将上面所有代码按顺序贴到同一个文件中运行,得到的输出变化曲线如下图所示,可以看到小车位移的稳态值为0.1,与预期值相符。

小车位置和摆杆角度变化曲线

参考资料

  1. Inverted Pendulum: State-Space Methods for Controller Design
  2. 现代控制工程
  3. Controller Design in Sate-Space - Engineering LibreTexts