MATLAB 状态空间设计 —— LQG/LQR 和极点配置算法

系列文章目录

前言

状态空间控制设计方法，如 LQG/LQR 和极点配置算法，适用于 MIMO 设计。

一、相关函数 —— LQG/LQR 和极点配置算法

1.1 LQR —— lqr 函数

1.1.1 函数用法

[K,S,P] = lqr(sys,Q,R,N)
[K,S,P] = lqr(A,B,Q,R,N)

说明

[K,S,P] = lqr(sys,Q,R,N) 计算连续时间或离散时间状态空间模型 sys 的最优增益矩阵
K、相关代数黎卡提方程的解 S 和闭环极点 P。Q 和 R 分别是状态和输入的权重矩阵。交叉项矩阵 N
在省略时设为零。

[K,S,P] = lqr(A,B,Q,R,N) 使用连续时间状态空间矩阵 A 和 B 计算最佳增益矩阵 K、相关黎卡提方程的解 S 以及闭环极点 P。对于离散时间模型，请使用 dlqr。

输入参数

sys - 动态系统模型，动态系统模型，以 ss 模型对象的形式指定。

A - 状态矩阵，状态矩阵，指定为 n x n 矩阵，其中 n 为状态数。

B - 输入到状态矩阵，输入 - 状态矩阵，指定为 n x m 的输入 - 状态矩阵，其中 m 为输入个数。

Q - 状态成本权重矩阵，状态-成本加权矩阵，指定为 n x n 矩阵，其中 n 为状态数。你可以使用 Bryson 规则来设置 Q 的初始值，其给定值为
Q i , i = 1 m a x i m u m   a c c e p t a b l e   v a l u e   o f   ( e r r o r s t a t e s ) 2 ,   i ∈ { 1 , 2 , . . . , n } Q_{i,i}=\dfrac{1}{\mathrm{maximum~acceptable~value~of~(error_{states})^{2}}},~\mathrm{i}\in\{1,2,...,n\} Qi,i​=maximum acceptable value of (errorstates​)21​, i∈{1,2,...,n}
$$Q=\left[\begin{array}{cccc}{Q}_{1,1}& 0& \cdots & 0\\ 0& Q_{2,2}& \cdots & 0\\ 0& 0& \ddots & ⋮\\ 0& 0& \cdots & Q_{n,n}\end{array}\right]$$

R - 输入成本权重矩阵，输入成本加权矩阵，指定为标量或与 D'D 相同大小的矩阵。这里，D 是馈通状态空间矩阵。可以使用 Bryson 规则设置 R 的初始值，其给定值为
R i , i = 1 m a x i m u m   a c c e p t a b l e   v a l u e   o f   ( e r r o r s t a t e s ) 2 ,   i ∈ { 1 , 2 , . . . , m } R_{i,i}=\dfrac{1}{\mathrm{maximum~acceptable~value~of~(error_{states})^{2}}},~\mathrm{i}\in\{1,2,...,m\} Ri,i​=maximum acceptable value of (errorstates​)21​, i∈{1,2,...,m}
$$R=\left[\begin{array}{cccc}{R}_{1,1}& 0& \cdots & 0\\ 0& R_{2,2}& \cdots & 0\\ 0& 0& \ddots & ⋮\\ 0& 0& \cdots & R_{n,n}\end{array}\right]$$
这里，m 是输入的个数。

N - 可选的交叉项矩阵，可选的交叉项矩阵，以矩阵形式指定。如果没有指定 N，lqr 默认将 N 设为 0。

输出参数

K - 最优增益，闭环系统的最优增益，以大小为 n 的行向量形式返回，其中 n 为状态数。

S - 相关代数黎卡提方程的解，相关代数黎卡提方程的解，以 n x n 矩阵形式返回，其中 n 为状态数。换句话说，S
的维度与状态空间矩阵 A 相同。更多信息，请参见 icare 和 idare。

P - 闭环系统的极点，闭环系统的极点，以大小为 n 的列向量形式返回，其中 n 为状态数。

局限性

输入数据必须满足以下条件：

• 一对矩阵 A 和 B 必须是可稳定的。

