CN101763087A - 一种基于非线性共轭梯度法的工业过程动态优化***及方法 - Google Patents

Abstract

一种基于非线性共轭梯度法的工业过程动态优化***，包括与工业过程对象连接的现场智能仪表、DCS***和上位机；上位机包括：约束处理模块，用于处理优化过程中的控制变量边界约束；初始化处理模块，用于初始参数的设置；ODE求解模块，用于求解动态优化问题的常微分方程组；迭代优化模块，用于搜寻使目标函数J最优的决策向量w；收敛判断模块，用于判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，若成立，将当前的最优决策向量w^*、最优目标值J^*和迭代次数k保存。以及提供一种基于非线性共轭梯度法的工业过程动态优化方法。本发明能同时满足在线动态优化求解高效率、高精度。

Description

一种基于非线性共轭梯度法的工业过程动态优化***及方法

技术领域

本发明涉及最优控制领域，尤其是一种基于非线性共轭梯度法的工业过程动态优化***及方法。

背景技术

工业过程动态优化是过程模拟技术的核心，是过程优化设计、操作和控制的一个重要环节。采用动态优化方法来解决过程优化控制中的瓶颈问题和挖潜增效，已经越来越受到国内外学术界和工业界的重视。

随着工业过程对在线最优控制的需求的不断增加，改进动态优化算法性能，提高其在线应用的计算效率和准确性，变得越来越重要。

工业过程动态优化的难点在于需要在动态模型的基础上寻求目标泛函的最优值(动态模型往往是由一组复杂的大规模非线性微分方程组成的)。尽管目前的一些动态优化方法，如同步策略法、控制变量参数化法、遗传算法、粒子群算法、模拟退火法等，已经能够找到工业过程动态优化问题的最优解，但是往往出现计算不准确或收敛缓慢的问题，很难既保证所得最优控制结果的较好的准确性又满足动态优化求解过程的高效性。

发明内容

为了克服现有的工业过程动态优化***和方法很难同时满足在线动态优化求解准确性和高效性要求的不足，本发明提供一种能同时满足在线动态优化求解高效率、高精度的基于非线性共轭梯度法的工业过程动态优化***及方法。

本发明解决其技术问题所采用的技术方案是：

一种基于非线性共轭梯度法的工业过程动态优化***，包括与工业过程对象连接的现场智能仪表、DCS***和上位机；所述的DCS***包括操作站和数据库；所述现场智能仪表与DCS***连接，所述DCS***与上位机连接，其特征在于：所述的上位机包括：

约束处理模块，用于处理优化过程中的控制变量边界约束，采用以下转换方程：

u(t)＝0.5(u_max-u_min)×{sin[w(t)]+1}+u_min    (1)

将带有边界约束u_min≤u(t)≤u_max的控制变量u(t)转化为不受边界约束的变量w(t)，其中下标min、max分别表示最小值和最大值，u_min、u_max分别对应控制变量的下界和上界，并将w(t)作为优化问题的决策变量进行求解；

初始化处理模块，用于初始参数的设置，决策变量w(t)的离散化和初始赋值，具体步骤如下：

[3.1]设置判断迭代优化是否终止的收敛精度值ξ(当优化目标值迭代误差小于ξ时，停止迭代)；取迭代次数k初值为0；

[3.2]将时间域[t₀，tf]平均分成n小段：[t₀，t₁]，[t₁，t₂]，...，[t_n-1，t_n]；每个时间段的长度为h＝(tf-t₀)/n；

[3.3]对决策变量w(t)在[3.2]所述时间分段上进行离散化，将w(t)替换为由n个分段常值组成的向量w(即决策向量)，并选取任意常数向量作为决策向量的初始值w⁰；

[3.4]设置一维搜索的初始步长α，并设定一维搜索的相关参数α1，α2，p，q；

ODE求解模块，用于求解动态优化问题的常微分方程组，为迭代优化过程的梯度计算和收敛条件判断做准备，采取以下步骤来完成：

[4.1]数值求解状态方程：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),w\left(t\right),t]$ 边界条件：x(t₀)＝x0；(2)

其中f表示微分函数向量，x(t)为由m个状态变量组成的向量，x₀为状态向量x在初始时刻t₀的值，采用龙格库塔法由初始值x(t₀)正向积分求出；

[4.2]数值求解协态方程：

$\frac{\mathrm{d\λ}\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂x}-\mathrm{\λ}\left(t\right)\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂x},$ 含边界条件：

其中φ、ψ分别是目标函数 $J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}$ 中的非积分项和定积分项，λ(t)为由m个协态变量组成的向量，λ(tf)为协态向量λ在终端时刻tf的值，采用龙格库塔法由终端值λ(tf)通过反向积分求出；

[4.3]由所得的状态向量和协态向量计算出目标函数值：

$J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}---\left(4\right)$

迭代优化模块，用于搜寻使目标函数J最优的决策向量w，采取以下步骤来完成，上标k均表示迭代次数，初始赋值为零，获取ODE求解模块输出的状态向量、协态向量及目标函数值，即当前目标值J^k，进行如下操作：

[6.1]计算当前迭代的梯度g^k，其中，x(t)和λ(t)分别为ODE求解模块得出的当前状态向量和协态向量，上标T表示转置：

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^k}---\left(5\right)$

[6.2]计算当前迭代的搜索方向d^k，d^k-1表示第k-1次迭代的搜索方向，β^k是中间参数：

$d^k=\left\{\begin{array}{cc}{-g}^k,& k=1;\\ {-g}^k+{\mathrm{\β}}^k{d}^{k-1},& k\≥2.\end{array}\right.$ 其中， ${\mathrm{\β}}^k=\frac{{\left(g^k\right)}^T(g^k-g^{k-1})}{{\left|\right|g^{k-1}\left|\right|}^2}---\left(6\right)$

[6.3]进行一维搜索，从当前的迭代点w^k出发，沿方向d^k作一维搜索，计算出最佳步长α^*，满足

$H(w^k+{\mathrm{\α}}^{*}\·d^k)=\underset{\mathrm{\α}\≥0}{\min }H(w^k+\mathrm{\α}\·d^k)---\left(7\right)$

其中H表示哈密顿函数，由下式计算出，

表示在[0，+∞)中寻找使H达到最小值的步长α；

H＝ψ[x(t)，w(t)，t]+λ(t)^Tf[x(t)，w(t)，t](8)

[6.4]取当前迭代的搜索步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{B\·\left|\right|d^k\left|\right|},{\mathrm{\α}}^{*}),$ B取区间[5，8]内的整数值；

