对于无稳态非线性系统,局部线性模型无法通过稳态点上数据辨识确定,需进行寻优.双操作变量LPV模型的线性权重形式已在文献[12]中进行推导并给出计算公式,在此不做赘述.线性权重中无待定参数,可直接由当前操作变量值w₁(t)、w₂(t)确定,所以,辨识无稳态对象的线性权重LPV模型仅需对局部线性模型进行寻优.在式(2)和(3)中,记θ=[a^1,1₁ …a^1,1_r,b^1,1,1₁ …b^n₁,n₂,m_r]为待定参数序列,N为采样数据长度,则θ可通过极小化LPV模型的输出误差ε_OE(t)来求解:

(＾overθ)=argmin_θ ∑^N_t=r+1ε_OE(t)²=

argmin_θ∑^N_t=r+1[y(t)-(^overy)(t)]²,(4)

其中,y(t)为系统实际输出,将式(1)和(2)代入(^overy)(t)可得

y(t)-(^overy)(t)=y(t)-∑n1i=1∑n2j=1αi,j(wi1,wj2)(^overy)i,j(t)=y(t)-∑n1i=1∑n2j=1αi,j(wi1,wj2)1/(Ai,j(q))

[Bi,j1(q)u1(t)+…+Bi,jm(q)um(t)]=y(t)-∑n1i=1∑n2j=1αi,j(wi1,wj2)

([(b^i,j,1₁ q^-1+…+b^i,j,1_r q^-r)u₁(t)+…+(b^i,j,m₁ q^-1+…+b^i,j,m_r q^-r)u_m(t)])/(1+a^i,j₁ q^-1+…+a^i,j_r q^-r).(5)

从式(5)可看出,由于线性权重没有待定参数,所有待定参数均集中在局部线性模型中,并且,a^i,j_k(k=1,2,…,r)与ε_OE呈非线性关系,必须采用非线性数值计算方法求解式(4)的优化命题.在此采用高斯牛顿法^[13]和最小二乘法相结合的方法进行迭代求解.在每一次迭代中,令式(4)中的ε_OE(t)在局部最优解附近进行近似泰勒展开,通过求偏导实现待定参数和输出误差之间关系的线性化,随后进行最小二乘解的计算确定下一组待定参数,直至输出误差小于某阈值.具体算法如下:假设(＾overθ)^k_Model 是第k次迭代获得的模型参数序列[a^1,1(k)₁ …a^1,1(k)_r,b^1,1,1(k)₁ …b^{n₁,n₂,m(k)}_r],且接近全局最优解.忽略泰勒展开式中的高次项,则输出误差可近似为

ε_OE(t,θ_Model)≈ε(t,(＾overθ)^k_Model)+

(ε(t,θ_Model))/(θ^T_Model)〖JB>2|〗_{θ=(＾overθ)^k}(θ_Model-(＾overθ)^k_Model)=

ε(t,(＾overθ)^k_Model)+φ^k(t)(θ_Model-(＾overθ)^k_Model)=

-[φ^k(t)(＾overθ)^k_Model-ε(t,(＾overθ)^k_Model)]+φ^k(t)θ_Model,(6)

其中,

φ^k(t)=(ε(t,θ_Model))/(θ^T_Model)〖JB>2|〗_{θModel=(＾overθ)^kModel}=

[(ε(t,θModel))/(a1,11 )〖JB>2|〗(^overa)1,1(k)1…(ε(t,θModel))/(an1,n2r )〖JB>2|〗(^overa)n1,n2(k)r(ε(t,θModel))/(b1,1,11 )〖JB>2|〗b1,1,1(k)1…(ε(t,θModel))/(bn1,n2,mr )〖JB>2|〗bn1,n2,m(k)r]=

[α^1,1(w¹₁,w¹₂)(B^1,1(k)₁(q)u₁(t-1)+…+B^1,1(k)_m(q)u_m(t-1))/(A^1,1(k)(q)²)



αn1,n2(wn11,wn22)(Bn1,n2(k)1(q)u1(t-r)+…+Bn1,n2(k)m(q)um(t-r))/(An1,n2(k)(q)2)

-α1,1(wn11,wn22)1/(A1,1(k)(q))u1(t-1)



-αn1,n2(wn11,wn22)1/(An1,n2(k)(q))um(t-r)]T.(7)

处理后,a^i,j_k(k=1,2,…,r)与ε_OE(t)呈线性关系,可直接用最小二乘法求解.将式(6)～(7)代入式(4),结合最小二乘解式,推导得待定参数序列的下一步估计值为

(＾overθ)^k+1_Model=

[∑^N_t=r+1[φ^k(t)]^Tφ^k(t)]^-1∑^N_t=r+1[φ^k(t)]^T

