系统辨识作业
系统辨识试卷A参考答案
襄樊学院2008-2009学年度上学期《系统辨识》试题A卷参考答案及评分标准一、选择题：（从下列各题的备选答案中选出一个或几个正确答案，并将其代号写在题干后面的括号内。
答案选错或未选全者，该题不得分。
每空2分，共12分）1、（C）2、（D）3、（ACD）4、（D）5、（A）6、（ABC）二、填空题：（每空2分，共14分）1、计算。
2、阶次和时滞3、极大似然法和预报误差法4、渐消记忆的最小二乘递推算法和限定记忆的最小二乘递推算法三、判断题（下列命题你认为正确的在题后括号内打“√”；错误的打“×”并改正；每小题2分，共20分）（注：正确的题目括号内打“√”得2分，打“×”得0分；错误的题目括号内打“×”得1分，改正正确再得1分，错误的题目括号内打“√”得0分；）1、（√）2、（×）参数型→非参数型3、（√）4、（×）没有→有5、（√）6、（×）考虑→基本不考虑7、（√）8、（√）9、（×）完全相同→不完全相同 10、（×）不需要→需要四、简答题：（回答要点，并简明扼要作解释，每小题6分，共18分）1、答：相关分析法的主要优点是由于M序列信号近似于白噪声，噪声功率均匀分布于整个频带，从而对系统的扰动甚微，保证系统能正常工作（1.5分）。
此外。
因为相关函数的计算是一种统计平均的方法，具有信息滤波的功能，因此，在有噪声污染下，仍可提取有用信息，准确地求出系统的脉冲响应（1.5分）。
相关辨识技术在工程中的应用、可归结为下述几个方面：（1）系统动态特性的在线测试。
包括机、炉、电等一次设备，风机、水泵等辅机以及二次自动控制系统；（1分）（2）对控制系统进行在线调试，使调节系统参数优化；（1分）（3）自适应控制中的非参数型模型辨识等。
（1分）2、答：计算中用一个数值来表示对观测数据的相对的“信任程度”，这就是权。
（2分）对于时变参数系统，其当前的观测数据最能反映被识对象当前的动态特性，数据愈“老”，它偏离当前对象特性的可能性愈大。
系统辨识大作业1201张青
《系统辨识》大作业学号：********班级：自动化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 系统模型与辨识系统模型表示了系统内部因果关系和输入输出关系，它是描述系统行为的数学方程。
系统辨识则是通过收集系统输入输出数据，根据这些数据建立模型，进而估计模型参数。
1.2 时域与频域方法在进行系统辨识时，可以采用时域方法或频域方法。
时域方法是指通过观察系统的时域响应，建立时间上的模型。
频域方法是指将系统输入输出的频谱进行分析，建立频域模型。
1.3 参数辨识与结构辨识参数辨识是指根据已知的系统输入输出数据，估计系统模型中的参数。
而结构辨识是指在已知系统输入输出数据的基础上，确定系统模型的结构或形式。
二、系统辨识方法2.1 线性系统辨识方法线性系统辨识是指对线性系统进行辨识，常用的方法包括最小二乘法、最大似然法、滑动模式控制等。
这些方法都基于线性系统的假设，且对噪声具有一定的假设条件。
2.2 非线性系统辨识方法非线性系统辨识是指对非线性系统进行辨识，因为非线性系统的行为较为复杂，因此常常需要更加复杂的模型和算法来进行辨识。
常见的方法包括神经网络、遗传算法等。
2.3 时间序列分析时间序列分析是指对系统输入输出数据在时间上的变化进行分析，用来建立系统的模型。
常用的方法包括自回归模型、移动平均模型等。
2.4 频域分析频域分析是指对系统输入输出数据的频谱进行分析，从而建立频域模型。
常用的方法包括傅里叶变换、功率谱估计等。
三、系统辨识实践练习在进行系统辨识实践练习时，首先需要明确辨识的目标和问题。
然后，收集系统的输入输出数据，并对数据进行预处理，如去噪、插值等。
接下来，选择合适的辨识方法，建立系统的数学模型，并进行参数估计。
最后，对辨识结果进行验证和评估。
系统辨识结课大作业
PDF 文件使用 "pdfFactory Pro" 试用版本创建
模型辨识
m3=pem(dat1e,‘p2zd’,‘td’,{‘max’,5}); 1+Tz*s G(s) = K * ------------------ * exp(-Td*s) (1+Tp1*s)(1+Tp2*s) Estimated using PEM from data set dat1e Loss function 0.00864915 and FPE 0.00894234
1）获取数据：
Page
3
PDF 文件使用 "pdfFactory Pro" 试用版本创建
2）数据预处理
Page
4
PDF 文件使用 "pdfFactory Pro" 试用版本创建
模型辨识 （ARX）
建立几个低阶且不同时延的ARX模型，并找出最合适 的模型阶次
结果分析
从仿真的结果来看，除了一阶无迟延模型外，其余几种辨识模型 与原吹风机的实际模型都有一定的可比性，也具有一定的跟踪能 力，在有效性校验的比较中，可以看出二阶迟延模型为88.13%， 显然二阶迟延模型的匹配度要高一点，而且二阶迟延模型的残差 均值近似为零，但仍需要进一步的优化。
Page
18
PDF 文件使用 "pdfFactory Pro" 试用版本创建
Page 13
PDF 文件使用 "pdfFactory Pro" 试用版本创建
辨识模型分析
用后300个样本值进行验证，对所取的6种模型的有效性进行比较。
二阶延迟模型的拟 合度最高
Page
系统辨识作业及答案
一． 问答题1. 介绍系统辨识的步骤。
答：（1）先验知识和建模目的的依据；（2）实验设计；（3）结构辨识；（4）参数估计；（5）模型适用性检验。
2. 考虑单输入单输出随机系统，状态空间模型[])()(11)()(11)(0201)1(k v k x k y k u k x k x +=⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡=+ 转换成ARMA 模型。
答：ARMA 模型的特点是u(k)=0,[])()(11)()(0201)1(k v k x k y k x k x +=⎥⎦⎤⎢⎣⎡=+3. 设有一个五级移位寄存器，反馈取自第2级和第3级输出的模2加法和。
试说明：（1） 其输出序列是什么？ （2） 是否是M 序列？（3） 它与反馈取自第4级与第3级输出模2加法和所得的序列有何不同？ （4） 其逆M 序列是什么？ 答：（1）设设输入序列1 1 1 1 1111018110107101006010015100114001113011112111111）（）（）（）（）（）（）（）（()()()()()()()01110161110115110101410100)13(010011210011110011110011109()()()()()()()001112401110)23(111012211010211010020010011910011180011117()()()()()()()()10011320011131011103000111291101028101002701001261001125 其输出序列为：1 1 1 1 1 0 0 1 0 1⑵不是M 序列⑶第4级与第3级模2相加结果100108001007010006100015000114001113011112111111）（）（）（）（）（）（）（）（()()()()()()()11110161110115110101410101)13(010111210110110110010110019()()()()()()()110012410010)23(001002201000211000120000111900111180111117()()()()()()()()01111321111031111013011010291010128010112710110260110025 不同点：第2级和第3级模二相加产生的序列，是从第4时刻开始，每隔7个时刻重复一次；第4级与第3级模2相加产生的,序列，是从第2时刻开始每隔15个时刻重复一次。
系统辩识作业题
系统辨识大作业
一.设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.比较所用三种方法的优劣及有效性；
五.给出由正弦输入求取系统开环频率响应特性曲线的辨识方法。
六.提出一种自己创造的辨识新方法，并用所给数据进行辨识验证。
注：闭卷考试时提交大作业报告。
系统辨识大作业
系统辨识大作业专业班级：自动化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的遗忘因子得到的参数辨识过程与残差的变化过程根据上面两种方法所得到的误差曲线和参数过渡过程曲线，我们可以看出来利用最小二乘法得到的参数最终趋于稳定，为利用带遗忘因子的最小二乘算法，曲线参数最终还是有小幅度震荡。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1、二阶反馈且存在反馈通道噪声%*******************************%闭环系统%最小二乘辨识的递推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=0.33*z(k)+0.033*z(k-1)-0.4*z(k)+w(k)%*******************************clearclc%**************产生均值为0，方差为1 的正态分布噪声************** v=random('Normal',0,1,1,402);w=random('Normal',0,1,1,402);%************产生观测序列z*******************z=zeros(402,1);u=zeros(402,1);z(1)=-1';z(2)=0;u(1)=-1';u(2)=0;for i=3:402z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=0.33*z(i)+0.033*z(i-1)-0.4*z(i-2)+w(i);end%************递推求解************P=100*eye(4); %估计方差Pstore=zeros(4,401);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,401); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=3:402h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end%***************************************************disp('参数a1、a2、b1、b2 估计结果：') Theta(:,401) i=1:401; figure(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)) title('待估参数过渡过程') gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程') Grid参数a1、a2、b1、b2 估计结果： ans =1.3842 0.4367 1.0082 0.6828500100015002000250000.511.522.533.54待估参数过渡过程a1a2b1b2二阶反馈且存在反馈通道噪声模型参数辨识过程2、二阶反馈但不加反馈噪声%*******************************%闭环系统%最小二乘辨识的递推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=0.33*z(k)+0.033*z(k-1)-0.4*z(k)%*******************************clearclc%**************产生均值为0，方差为1 的正态分布噪声************** v=random('Normal',0,1,1,2002);%w=random('Normal',0,1,1,402);%************产生观测序列z*******************z=zeros(2002,1);u=zeros(2002,1);z(1)=-1';z(2)=0;u(1)=-1';u(2)=0;for i=3:2002z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=0.33*z(i)+0.033*z(i-1)-0.4*z(i-2);end%************递推求解************P=100*eye(4); %估计方差Pstore=zeros(4,2001);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,2001); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=3:2002h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end%***************************************************disp('参数a1、a2、b1、b2 估计结果：') Theta(:,2001) i=1:2001; figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*') legend('a1','a2','b1','b2') title('待估参数过渡过程') gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程') grid参数a1、a2、b1、b2 估计结果： ans =1.3812 0.4188 0.9284 0.72315001000150020002500-0.50.511.522.533.54待估参数过渡过程a1a2b1b2二阶反馈但不加反馈噪声模型辨识过程3、一阶反馈并加反馈噪声%*******************************%闭环系统%最小二乘辨识的递推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=z(k)+0.2*z(k-1)+w(k)%*******************************clearclc%*************产生均值为0，方差为1 的正态分布噪声************* v=random('Normal',0,1,1,2002);w=random('Normal',0,1,1,2002);%************产生观测序列z*******************z=zeros(2002,1);u=zeros(2002,1);z(1)=-1';z(2)=0;u(1)=-1';u(2)=0;for i=3:2002z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=z(i)+0.2*z(i-1)+w(i);end%************递推求解************P=100*eye(4); %估计方差Pstore=zeros(4,2001);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,2001); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=3:2002h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end%***************************************************disp('参数a1、a2、b1、b2 估计结果：') Theta(:,2001) i=1:2001; figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*') legend('a1','a2','b1','b2') title('待估参数过渡过程') gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程') grid参数a1、a2、b1、b2 估计结果： ans =1.4047 0.4392 1.01270.65695001000150020002500-3-2-112345待估参数过渡过程a1a2b1b2一阶反馈并加反馈噪声模型参数辨识过程4、一阶反馈不加反馈噪声%*******************************%闭环系统%最小二乘辨识的递推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=z(k)+0.2*z(k-1)%*******************************clearclc%**************产生均值为0，方差为1 的正态分布噪声************** v=random('Normal',0,1,1,2002);%w=random('Normal',0,1,1,2002);%************产生观测序列z*******************z=zeros(2002,1);u=zeros(2002,1);z(1)=-1';z(2)=0;u(1)=-1';u(2)=0;for i=3:2002z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=z(i)+0.2*z(i-1);end%************递推求解************P=100*eye(4); %估计方差Pstore=zeros(4,2001);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,2001); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=3:2002h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end%***************************************************%*************************************************** disp('参数a1、a2、b1、b2 估计结果：') Theta(:,2001) i=1:2001; figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*') legend('a1','a2','b1','b2') title('待估参数过渡过程') gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程') grid参数a1、a2、b1、b2 估计结果： ans =8.4624 1.7891 8.05780.63455001000150020002500012345678910待估参数过渡过程a1a2b1b2一阶反馈并加反馈噪声模型参数辨识过程50010001500200025000102030405060708090100估计方差变化过程5、比例反馈加反馈噪声%******************************* %闭环系统%最小二乘辨识的递推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k) %u(k)=1.2*z(k)+w(k)%******************************* clear clc%===========产生均值为0，方差为1 的正态分布噪声========= v=random('Normal',0,1,1,2002); w=random('Normal',0,1,1,2002);%************产生观测序列z******************* z=zeros(2002,1); u=zeros(2002,1); z(1)=-1'; z(2)=0; u(1)=-1'; u(2)=0;for i=3:2002z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=1.2*z(i)+w(i);end%************递推求解************P=100*eye(4); %估计方差Pstore=zeros(4,2001);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,2001); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=3:2002h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end%*************************************************** %*************************************************** disp('参数a1、a2、b1、b2 估计结果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*') legend('a1','a2','b1','b2')title('待估参数过渡过程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('估计方差变化过程')grid参数a1、a2、b1、b2 估计结果：ans =1.45080.47301.04140.71355001000150020002500-5-4-3-2-11234待估参数过渡过程a1a2b1b2比例反馈加反馈噪声模型参数辨识过程6、比例反馈不加反馈噪声%******************************* %闭环系统%最小二乘辨识的递推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k) %u(k)=1.2*z(k)%******************************* clear clc%**************产生均值为0，方差为1 的正态分布噪声************** v=random('Normal',0,1,1,2002); w=random('Normal',0,1,1,2002);%************产生观测序列z******************* z=zeros(2002,1); u=zeros(2002,1); z(1)=-1'; z(2)=0; u(1)=-1'; u(2)=0;for i=3:2002z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=1.2*z(i);end%************递推求解************P=100*eye(4); %估计方差Pstore=zeros(4,2001);Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)];Theta=zeros(4,2001); %参数的估计值，存放中间过程估值Theta(:,1)=[3;3;3;3];% K=zeros(4,400); %增益矩阵K=[10;10;10;10];for i=3:2002h=[-z(i-1);-z(i-2);u(i-1);u(i-2)];K=P*h*inv(h'*P*h+1);Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));P=(eye(4)-K*h')*P;Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4)];end%*************************************************** disp('参数a1、a2、b1、b2 估计结果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*') legend('a1','a2','b1','b2')title('待估参数过渡过程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:))title('估计方差变化过程')grid参数a1、a2、b1、b2 估计结果：ans =3.33093.54102.60293.271350010001500200025000.511.522.533.544.5待估参数过渡过程a1a2b1b2比例反馈不加反馈噪声模型参数辨识过程5001000150020002500102030405060708090100估计方差变化过程闭环辨识结果反馈律反馈噪声参数估计值可辨识性a1(1.4) a2(0.45) b1(1.0) b2(0.7)2阶反馈：u(k)=0.33*z(k)+0.033*z (k-1)-0.4*z(k) w(k)≠0 1.3842 0.4367 1.0082 0.6828 唯一，一致w(k)=0 1.3812 0.4188 0.9284 0.7231 唯一，一致1阶反馈：u(k)=z(k)+0.2*z(k-1) w(k)≠0 1.4047 0.4392 1.01270.6569唯一，一致w(k)=0 8.4624 1.7891 8.0578 0.6345 不可辨识0阶反馈：u(k)=1.2*z(k) w(k)≠0 1.4508 0.4730 1.0414 0.7135 唯一，一致w(k)=0 3.3309 3.5410 2.6029 3.2713 不可辨识实验结果分析：1）当反馈通道为二阶方程时，无论反馈通道是否有噪声输入，闭环系统前向通道模型参数估计值有唯一解，且具有一致性。