系统辨识MATLAB仿真汇编

1、设一非线性系统如下所示y k 1 e yk 1 u k 1 k0.75 0.35 0.25, k 是零均值，方差为 0.01 的噪声序列(均匀白噪声)1)试设计一种激励信号能持续激励系统的各工作点(平衡点)2)用适当的方法辨识出系统的等价模型(用另一组数据来检验模型的泛化性)说明：下面讨论的都是离散系统，所以时间坐标均采用离散时间节点 k 。解：(1) 线性化处理寻找系统的合理输入信号 可以求得系统的平衡点为：y k 0.75 (1.1) 按题意要求最后系统必须收敛于平衡点附近，即：lim y k 0.75 (1.2) 为了找出系统的合理输入信号， 使得系统最终工作在平衡点附近， 这里可以将系

2、统线性 化处理，将上述非线性系统进行泰勒展开得：y k y k 121!2y k1 231!3y k13n1!ny k 1 nu k 1 k(1.3)因为 01 ，后面 y k 1 的高阶项都可以扔掉(只作为寻找输入信号使用) ，所以系统可以化为下式：y k y k 1 u k 1 k (1.4) 此时不妨设系统输出 y k 的最后的极限为 A, 从式 1.2 得 A 0.75。那么应该满足lkim y k lkim y k 1 A (1.5)从而有1 A u k 1 k(1.6)同时为了抵消系统的部分噪声，这里采用 MATLAB 软件编程产生另一服从同一分布 的均匀噪声 1 k ，将式 1.

3、6 变形得：11u k 1 1u0 k 1 1 k (1.7)式 1.7 中 u0 k 是一个最后收敛于系统平衡点 0.75 的基本信号，这里可以采用一阶线性系统的阶跃响应曲线作为基本信号u0 k ，同时考虑系统的平衡点，即设计为：(1.8)u0 k 0.75 e k/TT 是一阶线性系统对应的时间常数， 反应到输入基本信号 u0 k 上就是过零点作 u0 k 对应1曲线的切线的斜率 K 。T阶跃信号系统输入基本信号 u0 kTs 1(1.9)(1.10) u0 k(1.11)T尽图 2 噪声随机抵消情况下的原系统响应曲线图1 基本信号产生框图很显然有lim u0 k 0.75t从而最后的输入

4、信号设计为：11u ku0 k 1 1 k 1式 1.10 中 , , 题中已知， 1 k 满足： E 1 k0 ，D 1 k0.01 ，的产生参照图 1。代入式 1.8 得：u k 1 (0.75 e (k 1)/T ) 1 1 k 1产生 u0 k 过程中的时间常数 T 由仿真效果决定，如果从系统的性能指标而言， 可能小，仿真如下：图 3 噪声完全抵消情况下的原系统响应曲线仿真分析：图 2 与图 3 中都是原非线性系统的响应曲线图，两者的输入信号有差别，分别为u1 k ,u2 k ， 如下式表示：11u1 ku0 k 1 1 k 1 (1.12)11u2 ku0 k 1 k 1 (1.13

5、)u0 k 参照式 1.8。式 1.12 中 1 k 1 是 MATLAB 中产生的符合同一均匀分布的给定随机噪声，从而式1.12 中输入信号 u1 k 实现噪声随机抵消，结果系统的响应 y k 呈现出很大的噪声干扰， 这是符合对未知系统认识的逻辑过程的， 如果要得出理想的响应曲线， 也就是要求系统最终 在平衡点附近趋于稳定， 那么需要设计相应的滤波器， 这里不再讨论， 后续会尝试设计滤波 器滤波。式 1.13 中 k 1 是系统给定的噪声，那么设计的输入信号 u k 实现了对噪声的完全抵消， 得到的理想响应曲线如图 3所示，从图 3中可以看出， y k 很快的收敛于 0.75 附 近，且上升

6、时间以及超调量等性能指标均符合理想响应要求， 这里就没有具体计算， 显然符 合控制系统的快稳准的三要素。(2) 最小二乘法辨识模型通过 MATLAB 编程从原非线性系统中采集一些输入输出数据，利用最小二乘法辨识出系统近似的线性模型，实现框图如图4。时不变 SISO 系统动态过程的数学模型为11A z 1 y k B z 1 u k k其中 u k 是系统的输入， y k 是系统的输出， A z 1 ,B z 1 的展开式分别如下：A z1 1 a1z1 a2z2anaz naB z 1 b1z 1 b2z 2bnb z nb参数 a1,a2, ,ana ,b1,b2, ,bnb 待辨识。式 2

