系统辨识MATLAB仿真.doc

设一非线性系统如下所示: ,是零均值，方差为0.01的噪声序列（均匀白噪声）。（1） 试设计一种激励信号能持续激励系统的各工作点（平衡点）（2） 用适当的方法辨识出系统的等价模型（用另一组数据来检验模型的泛化性）说明：下面讨论的都是离散系统，所以时间坐标均采用离散时间节点。解：(1) 线性化处理寻找系统的合理输入信号可以求得系统的平衡点为： (1.1)按题意要求最后系统必须收敛于平衡点附近，即： (1.2) 为了找出系统的合理输入信号，使得系统最终工作在平衡点附近，这里可以将系统线性化处理，将上述非线性系统进行泰勒展开得：（1.3） 因为 ，后面的高阶项都可以扔掉（只作为寻找输入信号使用）， 所以系统可以化为下式：(1.4) 此时不妨设系统输出 的最后的极限为, 从式1.2得。 那么应该满足 (1.5)从而有 (1.6) 同时为了抵消系统的部分噪声，这里采用MATLAB软件编程产生另一服从同一分布的均匀噪声，将式1.6变形得： (1.7) 式1.7中是一个最后收敛于系统平衡点0.75的基本信号，这里可以采用一阶线性系统的阶跃响应曲线作为基本信号，同时考虑系统的平衡点，即设计为： (1.8)是一阶线性系统对应的时间常数，反应到输入基本信号上就是过零点作对应曲线的切线的斜率。系统输入基本信号阶跃信号 图1 基本信号产生框图 很显然有 (1.9)从而最后的输入信号设计为： (1.10)式1.10中题中已知，满足： ，的产生参照图1。代入式1.8得： (1.11) 产生 过程中的时间常数由仿真效果决定，如果从系统的性能指标而言，T尽可能小，仿真如下： 图2 噪声随机抵消情况下的原系统响应曲线图3噪声完全抵消情况下的原系统响应曲线仿真分析：图2与图3中都是原非线性系统的响应曲线图，两者的输入信号有差别，分别为，如下式表示： (1.12) (1.13) 参照式1.8。 式1.12中是MATLAB中产生的符合同一均匀分布的给定随机噪声，从而式1.12中输入信号实现噪声随机抵消，结果系统的响应呈现出很大的噪声干扰，这是符合对未知系统认识的逻辑过程的，如果要得出理想的响应曲线，也就是要求系统最终在平衡点附近趋于稳定，那么需要设计相应的滤波器，这里不再讨论，后续会尝试设计滤波器滤波。 式1.13中是系统给定的噪声，那么设计的输入信号实现了对噪声的完全抵消，得到的理想响应曲线如图3所示，从图3中可以看出，很快的收敛于0.75附近，且上升时间以及超调量等性能指标均符合理想响应要求，这里就没有具体计算，显然符合控制系统的快稳准的三要素。(2) 最小二乘法辨识模型 通过MATLAB编程从原非线性系统中采集一些输入输出数据，利用最小二乘法辨识出系统近似的线性模型，实现框图如图4。 +非线性系统MSE准则下 LS辨识-等价模型图4基于MSE的辨识过程实现框图时不变SISO系统动态过程的数学模型为 （2.1）其中是系统的输入，是系统的输出，的展开式分别如下： （2.2） （2.3）参数待辨识。式2.1写成最小二乘格式 （2.4）式2.4中： （2.5） （2.6）对于方程2.4构成一个线性方程组，可以把它写成式2.7， （2.7）其中 （2.8） （2.9） （2.10）由题意知为均匀分布的噪声所以 （2.11） （2.12）式中，是单位矩阵。这里认为噪声与输入是不相关的，从而 （2.13） 这里需要考虑的是记忆长度的选择，上述矩阵方程包含个未知数，一般要求，这样才可能确定一个“最优”的模型参数，而且为了保证辨识的精度，L必须充分大，由AIC定阶准则知同样决定着所辨识系统的阶次，下面会阐述相应理论。定义准则函数： （2.14） 其中是加权因子，它衡量的是采样数据的可信度，这里对数据的可信度难以肯定，就简单的将取。将准则函数写成二次型的形式如式2.15。 （2.15）设 使得 （2.16）利用下面的两个向量微分公式展开得 （2.17）其中是对称阵，得到正则方程 （2.18）当为正定矩阵时，有： （2.19）从而 （2.20）结合式2.16到2.20知，使，并且是唯一的。从而求得最小二乘估计值 （2.21）定阶： 系统的等价线性模型的阶次确定可以基于AIC准则来确定，首先将系统写成如式2.7的形式。输出变量在条件下的似然函数为： （2.22）对应的对数似然函数为： （2.23） 其中L是数据长度，根据极大似然原理得到模型参数及白噪声方差的极大似然估计值如下： （2.24） 将与回代式2.23可得： （2.25）AIC准则的标准式为： （2.26）将式2.25代入AIC准则的标准式2.26中得到： （2.27） 上面各式中表示的就是模型要估计的阶次，最后找到使的作为模型的阶次。 本次辨识采样的数据长度, 很显然足够小，当数据长度取定时，显然可以尽可能的小，暂且取=4。仿真分析： 图5 定4阶情况下系统的辨识效果图6 定4阶情况下更换输入为原来的情况下系统的辨识效果 图7 定4阶情况下更换输入为原来的4倍情况下系统的辨识效果图5与图6及图7中分别给出了定4阶情况下的系统利用最小二乘法辨识的结果，图5中的辨识效果在幅值上出现减小，有一定的误差，但基本可以与原系统同步，为了检验模型的泛化性，在MATLAB编程过程中，将输入信号减小为原来的以及增大为原来的4倍，图6中得到了较好的辨识效果，但图7却进一步增大幅度误差，从这里看出所辨识模型的泛化性并不好，即随着输入信号的增大，模型与系统的幅度差值增大。其次所辨识系统与原系统都有很大的噪声干扰，这在工程应用中需要设计滤波器进行滤波。综合分析： 从设计合理的输入信号到最后辨识系统的模型，设计输入信号显得尤为重要，可持续激励系统的各个工作点（平衡点）是系统可辨识的必要条件，从（2）中的辨识效果看，本文（1）中设计的输入信号并不是理想的，它显然满足持续激励系统的平衡点，但辨识的模型泛化性并不理想，随着输入信号的改变，模型与过程的吻合程度有变化。文中有很多不解的问题，比如非线性系统的平衡点理论（系统差分方程表示下），数据长度的选择，初衷是完整地从设计输入信号到辨识建模，这其中包括定阶问题，不同的辨识实现算法，最后到模型的检验问题，以及得到平滑的响应曲线需要加入滤波，但由于本人水平有限，有一些问题还没有弄懂，最后没能实现，希望恩师不吝赐教。附录：图2程序噪声随机抵消情况下原系统的响应程序（数据长度200）clear all；y(1,1)=0;k = 2;x = (0.01:0.01:2); %采样频率100HZ,200个点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 = 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;endy = y(1,1:200);plot(y,b);title(噪声随机抵消情况下的原系统响应曲线); %给图形添加标题xlabel(k,FontName,Times New Roman,FontSize,16);ylabel(y(k),FontName,Times New Roman,FontSize,16,Rotation,0);图3程序噪声完全抵消情况下原系统的响应程序（数据长度200）clear all；y(1,1)=0;k = 2;x = (0.01:0.01:2); %采样频率100HZ,200个点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 = 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;endy = y(1,1:200);plot(y,b);title(噪声完全抵消情况下的原系统响应曲线); %给图形添加标题xlabel(k,FontName,Times New Roman,FontSize,16);ylabel(y(k),FontName,Times New Roman,FontSize,16,Rotation,0);图5程序定4阶情况下系统的辨识效果（数据长度1000）%数据生成%clear all；y(1,1)=0;k = 2;x = (0.01:0.01:10); %采样频率100HZ,1000个点u1 = 2.2125-2.95*exp(-10*x); %产生频率为5HZ,幅度为20的sin函数u2 = -0.18+0.36*rand(1,1000); %产生-0.5到0.5的均匀白噪声，1000个点u3 = -0.18+0.36*rand(1,1000); %产生-0.5到0.5的均匀白噪声，1000个点作为噪音u = u1-4*u2; %输入信号u(k) while (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;endy = y(1,1:1000);plot(y,b);title(噪声随机抵消情况下的原系统响应曲线); %给图形添加标题xlabel(k,FontName,Times New Roman,FontSize,16);ylabel(y(k),FontName,Times New Roman,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_evalue(1,k) =Q(5,1)*u(1,k-1)+Q(6,1)*u(1,k-2)+Q(7,1)*u(1,k-3)+Q(8,1)*u(1,k-4)+u3(1,k)-(Q(1,1)*y(1,k-1)+Q(2,1)*y(1,k-2)+Q(3,1)*y(1,k-3)+Q(4,1)*y(1,k-4); k=k+1;endhold on;plot(y_evalue,r);xlabel(k,FontName,Times New Roman,FontSize,16);ylabel(y(k),FontName,Times New Roman,FontSize,16,Rotation,0);title(系统输出结果对比); %给图形添加标题%text(120,-1,红线代表估计系统输出);% 在图形的任意位置增加说明性文本信息%text(12