matlab 系统识别练习题
Problem1
Parametric estimation of a resistance from current-voltage data.
Consider a resistor described by the following current-voltage relation:
V=Ri+
where R is the resistance,i is the current in the resistor,V is the voltage across the resistor terminals, and is a white noise with zero mean value.Suppose that the following data are available:
[i1i2:::i N]T:N=20current measurements equispaced in the interval[1;10]Ampere;
[V1V2:::V N]T:N=20noise-corrupted voltage measurements corresponding to the above current values.
Problem:Compute the following estimates of the resistance R:
Sample Mean estimate b R N.
Least Squares estimate b R LS.
Gauss-Markov estimate b R GM.
Bayesian estimate b R B.
Evaluate the quality of the obtained estimates by comparing on a plot the measured and“predicted”voltage.
Main steps:
(1)Load the data from the…le resistor_data_1.mat(load Matlab command).
(2)Plot on a…gure the measured data(plot Matlab command).
(3)Compute the Sample Mean estimate b R N.Plot in the…gure generated at step(2)the estimated resistor characteristic b V N=b R N i.
(4)Compute the Least Squares estimate b R LS(mldivide Matlab command or“n”Matlab command). Plot in the…gure generated at step(2)the estimated resistor characteristic b V LS=b R LS i. (5)Compute the Gauss-Markov estimate b R GM considering the disturbance variance matrix = 2 I N,where =5V olt and I N is the identity matrix.Plot in the…gure generated at step(2)the estimated resistor characteristic b V GM=b R GM i.
(6)Compute the Bayesian estimate as
b R B=R+ RV 1V V(V V)
where:
R=E[R]is the mean value of R.This mean value is assumed a priori(try with a value around3). V=Ri:
RV=E[(R R)(V V)T]=E[(R R)(Ri+ Ri)T]=E[(R R)(Ri Ri)T]=:::= RR i T. RR=E[(R R)(R R)]=variance of R.This variance is assumed a priori(try with the following values:10;1;0:1;0:01).
V V=E[(V V)(V V)T]=:::= RR ii T+ ,where = 2 I N.
Plot in the…gure generated at step(2)the estimated resistor characteristic b V B=b R B i.
(7)Repeat the steps(1)-(6)using the data stored in the…le resistor_data_2.mat.In steps(5)and (6),use a variance matrix consistent with the data.
Problem2
Estimation of the parameters of an auto-regressive system from output data.
Consider the following auto-regressive system:
y(t)=a1y(t 1)+a2y(t 2)+b+ (t)
where (t)is a white noise with zero mean value and variance 2 =40.Suppose that a set of L=500 data has been generated by this system.
Problem:Estimate the parameters a1;a2and b by means of Least https://www.360docs.net/doc/fe12361587.html,pare on a plot the measured and predicted output.
Main steps:
(1)Load the data from the…le AR_data.mat.
(2)De…ne the output vector y and the matrix according to
y=264y(3)...y(L)37
5 =
26
4y(2)y(1)1
..
.
..
.
..
.
y(L 1)y(L 2)1
37
5:
The measurement equation is y= + where =[a1a2b]T.
(3)Compute the Least Squares estimate b LS of .
(4)Compute the predicted output as y pred= b LS.Plot on a…gure y and y pred.
系统辨识最小二乘参数估计matlab
最小二乘参数估计 摘要： 最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程()Z l T l l T l ΦΦΦ-∧=1θ。 最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。 关键词： 最小二乘（Least-squares ），系统辨识（System Identification ） 目录： 1.目的 (1) 2.设备 (1) 3引言 (1) 3.1 课题背景 (1) 4数学模型的结构辨识 (2) 5 程序 (3) 5.1 M 序列子函数 ................................................................................. 错误！未定义书签。 5.2主程序............................................................................................... 错误！未定义书签。 6实验结果： ................................................................................................................................... 3 7参考文献： ................................................................................................. 错误！未定义书签。 1.目的 1.1掌握系统辨识的理论、方法及应用 1.2熟练Matlab 下最小二乘法编程 1.3掌握M 序列产生方法 2.设备 PC 机1台（含Matlab 软件） 3引言 3.1 课题背景 最小二乘理论是有高斯（K.F.Gauss ）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。 最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最
用matlab实现最小二乘递推算法辨识系统参数
用matlab实现最小二乘递推算法辨识系统参 数 自动化系统仿真实验室指导教师： 学生姓名班级计082-2 班学号撰写时间： 全文结束》》-3-1 成绩评定： 一．设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作； 3、了解并掌握Matlab中一些函数的作用与使用；二．设计要求最小二乘递推算法辨识系统参数，利用matlab编程实现，设初始参数为零。z(k)-1、5*z(k-1)+0、7*z(k-2)=1*u(k-1)+0、5*u(k-2)+v(k); 选择如下形式的辨识模型：z(k)+a1*z(k- 1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k);三．实验程序 m=3;N=100;uk=rand(1,N);for i=1:Nuk(i)=uk(i)*(-1)^(i-1);endyk=zeros(1,N); for k=3:N yk(k)=1、5*yk(k-1)-0、 7*yk(k-2)+uk(k-1)+0、5*uk(k-2); end%j=100;kn=0;%y=yk(m:j);%psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j- 2)];%pn=inv(psi*psi);%theta=(inv(psi*psi)*psi*y);theta=[0 ;0;0;0];pn=10^6*eye(4);for t=3:Nps=([yk(t-1);yk(t-
2);uk(t-1);uk(t-2)]);pn=pn- pn*ps*ps*pn*(inv(1+ps*pn*ps));theta=theta+pn*ps*(yk(t)-ps*theta);thet=theta;a1=thet(1);a2=thet(2);b1=thet(3);b2= thet(4); a1t(t)=a1;a2t(t)=a2;b1t(t)=b1;b2t(t)=b2;endt=1:N;plot(t,a 1t(t),t,a2t(t),t,b1t(t),t,b2t(t));text(20,1、 47,a1);text(20,-0、67,a2);text(20,0、97,b1);text(20,0、47,b2);四．设计实验结果及分析实验结果图：仿真结果表明，大约递推到第步时，参数辨识的结果基本到稳态状态，即a1=1、5999，b1=1，c1=0、5，d1=-0、7。五、设计感受这周的课程设计告一段落了，时间短暂，意义重大。通过这次次练习的机会，重新把matlab课本看了一遍，另外学习了系统辨识的有关内容，收获颇丰。对matlab的使用更加纯熟，也锻炼了自己在课本中搜索信息和知识的能力。在设计过程中虽然遇到了一些问题，但经过一次又一次的思考，一遍又一遍的检查终于找出了原因所在，也暴露出了前期我在这方面的知识欠缺和经验不足。同时我也进一步认识了matlab软件强大的功能。在以后的学习和工作中必定有很大的用处。
Matlab系统辨识尝试之详细过程1
Matlab系统辨识尝试之详细过程1 前面介绍了Matlab系统辨识工具箱的一些用法，这里拿一个直观的例子来尝试工具箱的具体用法。比较长，给个简单目录吧： 1.辨识的准备 2.辨识数据结构的构造 3.GUI辨识 4.辨识效果 5.对固有频率的辨识 6.结构化辨识 7.灰箱辨识 8.加入kalman滤波的灰箱辨识 1.辨识的准备 在辨识前，首先要根据自己辨识的情况，确定要辨识的状态空间模型的一些特点，如连续还是离散的；有无直通 分量（即从输入直通到输出的分量）；输入延迟；初始状态等。了解了这些情况就可以更快速的配置辨识时的一些设 置选项。 2.辨识数据结构的构造 使用原始数据构造iddata结构： data=iddata(y,u,Ts); 这里以一个弹簧质量系统的仿真为例 代码如下，其中用到了函数MDOFSolve，这在之前的博文介绍过（https://www.360docs.net/doc/fe12361587.html,/?p=183），拿来用即可。如果发现运行有错误，可以将MDOFSolve函数开头的一句 omega2=real(eval(omega2)); 注释掉。 %弹簧质量系统建模 clc clear close all m=200; k=980*1000;
c=1.5*1000; m1=1*m; m2=1.5*m; k1=1*k; k2=2*k; k3=k1; %%由振动力学知识求固有频率 M=[m10;0m2]; K=[k1+k2-k2;-k2k3+k2]; [omega,phi,phin]=MDOFSolve(M,K); fprintf('固有频率:%fHz\n',subs(omega/2/pi)); %%转化到状态空间 innum=2; outnum=2; statenum=4; A=[0100; -(k1+k2)/m10k2/m10; 0001; k2/m20-(k3+k2)/m20]; B=[00; 1/m10; 00; 01/m2]; C=[1000; 0010]; D=zeros(outnum,innum); K=zeros(statenum,innum); mcon=idss(A,B,C,D,K,'Ts',0);%连续时间模型 figure impulse(mcon) %%信号仿真，构造数据供辨识 n=511;%输入信号长度 Ts=0.001; t=0:Ts:(n-1)*Ts; u1=idinput(n,'prbs');%输入1为伪随机信号 u2=zeros(n,1);%输入2为空 u=[u1u2]; simdat=iddata([],u,Ts);%形成输入数据对象 e=randn(n,2)*1e-7; simopt=simOptions('AddNoise',true,'NoiseData',e);%添加噪声yn=sim(mcon,simdat,simopt);%加噪声仿真 y=sim(mcon,simdat);%无噪声仿真
Matlab_系统辨识_应用例子
例1、考虑仿真对象 )()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+-- 其中，)(k v 是服从正态分布的白噪声N )1,0(。输入信号采用4阶M 序列，幅度为1。选择如下形式的辨识模型 )()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+ 设输入信号的取值是从k =1到k =16的M 序列，则待辨识参数LS θ?为LS θ?=(T T -ΦΦΦ1)z 。其中，被辨识参数LS θ?、观测矩阵Φ的表达式为: ????? ???????=2121?b b a a LS θ (3)(4)(16)z z z ??????=???? ??z (2)(1)(2)((3)(2)(3)(2)(15)(14)(15)(14)z z u u z z u u z z u u --????--??Φ=????--?? 程序框图如图1所示。Matlab 仿真程序如下： %二阶系统的最小二乘一次完成算法辨识程序，文件名：LS.m
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画输入信号u的径线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(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)] %给样本矩阵 赋值
系统辨识及其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:100 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=v1; x0=x1; v0=v1; end v2=v k1=k; %grapher k=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:N x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*f; x0=x1;
v0=v1; end v2=v k1=k; %grapher k=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==0 U(i)=-1; else U(i)=Y4; end end M=U %绘图 i1=i k=1:1:i1; plot(k,U,k,U,'rx') xlabel('k') ylabel('M序列') title('移位寄存器产生的M序列') 【4】二阶系统一次性完成最小二乘辨识程序 %FLch3LSeg1 u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画出输入信号u的经线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线
系统辨识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变形得： ()()()0111 1u k u k k αβ εγ γ --= - (1.7) 式1.7中()0u k 是一个最后收敛于系统平衡点0.75的基本信号，这里可以采用一阶线性系统的阶跃响应曲线作为基本信号()0u k ，同时考虑系统的平衡点，即设计为： ()/00.75k T u k e -=- (1.8) T 是一阶线性系统对应的时间常数， 反应到输入基本信号()0u k 上就是过零点作()0u k 对应()()()()11y k y k u k k αβγε--=-+01αβ<（1.3） (1.4)
M序列辨识 Matlab程序 系统辨识
M序列发生程序 function [out]=Mseries(order) %M-series function % by Chen Hao,2011/12/16 for i=1:order x(i)=1; end for n=1:252 y(1:2)=x(order-1:order); for i=order:-1:2 x(i)=x(i-1); end x(1)=xor(y(1),y(2)); out(n)=1-2*x(order); end stem(out); end 计算程序 function [out]=addfunction( A,B,order,n ) %addfunction % by Chen Hao,2011/12/9 order=n*order; for i=1:order out(i)=0; end for k=0:order-1 for i=0:order-1 out(k+1)=out(k+1)+A(i-k+order+1)*B(i+1); end end for k=1:order out(k)=(out(k)-out(order))/(order+1); end end 主程序 function main() %Main function % by Chen Hao,2011/12/9 out=Mseries(6); time=1:252; simin=[time' out']; assignin('base','simin',simin); sim('Mseries_distinguish'); z=simout.signals.values(2:253); out=[out out]; g=addfunction(out,z,63,4); g=g(1:63); plot(g); num=[1 2]; den=[3 1 1]; s=tf(num,den); t=0:62; h=impulse(s,t);
matlab系统辨识
（System Identification Tool）系统辨识工具箱 早听说matlab博大精深，神通广大了，于是乎我确定肯定有更简单、直观、强大的工具来完成这小儿科把戏。查资料琢磨之后，我做了个小实验，在simulink里验证了该种方法。该方法的大原则是：在确定了系统的输入输出数据（两个列向量N×1形式，如果是1×N，会提示出错！）之后，设计好一定的辨识原则（比如说是2阶？3阶？，传递函数是零极点形式，还是带阻尼形式，等等），然后就交给强大的matlab，得到辨识结果。Step by step，plz！ Step1、建立模型获取系统输入输出数据 图1 图1系统的输入是阶跃信号，用Scope1监视，并输出到workspace （这步不会的自己百度哦），采样周期是0.1s，得到输入变量u（101×1的矩阵）；本人在系统的阶跃响应上叠加了一白噪声，当然也可以不加噪声，加了噪声就是期望更真实的模拟实际情况，白噪声参数设置见图2
图2 同样在Scope2监视，也将结果输出到workspace，得到响应数据y（同样也是101×1的矩阵） Step 2、进入辨识工具箱&设置辨识规则 直接在command window 输入ident，回车，进入辨识工具箱图3
图3 点击import下拉菜单，选时域数据time domain data，见图4 图4 在下图5红色圈区域输入之前得到的系统输入和输出数据，u和y
图5 在下图6绿色圈输入数据的一些信息，因为之前模型中，阶跃起点我是放在0s处的，这里也设置0，如果前面模型仿真是1s，这里应该也是1s；采样时间是0.1s，根据实际情况设置统一哦
系统辨识matlab最小二乘法
一、 实验题目：最小二乘法在系统辨识中的应用 二、 实验目的 1.掌握系统辨识的理论、方法及应用 2.熟练Matlab 下最小二乘法编程 3.掌握M 序列产生方法 三、 实验设备 1、 硬件设备：计算机配置，P4、32位CPU 、512M 内存 2、 软件设备： windows xp 操作系统 、matlab6.5软件包 四、 实验原理 最小二乘理论是有高斯（K.F.Gauss ）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”。 单输入单输出离散时间动态系统差分方程为： )()()()k (1i 1i k e i k u i k Z Z bn i na i b a +-=-+∑∑==其中Z （k ）为输出变量，u(k)为输入变量，e(k)为偏差。 上式可以表示为)()()(-)k (1i 1i k e i k u i k Z Z bn i na i b a +--=∑∑==各参数用矩阵表示 ????????????--------------=)()1()()1()2()1()2()1()1() 0()1()0(a a b a b a n l u l u n l z l z n u u n z z n u u n z z H （1） T l z z z Z )](),2(),1([，?= （2） 其中l 为所需要采集的点数。 []T n e E )(........).........0(e = （3） []T bn anb a ........1.........1=θ （4） Z=H*θ+E ，E=Z-H*θ，根据最小二乘理论E 必须最小对上式进行求导，推出 Z H H H T T 1)(-=θ 根据表达式Z H H H T T 1)(-=θ带入（1）（2）（4）即可求出a1....an b1.......bn 。
系统辨识最小二乘参数估计matlab
最小二乘参数估计 摘要： 最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程()Z l T l l T l ΦΦΦ-∧=1θ。 最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。 关键词： 最小二乘（Least-squares ），系统辨识（System Identification ） 目录： 1.目的 (2) 2.设备 (2) 3引言 (2) 课题背景 (2) 4数学模型的结构辨识 (3) 5 程序 (4) M 序列子函数 (4) 主程序 (5) 6实验结果： (7) 7参考文献： (7) 1.目的 掌握系统辨识的理论、方法及应用 熟练Matlab 下最小二乘法编程 掌握M 序列产生方法 2.设备 PC 机1台（含Matlab 软件） 3引言 课题背景 最小二乘理论是有高斯（）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。 最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最好拟合的数学模型。递推最小二乘法是在最小二乘法得到的观测数据的基础上，用新引入的
数据对上一次估计的结果进行修正递推出下一个参数估计值，直到估计值达到满意的精确度为止。 4数学模型的结构辨识 根据汉格尔矩阵估计模型的阶次 设一个可观可控的SISO 过程的脉冲响应序列为{个g(1),g(2),……g(L)},可以通过汉格尔（Hankel ）矩阵的秩来确定系统的阶次。 令Hankel 阵为： ????? ???????-++-++++-++=)22()1()1()()2()1()1()1()(),(l k g k g l k g l k g k g k g l k g k g k g k l H ，其中l 决定),(k l H 阵地维数，k 可在1至()22+-l L 间任意选择。则有[]k n l n k l H rank ?≥=,,),(00。 如果0n l ≥（过程的真实阶次），那么Hankel 阵的秩等于0n 。因此可以利用Hankel 阵的奇异性来确定系统的阶次0n 。 根据残差平方和估计模型的阶次 SISO 过程的差分方程模型的输出残差为)(~k z ，数据长度L ，n H ?为n ?阶时的数据矩阵，n ??θ为n ?阶时的参数的估计量，n ?为模型阶次估计值，0n 为真实阶次，则残差平方和函数J ： )(~1)?()?(1~~1)?(1 2??????00k z L H z H z L z z L n J L n n k n n n T n n n n T n ∑++==--==θθ 残差平方和有这样的性质：当L 足够大时，随着n ?增加)?(n J 先是显著地下降，当n ?>0n 时，)?(n J 值显著下降的现象就终止。这就是损失函数法来定阶的原理。
系统辨识的Matlab实现方法(手把手)
最近在做一个项目的方案设计，应各位老总的要求，只有系统框图和器件选型可不行，为了凸显方案设计的高大上，必须上理论分析，炫一下“技术富”，至于具体有多大实际指导意义，那就不得而知了！本人也是网上一顿百度，再加几日探索，现在对用matlab 实现系统辨识有了一些初步的浅薄的经验，在此略做一小节。 必须要指出的是，本文研究对象是经典控制论理最简单最常用的线性时不变的siso 系统，而且是2阶的哦，也就是具有如下形式的传递函数： 1 21)(2 2++=Ts s T s G ξ 本文要做的就是，对于有这样传递函数的一个系统，要辨识得到其中的未知数T ， ξ!!这可是控制系统设计分析的基础哦，没有系统模型，啥理论、算法都是白扯，在实际工程中非常重要哦！ 经过总结研究，在得到系统阶跃响应实验数据之后（当然如果是其他响应，也有办法可以辨识，在此还是只讨论最简单的阶跃响应实验曲线，谁让你我是菜鸟呢），利用matlab 至少可以有两种方法实现实现（目前我只会两种，呵呵）！ 一、函数法 二、GUI 系统辨识工具箱 下面分别作详细介绍！
一、 函数法 看官别着急，先来做一段分析（请看下面两排红*之间部分），这段分析是网上找来的，看看活跃一下脑细胞吧，如果不研读一下，对于下面matlab 程序，恐怕真的就是一头雾水咯！ ******************************************************************************* G(s)可以分解为：) )((1)(212ωω++=s s T s G 其中， [] [] 1 1 1 1 2221--=-+=ξξωξξωT T 1ω、2ω都是实数且均大于零。 则有： 211 ωω=T ，2 12 12ωωωωξ+= 传递函数进一步化为： ) )(()(212 1ωωωω++= s s s G 因此，辨识传递函数就转化为求解1ω、2ω。 当输入为单位阶跃函数时，对上式进行拉普拉斯反变换，得系统时域下的单位阶跃响应为： t t e e t y 21 2111 221)(ωωωωωωωω---+ -- =
系统参数辨识+matlab+实现
4. 设某物理量Y 与X 满足关系式Y=aX 2+bX+c ，实验获得一批数据如下表，试辨识模型参数a ，b 和c 。（50分） X 1.01 2.03 3.02 4.01 5 6.02 7.03 8.04 9.03 10 Y 9.6 4.1 1.3 0.4 0.05 0.1 0.7 1.8 3.8 9.0 报告要求：要有问题描述、参数估计原理、程序流程图、程序清单，最后给出结果及分析。 (1)问题描述： 由题意知，这是一个已知模型为Y=aX 2+bX+c ，给出了10组实验输入输出数 据，要求对模型参数a ，b ，c 进行辨识。这里对该模型参数辨识采用递推最小二乘法。 (2)参数估计原理 对该模型参数辨识采用递推最小二乘法，即RLS （ recurisive least square ）， 它是一种能够对模型参数进行在线实时估计的辨识方法。 其基本思想可以概括为：新的估计值)(?k θ=旧的估计值)1(?-k θ+修正项 下面将批处理最小二乘法改写为递推形式即递推最小二乘参数估计的计算方法。 批处理最小二乘估计θ?为Y T T ΦΦΦ=-1 ) (?θ，设k 时刻的批处理最小二乘估计为： k T k k T k Y ΦΦΦ=-1 )(?θ令1 11 )] 1()()1([) ()(----+-=ΦΦ=k k k P k P T k T k ?? K 时刻的最小二乘估计可以表示为 k T k Y k P k Φ=)()(?θ=)]()()[(1 1k y k Y k P k T k ?+Φ-- =)]1(?)()()[()1(?--+-k k k y k K k T θ?θ ；式中)()()(k k P k K ?=，因为要 推导出P(k)和K(k)的递推方程，因此这里介绍一下矩阵求逆引理：设A 、（A+BC ）和（I +B CA 1 -）均为非奇异方阵，则1 11 11 1 )() (------+-=+CA B CA I B A A BC A 通过运用矩 阵求逆引理，把复杂的矩阵求逆转化为标量求倒数，大大减小了计算量。1+N P 与N P 间的递推关系。最终得到递推最小二乘参数递推估计公式如下： ??? ()(1)()[()()(1)]T k k K k y k k k θθ?θ=-+-- ① ()[() ()](T P k I K k k P k ?=-- ②
系统辨识及其matlab仿真(一些噪声和辨识算法)
1】随机序列产生程序 2】白噪声产生程序 【3】M序列产生程序 4】二阶系统一次性完成最小二乘辨识程序 5】实际压力系统的最小二乘辨识程序 6】递推的最小二乘辨识程序 7】增广的最小二乘辨识程序 8】梯度校正的最小二乘辨识程序 9】递推的极大似然辨识程序 10 】Bayes 辨识程序 【11】改进的神经网络MBP算法对噪声系统辨识程序 1 2 】多维非线性函数辨识程序的Matlab 程序 13】模糊神经网络解耦Matlab 程序 14】F- 检验法部分程序 【1】随机序列产生程序 A=6; x0=1;M=255; for k=1:100 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=v1; x0=x1; v0=v1; end v2=v k1=k; %grapher k=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:N x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*f; x0=x1; v0=v1; end
v2=v k1=k; %grapher k=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==0 U(i)=-1; else U(i)=Y4; end end M=U %绘图 i1=i k=1:1:i1; plot(k,U,k,U,'rx') xlabel('k') ylabel('M 序列 ') ti tle(' 移位寄存器产生的M序列') 【 4】二阶系统一次性完成最小二乘辨识程序 %FLch3LSeg1 u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; % 系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); % 定义输岀观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); % 用理想输岀值作为观测值 end subplot(3,1,1) % 画三行一列图形窗口中的第一个图形 stem(u) % 画岀输入信号 u 的经线图形 subplot(3,1,2) % 画三行一列图形窗口中的第二个图形
系统辨识MATLAB
系统辨识作业 方法一: 递推最小二乘法 [NUM]=xlsread('shuju','B4:B257'); u=NUM; [NUM1]=xlsread('shuju','C4:C257'); z=NUM1; N=length(u); c0=[0.001,0.001,0.001,0.001]'; %直接给出被辨识参数的初始值，取一个充分小的实向量 p0=10^7*eye(4,4); %初始状态P0也采用直接取方式，取一个充分大的实数单位矩阵 E=0.00000005; %相对误差E参考值 c=[c0,zeros(4,253)]; %被辨识参数矩阵的初始值及大小 e=zeros(4,254); %相对误差的初始值及大小 %开始递推运算 for k=3:N; h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]'; %求h1 x=h1'*p0*h1+1; x1=inv(x); k1=p0*h1*x1;%求k1 c1=c0+k1*(z(k)-h1'*c0); %求c e1=(c1-c0)./c0;%求参数的相对变化 e(:,k)=e1; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把当前所辨识参数的c1列向量加入辨识参数矩阵的最后一列 p1=p0-k1*h1'*p0; %求p1值 p0=p1; %把当前值给下次用 if norm(e1)<=E break; %若参数收敛满足要求，终止计算 end end %分离参数 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(1); %画第1个图形 i=1:254; %横坐标从1到254 plot(i,a1,'k',i,a2,'b',i,b1,'r',i,b2,'g') %画出a1，a2，b1，b2的各次辨识结果
Matlab_系统辨识_应用例子
例1、考虑仿真对象 z(k) 1.5z(k 1) 0.7z(k 2) u(k 1) 0.5u(k 2) v(k) 其中，v(k)是服从正态分布的白噪声N(0,1)。输入信号采用4阶M序列，幅度为1。选择如下形式的辨识模型 z(k) a1z(k 1) a2z(k 2) b1u(k 1) pu(k 2) v(k) 设输入信号的取值是从k =1到k =16的M序列，则待辨识参数?LS 为?LS=( T ) 1 T z。其中，被辨识参数？LS、观测矩阵的表达式为: a1 z(3) ? a2 址s . z(4) z L b1 b2 z( 16 ) z(2) z(1) u(2) u(1) z(3) z(2) u(3) u(2) L L z(15) z(14)u(15)u(14)程序框图如图1所示。Matlab仿真程序如下: %二阶系统的最小二乘一次完成算法辨识程序，文件名: 图1最小二乘一次完成算法程序框图
u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周 期的M 序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=*z(k-1)*z(k-2)+u(k-1)+*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画输入信号u 的径线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(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)
系统辨识试验(苍松书屋)
2.用普通最小二乘法（OLS ）法辨识对象数学模型 选择的仿真对象的数学模型如下 )()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+-- 其中，)(k v 是服从正态分布的白噪声N )1,0(。输入信号采用4阶M 序列，幅度为1。选择如下形式的辨识模型 )()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+ 设输入信号的取值是从k =1到k =16的M 序列，则待辨识参数LS θ?为LS θ?=L τL 1L τL z H )H H -(。其中，被辨识参数LS θ?、观测矩阵z L 、H L 的表达式为 ? ? ??? ? ??????=2121? b b a a LS θ , ????????????=)16()4()3(z z z L z ， ????????????------=)14()2()1()15()3()2()14()2()1()15()3()2(u u u u u u z z z z z z L H 程序框图如下所示： 赋输入信号初值u 定义输出观测值的长度并计算系统的输出值 画出输入和输出观测值的图形 给样本矩阵H L 和z L 赋值 根据公式计算参数LS θ? 从LS θ?中分离出并显示出被辨识参数a 1, a 2, b 1, b 2 停机 图2 最小二乘一次完成算法程序框图
参考程序： %ols u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M 序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画出输入信号u的经线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(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 PARAMETERS a1=c(1), a2=c(2), b1=c(3), b2=c(4) %从中分离出并显示a1 、a2、 b1、 b2
多变量系统辨识matlab程序
多变量系统辨识matlab程序 y(i)=0.05129*u1_1+0.0418；u1_3=u1_2;u1_2=u1_1;u1_1；u2_3=u2_2;u2_2=u2_1;u2_1；y_3=y_2;y_2=y_1;y_1=y(i)；r_3=r_2;r_2=r_1;r_1=r(i)；end；plot(time,y,'b')；holdon；xi=y';；savesub. y(i)=0.05129*u1_1+0.0418*u2_1+0.6386*y_1+0.06268*u1_2+0.0346*u2_2 -0.1179*y_2-0.004184*u1_3-0.00218*u2_3+0.006738*y_3+0.091*r_1-0.114*r _2+0.0509*r_3; u1_3=u1_2;u1_2=u1_1;u1_1=u1(i); u2_3=u2_2;u2_2=u2_1;u2_1=u2(i); y_3=y_2;y_2=y_1;y_1=y(i); r_3=r_2;r_2=r_1;r_1=r(i); end plot(time,y,'b') hold on xi=y'; save sub.txt xi –ascii 程序5 clear %CRA模型基于模型阶次递增的辨识。 clc close all
z=load('sub.txt'); u1=load('prbs1.txt'); u2=load('prbs2.txt'); for i=1:1:100 H(i,:)=[u1(20+i-1) u2(20+i-1) -1*z(20+i-1)]; end theta=(1e-3)*ones(3,1); P=(1e8)*eye(3); for i=1:1:100 K=P*H(i,:)'./(H(i,:)*P*H(i,:)'+1); theta=theta+K*(z(i+20)-H(i,:)*theta); P=(eye(3)-K*H(i,:))*P; end theta1=theta H1=H; J(1)=(z(21:120)-H1*theta1)'*(z(21:120)-H1*theta1); ZZ=inv(H1'*H1); %************************** for n=2:1:10 for i=1:1:100