• $[Q,N;N^{\prime },R]$ 必须是非负定。

• $R>0$ 并且 $Q-N{R}^{-1}{N}^T\ge 0.$

• $(Q-N{R}^{-1}{N}^T,A-B{R}^{-1}{N}^T)$ 在虚轴（或离散时间中的单位圆）上没有不可观测的模式。

提示： lqr 支持非奇异 E 的描述符模型。lqr 的输出 S 是等效显式状态空间模型的代数黎卡提方程的解： $$\frac{dx}{dt}=E^{-1}Ax+E^{-1}Bu$$

LQR 算法

对于连续时间系统，lqr 计算的状态反馈控制 u=-Kx 可使二次成本函数最小化
$$J(u)={\int }_0^{\infty }(x^{T}Qx+u^{T}Ru+2x^{T}Nu)dt$$ 系统动力学
$\dot{x}=Ax+Bu.$

除了状态反馈增益 K 之外，lqr 还返回相关代数黎卡提方程的解 S $$A^{T}S+SA-(SB+N)R^{-1}(B^{T}S+N^T)+Q=0$$ 和闭环极点 P = eig(A-BK)。增益矩阵 K 由 S 得出，计算公式为
$$K=R^{-1}(B^{T}S+N^T).$$ 对于离散时间系统，lqr 可计算状态反馈控制 $u_n=-K{x}_n$
最小化 J = ∑ n = 0 ∞ { x T Q x + u T R u + 2 x T N u } J=\sum_{n=0}^{\infty}\{x^{T}Q x+u^{T}R u+2x^{T}N u\} J=n=0∑∞​{xTQx+uTRu+2xTNu} 系统动力学
$x_{n+1}=A{x}_n+B{u}_n$

在任何情况下，如果省略交叉项矩阵 N，lqr 都会将 N 设为 0。

1.1.2 举例

1.1.2.1 倒摆模型的 LQR 控制

pendulumModelCart.mat 包含小车上倒立摆的状态空间模型，其输出为小车位移 x 和摆角 θ，控制输入 u 为小车上的水平力。
KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ \begin{bmatrix…
KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ \bm{y}=\begin{…
首先，将状态空间模型 sys 加载到工作区。

load('pendulumCartModel.mat','sys')

由于输出为 $x$ 和 $\theta$，而输入只有一个，因此使用 Bryson 规则确定 Q 和 R。

Q = [1,0,0,0;...
    0,0,0,0;...
    0,0,1,0;...
    0,0,0,0];
R = 1;

使用 lqr 求增益矩阵 K。由于没有指定 N，lqr 将 N 设为 0。

[K,S,P] = lqr(sys,Q,R)

K = 1×4

   -1.0000   -1.7559   16.9145    3.2274

S = 4×4

    1.5346    1.2127   -3.2274   -0.6851
    1.2127    1.5321   -4.5626   -0.9640
   -3.2274   -4.5626   26.5487    5.2079
   -0.6851   -0.9640    5.2079    1.0311

P = 4×1 complex

  -0.8684 + 0.8523i
  -0.8684 - 0.8523i
  -5.4941 + 0.4564i
  -5.4941 - 0.4564i

虽然 Bryson 规则通常能提供令人满意的结果，但它通常只是根据设计要求调整闭环系统响应的试错迭代设计程序的起点。

1.2 LQG —— lqg() 函数

1.2.1 函数用法

reg = lqg(sys,QXU,QWV)
reg = lqg(sys,QXU,QWV,QI)
reg = lqg(sys,QXU,QWV,QI,'1dof')
reg = lqg(sys,QXU,QWV,QI,'2dof')
reg = lqg(___,'current')
[reg,info] = lqg(___)

说明

reg = lqg(sys,QXU,QWV) 给定一个被控对象的状态空间模型 sys 以及权重矩阵 QXU 和
QWV，计算出一个最优线性二次高斯（LQG）调节器 reg。动态调节器 reg 利用测量值 y 生成控制信号 u，将
y 调节到零值附近。使用正反馈将该调节器与被控对象的输出 y 连接起来。

LQG 调节器使成本函数最小化
$$J=E\left\{*{lim}_{\tau \to \infty }\frac{1}{\tau }{\int }_0^{\tau }[x^T,u^T]Qxu\left[\begin{array}{c}x\\ u\end{array}\right]dt\right\}$$

受被控对象方程的约束

$$\begin{array}{c}dx/dt=Ax+Bu+w\\ y=Cx+Du+\nu \end{array}$$

其中，过程噪声 w 和测量噪声 v 均为具有协方差的高斯白噪声：

KaTeX parse error: Undefined control sequence: \cal at position 95: …rime}]\right)={\̲c̲a̲l̲ ̲Q}{\cal W}V

