当前位置:文档之家› 系统辨识MATLAB仿真

系统辨识MATLAB仿真

设一非线性系统如下所示:()()()()()111y k y k e u k k βαε--=-+-+0.75α= 0.35β=0.25γ=,()k ε是零均值,方差为0.01的噪声序列(均匀白噪声)。

(1)试设计一种激励信号能持续激励系统的各工作点(平衡点)(2)用适当的方法辨识出系统的等价模型(用另一组数据来检验模型的泛化性) 说明:下面讨论的都是离散系统,所以时间坐标均采用离散时间节点k 。

解:(1) 线性化处理寻找系统的合理输入信号 可以求得系统的平衡点为:()0.75y k α== (1.1)按题意要求最后系统必须收敛于平衡点附近,即:()lim 0.75k y k →∞= (1.2)为了找出系统的合理输入信号,使得系统最终工作在平衡点附近,这里可以将系统线性化处理,将上述非线性系统进行泰勒展开得:()()()()()()()2323111111112!3!!n n y k y k y k y k y k u k k n αβαβαβαβγε--+---++-=-+因为 ,后面()1y k -的高阶项都可以扔掉(只作为寻找输入信号使用), 所以系统可以化为下式:此时不妨设系统输出()y k 的最后的极限为A , 从式1.2得0.75A =。

那么应该满足()()lim lim 1k k y k y k A →∞→∞=-= (1.5)从而有 ()()()11A u k k αβγε-=-+ (1.6)同时为了抵消系统的部分噪声,这里采用MA TLAB 软件编程产生另一服从同一分布的均匀噪声()1k ε,将式1.6变形得:()()()01111u k u k k αβεγγ--=-(1.7)式1.7中()0u k 是一个最后收敛于系统平衡点0.75的基本信号,这里可以采用一阶线性系统的阶跃响应曲线作为基本信号()0u k ,同时考虑系统的平衡点,即设计为:()/00.75k Tu k e-=- (1.8)T 是一阶线性系统对应的时间常数,反应到输入基本信号()0u k 上就是过零点作()0u k 对应()()()()11y k y k u k k αβγε--=-+01αβ<(1.3)(1.4)曲线的切线的斜率1KT =。

图1 基本信号产生框图很显然有()lim0.75tu k→∞=(1.9)从而最后的输入信号设计为:()()()011111u k u k kαβεγγ-=+-+(1.10)式1.10中,,αβγ题中已知,()1kε满足:()1E kε=⎡⎤⎣⎦,()10.01D kε=⎡⎤⎣⎦,()u k 的产生参照图1。

代入式1.8得:()()(1)/111(0.75)1k Tu k e kαβεγγ-+-=--+(1.11)产生()u k过程中的时间常数T由仿真效果决定,如果从系统的性能指标而言,T尽可能小,仿真如下:图2 噪声随机抵消情况下的原系统响应曲线系统输入基本信号()u k阶跃信号11Ts+图3噪声完全抵消情况下的原系统响应曲线仿真分析:图2与图3中都是原非线性系统的响应曲线图,两者的输入信号有差别,分别为()()12,u k u k ,如下式表示:()()()1011111u k u k k αβεγγ-=+-+ (1.12)()()()201111u k u k k αβεγγ-=+-+ (1.13)()0u k 参照式1.8。

式1.12中()11k ε+是MA TLAB 中产生的符合同一均匀分布的给定随机噪声,从而式1.12中输入信号()1u k 实现噪声随机抵消,结果系统的响应()y k 呈现出很大的噪声干扰,这是符合对未知系统认识的逻辑过程的,如果要得出理想的响应曲线,也就是要求系统最终在平衡点附近趋于稳定,那么需要设计相应的滤波器,这里不再讨论,后续会尝试设计滤波器滤波。

式1.13中()1k ε+是系统给定的噪声,那么设计的输入信号()u k 实现了对噪声的完全抵消,得到的理想响应曲线如图3所示,从图3中可以看出,()y k 很快的收敛于0.75附近,且上升时间以及超调量等性能指标均符合理想响应要求,这里就没有具体计算,显然符合控制系统的快稳准的三要素。

(2) 最小二乘法辨识模型通过MATLAB 编程从原非线性系统中采集一些输入输出数据,利用最小二乘法辨识出系统近似的线性模型,实现框图如图4。

时不变SISO 系统动态过程的数学模型为()()()()()11A z y k B z u k k ε--=+ (2.1)其中()u k 是系统的输入,()y k 是系统的输出,()()11,B A z z --的展开式分别如下:()112121a a n n A z a z a z a z ----=+++ (2.2) ()11212b b n n B z b z b z b z ----=++ (2.3)参数1212,,,,,,,a b n n a a a b b b 待辨识。

