随着市场竞争的日益激烈，环境保护意识的逐渐增强，如何更好地利用资源与能源以达到最优操作的目标，已经成为过程工业研究的重要内容. 传统的过程设计方法不能兼顾过程优化与运行控制的集成问题，造成过程设计与运行中各个环节相互脱节、衔接不够，成为制约过程工业发展的瓶颈问题.

过程工业中的实时优化（real time optimization，RTO）是指结合工艺知识和现场操作数据，通过快速、高效的优化计算，对操作运行中的生产装置参数进行调整，增强其对环境变化、市场变化的适应能力，保持生产装置始终处于高效、低耗、安全的最优工作状态^[1]. 过程装置一旦投入运行，将始终被一系列影响因素所干扰. 影响因素可以分为外部不确定性和内部不确定性2种^[2]：外部不确定性包括原料变化、公用资源的供应约束、市场需求变化和气候变化等；内部不确定性包括过程设备特性漂移以及来自过程中其他单元的影响等. 实时优化的目的是在满足所有约束条件的前提下，不断调整最佳工作点，以克服外部和内部不确定性因素，从而保证过程始终能够达到最佳经济效益^[3]. 现有的工业生产系统大多能保持平稳的运行状态. 因此，在过程工业领域，进行稳态实时优化是合理的简化，避免动态建模和求解的困难. 本研究所提出的优化方法也属于稳态实时优化方法.

复杂工业过程的实时优化问题常采用递阶控制方法. 这种方法既包含优化层，又包含控制层，其核心思想是每个决策单元同时响应子过程优化，所有的优化层通过相互迭代达到最优解. 但是在实际应用中，迭代的发生使工业过程不能达到最优解，同时还可能违反某些约束. Brdys等^[4]提出，从实际过程中提取变量的稳态信息，反馈至上一层的决策单元中，并利用它修正最优解，使其接近真实最优解. 另外，递阶控制方法的求解精度依赖于初始点的选取，Roberts^[5]提出将系统优化与参数估计2个问题分开处理，使用修正乘子交替处理系统优化与参数估计2个问题，直至收敛到最优解. 这一方法称为系统优化与参数估计集成方法（integrated system optimization and parameter estimation，ISOPE）. Tatjewski等^[6-7]通过对迭代优化算法中设定点控制问题的研究进一步改进ISOPE方法. 传统的实时优化问题是根据过程的数学模型，在相应的约束条件下，优化所提出的目标函数. 近年来，随着生产规模越来越大，生产过程越来越复杂，传统基于模型的实时优化方法由于很难获得精确可靠的过程模型而变得具有局限性^[8]. 测量值优化概念的引入，极大克服了传统优化方法的弊端，提高了预期的性能指标^[9-10]. 最优性必要条件跟踪法^[11]和自优化控制法^[12]就属于这类优化策略. Araújo等^[13]给出测量值优化策略更精确的定义，并将其成功应用于田纳西-伊斯曼挑战问题^[14]. Halvorsen等^[15]指出最优操作点附近的二阶近似值可以被用于寻找测量值的组合方式. 最优输出灵敏度矩阵也可以被用于寻找测量值的组合^[16]. Srinivasan等^[17-18]完成了其他相关优化策略的研究，并且考虑了不确定性因素的影响. 上述测量值组合的计算都是基于离线方法得到的，不需要使用显式的过程模型^[19].

本研究针对工业过程的最优操作问题，提出基于梯度信息的实时优化与控制集成策略，通过对被控变量的合理选取，将控制层与优化层联系起来，形成优化控制集成的级联结构；给出求解被控变量的解析方法，提高计算效率，使优化效果具有全局性；利用算例验证所提方法的优越性.

完整的过程系统可以纵向地分解为不同的结构层次，不同的结构层次运行在不同的时间尺度上^[20-21]，如图1所示. 层级模型的最上层是计划和调度层，主要针对供应链决策，如产品质量、安全性和经济效益等因素. 这一层级的时间尺度一般为周或月，取决于过程的类型和产品的规模. 优化层位于计划和调度层以下，主要负责实现计划和调度层的指令. 优化层包括企业级优化和实时优化2个部分. 企业级优化面向1个或若干个产品生产线的组合，处理全流程的工艺参数匹配，其运行周期为天或周；实时优化针对具体的生产流程，对流程运行中的设定值进行参数优化.

近年来，随着优化技术的发展，在线优化技术被越来越多地用于寻找最优设定值. 不过，在线优化技术包括稳态检测、数据估计和数据调和等复杂的前序步骤，同时在线优化问题本身也是大规模的非线性优化问题. 鉴于以上原因，本研究提出基于梯度信息的稳态实时优化方法，这种方法不需要使用显式的过程模型，而是通过收集测量值，计算过程的梯度信息. 梯度信息可以被用来更新最优设定值，一旦设定值更新完成，将被传递到控制层作为目标去实现，设定值的更新大约需要几分钟至几小时. 优化层以下是控制层，控制层可以产生操纵变量直接作用于实际过程装置. 控制层包含2个部分，分别为多变量协调控制和基础回路调节. 多变量协调控制主要用于实现实时优化给出的最优设定值，其运行周期为秒或分钟. 目前过程工业主流的多变量协调控制算法是模型预测控制（model predictive control，MPC）. MPC控制在处理各类操作变量和输出变量约束方面，有显著优势^[22-23]. 基础回路调节主要用于消除快速干扰，实现单回路或串级回路的基本平稳运行，其运行周期一般等于集散控制系统的采样周期. 基础回路调节的主流控制算法是比例积分微分控制（proportion integral derivative, PID）. 为了集成实时优化与控制，被控变量的选取尤为重要. 被控变量选取方法的主要思想为，在最优设定值更新期间，使所选被控变量对扰动信息不敏感. 在设定值更新期间只需要将被控变量调节到期望的设定值，就可以达到满意的控制效果. 所提出的实时优化与控制集成的级联结构如图2所示.

