间歇过程的批间自优化控制

叶凌箭，宋执环，马修水

(1 浙江大学宁波理工学院, 浙江 宁波 315100；2 浙江大学控制科学与工程学系, 浙江 杭州 310027)

引 言

间歇过程的优化是一个动态优化问题，数学上求解相对于稳态过程难度更大。解析法如应用Pontryagin 极值原理最小化Hamiltonian 函数等，对于反应机理复杂、规模庞大的间歇过程几乎很难实现。应用数值法求解间歇过程的优化问题更为常见，如控制向量参数化法（control vector parameterization，CVP）等[1]。

对间歇生产过程进行操作优化时，如果仅在标称状态下进行一次离线优化，优化性能具有局限性。主要原因是工业过程在实际运行时受到各种不确定扰动的影响，引起了系统工作点的漂移。标称工作点的最优输入轨迹无法在真实工况下保持最优性，引起了目标函数损失。实时优化（real-time optimization，RTO）的任务即在不确定扰动影响下使过程系统重新运行于新工况下的最优点[2]。

根据RTO 策略和实现方法不同，Chachuat 等[3]将现有的RTO 方法主要分为3 大类：（1）参数估计-重优化两步法[4]；（2）修正适应法（modifier adaptation）[5]；（3）直接调节输入法。方法（1）利用过程的测量变量数据，首先基于参数化模型估计出不确定扰动，然后重新优化求解。方法（2）不对模型进行更新，而是对约束条件和目标函数进行合理修正，使得修正后的最优性条件（necessary conditions of optimality，NCO）和真实过程的NCO相匹配，从而找到最优输入。方法（3）以NCO 跟踪法[2,6]和自优化控制（self-optimizing control，SOC）[7]为代表，NCO 跟踪的主要特征是测量/估计出当前的NCO，设计控制律使得系统最终满足NCO。SOC则通过跟踪控制某些经过离线设计的被控变量，即使不在线调整设定值，系统也在反馈控制器作用下自动运行于最优点附近。以上优化策略中，大多需要在线测量/估计目标函数和测量变量的梯度信息，或者进行频繁的高负荷数学运算。相比而言，SOC的在线实现最为简便，控制结构中无须单独的优化层，而是通过常规控制作用实现RTO 效果。由于SOC 设计的被控变量（测量变量的组合）间接地获取了梯度信息，因此无须像有限差分法对系统进行激励以获取实验梯度，从而极大地加快了RTO 速度。

由于具备良好RTO 性能的同时保持了简洁的可实现性，SOC 自Skogestad[7]于2000年提出后就引起了相关学者的关注和研究[8-13]，大量平台上的成功应用不断报道。现有的SOC 主要应用于连续过程，在间歇过程上的应用较少。文献[14-16]基于SOC 思想对间歇过程进行了一些初步的工作，但具有如下两个方面的局限性：（1）满足最优性条件的Hamiltonian 函数经解析法推导，因此只适用于一些机理简单、规模微小的间歇过程；（2）要求实时测量一些关键的过程变量如物料浓度等，实际操作时通常很难满足这一点。如果过程变量只在批次运行结束后才被测量，只能通过批间方式实现RTO。目前，基于SOC 策略的间歇过程批间优化，仍未见 报道。

1 连续过程自优化控制

自优化控制的核心思想是对系统的控制结构进行干预，通过选择（设计）合理的被控变量实现RTO。扰动产生时，无须进行重新优化求解更新设定值，而是只在常规控制器下跟踪被控变量的恒定设定值，就能够使系统自动靠近最优点运行[7]。

在自优化控制方法发展的十多年中，文献[8-13]报道了不同方法求解最优被控变量，如使用单独的测量变量，或者它们的组合函数。本文采用的自优化控制法为作者提出的NCO 近似法[13]，并将其推广到间歇过程。

考虑如下连续过程的优化问题

式中，标量J 表示需要最小化的性能指标，u ∈ Rnu,d ∈Rnd分别是操纵变量和扰动变量。y ∈Rny为过程输出模型， G ∈Rng表示约束条件。假设uopt是最优输入，根据优化理论，在uopt点应该满足如下KKT 条件[17]

式中，Φ 为拉格朗日函数，μ ∈ Rng为拉格朗日乘子。满足上述KKT 条件的第一步是位于边界的积极约束Ga为0（对应非零拉格朗日乘子μi），Ga确定后可从式(2)约去μ 并得到[3]

