面向建模误差PDF 形状与趋势拟合优度的动态过程优化建模

周 平 赵向志

实际复杂工业过程如高炉炼铁过程、磨矿过程、造纸制浆过程、污水处理过程等,通常涉及复杂的物理化学反应,具有多变量、强耦合、非线性、大滞后等综合复杂动态特性,利用传统的物理化学等机理方法时,难以建立精确的数学模型[1−4].近年来,随着大数据和人工智能技术的发展,对于难以进行机理建模的工业过程,甚至缺乏机理模型的复杂工业系统,数据驱动建模常被看作是一种非常有效的建模方法[1,4−5].目前,数据驱动建模主要采用人工智能技术或多元统计分析技术来描述过程输入与输出之间的复杂未知动态关系,在此基础上建立具有一定结构和适当模型参数的过程数据模型.由于数据模型输出与实际过程输出之间存在一定的偏差,为了最优化模型性能,通常需要采用相关算法来优化关于建模误差的性能指标,如均方根误差(Root mean square error,RMSE)、均方差(MSE)及平均绝对误差(Mean absolute error,MAE)等,以此获得满意的数据模型参数[6].

在实际工业过程建模时,单纯的RMSE、MSE、MAE 等性能指标均以最小化一维的统计建模误差均值为目标,并不能完全描述和刻画动态系统建模过程的随机性和不确定性[6].在复杂工业过程中,外部不确定因素和随机动态干扰往往具有非高斯特性,将其假设为高斯分布特征时,难以获得满意的建模效果[7].另一方面,复杂工业过程中,其建模误差通常是一个未知的随机变量,因而建模误差在时空尺度上的二维概率密度函数(Probability density function,PDF)形状分布能够包含动态系统建模误差的所有分布和统计信息.因此,误差PDF 形状优化的思想受到了越来越广泛的关注,并逐渐用于复杂工业过程的数据建模与控制中,如以非高斯动态系统PDF 形状为目标的随机分布控制已成为解决随机动态系统控制的非常有效的方法,得到了广泛的应用[7−8].近年,面向建模误差PDF 形状优化的数据建模方法也逐渐引起重视.文献[6]中,作为有界随机分布系统建模与控制方法的扩展,Zhou等将输出PDF 控制或随机分布控制思想引入到过程建模中.通过优化建模误差PDF 形状,进而对模型参数求解,使得实际建模误差PDF 形状接近设定的期望PDF 形状[6].该方法不仅可以获得较为满意的建模效果,而且一定程度降低建模过程中的随机性和不确定性.此外,文献[9]通过优化建模误差PDF 形状,建立了选矿过程精矿品位的最小支持向量机模型.而文献[10]通过优化建模误差PDF 形状,间接对模型参数进行调节,建立了间歇过程的模糊神经网络模型.

上述方法均是期望实际模型的建模误差PDF形状更好地跟踪期望的高斯分布形状,以此建立具有最优参数的过程数据模型.然而,不管是常规建模方法的误差RMSE 指标,还是上述改进方法提到的建模误差PDF 指标,均仅仅体现过程模型输出与实际输出之间的误差大小情况,难以衡量模型输出与实际动态过程输出之间拟合趋势是否一致.实际上,实际工业动态系统中,过程输出变化趋势的估计和预测,对于基于模型的预测控制、生产过程运行态势的把握与调控等诸多工程应用,都具有十分重要的作用.因此,在动态系统建模时,除了需要优化建模误差的PDF 形状,同时也需要考虑建模输出与样本数据之间拟合趋势最接近,即曲线拟合动态变化趋势的相似度最大[11].

针对上述动态系统建模的实际需求和现有方法的不足,本文以小波神经网络(Wavelet neural network,WNN)[12−13]数据建模为例,提出一种新型的面向建模误差PDF 形状与趋势拟合优度(相似度)多目标优化的动态系统数据建模方法.所提方法不仅引入二维尺度的PDF指标来对动态建模误差在时间和空间进行全面刻画,同时引入拟合优度(相似度)指标[11,14]刻画动态系统数据建模的拟合趋势.通过采用核密度估计(Kernel density estimation,KDE)[15−17]技术对实际建模误差PDF 形状进行估计,以及采用NSGA-II 算法[18]对建模误差PDF 形状的偏差以及拟合优度指标进行多目标优化,从而建立具有最优模型参数的WNN 模型.数值仿真以及污水处理过程[19−20]数据验证表明所提方法的实际建模误差PDF 能够更好地逼近设定的期望PDF,并且模型输出与样本数据拟合趋势接近.

