模型预测控制快速求解算法

1、模型预测控制快速求解算法模型预测控制(ModelPredictiveControl,MPC)是一种基于在线计算的控制优化算法，能够统一处理带约束的多参数优化控制问题。当被控对象结构和环境相对复杂时，模型预测控制需选择较大的预测时域和控制时域，因此大大增加了在线求解的计算时间，同时降低了控制效果。从现有的算法来看，模型预测控制通常只适用于采样时间较大、动态过程变化较慢的系统中。因此，研究快速模型预测控制算法具有一定的理论意义和应用价值。虽然MPC方法为适应当今复杂的工业环境已经发展出各种智能预测控制方法，在工业领域中也得到了一定应用，但是算法的理论分析和实际应用之间仍然存在着一定差距，尤其在多输

2、入多输出系统、非线性特性及参数时变的系统和结果不确定的系统中。预测控制方法发展至今，仍然存在一些问题，具体如下：模型难以建立。模型是预测控制方法的基础，因此建立的模型越精确，预测控制效果越好。尽管模型辨识技术已经在预测控制方法的建模过程中得以应用，但是仍无法建立非常精确的系统模型。在线计算过程不够优化。预测控制方法的一大特征是在线优化，即根据系统当前状态、性能指标和约束条件进行在线计算得到当前状态的控制律。在在线优化过程中，当前的优化算法主要有线性规划、二次规划和非线性规划等。在线性系统中，预测控制的在线计算过程大多数采用二次规划方法进行求解，但若被控对象的输入输出个数较多或预测时域较大时，该

3、优化方法的在线计算效率也会无法满足系统快速性需求。而在非线性系统中，在线优化过程通常采用序列二次优化算法，但该方法的在线计算成本相对较高且不能完全保证系统稳定，因此也需要不断改进。误差问题。由于系统建模往往不够精确，且被控系统中往往存在各种干扰，预测控制方法的预测值和实际值之间一定会产生误差。虽然建模误差可以通过补偿进行校正，干扰误差可以通过反馈进行校正，但是当系统更复杂时，上述两种校正结合起来也无法将误差控制在一定范围内。模型预测控制区别于其它算法的最大特征是处理多变量多约束线性系统的能力，但随着被控对象的输入输出个数的增多，预测控制方法为保证控制输出的精确性，往往会选取较大的预测步长和控制

4、步长，但这样会大大增加在线优化过程的计算量，从而需要更多的计算时间。因此，预测控制方法只能适用于采样周期较大且动态变化过程较慢的系统中。为使预测控制方法能在更多场合中应用，快速模型预测控制算法成为了一个新的研究方向。国内外研究现状近年来人们对预测控制算法的不足有了越来越清晰的认识，为了将该算法应用到更多领域，越来越多的学者对其进行不断研究和改进。阅读近些年国内外核心期刊的文献可知，人们对预测控制方法能够在更大更复杂的系统中应用寄予了很高期望，同时也在其不足方面做了很多探究和尝试，发展出了多种智能预测控制方法。快速模型预测控制算法作为目前智能预测控制方法之一，其研究方向主要有以下几个方面：（1）

5、显式模型预测控制（ExplicitModelPredictiveControl,EMP。2002年AlbertoBemporad等学者1提出了显式模型预测控制，该方法在预测控制基础之上，在线性时不变系统优化求解过程中引入多参数二次规划理论，对系统的状态区域进行凸划分，根据最优控制问题计算得到状态分区上相应的控制律。EMPC将模型预测控制的在线计算过程转化为离线和在线计算相结合的过程，大大减少了算法的在线计算时间，弥补了MPC方法反复在线计算的不足，EMPC也在电力电子2、电机控制3等领域得到了很好应用。但随着被控对象问题规模（如输入、状态维数、约束等）的增大，EMPC算法在离线计算过程中所求的

6、状态分区数会呈指数倍增加，而状态分区数的增加不仅会导致存储状态分区所需的存储空间增加，还会导致EMPC算法在线查找最优解所需的计算时间增加，因此该算法很难适用于状态约束较多（状态变量往往不超过5）、预测步长较大的复杂系统。基于以上原因，许多学者提出用近似的显式模型预测控制方法4,5来代替精确的EMPC算法，即通过牺牲一定的控制精度来降低计算过程中的复杂度，从而简化整个求解过程。如文献6提出了显式模型预测控制的多尺度近似方法，通过引入分段线性插入法和自适应分层函数近似法，运用重心插值理论对EMPC离线计算出的状态空间进行网格划分，得到近似的状态分区和近似控制律。基于同样的思想，文献7提出了显式模

