系统辨识广义最小二乘
系统辨识相关分析最小二乘
相关分析法辨识系统单位脉冲响应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 ≥。
系统辨识之最小二乘法
系统辨识之最小二乘法方法一、最小二乘一次性算法：首先对最小二乘法的一次性辨识算法做简要介绍如下：过程的黑箱模型如图所示：其中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 文件，实现上述算法。
基于广义最小二乘法的系统参数辨识与仿真
基于广义最小二乘法的系统参数辨识与仿真摘要：系统辨识是自动控制学科的一个重要分支，由于其特殊作用，已经广泛应用于各种领域，尤其是复杂系统或参数不容易确定的系统的建模。
过去，系统辨识主要用于线性系统的建模，经过多年的研究，已经形成成熟的理论。
但随着社会、科学的发展，非线性系统越来越受到人们的关注，其控制与模型之间的矛盾越来越明显，因而非线性系统的辨识问题也越来越受到重视，其便是理论不断发展和完善。
本文重点介绍了系统参数辨识中广义最小二乘法的基本原理，具体说明了基于广义最小二乘法参数辨识在Matlab中的实现方法，结合实例给出相应的仿真程序及结果分析。
关键词：系统辨识；参数辨识；广义最小二乘法；Matlab1.引言所谓辨识就是通过测取研究对象在人为输入作用下的输出响应，或正常运行时的输入输出数据记录，加以必要的数据处理和数学计算，估计出对象的数学模型。
这是因为对象的动态特性被认为必然表现在它的变化着的输入输出数据之中，辨识只不过是利用数学的方法从数据序列中提炼出对象的数学模型而已。
在系统辨识领域中，最小二乘法是最基本最常用的方法。
可用于动态、静态、线性、非线性系统。
这种方法只适用于噪声是不相关随机序列时才是无偏估计，但大多数情况下噪声却是相关随机序列。
所以本文讨论克服最小二乘法有偏估计的一种方法—广义最小二乘法。
2.系统辨识一般而言，建立系统的数学模型有两种方法：激励分析法和系统辨识法。
前者是按照系统所遵循的物化（或社会、经济等）规律分析推导出模型。
后者则是从实际系统运行和实验数据处理获得模型。
如图1所示，系统辨识就是从系统的输入输出数据测算系统数学模型的理论和方法。
更进一步的定义是L.A.Zadeh曾经与1962年给出的，即“系统辨识是在输入和输出的基础上，从系统的一类系统范围内，确立一个与所实验系统等价的系统。
”另外，系统辨识还应该具有3个基本要素，即模型类、数据和准则。
被辨识系统模型根据模型形式可分为参数模型和非参数模型两大类。
第六章 最小二乘类参数辨识方法
《系统辨识基础》第20讲要点第6章 最小二乘类参数辨识方法6.1 引言最小二乘法是一种最基本的辨识方法，但如果模型的噪声不是白噪声，最小二乘法不一定能给出无偏、一致估计。
以下着重讨论模型噪声是有色噪声时的各种最小二乘辨识方法。
6.2 增广最小二乘法 6.2.1 增广最小二乘原理 考虑如下模型A z z kB z u k N z v k ()()()()()()---=+111式中u (k )和z (k ) 分别为模型输入和输出变量；v (k ) 是均值为零、方差为σv 2的不相关随机噪声或称白噪声；N z ()-1为噪声模型；A z ()-1 和B z ()-1为迟延算子多项式，记作A z a z a z a zB z b z b z b z n n n n a ab b()()--------=++++=+++⎧⎨⎪⎩⎪11122111221 其中n a 和n b 为模型阶次。
为了运用最小二乘原理来辨识这种模型的参数，需要把模型(4-1)式写成最小二乘格式)()()(k v k k z +=θτh这样就必须把噪声模型的参数包含在参数向量θ 中，从而引出增广概念，用来构造模型的参数向量θ和数据向量h(k )，具体的构成形式会因噪声模型的结构不同而不同。
下面是三种不同噪声模型的向量构成方法：① 若N z D z d z d z d z n n dd()()==++++----111221 ，可按下式构成参数向量和数据向量：⎪⎪⎩⎪⎪⎨⎧--------+==--------=∑∑∑===)(ˆ)1(ˆ)()1(ˆ)()1(ˆ)()(ˆ],,,,,,,,[)](ˆ,),1(ˆ),(,),1(),(,),1([)(111111i k v k d i k u k b i k z k a k z k v d d b b a a n k v k v n k u k u n k z k z k db a d b a n i i n i i n i i n n n d b a ττθ h② 若N z C zc zc zc zn n c c()()==++++----11111122，参数向量和数据向量的构成形式为：⎪⎪⎩⎪⎪⎨⎧-----+==----------=∑∑==)()1(ˆ)()1(ˆ)()(ˆ],,,,,,,,[)](ˆ,),1(ˆ),(,),1(),(,),1([)(11111i k u k b i k z k a k z k e c c b b a a n k e k e n k u k u n k z k z k ba cb a n i i n i i n n nc ba ττθ h③ 若N z D z C zd z d z d z c zc zc zn n n n d dc c()()()==++++++++--------111122112211 ，参数向量和数据向量的构成形式为：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧-----+=-----+==------------=∑∑∑∑====)()1(ˆ)()1(ˆ)()(ˆ)(ˆ)1(ˆ)(ˆ)1(ˆ)(ˆ)(ˆ],,,,,,,,,,,[)](ˆ,),1(ˆ),(ˆ,),1(ˆ),(,),1(),(,),1([)(11111111i k u k b i k z k a k z k e i k v k d i k e k c k e k v d d c c b b a a n k v k v n k e k e n k u k u n k z k z k b a dc cc b a n i i n i i n i i n i i n n n nd c b a ττθ h以上这种构成参数向量和数据向量的思想就是所谓的增广原理，它是增广最小二乘法的根本。
系统辨识最小二乘法大作业 (2)
系统辨识大作业最小二乘法及其相关估值方法应用学院：自动化学院学号：姓名：日期：基于最小二乘法的多种系统辨识方法研究一、实验原理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) 为最小来确定估值。
广义最小二乘法
4.5 广义最小二乘法（GLS ） GLS----Generalized Least Squares 1. 基本原理广义最小二乘法的基本思想在于引入一个所谓成形滤波器（白化滤波器），把相关噪声)(k ξ转化成白噪声)(k ε。
由方程（4-4）、（4-5），系统的差分方程可以表示为)()()()()(11k k u z b k y z a ξ+=-- （4-114）式中n n z a z a z a z a ----++++=ΛΛ221111)(nn z b z b z b b z b ----++++=ΛΛ221101)(如果知道有色噪声序列)(k ξ的相关性，则可以把)(k ξ看成白噪声通过线性系统后所得的结果。
这种线性系统通常称为成形滤波器，其差分方程为)()()()(11_k z d k zc εξ---= （4-115）式中)(k ε是均值为零的白噪声序列，)()(11_---z d 、z c 是1-z 的多项式。
令 _111212_1()()1()m m c z f z f z f z f z d z ------==+++L L （4-116）有 )()(1)()()()(11k z f k k k z f εξεξ--==或 （4-117）即1212(1)()()m m f z f z f z k k ξε---++++=L L （4-118）或)()()2()1()(21k m k f k f k f k m εξξξξ+-------=ΛΛ ()1,,n k n N =++L L（4-119）这一噪声模型（自回归模型）的阶m ，一般事先是不知道的，实际经验表明，若指定m为2或3，就可以获得令人满意的描述)(k ξ的模型。
把方程（4-119）看作输入为零的差分方程，并由此式来写出N 个方程。
⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧++-+---+--+-=+++-+---+-=+++-+-----=+)()()2()1()()2()2()()1()2()1()1()1()()1(212121N n m N n f N n f N n f N n n m n f n f n f n n m n f n f n f n m m m εξξξξεξξξξεξξξξΛΛM ΛΛΛΛ写成向量矩阵形式为εξ+Ω=f （4-120）其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡++=)()1(N n n ξξξM ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=m f f f M 1，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡++=)()1(N n n εεεM ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-+--+--+--+--+--+----=Ω)()2()1()2()()1()1()1()(m N n N n N n m n n n m n n n ξξξξξξξξξM Λ(4-120)式所示的线性组合关系是辨识问题的基本表达形式,称作最小二乘格式。
系统辨识第5-7讲
《系统辨识》第5讲要点第5章 线性动态模型参数辨识－最小二乘法5.1 辨识方法分类根据不同的辨识原理，参数模型辨识方法可归纳成三类：① 最小二乘类参数辨识方法，其基本思想是通过极小化如下准则函数来估计模型参数：其中代表模型输出与系统输出的偏差。
典型的方法有最小二乘法、增广最小二乘法、辅助变量法、广义最小二乘法等。
② 梯度校正参数辨识方法，其基本思想是沿着准则函数负梯度方向逐步修正模型参数，使准则函数达到最小，如随机逼近法。
③ 概率密度逼近参数辨识方法，其基本思想是使输出z 的条件概率密度最大限度地逼近条件下的概率密度，即。
典型的方法是极大似然法。
5.2 最小二乘法的基本概念● 两种算法形式① 批处理算法：利用一批观测数据，一次计算或经反复迭代，以获得模型参数的估计值。
② 递推算法：在上次模型参数估计值的基础上，根据当前获得的数据提出修正，进而获得本次模型参数估计值，广泛采用的递推算法形式为其中表示k 时刻的模型参数估计值，K(k)为算法的增益，h(k-d)是由观测数据组成的输入数据向量，d为整数，表示新息。
● 最小二乘原理定义：设一个随机序列的均值是参数 的线性函数其中h(k)是可测的数据向量，那么利用随机序列的一个实现，使准则函数达到极小的参数估计值称作的最小二乘估计。
● 最小二乘原理表明，未知参数估计问题，就是求参数估计值，使序列的估计值尽可能地接近实际序列，两者的接近程度用实际序列与序列估计值之差的平方和来度量。
● 如果系统的输入输出关系可以描述成如下的最小二乘格式式中z(k)为模型输出变量，h(k)为输入数据向量，为模型参数向量，n(k)为零均值随机噪声。
为了求此模型的参数估计值，可以利用上述最小二乘原理。
根据观测到的已知数据序列和，极小化下列准则函数即可求得模型参数的最小二乘估计值。
● 最小二乘估计值应在观测值与估计值之累次误差的平方和达到最小值处，所得到的模型输出能最好地逼近实际系统的输出。
系统辨识—最小二乘法_3
---------------------------------------------------------------最新资料推荐------------------------------------------------------系统辨识—最小二乘法最小二乘法参数辨识 1 引言系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。
现代控制理论中的一个分支。
通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。
对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。
对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。
而系统辨识所研究的问题恰好是这些问题的逆问题。
通常，预先给定一个模型类={M}（即给定一类已知结构的模型），一类输入信号 u 和等价准则 J=L(y，yM)(一般情况下，J 是误差函数，是过程输出 y 和模型输出 yM 的一个泛函)；然后选择使误差函数J 达到最小的模型，作为辨识所要求的结果。
系统辨识包括两个方面：结构辨识和参数估计。
在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
2 系统辨识的目的在提出和解决一个辨识问题时，明确最终使1 / 17用模型的目的是至关重要的。
它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。
通过辨识建立数学模型通常有四个目的。
①估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。
这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。
②仿真仿真的核心是要建立一个能模仿真实系统行为的模型。
用于系统分析的仿真模型要求能真实反映系统的特性。
系统辨识各类最小二乘法汇总
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 =
系统辨识之最小二乘
写成矩阵的形式为：
Y = Fq + e
3.6
é - y(n)
- y(n -1) … - y(1) u(n) … u(1) ù
F
=
ê ê ê
-
y(n "
-1)
- y(n) "
… - y(2)
u(n +1)
…
u(2)
ú ú
…"
" # "ú
êë- y(n + N ) - y(n + N -1) ! - y(N ) u(n + N ) ! u(N )úû
q = F -1y
如果噪声x ¹ 0，则 q = F-1y - F-1x
从上式可以看出噪声 x 对参数估计有影响，为了尽量减小噪声 x 对q 估值的影响，应取 N>(2n+1), 即方程数大于未知数数目。在这种情况下，不能用解方程的办法来求q ，而要采用数理统计的办法， 以便减小噪声对q 估值的影响。在给定输出向量 y 和测量矩阵 F 的条件下求系统参数q 的估值，这 就是系统辨识问题。可用最小二乘法来就q 的估值。
在科学研究中，为了揭示某些相关量之间的关系，找出其规律，往往需要求解其函数解析式。
一种方法是采用插值逼近法，即所构造的近似函数 f (x) 在已知节点 xi 上必须满足 f(xi) = yi 要求逼 近函数 f (xi) 与被逼近函数 f (x) 在各已知点 xi 处的误差为零，即要求 f (x) 的曲线必须通过所有的
取泛函 J (q )为
N
å å J (q ) = (Y - Fq )2 = e2 (n + i) = eT • e = (Y - Fq )T (Y - Fq ) i =1
广义递推最小二乘辨识
广义递推最小二乘辨识一、实验目的1 通过实验掌握广义最小二乘辨识算法；2 运用MATLAB编程，掌握算法实现方法。
二、实验原理广义最小二乘法的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。
如果滤波模型选择得合适，对数据进行了较好的白色化处理，那么直接利用普通最小二乘法就能获得无偏一致估计。
广义最小二乘法所用的滤波模型实际上就是一种动态模型，在整个迭代过程中不断靠偏差信息来调整这个滤波模型，使它逐渐逼近于一个较好的滤波模型，以便对数据进行较好的白色化处理，使模型参数估计称为无偏一致估计。
理论上说，广义最小二乘法所用的动态模型经过几次迭代调整后，便可对数据进行较好的白化处理，但是，当过程的输出噪信比比较大或模型参数比较多时，这种数据白色化处理的可靠性就会下降。
此时，准则函数可能出现多个局部收敛点，因而辨识结果可能使准则函数收敛于局部极小点上而不是全局极小点上。
这样，最终的辨识结果往往也会是有偏的。
其收敛速度比较慢，需要经过多次迭代计算，才能得到较准确的参数估计值。
一般情况下，经过多次迭代后，估计值便会收敛到稳态值。
但在某些情况下（如噪声比较低时）存在局部极小值，估计值不一定收敛到准则函数的全局极小值上。
为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。
如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。
三、实验内容<1> 数据获取：实验数据按照表9-1，为二阶线性离散系统的输入输出数据<2> 数据处理：为了提高辨识精度，实验者必须对原始数据进行剔除坏数据、零均值化、工频滤波等处理。
实验进行了白化滤波处理。
<3> 辨识算法：利用处理过的数据（取适当的数据长度），选择某种辨识方法（如RLS递推最小二乘法、RELS、RIV或RML等参数估计算法及F-检验或AIC定阶法），估计出模型参数和阶次，同时分析辨识结果。
系统辨识-最小二乘法
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)。
二阶系统的最小二乘辨识和广义预测控制仿真
二阶系统的最小二乘辨识和广义预测控制仿真1. 二阶系统模型1.1 系统模型及输入输出变量设传递函数为()G s =25+31s s +，()1H s =，且为负反馈，系统结构如图1. 1所示。
图1. 1 系统结构图系统的闭环传递函数为式(1-1)（程序如附件I ）。
2()5()()36C s s R s s s φ==++ (1-1) 取采样周期为0.1，在MATLAB 中采用零阶保持器离散化得到脉冲传递函数如式(1-2)（程序如附件II ）。
1112(0.0226+0.0204)()1 1.68920.7408z z G z z z----=-+ (1-2) 根据极点分布如图1. 2，判断系统稳定（判别程序如附件III ）。
对式(1-2)进行差分处理，则可以得到式(1-3)。
() 1.6892(1)0.7408(2)0.0226(1)0.0204(2)y k y k y k u k u k =---+-+- (1-3)在实际的工作过程中，存在一定的误差或者环境会造成参数有所偏差，会导致实际测量到的输出有噪声干扰，因此输出要考虑噪声对过程的干扰。
差分方程为式(1-4)。
() 1.6892(1)0.7408(2)0.0226(1)0.0204(2)+()y k y k y k u k u k w k =---+-+- (1-4) 由式(1-4)可知，(1)[(1),(2),(1),(2)]T k y k y k u k u k ϕ-=----，[-1.6892,0.7408,0.0226,0.0204]T θ=，最小二乘法的自回归方程如式(1-5)。
()(1)()T y k k w k ϕθ=-+ (1-5)其中()y k 为过程输出，()k ϕ 为观测数据向量，θ为回归参数向量，()w k 为统计噪声或误差。
图1. 2 零极点分布图1.2 辨识准则函数及辨识方法由于成批型最小二乘(Batch Least-square , BLS )法不仅计算量大，占用内存多，而且不能很好适用于在线辨识。
第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 序列，于是可得
系统辨识第5章 线性动态模型参数辨识 最小二乘法
度函数
，则称uS(uk()为) “持续激励”信号。
● 定义4 一个具有谱密度 Fn (为z 1的) 平f1z稳1 信f2号z 2u(k)称fn为z nn 阶
“持续激励”Fn信(e号j )，2 S若u (对) 一0 切形如 Fn (e j ) 0
的滤波器，关系式
，意味着
。
● 定理2 设输入信号u(kR)u是(0)平稳R随u (1机) 信号，Ru (如n 果1)相关函数矩阵
式中
zL H L nL
nzHLLL[[zn(h(hh11TT)T),((,(zL12n())()22)),,,,znz(((LzLzL)(()]10]))1)
z(1 na ) z(2 na )
z(L na )
u(0) u(1)
u(L 1)
u(1 nb )
u(2
nb
)
u(L nb )
5.2 最小二乘法的基本概念
● 两种算法形式
① 批处理算法：利用一批观测数据，一次计算或经反复迭代，
以获得模型参数的估计值。
②
递推算法：在上次模型参数估计值
ˆ
(k
1)的基础上，根据当
前获得的数据提出修正，进而获得本次模型参数估计值ˆ (k )，
广泛采用的递推算法形式为
(k ) (k 1) K (k )h(k d )~z (k )
z(k ) h (k ) n(k )
式中z(k)为模型输出变量，h(k)为输入数据向量， 为模型参
数向量，n(k)为零均值随机噪声。为了求此模型的参数估计值， 可以利用上述最小二乘原理。根据观测到的已知数据序列
和{z(k)} ，{h极(k小)} 化下列准则函数
L
J ( ) [z(k ) h (k ) ]2
系统辨识之最小二乘
设
n
x (k) = n(k) + å ain(k - i) i =1
（4）
则式（3）可以写成
y(k) = -a1y(k -1) -!- an y(k - n) + b0u(k) + b1u(k -1) +!+
bnu(k - n) + x (k)
(5)
在测量 u(k) 时也有测量误差，系统内部也可能有噪声，应当考虑它们的影响。因此假定
ê ê
y(n
+
2)
ú ú
=
ê ê
- y(n + 1)
ê!úê
"
! - y(2)
u(n + 2)
u(2)
ú ú
´
"
"
"ú
ê ë
y(n
+
N
)úû
êë- y(n + N -1) ! - y(N ) u(n + N ) u(N )úû
éa1 ù
ê ê
!
êêêbb0n
ê! ê
ú ú ú ú ú ú ú
+
é x (n + 1) ù
个方程，即
y(n +1) = -a1 y(n) - a2 y(n -1) -!- an y(1) + b0u(n +1) + b1u(n) + !+ bnu(1) + x (n +1) y(n + 2) = -a1 y(n +1) - a2 y(n) -!- an y(2) +
b0u(n + 2) + b1u(n +1) +!+ bnu(2) + x (n + 2)
二阶系统的最小二乘辨识和广义预测控制仿真
二阶系统的最小二乘辨识和广义预测控制仿真二阶系统的最小二乘辨识和广义预测控制仿真1. 二阶系统模型1.1 系统模型及输入输出变量设传递函数为()G s =25+31s s +，()1H s =，且为负反馈，系统结构如图1. 1所示。
图1. 1 系统结构图系统的闭环传递函数为式(1-1)（程序如附件I ）。
2()5()()36C s s R s s s φ==++ (1-1) 取采样周期为0.1，在MATLAB 中采用零阶保持器离散化得到脉冲传递函数如式(1-2)（程序如附件II ）。
1112(0.0226+0.0204)()1 1.68920.7408z z G z z z----=-+ (1-2) 根据极点分布如图1. 2，判断系统稳定（判别程序如附件III ）。
对式(1-2)进行差分处理，则可以得到式(1-3)。
() 1.6892(1)0.7408(2)0.0226(1)0.0204(2)y k y k y k u k u k =---+-+- (1-3)在实际的工作过程中，存在一定的误差或者环境会造成参数有所偏差，会导致实际测量到的输出有噪声干扰，因此输出要考虑噪声对过程的干扰。
差分方程为式(1-4)。
() 1.6892(1)0.7408(2)0.0226(1)0.0204(2)+()y k y k y k u k u kw k =---+-+- (1-4)由式(1-4)可知，(1)[(1),(2),(1),(2)]Tk y k y k u k u k ?-=----，[-1.6892,0.7408,0.0226,0.0204]T θ=，最小二乘法的自回归方程如式(1-5)。
()(1)()T y k k w k ?θ=-+ (1-5)其中()y k 为过程输出，()k ? 为观测数据向量，θ为回归参数向量，()w k 为统计噪声或误差。
图1. 2 零极点分布图1.2 辨识准则函数及辨识方法由于成批型最小二乘(Batch Least-square , BLS )法不仅计算量大，占用内存多，而且不能很好适用于在线辨识。
系统辨识—最小二乘法
最小二乘法参数辨识1 引言系统辨识是根据系统的输入输出时间函数来确定描述系统行为的数学模型。
现代控制理论中的一个分支。
通过辨识建立数学模型的目的是估计表征系统行为的重要参数，建立一个能模仿真实系统行为的模型，用当前可测量的系统的输入和输出预测系统输出的未来演变，以及设计控制器。
对系统进行分析的主要问题是根据输入时间函数和系统的特性来确定输出信号。
对系统进行控制的主要问题是根据系统的特性设计控制输入，使输出满足预先规定的要求。
而系统辨识所研究的问题恰好是这些问题的逆问题。
通常，预先给定一个模型类μ={M}（即给定一类已知结构的模型），一类输入信号u和等价准则J=L(y，yM)(一般情况下，J是误差函数，是过程输出y和模型输出yM的一个泛函)；然后选择使误差函数J达到最小的模型，作为辨识所要求的结果。
系统辨识包括两个方面：结构辨识和参数估计。
在实际的辨识过程中，随着使用的方法不同，结构辨识和参数估计这两个方面并不是截然分开的，而是可以交织在一起进行的。
2 系统辨识的目的在提出和解决一个辨识问题时，明确最终使用模型的目的是至关重要的。
它对模型类（模型结构）、输入信号和等价准则的选择都有很大的影响。
通过辨识建立数学模型通常有四个目的。
①估计具有特定物理意义的参数有些表征系统行为的重要参数是难以直接测量的，例如在生理、生态、环境、经济等系统中就常有这种情况。
这就需要通过能观测到的输入输出数据，用辨识的方法去估计那些参数。
②仿真仿真的核心是要建立一个能模仿真实系统行为的模型。
用于系统分析的仿真模型要求能真实反映系统的特性。
用于系统设计的仿真，则强调设计参数能正确地符合它本身的物理意义。
③预测这是辨识的一个重要应用方面，其目的是用迄今为止系统的可测量的输入和输出去预测系统输出的未来的演变。
例如最常见的气象预报，洪水预报，其他如太阳黑子预报，市场价格的预测，河流污染物含量的预测等。
预测模型辨识的等价准则主要是使预测误差平方和最小。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
目录
1.最小二乘法简介
2.广义最小二乘法简介 3.用最小二乘法与广义最小二乘法辨识系统 4.总结
抓服务, 促教学, 兴科研!
2
电气工程与控制科学学院
1.最小二乘法介绍
抓服务, 促教学, 兴科研!
3
电气工程与控制科学学院
最小二乘估值算法
抓服务, 促教学, 兴科研!
4
电气工程与控制科学学院
抓服务, 促教学, 兴科研!
7
电气工程与控制科学学院
抓服务, 促教学, 兴科研!
8
电气工程与控制科学学院
最小二乘辨识结果：
最小二乘法
方差 0.0001 0.001 0.01 1.5954 1.5858 1.5227 0.7920 0.7887 0.7686 0.3991 0.3914 0.3773 0.5058 0.4959 0.4196
抓服务, 促教学, 兴科研!
13
电气工程与控制科学学院
College of Electrical Engineering and Control Science
智能电网 - 智能机器 - 智能建筑 - 智能制造
谢谢!
抓服务, 促教学, 兴科研!
0.1
真实值
0.8231
1.6
0.1593
0.8
0.4314
0.4
0.1495
0.5
抓服务, 促教学, 兴科研!
9
电气工程与控制Βιβλιοθήκη 学学院（2）广义最小二乘对系统的参数估计
①求出系统最小二乘估计
uk=[1.147 0.201 -0.787 -1.159 -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]; %输入 mm=0.01; %白噪声的方差 nn=normrnd(0,sqrt(mm),1,31); %产生个均值为0方差为mm的1x31的白噪声序列 yk(1)=0; %得出测量输出yk yk(2)=0; for i=3:31 yk(i)=-1.6*yk(i-1)-0.8*yk(i-2)+0.4*uk(i-1)+0.5*uk(i-2)+nn(i)+1.2*nn(i-1)+0.9*nn(i-2); end for i=1:29 %求矩阵������用A表示 A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)]; end
sita=inv(A‘*A)*A’*yk(3:31)‘;
抓服务, 促教学, 兴科研! 10
电气工程与控制科学学院
②迭代过程
for j=1:60 for i=3:31 yk(i)=-1.6*yk(i-1)-0.8*yk(i-2)+0.4*uk(i1)+0.5*uk(i-2)+nn(i)+1.6*nn(i-1)+0.8*nn(i-2); end for i=1:29 A(i,:)=[-yk(i+1) -yk(i) uk(i+1) uk(i)]; end sita=(inv(A'*A))*A'*yk(3:31)'; e(1)=yk(1); %计算残差 e(2)=yk(2)+sita(1)*yk(1)-sita(3)*uk(1); for i=3:31 e(i)=yk(i)+sita(1)*yk(i-1)+sita(2)*yk(i-2)sita(3)*uk(i-1)-sita(4)*uk(i-2); end
2.广义最小二乘法
抓服务, 促教学, 兴科研!
5
电气工程与控制科学学院
广义最小二乘算法具体步骤
抓服务, 促教学, 兴科研!
6
电气工程与控制科学学院
3.系统辨识实例：
设单输入-单输出系统的差分方程为： ������ ������ = −������1 ������ ������ − 1 − ������2 ������ ������ − 2 + ������1 ������ ������ − 1 + ������2 ������ ������ − 2 + ������ (������) ������ ������ = ������ ������ + ������1 ������ ������ − 1 + ������2 ������ ������ − 2 取真实值������������ = ������������ ������������ ������������ ������������ = 1.6 0.8 0.4 0.5 ，������1 = 1.2，������2 = 0.9。 输入数据如下所示： k u(k) k u(k) k 1 2 3 4 5 6 7 8 9 10 1.147 0.201 -0.787 -1.159 -1.052 0.866 1.152 1.573 0.626 0.433 11 12 13 14 15 16 17 18 19 20 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 -1.099 1.450 1.151 21 22 23 24 25 26 27 28 29 30
u(k) 0.485 1.633 0.043 1.326 1.706 -0.340 0.890 1.144 1.177 -0.390
用������的真实值，利用差分方程求出������ ������ 作为测量值，������ ������ 为均值为 0，方差为 0.0001， 0.001，0.01，0.1 的不相关随机序列。 (1) 用最小二乘法估计参数������������ = ������������ ������������ ������������ ������������ 。 (2) 用广义最小二乘法估计参数������������ = ������������ ������������ ������������ ������������ 。
电气工程与控制科学学院
College of Electrical Engineering and Control Science
智能电网 - 智能机器 - 智能建筑 - 智能制造
最小二乘与广义最小 二乘法系统辨识
组别：第2组
组员：史玄玄 薛竹韵 岳壮壮
抓服务, 促教学, 兴科研!
电气工程与控制科学学院
抓服务, 促教学, 兴科研!
11
电气工程与控制科学学院
广义最小二乘仿真结果
抓服务, 促教学, 兴科研!
12
电气工程与控制科学学院
4.总结
通过编程计算，发现在噪声方差比较小的情况下，各种方法所获得的估 值比较理想.
但随着噪声方差的增大，用最小二乘估值的偏差随之增大，但用广义最 小二乘法通过有限次的迭代运算能够更好地还原参数值.