CN101763082B - 一种有效的基于两点步长梯度法的工业过程动态优化系统及方法 - Google Patents

Abstract

一种有效的基于两点步长梯度法的工业过程动态优化系统，包括与工业过程对象连接的现场智能仪表、DCS系统和上位机，上位机包括：初始化模块，用于初始参数的设置，决策变量z(t)的离散化和初始赋值；约束转换模块，用于通过中间变量处理优化过程中的控制变量边界约束；ODE求解模块，用于求解工业过程动态优化问题的常微分方程组；迭代寻优模块，用于根据两点步长梯度法搜寻使目标函数J最优的决策向量z^＊，迭代控制模块，用于控制迭代寻优模块的功能状态，并保存迭代优化的结果，根据收敛条件停止迭代优化。以及提供一种基于两点步长梯度法的工业过程动态优化方法。本发明能够准确找到全局最优解且求解高效、适用性广。

Description

一种有效的基于两点步长梯度法的工业过程动态优化系统及方法

技术领域

本发明涉及最优控制领域，尤其是一种基于两点步长梯度法的工业过程动态优化系统。 

背景技术

工业过程从严格意义上说均是一个动态的过程，即描述过程的状态变量(如流量、温度、压力、液位等)随时间的演进、空间的转移而发生改变。动态过程由微分方程或差分方程描述，称为动态模型。动态优化就是对动态模型中的操作变量实施控制，使得过程的性能指标达到最优。 

从20世纪60年代至今，动态优化(也常称为最优控制)在理论研究和实际应用领域的发展都十分令人瞩目，它已经由一种学术界乐于探究的方法论演变为一项正在并将继续对过程工业产生重大影响的技术。随着ASPEN CUSTOM MODELER，FEMLAB，gPROMS等强大动态仿真商业软件的开发和应用，工业过程动态建模已得到越来越广泛的发展，而动态优化将在这些动态仿真工具的基础上实现工程决策的自动化。 

动态优化的应用领域非常广泛，一个典型的工程在线应用是求解非线性模型预测控制(Nonlinear Model Predictive Control，NMPC)中的优化问题。NMPC的本质特征是在从现时到未来的有限时域内滚动优化开环性能目标，从这个意义上说，能够高效地求解优化问题以实时获得最优控制轨迹的动态优化算法是实行NMPC的一个关键因素。由于实际工业应用要求NMPC求解的动态过程模型日趋大规模和精细化，所以发展高效的动态优化算法的重要性已经越来越突出。 

目前的一些动态优化方法，如迭代动态规划法、随机算法、同步策略法、控制变量参数化法等，虽然能够找到工业过程动态优化问题的最优解，但是往往出现收敛缓慢和不稳定问题，还可能陷入局部最优，很难既保证所得解的全局最优性，又使得动态优化问题的求解稳定、快速。 

发明内容

为了克服已有的工业过程动态优化系统和方法很难既准确又快速地找到最优解、适用性差的不足，本发明提供了一种能够准确找到全局最优解且求解高效、适用性广的基于两点梯度法的工业过程动态优化系统及方法。 

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

一种有效的基于两点步长梯度法的工业过程动态优化系统，包括与工业过程对象连接的现场智能仪表、DCS系统和上位机，所述的DCS系统包括数据库和控制站，所述现场智能仪表与DCS系统连接，所述DCS系统与上位机连接，所述的上位机包括： 

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

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

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

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

(3.4)设置迭代搜索的初始步长α0(通常为(0，1]区间的值)； 

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

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

将带有边界约束u_min≤u(t)≤u_max的控制变量u(t)替换为不受边界约束的中间变量z(t)的三角函数表达式，其中，下标min、max分别表示最小值和最大值，u_min、u_max分别对应控制变量的下界和上界；并将z(t)作为动态优化问题的决策变量进行求解； 