[φ^k(t)(＾overθ)^k_Model-ε(t,(＾overθ)^k_Model)]=(＾overθ)^k_Model-

[∑^N_t=r+1[φ^k(t)]^Tφ^k(t)]^-1∑^N_t=r+1[φ^k(t)]^Tε(t,(＾overθ)^k_Model).(8)

数次迭代之后,当输出误差小于某个设定的阈值时,即可停止迭代.

上述算法的有效性很大程度取决于待定参数的初值.若初值能够尽可能接近局部最优解,则式(6)的近似误差较小,能够在有限次数内结束迭代.若初值与最优解偏差较大,则难以保证算法的收敛性.在实际应用时,为提高初始局部模型的精度,在无稳态阶段的情况下,应尽量截取所选工作点附近的数据片段进行辨识.常规的线性模型辨识手段均可用于局部模型的初始建模.

双操作变量的高斯权重可记为

α^i,j(wⁱ₁,w^j₂)=exp[-1/2(((w₁(t)-wⁱ₁)/(σ¹_i,j ))²+

((w₂(t)-w^j₂)/(σ²_i,j ))²)],(9)

其中σ¹_i,j 和σ²_i,j是分别对应2个工作点变量的待定分离系数.

与线性权重不同,高斯权重中含有待定参数,故若继续采用式(4)进行参数寻优,则待定参数包含2部分,一部分在权重函数中,另外一部分在局部模型中,且高斯权重的待定参数在指数函数的分母上,显著加大了求解难度.针对此问题,可采用Narendra-Gallman算法^[14]求解.该方法的基本思想是,根据待定参数与优化命题的关系,将待定参数分为线性部分和非线性部分,每次迭代分为2步:先固定非线性部分,优化线性部分; 再固定线性部分,优化非线性部分.已有文献证明该方法是处理复杂线性、非线性混合优化命题的有效方法,且具有鲁棒性强、收敛性强的优点.本研究将其应用于无稳态系统的高斯权重LPV模型辨识中,具体步骤如下:

1)确定待定参数初值、迭代结束阈值.

待定参数包括各工作点处的高斯函数分离系数,以及局部线性模型的分子分母系数.其中,根据高斯函数图像特性以及LPV模型的结构特点,分离系数初值设为相邻操作变量值的1/3,如σ¹_1,1 的初值设为(w²₁-w¹₁)/3; σ²_1,1 的初值设为(w²₂-w¹₂)/3.对各组局部线性模型的初值,则选择对应工作点附近的数据进行辨识.设置足够小的正数ε作为迭代结束阈值,当LPV模型输出误差小于ε时停止迭代.设置迭代计数器初值n_Count=1.

2)固定优化命题的非线性部分,优化线性部分.

同样采用式(4)作为优化命题.所有待定参数中,局部模型分子上的参数[b^1,1,1₁ …b^n₁,n₂,m_r]与ε_OE呈线性关系,权重参数[σ¹_1,1,σ²_1,1,…σ¹_n1,n2,σ²_n1,n2]以及局部模型分母参数[a^1,1₁ …a^1,1_r]与ε_OE呈非线性关系.按照Narendra-Gallman算法思路,保持非线性部分不变,仅优化线性部分.记当前局部模型分子参数序列为(＾overθ)^(k)_B=[b^1,1,1(k)₁,…,b^{n₁,n₂,m(k)}_r],直接采用最小二乘法进行优化,得:

(＾overθ)^(k+1)_B=[Φ^TΦ]^-1Φ^TY,(10)

其中,

Y=[y(r+1)

y(r+2)



y(N)],Φ=[(α1,1(w11,w12)u1(r))/(A1,1(k))(α1,1(w11,w12)u1(r-1))/(A1,1(k))…(αn1,n2(wn11,wn22)um(1))/(An1,n2(k))

 

(α1,1(w11,w12)u1(N-1))/(A1,1(k))(α1,1(w11,w12)u1(N-2))/(A1,1(k))…(αn1,n2(wn11,wn22)um(N-r))/(An1,n2(k))].

3)固定优化命题的线性部分,优化非线性部分.

从理论角度,应保持2)计算出的(＾overθ)^(k+1)_B不变,同时优化剩余的高斯权重参数和局部模型分母参数.然而,经过多次实验表明,由于没有对局部模型分母参数做约束,寻优过程中难以保证局部模型的稳定性,从而影响了插值后LPV模型的稳定性.为确保各局部模型的稳定性,在此仅优化权重参数.很多数值计算方法都可以用来优化权重参数,如最速下降法等.在本步骤中进行内层迭代,每次迭代通过对权重参数序列求偏导,计算负梯度向量并由此确定搜索方向.