考虑如下带有约束的稳态实时优化问题：

(1) $\left.{\begin{array}{*{20}{c}} {\mathop {{\rm{min}}}\limits_{ v} }\;\;\;\;{{ Q}({ v,{ d}});}\\ {{\rm{s.t.}}}\;\;\;\; {{{ f}({ v,{ d}}) = {\bf 0},}\;\; {{ g}({ v,{ d}}) \leqslant {\bf 0},}\;\; {{ y} = {{ f}_y}({ v},{ d}).}} \end{array}}\right \} $

式中： ${ Q} ({ v},{ d})$ 为稳态实时优化问题的目标函数， ${ v} \in {{\bf{R}}^{{n_v}}}$ 为过程设定值， ${ d} \in {{\bf{R}}^{{n_d}}}$ 为扰动变量， ${ f}$ 为与系统模型相关的等式约束， ${ g} $ 为与操作相关的不等式约束（温度限制、流量限制、安全性指标等）， ${ y}$ 为测量值. 在求解此稳态实时优化问题时，发现某些约束是积极的，对应 ${{g} _i}( v, d) = 0.$ 为了达到最优操作，假设已经控制了所有的积极约束，并且这些积极约束不会发生变化. 控制积极约束会消耗一定的自由度，则剩余自由度对应的无约束简约空间优化问题为

(2) $\mathop {\min\; }\limits_{{w}} { P} ({{w}},{ d}).$

式中：w 为除去控制积极约束消耗的自由度后，剩余自由度对应的过程设定值向量.

式（2）中的等式约束（包括模型方程和积极约束）隐式地包含在目标函数 ${ P} $ 中，所以 ${ P} $ 不是 ${{w}}$ 和 ${ d}$ 的简单函数关系. 将目标函数 ${ P} $ 在当前设定点 ${{{w}}^k}$ 处，按二阶泰勒公式展开为

(3) ${ P} ({{{w}}^k} + {{\delta }},{ d}) = { P} ({{{w}}^k},{ d}) + {{G}}{({{{w}}^k})^{\rm{T}}}{{\delta }} + {{{\delta }}^{\rm{T}}}{{H}}({{{w}}^k}){{\delta }}/2.$

式中： ${{G}}({{{w}}^k})$ 为 ${ P} $ 在 ${{{w}}^k}$ 处的梯度向量， ${{H}}({{{w}}^k})$ 为 ${ P} $ 在 ${{{w}}^k}$ 处的海森矩阵， ${{\delta }} $为迭代步长. 当 ${{\delta }} $充分小时式(3)成立.

设新的迭代点为 ${{{w}}^{k + 1}} = {{{w}}^k} + {{\delta }}$，则 ${ P} $ 在 ${{{w}}^{k + 1}}$ 处取极小值的必要条件为

(4) ${\nabla _{{\delta }} }{ P} ({{{w}}^k}{{ + \delta }},d) = {{G}}({{{w}}^k}) + {{H}}({{{w}}^k}){{\delta }} = {{0}}.$

式中： ${\nabla _{{\delta}} }$ 为在新的迭代点处，P对 ${{\delta }}$ 的一阶微分.

假设海森矩阵 ${{H}}({{{w}}^{{k}}})$ 非奇异，可以得到

(5) ${{\delta = }} - {{{H}}^{{{ - 1}}}}{{(}}{{{w}}^k}{{)G(}}{{{w}}^k}{{)}}.$

上述方法可以在一步迭代之后达到一阶最优性必要条件. G 、 H都与梯度信息有关，因此将此方法命名为基于梯度信息的稳态实时优化方法. 在通常情况下，由于可行性问题，需要在增量 ${{\delta }}$ 前乘调节因数 $\alpha $， $\alpha \in [0,1.0]$，最终得到更新后的设定点为

(6) ${{{w}}^{k + 1}} = {{{w}}^k} + \alpha {{\delta }}.$

为了计算增量 ${{\delta }}$，需要分别计算梯度向量 ${{G}}$ 和海森矩阵逆矩阵 ${{{H}}^{ - 1}}$. 针对 ${{{H}}^{ - 1}}$ 的计算方法，已经得到广泛研究，其中最著名的方法是由Broyden、Fletcher、Goldfarb、Shannon分别但几乎同时提出的BFGS方法^[24].

