网站漏洞 在线扫描,贵州网站建设营销公司,推广的方式,企业查询天眼查免费目录 前言【基础】二次型是凸函数【基础】常见的凸约束 一、约束MPC的求解1.1 等式约束MPC1.2 不等式约束MPC 二、稳定性分析2.1 终端等式约束2.1.1 稳定性证明2.1.2 迭代可行性分析2.1.3 MATLAB应用实例2.1.4 总结 2.2 终端不等式约束#xff08;优化无限时域#xff09;2.2… 目录 前言【基础】二次型是凸函数【基础】常见的凸约束 一、约束MPC的求解1.1 等式约束MPC1.2 不等式约束MPC 二、稳定性分析2.1 终端等式约束2.1.1 稳定性证明2.1.2 迭代可行性分析2.1.3 MATLAB应用实例2.1.4 总结 2.2 终端不等式约束优化无限时域2.2.1 求解步骤与稳定性保证2.2.2 迭代可行性分析2.2.3 MATLAB应用实例2.2.4 总结 附录1附录2 前言
由上一节可以看到MPC是通过优化代价函数来求解最优控制序列这一般是通过凸优化来进行的。 凸优化问题已有高效的算法来进行求解我们确保构建的优化问题是凸优化即可。 致谢 【模型预测控制2022春lecture 2-1 Constrained MPC】、【模型预测控制2022春lecture 2-2 Constrained MPC】
【基础】二次型是凸函数
凸函数定义设 C C C 是非空凸集 f f f 是定义在 D D D 上的函数对任意 x 1 , x 2 ∈ C x_1, x_2 \in C x1,x2∈C λ ∈ ( 0 , 1 ) \lambda \in (0,~1) λ∈(0, 1) 均有 f ( λ x 1 ( 1 − λ ) x 2 ) ≤ λ f ( x 1 ) ( 1 − λ ) f ( x 2 ) (1) f(\lambda x_1 (1-\lambda) x_2) \le \lambda f(x_1)(1-\lambda)f(x_2) \tag{1} f(λx1(1−λ)x2)≤λf(x1)(1−λ)f(x2)(1)则称 f f f 为 C C C 上的凸函数。
证明二次型 f ( x ) 1 2 x T H x f T x f(x) \frac{1}{2}x^THxf^Tx f(x)21xTHxfTx是凸函数 其中 x ∈ R n x \in \mathbb{R}^n x∈Rn
任取 x 1 , x 2 ∈ R n x_1, x_2 \in \mathbb{R}^n x1,x2∈Rn 计算式 (1) 左侧 f ( λ x 1 ( 1 − λ ) x 2 ) 1 2 [ λ 2 x 1 T H x 1 ( 1 − λ ) 2 x 2 T H x 2 2 λ ( 1 − λ ) x 1 T H x 2 ] f T [ λ x 1 ( 1 − λ ) x 2 ] \begin{align*} f(\lambda x_1 (1-\lambda) x_2) \frac{1}{2} \left[ \lambda ^2 x_1^THx_1 (1-\lambda) ^2 x_2^THx_2 2\lambda (1-\lambda) x_1^THx_2\right] f^T \left[ \lambda x_1 (1-\lambda) x_2\right] \end{align*} f(λx1(1−λ)x2)21[λ2x1THx1(1−λ)2x2THx22λ(1−λ)x1THx2]fT[λx1(1−λ)x2] 计算式 (1) 右侧 λ f ( x 1 ) ( 1 − λ ) f ( x 2 ) 1 2 λ x 1 T H x 1 1 2 ( 1 − λ ) x 2 T H x 2 λ f T x 1 ( 1 − λ ) f T x 2 \begin{align*} \lambda f(x_1) (1-\lambda) f(x_2) \frac{1}{2} \lambda x_1^THx_1 \frac{1}{2} (1-\lambda) x_2^THx_2 \lambda f^Tx_1 (1-\lambda) f^Tx_2 \end{align*} λf(x1)(1−λ)f(x2)21λx1THx121(1−λ)x2THx2λfTx1(1−λ)fTx2 左侧减右侧得 f ( λ x 1 ( 1 − λ ) x 2 ) − ( λ f ( x 1 ) ( 1 − λ ) f ( x 2 ) ) − 1 2 λ ( 1 − λ ) x 1 T H x 1 − 1 2 λ ( 1 − λ ) x 2 T H x 2 λ ( 1 − λ ) x 1 T H x 2 1 2 λ ( 1 − λ ) [ x 1 T H x 2 − x 1 T H x 1 x 1 T H x 2 − x 2 T H x 2 ] 1 2 λ ( 1 − λ ) [ x 1 T H ( x 2 − x 1 ) ( x 1 − x 2 ) T H x 2 ] 1 2 λ ( 1 − λ ) [ x 1 T H ( x 2 − x 1 ) − x 2 T H ( x 2 − x 1 ) ] − 1 2 λ ( 1 − λ ) ( x 2 − x 1 ) T H ( x 2 − x 1 ) ≤ 0 (当 H 半正定) \begin{align*} \hspace{0.5cm} f(\lambda x_1 (1-\lambda) x_2) - \left(\lambda f(x_1) (1-\lambda) f(x_2) \right) \\ -\frac{1}{2} \lambda (1 - \lambda )x_1^THx_1 -\frac{1}{2}\lambda(1-\lambda) x_2^THx_2 \lambda (1-\lambda) x_1^THx_2 \\ \frac{1}{2}\lambda(1-\lambda) \left[ x_1^THx_2 - x_1^THx_1 x_1^THx_2 - x_2^THx_2 \right] \\ \frac{1}{2}\lambda(1-\lambda) \left[ x_1^TH(x_2 - x_1) (x_1 - x_2)^THx_2 \right] \\ \frac{1}{2}\lambda(1-\lambda) \left[ x_1^TH(x_2 - x_1) - x_2^TH(x_2 - x_1) \right] \\ -\frac{1}{2}\lambda(1-\lambda) (x_2-x_1)^TH(x_2 - x_1) \le 0 \quad \text{(当 $H$ 半正定)} \end{align*} f(λx1(1−λ)x2)−(λf(x1)(1−λ)f(x2))−21λ(1−λ)x1THx1−21λ(1−λ)x2THx2λ(1−λ)x1THx221λ(1−λ)[x1THx2−x1THx1x1THx2−x2THx2]21λ(1−λ)[x1TH(x2−x1)(x1−x2)THx2]21λ(1−λ)[x1TH(x2−x1)−x2TH(x2−x1)]−21λ(1−λ)(x2−x1)TH(x2−x1)≤0(当 H 半正定) 即 f ( λ x 1 ( 1 − λ ) x 2 ) ≤ λ f ( x 1 ) ( 1 − λ ) f ( x 2 ) f(\lambda x_1 (1-\lambda) x_2) \le \lambda f(x_1) (1-\lambda) f(x_2) f(λx1(1−λ)x2)≤λf(x1)(1−λ)f(x2)
【基础】常见的凸约束
约束条件构成凸集即为凸约束 凸集定义如果 C C C 中任意两点间的线段仍然在 C C C 中则该集合为凸集即对于任意 x 1 , x 2 ∈ C x_1, x_2 \in C x1,x2∈C λ ∈ [ 0 , 1 ] \lambda \in [0,~1] λ∈[0, 1]都有 λ x 1 ( 1 − λ ) x 2 ∈ C \lambda x_1 (1-\lambda) x_2 \in C λx1(1−λ)x2∈C. 来源凸优化(Convex Optimization) (Stephen Boyd, Lieven Vandenberghe, 王书宁译) 常见凸约束 (1) 等式约束 A e q x b e q A_{eq}xb_{eq} Aeqxbeq (2) 不等式约束 A x ≤ b Ax \le b Ax≤b
证明 1任取 x 1 , x 2 ∈ { x ∣ A e q x b e q } x_1, x_2 \in \{x |A_{eq}xb_{eq}\} x1,x2∈{x∣Aeqxbeq} λ ∈ [ 0 , 1 ] \lambda \in [0,~1] λ∈[0, 1]有 A e q [ λ x 1 ( 1 − λ ) x 2 ] λ A e q x 1 ( 1 − λ ) A e q x 2 λ b e q ( 1 − λ ) b e q b e q \begin{align*} \hspace{0.5cm} A_{eq}\left[\lambda x_1 (1-\lambda) x_2 \right] \\ \lambda A_{eq}x_1 (1-\lambda) A_{eq}x_2 \\ \lambda b_{eq} (1-\lambda) b_{eq} \\ b_{eq} \end{align*} Aeq[λx1(1−λ)x2]λAeqx1(1−λ)Aeqx2λbeq(1−λ)beqbeq故 λ x 1 ( 1 − λ ) x 2 ∈ { x ∣ A e q x b e q } \lambda x_1 (1-\lambda) x_2 \in \{x |A_{eq}xb_{eq}\} λx1(1−λ)x2∈{x∣Aeqxbeq}.
2任取 x 1 , x 2 ∈ { x ∣ A x ≤ b } x_1, x_2 \in \{x |Ax \le b\} x1,x2∈{x∣Ax≤b} λ ∈ [ 0 , 1 ] \lambda \in [0,~1] λ∈[0, 1]有 A [ λ x 1 ( 1 − λ ) x 2 ] ≤ λ A x 1 ( 1 − λ ) A x 2 λ b ( 1 − λ ) b b \begin{align*} \hspace{0.5cm} A\left[\lambda x_1 (1-\lambda) x_2 \right] \\ \le \lambda Ax_1 (1-\lambda) Ax_2 \\ \lambda b (1-\lambda) b \\ b \end{align*} A[λx1(1−λ)x2]≤λAx1(1−λ)Ax2λb(1−λ)bb故 λ x 1 ( 1 − λ ) x 2 ∈ { x ∣ A x b } \lambda x_1 (1-\lambda) x_2 \in \{x |Axb\} λx1(1−λ)x2∈{x∣Axb}. 一、约束MPC的求解
约束MPC的求解相似于上一节【MPC】模型预测控制笔记 (1)无约束MPC同样是构建二次规划问题来求解最优控制序列 U ∗ U^* U∗即 U ∗ a r g min U 1 2 U T H U f T U U^* \mathrm{arg} \min_U \frac{1}{2}U^THUf^TU U∗argUmin21UTHUfTU 符号说明 arg min (argument of the minimum): 使得目标函数取得最小值的自变量. s. t. (subject to): 满足以下条件 1.1 等式约束MPC
在等式约束下优化问题可表述为 U ∗ a r g min U 1 2 U T H U f T U s . t . A e q U b e q \begin{align*} U^* \mathrm{arg} \min_U \frac{1}{2}U^THUf^TU \\ \mathrm{s. t.} \quad A_{eq}U b_{eq} \end{align*} U∗argUmin21UTHUfTUs.t.AeqUbeq 可通过拉格朗日乘子法将带约束优化问题转换为无约束优化问题。 引入未知的拉格朗日乘子 λ \lambda λ构造拉格朗日函数为 L 1 2 U T H U f T U λ ( A e q U − b e q ) \mathcal{L} \frac{1}{2}U^THUf^TU \lambda (A_{eq}U - b_{eq}) L21UTHUfTUλ(AeqU−beq) 通过令一阶导数为 0 求解 U U U 和 λ \lambda λ: { ∂ L ∂ U H U f λ A e q T 0 ∂ L ∂ λ A e q U − b e q 0 \left\{ \begin{aligned} \frac{\partial \mathcal{L}}{\partial U} HU f \lambda A_{eq}^T 0 \\ \frac{\partial \mathcal{L}}{\partial \lambda} A_{eq}U - b_{eq} 0 \end{aligned} \right. ⎩ ⎨ ⎧∂U∂L∂λ∂LHUfλAeqT0AeqU−beq0
1.2 不等式约束MPC
在不等式约束下优化问题可表述为 U ∗ a r g min U 1 2 U T H U f T U s . t . A U ≤ b \begin{align*} U^* \mathrm{arg} \min_U \frac{1}{2}U^THUf^TU \\ \mathrm{s. t.} \quad AU \le b \end{align*} U∗argUmin21UTHUfTUs.t.AU≤b 在 active-set 求解方法中首先忽略约束求解 U u c ∗ U^*_{uc} Uuc∗若 U u c ∗ U^*_{uc} Uuc∗ 超出了不等式的约束范围说明 U ∗ U^* U∗ 落在了不等式约束的某一边界上即满足 A ′ U b ′ A^\prime Ub^\prime A′Ub′此时可根据被激活的约束构建等式约束优化问题来求解。
二、稳定性分析
最优不能保证系统稳定故需要一些措施来保证系统是渐近稳定的。 以下所有内容均针对线性定常系统 x k 1 A x x B u k x_{k1} Ax_{x} Bu_k xk1AxxBuk
2.1 终端等式约束
终端等式约束即强制令MPC最后一个状态优化为0 即在原本的优化问题中增加约束 x ( N ∣ k ) 0 x_{(N|k)} 0 x(N∣k)0.
2.1.1 稳定性证明
选取最优的代价函数作为李雅普诺夫函数 V ( x k ) V(x_k) V(xk) V ( x k ) J k ∗ ∑ i 1 N ( x ( i ∣ k ) T Q x ( i ∣ k ) u ( i ∣ k ) T R u ( i − 1 ∣ k ) ) V(x_k) J^*_k\sum^N_{i1} \left(x^T_{(i|k)}Qx_{(i|k)} u^T_{(i|k)}Ru_{(i-1|k)} \right) V(xk)Jk∗i1∑N(x(i∣k)TQx(i∣k)u(i∣k)TRu(i−1∣k)) 显然满足 V ( x k ) 0 V(x_k) 0 V(xk)0. 设 J k 1 J_{k1} Jk1 为完全按 U k ∗ U_k^* Uk∗ 序列执行下的代价且令 u ( N − 1 ∣ k 1 ) 0 u_{(N-1|k1)} 0 u(N−1∣k1)0 通过终端约束 x ( N ∣ k ) 0 x_{(N|k)} 0 x(N∣k)0 可知在 U k ∗ U_k^* Uk∗ 序列下有 x ( N ∣ k 1 ) A x ( N − 1 ∣ k 1 ) B u ( N − 1 ∣ k 1 ) A x ( N ∣ k ) 0 0 x_{(N|k1)} Ax_{(N-1|k1)} Bu_{(N-1|k1)} Ax_{(N|k)} 0 0 x(N∣k1)Ax(N−1∣k1)Bu(N−1∣k1)Ax(N∣k)00. 故有 J k 1 ∗ − J k ∗ ≤ J k 1 − J k ∗ ∑ i 1 N ( x ( i ∣ k 1 ) T Q x ( i ∣ k 1 ) u ( i − 1 ∣ k 1 ) T R u ( i − 1 ∣ k 1 ) ) − ∑ i 1 N ( x ( i ∣ k ) T Q x ( i ∣ k ) u ( i − 1 ∣ k ) T R u ( i − 1 ∣ k ) ) [ ∑ i 1 N − 1 ( x ( i 1 ∣ k ) T Q x ( i 1 ∣ k ) u ( i ∣ k ) T R u ( i − 1 ∣ k ) ) 0 0 ] − ∑ i 1 N ( x ( i ∣ k ) T Q x ( i ∣ k ) u ( i − 1 ∣ k ) T R u ( i − 1 ∣ k ) ) − ( x ( 1 ∣ k ) T Q x ( 1 ∣ k ) u ( 0 ∣ k ) T R u ( 0 ∣ k ) ) ≤ 0 \begin{align*} J_{k1}^* - J_{k}^* \le J_{k1} - J_{k}^* \\ \sum^{N}_{i1} \left( x^T_{(i|k1)}Qx_{(i|k1)} u^T_{(i-1|k1)}Ru_{(i-1|k1)} \right) - \sum^{N}_{i1} \left( x^T_{(i|k)}Qx_{(i|k)} u^T_{(i-1|k)}Ru_{(i-1|k)} \right) \\ \left[ \sum^{N-1}_{i1} \left( x^T_{(i1|k)}Qx_{(i1|k)} u^T_{(i|k)}Ru_{(i-1|k)} \right) 0 0 \right] - \sum^{N}_{i1} \left( x^T_{(i|k)}Qx_{(i|k)} u^T_{(i-1|k)}Ru_{(i-1|k)} \right) \\ - \left( x^T_{(1|k)}Qx_{(1|k)} u^T_{(0|k)}Ru_{(0|k)} \right) \\ \le 0 \end{align*} Jk1∗−Jk∗≤Jk1−Jk∗i1∑N(x(i∣k1)TQx(i∣k1)u(i−1∣k1)TRu(i−1∣k1))−i1∑N(x(i∣k)TQx(i∣k)u(i−1∣k)TRu(i−1∣k))[i1∑N−1(x(i1∣k)TQx(i1∣k)u(i∣k)TRu(i−1∣k))00]−i1∑N(x(i∣k)TQx(i∣k)u(i−1∣k)TRu(i−1∣k))−(x(1∣k)TQx(1∣k)u(0∣k)TRu(0∣k))≤0 其中 Q Q Q、 R R R 正定。故有 Δ V ( x ) J k 1 ∗ − J k ∗ ≤ 0 \Delta V(x) J_{k1}^* - J_{k}^* \le 0 ΔV(x)Jk1∗−Jk∗≤0当且仅当 x ( 1 ∣ k ) T 0 x^T_{(1|k)}0 x(1∣k)T0 u ( 0 ∣ k ) T 0 u^T_{(0|k)}0 u(0∣k)T0时等号成立。 故系统渐近稳定。
2.1.2 迭代可行性分析
终端约束 x ( N ∣ k ) 0 x_{(N|k)} 0 x(N∣k)0 其实是难以保证可达到的 可以证明的是当在初始状态下可达则未来迭代中均可达 所谓迭代可行性。
在无约束系统中 一阶系统可以在 1 步的跌代中使状态归零 如系统 v ˙ F / m \dot{v} F/m v˙F/m无论初值 v 0 v_0 v0 为多少都可施加适当的作用力使速度在一个迭代中归零即令 v v 0 ( F ∗ / m ) Δ t 0 v v_0 (F^*/m) \Delta t 0 vv0(F∗/m)Δt0. 二阶系统可以在 2 步的跌代中使状态归零 如系统 { s ˙ v v ˙ F / m \left\{ \begin{aligned} \dot{s} v\\ \dot{v} F/m \end{aligned} \right. {s˙v˙vF/m 在第一步迭代中 s 1 s 0 v 0 Δ t s_1 s_0 v_0 \Delta t s1s0v0Δt使 v 1 s 1 / Δ t v_1 s_1/\Delta t v1s1/Δt则 F 0 m ( v 1 − v 0 ) / Δ t F_0 m(v_1-v_0)/ \Delta t F0m(v1−v0)/Δt. 在第二步迭代中 s 2 s 1 v 1 Δ t 0 s_2 s_1 v_1 \Delta t 0 s2s1v1Δt0使 v 2 v 1 ( F 1 / m ) Δ t 0 v_2 v_1 (F_1/m)\Delta t 0 v2v1(F1/m)Δt0 即完成系统状态归零. 但在约束系统中这是难以保证的。
我们可分析初始状态下终端约束 x ( N ∣ k ) 0 x_{(N|k)} 0 x(N∣k)0 是否可达来保证终端约束可达。
迭代可行性终端约束 x ( N ∣ k ) 0 x_{(N|k)} 0 x(N∣k)0 当在初始状态下可达则未来迭代中均可达
证明 若在初始状态下有可行解 U 0 [ u ( 0 ∣ 0 ) , u ( 1 ∣ 0 ) , u ( 2 ∣ 0 ) , ⋯ , u ( N − 1 ∣ 0 ) ] U_0 \left[ u_{(0|0)},~ u_{(1|0)}, ~u_{(2|0)},~\cdots,~ u_{(N-1|0)} \right] U0[u(0∣0), u(1∣0), u(2∣0), ⋯, u(N−1∣0)]使得 x ( N ∣ 0 ) 0 x_{(N|0)} 0 x(N∣0)0. 则在执行 u ( 0 ∣ 0 ) u_{(0|0)} u(0∣0) 后进入下一状态必然存在可行解 U 1 [ u ( 1 ∣ 0 ) , u ( 2 ∣ 0 ) , ⋯ , u ( N − 1 ∣ 0 ) , u ( N − 1 ∣ 1 ) ] U_1 \left[ u_{(1|0)}, ~u_{(2|0)},~\cdots,~ u_{(N-1|0)},~u_{(N-1|1)} \right] U1[u(1∣0), u(2∣0), ⋯, u(N−1∣0), u(N−1∣1)]. 其中 u ( N − 1 ∣ 1 ) 0 u_{(N-1|1)} 0 u(N−1∣1)0. 因 x ( N − 1 ∣ 1 ) x ( N ∣ 0 ) 0 x_{(N-1|1)} x_{(N|0)} 0 x(N−1∣1)x(N∣0)0仅当 u ( N − 1 ∣ 1 ) 0 u_{(N-1|1)} 0 u(N−1∣1)0 时有 x ( N ∣ 1 ) A x ( N − 1 ∣ 1 ) B u ( N − 1 ∣ 1 ) 0 x_{(N|1)} Ax_{(N-1|1)} Bu_{(N-1|1)} 0 x(N∣1)Ax(N−1∣1)Bu(N−1∣1)0. 由此类推 当在初始状态下有可行解 U 0 [ u ( 0 ∣ 0 ) , u ( 1 ∣ 0 ) , u ( 2 ∣ 0 ) , ⋯ , u ( N − 1 ∣ 0 ) ] U_0 \left[ u_{(0|0)},~ u_{(1|0)}, ~u_{(2|0)},~\cdots,~ u_{(N-1|0)} \right] U0[u(0∣0), u(1∣0), u(2∣0), ⋯, u(N−1∣0)] 时 系统必然有可行解 [ u ( 0 ∣ 0 ) , u ( 1 ∣ 0 ) , u ( 2 ∣ 0 ) , ⋯ , u ( N − 1 ∣ 0 ) , 0 , 0 , ⋯ ] \left[ u_{(0|0)},~ u_{(1|0)}, ~u_{(2|0)},~\cdots,~ u_{(N-1|0)},~0,~ 0,~ \cdots \right] [u(0∣0), u(1∣0), u(2∣0), ⋯, u(N−1∣0), 0, 0, ⋯] 可满足终端约束.
2.1.3 MATLAB应用实例
针对系统 x k 1 A x k B u k x_{k1} Ax_k Bu_k xk1AxkBuk 其中 x ∈ R 2 x \in \mathbb{R}^2 x∈R2 u ∈ R 1 u \in \mathbb{R}^1 u∈R1 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]且输入需要满足约束 − 2 ≤ u ≤ 2 -2 \le u \le 2 −2≤u≤2.
1 N N N 步预测空间的状态可写为 X G x ( 0 ∣ k ) H U (1) X\mathcal{G}x_{(0|k)} \mathcal{H}U \tag{1} XGx(0∣k)HU(1) 其中 X [ x ( 1 ∣ k ) x ( 2 ∣ k ) ⋯ x ( N ∣ k ) ] T U [ u ( 0 ∣ k ) u ( 1 ∣ k ) ⋯ u ( N − 1 ∣ k ) ] T G [ A A 2 ⋯ A N ] T H [ B 0 0 ⋯ 0 A B B 0 ⋯ 0 A 2 B A B B ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ A N − 1 B A 2 B A B ⋯ B ] \begin{align*} X [x_{(1|k)} ~ x_{(2|k)} ~ \cdots~x_{(N|k)}]^T \\ U [u_{(0|k)} ~ u_{(1|k)} ~ \cdots~u_{(N-1|k)}]^T \\ \mathcal{G} \left[ A ~ A^2 ~\cdots ~ A^N \right]^T \\ \mathcal{H} \begin{bmatrix} B 0 0 \cdots 0\\ AB B 0 \cdots 0\\ A^2B AB B \cdots 0 \\ \vdots \vdots \vdots \ddots \vdots \\ A^{N-1}B A^2B AB \cdots B \end{bmatrix} \end{align*} XUGH[x(1∣k) x(2∣k) ⋯ x(N∣k)]T[u(0∣k) u(1∣k) ⋯ u(N−1∣k)]T[A A2 ⋯ AN]T BABA2B⋮AN−1B0BAB⋮A2B00B⋮AB⋯⋯⋯⋱⋯000⋮B 2定义代价函数为 J X T Q X U T R U (2) JX^T\mathcal{Q}X U^T\mathcal{R}U \tag{2} JXTQXUTRU(2) 3将式 (1) 代入式 (2)并考虑约束条件增加终端约束可构建二次规划问题 U k ∗ a r g min U k [ ( G x ( 0 ∣ k ) ) T Q ′ G x ( 0 ∣ k ) 2 x ( 0 ∣ k ) T G T Q ′ H U k U k T ( H T Q ′ H R ) U k ] s . t . − 2 N × 1 ≤ U k ≤ 2 N × 1 x ( N ∣ k ) 0 \begin{align*} U_{k}^* \mathrm{arg} \min_{U_k} \left[ (\mathcal{G}x_{(0|k)})^T \mathcal{Q}^\prime \mathcal{G}x_{(0|k)} 2x_{(0|k)}^T\mathcal{G}^T \mathcal{Q}^\prime \mathcal{H} U_k U_k^T (\mathcal{H}^T \mathcal{Q}^\prime \mathcal{H} \mathcal{R})U_k \right] \\ \hspace{3cm} \mathrm{s. t.} \quad \mathbf{-2}_{N \times1} \le U_k \le \mathbf{2}_{N \times1} \\ \hspace{4cm} x_{(N|k)} 0 \end{align*} Uk∗argUkmin[(Gx(0∣k))TQ′Gx(0∣k)2x(0∣k)TGTQ′HUkUkT(HTQ′HR)Uk]s.t.−2N×1≤Uk≤2N×1x(N∣k)0 其中 x ( N ∣ k ) [ 0 , 0 , ⋯ , I 2 × 2 ] X k [ 0 , 0 , ⋯ , I 2 × 2 ] [ G x ( 0 ∣ k ) H U k ] 0 x_{(N|k)} [0,~0,~\cdots, ~I_{2 \times 2}]X_k [0,~0,~\cdots, ~I_{2 \times 2}] \left[ \mathcal{G}x_{(0|k)} \mathcal{H}U_k \right] 0 x(N∣k)[0, 0, ⋯, I2×2]Xk[0, 0, ⋯, I2×2][Gx(0∣k)HUk]0. 可以改写为 A e q U k b e q A_{eq}U_k b_{eq} AeqUkbeq其中 A e q [ 0 , 0 , ⋯ , I 2 × 2 ] H A_{eq} [0,~0,~\cdots, ~I_{2 \times 2}]\mathcal{H} Aeq[0, 0, ⋯, I2×2]H b e q − [ 0 , 0 , ⋯ , I 2 × 2 ] G x ( 0 ∣ k ) b_{eq} -[0,~0,~\cdots, ~I_{2 \times 2}]\mathcal{G}x_{(0|k)} beq−[0, 0, ⋯, I2×2]Gx(0∣k). 忽略与输入无关的项优化问题可重新写为 U k ∗ a r g min U k [ 1 2 U k T H U k f T U k ] s . t . − 2 N × 1 ≤ U k ≤ 2 N × 1 A e q U k b e q \begin{align*} U_{k}^* \mathrm{arg} \min_{U_k} \left[ \frac{1}{2}U_k^THU_k f^TU_k \right] \\ \mathrm{s. t.} \quad \mathbf{-2}_{N \times1} \le U_k \le \mathbf{2}_{N \times1} \\ \hspace{1cm} A_{eq}U_k b_{eq} \end{align*} Uk∗argUkmin[21UkTHUkfTUk]s.t.−2N×1≤Uk≤2N×1AeqUkbeq其中 H 2 ( H T Q ′ H R ) H 2(\mathcal{H}^T \mathcal{Q}^\prime \mathcal{H} \mathcal{R}) H2(HTQ′HR) f T 2 x ( 0 ∣ k ) T G T Q ′ H f^T 2x_{(0|k)}^T\mathcal{G}^T \mathcal{Q}^\prime \mathcal{H} fT2x(0∣k)TGTQ′H. 以上形式可使用MATLAB的 quadprog 函数求解MATLAB代码见附录1. 设置控制时域 N 20 N20 N20得到结果如下
2.1.4 总结
终端等式约束MPC特点 简单 可保证系统稳定性但需确保初始可行- 需要足够长的控制时域 N N N 来保证初始可行- 增加终端等式约束会削弱求解可行性原来的优化问题是凸优化加上等式约束后可能使原本的凸约束集变为空集无法求解
2.2 终端不等式约束优化无限时域
与上一节 【MPC】模型预测控制笔记 (1)无约束MPC 一致优化无限步预测空间时可通过选取李雅普诺夫直接法证明系统的稳定性。 但找到一个反馈增益 K K K 使系统 x k 1 ( A − B K ) x k x_{k1} (A-BK)x_{k} xk1(A−BK)xk 渐近稳定时显然不能在任意状态下使所有 u k − K x k u_k -Kx_k uk−Kxk 满足约束。
设集合 U \mathcal{U} U 为满足约束条件的输入集合。 我们可以寻找一个不变集 Ω \Omega Ω : 若 x k ∈ Ω x_k \in \Omega xk∈Ω则 x k 1 ∈ Ω x_{k1} \in \Omega xk1∈Ω. 当找到一个反馈增益 K K K 可使系统 x k 1 ( A − B K ) x k x_{k1} (A-BK)x_{k} xk1(A−BK)xk 渐近稳定 且满足 x k ∈ Ω x_k \in \Omega xk∈Ω u k ∈ − K Ω ⊂ U u_k \in -K\Omega \subset \mathcal{U} uk∈−KΩ⊂U 则此时系统未来时刻始终满足约束条件可使用上一节中的无约束MPC的方法求解。
2.2.1 求解步骤与稳定性保证
1增加终端不等式约束 A x ( N , k ) ≤ b Ax_{(N,k)} \le b Ax(N,k)≤b 保证优化的终端状态进入不变集 x ( N , 0 ) ∈ Ω x_{(N,0)} \in \Omega x(N,0)∈Ω 2选择一个反馈增益 K K K 使系统 x k 1 ( A − B K ) x k x_{k1} (A-BK)x_{k} xk1(A−BK)xk 渐近稳定且满足 u k ∈ − K Ω ⊂ U u_k \in -K\Omega \subset \mathcal{U} uk∈−KΩ⊂U 3通过优化无限时域来保证稳定性证明参考【MPC】模型预测控制笔记 (1)无约束MPC 终端不等式约束保证了在 N N N 到 ∞ \infty ∞ 的时域中可使 u k − K x k u_k -Kx_k uk−Kxk满足 u k ∈ U u_k \in \mathcal{U} uk∈U.
2.2.2 迭代可行性分析
当初始状态可满足终端不等式约束 A x ( N , k ) ≤ b Ax_{(N,k)} \le b Ax(N,k)≤b 未来所有状态都可以满足
当在初始状态下有可行解 U 0 [ u ( 0 ∣ 0 ) , u ( 1 ∣ 0 ) , u ( 2 ∣ 0 ) , ⋯ , u ( N − 1 ∣ 0 ) ] U_0 \left[ u_{(0|0)},~ u_{(1|0)}, ~u_{(2|0)},~\cdots,~ u_{(N-1|0)} \right] U0[u(0∣0), u(1∣0), u(2∣0), ⋯, u(N−1∣0)] 时 系统必然有可行解 [ u ( 0 ∣ 0 ) , u ( 1 ∣ 0 ) , u ( 2 ∣ 0 ) , ⋯ , u ( N − 1 ∣ 0 ) , − K x N 1 , − K x N 2 , ⋯ ] \left[ u_{(0|0)},~ u_{(1|0)}, ~u_{(2|0)},~\cdots,~ u_{(N-1|0)},~-Kx_{N1},~ -Kx_{N2},~ \cdots \right] [u(0∣0), u(1∣0), u(2∣0), ⋯, u(N−1∣0), −KxN1, −KxN2, ⋯]使系统在 N N N 步状态后始终满足终端不等式约束.
2.2.3 MATLAB应用实例
同样针对系统 x k 1 A x k B u k x_{k1} Ax_k Bu_k xk1AxkBuk 其中 x ∈ R 2 x \in \mathbb{R}^2 x∈R2 u ∈ R 1 u \in \mathbb{R}^1 u∈R1 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]且输入需要满足约束 − 2 ≤ u ≤ 2 -2 \le u \le 2 −2≤u≤2.
1 N N N 步预测空间的状态可写为 X G x ( 0 ∣ k ) H U X\mathcal{G}x_{(0|k)} \mathcal{H}U XGx(0∣k)HU 其中 X [ x ( 1 ∣ k ) x ( 2 ∣ k ) ⋯ x ( N ∣ k ) ] T U [ u ( 0 ∣ k ) u ( 1 ∣ k ) ⋯ u ( N − 1 ∣ k ) ] T G [ A A 2 ⋯ A N ] T H [ B 0 0 ⋯ 0 A B B 0 ⋯ 0 A 2 B A B B ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ A N − 1 B A 2 B A B ⋯ B ] \begin{align*} X [x_{(1|k)} ~ x_{(2|k)} ~ \cdots~x_{(N|k)}]^T \\ U [u_{(0|k)} ~ u_{(1|k)} ~ \cdots~u_{(N-1|k)}]^T \\ \mathcal{G} \left[ A ~ A^2 ~\cdots ~ A^N \right]^T \\ \mathcal{H} \begin{bmatrix} B 0 0 \cdots 0\\ AB B 0 \cdots 0\\ A^2B AB B \cdots 0 \\ \vdots \vdots \vdots \ddots \vdots \\ A^{N-1}B A^2B AB \cdots B \end{bmatrix} \end{align*} XUGH[x(1∣k) x(2∣k) ⋯ x(N∣k)]T[u(0∣k) u(1∣k) ⋯ u(N−1∣k)]T[A A2 ⋯ AN]T BABA2B⋮AN−1B0BAB⋮A2B00B⋮AB⋯⋯⋯⋱⋯000⋮B 2增加终端约束 x ( N ∣ k ) ∈ Ω x_{(N|k)} \in \Omega x(N∣k)∈Ω其中 x ( N ∣ k ) [ 0 , 0 , ⋯ , I 2 × 2 ] X k [ 0 , 0 , ⋯ , I 2 × 2 ] [ G x ( 0 ∣ k ) H U k ] x_{(N|k)} [0,~0,~\cdots, ~I_{2 \times 2}]X_k [0,~0,~\cdots, ~I_{2 \times 2}] \left[ \mathcal{G}x_{(0|k)} \mathcal{H}U_k \right] x(N∣k)[0, 0, ⋯, I2×2]Xk[0, 0, ⋯, I2×2][Gx(0∣k)HUk] Ω \Omega Ω 为不变集. 不变集的选取待补充目前笔者也不清楚 可以选取足够小的集合 X \mathcal{X} X 来替代不变集但集合范围越大可行性越强。 此处选取 X [ − 0.2 , 0.2 ] × [ − 0.1 , 0.1 ] ( 笛卡尔积形式等价于 { ( x 1 , x 2 ) ∈ R 2 ∣ x 1 ∈ [ − 0.2 , 0.2 ] , x 2 ∈ [ − 0.1 , 0.1 ] } ) \mathcal{X} [-0.2,~0.2] \times [-0.1,~0.1] ~~(\text{笛卡尔积形式等价于}\{(x_1, x_2) \in \mathbb{R}^2 |~ x_1 \in [-0.2,~0.2], x_2 \in [-0.1,~0.1]\}) X[−0.2, 0.2]×[−0.1, 0.1] (笛卡尔积形式等价于{(x1,x2)∈R2∣ x1∈[−0.2, 0.2],x2∈[−0.1, 0.1]}) 约束 x ( N ∣ k ) ∈ X x_{(N|k)} \in \mathcal{X} x(N∣k)∈X 可以改写为 A i n U k ≤ b i n A_{in}U_k \le b_{in} AinUk≤bin其中 A i n [ 1 0 − 1 0 0 1 0 − 1 ] [ 0 , 0 , ⋯ , I 2 × 2 ] H b i n [ 0.2 0.2 0.1 0.1 ] − [ 1 0 − 1 0 0 1 0 − 1 ] [ 0 , 0 , ⋯ , I 2 × 2 ] G x ( 0 ∣ k ) \begin{align*} A_{in} \begin{bmatrix} 1 0 \\ -1 0 \\ 0 1 \\ 0 -1 \end{bmatrix} [0,~0,~\cdots, ~I_{2 \times 2}] \mathcal{H} \\ b_{in} \begin{bmatrix} 0.2 \\ 0.2 \\ 0.1 \\ 0.1 \end{bmatrix}-\begin{bmatrix} 1 0 \\ -1 0 \\ 0 1 \\ 0 -1 \end{bmatrix}[0,~0,~\cdots, ~I_{2 \times 2}] \mathcal{G} x_{(0|k)} \end{align*} Ainbin 1−100001−1 [0, 0, ⋯, I2×2]H 0.20.20.10.1 − 1−100001−1 [0, 0, ⋯, I2×2]Gx(0∣k). 3选取反馈增益 K [ 1.4 5.76 ] K[1.4 \quad 5.76] K[1.45.76]可验证其满足 ∣ e i g ( A − B K ) ∣ 1 |\mathrm{eig}(A-BK)|1 ∣eig(A−BK)∣1 − 2 ≤ − K X ≤ 2 -2 \le -K\mathcal{X} \le 2 −2≤−KX≤2 4通过 P − ( A − B K ) T P ( A − B K ) Q K T R K P-(A-BK)^TP(A-BK) Q K^TRK P−(A−BK)TP(A−BK)QKTRK 求解 P 5定义代价函数为 J X T Q p X U T R p U JX^T\mathcal{Q_p}X U^T\mathcal{R_p}U JXTQpXUTRpU其中 Q p d i a g ( Q , Q , ⋯ , P ) Q_p \mathrm{diag}(Q, Q, \cdots, P) Qpdiag(Q,Q,⋯,P) R p d i a g ( R , R , ⋯ , R ) R_p \mathrm{diag}(R, R, \cdots, R) Rpdiag(R,R,⋯,R). 6构建二次规划问题 U k ∗ a r g min U k [ 1 2 U k T H U k f T U k ] s . t . − 2 N × 1 ≤ U k ≤ 2 N × 1 A i n U k ≤ b i n \begin{align*} U_{k}^* \mathrm{arg} \min_{U_k} \left[ \frac{1}{2}U_k^THU_k f^TU_k \right] \\ \mathrm{s. t.} \quad \mathbf{-2}_{N \times1} \le U_k \le \mathbf{2}_{N \times1} \\ \hspace{1cm} A_{in}U_k \le b_{in} \end{align*} Uk∗argUkmin[21UkTHUkfTUk]s.t.−2N×1≤Uk≤2N×1AinUk≤bin其中 H 2 ( H T Q ′ H R ) H 2(\mathcal{H}^T \mathcal{Q}^\prime \mathcal{H} \mathcal{R}) H2(HTQ′HR) f T 2 x ( 0 ∣ k ) T G T Q ′ H f^T 2x_{(0|k)}^T\mathcal{G}^T \mathcal{Q}^\prime \mathcal{H} fT2x(0∣k)TGTQ′H. MATLAB代码见附录2. 设置控制时域 N 20 N20 N20得到结果如下
2.2.4 总结
终端不等式约束相比终端等式约束 可行性更强但仍需确保初始可行- 更复杂 附录1
%%
A [1.1 2;0 0.95];
B [0; 0.079];
Q eye(2);
R 0.1;% 预测空间
N 20; % 当控制时域过小时如 N10将无法满足终端约束条件求解会导致错误[G, H] getGH(N, A, B);
[Qp, Rp] getQR(N, Q, Q, R);xCur [1;1]; % 设初始状态为[1;1]
%% 约束条件
lb -2 * ones(N, 1);
ub 2 * ones(N, 1);n size(A, 2);
tmp kron(ones(1, N-1), zeros(n));
tmp [tmp, eye(n)];
% Aeq tmp * H;
% beq -tmp * G * xCur;
%% 效果演示
xCur [1; 1]; % 设初始状态为[1;1]
xLog xCur;
uLog [];options optimoptions(quadprog, MaxIterations, 200, Display,none);step 0:50;
for i stepHp 2 * (H * Qp * H Rp);fp 2 * xCur * G * Qp * H;Hp 0.5 * (Hp Hp);Aeq tmp * H;beq -tmp * G * xCur;U quadprog(Hp, fp, [], [], Aeq, beq, lb, ub, zeros(N,1), options);u U(1);xCur A*xCur B*u;xLog [xLog, xCur];uLog [uLog, u];
endfigure(1)
subplot(3,1,1)
plot(step, xLog(1,1:end-1))
title(x1)
grid on
subplot(3,1,2)
plot(step, xLog(2,1:end-1))
title(x2)
grid on
subplot(3,1,3)
plot(step, uLog)
title(u)
grid on
%%
function [Qp, Rp] getQR(N, Q, P, R)Qp eye(N);Qp(end) 0;Qp kron(Qp, Q) kron(eye(N)-Qp, P);Rp eye(N);Rp kron(Rp, R);
endfunction [G, H] getGH(N, A, B) % N1tmp A;G tmp;for i2:Ntmp A*tmp;G [G; tmp];endr size(B, 1);c size(B, 2);H zeros(r * N, c * N);tmp B;for j N:-1:1H( (j-1)*r1:j*r, (j-1)*c1:j*c ) tmp;endfor i 2:Ntmp A*tmp;for j i:NH( (j-1)*r1:j*r, (j-i)*c1:(j-i1)*c ) tmp;endend
end附录2
%%
A [1.1 2;0 0.95];
B [0; 0.079];
Q eye(2);
R 0.1;% 预测空间
N 20;% 求解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] getGH(N, A, B);
[Qp, Rp] getQR(N, Q, Psol, R);xCur [1;1]; % 设初始状态为[1;1]
%% 约束条件
lb -2 * ones(N, 1);
ub 2 * ones(N, 1);n size(A, 2);
tmpReshape kron(ones(1, N-1), zeros(n));
tmpReshape [tmpReshape, eye(n)];
tmpAin [1 0;-1 0;0 1;0 -1];
tmpbin [0.2; 0.2; 0.1; 0.1];
% Ain tmpAin * tmpReshape * H;
% bin tmpbin - tmpAin * tmpReshape * G * xCur;
%% 效果演示
xCur [1; 1]; % 设初始状态为[1;1]
xLog xCur;
uLog [];options optimoptions(quadprog, MaxIterations, 200, Display,none);step 0:50;
u zeros(N, 1);
for i stepHp 2 * (H * Qp * H Rp);fp 2 * xCur * G * Qp * H;Hp 0.5 * (Hp Hp);Ain tmpAin * tmpReshape * H;bin tmpbin - tmpAin * tmpReshape * G * xCur;U quadprog(Hp, fp, Ain, bin, [], [], lb, ub, u, options);u U(1);xCur A*xCur B*u;xLog [xLog, xCur];uLog [uLog, u];
endfigure(1)
subplot(3,1,1)
plot(step, xLog(1,1:end-1))
title(x1)
grid on
subplot(3,1,2)
plot(step, xLog(2,1:end-1))
title(x2)
grid on
subplot(3,1,3)
plot(step, uLog)
title(u)
grid on
%%
function [Qp, Rp] getQR(N, Q, P, R)Qp eye(N);Qp(end) 0;Qp kron(Qp, Q) kron(eye(N)-Qp, P);Rp eye(N);Rp kron(Rp, R);
endfunction [G, H] getGH(N, A, B) % N1tmp A;G tmp;for i2:Ntmp A*tmp;G [G; tmp];endr size(B, 1);c size(B, 2);H zeros(r * N, c * N);tmp B;for j N:-1:1H( (j-1)*r1:j*r, (j-1)*c1:j*c ) tmp;endfor i 2:Ntmp A*tmp;for j i:NH( (j-1)*r1:j*r, (j-i)*c1:(j-i1)*c ) tmp;endend
end
