模型预测控制(MPC)的公式推导与理解（转）

转自：https://blog.csdn.net/qq_40241332/article/details/105023605

（公式：预测域P和控制域M大小不同，且有常干扰）

在本文中，主要是针对线性无约束系统，设计模型预测控制算法。首先给出一个离散的数学模型，再根据模型预测控制“三步走”战略，实现控制器的设计

（相比原文，修改了一个小小的错）

“三步走”策略

预测系统未来动态
求解优化问题
解第一个元素作用于系统

模型：

我们引入离散时间的状态空间模型，如下：

其中 x(k) 为系统内部状态变量；A 为系统矩阵；Bu 为控制输入矩阵；u(k) 为控制输入变量；Bd 为外部干扰输入矩阵；d(k) 为可测外部干扰变量。

预测方程：

将模型改成增量模型

两式相减得：
其中 Δ x ( k ) = x ( k ) − x ( k − 1 ) Δx(k)=x(k)-x(k-1) Δx(k)=x(k)−x(k−1) , Δ u ( k ) = u ( k ) − u ( k − 1 ) Δu(k)=u(k)-u(k-1) Δu(k)=u(k)−u(k−1) , Δ d ( k ) = d ( k ) − d ( k − 1 ) Δd(k)=d(k)-d(k-1) Δd(k)=d(k)−d(k−1) 。

那输出y可由:

合并得：

以最新测量值为初始条件，预测时域为p，控制时域为m，并且有以下两个假设：

则：

其中 Δ x ( k + 1 ∣ k ) Δx(k+1|k) Δx(k+1∣k) 表示第k个时刻对第k+1个时刻的状态预测，那么接下来预测域 p p p长个状态表示为：。

进一步输出方程可以为

定义以下向量：

那么，对系统未来p步预测的输出：

其中：

这里 I I I中ne指的是： C ∗ A ∗ B C*A*B C∗A∗B矩阵的维度，是为了能够将 y ( k ) y(k) y(k)和其他项能够合并。

优化问题描述：

对于此优化问题考虑两个方面，一是输出跟踪上参考输入；二是尽可能减少控制幅度。目标函数选择如下：

其中 Γ y , i = def d i a g ( Γ y 1 , i \Gamma_{y,i}\stackrel{\text { def }}{=}diag(\Gamma_{y_1,i} Γy,i= def diag(Γy1,i, Γ y 2 , i , ⋯ , Γ y n c , i ) \Gamma_{y_2,i},\cdots,\Gamma_{y_{n_c},i}) Γy2,i,⋯,Γync,i)， Γ u , i = def d i a g ( Γ u 1 , i \Gamma_{u,i}\stackrel{\text { def }}{=}diag(\Gamma_{u_1,i} Γu,i= def diag(Γu1,i, Γ u 2 , i \Gamma_{u_2,i} Γu2,i, ⋯ , Γ u n u , i ) \cdots,\Gamma_{u_{n_u},i}) ⋯,Γunu,i)。 Γ y j , i \Gamma_{y_j,i} Γyj,i是第i步预测的第j个输出的误差的加权因子， Γ u j , i \Gamma_{u_j,i} Γuj,i是第i步预测对控制增量第j个分量的加权因子。改写为矩阵形式：

其中： Γ y = d i a g ( Γ y , 1 \Gamma_y=diag(\Gamma_{y,1} Γy=diag(Γy,1, Γ y , 2 \Gamma_{y,2} Γy,2, ⋯ , Γ y , p ) \cdots,\Gamma_{y,p}) ⋯,Γy,p)， Γ u = d i a g ( Γ u , 1 \Gamma_u=diag(\Gamma_{u,1} Γu=diag(Γu,1, Γ u , 2 \Gamma_{u,2} Γu,2, ⋯ , Γ u , p ) \cdots,\Gamma_{u,p}) ⋯,Γu,p)，参考输入序列为

问题可以被描述为： min ⁡ Δ U ( k ) J ( x ( k ) , Δ U ( k ) , m , p ) \min _{\Delta U(k)} J(x(k), \Delta U(k), m, p) minΔU(k)J(x(k),ΔU(k),m,p)
满足动力学模型： Y p ( k + 1 ∣ k ) = S x Δ x ( k ) + I y c ( k ) + S d Δ d ( k ) + S u Δ U ( k ) Y_p(k+1|k)=S_x\Delta x(k)+\mathcal {I} y_c(k)+S_d\Delta d(k)+S_u\Delta U(k) Yp(k+1∣k)=SxΔx(k)+Iyc(k)+SdΔd(k)+SuΔU(k)

定义一个辅助变量

将预测方程带入可得

其中：