4)用式(4)计算LPV模型输出误差V^N_OE.若V^N_OE<ε则停止迭代; 反之令n_Count=n_Count+1,返回2).

值得注意的是,若系统仅含有一个操作变量,1.1节和1.2节的所有算法和步骤依然适用.辨识时仅需删去权重函数的其中一维,减少对应的待定参数即可.

为方便建模效果的对比,辨识所用原始数据以及选择的输入和输出变量均与文献[15]相同.即在不影响正常工业生产的前提下,在系统输入端分别叠加小幅值广义二进制噪声进行4.5 d连续不间断的数据采集,采样周期为1 min.通过对运行数据的分析以及与现场工程师交流得知,出口蒸汽流量的波动源于白天和夜晚产量的不同,并且是造成系统非线性的主要因素,系统的其他变量均随着蒸汽流量的变化而做出相应调整,将其作为唯一的操作变量.图2为蒸汽流量原始采样数据,由于数据历时4.5 d,从图中可看出处于白天的高负荷和处于夜晚的低负荷阶段交替进行,各出现5次.根据低、中、高负荷选择如下3个典型工作点:

工作点1:蒸汽流量=20 m³/h;

工作点2:蒸汽流量=30 m³/h;

工作点3:蒸汽流量=40 m³/h.

图2 蒸汽流量采样数据
Fig.2 Sampling data of the steam flow

以2.1节辨识获得的线性模型,作为式(1)中LPV模型中局部模型初值分别用1.1和1.2节所提方法,构建基于线性权重和高斯权重的LPV模型.文献[15]的结果中,不管采用何种权重,5个输出中蒸汽压力、蒸汽温度和炉膛温度的拟合效果很差,在此着重对这3个输出进行对比.由于2种权重值略有不同,寻优后的局部线性模型也不尽相同,但增益大小总体相差不大.基于线性权重的LPV模型中,所有待定参数均集中在局部模型中,由于共有5个输入且为3阶模型,故每个输出对应的LPV模型待定参数总数=输入个数×模型分子参数×工作点数=5×3×3=45; 以蒸汽压力为例,寻优共经过132次迭代结束.基于高斯权重的LPV模型中,待定参数分布于权重部分和局部模型中,故每个输出对应的LPV模型待定参数总数比基于线性权重的LPV模型多3个,共48个; 蒸汽压力的寻优共经过356次迭代结束.虽然待定参数仅相差3,但由于Narendra-Gallman算法在每次迭代中需要进行线性部分和非线性部分的交替寻优,在优化非线性部分时还引入了最速下降法,最终寻优时间远大于基于线性权重的LPV模型.图3为2种LPV模型在3个输出的拟合曲线.从图中可以看出,2种权重的LPV模型输出均能够较好地拟合实际值.为使比较结果更为直观,表1对比了文献[15]中带稳态LPV模型和本文中无稳态LPV模型的最优匹配率R_BFT:

R_BFT=100%×max(1-(=y(t)-(^overy)(t)=₂)/(=y(t)-(-overy)=₂),0),(11)

其中,(-overy)为系统输出均值.

综合考虑图3和表1的对比结果,用无稳态LPV模型对循环硫化床锅炉进行建模,获得了比带稳态LPV模型更好的拟合效果和最优匹配率,3个输出的最优匹配率均有了显著的提高.这是因为该锅炉运行过程中,造成非线性的主要因素,即蒸汽流量的波动无明显稳态,尤其在高负荷阶段,蒸汽流量在35～50 m³/h范围内持续并快速变化(如图2所示),系统运行模型也随之发生改变.采用带稳态LPV模型对其建模时,分别截取了高、低负荷附近的数据辨识局部模型并直接参与插值,局部模型本身带有较大的误差,后续虽调整模型权重仍不足以消除这些误差,导致最终拟合效果不尽如人意; 而采用无稳态LPV模型时,辨识过程中不仅调整了局部模型权重,同时也对局部线性模型参数做了优化,使局部模型更加精确,使得整体LPV模型获得更好的拟合效果.在均采用无稳态LPV模型的前提下,由于高斯函数和线性函数的图像形式和权重值均相似,采用高斯权重并没有比采用线性权重获得明显更好的效果.由于在辨识过程中,高斯权重LPV模型的寻优难度和耗费时间远大于线性权重LPV模型,故针对此系统,采用基于线性权重的无稳态LPV模型建模是较为合适的选择.

图3 实际系统与LPV模型输出对比
Fig.3 Comparison of outputs of real process and LPV model

表1 带稳态与无稳态LPV模型最优匹配率对比
Tab.1 Comparison the optimal matching rates of LPV model between steady state and non-steady state %