ODE求解模块，用于求解工业过程动态优化问题的常微分方程组 (Ordinary Differential Equations，ODE)，为迭代寻优模块的梯度计算提供状态变量和协态变量信息，也为迭代控制模块的收敛条件判断提供目标函数信息，采取以下步骤来完成： 

(4.1)数值求解状态方程组： 

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],$ x_i(t₀)＝x_i0，其中i＝1，2，...，m                            (2) 

其中，f表示微分函数向量，x(t)为m个状态变量组成的向量，x_i(t)表示第i个状态变量，x_i0为状态变量x_i在初始时刻t₀的值，采用龙格库塔法由状态变量初始值x_i0通过正向积分求出状态变量在每个离散时刻的值x_i(t_i)，i＝1，2，...，m，j＝1，2，...，n； 

(4.2)数值求解协态方程组： 

$\frac{d{\mathrm{\λ}}_i\left(t\right)}{\mathrm{dt}}=-\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{{\∂x}_i}-{\mathrm{\λ}}_i\left(t\right)\·\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂x_i},$

其中， 

ψ分别是给定的目标函数 

ψ[x(t)，z(t)，t]dt的非积分项和定积分项，λ_i(t)为第i个协态变量，λ(t)为m个协态变量组成的向量，λ_i(tf)为协态变量λ_i在终端时刻tf的值，采用龙格库塔法由协态变量终端值λ_i(tf)通过反向积分求出协态变量在每个离散时刻的值λ_i(t_j))，i＝1，2，...，m；j＝n-1，n-2，...，1，0； 

(4.3)由所得的状态向量和决策向量计算出目标函数值： 

迭代寻优模块，用于调用ODE求解模块，保存所得的状态向量、协态向量及目标函数值，所述目标函数值为当前目标值J^k；根据两点步长梯度法搜寻使目标函数J最优的决策向量z^*，采取以下步骤来完成，上标k均表示迭代次数，初始赋值为零： 

(5.1)计算当前梯度g^k，上标T表示向量或矩阵的转置： 

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}\}{|}_{z\left(t\right)=z^k,t=t_j},j=\mathrm{0,1,2},...,n---\left(5\right)$

(5.2)保存当前迭代点z^k及梯度信息g^k； 

(5.3)若k＝0，则搜索步长α^k取为初始值，即α^k＝α0，转步骤(5.5)； 

否则，依据两点步长梯度法，利用迭代当前点和前一点的信息来确定步长因子l^k： 

$l^k=\frac{{\left(s^{k-1}\right)}^T\·y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(6\right)$

其中，s^k-1表示当前迭代点与前一迭代点的误差，计算式为： 

s^k-1＝z^k-z^k-1        (7) 

y^k-1表示当前迭代点与前一迭代点的梯度误差，计算式为： 

y^k-1＝g^k-g^k-1        (8) 

(5.4)取最佳步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{D\·\max \left(g^k\right)},l^k)---\left(9\right)$

其中，D取区间[5，8]内的整数值； 

(5.5)计算下一个迭代点： 

z^k+1＝z^k-α^k·g^k    (10) 

把新的迭代点z^k+1传给ODE求解模块以计算新的目标函数值J^k+1，然后进入迭代控制模块； 

迭代控制模块，用于控制迭代寻优模块的功能状态，并保存迭代优化的结果；首先判断是否满足收敛条件： 

|J^k-J^k+1|≤ζ        (11) 

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；若上式(11)成立，表明当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值不超过设定的收敛精度ζ，则停止迭代优 化计算，z^k+1就是最优决策向量z^*，J^k+1就是最优目标值J^*，将z^*、J^*以及相应迭代次数(k+1)保存到结果输出模块；若上式(11)不成立，则保存目标值J^k+1，取k＝k+1，然后返回迭代寻优模块进行新一轮的迭代求解。 

作为优选的一种方案：所述的上位机还包括：信息采集模块，用于设定采样时间，采集由现场智能仪表上传的工业过程对象的动态信息。 

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

