一种高阶线性时不变系统的辨识算法仿真

吴相甫

（航空工业中国飞机强度研究所，西安710065）

在现代控制理论系统辨识的研究领域里，线性系统建模是控制学科的重要课题。在古典控制中，常常通过测试被控对象的单位阶跃响应来求其传递函数[1-3]。利用阶跃响应曲线来确定典型工业过程传递函数的方法很多，常用的有近似法、半对数法、切线法、两点法和面积法等。当阶跃响应曲线比较规则时，近似法、切线法就、半对数法和两点法都能比较有效地导出传递函数[4-6]。而在现代控制中，伴随着自适应控制算法的完善和发展，以系统阶跃响应或脉冲响应为依据而进行建模的算法也应运而生，如动态矩阵控制算法，它的模型是建立在对象的单位阶跃响应的基础之上。目前，工程上用得最多的是用阶跃扰动和方波扰动来测取阶跃响应。

从某种意义上讲，阶跃响应建模法提供了过程数学模型建立的一种解决方案。阶跃响应曲线形象地表示了对象的动态特性，但对于进一步地分析和研究系统却很不方便，因此需要从阶跃响应曲线求出对象的传递函数，过去常用的方法是图表法和积分法。图表法包括切线法和两点法，但都有它的局限性对于有自平衡能力的对象，要在阶跃响应曲线上找出拐点，在其拐点处作切线，由于拐点位置不易选准，且切线的方向斜率也难于确定，以致传递函数的特征参量不能准确确定。再者，根据曲线所得的数据使传递函数不便于运算和模拟，给传递函数的确定带来了一定的误差，同样，数值积分法虽是一种从阶跃响应曲线求出传递函数的通用数学方法，可以适用于各种形式的传递函数，但是由于数值计算的误差，实际上只能有效地定出少量的常数[7-8]。

本文利用系统阶跃响应的曲线特征来确定传递函数，思路是利用传递函数可以分解为一二阶典型环节的并联原理，假设模型的阶数后，调节典型环节并联模型中的参数值，用拟合曲线逼近实际的阶跃响应曲线，从而得到系统的传递函数[9]。研究时借助Matlab 平台开发辨识软件，通过曲线拟合确定典型环节并联模型的参数值，可以在一定精度内获得系统的传递函数。

1 高阶线性时不变系统模型分析

一般情况下，高阶线性时不变系统的传递函数分子分母都是s 的多项式，可以写为

式（1）可表示成如下形式：

式中：m＜n，q+2r=n

由式（2）可知，任何高阶线性时不变系统都可以分解为若干个一阶系统和二阶系统的典型环节之和。

另外，式（1）可表示为如下因式的乘积形式：

在输入为单位阶跃函数时，输出量可表示为

式中：q 为实数极点的个数；r 为共轭复数极点的对数。设0<ζk<1，将式（5）展成部分分式

设初始条件全部为零，将式（6）进行拉氏反变换，可得高阶系统的单位阶跃响应：

式（7）表明，高阶系统的时间响应是由一阶系统和二阶系统的时间响应叠加组成的。

综上所述，任何线性高阶系统都可以分解为若干个一阶和二阶典型环节的并联，且其响应也可以分解成若干个一阶和二阶典型环节的阶跃响应之和。更进一步，通过一个一阶系统和一个二阶系统的典型环节的并联模型的阶跃响应曲线逼近实际系统的阶跃响应曲线，通过调节参数

改变曲线形状，不断地逼近高阶系统的阶跃响应曲线，能够反映原系统的整体特性[10]，即可用一个一阶系统和一个二阶系统典型环节的并联模型去等效高阶线性系统的传递函数。

2 阶跃响应曲线特征分析

传递函数的参数值对系统的动态性能和稳态性能有着直接的影响，各参数的取值情况直接决定系统的阶跃响应曲线特征。一阶系统与二阶系统典型环节并联模型可表示为

各参数对阶跃响应曲线形状影响的规律如下：

增益值的影响系统阶跃响应稳态值为增益值K1、K2之和，K1、K2对上升时间、峰值时间、调节时间无影响；K1、K2之和不变时，K2越大，峰值越大，即超调量越大。

时间常数的影响T1、T2由小变大时，峰值时间和调节时间变大，这是因为一阶环节响应速度变慢；T2相对于T1会使响应的速度变慢的更快。按照T1、T2的比例关系将该结构的三阶系统分为3 类：T1≥5T2、T2＜T1＜5T2、T1＜T2，时间常数对其阶跃响应曲线形状的影响如下[11]：

（1）T1≥5T2： 阶跃响应曲线的第一个波峰值要小于稳态值。二阶系统在阶跃响应曲线前期产生作用，一阶系统在后期产生作用。

（2）T2＜T1＜5T2： 阶跃响应曲线的第一个波峰值大于稳态值。动态部分上升速度慢，超调小甚至没有超调，但曲线上有拐点。

（3）T1＜T2： 阶跃响应曲线的第一个波峰值大于稳态值。一阶环节发生作用比二阶环节早，一阶环节在上升段初期发生作用，二阶环节产生影响的时间靠后。

