用ode45做控制系统仿真

发表于 分类于 Waline： 本文字数： 1.1k 阅读时长 ≈ 4 分钟

本文是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

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

倒立摆模型仿真结果

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