求解MPC: 在滚动时间窗内建立并求解QP问题

该小节旨在根据上一节得到的离散车辆横向动力学模型,在滚动时间窗内建立并求解QP问题,实现MPC横向控制.
在开始求解之前,我们需要回答以下两个问题:
1.什么是QP问题?其标准形式的什么样的?2.怎么把MPC中的优化问题转化为求解QP问题?\begin{aligned} &1.什么是QP问题?其标准形式的什么样的? \\ &2.怎么把MPC中的优化问题转化为求解QP问题? \end{aligned} 1.什么是QP问题?其标准形式的什么样的?2.怎么把MPC中的优化问题转化为求解QP问题?
简要回答第一个问题: QP(Quadratic programming),即二次规划.旨在求解最小化一个带线性约束的多元二次型代价函数的最优变量值.其标准形式为
minxJ=minx12xTQx+xTgs.t.Ax≤b\begin{aligned} &\mathop{min}\limits_{x}\quad J \quad=\quad \mathop{min}\limits_{x} &&\frac{1}{2}x^TQx+x^Tg\\ &s.t. && Ax \leq b\\ \end{aligned} xminJ=xmins.t.21xTQx+xTgAx≤b
其中,当且仅当 Q⪰0Q \succeq 0Q⪰0,即 QQQ 为半正定矩阵时,QP问题为一个凸优化问题.

目前有很多开源的QP求解库,如 OSQP/qpOASESOSQP / qpOASESOSQP/qpOASES,以及 MatlabMatlabMatlab 中的 quadprogquadprogquadprog 都为以上QP问题提供了求解的API.

接下来回答第二个问题,如何将MPC转化为求解QP问题?
通常的,MPC控制算法包含以下四个基本步骤:
1.根据已知的系统模型,控制序列与当前的系统状态变量,对预测时域内各状态变量作出预测2.根据第一步的预测结果,计算预测时域内,满足约束的条件,最小化损失函数 J的控制序列3.将第二步中得到的控制序列的第一个控制量输入到系统中4.在下一次采样时间,将时间往前推动一步,更新第一步中的当前的系统状态变量,重复步骤1-3\begin{aligned} &\text{1.根据已知的系统模型,控制序列与当前的系统状态变量,对预测时域内各状态变量作出预测}\\ &\text{2.根据第一步的预测结果,计算预测时域内,满足约束的条件,最小化损失函数 $J$ 的控制序列}\\ &\text{3.将第二步中得到的控制序列的第一个控制量输入到系统中}\\ &\text{4.在下一次采样时间,将时间往前推动一步,更新第一步中的当前的系统状态变量,重复步骤1-3} \end{aligned} 1.根据已知的系统模型,控制序列与当前的系统状态变量,对预测时域内各状态变量作出预测2.根据第一步的预测结果,计算预测时域内,满足约束的条件,最小化损失函数 J 的控制序列3.将第二步中得到的控制序列的第一个控制量输入到系统中4.在下一次采样时间,将时间往前推动一步,更新第一步中的当前的系统状态变量,重复步骤1-3
步骤一中,我们需要根据当前系统状态和已知的控制序列对整个时间窗内的状态变量值作出预测. 我们可以通过迭代系统的离散模型得到 t+kt+kt+k 时刻的状态变量与 ttt 时刻(当前时刻)的状态变量和控制序列之间的关系,即系统的预测方程.

步骤二中,我们将时间窗内的系统性能期望构造为一个二次型损失函数,并将控制过程中所涉及到的约束条件都改写为 Ax⪯bAx \preceq bAx⪯b 的形式.至此,我们便构造了一个标准QP问题,接下来就是采用合适的优化算法(GD,Newton,ADMM等)对问题进行求解,得到时间窗内的最优控制序列.

步骤三中,我们只将控制序列中的第一个控制量作为系统的真实输入.采用这种策略,我认为原因有二.
原因一是因为步骤一和步骤二中的预测以及优化过程均为开环.在系统不确定度的影响下,这种计算是粗略的.且预测时域越长,与实际系统的偏差越大.因此,我们需要引入真实系统的反馈去修正预测.那么在当前时刻,新的系统反馈抵达之前,采用控制序列中的第一个控制量作为系统的真实输入是合适的.

略作延伸, 我们可以把MPC算法与Kalman滤波算法作对比.Kalman滤波算法中也有相似的预测过程,同理,若想要预测结果收敛于真值,单有预测是不够的.我们需要持续将最新的观测纳入考虑,利用其中蕴含的信息对预测过程中的相关量进行修正.Kalman的做法是通过贝叶斯定理,利用观测量对协方差矩阵进行修正.

那么MPC是如何利用最新的观测的呢? MPC的策略是在进入下一次采样时间后, 将最新的观测量转化为预测步骤的初始条件,然后在新的时间窗内进行预测,优化,即步骤四.这种方法称为滚动时域优化,因此MPC也称作滚动时域控制.
由此,我们便引出了原因二.通过滚动时域的方法,利用真实系统的反馈不断修正控制量,保障步骤三的操作可以满足最终的控制目标.而每一个时间片的控制过程又符合优化损失函数的预期.两者相互配合,形成闭环.

值得一提的是,滚动时域的策略使MPC的优化求解都是在一个有限时间的窗口下进行的,因此是无法保障得到整个时域的全局最优.但是因为持续优化,最终也能获得较好的结果. 虽然MPC策略损失了全局最优性,但其获得了更大灵活性.M因此相较于如LQR/LQG等最优控制器,MPC对于复杂非线性约束的对应更加自如,以及模型误差容忍度更高.

综上,回答第二个问题:通过滚动时域优化的方法,将MPC算法转化为在每一个时间窗下求解一个QP问题.

预测方程

我们将离散误差动力学模型改写为以控制量增量 Δδ\Delta\deltaΔδ 的输入的增量式模型.当然了你也可以不用增量式模型,直接用原有模型也是可以的.这并不影响预测方程的形式,只不过我选择了增量式模型作为例子,MPC的模型和损失函数的构造都是很灵活的.回到主题,增量式模型如下:
x(k+1)=Adx(k)+Bdδ(k)+Bcdψ˙des(k)=Adx(k)+Bd[δ(k−1)+Δδ(k)]+Bcdψ˙des(k)\begin{aligned} x(k+1) &= A_dx(k)+B_d\delta(k)+B_{cd}\dot{\psi}_{des}(k) \\ &= A_dx(k)+B_d[\delta(k-1)+\Delta\delta(k)]+B_{cd}\dot{\psi}_{des}(k) \end{aligned} x(k+1)=Adx(k)+Bdδ(k)+Bcdψ˙des(k)=Adx(k)+Bd[δ(k−1)+Δδ(k)]+Bcdψ˙des(k)

