系统辨识作业1

系统辨识大作业
已知一系统为两输入单输出系统，观测数据受有色噪声污染，噪信比为N/S=0.1。

系统经2000次采样，存放于文件T3T.TXT中。

系统输入u1为7级M序列，u2为u1的63步移位序列。

模型类可选为：A(q-1)y(k)=B1(q)u1(k)+ B2(q)u2(k)+w(k)/C(q-1)。

要求编制程序，辨识出该模型的结构及参数。

作业文档要求：
描述问题；
选择辨识方法并简单说明所选方法中的结构辨识原理和参数估计原理；
程序流程图及程序清单；
说明程序中用到的一些技术，如数据标准化、UD分解、稳定性判断等；
结构搜索路线及各结构下的参数、残差；
给出最终结果：A(q-1)=
B1(q)=
B2(q)=
C(q-1)=
并给出选择此最终结果的理由；
用你的辨识结果来预报系统输出误差e(k)=y(k)-y(k)，并画出e(101)-e(400)的曲线图。

《系统辨识》大作业学号：********班级：自动化1班姓名：**信息与控制工程学院自动化系2015-07-11第一题模仿index2，搭建对象，由相关分析法，获得脉冲响应序列ˆ()g k，由ˆ()g k，参照讲义，获得系统的脉冲传递函数()G z和传递函数()G s;应用最小二乘辨识,获得脉冲响应序列ˆ()g k;同图显示两种方法的辨识效果图；应用相关最小二乘法，拟合对象的差分方程模型；构建时变对象，用最小二乘法和带遗忘因子的最小二乘法，(可以用辨识工具箱) 辨识模型的参数，比较两种方法的辨识效果差异;答：根据index2搭建结构框图：相关分析法：利用结构框图得到UY 和tout其中的U就是题目中要求得出的M序列，根据结构框图可知序列的周期是1512124=-=-=npN。

在command window中输入下列指令，既可以得到脉冲相应序列()g k：aa=5;NNPP=15;ts=2; RR=ones(15)+eye(15); for i=15:-1:1UU(16-i,:)=UY(16+i:30+i,1)'; endYY=[UY(31:45,2)];GG=RR*UU*YY/[aa*aa*(NNPP+1)*ts]; plot(0:2:29,GG) hold onstem(0:2:29,GG,'filled') Grid;title('脉冲序列g(τ)')最小二乘法建模的响应序列由于是二阶水箱系统，可以假设系统的传递函数为221101)(sa s a sb b s G +++=，已知)(τg ，求2110,,,a a b b已知G （s ）的结构，用长除法求得G(s)的s 展开式，其系数等于从 )( g 求得的各阶矩，然后求G(s)的参数。

得到结果： a1 =-1.1561 a2 =0.4283 b0 =-0.0028 b1=0.2961在command window 中输入下列指令得到传递函数：最小二乘一次算法相关参数%最小二乘法一次完成算法 M=UY(:,1); z=UY(:,2); H=zeros(100,4); for i=1:100 H(i,1)=-z(i+1); H(i,2)=-z(i); H(i,3)=M(i+1); H(i,4)=M(i); endEstimate=inv(H'*H)*H'*(z(3:102)) %结束得到相关系数为：Estimate =-0.7866 0.1388 0.5707 0.3115带遗忘因子最小二乘法：%带遗忘因子最小二乘法程序M=UY(:,1);z=UY(:,2);P=1000*eye(5); %Theta=zeros(5,200); %Theta(:,1)=[0;0;0;0;0];K=zeros(4,400); %K=[10;10;10;10;10];lamda=0.99;%遗忘因数for i=3:201h=[-z(i-1);-z(i-2);M(i);M(i-1);M(i-2)];K=P*h*inv(h'*P*h+lamda);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(5)-K*h')*P/lamda;endi=1:200;figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:) )title('带遗忘因子最小二乘法')grid%结束Estimate 可由仿真图得出，可知两种方法参数确定十分接近。

论文系统辨识姿态角控制1.系统辨识概述辨识、状态估计和控制理论是现代控制理论三个相互渗透的领域。

辨识和状态估计离不开控制理论的支持，控制理论的应用又几乎不能没有辨识和状态估计技术。

随着控制过程复杂性的提高，控制理论的应用日益广泛，但其实际应用不能脱离被控对象的数学模型。

然而在大多数情况下，被控对象的数学模型是不知道的，或者在正常运行期间模型的参数可能发生变化，因此利用控制理论去解决实际问题时，首先需要建立被控对象的数学模型。

系统辨识正是适应这一需要而形成的，他是现代控制理论中一个很活跃的分支。

