系统辨识作业1

《系统辨识》

作业 1

题目最小二乘参数估计仿真学院信息科学与工程学院

专业控制科学与工程

姓名王瑞荣

学号 030120694

教师顾幸生

2012年12月26日

一：一般最小二乘法 考虑仿真系统：

式中，)(k v 为方差为1的白噪声。

选用幅值为1的逆M 序列作为输入信号)(k u ，利用LS 算法进行参数估计，仿真结果如下表所示。

表1.1 批处理最小二乘法的参数估计结果

仿真程序：

%批处理最小二乘参数估计（LS ） clear all;

a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数

na=length(a)-1; nb=length(b)-1; %na 、nb 为A 、B 阶次

L=500; %数据长度

uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i) yk=zeros(na,1); %输出初值

x1=1; x2=1; x3=1; x4=0; S=1; %移位寄存器初值、方波初值 xi=randn(L,1); %白噪声序列

theta=[a(2:na+1);b]; %对象参数真值 for k=1:L

phi(k,:)=[-yk;uk(d:d+nb)]'; %此处phi(k,:)为行向量，便于组成phi 矩阵 y(k)=phi(k,:)*theta+xi(k); %采集输出数据

IM=xor(S,x4); %产生逆M 序列 if IM==0 u(k)=-1; else

u(k)=1; end

S=not(S); M=xor(x3,x4); %产生M 序列

%更新数据

)

()4(5.0)3()2(7.0)1(5.1)(k v k u k u k y k y k y +-+-=-+--

x4=x3; x3=x2; x2=x1; x1=M;

for i=d+nb:-1:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=na:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

end

thetae=inv(phi'*phi)*phi'*y' %计算参数估计值thetae

二：递推最小二乘算法 算法：

已知：n a 、n b 和d 。

Step1 设置初值)0(P 和)0(∧

θ，输入初始数据； Step2 采样当前输出)(k y 和输入)(k u ；

Step3 利用递推公式，计算K(k)、)(k ∧

θ和)(k P ； Step4 1+→k k ,返回Step2,继续循环。

考虑仿真系统：

式中，)(k v 为方差为1的白噪声。

取初值I P 6

10)0(=、0)0(=∧

θ。选择方差为1的白噪声作为输入信号)(k u ，

采用RLS 算法进行参数估计，仿真结果如图1.2所示。

51015

20253035

-2.5-2-1.5-1-0.500.5

11.52k

参数估计a 、b

图1.2 递推最小二乘法的参数估计结果

当k=35时，参数估计值为0.48780.9276,0.6389,-1.4399,1021====∧

∧

∧

∧

b b a a 。

)

()4(5.0)3()2(7.0)1(5.1)(k v k u k u k y k y k y +-+-=-+--

仿真程序：

%递推最小二乘参数估计（RLS）

clear all; close all;

a=[1 -1.5 0.7]'; b=[1 0.5]'; d=3; %对象参数

na=length(a)-1; nb=length(b)-1; %na、nb为A、B阶次L=35; %仿真长度

uk=zeros(d+nb,1); %输入初值：uk(i)表示u(k-i)

yk=zeros(na,1); %输出初值

u=randn(L,1); %输入采用白噪声序列

xi=sqrt(0.1)*randn(L,1); %白噪声序列

theta=[a(2:na+1);b]; %对象参数真值

thetae_1=zeros(na+nb+1,1); %thetae初值

P=10^6*eye(na+nb+1);

for k=1:L

phi=[-yk;uk(d:d+nb)]; %此处phi为列向量

y(k)=phi'*theta+xi(k); %采集输出数据

%递推最小二乘法

K=P*phi/(1+phi'*P*phi);

thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);

P=(eye(na+nb+1)-K*phi')*P;

%更新数据

thetae_1=thetae(:,k);

for i=d+nb:-1:2

uk(i)=uk(i-1);

end

uk(1)=u(k);

for i=na:-1:2

yk(i)=yk(i-1);

end

yk(1)=y(k);

end

plot([1:L],thetae); %line([1,L],[theta,theta]);

xlabel('k'); ylabel('参数估计a、b');

legend('a_1','a_2','b_0','b_1');

%axis([0 L -2 2]);