系统辨识最小二乘法大作业 (2)

相关分析法辨识系统单位脉冲响应1辨识原理对于下图示的单输入单输出线性系统，其输入输出的因果关系可用卷积公式描述。

公式为：0()()()y t g x t d λλλ∞=-⎰把变量t 换成t +τ，上式两边同乘以x (t )，取时间的平均值，得11lim()(+)(){lim()(+)}22TTTTT T x t y t dt g x t x t dt d TTτλτλλ∞--→∞→∞=-⎰⎰⎰即 0()()()x y x R g R d τστλλ∞=-⎰上式即为维纳-霍夫方程，其给出了输入的自相关函数，输入、输出的互相关函数及脉冲响应函数三者之间的关系。

令x (t )为白噪声信号，则其自相关函数为：()(), ()()x x R k R k τδττλδτλ=-=-代入维纳-霍夫方程得：()()()()xy x R g R d kg τλτλλτ∞=-=⎰则有：()()xy R g kττ=这样，只要记录x(t)、y(t)的值，并计算它们的互相关函数，即可求得脉冲响应函数g(τ)。

在系统有正常输入的情形下，辨识脉冲响应的原理图如下图所示。

2辨识过程2.1预备实验以二阶系统22()2G s s s ++=作为辨识对象。

在实验前首先要进行预备实验，以了解系统特性。

通过简单阶跃响应确定系统过度过程时间T s 大约为11s ，如下图所示。

给系统施加不同周期的正弦信号，系统输出为输入的0.707倍时，确定截止频率f M 大约为0.318Hz 。

2.2选择二位式伪随机序列的参数（1）选择t 和N 由0.3Mt f ∆≤，得0.94t s ∆≤。

因为系统的时间常数1nT s ζω=，根据时间常数可按照0.050.1t T ∆= （）选择t ∆。

由二位式伪随机序列周期要大于系统过渡过程时间，若t ∆选择0.94s ，则由(1)s N t T -⨯∆≥，得12.7021N ≥；若t ∆选择0.195s ，则由(1)s N t T -⨯∆≥，得57.4103N ≥。

1、利用cftool拟合曲线得Coefficients (with 95% confidence bounds):a = 0.4601 (0.3711, 0.5492)b = -5.1 (-6.108, -4.092)c = 13.42 (11, 15.83)2、最小二乘法：最小二乘的思想就是寻找一个θ 的估计值θ ，使得各次测量的Z i (i = 1,⋯ m)与由估计θ 确定的量测估计Z i = H iθ 之差的平方和最小，由于此方法兼顾了所有方程的近似程度，使整体误差达到最小，因而对抑制误差是有利的。

但是最小二乘一般方法的估计精度不够高，这是由于对各个测量数据同等对待，而各次测量数据一般不会在相同的条件下获得，造成测量数据的置信度不变较大，当新数据源源而来时，将出现以下问题：1.数据增加，要求计算机的存储空间增加2、每增加一组数据，即作一次求逆，导致计算量增加，难以用于在线辨识。

只有当模型的噪声项是独立的随机变量时，普通最小二乘法才能得到真实参数的无偏估计值，否则所得到的估计值是有偏的。

因此，最小二乘法有很多改进算法，虽然没有一个是完美的，但是能够适应不同的情况、条件，对应选择不同的算法，其各自的性能及优缺点如下：广义最小二乘法的优点是计算精度高，估计的效果比较好，是无偏估计，但广义最小二乘法缺点是计算量大，其收敛是比较缓慢的，为了得到准确的参数估值，往往需要进行多次迭代计算，另外，对于循环的循环性没有给出证明，并非总是收敛于最优估值上。

由于一般情况下，系统信噪比比较低，准则函数为非单值函数（即存在多个局部极小值），如果初值给的不合理，用GLS方法得到的将是局部极小值，若想得到总体最优解，初值应接近该最优值。

递推算法实现了实时控制，减少了计算量和存储量，但未解决最小二乘法的递推算法有偏估计问题。

矩形窗（限定记忆）RLS 方法需要保留一定的数据存储量，此存储量大小取决于矩形窗宽度，因而在应用范围上有一定程度的限制。

