递推最小二乘估计及模型阶次辨识
递推最小二乘估计及模型阶次辨识
实验二 递推最小二乘估计(RLS)及模型阶次辨识(F-Test ）1 实验方案设计1.1 生成输入数据和噪声用M 序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。
生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。
1.2 过程仿真辨识模型的形式取)()()()()(11k e k u z B k z z A +=--，为方便起见，取n n n b a == 即nn n n zb a b z b z B z a a a z a z A ------++++=++++=...1)(...1)(22112211用M 序列作为辨识的输入信号。
1.3 递推遗忘因子法数据长度L 取534，初值⎪⎪⎪⎩⎪⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡==1000010000100001)0(001.0)0(P θ 1.4 计算损失函数、噪声标准差损失函数⎥⎦⎤⎢⎣⎡+---+-=μθμττ)()1()()]1(ˆ)()([)1()(2k h k P k h k k h k z k J k J噪声标准差θλdim )(ˆ-=L L J1.6 F-Test 定阶法计算模型阶次统计量t)22,2(~222)1()1()()1,(----++-=+n L F n L n J n J n J n n t其中，)(∙J 为相应阶次下的损失函数值，L 为所用的数据长度，n 为模型的估计阶次。
若a t n n t >+)1,(，拒绝00:n n H >，若a t n n t <+)1,(，接受00:n n H >，其中αt 为风险水平α下的阀值。
这时模型的阶次估计值可取1+n 。
1.6 计算噪信比和性能指标噪信比22ye σση= 参数估计平方相对偏差i i i ni i i θθθθθδˆ~,~1221-=⎪⎪⎭⎫ ⎝⎛=∑= 参数估计平方根偏差ii i n i ini iθθθθθδˆ~,)()~(2122122-==∑∑== 2 编程说明M 序列中，M 序列循环周期取15124=-=p N ，时钟节拍t ∆=1Sec ，幅度1=a ，特征多项式为1)(56⊕⊕=s s s F 。
第六章 最小二乘类参数辨识方法
《系统辨识基础》第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以上这种构成参数向量和数据向量的思想就是所谓的增广原理，它是增广最小二乘法的根本。
递推最小二乘估计
当C=I 时， [A+BD]-1 = A-1 -A-1B[I +DA-1B]-1DA-1
理论。 • 高斯仅用1小时就算出了谷神星的
轨道形状，并进行了预测
•1794年，高斯提出了最小二乘的思想。
1：引言
最小二乘法(Least Square)Gauss 1795年提出 在预测卫星和彗星运动的轨道时，需要处理由望远镜获得的观测数 据，以估计天体运动的六个参数。
Gauss在《天体运动理论》一书中写道：“未知量的最大概值是这 样的数值，它使各实测值与计算值之差的平方乘以度量其精度的数 值后，所得的和值达到最小。” ——著名的最小二乘思想 在系统辨识中，LS已成功应用于系统参数估计。 在自校正控制中，LS是应用最广泛的算法之一。
2：原理
2：原理
2：原理
2：原理
2：原理
3：特点
3：特点
（1）：无需存储全部数据，取得一组观测数据便 可估计一次参数，而且都能在一个采样周期中完 成，所需计算量小，占用的存储空间小。
（2）：具有一定的实时处理能力
谢谢
递推最小二乘估计（RLS）
董博南
一：最小二乘法回顾
二：递推最小二乘估计
一：最小二乘法回顾
1：引言 2：原理
3：特点
1：引言
• 1801年初，天文学家皮亚齐发现了谷神星。
•1801年末，天文爱好者奥博斯，在高斯预 言的时间里，再次发现谷神星。 •1802年又成功地预测了智神星的轨道。
有色噪声干扰输出误差系统的偏差补偿递推最小二乘辨识方法
y (t) = ϕ T (t)θ + ϕ T e (t)θ + e(t)
极小化输出误差准则函数
t
R 和向量 p 中 0 表示适当维数的零矩阵, 则利用式 (2), 由式 (4) 可得 1 −1 ˆ P (t)[θ LS (t) − θ ] = t→∞ t t 1 1 lim [ ϕ (i)ϕ T e (i)]θ + lim t→∞ t t→∞ t i=1 lim θ + p) − (Rθ
A(z ) = 1 + a1 z −1 + a2 z −2 + · · · + ana z −na B (z ) = b1 z −1 + b2 z −2 + · · · + bnb z −nb
不妨设 t ≤ 0 时, u(t) = 0, y (t) = 0, e(t) = 0. 假设
收稿日期 2006-9-12 收修改稿日期 2007-2-1 Received September 12, 2006; in revised form February 1, 2007 国家自然科学基金 (60674092, 60574051), 江苏省高技术研究 (工业) (BG2006010) 资助 Supported by National Natural Science Foundation of China (60674092, 60574051), High-tech R&D Program of Jiangsu (Industry) (BG2006010) 1. 江南大学通信与控制工程学院 无锡 214122 1. School of Communication and Control Engineering, Jiangnan University, Wuxi 214122 DOI: 10.1360/aas-007-1053
系统辩识作业题
系统辨识大作业
一.设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.比较所用三种方法的优劣及有效性；
五.给出由正弦输入求取系统开环频率响应特性曲线的辨识方法。
六.提出一种自己创造的辨识新方法，并用所给数据进行辨识验证。
注：闭卷考试时提交大作业报告。
使用最小二乘法法进行系统辨识的两种方法
递推最小二乘法辨识与仿真现在有如下的辨识仿真对象：图中， )(k v 是服从N )1,0(分布的不相关随机噪声。
且)(1-zG )()(11--=z A z B ，)(1-z N )()(11--=zC zD ， （1）⎪⎪⎩⎪⎪⎨⎧=+==+-=--------1)(5.00.1)()(7.05.11)(121112111z D z zz B z C z z a z A选择上图所示的辨识模型。
仿真对象选择如下的模型结构：)()2()1()2()1()(2121k w k u b k u b k y a k y a k y +-+-=-+-+可得系统模型为：)()2(5.0)1()2(7.0)1(5.1)(k w k u k u k y k y k y +-+-=-+-- 递推最小二乘法的推导公式如下：)]1(ˆ)()()[()1(ˆ)(ˆ--+-=k k k z k k k θh K θθτ )1()]()([)]()()1([)(11--=+-=--k k k k k k k τP h K I h h PP τ1]1)()1()()[()1()(-+--=k k k k k k h P h h P K τ相关程序如下：% exp053 %%递推最小二乘法程序%clear%清理工作间变量L=55;% M序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值for i=1:L;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，if y(i)>0.5,u(i)=-0.03;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”else u(i)=0.03;%当M序列的值为"0"时,辨识的输入信号取“0.03”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号uw=normrnd(0, sqrt(0.1), 1, 55);%加入白噪声figure(1);%第1个图形，伪随机序列stem(u),grid on%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零for k=3:55;%循环变量从3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k);%给出理想的辨识输出采样信号endc0=[0.001 0.001 0.001 0.001,0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(5,5);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=[c0,zeros(5,54)];%被辨识参数矩阵的初始值及大小e=zeros(5,55);%相对误差的初始值及大小for k=3:55; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2),w(k)]'; 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<=E break;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c;%显示被辨识参数e;%显示辨识结果的收敛情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);d1=c(5,:);ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);figure(2);%第2个图形i=1:55;%横坐标从1到55plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':',i,b1,'k') %画出a1，a2，b1，b2的各次辨识结果title('系统辨识结果')%图形标题如果系统中的参数为：⎪⎪⎩⎪⎪⎨⎧+-=+==+-=----------211211121112.01)(5.00.1)()(7.05.11)(z z z D z zz B z C z z a z A 那么系统模型机构为：)2(2.0)1()()2(5.0)1()2(7.0)1(5.1)(-+--+-+-=-+--k w k w k w k u k u k y k y k y 相关程序如下： % exp054.m % %递推最小二乘法编程%clear%清理工作间变量 L=55;% M 序列的周期y1=1;y2=1;y3=1;y4=0;%四个移位寄存器的输出初始值 for i=1:L;%开始循环，长度为Lx1=xor(y3,y4);%第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或” x2=y1;%第二个移位积存器的输入是第3个移位积存器的输出 x3=y2;%第三个移位积存器的输入是第2个移位积存器的输出 x4=y3;%第四个移位积存器的输入是第3个移位积存器的输出 y(i)=y4;%取出第四个移位积存器幅值为"0"和"1"的输出信号，if y(i)>0.5,u(i)=-0.03;%如果M序列的值为"1"时,辨识的输入信号取“-0.03”else u(i)=0.03;%当M序列的值为"0"时,辨识的输入信号取“0.03”end%小循环结束y1=x1;y2=x2;y3=x3;y4=x4;%为下一次的输入信号做准备end%大循环结束，产生输入信号uw=normrnd(0, sqrt(0.1), 1, 55);figure(1);%第1个图形，伪随机序列stem(u),grid on%以径的形式显示出输入信号并给图形加上网格z(2)=0;z(1)=0;%取z的前两个初始值为零for k=3:55;%循环变量从3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k)-w(k-1)+0.2*w(k-2);%给出理想的辨识输出采样信号endc0=[0.001 0.001 0.001 0.001,0.001]';%直接给出被辨识参数的初始值,即一个充分小的实向量p0=10^6*eye(5,5);%直接给出初始状态P0，即一个充分大的实数单位矩阵E=0.000000005;%相对误差E=0.000000005c=[c0,zeros(5,54)];%被辨识参数矩阵的初始值及大小e=zeros(5,55);%相对误差的初始值及大小for k=3:55; %开始求Kh1=[-z(k-1),-z(k-2),u(k-1),u(k-2),w(k),w(k-1),w(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<=E break;%若参数收敛满足要求，终止计算end%小循环结束end%大循环结束c;%显示被辨识参数e;%显示辨识结果的收敛情况%分离参数a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);d1=c(5,:);d2=c(6,:);d3=c(7,:); figure(2);%第2个图形i=1:55;%横坐标从1到55plot(i,a1,'r',i,a2,':',i,b1,'g',i,b2,':',i,d1,'k',i,d2,':',i,d3,'*') %画出a1，a2，b1，b2的各次辨识结果title('系统辨识结果')%图形标题。
渐消记忆递推最小二乘算法
是相互独立，服从N（0，）分布的随机变量。 1
M序列产生方法
X1
移位脉冲 XOR
X2
X3
X4
输出
1111 0111 0010 1001 1010 1101
0011 1100 1110
0001 0110 1111
14
1000 0100 1011 0101
111100010011010
M序列的性质
（1）由n级移位寄存器产生的M序列是确定的周期性序 列，它的周期长度为N=2n-1。 （2）n级移位寄存器中必须避免全部为“0”的状态，M 序列一个周期内状态“0”出现的次数比状态“1” 少1。 （3）M序列中，状态“0”或“1”连续出现的段称为游程。 游程中“0”或“1”的个数称为游程长度。
y0(t)+yw(t)
白噪声
Xw(t)
K g(τ)
延迟τ
Xw(t-τ)
乘法器
积分器
具有正常输入时的系统辨识原理图
23
用M序列辨识系统脉冲响应
（1） 计算方法 M序列的循环周期为N，移位脉冲的周期为t
a 2 , Rx ( ) a 2 , N
回顾
0, N ,2 N , 0, N ,2 N ,
1 N 1 RyM (k ) M ( j k )y( j ) N j 0
工程上：
a 2 t N 1 ˆ ( j) c g N j 0
25
c RyM ( NP 1)
（2）通过脉冲响应求传递函数
G(s)为系统的传递函数
G(s) C (s) M (s) R( s) N ( s)
阶跃响应法
输入 x(t) (a)
x(t)
递推最小二乘辨识
实验4 递推最小二乘法的实现实验报告哈尔滨工业大学航天学院控制科学与工程系专业：自动化班级：1040101姓名：日期：2013 年10 月23 日1．实验题目： 递推最小二乘法的实现 2．实验目的：(n a y k -()(n b u k n ξ+-+,,n a ，0,n b b 为待辨识的未知参数，(k ξ为系统的输出，为系统的输入。
分别测出n N +个输出、(3),(),(1),(2),()y y n N u u u n N ++，则可写出矩阵形式，有0()(1)(1)(1)((1)(2)(2)(2)((1)()()()(n n y n y u n u n a y n y u n u n b y n N y N u n N u N n b ξξξ⎢⎥--+⎤⎡⎤⎡⎢⎥⎢⎥⎢⎢⎥-+-+⎢⎥⎢+⎢⎥⎢⎥⎢⎢⎥⎢⎥⎢⎢⎥-+--+⎣⎦⎣⎢⎥⎢⎥⎣⎦（2）101)(1),,()()n n a n a b y n N n N b ξξθξξ⎡⎤⎢⎥++⎤⎡⎤⎢⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥++⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦，()(1)(1)(1)(1)(2)(2)(2)1)()()()y n y u n u n y u n u N y N u n N u N -+⎤⎥+-+⎥⎥⎥+--+⎦y θξ=Φ+ （维输出向量；ξ为N 维噪声向量；θ为21n +维参数向量；Φ测量矩阵。
为了尽量减小噪声ξ对θ估值的影响，应取2N n >+()T θ=ΦΦ为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨设已获得的观测数据长度为N ，将式N N Y θ=ΦN θ表示θ的最小二乘估计，则(TN Nθ=ΦΦ(T N N P =ΦΦT N N N P θ=Φ如果再获得一组新的观测值(u n N +11TN ψθ++=+(1),n N ++](1)(1)(1)y N u n N u N -++++）合并，并写成分块矩阵形式，可得T 111N N N N N Y y ξθξψ+++⎡⎤Φ⎡⎤⎤⎢⎥⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦（9）于是，类似地可得到新的参数估值1T T1T T T 1111N NN N N N N N N Y y θψψψ-+++++⎧⎫⎡⎤⎡⎤⎡⎤ΦΦΦ⎡⎪⎪⎢⎥⎢⎥⎢⎥⎪⎪⎢⎥=⎨⎬⎢⎥⎢⎥⎢⎥⎢⎥⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎪⎪⎩⎭ )T1T 1111N N N N N T N N N N Y P y Y y ψψ+++++⎡⎤Φ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦Φ+)T T 111T 11N NN N N N N ψψψψ++-++⎪⎢⎥⎢⎥⎪⎬⎢⎥⎢⎥⎪⎢⎥⎢⎥⎣⎦⎣⎦⎪⎭+ )1T 11N N ψ-++ 与P 的递推关系式出发，导，最终可求得递推最小二乘法辨识公式：()T1111N N N N N N K y θθψθ++++=+- ()1T 11111N N N N N N K P P ψψψ-++++=+ ()1T T 11111N N N N N N N N P P P P ψψψψ-++++=-+ 为了进行递推计算，需要给出N P 和N θ的初值0P 和0θ。
第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 )
广义递推最小二乘辨识
广义递推最小二乘辨识一、实验目的1 通过实验掌握广义最小二乘辨识算法；2 运用MATLAB编程，掌握算法实现方法。
二、实验原理广义最小二乘法的基本思想是基于对数据先进行一次滤波预处理，然后利用普通最小二乘法对滤波后的数据进行辨识。
如果滤波模型选择得合适，对数据进行了较好的白色化处理，那么直接利用普通最小二乘法就能获得无偏一致估计。
广义最小二乘法所用的滤波模型实际上就是一种动态模型，在整个迭代过程中不断靠偏差信息来调整这个滤波模型，使它逐渐逼近于一个较好的滤波模型，以便对数据进行较好的白色化处理，使模型参数估计称为无偏一致估计。
理论上说，广义最小二乘法所用的动态模型经过几次迭代调整后，便可对数据进行较好的白化处理，但是，当过程的输出噪信比比较大或模型参数比较多时，这种数据白色化处理的可靠性就会下降。
此时，准则函数可能出现多个局部收敛点，因而辨识结果可能使准则函数收敛于局部极小点上而不是全局极小点上。
这样，最终的辨识结果往往也会是有偏的。
其收敛速度比较慢，需要经过多次迭代计算，才能得到较准确的参数估计值。
一般情况下，经过多次迭代后，估计值便会收敛到稳态值。
但在某些情况下（如噪声比较低时）存在局部极小值，估计值不一定收敛到准则函数的全局极小值上。
为了防止参数估计值收敛到局部极小值，最好选定初值接近最优解，一般可以用最小二乘法的批处理估计值作为初值。
如果系统是时变的，或为了克服数据饱和现象，可以在两次RLS算法中分别引进遗忘因子。
三、实验内容<1> 数据获取：实验数据按照表9-1，为二阶线性离散系统的输入输出数据<2> 数据处理：为了提高辨识精度，实验者必须对原始数据进行剔除坏数据、零均值化、工频滤波等处理。
实验进行了白化滤波处理。
<3> 辨识算法：利用处理过的数据（取适当的数据长度），选择某种辨识方法（如RLS递推最小二乘法、RELS、RIV或RML等参数估计算法及F-检验或AIC定阶法），估计出模型参数和阶次，同时分析辨识结果。
基于递推最小二乘估计的动力学模型参数识别
基 于 递 推 最 小 二 乘 估 计 的 动 力 学 模 型 参 数 识 别
李 光 ， 肖国伟 ， 志 浩 艾
（ 南 工 业 大 学 机 械工 程 学 院 ， 南 株 洲 湖 湖 ４２０ ） １ ０ ８
摘
要 ： 力 学建 模 过 程 中 ， 型 参 数 的 不 确 定 性 和 参 数 的 未 知 变化 使 得 建 模 精 度 大 为 降 低 。 通 过 对 动 模
非 线 性 最 小二 乘 估 计 理 论 的 研 究 ， 出动 力 学模 型参 数 的 最 小二 乘 估 计 ， 设 计 出相 应 的估 计 器 ， 于 现 场 提 并 基 数 据 给 予模 型较 为 准 确 的 定 论 。 应 用 于 曲柄 滑 块 机 构 中滑 块 与 滑 动 面 间 的 摩 擦 系数 的 估 计 ， 过 较 少 的 迭 经 代 次数 获 得 满 意 的 估 计 效 果 。 为 动 力 学 模 型 参 数 识 别 ， 立更 精 确 模 型 提 供 一 个 有 效 的途 径 。 建
ＬＩＧｕ ｎ ＸＩ ａ ｇ， ＡＯ ｏ ｗ ，ＡＩＺｈ — ｏ Ｇｕ — ｉ ｉｈａ
（ ｎ ｎＵｎｖ ｒ ｔ ｏ ｅ ｈ ｏ ｇ ， ｈ ｚ ｏ ４ ２ ０ ， ｈｎ ） Ｈｕ ａ ｉ ｓ ｙ ｆ ｃ ｎ ｌ ｙ Ｚ ｕｈ ｕ １ ０ ８ Ｃ ｉ ｅ ｉ Ｔ ｏ ａ
Ａｂ ｔ ａ ｔ Ｉ ｓ ｈ ｒ ｌ Ｏ ｏ ｔｉ ｃ ｕ ａ ｅ ｍｏ ｅ ｎ ｔ ｅ ｃ ｕ ｓ ｆｄ ｎ ｍｉ ｍｏ ｅｉｇ ｆ ｒ ｔ ｅ ｕ ｃ ｒ ａｎ ｙ ａ ｄ ｔ ｅ ｕ ｋ ｏ ｓｒ ｃ ：ｔ ｉ ａ ｄｙ ｔ ｂ ａ ｎ ａ ａ ｃ ｒ ｔ ｄ ｌｉ ｈ ｏ ｒｅ ｏ ｙ ａ ｃ ｄ ｌ ｏ ｈ ｎ ｅ ｔｉ ｔ ｎ ｈ ｎ ｎ ｗｎ ｎ ｃ ａ ｇ ｓｏ ｄ ｌ ａ ａ ｔ ｒ ．Ｉ ｈ ｓｐ ｐ ｒ ｈ ｅ ｕ ｓｖ ａ ｔ ｑ ａ ｅ ｅ ｔｍａ ｉ ｎ ｆ ｒ ｄ ｎ ｍｉ ｍｏ ｅ ａ ｅ ｎ ｐ ｏ ｏ ｅ ｎ ｈ ｈ ｎ ｅ ｆｍｏ ｅ ｒｍｅ ｅ ｓ ｎ ｔ ｉ ａ ｅ ，ｔ ｅ ｒ ｃ ｒｉ ｅ ｌｓ —ｓ ｕ ｒ ｓｉ ｔ ｏ ｙ ａ ｃ ｐ ｏ ｄ ｌｈ ｓ ｂ ｅ ｒ ｐ ｓ ｄ ａ ｄ ｔ ｅ ｃ ｒｅ ｐ ｎ ｉｇ ｅ ｔ ｔ ｒｈ ｓａｓ ｅ ｎ ｄ ｓｇ ｅ ｆｅ ｔ ｄｎ ｎ ｔ ｅ ｎ ｎ— ｌ ｅ ｒｌａ ｔ ｑ ａ ｅ ｓ ｉ ｔ ｎｔ ｅ ｒ ．Ｔｈ ｎ ａ ｍｏ ｅａ Ｃ — ｏ ｒ ｓ ｏ ｄ ｎ ｓ ｉ ｏ ａ ｌ ｂ ｅ ｅ ｉ ｎ ｄ ａ ｔｒｓｕ ｉ ｇｏ ｈ ｏ ｍａ ｏ ｉ ａ ｅ ｓ ｕ ｒｓｅ ｔ ｎ ｓ ｍａ ｉ ｈ ｏ ｙ ｏ ｅ ｒ Ｃ Ｈ ｒ ｔ ｏ ｃｕ ｉｎ ｆ ｒｍｏ ｅ ｐ ｒ ｍｅ ｅ ａ ｅ ａ ａｎ ｄ ｂ ｓ ｇｏ ａ ｅ ｃ ｎ ｌｓｏ ｏ ｄ ｌ ａ ａ ｔ ｒｃ ｎ ｂ ｂ ｉ ｅ ａ ｉ ｎ—ｌ ｅｄ ｔ ．Ｆ ｎ ｌ ｎ ｉ ａａ ｉ ａｌ ｔｈ ｓｂ ｅｓ ｃ ｅｓ ａｙａ ｐ ｉｄ ｔ ｔｍａ ｅｔ ｅｆｉ— ｎ ｙ，ｉ ａ ｅ ｕ ｃ ｓｆ ｌ ｐ ｌ Ｏｅ ｉ ｔ ｈ ｒ ｅ ｓ ｃ ｔ ｎ ｃ ｅｆ ｉｎ ｅ ｗｅ ｎ ｔ ｅ ｓ ｄ ｒａ ｄ ｔ ｅｓｉｉ ｇ ｓ ｒａ ｅｏ ｌ ｅ —ｃ ａ ｋ ｍｅ ｈ ｎ ｓ ａ ｄ ａ ｓ ｔ ｆ ｃｏ ｙ ｖ ｌｅ ｈ ｓ ｂ ｅ ａｎ ｄ ａｔ ｒ ｉ ｏ ｆｉ ｅ ｔｂ ｔ ｅ ｈ ｌ ｅ ｎ ｈ ｌ ｎ ｕ ￣ ｃ ｆ ｉ ｒ ｒ ｎ ｃ ａ ｉ ｏ ｃ ｉ ｄ ｆ ｓｄ ｍ ｎ ａ ｉ ａ ｔ ｒ ａｕ ａ ｅ ｎ ｇ ｉ ｅ ｆｅ ｓ ｌｓ ｕ ｒ ｕ ｔｒ ｔ ｎ ｔ ｓ Ａｃ ｕ ｌ ｔ ｌ ｐ ｏ ｉｅ ａ ｆｅ ｔ ｅｗａ ｅ ｔ ｙ ｄ ｎ ｍｉ ｍｏ ｅ ａ ａ ｔｒ ｎ ｕ ｌ ｒ ｃ ｕ ｅ ｓ ｎ ｍｅｏ ｓｉ ａ ｉ ｉ ｅ ｏ ｍｅ ． ｔ ａｌ ｉ ｗｉ ｒ ｖ ｄ ｎ ｅ ｆｃ ｉ ｙ ｔ ｉ ｎ ｉ ｙ ａ ｃ ｙ ｌ ｖ ｏｄ ｆ ｄ ｌ ｒ ｍｅ ｅｓａ ｄ ｂ ｉ ａ ｍｏ ｅａ ｃ — ｐ ｄ ｒｔ ｄ １ ａ ｅ ｍｏ ｅ ． Ｋｅ ｒ ｓ ｒ ｃ ｒｉ ｅ ｌ ｓ —ｓ ｕ ｒ ｓｉ ｔ ｎ；ｐ ｒ ｍｅｅ ｅ ｏ ｎ ｔ ｎ；ｄ ｎ ｍｉ ｍｏ ｅ ｙ ｗｏ ｄ ：ｅ ｕ ｓｖ ｅ ｔ ｑ ａ ｅ ｅ ｔ ａ ｍａ ｉ ｏ ａ ａ ｔｒｒ ｃ ｇ ｉ ｏ ｉ ｙａ ｃ ｄｌ
递推最小二乘估计及模型阶次辨识
《糸统辨枳基础》第19讲要点卖验二盪推最小二乘估计(RLS)及模型阶次辨枳(F・Test)、卖验目的①通过卖验，•拿握透推最小二乘参数辨识方法②通过卖验，•拿握F・Test模型阶次辨识方冻、卖验内彖1、仿真模型卖验所用的仿真模型如下：框图表示z(k) 1.5z(k 1) 0.7z(k 2) u(k 1) 0.5u(k) v(k)其中u(k)和z(k)分别为模型的输入和输出变量；v(k)为零均值、方差为1、服从正态分布的自噪步；为噪＞的标准差(卖验肘，可取0.0、0.1 . 0.5、1.0)；输入变量u( k)采用M序列，其特征多项式取F(s) s4 s 1,幅度取1.0o2、辨识模型辨识模型的形式取A(z 1)z(k) B(z 1)u(k) e(k)为方便起见，取n a nbn ,即A(z) 1 aiz 1 a2z 2a n zB(z) biz 1 b2z 2b n Z n根据仿真模型生成的数据u(k),k 1, ,L和z(k),k 1, ,L ,辨帜模型的参＜ai,a2, ,a n和bi,b2, , bn ；并确定模型阶次n,同肘估计出模型誤案e(k)的方差(应近似等于模型噪＞v(k)的方差，即为2)和模型的静态增益Ko3、辨识算法①采用适推遗忘因子冻：(k) (k1)K(k) z(k) h (k) (k1)K(k) P(k1)h(k)h(k)P(k1)h(k)|i1P(k)dK(k)h (k) P(k1)p其中，遗忘因子0 1 (具体值根据情况自己确定)；数据长度L可取100、300、500；初始(0) P(0) a2 I②赖夾函数的透推计算：J(k)|iJ(k1)z(k)h(k) (k1)2h (k) P(k1)h(k) |44、 F-Test 走阶法统计量tJ(n 1)2t(n,n J J(n)J(n"L2n2~F(2,L2n 2)其中，J （）为相应阶次下的赖夾函数值，L 为所用的数据长度，n 为模型的估计阶次。
最小二乘参数辨识方法及原理PPT学习课件
Gauss（1777-1855）
m
使 w(k ) | z(k ) y最(k小) |2 k 1
1、问题的提出
1795年，高斯提出的最小二乘的基本原理是
未知量的最可能值是使各项实际观测值和计 算值之间差的平方乘以其精确度的数值以后的和 为最小。
z(k) y(k) v(k)
Gauss（1777-1855）
Z m H m Vm
2.2 一般最小二乘法原理及算法
最小二乘的思想就是寻找一个 的估计值ˆ ，使得各次测量 的 Z i (i 1, m) 与由估计ˆ 确定的量测估计 Zˆi Hiˆ 之差的平方
和最小，即
J (ˆ) (Zm Hmˆ)T (Zm Hmˆ) min
J
ˆ
2H
T m
(Z
m
H mˆ)
i1
2.2 一般最小二乘法原理及算法
u(k )
y(k )
G(k )
v(k ) z(k )
图 3.4 SISO 系统的“灰箱”结构
G(z)
y(z) u(z)
b1z 1 b2 z 2 1 a1z 1 a2 z 2
bn z n an z
n
n
n
y(k ) ai y(k i) biu(k i)
i 1
0 0
J a J b
a b bˆ
aˆ
N
2 (Ri a bti )
i 1 N
2 (Ri a bti )ti
i 1
0 0
Naˆ
N
bˆ
N i 1
ti
N
N i 1
Ri
N
aˆ
i 1
ti
bˆ
t
2 i
i 1
递推最小二乘法相关辨识方法(97-2003)
定义 式中
k 1 T T P k H k H k h i h i i 1 k 1 P 1 k 1 H T H h i hT i k 1 k 1 i 1
hT 1 T h 2 Hk T h k
把P(k)写成
利用矩阵反演公式 可得
P k P
1
1
k 1 h k h k
T
1
1
A CCT A1 A1C I CT A1C CT A1
1
P k P k 1 P k 1 h k hT k P k 1 hT k P k 1 h k 1 P k 1 h k hT k I T P k 1 h k P k 1 h k 1
z k a1z k 1 ana z k na bu k 1 bnb u k nb e k 1
令k=1,2,3,…,L，可构成线性方程组 zL H L eL 式中
z 1 e 1 z 2 e 2 zL , eL z L e L u 0 u 1 nb z 0 z 1 na z 1 z 2 na u 1 u 2 nb HL z L 1 z L na u L 1 u L nb
hT 1 T h 2 H k 1 T h k 1
式中，h(i)是一个列向量，也就是HL 的第i行的倒置， P(k)是一个方阵，它的维数取决于未知参数的个数。
递推最小二乘辨识讲解
（2）每增加一次观测量，都必须重新计算φ，（）T -1。
（3）如果出现φ列相关，既不满秩的情况， T 为病态矩阵，则不能得到最 小二乘估计值。 *递推最小二乘参数辨识，就是当被辨识的系统在运行时，每取得一次新的观
测数据后，就在前一次估计结果的基础上，利用新引入的观测数据对前次 估计的结果，根据递推算法进行修正，从而递推地得出新的参数估计值。 这样，随着新的观测数据的逐次引入，一次接一次的进行参数计算，直到 参数估计值达到满意的精确程度为止。
P(k -1)(k -1) (k -1)
P
(k
)
I
-
1
(k
- 1) P
(k
-1)
(k
-1)
P
(k
-1)
其中的计算顺序为先计算P(k),然后再计算 θˆ(k) .
(7) (8)
有时,为计算方便并便于理解,上述RLS估计算法又可表示 为
常规最小二乘辨识的递推算法
主要内容
1.思想及原理 2.实例与MATLAB仿真 3.应用
1.递推最小二乘法的思想及原理
1.1递推最小二乘法的引入 *最小二乘法的缺陷 （1）数据量越多，系统参数估计的精度就越高，为了获得满意的辨识结果，
矩阵的阶数 T 常常取得相当大。这样矩阵求逆的计算量很大，存储量
1.2递推算法的思想 * 递推辨识算法的思想可以概括成
新的参数估计值=旧的参数估计值+修正项 即新的递推参数估计值是在旧的递推估计值的基础上而成，
这就是递推的概念. 递推算法不仅可减少计算量和存储量,而且能实现在线
实时辨识. * 递推算法是依时间顺序,每获得一次新的观测数据就修
正一次参数估计值,随着时间的推移,便能获得满意的辨 识结果. RLS法即为成批型LS算法的递推化，即将成批型LS算法
系统辨识最小二乘模型阶次大全
A:仿真对象为：z(k)-1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k)B：最小二乘一次算法和递推算法程序：clearclcNp=400;u=[1 0 0 1 1 0 1 0 1 1];for i=10:(Np-1)u(i)=xor(u(i-4),u(i-9));end%²úÉúMÐòÁÐfor j=10:(Np-1)if mod(j,2)==0y(j)=1;elsey(j)=0;endu(i)=xor(u(i),y(j));end%²úÉúÄæMÐòÁÐv=randn(3,400);z=[1,398]';for i=3:400z(i)=1.5*z(i-1)-0.7*z(i-2)+1*u(i-1)+0.5*u(i-2)+v(i);endH=zeros(400,4);for i=3:400H(i,1)=-z(i-1);H(i,2)=-z(i-2);H(i,3)=u(i-1);H(i,4)=u(i-2);endtheta1=inv(H'*H)*H'*zi=1:400;subplot(1,2,1)plot(i,theta1(1,:),'r',i,theta1(2,:),'k',i,theta1(3,:),'g',i,theta1(4 ,:),'b')title('惠一龙，一次最小二乘算法')P=10^6*eye(4);theta2=0.001*eye(4,1);I=eye(4);h=H';for i=3:400K=P*h(:,i)*[H(i,:)*P*h(:,i)+1]^(-1);P=[I-K*H(i,:)]*P;theta2=theta2+K*[z(i)-H(i,:)*theta2];end disp('') theta2 i=1:400; subplot(1,2,2)plot(i,theta2(1,:),'r',i,theta2(2,:),'k',i,theta2(3,:),'g',i,theta2(4,:),'b')title('惠一龙，递推算法')C ：结果示意图：0100200300400-1.5-0.5惠一龙，一次最小二乘算法0100200300400惠一龙，递推算法D ：结果分析：辨识结果误差较小。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2
A( z 1 ) a 0 a1 z 1 a 2 z 2 a n z n , a 0 0 B( z 1 ) b0 b1 z 1 b2 z 2 bn z n ,
则有
2 y
1 1 B( z ) B( z 1 ) dz 1 dz G ( z ) G ( z ) ＝ 2j l z 2j l A( z ) A( z 1 ) z
6、计算性能指标 ① 参数估计平方相对偏差
~ i 1 i 1 i
2n
~ , i i i
2
② 参数估计平方根偏差
2
③ 静态增益估计相对偏差
i 1
2n
~2
i
பைடு நூலகம்
i 1
2n
2 i
~ , i i i
~ K K K, K K K
4
平 下的阀值。这时模型的阶次估计值可取 n 1 。 注：F 分布值表（风险水平 ＝ 0。 05 ） 阀值 t 自由度 2 100 300 500 5、噪信比计算 ● 噪信比定义 噪信比
2 e 2 y
自由度 1
2
3.09 3.03 3.01
u(k)
e(k) G(z-1) y(k) z(k)
2 其中， e2 为噪声方差， y 为过程输出方差。
2 ● 过程输出方差 y 的计算
2 y
1 dz G( z )G( z 1 ) l 2j z
其中，积分围线 l 是 z 平面内沿逆时间方向的单位圆圆周。若定义 B( z 1 ) 1 G( z ) A( z 1 ) 式中
1 ＝ a0
i 0
n
b
i
i 2
i a0
k 1 k 1 k 1 k 1 k a0 ai ak 1 a k 1 i a i k 1 a0 k 1 k 1 k a b b k 1 a k 1 k 1 k 1 i b 0 i k 1 i a0 k n 1, n 2, ,1, 0 ; i 0, 1, , k n n a a , b bi i i i
1
2 z (k ) h (k ) (k 1) J (k ) μ J (k 1) h (k ) P (k 1)h(k ) μ
③ 噪声标准差的估计
④ 模型静态增益估计
＝
J ( L) L dim
i K＝ i 1 n 1 ai i 1
A( z 1 ) z(k ) B( z 1 )u(k ) e(k )
为方便起见，取 n a n b n ，即
A( z ) 1 a1 z 1 a 2 z 2 a n z n B( z ) b1 z 1 b2 z 2 bn z n 根 据 仿 真 模 型 生 成 的 数 据 u(k ), k 1,, L 和 z(k ), k 1,, L ， 辨 识 模 型 的 参 数 a1 , a 2 ,, a n 和b1 , b2 ,, bn ；并确定模型阶次 n ，同时估计出模型误差 e(k ) 的方差（应近似等
《系统辨识基础》第 19 讲要点
实验二 递推最小二乘估计(RLS)及模型阶次辨识(F-Test） 一、实验目的 ① 通过实验，掌握递推最小二乘参数辨识方法 ② 通过实验，掌握 F-Test 模型阶次辨识方法 二、实验内容 1、仿真模型 实验所用的仿真模型如下： 框图表示 v(k)
1 1.5z 1
1
0.7 z 2
e(k) y(k) + + z(k)
u(k) 模型表示
z 1 0.5z 2 1 1.5z 1 0.7 z 2
z(k ) 1.5z(k 1) 0.7 z(k 2) u(k 1) 0.5u(k ) v(k )
其中 u(k)和 z(k)分别为模型的输入和输出变量；v(k)为零均值、方差为 1、服从正态分布的 白噪声； 为噪声的标准差（实验时，可取 0.0、0.1、0.5、1.0）；输入变量 u(k)采用 M 序 列，其特征多项式取 F ( s) s 4 s 1 ,幅度取 1.0。 2、辨识模型 辨识模型的形式取
i K＝ i 1 n , K＝ i 1 n 1 ai 1 ai i 1 i 1
~
bi
n
b
n
三、程序流程（供参考） 启动 定维
输入数据 u(534)，输出数据 z(534)，M 序列 M(5) 参数估计向量 THETA(8)， 数据向量 h(8)， 协方差矩阵 P(8,8) 损失函数 J(4)，噪声标准差 LAMBDA
b
n
4、F-Test 定阶法 统计量 t
J (n) J (n 1) L 2n 2 ~ F (2, L 2n 2) J (n 1) 2 其中， J () 为相应阶次下的损失函数值， L 为所用的数据长度， n 为模型的估计阶次。 t (n, n 1)
t ,接受 H 0 : n 1 n0 ，其中 t 为风险水 t ,拒绝 H 0 : n n 0 ,若 t (n, N 1〕 若 t (n, n 1〕
于模型噪声 v(k ) 的方差，即为 2 ）和模型的静态增益 K。 3、辨识算法 ① 采用递推遗忘因子法： ( k ) ( k 1) K ( k ) z ( k ) h ( k ) ( k 1) 1 K ( k ) P ( k 1) h( k ) h ( k ) P ( k 1) h( k ) μ 1 P ( k ) I K ( k ) h ( k ) P ( k 1) μ 其中，遗忘因子 0 1 (具体值根据情况自已确定)；数据长度 L 可取 100、300、500；初始 (0) 值 。 2 P (0) a I ② 损失函数的递推计算：
3
赋初值
生成 M 序列参数：a=1, P＝4, M(0),…,M(5)不能全为 0 生成白噪声参数：M＝32768, A=179, x0=11
人机对话
噪声标准差：Lambda；数据长度：L；遗忘因子： Model Order: from Nbeg to Nend
过程仿真
生成 M 序列；生成白噪声；生成过程输入和输出数据
设定模型阶次，从 Nbeg 到 Nend
模型参数估计：RLS 算法 损失函数计算
模型阶次辨识 计算性能指标 打印实验结果及性能指标 四、实验步骤 (1) 掌握最小二乘递推算法和 F-Test 模型阶次辨识的基本原理。 (2) 设计实验方案。 (3) 编制实验程序。 (4) 调试程序，研究实验问题，记录数据。 (5) 分析实验结果，完成实验报告。 五、实验报告 实验报告包括实验方案设计、编程说明、源程序清单、数据记录、结果分析、误差计算、 数据列表、曲线打印、实验体会等。 注意：为实验二安排使用计算机时间如下： 时间：2001 年 12 月 17－19 日晚上 7:00-9:45 地点：中央主楼 5 楼计算机实验室