社会科学和自然科学领域已经投入相当多的人力去观察、研究有关的系统辨识问题。

系统辨识是建模的一种方法，不同的学科领域，对应着不同的数学模型。

从某种意义上来说，不同学科的发展过程就是建立他的数学模型的过程。

辨识问题可以归结为用一个模型来表示可观系统（或将要改造的系统）本质特征的一种演算，并用这个模型吧对客观系统的理解表示成有用的形式。

当然可以刻有另外的描述，辨识有三个要素：数据，模型类和准则。

辨识就是按照一个准则在一组模型类中选择一个与数据拟合得最好的模型。

总而言之，辨识的实质就是从一组模型类中选择一个模型，按照某种准则，使之能最好地拟合所关心的实际过程的静态或动态特性。

通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。

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

对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。

而系统辨识所研究的问题恰好是这些问题的逆问题。

通常，预先给定一个模型类{}M（即给定一类已知结构的模型），一类输入信号u和等价准则(,)JLyyM（一般情况下，J是误差函数，是过程输出y和模型输出yM的一个泛函）；然后选择是误差函数J达到最小的模型，作为辨识所要求的结果。

1.1讨论系统方程为：（信噪比）至少30%)()1()()1()(k e k u k b k ay k y +-=-+)(k e 为零均值白噪声， =)(k a =)(k b要求对系统参数辨识进行讨论（1） 定常系统 a=0.8,b=0.5参数递推估计 （2） 时变系统取不同值是的不同结果并讨论。

（1）取初值P (0)=I 、∧θ（0）=0选择方差为1的白噪声作为输入信号u （k ），L=300，采用RLS 算法进行参数估计，代码及仿真结果图如下：clear all ;close all ;a=[1 0.8]';b=[0.5]';d=1; %对象参数na=length(a)-1;nb=length(b)-1; %na 、nb 为A 、B 阶次L=300; %仿真长度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:2 uk(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','b_0');axis([0 L -2 2]);图1-1 递推最小二乘法的参数估计结果（2））取初值P(0)=I、∧θ（0）=0选择方差为1的白噪声作为输入信号u（k），取分别为0.91,0.95,0.98,1.00时，L=600，采用FFRLS算法进行参数估计，代码及仿真图如下所示：clear all;close all;a=[1 0.8]';b=[0.5]';d=1; %对象参数na=length(a)-1;nb=length(b)-1; %na、nb为A、B阶次L=600; %仿真长度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); %白噪声序列thetae_1=zeros(na+nb+1,1); %thetae初值P=10^6*eye(na+nb+1);lambda=0.98; %遗忘因子范围[0.9 1]for k=1:Lif k==301a=[1 0.6]';b=[0.3]'; %对象参数突变endtheta( : ,k)=[a(2 : na+1);b]; %对象参数真值phi = [-yk;uk(d:d+nb)]; %´此处phi为列向量y(k)=phi'*theta( : ,k) + xi(k); %采集输出数据%遗忘因子递推最小二乘法·¨K=P*phi/(lambda+phi'*P*phi);thetae( : ,k)=thetae_1 + K*(y(k)-phi'*thetae_1);P=(eye(na+nb+1)-K*phi')*P/lambda;%更新数据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);endsubplot(1,2,1)plot([1:L],thetae(1:na,:));hold on;plot([1:L],theta(1:na,:),'k:'); xlabel('k');ylabel('参数估计a');legend('a_1');axis([0 L -2 2]);subplot(1,2,2)plot([1:L],thetae(na+1:na+nb+1,:));holdon;plot([1:L],theta(na+1:na+nb+1,:),'k:');xlabel('k');ylabel('参数估计b');legend('b_0');axis([0 L -0.5 2]);图1-2-1 遗忘因子递推最小二乘法的参数估计结果（λ=0.91）图1-2-2 遗忘因子递推最小二乘法的参数估计结果（λ=0.95）图1-2-3遗忘因子递推最小二乘法的参数估计结果（ =0.98）图1-2-4 遗忘因子递推最小二乘法的参数估计结果（λ=1）由以上可以看出，技术对于参数突变的系统，采用FFRLS算法也能够有效地进行参数估计。

系统辨识大作业
一.设SlSO系统差分方程为
y(k)=—α1y(k-1)-a2y(k-2)+bλu(k-1)+b2u(k-2)+ξ{k)
辨识参数向量为θ=[q a2b l b2]r,输入输出数据详见数据文件UyLtXt—uy3.txtoξ(k)为噪声方差各异的白噪声或有色噪声。

试求解：
1)用n元一次方程解析法，再求其平均值方法估计。

2)用最小二乘及递推最小二乘法估计。

；
3)用辅助变量法及其递推算法估计
4)用广义最小二乘法及其递推算法估计
5)用夏氏偏差修正法、夏氏改良法及其递推算法估计
6)用增广矩阵法估计
7)分析噪声父攵)特性；
二.用极大似然法估计6。