学号：XXXX大学系统辨识实验报告最小二乘辨识的应用院（系）计算机与信息工程学院专业控制理论与控制工程学生姓名XXXXX成绩指导教师2013年6月摘要：系统辨识在工程中的应用非常广泛，系统辨识的方法有很多种，最小二乘法是一种应用极其广泛的系统辨识方法。

本文阐述了动态系统模型的建立及其最小二乘法在系统辨识中的应用，并通过实例分析说明了最小二乘法应用于系统辨识，比如模糊最小二乘辨识应用于等精度测量，最小二乘法辨识算法在故障检测中的应用，最小二乘辨识模型的控制器误差修正。

关键词：最小二乘法；系统辨识；辨识精度随着科学技术的不断发展，人们认识自然、利用自然的能力越来越强，对于未知对象的探索也越来越深入。

我们所研究的对象，可以依据对其了解的程度分为三种类型：白箱、灰箱和黑箱。

如果我们对于研究对象的内部结构、内部机制了解很深入的话，这样的研究对象通常称之为“白箱”；而有的研究对象，我们对于其内部结构、机制只了解一部分，对于其内部运行规律并不十分清楚，这样的研究对象通常称之为“灰箱”；如果我们对于研究对象的内部结构、内部机制及运行规律均一无所知的话，则把这样的研究对象称之为“黑箱”。

研究灰箱和黑箱时，将研究的对象看作是一个系统，通过建立该系统的模型，对模型参数进行辨识来确定该系统的运行规律。

对于动态系统辨识的方法有很多，但其中应用最广泛，辨识效果良好的就是最小二乘辨识方法，研究最小二乘法在系统辨识中的应用具有现实的、广泛的意义。

1 动态系统模型的建立要想了解动态系统的特性，首先就要根据先验知识和对于该系统了解的程度建立该动态系统的模型。

动态系统的建模问题是进行系统辨识首先要解决的问题。

建立被控系统的数学模型通常有理论分析和试验分析两种方法。

理论分析方法就是在已知系统内部规律的基础上(此时该系统即为“白箱”)，推导出系统的动态方程和输入输出关系。

然而人们对一些复杂被控对象(“灰箱”或“黑箱”)的内部规律及其诸多重要参数认识得并不十分清楚。

