系统辨识最小二乘法大作业 (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的时候，加权最小二乘法就退化成最小二乘法。
第七章—最小二乘参数辨识(II)
h ( k ) = [ − z ( k − 1),
aБайду номын сангаас
1、一次完成算法 ⎧ A( z
, − z ( k − na ), u ( k − 1),
b
A( z −1 ) z (k ) = B( z −1 )u (k ) + v(k )
−1
, u ( k − nb )]
T
θ = [a1 , , an , b1 , , bn ]T
限定记忆法的递推算法（RFM）
ˆ ˆ θ (k + 1, k + L) = θ (k , k + L) − K (k + 1, k + L)[ z (k ) − hT (k )θˆ(k , k + L)] K (k + 1, k + L) = P(k , k + L)h(k )[1 − hT (k ) P(k , k + L)h(k )]−1 P(k + 1, k + L) = ⎡ I + K (k + 1, k + L)hT (k ) ⎤ P(k , k + L) ⎣ ⎦ ˆ ˆ ˆ θ (k , k + L) = θ (k , k + L − 1) + K (k , k + L)[ z (k + L) − hT (k + L)θ (k , k + L − 1)] K (k , k + L) = P(k , k + L − 1)h(k + L)[1 + hT (k + L) P(k , k + L − 1)h(k + L)]−1 P(k , k + L) = ⎡ I − K (k , k + L)hT (k + L) ⎤ P(k , k + L − 1) ⎣ ⎦
第五章 最小二乘法辨识2
e(2) (k)
，按步骤
（3）重新估计f，得到估计值
^
f
(2)
。再按步骤（4）
计算
(2)
y (k)
和
u
(
2
)
(k
)
，按步骤（5）求
的第3次估计
值。
重复上述循环步骤，直到
的估计值
^
(i)
收敛为止。
上述循环的收敛性可用下式判断，即
^ (i)
lim f (z 1 ) 1
i
即当
i 比较大时，如果
5、辅助变量法
对于原辨识问题
y
（1）
当 (k) 是不相关随机序列时，最小二乘法可得到 参数向量 的一致性无偏估计。
但实际应用中，(k) 往往是相关随机序列。
假定存在一个 (2n 1) N 的矩阵 （Z 与 同阶数）
满足约束条件
Nlim
1 N
ZT
E
ZT
0
Nlim
1 N
ZT
估计值
^ (1)
f
((1) )T (1)
1 ((1) )T e(1)
式中
^ (1)
f1
^ (1)
^ (1)
f f2
^
(1)
f m
e(1) (n 1)
e (1)
e
(1)
(n
2)
e
(1)
(n
N
)
e(1) (n)
e(1) (n 1)
(1)
e(1) (n 1)
)
是持续激励信号时，必有
E
^
T k
是非奇异矩
阵。又因为
^
系统辨识作业_梯度法最小二乘法
系统辨识一、最小二乘算法仿真对象为：)()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。
系统辨识各类最小二乘法汇总
yk(k)=1.5*yk(k-1)-0.7*yk(k-2)+uk(k-1)+0.5*uk(k-2)+y1(k); end figure(3); plot(yk); title('对应输出曲线');
theta=[0;0;0;0]; p=10^6*eye(4);
9
for t=3:N h=([-yk(t-1);-yk(t-2);uk(t-1);uk(t-2)]); x=1+h'*p*h; p=(p-p*h*1/x*h'*p); theta=theta+p*h*(yk(t)-h'*theta);
12
p=(p-p*h*1/x*h'*p); theta=theta+p*h*(yk(t)-h'*theta);
a1t(t)=theta(1); a2t(t)=theta(2); b1t(t)=theta(3); b2t(t)=theta(4); d1t(t)=theta(5); d2t(t)=theta(6);
end 5、RGLS 试验程序（部分） for t=3:N
he=([-e(t-1);-e(t-2)]); xe=1+he'*pe*he; pe=(pe-pe*he*1/xe*he'*pe); thete=thete+pe*he*(e(t)-he'*thete);
c1t(t)=thete(1); c2t(t)=thete(2);
7
RELS： 当噪声模型： e k = D Z −1 ∗ v k ( v(k) 为白噪声 )时，我们采用增广最 小二乘方法。能辨识出参数（包括噪声参数）的无偏估计。 RGLS： 当噪声模型： e 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的遗忘因子得到的参数辨识过程与残差的变化过程根据上面两种方法所得到的误差曲线和参数过渡过程曲线，我们可以看出来利用最小二乘法得到的参数最终趋于稳定，为利用带遗忘因子的最小二乘算法，曲线参数最终还是有小幅度震荡。
第4章最小二乘类辨识算法 (2)
另外
设模型的噪声 n(k) 特征为
E{n(1)} E{n(2)} 0 E{n L } E{n( L)}
E{n 2 (1)} T E{n(2)n(1)} cov{nL } E{nL nL } E{n( L)n(1)} E{n(1)n(2)} E{n 2 (2)} E{n( L)n(2)}
对 k 1,2,, L
(4.1.5)式构成一个线性方程组 可以写成
z L H L nL
（4.2.4）
22
z L [ z (1), z (2), z ( L)]
hT (1) z (0) T h (2) z (1) HL T h ( L) z( L 1)
最后，假设数据长度 L (na nb )
k, l
（4.2.8）
25
(4.2.4)式有L个方程，包括 na nb 个未知数。 如果 L na nb ，方程的个数少于未知数的 个数，模型参数 不是唯一确定。 如果 L na nb ，则只有当 n 0 时， L 才有唯一确定解。 当 nL 0 时，只有取 L na nb ，才有可 能确定一个最优的模型参数 ，而且为了 保证辨识的精度，L必须充分大。
（4.0.5）
10
各种方法所用的辨识模型结构略有不同
最小二乘法（受控自回归 CAR模型）
A( z 1 ) z (k ) B( z 1 )u (k ) v(k )
（4.0.7）
增广最小二乘法(受控自回归滑动平均 CARMA模型)
A( z 1 ) z (k ) B( z 1 )u (k ) D( z 1 )v(k )
系统辨识大作业.
一、 问题描述考虑仿真对象：()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的不相关随机噪声。
(完整word版)2003版系统辨识最小二乘法大作业
西北工业大学系统辩识大作业题目：最小二乘法系统辨识一、 问题重述：用递推最小二乘法、加权最小二乘法、遗忘因子法、增广最小二乘法、广义最小二乘法、辅助变量法辨识如下模型的参数离散化有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 是一套高性能数字计算和可视化软件 ,它集成概念设计 ,算法开发 ,建模仿真 ,实时实现于一体 ,构成了一个使用方便、界面友好的用户环境 ,其强大的扩展功能为各领域的应用提供了基础。
对于一个简单的系统 ,可以通过分析其过程的运动规律 ,应用一些已知的定理和原理,建立数学模型 ,即所谓的“白箱建模 ”。
但对于比较复杂的生产过程 ,该建模方法有很大的局限性。
由于过程的输入输出信号一般总是可以测量的 ,而且过程的动态特性必然表现在这些输入输出数据中 ,那么就可以利用输入输出数据所提供的信息来建立过程的数学模型。
系统辨识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六、实验结果分析通过实验结果可知所得的实验结果与待辨识的系统传递函数的系数很接近了。
MIMO系统最小二乘法参数辨识
转化为最小二乘形式
转换为标准最小二乘写法
式中，
转换最小二乘形式
2. 根据最小二乘法，求出对应的第i行对应的元素 3. 依次对每个输出yi，重复1-2 4. 将结果组合，得到相应的A和R 5. 由A和R，反求原来形式的矩阵，得到解。
最小二乘法只能处理如下的形式 其中，Y为Nx1的列向量，H为矩阵，r为参数列向量
最小二乘法特点
因为最小二乘发只能处理写成最小二乘表达式，
因此，最小二乘法很适合处理SISO的辨识。过程如下： 1. 将SISO系统写成差分方程表达形式，即连续传递函数离散化
最小二乘法
最小二乘法只能处理最小二乘形式，对系统辨识而言，只能直接处 理如下格式
s^2 + 2.545 s + 3.24 为了计算机便于计算，需要离散化后，得到离散格式的传递函数 dsys，对得到的dsys使用扫频信号u激励，得到输出y。根据得到 的u和y，辨识得到离散的传递函数dsysL，然后得到连续传递函数 sysL。
过程数据
辨识得到的结果
响应曲线验证
最小二乘法特点
最小二乘法特点
2. 上式中生成的a_i，b_i当做未知参数theta的分量 3. 写成最小二乘表达法
4. 求出a_i,b_i,得到离散传递函数（1步） 5. 由离散传递函数求得连续传递函数。
MIMO特点
状态空间表达式一般为
要想用最小二乘法辨识，有如下两个问题 1. 状态变量X需要消去
MIMO 最小二乘辨识
李曹磊
辨识目的
某些小型飞机不知道转动惯量、气动导数，又无法由传统方法获得 （或经费所限），可以通过辨识来求。从而为飞控系统设计提供自 然飞机模型。 对于某些飞机，可以验证传统方法得到的气动导数。
第4章最小二乘类辨识算法2
k 1
L
(k )[ z (k ) h T (k ) ]2
(k ) - 加权因子，对 k , (k ) 0
如
(k ) Lk
L 1
0 1
K = 1 时 (1) K = L 时
1
19
( L) 1
准则函数 J ( ) 可写成二次型形式
ˆ 所以 WLS 是无偏估计。
35
无偏性并不要求噪声 nL 一定是白噪声， HL 只要求它与 统计独立即可。如果 nL HL nL 是白噪声，则 与 一定统计独立。 另外，定理１所给出是条件是 ˆ θ WLS 为无偏估计的充分条件，并不是必要条 件。它的必要条件应是
1 T E{(HT Λ H ) H L Λ L nL } 0 L L L
J ( ) ( z L H L )T L ( z L H L )
L - 加权矩，一般为正定的对角矩阵
(1) 0 L 0 (2) 0 0 0 0 ( L)
20
设
使
WLS
J ( )
WLS
min
则有
6
各种方法所用的辨识模型结构略有不同
最小二乘法（受控自回归 CAR模型）
A( z 1 ) z(k ) B( z 1 )u(k ) v(k )
增广最小二乘法(受控自回归滑动平均 CARMA模型)
A( z 1 ) z(k ) B( z 1 )u(k ) D( z 1 )v(k )
4 阶 M 序 列
输 出 信 号
解：由于输入信号为 4 阶 M 序列，所以 M 序列的循 环长度为 L 2 4 1 15 。因此设输入信号的取值从 k 1 到 k 16 的 M 序列，于是可得
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院学号：姓名：日期：基于最小二乘法的多种系统辨识方法研究一、实验原理1.最小二乘法在系统辨识中用得最广泛的估计方法是最小二乘法(LS)。
设单输入-单输出线性定长系统的差分方程为(5.1.1)式中:为随机干扰；为理论上的输出值。
只有通过观测才能得到，在观测过程中往往附加有随机干扰。
的观测值可表示为(5.1.2)式中:为随机干扰。
由式(5.1.2)得(5.1.3)将式(5.1.3)带入式(5.1.1)得(5.1.4)我们可能不知道的统计特性，在这种情况下，往往把看做均值为0的白噪声。
设(5.1.5)则式(5.1.4)可写成(5.1.6)在观测时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。
因此假定不仅包含了的测量误差，而且包含了的测量误差和系统内部噪声。
假定是不相关随机序列(实际上是相关随机序列)。
现分别测出个随机输入值，则可写成个方程，即上述个方程可写成向量-矩阵形式(5.1.7) 设则式(5.1.7)可写为(5.1.8)式中：为维输出向量；为维噪声向量；为维参数向量；为测量矩阵。
因此式(5.1.8)是一个含有个未知参数，由个方程组成的联立方程组。
如果，方程数少于未知数数目，则方程组的解是不定的，不能唯一地确定参数向量。
如果，方程组正好与未知数数目相等，当噪声时，就能准确地解出(5.1.9)如果噪声，则(5.1.10)从上式可以看出噪声对参数估计是有影响的，为了尽量较小噪声对估值的影响。
在给定输出向量和测量矩阵的条件下求系统参数的估值，这就是系统辨识问题。
可用最小二乘法来求的估值，以下讨论最小二乘法估计。
2.最小二乘法估计算法设表示的最优估值，表示的最优估值，则有(5.1.11)写出式(5.1.11)的某一行，则有(5.1.12) 设表示与之差，即-(5.1.13)式中成为残差。
把分别代入式(5.1.13)可得残差。
设则有(5.1.14) 最小二乘估计要求残差的平方和为最小，即按照指数函数(5.1.15) 为最小来确定估值。
求对的偏导数并令其等于0可得(5.1.16)(5.1.17)由式(5.1.17)可得的最小二乘估计(5.1.18) 3.递推最小二乘法为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。
设已获得的观测数据长度为，将式(5.1.8)中的和分别用来代替，即(5.3.1) 用的最小二乘估计，则(5.3.2) 设(5.3.5)于是(5.3.6)如果再获得1组新的观测值和，则又增加1个方程(5.3.7) 式中将式(5.3.1)和式(5.3.7)合并，并写成分块矩阵形式，可得(5.3.8) 根据上式可得到新的参数估值(5.3.9) 式中根据矩阵求逆引理可以求得递推最小二乘法辨识公式(5.3.19)(5.3.20)(5.3.21) 由于进行递推计算需要给出和初值和，通过计算证明，可以取初值：，，c是充分大的常数，为单位矩阵，则经过若干次递推之后能够得到较好的参数估计。
3.辅助变量法辅助变量法是一种可克服最小二乘有偏估计的一种方法，对于原辨识方程(5.4.1)当是不相关随机序列时，最小二乘法可以得到参数向量的一致无偏估计。
但是，在实际应用中往往是相关随机序列。
假定存在着一个的矩阵满足约束条件(5.4.2)式中是非奇异的。
用乘以式(5.4.1)等号两边得(5.4.3)由上式得(5.4.4)如果取(5.4.5)作为估值，则称为辅助变量估值，矩阵成为辅助变量矩阵，中的元素称为辅助变量。
常用的辅助变量法有递推辅助变量参数估计法，自适应滤波法，纯滞后等。
5.广义最小二乘法广义最小二乘法是能克服最小二乘法有偏估计的另一种方法，这种方法计算比较复杂但效果比较好。
下面直接介绍广义最小二乘法的计算步骤：(1)应用得到的输入和输出数据和，按模型求出的最小二乘估计(2)计算残差(3)用残差代替，计算(4)计算和(5)应用得到的和按模型用最小二乘法重新估计，得到的第2次估值。
然后按步骤(2)计算残差，按步骤(3)重新估计，得到估值。
再按照步骤(4)计算和，按照步骤(5)求的第3次估值。
重复上述循环，之道的估值收敛为止。
二、实验要求设但输入-单输出系统的差分方程为取真实值，输入数据：X1=1;X2=1;X3=1;X4=1; %移位寄存器输入Xi初T态（1111），Yi为移位寄存器各级输出n=2000; %置M序列总长度for i=1:nY4=X4; Y3=X3; Y2=X2; Y1=X1;X4=Y3; X3=Y2; X2=Y1;X1=xor(Y3,Y4); %异或运算U(i)=Y4;if Y4==0U(i)=1;elseU(i)=Y4*(-1);endendu=U;用的真实值利用查分方程求出作为测量值，为均值为0，方差为0.1，0.5的不相关随机序列。
(1)用最小二乘法估计参数。
(2)用递推最小二乘法估计。
(3)用辅助变量法估计参数。
(4)设，用广义最小二乘法估计参数。
(5)详细分析和比较所获得的参数辨识结果，并说明上述参数便是方法的优缺点。
根据题目要求的解法，利用Matlab编程实现系统辨识的估值三、实验结果实验方法方差试验次数误差最小二乘0.5 1 0.1708 0.2070 0.3920 0.2089 1.81482 0.1680 0.1983 0.4008 0.2119 1.77483 0.1820 0.1724 0.3843 0.2347 1.82074 0.2230 0.1969 0.3857 0.2234 1.82390.1 1 0.8281 0.1373 0.3908 0.0335 1.04712 0.8005 0.1191 0.3866 0.0207 1.08243 0.7963 0.1147 0.3899 0.0179 1.08904 0.7828 0.1030 0.3887 0.0160 1.1065 0.01 1 1.6217 0.7008 0.3901 0.3418 0.02612 1.6224 0.7007 0.3893 0.3424 0.02543 1.6214 0.7008 0.3894 0.3410 0.02674 1.6225 0.7007 0.3898 0.3426 0.0253 0.001 1 1.6419 0.7149 0.3900 0.3499 0.00022 1.6417 0.7148 0.3899 0.3499 0.00033 1.6417 0.7149 0.3899 0.3498 0.00044 1.6418 0.7148 0.3900 0.3500 0.0002实验方法方差试验次数误差递推最小二乘0.5 1 0.2015 0.2059 0.3951 0.2111 1.86602 0.1634 0.1894 0.3997 0.2410 1.86603 0.2259 0.1747 0.3867 0.2027 1.86604 0.1797 0.1980 0.4151 0.2148 1.8660 0.1 1 0.8272 0.1364 0.3906 0.0300 1.68602 0.7653 0.0892 0.3940 0.0146 1.86603 0.8236 0.1311 0.3938 0.0340 1.86604 0.7829 0.1065 0.3892 0.0113 1.6800 0.01 1 1.6223 0.7007 0.3912 0.3430 0.02302 1.6236 0.7014 0.3908 0.3439 0.02453 1.6205 0.6999 0.3906 0.3419 0.02344 1.6224 0.7006 0.3905 0.3429 0.0236 0.001 1 1.6418 0.7148 0.3900 0.3500 0.0012 1.6419 0.7149 0.3900 0.3500 0.0023 1.6418 0.7149 0.3900 0.3500 0.00214 1.6417 0.7148 0.3899 0.3498 0.0031实验方法方差试验次数误差辅助变量0.5 1 0.7238 0.0910 0.3711 0.0015 1.86522 0.0015 0.1831 0.4055 0.1476 1.86603 1.2077 0.4235 0.4109 0.1793 1.86604 0.7185 0.1155 0.3840 0.0381 1.8593 0.1 1 1.6352 0.7062 0.3961 0.3557 1.86602 1.5911 0.6853 0.3955 0.3284 1.86603 1.6127 0.6959 0.3962 0.3411 1.86604 1.6318 0.7066 0.3905 0.3467 1.86600.01 1 1.6421 0.7150 0.3894 0.3500 0.23542 1.6420 0.7148 0.3903 0.3505 0.26123 1.6401 0.7144 0.3900 0.3487 0.24314 1.6423 0.7149 0.3892 0.3500 0.23450.001 1 1.6421 0.7150 0.3900 0.3501 0.00012 1.6420 0.7149 0.3900 0.3500 0.00023 1.6421 0.7150 0.3900 0.3500 0.00014 1.6421 0.7150 0.3900 0.3501 0.0001实验方法方差试验次数误差广义最小二乘0.5 1 0.9019 0.2338 0.3985 0.0448 0.97502 0.9713 0.2798 0.3953 0.0725 0.90233 0.8912 0.2437 0.4014 0.0229 0.78144 0.9983 0.3069 0.4035 0.0693 0.9341 0.1 1 1.6158 0.6994 0.3915 0.3391 0.03242 1.6140 0.6980 0.3931 0.3403 0.03433 1.6141 0.6977 0.3905 0.3383 0.03494 1.6194 0.7014 0.3927 0.3408 0.0281 0.01 1 1.6437 0.7162 0.3901 0.3509 0.00232 1.6444 0.7166 0.3898 0.3508 0.00303 1.6439 0.7164 0.3900 0.3505 0.00244 1.6441 0.7162 0.3900 0.3510 0.0026 0.001 1 1.6420 0.7150 0.3900 0.3500 0.00002 1.6419 0.7150 0.3900 0.3499 0.00013 1.6420 0.7150 0.3900 0.3500 0.00014 1.6419 0.7150 0.3900 0.3500 0.00011）输入信号图像：2）递推最小二乘图像3）辅助变量法最小二乘四、实验结论1)在噪声方差比较小的情况下，各种方法所获得的估值比较理想，但随着噪声方差的增大，估值的偏差随之增大。