1 建模思路与策略

1.1 WNN 简介

小波神经网络是结合小波分析与神经网络的一种前馈型网络.WNN 用小波函数代替传统Sigmoid 函数作为激励函数,通过仿射变换建立起小波变换与网络参数之间的连接,能以任意精度对函数进行逼近[12−13].如图1 所示,WNN 通常采用三层网络结构,其中:x1,x2,···,xM为WNN 的输入变量,y1,y2,···,yN为输出变量,ωI,ji和ωH,lj分别是输入层到隐含层的连接权值以及隐含层到输出层的连接权值,θH,j,θO,l分别为隐含层和输出层节点阈值,aj和bj分别为小波基函数的伸缩因子和平移因子.

图1 WNN 结构图Fig.1 Structure diagram of WNN

WNN 隐层激励函数ϕ(t) 通常采用如下的Morlet母小波函数[12]:

输出层激励函数f(t) 则采用Sigmoid 函数,即:

此外,定义三层WNN 的输入层、隐含层、输出层节点数分别为M,n和N,则隐含层第j个节点的输入netj和输出zj分别为:

式中,bj为小波基函数φ(t)的平移因子,aj为小波基函数的伸缩因子,φ(t) 为Morlet 母小波函数,θH,j为隐含层节点的阈值,则WNN 的最终输出为:

式中,θO,l为输出层节点的阈值,函数f(·) 为输出层 的激励函数.

1.2 建模策略

对于常规WNN 等现有多数建模方法,通常采用如下均方根误差(RMSE)、均方差(MSE)、平均绝对误差(MAE)等单一的误差性能指标,通过性能指标数值大小来评价建模精度.

然而,式(6)所示传统性能指标是从建模误差的均值角度评价模型精度,并不能全面描述动态系统建模误差在时空尺度上的随机特性.此外,对于时序相关动态系统建模,运行数据拟合趋势的估计对于建模效果有很大影响,并且更有实际意义.而式(6)所示常规建模性能评价指标仅希望建模输出与实际数据之间偏差最小,却难以描述动态系统的拟合趋势好坏.

为了解决上述问题,本文以WNN 智能建模为基础,通过引入建模误差概率密度函数(PDF)指标从时空二维角度对建模误差进行全面刻画,以及引入拟合优度指标对动态系统数据建模的拟合趋势进行相似性评估,从而提出图2 所示的面向建模误差PDF 形状与趋势拟合优度的动态系统优化建模方法,具体如下:

图2 面向建模误差PDF 形状与趋势拟合优度的优化建模策略Fig.2 Optimized modeling strategy towards modeling error PDF shape and goodness of fit

1)首先,构建动态系统数据建模的实际建模误差PDF 与期望建模误差PDF 的偏差平方积分作为多目标优化计算的第一个评价指标,如下所示:

式中,Γ(e) 和 Γtarget(e) 分别为实际建模误差PDF和期望建模误差PDF.本文实际建模误差PDF 是采用核密度估计技术对所建立数据模型的建模误差序列进行求解获得,而期望建模误差PDF 是设置的一个较为理想的(即均值为0、方差尽量小)高斯分布形状的PDF,e为建模误差PDF 的自变量.

2)其次,引入式(8)所示的拟合优度指标ρAB[11,14]对动态系统数据建模的动态拟合趋势进行相似性评估,然后构建式(9)所示关于拟合优度的性能指标作为第二个评价指标.

式中,A,B为两个数据矩阵,分别为数据矩阵A,B的均值.事实上,ρAB是衡量数据矩阵A 和B之间近似程度的量,|ρAB|→1 表示数据矩阵A 和B 之间相关性很强,而|ρAB|→0 意味着数据矩阵A 和B 之间相关性较弱.由于本文要衡量建模输出与实际输出的时序相关数据之间的动态拟合趋势,所以本文式(8)中A 和B 分别表示小波神经网络模型输出和实际过程输出所构成的时序相关数据矩阵.

3)最后,分别将式(8)和式(9)作为数据建模的综合性能评价指标的适应度函数,采用运算速度快、解集收敛性好的NSGA-II 算法[11]来获得WNN模型的最优参数集 [ωI,ji,ωH,lj,aj,bj,θH,j,θO,l].

2 建模算法

