直接radau配置的多阶非线性模型预测聚合反应控制

孙京诰 ，占子赫

1.华东理工大学 信息科学与工程学院，上海 200237

2.华东理工大学 化工过程先进控制和优化技术教育部重点实验室，上海 200237

1 引言

半间歇聚合反应的温度控制一直是过程控制研究的热点问题[1-2]。非线性模型预测控制（Nonlinear Model Predictive Control，NMPC）具有显示处理系统约束和在不确定环境下优化控制的共性机理，使其成为继PID控制之后在控制工程中被广泛应用和研究的先进控制技术之一[3]。NMPC运用更加合理的目标函数，比如最小化批次时间，使得反应过程更加注重经济效益而不是过度地追求与目标温度的温差最小化[4]。然而在实际的化工过程中，往往存在实际参数与模型不匹配的情况，这使得NMPC的性能极大地下降。一类基于不确定预测模型鲁棒模型预测控制旨在弥补传统NMPC的不足，近年被众多学者广泛研究。最早有学者采用开环的min-max策略，即在通过最小化最差情况下的目标函数得到控制输入[5]。这种方法被证实不仅具有一定程度的保守性，而且由于基于开环控制的假设与真实的闭环系统不匹配，导致其并不具备鲁棒性。为了攻克这个问题，有学者提出了闭环min-max策略[6]，但是由于内在的计算复杂度，该算法在实际应用中的成功率有限。Tube-based mpc运用副控制器保证在不确定因素存在的情况下满足限制条件并保证控制的稳定性，但是该方法不能保证是最优控制[7]。近年来，多阶非线性模型预测控制（Multi-stage Nonlinear Model Predict Control，Multistage NMPC）在基于闭环鲁棒MPC和运用场景树对不确定量的建模的思想上产生[8-9]。在这个策略中，不确定量表示成场景树的形式，未来的控制变量作为辅助变量来优化，以此来降低该策略的保守性，Multi-stage NMPC已经被证实具备很好的解决传统NMPC面临的模型参数常量与实际不匹配问题的能力。

NMPC是一类基于在线优化的控制算法，其在线优化计算速度直接关系到能否成功应用。相比于线性模型预测控制，NMPC计算的最大困难在于要在有限的采样时间内计算一个非凸的非线性规划问题（Nonlinear Programming，NLP）[3]。求解该问题的计算量随决策变量的维数呈指数增长。故对于规模更大的Multi-stage NMPC，快速高效的求解方法是其能否成功应用的关键。现有的NMPC数值算法大致分为序贯法和联立法两种[10]，联立法也称为直接转化法，它通过同时离散状态变量和控制变量将问题转化为NLP。这样将DAE或ODE系统的求解和优化问题结合在一起，避免了中间解的寻找，并且具有更灵活的处理路径不等式的能力和处理不稳定和病态条件问题的能力[11]。有限元正交配置（Orthogonal Collocation on Finite Element，OCFE）是联立法常用的局部离散方法，具有精度高，计算量小的特点，在动态优化问题中有广泛的应用。此外，配置点的选择合理与否，不仅关系着解的精确性，也关系着运算的复杂度，正交配置中常用的配置点有radau点、legendre点和lobatto点[12]。

本文提出了多阶非线性模型预测半间歇聚合反应控制器设计方法，将NMPC滚动优化问题描述成最优控制问题，其目标函数综合考虑了温度偏移和批次反应时间。为了提高动态优化问题求解效率，采用基于有限元正交配置的联立法求解，通过实验选取合适的radau配置点，将状态变量和控制变量同时离散化，推导出基于multi-stage NMPC的可求解 NLP问题，然后借助IPOPT工具箱调用内点法求解[13]。

2 Multi-stage NMPC

Multi-stage NMPC最重要的理论基础是不确定量可以使用离散场景树的形式表示出来，如图1所示，在每一个节点的每一个分支都代表了一种未知的不确定因素和选择的控制变量对系统的影响，这些不确定因素可能是参数误差或者外部扰动。通过这个方程，在未来采样时刻将得到新的信息，因此最优化的控制变量可以当作抵消各种不确定参数的辅助变量，这种方法降低该策略的保守性。尽管在每个阶段场景树都有不同分支，但是这不表示在每个采样点不确定量是变化的。这只是表示如果不确定量如果在第K阶段是未知的，那么在第K+1阶段也是未知的。因此在下一个采样点，新的场景树将被生成以求出下一个时刻的控制量。

图1 Multi-stage NMPC的场景树示意图

由于场景树的离散设定，在场景树的设定下一个带有不确定参数的不确定非线性系统表述如式（1）：