一种用所述的基于两点步长梯度法的工业过程动态优化系统实现的动态优化方法，所述的动态优化方法包括以下步骤： 

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

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

u(t)＝0.5(u_max-u_min)×{cos[z(t)]+1}+u_min；        (1) 

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

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

(3.1)在上位机的迭代控制模块，设置收敛精度ζ值，一般取为10^-6即可；设置优化迭代次数k初始计数为0； 

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

(3.3)对决策变量z(t)在步骤(3.2)所述时间分段上进行离散化，使用由n个分段常值组成的决策向量z来表示z(t)，并选取决策向量的初始值z⁰，取为简单的常数向量； 

(3.4)设置设置迭代搜索的初始步长α0＞0(通常为(0，1]区间的值)； 

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

(4.1)数值求解状态方程组： 

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],x_i\left(t_0\right)=x_{i0}(i=\mathrm{1,2},...,m)---\left(2\right)$

其中，f表示微分函数向量，x(t)为m个状态变量组成的向量，x_i(t)表示第i个状态变量，x_i0为状态变量x_i在初始时刻t₀的值，采用龙格库塔法由状态变量初始值x_i0通过正向积分求出状态变量在每个离散时刻的值x_i(t_j)，i＝1，2，...，m，j＝1，2，...，n； 

(4.2)数值求解协态方程组： 

$\frac{d{\mathrm{\λ}}_i\left(t\right)}{\mathrm{dt}}=\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂x_i}-{\mathrm{\λ}}_i\left(t\right)\·\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂x_i},$

其中， ψ分别是给定的目标函数 

的非积分项和定积分项，λ_i(t)为第i个协态变量，λ(t)为m个协态变量组成的向量，λ_i(tf)为协态变量λ_i在终端时刻tf的值，采用龙格库塔法由协态变量终端值λ_i(tf)通过反向积分求出协态变量在每个离散时刻的值λ_i(t_j)，i＝1，2，...，m；j＝n-1，n-2，...，1，0； 

(4.3)由所得的状态向量和决策向量计算出目标函数值： 

5)由ODE求解模块得出的状态向量和协态变量信息计算出迭代优化的搜索方向和步长，求解使目标函数J更接近最优的决策向量 

实施迭代寻优的步骤如下，上标k均表示迭代次数，初始赋值为零： 

(5.1)计算当前梯度g^k，即迭代优化的搜索方向，其中，x(t)和λ(t)分别为求解ODE得出的当前状态向量和协态向量，上标T表示向量或矩阵的转置： 

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}\}{|}_{z\left(t\right)=z^k,t=t_j},j=\mathrm{0,1,2},...,n---\left(5\right)$

(5.2)保存当前迭代点z^k及梯度信息g^k； 

(5.3)若k＝0，则搜索步长α^k取为初始值，即α^k＝α0，转步骤(5.5)；否则，依据两点步长梯度法，利用迭代当前点和前一点的信息来确定步长因子l^k： 