[6.5]计算下一个迭代点w^k+1＝w^k+α^k·d^k，增加迭代次数：k＝k+1，并把新的迭代点传给ODE求解模块以计算新的状态向量、协态向量及目标函数值J^k+1；

收敛判断模块，用于判定是否终止迭代计算，即判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，即收敛条件为：

|J^k-J^k+1|≤ξ    (14)

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；

若上式(14)成立，则将当前的最优决策向量w^*、最优目标值J^*以及迭代次数k保存到结果输出模块；若上式(14)不成立，则保存当前目标值J^k，然后转入迭代优化模块进行下一次迭代求解。

作为优选的的一种方案，所述上位机还包括：一维搜索模块，用于寻找使哈密顿函数H达到最小值的步长α^k，返回给迭代优化模块；依据强Wolfe-Powell准则进行一维搜索，按照以下步骤实现：

6.3.1)选取初始数据：在搜索区间[0，+∞)中选取步长的初始值α＞0，取α1＝0，α2＝∞，p∈(0，0.5)，q∈(p，1)；计算步长α^k为零时的哈密顿函数H₁＝H(w^k)，及其导数H₁’＝g^kd^k；

6.3.2)计算哈密顿函数H＝H(w^k+α^kd^k)，检验下式是否成立：

H≤H₁+pα^k(g^k)^Td_k    (9)

若成立，则转步骤6.3.4)；否则，利用二次插值公式计算试探点：

${\mathrm{\α}}_S=\mathrm{\α}1+\frac{1}{2}\·\frac{\mathrm{\α}-\mathrm{\α}1}{1+\frac{H_1-H}{(\mathrm{\α}-\mathrm{\α}1)\·H_1^{\′}}}---\left(10\right)$

6.3.3)令α2＝α，α＝α_s，转步骤6.3.2)；

6.3.4)计算哈密顿函数导数H’＝g^k+1d^k以及梯度向量

$g^{k+1}=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^{k+1}}---\left(11\right)$

6.3.5)检验下式是否成立：

|(g^k+1)^Td^k|≤-q·(g^k)^Td_k    (12)

若成立，则停止一维搜索，令α^k＝α，并输出α^k给迭代优化模块；否则，利用二次插值公式计算试探点α_s：

${\mathrm{\α}}_S=\mathrm{\α}+\frac{1}{2}\·\frac{(\mathrm{\α}-\mathrm{\α}1)H^{\′}}{(H_1^{\′}-H^{\′})}---\left(13\right)$

6.3.6)令α1＝α，H₁＝H，H₁’＝H’，α＝α_s，转步骤6.3.2)。

进一步，所述的上位机还包括：信号采集模块，用于设定采样时间，采集现场智能仪表的信号。

再进一步，所述的上位机还包括：结果输出模块，用于将迭代优化模块所获得的最优决策轨线为w^*(t)通过式(1)转化为最优控制轨线u^*(t)，然后将此最优控制信号传输给DCS***，并在DCS***中显示所得到的优化结果信息。

一种用所述的基于非线性共轭梯度法的工业过程动态优化***实现的动态优化方法，所述的动态优化方法包括以下步骤：

1)在DCS***中指定动态优化的状态变量和控制变量，根据实际生产环境的条件和操作限制的条件设定控制变量的上下边界u_max、u_min和DCS的采样周期，并将DCS数据库中相应各变量的历史数据，控制变量上下边界值u_max、u_min传送给上位机；

2)通过三角函数代换，对受边界约束的控制变量u(t)进行无约束转换，将其替换为另一无约束变量w(t)的函数表达式，即：

u(t)＝0.5(u_max-u_min)×{sin[w(t)]+1}+u_min    (1)

通过上式，控制变量u(t)∈[u_min，u_max]转化为无约束变量w(t)∈R的三角函数表达式；然后，将w(t)作为决策变量进行优化求解，最终求得的w(t)通过代人(1)式，即得到相应的u(t)；

3)对模块初始参数进行设置，并对DCS***输入的数据进行初始化处理，按照以下步骤完成：

[3.1]设置判断迭代优化是否终止的收敛精度值ξ(当优化目标值迭代误差小于ξ时，停止迭代)；取迭代次数k初值为0；

[3.2]将时间域[t₀，tf]平均分成n小段：[t₀，t₁]，[t₁，t₂]，...，[t_n-1，t_n]；每个时间段的长度为h＝(tf-t₀)/n；

[3.3]对决策变量w(t)在[3.2]所述时间分段上进行离散化，将w(t)替换为由n个分段常值组成的向量w，并选取任意常数向量作为决策向量的初始值w⁰；

[3.4]设置一维搜索的初始步长α，并设定一维搜索的相关参数α1，α2，p，q；

4)当前迭代的决策变量w＝w^k，迭代次数k为0时，w＝w⁰，采用数值积分方法求出当前状态变量、协态变量，并得出相应的当前目标值J^k。实现步骤如下：

[4.1]数值求解状态方程：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),w\left(t\right),t]$ 含边界条件：x(t₀)＝x0；(2)

其中f表示微分函数向量，x(t)为由m个状态变量组成的向量，x₀为状态向量x在初始时刻t₀的值，采用龙格库塔法由初始值x(t₀)正向积分求出；

[4.2]数值求解协态方程：

$\frac{\mathrm{d\λ}\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂x}-\mathrm{\λ}\left(t\right)\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂x}$ 含边界条件：

其中φ、ψ分别是目标函数 $J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}$ 中的非积分项和定积分项，λ(t)为由m个协态变量组成的向量，λ(tf)为协态向量λ在终端时刻tf的值，采用龙格库塔法由终端值λ(tf)通过反向积分求出；

[4.3]由[4.1]、[4.2]所得的状态向量和协态向量计算出当前目标函数值：

$J^k=J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}---\left(4\right)$

5)判断收敛条件|J^k-J^k+1|≤ξ是否满足，其中J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值，若不满足，则保存当前目标值J^k，再转入下一步骤进行下一次迭代优化；若满足，则终止迭代计算，当前决策向量就是最优决策向量w^*，当前目标值就是最优目标值J^*，保存w^*、J^*及迭代次数k到结果输出模块；

6)由状态向量和协态变量值计算出进一步迭代优化的搜索方向和步长，求解使目标函数J更接近最优的决策向量w，然后返回步骤4)求解ODE；实施迭代优化的步骤如下，上标k均表示迭代次数，初始赋值为零：

[6.1]计算当前迭代的梯度g^k，其中，x(t)和λ(t)分别为步骤4)得出的当前状态向量和协态向量，上标T表示转置：

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^k}---\left(5\right)$