式中，(·)+表示矩阵的Pseudo 伪逆矩阵。式(3)包含两部分积极约束Ga和简约梯度∇rJ，称为连续过程优化问题的NCO 条件。

NCO 近似法的基本想法是利用系统处于最优状态时必须满足NCO 条件（Ga和∇rJ）为0 这一特征，将其直接作为被控变量（设定值为0）实现理想的自优化控制。

假设NCO 条件（3）中的Ga部分可直接测量，则Ga直接作为被控变量进行卡边控制。而简约梯度∇rJ 往往依赖于不可测扰动d，因此不能在线测量。因此考虑建立∇rJ 的回归模型Z(y)，使用可测变量y估计∇rJ，进而设计被控变量c=Z(y)对其进行间接控制，其中Z 为可测变量y 的线性或非线性函数。

以线性最小二乘拟合为例，被控变量记为c = θ0+θTy ( tf)。将过程M 个不同操作点的测量变量和简约梯度排列为矩阵形式

式中，下角标i 表示第i 个操作点下的对应数据，γi表示梯度值。为使自优化控制效果在整个扰动空间都有效，M 个操作点应涵盖所有扰动工况。最小二乘拟合的结果为

除线性最小二乘回归外，多项式回归也可视拟合精度效果选择性采用，也可以使用更加复杂的人工神经网络、支持向量机回归等。对于具有多个决策变量的系统，应对每个简约梯度逐个拟合获取多个被控变量。文献[18]将自优化控制问题归结为一个软测量问题，这里不再赘述。

2 间歇过程的批间自优化控制

2.1 间歇过程优化

本文研究一类参数不确定型间歇过程，将自优化控制方法推广到间歇过程的批间实时优化。操作目标表示为如下最优化问题

式中，tf为批次运行时间，φ为标量目标函数，u(t)为连续的操纵输入轨迹，x 为nx维状态向量（初态x0）。J、y、d 的定义如前，其中y 测量值在批次终端获得，即y(tf)。F、H 和G 分别为nF维动态模型、nh维输出方程和ng维终端约束。

最优化问题式(6)是一个动态优化问题，可使用数值求解法如控制向量参数化法（CVP）求解得到最优输入轨迹uopt(t)。但由于系统受到不确定扰动d的影响，因此使用标称模型求解得到的最优输入轨迹并非适用于所有工况，需根据实际操作工况进行调整。对此，参数估计-重复优化法[4]采用了一种直接的优化思路，即首先根据当前测量变量y 估计出不确定扰动，然后更新模型并对其重新优化求解。缺点是此方法的计算负荷较大，另外，由于参数估计和优化两者之间存在冲突，需要对系统进行充 分激励才能准确估计出扰动，因此需要折中考虑 两者[4]。

2.2 输入轨迹参数化

动态优化问题式(1)可等价转化为如下静态优化问题[6]

操纵输入轨迹u(t)的参数化方案对求解式(7)具有重要影响，通常做法是将批次运行时间区间[0, tf]等间距分割成N 个时间间隔，在单个时间间隔内u(t)为常数或多项式曲线，其中视为常数时即CVP法。可以证明，当N→∞时，参数化后的近似曲线可逼近任意形状的曲线。如图1 所示，经CVP 法参数化后新的决策变量为u1,…,uN。

图1 基于CVP 法的输入轨迹参数化Fig.1 Input parameterization based on CVP (N=12)

提高CVP 法的算法精度可以通过增加N 来调节轨迹形状，也增加了优化问题中的决策变量个数。因此CVP 法优化性能的提高以增加计算复杂度为代价。这种方法的缺点在于，间歇过程中输入轨迹的变化趋势往往具有阶段性，不同时间段内的变化特性不同。对于较为平滑的部分，只需少量参数即可描述，而对较为“陡”的部分，才需较多的参数变量。对于CVP 法来说，由于整个时间区间是等间距划分的，因此需要大幅度增加N，以满足最“陡”的那部分输入轨迹曲线的需求，造成了计算资源的浪费。

本文考虑另一种更为灵活的输入轨迹参数化方法，下面结合图2 进行说明。图2 中均使用了非等间距的分隔方式对整个输入轨迹曲线进行划分，两者参数化后的决策变量均为u=[u1,u2,u3,t1,t2]T，图2(a)中单个间隔内u 均保持为常数，图2 (b)中后两段u 为斜坡形曲线。对于更为一般化的情况，也可以采用多项式曲线，n 阶多项式由n+1 个参数确定。注意，这里3 段曲线的切换时间t1和t2也作为自由决策变量考虑，这是因为间歇过程在受到外界不确定扰动的影响时，分阶段轨迹曲线往往呈现出一种在时间尺度上的拉伸或压缩的特性，此时使用时间参数作为自由变量更合理。