$l^k=\frac{{\left(s^{k-1}\right)}^T\·y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(6\right)$

其中，s^k-1表示当前迭代点与前一迭代点的误差，计算式为： 

s^k-1＝z^k-z^k-1        (7) 

y^k-1表示当前迭代点与前一迭代点的梯度误差，计算式为： 

y^k-1＝g^k-g^k-1        (8) 

(5.4)取最佳步长 ${\mathrm{\α}}^k=\min [\frac{\mathrm{\π}}{D\·\max \left(g^k\right)},l^k],---\left(9\right)$

其中D取区间[5，8]内的整数值； 

(5.5)计算下一个迭代点： 

z^k+1＝z^k-α^k·g^k    (10) 

并将z^k+1传给ODE求解模块计算出相应的目标函数值J^k+1； 

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

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

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

本发明的有益效果主要表现在：1、能够寻找到工业过程非线性系统动态优化的全局最优解；2、具有很高的求解效率，收敛稳定性好。因此，在工业过程动态仿真和最优控制的各个领域都具有广泛的应用前景。 

附图说明

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

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

具体实施方式

下面根据附图具体说明本发明。 

参照图1、图2，一种有效的基于两点步长梯度法的工业过程动态优化系统，包括与工业过程对象1连接的现场智能仪表2、DCS系统以及上位机6，所述的DCS系统由总线接口3、控制站4、数据库5构成；现场智能仪表2与现场总线连接，所述现场总线与总线接口3连接，所述总线接口分别与控制站4、数据库5和上位机6连接，所述的上位机包括： 

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

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

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

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

(3.4)设置迭代搜索的初始步长α0(通常为(0，1]区间的值)；约束转换模块8：通过中间变量处理优化过程中的控制变量边界约束，采用以下转换方程： 

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

将带有边界约束u_min≤u(t)≤u_max的控制变量u(t)替换为不受边界约束的中间变量z(t)的三角函数表达式，，(其中下标min、max分别表示最小值和最大值，u_min、u_max分别对应控制变量的下界和上界；并将z(t)作为动态优化问题的决策变量进行求解。 

ODE求解模块10：用于求解工业过程动态优化问题的常微分方程组，为迭代寻优模块11的梯度计算提供状态变量和协态变量信息，也为迭代控制模块12的收敛条件判断提供目标函数信息，采取以下步骤来完成： 

(4.1)数值求解状态方程组： 

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],x_i\left(t_0\right)=x_{i0}(i=\mathrm{1,2},...,m)---\left(2\right)$

其中，f表示微分函数向量，x(t)为m个状态变量组成的向量，x_i(t)表示第i个状态变量，x_i0为状态变量x_i在初始时刻t₀的值，采用龙格库塔法由状态变量初始值x_i0通过正向积分求出状态变量在每个离散时刻的值x_i(t_j)，i＝1，2，...，m，j＝1，2，...，n； 

(4.2)数值求解协态方程组： 

$\frac{d{\mathrm{\λ}}_i\left(t\right)}{\mathrm{dt}}=\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂x_i}-{\mathrm{\λ}}_i\left(t\right)\·\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂x_i},$

其中， 

ψ分别是给定的目标函数 

的非积分项和定积分项，λ_i(t)为第i个协态变量，λ(t)为m个协态变量组成的向量，λ_i(tf)为协态变量λ_i在终端时刻tf的值，采用龙格库塔法由协态变量终端值λ_i(tf)通过反向积分求出协态变量在每个离散时刻的值λ_i(t_j)，i＝1，2，...，m；j＝n-1，n-2，...，1，0； 

(4.3)由所得的状态向量和决策向量计算出目标函数值： 

迭代寻优模块11：用于调用ODE求解模块10，保存所得的状态向量、协态向量及目标函数值(即当前目标值J^k)，根据两点步长梯度法搜寻使目标函数J最优的决策向量z^*，采取以下步骤来完成，上标k均表示迭代次数，初始赋值为零： 

(5.1)计算当前梯度g^k(上标T表示向量或矩阵的转置)，即迭代优化的搜索方向： 

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}\}{|}_{z\left(t\right)=z^k,t=t_j},j=\mathrm{0,1,2},...,n---\left(5\right)$

(5.2)保存当前迭代点z^k及梯度信息g^k； 

(5.3)若k＝0，则搜索步长α^k取为初始值，即α^k＝α0，转步骤(5.5)；否则，依据两点步长梯度法，利用迭代当前点和前一点的信息来确定步长因子l^k： 

$l^k=\frac{{\left(s^{k-1}\right)}^T\·y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(6\right)$

其中，s^k-1表示当前迭代点与前一迭代点的误差，计算式为： 

s^k-1＝z^k-z^k-1        (7) 

y^k-1表示当前迭代点与前一迭代点的梯度误差，计算式为： 

y^k-1＝g^k-gk^-1        (8) 

(5.4)取最佳步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{D\·\max \left(g^k\right)},l^k)---\left(9\right)$