阻尼比的影响阻尼比ζ 增大时，超调减小，振荡次数减少。峰值时间tp由其中的二阶系统决定，根据二阶系统峰值时间的计算方法，T2与ζ 满足关系式

该式在K1、K2之和不变时，K1、K2的变化不影响该结论。

3 辨识方法

根据上一节的分析，阶跃响应稳态值由K1、K2之和决定，T1、T2决定曲线的响应速度，阻尼比决定曲线的超调。按照T1、T2的比例关系对阶跃响应曲线特征的影响将分3 类：T1≥5T2、T2＜T1＜5T2、T1＜T2，分析辨识方法。

T1≥5T2这种情况下，一阶环节的时间常数比二阶环节时间常数大得多，因此一阶环节发生作用比二阶环节晚。此时，阶跃响应曲线的第一个波峰值要小于稳态值，该特征可用来判断此系统属于T1≥5T2这一类。

具体参数辨识算法如下：

（1）一阶时间常数T1：由于一阶环节发生作用要比二阶环节要晚，曲线会先发生振荡，然后曲线会平缓上升直到达到稳态值，因此可利用一阶系统求时间常数的方法求出T1，即在阶跃响应曲线上找到输出量变化至终值95%时的坐标点，它所对应的时刻为调节时间ts，则T1=ts/3。

（2）ζ、T2、K1、K2：这4 个参数要由循环搜索逼近的方法确定。

算法包括两个循环，外层循环调节ζ、T2，内层循环调节K1、K2：首先，求出峰值时间tp，令阻尼比ζ=0.05，由式（10）可以算出T2，并令K1=K2=K/2，在该组数值下求出阶跃响应数据，算出第一个波峰差p=y1-y1′（实际阶跃曲线值减去拟合曲线值），若p>0，则令K2=K2+0.1，K1=K1-0.1（若p<0，则K2=K2-0.1，K1=K1+0.1），此为内循环（用来调节K1、K2）；当p 的绝对值小于0.1 时，在当前参数值之下算出第一个波谷差bg=y2-y2′（实际阶跃曲线值减去拟合曲线值），若bg>0.1，则说明阻尼比ζ 的值太小（振荡太大），此时令ζ=ζ+0.05，重新计算T2，此为外循环（用来调节ζ、T2），并令K1=K2=K/2，再根据波峰差p 调节K1、K2。参数辨识流程见图1。

图1 参数辨识流程（T1≥5T2）Fig.1 Flow chart of parameter identification（T1≥5T2）

T1＜T2这种情况下，一阶环节的时间常数比二阶环节时间常数小，因此一阶环节发生作用比二阶环节早。此时，阶跃响应曲线的第一个波峰值大于稳态值，该特征可用来判断此系统不属于T1≥5T2类。

在T1＜T2以及T2＜T1＜5T2时都首先假设该系统为二阶系统，求出阻尼比ζ、时间常数T2、增益K，在这3 个参数下求出这个拟合的阶跃响应拟合值y′，用实际阶跃响应曲线值y 减去y′，求出动态误差d，这样得到的动态误差d 的值会很大。经过多组实验仿真，当d>0.4 时，将待辨识系统归为T1＜T2类；当d<0.4 时，将待辨识系统归为T2＜T1＜5T2类。

对 于T1＜T2类，ζ、T1、T2、K1、K2这5个参数由3个循环搜索逼近的方法确定，具体参数辨识算法如下：

（1）ζ、T2： 这2 个参数由最外层循环搜索确定。由于含1 个一阶环节，该环节会减缓系统的振荡，因此由假设的二阶系统求出的阻尼比ζ 是偏大的，因此要向下调整，T2则由式（10）计算得出。

（2）T1、K1、K2：这3 个参数由中层和最内层循环搜索确定。第一步确定了阻尼比ζ 和T2的值后，令T1=T2，K1=K2=K/2。由于T1＜T2，T1应向下调整，每次减小0.1（此值可调整），约束条件为T1>0.1。每赋T1一个值时，K1、K2在最内层循环从K/2 开始根据峰值差p 调整，算出峰值差p=y1-y1′（实际阶跃曲线值减去拟合曲线值），若p>0，则令K2=K2+0.1，K1=K1-0.1；若p<0，则令K2=K2-0.1，K1=K1+0.1。

（3）当T1减小至接近0.1 时，若误差仍然有d>0.05，则另阻尼比ζ 减小0.05，重新计算，然后再进入第二步，如此循环搜索直至d<0.05。

参数辨识流程见图2。

图2 参数辨识流程图（T1＜T2）Fig.2 Flow chart of parameter identification（T1＜T2）

T2＜T1＜5T2这种情况与T1＜T2相同，阶跃响应曲线的第一个波峰值大于稳态值，该特征可用来判断此系统不属于T1≥5T2类。建模思想同T1＜T2相同，首先认为待辨识系统为二阶系统，并算出动态误差d，当d<0.4 时，将待辨识系统归为T2＜T1＜5T2类。