7、型预测控制多胞体近似方法，主要利用双描述法对最优控制问题进行近似处理，再通过重心插值得到近似控制律。文献8基于小波的多分辨率分析提出了近似EMPC通过二次插值和网格划分得到低复杂度且可保证系统稳定的近似控制律。上述文献提出的方法均能在误差允许的范围内保证系统的控制性能，在一定程度上解决了EMPC随着问题规模增大而带来的复杂度和存储空间增大的问题。(2)模型预测控制的简化算法MPC算法采用的是在线滚动优化的控制策略，但随着工业模型和环境越来越复杂，其在线计算量越来越大，所以限制了MPC算法在动态变化较快系统中的应用。为减少MPC在线优化求解的计算量，有学者考虑对参数进行优化，提出了预测控制的简化

8、算法，如将参数分块化的blocking技术9,其思想是将越远离当前时刻的控制输入越粗略计算，从而减少在线计算量。在此基础上，又有学者提出了移动blocking方法，其核心是限制系统优化变量个数同时增大系统的有效输入步长10。相较于blocking技术，文献11提出的简化方法主要把约束分为起作用约束集、不起作用约束集和不确定约束集，对不起作用的约束集进行忽略操作，对起作用约束集和不确定约束集进行优化计算，从而降低在线计算量，加快控制进程。(3)改进的在线优化算法模型预测控制的核心是采用反复的迭代优化进行在线求解，选取适当的在线优化方法可以提高在线计算速度。近年来，有学者试图对模型预测控制的标准形

9、式做适当变形或者近似处理，继而降低预测控制方法的在线计算量。如文献12提出了一种用扩展的牛顿拉夫逊(Newton-Raphson)方法来代替传统的二次规划方法，当问题规模增大时，不仅能够保证优化问题总是收敛，还可以有效解决MPC在线计算量过大的问题。文献13提出了将表存储和在线优化相结合的部分枚举法(PartialEnumeration,PE),对于规模较大的线性模型不仅有很好的控制效果，而且求解速度是传统MPC算法的5倍以上。文献14提出了基于降精度求解准则的序列二次规划法，主要通过牺牲一定的精确度来降低在线优化所需的迭代次数，从而提高在线求解速度。文献15提出了可用于预测控制在线优化过程的

10、有效内点法，该内点法不仅可以减少代码量，而且可以使在线计算时间提高2-5倍。理论基础凸集与凸函数首先给出凸优化理论中凸集、凸函数和仿射的定义。凸集：假设C为n维实数空间Rn中的集合，若C中任意两点之间的线段仍然在集合C中，即对于任意xi,x2C且对任意实数0,1,都满足：xi(1)X2C(3-1)则称集合C为凸集，反之则为非凸集。满足凸集的集合具有以下性质：(1)若集合CC2都是凸集则集合C1C2x|xy1N2MC1,y2C2也是凸集。(2)两个或多个凸集的交集仍为凸集。(3)空集也是凸集。凸函数：设fx是定义在非空凸集CRn上的函数，对于凸集C中任意的x1和x2以及任意实数0,1均满足f(x

11、1(1)x2)f(x1)(1)f(x2)(3-2)则称f(x)为凸函数，反之则为非凸函数(或称凹函数)。若对x1x2和任意实数0,1均满足f(x1(1)x2)f(x1)(1)f(x2)(3-3)则称f(x)为严格凸函数。同样，根据凸函数的定义可以推导出凸函数的几个性质，具体如下：(1)若fi(X),L,fm(x)是定义在凸集C上的凸函数,则1f1(x)Lmfm(x)也是凸集C上的凸函数，其中为任意实数。(2)若函数f(x)是定义在凸集C上的凸函数，则对于任意实数，函数f(x)也是凸集C上的凸函数。(3)若g(x)是单调递增的凸函数，f(x)也是凸函数，则复合函数gf(x)也是凸函数。(4)若f

12、(x)和g(x)都是定义在凸集C中非单调递增的凸函数，则f(x)g(x)也是定义在凸集C中的凸函数。仿射：对于集合C中任意两点xi和x2且任意实数R,若满足x1(1)x2c(3-4)则称集合C是仿射的。也就是说，通过集合CRn中任意两点的直线仍在集合C中。凸优化问题凸优化问题是指性能指标函数和不等式约束函数均为凸函数且由约束条件所得到的集合是凸集的最优化问题。由于凸优化问题的局部最优解就是全局最优解，因此对于一个实际问题，如果能表示成凸优化问题，则意味着该问题可以得到彻底解决，对实际应用有着十分重要的意义。凸优化问题有多种形式，其基本形式可表示为：minf(x)s.t.g,(x)0i1,L,m