[6.2]计算当前迭代的搜索方向d^k，d^k-1表示第k-1次迭代的搜索方向，β^k是中间参数：

$d^k=\left\{\begin{array}{cc}{-g}^k,& k=1;\\ {-g}^k+{\mathrm{\β}}^k{d}^{k-1},& k\≥2.\end{array}\right.$ 其中， ${\mathrm{\β}}^k=\frac{{\left(g^k\right)}^T(g^k-g^{k-1})}{{\left|\right|g^{k-1}\left|\right|}^2}---\left(6\right)$

[6.3]进行一维搜索，从当前的迭代点w^k出发，沿方向d^k作一维搜索，计算出最佳步长α^*，满足

$H(w^k+{\mathrm{\α}}^{*}\·d^k)=\underset{\mathrm{\α}\≥0}{\min }H(w^k+\mathrm{\α}\·d^k)---\left(7\right)$

其中H表示哈密顿函数，由下式计算出，

表示在[0，+∞)中寻找使H达到最小值的步长α；

H＝ψ[x(t)，w(t)，t]+λ(t)^Tf[x(t)，w(t)，t]    (8)

[6.4]取当前迭代的搜索步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{B\·\left|\right|d^k\left|\right|},{\mathrm{\α}}^{*}),$ B取区间[5，8]内的整数值；

[6.5]计算下一个迭代点w^k+1＝w^k+α^k·d^k，增加迭代次数：k＝k+1，并把新的迭代点传给步骤4)以计算新的状态向量、协态向量及目标函数值J^k+1；

7)判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，即收敛条件为：

|J^k-J^k+1|≤ξ    (14)

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；

若上式(14)成立，则将当前的最优决策向量w^*、最优目标值J^*以及迭代次数k保存到结果输出模块；若上式(14)不成立，则保存当前目标值J^k，然后转入迭代优化模块进行下一次迭代求解。

作为优选的一种方案：上述步骤[6.3]是依据强Wolfe-Powell准则进行一维搜索，按照以下步骤实现：

6.3.1)选取初始数据：在搜索区间[0，+∞)中选取步长的初始值α＞0，取α1＝0，α2＝∞，p∈(0，0.5)，q∈(p，1)；计算步长α^k为零时的哈密顿函数H₁＝H(w^k)，及其导数H₁’＝g^kd^k；

6.3.2)计算哈密顿函数H＝H(w^k+α^kd^k)，检验下式是否成立：

H≤H₁+pα^k(g^k)^Td_k    (9)

若成立，则转步骤6.3.4)；否则，利用二次插值公式计算试探点：

${\mathrm{\α}}_S=\mathrm{\α}1+\frac{1}{2}\·\frac{\mathrm{\α}-\mathrm{\α}1}{1+\frac{H_1-H}{(\mathrm{\α}-\mathrm{\α}1)\·H_1^{\′}}}---\left(10\right)$

6.3.3)令α2＝α，α＝α_s，转步骤6.3.2)；

6.3.4)计算哈密顿函数导数H’＝g^k+1d^k以及梯度向量

$g^{k+1}=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^{k+1}}---\left(11\right)$

6.3.5)检验下式是否成立：

|(g^k+1)^Td^k|≤-q·(g^k)^Td_k    (12)

若成立，则停止一维搜索，令α^k＝α，并输出α^k给迭代优化模块；否则，利用二次插值公式计算试探点α_s：

${\mathrm{\α}}_S=\mathrm{\α}+\frac{1}{2}\·\frac{(\mathrm{\α}-\mathrm{\α}1)H^{\′}}{(H_1^{\′}-H^{\′})}---\left(13\right)$

6.3.6)令α1＝α，H₁＝H，H₁’＝H’，α＝α_s，转步骤6.3.2)。

进一步，所述的动态优化方法还包括：将现场智能仪表所采集的工业过程对象的数据传送到DCS***的实时数据库中，在每个采样周期从DCS***的数据库得到的最新数据输出到上位机，并在上位机的初始化处理模块进行初始化处理。

再进一步，所述动态优化方法还包括：在所述步骤5)中得到的最优决策轨线w^*将通过结果输出模块转换为最优控制曲线u^*(t)，并在上位机的人机界面上显示u^*(t)和最优目标值J^*；同时，最优控制曲线u^*(t)将通过通讯接口传给DCS***，并在DCS***中显示所得到的优化结果信息。

本发明的技术构思为：工业过程动态优化，往往需要在求解大规模复杂非线性微分方程***的基础上寻求目标泛函的最优值。通常使用的同步策略法、控制变量参数化法、遗传算法、粒子群算法、模拟退火法等，在处理大规模复杂非线性微分方程***时，往往出现计算结果不准确或收敛缓慢的问题。本发明的动态优化方法，采用非线性共轭梯度法，结合强Wolfe-Powell线搜索进行迭代寻优，并对具有边界约束的控制变量实施无约束转化，使得动态优化问题的求解更加简便、高效，同时也提高了算法收敛的稳定性和准确性。

本发明的有益效果主要表现在：1、能够寻找到工业过程非线性***动态优化的最优解；2、寻优过程准确、高效、稳定。因此，可以广泛地应用于工业过程动态优化的各个领域。

附图说明

图1是本发明所提供的工业过程动态优化***的硬件结构图；

图2是本发明上位机实现动态优化方法的功能结构图。

具体实施方式

下面根据附图具体说明本发明。本发明实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权力要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。

实施例1

参照图1、图2，一种基于非线性共轭梯度法的工业过程动态优化***，包括与工业过程对象1连接的现场智能仪表2、DCS***以及上位机6，所述的DCS***由数据接口3、操作站4、数据库5构成；现场智能仪表2与通讯网络连接，所述通讯网络与数据接口3连接，所述数据接口3与现场总线连接，所述现场总线与操作站4、数据库5和上位机6连接，所述的上位机包括：

约束处理模块8，用于处理优化过程中的控制变量边界约束，采用以下转换方程：

u(t)＝0.5(u_max-u_min)×{sin[w(t)]+1}+u_min    (1)

将带有边界约束u_min≤u(t)≤u_max的控制变量u(t)转化为不受边界约束的变量w(t)，其中下标min、max分别表示最小值和最大值，u_min、u_max分别对应控制变量的下界和上界，并将w(t)作为优化问题的决策变量进行求解；

初始化处理模块9，用于初始参数的设置，决策变量w(t)的离散化和初始赋值，具体步骤如下：