reg = lqg(sys,QXU,QWV,QI) 使用设定点指令 r 和测量值 y 来生成控制信号 u。

LQG 伺服控制器可使成本函数最小化

$$J=E\left\{*{lim}_{\tau \to \infty }\frac{1}{\tau }{\int }_0^{\tau }\left([x^T,u^T]Qxu{|}_u^x|+x_i^T{Q}_i{x}_i\right)dt\right\}$$

对于多输入多输出系统，r、y 和 $x_i$ 的长度必须相同。

reg = lqg(sys,QXU,QWV,QI,'1dof')计算以 e = r - y 而不是 [r ; y]
作为输入的一自由度伺服控制器。

reg = lqg(sys,QXU,QWV,QI,'2dof') 等价于
LQG(sys,QXU,QWV,QI)，并产生前面显示的二自由度伺服控制器。

reg = lqg(____,'current')使用 "current "卡尔曼估计器，该估计器在计算离散时间系统的 LQG
调节器时使用 x[n|n] 作为状态估计值。

[reg,info] = lqg(____)返回前面任何语法结构 info
中的控制器和估计器增益矩阵。例如，您可以使用控制器和估计器增益以观测器形式实现控制器。更多信息，请参阅算法。

1.2.2 举例

线性-二次高斯 (LQG) 调节器和伺服控制器设计

本例介绍如何为以下系统设计线性二次高斯（LQG）调节器、一自由度 LQG 伺服控制器和二自由度 LQG 伺服控制器。

被控对象有三个状态 (x)、两个控制输入 (u)、三个随机输入 (w)、一个输出 (y)、输出的测量噪声 (v)，以及以下状态方程和测量方程。

$$\begin{gathered} \frac{dx}{dt}=Ax+Bu+w \\ y=Cx+Du+\nu \end{gathered}$$

$$\begin{array}{rl}A=\left[\begin{array}{lll}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right]& B=\left[\begin{array}{ll}0.3& 1\\ 0& 1\\ -0.3& 0.9\end{array}\right]\end{array}$$$$C=\left[1.91.31\right]D=\left[0.53-0.61\right]$$

系统的噪声协方差数据如下

$Qn=E(w{w}^T)=\left[\begin{array}{ccc}4& 2& 0\\ 2& 1& 0\\ 0& 0& 1\end{array}\right]$

$R_n=E(\nu {\nu }^T)=0.7$

对于调节器，使用以下成本函数来定义调节性能和控制输入之间的权衡：

$J(u)={\int }_0^{\infty }(0.1x^{T}x+u^T\left[\begin{array}{cc}1& 0\\ 0& 2\end{array}\right]u)dt$

对于伺服控制器，使用以下成本函数来定义跟踪器性能和控制力之间的权衡：

$J(u)={\int }_0^{\infty }(0.1x^{T}x+x_i^2+u^T\left[\begin{array}{cc}1& 0\\ 0& 2\end{array}\right]u)dt$

为该系统设计 LQG 控制器：

1. 在 MATLAB 命令窗口中输入以下内容，创建状态空间系统：

A = [0 1 0;0 0 1;1 0 0];    
B = [0.3 1;0 1;-0.3 0.9];
C = [1.9 1.3 1];  
D = [0.53 -0.61];
sys = ss(A,B,C,D);

2. 输入以下命令，定义噪声协方差数据和加权矩阵：