其中D取区间[5，8]内的整数值； 

(5.5)计算下一个迭代点： 

z^k+1＝z^k-α^k·g^k    (10) 

(5.6)把新的迭代点z^k+1传给ODE求解模块以计算新的目标函数值J^k+1，然后进入迭代控制模块12； 

迭代控制模块12：用于控制迭代寻优模块11的功能状态，并保存迭代优化的结果。首先判断是否满足收敛条件： 

|J^k-J^k+1|≤ζ 

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值。若上式(11)成立，表明当前迭代得到的目标值与前一次迭代所得目标值的误差绝对值不超过设定的收敛精度ζ，则停止迭代优化计算，z^k+1就是最优决策向量z^*(＝z^k+1)，J^k+1就是最优目标值J^*(＝J^k+1)，将z^*、J^*以及相应迭代次数(k+1)保存到结果输出模块13；若上式(11)不成立，则保存目标值J^k+1，取k＝k+1，然后返回迭代寻优模块11进行新一轮的迭代求解。 

所述的上位机6还包括信息采集模块7，用于设定采样时间，采集由现场智能仪表上传的工业过程对象的动态信息。 

所述的上位机6还包括结果输出模块13，用于将迭代寻优模块得到的最优决策轨线z^*(t)(由向量z^*分时段表示)通过式(1)转化为最优控制轨线u^*(t)，然后将u^*(t)传输给DCS系统，并在DCS系统中显示所得到的优化结果信息。 

本实施案例的系统硬件结构图如附图1所示，所述的动态优化系统核心包括带人机界面的上位机6中的约束转换模块8、初始化模块9、 ODE求解模块10、迭代寻优模块11、迭代控制模块12等5大功能模块，此外还包括：现场智能仪表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)∈[u_min，u_max]转化为另一无约束变量z(t)的函数表达式，即： 

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

然后，将z(t)作为决策变量进行优化求解，最终求得的z(t)代入(1)式，即得到相应的u(t)； 

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

(3.1)在上位机的迭代控制模块12，设置收敛精度ζ值，一般取为10^-6即可；设置优化迭代次数k初始计数为0； 

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

(3.3)对决策变量z(t)在(2.2)所述时间分段上进行离散化，使用由n个分段常值组成的向量z(即决策向量)来表示z(t)，并选取决策向量的初始值z⁰，可取为简单的常数向量； 

(3.4)设置迭代搜索的初始步长α0＞0(通常为(0，1]区间的值)； 

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

(4.1)数值求解状态方程组： 

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],x_i\left(t_0\right)=x_{i0}(i=\mathrm{1,2},...,m)---\left(2\right)$

其中，f表示微分函数向量，x(t)为m个状态变量组成的向量，x_i(t)表示第i个状态变量，x_i0为状态变量x_i在初始时刻t₀的值，采用龙格库塔法由状态变量初始值x_i0通过正向积分求出状态变量在每个离散时刻的值x_i(t_j)，i＝1，2，...，m，j＝1，2，...，n； 

(4.2)数值求解协态方程组： 

$\frac{d{\mathrm{\λ}}_i\left(t\right)}{\mathrm{dt}}=\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂x_i}-{\mathrm{\λ}}_i\left(t\right)\·\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂x_i},$

其中， 

ψ分别是给定的目标函数 

的非积分项和定积分项，λ_i(t)为第i个协态变量，λ(t)为m个协态变量组成的向量，λ_i(tf)为协态变量λ_i在终端时刻tf的值，采用龙格库塔法由协态变量终端值λ_i(tf)通过反向积分求出协态变量在每个离散时刻的值λ_i(t_j)，i＝1，2，...，m；j＝n-1，n-2，...，1，0； 