进一步的, 将以上模型改写为矩阵形式,可得
[x(k+1)δ(k)]=[AdBd0I][x(k)δ(k−1)]+[BdI]Δδ(k)+[Bcd0]ψ˙des(k)y(k)=[Cd0][x(k)δ(k−1)]\begin{aligned} &\begin{bmatrix} x(k+1) \\ \delta(k) \end{bmatrix} = \begin{bmatrix} A_d & B_d\\ 0 & I \end{bmatrix} \begin{bmatrix} x(k) \\ \delta(k-1) \end{bmatrix} + \begin{bmatrix} B_d \\ I \end{bmatrix}\Delta\delta(k)+ \begin{bmatrix} B_{cd} \\ 0 \end{bmatrix}\dot{\psi}_{des}(k)\\ \\ &y(k) = \begin{bmatrix} C_d & 0 \end{bmatrix} \begin{bmatrix} x(k) \\ \delta(k-1) \end{bmatrix} \end{aligned} [x(k+1)δ(k)]=[Ad0BdI][x(k)δ(k−1)]+[BdI]Δδ(k)+[Bcd0]ψ˙des(k)y(k)=[Cd0][x(k)δ(k−1)]

定义
ξ(k∣t)=[x(k)δ(k−1)]A~d=[AdBd0I]B~d=[BdI]B~cd=[Bcd0]C~d=[Cd0]\begin{aligned} &\xi(k|t) = \begin{bmatrix} x(k) \\ \delta(k-1) \end{bmatrix} \quad \widetilde{A}_d = \begin{bmatrix} A_d & B_d\\ 0 & I \end{bmatrix} \\ &\widetilde{B}_d = \begin{bmatrix} B_d \\ I \end{bmatrix} \quad \widetilde{B}_{cd} = \begin{bmatrix} B_{cd} \\ 0 \end{bmatrix} \quad \widetilde{C}_d = \begin{bmatrix} C_d & 0 \end{bmatrix} \end{aligned} ξ(k∣t)=[x(k)δ(k−1)]Ad=[Ad0BdI]Bd=[BdI]Bcd=[Bcd0]Cd=[Cd0]

其中, ξ(k∣t)\xi(k|t)ξ(k∣t) 代表在ttt时刻,预测kkk个周期得到的状态.为了简便, 在后续的使用中 ξ(k)\xi(k)ξ(k) 与 ξ(k∣t)\xi(k|t)ξ(k∣t) 通用.

综上,增量式离散误差动力学模型为
ξ(k+1)=A~dξ(k)+B~dΔδ(k)+B~cdψ˙des(k)y(k)=C~dξ(k)\begin{aligned} &\xi(k+1) = \widetilde{A}_d\xi(k)+\widetilde{B}_d\Delta\delta(k)+\widetilde{B}_{cd}\dot{\psi}_{des}(k)\\ &y(k) = \widetilde{C}_d\xi(k) \end{aligned} ξ(k+1)=Adξ(k)+BdΔδ(k)+Bcdψ˙des(k)y(k)=Cdξ(k)

假设预测步数为kkk,即MPC的预测时域为kkk个周期,则预测方程有
ξ(1)=A~dξ(0)+B~dΔδ(0)+B~cdψ˙des(0)ξ(2)=A~dξ(1)+B~dΔδ(1)+B~cdψ˙des(1)=A~d[A~dξ(0)+B~dΔδ(0)+B~cdψ˙des(0)]+B~dΔδ(1)+B~cdψ˙des(1)=A~d2ξ(0)+A~dB~dΔδ(0)+B~dΔδ(1)+A~dB~cdψ˙des(0)+B~cdψ˙des(1)ξ(3)=A~dξ(2)+B~dΔδ(2)+B~cdψ˙des(2)=A~d[A~d2ξ(0)+A~dB~dΔδ(0)+B~dΔδ(1)+A~dB~cdψ˙des(0)+B~cdψ˙des(1)]+B~dΔδ(2)+B~cdψ˙des(2)=A~d3ξ(0)+A~d2B~dΔδ(0)+A~dB~dΔδ(1)+B~dΔδ(2)+A~d2B~cdψ˙des(0)+A~dB~cdψ˙des(1)+B~cdψ˙des(2)...ξ(k)=A~d(k)ξ(0)+∑i=0k−1A~d(i)B~dΔδ(k−i−1)+∑j=0k−1A~d(j)B~cdψ˙des(k−j−1)y(k)=C~dξ(k)=C~dA~d(k)ξ(0)+C~d∑i=0k−1A~d(i)B~dΔδ(k−i−1)+C~d∑j=0k−1A~d(j)B~cdψ˙des(k−j−1)\begin{aligned} &\xi(1) &&=\widetilde{A}_d\xi(0)+\widetilde{B}_d\Delta\delta(0)+\widetilde{B}_{cd}\dot{\psi}_{des}(0) \\ \\ &\xi(2) &&=\widetilde{A}_d\xi(1)+\widetilde{B}_d\Delta\delta(1)+\widetilde{B}_{cd}\dot{\psi}_{des}(1) \\ & &&=\widetilde{A}_d[\widetilde{A}_d\xi(0)+\widetilde{B}_d\Delta\delta(0)+\widetilde{B}_{cd}\dot{\psi}_{des}(0)] +\widetilde{B}_d\Delta\delta(1)+\widetilde{B}_{cd}\dot{\psi}_{des}(1) \\ & &&=\widetilde{A}^2_d\xi(0)+\widetilde{A}_d\widetilde{B}_d\Delta\delta(0)+\widetilde{B}_d\Delta\delta(1)+ \widetilde{A}_d\widetilde{B}_{cd}\dot{\psi}_{des}(0)+\widetilde{B}_{cd}\dot{\psi}_{des}(1) \\ \\ &\xi(3) &&=\widetilde{A}_d\xi(2)+\widetilde{B}_d\Delta\delta(2)+\widetilde{B}_{cd}\dot{\psi}_{des}(2) \\ & &&=\widetilde{A}_d[\widetilde{A}^2_d\xi(0)+\widetilde{A}_d\widetilde{B}_d\Delta\delta(0)+\widetilde{B}_d\Delta\delta(1)+ \widetilde{A}_d\widetilde{B}_{cd}\dot{\psi}_{des}(0)+\widetilde{B}_{cd}\dot{\psi}_{des}(1)]+\widetilde{B}_d\Delta\delta(2)+\widetilde{B}_{cd}\dot{\psi}_{des}(2)\\ & &&=\widetilde{A}^3_d\xi(0)+\widetilde{A}^2_d\widetilde{B}_d\Delta\delta(0)+\widetilde{A}_d\widetilde{B}_d\Delta\delta(1)+\widetilde{B}_d\Delta\delta(2)+ \widetilde{A}^2_d\widetilde{B}_{cd}\dot{\psi}_{des}(0)+\widetilde{A}_d\widetilde{B}_{cd}\dot{\psi}_{des}(1)+\widetilde{B}_{cd}\dot{\psi}_{des}(2) \\ \\ & \ ... \\ \\ &\xi(k) &&= \widetilde{A}^{(k)}_d\xi(0)+\sum_{i=0}^{k-1}\widetilde{A}_d^{(i)}\widetilde{B}_d\Delta\delta(k-i-1)+ \sum_{j=0}^{k-1}\widetilde{A}_d^{(j)}\widetilde{B}_{cd}\dot{\psi}_{des}(k-j-1) \\ &y(k) &&=\widetilde{C}_d\xi(k)=\widetilde{C}_d\widetilde{A}^{(k)}_d\xi(0)+\widetilde{C}_d\sum_{i=0}^{k-1}\widetilde{A}_d^{(i)}\widetilde{B}_d\Delta\delta(k-i-1)+ \widetilde{C}_d\sum_{j=0}^{k-1}\widetilde{A}_d^{(j)}\widetilde{B}_{cd}\dot{\psi}_{des}(k-j-1) \end{aligned} ξ(1)ξ(2)ξ(3) ...ξ(k)y(k)=Adξ(0)+BdΔδ(0)+Bcdψ˙des(0)=Adξ(1)+BdΔδ(1)+Bcdψ˙des(1)=Ad[Adξ(0)+BdΔδ(0)+Bcdψ˙des(0)]+BdΔδ(1)+Bcdψ˙des(1)=Ad2ξ(0)+AdBdΔδ(0)+BdΔδ(1)+AdBcdψ˙des(0)+Bcdψ˙des(1)=Adξ(2)+BdΔδ(2)+Bcdψ˙des(2)=Ad[Ad2ξ(0)+AdBdΔδ(0)+BdΔδ(1)+AdBcdψ˙des(0)+Bcdψ˙des(1)]+BdΔδ(2)+Bcdψ˙des(2)=Ad3ξ(0)+Ad2BdΔδ(0)+AdBdΔδ(1)+BdΔδ(2)+Ad2Bcdψ˙des(0)+AdBcdψ˙des(1)+Bcdψ˙des(2)=Ad(k)ξ(0)+i=0∑k−1Ad(i)BdΔδ(k−i−1)+j=0∑k−1Ad(j)Bcdψ˙des(k−j−1)=Cdξ(k)=CdAd(k)ξ(0)+Cdi=0∑k−1Ad(i)BdΔδ(k−i−1)+Cdj=0∑k−1Ad(j)Bcdψ˙des(k−j−1)