提出利用最小二乘思想求解梯度向量 ${{G}}$ 的方法. 首先假设在第 $i$ 次迭代时，已经存在前 $m$ 次的迭代点 ${{{w}}^i},{{{w}}^{i - 1}}, \cdots ,{{{w}}^{i - m}}$，且 $\Delta {{{w}}^{ik}} = {{{w}}^i} - {{{w}}^{i - k}}(k = 1,$ $2, \cdots ,m)$线性无关. $\Delta {{ P} ^{ik}} = {{ P} ^i} - {{ P} ^{i - k}}$ 为当前目标函数值与第 $k$ 次迭代前的目标函数值的差值.

设

(7) ${{{S}}^i} = {[\Delta {{{w}}^{i1}},\Delta {{{w}}^{i2}}, \cdots, \Delta {{{w}}^{ik}}]^{\rm{T}}}.$

则在设定点 ${{{w}}^i}$ 处，梯度向量满足

(8) ${{{S}}^i} \times {{G}}\left( {{{{w}}^i}} \right) = {[\Delta {{ P}^{i1}},\Delta {{ P} ^{i2}}, \cdots, \Delta {{ P} ^{ik}}{\rm{]}}^{\rm{T}}}.$

令 ${{{S}}^i} = {{A}},$ ${{G}}\left( {{{{w}}^i}} \right) = {{X}}, $ ${[\Delta {{ P} ^{i1}},\Delta {{ P} ^{i2}}, \cdots, \Delta {{ P} ^{ik}}]^{\rm{T}}} = {{B}}.$ 则式（8）可以转换为如下矩阵方程问题：

(9) ${{AX}} - {{B}} = {\bf 0}.$

设

(10) ${{X}} = \left[ {\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1n}}} \\ {{x_{21}}}&{{x_{22}}}& \cdots &{{x_{2n}}} \\ \vdots & \vdots & & \vdots \\ {{x_{m1}}}&{{x_{m2}}}& \cdots &{{x_{mn}}} \end{array}} \right],$

则有

(11) $\begin{gathered} {{AX}} = \left[ {\begin{array}{*{20}{c}} \!\!\!\! {\Delta w_1^{i1}{x_{11}} + \Delta w_2^{i1}{x_{21}} + \cdots + \Delta w_m^{i1}{x_{m1}}} \!\!&\!\! {\Delta w_1^{i1}{x_{12}} + \Delta w_2^{i1}{x_{22}} + } \\ \!\!\! {\Delta w_1^{i2}{x_{11}} + \Delta w_2^{i2}{x_{21}} + \cdots + \Delta w_m^{i2}{x_{m1}}} \!\!&\!\! {\Delta w_1^{i2}{x_{12}} + \Delta w_2^{i2}{x_{22}} + } \\ \vdots \!\!&\qquad\quad\qquad\qquad\!\! \vdots \\ \!\!\!\! {\Delta w_1^{im}{x_{11}} + \Delta w_2^{im}{x_{21}} + \cdots + \Delta w_m^{im}{x_{m1}}} \!\!&\!\! {\Delta w_1^{im}{x_{12}} + \Delta w_2^{im}{x_{22}} + } \end{array}} \right. \left. {\begin{array}{*{20}{c}} \!\!\!\!\!\!\!\!\!\!\!{ \cdots + \Delta w_m^{i1}{x_{m2}}}& \cdots &{\Delta w_1^{i1}{x_{1n}} + \Delta w_2^{i1}{x_{2n}} + \cdots + \Delta w_m^{i1}{x_{mn}}} \!\!\!\! \\ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! { \cdots + \Delta w_m^{i2}{x_{m2}}}& \cdots &{\Delta w_1^{i2}{x_{1n}} + \Delta w_2^{i2}{x_{2n}} + \cdots + \Delta w_m^{i2}{x_{mn}}} \!\!\!\! \\ \!\!\!\! & & \qquad \vdots \!\!\!\! \\ \!\!\!\!\!\!\!\!\!\!\!\! { \cdots + \Delta w_m^{im}{x_{m2}}}& \cdots &{\Delta w_1^{im}{x_{1n}} + \Delta w_2^{im}{x_{2n}} + \cdots + \Delta w_m^{im}{x_{mn}}} \!\!\!\! \end{array}} \right]. \\ \end{gathered} $

又因

(12) ${{B}} = \left[ {\begin{array}{*{20}{c}} {\Delta {{P}} _1^{i1}}&{\Delta {{P}} _2^{i1}}& \cdots &{\Delta {{P}} _n^{i1}} \\ {\Delta {{P}} _1^{i2}}&{\Delta {{P}} _2^{i2}}& \cdots &{\Delta {{P}} _n^{i2}} \\ \vdots & \vdots & & \vdots \\ {\Delta {{P}} _1^{im}}&{\Delta {{P}} _2^{im}}& \cdots &{\Delta {{P}} _n^{im}} \end{array}} \right].$

根据矩阵 ${{AX}}$ 和 ${{B}}$ 对应元素相等可得到如下线性方程组：

