系统辨识理论及matlab仿真

系统辨识和最小二乘参数估计Matlab仿真一、系统辨识在控制系统的分析中，首先要建立系统的数学模型，控制系统的数学模型是定量描述系统或过程内部物理量（或变量）之间关系的数学表达式。

一般来说，建立控制系统数学模型有两种基本方法：（1）机理建模（白箱模型）：即根据系统内在运行机制、物料和能量守恒等物理学、化学规律建立系统的数学模型，一般步骤如下：Step1：根据系统工作原理及其在控制系统中的作用，确定输入和输出；Step2：根据物料和能量守恒等关系列写基本方程式；Step3：消去中间量；Step4：获得系统模型；（2）实验法建模（黑箱模型）：即对于机理尚不清楚或机理过于复杂的系统，可以人为的对其施加某种测试信号，并记录其输出响应，或者记录正常运行时的输入输出数据，然后利用这些输入输出数据确定系统模型结构和参数。

多年来，系统辨识已经发展为一门独立学科分支，通过系统辨识建立一个对象的数学模型，通常包括两方面的工作：一是模型结构的确定（模型的类型、阶次），二是模型参数估计。

根据时间是否连续，参数模型又可以分为连续时间系统模型和离散时间系统参数模型，这两类模型均可采用输入输出模型和状态空间模型描述，离散系统采用差分方程描述，以单输入单输出（SISO）离散系统参数模型为例。

1.确定性模型SISO系统确定性模型可表示为：u(k)和y(k)分别为输入和输出，d为纯延时。

2.随机性模型如果受到随机扰动，则式子可写为：为系统随机扰动，其结构如图：系统辨识的一般步骤如图：从图中可以看出，利用辨识的方法建立系统数学模型，从实验设计到模型获得，需要这些步骤。

二、最小二乘参数估计1.批处理最小二乘考虑以下CAR模型：为白噪声，结构参数na、nb和d已知，参数估计的任务就是根据可测量的输入输出，确定如下参数：仿真实例：式中，为方差为1的白噪声，选用幅值为1的逆M序列作为输入，LS算法进行参数估计，仿真结果如图：仿真程序(Matlab):%批处理最小二乘参数估计（LS）clear all;a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=500; %数据长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值x1=1; x2=1; x3=1; x4=0; S=1; %移位寄存器初值、方波初值xi=randn(L,1); %白噪声序列theta=[a(2:na+1);b]; %对象参数真值for k=1:Lphi(k,:)=[-yk;uk(d:d+nb)]'; %此处phi(k,:)为行向量，便于组成phi矩阵y(k)=phi(k,:)*theta+xi(k); %采集输出数据IM=xor(S,x4); %产生逆M序列if IM==0u(k)=-1;elseu(k)=1;endS=not(S); M=xor(x3,x4); %产生M序列%更新数据x4=x3; x3=x2; x2=x1; x1=M;for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endthetae=inv(phi'*phi)*phi'*y' %计算参数估计值thetae2.递推最小二乘在批处理最小二乘法时，由于每次处理的数据量较大，而且常常要求对象参数能够在线实时估计，解决的方法就是将其化成递推算法，其基本思想为：算法介绍：仿真实例:式中，为方差为0.1的白噪声，取初值，选择方差为1的白噪声作为输入信号u(k)，采用RLS算法进行参数估计，仿真结果如图：仿真程序（Matlab）：%递推最小二乘参数估计（RLS）clear all; close all;a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=400; %仿真长度uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)yk=zeros(na,1); %输出初值u=randn(L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); %白噪声序列theta=[a(2:na+1);b]; %对象参数真值thetae_1=zeros(na+nb+1,1); %thetae初值P=10^6*eye(na+nb+1);for k=1:Lphi=[-yk;uk(d:d+nb)]; %此处phi为列向量y(k)=phi'*theta+xi(k); %采集输出数据%递推最小二乘法K=P*phi/(1+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P;%更新数据thetae_1=thetae(:,k);for i=d+nb:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=na:-1:2yk(i)=yk(i-1);endyk(1)=y(k);endplot([1:L],thetae); %line([1,L],[theta,theta]); xlabel('k'); ylabel('参数估计a、b');legend('a_1','a_2','b_0','b_1'); axis([0 L -2 2]);。