nx = 3;    %Number of states
ny = 1;    %Number of outputs
Qn = [4 2 0; 2 1 0; 0 0 1];
Rn = 0.7;
R = [1 0;0 2]
QXU = blkdiag(0.1*eye(nx),R);
QWV = blkdiag(Qn,Rn);
QI = eye(ny);

3. 键入以下命令，组建 LQG 调节器：

KLQG = lqg(sys,QXU,QWV)

This command returns the following LQG regulator:
A = 
           x1_e    x2_e    x3_e
   x1_e  -6.212  -3.814  -4.136
   x2_e  -4.038  -3.196  -1.791
   x3_e  -1.418  -1.973  -1.766
 
B = 
             y1
   x1_e   2.365
   x2_e   1.432
   x3_e  0.7684
 
C = 
            x1_e       x2_e       x3_e
   u1   -0.02904  0.0008272     0.0303
   u2    -0.7147    -0.7115    -0.7132
 
D = 
       y1
   u1   0
   u2   0
 
Input groups:              
       Name        Channels
    Measurement       1    
                           
Output groups:             
      Name      Channels   
    Controls      1,2      
                           
Continuous-time model.

4. 键入以下命令，形成单自由度 LQG 伺服控制器：

KLQG1 = lqg(sys,QXU,QWV,QI,'1dof')

This command returns the following LQG servo controller:
A = 
           x1_e    x2_e    x3_e     xi1
   x1_e  -7.626  -5.068  -4.891  0.9018
   x2_e  -5.108  -4.146  -2.362  0.6762
   x3_e  -2.121  -2.604  -2.141  0.4088
   xi1        0       0       0       0
 
B = 
              e1
   x1_e   -2.365
   x2_e   -1.432
   x3_e  -0.7684
   xi1         1
 
C = 
          x1_e     x2_e     x3_e      xi1
   u1  -0.5388  -0.4173  -0.2481   0.5578
   u2   -1.492   -1.388   -1.131   0.5869
 
D = 
       e1
   u1   0
   u2   0
 
Input groups:           
    Name     Channels   
    Error       1       
                        
Output groups:          
      Name      Channels
    Controls      1,2   
                        
Continuous-time model.

5. 输入以下命令，生成二自由度 LQG 伺服控制器：

KLQG2 = lqg(sys,QXU,QWV,QI,'2dof')

This command returns the following LQG servo controller:
A = 
           x1_e    x2_e    x3_e     xi1
   x1_e  -7.626  -5.068  -4.891  0.9018
   x2_e  -5.108  -4.146  -2.362  0.6762
   x3_e  -2.121  -2.604  -2.141  0.4088
   xi1        0       0       0       0
 
B = 
             r1      y1
   x1_e       0   2.365
   x2_e       0   1.432
   x3_e       0  0.7684
   xi1        1      -1
 
C = 
          x1_e     x2_e     x3_e      xi1
   u1  -0.5388  -0.4173  -0.2481   0.5578
   u2   -1.492   -1.388   -1.131   0.5869
 
D = 
       r1  y1
   u1   0   0
   u2   0   0
 
Input groups:              
       Name        Channels
     Setpoint         1    
    Measurement       2    
                           
Output groups:             
      Name      Channels   
    Controls      1,2      
                           
Continuous-time model.

小贴士

• lqg 可用于连续时间和离散时间被控对象。在离散时间情况下，lqg 默认使用 x[n|n-1] 作为状态估计值。要使用 x[n|n] 作为状态估计并计算最优 LQG 控制器，请使用 "current "输入参数。有关状态估计器的详细信息，请参见
  kalman。

• 计算 LQG调节器时，lqg 使用 lqr 和 kalman 命令。要计算伺服控制器，lqg 使用 lqi 和 kalman 命令。

• 如果希望更灵活地设计调节器，可以使用 lqr、kalman 和 lqgreg 命令。在设计伺服控制器时，如果需要更大的灵活性，可以使用 lqi、kalman 和 lqgtrack
  命令。有关使用这些命令以及如何决定何时使用这些命令的更多信息，请参阅线性-二次方-高斯 (LQG)
  调节设计和带积分动作的伺服控制器的线性-二次方-高斯 (LQG) 设计。

LQG 算法

控制器方程为