基于以上的预测方程,我们可以通过ttt时刻的状态量初值 ξ(0)\xi(0)ξ(0),建立状态量 ξ\xiξ 在 kkk 个时域预测周期内与各周期控制增量 Δδ\Delta\deltaΔδ 之间的联系,即
[ξ(t+1∣t)ξ(t+2∣t)ξ(t+3∣t)⋮ξ(t+k∣t)]=[A~dA~d2A~d3⋮A~dk]ξ(0)+[B~d00…0A~dB~dB~d0…0A~d2B~dA~dB~dB~d…0⋮⋮⋮⋱⋮A~dk−1B~dA~dk−2B~dA~dk−3B~d…B~d][Δδ(0)Δδ(1)Δδ(2)⋮Δδ(k−1)]+[B~cd00…0A~dB~cdB~cd0…0A~d2B~cdA~dB~cdB~cd…0⋮⋮⋮⋱⋮A~dk−1B~cdA~dk−2B~cdA~dk−3B~cd…B~cd][ψ˙des(0)ψ˙des(1)ψ˙des(2)⋮ψ˙des(k−1)]\begin{bmatrix} \xi(t+1 | t) \\ \xi(t+2 | t) \\ \xi(t+3 | t) \\ \vdots \\ \xi(t+k | t)\end{bmatrix} = \begin{bmatrix} \widetilde{A}_d \\ \widetilde{A}^2_d \\ \widetilde{A}^3_d \\ \vdots \\ \widetilde{A}^k_d\end{bmatrix}\xi(0)+ \begin{bmatrix} \widetilde{B}_d & 0 & 0 & \dots & 0 \\ \widetilde{A}_d\widetilde{B}_d & \widetilde{B}_d & 0 & \dots & 0 \\ \widetilde{A}_d^2\widetilde{B}_d & \widetilde{A}_d\widetilde{B}_d & \widetilde{B}_d & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \widetilde{A}_d^{k-1}\widetilde{B}_d & \widetilde{A}_d^{k-2}\widetilde{B}_d &\widetilde{A}_d^{k-3}\widetilde{B}_d & \dots & \widetilde{B}_d \end{bmatrix} \begin{bmatrix} \Delta\delta(0) \\ \Delta\delta(1) \\ \Delta\delta(2) \\ \vdots \\ \Delta\delta(k-1) \end{bmatrix} + \begin{bmatrix} \widetilde{B}_{cd} & 0 & 0 & \dots & 0\\ \widetilde{A}_d\widetilde{B}_{cd} & \widetilde{B}_{cd} & 0 & \dots & 0 \\ \widetilde{A}_d^2\widetilde{B}_{cd} & \widetilde{A}_d\widetilde{B}_{cd} & \widetilde{B}_{cd} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \widetilde{A}_d^{k-1}\widetilde{B}_{cd} & \widetilde{A}_d^{k-2}\widetilde{B}_{cd} &\widetilde{A}_d^{k-3}\widetilde{B}_{cd} & \dots & \widetilde{B}_{cd} \end{bmatrix} \begin{bmatrix} \dot{\psi}_{des}(0) \\ \dot{\psi}_{des}(1) \\ \dot{\psi}_{des}(2) \\ \vdots \\ \dot{\psi}_{des}(k-1) \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎡ξ(t+1∣t)ξ(t+2∣t)ξ(t+3∣t)⋮ξ(t+k∣t)⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡AdAd2Ad3⋮Adk⎦⎥⎥⎥⎥⎥⎤ξ(0)+⎣⎢⎢⎢⎢⎢⎡BdAdBdAd2Bd⋮Adk−1Bd0BdAdBd⋮Adk−2Bd00Bd⋮Adk−3Bd………⋱…000⋮Bd⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡Δδ(0)Δδ(1)Δδ(2)⋮Δδ(k−1)⎦⎥⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎢⎡BcdAdBcdAd2Bcd⋮Adk−1Bcd0BcdAdBcd⋮Adk−2Bcd00Bcd⋮Adk−3Bcd………⋱…000⋮Bcd⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡ψ˙des(0)ψ˙des(1)ψ˙des(2)⋮ψ˙des(k−1)⎦⎥⎥⎥⎥⎥⎤