基于Matlab系统辨识的参数辨识与仿真【摘要】论述了系统辨识的基本理论，分别用最小二乘法参数辨识和辅助变量法参数辨识。

根据Matlab系统辨识工具箱中的一些基本函数，结合实例来熟悉基于系统辨识工具箱的建模方法。

【关键词】Matlab;参数辨识;最小二乘法;辅助变量法1.系统辨识的基本理论系统辨识是根据系统的输入输出的时间函数来确定描述系统行为的数学模型，是现代控制理论中的一个分支。

对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。

它包括确定系统数学模型结构和估计其参数的方法。

系统辨识的流程如图1所示。

图1 系统辨识过程流程图2.模型参数辨识的方法系统辨识包括模型阶次辨识和参数辨识。

经典参数辨识的方法主要有他包括脉冲响应法、阶跃响应法、频率响应法、最小二乘法、相关分析法、谱分析法和极大似然法等，其中最小二乘法是最基本和最经典的，也是其他方法基本的思想的来源。

比如辅助变量法。

2.1 最小二乘法辨识考虑如下CAR模型：（1）参数估计的任务是根据可测量的输入和输出，确定如下个参数：对象（1）可以写成如下最小二乘形式：（2）现有L组输入输出观测数据：利用最小二乘法得到系统参数的估计值为：（3）2.2 辅助变量法辨识当为有色噪声时，利用最小二乘法进行参数辨识时往往得不到无偏一致的参数估计量。

在这个时候可以引入变量，然后利用最小二乘法进行辨识就可得到无偏一致的参数估计量。

因此，对于线性或本质线性系统，其过程的模型都可以化成最小二乘形式，考虑如下所示的模型方程：（4）将上式写成最小二乘格式，则得：假定存在一个辅助变量矩阵，维数与H相同，它满足以下极限特性：式中Q是非奇异矩阵。

如果辅助变量满足上述条件，则有：（5）图2 系统仿真图3.建模实例3.1 非参数模型辨识某被控对象的数学模型可以表示为：，式中：;为白噪声，编制MATLAB程序，分别对上述对象进行ARX建模和辅助变量法建模，并比较两种方法得到的脉冲响应。

系统辨识理论及matlab仿真系统辨识是一门评估和改进系统性能的研究领域，它利用外部观测数据对系统进行建模，使用这种模型来识别系统行为以及建议改进措施。

它结合了计算机科学、工程学、统计学和应用数学等多种学科的方法。

系统辨识的方法可以用于分析机器人系统、有限元素模型、分散系统和非线性系统等。

系统辨识理论是现代工程中最重要的技术，它能够有效地分析和模拟系统，以解决重大工程挑战。

系统辨识理论的基础可以回溯到20世纪40年代，当时大量研究的重点在于控制系统的设计与形成。

随着进步，它拓展到多种应用，其中包括一般系统建模和优化，以及系统健康监测，智能控制系统，复杂过程建模等等。

Matlab是一种基于矩阵的编程语言，它可以提供强大的工具来支持系统辨识理论的应用。

利用Matlab，可以非常方便地实现模型建模、数据处理、数值求解以及可视化等功能。

此外，它还提供了大量预定义函数，可以极大地简化系统辨识理论中所需要实现的功能，如参数估计、优化、数据验证等。

系统辨识理论和Matlab仿真在工程实践中应用十分广泛，其中最值得一提的有：（1）机器人控制：利用系统辨识理论和Matlab仿真，可以对机器人运动特性进行模型建模，并实现复杂的运动控制；（2）流体动力学：利用系统辨识理论和Matlab仿真，可以对流体动力学进行建模和模拟，以便改善和优化系统性能;（3）模糊控制：利用系统辨识理论和Matlab仿真，可以实现模糊控制系统的建模和仿真，用于实现智能控制等。

