离散状态观测器设计与仿真

发表于 2023-12-07 分类于大学课程 Waline：本文字数： 1.3k 阅读时长 ≈ 5 分钟

系统离散化

由于计算机控制在实际中占绝对主导地位，所以离散系统非常重要。当采用零阶保持器时，连续系统到离散系统的转换方法如下图所示：

可以采用MATLAB的符号计算功能得到离散模型的参数，也可以采用c2d函数完成转换。以倒立摆为例，从连续模型得到离散模型的代码如下：

% 连续系统的参数
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_continuous = ss(A, B, C, D);
% 得到采样间隔为0.01s时的离散模型
ts = 0.01;
sys_discrete = c2d(sys_continuous, ts, 'zoh');
% 离散系统的参数
A = sys_discrete.a; B = sys_discrete.b;
C = sys_discrete.c; D = sys_discrete.d;

转换得到的离散系统由差分方程描述，标准形式如下：

\[ x(k+1)=Ax(k)+Bu(k) \tag{1} \] \[ y(k)=Cx(k)+Du(k) \tag{2} \]

接下来我们设计状态观测器，并基于状态估计值设计反馈控制器。倒立摆的状态观测器设计与仿真一文介绍了连续状态观测器的设计过程，本文设计离散状态观测器的过程和它类似，因此两篇文章可以对照着看。

能控性和能观性

在设计控制器和观测器之前，先判断系统的能控性和能观性。运行如下命令：

co = ctrb(sys_discrete);
controllability = rank(co);
ob = obsv(sys_discrete);
observability = rank(ob);

两个矩阵的秩都是4，所以系统是能控能观的。

离散观测器

本文采用如下形式的离散Luenberger状态观测器

\[ \hat{x}(k+1)=A\hat{x}(k)+Bu(k)+L\left[y(k)-\hat{y}(k)\right] \tag{3} \] \[ \hat{y}(k)=C\hat{x}(k)+Du(k) \tag{4} \]

状态估计误差\(e(k)=x(k)-\hat{x}(k)\)，所以有

\[ e(k+1)=x(k+1)-\hat{x}(k+1) \tag{5} \]

将\((1)(2)(3)(4)\)代入\((5)\)，整理可得

\[ e(k+1)=(A-LC)e(k) \tag{6} \]

如果矩阵\(A-LC\)的特征值都在单位圆之内，则\(e(k)\)会趋于零，通过配置\(A-LC\)的特征值可以决定\(e(k)\)的收敛速度。

反馈控制器

现在基于状态估计\(\hat{x}(k)\)对倒立摆进行反馈控制，系统输入信号为\(r(k)\)，采用如下形式的反馈控制律

\[ u(k)=r(k)-K\hat{x}(k) \tag{7} \]

系统状态方程为

\[ \begin{align} x(k+1)&=Ax(k)+B\left[r(k)-K\hat{x}(k)\right] \\ &=(A-BK)x(k)+BK\left[x(k)-\hat{x}(k)\right]+Br(k) \\ \end{align} \tag{8} \]

即

\[ x(k+1)=(A-BK)x(k)+BKe(k)+Br(k) \tag{9} \]

将\((6)\)和\((9)\)合并，得到如下矩阵方程

\[ \left[ \begin{array}{c} x(k+1)\\ e(k+1)\\ \end{array} \right] =\left[ \begin{matrix} A-BK& BK\\ 0& A-LC\\ \end{matrix} \right] \left[ \begin{array}{c} x(k)\\ e(k)\\ \end{array} \right] +\left[ \begin{array}{c} B\\ 0\\ \end{array} \right] r(k) \tag{10} \]

系统输出

\[ y(k)=\left[ \begin{matrix} C& 0\\ \end{matrix} \right] \left[ \begin{array}{c} x(k)\\ e(k)\\ \end{array} \right] \tag{11} \]

稳态值分析

接下来仿照倒立摆的状态反馈控制中的方法，分析系统状态\(x(k)\)和系统输出\(y(k)\)的稳态值。当\(k\to\infty\)时，有

\[ \begin{align} x(k+1)&=x(k)=x_{ss} \\ e(k)&=0 \\ \end{align} \tag{12} \]

如果输入为定值，即\(r(k)=r_{ss}\)，将\((12)\)代入到\((9)\)中可得：

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

所以

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

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

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

本例中\(D=0\)，计算可得\(C(I-A+BK)^{-1}B\)的值为

>> C*inv(eye(4)-A+B*K)*B
ans = 
   -0.0111
    0.0000

第一个分量不为零，所以，通过设置输入信号\(r(k)=r_{ss}\)的值可以决定小车位移的稳态值，摆杆角度则必然趋于零。

模型求解

通过极点配置，得到\(K\)和\(L\)矩阵。建立\((10)(11)\)所描述系统的模型并求解，代码如下：

% 极点配置，得到K和L
poles = [0.95; 0.96; 0.90+0.10j; 0.90-0.10j];
K = place(A, B, poles);
poles = [0.90; 0.91; 0.93; 0.94];
L = place(A', C', poles)';

% (10)(11)所描述的系统
Aco = [A-B*K, B*K;
       zeros(size(A)), A-L*C];
Bco = [B; zeros(size(B))];
Cco = [C, zeros(size(C))]; Dco = 0;
sys_co = ss(Aco, Bco, Cco, Dco, ts);

final_time = 2.5;     % 仿真时间
t = 0:ts:final_time;  % 时间轴
% 配置输入信号，使位移x的稳态值为0.1
rss = 0.1/(-0.0111);
r = rss*ones(size(t));
% 初始条件
x0 = [-0.2; 0; 0.1; 0]; e0 = -x0;
xe0 = [x0; e0];

% 用lsim仿真，输出保存到y，状态保存到xe
[y, t, xe] = lsim(sys_co, r, t, xe0);

% 分别绘制输出y和状态估计误差e的图像
figure(1); plot(t, y(:, 1:2), 'LineWidth', 1.5); grid on 
legend('Cart Position (m)', 'Pendulum Angle (rad)');
title('Step Response with Observer Feedback Control')
figure(2); plot(t, xe(:, 5:8), 'LineWidth', 1.5); grid on
title('State Estimate Error')

将上面所有代码按顺序贴到同一个文件中运行，可得到如下仿真结果图，控制器和观测器均满足设计要求。