13、(3-5)hi(x)0i1,L,p其中f(x)为性能指标函数，gi(x)0为不等式约束条件，hi(x)0为等式约束条件，Cxgi(x)0i1,L,m为可行域(也称可行集或约束集)。%(x)0i1,L,p如果上述问题中性能指标函数f(x)和不等式约束gi(x)都是仿射函数，则该(3-6)问题又称线性规划问题(LinearProgrammingLP),其形式可表示为：mincTxds.t.GxhAxb其中GRmn,ARpn,对于线性规划问题而言，通常将性能指标函数中常数项d忽略，因为其不影响最优解的集合。该问题的几何解释可由图表示，图中阴影多面体为可行域P,若性能指标函数中cTx为线性的，则其等位

14、曲线是和c正父的超平面(图中的虚线)，点x对即为最优点。图线性规划的几何解释如果性能指标函数f(x)是凸二次型且约束条件是仿射函数，则该问题又称二次规划问题(QuadraticProgramming,QP)。对于二次规划问题，其标准形式可表小为：.1TTmin-xQxqxr2(3-7)s.t.AxbGxh其中QSn,qRn,ARmn,bRm,GRpn,hRp。对于该问题，最优解的求解即是在多面体上极小化一个凸二次函数，其几何解释如图所示，该问题的多面体可行域P如图阴影所示，性能指标函数的等位线如图虚线所示，则点x*即为最优点图二次规划问题的几何解释对于凸优化问题的求解，理论上很容易得到其最优解

15、，因为其局部极小点就是全局极小点，如果求解的凸优化问题中性能指标函数是严格凸函数，则求得的极小点是唯一极小点。基于凸优化的快速模型预测控制原始对偶内点法求解QP问题对于标准二次规划问题（3-7）,本算法采用原始对偶内点法进行求解。相较于有效集等方法，原始对偶内点法无需热启动，且通过适当优化即可保证算法在5-25步内迭代便可精确求解。首先引入松弛变量sRp,将标准二次规划问题转化为：n1TTmin-xQxqx2s.t.Axb（3-8）Gxshs0对于式（3-8）,定义双变量yRm且满足等式约束条件，zRp满足不等式约束，求导后得到其最优化条件（KKT条件）：AxbGxshs0,z0（3-9）Qx

16、qGTzATy0ZS0,i1L,p其中ZiSi0为互补松弛条件模型预测控制问题模型预测控制算法需要依靠计算机实现求解，其处理的问题形式一般是离散控制问题而非连续问题。对任意线性时不变随机控制问题，其动态方程可表示为:(3-10)x(t1)Ax(t)Bu(t)(t)其中t0,1,L表示时间，x(t)Rn为状态变量，u(t)Rm为控制输入,(t)Rn为扰动项，ARnn和BRnm分别为状态和输入矩阵。预测控制方法的核心任务往往是求解系统的性能指标函数，在优化求解过程中，当该函数的值达到极大或极小时，表明该情况下的系统处于最优状态，而此时的控制输入即为系统的最优控制律，因此系统最优控制律的选择往往取决

17、于系统的性能指标函数。可以定义上述问题中的性能指标函数为：TxQxtTuRux(tT)TQfinaix(tT)(3-11)其中Qdiagq1,q2,L,qn为状态权重矩阵，Rdiagri,r2,L,rm为输入权重矩阵，Qfinal为状态终端权重矩阵，且均为对角阵；T为预测步长。由于实际被控系统的输入输出往往存在局限性，需考虑状态变量和控制变量的约束条件，该问题的模型预测控制形式可表示为：minJs.t.x(1)AxBuuminuumaxxminxxmaxt,L,tT1(3-12)其中x(t1),L,x(tT)为状态变量，u(t1),L,u(tT)为控制变量。对于该次规划问题，求解后的最优解可表

18、示为u*(t),L,u*(tT1)和x*(t),L,x*(tT1),预测控制方法仅取最优控制序列的第一项u*()作为控制输入。障碍内点法为更快地求解上式(3-12)中的最优控制问题，利用其特殊结构，采用改进后的原障碍内点法可以降低问题求解的复杂度，从而达到快速求解的目的。首先定义全局优化变量zu(t),x(t1),L,u(tT1),x(tT)对式子(3-12)进行重写，代入化简后可转化为：(3-13)minzTHzgTzst.CzbPzh00-0尺00n。0000UU。010-0a巾1A0077-000+9VH+-i*i+9-+,1*+r=000-1b0-000。U0。00-1000000nR

19、11000。000一月/00000xm0-打/。00.000-A-B00C-1hiEiax土«+p1喝»44V-nuib000000r"ikllll00000*AB.r矩P$H是由各权重矩阵Q、R和Qfinal组成的块对角矩阵，C是由状态矩阵A、B和单位矩阵I组成的稀疏矩阵，矩阵h由状态和输入的最大值和最小值组成。利用原障碍内点法对重组后的QP问题进行求解。原障碍内点法的实质是通过引入障碍参数，将含不等式约束的最优控制问题转化为等式约束的凸优化问题。因此，令障碍参数>0,替换(3-13)中的不等式约束后得到其近似问题为：一TT(3-14)minzHzgzzs