定义
X=[ξ(t+1∣t)ξ(t+2∣t)⋮ξ(t+k∣t)]ΔU=[Δδ(0)Δδ(1)⋮Δδ(k−1)]Λ=[A~dA~d2⋮A~dk]Ψ=[ψ˙des(0)ψ˙des(1)ψ˙des(2)⋮ψ˙des(k−1)]Γ=[B~d00…0A~dB~dB~d0…0A~d2B~dA~dB~dB~d…0⋮⋮⋮⋱⋮A~dk−1B~dA~dk−2B~dA~dk−3B~d…B~d]Υ=[B~cd00…0A~dB~cdB~cd0…0A~d2B~cdA~dB~cdB~cd…0⋮⋮⋮⋱⋮A~dk−1B~cdA~dk−2B~cdA~dk−3B~cd…B~cd]\begin{aligned} &X = \begin{bmatrix} \xi(t+1 | t) \\ \xi(t+2 | t) \\ \vdots \\ \xi(t+k | t)\end{bmatrix} \quad \Delta U = \begin{bmatrix} \Delta\delta(0) \\ \Delta\delta(1) \\ \vdots \\ \Delta\delta(k-1) \end{bmatrix} \quad \Lambda = \begin{bmatrix} \widetilde{A}_d \\ \widetilde{A}^2_d \\ \vdots \\ \widetilde{A}^k_d\end{bmatrix} \quad \Psi = \begin{bmatrix} \dot{\psi}_{des}(0) \\ \dot{\psi}_{des}(1) \\ \dot{\psi}_{des}(2) \\ \vdots \\ \dot{\psi}_{des}(k-1) \end{bmatrix}\\ \\ &\Gamma = \begin{bmatrix} \widetilde{B}_d & 0 & 0 & \dots & 0 \\ \widetilde{A}_d\widetilde{B}_d & \widetilde{B}_d & 0 & \dots & 0 \\ \widetilde{A}_d^2\widetilde{B}_d & \widetilde{A}_d\widetilde{B}_d & \widetilde{B}_d & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \widetilde{A}_d^{k-1}\widetilde{B}_d & \widetilde{A}_d^{k-2}\widetilde{B}_d &\widetilde{A}_d^{k-3}\widetilde{B}_d & \dots & \widetilde{B}_d \end{bmatrix}\\ \\ &\Upsilon = \begin{bmatrix} \widetilde{B}_{cd} & 0 & 0 & \dots & 0\\ \widetilde{A}_d\widetilde{B}_{cd} & \widetilde{B}_{cd} & 0 & \dots & 0 \\ \widetilde{A}_d^2\widetilde{B}_{cd} & \widetilde{A}_d\widetilde{B}_{cd} & \widetilde{B}_{cd} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \widetilde{A}_d^{k-1}\widetilde{B}_{cd} & \widetilde{A}_d^{k-2}\widetilde{B}_{cd} &\widetilde{A}_d^{k-3}\widetilde{B}_{cd} & \dots & \widetilde{B}_{cd} \end{bmatrix} \end{aligned} X=⎣⎢⎢⎢⎡ξ(t+1∣t)ξ(t+2∣t)⋮ξ(t+k∣t)⎦⎥⎥⎥⎤ΔU=⎣⎢⎢⎢⎡Δδ(0)Δδ(1)⋮Δδ(k−1)⎦⎥⎥⎥⎤Λ=⎣⎢⎢⎢⎡AdAd2⋮Adk⎦⎥⎥⎥⎤Ψ=⎣⎢⎢⎢⎢⎢⎡ψ˙des(0)ψ˙des(1)ψ˙des(2)⋮ψ˙des(k−1)⎦⎥⎥⎥⎥⎥⎤Γ=⎣⎢⎢⎢⎢⎢⎡BdAdBdAd2Bd⋮Adk−1Bd0BdAdBd⋮Adk−2Bd00Bd⋮Adk−3Bd………⋱…000⋮Bd⎦⎥⎥⎥⎥⎥⎤Υ=⎣⎢⎢⎢⎢⎢⎡BcdAdBcdAd2Bcd⋮Adk−1Bcd0BcdAdBcd⋮Adk−2Bcd00Bcd⋮Adk−3Bcd………⋱…000⋮Bcd⎦⎥⎥⎥⎥⎥⎤

进一步的,定义
Y=[y(t+1∣t)y(t+2∣t)⋮y(t+k∣t)]Ω=[C~dA~dC~dA~d2⋮C~dA~dk]Θ=[C~dB~d00…0C~dA~dB~dC~dB~d0…0C~dA~d2B~dC~dA~dB~dC~dB~d…0⋮⋮⋮⋱⋮C~dA~dk−1B~dC~dA~dk−2B~dC~dA~dk−3B~d…C~dB~d]Φ=[C~dB~cd00…0C~dA~dB~cdC~dB~cd0…0C~dA~d2B~cdC~dA~dB~cdC~dB~cd…0⋮⋮⋮⋱⋮C~dA~dk−1B~cdC~dA~dk−2B~cdC~dA~dk−3B~cd…C~dB~cd]\begin{aligned} &Y = \begin{bmatrix} y(t+1 | t) \\ y(t+2 | t) \\ \vdots \\ y(t+k | t)\end{bmatrix} \quad \Omega = \begin{bmatrix} \widetilde{C}_d\widetilde{A}_d \\ \widetilde{C}_d\widetilde{A}^2_d \\ \vdots \\ \widetilde{C}_d\widetilde{A}^k_d\end{bmatrix} \\ \\ &\Theta = \begin{bmatrix} \widetilde{C}_d\widetilde{B}_d & 0 & 0 & \dots & 0 \\ \widetilde{C}_d\widetilde{A}_d\widetilde{B}_d & \widetilde{C}_d\widetilde{B}_d & 0 & \dots & 0 \\ \widetilde{C}_d\widetilde{A}_d^2\widetilde{B}_d & \widetilde{C}_d\widetilde{A}_d\widetilde{B}_d & \widetilde{C}_d\widetilde{B}_d & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \widetilde{C}_d\widetilde{A}_d^{k-1}\widetilde{B}_d & \widetilde{C}_d\widetilde{A}_d^{k-2}\widetilde{B}_d &\widetilde{C}_d\widetilde{A}_d^{k-3}\widetilde{B}_d & \dots & \widetilde{C}_d\widetilde{B}_d \end{bmatrix}\\ \\ &\Phi = \begin{bmatrix} \widetilde{C}_d\widetilde{B}_{cd} & 0 & 0 & \dots & 0\\ \widetilde{C}_d\widetilde{A}_d\widetilde{B}_{cd} & \widetilde{C}_d\widetilde{B}_{cd} & 0 & \dots & 0 \\ \widetilde{C}_d\widetilde{A}_d^2\widetilde{B}_{cd} & \widetilde{C}_d\widetilde{A}_d\widetilde{B}_{cd} & \widetilde{C}_d\widetilde{B}_{cd} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \widetilde{C}_d\widetilde{A}_d^{k-1}\widetilde{B}_{cd} & \widetilde{C}_d\widetilde{A}_d^{k-2}\widetilde{B}_{cd} &\widetilde{C}_d\widetilde{A}_d^{k-3}\widetilde{B}_{cd} & \dots & \widetilde{C}_d\widetilde{B}_{cd} \end{bmatrix} \end{aligned} Y=⎣⎢⎢⎢⎡y(t+1∣t)y(t+2∣t)⋮y(t+k∣t)⎦⎥⎥⎥⎤Ω=⎣⎢⎢⎢⎡CdAdCdAd2⋮CdAdk⎦⎥⎥⎥⎤Θ=⎣⎢⎢⎢⎢⎢⎡CdBdCdAdBdCdAd2Bd⋮CdAdk−1Bd0CdBdCdAdBd⋮CdAdk−2Bd00CdBd⋮CdAdk−3Bd………⋱…000⋮CdBd⎦⎥⎥⎥⎥⎥⎤Φ=⎣⎢⎢⎢⎢⎢⎡CdBcdCdAdBcdCdAd2Bcd⋮CdAdk−1Bcd0CdBcdCdAdBcd⋮CdAdk−2Bcd00CdBcd⋮CdAdk−3Bcd………⋱…000⋮CdBcd⎦⎥⎥⎥⎥⎥⎤