图2 输入轨迹参数化Fig.2 Input trajectory parameterization

图3 使用二段2 次多项式曲线对图1 所示的输入曲线进行了近似处理，此时输入轨迹被参数化为7 个变量u=[a0,a1,a2,b0,b1,b2,t1]T，相比原来N=12可大大降低求解复杂度。

图3 使用二段2 次多项式曲线进行输入轨迹参数化Fig.3 Input parameterization using two-stage second order polynomial functions

参数化法的性能好坏，很大程度上取决于参数化后的输入轨迹是否吻合最优输入轨迹的形态特征，以尽量少的参数准确描述。如果缺乏先验知识，仍然需要先借助如CVP 法求取最优轨迹，然后再对其进行合理的参数化。这是因为，在系统受到不确定扰动影响时，最优输入轨迹主要体现为这些参数的变化，而主要的形态特征往往不会发生较大改变，因此参数化法能大大降低优化求解难度。本质上来说，参数化就是对输入轨迹曲线施加额外的约束条件。一般来说，如果参数化后的决策变量过少，导致新的轨迹曲线和目标轨迹有较大差距，就可能引起一定程度的性能损失。此时需要设计人员进行折中考虑，在性能损失可以承受的范围内，尽量减少决策变量数目。间歇过程的特点是决策变量数量庞大，各个输入输出通道之间相互耦合，对实现反馈控制是不利的。输入轨迹参数化方案能够大幅度减少决策变量数目，极大地利于反馈控制的实施，因此在本文背景下具有重要意义。

2.3 间歇过程批间自优化

从前文分析可以看到，经过稳态转化和输入轨迹参数化后，最优化问题式(7)和稳态自优化控制问题具有相同的形式。因此分析式(7)的最优性条件，然后应用NCO 近似法，可以实现间歇过程的批间自优化控制。具体来说，是建立y(tf)和最优性条件之间的回归模型c=Z(y(tf))，设定值恒为0，然后对其进行反馈控制。

和连续过程自优化的不同点在于：（1）连续过程的操纵变量往往是确定不变的，而间歇过程的决策变量经过轨迹参数化得到，可以根据具体问题进行抉择；（2）连续过程的控制器工作于连续或者高频状态，而间歇过程的控制器工作频率为每批次进行，以最简单的PI 控制器为例，第k+1 批次的操纵输入uk+1的更新律为

式中，Kp和TI分别为比例系数和积分时间。

综上，本文提出的间歇过程批间自优化控制方法，具体实施步骤如下：

（1）使用传统的数值优化算法，在标称工作状态下对间歇过程进行优化求解，获得最优输入轨迹 ( )u t∗ ；

（2）观察最优输入轨迹的形态特征，将其参数化为若干少量的决策变量u。然后根据参数化方案产生的输入轨迹曲线u(t)，计算新的目标函数J(u(t))。衡量J(u(t))相对于 ( ( ))J u t∗ 的损失程度，如果损失可以接受，则进行下一步，否则重新调整参数化方案；

图4 本文方法的控制方框图Fig.4 Control diagram of proposed scheme

（3）基于间歇过程的动态模型，以计算机仿真的方式在整个操作空间内（包括扰动空间和决策变量空间）进行开环实验，得到测量变量y(tf)和目标函数J 相对于u 的梯度值∂J/∂u（可通过有限差分法得到）；

（4）以y(tf)为输入，∂J/∂u 为输出，建立回归模型c=Z(y(tf))；

（5）设计控制器，被控变量为c，设定值0，控制输入为u。

控制系统的在线运行如图4 所示，控制系统对第k 批次的终端输出值yk(tf)进行采样，计算NCO的估计值并作为反馈控制器的被控变量ck，以批间方式跟踪控制其设定值0。最后，使用控制输出uk产生轨迹曲线uk(t)使用于下一批次。

3 仿真实例

3.1 过程描述

研究一个间歇反应器[3]，原料A 和B 反应生成产物C，同时B 自身合成生成副产物D，如图5 所示。A 在反应初始时刻一次性投放，B 在反应过程中实时投放，流率为操纵变量F(t)

体系的模型方程为

图5 间歇反应器示意图Fig.5 Schematic diagram of batch reactor