(4.3)由所得的状态向量和决策向量计算出目标函数值： 

5)由ODE求解模块10得出的状态向量和协态变量信息，计算出迭代优化的搜索方向和步长，求解使目标函数J更接近最优的决策向量z^*，实施迭代寻优的步骤如下，上标k均表示迭代次数，初始赋值为零： 

(5.1)计算当前梯度g^k，即迭代优化的搜索方向，其中，x(t)和λ(t)分别为求解ODE得出的当前状态向量和协态向量，上标T表示向量或矩阵的转置： 

$g^k=\{\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}\}{|}_{z\left(t\right)=z^k,t=t_j},j=\mathrm{0,1,2},...,n---\left(5\right)$

(5.2)保存当前迭代点z^k及梯度信息g^k； 

(5.3)若k＝0，则搜索步长α^k取为初始值，即α^k＝α0，转步骤(5.5)；否则，依据两点步长梯度法，利用迭代当前点和前一点的信息来确定步长因子l^k： 

$l^k=\frac{{\left(s^{k-1}\right)}^T\·y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(6\right)$

其中，s^k-1表示当前迭代点与前一迭代点的误差，计算式为： 

s^k-1＝z^k-z^k-1        (7) 

y^k-1表示当前迭代点与前一迭代点的梯度误差，计算式为： 

y^k-1＝g^k-g^k-1        (8) 

(5.4)取最佳步长 ${\mathrm{\α}}^k=\min [\frac{\mathrm{\π}}{D\·\max \left(g^k\right)},l^k]---\left(9\right)$

其中，D取区间[5，8]内的整数值； 

(5.5)计算下一个迭代点： 

z^k+1＝z^k-α^k·g^k   (10) 

并将z^k+1传给ODE求解模块10计算出相应的目标函数值J^k+1； 

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

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

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

本实施例中，系统开始投运，具体的过程为： 

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

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

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

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

5)上位机6的ODE求解模块10，根据迭代寻优模块11输入的初始决策向量或者迭代决策向量求解ODE模型，所得的状态向量、协态向量和目标值再传回迭代寻优模块11； 

6)上位机6的迭代寻优模块11，依据约束转换模块8的变量代换关系和ODE求解模块10得出的变量信息进行梯度计算，并根据迭代控制模块12的判定结果实施迭代优化，优化的结果传给结果输出模块13； 

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

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

上述实施例用来解释说明本发明，而不是对本发明进行限制，在 本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。 

Claims (5)

1.一种有效的基于两点步长梯度法的工业过程动态优化系统，包括与工业过程对象连接的现场智能仪表、DCS系统和上位机，所述的DCS系统包括数据库和控制站，所述现场智能仪表与DCS系统连接，所述DCS系统与上位机连接，其特征在于：所述的上位机包括：

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

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

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

(3.3)设置判断迭代优化是否终止的收敛精度值ζ，当迭代优化的目标函数值迭代的误差绝对值小于等于ζ时，停止迭代；取迭代次数k初始值为0；

(3.4)设置迭代搜索的初始步长α0；

约束转换模块，用于通过决策变量z(t)处理优化过程中的控制变量u(t)边界约束，采用以下转换方程：

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

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

ODE求解模块，用于求解工业过程动态优化问题的常微分方程组，采取以下步骤来完成：

(4.1)数值求解状态方程组：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],x_i\left(t_0\right)=x_{i0},$ 其中i＝1，2，...，m                             (2)

其中，f表示微分函数向量，x(t)为m个状态变量组成的状态向量，x_i(t)表示第i个状态变量，x_i0为状态变量x_i(t)在初始时刻t₀的初始值，采用龙格库塔法由初始值x_i0通过正向积分求出状态变量x_i(t)在每个离散时刻的值x_i(t_j)，其中，i＝1，2，...，m，j＝1，2，...，n；