[3.1]设置判断迭代优化是否终止的收敛精度值ξ，当优化目标值迭代误差小于ξ时，停止迭代；取迭代次数k初值为0；

[3.2]将时间域[t₀，tf]平均分成n小段：[t₀，t₁]，[t₁，t₂]，...，[t_n-1，t_n]；每个时间段的长度为h＝(tf-t₀)/n；

[3.3]对决策变量w(t)在[3.2]所述时间分段上进行离散化，将w(t)替换为由n个分段常值组成的向量w，并选取任意常数向量作为决策向量的初始值w⁰；

[3.4]设置一维搜索的初始步长α，并设定一维搜索的相关参数α1，α2，p，q；

ODE求解模块10：用于求解动态优化问题的常微分方程组，为迭代优化过程的梯度计算和收敛条件判断做准备。采取以下步骤来完成：

[4.1]数值求解状态方程：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),w\left(t\right),t]$ 边界条件：x(t₀)＝x0；(2)

其中f表示微分函数向量，x(t)为由m个状态变量组成的向量，x₀为状态向量x在初始时刻t₀的值，采用龙格库塔法由初始值x(t₀)正向积分求出；

[4.2]数值求解协态方程：

$\frac{\mathrm{d\λ}\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂x}-\mathrm{\λ}\left(t\right)\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂x},$ 含边界条件：

其中φ、ψ分别是目标函数 $J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}$ 中的非积分项和定积分项，λ(t)为由m个协态变量组成的向量，λ(tf)为协态向量λ在终端时刻tf的值，采用龙格库塔法由终端值λ(tf)通过反向积分求出；

[4.3]由[4.1]、[4.2]所得的状态向量和协态向量计算出目标函数值：

$J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}---\left(4\right)$

迭代优化模块11：用于搜寻使目标函数J最优的决策向量w，采取以下步骤来完成，上标k均表示迭代次数，初始赋值为零，获取ODE求解模块输出的状态向量、协态向量及目标函数值，即当前目标值J^k，进行如下操作：

[6.1]计算当前迭代的梯度g^k，其中，x(t)和λ(t)分别为ODE求解模块10得出的当前状态向量和协态向量，上标T表示转置：

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^k}---\left(5\right)$

[6.2]计算当前迭代的搜索方向d^k，d^k-1表示第k-1次迭代的搜索方向，β^k是中间参数：

$d^k=\left\{\begin{array}{cc}{-g}^k,& k=1;\\ {-g}^k+{\mathrm{\β}}^k{d}^{k-1},& k\≥2.\end{array}\right.$ 其中， ${\mathrm{\β}}^k=\frac{{\left(g^k\right)}^T(g^k-g^{k-1})}{{\left|\right|g^{k-1}\left|\right|}^2}---\left(6\right)$

[6.3]进行一维搜索，从当前的迭代点w^k出发，沿方向d^k作一维搜索，计算出最佳步长α^*，满足

$H(w^k+{\mathrm{\α}}^{*}\·d^k)=\underset{\mathrm{\α}\≥0}{\min }H(w^k+\mathrm{\α}\·d^k)---\left(7\right)$

其中H表示哈密顿函数，由下式计算出，

表示在[0，+∞)中寻找使H达到最小值的步长α；

H＝ψ[x(t)，w(t)，t]+λ(t)^Tf[x(t)，w(t)，t](8)

[6.4]取当前迭代的搜索步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{B\·\left|\right|d^k\left|\right|},{\mathrm{\α}}^{*}),$ B取区间[5，8]内的整数值；

[6.5]计算下一个迭代点w^k+1＝w^k+α^k·d^k，增加迭代次数：k＝k+1，并把新的迭代点传给ODE求解模块以计算新的状态向量、协态向量及目标函数值J^k+1；

一维搜索模块12：用于寻找使哈密顿函数H达到最小值的步长α^k，返回给迭代优化模块。模块采用强Wolfe-Powell准则进行一维搜索，实现的步骤如下：

6.3.1)选取初始数据：在搜索区间[0，+∞)中选取步长的初始值α＞0，取α1＝0，α2＝∞，p∈(0，0.5)，q∈(p，1)；计算步长α^k为零时的哈密顿函数H₁＝H(w^k)，及其导数H₁’＝g^kd^k；

6.3.2)计算哈密顿函数H＝H(w^k+α^kd^k)，检验下式是否成立：

H≤H₁+pα^k(g^k)^Td_k    (9)

若成立，则转步骤6.3.4)；否则，利用二次插值公式计算试探点：

${\mathrm{\α}}_S=\mathrm{\α}1+\frac{1}{2}\·\frac{\mathrm{\α}-\mathrm{\α}1}{1+\frac{H_1-H}{(\mathrm{\α}-\mathrm{\α}1)\·H_1^{\′}}}---\left(10\right)$

6.3.3)令α2＝α，α＝α_s，转步骤6.3.2)；

6.3.4)计算哈密顿函数导数H’＝g^k+1d^k以及梯度向量

$g^{k+1}=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^{k+1}}---\left(11\right)$

6.3.5)检验下式是否成立：

|(g^k+1)^Td^k|≤-q·(g^k)^Td_k    (12)

若成立，则停止一维搜索，令α^k＝α，并输出α^k给迭代优化模块；否则，利用二次插值公式计算试探点α_s：

${\mathrm{\α}}_S=\mathrm{\α}+\frac{1}{2}\·\frac{(\mathrm{\α}-\mathrm{\α}1)H^{\′}}{(H_1^{\′}-H^{\′})}---\left(13\right)$

6.3.6)令α1＝α，H₁＝H，H₁’＝H’，α＝α_s，转步骤6.3.2)。

收敛判断模块13：用于判定是否终止迭代计算，也就是判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，即收敛条件为：

|J^k-J^k+1|≤ξ    (14)

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；

若上式(14)成立，则将当前的最优决策向量w^*、最优目标值J^*以及迭代次数k保存到结果输出模块；若上式(14)不成立，则保存当前目标值J^k，然后转入迭代优化模块进行下一次迭代求解。

所述的上位机还包括：信号采集模块7，用于设定采样时间，采集现场智能仪表的信号；以及结果输出模块14，用于将迭代优化模块所获得的最优决策轨线为w^*(t)通过式(1)转化为最优控制轨线u^*(t)，然后将此最优控制信号传输给DCS***，并在DCS***中显示所得到的优化结果信息，同时，结果输出模块14还通过DCS***和现场总线将所得到的优化结果信息传递到现场工作站进行显示，并由现场操作站来执行最优操作。

