实验4 递推最小二乘法的实现

y (1) u (n 1) u (1) y ( n) y (n 1) y (2) u (n 2) u (2) y (n N 1) y ( N ) u (n N ) u ( N )
T T N N N N 1 1 T PN1 N 1 N 1 1
1
1
（11）
应用矩阵求逆引理，从求得 PN 1 与 PN 的递推关系式出发，经过一系列的推导，最终 可求得递推最小二乘法辨识公式：
T N 1 N K N 1 yN 1 N 1 N T K N 1 PN N 1 1 N 1 P N N 1 1
用户输入所需递推次数rec_num 将θN+1存入final_data θN+1与实际值得相对误差存入 final_percent
生成零均值白噪声ξ(k)
根据状态反馈控制律 计算增益系数矩阵K及控制U
满足递推次数？ 打印结果并绘制相关图形 是 结束
否
4
6．程序代码 function [] = Lab_4() %--------------------------------------------------------实验题目及初始化定义 disp(' ') disp('实验 4 递推最小二乘法的实现') disp(' ') disp('** 本实验对象系统的辨识模型为：') disp(' y(k)=-a1*y(k-1)-a2*y(k-2)+b0*u(k)+b1*u(k-1)+b2*u(k-2)+ξ(k).') disp(' 假定实际参数：a1=2, a2=1.3, b0=0.4，b1=0.88, b2=2.2. ') disp(' ') disp('** 本实验系统输入信号取为 u(k)=1.5*sin(0.2*k); 输出受到零均值白噪声ξ (k) ∈ [-0.1,0.1]污染.') disp(' ') a0=65539;M=2147483647;x=123456;b=10000; Y(1)=0;Y(2)=0;U(1)=1.5*sin(0.2*1);U(2)=1.5*sin(0.2*2);U(3)=1.5*sin(0.2*3); V=[];k_num=2;T_NT=[2;1.3;0.4;0.88;2.2]; c = input('** 为了进行推算，初值设定时，假定 P0=c*c*I，请输入充分大的常数 c = '); P_N0=c*c*eye(5); T_N0=zeros(5,1); rec_num = input(' 请输入所需递推步数 rec_num = '); disp(' ') rec_NUM = rec_num; final_data=zeros(rec_num,5); final_percent=zeros(rec_num,5); %--------------------------------------------------产生白噪声作为干扰信号ξ(k) for i=1:rec_num+2 x=mod(a0*x+b,M); V(i)=0.2*(x/M-0.5); end %------------------------------------------------------递推最小二乘法在线计算 while (rec_num~=0) HL(1)=-Y(k_num);HL(2)=-Y(k_num-1); HL(3)=U(k_num+1);HL(4)=U(k_num);HL(5)=U(k_num-1); Y(k_num+1)=HL*[2;1.3;0.4;0.88;2.2]+V(k_num+1); K_NN=P_N0*HL'*(inv(1+HL*P_N0*HL')); P_NN=P_N0-K_NN*HL*P_N0; T_NN=T_N0+K_NN*(Y(k_num+1)-HL*T_N0); final_data(k_num-1,:)=T_NN; T_ND1=abs(T_NN-T_NT); T_ND2=[T_ND1(1,1)/T_NT(1,1) T_ND1(2,1)/T_NT(2,1) T_ND1(3,1)/T_NT(3,1) T_ND1(4,1)/T_NT(4,1) T_ND1(5,1)/T_NT(5,1)]; final_percent(k_num-1,:)=T_ND2; P_N0=P_NN; T_N0=T_NN; k_num=k_num+1; K=[T_N0(2)-16 T_N0(1)-4]; U(k_num+1)=K*[Y(k_num);Y(k_num-1)]+[T_N0(3) T_N0(4) [1.5*sin(0.2*(k_num+1));1.5*sin(0.2*k_num);1.5*sin(0.2*(k_num-1))] T_N0(5)] *
5
rec_num=rec_num-1; end %-----------------------------------------------显示识别结果并与实际值进行比较 disp('** 运用递推最小二乘法对这一系统参数进行辨识，计算得到:') disp(' ') disp('相对于理论值[a1 a2 b0 b1 b2]=[2 1.3 0.4 0.88 2.2]，各次递推结果对应如 final_data 矩 阵所示：') final_data disp('实际值相对于理论值的误差如 final_percent 矩阵所示:') final_percent %-------------------------------------------------------------------相关绘图 k=0:1:rec_NUM-1; subplot(3,2,1) stairs(k,final_data(:,1)','k-');hold on line([0 rec_NUM-1],[T_NT(1,1) T_NT(1,1)]) axis([0 rec_NUM-1 -0.5 2.5]) title('a_1 收敛情况图形') subplot(3,2,3) stairs(k,final_data(:,2)','k-');hold on line([0 rec_NUM-1],[T_NT(2,1) T_NT(2,1)]) axis([0 rec_NUM-1 -0.2 1.5]) title('a_2 收敛情况图形') subplot(3,2,2) stairs(k,final_data(:,3)','k-');hold on line([0 rec_NUM-1],[T_NT(3,1) T_NT(3,1)]) axis([0 rec_NUM-1 0.1 1.4]) title('b_0 收敛情况图形') subplot(3,2,4) stairs(k,final_data(:,4)','k-');hold on line([0 rec_NUM-1],[T_NT(4,1) T_NT(4,1)]) axis([0 rec_NUM-1 0.6 1.3]) title('b_1 收敛情况图形') subplot(3,2,6) stairs(k,final_data(:,5)','k-');hold on line([0 rec_NUM-1],[T_NT(5,1) T_NT(5,1)]) axis([0 rec_NUM-1 -0.5 2.5]) title('b_2 收敛情况图形') end
（1）
其中 a1，a2…an，b0，b1，b2…bn 为待辨识的未知参数， ( k ) 是不相关随机序列。y 为系 统的输出， u 为系统的输入。 分别测出 n+N 个输出、 n+N 输入值 y(1)， y(2)…y(n+N)， u(1)， u(2)…u(n+N)，则可写出 N 个方程，具体写成矩阵形式，有
a1 y (1) u (n 1) u (1) (n 1) y (n 1) y (n) y (n 2) y (n 1) a y (2) u (n 2) u (2) n (n 2) （2） b0 y (n N ) y (n N 1) y ( N ) u (n N ) u ( N ) (n N ) bn
T N y ( N 1) u (n N 1) u ( N 1) 1 y ( n N )
将式（5）和式（8）合并，并写成分块矩阵形式，可得
YN N T y N 1 N 1 N 1
实验 4
递推最小二乘法的实现
实验报告
专业：
U
自动化
班级： 姓名： 日期： 年 月 日
1．实验题目：
U
递推最小二乘法的实现
2．实验目的 熟悉并掌握递推最小二乘法的算法原理。
3．实验主要原理 给定系统
y (k ) a1 y (k 1) a2 y (k 2) an y (k n) b0u (k ) b1u ( k 1) bnu (k n) (k )

（4）
为了实现实时控制，必须采用递推算法，这种辨识方法主要用于在线辨识。设已获 得的观测数据长度为 N，将式（3）中的 y 、 和 分别用 YN , N , N 来代替，即
YN N N
用 N 表示 的最小二乘估计，则
T N T N N N YN 1



（12） （13） （14）
T T PN 1 PN PN N 1 1 N N N 1 N 1 P N 1 P 1
为了进行递推计算，需要给出 PN 和 N 的初值 P0 和 0 。 推荐取值方法为： 假定 0 0, P0 c I ， c 是充分大的常数， I 为 (2n 1) (2n 1) 单
2



位矩阵，则经过若干次递推之后能得到较好的参数估计。
3
4．实验对象或参数 给定系统
y (k ) a1 y (k 1) a2 y ( k 2) b0u ( k ) b1u (k 1) b2u (k 2) (k ) （15）
即 n 2 。假设实际系统的参数为 a1 2 ， a2 1.3 ， b0 0.4 ， b1 0.88 ， b2 2.2 ， 但是不已知，即不可测。取 (k ) [0.1, 0.1] 的零均值白噪声。输入信号取为
u (k ) 1.5sin 0.2k
（16）
要求编制 MATLAB 程序，运用递推最小二乘法对这一系统的参数进行在线辨识，并 将辨识结果与实际参数进行对比。
5．程序框图
开始
程序初始化与赋初值
由先前输入输出生成 ΨN+1
用户输入初值设定常数c
由系统模型计算获得yN+1
设定初值θ0和P0
由递推最小二乘辨识公式得到 KN+1、PN+1、θN+1
（5）


（6）
令 PN N N
T


1
，则
N PN T N YN
如果再获得一组新的观测值 u (n N 1) 和 y (n N 1) ，则又增加一个方程
T yN 1 N 1 N 1

（7）
（8）
式中
yN 1 y (n N 1), N 1 (n N 1)
设
a1 y (n 1) (n 1) y (n 2) a , n , (n 2) ， y b0 (n N ) y (n N ) bn
T
1
N T N 1
T
YN y N 1
PN 1 T N YN N 1 y N 1
式中
YN y N 1
（10）
Βιβλιοθήκη Baidu
T N N PN 1 T T N 1 N 1
于是，类似地可得到新的参数估值
（9）
2
T N N N 1 T T N 1 N 1 N PN 1 T N 1
则式（2）可写为
y
（3）
1
式中：y 为 N 维输出向量； 为 N 维噪声向量； 为 2n+1 维参数向量； 为 N×(2n+1) 测量矩阵。为了尽量减小噪声 对 估值的影响，应取 N＞2n+1，即方程数目大于未知数 数目。
的最小二乘估计为
(T ) 1 T y