所提方法的具体建模算法包括如下几个过程:首先,采用式 (1)∼(5) 所示算法构建初始的WNN数据模型,通过比较WNN 数据模型输出与过程输出或者相应的实际值,可以得到特定时间内的建模误差序列.然后,采用第3.1 节的核密度估计技术对实际建模误差PDF 进行计算.最后,采用NSGAII 算法优化式(8)和式(9)所示的多目标性能指标,获得同时具有较好建模误差PDF 形状与拟合优度值的多组WNN 模型参数集解.

2.1 建模误差PDF 核密度估计

核密度估计(Kernel density estimation,KDE)是由Parzen 提出的一种非参数估计方法[15−17],用于求解给定随机变量数据集合分布的概率密度函数.假设xi∈R,i=1,···,n为独立同分布的随机变量数据集,其所服从的分布密度函数为f(x),x ∈R,则f(x)的核密度估计定义如下:

式中,窗宽hp是一个给定的正数,ϕ(x) 为核函数,n为样本数.

对于所提建模方法,采用KDE 对建模误差PDF进行估计,可以得到估计的建模误差概率密度函数Γe为:

式中,K为建模误差样本数目,通过设置期望的建模误差PDF,可以构造式(7)所示的性能指标.式(11)所示基于KDE 的WNN 实际建模误差PDF估计求解步骤如下:

1)选择核函数:在估计随机变量未知概率密度函数时,常用的核函数有高斯核函数、矩形窗核函数、Epanechnikov 核函数等.核函数的不同选择在KDE 中不敏感,当样本数据很大时,对核函数密度估计的结果影响不大.本文选取高斯核函数,其表式如下:

2)选择窗宽:窗宽hp的选择对核函数的密度估计起着局部光滑的作用,如果hp过大会使模型误差PDF 形状很光滑,使其主要部分的某些特征(如多峰性)被掩盖起来,从而增加估计量的偏差;而若hp过小,则整个密度函数表现粗糙.本文基于正态参照规则方法[17]进行窗宽选择,假设建模误差服从正态分布,则窗宽hp设置为hp=1.06σK−1/5其中σ由 min{s,Q/1.34}估计,s 表示样本标准差,Q为四分位数间距.

3)求解模型误差PDF:根据步骤1)和2)选择合适窗函数和窗宽参数,然后代入式(11),可以得到WNN 建模误差PDF 函数的估计值为:

2.2 模型参数多目标优化求解

式(1)～(5)所示基本WNN 数据模型的参数主要包括:输入层连接权值ωI,ji、隐含层连接权值ωH,lj、隐含层阈值θH,j、输出层阈值θO,l、小波基函数的伸缩因子aj以及平移因子bj.这些参数的取值直接决定了WNN 数据模型的性能,因而基于前述构建的多目标建模性能指标,采用NSGA-II 算法对模型参数进行优化,步骤如下:

1)网络参数的编码.将WNN 模型参数集θWNN=[ωI,ji,ωH,lj,aj,bj,θH,j,θO,l]与每条染色体相对应,即对WNN 模型参数进行如下形式的编码:

式中,染色体基因数为S=(M+3)n+(n+1)N,P=[S1,S2,···,Si,···,SQ]T表示包含Q条染色体的初始种群.

2)个体适应度计算.每条染色体中的各个基因分别代表WNN 的各个参数,将第t代种群中第h条染色体上的各个基因代入下式的第t代、第h个个体多适应度函数中:

3)选择算子.根据非支配排序结果,选择非支配排序中支配层较低的个体.如果有多个个体在同一支配层,从种群多样性角度考虑,选择拥挤度距离较大的个体.

4) 模拟二进制交叉.基于实数编码,交叉后代为父代的线性组合,即

5) 多项式变异.二进制交叉后,进行多项式变异,变异后的个体为:

采用NSGA-II 算法优化WNN 模型参数时,每个待优化的参数对应染色体上的一个基因.在遗传算法中,适应度函数的选择决定着遗传优化的精度和收敛速度.描述个体性能的指标主要通过适应度函数值体现,依据适应度值的大小对个体进行优胜劣汰.本文多目标适应度函数为实际建模误差PDF与期望建模误差PDF 之间的二维偏差平方和以及趋势拟合优度的倒数,并通过基因之间的选择、二进制交叉、变异产生最优个体即最优模型参数.

3 数值仿真及工业数据验证

3.1 数值仿真