本实施案例的***硬件结构图如附图1所示，所述的动态优化***核心包括带人机界面的上位机6中的约束处理模块8、初始化处理模块9、ODE求解模块10、迭代优化模块11、一维搜索模块12和收敛判断模块13等6大功能模块，此外还包括：现场智能仪表2、DCS***和现场总线。所述的DCS***由数据接口3、操作站4、数据库5组成；工业过程对象1、现场智能仪表2、DCS***、上位机6通过现场总线依次相连，实现信息流的上传和下达，上位机与底层***及时进行信息交换，实现***的在线优化。

实施例2

参照图1和图2，一种基于非线性共轭梯度法的工业过程动态优化，所述的动态优化方法按照以下步骤实施：

1)在DCS***中指定动态优化的状态变量和控制变量，根据实际生产环境的条件和操作限制的条件设定控制变量的上下边界u_max、u_min和DCS的采样周期，并将DCS数据库5中相应各变量的历史数据，控制变量上下边界值u_max、u_min传送给上位机6。

2)在上位机的约束处理模块8中，通过三角函数代换，对受边界约束的控制变量u(t)进行无约束转换，将其替换为中间变量w(t)的函数表达式，即：

u(t)＝0.5(u_max-u_min)×{sin[w(t)]+1}+u_min    (1)

通过上式，控制变量u(t)∈[u_min，u_max]转化为无约束变量w(t)∈R的三角函数表达式。然后，将w(t)作为决策变量进行优化求解，最终求得的w(t)代入(1)式，即得到相应的u(t)。

3)在上位机的初始化处理模块9中，对上位机各模块初始参数进行设置，并对DCS***输入的数据进行初始化处理，按照以下步骤完成：

[3.1]设置判断迭代优化是否终止的收敛精度值ξ，当优化目标值迭代误差小于ξ时，停止迭代；取迭代次数k初值为0；

[3.2]将时间域[t₀，tf]平均分成n小段：[t₀，t₁]，[t₁，t₂]，...，[t_n-1，t_n]；每个时间段的长度为h＝(tf-t₀)/n；

[3.3]对决策变量w(t)在[3.2]所述时间分段上进行离散化，将w(t)替换为由n个分段常值组成的向量w，并选取任意常数向量作为决策向量的初始值w⁰；

[3.4]设置一维搜索的初始步长α，并设定一维搜索的相关参数α1，α2，p，q；

4)将当前迭代的决策变量w＝w^k(迭代次数k为0时，w＝w⁰)代入ODE求解模块10，采用数值积分方法求出当前状态变量、协态变量，并得出相应的当前目标值J^k。实现步骤如下：

[4.1]数值求解状态方程：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),w\left(t\right),t]$ 边界条件：x(t₀)＝x0；(2)

其中f表示微分函数向量，x(t)为由m个状态变量组成的向量，x₀为状态向量x在初始时刻t₀的值，采用龙格库塔法由初始值x(t₀)正向积分求出；

[4.2]数值求解协态方程：

$\frac{\mathrm{d\λ}\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂x}-\mathrm{\λ}\left(t\right)\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂x}$ 边界条件：

其中φ、ψ分别是目标函数 $J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}$ 中的非积分项和定积分项，λ(t)为由m个协态变量组成的向量，λ(tf)为协态向量λ在终端时刻tf的值，采用龙格库塔法由终端值λ(tf)通过反向积分求出；

[4.3]由[4.1]、[4.2]所得的状态向量和协态向量计算出当前目标函数值：

$J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}---\left(4\right)$

5)收敛判断模块13判断收敛条件|J^k-J^k+1|≤ξ(其中J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值)是否满足，若不满足，则保存当前目标值J^k，再转入下一步骤进行下一次迭代优化；若满足，则终止迭代计算，当前决策变量就是最优决策变量(表示为w^*)，当前目标值就是最优目标值(表示为J^*)，保存w^*、J^*及迭代次数k到结果输出模块；

6)由ODE求解模块10得出的状态向量和协态变量值，计算出进一步迭代优化的搜索方向和步长，求解使目标函数J更接近最优的决策向量w，然后返回步骤4)求解ODE。实施迭代优化的步骤如下(上标k均表示迭代次数，初始赋值为零)：

[6.1]计算当前迭代的梯度g^k，其中，x(t)和λ(t)分别为ODE求解模块10得出的当前状态向量和协态向量，上标T表示转置：

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^k}---\left(5\right)$

[6.2]计算当前迭代的搜索方向d^k，d^k-1表示第k-1次迭代的搜索方向，β^k是中间参数：

$d^k=\left\{\begin{array}{cc}{-g}^k,& k=1;\\ {-g}^k+{\mathrm{\β}}^k{d}^{k-1},& k\≥2.\end{array}\right.$ 其中， ${\mathrm{\β}}^k=\frac{{\left(g^k\right)}^T(g^k-g^{k-1})}{{\left|\right|g^{k-1}\left|\right|}^2}---\left(6\right)$

[6.3]进行一维搜索，从当前的迭代点w^k出发，沿方向d^k作一维搜索，计算出最佳步长α^*，满足

$H(w^k+{\mathrm{\α}}^{*}\·d^k)=\underset{\mathrm{\α}\≥0}{\min }H(w^k+\mathrm{\α}\·d^k)---\left(7\right)$

其中H表示哈密顿函数，由下式计算出，

表示在[0，+∞)中寻找使H达到最小值的步长α；

H＝ψ[x(t)，w(t)，t]+λ(t)^Tf[x(t)，w(t)，t](8)

[6.4]取当前迭代的搜索步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{B\·\left|\right|d^k\left|\right|},{\mathrm{\α}}^{*}),$ B取区间[5，8]内的整数值；

[6.5]计算下一个迭代点w^k+1＝w^k+α^k·d^k，增加迭代次数：k＝k+1，并把新的迭代点传给步骤4)的ODE求解模块以计算新的状态向量、协态向量及目标函数值J^k+1；

需要注意的是：上述第[6.3]步骤是依据强Wolfe-Powell准则进行一维搜索，按照以下步骤实现：

6.3.1)选取初始数据：在搜索区间[0，+∞)中选取步长的初始值α＞0，取α1＝0，α2＝∞，p∈(0，0.5)，q∈(p，1)；计算步长α^k为零时的哈密顿函数H₁＝H(w^k)，及其导数H₁’＝g^kd^k；

6.3.2)计算哈密顿函数H＝H(w^k+α^kd^k)，检验下式是否成立：

H≤H₁+pα^k(g^k)^Td_k    (9)