综上有
X=Λξ(0)+ΓΔU+ΥΨY=Ωξ(0)+ΘΔU+ΦΨ\begin{aligned} &X=\Lambda\xi(0)+\Gamma\Delta U+\Upsilon\Psi\\ &Y=\Omega\xi(0)+\Theta\Delta U +\Phi\Psi \end{aligned} X=Λξ(0)+ΓΔU+ΥΨY=Ωξ(0)+ΘΔU+ΦΨ

代价函数

我构造的代价函数惩罚输出量以及控制增量,因此
J=∣∣Y−Ydes∣∣Qc+∣∣ΔU∣∣Rc=(Y−Ydes)TQc(Y−Ydes)+ΔUTRcΔU=(Ωξ(0)+ΘΔU+ΦΨ−Ydes)TQc(Ωξ(0)+ΘΔU+ΦΨ−Ydes)+ΔUTRcΔU\begin{aligned} J &= ||Y-Y_{des}||_{Q_c}+||\Delta U||_{R_c}\\ &=(Y-Y_{des})^TQ_c(Y-Y_{des})+\Delta U^TR_c\Delta U\\ &=(\Omega\xi(0)+\Theta\Delta U+\Phi\Psi-Y_{des})^TQ_c(\Omega\xi(0)+\Theta\Delta U+\Phi\Psi-Y_{des})+\Delta U^TR_c\Delta U \end{aligned} J=∣∣Y−Ydes∣∣Qc+∣∣ΔU∣∣Rc=(Y−Ydes)TQc(Y−Ydes)+ΔUTRcΔU=(Ωξ(0)+ΘΔU+ΦΨ−Ydes)TQc(Ωξ(0)+ΘΔU+ΦΨ−Ydes)+ΔUTRcΔU

我们优化的对象参数是 ΔU\Delta UΔU,因此将代价函数中不含对象参数的部分定义为 NNN,则
N=Ωξ(0)+ΦΨ−YdesN = \Omega\xi(0)+\Phi\Psi-Y_{des} N=Ωξ(0)+ΦΨ−Ydes

将 NNN 代入到代价函数中可得
J=(N+ΘΔU)TQc(N+ΘΔU)+ΔUTRcΔU=NTQcN+2ΔUTΘTQcN+ΔUTΘtQcΘΔU+ΔUTRcΔU=NTQcN+2ΔUTΘTQcN+ΔUT(ΘtQcΘ+Rc)ΔU\begin{aligned}J &= (N+\Theta\Delta U)^TQ_c(N+\Theta\Delta U)+\Delta U^TR_c \Delta U\\ &= N^TQ_cN + 2\Delta U^T\Theta^TQ_cN + \Delta U^T \Theta^tQ_c\Theta\Delta U + \Delta U^TR_c\Delta U\\ &= N^TQ_cN + 2\Delta U^T\Theta^TQ_cN + \Delta U^T(\Theta^tQ_c\Theta+R_c)\Delta U \end{aligned} J=(N+ΘΔU)TQc(N+ΘΔU)+ΔUTRcΔU=NTQcN+2ΔUTΘTQcN+ΔUTΘtQcΘΔU+ΔUTRcΔU=NTQcN+2ΔUTΘTQcN+ΔUT(ΘtQcΘ+Rc)ΔU

针对该优化问题,我们可以进一步简化代价函数JJJ,即消去代价函数中不含优化对象参数 ΔU\Delta UΔU的项
argminΔU(J)=argminΔU(NTQcN+2ΔUTΘTQcN+ΔUT(ΘtQcΘ+Rc)ΔU)=argminΔU(2ΔUTΘTQcN+ΔUT(ΘtQcΘ+Rc)ΔU)=argminΔU(J∗)\begin{aligned} \mathop{argmin}\limits_{\Delta U}(J) &=\mathop{argmin}\limits_{\Delta U}(N^TQ_cN+2\Delta U^T\Theta^TQ_cN + \Delta U^T(\Theta^tQ_c\Theta+R_c)\Delta U)\\ &=\mathop{argmin}\limits_{\Delta U}(2\Delta U^T\Theta^TQ_cN + \Delta U^T(\Theta^tQ_c\Theta+R_c)\Delta U)\\ &=\mathop{argmin}\limits_{\Delta U}(J^*) \end{aligned} ΔUargmin(J)=ΔUargmin(NTQcN+2ΔUTΘTQcN+ΔUT(ΘtQcΘ+Rc)ΔU)=ΔUargmin(2ΔUTΘTQcN+ΔUT(ΘtQcΘ+Rc)ΔU)=ΔUargmin(J∗)
其中, J∗J^*J∗ 为简化后的代价函数

定义
H=ΘtQcΘ+RcK=2ΘTQcN\begin{aligned} &H = \Theta^tQ_c\Theta+R_c \\ &K = 2\Theta^TQ_cN \end{aligned} H=ΘtQcΘ+RcK=2ΘTQcN

则QP问题为
minΔU(J∗)=minΔU(ΔUTHΔU+ΔUTK)\mathop{min}\limits_{\Delta U}(J^*) =\mathop{min}\limits_{\Delta U}(\Delta U^TH\Delta U+ \Delta U^TK) ΔUmin(J∗)=ΔUmin(ΔUTHΔU+ΔUTK)