系统辨识理论和Matlab仿真的研究成果不仅可以改善系统性能，还可以用于展示系统原理及其实现。

以上是我关于系统辨识理论及Matlab仿真的研究结果及其应用。

经过以上介绍，可以看出系统辨识理论及Matlab仿真在工程应用中有重要的作用，它们为解决复杂工程问题提供了可行的解决方案。

然而，在实际应用过程中，系统辨识理论及Matlab仿真也存在一定的不足，比如误差控制方面的困难，这就要求我们不断改进和完善理论和技术，以提高系统的有效性和可靠性。

【1】随机序列产生程序【2】白噪声产生程序【3】M序列产生程序【4】二阶系统一次性完成最小二乘辨识程序【5】实际压力系统的最小二乘辨识程序【6】递推的最小二乘辨识程序【7】增广的最小二乘辨识程序【8】梯度校正的最小二乘辨识程序【9】递推的极大似然辨识程序【10】Bayes辨识程序【11】改进的神经网络MBP算法对噪声系统辨识程序【12】多维非线性函数辨识程序的Matlab程序【13】模糊神经网络解耦Matlab程序【14】F-检验法部分程序【1】随机序列产生程序A=6;x0=1;M=255;for k=1:100x2=A*x0;x1=mod (x2,M);v1=x1/256;v(:,k)=v1;x0=x1;v0=v1;endv2=vk1=k;%grapherk=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');title('(0,1)均匀分布的随机序列') 【2】白噪声产生程序A=6; x0=1; M=255; f=2; N=100;for k=1:Nx2=A*x0;x1=mod (x2,M);v1=x1/256;v(:,k)=(v1-0.5)*f;x0=x1;v0=v1;endv2=vk1=k;%grapherk=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');title('(-1,+1)均匀分布的白噪声')【3】M序列产生程序X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初T态（0101），Yi为移位寄存器各级输出m=60; %置M序列总长度for i=1:m %1#Y4=X4; Y3=X3; Y2=X2; Y1=X1;X4=Y3; X3=Y2; X2=Y1;X1=xor(Y3,Y4); %异或运算if Y4==0U(i)=-1;elseU(i)=Y4;endendM=U%绘图i1=ik=1:1:i1;plot(k,U,k,U,'rx')xlabel('k')ylabel('M序列')title('移位寄存器产生的M序列')【4】二阶系统一次性完成最小二乘辨识程序%FLch3LSeg1u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列z=zeros(1,16); %定义输出观测值的长度for k=3:16z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值endsubplot(3,1,1) %画三行一列图形窗口中的第一个图形stem(u) %画出输入信号u的经线图形subplot(3,1,2) %画三行一列图形窗口中的第二个图形i=1:1:16; %横坐标范围是1到16，步长为1plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线subplot(3,1,3) %画三行一列图形窗口中的第三个图形stem(z),grid on%画出输出观测值z的经线图形，并显示坐标网格u,z%显示输入信号和输出观测信号%L=14%数据长度HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10)u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13)u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵HL赋值ZL=[z(3);z(4);z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)]% 给样本矩阵zL 赋值%calculating parameters%计算参数c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 %计算并显示%DISPLAY PARAMETERSa1=c(1), a2=c(2), b1=c(3), b2=c(4) %从中分离出并显示a1 、a2、b1、b2%End【5】实际压力系统的最小二乘辨识程序%FLch3LSeg2clear%工作间清零V=[54.3,61.8,72.4,88.7,118.6,194.0]'%赋初值V，并显示P=[61.2,49.5,37.6,28.4,19.2,10.1]'%赋初值P，并显示%logP=-alpha*logV+logbeita=[-logV,1][alpha,log(beita)]'=HL*sita%注释P、V之间的关系for i=1:6;%循环变量的取值为从1到6Z(i)=log(P(i));%赋系统的输出采样值end%循环结束ZL=Z'%给zL赋值HL=[-log(V(1)),1;-log(V(2)),1;-log(V(3)),1;-log(V(4)),1;-log(V(5)),1;-log(V(6)),1] %给HL赋值%calculating parameters%计算参数c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c4=c2*c3%计算%Separation of Parameters%分离变量alpha=c4(1) % 为c4的第1个元素beita=exp(c4(2)) % 为以自然数为底的c4的第2个元素的指数【6】递推的最小二乘辨识程序%FLch3RLSeg3clear%清理工作间变量L=15;% M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值for i=1;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，if y(i)>0.5,u(i)=-0.03;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”else u(i)=0.03;%当M序列的值为"0"时,辨识的输入信号取“0.03”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号ufigure(1);%第1个图形stem(u),grid on%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零for k=3:15;%循环变量从3到15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%给出理想的辨识输出采样信号end%RLS递推最小二乘辨识c0=[0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=[c0,zeros(4,14)];%被辨识参数矩阵的初始值及大小e=zeros(4,15);%相对误差的初始值及大小for k=3:15; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; x=h1'*p0*h1+1; x1=inv(x); %开始求K(k)k1=p0*h1*x1;%求出K的值d1=z(k)-h1'*c0; c1=c0+k1*d1;%求被辨识参数ce1=c1-c0;%求参数当前值与上一次的值的差值e2=e1./c0;%求参数的相对变化e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列c0=c1;%新获得的参数作为下一次递推的旧参数c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵的最后一列p1=p0-k1*k1'*[h1'*p0*h1+1];%求出p(k)的值p0=p1;%给下次用if e2<=E break;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c%显示被辨识参数e%显示辨识结果的收敛情况%分离参数a1=c(1,; a2=c(2,; b1=c(3,; b2=c(4,; ea1=e(1,; ea2=e(2,; eb1=e(3,; eb2=e(4,; figure(2);%第2个图形i=1:15;%横坐标从1到15plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':') %画出a1，a2，b1，b2的各次辨识结果title('Parameter Identification with Recursive Least Squares Method')%图形标题figure(3); %第3个图形i=1:15; %横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题【7】增广的最小二乘辨识程序%FLch3ELSeg4clearL=60;%四位移位积存器产生的M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位积存器的输出初始值for i=1;x1=xor(y3,y4);%第一个移位积存器的输入信号x2=y1;%第二个移位积存器的输入信号x3=y2;%第三个移位积存器的输入信号x4=y3;%第四个移位积存器的输入信号y(i)=y4;%第四个移位积存器的输出信号，幅值"0"和"1"if y(i)>0.5,u(i)=-1;%M序列的值为"1"时,辨识的输入信号取“-1”else u(i)=1;%M序列的值为"0"时,辨识的输入信号取“1”endy1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号作准备endfigure(1);%画第一个图形subplot(2,1,1); %画第一个图形的第一个子图stem(u),grid on%画出M序列输入信号v=randn(1,60); %产生一组60个正态分布的随机噪声subplot(2,1,2); %画第一个图形的第二个子图plot(v),grid on;%画出随机噪声信号R=corrcoef(u,v);%计算输入信号与随机噪声信号的相关系数r=R(1,2);%取出互相关系数u%显示输入型号v%显示噪声型号z=zeros(7,60);zs=zeros(7,60);zm=zeros(7,60);zmd=zeros(7,60);%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出矩阵的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%给输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出赋初值%增广递推最小二乘辨识c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(7,7);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=5.0e-15;%取相对误差Ec=[c0,zeros(7,59)];%被辨识参数矩阵的初始值及大小e=zeros(7,60);%相对误差的初始值及大小for k=3:60; %开始求Kz(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);%系统在M序列输入下的输出采样信号h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';%为求K(k)作准备x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1; %Kd1=z(k)-h1'*c0; c1=c0+k1*d1;%辨识参数czs(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);%系统在M序列的输入下不考虑扰动时的输出响应zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)];%模型在M序列的输入下不考虑扰动时的的输出响应zmd(k)=h1'*c1;%模型在M序列的输入下的的输出响应e1=c1-c0;e2=e1./c0; %求参数的相对变化e(:,k)=e2;c0=c1;%给下一次用c(:,k)=c1;%把辨识参数c 列向量加入辨识参数矩阵p1=p0-k1*k1'*[h1'*p0*h1+1];%find p(k)p0=p1;%给下次用if e2<=E break;%若收敛情况满足要求，终止计算end%判断结束end%循环结束c, e, %显示被辨识参数及参数收敛情况z,zmd %显示输出采样值、系统实际输出值、模型输出值%分离变量a1=c(1,; a2=c(2,; b1=c(3,; b2=c(4,;%分离出a1、a2、b1、b2d1=c(5,; d2=c(6,; d3=c(7,; %分离出d1、d2、d3ea1=e(1,; ea2=e(2,; eb1=e(3,; eb2=e(4,; %分离出a1、a2、b1、b2的收敛情况ed1=e(5,; ed2=e(6,; ed3=e(7,; %分离出d1 、d2 、d3的收敛情况figure(2);%画第二个图形i=1:60;plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+')%画出各个被辨识参数title('Parameter Identification with Recursive Least Squares Method')%标题figure(3);%画出第三个图形i=1:60;plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed2,'r+')%画出各个参数收敛情况title('Identification Error')%标题%response%响应figure(4);%画出第四个图形subplot(4,1,1); %画出第四个图形中的四个子图的第一个子图i=1:60;plot(i,zs(i),'r')%画出被辨识系统在没有噪声情况下的实际输出响应subplot(4,1,2); %画出第四个图形中的四个子图的第二个子图i=1:60;plot(i,z(i),'g')%画出被辨识系统的采样输出响应subplot(4,1,3); %画出第四个图形中的四个子图的第三个子图i=1:60;plot(i,zmd(i),'b')%画出模型含有噪声的输出响应subplot(4,1,4); %画出第四个图形中的四个子图的第四个子图i=1:60;plot(i,zm(i),'b')%画出模型去除噪声后的输出响应【8】梯度校正的最小二乘辨识程序%FLch4GAeg1clearu=[-1,-1,-1,-1,1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,-1,-1,1];y=[0,-2,-6,-7,-7,-3,5,7,3,-1,5,3,-5,-3,1,-1,1,-5,-7,-7];%画出u和y图形figure(1), subplot(2,1,1), stem(u), subplot(2,1,2), stem(y), hold onk=1:20plot(k,y)%给出初始值h1=[-1,0,0]';h2=[-1,-1,0]';g=[0,0,0]';I=[1,0,0;0,1/2,0;0,0,1/4];h=[h1,h2,zeros(3,16)];%计算输入样本数据h(k)for k=3:18h(:,k)=[u(k),u(k-1),u(k-2)]';end%计算出权矩阵R(k)和估计值gfor k=1:18a=h(1,k)^2+(h(2,k)^2)/2+(h(3,k)^2)/4;%按照式（4.45)开始计算权矩阵a1=1/a;R=a1*I;%按照式（4.45）计算出权矩阵g(:,k+1)=g(:,k)+R*h(:,k)*(y(k+1)-h(:,k)'*g(:,k));%按照式（4.44)计算脉冲响应估计值end%画出图形g1=g(1,; g2=g(2,; g3=g(3,;figure(2)k=1:19;subplot(121)plot(k,g1,'r',k,g2,'g',k,g3,'b'),grid on%计算模型输出值ym及系统输出与模型输出之间的误差Eyfor k=1:18ym(k)=h(:,k)'*g(:,k); Ey(k)=y(k+1)-ym(k);endk=1:18;subplot(122)plot(k,Ey),grid ong, ym, Ey %显示脉冲响应估计值、模型输出值及系统输出与模型输出之间的误差figure(3)%画出脉冲响应曲线x=0:1:3;y=[0,g(1,18),g(2,18),g(3,18)];xi=linspace(0,3);yi=interp1(x,y,xi,'cubic');%三次插值plot(x,y,'o',xi,yi,'m'),grid on%画出脉冲响应估计值及其三次插值曲线【9】递推的极大似然辨识程序clear %清零a(1)=1;b(1)=0;c(1)=1;d(1)=0;u(1)=d(1);z(1)=0;z(2)=0; %初始化for i=2:1200 %产生m序列u(i)a(i)=xor(c(i-1),d(i-1));b(i)=a(i-1);c(i)=b(i-1);d(i)=c(i-1);u(i)=d(i);endu; %若取去‘；’可以在程序运行中观测到m序列v=randn(1200,1); %产生正态分布随机数V=0; %计算噪声方差for i=1:1200V=V+v(i)*v(i);endV1=V/1200;for k=3:1200 %根据v和u计算zz(k)=1.2*z(k-1)-0.6*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);endo1=0.001*ones(6,1);p0=eye(6,6); %赋初值zf(1)=0.1;zf(2)=0.1;vf(2)=0.1;vf(1)=0.1;uf(2)=0.1;uf(1)=0.1;%迭代计算参数值和误差值for k=3:1200h=[-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2)];hf=h;K=p0*hf*inv(hf'*p0*hf+1);p=[eye(6,6)-K*hf']*p0;v(k)=z(k)-h'*o1;o=o1+K*v(k) ;p0=p;o1=o;a1(k)=o(1);a2(k)=o(2);b1(k)=o(3);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=[-zf(k-1);-zf(k-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2)];endo1 %若取去‘；’可以在程序运行中观测到参数V1%绘图subplot(4,1,1)k=1:1200;plot(k,a1,'k:',k,a2,'b',k,b1,'r',k,b2,'m:',k,d1,'g',k,d2,'k');xlabel('k')ylabel('parameter')legend('a1=-1.2,','a2=0.6','b1=1.0','b2=0.5','d1=-1.0','d2=0.2'); %图标炷title('The parameter idendification of the RML');endsubplot(4,1,2)k=1:1200;plot(k,e1,'k',k,e2,'b',k,e3,'r',k,e4,'m',k,e5,'g',k,e6,'k');xlabel('k'); ylabel('error'); %title('误差曲线')endsubplot(4,1,3)k=1:1200;plot(k,u);xlabel('k'); ylabel('input'); %title('系统输入信号')endsubplot(4,1,4)k=1:1200;plot(k,v);xlabel('k'); ylabel('random noise'); %title('系统所加的随机噪声') end。