若成立，则转步骤6.3.4)；否则，利用二次插值公式计算试探点：

${\mathrm{\α}}_S=\mathrm{\α}1+\frac{1}{2}\·\frac{\mathrm{\α}-\mathrm{\α}1}{1+\frac{H_1-H}{(\mathrm{\α}-\mathrm{\α}1)\·H_1^{\′}}}---\left(10\right)$

6.3.3)令α2＝α，α＝α_s，转步骤6.3.2)；

6.3.4)计算哈密顿函数导数H’＝g^k+1d^k以及梯度向量

$g^{k+1}=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^{k+1}}---\left(11\right)$

6.3.5)检验下式是否成立：

|(g^k+1)^Td^k|≤-q·(g^k)^Td_k    (12)

若成立，则停止一维搜索，令α^k＝α，并输出α^k给迭代优化模块；否则，利用二次插值公式计算试探点α_s：

${\mathrm{\α}}_S=\mathrm{\α}+\frac{1}{2}\·\frac{(\mathrm{\α}-\mathrm{\α}1)H^{\′}}{(H_1^{\′}-H^{\′})}---\left(13\right)$

6.3.6)令α1＝α，H₁＝H，H₁’＝H’，α＝α_s，转步骤6.3.2)。

7)判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，即收敛条件为：

|J^k-J^k+1|≤ξ    (14)

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；

若上式(14)成立，则将当前的最优决策向量w^*、最优目标值J^*以及迭代次数k保存到结果输出模块；若上式(14)不成立，则保存当前目标值J^k，然后转入迭代优化模块进行下一次迭代求解。

本实施例中，***开始投运，采用以下步骤完成：

1)利用定时器，设置好每次数据检测和采集的时间间隔；

2)现场智能仪表2检测工业过程对象1的数据并传送至DCS***的实时数据库5中，得到最新的变量数据；

3)在上位机6的约束处理模块8中，对控制变量边界约束进行处理，将处理的结果作为初始化处理模块9的输入；

4)在上位机6的初始化处理模块9中，根据实际生产需求和操作限制条件对各模块相关参数和变量进行初始化处理，将处理的结果作为ODE求解模块10和迭代优化模块11的输入；

5)上位机6的ODE求解模块10，根据初始化处理模块9输入的初始决策向量或者迭代优化模块11输入的迭代决策向量求解ODE，所得的状态向量、协态向量和目标值作为迭代优化模块11的输入；

6)上位机6的迭代优化模块11，利用ODE求解模块10得出的变量信息进行梯度计算，并调用一维搜索模块12求取最佳步长，同时结合收敛判断模块13的判定信息实施迭代优化，优化的结果传给结果输出模块14；

7)上位机6的收敛判断模块13，根据收敛条件判定是否终止迭代优化，所得的结果传给迭代优化模块11和结果输出模块14；

8)上位机6的人机界面上显示工业过程最优控制的结果信息，上位机6将所得到的最优控制曲线传给DCS***，并在DCS***的操作站4显示所得到的优化结果信息，同时通过DCS***和现场总线将所得到的优化结果信息传输到现场工作站进行显示，并由现场操作站来执行最优操作。

Claims (8)

1.一种基于非线性共轭梯度法的工业过程动态优化***，包括与工业过程对象连接的现场智能仪表、DCS***和上位机；所述的DCS***包括操作站和数据库；所述现场智能仪表与DCS***连接，所述DCS***与上位机连接，其特征在于：所述的上位机包括：

约束处理模块，用于处理优化过程中的控制变量边界约束，采用以下转换方程：

u(t)＝0.5(u_max-u_min)×{sin[w(t)]+1}+u_min    (1)

将带有边界约束u_min≤u(t)≤u_max的控制变量u(t)转化为不受边界约束的变量w(t)，其中下标min、max分别表示最小值和最大值，u_min、u_max分别对应控制变量的下界和上界，并将w(t)作为优化问题的决策变量进行求解；

初始化处理模块，用于初始参数的设置，决策变量w(t)的离散化和初始赋值，具体步骤如下：

[3.1]设置判断迭代优化是否终止的收敛精度值ξ，当优化目标值迭代误差小于ξ时，停止迭代；取迭代次数k初值为0；

[3.2]将时间域[t₀，tf]平均分成n小段：[t₀，t₁]，[t₁，t₂]，…，[t_n-1，t_n]；每个时间段的长度为h＝(tf-t₀)/n；

[3.3]对决策变量w(t)在[3.2]所述时间分段上进行离散化，将w(t)替换为由n个分段常值组成的向量w，并选取任意常数向量作为决策向量的初始值w⁰；

[3.4]设置一维搜索的初始步长α，并设定一维搜索的相关参数α1，α2，p，q；

ODE求解模块，用于求解动态优化问题的常微分方程组，为迭代优化过程的梯度计算和收敛条件判断做准备，采取以下步骤来完成：

[4.1]数值求解状态方程：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),w\left(t\right),t]$ 边界条件： $x\left(t_0\right)=x0;---\left(2\right)$

其中f表示微分函数向量，x(t)为由m个状态变量组成的向量，x₀为状态向量x在初始时刻t₀的值，采用龙格库塔法由初始值x(t₀)正向积分求出；

[4.2]数值求解协态方程：

$\frac{\mathrm{d\λ}\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂x}-\mathrm{\λ}\left(t\right)\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂x},$ 含边界条件：

其中φ、ψ分别是目标函数 $J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}$ 中的非积分项和定积分项，λ(t)为由m个协态变量组成的向量，λ(tf)为协态向量λ在终端时刻tf的值，采用龙格库塔法由终端值λ(tf)通过反向积分求出；

[4.3]由[4.1]、[4.2]所得的状态向量和协态向量计算出目标函数值：

$J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}---\left(4\right)$

迭代优化模块，用于搜寻使目标函数J最优的决策向量w，采取以下步骤来完成，上标k均表示迭代次数，初始赋值为零，获取ODE求解模块输出的状态向量、协态向量及目标函数值，即当前目标值J^k，进行如下操作：

[6.1]计算当前迭代的梯度g^k，其中，x(t)和λ(t)分别为ODE求解模块得出的当前状态向量和协态向量，上标T表示转置：

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^k}---\left(5\right)$

[6.2]计算当前迭代的搜索方向d^k，d^k-1表示第k-1次迭代的搜索方向，β^k是中间参数：