约束条件

控制增量约束

在增量式模型中,我们优化的对象参数是控制量量增量 ΔU\Delta UΔU,因此对控制量增量直接约束为
ΔUmin⪯ΔU⪯ΔUmax\Delta U^{min} \preceq \Delta U \preceq \Delta U^{max} ΔUmin⪯ΔU⪯ΔUmax
其中
ΔUmin=[Δδ(0)minΔδ(1)min…Δδ(k−1)min]TΔUmax=[Δδ(0)maxΔδ(1)max…Δδ(k−1)max]T\begin{aligned} &\Delta U^{min} = \begin{bmatrix} \Delta\delta(0)^{min} & \Delta\delta(1)^{min} & \dots & \Delta\delta(k-1)^{min} \end{bmatrix} ^T \\ &\Delta U^{max} = \begin{bmatrix} \Delta\delta(0)^{max} & \Delta\delta(1)^{max} & \dots & \Delta\delta(k-1)^{max} \end{bmatrix} ^T \end{aligned} ΔUmin=[Δδ(0)minΔδ(1)min…Δδ(k−1)min]TΔUmax=[Δδ(0)maxΔδ(1)max…Δδ(k−1)max]T
以上约束条件可改写成

[−10…000−1…00⋮⋮⋱⋮⋮00…−1000…0−110…0001…00⋮⋮⋱⋮⋮00…1000…01][Δδ(0)Δδ(1)⋮Δδ(k−2)Δδ(k−1)]⪯[−Δδ(0)min−Δδ(1)min⋮−Δδ(k−2)min−Δδ(k−1)minΔδ(0)maxΔδ(1)max⋮Δδ(k−2)maxΔδ(k−1)max]\left[\begin{array}{rrrrr} -1 & 0 & \dots & 0 & 0\\ 0 & -1 & \dots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \dots & -1 & 0\\ 0 & 0 & \dots & 0 & -1\\ 1 & 0 & \dots & 0 & 0\\ 0 & 1 & \dots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \dots & 1 & 0\\ 0 & 0 & \dots & 0 & 1 \end{array}\right] \begin{bmatrix} \Delta\delta(0)\\ \Delta\delta(1)\\ \vdots \\ \Delta\delta(k-2)\\ \Delta\delta(k-1) \end{bmatrix} \preceq \begin{bmatrix} -\Delta\delta(0)^{min}\\ -\Delta\delta(1)^{min}\\ \vdots\\ -\Delta\delta(k-2)^{min}\\ -\Delta\delta(k-1)^{min}\\ \Delta\delta(0)^{max}\\ \Delta\delta(1)^{max}\\ \vdots\\ \Delta\delta(k-2)^{max}\\ \Delta\delta(k-1)^{max}\\ \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡−10⋮0010⋮000−1⋮0001⋮00……⋱…………⋱……00⋮−1000⋮1000⋮0−100⋮01⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡Δδ(0)Δδ(1)⋮Δδ(k−2)Δδ(k−1)⎦⎥⎥⎥⎥⎥⎤⪯⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡−Δδ(0)min−Δδ(1)min⋮−Δδ(k−2)min−Δδ(k−1)minΔδ(0)maxΔδ(1)max⋮Δδ(k−2)maxΔδ(k−1)max⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

定义
Sdul=[−10…000−1…00⋮⋮⋱⋮⋮00…−1000…0−110…0001…00⋮⋮⋱⋮⋮00…1000…01]=[−II]Sdur=[−Δδ(0)min−Δδ(1)min⋮−Δδ(k−2)min−Δδ(k−1)minΔδ(0)maxΔδ(1)max⋮Δδ(k−2)maxΔδ(k−1)max]=[−ΔUminΔUmax]\begin{aligned} &S^{l}_{du} = \left[\begin{array}{rrrrr} -1 & 0 & \dots & 0 & 0\\ 0 & -1 & \dots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \dots & -1 & 0\\ 0 & 0 & \dots & 0 & -1\\ 1 & 0 & \dots & 0 & 0\\ 0 & 1 & \dots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \dots & 1 & 0\\ 0 & 0 & \dots & 0 & 1 \end{array}\right] = \begin{bmatrix} -I \\ \quad I \end{bmatrix}\\ &S^{r}_{du} = \begin{bmatrix} -\Delta\delta(0)^{min}\\ -\Delta\delta(1)^{min}\\ \vdots\\ -\Delta\delta(k-2)^{min}\\ -\Delta\delta(k-1)^{min}\\ \Delta\delta(0)^{max}\\ \Delta\delta(1)^{max}\\ \vdots\\ \Delta\delta(k-2)^{max}\\ \Delta\delta(k-1)^{max}\\ \end{bmatrix} = \begin{bmatrix} -\Delta U^{min} \\ \quad\Delta U^{max} \end{bmatrix} \end{aligned} Sdul=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡−10⋮0010⋮000−1⋮0001⋮00……⋱…………⋱……00⋮−1000⋮1000⋮0−100⋮01⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=[−II]Sdur=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡−Δδ(0)min−Δδ(1)min⋮−Δδ(k−2)min−Δδ(k−1)minΔδ(0)maxΔδ(1)max⋮Δδ(k−2)maxΔδ(k−1)max⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=[−ΔUminΔUmax]
则控制增量的约束可表示为
SdulΔU⪯SdurS^{l}_{du} \Delta U \preceq S^{r}_{du} SdulΔU⪯Sdur

控制量约束

首先我们要将控制量矩阵 UUU 以控制增量 ΔU\Delta UΔU 的形式表达
U=[δ(0)δ(1)⋮δ(k−2)δ(k−1)]=[11⋮11]δ(−1)+[100…00110…00111…00⋮⋮⋮⋱⋮⋮111…10111…11][Δδ(0)Δδ(1)⋮Δδ(k−2)Δδ(k−1)]U = \begin{bmatrix} \delta(0)\\ \delta(1)\\ \vdots \\ \delta(k-2)\\ \delta(k-1) \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \\ 1 \end{bmatrix} \delta(-1)+ \begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0\\ 1 & 1 & 0 & \dots & 0 & 0\\ 1 & 1 & 1 & \dots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots &\vdots & \vdots\\ 1 & 1 & 1 & \dots & 1 & 0\\ 1 & 1 & 1 & \dots & 1 & 1 \end{bmatrix} \begin{bmatrix} \Delta\delta(0)\\ \Delta\delta(1)\\ \vdots \\ \Delta\delta(k-2)\\ \Delta\delta(k-1) \end{bmatrix} U=⎣⎢⎢⎢⎢⎢⎡δ(0)δ(1)⋮δ(k−2)δ(k−1)⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡11⋮11⎦⎥⎥⎥⎥⎥⎤δ(−1)+⎣⎢⎢⎢⎢⎢⎢⎢⎡111⋮11011⋮11001⋮11………⋱……000⋮11000⋮01⎦⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡Δδ(0)Δδ(1)⋮Δδ(k−2)Δδ(k−1)⎦⎥⎥⎥⎥⎥⎤
其中 δ(−1)\delta(-1)δ(−1) 代表 t−1t-1t−1 时刻的控制量,即上一周期的控制量