上式表明k时刻的状态取决于前一个状态，相应的控制输入以及不确定参数简洁起见，所有的不确定参数和场景树都设置为统一的形式，在所有的节点上都有相同数目的分支。本文约定在一个节点处，一个不确定量用三个分支来表示，假设有m个不确定量，易知在每个节点需要有3m个分支，这样随着预测时域的增大，优化问题的规模也级数倍增长。为了避免场景树的规模过大而带来不必要的运算，设置鲁棒时域NR，在鲁棒时域后场景树将不再分支，如图1中，鲁棒时域设置为2。场景树中出现的(j,k)统一用符号S表示。每条从节点xk到叶子节点的路径称为一个场景。

在上述符号表述的基础上，基于场景树的最优化问题可如式（2）描述：

公式中，N表示场景的数目，Ji和wi分别表示每个场景Si的目标方程和其相应权重。X、U分别表示状态变量和控制变量的取值范围，xsoft和usoft分别表示对状态变量和控制变量的软化约束（Soft constraints）。式（6）称为非预期约束条件，表示在相同的分支处，设定控制量相同。对于一般的间歇反应模型，单个场景的目标方程有如下形式：

其中NP表示预测时域，Q和R为对角加权矩阵，目标方程的第一项表示状态变量跟踪设定值的惩罚项。第二项表示对控制变量变化量惩罚项，旨在防止控制变量剧烈变化。此外，限制条件包括对状态变量和控制变量取值范围的限制，同时，动态模型也作为最优化问题的约束条件。

3 将动态优化问题转化成NLP问题

使用联立法求解动态优化问题需要精确的离散方法以保证对连续的轨迹的逼近精度。为了使得到的NLP问题规模不至于过大而增加求解难度，一般情况下不采用全局的离散化方法。OCFE是一种局部的离散化方法，属于隐式Runge-Kutta法，可以提供高阶的精度和良好的稳定性。在一段采样区间内，定义有限元个数为nκ，有限元长度hκ。在每个有限元上通过Lagrange插值对状态变量，控制变量进行逼近：

其中(t)表示在场景树的(j,k)位置上，处于第κ个有限元区间内的状态变量值，(t)亦然。 τ∈[0,1]表示配置点。设nθ为插值的阶次，则一个有限元内的配置点 用 τ0,…,τnθ表示，本文选用Radau配置点，具体配置点信息将在第 5章给出。 lθ(τ)和 lˉθ(τ)分别表示状态变量和控制变量的Lagrange插值多项式，可用如下形式表示：

Lagrange多项式插值具有如下性质：

定义tk,κ,θ为第k阶段第κ个有限元第θ个插值点的时间，为第k阶段第κ个有限元的终端时间，该时间仅仅和有限元的个数nκ和插值阶次nθ相关。于是，在场景树中每个正交配置点上，满足

于是可以得到离散化后的NLP问题：

其中，式（12）通过将每个有限元微分得到，其中微分变量需要保持状态的连续性，式（13）保证了每个相邻有限元之间的连续性，式（14）保证了相邻控制时域内的连续性。状态变量约束和控制变量约束条件由式（15）给出。

4 半间歇聚合反应模型

为了验证提出的算法的可行性，本文对一种间歇性化工聚合反应进行研究，该模型来自于一个实际的化工间歇聚合反应模型，它的机理模型由BASF SE公司提供。

4.1 机理模型

该系统主要由一个反应釜，夹套以及一个外部热转换器（External Heat Exchanger，EHE）组成，夹套和EHE可以用来控制反应釜内的温度。其工艺图如图2。在一个批次的反应中，开始阶段将反应所需的引发剂，水等反应所需物质加入反应釜内，将釜温加热到设定温度后进入进料阶段，即开始持续向釜内加入单体，当加入的单体达到一定质量时停止进料，进入保持阶段，直到产生到规定质量的聚合物后反应停止。

图2 半间歇式聚合反应工艺图

根据物料守恒和能量守恒，其机理模型由主要由8个常微分方程构成，如式（17）～（24）所示：

该机理模型包括水mW，单体mA和聚合产物mP的物料守恒方程以及反应釜内温TR，管道温度TS，夹套温度TM，外部热转换器温度TEK和冷却剂在热转换器的出口温度TAWT的能量守恒方程。反应过程中有三个控制变量，分别是进料速率mF，外部热交换器入口温度夹套入口冷却剂温度TINM。

U表示反应釜内聚合物与单体的比例：

mges表示釜内物质的总质量：

kR1、kR2分别表示表示釜内的反应速率和外部热转换器的反应速率，如式（27）、（28）所示：

kK表示反应釜内的混合物总传热系数，由式（29）给出：

