2.4m风洞是中国空气动力研究与发展中心所属的一座半回流暂冲式跨声速风洞[1].对风洞系统建模的主要目的是为风洞控制系统的设计提供精确可靠的数学模型,由于风洞内部结构复杂,各部分相互耦合,并涉及三维非定常流动的复杂气动问题,很难建立精确的机理模型.
近年来,利用运行数据对风洞系统建模的研究取得了一些成果,主要方法有:基于特征子集的神经网络融合(ENN)模型[2],基于BP网络的NARMAX(BP-NARMAX)模型[3],以及根据输入输出数据获得的线性降阶模型[4].这几种方法实现了对风洞系统的建模,但BP神经网络模型在训练过程中存在收敛速度慢、容易陷入局部最优等固有缺陷,且模型精度容易受到训练样本的影响,很难保证风洞控制系统的实时性和稳定性要求;线性降阶模型虽然能为控制器的设计提供便利,但是降阶模型无法精确描述风洞这样复杂系统的动态特性.另外,在某些特种试验中,采用运行数据建模难以满足充分激励条件[5],现有智能及现代辨识法难以有效应用.
本文针对样本数较少,无法满足持续激励的风洞试验工况,将风洞流场近似为线性系统[6-7],通过实验测取风洞系统的阶跃响应.采用三次样条插值法对数据进行滤波,再从耦合的数据中得到单个通道的阶跃响应数据,采用改进的面积法辨识出各个通道的传递函数模型,并通过对不同阶次模型输出均方根误差(RMSE)的比较,实现了各个通道最佳模型阶次的选择.
1 风洞系统模型风洞试验阶段是通过调节主排气阀开度和栅指位移控制j1总压(y1)和静压(y2).主排气阀位移和栅指开度与输出总压和静压之间存在着时滞和耦合.检验流场品质的两个重要指标是稳定段总压和试验段马赫数(Ma).马赫数是由总压和静压通过式(1)得到.
(1) |
风洞系统模型结构可简化成图 1所示形式.
其中:u1,u2分别为风洞系统主排气阀和栅指的控制量输入;y1,y2分别为输出总压和静压;g11(s),g21(s),g12(s)和g22(s)分别为各个通道的传递函数.由图 1可得到风洞流场模型为
(2) |
(3) |
各个通道的传递函数定义为
(4) |
(5) |
(5) |
(7) |
其中:K11,K12,K21,K22为各个通道的静态增益;L11,L12,L21,L22为各个通道的延迟时间;a1,1,…,a1,n1-1,a1,n1;a2,1,…,a2,n2-1,a2,a2;a3,1,…,a3,n3-1,a3,n3;a4,1,…,a4,n4-1,a4,n4;b1,m1-1,b1,m1;b2,1,…,b2,m2-1,b2,m2;b3,1,…,b3,m3-1,b3,m3;b4,1,…,b4,m4-1,b4,m4分别为模型参数.本文主要工作就是通过阶跃响应数据,确定各个通道传递函数最佳阶次及相应的模型参数.
2 基于面积法的模型参数辨识利用阶跃响应曲线确定典型工业过程传递函数的方法很多[8-10],其中面积法有广泛的工程应用价值,更适用于风洞系统的建模.经典面积法辨识模型参数的过程如下,假设过程的传递函数为
(8) |
其中:m,n(n≥m)为模型阶次;b1,…,bm-1,bm;a1,…,an-1,an为模型参数;L为延迟时间;K为增益.
首先,将得到的阶跃响应数据转化成无因次的形式,如图 2所示.
其中:u(t)为控制量阶跃值;h(t)为阶跃响应数据;u*(t),h*(t)为经过无因次转化后的控制量阶跃值和阶跃响应数据.则有
(9) |
令
(10) |
式中,
(11) |
则[1-h*(t)]的Laplace变换为
(12) |
则图 2中一阶面积A1为
(13) |
再令
(14) |
并定义二阶面积A2为
(15) |
依此类推,i阶面积Ai为
(16) |
其中
(17) |
利用
(18) |
可得
(19) |
其中,
(20) |
由式(12),式(19)及ci=Ai(i=1,2,…),得
(21) |
则有
(22) |
根据式(8)及ci=Ai(i=1,2,…),可得
(23) |
比较式(23)两边各次幂的系数可得到以下模型参数对应关系矩阵(24)和(25).
(24) |
(25) |
为获得所需的阶跃响应数据,在风洞流场达到稳态后,先保持栅指位移不变,对主排气阀给定阶跃值,待系统进入新的稳态后,保持主排气阀开度不变,对栅指给定阶跃值,并采集总压和静压数据,经过相应计算获得单个控制量变化分别对总压和静压的影响.本文选取总压110kPa,马赫数为0.578工况下的阶跃响应试验数据进行研究,获得的试验数据如图 3所示.
其中:Steady1为主排气阀阶跃前总压和静压稳态区间,Steady2为主排气阀阶跃后及栅指阶跃前总压和静压达到的新的稳态区间,Steady3为栅指阶跃响应后总压和静压达到的稳态区间;Dynamic1,Dynamic2分别为主排气阀和栅指阶跃响应动态区间.定义y1_Steady1,y2_Steady1分别为Steady1区间内总压和静压的稳态值,即主排气阀阶跃前总压和静压的初始值;y1_Steady2,y2_Steady2分别为Steady2区间内总压和静压的稳态值,即主排气阀阶跃后,栅指阶跃前的稳态值;y1_Steady3,y2_Steady3分别为Steady3区间内总压和静压的稳态值,而y1_Steady3,y2_Steady3是由主排气阀和栅指阶跃后得到的总压和静压的稳态值,各稳态值通过多点平均法得到.
3.2 数据预处理采集到的试验数据受噪声及测量误差影响,存在不同程度的波动,数据需预处理来提高数据质量.分别选取图 4中Dynamic1区间的总压,静压数据和Dynamic2区间的总压,静压数据作为待处理数据,应用三次样条插值法对数据进行滤波处理.图 4表示Dynamic1区间总压和静压相对主排气阀单位阶跃响应变化曲线,图 5为Dynamic2区间总压和静压相对栅指单位阶跃响应变化曲线.
由采集到的阶跃响应数据,确定各个通道的时间延迟L11,L12,L21,L22;并计算式(4)~式(7)中各通道传递函数的静态增益,计算公式为
(26) |
(27) |
(28) |
(29) |
其中:u1_step1,u2_step1分别为主排气阀和栅指阶跃前控制量输出;u1_step2,u2_step2分别为主排气阀和栅指阶跃后控制量输出.
得到的阶跃响应数据存在着耦合,需要计算出单个控制量阶跃分别对总压和静压的阶跃响应数据,才能应用面积法得到相应的传递函数.在Dynamic1区间,栅指位移保持不变,总压和静压的变化是由主排气阀阶跃引起的,因此有u2(t)=u2_step1,主排气阀的阶跃值Δu1=u1_step2-u1_step1对应的总压和静压阶跃响应数据分别由式(30)、式(31)得到:
(30) |
(31) |
其中,y1(t),y2(t)分别为总压和静压测量值,从而可用面积法得到传递函数g11(s)和g21(s)的模型参数.
同理,在Dynamic2区间,主排气阀开度保持不变,有u1(t)=u1_step2,总压和静压的变化是由栅指阶跃引起的,栅指的阶跃值为Δu2=u2_step2-u2_step1,通过式(32)、式(33)得到栅指阶跃引起的总压和静压的阶跃响应数据:
(32) |
(33) |
再通过面积法得到g12(s)和g22(s)的参数.
使用面积法辨识各个通道的参数时,对每个通道的传递函数模型给定不同的阶次,从而得到各个通道的多个传递函数模型.各个通道不同阶次的传递函数模型对应的RMSE见表 1,从表中可看出模型的阶次是影响模型预报精度的主要因素.在工业应用中,应同时考虑模型复杂程度和模型预报精度,折中选取模型阶次.
为检验本文所建模型的有效性,将得到的模型在风洞现场的控制系统平台上进行验证,将模型估计值与采集到的实时试验数据进行对比,总压和马赫数的估计值与测量值对比曲线如图 6,图 7所示.
将本文所建模型与文献[3]提出的BP-NARMAX模型进行比较,选取试验工况为总压110kPa,马赫数0.578的13组测量数据(数据采样周期为10ms)作为样本,以模型估计值与测量值的RMSE作为精度指标对两种建模方法进行比较,得到的RMSE的平均值及其范围的统计结果如表 2所示.
结果显示,本文提出的建模方法有较高的预报精度,在试验样本数据有限的情况下,RMSE模型的预报精度优于BP-NARMAX模型.
由图 6,图 7可见,估计值与测量值基本一致,表明所建模型的有效性.
5 结论本文采用多变量阶跃响应方法辨识风洞系统模型,用三次样条插值法对阶跃响应数据进行滤波处理.得到的阶跃响应数据存在着耦合,通过各阶段稳态值求取各个通道的静态增益,进而得到单个通道的阶跃响应数据,再使用面积法获得各个通道的传递函数模型.通过对不同通道得到的不同阶次传递函数模型估计的RMSE的比较,实现了模型最佳阶次的选择,从而得到风洞系统的两输入-两输出耦合多变量传递函数模型.仿真和风洞现场测试结果验证了所建模型的有效性.
[1] |
施洪昌.
高低速风洞测量与控制系统设计[M]. 北京: 国防工业出版社, 2001 .
( Shi Hong-chang. Measurement and control system design in high and low speed wind tunnel[M]. Beijing: National Defense Industry Press, 2001 . ) |
[2] | Jin Z W,Li Z,Rao Z Z.Mach number prediction models based on ensemble neural networks for wind tunnel testing[C]//26th Chinese Control and Decision Conference.Changsha,2014:1197-1200. |
[3] | Rui W,Qin J H,Ma Y Y.A novel approach for modelling of an injector powered transonic wind tunnel[C]//26th Chinese Control and Decision Conference.Changsha,2014:1637-1640. |
[4] | Wang X, Li S Y, Chai T Y. Multi-model direct adaptive decoupling control with application to the wind tunnel system[J]. ISA Transactions , 2005, 44 (11) : 131–143. |
[5] | Ljung L. System identification:theory for the user[M]. Jersey City: Prentice Hall, 1999 . |
[6] | Yu W,Zhang G J.Modelling and controller design for 2.4 m injector powered transonic wind tunnel[C]// Proceedings of the American Control Conference.Albuquerque,1997:1544-1545. |
[7] | Soeterboek A R M, Pels A F, Verbruggen H B, et al. A predictive controller for the mach number in a transonic wind tunnel[J]. IEEE Control Systems , 1991, 11 (1) : 63–72. DOI:10.1109/37.103359 |
[8] | Mei H, Li S Y. Decentralized identification for multivariable integrating processes with time delays from closed-loop step tests[J]. ISA Transactions , 2007, 46 (2) : 189–198. DOI:10.1016/j.isatra.2006.11.001 |
[9] | Chen L, Li J H, Ding R F. Identification for the second-order systems based on the step response[J]. Mathematical and Computer Modelling , 2011, 53 (5/6) : 1074–1083. |
[10] |
方崇智, 萧德云.
过程辨识[M]. 北京: 清华大学出版社, 2000 .
( Fang Chong-zhi, Xiao De-yun. Process identification[M]. Beijing: Tsinghua University Press, 2000 . ) |