$d^k=\left\{\begin{array}{c}-g^k,k=1;\\ -g^k+{\mathrm{\β}}^k{d}^{k-1},k\≥2.\end{array}\right.$ 其中， ${\mathrm{\β}}^k=\frac{{\left(g^k\right)}^T(g^k-g^{k-1})}{{\left|\right|g^{k-1}\left|\right|}^2}---\left(6\right)$

[6.3]进行一维搜索，从当前的迭代点w^k出发，沿方向d^k作一维搜索，计算出最佳步长α^*，满足

$H(w^k+{\mathrm{\α}}^{*}\·d^k)=\underset{\mathrm{\α}\≥0}{\min }H(w^k+\mathrm{\α}\·d^k)---\left(7\right)$

其中H表示哈密顿函数，由下式计算出， $\underset{\mathrm{\α}\≥0}{\min }$ 表示在[0，+∞)中寻找使H达到最小值的步长α；

H＝ψ[x(t)，w(t)，t]+λ(t)^Tf[x(t)，w(t)，t]    (8)

[6.4]取当前迭代的搜索步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{B\·\left|\right|d^k\left|\right|},{\mathrm{\α}}^{*}),$ B取区间[5，8]内的整数值；

[6.5]计算下一个迭代点w^k+1＝w^k+α^k·d^k，增加迭代次数：k＝k+1，并把新的迭代点传给ODE求解模块以计算新的状态向量、协

态向量及目标函数值J^k+1；

收敛判断模块，用于判定是否终止迭代计算，即判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，即收敛条件为：

$J^k-J^{k+1}|\≤\mathrm{\ξ}---\left(14\right)$

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；

若上式(14)成立，则将当前的最优决策向量w^*、最优目标值J^*以及迭代次数k保存到结果输出模块；若上式(14)不成立，则保存当前目标值J^k，然后转入迭代优化模块进行下一次迭代求解。

2.如权利要求1所述的基于非线性共轭梯度法的工业过程动态优化***，其特征在于：所述上位机还包括：

一维搜索模块，用于寻找使哈密顿函数H达到最小值的步长α^k，返回给迭代优化模块；依据强Wolfe-Powell准则进行一维搜索，按照以下步骤实现：

6.3.1)选取初始数据：在搜索区间[0，+∞)中选取步长的初始值α＞0，取α1＝0，α2＝∞，p∈(0，0.5)，q∈(p，1)；计算步长α^k为零时的哈密顿函数H₁＝H(w^k)，及其导数H₁’g^kd^k；

6.3.2)计算哈密顿函数H＝H(w^k+α^kd^k)，检验下式是否成立：

H≤H₁+pα^k(g^k)Td^k(9)

若成立，则转步骤6.3.4)；否则，利用二次插值公式计算试探点：

${\mathrm{\α}}_s=\mathrm{\α}1+\frac{1}{2}\·\frac{\mathrm{\α}-\mathrm{\α}1}{1+\frac{H_1-H}{(\mathrm{\α}-\mathrm{\α}1)\·{H_{}}_1^{\text{'}}}}---\left(10\right)$

6.3.3)令α2＝α，α＝α_s，转步骤6.3.2)；

6.3.4)计算哈密顿函数导数H’＝g^k+1d^k以及梯度向量

$g^{k+1}=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^{k+1}}---\left(11\right)$

6.3.5)检验下式是否成立：

${\left(g^{k+1}\right)}^T{d}^k{\left(g^k\right)}^T{d}_k---\left(12\right)$

若成立，则停止一维搜索，令α^k＝α，并输出α^k给迭代优化模块；否则，利用二次插值公式计算试探点α_s：

${\mathrm{\α}}_S=\mathrm{\α}+\frac{1}{2}\·\frac{(\mathrm{\α}-\mathrm{\α}1)H^{\text{'}}}{(H_1^{\text{'}}-H^{\text{'}})}---\left(13\right)$

6.3.6)令α1＝α，H₁＝H，H₁’＝H’，α＝α_s，转步骤6.3.2)。

3.如权利要求1或2所述的基于非线性共轭梯度法的工业过程动态优化***，其特征在于：所述的上位机还包括：信号采集模块，用于设定采样时间，采集现场智能仪表的信号。

4.如权利要求1或2所述的基于非线性共轭梯度法的工业过程动态优化***，其特征在于：所述的上位机还包括：结果输出模块，用于将迭代优化模块所获得的最优决策轨线为w^*(t)通过式(1)转化为最优控制轨线u^*(t)，然后将此最优控制信号传输给DCS***，并在DCS***中显示所得到的优化结果信息。

5.一种用如权利要求1所述的基于非线性共轭梯度法的工业过程动态优化***实现的动态优化方法，其特征在于：所述的动态优化方法包括以下步骤：

1)在DCS***中指定动态优化的状态变量和控制变量，根据实际生产环境的条件和操作限制的条件设定控制变量的上下边界u_max、u_min和DCS的采样周期，并将DCS数据库中相应各变量的历史数据，控制变量上下边界值u_max、u_min传送给上位机；

2)通过三角函数代换，对受边界约束的控制变量u(t)进行无约束转换，将其替换为另一无约束变量w(t)的函数表达式，即：

u(t)＝0.5(u_max-u_min)×{sin[w(t)]+1}+u_min    (1)

通过上式，控制变量u(t)∈[u_min，u_max]转化为无约束变量w(t)∈R的三角函数表达式；然后，将w(t)作为决策变量进行优化求解，最终求得的w(t)通过代入(1)式，即得到相应的u(t)；

3)对模块初始参数进行设置，并对DCS***输入的数据进行初始化处理，按照以下步骤完成：

[3.1]设置判断迭代优化是否终止的收敛精度值ξ，当优化目标值

迭代误差小于ξ时，停止迭代；取迭代次数k初值为0；

[3.2]将时间域[t₀，tf]平均分成n小段：[t₀，t₁]，[t₁，t₂]，…，[t_n-1，t_n]；每个时间段的长度为h＝(tf-t₀)/n；

[3.3]对决策变量w(t)在[3.2]所述时间分段上进行离散化，将w(t)替换为由n个分段常值组成的向量w，并选取任意常数向量作为决策向量的初始值w⁰；

[3.4]设置一维搜索的初始步长α，并设定一维搜索的相关参数α1，α2，p，q；

4)当前迭代的决策变量w＝w^k，迭代次数k为0时，w＝w⁰，采用数值积分方法求出当前状态变量、协态变量，并得出相应的当前目标值J^k。实现步骤如下：

[4.1]数值求解状态方程：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),w\left(t\right),t]$ 边界条件：x(t₀)＝x0；(2)