式中，cX表示物料X 的浓度，其他各符号含义参见表1。

表1 反应过程参数Table 1 Parameters for reactor process

该过程的扰动变量为主反应常数k1，可能的波动范围为[0.02, 0.08]L·mol·min-1，过程运行时k1不可测。反应器过程的目标函数表示为

式中，第1 项为产物C 的产量，第2 项代表了操作成本，ω 为2500 L2·min-1·mol-1。

3.2 输入轨迹参数化

首先使用数值法求解标称工作点的最优输入轨迹，取CVP 法的N=15，使用Dynopt 工具箱[19]得到最优输入轨迹如图6 左上角所示，此时最优目标函数为= 0.2412。从图中可以观察到，最优输入轨迹形状的规律为：反应初期反应物B 的进料流率F 较大，以利于反应进行得到更多产物C；反应后期则减少B 的进料，这是因为反应生成C 需要一定时间，后续进料来不及在反应结束前产生效用，因此应该少量进料。

图6 CVP 法（N=15）和输入参数化法对比Fig.6 Comparison between CVP and input parameterization

根据上述规律，考虑如下一种简单的参数化方案：将整个反应时间段分为两个阶段，第一阶段[0, ts)实行恒定进料F(t)=FB，第二阶段[ts, tf]停止进料F(t)=0。这样整条输入轨迹只需两个参数[FB, ts]就完全确定。在此基础上对系统进行重新优化求解，最优目标函数Jopt=0.2395，FB和ts分别为3.98×10-4L·min-1和234.6 min。可以看到，使用参数化方案将决策变量从15 减少至2 个后，只引起了目标函数的少量损失，是非常值得的。两者的最优输入轨迹和状态变量对比如图6 所示，可以看到cA、V、cC、cD的变化曲线几乎是重叠的，cB略有差异，这是因为物质B 本身就是进料。

3.3 批间自优化控制

最优决策变量的分布情况是假设扰动变量已知，通过离线优化求解得到的。而系统在线运行时的真实扰动值不可测，因此只能通过间接的方式实现优化。

由于该间歇反应器不存在积极约束部分，因此只需控制两个梯度∂J/∂FB和∂J/∂ts为0即可保持系统最优性。在扰动变量k1∈[0.02, 0.08]的范围内，分别优化求解两个决策变量，得到它们的分布情况如图7 所示。可以看到，FB的最优值受扰动k1变化较为敏感，而ts几乎保持不变，变化范围为标称工作点上下±2 min 内波动。可以验证，ts在小范围内变化时对目标函数的影响不大，可以考虑取恒值ts=234.6 min 进一步简化操作。

最终，使用NCO 近似法对该间歇反应器进行批间自优化控制，只需得到可测变量y(tf)估计梯度∂J/∂FB的回归模型作为被控变量。开环状态下，随机变化决策变量和扰动变量，仿真获取100 组样本数据，每组样本数据包含可测变量 y(tf)和梯度∂J/∂FB。使用线性最小二乘拟合得到如下被控变量

图7 最优决策变量在整个扰动范围内的分布Fig.7 Distributions of optimal decision variables in whole disturbance range

使用PI 控制器，设定c 的参考值为0 进行批间自优化控制。安排实验如下：第1～8 批次k1=0.02，第9～18 批次k1=0.08，第19～24 批次k1=0.053，实验结果显示于图8 和表2。从图表中可以看到，方法具有良好的自优化控制效果。扰动发生后，由于被控变量c 在反馈控制器作用下调整Fs跟踪设定值0，使反应过程自动向最优工作点逼近，并且在若干批次内收敛到各自的最终工作点。3 种扰动情形中，k1= 0.02 时系统的目标函数损失最大为1.18%。这部分损失主要来自两方面：（1）从图7中可以验证，k1=0.02 时ts的最优值约为231.5 min，相对而言偏离固定值ts=234.6 min 最远，引起了一定程度损失；（2）对梯度∂J/∂FB拟合时引入了一定的拟合误差，如需进一步提高自优化控制效果，可以考虑使用多项式拟合提高精度。k1= 0.08 和k1=0.053 时系统的目标函数损失都很小（0.13%、0），可认为已近似达到了最优控制。

表2 批间自优化控制结果和真实最优值对比Table 2 Batch-to-batch self-optimizing control performances against true optimums

4 结 论