7、.1 写成最小二乘格式y k hT k k式 2.4 中：h k y k 1, , y k na ,u k 1, ,u k nb TTa1,a2, ,ana ,b1,b2, ,bnb对于 k 1,2, ,L,方程 2.4构成一个线性方程组，可以把它写成式2.7，YL H LEL其中YLy 1 ,y 2 , ,y L TEL 1 , 2 , , L THLhT 1hT 2hT Ly 0y 1nau 0u 1nby 1y 2nau 1u 2nby L 1 y L na u L 1 u L nb（2.1）（2.2）（2.3）（2.4）（2.5）（2.6）（2.7）（2.8）（2.9）2.10）由题意

8、知 k 为均匀分布的噪声所以E ELE1E22.11)Cov EL E ELE LTE2 1E 1 2E1LE21E 2 2E2LE L 1 E L 2 E 2 LE22I2.12)2式中 2 0.01 ， I 是单位矩阵。这里认为噪声 k 与输入 u k 是不相关的，从而E k u k l 0, k, l2.13)这里需要考虑的是记忆长度的选择，上述矩阵方程包含（na nb） 个未知数，一般要求L na nb ，这样才可能确定一个“最优”的模型参数，而且为了保证辨识的精度， L 必须充分大，由 AIC 定阶准则知 L 同样决定着所辨识系统的阶次，下面会阐述相应理论。 定义准则函数：2.14)

9、LJ k y k hT k1其中 k 是加权因子，它衡量的是采样数据的可信度，这里对数据的可信度难以肯 定，就简单的将取 k 1。将准则函数写成二次型的形式如式2.15。J YL HL T Y L HL（2.15）设 wls 使得 J | minJ T T| YLHLYLH L| 0T（2.16）wls LL L L wls利用下面的两个向量微分公式展开得TTxxxT Ax 2xT Ax（2.17）其中 A 是对称阵，得到正则方程TTH LT H L wls H LTYL（2.18）当 HLTHL 为正定矩阵时，有：T 1 T wlsHLTH LHLTYL（2.19）从而J 2 |2HLT H

10、L 02 wls L L（ 2.20 ）结合式 2.16 到 2.20 知，wls 使J | min ，并且 wls是唯一的。从而求得最小二乘估计值T 1 T LSH LTH LH LTYL（ 2.21）定阶：系统的等价线性模型的阶次确定可以基于 AIC 准则来确定，首先将系统写成如式 2.7 的形式。输出变量在 n 条件下的似然函数为：L2.22）H n n （2.23）2 1 Texp 2 Yn H n nYn H n n2对应的对数似然函数为：2L L 2 1 Tln L n l n 2 ln 2 2 ln 2 Yn H n n Yn其中 L 是数据长度，根据极大似然原理得到模型参数 及

11、白噪声 k 方差的极大似然估计值如下：2.24）Yn H n MLT 1 TMLHnT HnHnTYn2 1 T1L Yn H n ML2将 ML 与 回代式 2.23 可得：NN2ln LML2ln LMLAIC2NAIC 准则的标准式为：l MLL2const ln22.25)2.26)将式 2.25 代入 AIC 准则的标准式 2.26 中得到：2AIC n Lln 4n2.27)面各式中 n表示的就是模型要估计的阶次， 最后找到使 AIC n min 的 n作为模型的阶次。本次辨识采样的数据长度L 1000, 很显然 足够小，当数据长度取定时，显然n可以尽可能的小，暂且取 n=4。仿真

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

13、辨识系统与原系 统都有很大的噪声干扰，这在工程应用中需要设计滤波器进行滤波。综合分析： 从设计合理的输入信号到最后辨识系统的模型，设计输入信号显得尤为重要，可持续激 励系统的各个工作点（平衡点）是系统可辨识的必要条件，从（2）中的辨识效果看，本文（1）中设计的输入信号 u1 k 并不是理想的，它显然满足持续激励系统的平衡点，但辨识 的模型泛化性并不理想， 随着输入信号的改变， 模型与过程的吻合程度有变化。 文中有很多 不解的问题，比如非线性系统的平衡点理论（系统差分方程表示下） ，数据长度的选择，初 衷是完整地从设计输入信号到辨识建模， 这其中包括定阶问题， 不同的辨识实现算法， 最后 到模型

