系统辨识—最小二乘法汇总

航空航天飞行器动力学模型参数辨识方法对比航空航天飞行器的动力学模型参数辨识是飞行器设计、控制和仿真过程中的重要环节。

通过对飞行器的动力学特性进行建模和参数辨识，可以更好地理解和预测其飞行行为，并为飞行器的性能优化和控制提供可靠的依据。

在实际应用中，有多种方法可以用于航空航天飞行器动力学模型参数的辨识，本文将就常见的几种方法进行对比分析。

1. 最小二乘法（Least Squares Method）最小二乘法是一种经典的参数估计方法，在航空航天领域被广泛应用于动力学模型参数的辨识。

该方法通过最小化观测值与模型预测值之间的误差平方和，来确定参数的最佳估计值。

最小二乘法简单易用，计算效率高，能够处理噪声和测量误差的影响。

然而，该方法对数据的要求较高，需要大量高质量的测量数据。

此外，最小二乘法对数据的分布情况较为敏感，当数据存在离群点或者异常值时，估计结果可能不准确。

2. 基于频域的方法（Frequency-domain Methods）基于频域的方法是一种常用的动力学模型参数辨识方法，采用系统的频率响应数据进行模型参数的估计。

该方法通过对输入-输出频率特性的分析，可以得到系统的衰减特性、共振频率等重要参数。

基于频域的方法在辨识非线性系统和具有频率特性的系统上表现出良好的性能。

然而，该方法需要较为复杂的实验设置和信号处理技术，对测量设备和环境条件要求较高。

3. 系统辨识方法（System Identification Methods）系统辨识方法是一种通过对系统的输入-输出数据进行分析，确定系统动态行为和特性的方法。

该方法可以采用多种数学模型和算法，如ARMA模型、ARMAX模型、神经网络模型等，来描述系统的动力学特性。

系统辨识方法具有较高的灵活性和适用性，能够处理非线性系统和时变系统的参数辨识问题。

然而，系统辨识方法的参数辨识过程和计算复杂度较高，需要对模型结构和参数进行合理选择，以获得准确的辨识结果。