• 对于连续时间： $$\begin{gathered} d{x}_e=A{x}_e+Bu+L(y-C{x}_e-Du) \\ u=-K{x}_e-K_i{x}_i \end{gathered}$$
• 对于离散时间： $x[n+1|n]=Ax[n|n-1]+Bu[n]+L(y[n]-Cx[n]n-1]-Du[n])$
  • 延迟估计器： $u[n]=-K_{x}x[n]n-1]-K_i{x}_i[n]$
  • 目前的估算器： $u[n]=-K_{x}x[n]n]-K_i{x}_i[n]-K_{w}w[n]=-K_{x}x[n]n-1]-K_i{x}_i[n]-(K_x{M}_x+K_w{M}_w)y_{inn}[n]$ $y_{inn}[n]=y[n]-Cx[n]n-1]-Du[n]$

这里

• A、B、C 和 D 是 LQG 调节器的状态空间矩阵，reg.

• $x_i$ 是跟踪误差 r - y 的积分。

• $K_x$、$K_w$、$K_i$、L、$M_x$ 和 $M_w$ 是信息中返回的控制器和估计器增益矩阵。

1.3 极点配置 —— place() 函数

1.3.1 函数用法

K = place(A,B,p)
[K,prec] = place(A,B,p)

说明

极点配置是一种计算最优增益矩阵的方法，用于将闭环极点分配给指定位置，从而确保系统稳定性。闭环极点位置会直接影响上升时间、稳定时间和瞬变振荡等时间响应特性。有关详细信息，请参阅极点配置。

从图中，假设有以下状态空间形式的线性动态系统：
$$\begin{array}{c}\dot{x}=Ax+Bu\\ \\ y=Cx+Du\end{array}$$
对于期望的自共轭闭环极点位置的给定向量 p，place 计算增益矩阵 K，使得状态反馈 u = –Kx 将极点配置在位置 p。换句话说，A - BK 的特征值将匹配 p 的条目（取决于排序）。

K = place(A,B,p) 通过计算状态反馈增益矩阵 K，配置所需的闭环极点 p。被控对象的所有输入都假定为控制输入。
place也适用于多输入系统，并且基于 [1] 中的算法。此算法使用额外的自由度来求一个解，以使闭环极点对于 A 或 B 中的扰动具有最小的敏感度。

[K,prec] = place(A,B,p) 还返回 prec，用于精确估计 A - BK 的特征值与指定位置 p 的匹配程度（prec 可计算实际闭环极点中的精确小数位数）。如果某个非零闭环极点偏离期望位置超出 10%，则系统会发出警告。

输入参数

A — 状态矩阵 状态矩阵，指定为一个 Nx×Nx 矩阵，其中 Nx 是状态数。

B — 输入-状态矩阵 输入-状态矩阵，指定为 Nx×Nu 矩阵，其中 Nx 是状态数，Nu 是输入数。

p — 闭环极点位置 闭环极点位置，指定为长度为 Nx 的向量，其中 Nx 是状态数。换句话说，p 的长度必须与 A 的行大小匹配。闭环极点位置会直接影响上升时间、稳定时间和瞬变振荡等时间响应特性。有关选择极点的示例，请参阅二阶系统的极点配置设计。

如果 p 中某些极点的重数大于 rank(B)，则 place 返回错误。

在高阶问题中，选择某些极点位置会导致增益非常大。大增益会带来敏感性问题，这表明在使用极点配置方法时要小心。有关数值测试的结果，请参阅 [2]。

输出参数

K — 最优增益 最优增益或全状态反馈增益，以 Ny×Nx 矩阵形式返回，其中 Nx 是状态数，Ny 是输出数。place 计算增益矩阵
K，使得状态反馈 u = -Kx 将闭环极点配置于位置 p。

当矩阵 A 和 B 为实数时，则 K 为：

• 实数，前提是 p 具有自共轭性。

• 复数，前提是极点位置不具有复共轭性。

prec — 指定极点的准确性估计值 指定极点的准确性估计值，以标量形式返回。prec 对比 p
中指定的极点位置来计算实际闭环极点的精确小数位数。

1.3.2 示例