为了验证所提方法的有效性和优越性,首先使用下述两输入一输出非线性动态系统进行数值验证:

式中,y(0)=0.1,u(k) 为在区间(0,1)内服从均匀分布的随机序列,ω(k)为通过参数σ描述的、服从瑞利分布的非高斯随机干扰序列.针对以上非线性系统,利用提出的建模方法进行建模,所要建立的WNN 数据模型可以表示为:

假设ω为随机产生、服从瑞利分布且参数为0.2的非高斯干扰.WNN 隐层节点数选择为6,迭代优化步长r为0.003.采用NSGA-II 算法对WNN 模型参数进行寻优时,交叉分布指数ηc=20,变异分布指数ηm=20,优化变量的上限与下限分别设定为1 和−1,交叉率和突变率分别设为0.9 和0.1.

建模后,得到60 组Pareto 前沿解进化过程如图3 所示,而图4 为60 组多目标优化解对应的拟合优度变化曲线,这里将所提方法与常规WNN 方法以及近年文献[6]中提出的面向建模误差PDF优化的WNN 方法进行比较.由于文献[6]是采用梯度下降方法来优化WNN 模型的建模误差PDF,因而本文将其称为GD-WNN.从图4 可以看出,采用所提建模方法可以获得具有较大动态变化趋势拟合优度的一组解集,这些解对应的拟合优度均远好于常规WNN 方法以及GD-WNN 方法.所提方法得到的解对应的拟合优度指标最高达到0.96,而文献[6]中方法的拟合优度指标仅为0.83,以及常规WNN 方法的拟合优度指标甚至仅为0.75.并且图4还可以看出所提方法有39 组解的拟合优化度指标好于文献[6]中方法得到的拟合优度值.图5 是图3中1 到60 号解对应的建模误差PDF 曲线变化图,可以看出60 号解对应模型的建模误差PDF 最好.图6 为多目标优化后30 号解与60 号解的建模误差PDF 曲线与其他两种现有方法的建模误差PDF曲线的对比图,本文设置的期望PDF 为均值为0方差为0.25 的高斯型概率密度函数.

图3 Pareto 前沿解进化过程Fig.3 Evolution process of Pareto front

图4 不同优化解对应的拟合优度值变化曲线Fig.4 Change curve of goodness of fit corresponding to different optimization solutions

图5 不同优化解对应的建模误差PDF 变化曲面Fig.5 PDF changing surface corresponding to different optimization solutions

图6 不同方法建模误差PDF 比较Fig.6 Comparison of modeling error PDF with different methods

可以看出,本文方法获得的非最优30 号解对应的建模误差PDF 曲线也要远好于常规WNN 方法和GD-WNN 方法.所提方法得到的建模误差PDF 较高且较窄,与其他方法相比方差更小,即模型的随机性和不确定性更小,这也表明所提方法的有效性和优越性.图7 和图8 分别是本文方法中非最优的30 号解对应的建模效果和新样本测试效果,可以看出所提方法得到非最优解对应的模型不论是建模和新样本测试均好于其他两种方法.

图7 所提方法30 号优化解对应的建模效果Fig.7 Modeling result corresponding to the 30th optimization solution of the proposed method

图8 所提方法30 号优化解对应的测试效果Fig.8 Testing result corresponding to the 30th optimization solution of the proposed method

3.2 活性污泥污水处理过程生化反应池出水COD含量建模

目前,城市污水处理广泛采用活性污泥法[3,19−22].图9 为典型活性污泥污水处理的工艺流程图,主要包括三个级别的处理过程.一级处理主要进行物理反应,除去原生污水中的悬浮固体.二级处理过程包括曝气池及二沉池处理两个部分.曝气池是污水处理的核心部分,主要进行微生物自身的代谢活动,从而达到对污水中有机污染物如氮、磷的去除以及有氧生物的降解.

图9 典型活性污泥法污水处理过程工艺流程图Fig.9 Flow chart of a typical activated sludge wastewater treatment process