(13) $ \left. {\begin{array}{*{20}{c}} {\Delta w_1^{i1}{x_{11}} + \Delta w_2^{i1}{x_{21}} + \cdots + \Delta w_m^{i1}{x_{m1}} = \Delta {{P}} _1^{i1}} ,\\ \vdots \\ {\Delta w_1^{i1}{x_{1n}} + \Delta w_2^{i1}{x_{2n}} + \cdots + \Delta w_m^{i1}{x_{mn}} = \Delta {{P}} _n^{i1}} ,\\ {\Delta w_1^{i2}{x_{11}} + \Delta w_2^{i2}{x_{21}} + \cdots + \Delta w_m^{i2}{x_{m1}} = \Delta {{P}} _1^{i2}} ,\\ \vdots \\ {\Delta w_1^{i2}{x_{1n}} + \Delta w_2^{i2}{x_{2n}} + \cdots + \Delta w_m^{i2}{x_{mn}} = \Delta {{P}} _n^{i2}} ,\\ \vdots \\ {\Delta w_1^{im}{x_{11}} + \Delta w_2^{im}{x_{21}} + \cdots + \Delta w_m^{im}{x_{m1}} = \Delta {{P}} _1^{im}} ,\\ \vdots \\ {\Delta w_1^{im}{x_{1n}} + \Delta w_2^{im}{x_{2n}} + \cdots + \Delta w_m^{im}{x_{mn}} = \Delta {{P}} _n^{im}}. \end{array}}\right\} $

根据最小二乘的思想，求解线性矩阵方程组式（13）的解，即求解如下最小二乘问题：

(14) $\mathop {\min }\limits_{\tilde{{ X}}} \;{\left\| {\tilde{{ A}}\tilde {{X}} - \tilde {{B}}} \right\|^2}.$

式中： $\tilde{{ A}}$ 为线性矩阵方程组式（13）的系数矩阵，

(15) $\tilde{{ B}}\!\! =\!\! {[\Delta {{P}} _1^{i1}, \cdots,\, \Delta {{P}} _n^{i1},\,\Delta {{P}} _1^{i2},\, \cdots ,\,\Delta {{P}} _n^{i2},\, \cdots ,\,\Delta {{P}} _1^{im},\, \cdots ,\,\Delta {{P}} _n^{im}]^{\rm{T}}},$

(16) $\tilde{{ X}} = {[{x_{11}},\; \cdots ,\;{x_{m1}},\;{x_{12}}, \;\cdots ,\;{x_{m2}}, \;\cdots ,\;{x_{1n}}, \;\cdots ,\;{x_{mn}}]^{\rm{T}}}.$

最小二乘问题式（14）的最优解 $\tilde{{ X}}$ 即为梯度向量在 $m$ 次迭代过程中完整的表达形式.

真实的工业过程具有不同程度的非线性，如果将已有的线性优化控制方法应用于非线性过程，控制效果只能在标称工作点附近得到保证，具有较大的局限性. 考虑如下过程的最优化操作问题：

(17) $ \left. \begin{gathered} \mathop {\min } \limits_{ u} { J} ({ u},{ d}'); \\ {\rm s.t.}\;\; {{ G} ({ u,{ d}}') \leqslant {\bf 0},} \;\; {{ y} = { f} ({ u,{ d}}').} \end{gathered} \right \} $

式中： ${ J} $ 为需要优化的目标函数， ${ u} \in {{\bf{R}}^{{n_u}}}$ 为操纵变量， ${ d}' \in {{\bf{R}}^{{n_{d'}}}}$ 为扰动变量， ${ G} \in {{{\bf{R}}}^{{n_g}}}$ 为约束条件， ${ y} = { f} ({ u},{ d}')$ 为非线性过程的输入输出模型. 非线性过程被控变量的选取目标如下：选取最优组合矩阵 ${{H}}$，使得当 ${{c}} = {{{H}}^{\rm{T}}}{{y}}\; ({{c}} \in {{\bf{R}}^{{n_c}}},{n_c} = {n_u},\;{{y}} $为测量值所组成的列向量)为被控变量时，系统能够在跟踪设定值的同时，自动接近最优工作点.

定义损失函数表达式为

(18) ${ L} = { J} ({ u},{ d}') - { J} ({{ u}^{{\rm{opt}}}},{ d}').$

式中： ${{{{ u}^{\rm opt}}}}$ 为操纵变量最优值.

将式（18）相对于 ${{c}}$ 按二阶泰勒公式展开为

(19) ${ L} = {({{c}} - {{{c}}^{{\rm{opt}}}})^{\rm{T}}}{{{J}}_{{\rm{cc}}}}({{c}} - {{{c}}^{{\rm{opt}}}})/2.$

式中： ${{{J}}_{{{{\rm{cc}}}}}}$ 为 ${ J} $ 相对于 ${{c}}$ 的海森矩阵， ${{c}}^{{\rm{opt}}}$为c的最优值.

关于被控变量选取的推论如下：任意被控变量 ${{c}} = {{{H}}^{\rm{T}}}{{y}}$ 的非0设定值 ${{{c}}_{\rm{s}}}$，可以通过以下方式变换为0: 将恒为1的测量变量 ${{{y}}_{{0}}}$ 加入到原测量变量 ${{y}}$ 中形成增广变量集 $\tilde{{ y}} = $ $ [{{{y}}_{{0}}};{{y}}]$，相应的增广最优组合矩阵为 $\widetilde{{ H}} = [ - {{c}}_{\rm{s}}^{\rm{T}};{{H}}]$.

上述推论的证明如下：

(20) $\tilde{{ c}} = {\tilde{{ H}}^{\rm{T}}}\tilde {{y}},$

(21) $\tilde{{ c}} = [{{ - }}{{{c}}_{{{\rm{s}}}}}{{,}}\;{{{{ H}}}^{{{\rm{T}}}}}] \times \left[ {\begin{array}{*{20}{c}} {{{{y}}_{{0}}}} \\ {{y}} \end{array}} \right] = - {{ c}_{{{\rm{s}}}}}{{{y}}_{{0}}} + {{{H}}^{\rm{T}}}{{y}} = {{ - }}{{{c}}_{{{\rm{s}}}}} + {{{c}}_{{{\rm{s}}}}} = {{{0}}}.$

在不引起歧义的前提下，为了表述方便，后文直接用 ${{c}}$、 ${{y}}$、 ${{H}}$ 代替 $\tilde{{ c}}$、 $\tilde {{y}}$、 $\tilde {{H}}.$ 根据上述推论，式（19）可以简化为

(22) ${ L} = {({{{c}}^{{\rm{opt}}}})^{\rm{T}}}{{{J}}_{{\rm{cc}}}}{{{c}}^{{\rm{opt}}}}/2.$

针对连续的扰动空间 ${ d}'$，将其离散化为 $N$ 种独立的扰动情况 ${{ d}'_i}\;(i = 1,2,\cdots,N)$，并定义最小化全局平均损失函数为