1.3.2.1 二阶系统的极点配置设计

对于此示例，假设有一个具有以下状态空间矩阵的简单二阶系统：
$$A=\left[\begin{array}{ll}-1& -2\\ 1& 0\end{array}\right]B=\left[\begin{array}{l}2\\ 0\end{array}\right]C=\left[01\right]D=0$$
输入矩阵并创建状态空间系统。

A = [-1,-2;1,0];
B = [2;0];
C = [0,1];
D = 0;
sys = ss(A,B,C,D);

计算开环极点并检查开环系统的阶跃响应。

Pol  = pole(sys)

Pol = 2×1 complex

  -0.5000 + 1.3229i
  -0.5000 - 1.3229i

figure(1)
step(sys)
hold on;

请注意，生成的系统为欠阻尼系统。因此，选择复平面左半部分的实极点来消除振荡。

p = [-1,-2];

使用极点配置求增益矩阵 K，并检查 syscl 的闭环极点。

K = place(A,B,p);
Acl = A-B*K;
syscl = ss(Acl,B,C,D);
Pcl = pole(syscl)

Pcl = 2×1

   -2.0000
   -1.0000

现在，比较闭环系统的阶跃响应。

figure(1)
step(syscl)

因此，使用极点配置获得的闭环系统是稳定的，具有良好的稳态响应。

请注意，选择远离虚轴的极点可以缩短响应时间，但会减小系统的稳态增益。例如，假设上述系统使用极点 [-2,-3]。

p = [-2, -3];
K2 = place(A,B,p);
syscl2 = ss(A-B*K2,B,C,D);
figure(1);
step(syscl2);

stepinfo(syscl)

ans = struct with fields:
         RiseTime: 2.5901
    TransientTime: 4.6002
     SettlingTime: 4.6002
      SettlingMin: 0.9023
      SettlingMax: 0.9992
        Overshoot: 0
       Undershoot: 0
             Peak: 0.9992
         PeakTime: 7.7827

stepinfo(syscl2)

ans = struct with fields:
         RiseTime: 1.4130
    TransientTime: 2.4766
     SettlingTime: 2.4766
      SettlingMin: 0.3003
      SettlingMax: 0.3331
        Overshoot: 0
       Undershoot: 0
             Peak: 0.3331
         PeakTime: 4.1216

1.3.2.2 极点配置观测器设计

对于此示例，假设有以下 SISO 状态空间模型：

$A=\left[\begin{array}{ll}-1& -0.75\\ 1& 0\end{array}\right]B=\left[\begin{array}{l}1\\ 0\end{array}\right]C=\left[11\right]D=0$

创建由以下状态空间矩阵定义的 SISO 状态空间模型：

创建由以下状态空间矩阵定义的 SISO 状态空间模型：

现在，向被控对象提供一个脉冲，并使用 lsim 对其进行仿真。绘制输出。

N = 250;
t = linspace(0,25,N);
u = [ones(N/2,1); zeros(N/2,1)];
x0 = [1;2];
[y,t,x] = lsim(Plant,u,t,x0);

figure
plot(t,y);
title('Output');

对于此示例，假设所有状态变量都无法测量，只有输出才能测量。因此，使用这种测量方法设计一个观测器。使用 place 来计算估计器增益，方法是转置 A 矩阵，并用 C’ 代换矩阵 B。对于此实例，将所需的极点位置选为 -2 和 -3。

L = place(A',C',[-2,-3])';

使用估计器增益，依据对偶/分离原理代换状态矩阵，并创建估计的状态空间模型。

At = A-L*C;
Bt = [B,L];
Ct = [C;eye(2)];
sysObserver = ss(At,Bt,Ct,0);

使用相同的脉冲输入对系统的时间响应进行仿真。

[observerOutput,t] = lsim(sysObserver,[u,y],t);
yHat = observerOutput(:,1);
xHat = observerOutput(:,[2 3]);

比较实际系统和估计系统的响应。

figure;
plot(t,x);
hold on;
plot(t,xHat,'--');
legend('x_1','x_2','xHat_1','xHat_2')
title('Comparison - Actual vs. Estimated');