14、的检验问题， 以及得到平滑的响应曲线需要加入滤波， 但由于本人水平有限， 有一些 问题还没有弄懂，最后没能实现，希望恩师不吝赐教。附录： 图2 程序噪声随机抵消情况下原系统的响应程序（数据长度 200）clear all y(1,1)=0;k = 2;x = (0.01:0.01:2);%采样频率 100HZ,200 个点%产生频率为 5HZ, 幅度为 20 的 sin 函数%产生 -0.5 到 0.5 的均匀白噪声， 200 个点%产生 -0.5 到 0.5 的均匀白噪声， 200 个点作为噪音u1 = 2.2125-2.95*exp(-10*x);u2 = -0.18+0.36*rand(

15、1,200);u3 = -0.18+0.36*rand(1,200);u = u1-4*u2;%输入信号 u(k)while (k = 201) y(1,k) = 0.75*(1-exp(-0.35*y(1,k-1)+0.25*u(1,k-1)+u3(1,k-1);k=k+1;end y = y(1,1:200);plot(y, b );title( 噪声随机抵消情况下的原系统响应曲线 ); %给图形添加标题 xlabel( k , FontName , Times New Roman, FontSize ,16);ylabel( y(k) , FontName , Times NewRoma

16、n , FontSize ,16, Rotation ,0);图3 程序噪声完全抵消情况下原系统的响应程序（数据长度 200）clear all ； y(1,1)=0;k = 2;%采样频率 100HZ,200 个点x = (0.01:0.01:2);u1 = 2.2125-2.95*exp(-10*x);%产生频率为 5HZ, 幅度为 20 的 sin 函数u2 = -0.18+0.36*rand(1,200);%产生 -0.5 到 0.5 的均匀白噪声， 200 个点u3 = -0.18+0.36*rand(1,200);%产生 -0.5 到 0.5 的均匀白噪声， 200 个点作为噪音u

17、 = u1-4*u3;%输入信号 u(k)while (k = 201)y(1,k) = 0.75*(1-exp(-0.35*y(1,k-1)+0.25*u(1,k-1)+u3(1,k-1);k=k+1;end y = y(1,1:200);plot(y, b );title( 噪声完全抵消情况下的原系统响应曲线 ); %给图形添加标题 xlabel( k , FontName , Times New Roman, FontSize ,16);ylabel( y(k) , FontName , Times NewRoman , FontSize ,16, Rotation ,0); 图5 程序

18、定 4阶情况下系统的辨识效果(数据长度1000)%数%据%生%成% clear ally(1,1)=0;k = 2;x = (0.01:0.01:10);%采样频率 100HZ,1000 个点u1 = 2.2125-2.95*exp(-10*x);u2 = -0.18+0.36*rand(1,1000);u3 = -0.18+0.36*rand(1,1000);%产生频率为 5HZ, 幅度为 20 的 sin 函数%产生 -0.5 到 0.5 的均匀白噪声，1000 个点作为噪%产生 -0.5 到 0.5 的均匀白噪声， 1000 个点u = u1-4*u2;%输入信号 u(k)while (

19、k = 1001)y(1,k) = 0.75*(1-exp(-0.35*y(1,k-1)+0.25*u(1,k-1)+u3(1,k-1); k=k+1;end y = y(1,1:1000);plot(y, b );title( 噪声随机抵消情况下的原系统响应曲线 ); %给图形添加标题 xlabel( k , FontName , Times New Roman, FontSize ,16);ylabel( y(k) , FontName , Times NewRoman , FontSize ,16, Rotation ,0);%最%小%二%乘法求解系统模型 % l = 990;%共 1000 个数据，从 -9 到 990n = 4; %定阶为 4yl = (y(1,11:1000);u3l = (u3(1,11:1000);i=1;while (i=990) hl(i,:) = y(1,i+9),y(1,i+8),y(1,i+7),y(1,i+6),u(1,i+9),u(1,i+8),u(1,i+7),u(1,i+6);i = i+1;endQ = inv(hl*hl)*hl*yl;k = 5;while (k = 1000)y_