免费建站团队,wordpress 图片 宽 高,电子商务网站建设作业代码,河南网站建设报价目录 前言输出反馈MPC将要用到的基础概念 参考#xff1a;《现代控制理论#xff08;第3版#xff09; (刘豹,唐万生) 》【基础】系统能控性【基础】系统能观性【基础】能控性和能观性的对偶关系 一、观测器设计1.1 线性观测器设计1.2 误差#xff08;稳定性#xff09;分… 目录 前言输出反馈MPC将要用到的基础概念 参考《现代控制理论第3版 (刘豹,唐万生) 》【基础】系统能控性【基础】系统能观性【基础】能控性和能观性的对偶关系 一、观测器设计1.1 线性观测器设计1.2 误差稳定性分析1.3 MATLAB实例1.3.1 极点配置1.3.2 LQR设计最优增益 二、无约束输出反馈MPC2.1 稳定性分析2.2 MATLAB实例 前言
致谢 【模型预测控制2022春lecture 3-1 Output feedback MPC】
输出反馈MPC
针对更通用的系统 x k 1 A x k B u k y k C x k (1) \begin{align*} x_{k1} Ax_k Bu_k \\ y_k Cx_k \end{align*} \tag{1} xk1ykAxkBukCxk(1) 其中 x k ∈ R n x_k \in \mathbb{R}^n xk∈Rn 是系统的全状态 u k ∈ R p u_k \in \mathbb{R}^p uk∈Rp 是系统的输入 y k ∈ R m y_k \in \mathbb{R}^m yk∈Rm 是可知的系统输出。 当 m n mn mn 时说明部分状态不可测量需要增加观测器来观测系统的全状态此前提是系统是可观测的能观性。 由此即可通过系统的输出反馈来设计 MPC 控制器。
将要用到的基础概念 参考《现代控制理论第3版 (刘豹,唐万生) 》
【基础】系统能控性
定义若存在连续的输入 u ( t ) u(t) u(t)可在有限时间 [ t 0 , t f ] [t_0, ~ t_f] [t0, tf] 内使系统从状态 x ( t 0 ) x(t_0) x(t0)到达任一终端状态 x ( t f ) x(t_f) x(tf)则称该状态是能控的。若系统每一状态都能控则称系统是能观的。
简单理解就是系统中每个子状态都独立地直接或间接地受输入影响则可控。
【基础】系统能观性
定义若系统在有限观测时间 t f t 0 t_f t_0 tft0 内可根据 [ t 0 , t f ] [t_0, ~ t_f] [t0, tf] 期间输出的 y ( t ) y(t) y(t)唯一地确定系统状态 x ( t 0 ) x(t_0) x(t0)则称该状态是能观的。若系统每一状态都能观测则称系统是能观的。
针对离散系统的简单理解就是 将系统 k k k 时刻的全状态 x k x_k xk 分为两块记 x k x k ∪ x k − x_k x_k^ \cup x_k^- xkxk∪xk− 系统输出 y k y_k yk 可能只与部分的系统状态 x k x_k^ xk 存在联系若系统中与 y k y_k yk没有直接关系的状态 x k − x_k^- xk−其与 x k x_k^ xk 存在联系 就可通过有限的多时刻的系统输出 y k , y k 1 , ⋯ y_k,y_{k1},\cdots yk,yk1,⋯建立足够多的独立方程求解出系统的全状态 x k x_k xk. 故 x k − x_k^- xk−与 x k x_k^ xk 是否存在联系决定了系统的能观性。
《现代控制理论第3版 (刘豹,唐万生) 》中给出了许多系统能控/观性的判定方法。
【基础】能控性和能观性的对偶关系 对偶关系定义 系统 ∑ 1 \sum _1 ∑1 x ˙ 1 A 1 x 1 B 1 u 1 y ˙ 1 C 1 x 1 \begin{align*} \dot{x}_1 A_1 x_1 B_1 u_1 \\ \dot{y}_1 C_1 x_1 \end{align*} x˙1y˙1A1x1B1u1C1x1 系统 ∑ 2 \sum _2 ∑2 x ˙ 2 A 2 x 2 B 2 u 2 y ˙ 2 C 2 x 2 \begin{align*} \dot{x}_2 A_2 x_2 B_2 u_2 \\ \dot{y}_2 C_2 x_2 \end{align*} x˙2y˙2A2x2B2u2C2x2 若 A 2 A 1 T , B 2 C 1 T , C 2 B 1 T A_2 A_1^T, \quad B_2 C_1^T, \quad C_2 B_1^T A2A1T,B2C1T,C2B1T 则称系统 ∑ 1 \sum _1 ∑1 和 ∑ 2 \sum _2 ∑2 是互为对偶的。
对偶原理 系统 ∑ 1 \sum _1 ∑1 和 ∑ 2 \sum _2 ∑2 是互为对偶则 ∑ 1 \sum _1 ∑1 的能控性等价于 ∑ 2 \sum _2 ∑2 的能观性 ∑ 1 \sum _1 ∑1 的能观性等价于 ∑ 2 \sum _2 ∑2 的能控性。
应用 在设计观测器的观测增益 L L L 时可通过对偶的控制系统使用控制器设计方法设计最优反馈增益 K K K有 L K T L K^T LKT. 一、观测器设计
1.1 线性观测器设计
针对系统 (1)设 k k k 时刻观测的状态为 x ^ k \hat{x}_{k} x^k线性观测器可设计为 x ^ k 1 A x ^ k B u k H ( y k − C x ^ k ) (2) \hat{x}_{k1} A\hat{x}_{k} Bu_k H(y_k - C\hat{x}_{k}) \tag{2} x^k1Ax^kBukH(yk−Cx^k)(2) 上式可称为 观测器状态方程其中 H H H 为观测增益.
1.2 误差稳定性分析
定义误差为 x ~ k x k − x ^ k \tilde{x}_k x_k - \hat{x}_{k} x~kxk−x^k由系统状态方程-观测器状态方程得 x ~ k 1 A x ~ k − L ( y k − C x ^ k ) \tilde{x}_{k1} A\tilde{x}_k - L(y_k - C\hat{x}_{k}) x~k1Ax~k−L(yk−Cx^k) 将输出方程代入上式可整理得 x ~ k 1 ( A − L C ) x ~ k (3) \tilde{x}_{k1} (A-LC)\tilde{x}_k \tag{3} x~k1(A−LC)x~k(3) 由李雅普诺夫稳定性间接法可知当 ∣ e i g ( A − L C ) ∣ 1 |\mathrm{eig}(A-LC)|1 ∣eig(A−LC)∣1 时误差将趋于 0. L L L 可通过其对偶系统来设计。 式 (3) 可以写为 x ~ k 1 A ′ x ~ k \tilde{x}_{k1} A^\prime \tilde{x}_k x~k1A′x~k其中 A ′ ( A − L C ) A^\prime (A-LC) A′(A−LC). 对偶系统为 x ~ k 1 ( A − L C ) T x ~ k ⇒ x ~ k 1 A T x ~ k − C T L T x ~ k \begin{align*} \tilde{x}_{k1} (A-LC)^T \tilde{x}_k\\ \Rightarrow \tilde{x}_{k1} A^T\tilde{x}_k - C^T L^T \tilde{x}_k \end{align*} x~k1⇒x~k1(A−LC)Tx~kATx~k−CTLTx~k 以上形式就转为了我们熟悉的线性反馈控制系统的形式 x k 1 A x k − B K x k x_{k1} Ax_k - BKx_k xk1Axk−BKxk. 我们可用极点配置或最优控制来设计增益使系统渐近稳定. 因为 e i g ( A ′ ) e i g ( ( A ′ ) T ) \mathrm{eig}(A^\prime) \mathrm{eig}((A^\prime)^T) eig(A′)eig((A′)T)所以当控制系统渐近稳定观测器误差也将收敛为0.
1.3 MATLAB实例
针对式 (2) 形式的观测器需要使式 (3) 系统的状态稳定为 0 \mathbf{0} 0. 设系统中 A [ 1.1 2 0 0.95 ] A \begin{bmatrix} 1.1 2 \\ 0 0.95 \end{bmatrix} A[1.1020.95] B [ 0 0.079 ] B \begin{bmatrix} 0 \\ 0.079 \end{bmatrix} B[00.079] C [ 1 0 ] C \begin{bmatrix} 1 0 \end{bmatrix} C[10]. 设计观测增益 L L L使 ∣ e i g ( A − L C ) ∣ 1 |\mathrm{eig}(A-LC)|1 ∣eig(A−LC)∣1即可使观测器误差收敛为0.
下面将根据对偶原理使用极点配置法和基于对偶原理的LQR方法来设计. 式 (3) 的对偶系统为 x ~ k 1 A T x ~ k − C T L T x ~ k \tilde{x}_{k1} A^T\tilde{x}_k - C^T L^T \tilde{x}_k x~k1ATx~k−CTLTx~k 对应控制系统 x ~ k 1 ′ A ′ x ~ k ′ − B ′ K ′ x ~ k ′ \tilde{x}_{k1}^\prime A^\prime\tilde{x}_k^\prime - B^\prime K^\prime \tilde{x}_k^\prime x~k1′A′x~k′−B′K′x~k′ 其中 A ′ A T A^\prime A^T A′AT B ′ C T B^\prime C^T B′CT K ′ L T K^\prime L^T K′LT.
1.3.1 极点配置
A [1.1 2;0 0.95];
B [0; 0.079];
C [1, 0];Ap A;
Bp C;%% 极点配置
p [0.5 0.5j, 0.5 - 0.5j]; % 适当选取极点
K place(Ap, Bp, p);L K;
disp(eig(A - L*C))
disp(abs(eig(A - L*C)))得 L [ 1.0500 0.2263 ] L \begin{bmatrix} 1.0500 \\ 0.2263 \end{bmatrix} L[1.05000.2263].
1.3.2 LQR设计最优增益
LQR相关内容可参考 离散LQR原理 .
A [1.1 2;0 0.95];
B [0; 0.079];
C [1, 0];Ap A;
Bp C;%% LQR 最优增益
K LQR(Ap, Bp, eye(2), 0.1, 500, 1e-6);
L K;
disp(eig(A - L*C))
disp(abs(eig(A - L*C)))
%%
function K LQR(A, B, Q, R, maxIter, eps)
% A、B分别为系统矩阵和输入矩阵Q和R分别为状态误差和输入的对角权重矩阵
% maxIter为最大迭代步数Neps为迭代精度Ci 1; P Q; delta 1e9;while i maxIter delta epsPn Q A * (P - P*B* inv(RB*P*B) *B*P) * A;delta max(abs(Pn-P), [], all);P Pn;i i1;endK inv(R B * P * B) * B * P * A;
end得 L [ 1.8342 0.3571 ] L \begin{bmatrix} 1.8342 \\ 0.3571 \end{bmatrix} L[1.83420.3571].
二、无约束输出反馈MPC
针对系统 x k 1 A x k B u k y k C x k \begin{align*} x_{k1} Ax_k Bu_k \\ y_k Cx_k \end{align*} xk1ykAxkBukCxk 由 【MPC】模型预测控制笔记 (1)无约束MPC 可知 无约束MPC的输入可写为状态反馈形式 u k − K x k u_k -Kx_k uk−Kxk 其中 K [ I p × p 0 0 ⋯ 0 ] ( H T Q ′ H R ′ ) − 1 H T Q ′ G K [I_{p\times p} ~0 ~0 ~\cdots ~0](\mathcal{H}^T \mathcal{Q}^\prime \mathcal{H} \mathcal{R}^\prime)^{-1}\mathcal{H}^T \mathcal{Q}^\prime \mathcal{G} K[Ip×p 0 0 ⋯ 0](HTQ′HR′)−1HTQ′G. 无约束输出反馈MPC可直接使用观测器观测状态代替真实状态 u k − K x ^ k u_k -K\hat{x}_k uk−Kx^k 其中 x ^ k 1 A x ^ k B u k L ( y k − C x ^ k ) \hat{x}_{k1} A\hat{x}_{k} Bu_k L(y_k - C\hat{x}_{k}) x^k1Ax^kBukL(yk−Cx^k)记观测误差为 x ~ k x k − x ^ k \tilde{x}_k x_k - \hat{x}_{k} x~kxk−x^k. 系统可改写为 x k 1 A x k − B K x ^ k A x k − B K ( x k − x ~ k ) (4) \begin{align*} x_{k1} Ax_k - BK\hat{x}_k \\ Ax_k - BK(x_k - \tilde{x}_k) \end{align*} \tag{4} xk1Axk−BKx^kAxk−BK(xk−x~k)(4)
观测误差可写为 x ~ k 1 ( A − L C ) x ~ k (5) \tilde{x}_{k1} (A-LC)\tilde{x}_k \tag{5} x~k1(A−LC)x~k(5)
2.1 稳定性分析
构建增广系统 [ x k 1 x ~ k 1 ] [ A − B K B K 0 A − L C ] [ x k x ~ k ] \begin{align*} \begin{bmatrix} x_{k1} \\ \tilde{x}_{k1} \end{bmatrix} \begin{bmatrix} A-BK BK \\ 0 A - LC \end{bmatrix} \begin{bmatrix} x_{k} \\ \tilde{x}_{k} \end{bmatrix} \end{align*} [xk1x~k1][A−BK0BKA−LC][xkx~k] 系统的的稳定性取决于 A a u g [ A − B K B K 0 A − L C ] \begin{align*} A_{aug} \begin{bmatrix} A-BK BK \\ 0 A - LC \end{bmatrix} \end{align*} Aaug[A−BK0BKA−LC] 的特征值。 当 ∣ e i g ( A a u g ) ∣ 1 |\mathrm{eig}(A_{aug})|1 ∣eig(Aaug)∣1 时系统渐近稳定其中 e i g ( A a u g ) e i g ( A − B K ) ∪ e i g ( A − L C ) \mathrm{eig}(A_{aug}) \mathrm{eig}(A-BK) \cup \mathrm{eig}(A-LC) eig(Aaug)eig(A−BK)∪eig(A−LC). 由此可见整体系统的稳定性取决了控制器和观测器这两部分可独立设计 但需要注意的时观测器应收敛得比控制器更快否则控制器会根据不准确的观测结果将系统控歪。 线性代数基础三角分块矩阵的特征值为其对角子块上的特征值集合。 证明 首先证明三角分块矩阵的行列式等于其对角子块行列式的乘积参考分块矩阵行列式的性质证明 性质1从一行列中减去另一行列的倍数行列式保持不变 性质2三角矩阵的行列式为其主对角线上元素的乘积 参考《线性代数导论第五版 (美Gilbert Strang著 海昕, 文军, 屈龙江, 钱旭译) 》 设 ∣ D ∣ ∣ a 11 ⋯ a 1 k ⋮ ⋮ 0 a k 1 ⋯ a k k c 11 ⋯ c 1 k b 11 ⋯ b 1 n ⋮ ⋮ ⋮ ⋮ c n 1 ⋯ c n k b n 1 ⋯ b n n ∣ , ∣ D 1 ∣ ∣ a 11 ⋯ a 1 k ⋮ ⋮ a k 1 ⋯ a k k ∣ , ∣ D 2 ∣ ∣ b 11 ⋯ b 1 n ⋮ ⋮ b n 1 ⋯ b n n ∣ \begin{align*} |D| \begin{vmatrix} a_{11} \cdots a_{1k} \\ \vdots \vdots 0 \\ a_{k1} \cdots a_{kk} \\ c_{11} \cdots c_{1k} b_{11} \cdots b_{1n} \\ \vdots \vdots \vdots \vdots \\ c_{n1} \cdots c_{nk} b_{n1} \cdots b_{nn} \end{vmatrix}, |D_{1}| \begin{vmatrix} a_{11} \cdots a_{1k} \\ \vdots \vdots \\ a_{k1} \cdots a_{kk} \end{vmatrix}, |D_{2}| \begin{vmatrix} b_{11} \cdots b_{1n} \\ \vdots \vdots \\ b_{n1} \cdots b_{nn} \end{vmatrix} \end{align*} ∣D∣ a11⋮ak1c11⋮cn1⋯⋯⋯⋯a1k⋮akkc1k⋮cnkb11⋮bn10⋯⋯b1n⋮bnn ,∣D1∣ a11⋮ak1⋯⋯a1k⋮akk ,∣D2∣ b11⋮bn1⋯⋯b1n⋮bnn 分别通过行变换和列变换可得 ∣ D 1 ∣ ∣ p 11 ⋯ 0 ⋮ ⋱ ⋮ p k 1 ⋯ p k k ∣ p 11 p 22 ⋯ p k k , ∣ D 2 ∣ ∣ q 11 ⋯ 0 ⋮ ⋱ ⋮ q n 1 ⋯ q n n ∣ p 11 p 22 ⋯ p n n \begin{align*} |D_{1}| \begin{vmatrix} p_{11} \cdots 0 \\ \vdots \ddots \vdots \\ p_{k1} \cdots p_{kk} \end{vmatrix} \mathrm{p}_{11} \mathrm{p}_{22} \cdots \mathrm{p}_{\mathrm{kk}}, |D_{2}| \begin{vmatrix} q_{11} \cdots 0 \\ \vdots \ddots \vdots \\ q_{n1} \cdots q_{nn} \end{vmatrix} \mathrm{p}_{11} \mathrm{p}_{22} \cdots \mathrm{p}_{\mathrm{nn}} \end{align*} ∣D1∣ p11⋮pk1⋯⋱⋯0⋮pkk p11p22⋯pkk,∣D2∣ q11⋮qn1⋯⋱⋯0⋮qnn p11p22⋯pnn 有 ∣ D ∣ ∣ p 11 ⋯ 0 ⋮ ⋮ 0 p k 1 ⋯ p k k c 11 ⋯ c 1 k q 11 ⋯ 0 ⋮ ⋮ ⋮ ⋮ c n 1 ⋯ c n k q n 1 ⋯ q n n ∣ p 11 p 22 ⋯ p k k ⋅ p 11 p 22 ⋯ p n n ∣ D 1 ∣ ∣ D 2 ∣ \begin{align*} |D| \begin{vmatrix} p_{11} \cdots 0 \\ \vdots \vdots 0 \\ p_{k1} \cdots p_{kk} \\ c_{11} \cdots c_{1k} q_{11} \cdots 0 \\ \vdots \vdots \vdots \vdots \\ c_{n1} \cdots c_{nk} q_{n1} \cdots q_{nn} \end{vmatrix} \end{align*} \mathrm{p}_{11} \mathrm{p}_{22} \cdots \mathrm{p}_{\mathrm{kk}} \cdot \mathrm{p}_{11} \mathrm{p}_{22} \cdots \mathrm{p}_{\mathrm{nn}} |D_1||D_2| ∣D∣ p11⋮pk1c11⋮cn1⋯⋯⋯⋯0⋮pkkc1k⋮cnkq11⋮qn10⋯⋯0⋮qnn p11p22⋯pkk⋅p11p22⋯pnn∣D1∣∣D2∣ 在计算特征值时有 ∣ D − λ I ∣ ∣ D 1 − λ I ∣ ∣ D 2 − λ I ∣ 0 |D-\lambda I| |D_1-\lambda I| |D_2-\lambda I| 0 ∣D−λI∣∣D1−λI∣∣D2−λI∣0 故 三角分块矩阵的特征值为其对角子块上的特征值集合 得证. 2.2 MATLAB实例
针对系统 (1)设系统中 A [ 1.1 2 0 0.95 ] A \begin{bmatrix} 1.1 2 \\ 0 0.95 \end{bmatrix} A[1.1020.95] B [ 0 0.079 ] B \begin{bmatrix} 0 \\ 0.079 \end{bmatrix} B[00.079] C [ 1 0 ] C \begin{bmatrix} 1 0 \end{bmatrix} C[10]. 由 【MPC】模型预测控制笔记 (1)无约束MPC 可知控制器反馈增益 K [ I p × p 0 0 ⋯ 0 ] ( H T Q ′ H R ′ ) − 1 H T Q ′ G K [I_{p\times p} ~0 ~0 ~\cdots ~0](\mathcal{H}^T \mathcal{Q}^\prime \mathcal{H} \mathcal{R}^\prime)^{-1}\mathcal{H}^T \mathcal{Q}^\prime \mathcal{G} K[Ip×p 0 0 ⋯ 0](HTQ′HR′)−1HTQ′G 计算可得 K [ 2.6167 12.9286 ] K \begin{bmatrix} 2.6167 \\ 12.9286 \end{bmatrix} K[2.616712.9286]MATLAB代码如下
%% 选取K并检验系统是否稳定
A [1.1 2;0 0.95];
B [0; 0.079];
K [1.4 5.76];eigSys eig(A - B*K);
disp(abs(eigSys))
%% 求解P
K [1.4 5.76];
Q eye(2);
R 0.1;
syms P [2 2] % P 为2*2的矩阵
equ P - (A - B*K) * P * (A - B*K) Q K*R*K;
Psol solve(equ, P);
Psol [Psol.P1_1, Psol.P2_1; Psol.P2_1, Psol.P2_2];
Psol double(Psol);
disp(Psol)
%% 计算G、H
N 4;
[G, H] getGH(N, A, B);
%% 计算Q、R、K
[Qp, Rp] getQR(N, Q, Psol, R);
Kp [eye(1), kron(ones(1, N-1), zeros(1))] * inv(H*Qp*H Rp)*H*Qp*G;由 1.3 中可知观测增益 . 最终得到的系统的动态如下 ( L [ 1.8342 0.3571 ] L \begin{bmatrix} 1.8342 \\ 0.3571 \end{bmatrix} L[1.83420.3571]) ( L [ 1.0500 0.2263 ] L \begin{bmatrix} 1.0500 \\ 0.2263 \end{bmatrix} L[1.05000.2263])
MATLAB演示代码
A [1.1 2;0 0.95];
B [0; 0.079];
C [1, 0];
K [2.6167, 12.9286];
% L [1.05; 0.2263];
L [1.8342; 0.3571];xCur [1.2;-0.7]; % 设初始状态为[1;1]
xLog xCur;
xHatCur [0;0]; % 设初始状态为[0;0]
xHatLog xHatCur;
uLog [];step 0:50;
u 0;
for i stepu -K * xHatCur;yCur C * xCur; % y_kxHatCur A * xHatCur B*u L * (yCur - C * xHatCur); % xHat_k1xCur A*xCur B*u; % x_k1xHatLog [xHatLog, xHatCur];xLog [xLog, xCur];uLog [uLog, u];
endfigure(1)
subplot(3,1,1)
hold on
plot(step, xLog(1,1:end-1))
plot(step, xHatLog(1,1:end-1))
title(x1)
grid on
subplot(3,1,2)
hold on
plot(step, xLog(2,1:end-1))
plot(step, xHatLog(2,1:end-1))
title(x2)
grid on
subplot(3,1,3)
plot(step, uLog)
title(u)
grid on