定义
E=[11⋮11]F=[100…00110…00111…00⋮⋮⋮⋱⋮⋮111…10111…11]E = \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \\ 1 \end{bmatrix} \quad F = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0\\ 1 & 1 & 0 & \dots & 0 & 0\\ 1 & 1 & 1 & \dots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots &\vdots & \vdots\\ 1 & 1 & 1 & \dots & 1 & 0\\ 1 & 1 & 1 & \dots & 1 & 1 \end{bmatrix} E=⎣⎢⎢⎢⎢⎢⎡11⋮11⎦⎥⎥⎥⎥⎥⎤F=⎣⎢⎢⎢⎢⎢⎢⎢⎡111⋮11011⋮11001⋮11………⋱……000⋮11000⋮01⎦⎥⎥⎥⎥⎥⎥⎥⎤

则上式可改写为
U=Eδ(−1)+FΔUU = E\delta(-1)+F\Delta U U=Eδ(−1)+FΔU

控制量 UUU 需满足的约束条件为
Umin⪯U⪯UmaxU^{min} \preceq U \preceq U^{max} Umin⪯U⪯Umax
其中
Umin=[δ(0)minδ(1)min…δ(k−1)min]TUmax=[δ(0)maxδ(1)max…δ(k−1)max]T\begin{aligned} &U^{min} = \begin{bmatrix} \delta(0)^{min} & \delta(1)^{min} & \dots & \delta(k-1)^{min} \end{bmatrix} ^T \\ & U^{max} = \begin{bmatrix} \delta(0)^{max} & \delta(1)^{max} & \dots & \delta(k-1)^{max} \end{bmatrix} ^T \end{aligned} Umin=[δ(0)minδ(1)min…δ(k−1)min]TUmax=[δ(0)maxδ(1)max…δ(k−1)max]T
根据以上关系,改写控制量约束
[−FF]ΔU⪯[−Umin+Eδ(−1)Umax−Eδ(−1)]\begin{bmatrix} -F\\ \quad F \end{bmatrix} \Delta U \preceq \begin{bmatrix} -U^{min}+E\delta(-1)\\ \quad U^{max}-E\delta(-1) \end{bmatrix} [−FF]ΔU⪯[−Umin+Eδ(−1)Umax−Eδ(−1)]

定义
Sul=[−FF]Sur=[−Umin+Eδ(−1)Umax−Eδ(−1)]\begin{aligned} &S^{l}_{u} = \begin{bmatrix} -F\\ \quad F \end{bmatrix}\\ &S^{r}_{u} = \begin{bmatrix} -U^{min}+E\delta(-1)\\ \quad U^{max}-E\delta(-1) \end{bmatrix} \end{aligned} Sul=[−FF]Sur=[−Umin+Eδ(−1)Umax−Eδ(−1)]

则控制增量的约束可表示为

SulΔU⪯SurS^{l}_{u} \Delta U \preceq S^{r}_{u} SulΔU⪯Sur

状态量约束与输出量约束

状态量约束和输出量约束我们可以直接根据增量式动力学模型推导,步骤如下
X=Λξ(0)+ΓΔU+ΥΨY=Ωξ(0)+ΘΔU+ΦΨXlb−Λξ(0)−ΥΨ⪯ΓΔU⪯Xub−Λξ(0)−ΥΨYlb−Ωξ(0)−ΦΨ⪯ΘΔU⪯Yub−Ωξ(0)−ΦΨ[−ΓΓ]ΔU⪯[−Xlb+Λξ(0)+ΥΨXub−Λξ(0)−ΥΨ][−ΘΘ]ΔU⪯[−Ylb+Ωξ(0)+ΦΨYub−Ωξ(0)−ΦΨ]\begin{aligned} &X=\Lambda\xi(0)+\Gamma\Delta U+\Upsilon\Psi\\ &Y=\Omega\xi(0)+\Theta\Delta U +\Phi\Psi \\ \\ &X_{lb}-\Lambda\xi(0)-\Upsilon\Psi \preceq \Gamma\Delta U \preceq X_{ub}-\Lambda\xi(0)-\Upsilon\Psi\\ &Y_{lb}-\Omega\xi(0)-\Phi\Psi \preceq \Theta\Delta U \preceq Y_{ub}-\Omega\xi(0)-\Phi\Psi\\\\ &\begin{bmatrix} -\Gamma\\ \quad\Gamma \end{bmatrix} \Delta U \preceq \begin{bmatrix} -X_{lb}+\Lambda\xi(0)+\Upsilon\Psi \\ \quad X_{ub}-\Lambda\xi(0)-\Upsilon\Psi \end{bmatrix}\\ &\begin{bmatrix} -\Theta\\ \quad\Theta \end{bmatrix} \Delta U \preceq \begin{bmatrix} -Y_{lb}+\Omega\xi(0)+\Phi\Psi \\ \quad Y_{ub}-\Omega\xi(0)-\Phi\Psi \end{bmatrix} \end{aligned} X=Λξ(0)+ΓΔU+ΥΨY=Ωξ(0)+ΘΔU+ΦΨXlb−Λξ(0)−ΥΨ⪯ΓΔU⪯Xub−Λξ(0)−ΥΨYlb−Ωξ(0)−ΦΨ⪯ΘΔU⪯Yub−Ωξ(0)−ΦΨ[−ΓΓ]ΔU⪯[−Xlb+Λξ(0)+ΥΨXub−Λξ(0)−ΥΨ][−ΘΘ]ΔU⪯[−Ylb+Ωξ(0)+ΦΨYub−Ωξ(0)−ΦΨ]

定义

SXl=[−ΓΓ]SXr=[−Xlb+Λξ(0)+ΥΨXub−Λξ(0)−ΥΨ]SYl=[−ΘΘ]SYr=[−Ylb+Ωξ(0)+ΦΨYub−Ωξ(0)−ΦΨ]\begin{aligned} &S^{l}_{X} = \begin{bmatrix} -\Gamma\\ \quad \Gamma \end{bmatrix}\quad S^{r}_{X} = \begin{bmatrix} -X_{lb}+\Lambda\xi(0)+\Upsilon\Psi \\ \quad X_{ub}-\Lambda\xi(0)-\Upsilon\Psi \end{bmatrix}\\ &S^{l}_{Y}=\begin{bmatrix} -\Theta\\ \quad\Theta \end{bmatrix}\quad S^{r}_{Y} = \begin{bmatrix} -Y_{lb}+\Omega\xi(0)+\Phi\Psi \\ \quad Y_{ub}-\Omega\xi(0)-\Phi\Psi \end{bmatrix} \end{aligned} SXl=[−ΓΓ]SXr=[−Xlb+Λξ(0)+ΥΨXub−Λξ(0)−ΥΨ]SYl=[−ΘΘ]SYr=[−Ylb+Ωξ(0)+ΦΨYub−Ωξ(0)−ΦΨ]

