1、第5章 控制系统的仿真建模 5.1 问题的描述与模型的定义问题的描述与模型的定义 5.2 控制系统模型的建立控制系统模型的建立5.3 控制系统模型的确认和修改控制系统模型的确认和修改 5.4 控制系统仿真建模的实例控制系统仿真建模的实例 小结小结 5.1 问题的描述与模型的定义控制系统的分析、设计、仿真等都是以数学模型数学模型为基础的。上一页下一页返回5.2 控制系统模型的建立5.2.1 建模要求和原则建模要求和原则5.2.2 建模过程的信息源建模过程的信息源5.2.3 建模方法建模方法5.2.4 最小二乘参数估计最小二乘参数估计5.2.5 模型的阶次辨识模型的阶次辨识5.2.6 控制系统仿真
2、建模的步骤控制系统仿真建模的步骤上一页下一页返回5.2.1 建模要求和原则模型与真实世界之间最重要的关系就是抽象抽象和映射映射。在建立一个数学描述时,首先需要建立几个抽象,即定义以下几个集合:输入集、输出集、状态变量集输入集、输出集、状态变量集。定义了上述集合之后,再在这些抽象的基础上,建立复合的集合结构,包括一些特定的函数关系,通常称这个过程为理论构造。实现抽象模型结构与实际系统之间的联系的过程,称为映射。合理的模型应该是这样的:模型的复杂程度能适度描述一个给定的系统,即在实用前提下的最优(参数最少)。上一页下一页返回模型过于复杂如何?模型过于复杂如何?模型过于简单呢?模型过于简单呢?5.2
3、.2 建模过程的信息源 建模过程涉及许多信息源,其中主要的是三类:建模目的、先验知识和试验数据,它们的关系如图5.1所示。图5.1 1建模的目的 数学模型是对实际系统的一种相似描述。目的的不同、选择实体的不同将导致建模过程沿不同方向进行。2先验知识 建模过程应当尽量利用以往的知识源。3试验数据 合适的定量观测是解决建模的另一途径。上一页下一页返回5.2.3 建模方法 1机理建模法(机理建模法(质量质量-弹簧弹簧-阻尼器系统的建模阻尼器系统的建模)运用先验信息,通过数学上的逻辑推导和演绎推理,从理论上建立描述系统中各部分的数学表达式或逻辑表达式。是从一般到特殊的方法,试验数据只被用来证实或否定原
4、始的假设或原理。也称为演绎法或理论建模法。2试验建模法试验建模法 根据观测到的系统行为结果导出与观测结果相符合的模型,是一个从特殊到一般的方法。也称为归纳法或系统辨识法。系统辨识就是按照一个准则在一组模型中选取一个与数据拟合得最好的模型。上一页下一页返回 3综合(混合)建模法综合(混合)建模法 对于那些内部结构和特性有些了解但又不十分清楚的。“灰色”系统(灰箱问题),只能采用综合建模法(机理法、辨识法以及其他一些方法)。上一页下一页返回5.2.4 最小二乘参数估计1.最小二乘法原理最小二乘法原理 2.最小二乘估计的批处理算法最小二乘估计的批处理算法 3.最小二乘估计量的统计性质最小二乘估计量的
5、统计性质 4.参数个数的递推参数个数的递推 上一页下一页返回1.最小二乘法原理 假设z是一根金属轴的长度,t是该金属轴的温度,希望确定轴长z和温度t之间的关系。具体方法是:首先在不同温度t下对变量z进行观测,得到试验数据;然后根据试验数据,寻找一个函数去 拟合它们,同时要确定该函数关系式中的未知参数的值。假设已经确定了模型的类型和结构,即轴长z和温度t之间有如下的线性关系(先验知识)z=z0(1+t)式中,z0是0时金属轴的长度,是膨胀系数。上一页下一页返回 如果令z0=a,z0=b,则上述函数关系可写成 z=a+bt (5.1)每次观测中总带有测量误差。因此,每次观测所得的轴长并不是真正的轴
6、长zi,而是yi。yi可写成 yi=zi(真值)+vi(随机观测误差)或 yi=a+bti+vi (5.2)式中,yi是可观测的随机变量,ti是可观测的独立变量(非随机变量),vi是不可预测的随机性观测噪声,a和b是待估计的未知参数。根据N次观测数据ti,yi,i=1,N,来估计出模型(5.1)中的未知参数a和b的值。上一页下一页返回两个含义:两个含义:想求得a,b的真实值是不可能的,只有最优估计值 如何让估计出的值与真实值最接近?通常采用各次误差的平方和作为总误差 (5.4)这个误差平方和函数误差平方和函数就是在估计参数时所采用的准则函数(或称为性能指标)。由于平方运算也称为“二乘”运算,因
7、此,按照这种原则来估计参数a和b的值的方法称为最小二乘估计法(LS法)。要使使J达到极小值达到极小值,只需分别对a和b求偏导数,并令它们等于零。a和b的估计值应满足下列条件NiiiNjibtayvJ1212)(Niiibbaat bayaJ10)(2NiiiibbaattbaybJ10)(2上一页下一页返回而 和 由下列方程组确定 (5.5)得 (5.6)通常将 和 称为最小二乘估计量(LSE)。NiNiNiiiiiNiNiiitytatbyNatb1112112112111211211112)()(NiiNiiNiNiiNiiiiNiiNiiiNiiNiiNiiNiittNytytNbttN
8、yttytaa a bb上一页下一页返回2.最小二乘估计的批处理算法 把上述最小二乘法原理应用到 离散模型(5.7)(5.7)这里假设待辨识系统的阶次为n。由于观测数据有误差,实际上观测到的系统输出不是z(i),而是y(i),如图5.2所示。niiniiikubikzakz11)()()(Nk,2,1图5.2上一页下一页返回将y(i)代入(5.7)式,得 (5.8)式中,e(k)称为残差,表示用观测值y(i)取代模型(5.7)式的计算结果产生的误差。取 (5.9)作为模型参数估计的准则函数,则使J为极小的参数估计就是最小二乘估计。令 (5.10)niniiikeikubikyaky11)()(
9、)()(NiieJ12)(nnbbaa11)()2()1()(NyyyN)()2()1()(NeeeN上一页下一页返回式中,是待估计的2n维参数向量;(N)是N维观测向量;(N)是N维误差向量;(N)是N2n维数据矩阵。(5.8)式和(5.9)式可以重新表示为 (5.12)(5.13))()1()()1()2()1()2()1()1()0()1()0()(nNuNunNyNynuunyynuunyyN(5.11))()()(NNN)()(NNJT上一页下一页返回将(5.12)式代入(5.13)式,可得 (5.14)将J对求导,得 (5.15)在(5.15)式中令 ,由可求得的最小二乘估计 (5
10、.16)()()()(NNNNJT)()()()()()()()(NNNNNNNNTTTTTT)(2)()(2NNNJTT0J)(N)()()()(NNNNTT上一页下一页返回 这个方程叫做正则方程。由此方程就可以得到最小二乘估计值为 (5.17)()()()(1NNNNTT上一页下一页返回 (5.17)式所示的参数向量 是根据N组实际观测数据(N)和y对参数的估计值。选用的准则是最小二乘准则,所以称它为最小二乘估计量。为了与其它估计量区别起见,用符号 表示。显然,估计量 是被观测数据y的线性函数,所以最小二乘估计是一种线性估计。(5.9)式表示的准则函数中,假定了每次观测都具有相同的重要性。
11、若考虑到各次观测的不同情况,则可以在准则函数中引入对每次观测的重视程度因子(称为权重)w(i),得LSLS)()()2()2()1()1(222NeNwewewJwNiieiw12)()(上一页下一页返回用矩阵表示,就是 (5.18)式中,权矩阵W为N维对角线矩阵,并规定它是个对称正定阵。加权最小二乘估计量 应满足的正则方程为 (5.19)当 非奇异时,可求得加权最小二乘估计量为 (5.20)当W=I(单位矩阵),加权最小二乘估计就成为最小二乘估计。所以,最小二乘估计是加权最小二乘估计的特例。)()(NNJTWW)()()()(LSWNNNNTTWW)()(NNTW)()()()(1LSWNN
12、NNTTWWLSW上一页下一页返回3.最小二乘估计量的统计性质 可以证明:若观测噪声等引起的偏差(N)是白噪声,即:偏差序列e(i)是独立的序列,其数学期望为0,方差为2,用数学公式表示就是 LSE是无偏的、有效的和一致的。一般观测总是在尽可能相同的条件下进行的,每次观测不受其它观测的影响。因此,观测的噪声干扰是近似具有白噪声特性的。以上结论表明,LS法是一种比较好的参数估计法。0)(NE I2)()(NNET上一页下一页返回4.参数个数的递推 前面介绍的算法,是在假设系统阶次n已知的情况下,利用已得到观测数据来计算出参数估计量。但是,在实际中有时系统的阶次也不知道。这时,不仅要估计系统的参数
13、,而且要估计系统的阶次。即要用不同阶次的数学模型反复地进行参数估计,直到得到较理想的结果为止。要进行这样的工作就要用同一批数据,在选择不同的阶次n的情况下,反复去解方程(5.17)式或(5.19)式。假设系统阶次为n时,已求出2n个参数估计值后,在假设系统阶次为n+m时,看看这(n+m)个参数估计值如何从原来的2n个参数中递推出来。上一页下一页返回 先介绍一个分块矩阵求逆公式。一般的n+m阶非奇异分块矩阵 中,A A11和A A22分别为n阶和m阶可逆方阵,A A12和A A21分别为nm阶和mn阶矩阵,若D D=A A22-A A21A A11-1A A12为非奇异矩阵,则有 (5.22)设
14、新的参数估计值向量为 (5.23)式中,的维数为2n;的维数为2m,即有 22211211AAAAA1111211112111111211121111)(DAADDAAAADAIAA21*12上一页下一页返回*1*11nnbbaa*12mnmnmnnbbaa(5.24)同样,新的数据矩阵 也是在原有的基础上增维而成的,即有 (5.25)式中,由(5.11)式给定,而 (5.26)2*,)(N)(N)()1()()1()2()12()2()12()1()11()1()11()(22mnNunNumnNynNymnunumnynymnunumnynyN上一页下一页返回由正则方程 ,有 (5.27)
15、利用(5.22)式可得 (5.28)TT*TTTT22122,TTTTTT2212222TTTTTT21222221233123112222PPPPPPTTTTTT上一页下一页返回式中 (5.29)(5.30)(5.31)PT11212222PPTTT2231PPPT上一页下一页返回把(5.28)式代入(5.27)式,得 即 (5.32)于是,新的阶次下的参数估计量 就可以由原来阶次下的估计量 递推计算出来。PPPPPPTTTT2233123121*PPPPPPTTTTTTT223231231PP2223*TT*上一页下一页返回5.2.5 模型的阶次辨识 1模型的阶次与阶次辨识 建立一个系统的
16、数学模型时,可能会遇到下列三种情况 系统的运动规律是完全了解的,而且系统本身是线性的,可用常微分方程来描述的。这样的系统可以用一个相当精确的线性差分方程来描述。模型的阶次是已知的,建模工作只是估计模型中的参数,无需辨识其阶次。实际系统的阶数是客观存在的,但对系统的先验知识的了解不足以知道其确切值。于是在建模时既要进行阶次辨识,又要进行参数估计。上一页下一页返回 实际系统是分布参数类型的,需要用偏微分方程来描述,而且系统往往是非线性的,甚至其运动规律都还不很清楚。对于这类系统,即使建立了机理模型,要想辨识其参数也是极其困难的,并且实际应用起来也会十分不方便。为此,总希望能用一个线性的差分方程来近
17、似地描述之。这时,数学模型不再具有明显的物理意义,而只是一个抽象的模型,是人为地拼凑出来的。模型的阶次与原来系统就没有什么直接关系了,只是“拼凑”出来的。很多阶次辨识方法都是用参数估计方法作为工具。一般是假定系统阶次是n*=1,2,p,然后在各个假定的阶次下进行参数估计,同时确定一个准则或是一个标准,把不同阶次下得到的结果进行仿真比较,从中找出最合适的阶次以及相应的参数估计量。上一页下一页返回 2损失函数检测法(1)损失函数曲线 残差平方和是参数估计量和阶次的函数,又称为损失函数,是对模型精度的一种度量。一般而言,损失函数是随着阶次的增加而递减的,并不一定存在最小值。一般损失函数J与阶次的关系
18、如图5.3所示。起初,随着阶次的增高损失函数下降得很快,然后就基本上就保持不变,有时还会略微有点回升。图5.3上一页下一页返回 因此,阶次的估计值一般选择损失函数平坦部分左端点的阶次值,如图5.3中可选3或4。这时损失函数值较小,阶次也不太高。为此,需要检查J(n*+1)和J(n*),如果它们没有显著的差异,则n*就是阶次的估计值 。(2)显著性检验 为了确定什么是损失函数的“显著”变化,通常采用统计学中的F检验法。在进行最小二乘估计时,采用的观测方法应该是相同和独立的。得到的观测数据中包含的系统的噪声是正态分布的。这时,最小二乘法的残差也成正态分布。可以证明,随机变量J(n2*)与J(n1*
19、)-J(n2*)是独立的随机变量并且有x2分布的特性。如果阶次由n1*增加到n2*,需要辨识的参数个数便由2n1*增加到2n2*。n 上一页下一页返回 基于N 组观测数据的统计量是为 (5.33)如果N 很大(大样本),则 f 是 F(f1,f2)渐近分布。它有二个自由度 f1=2(n2*-n1*)和 f2=N-2n2*。由于在进行阶次辨识时,每次阶次增加1,n2*-n1*=1,因此 f1=2。)(22)()()(*1*2*2*2*2*1nnnNnJnJnJf上一页下一页返回f21020304060801001200.12.922.592.492.442.392.372.362.352.300
20、.054.103.493.323.233.153.113.093.073.000.017.565.855.395.14.984.884.854.794.61 表5.1给出了在 f1=2 情况下F分布的临界值表。其中为风险水平(1-为置信度)。表5.1 F分布表(f1=2)上一页下一页返回(3)置信度 在进行F检验时,选择置信度(即风险水平)。置信度的选择完全是人为的。当置信度增加时,风险水平就降低了。此时,临界值增大,按(5.33)式计算出来的f值,如果小于该临界值,就能以更大的把握断定J(n2*)与J(n1*)有显著的变化。一般统计检验的置信度在0.90与0.99之间(即=0.100.01)
21、,工程中常选置信度为0.95。从表5.1可以看出,当=0.05时,在大样本的情况下(f260),F分布的临界值的变化就很慢了。因此,在数据较多的情况下,可以使用 f2 的临界值。工程中进行系统辨识时,如果 f2150,在风险水平=0.05 的情况下,常常取F检验的临界值 f*=3.00。上一页下一页返回5.2.6 控制系统仿真建模的步骤假定N-2p150(p为待辨识模型可能的最大阶次),并且取风险水平=0.05。1采取数据u(i),y(i),置2按(5.10)式构造=(N),按(5.11)式构造=(N)3用(5.17)式求得参数估计量 ,并计算残差平方和。即 niniiiikubikyakyk
22、e11)()()()(Nk,2,1NkkeJ11)(00.3*f1n上一页下一页返回 4置 ,按(5.26)式构造2(m=1)5用(5.32)式,计算得 ,并计算残差平方和:6计算 7若ff*,则输出n及 ,停机;否则,转8 8判np,若成立,则输出辨识失败信息,停机;否则,置 并在新的水平下,重新构造=(N),转4 1*nn*NkniniiikeJikubikyakyke1211*)()()()()(22*221nNJJJf*nn*21JJ上一页下一页返回5.3 控制系统模型的确认和修改 确认模型是否合适的最简单方法是,将施加到实际系统上的输入同时施加到模型上,然后比较实际系统输出和模型输出
23、之间的一致性。还可以采用下面几种方法。1与验前信息的一致性 判断模型是否符合对所建模型提出的要求,与系统已有知识或直觉的理解相比较。2交叉验证 交叉验证就是将由系统的一组观测数据产生的模型与同一系统在相同条件下得到的另一组独立观测数据产生的模型进行比较。上一页下一页返回 3原始参数的核验 在系统的参数模型中,包含有各个原始参数i。这些参数往往具有一定的物理意义。如果出现原始参数的含义并未得到保持,即所求出的i的数值与其应有的物理含义并不一致的情况,所建模型是不能被接受的。4模型的满意度 对模型的满意度可以根据建模的目的来衡量。如果模型是用来设计调节器的,当根据模型所设计的调节器能够满意地运行,
24、该模型就是一个有用的模型;如果模型是用于预测的,当模型的预测结果与系统的实际结果相差不大时,该模型就是一个较好的预测模型。上一页下一页返回5.4 控制系统仿真建模的实例5.4.1 机理建模法机理建模法5.4.2 试验方法试验方法5.4.3 基于基于MATLAB的系统模型的估计的系统模型的估计上一页下一页返回5.4.1 机理建模法【例5.1】控制系统如图5.4所示,该系统是由自整角机和伺服结构组成的随动系统。理想位置和实际位置的偏差经过自整角机转变为电压信号,该信号与速度反馈信号之差经放大去控制电机,再经减速器带动负载。通过控制使系统的实际位置跟随理想位置。试运用机理模型法建立系统的数学模型。图
25、5.4上一页下一页返回 【解】由图5.4,系统为一位置伺服闭环控制系统,将其分解为基本元件或部件,按工作机理分别列写输入-输出动态方程,并按各元件、部件之间的关系,画出系统结构图,最后根据结构图求出系统的总传递函数,从而建立起系统的数学模型。(1)同步误差检测器 设其输入为给定角位移r与实际角位移c之差,输出为位移误差电压u1,且位移-电压转换系数为k1,有u1=k1(rc)上一页下一页返回 (2)放大器(A)设其输入为位移误差电压u1与测速发电机反馈电压u2之差,输出为直流电动机端电压u,电压放大系数k2,有 u=k2(u1-u2)上一页下一页返回 (3)直流电动机(M)设其输入为u,输出为
26、电动机电角度,R为电枢回路等效电阻,L为电枢回路等效电感,km为电磁转矩系数,J为电动机转动惯量,忽略反电动势和负载转矩影响,则由电动机电压平衡方程和力矩平衡方程,有 式中,ia为电枢电流。消去了中间变量电枢电流ia,有 式中,T为电动机电磁时间常数,T=L/R;k3为电压/速度转换系数,k3=km/RJ。aaRitiLuddtJikamddukttT322dddd上一页下一页返回(4)测速发电机 设其输入为电动机角速度,输出为测速电压值u2,速度-电压转换系数为k4,有 u2=k4(5)负载输出 设输入为电动机角速度,输出为负载角位移c,传动比 n=N1/N21,有(6)构成系统各部分结构图
27、 将上述各式进行拉氏变换,按输入输出关系表示出各环节传递函数,并据此画出各部分的结构如图5.5所示。ntcdd上一页下一页返回图5.5上一页下一页返回(7)求出系统模型 按相互之间的作用关系,构成系统总结构图(如图5.6所示)。化简后,求出该系统总传递函数GB(s)即为所需的系统数学模型。图5.6nkkkskkksTsnkkksssGrcB32143223321)()()(上一页下一页返回5.4.2 试验方法【例5.2】用试验方法测得某系统的开环频率响应数据如表5.2所示,试由表中数据建立该系统开环传递函数模型G(s)。0.100.140.230.370.600.951.532.443.916
28、.2510.0-0.049-0.102-0.258-0.638-1.507-3.270-6.315-10.81-16.69-23.65-31.27-9.72-14.12-22.45-35.35-54.56-81.25-115.5-157.2-207.8-271.7-358.9表5.2 系统的开环频率响应实测数据表)srad/(1dB/)(L)/()(o注:为输入信号角频率;L()为输出信号对数幅频特性值;()为输出信号对数相频特性值。上一页下一页返回【解】该问题可分为以下几步进行研究。由已知数据绘制该系统开环频率响应,如图5.7所示。图5.7上一页下一页返回 用20dB/dec及其倍数的折线逼
29、近幅频特性,得两个转折频率,即1=1rad/s,2=2.85rad/s 求出相应惯性环节的时间常数为T1=1/1=1s,T2=1/2=0.35s 由低频段幅频特性可知,L()|0=0,所以 K=1。由高频段相频特性知,相位滞后已超过-180,且随增大,滞后愈加严重,显然该系统存在纯滞后环节 e-s,为非最小相位系统。因此,系统开环传递函数应为以下形式sses.ssTsTKesG)1350)(1(1)1)(1()(21上一页下一页返回 设法确定纯滞后时间值。查图中=1=1rad/s时,(1)=86,而按所求得的传递函数,应有 (1)=arctan1 arctan0.35 1180/=86 易解得
30、1=0.37s当=2=2.85rad/s时,(2)=169,同样从(2)=arctan2.85 arctan(0.352.85)2.851180/=169 解得 2=0.33s取两次平均值得 =(1+2)/2=0.35s 最终求得该系统开环传递函数模型为s.-ses(ssTsTKesG35021)135.0)(11)1)(1()(上一页下一页返回5.4.3 基于MATLAB的系统模型的估计 1AR模型参数估计 AR(自回归)模型具有以下形式 (5.34)式中na为A(q-1)的阶次。)()()(1kekyqAnanaqaqaqaaqA221101)(上一页下一页返回 MATLAB中提供了ar函
31、数,用于进行AR模型的估计,基本调用格式为 thar(y,n)thar(y,n,approach)其中,y为采集的时间序列辨识信号(列向量);n为AR模型的阶次;approach为所给出的方法,选项有以下几种:fb 前向反馈方法(缺省方式);1s 最小二乘法;yw Yule-Walker方法;burg Burg方法。上一页下一页返回 2ARX模型参数估计 ARX(自回归滑动平均)模型具有以下形式 (5.35)式中 )()()()()(11kenkuqBkyqAbnanaqaqaqaaqA221101)(nbnbqbqbqbqB22111)(上一页下一页返回 MATLAB中提供了基于最小二乘法的
32、ARX模型估计函数arx与基于辅助变量法的iv4函数,基本调用格式为th=arx(Z,NN)thiv4(Z,NN)其中,Zy,u为输出输入数据集,y与u为列向量。对于多变量系统,Z=yl y2yp u1 u2um,也可以仅仅是时间序列;NN=na nb nk为该模型的阶次和延迟。th变量返回ARX的估计模型,包括估计协方差与结构信息,它的准确格式可以由theta(为系统初始估计模型的返回变量)加以显示。上一页下一页返回【例5.3】已知系统输入输出数据集dryer2,共1000组数据,图5.8给出了部分输入(功率)、输出(温度)数据图,试估计ARX模型的参数并对估计模型进行检验。图5.8上一页下
33、一页返回【解】采用MATLAB编程实现。首先调入相关数据并对其进行预处理,再应用arx函数、iv4函数对其进行参数估计,得到参数估计值。文件名为exam5_3.m。程序如下:%这是例5.3的程序load dryer2%从文件加载数据到workspacedry=iddata(y2,u2,0.08);%构造iddata object,y2为输出,u2为输入,%0.08 为采样间隔%给输入和输出命名dry.InputName=Power;dry.OutputName=Temperature;%选择前300个数据建模ze=dry(1:300);%将200到300的数据绘制成图pause,plot(ze
34、(200:300),pause%按任意键得到图 上一页下一页返回 ze=dtrend(ze);%应用dtrend函数对信号进行数据预处理,将其零均值化 m2=arx(ze,2 2 3);%估计ARX模型的参数%查看模型参数 pause,m2%按任意键得到模型参数 pause%按任意键继续运行%利用有效数据验证m2模型的特性 compare(ze,m2,*);pause 运行exam5_3.m。下面是程序运行过程中MATLAB Command Window显示的内容:上一页下一页返回Discrete-time IDPOLY model:A(q)y(t)=B(q)u(t)+e(t)A(q)=1-1
35、.274 q-1+0.3935 q-2 B(q)=0.06662 q-3+0.04448 q-4 Estimated using ARX from data set ze Loss function 0.00166284 and FPE 0.00170778 Sampling interval:0.08图5.9 并且得到如图5.9所示的结果图形。上一页下一页返回 由运行过程中显示的内容可以看出,采用最小二乘法得到ARX估计模型为即有 从图5.9可以看出,测量输出和模型输出的吻合率达到了88.93%。因此,所求得的模型是可用的。)()()04448.006662.0()()3935.0274.11(4321kekuzzkyzz)()4(04448.0)3(06662.0)2(3935.0)1(274.1)(kekukukykyky上一页下一页返回小结本章主要介绍控制系统仿真建模的定义、建模的要求和方法。对最小二乘参数估计基本原理、最小二乘估计的批处理算法和模型阶次辨识进行了详细介绍,并给出了控制系统仿真建模的步骤和模型的确认与修改方法。最后通过具体例子来说明机理建模法、试验方法和最小二乘参数估计法的实际应用。上一页返回