第二部分程序注释及运行结果读者须知:为了便于读者理解，现将光盘上第一部分可直接在MATLAB6.I 下运行的MATLAB程序的编号和书本上的内容对应如下，每个程序题目括号内的file.m是对应书本上的内容在光盘上第一部分的程序编号。

第二章的随机序列产生程序例2.1 用乘同余法产生随机数(见光盘FLch2sjxleg1.m)①编程如下：A=6; N=100; %初始化；x0=1; M=255;for k=1:N %乘同余法递推100次；x2=A*x0; %x2和x0分别表示x i和x i-1；x1=mod (x2,M); %将x2存储器的数除以M，取余数放x1（x i）中；v1=x1/256; %将x1存储器的数除以256得到小于1的随机数放v1中；）存放在矩阵存储器v的第k列中，v(:,k) v(:,k)=v1; % 将v1中的数（i%表示行不变、列随递推循环次数变化；x0=x1; %x i-1= x i；v0=v1;end %递推100次结束；v2=v %该语句末无‘；’，实现矩阵存储器v中随机数放在v2中，%且可直接显示在MA TLAB的window中；k1=k;%grapher %以下是绘图程序；k=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');tktle('(0-1)均匀分布的随机序列')②程序运行结果如图2.5所示。

图2.5 采用MA TLAB产生的(0,1)均匀分布的随机序列图③产生的(0-1)均匀分布的随机序列在程序运行结束后，产生的(0，1)均匀分布的随机序列，直接从MATLAB的window 界面中copy出来如下（v2中每行存6个随机数）：v2 =0.0234 0.1406 0.8438 0.0820 0.4922 0.96090.7852 0.7266 0.3750 0.2578 0.5508 0.31640.9023 0.4336 0.6094 0.6680 0.0234 0.14060.8438 0.0820 0.4922 0.9609 0.7852 0.72660.3750 0.2578 0.5508 0.3164 0.9023 0.43360.6094 0.6680 0.0234 0.1406 0.8438 0.08200.4922 0.9609 0.7852 0.7266 0.3750 0.25780.5508 0.3164 0.9023 0.4336 0.6094 0.66800.0234 0.1406 0.8438 0.0820 0.4922 0.96090.7852 0.7266 0.3750 0.2578 0.5508 0.31640.9023 0.4336 0.6094 0.6680 0.0234 0.14060.8438 0.0820 0.4922 0.9609 0.7852 0.72660.3750 0.2578 0.5508 0.3164 0.9023 0.43360.6094 0.6680 0.0234 0.1406 0.8438 0.08200.4922 0.9609 0.7852 0.7266 0.3750 0.25780.5508 0.3164 0.9023 0.4336 0.6094 0.66800.0234 0.1406 0.8438 0.0820第二章的白噪声产生程序例2.2 用乘同余法产生(见光盘FLch2bzsheg2.m)①编程如下：A=6; x0=1; M=255; f=2; N=100；%初始化；x0=1; M=255;for k=1: N %乘同余法递推100次；x2=A*x0; %分别用x2和x0表示x i+1和x i-1；x1=mod (x2,M); %取x2存储器的数除以M的余数放x1（x i）中；v1=x1/256; %将x1存储器中的数除以256得到小于1的随机数放v1中；）减去0.5再乘以存储器f中的系数，存放v(:,k)=(v1-0.5 )*f; %将v1中的数（i在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化；x0=x1; % x i-1= x i；v0=v1;end %递推100次结束；v2=v %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MA TLAB的window中；k1=k;%grapher %以下是绘图程序；k=1:k1;plot(k,v,k,v,'r');xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声')②程序运行结果如图2.6所示。