其中f表示微分函数向量，x(t)为由m个状态变量组成的向量，x₀为状态向量x在初始时刻t₀的值，采用龙格库塔法由初始值x(t₀)正向积分求出；

[4.2]数值求解协态方程：

$\frac{\mathrm{d\λ}\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂x}-\mathrm{\λ}\left(t\right)\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂x}$ 边界条件：

其中φ、ψ分别是目标函数 $J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}$ 中的非积分项和定积分项，λ(t)为由m个协态变量组成的向量，λ(tf)为协态向量λ在终端时刻tf的值，采用龙格库塔法由终端值λ(tf)通过反向积分求出；

[4.3]由[4.1]、[4.2]所得的状态向量和协态向量计算出当前目标函数值：

$J^k=J=\mathrm{\φ}\left[x\left(\mathrm{tf}\right)\right]+{\∫}_0^{\mathrm{tf}}\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]\mathrm{dt}---\left(4\right)$

5)判断收敛条件|J^k-J^k+1|≤ξ是否满足，其中J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值，若不满足，则保存当前目标值J^k，再转入下一步骤进行下一次迭代优化；若满足，则终止迭代计算，当前决策向量就是最优决策向量w^*，当前目标值就是最优目标值J^*，保存w^*、J^*及迭代次数k到结果输出模块；

6)由状态向量和协态变量值计算出进一步迭代优化的搜索方向和步长，求解使目标函数J更接近最优的决策向量w，然后返回步骤4)求解ODE；实施迭代优化的步骤如下，上标k均表示迭代次数，初始赋值为零：

[6.1]计算当前迭代的梯度g^k，其中，x(t)和λ(t)分别为步骤4)得出的当前状态向量和协态向量，上标T表示转置：

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}\left(t\right)}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^k}---\left(5\right)$

[6.2]计算当前迭代的搜索方向d^k，d^k-1表示第k-1次迭代的搜索方向，β^k是中间参数：

$d^k=\left\{\begin{array}{c}-g^k,k=1;\\ -g^k+{\mathrm{\β}}^k{d}^{k-1},k\≥2.\end{array}\right.$ 其中， ${\mathrm{\β}}^k=\frac{{\left(g^k\right)}^T(g^k-g^{k-1})}{{\left|\right|g^{k-1}\left|\right|}^2}---\left(6\right)$

[6.3]进行一维搜索，从当前的迭代点w^k出发，沿方向d^k作一维搜索，计算出最佳步长α^*，满足

$H(w^k+{\mathrm{\α}}^{*}\·d^k)=\underset{\mathrm{\α}\≥0}{\min }H(w^k+\mathrm{\α}\·d^k)---\left(7\right)$

其中H表示哈密顿函数，由下式计算出， $\underset{\mathrm{\α}\≥0}{\min }$ 表示在[0，+∞)中寻找使H达到最小值的步长α；

H＝ψ[x(t)，w(t)，t]+λ(t)^Tf[x(t)，w(t)，t]    (8)

[6.4]取当前迭代的搜索步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{B\·\left|\right|d^k\left|\right|},{\mathrm{\α}}^{*}),$ B取区间[5，8]内的整数值；

[6.5]计算下一个迭代点w^k+1＝w^k+α^k·d^k，增加迭代次数：k＝k+1，并把新的迭代点传给步骤4)以计算新的状态向量、协态向量及目标函数值J^k+1；

7)判断当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值是否小于等于设定的收敛精度ξ，即收敛条件为：

|J^k-J^k+1|≤ξ                         (14)

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；

若上式(14)成立，则将当前的最优决策向量w^*、最优目标值J^*以及迭代次数k保存到结果输出模块；若上式(14)不成立，则保存当前目标值J^k，然后转入迭代优化模块进行下一次迭代求解。

6.如权利要求5所述的动态优化方法，其特征在于：上述步骤[6.3]是依据强Wolfe-Powell准则进行一维搜索，按照以下步骤实现：

6.3.1)选取初始数据：在搜索区间[O，+∞)中选取步长的初始值α＞O，取α1＝0，α2＝∞，p∈(0，0.5)，q∈(p，1)；计算步长α^k为零时的哈密顿函数H₁＝H(w^k)，及其导数H₁’＝g^kd^k；

6.3.2)计算哈密顿函数H＝H(w^k+α^kd^k)，检验下式是否成立：

H≤H₁+pα^k(g^k)^Td^k(9)

若成立，则转步骤6.3.4)；否则，利用二次插值公式计算试探点：

${\mathrm{\α}}_s=\mathrm{\α}1+\frac{1}{2}\·\frac{\mathrm{\α}-\mathrm{\α}1}{1+\frac{H_1-H}{(\mathrm{\α}-\mathrm{\α}1)\·{H_{}}_1^{\text{'}}}}---\left(10\right)$

6.3.3)令α2＝α，α＝α_s，转步骤6.3.2)；

6.3.4)计算哈密顿函数导数H’＝g^k+1d^k以及梯度向量

$g^{k+1}=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}+{\mathrm{\λ}}^T\·\frac{\∂f[x\left(t\right),w\left(t\right),t]}{\∂w\left(t\right)}\}{|}_{w\left(t\right)=w^{k+1}}---\left(11\right)$

6.3.5)检验下式是否成立：

|(g^k+1)^Td^k|≤-q·(g^k)^Td^k(12)

若成立，则停止一维搜索，令α^k＝α，并输出α^k给迭代优化模块；否则，利用二次插值公式计算试探点α_s：

${\mathrm{\α}}_S=\mathrm{\α}+\frac{1}{2}\·\frac{(\mathrm{\α}-\mathrm{\α}1)H^{\text{'}}}{(H_1^{\text{'}}-H^{\text{'}})}---\left(13\right)$

6.3.6)令α1＝α，H₁＝H，H₁’＝H’，α＝α_s，转步骤6.3.2)。

7.如权利要求5或6所述的动态优化方法，其特征在于：所述的动态优化方法还包括：将现场智能仪表所采集的工业过程对象的数据传送到DCS***的实时数据库中，在每个采样周期从DCS***的数据库得到的最新数据输出到上位机，并在上位机的初始化处理模块进行初始化处理。

8.如权利要求5或6所述的动态优化方法，其特征在于：所述动态优化方法还包括：在所述步骤5)中得到的最优决策轨线w^*将通过结果输出模块转换为最优控制曲线u^*(t)，并在上位机的人机界面上显示u^*(t)和最优目标值J^*；同时，最优控制曲线u^*(t)将通过通讯接口传给DCS***，并在DCS***中显示所得到的优化结果信息。

Images