(23) $\mathop {\min }\limits_{{H}} \;{{ L} _{{\rm{av}}}} = \frac{1}{{2N}}\sum\limits_{i = 1}^N {({ y}_i^{{\rm{opt}}}} {)^{\rm{T}}}{{H(}}{{{J}}_{{{{\rm{cc}},{i}}}}}{{)}}{{{H}}^{\rm{T}}}{ y}_i^{{\rm{opt}}}.$

将式（23）写成矩阵的表达形式，并定义

(24) ${{Y}} = \left[ {\begin{array}{*{20}{c}} {{ y}_1^{{\rm{opt}}}} \\ \vdots \\ {{ y}_N^{{\rm{opt}}}} \end{array}} \right],\;{{W}} = \left[ {\begin{array}{*{20}{c}} {{{{J}}_{{{{\rm{cc}},}}1}}}& \ldots &{\bf 0} \\ \vdots & & \vdots \\ {\bf 0}& \cdots &{{{{J}}_{{{{\rm{cc}},}}N}}} \end{array}} \right].$

式中： ${{Y}}$ 的行向量为 $N$ 个扰动对应的最优测量变量 ${{y}}_i^{{\rm{opt}}}$， ${{W}}$ 为对角矩阵，对角线元素为 $N$ 个扰动对应的海森矩阵 ${{{H}}}{({{J}}_{{{{\rm{cc}},}}i})}$. 最小化全局平均损失矩阵表达形式可以改写为

(25) $\mathop {\min }\limits_{{H}}\; {{ L} _{{\rm{av}}}} = \frac{1}{{2N}}{{{Y}}^{\rm{T}}}{{HW}}{{{H}}^{\rm{T}}}{{Y}}.$

由于矩阵 ${{W}}$ 中每一项都是决策变量 ${{H}}$ 的复杂非线性函数，不容易直接求解，一般取矩阵 ${{W}}$ 为单位矩阵 ${{I}}$，即忽略 ${{W}}$ 对全局平均损失的影响. 则式（25）可以简化为

(26) $\mathop {\min }\limits_{{H}} \;({{{Y}}^{\rm{T}}}{{H}}{{{H}}^{\rm{T}}}{{Y}}).$

式（26）的最优解不一定存在，需要对决策变量 ${{H}}$ 施加一个约束条件^[25]，即 ${{{H}}^{\rm{T}}}{{H}}{{ = }}{{I}}.$ 则可将式（26）进一步简化为

(27) $\mathop {\min }\limits_{{H}} \;({{{Y}}^{\rm{T}}}{{Y}}).$

其最优解 ${{H}}$ 的构造形式如下：首先对 ${{{Y}}^{\rm{T}}}{{Y}}$ 进行特征值分解，然后取 ${{H}}$ 的列向量为最小的 ${n_c}$ 个特征值对应的特征向量.