三.以上题的结果为例，进行：
1.分析比较各种方法估计的精度；
2.分析其计算量；
3.分析噪声方差的影响；
4.比较白噪声和有色噪声对辨识的影响。

四.系统模型阶次的辨识：
1.用三种方法确定系统的阶次并辨识；
2.分析噪声对定阶的影响；
3.比较所用三种方法的优劣及有效性；
五.给出由正弦输入求取系统开环频率响应特性曲线的辨识方法。

六.提出一种自己创造的辨识新方法，并用所给数据进行辨识验证。

注：闭卷考试时提交大作业报告。

系 统 辨 识大 作 业系统辨识大作业1.考虑如下系统y() 1.5(1)0.7(2)(3)0.5(4)()k y k y k u k u k k ξ--+-=-+-+式中，()k ξ为白噪声。

取初值6(0)10P I =，ˆ(0)0θ=。

分别选择M 序列和方差为1的正态分布白噪声作为输入信号()u k ，采用递推最小二乘算法进行参数估计，迭代L=400步停止计算。

要求i) 给出基本迭代公式；ii) 画出程序流程框图；iii) 画出输入输出数据曲线、参数估计曲线、误差曲线；(提示：产生长度为L 方差为1的正态分布白噪声的MATLAB 命令为 randn （L,1）。

)答：选择如下形式的辨识模型)()4()3()2()1()(2121k k u b k u b k z a k z a k z ξ+-+-=-+-+设输入信号的取值是1．从1=k 到16=k 的M 序列，采用4位移位寄存器产生，幅度为0.03；2. 方差为1、长度为L 的正态分布白噪声，用randn(L,1)命令产生；i) 基本的迭代公式为1ˆˆˆ()(1)()[()()(1)]1()(1)()()(1)()()()[()()](1)T T T k k k z k k k k k k k k k k k I k k k θθθ-⎧=-+--⎪⎪⎡⎤=--+⎨⎢⎥Λ⎣⎦⎪⎪=--⎩K h K P h h P h P K h Pii) 程序流程框图iii)源代码与输入输出曲线、参数估计曲线、误差曲线%最小二乘递推算法辨识程序(输入信号采用M序列)clear %清理工作间变量L1=15; % M序列的周期L2=400; %迭代L2步后停止y1=1;y2=1;y3=1;y4=0; %四个移位寄存器的输出初始值%用4位移位寄存器生成M序列for i=1:L1; %开始循环，长度为L1x1=xor(y3,y4); %第一个移位寄存器的输入是第三个与第四个移位寄存器的输出的“或”x2=y1; %第二个移位寄存器的输入是第一个移位寄存器的输出x3=y2; %第三个移位寄存器的输入是第二个移位寄存器的输出x4=y3; %第四个移位寄存器的输入是第三个移位寄存器的输出y(i)=y4; %取出第四个移位寄存器的幅值为"0"和"1"的输出信号，即M序列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); %输出采样信号endfigure(2); %第2个图形stem(z),grid on %显示出输出信号径线图并给图形加上网格%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); %相对误差的初始值及大小num=0; %迭代次数初始化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; %给下次用num=num+1; %迭代次数加1if num==L2break; %如果参数收敛情况满足要求，终止计算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(3); %第3个图形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(4); %第4个图形i=1:15; %横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题输入输出曲线：参数估计曲线和误差曲线：%最小二乘递推算法辨识程序(输入信号采用正态白噪声)clear %清理工作间变量L2=400; %迭代L2步后停止%产生白噪声输入信号u=randn(L2,1); %产生长度为L2、方差为1的正态白噪声作为输入figure(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); %输出采样信号endfigure(2); %绘制输出信号图形stem(z),grid on %显示出输出信号径线图并给图形加上网格%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); %相对误差的初始值及大小num=0; %迭代次数初始化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; %给下次用num=num+1; %迭代次数加1if num==L2break; %如果参数收敛情况满足要求，终止计算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(3); %第3个图形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(4); %第4个图形i=1:15; %横坐标从1到15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') %画出a1，a2，b1，b2的各次辨识结果的收敛情况title('Identification Precision') %图形标题输入输出曲线：参数估计曲线：2.考虑如下确定性系统() 1.5(1)0.7(2)(3)0.5(4)y k y k y k u k u k --+-=-+-采用梯度校正算法进行参数估计，通过相对误差准则停止计算。