(4.2)数值求解协态方程组：

其中，

分别是给定的目标函数值的非积分项和定积分项；λ_i(t)为第i个协态变量，λ(t)为m个协态变量组成的协态向量，λ_i(tf)为协态变量λ_i(t)在终端时刻tf的终端值，采用龙格库塔法由终端值λ_i(tf)通过反向积分求出协态变量λ_i(t)在每个离散时刻的值λ_i(t_j)，其中，i＝1，2，...，m；j＝n-1，n-2，...，1，0；

(4.3)由所得的状态向量和决策向量计算出目标函数值：

迭代寻优模块，用于调用ODE求解模块，保存所得的状态向量、协态向量及目标函数值J，所述目标函数值J为当前目标函数值J^k；根据两点步长梯度法搜寻使目标函数值J最优的最优决策向量z^*，采取以下步骤来完成，上标k均表示迭代次数，初始赋值为零：

(5.1)计算当前梯度g^k，上标T表示向量或矩阵的转置：

$g^k={\{\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}\}|}_{z\left(t\right)=z^k,t=t_j},j=\mathrm{0,1,2},...,n---\left(5\right)$

(5.2)保存当前迭代点z^k及梯度信息g^k；

(5.3)若k＝0，则搜索步长α^k取为初始值，即α^k＝α0，转步骤(5.5)；否则，依据两点步长梯度法，利用当前迭代点z^k和前一迭代点z^k-1的信息来确定步长因子l^k：

$l^k=\frac{{\left(s^{k-1}\right)}^T\·y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(6\right)$

其中，s^k-1表示当前迭代点z^k与前一迭代点z^k-1的误差，计算式为：

s^k-1＝z^k-z^k-1                     (7)

y^k-1表示当前迭代点z^k与前一迭代点z^k-1的梯度误差，计算式为：

y^k-1＝g^k-g^k-1                     (8)

(5.4)取最佳步长 ${\mathrm{\α}}^k=\min (\frac{\mathrm{\π}}{D\·\max \left(g^k\right)},l^k)---\left(9\right)$

其中，D取区间[5，8]内的整数值；

(5.5)计算下一个迭代点：

z^k+1＝z^k-α^k·g^k                    (10)

把新的迭代点z^k+1传给ODE求解模块以计算新的目标函数值J^k+1，然后进入迭代控制模块；

迭代控制模块，用于控制迭代寻优模块的功能状态，并保存迭代优化的结果；首先判断是否满足收敛条件：

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

其中，J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值；若上式(11)成立，表明第k次迭代得到的目标函数值J^k与第k+1次迭代所得目标函数值J^k+1的误差绝对值不超过设定的收敛精度值ζ，则停止迭代优化计算，z^k+1就是最优决策向量z^*，J^k+1就是最优目标函数值J^*，将z^*、J^*以及相应迭代次数k+1保存到结果输出模块；若上式(11)不成立，则保存目标函数值J^k+1，取k＝k+1，然后返回迭代寻优模块进行新一轮的迭代求解。

2.如权利要求1所述的一种基于两点步长梯度法的工业过程动态优化系统，其特征在于：所述的上位机还包括：信息采集模块，用于设定采样时间，采集由现场智能仪表上传的工业过程对象的动态信息。

3.一种用如权利要求1所述的基于两点步长梯度法的工业过程动态优化系统实现的动态优化方法，其特征在于：所述的动态优化方法包括以下步骤：

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

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

u(t)＝0.5(u_max-u_min)×{cos[z(t)]+1}+u_min；                (1)

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

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

(3.1)在上位机的初始化模块，设置收敛精度值ζ，设置迭代优化的迭代次数k初始计数为0；

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

(3.3)对决策变量z(t)在步骤(3.2)所述各时间段上进行离散化，使用由n个分段常值组成的决策向量z来表示z(t)，并选取决策向量z的初始值z⁰，取为简单的常数向量；

(3.4)设置迭代搜索的初始步长α0＞0；

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

(4.1)数值求解状态方程组：

$\frac{\mathrm{dx}\left(t\right)}{\mathrm{dt}}=f[x\left(t\right),z\left(t\right),t],x_i\left(t_0\right)=x_{i0},(i=\mathrm{1,2},...,m)---\left(2\right)$

其中，f表示微分函数向量，x(t)为m个状态变量组成的状态向量，x_i(t)表示第i个状态变量，x_i0为状态变量x_i(t)在初始时刻t₀的初始值，采用龙格库塔法由初始值x_i通过正向积分求出状态变量x_i(t)在每个离散时刻的值x_i(t_j)，其中，i＝1，2，...，m，j＝1，2，...，n；

(4.2)数值求解协态方程组：

其中，

分别是给定的目标函数的非积分项和定积分项，λ_i(t)为第i个协态变量，λ(t)为m个协态变量组成的协态向量，λ_i(ft)为协态变量λ_i(t)在终端时刻tf的终端值，采用龙格库塔法由终端值λ_i(tf)通过反向积分求出协态变量λ_i(t)在每个离散时刻的值λ_i(t_j)，其中，i＝1，2，...，m；j＝n-1，n-2，...，1，0；

(4.3)由所得的状态向量和决策向量计算出目标函数值：

5)由状态向量和协态向量信息计算出迭代优化的搜索方向和步长，求解使目标函数值J更接近最优决策向量z^*；；实施迭代寻优的步骤如下，上标k均表示迭代次数，初始赋值为零：