综上所述，非线性过程被控变量的选取步骤如下. 1）将连续的扰动空间 ${ d}'$ 离散化，得到 $N$ 种扰动情况 ${{ d}'_i}(i = 1,2,\cdots,N)$；2）针对每一种扰动情况 ${{ d}'_i}$ 分别进行离线优化，得到最优的测量变量 ${ y}_i^{{\rm{opt}}}$，并构造矩阵 ${{Y}}$；3）对 ${{{Y}}^{\rm{T}}}{{Y}}$ 进行特征值分解，取 ${{H}}$ 的列向量为最小的 ${n_c}$ 个特征值对应的特征向量；4）由 ${{c}} = {{{H}}^{\rm{T}}}{{y}}$，确定非线性过程的被控变量.

综合第2、3节的内容可以得到非线性过程实时优化与控制集成策略，具体表述如下：

在每一采样时刻 $k$，

1）根据最小二乘思想，计算梯度向量 ${{G}}({{{w}}^k})$；

2）根据BFGS方法，计算海森矩阵的逆矩阵 ${{{H}}^{ - 1}}({{{w}}^k})$；

3）更新设定点 ${{{w}}^k}$，

$\begin{gathered} {{\delta }} = - {{{H}}^{ - 1}}({{{w}}^k}){{G}}({{{w}}^k}), \,\, {{{w}}^{k + 1}} = {{{w}}^k} + \alpha {{\delta }}. \\ \end{gathered} $

4）通过第3节方法选取非线性过程的被控变量，即 ${{c}} = {{{H}}^{\rm{T}}}{{y}}$；

5）通过相应的控制算法，将选取的被控变量 ${{c}}$ 调节到新的设定点处，即 ${{c}} = {{{w}}^{k + 1}}$，并产生控制增量 $\Delta {{u}}.$

6）更新操纵变量 ${{{u}}^k}$，

$\begin{gathered} {{{u}}^{k + 1}} = {{{u}}^k} + \Delta {{u}}, \,\, k = k + 1. \end{gathered} $

考虑如下单目标、单变量数值优化问题：

(28) $J= {(u - d)^2}/4 + 2.$

y 的表达式为

(29) $ {\begin{array}{*{20}{c}} {{y_1} = u + d,} \;\; {{y_2} = {u^2} - d.} \end{array}} $

其标称工作点为

(30) ${u^*} = {d^*} = 0,\;{{ y}^*} = [ 0,\; 0]^{\rm{T}}.$

$d$ 为连续的扰动空间，为 $[ - 2,\;2].$ 按照本研究选取被控变量的步骤，首先将连续的扰动空间离散化，取 $N = 5.$ 离散化后的扰动情况为 ${d_i} = - 2, - 1,$ 0, 1, 2. 通过离线优化得到 ${u^{{\rm{opt}}}} = d$，构造矩阵 ${{Y}}$ 为

(31) ${{Y}} = \left[ {\begin{array}{*{20}{c}} 1&{ - 4}&6 \\ 1&{ - 2}&2 \\ 1&0&0 \\ 1&2&0 \\ 1&4&2 \end{array}} \right].$

对 ${{{Y}}^{\rm{T}}}{{Y}}$ 进行特征值分解，取 ${{H}}$ 的列向量为最小的1个特征值对应的特征向量，可得

(32) ${{H}} = {[ {0.941\;9,}}\;\;\;\;{ - 0.155\;9,}\;\;\;\;{ - 0.297\;4} ]^{{\rm{T}}}.$

由 ${{c}} = {{{H}}^{\rm{T}}}{{y}}$ 确定非线性过程的被控变量为

(33) ${c_1} = 0.941\;9 - 0.155\;9{y_1} - 0.297\;4{y_2}.$

为了对比效果，使用文献[26]提出的基于零空间方法所选择的被控变量：

(34) ${c_2} = 0.447\;2{y_1} + 0.894\;4{y_2}.$

在扰动空间 $d \in [ - 2,2]$ 内，对比本研究所提出方法与文献[26]所提出方法的性能损失，如图3所示. 由图可以看出，文献[26]提出的方法只有在标称点 ${d^*} = 0$ 附近性能损失较小，在远离标称点的区域性能损失较大，且距离越远，性能损失越剧烈. 相比较而言，本研究所提出的方法在整个扰动空间内的性能损失均较小. 就全局平均损失而言，所提出的方法的全局平均损失 ${{L} _{{\rm{av1}}}} = 0.357\;7$，文献[26]提出的零空间方法的全局平均损失 ${{L} _{{\rm{av2}}}} = 2.689\;1$. 可以发现，该方法所选择的被控变量使全局平均损失得到了大幅度的降低，仅为同类算法得到的平均损失的 13.3%.

如图4所示为强制循环的蒸发器,稀释液体在高温蒸汽中蒸发,然后在分离器内分离为气液两相, 浓缩后的溶液一部分作为产品,另一部分与进料混合后重新进入换热器. 该蒸发过程的机理模型可参见文献[27]. 过程涉及到的变量物理含义和标称工作点如表1所示.

该蒸发过程的 ${{u}}$、 ${{d}}$、 ${{y}}$ 分别为

(35) $\left. \!\! \begin{gathered} {{u}} = {[{q_{m1}},{q_{m2}},{q_{m3}},{q_{{m200}}},{p_{100}}]^{\rm{T}}}, \\ {{d}} = {[{x_1},{\theta_1},{\theta_{200}}]^{\rm{T}}}, \\ \!\!\!\!\!{{y}} \!=\! {[{q_{m1}},{q_{m2}},{q_{m3}},{q_{m5}},{q_{{m100}}},{q_{{m200}}},{p_2},{\theta_2},{\theta_3},{\theta_{201}}]^{\rm{T}}}. \end{gathered} \right\}$

扰动变量 ${{d}}$ 的变化范围为