式2.1写成最小二乘格式()()()Ty k hk k θε=+ (2.4)式2.4中:()()()()()[y 1,,y ,u 1,,u ]T a b h k k k n k k n =------ (2.5)1212,,,,,,,a b Tn n a a a b b b θ⎡⎤=⎣⎦ (2.6)对于1,2,,,k L =方程2.4构成一个线性方程组,可以把它写成式2.7,L L L Y H E θ=+ (2.7)其中()()()1,y 2,,y TL Y y L =⎡⎤⎣⎦ (2.8) ()()()1,2,,TL E L εεε=⎡⎤⎣⎦ (2.9)()()()()()()()()()()()()()()()101011212211T a b T a b L T a b h y y n u u n y y n u u n h H y L y L n u L u L n h L ⎡⎤----⎡⎤⎢⎥⎢⎥----⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥------⎣⎦⎣⎦ (2.10)由题意知()k ε为均匀分布的噪声所以{}(){}(){}(){}120L E E E E E L εεε⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦(2.11){}{}(){}()(){}()(){}()(){}(){}()(){}()(){}()(){}(){}22221121212212T L L L E E E L E E E L Cov E E E E I E L E L E L εεεεεεεεεεεσεεεεε⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦(2.12)式中20.01εσ=,I 是单位矩阵。

这里认为噪声()k ε与输入()u k 是不相关的,从而()(){}0,,E k u k l k l ε-=∀ (2.13)这里需要考虑的是记忆长度的选择,上述矩阵方程包含()a b n n +个未知数,一般要求a b L n n >+,这样才可能确定一个“最优”的模型参数θ,而且为了保证辨识的精度,L 必须充分大,由AIC 定阶准则知L 同样决定着所辨识系统的阶次,下面会阐述相应理论。

定义准则函数:()()()()21LTk J k y k h k θθ=⎡⎤=Λ-⎣⎦∑ (2.14)其中()k Λ是加权因子,它衡量的是采样数据的可信度,这里对数据的可信度难以肯定,就简单的将取()1k Λ=。

将准则函数写成二次型的形式如式2.15。

()()()TL L L L J Y H Y H θθθ=-- (2.15)设wls θ∧使得()|min wls J θθ=()()()||0T T L L L L wls wls J Y H Y H θθθθθθθ∧∧∂∂=--=∂∂ (2.16) 利用下面的两个向量微分公式展开得()()2T TT T x xx Ax x A x αα∂⎧=⎪⎪∂⎨∂⎪=⎪∂⎩(2.17) 其中A 是对称阵,得到正则方程()TTwls LL L L HH H Y θ∧= (2.18)当TL L H H 为正定矩阵时,有:()1T Twls L L L L H H H Y θ∧-= (2.19)从而()22|20wlsTL L J H H θθθ∧∂=>∂ (2.20) 结合式2.16到2.20知,wls θ∧使()|min wls J θθ=,并且wls θ∧是唯一的。

从而求得最小二乘估计值 ()1T T LS LLL L H H H Y θ∧-= (2.21)定阶:系统的等价线性模型的阶次确定可以基于AIC 准则来确定,首先将系统写成如式2.7的形式。

输出变量在n θ条件下的似然函数为:()()()22212exp 2L T n n n n n n n L Y H Y H θπσθθσ-∧∧⎧⎫⎛⎫⎪⎪=---⎨⎬ ⎪⎝⎭⎪⎪⎩⎭(2.22) 对应的对数似然函数为:()()()()221ln ln 2ln 222Tn n n n n n n n L L L l Y H Y H θθπσθθσ∧∧==----- (2.23) 其中L 是数据长度,根据极大似然原理得到模型参数θ及白噪声()k ε方差的极大似然估计值如下:()121T TML n n n n T ML ML n n n n H H H Y Y H Y H L θσθθ∧-∧∧∧⎧=⎪⎪⎨⎛⎫⎛⎫⎪=-- ⎪ ⎪⎪⎝⎭⎝⎭⎩(2.24) 将ML θ∧与2σ∧回代式2.23可得:2ln 2ML L l const θσ∧∧⎛⎫=- ⎪⎝⎭(2.25)AIC 准则的标准式为:2ln 2ML AIC N L N θ∧∧∧⎛⎫⎛⎫=-+ ⎪ ⎪⎝⎭⎝⎭(2.26)将式2.25代入AIC 准则的标准式2.26中得到:2ln 4AIC n L n σ∧∧∧⎛⎫=+ ⎪⎝⎭(2.27)上面各式中n ∧表示的就是模型要估计的阶次,最后找到使min AIC n ∧⎛⎫= ⎪⎝⎭的n ∧作为模型的阶次。

本次辨识采样的数据长度1000L =, 很显然2σ∧足够小,当数据长度取定时,显然n ∧可以尽可能的小,暂且取n ∧=4。

仿真分析:图5 定4阶情况下系统的辨识效果图6 定4阶情况下更换输入为原来的14情况下系统的辨识效果图7 定4阶情况下更换输入为原来的4倍情况下系统的辨识效果图5与图6及图7中分别给出了定4阶情况下的系统利用最小二乘法辨识的结果,图5中的辨识效果在幅值上出现减小,有一定的误差,但基本可以与原系统同步,为了检验模型的泛化性,在MATLAB编程过程中,将输入信号减小为原来的14以及增大为原来的4倍,图6中得到了较好的辨识效果,但图7却进一步增大幅度误差,从这里看出所辨识模型的泛化性并不好,即随着输入信号的增大,模型与系统的幅度差值增大。

相关主题