《系统辨识》实验二要点实验二 递推最小二乘估计(RLS)及模型阶次辨识(F-Test ）一、实验目的① 通过实验，掌握递推最小二乘参数辨识方法 ② 通过实验，掌握F-Test 模型阶次辨识方法二、实验内容1、仿真模型实验所用的仿真模型如下： 框图表示模型表示)()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z λ+-+-=-+-- 其中u (k )和z (k )分别为模型的输入和输出变量；v (k )为零均值、方差为1、服从正态分布的白噪声；λ为噪声的标准差（实验时，可取0.0、0.1、0.5、1.0）；输入变量u (k )采用M 序列，其特征多项式取1)(4⊕⊕=s s s F ,幅度取1.0。

2、辨识模型辨识模型的形式取)()()()()(11k e k u z B k z z A +=--为方便起见，取n n n b a ==，即nn nn zb z b z b z B z a z a z a z A ------+++=++++= 22112211)(1)(根据仿真模型生成的数据{}L k k u ,,1),( =和{}L k k z ,,1),( =，辨识模型的参数n n b b b a a a ,,,,,,2121 和；并确定模型阶次n ，同时估计出模型误差)(k e 的方差（应近似等于模型噪声)(k v 的方差，即为2λ）和模型的静态增益K 。

3、辨识算法① 采用递推遗忘因子法：[][][]⎪⎪⎩⎪⎪⎨⎧--=+--=--+-=-)1()()(1)()()1()()()1()()1()()()()1()(1k k k μk μk k k k k k k k k z k k k P h K I P h P h h P K h K τττθθθ 其中，遗忘因子10≤<μ(具体值根据情况自已确定)；数据长度L 可取100、300、500；初始值⎩⎨⎧==IP 2)0()0(a εθ。

作业四宋家亮15030024一、实验题目用MatLAB辨识系统，系统输出分别叠加两种不同类型的噪声，针对每一种情况使用最小二乘整批算法、递推算法和广义最小二乘法实现辨识。

图1 图2二、实验思路及目的对图1、图2分别通过构造的系统给定输入（白噪声和M序列）并叠加噪声（白噪声和有色噪声），测得输出数据并显示，对未知系统定阶并用整批、递推、广义最小二乘算法进行辨识，显示辨识结果；比较三种算法对不同噪声模型的辨识精度（λ取相同的值），显示辨识结果的脉冲响应图像并于理想系统响应对比，对结果给予合理的解释；改变 值(0.01,0.02,0.03等)，比较辨识结果的精度，说明信噪比对各种算法辨识精度的影响。

三、实验过程1、设计界面利用guide编辑器针对所要实现的功能添加控件进行界面初步的设计，设计结果如图1.1所示图1.12、设置控件属性对界面的控件设置属性，形成实验所需要的最终界面，如图2.1所示：图2.1其中系统噪声和输入信号设置成可选择的属性，其value值作为后面函数中的控制条件，这个刚开始走了弯路，直接用string属性做控制条件，结果发现根本实现不了，请教大神后选择value值作为控制条件，因为value值默认为1，当通过下拉选项选择下一个信号时其value值会自动加1。

如图2.2所示。

图2.23、设计guide回调函数（1）噪声参数输入的回调函数本实验设计成噪声参数可输入的程序，所以为达到这一目的，添加了参数输入控件，通过如图3.1所示打开回调函数，添加如下程序，可以将输入的容转变成数字用于后续函数的运算。

global zaoshengxishuzaoshengxishu=get(handles.zaoshengxishu,'String'); zaoshengxishu=str2num(zaoshengxishu);%将输入的字符串λ变成数值图3.1（2）开始辨识按钮的回调函数（主体程序）如图3.2打开开始辨识按钮回调函数，添加主程序（见附录），这样通过点击该按钮可实现本实验的所有功能。

2012系统辨识实验姓名：周自飞学号：P4*******年级专业：09自动化线性系统参数估计的最小二乘法一、实验目的：1、掌握线性离散系统的数学模型；2、掌握线性离散系统在无噪声时的单位脉冲响应；3、掌握线性系统参数估计的最小二乘法。

二、实验原理：1、基本最小二乘算法2、递推最小二乘算法3、渐消记忆递推最小二乘算法4、辅助变量递推最小二乘算法5、增广递推最小二乘算法三、实验内容：设带有噪声的离散系统模型为：()()()()()()k213.0-13.0-+.028.11-412+ukkukyk--yξky=现给定系统的一列长度为50的输入序列：u(k)= 1,0.3,-0.5,0.9,-0.5,-0.3,-0.2,0.4,0.8,-0.6,-0.1,0,0.1,0.5,-0.6,-0.2,0.3,0.9,0.5,0.2,-0.6,-0.3,-0.1,0.7,0,0.3,-0.6,0.3,0.1,0.5,-0.7,-0.4,-0.9,-0.6,-0.2,-0.4,-0.2,0.1,-0.1,0.1,0.9,0.5,0.3,0.7,0.4,-0.2,-0.7,-0.2,0.1,01. 构造离散系统模型函数，给出无噪声时系统在单位脉冲输入下的响应序列；%zhouzifei.mb=[0.3 -0.213];a=[1 -1.28 0.41];impz(b,a,30);title(‘系统单位脉冲响应’);axis([-1 30 -0.2 0.5]);2. 给出系统在给定输入u(k)和噪声()kξ下的输出信号y(k)；%zifei2.m3. 利用已有输入信号u(k)和在第2步中得出的输出信号y(k)，使用下列方法辨识系统的参数：(1) 基本最小二乘算法%LS.m参数：a1 =-1.2803a2 =0.4138b0 = 0.2867b1 = -0.2104(2) 递推最小二乘算法%RLS.m参数：a1 =-1.2800 a2 =0.4100 b0 = 0.3000 b1 = -0.2130(3) 渐消记忆递推最小二乘算法%FMRLS参数：a1 =-1.2800 a2 =0.4100 b0 = 0.3000 b1 = -0.2130(4) 辅助变量递推最小二乘算法%IVLS.m参数：a1 =-1.2800 a2 =0.4100b0 = 0.3000 b1 = -0.2130(5) 增广递推最小二乘算法（可选做，有加分）%ARLS.m参数：a1 =-1.2800 a2 =0.4100 b0 = 0.3000 b1 = -0.2130四、实验结果：最小二乘法思想是使各次观测值和计算值之间差值的平方乘以度量其精确度的数值后的和为最小。

系统辨识与参数估计大作业第一题 递推最小二乘估计参数考虑如上图所示的仿真对象，选择模型结构为：)()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+，其中)(k v 是服从)1,0(N 正态分布的不相关随机噪声；输入信号)(k u 采用4阶逆M 序列，特征多项式取41)(s s s F ⊕⊕=，幅度为1，循环周期为bit N p 62=；控制λ值，使数据的噪信比分别为10％，73％，100％三种情况。

加权因子1)(=Λk ；数据长度L=500;初始条件取001.0)0(ˆ,10)0(6==θI P ， （1） 利用递推最小二乘算法在线估计参数，（2） 利用模型阶次辨识方法（AIC 准则），确定模型的阶次。

（3） 估计噪声)(k v 的方差和模型静态增益K （4） 作出参数估计值随时间的变化图 答：设过程的输入输出关系可以描述成()()()Tz k h k n k θ=+()z k 是输出量，()h k 是可观测的数据向量，n (k)是均值为0的随机噪声[]()(1),(2),(1),(2)Th k z k z k u k u k =------[]1212,,,Ta ab b θ=选取的模型为结构是1212()(1)(2)(1)(2)z k a z k a z k bu k bu k =----+-+- 12121.5,0.7, 1.0,0.5a a b b =-===加权最小二乘参数估计递推算法RWLS 的公式如下，11()(1)()()(1)()()()(1)()()()(1)()()()(1)T T TK k p k h k h k p k h k k k k K k z k h k k p k I K k h k p k θθθ-⎡⎤=--+⎢⎥Λ⎣⎦⎡⎤=-+--⎣⎦⎡⎤=--⎣⎦为了把p(k)的对称性，可以把p(k)写成1()(1)()()()(1)()()T T p k p k K k K k h k p k h k k ⎡⎤=---+⎢⎥Λ⎣⎦如果把()k Λ设成1的时候，加权最小二乘法就退化成最小二乘法。

系统辨识一、最小二乘算法仿真对象为：)()2(5.0)1()2(7.0)1(5.1)(k v k k u k k z k z +-+-=-+-+，)(k v 为服从N(0,1)分布的白噪声，输入信号采用4阶M 序列。

辨识模型为：)()2()1()2()1()(2121k v k z b k z b k z a k z a k z +-+-+----=，观测数据长度取400=L ，加权矩阵取I =Λ。

则matlab 程序如下:% M array constructionNp=15;r=4;X1=1;X2=1;X3=1;X4=1;m_length = r*Np;a=1;for i=1:1:m_lengthY4=X4;Y3=X3;Y2=X2;Y1=X1;X4=Y3;X3=Y2;X2=Y1;X1=xor(Y3,Y4);if Y4==0M(i)=-a;elseM(i)=a;endendfigure;i=1:1:m_length;plot(i,M);% 白噪声noise = zeros(1,m_length);for i=1:1:m_lengthtemp = noise + 0.5*rands(1,m_length);noise = temp;endnoise = noise/12;%noise = temp;% parameter of systemn=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;S_U0=0.2;S_Y0=0.2;% generate u,yu0=ones(1,m_length)*S_U0;U = M + u0 + noise;figure;i=1:1:m_length;plot(i,U);%formulationy(1)=0;y(2)=0;y(3)=0;Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3); for k=4:m_lengthy(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2);Y(k)=S_Y0+y(k)+noise(k);endfigure;i=1:1:m_length;plot(i,Y);%辨识% premitive valuec=2;P = (c^3)*eye(3*n+d);sita(:,3) = [0,0,0,0,0,0,0]';alf = 0.995;% M%compute U0,Y0sum_U = 0;sum_Y = 0;for k=1:1:m_lengthsum_U = sum_U + U(k);sum_Y = sum_Y + Y(k);endt_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length;for k=1:1:m_lengtht_u(k) = U(k) - t_U0;t_y(k) = Y(k) - t_Y0;endfigure;i=1:1:m_length;plot(i,t_u);figure;i=1:1:m_length;plot(i,t_y);v(1)=0;v(2)=0;v(3)=0;for k=4:1:m_lengthfai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d),v(k-1),v(k-2),v(k-3)]';v(k) = t_y(k) - fai'*sita(:,k-1);W = P*fai;dan = 1/(alf + fai'*W);sita(:,k) = sita(:,k-1) + dan*W*v(k);P = ( P - dan * W*W')/alf;end仿真结果如下图所示：则各个参数的估计值为：⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛∧∧∧∧5048.00062.17003.04961.1121bbaa二、梯度校正辨识方法仿真对象为)()2(5.0)1()2(7.0)1(5.1)(kvkkukkzkz+-+-=-+--，其中)(kv为服从N(0,1)分布的白噪声，取初值0)0(=∧θ，1=α，1.0=c。

系统辨识大作业专业班级：自动化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的不相关随机噪声。

西北工业大学系统辩识大作业题目：最小二乘法系统辨识一、 问题重述：用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数离散化有z^4 - 3.935 z^3 + 5.806 z^2 - 3.807 z + 0.9362---------------------------------------------- =z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187噪声的成形滤波器离散化有4.004e-010 z^3 + 4.232e-009 z^2 + 4.066e-009 z + 3.551e-010-----------------------------------------------------------------------------=z^4 - 3.808 z^3 + 5.434 z^2 - 3.445 z + 0.8187采样时间0.01s要求：1.用Matlab 写出程序代码；2.画出实际模型和辨识得到模型的误差曲线；3.画出递推算法迭代时各辨识参数的变化曲线；最小二乘法:在系统辨识领域中 ,最小二乘法是一种得到广泛应用的估计方法 ,可用于动态 ,静态 , 线性 ,非线性系统。

在使用最小二乘法进行参数估计时 ,为了实现实时控制 ,必须优化成参数递推算法 ,即最小二乘递推算法。

这种辨识方法主要用于在线辨识。

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 bni na i b a +-=-+∑∑==其中Z （k ）为输出变量，u(k)为输入变量，e(k)为偏差。

上式可以表示为)()()(-)k (1i 1i k e i k u i k Z Z bni 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 为所需要采集的点数。

[]Tn 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 。

五、实验代码以及实验结果m=20; %置M序列总长度y1=1;y2=1;y3=1;y4=0;for i=1:mx1=xor(y3,y4);%异或运算x2=y1;x3=y2;x4=y3;if y4==0;u(i)=1;elseu(i)=-1;endy1=x1;y2=x2;y3=x3;y4=x4;endz=zeros(21,1);%定义输出观测值的长度21行*1列的0矩阵ZL=zeros(19,1);%定义输出观测值的长度19行*1列的0矩阵for k=3:21z(k)=-1.5*z(k-1)-z(k-2)+u(k-1)+3*u(k-2) ;%用理想输出值作为观测值ZL(k-2)=z(k);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 %画出输出观测值z的图形，并显示坐标网格for m=2:20HL=[-z(m) -z(m-1) u(m) u(m-1)];for n=1:4;H((m-1),n)=HL(n);endendc1=H'*H;c2=inv(c1);c3=H'*ZL;c=c2*c3a1=c(1), a2=c(2), b1=c(3), b2=c(4)实验结果：c =1.50001.00001.00003.0000a1 =1.5000a2 =1.0000b1 =1b2 =3.0000六、实验结果分析通过实验结果可知所得的实验结果与待辨识的系统传递函数的系数很接近了。