系统辨识大作业专业班级：自动化09-3学号：09051325姓名：吴恩作业一:设某物理量Y与X满足关系式2=++，实验获得一批数据Y aX bX c如下表，试辨识模型参数,,a b c。

（15分）解答：问题描述：由题意知，这是一个已知模型为Y=aX2+bX+c，给出了10组实验输入输出数据，要求对模型参数a，b，c进行辨识。

问题求解：这里对该模型参数辨识采用最小二乘法的一次算法（LS）求解。

2=++可以写成矩阵形式Y=AE+e;其中A=[X^2,X,1]构成, Y aX bX c利用matlab不难求解出结果。

运行结果：利用所求的的参数，求出给定的X对应的YE值，列表如下做出上表的图形如下12345678910xyy=ax 2+bx+c 参数求解结果分析：根据运行结果可以看出，拟合的曲线与真是观测的数据有误差，有出入，但是误差较小，可以接受。

出现误差的原因，一方面是由于给出的数据只有十个点，数据量太少，难以真正的充分的计算出其参数，另外，该问题求解采用的是LS 一次算法，因此计算方法本身也会造成相应的误差。

作业二：模仿实验二，搭建对象，由相关分析法，获得脉冲相应序列()g k，由()G z；和传递函数g k，参照讲义，获得系统的脉冲传递函数()G s及应用相关最小二乘法，拟合对象的差分方程模型；加阶跃()扰动，用最小二乘法和带遗忘因子的最小二乘法，辨识二阶差分方程的参数，比较两种方法的辨识差异；采用不少于两种定阶方法，确定对象的阶次。

对象模型如图：利用相关分析法，得到对象的脉冲相应序列。

如下图：（1）.由脉冲相应序列，求解系统的脉冲传递函数G（z）Transfer function:0.006072 z^2 + 0.288 z + 0.1671-------------------------------z^2 + 0.1018 z - 0.7509Sampling time: 2(2).由脉冲相应序列求解系统的传递函数G(s)Transfer function:(0.04849+2.494e-018i)-----------------------s^2 + 0.1315 s + 0.6048(3).利用相关最小二乘法拟合系统的差分方程模型如下：(4).在t=100，加入一个0.5的阶跃扰动，，利用RLS求解差分方程模型：RLS加入遗忘因子之后与未加之前的曲线情况如下：未加遗忘因子之前参数以及残差的计算过程加入0.99的遗忘因子得到的参数辨识过程与残差的变化过程根据上面两种方法所得到的误差曲线和参数过渡过程曲线，我们可以看出来利用最小二乘法得到的参数最终趋于稳定，为利用带遗忘因子的最小二乘算法，曲线参数最终还是有小幅度震荡。

一、 问题描述考虑仿真对象：()0.9(1)0.15(2)0.02(3)0.7(1)0.15(2)()z k z k z k z k u k u k e k +-+-+-=---+ e() 1.0e(1)0.41e(2)(),~(0,1)k k k v k v N λ+-+-=式中，u(k)和z(k)是输入输出数据，v(k)是零均值、方差为1的不相关的随机噪声；u(k)采用与e(k)不相关的随机序列。

1. 设计实验，产生输入输出数据；2. 使用基本最小二乘法估计参数；3. 考虑其他适用于有色噪声的辨识方法估计参数；4. 模型验证。

二、 问题分析对于单输入单输出系统（Single Input Single Output, SISO ），如图 1所示，将待辨识的系统看作“灰箱”，它只考虑系统的输入输出特性，而不强调系统的内部机理。

图 1中，输入u(k)和输出z(k)是可以测量的，1()G z -是系统模型，用来描述系统的输入输出特性，y(k)是系统的实际输出。

1()N z -是噪声模型，v(k)是均值为零的不相关随机噪声，e(k)是有色噪声。

图 1 SISO 系统的“灰箱”结构对于SISO 随机系统，被辨识模型()G z 为：12121212()()()1n n nn b z b z b z y z G z u z a z a z a z ------+++==++++ 其相应的差分方程为11()()()n ni i i i y k a y k i b u k i ===--+-∑∑若考虑被辨识系统或观测信息中含有噪声，被辨识模型可改写为11()()()()n ni i i i z k a y k i b u k i v k ===--+-+∑∑式中，z(k)为系统输出量的第k 次观测值；y(k)为系统输出量的第k 次真值，y(k-1)为系统输出量的第k-1次真值，以此类推；u(k)为系统的第k 个输入值，u(k-1)为系统的第k-1个输入值；v(k)为均值为0的不相关随机噪声。