对于T2＜T1＜5T2类，ζ、T1、T2、K1、K2这5 个 参 数同样由3 个循环搜索逼近的方法确定，具体参数辨识算法如下：

（1）ζ、T2：同T1＜T2类辨识算法第一步；

（2）T1、K1、K2：这3 个参数由中层和最内层循环搜索确定。第一步确定了阻尼比ζ 和T2的值后，令T1=T2，K1=K2=K/2。由于T2＜T1＜5T2，T1应向上调整，每次增加0.1（此值可调整）。每赋T1一个值时，K1、K2在最内层循环从K/2 开始根据峰值差p 调整，算出峰值差p=y1-y1′（实际阶跃曲线值减去拟合曲线值），若p>0，则令K2=K2+0.1，K1=K1-0.1；若p<0，则令K2=K2-0.1，K1=K1+0.1。当K1、K2的调整使拟合模型的峰值满足要求时，再计算已知阶跃响应和拟合模型的第一个波谷差bg，判断当前阻尼比是否正确。

（3）当T1调整到使估计模型的峰值与稳态值大约相等时，说明T1已足够大，此时若仍有动态误差d>0.05、bg 的绝对值大于0.15 则进入第一步重新开始搜索。

参数辨识流程见图3。

图3 参数辨识流程图（T2＜T1＜5T2）Fig.3 Flow chart of parameter Identification（T2＜T1＜5T2）

4 辨识仿真

为了验证辨识算法，利用Matlab 图形用户界面GUI 来设计辨识软件，将辨识算法嵌入到软件操作界面中。基于某实验平台搭建模拟电路，如图4 所示，输入方波信号，用示波器观察响应曲线，当该电路的阶跃响应进入稳态值后，点击示波器的停止按钮，保存并导出阶跃响应数据，保存为Excel 格式文件，并导入到辨识软件中，计算出模型中的参数。

图4 模拟电路图Fig.4 Analog circuit

该电路由1 个一阶系统与1 个二阶系统并联而成，当电阻、电容参数值不相同时，这类结构的三阶系统可分为3 种情况。

T1≥5T2模拟电路各参数取值：R1=500 kΩ，C1=2 μF，R01=250 kΩ；R2=50 kΩ，C2=1 μF；R3=200 kΩ，C3=1 μF，R02=200 kΩ。

一阶环节的传递函数为

式中：时间常数T1为1，K1为2。

二阶环节的开环传递函数为

那么该二阶环节的闭环传递函数为

式中：K2为1，T2为0.1，阻尼比ζ 为0.25。

故该三阶系统的传递函数为

辨识结果如图5 所示，根据辨识算法计算出的参数值与原电路模型参数值基本接近。

T1＜T2模拟电路各参数取值：R1=100 kΩ，C1=1 μF，R01=50 kΩ；R2=500 kΩ，C2=1 μF；R3=500 kΩ，C3=1 μF，R02=500 kΩ。

一阶环节的传递函数为

图5 辨识结果（T1≥5T2）Fig.5 Identification results（T1≥5T2）

式中：时间常数T1为0.1，K1为2。

二阶环节的开环传递函数为

那么该二阶环节的闭环传递函数为

式中：K2为1，T2为0.5，阻尼比ζ 为0.5。

故该三阶系统的传递函数为

辨识结果如图6 所示，根据辨识算法计算出的参数值与原电路模型参数值基本接近。

图6 辨识结果（T1＜T2）Fig.6 Dentification results（T1＜T2）

T2＜T1＜5T2模拟电路各参数取值：R1=200 kΩ，C1=5 μF，R01=200 kΩ；R2=50 kΩ，C2=10 μF；R3=500 kΩ，C3=1 μF，R02=500 kΩ。

一阶环节的传递函数为

式中：时间常数T1为1，K1为1。

二阶环节的开环传递函数为

那么该二阶环节的闭环传递函数为

式中：K2为2，T2为0.5，阻尼比ζ 为0.5。

故该三阶系统的传递函数为

辨识结果如图7 所示，根据辨识算法计算出的参数值与原电路模型参数值基本接近。

图7 辨识结果（T2＜T1＜5T2）Fig.7 Dentification results（T2＜T1＜5T2）

5 结语

由辨识仿真结果看出，根据辨识算法计算出的参数值与原电路模型参数值基本接近，并且拟合曲线和原曲线十分接近，动、静态误差也满足辨识要求。辨识结果存在误差主要有两个原因：

（1）实验箱电阻、电容的标识值与实际值可能有偏差：阶跃响应数据是基于实验箱模拟电路得到的，若电阻、电容的标识值与实际值有偏差，那么从实验箱获取的阶跃响应数据也就和事先设计的系统参数不同，根据这样的阶跃响应数据辨识出来的参数也就必然与事先设计的系统不同。

（2）算法需深入研究：搜索逼近参数时设定的步长在算法中是固定的，会影响辨识精度和软件搜索时间；另外，用于判断终止搜索条件的动态误差计算方法也会影响辨识的结果。