(5.1)计算当前梯度g^k，即迭代优化的搜索方向，其中，x(t)和λ(t)分别为求解ODE得出的当前状态向量和协态向量，上标T表示向量或矩阵的转置：

$g^k={\{\frac{\∂\mathrm{\ψ}[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}+\mathrm{\λ}{\left(t\right)}^T\·\frac{\∂f[x\left(t\right),z\left(t\right),t]}{\∂z\left(t\right)}\}|}_{z\left(t\right)=z^k,t=t_j},j=\mathrm{0,1,2},...,n---\left(5\right)$

(5.2)保存当前迭代点z^k及梯度信息g^k；

(5.3)若k＝0，则搜索步长α^k取为初始值，即α^k＝α0，转步骤(5.5)；否则，依据两点步长梯度法，利用当前迭代点z^k和前一迭代点z^k-1的信息来确定步长因子l^k：

$l^k=\frac{{\left(s^{k-1}\right)}^T\·y^{k-1}}{{\left|\right|y^{k-1}\left|\right|}^2}---\left(6\right)$

其中，s^k-1表示当前迭代点z^k与前一迭代点z^k-1的误差，计算式为：

s^k-1＝z^k-z^k-1                    (7)

y^k-1表示当前迭代点z^k与前一迭代点z^k-1的梯度误差，计算式为：

y^k-1＝g^k-g^k-1                    (8)

(5.4)取最佳步长 ${\mathrm{\α}}^k=\min [\frac{\mathrm{\π}}{D\·\max \left(g^k\right)},l^k],---\left(9\right)$

其中D取区间[5，8]内的整数值；

(5.5)计算下一个迭代点：

z^k+1＝z^k-α^k·g^k                    (10)

并将z^k+1传给ODE求解模块计算出相应的目标函数值J^k+1；

6)判断收敛条件

是否满足，其中J^k和J^k+1分别表示第k次和第k+1次迭代计算得到的目标函数值，若不满足，则保存目标函数值J^k+1，取k＝k+1，再转入步骤4)，进行新一轮的迭代寻优；若满足，则终止迭代计算，z^k+1就是最优决策向量z^*，J^k+1就是最优目标函数值J^*，保存最优决策向量z^*、最优目标函数值J^*及迭代次数k+1到结果输出模块。

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

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

Images