系统辨识之最小二乘法方法一、最小二乘一次性算法：首先对最小二乘法的一次性辨识算法做简要介绍如下：过程的黑箱模型如图所示：其中u(k)和z(k)分别是过程的输入输出，)(1-z G 描述输入输出关系的模型，成为过程模型。

过程的输入输出关系可以描述成以下最小二乘格式：)()()(k n k h k z T +=θ （1）其中z(k)为系统输出，θ是待辨识的参数，h(k)是观测数据向量，n(k)是均值为0的随机噪声。

利用数据序列{z （k ）}和{h （k ）}极小化下列准则函数：∑=-=Lk T k h k z J 12])()([)(θθ （2）使J 最小的θ的估计值^θ，成为最小二乘估计值。

具体的对于时不变SISO 动态过程的数学模型为 )()()()()(11k n k u z B k z z A +=-- （3）应该利用过程的输入、输出数据确定)(1-z A 和)(1-Z B 的系数。

对于求解θ的估计值^θ，一般对模型的阶次a n ，b n 已定，且b a n n >；其次将（3）模型写成最小二乘格式)()()(k n k h k z T +=θ （4）式中=------=T n n T b a b a b b b a a a n k u k u n k z k z k h ],,,,,,,[)](,),1(),(,),1([)(2121 θ （5）L k ,,2,1 =因此结合式（4）（5）可以得到一个线性方程组L L L n H Z +=θ （6）其中==T L TL L n n n n L z z z z )](),2(),1([)](),2(),1([ （7）对此可以分析得出，L H 矩阵的行数为),max(b a n n L -，列数b a n n +。

在过程的输入为2n 阶次，噪声为方差为1，均值为0的随机序列，数据长度)(b a n n L +>的情况下，取加权矩阵L Λ为正定的单位矩阵I ，可以得出：L T L L T L z H H H 1^)(-=θ （8）其次，利用在Matlab 中编写M 文件，实现上述算法。

控制系统设计中的模型鉴别方法综述在控制系统设计中，模型鉴别方法是一项关键性工作。

模型鉴别方法可以帮助工程师准确地识别出待控系统的数学模型，为后续的控制器设计和性能优化提供基础。

本文将对控制系统设计中常用的模型鉴别方法进行综述。

一、最小二乘法最小二乘法是一种常见的模型鉴别方法，它通过最小化误差的平方和来拟合实际测量数据和理论模型之间的差异。

最小二乘法可以用于线性和非线性模型的鉴别。

对于线性模型，最小二乘法可以通过矩阵运算求解最优解。

而对于非线性模型，最小二乘法可以通过迭代优化算法求解。

二、频域方法频域方法是一种将系统响应与频率特性相关联的模型鉴别方法。

它通常基于输入和输出信号的频谱分析，可以用于连续时间和离散时间系统。

频域方法可以采用傅里叶变换、拉普拉斯变换等数学工具，通过求解传递函数或频率响应函数来获得系统模型。

频域方法适用于具有周期性输入和输出信号的系统。

三、时域方法时域方法是一种将系统响应与时间域特性相关联的模型鉴别方法。

它通常基于实际采集到的离散时间数据，通过插值、拟合等技术来获得离散时间系统的模型。

时域方法可以采用多项式插值、曲线拟合等数学工具，通过建立系统差分方程或状态空间模型来进行模型鉴别。

时域方法适用于实际工程中获得的离散时间数据。

四、系统辨识方法系统辨识方法是一种通过试验数据来识别系统动态特性的模型鉴别方法。

它可以通过对系统施加特定的输入信号，观测系统输出响应来获得系统模型。

系统辨识方法可以分为参数辨识和非参数辨识两种方法。

参数辨识方法假设系统具有某种结构，通过最小化残差的平方和来确定模型参数。

非参数辨识方法不对系统结构进行假设，通过直接拟合试验数据来获得系统模型。

五、神经网络方法神经网络方法是一种基于人工神经网络的模型鉴别方法。

它可以通过输入输出数据训练神经网络，从而获得系统的模型。

神经网络方法可以适用于非线性系统的建模和鉴别。

神经网络方法具有较强的自适应能力和非线性拟合能力，但对于网络结构和训练样本的选择具有一定的要求。

——第三章辨识方法
2 最小二乘法
1
清华大学电机系
辨识技术
辨识技术
()
*
lim k k α
α
→∞
=辨识技术
远离最小值点
接近最小值点
最速下降法
高斯牛顿法
清华大学电机系
辨识技术
42
2
三种算法的比较
y
香蕉函数最小化问题：
f (x ) = 100 ( X 2–X 12)2+ ( 1 –X 1)2
高斯-牛顿法：11次迭代
最速下降法：60次迭代
阻尼最小二乘法：18次迭代
清华大学电机系
辨识技术
43
例：非线性函数的参数辨识
x
C e
C x f ⋅⋅=21)(高斯-牛顿法：80次迭代
最速下降法：401次依然未收敛
阻尼最小二乘法：10次迭代
真值C * = [5.420136187 -0.25436189]
当x m = [ 1 2 3 4 5 ]时
测量值：f m (x ) = [ 4.20 3.25 2.52 1.95 1.51 ]f *(x ) = [ 4.202834 3.258924 2.527006 1.959469 1.519394 ]
真值
清华大学电机系
辨识技术
清华大学电机系辨识技术
清华大学电机系辨识技术
辨识技术
稳态电路。

《系统辨识基础》第17讲要点第5章 最小二乘参数辨识方法5.9 最小二乘递推算法的逆问题辨识是在状态可测的情况下讨论模型的参数估计问题，滤波是在模型参数已知的情况下讨论状态估计问题，两者互为逆问题。

5.10 最小二乘递推算法的几种变形最小二乘递推算法有多种不同的变形，常用的有七种情况：① 基于数据所含的信息内容不同，对数据进行有选择性的加权； ② 在认为新近的数据更有价值的假设下，逐步丢弃过去的数据； ③ 只用有限长度的数据；④ 加权方式既考虑平均特性又考虑跟综能力； ⑤ 在不同的时刻，重调协方差阵P (k )； ⑥ 设法防止协方差阵P (k )趋于零； 5.10.1 选择性加权最小二乘法 把加权最小二乘递推算法改写成[]⎪⎪⎩⎪⎪⎨⎧--=+--=--+-=-)1()]()([)(1)()1()()()()1()()()]1(ˆ)()()[()1(ˆ)(ˆ1k k k k k k k k k k k k k k k z k k k P h K P h P h h P K h K τττθθθI ΛΛ算法中引进加权因子，其目的是便于考虑观测数据的可信度．选择不同的加权方式对算法的性质会有影响，下面是几种特殊的选择：① 一种有趣的情况是Λ()k 取得很大，在极限情况下，算法就退化成正交投影算法。

也就是说，当选择⎩⎨⎧=-≠-∞=0)()1()(,00)()1()(,)(k k k k k k k h P h h P h ττΛ 构成了正交投影算法⎪⎪⎩⎪⎪⎨⎧--=--=--+-=)1()]()([)()()1()()()1()()]1(ˆ)()()[()1(ˆ)(ˆk k k k k k k k k k k k k z k k k P h K P h P h h P K h K τττθθθI 算法初始值取P ()0=I 及 ()θε0=(任定值)，且当0)()1()(=-k k k h P h τ时，令K ()k =0。

---------------------------------------------------------------最新资料推荐------------------------------------------------------系统辨识—最小二乘法最小二乘法参数辨识 1 引言系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。

现代控制理论中的一个分支。

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

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

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

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

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

系统辨识包括两个方面：结构辨识和参数估计。

在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。

2 系统辨识的目的在提出和解决一个辨识问题时，明确最终使1 / 17用模型的目的是至关重要的。

它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。

通过辨识建立数学模型通常有四个目的。

①估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。

这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。

②仿真仿真的核心是要建立一个能模仿真实系统行为的模型。

用于系统分析的仿真模型要求能真实反映系统的特性。

最小二乘法参数辨识201403027摘要:系统辨识在工程中的应用非常广泛,系统辨识的方法有很多种,最小二乘法是一种应用极其广泛的系统辨识方法.阐述了动态系统模型的建立及其最小二乘法在系统辨识中的应用,并通过实例分析说明了最小二乘法应用于系统辨识中的重要意义.关键词:最小二乘法;系统辨识;动态系统Abstract: System identification in engineering is widely used, system identification methods there are many ways， least squares method is a very wide range of application of system identification method and the least squares method elaborated establish a dynamic system models in System Identification applications and examples analyzed by the least squares method is applied to illustrate the importance of system identification.Keywords: Least Squares; system identification; dynamic system引言随着科学技术的不断发展,人们认识自然、利用自然的能力越来越强,对于未知对象的探索也越来越深入.我们所研究的对象,可以依据对其了解的程度分为三种类型:白箱、灰箱和黑箱.如果我们对于研究对象的内部结构、内部机制了解很深入的话,这样的研究对象通常称之为“白箱”;而有的研究对象,我们对于其内部结构、机制只了解一部分,对于其内部运行规律并不十分清楚,这样的研究对象通常称之为“灰箱”;如果我们对于研究对象的内部结构、内部机制及运行规律均一无所知的话,则把这样的研究对象称之为“黑箱”.研究灰箱和黑箱时,将研究的对象看作是一个系统,通过建立该系统的模型,对模型参数进行辨识来确定该系统的运行规律.对于动态系统辨识的方法有很多,但其中应用最广泛,辨识效果良好的就是最小二乘辨识方法,研究最小二乘法在系统辨识中的应用具有现实的、广泛的意义.1．1 系统辨识简介系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。

最小二乘法一次完成算法的MATLAB仿真一、二阶系统的最小二乘法一次完成算法的辨识程序及结果程序源代码：clear % 工作间清零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.8*z(k-1) - 0.8*z(k-2) + 1.1*u(k-1) + 0.4*u(k-2);ZL( k, : ) = z(k); % 给样本矩阵ZL赋值Hl = [ -z(k-1), -z(k-2), u(k-1), u(k-2) ]; % 用理想输出值作为观测值HL( k, : ) = Hl;endsubplot(3,1,1) % 画三行一列图形窗口中的第一个图形stem(u) % 画出输入信号u的经线图形subplot(3,1,2) % 画三行一列图形窗口中的第二个图形i=1:1:16; % 横坐标范围是1到16，步长为1% 图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线plot(i,z)subplot(3,1,3) % 画三行一列图形窗口中的第三个图形stem(z),grid on % 画出输出观测值z的经线图形，并显示坐标网格u,z % 显示输入信号和输出观测信号c1=HL'*HL; % 计算参数并显示c2=inv(c1);c3=HL'*ZL;c=c2*c3a1=c(1), a2=c(2), b1=c(3), b2=c(4) % 从中分离出并显示a1 、a2、b1、b2运算结果如下：u =-1 1 -1 1 1 1 1 -1 -1 -1 1 -1 -1 1 1z =Columns 1 through 90 0 0.7000 0.5600 1.1480 3.1184 6.1947 10.1558 12.6246 Columns 10 through 1613.0997 11.9798 11.7838 10.9270 8.7416 7.6933 8.3546c =-1.80000.80001.10000.4000a1 =-1.8000a2 =0.8000b1 =1.1000b2 =0.4000所画图形如下：二、实际压力系统的最小二乘辨识程序及结果程序源代码：clear % 工作间清零V = [ 54.3, 61.8, 72.4, 88.7, 118.6, 194.0]' % 赋初值V，并显示P = [ 61.2, 49.5, 37.6, 28.4, 19.2, 10.1]' % 赋初值P，并显示% P、V之间的关系:logP=-alpha*logV% +logbeita=[-logV,1][alpha,log(beita)]'=HL*sitafor i = 1 : 6; % 循环变量的取值为从1到6 Z(i) = log( P (i) ); % 赋系统的输出采样值ZL( i, : ) = Z(i); % 给样本矩阵ZL赋值Hl = [ log( V(i) ),1 ]; % 给HL赋值HL( i, : ) = Hlend % 循环结束c1=HL'*HL; % 计算参数并显示c2=inv(c1);c3=HL'*ZL;c4=c2*c3alpha = c4(1) % 为c4的第1个元素beita = exp(c4(2)) % 为以自然数为底的c4的第2个元素的指数运算结果如下：V =54.300061.800072.400088.7000118.6000194.0000P =61.200049.500037.600028.400019.200010.1000HL =3.9945 1.0000HL =3.9945 1.00004.1239 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.00004.7758 1.0000HL =3.9945 1.00004.1239 1.00004.2822 1.00004.4853 1.00004.7758 1.00005.2679 1.0000 c4 =-1.40429.6786alpha =-1.4042beita =1.5972e+004。

5.3 Matlab源程序：%最小二乘估计clearu=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1z=zeros(1,30);for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endh0=[-z(2) -z(1) u(2) u(1)]';HLT=[h0,zeros(4,28)];for k=3:30h1=[-z(k) -z(k-1) u(k) u(k-1)]';HLT(:,k-1)=h1;endHL=HLT';y=[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) ;z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29); z(30);z(31)];%求出FAIc1=HL'*HL;c2=inv(c1);c3=HL'*y;c=c2*c3;%display('方差=0.1时，最小二乘法估计辨识参数θ如下：');a1=c(1)a2=c(2)b1=c(3)b2=c(4)%递推最小二乘法估计clearu=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];z(2)=0;z(1)=0;n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endc0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(4,30)]; %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求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; %给下次用if e2<=Ebreak; %如果参数收敛情况满足要求，终止计算endend%display('方差为0. 1递推最小二乘法辨识后的结果是：');a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');for i=3:31;if(c(1,i)==0)q1=c(1,i-1);break;endendfor i=3:31;if(c(2,i)==0)q2=c(2,i-1);break;endendfor i=3:31;if(c(3,i)==0)q3=c(3,i-1);break;endendfor i=3:31;if(c(4,i)==0)q4=c(4,i-1);break;endenda1=q1a2=q2b1=q3b2=q4%辅助变量递推最小二乘法估计clearna=2;nb=2;siitt=[1.642 0.715 0.39 0.35]';siit0=0.001*eye(na+nb,1);p=10^6*eye(na+nb);siit(:,1)=siit0;y(2)=0;y(1)=0;x(1)=0;x(2)=0;j=0;u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1for k=3:31;h=[-y(k-1),-y(k-2),u(k-1),u(k-2)]';y(k)=h'*siitt+n(k)+1.642*n(k-1)+0.715*n(k-2);hx=[-x(k-1),-x(k-2),u(k-1),u(k-2)]';kk=p*hx/(h'*p*hx+1);p=[eye(na+nb)-kk*h']*p;siit(:,k-1)=siit0+kk*[y(k)-h'*siit0];x(k)=hx'*siit(:,k-1);j=j+(y(k)-h'*[1.642 0.715 0.39 0.35]')^2;e=max(abs((siit(:,k-1)-siit0)./siit0));ess(:,k-2)=siit(:,k-1)-siitt;siit0=siit(:,k-1);enda1=siit0(1)a2=siit0(2)b1=siit0(3)b2=siit0(4)clear%广义最小二乘估计clear;nn = normrnd(0,sqrt(0.1),1,31)'; %方差为0.1uk=[1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.626 0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 1.144 1.177 -0.390];yk(1)=0;yk(2)=0;for i=1:29;yk(i+2)=-1.642*yk(i+1)-0.715*yk(i)+0.39*uk(i+1)+0.35*uk(i)+nn(i+2)+1.642*nn(i+1 )+0.715*nn(i);end;for i=1:29;A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];endsiit=inv(A'*A)*A'*(yk(3:31)+nn(2:30)')';e(1)=yk(1);e(2)=yk(2)+siit(1)*yk(1)-siit(3)*uk(1);for i=3:31;e(i)=yk(i)+siit(1)*yk(i-1)+siit(2)*yk(i-2)-siit(3)*uk(i-1)-siit(4)*uk(i-2); endfor i=1:29;fai(i,:)=[-e(i+1) -e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';for i=3:31;yk(i)=yk(i)+f(1)*yk(i-1)+f(2)*yk(i-2);endyk(2)=yk(2)+f(1)*yk(1);for i=3:30;uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);for j=1:30for i=1:29;A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)];endsiit=inv(A'*A)*A'*yk(3:31)';e(1)=yk(1);e(2)=yk(2)+siit(1)*(yk(1))-siit(3)*uk(1);for i=3:31;e(i)=yk(i)+siit(1)*(yk(i-1))+siit(2)*(yk(i-2))-siit(3)*uk(i-1)-siit(4)*uk(i-2); endfor i=1:29;fai(i,:)=[-e(i+1) -e(i)];endf=inv(fai'*fai)*fai'*e(3:31)';k1(j)=f(1);k2(j)=f(2);for i=3:31;yk(i)=yk(i)+f(1)*(yk(i-1))+f(2)*(yk(i-2));endyk(2)=yk(2)+f(1)*yk(1);for i=3:30uk(i)=uk(i)+f(1)*uk(i-1)+f(2)*uk(i-2);enduk(2)=uk(2)+f(1)*uk(1);endsiit'%增广矩阵估计参数clearu=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1z=zeros(7,30);zs=zeros(7,30);zm=zeros(7,30);zmd=zeros(7,30);%输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出矩阵的大小z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;zmd(2)=0;zmd(1)=0;%给输出采样、不考虑噪声时系统输出、不考虑噪声时模型输出、模型输出赋初值for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endc0=[0.001 0.001 0.001 0.001]'; %直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(4,4); %直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005; %取相对误差E=0.000000005c=[c0,zeros(4,30)]; %被辨识参数矩阵的初始值及大小e=zeros(4,30); %相对误差的初始值及大小for k=3:30; %开始求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; %求被辨识参数czs(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2);%系统在M序列的输入下不考虑扰动时的输出响应zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)];%模型在输入下不考虑扰动时的的输出响应zmd(k)=h1'*c1;%模型在输入下的的输出响应e1=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; %给下次用if e2<=Ebreak; %如果参数收敛情况满足要求，终止计算endend%display('方差为0. 1递推最小二乘法辨识后的结果是：');a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);%display('a1,a2,b1,b2经过递推最小二乘法辨识的结果是：');for i=3:31;if(c(1,i)==0)q1=c(1,i-1);break;endendfor i=3:31;if(c(2,i)==0)q2=c(2,i-1);break;endendfor i=3:31;if(c(3,i)==0)q3=c(3,i-1);break;endendfor i=3:31;if(c(4,i)==0)q4=c(4,i-1);break;endenda1=q1a2=q2b1=q3b2=q4%用夏氏偏差修正法估计参数clear;u=[ 1.147 0.201 -0.787 -1.589 -1.052 0.866 1.152 1.573 0.6260.433 -0.985 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.4501.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 0.144 1.177-0.390];n=normrnd(0, sqrt(0.1), 1, 31); %方差为0.1z=zeros(1,30);for k=3:31z(k)=-1.642*z(k-1)-0.715*z(k-2)+0.39*u(k-1)+0.35*u(k-2)+n(k)+1.642*n(k-1)+0.715 *n(k-2);endh0=[-z(2) -z(1) u(2) u(1)]';HLT=[h0,zeros(4,28)];for k=3:30h1=[-z(k) -z(k-1) u(k) u(k-1)]';HLT(:,k-1)=h1;endHL=HLT';y=[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) ;z(17);z(18);z(19);z(20);z(21);z(22);z(23);z(24);z(25);z(26);z(27);z(28);z(29); z(30);z(31)];%求出FAIc1=HL'*HL;c2=inv(c1);c3=HL'*y;c=c2*c3;fai1_xiashi=HL;siit1_ls_guzhi=c;ka_ma=inv(fai1_xiashi'*fai1_xiashi)*fai1_xiashi';M_xiashi=eye(29)-fai1_xiashi*ka_ma;siit_N0_xshi=siit1_ls_guzhi;siit_bxs=[10,10,10,10]';for n=1:30e=zeros(29,1);e(3)=y(5)+siit_N0_xshi(1,1)*y(4)+siit_N0_xshi(2,1)*y(3)-siit_N0_xshi(3,1)*u(4)-siit_N0_xshi(4,1)*u(3);e(4)=y(6)+siit_N0_xshi(1,1)*y(5)+siit_N0_xshi(2,1)*y(4)-siit_N0_xshi(3,1)*u(5)-siit_N0_xshi(4,1)*u(4);e(5:29) = y(5:29) - fai1_xiashi(5:29,:) * siit_N0_xshi;ou_mei_ka=zeros(29,1);for i=1:29ou_mei_ka(i,:)=-e(i);endD_xiashi=ou_mei_ka'*M_xiashi*ou_mei_ka;f_xiashi=inv(D_xiashi)*ou_mei_ka'*M_xiashi*y(1:29);siit_b=ka_ma*ou_mei_ka*f_xiashi;siit_xiashi=siit1_ls_guzhi-siit_b;if norm(siit_b-siit_bxs)/norm(siit_bxs)<0.0001break ;endsiit_N0_xshi=siit_xiashi;siit_bxs=siit_b;endsiit1_xiashi=siit_xiashi;%display('利用夏氏修正法得出的辨识结果：');a1=siit1_xiashi(1)a2=siit1_xiashi(2)b1=siit1_xiashi(3)b2=siit1_xiashi(4)。

matlab系统辨识工具箱使用的算法MATLAB的系统辨识工具箱使用多种算法来进行系统辨识。

这些算法通常包括以下几种：
1. 最小二乘法（Least Squares）：这是最常用的系统辨识方法。

最小二乘法试图找到一组参数，使得实际数据和模型预测之间的误差平方和最小。

2. 极大似然估计（Maximum Likelihood Estimation）：这种方法基于数据生成的模型概率密度函数，通过最大化似然函数来估计模型参数。

3. 递归最小二乘法（Recursive Least Squares）：这是一种在线算法，可以在数据流中实时更新模型参数。

4. 扩展最小二乘法（Extended Least Squares）：这种方法可以处理包含噪声和异常值的数据，通过引入权重来调整误差平方和。

5. 非线性最小二乘法（Nonlinear Least Squares）：对于非线性系统，需要使用非线性最小二乘法来估计参数。

6. 遗传算法（Genetic Algorithms）：这是一种启发式搜索算法，通过模拟自然选择和遗传过程来寻找最优解。

7. 粒子群优化算法（Particle Swarm Optimization）：这是一种基于群体智能的优化算法，通过模拟鸟群、鱼群等生物群体的行为来寻找最优解。

以上这些算法都是MATLAB系统辨识工具箱中常用的算法，根据具体的问题和数据，可以选择最适合的算法来进行系统辨识。