系统辨识最小二乘参数估计matlab

最小二乘参数估计

摘要：

最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程()Z l T l l T

l ΦΦΦ-∧=1θ。

最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。

关键词：

最小二乘（Least-squares ），系统辨识（System Identification ） 目录：

1.目的 (2)

2.设备 (2)

3引言 (2)

课题背景 (2)

4数学模型的结构辨识 (3)

5 程序 (4)

M 序列子函数 (4)

主程序 (5)

6实验结果： (7)

7参考文献： (7)

1.目的

掌握系统辨识的理论、方法及应用

熟练Matlab 下最小二乘法编程

掌握M 序列产生方法

2.设备

PC 机1台（含Matlab 软件）

3引言

课题背景

最小二乘理论是有高斯（）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。

最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最好拟合的数学模型。递推最小二乘法是在最小二乘法得到的观测数据的基础上，用新引入的

数据对上一次估计的结果进行修正递推出下一个参数估计值，直到估计值达到满意的精确度为止。

4数学模型的结构辨识

根据汉格尔矩阵估计模型的阶次

设一个可观可控的SISO 过程的脉冲响应序列为{个g(1),g(2),……g(L)},可以通过汉格尔（Hankel ）矩阵的秩来确定系统的阶次。

令Hankel 阵为：

?????

???????-++-++++-++=)22()1()1()()2()1()1()1()(),(l k g k g l k g l k g k g k g l k g k g k g k l H ，其中l 决定),(k l H 阵地维数，k 可在1至()22+-l L 间任意选择。则有[]k n l n k l H rank ?≥=,,),(00。

如果0n l ≥（过程的真实阶次），那么Hankel 阵的秩等于0n 。因此可以利用Hankel 阵的奇异性来确定系统的阶次0n 。

根据残差平方和估计模型的阶次

SISO 过程的差分方程模型的输出残差为)(~k z ，数据长度L ，n H ?为n

?阶时的数据矩阵，n ??θ为n

?阶时的参数的估计量，n ?为模型阶次估计值，0n 为真实阶次，则残差平方和函数J ： )(~1)?()?(1~~1)?(1

2??????00k z L H z H z L z z L n J L n n k n n n T n n n n T n ∑++==--==θθ 残差平方和有这样的性质：当L 足够大时，随着n

?增加)?(n J 先是显著地下降，当n ?>0n 时，)?(n

J 值显著下降的现象就终止。这就是损失函数法来定阶的原理。

??

???--==-)?()?(1?)(?21ML n n T ML n n v n T n n T n ML H z H z L z H H H θ?σθ n L n

AIC v ?4?log )?(2+=σ 具体的定阶用法是：对不同阶次首先用极大似然法估计参数，然后计算似然函数值及)?(n

AIC 值，找到使min )?(=n AIC 的n ?作为0n 。

5 程序

%待辨识系统 z(k)=*z(k-1)*z(k-2)+*z(k-3)+u(k-1)+*u(k-2)*u(k-3)+v(k)/800%

clc

clear %清理工作间变量

L=300; % M 序列的周期

x1=1;x2=1;x3=1;x4=0;x5=1;x6=0; %四个移位积存器的输出初始值

for k=1:L; %开始循环，长度为L

u(k)=xor(x3,x4); %第一个移位积存器的输入是第3个与第4个移位积存器的输出的“或”

x6=x5;x5=x4;x4=x3;x3=x2;x2=x1;x1=u(k);

end %大循环结束，产生输入信号u

plot(u) %绘图M 序列

v=randn(300,1); %随机误差干扰

z=zeros(1,300);

for k=4:300

z(k)=*z(k-1)+*z(k-2)+*z(k-3)+*u(k-1)+*u(k-2)*u(k-3)+v(k)/400; %用理想输出值

end

H=zeros(300,6); %定义一个H“0”矩阵

for i=4:300

H(i,:)=[-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3)];%用循环产生H矩阵 z1(i,:)=[z(i)]; %用循环产生z矩阵

end

%计算参数%

c=inv(H'*H)*H'*z1%带入公式书上辨识出参数

%系统阶次辨识AIC算法%

bb=zeros(5,1);

n=1; %假设为1阶

for i=2:300

H1(i,:)=[-z(i-1) u(i-1)];

zz1(i,:)=[z(i)];

end

aa1=inv(H1'*H1)*H1'*zz1

bb(1)=(zz1-H1*aa1)'*(zz1-H1*aa1)/L;

AIC(1)=L*log(bb(1))+4*n;

n=2; %假设为2阶

for i=3:300

H2(i,:)=[-z(i-1) -z(i-2) u(i-1) u(i-2)];

zz2(i,:)=[z(i)];

end

aa2=inv(H2'*H2)*H2'*zz2

bb(2)=(zz2-H2*aa2)'*(zz2-H2*aa2)/L;

AIC(2)=L*log(bb(2))+4*n;

n=3; %假设为3阶

for i=4:300

H3(i,:)=[-z(i-1) -z(i-2) -z(i-3) u(i-1) u(i-2) u(i-3)];

zz3(i,:)=[z(i)];

end

aa3=inv(H3'*H3)*H3'*zz3

bb(3)=(zz3-H3*aa3)'*(zz3-H3*aa3)/L;

AIC(3)=L*log(bb(3))+4*n;

n=4; %假设为4阶

for i=5:300

H4(i,:)=[-z(i-1) -z(i-2) -z(i-3) -z(i-4) u(i-1) u(i-2) u(i-3) u(i-4)]; zz4(i,:)=[z(i)];

end

aa4=inv(H4'*H4)*H4'*zz4

bb(4)=(zz4-H4*aa4)'*(zz4-H4*aa4)/L;

AIC(4)=L*log(bb(4))+4*n;

n=5; %假设为5阶

H5(i,:)=[-z(i-1) -z(i-2) -z(i-3) -z(i-4) -z(i-5) u(i-1) u(i-2) u(i-3) u(i-4) u(i-5)];

zz5(i,:)=[z(i)];

end

aa5=inv(H5'*H5)*H5'*zz5

bb(5)=(zz5-H5*aa5)'*(zz5-H5*aa5)/L;

AIC(5)=L*log(bb(5))+4*n;

x=min(AIC)

for i=1:5

if(AIC(i)==x)

N=i %所辨识出的阶次N

end

end

plot(1:5,[AIC(1) AIC(2) AIC(3) AIC(4) AIC(5)])

6实验结果：

c =

a1 =

a2 =

a3 =

b1 =

b2 = b3 =

3