mA,R表示釜内的单体质量，由式（30）得出：

4.2 控制目标

在一个批次反应过程中，在满足安全条件限制和产物质量的前提下，要在尽可能短的时间内通过调整mF、TINAWT、TINM使反应安全稳定迅速地进行。温度的控制在聚合反应中占据了很重要的位置，它和聚合产物的质量息息相关。在本问题中，反应釜的温度应该被控制在设定温度Tset=(363.15±2.0)K范围内。在真实的反应过程中，出于对安全性和反应器材耐受能力的考虑，反应过程温度还需要满足安全限制。在本模型中，为了避免冷却失败的情况，反应釜最高温度被限制在382.15 K。在任意时刻反应釜的安全温度由式（31）给出：

此外，加入反应釜的单体的最大质量为30 000 kg。进料阶段投入单体，当添加完所需的所有的单体，反应进入保持阶段，如果聚合物的质量达到20 680 kg，则反应停止。若在给定的单体范围内没有生产出规定质量的聚合物则反应失败。在真实的系统中，通常模型的参数不能被准确地测量。在本模型中，有两个关键的参数：反应热焓ΔHR和反应速率k0，这两个参数通常会在其给定值的±30%内变化，且在一个批次的反应过程中值不变。模型的状态变量的初始值与取值范围和控制变量的取值范围分别如表1，表2所示。

表1 状态变量初始值和取值范围

表2 控制变量取值范围

5 仿真分析

所有的计算都运行在Windows环境台式电脑（Inter Core i7-6700K@4.00 GHz 16 GB DDR4 2 801 MHz）。Multi-stage NMPC所产生的NLP问题都采用最优化工具包IPOPT来求解，IPOPT工具包内部使用内点法对最优化问题进行求解。设定Multi-stage NMPC和NMPC的采样时间tstep=50 s，预测时域Np=10，设定Multi-stage NMPC的鲁棒时域NR=1。反应模型中两个不确定反应热焓ΔHR和反应速率k0在其给定值的±30%内变化，其真实值分表由，表示。由于有两个不确定量，在鲁棒时域内，场景树的每个节点都有9个分支，分支的不确定参数以最坏情况和标准情况的组合，即=(0.7,1.0,1.3)×ΔHR分别与=(0.7,1.0,1.3)×k0组合得到的9种可能。优化问题目标方程综合考虑了温度跟踪和反应的速率，设定式（6）中-1,100,0,0,0,0,0,0],R=diag[0.2,0.03,0.02]。此外，仿真过程中假定所需要的变量都可以被准确地测量。

5.1 正交配置点的选取

配置点τ的选择对求解的精度，速度以及解的存在与否都有很大的影响[15]。若采用等间隔方式选取配置点，当点过于密集时可能产生龙格（Runge）现象[14]。本文选择radau点作为配置点，使用radau配点的OCFE可以有效避免龙格现象，并获得平滑的逼近轨迹，且相比于legendre点，radau点可以提供更佳的稳定性。radau配置点取值如表3。

表3 radau配置点取值

设定有限元个数nκ=2，配置点阶次nθ分别取2、3、4进行实验仿真，由于模型参数的不确定性，在ΔHrealR=(0.7,1.0,1.3)×ΔHR分别与 kreal0=(0.7,1.0,1.3)×k0组合得到的9种可能情况下进行实验，分别得到批次反应时间和求解一次优化问题的平均时间，并以箱线图的形式画出，如图3、4。从图3中可以看出，随着nθ的增加，批次反应时间有减少的趋势，这符合目标方程的设定。可见nθ在一定范围内增加可以通过更好地逼近真实曲线而带来更好的控制效果。而从图4反应出，随着nθ的增加，在反应过程中，在每个采样时间求解最优化问题所需的平均运算时间也大幅度地增加，以中位数为例，nθ=2时，中位数为2.1 s，而在 nθ=3,4 时，这一数值分别达到了6.5 s和26.2 s。故虽然配点阶次增加可以带来控制效果的略微提升，但由于过度地增加了优化问题的计算时间使得不符合实际过程控制的实时性要求。

图3 radau配置点阶次为2、3、4时，一个批次反应所需时间

图4 radau配置点阶次为2、3、4时，求解一次优化问题平均时间

5.2 控制效果分析