20、.t.CzblTk其中zloghiPiTz,Pi表示矩阵P的行。此时，式(3-13)已经转化为由确定目标和线性等式约束所组成的凸优化问题。同时由式（3-14）易知，障碍参数越接近于零，近似性能指标函数越接近于精确的性能指标函数，求得的最优解也越精确。对于式（3-14）,采用不可行初始点牛顿法进行求解。定义双变量vRTn且满足上述等式约束条件，求导后得到最优化条件：(3-15)Czb0JT_T2HzgPdCv0z其中diKt-pT表示矩阵P的行，PTd是z的梯度。此时，mPiz）优化问题（3-14）的求解可等价为对优化条件（3-15）的求解。对于当前优化条件，假设当前时刻全局优化变量z不可行（z

21、满足不等式约束，但不满足等式约束条件），则需找到一个方向z和v分别使得zz和vv满足（至少近似满足）最优化条件，即zzz*、vvv*。然后在优化条件中用zz代替z*、.*.一，、，,vv代替v,并利用泰勒公式的一阶近似，将上述问题（3-15）进行整理后可得关于梯度z和v的一组线性方程：(3-16)2HPTdiag（d）2PCTz2HzgPTdCTvC0vCzb其中PTdiag（d）2P为z的Hessian阵，同时将上述方程等式右边变量分别记为对偶残差rd2HzgPTdCTv,原残差rpCzb,残差向量r（z,v）（Q（x,v）,rp（x,v）。对式（3-16）进行计算后可得到梯度z和v的值，再

22、利用回溯直线搜索法得到步长s0,1,同时可保证更新后的优化点满足不等式Pzh0最后利用步长s更新原始和对偶变量z:zsz,v:vsv0该迭代修正过程不断重复进行直到残余向量的二范数小于设定的误差阈值。上述问题若严格可行，则需要有限步数的计算便能获得一定的原始可行性，而一旦e=0,那么在以后的迭代过程中残余向量会一直等于0,且变量z和v将趋于优化点。为加快求解，本算法采用固定的障碍参数和固定的最大牛顿迭代步Kmax来简化计算。在求解优化问题(3-16)时，原障碍内点法分别对一系列不同障碍参数X的优化问题进行依次求解，再分别比较不同障碍参数下的控制性能得出最优解，因此会花费较多的求解时间。大量仿真

23、表明只需使得性能指标函数达到最小就能得到很好的控制律，而非精确求解QP问题，因此对于优化问题(3-16),障碍参数越接近于零，近似性能指标函数越接近于精确的性能指标函数J,求得的最优解也越精确。固定障碍参数X能够减少预测控制问题次迭代求解所需的牛顿迭代次数(最大牛顿迭代步)，从而加快凸优化问题的求解。在标准的牛顿法中，只有当残差达到足够小或者迭代达到上限时则停止迭代寻优。若问题不是发散的，往往最大迭代步Kmax在3至IJ10之间就可以求出最优解。因此，固定的迭代步Kmax可能会使残差未能达到足够小，优化变量没有达到最优，但这并不影响最终控制效果，还能提高在线求解速度。1 BemporadA,M

24、orariM,DuaV,PistikopoulosEN.TheExplicitLinearQuadraticRegulatorforConstrainedSystemsJ.Automatica,2002,38(1):3-20.2 BeccutiAGPapafotiouG,RobertoFrasca,MarariM.ExplicitHybridModelPredictiveControloftheBoostDC-DCConverterA.IEEEPowerPowerElectronicsSpecialistsConference,Orlando,Florida,USA,2007:2503-250

25、9.3 MariethozS,MorariM.ExplicitModel-PredictiveControlofaPWMImverterWithanLCLFilterJ.IEEETransactionsonIndustrialElectronics,2009,56(2):389-399.4NguyenH,Olarus,HovdpatchyapproximationofexplicitmodelpredictivecontrolJ.IntermationalJourmalofControl,2012,85(12):1929-1941.5 BemporadA,FilippiC.Analgorith

26、mforapproximatemultiparametricconvexprogrammingJ.ComputationalOptimizationandApplications,2006,35(1):87-108.6 SummersS,JonesCN,LygerosJ,Morarimultiscaleapproximationschemeforexplicitmodelpredictivecontrolwithstability,feasibility,andperformanceguaranteesC.48thEEEConferenceonDecisionandControland28thChineseControlConferenceC,2009:6327-6332.7 JonesCN,MorariM.PolytopicApproximationofExplicitModelPredictiveControllersJ.IEEETransactionsonAutomaticControl,2010,55(11):2542-2553.8 SummersS,JonesCN,Morarimultiresolutionapproximationmethodfo