如何使用MATLAB进行系统辨识与模型建模引言：近年来，随着科学技术的飞速发展，各行各业都在努力寻求更高效、更智能的解决方案。

系统辨识与模型建模作为一种重要方法和工具，被广泛应用于控制系统、信号处理、机器学习等领域。

在这些领域中，MATLAB作为一款功能强大的数值计算软件，为我们提供了丰富的工具和函数，可用于进行系统辨识与模型建模的分析和实现。

本文将详细介绍如何使用MATLAB进行系统辨识与模型建模，并探讨其在实际应用中的意义和局限性。

一、系统辨识的基本原理1.1 系统辨识的概念及意义系统辨识是指通过对已有数据的分析和处理，建立描述该系统行为的数学模型的过程。

在实际应用中，系统辨识可以帮助我们了解系统的结构和特性，预测系统的行为，并为系统控制、优化提供依据。

1.2 系统辨识的方法系统辨识的方法主要包括参数辨识和结构辨识两种。

参数辨识是指通过拟合已知数据，确定数学模型中的参数值的过程。

常用的参数辨识方法有最小二乘法、极大似然估计法等。

结构辨识是指通过选择适当的模型结构和参数化形式，使用已知数据确定模型结构的过程。

常用的结构辨识方法有ARX模型、ARMA模型等。

二、MATLAB在系统辨识中的应用2.1 数据准备与预处理在进行系统辨识之前，我们首先需要准备好相关的数据。

数据的质量和数量对系统辨识的结果有着重要的影响，因此在数据准备阶段应尽量确保数据的准确性和完整性。

MATLAB提供了丰富的数据处理和分析函数，可用于数据预处理、数据清洗、数据归一化等操作，以提高数据的质量和可用性。

2.2 参数辨识的实现参数辨识是系统辨识的重要步骤之一，其主要目标是通过适当的数学模型拟合已知数据，确定模型中的参数值。

在MATLAB中，我们可以使用curve fitting工具箱中的函数，如fit、cftool等，来进行参数辨识的实现。

同时，MATLAB还提供了最小二乘法等常用的参数辨识算法，方便我们根据实际需求进行选择和应用。