经过曝气池处理后的污水流入二沉池进行固液分离,上层是澄清的液体,下层的污泥一部分回流至曝气池,以维持曝气池内的污泥浓度,另一部分污泥排出系统.经过二级处理后的污水进入三级处理过程,通过加入药剂进而得到达标的出水.判断水质是否达标主要通过水质参数进行衡量.在众多水质参数指标中,出水COD 含量不仅代表污水中含有的有机物的量,同时还包括污水中还原性无机物被氧化时所消耗的氧气量.COD 数值越小,在氧化过程中氧气的消耗量就越少,即水体中有机物的量越少.因此,该指标能够反映有机污染物受纳的程度,是非常重要的出水水质指标.虽然目前有许多COD 含量的在线检测仪,但是都存在检测周期长、价格昂贵的问题.所以,通过基于数据的智能建模技术来预测曝气池出水COD 含量对于判断出水水质是否达标具有重要意义.

影响污水出水COD 含量的参数较多,包括:入水流量(Q)、入水化学需氧量(COD)、溶解氧浓度(DO)、污泥浓度(MLSS)、悬浮固体浓度(SS)、出水PH 值等.为此,根据过程机理分析,确定出水COD 预测建模的输入变量为:入水COD、入水流量(Q)和污泥浓度(MLSS).为消除变量间的量纲影响,建模所用训练与测试数据都归一化处理.设定WNN 隐层节点数为6,迭代优化步长为0.003,得到的100 组多目标优化Pareto 前沿解如图10 所示.图11 和图12 分别为100 组多目标优化解对应的拟合优度变化曲线和建模误差PDF 变化曲线.同样,这里也将所提方法与常规WNN 方法以及文献[6]中的GD-WNN 方法进行比较.可以看出,采用所提建模方法可以获得具有较大动态趋势拟合优度的70 余组解,这些解对应的拟合优度均远好于常规WNN方法以及文献[6]中的GD-WNN 方法,并且所提方法所得解的最优拟合优度已非常接近1.从图12所有多目标优化解的建模误差PDF 变化曲线可以看出,从第1 号解到第100 号解,建模误差PDF 的形状越来越窄而尖,并且越来越接近设定的理想PDF 形状,即模型的随机性和不确定性很小.

图10 COD 含量建模Pareto 前沿进化过程Fig.10 Pareto front evolution process of COD content modeling

图11 不同优化解对应的拟合优度值变化曲线Fig.11 Change curve of goodness of fit corresponding to different optimization solutions

图12 不同优化解对应的COD 含量建模误差PDF 变化曲面Fig.12 PDF changing surface of COD content modeling error corresponding to different optimization solutions

图13 为50 号解与100 号解的建模误差PDF曲线与其它方法的建模误差PDF 曲线对比图,这里设置的期望PDF 为均值为0 方差为0.3 的高斯型概率密度函数.图中可以看出,本文方法获得的非最优50 号解对应的建模误差PDF 形状为窄而高的形状,仍然远好于常规WNN 方法和GDWNN 方法.图14 和图15 分别是本文方法所得非最优50 号解的建模效果以及对新样本的测试效果,虽然是选取的本文方法的非最优解,但是从图中可以看出非最优解的建模和新样本测试效果均好于其他两种对比方法.

图13 不同方法COD 含量建模误差PDF 比较Fig.13 PDF comparison of COD content modeling error with different methods

图14 所提方法50 号优化解对应的COD 含量建模效果Fig.14 Modeling result of COD content corresponding to the 50th optimization solution of the proposed method

图15 所提方法50 号优化解对应的COD 含量测试效果Fig.15 Testing result of COD content corresponding to the 50th optimization solution of the proposed method

4 结语

基于误差最小的数据驱动工业系统建模时,通常基于单一的RMSE 等一维性能指标.但是RMSE等时间维度的一维性能指标并不能充分体现动态系统建模的随机性和不确定性.同时,对于传统动态系统建模方法,并没有考虑模型输出和动态系统实际输出之间的拟合趋势.为此,本文基于数据驱动小波神经网络智能建模、多目标参数优化以及核密度估计技术,提出综合考虑建模误差PDF 形状与趋势拟合优度的动态系统优化建模方法.其中多目标参数优化的性能指标分别为实际建模误差PDF与期望建模误差PDF 之间二维偏差平方、趋势拟合优度.仿真实验以及污水处理过程数据验证表明:相比于对比的两种现有建模方法,所提方法不仅具有更好的建模精度和泛化能力,还可控制建模误差的空间分布状态,使得所提方法的建模误差PDF比传统建模方法的建模误差PDF 更高、更窄,即模型中含有的随机性和不确定性更小.此外,所提方法还可以获得一大类具有建模误差PDF 形状接近期望分布形状,且模型输出与实际输出的趋势拟合优度值较大的数据模型参数解,因而具有更好的实用性.