E p ( k + 1 ∣ k ) = R ( k + 1 ) − S x Δ x ( k ) − I y c ( k ) − S d Δ d ( k ) E_{p}(k+1 | k)=R(k+1)-\mathcal{S}_{x} \Delta x(k)-\mathcal{I} y_{c}(k)-\mathcal{S}_{d} \Delta d(k) Ep(k+1∣k)=R(k+1)−SxΔx(k)−Iyc(k)−SdΔd(k)
因此 min ⁡ z ρ T ρ = ( A z − b ) T ( A z − b ) \min _{z} \rho^{\mathrm{T}} \rho=(\mathcal{A}z-b)^T(\mathcal {A}z-b) minzρTρ=(Az−b)T(Az−b)的极值解为 z ∗ = ( A T A ) − 1 A T b z^*=(\mathcal A^T\mathcal A)^{-1}\mathcal A^Tb z∗=(ATA)−1ATb。（充要条件自证）

最优控制序列为：

解第一个元素

Δ u ( k ) = [ I n u × n u 0 ⋯ 0 ] 1 × m Δ U ∗ ( k ) \Delta u(k)=\left[\begin{array}{cccc} I_{n_{u} \times n_{u}} & \mathbf{0} & \cdots & \mathbf{0} \end{array}\right]_{1 \times m} \Delta U^{*}(k) Δu(k)=[Inu×nu0⋯0]1×mΔU∗(k)
定义 K m p c = [ I n u × n u 0 ⋯ 0 ] 1 × m ( S u T Γ y T Γ y S u + Γ u T Γ u ) − 1 S u T Γ y T Γ y K_{\mathrm{mpc}}=\left[\begin{array}{cccc} I_{n_{u} \times n_{u}} & \mathbf{0} & \cdots & \mathbf{0} \end{array}\right]_{1 \times m}\left(\mathcal{S}_{u}^{\mathrm{T}} \Gamma_{y}^{\mathrm{T}} \Gamma_{y} \mathcal{S}_{u}+\Gamma_{u}^{\mathrm{T}} \Gamma_{u}\right)^{-1} \mathcal{S}_{u}^{\mathrm{T}} \Gamma_{y}^{\mathrm{T}} \Gamma_{y} Kmpc=[Inu×nu0⋯0]1×m(SuTΓyTΓySu+ΓuTΓu)−1SuTΓyTΓy，因此控制增量可以写为

Δ u ( k ) = K m p c E p ( k + 1 ∣ k ) \Delta u(k)=K_{\mathrm{mpc}} E_{p}(k+1 | k) Δu(k)=KmpcEp(k+1∣k)

理解：

将 E p ( k + 1 ∣ k ) = R ( k + 1 ) − S x Δ x ( k ) − I y c ( k ) − S d Δ d ( k ) E_{p}(k+1 | k)=R(k+1)-\mathcal{S}_{x} \Delta x(k)-\mathcal{I} y_{c}(k)-\mathcal{S}_{d} \Delta d(k) Ep(k+1∣k)=R(k+1)−SxΔx(k)−Iyc(k)−SdΔd(k)带入 Δ u ( k ) = K m p c E p ( k + 1 ∣ k ) \Delta u(k)=K_{\mathrm{mpc}} E_{p}(k+1 | k) Δu(k)=KmpcEp(k+1∣k)，可得

Δ u ( k ) = K m p c R ( k + 1 ) − K m p c ( S x + I C c ) Δ x ^ ( k ) − K m p c I y ^ c ( k − 1 ) − K m p c S d Δ d ( k ) \begin{aligned} \Delta u(k)=& K_{\mathrm{mpc}} R(k+1)-K_{\mathrm{mpc}}\left(\mathcal{S}_{x}+\mathcal{I} C_{c}\right) \Delta \hat{x}(k) -K_{\mathrm{mpc}} \mathcal{I} \hat{y}_{c}(k-1)-K_{\mathrm{mpc}} \mathcal{S}_{d} \Delta d(k) \end{aligned} Δu(k)=KmpcR(k+1)−Kmpc(Sx+ICc)Δx^(k)−KmpcIy^c(k−1)−KmpcSdΔd(k)

K m p c R ( k + 1 ) K_{\mathrm{mpc}} R(k+1) KmpcR(k+1)相当于未来参考输入的前馈补偿；
− K m p c S d Δ d ( k ) -K_{\mathrm{mpc}} \mathcal{S}_{d} \Delta d(k) −KmpcSdΔd(k)相当于可预测扰动的前馈补偿；
− K m p c ( S x + I C c ) Δ x ^ ( k ) − K m p c I y ^ c ( k − 1 ) -K_{\mathrm{mpc}}\left(\mathcal{S}_{x}+\mathcal{I} C_{c}\right) \Delta \hat{x}(k) -K_{\mathrm{mpc}} \mathcal{I} \hat{y}_{c}(k-1) −Kmpc(Sx+ICc)Δx^(k)−KmpcIy^c(k−1)相当于状态反馈补偿。

初学，转发保存。