本文提出了一种基于NCO 近似法的间歇过程的批间自优化控制策略。使用批次运行终端的测量变量y(tf)，设计被控变量c=Z(y(tf))使系统能够在后续批次中跟踪设定值的同时，自动运行于（近似）最优点。应用本文方法时，一个重要步骤是对输入轨迹曲线进行准确的参数化。从间歇反应器研究的例子中可以看到，合理的参数化方案能够在维持系统优化性能的同时大大降低问题复杂度。间歇反应器的仿真结果表明，本文方法能有效实现间歇过程的批间自优化控制。

图8 3 种典型扰动情形下的控制效果Fig.8 Control performances under 3 typical disturbance scenarios

[1]Srinivasan B, Palanki S, Bonvin D.Dynamic optimization of batch processes (Ⅰ): Characterization of the nominal solution [J].Comput.Chem.Eng., 2003, 27 (1): 1-26.

[2]Srinivasan B, Bonvin D, Visser E, Palanki S.Dynamic optimization of batch processes (Ⅱ): Role of measurements in handling uncertainty [J].Comput.Chem.Eng., 2003, 27 (1): 27-44.

[3]Chachuat B, Srinivasan B, Bonvin D.Adaptation strategies for real-time optimization [J].Comput.Chem.Eng., 2009, 33 (10): 1557-1567.

[4]Chen T, Liu Y, Chen J.An integrated approach to active model adaptation and on-line dynamic optimisation of batch processes [J].J.Proc.Control, 2013, 23 (10): 1350-1359.

[5]Marchetti A, Chachuat B, Bonvin D.Modifier-adaptation methodology for real-time optimization [J].Ind.Eng.Chem.Res., 2009, 48 (13): 6022-6033.

[6]Francois G, Srinivasan B, Bonvin D.Use of measurements for enforcing the necessary conditions of optimality in the presence of constraints and uncertainty [J].J.Proc.Control, 2005, 15 (6): 701-712.

[7]Skogestad S.Plantwide control: the search for the self-optimizing control structure [J].J.Proc.Control, 2000, 10 (5): 487-507.

[8]Halvorsen I J, Skogestad S, Morud J C, Alstad V.Optimal selection of controlled variables [J].Ind.Eng.Chem.Res., 2003, 42 (14): 3273-3284.

[9]Alstad V, Skogestad S.Null space method for selecting optimal measurement combinations as controlled variables [J].Ind.Eng.Chem.Res., 2007, 46 (3): 846-853.

[10]Kariwala V.Optimal measurement combination for local self-optimizing control [J].Ind.Eng.Chem.Res., 2007, 46 (11): 3629-3634.

[11]Kariwala V, Cao Y, Janardhanan S.Local self-optimizing control with average loss minimization [J].Ind.Eng.Chem.Res., 2008, 47 (4): 1150-1158.

[12]Alstad V, Skogestad S, Hori E S.Optimal measurement combinations as controlled variables [J].J.Proc.Control, 2009, 19 (1): 138-148.

[13]Ye L, Cao Y, Li Y, Song Z.Approximating necessary conditions of optimality as controlled variables [J].Ind.Eng.Chem.Res., 2013, 52 (2): 798-808.

[14]Ye L, Kariwala V, Cao Y.Dynamic optimization for batch processes with uncertainties via approximating invariant//8th IEEE Conference on Industrial Electronics and Applications [C].Melbourne, Australia, 2013: 1786-1791.

[15]Hu W, Mao J, Xiao G, Kariwala V.Selection of self-optimizing controlled variables for dynamic processes//Proceedings of 8th International Symposium on the Advanced Control of Chemical Processes [C].Singapore, 2012: 774-779.

[16]Jaschke J, Fikar M, Skogestad S.Self-optimizing invariants in dynamic optimization//50th IEEE Conference on Decision and Control and European Control Conference [C].Orlando, FL, 2011: 7753-7758.

[17]Edgar T F, Himmelblau D M, Lasdon L S.Optimization of Chemical Processes [M].2nd Ed.New York: McGraw-Hill, 2001.

[18]Ye Lingjian (叶凌箭), Ma Xiushui (马修水).Optimal control strategy for chemical processes based on soft-sensoring technique [J].Journal of Zhejiang University: Engineering Science (浙江大学学报: 工学版), 2013, (7): 1253-1257, 1266.

[19]Cizniar M, Fikar M, Latifi M A.A Matlab package for dynamic optimisation of processes//Proceedings of the 7th International Scientific-Technical Conference-PROCESS CONTROL[C].Kouty nad Desnou, Czech Republic, 2006.