根据5.1节实验结果，设定nκ=2，nθ=2，分别采用直接radau配置的Multi-stage NMPC（以下简称为RMSNMPC）和传统NMPC对模型进行仿真，模型参数设置如5.1节。在9种参数组合中有3种情况NMPC无法完成一个批次的反应，故不作图分析。其他仿真结果如图5～9，图中的红线表示温限。从图中可以看出，在任意一个情形下，RMSNMPC均可以较好地完成一个批次的反应，釜温不仅被控制在安全温度以内，也被控制在(363.15±2.0)K精度范围内。而NMPC在多数情形下，如图6、7、9，反应釜釜温都没有达到要求的控制精度要求。更为严重的是，在 ΔHrealR=1.0×ΔHR，kreal0=1.3×k0的情况下，如图7，反应进行到60 min时，处于进料阶段，釜温一度超过了安全温度，这会给实际生产过程带来严重的安全隐患。此外，虽然NMPC在标称情况和 ΔHrealR=0.7×ΔHR，kreal0=1.3×k0的情况下可以达到温度的精度和安全要求，但是进料量都有剧烈震荡的现象发生，这会给实际控制带来操作难度。相比之下，RMSNMPC在各个情形下都可以带来更加平滑的控制量。可见相比于传统的NMPC，RMSNMPC利用场景树模拟不确定量的优势，可以很好地处理具有不确定参数的模型。

图5 =1.0×ΔHR，=1.0×k0时，釜温、安全温度、进料量、夹套入口冷却剂温度曲线图

图6 =1.0×ΔHR，=1.3×k0时，釜温、安全温度、进料量、夹套入口冷却剂温度曲线图

图7 =1.3×ΔHR，=1.0×k0时，釜温、安全温度、进料量、夹套入口冷却剂温度曲线图

图8 =0.7×ΔHR，=1.3×k0时，釜温、安全温度、进料量、夹套入口冷却剂温度曲线图

图9 =1.3×ΔHR=0.7×k0时，釜温、安全温度、进料量、夹套入口冷却剂温度曲线图

6 结束语

针对半间歇聚合反应模型中存在不确定参数的情况，基于场景树模拟不确定量的Multi-stage NMPC方法，提出使用基于有限元正交配置的联立法对优化问题进行求解。通过大量实验选取最佳的radau配置点，形成有效的RMSNMPC方法。仿真结果表明在标称模型和具有不确定参数的模型下，所采用的算法不仅能使反应满足被控的温度精度要求和安全要求，还具备良好实时求解能力和鲁棒性。相比于传统的NMPC算法，更适用于复杂的半间歇聚合反应的控制。

[1]Wang J，Wang Y，Cao L，et al.Adaptive interative learning control based on unfalsified strategy for Chylla-Haase reactor[J].IEEE/CAA Journal of Automatica Sinica，2014，1（4）：347-360.

[2]Jang H，Lee J H，Biegler L T.A robust NMPC scheme for semi-batch polymerization reactors[J].IFAC-Papers on Line，2016，49（7）：37-42.

[3]何德峰，丁宝苍，于树友.非线性系统模型预测控制若干基本特点与主题回顾[J].控制理论与应用，2013，30（3）：273-287.

[4]Idris E A N，Engell S.Economics-based NMPC strategies for the operation and control of a continuous catalytic distillation process[J].Journal of Process Control，2012，22（10）：1832-1843.

[5]Campo P J，Morari M.Robust model predictive control[C]//American Control Conference，1987：1021-1026.

[6]Lee J H，Yu Z.Worst-case formulations of model predictive control for systems with bounded parameters[J].Automatica，1997，33（5）：763-781.

[7]Mayne D Q，Kerrigan E C，Van Wyk E J，et al.Tube based robust nonlinear model predictive control[J].International Journal of Robust and Nonlinear Control，2011，21（11）：1341-1353.

[8]Lucia S，Finkler T，Engell S.Multi-stage nonlinear model predictive control applied to a semi-batch polymerization reactor under uncertainty[J].Journal of Process Control，2013，23（9）：1306-1319.

[9]Lucia S，Engell S.Multi-stage and two-stage robust nonlinear model predictive control[C]//Nonlinear Model Predictive Control，2012，4（1）：181-186.

[10]陈杨，邵之江，钱积新，等.联立法中全局和局部正交配置算法[J].化工学报，2010（2）：384-391.

[11]黄森.动态优化问题的移动有限元算法研究[D].杭州：浙江大学，2012.

[12]Canuto C，Hussaini M Y，Quarteroni A M，et al.Spectral methods in fluid dynamics[M].New York：Springerverlay，2012.

[13]Biegler L T.An overview of simultaneous strategies for dynamic optimization[J].Chemical Engineering and Processing：Process Intensification，2007，46（11）：1043-1053.

[14]Ascher U M，Petzold L R.Computer methods for ordinary differential equations and differential-algebraic equations[M].Philadephia，PA，USA：SIAM，1998.

[15]陈杨.基于微分—代数混合方程机理模型的非线性预测控制[D].杭州：浙江大学，2011.