(36) $\left. \begin{gathered} 4.75{\text{% }} \leqslant {x_1} \leqslant 5.25\text{%} , \\ 32 \leqslant {\theta_1} \leqslant 48, \\ 20 \leqslant {\theta_{200}} \leqslant 30. \\ \end{gathered} \right\}$

该过程的目标函数为最小化操作成本 $J$，包括蒸汽、冷却水和泵的功率、原料成本和产品价值：

(37) $\begin{split} {J} = 600{q_{m100}} + 0.6{q_{m200}} + 1.009({q_{m2}} + {q_{m3}}) + \\ {0.2{q_{m1}}} - 4\;800{q_{m2}}. \qquad \qquad \end{split} $

该过程与产品质量、安全性、操作限制相关的约束条件如下：

(37) $\left. \begin{gathered} {x_2} \geqslant 35.5\text% , \,\, 40 \leqslant {p_2} \leqslant 80 , \,\, {p_{100}} \leqslant 400,\\ 0 \leqslant {q_{m1}} \leqslant 20 , \,\, 0 \leqslant {q_{m3}} \leqslant 100, \,\, 0 \leqslant {q_{m200}} \leqslant 400. \end{gathered} \right\}$

式（38）中有2个约束在整个扰动空间中都是积极约束，即 ${x_2} = 35.5\text{%} $ 和 ${p_{100}} = 400.$ 这2个积极约束消耗了2个自由度，另外分离器的液位 ${L_2}$ 需要稳定控制在其标称工作点处，也消耗了1个自由度. 由于该蒸发过程测量变量数目较多，不需要使用全部测量变量去构造被控变量. 考虑如下的测量变量子集（ ${{{y}}_1} \sim {{{y}}_4}$），来对比本研究提出的算法和文献[26]提出的算法的性能. 2种算法均在Windows 7的Matlab 2014a环境下运行，硬件配置为Intel Core i7 @2.6 GHz处理器，运行内存为16 GB.

(39) $\left. \begin{gathered} {{{y}}_1} = {[{q_{m3}},{q_{m200}}]^{\rm{T}}}, \\ {{{y}}_2} = {[{q_{m2}},{q_{m100}},{q_{m200}}]^{\rm{T}}}, \\ {{{y}}_3} = {[{q_{m2}},{q_{m3}},{q_{m200}},{\theta_{201}}]^{\rm{T}}}, \\ {{{y}}_4} = {[{q_{m1}},{q_{m2}},{q_{m100}},{q_{m200}},{P_2}]^{\rm{T}}}. \\ \end{gathered} \right\}$

为了对比效果，取 $N = 500$，对扰动空间进行随机采样，产生100组随机扰动，计算平均损失 ${{L} _{{\rm{av}}}}$、最大损失 ${{L} _{{\rm{max}}}}$ 和运算时间 $t$，结果如表2所示. 由表2可知，利用该方法大幅度提高了优化效果，4种测量变量所产生的平均损失分别降低了56.4%、30.2%、39.1%、31.1%，运算时间分别降低了55.1%、72.1%、77.4%、81.3%，同时最大损失得到了相应的降低. 除了上述测量子集外，该研究对其他子集进行了验证，得到了相似的结果.

连续搅拌釜式反应器（continuous stirred tank reactor，CSTR）包含可逆放热反应过程 $A \rightleftharpoons B$，反应原理如图5所示.

根据质量守恒定律和能量守恒定律，可以得到系统的微分方程描述形式：

(40) $\left. \begin{gathered} {{{\rm{d}}{c_{\rm{o}}}{{(A)}}}}/{{{\rm{d}}t}} = ({c_{\rm{i}}}{{(A)}} - {c_{\rm{o}}}{{(A)}})/ \;{\tau } - r, \\ {{{\rm{d}}{c_{\rm{o}}}{{(B)}}}}/{{{\rm{d}}t}} = ({c_{\rm{i}}}{{(B)}} - {c_{\rm{o}}}{{(B)}})/\;{\tau } + r, \\ {{{\rm{d}}{T_{\rm{o}}}}}/{{{\rm{d}}t}} = ({T_{\rm{i}}} - {T_{\rm{o}}})/\;{\tau } - 5r. \\ \end{gathered} \right\}$

式中： ${c_{\rm{o}}}{{(A)}}$、 ${c_{\rm{o}}}{{(B)}}$、 ${T_{\rm{o}}}$ 分别为出口处物料 ${{A}}$、 ${{B}}$ 的浓度和出口处温度； ${c_{\rm{i}}}{{(A)}}$、 ${c_{\rm{i}}}{{(B)}}$、 ${T_{{\rm{i}}}}$ 分别为入口处物料 ${{A}}$、 ${{B}}$ 的浓度和入口处温度； $\tau $ 为停留时间； $r$ 为反应速率. $r$ 的定义为

(41) $r = 5\;000{{\rm exp}\;{\left( \frac{ - 10\;000}{{1.987{T_{\rm{o}}}}} \right)} }{c_{\rm{o}}}{{(A)}} - {10^6}{{\rm exp}\;{\left(\frac{{ - 15\;000}}{{1.987{T_{\rm{o}}}}} \right)}}{c_{\rm{o}}}{{(B)}}.$