则状态量与输出量的约束可表示为

SXlΔU⪯SXrSYlΔU⪯SYr\begin{aligned} &S^{l}_{X} \Delta U \preceq S^{r}_{X} \\ &S^{l}_{Y} \Delta U \preceq S^{r}_{Y} \end{aligned} SXlΔU⪯SXrSYlΔU⪯SYr

值得注意的是,因为本文中代价函数惩罚的是输出量,因此我们选择输出量约束作为QR问题的约束条件之一.

总结

综上,在每一个时间窗内MPC需求解的QR问题为:
minΔU(J∗)=minΔU(ΔUTHΔU+ΔUTK)s.t.SdulΔU⪯SdurSulΔU⪯SurSYlΔU⪯SYr\begin{aligned} &\mathop{min}\limits_{\Delta U}(J^*) =\mathop{min}\limits_{\Delta U}(\Delta U^TH\Delta U+ \Delta U^TK)\\ \\ &s.t. \qquad \qquad \qquad \begin{aligned} &S^{l}_{du} &&\Delta U &&\preceq S^{r}_{du}\\ &S^{l}_{u} &&\Delta U &&\preceq S^{r}_{u}\\ &S^{l}_{Y} &&\Delta U &&\preceq S^{r}_{Y} \end{aligned} \end{aligned} ΔUmin(J∗)=ΔUmin(ΔUTHΔU+ΔUTK)s.t.SdulSulSYlΔUΔUΔU⪯Sdur⪯Sur⪯SYr
其中
H=ΘtQcΘ+RcK=2ΘTQcNSdul=[−II]Sdur=[−ΔUminΔUmax]Sul=[−FF]Sur=[−Umin+Eδ(−1)Umax−Eδ(−1)]SYl=[−ΘΘ]SYr=[−Ylb+Ωξ(0)+ΦΨYub−Ωξ(0)−ΦΨ]\begin{aligned} &H = \Theta^tQ_c\Theta+R_c \\ &K = 2\Theta^TQ_cN \\ \\ &S^{l}_{du} = \begin{bmatrix} -I \\ \quad I \end{bmatrix} &&S^{r}_{du} = \begin{bmatrix} -\Delta U^{min} \\ \quad\Delta U^{max} \end{bmatrix} \\ \\ &S^{l}_{u} = \begin{bmatrix} -F\\ \quad F \end{bmatrix} &&S^{r}_{u} = \begin{bmatrix} -U^{min}+E\delta(-1)\\ \quad U^{max}-E\delta(-1) \end{bmatrix}\\\\ &S^{l}_{Y}=\begin{bmatrix} -\Theta\\ \quad\Theta \end{bmatrix} &&S^{r}_{Y} = \begin{bmatrix} -Y_{lb}+\Omega\xi(0)+\Phi\Psi \\ \quad Y_{ub}-\Omega\xi(0)-\Phi\Psi \end{bmatrix} \end{aligned} H=ΘtQcΘ+RcK=2ΘTQcNSdul=[−II]Sul=[−FF]SYl=[−ΘΘ]Sdur=[−ΔUminΔUmax]Sur=[−Umin+Eδ(−1)Umax−Eδ(−1)]SYr=[−Ylb+Ωξ(0)+ΦΨYub−Ωξ(0)−ΦΨ]

_{第二章完结撒花}

模型预测控制器(MPC)系列: 2.求解车辆横向控制中的MPC相关推荐

模型预测控制器(MPC)系列: 1.建立车辆横向动力学模型
勘误 Update 02/23/2021 之前的文章中有不严谨的地方,这里做一个勘误.错误就在下面描述坐标系的图中.<更正后的图已覆盖到坐标系小节下> 在这个图中,我指出ENU坐标下,车自 ...
基于arx模型的MPC预测控制器simulink建模与仿真实现
目录一.理论基础二.核心程序三.测试结果一.理论基础 MPC的优点模型预测控制善于处理多输入多输出系统对于MIMO系统,PID需要为每个子系统单独设计PID控制器,由于存在耦 ...
模型预测控制算法基础与车辆纵向控制仿真分析
模型预测控制算法基础与车辆纵向控制仿真分析第三章模型预测控制算法基础与控制仿真分析模型预测控制算法基础模型预测控制的基本思想就是利用已有的模型.系统当前的状态和未来的控制量去预测系统未来的输出 ...
如何理解MPC模型预测控制理论
MPC模型预测控制理论分析详解 MPC的基本原理与特点基本原理特点 1.基于模型 2.滚动优化 3.前馈-反馈结构从无约束模型预测控制出发理解推导过程 1.使用状态空间模型 2.模型的预测功能 ...
预测电流FCS-MPC模型预测电流控制
预测电流FCS-MPC模型预测电流控制级联(链式):5个单元下载地址[https://download.csdn.net/download/a_zxswer/20009112](https://d ...
滚动时域控制 matlab,在 Simulink 中设计神经网络预测控制器
在 Simulink 中设计神经网络预测控制器在 Deep Learning Toolbox™ 软件中实现的神经网络预测控制器使用非线性被控对象的神经网络模型来预测被控对象将来的性能.然后,控制器计 ...
异步电机模型预测磁链控制（MPFC）
导读:本期主要介绍异步电机模型预测磁链控制(MPFC),与模型预测转矩控制(MPTC)进行对比. 需要文中的MPFC仿真模型的学友,可关注微信公众号:浅谈电机控制,获取. 图1 异步电机模型预测转矩控 ...
2021全球与中国车辆线控转向系统市场现状及未来发展趋势
全球车辆线控转向系统(Vehicle Steer-by-wire System)核心企业主要分布在日本.北美.欧洲以及中国地区.其中头部企业主要有KYB Corporation.Nexteer Aut ...
ADRC，自抗扰控制器，扩张状态观测器，ESO，模型预测控制算法MPC
ADRC,自抗扰控制器,扩张状态观测器,ESO,模型预测控制算法MPC,自适应模型预测控制算法,时变模型预测控制算法,H无穷算法,混合灵敏度,鲁棒控制算法,四旋翼,直升机,控制算法设计,仿真模型,算法 ...

模型预测控制器(MPC)系列: 2.求解车辆横向控制中的MPC