用ode45做控制系统仿真

本文是reddit r/ControlTheory一篇帖子的总结。这篇帖子以倒立摆为例,介绍如何在MATLAB中用解微分方程的方式做控制系统仿真,并与常规仿真方法做比较。

所用模型是倒立摆在平衡点附近的近似线性模型,这个模型很常见,从代码中也能很容易看出来,所以我就不列出。采用反馈控制方式,将测量输出\(y\)(本例中就是系统的状态,\(y=Cx\)\(C\)矩阵为单位阵)与参考输入\(r\)做比较,得到差值,将差值乘以反馈矩阵\(K\)后,作为倒立摆模型的输入,即 \[ u = K(r-y) \]

控制本质上要求解微分方程,MATLAB解微分方程最常用的是ode45函数,它求解的微分方程的标准形式为\(y'=f(t,y)\),右边的表达式是时间\(t\)和变量\(y\)的函数。我以前被这个表达式迷惑了,以为它只能求解没有外部输入的非强迫系统(unforced system)。其实不然,右边的表达式可以非常灵活,能够求解有外部输入信号的控制系统。全部代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
% 倒立摆近似线性模型的参数
m = 1; M = 5; L = 2; g = 10;
A = [0, 1, 0, 0;
0, -1/M, m*g/M, 0;
0, 0, 0, 1;
0, -1/(M*L), (m+M)*g/(M*L), 0];
B = [0; 1/M; 0; 1/(M*L)];

x0 = [0; 0; pi+0.1; 0.5]; % 初始条件
tspan = 0:0.001:10; % 仿真时间
r = [0; 0; pi; 0]; % 参考输入

% 配置极点,得到状态反馈矩阵K
p = [-1.5; -1.6; -1.7; -1.8];
K = place(A, B, p);

% 输入给系统的控制信号为u = K*(r-y),将它作为模型的参数
[t, y] = ode45(@(t, y) model(y, m, M, L, g, K*(r-y)), tspan, x0);

% 绘制各个状态值随时间变化的曲线
x = y(:, 1); x_dot = y(:, 2);
theta = y(:, 3); theta_dot = y(:, 4);
figure(1); hold on
plot(t, x, 'r', t, x_dot, 'g'); grid on
plot(t, theta, 'm', t, theta_dot, 'k'); grid on

% 倒立摆的微分方程模型,形式可以非常灵活
function dydt = model(y, m, M, L, g, u)
dydt = zeros(4, 1);
dydt(1) = y(2);
dydt(2) = (-1/M)*y(2) + (m*g/M)*(y(3)-pi) + (1/M)*u;
dydt(3) = y(4);
dydt(4) = (-1/(M*L))*y(2) + ((m+M)*g/(M*L))*(y(3)-pi) + (1/(M*L))*u;
end

最后几行代码是倒立摆微分方程模型的描述,将输入信号u作为model函数的参数。第18行调用ode45函数求解微分方程,在这里将输入信号的计算式u = K*(r-y)传递给model函数。

MATLAB中对这样的系统做仿真一般用lsim函数,因为它可以自定义输入信号,还可以包含初始条件。由于该模型是平衡点附近的近似线性模型,控制目标是将系统输出(或系统状态)保持在参考输入\(r\)附近,它不是零点,因此模型变成如下形式: \[ \dot{x} = A(x-r) + Bu \]

\(u = K(r-y)=K(r-Cx)\)代入得: \[ \dot{x} = (A-BKC)x+(BK-A)r \]

本文使用lsim函数的如下接口:

1
[y, tOut, x] = lsim(sys, u, t, x0)

\(r\)即为lsim的参数u\(BK-A\)是模型中与输入信号相乘的\(B\)矩阵。因此,得到如下仿真代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
% 倒立摆近似线性模型的参数
m = 1; M = 5; L = 2; g = 10;
A = [0, 1, 0, 0;
0, -1/M, m*g/M, 0;
0, 0, 0, 1;
0, -1/(M*L), (m+M)*g/(M*L), 0];
B = [0; 1/M; 0; 1/(M*L)];
C = eye(size(A));

% 配置极点,得到状态反馈矩阵K
p = [-1.5; -1.6; -1.7; -1.8];
K = place(A, B, p);

x0 = [0; 0; pi+0.1; 0.5]; % 初始条件
tspan = 0:0.001:10; % 仿真时间
r = [0; 0; pi; 0]; % 参考输入

% 输入信号的维数必须满足要求
r_t = r * ones(size(tspan));
% 闭环系统的模型
sys_cl = ss(A-B*K*C, B*K-A, C, 0);
[y, tspan, x] = lsim(sys_cl, r_t, tspan, x0);

% 绘制系统状态变化曲线
figure(2); hold on
plot(tspan, y(:, 1), 'r', tspan, y(:, 2), 'g'); grid on
plot(tspan, y(:, 3), 'm', tspan, y(:, 4), 'k'); grid on

可以验证,用ode45lsim得到的结果完全一样,系统状态随时间变化的曲线如下图所示,达到了控制目标。

倒立摆模型仿真结果

对比来看,lsim是MATLAB中专门用于控制系统仿真的函数,使用更方便,ode45则更通用,可以做更复杂更灵活的仿真,如事件触发控制(event-triggered control)。