该放热反应过程的 ${{u}}$、 ${{d}}$、 ${{y}}$ 分别为

(42) $\left. \begin{gathered} {{{u}}} = \left[ {{T_{\rm{i}}}} \right], \\ {{{d}}} = {\left[ \!\!\!\!\!{\begin{array}{*{20}{c}} {{d_1}},\;\;{{d_2}} \end{array}} \!\!\!\!\! \right]^{\rm{T}}} = {\left[\!\!\!\!\! {\begin{array}{*{20}{c}} {{c_{\rm{i}}}{\rm{(A)}}},\;\;{{c_{\rm{i}}}{\rm{(B)}}} \end{array}} \!\!\!\!\!\right]^{\rm{T}}}, \\ {{y}} = {\left[\!\!\!\!\! {\begin{array}{*{20}{c}} {{y_1}},\;\;{{y_2}},\;\;{{y_3}} \end{array}}\!\!\!\!\! \right]^{\rm{T}}} = {\left[ \!\!\!\!\!{\begin{array}{*{20}{c}} {{c_{\rm{o}}}{\rm{(A)}}},\;\;{{c_{\rm{o}}}{\rm{(B)}}},\;\;{{T_{\rm{o}}}} \end{array}}\!\!\!\!\! \right]^{\rm{T}}}. \\ \end{gathered} \right\}$

扰动变量 ${{d}}$ 的变化范围为

(43) $ \begin{gathered} 0.5 \leqslant {c_{\rm{i}}}{{(A)}} \leqslant 1.5 , \;0 \leqslant {c_{\rm{i}}}{{(B)}} \leqslant 0.5 . \end{gathered} $

目标函数定义如下：

(44) ${J} = {(0.001\;657{T_{\rm{i}}})^2} - 2.009{c_{\rm{o}}}{{(B)}}.$

CSTR模型的标称工作点如表3所示. 扰动变量的变化轨迹如图6所示，在t=1 000 s时， ${c_{\rm{i}}}{{(A)}}$ 产生0.3 mol/L的阶跃变化；在t=1 500 s时， ${c_{\rm{i}}}{{(B)}}$ 产生0.2 mol/L的阶跃变化；在t=2 000 s时， ${c_{\rm{i}}}{{(A)}}$ 产生−0.3 mol/L的阶跃变化；在t=3 000 s时， ${c_{\rm{i}}}{{(A)}}$ 产生−0.3 mol/L的阶跃变化；在t=4 000 s时， ${c_{\rm{i}}}{{(A)}}$ 产生0.3 mol/L的阶跃变化.

为了对比控制效果，选择文献[28]提出的最优条件跟踪法. 文献[28]提出的方法得到的测量变量 ${c_{\rm{o}}}{{(A)}}$、 ${c_{\rm{o}}}{{(B)}}$、 ${T_{\rm{o}}}$ 的变化轨迹如图7所示，本研究提出的方法得到的测量变量 ${c_{\rm{o}}}{{(A)}}$、 ${c_{\rm{o}}}{{(B)}}$、 ${T_{\rm{o}}}$ 的变化轨迹如图8所示，本研究提出的方法与文献[28]提出的方法得到的操纵变量 ${T_{\rm{i}}}$ 的变化轨迹对比如图9所示. 通过对比可以发现，使用本研究提出的方法可以使测量变量轨迹和操纵变量轨迹较平滑，意味着系统的鲁棒性得到了一定程度的提升. 本研究提出的方法的另一个优势是，其不需要使用显式的过程模型，可以有效抑制模型失配的影响. 假设在t=1 500 s时发生模型失配的情况，比较本研究方法和文献[29]方法的目标函数值变化，如图10所示. 可以发现，基于梯度信息的稳态实时优化方法可以有效抑制模型失配对优化目标的影响.

（1）提出实时优化与控制集成的级联结构，对于优化层采用基于梯度信息的稳态实时优化方法，控制层中的多变量协调控制采用MPC控制算法，基础回路调节采用PID控制算法，通过对被控变量的合理选择，将优化层与控制层联系起来.

（2）提出基于梯度信息的稳态实时优化方法，通过对过程测量值的在线采集，估计过程的梯度信息，进而更新设定值. 由于不需要使用显式的过程模型，可以有效抑制模型失配对优化目标的影响.

（3）提出利用最小二乘思想求解梯度向量的方法，降低了计算成本，可以应用于大规模工业过程的稳态实时优化.

（4）提出非线性过程被控变量的选取方法，利用非线性模型计算平均损失，优化效果具有全局性. 为了快速求解非线性规划问题，对某些条件进行合理假设，从而获得次优解. 同时给出求解次优解被控变量的解析方法，由于只需要对构造出的矩阵进行特征值分解，不需要计算一阶或二阶灵敏度矩阵，提高了计算效率.

通过对数值算例、蒸发过程和放热反应过程的研究，验证了本研究提出的方法的有效性. 下一步研究的重点是考虑更加真实的过程对象，例如包含随机噪声的过程系统，或者过程的约束条件存在病态的情况. 另外，如何选择测量变量子集也是值得研究的方向.