系统辨识之最小二乘法

方法一、最小二乘一次性算法：

首先对最小二乘法的一次性辨识算法做简要介绍如下： 过程的黑箱模型如图所示：

其中u(k)和z(k)分别是过程的输入输出，G(z)描述输入输出关系的模型，成为过程模型。

过程的输入输出关系可以描述成以下最小二乘格式：

1z(k)hT(k)n(k) （1）

其中z(k)为系统输出，是待辨识的参数，h(k)是观测数据向量，n(k) 是均值为0的随机噪声。

利用数据序列{z（k）}和{h（k）}极小化下列准则函数：

J()[z(k)hT(k)]2 （2）

k1L使J最小的的估计值，成为最小二乘估计值。

具体的对于时不变SISO动态过程的数学模型为

^A(z1)z(k)B(z1)u(k)n(k) （3）

11A(z)B(Z)的系数。 应该利用过程的输入、输出数据确定和

对于求解的估计值，一般对模型的阶次型写成最小二乘格式

^na，

nb已定，且

nanb；其次将（3）模

z(k)hT(k)n(k) （4）

式中

Th(k)[z(k1),,z(kna),u(k1),,u(knb)]T[a,a,,a,b,b,,b]12n12nab （5）

k1,2,,L

因此结合式（4）（5）可以得到一个线性方程组

ZLHLnL （6）

其中

zL[z(1),z(2),z(L)]TTnL[n(1),n(2),n(L)] （7）

对此可以分析得出，HL矩阵的行数为Lmax(na,nb)，列数nanb。

在过程的输入为2n阶次，噪声为方差为1，均值为0的随机序列，数据长度L(nanb)的情况下，取加权矩阵L为正定的单位矩阵I，可以得出：

TT(HLHL)1HLzL （8） ^ 其次，利用在Matlab中编写M文件，实现上述算法。

此次算法的实现，采用6阶M序作为过程黑箱的输入；噪声采用方差为1，均值为0的随机数序列；黑箱模型假设为：y(k)-1.5y(k-1)+0.7y(k-2)=2u(k-1)+0.5u(k-2)，则系统输出为Z(k)-1.5Z(k-1)+0.7Z(k-2)=2U(k-1)+0.5U(k-2)+n（k）；模型的阶次na2,nb2；数据长度取L=200。

程序清单如下见附录：最小二乘一次性算法Matlab程序 运行结果如下：

图1 最小二乘一次性算法参数真值与估计值

其中re为真值，ans为估计值

结果发现辨识出的参数与真值之间存在细微误差，这是由于系统噪声以及数据长度L的限制引起的，最小二乘辨识法是一种无偏估计方法。

^方法二、最小二乘递推算法：

最小二乘一次性算法计算量大，并且浪费存储空间，不利于在线应用，由此引出最小

二乘参数估计的递推算法，其基本思想是按观测次序一步一修正，即：

(k1)K(k)[z(k)h(k)(k1)] K(k)P(k1)h(k)[h'(k)P(k1)h(k)(1k)]1

P(k)[IK(k)h'(k)]P(k1)

这里选取估计参数的初值为全0矩阵，P矩阵设置参数为100的对角阵，增益阵K初值设置为10。程序清单如下见附录：最小二乘递推算法Matlab程序。

程序运行结果如下：

^^'^

图2 最小二乘递推算法参数真值与估计值

图3 最小二乘递推算法待估参数过渡过程

附录：

1. 最小二乘一次性算法Matlab程序

%最小二乘一次性算法

%Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+2*u(k-1)+0.5u(k-2)+v(k) % na=2,nb=2 clear clc

tic %计时开始

L=200; %数据长度

%生成M序列 6阶 a1=1; a2=0; a3=1; a4=0; a5=1;

a6=0; M=zeros(L,1); for i=1:1:L M(i)=a1; a1=a2; a2=a3; a3=a4; a4=a5; a5=a6;

a6=xor(M(i),a1); end

%噪声

v=randn(1,L);

% na=2, 设置测量输出z(1)到z(na)的值为0 z=[]; z(1)=0; z(2)=0;

%产生输出序列z 从z(na+1)开始,即z(3)开始 for i=3:L

z(i)=1.5*z(i-1)-0.7*z(i-2)+2*M(i-1)+0.5*M(i-2)+v(i-2); end

%构造HL矩阵

na=2;nb=2; %模型阶次，并找出最大者赋给max if na>nb

max=na; else

max=nb; end

%HL矩阵的的维数 hang=L-max; lie=na+nb;

%构建HL矩阵 用z和M填充HL矩阵与na,nb的最大数max有关 HL=zeros(hang,lie);

a=na; for j=1:1:na

for i=1:1:hang

HL(i,j)=-z(i+a-1); end a=a-1; end b=nb;

for j=na+1:1:lie for i=1:1:hang

HL(i,j)=M(i+b-1); end b=b-1; end

%构建ZL矩阵 ZL=zeros(hang,1); for i=1:1:hang; ZL(i)=z(i+max); end

%参数估计值 th=zeros(na+nb,1);

th=inv((HL'*HL))*HL'*ZL; th;

%比较估计误差 re为真值，th为估计值 re=[-1.5,0.7,2,0.5] th'

toc %计时结束

2. 最小二乘递推算法Matlab程序 %最小二乘递推算法

%Z(k)=1.5*Z(k-1)-0.7*Z(k-2)+2*u(k-1)+0.5u(k-2)+v(k) % na=2,nb=2

clear clc

tic %计时开始

L=200; %数据长度

%生成M序列 6阶 a1=1; a2=0; a3=1; a4=0; a5=1;

a6=0; M=zeros(L,1); for i=1:1:L M(i)=a1; a1=a2; a2=a3; a3=a4; a4=a5; a5=a6;

a6=xor(M(i),a1); end

%噪声

v=randn(1,L);

% na=2, 设置测量输出z(1)到z(na)的值为0 z=[]; z(1)=0; z(2)=0;

%产生输出序列z 从z(na+1)开始,即z(3)开始 for i=3:L

z(i)=1.5*z(i-1)-0.7*z(i-2)+2*M(i-1)+0.5*M(i-2)+v(i-2); end

na=2;nb=2; %模型阶次，并找出最大者赋给max if na>nb

max=na; else

max=nb; end

%HL矩阵的的维数 hang=L-max; lie=na+nb;

%构建HL矩阵 用z和M填充HL矩阵与na,nb的最大数max有关 HL=zeros(hang,lie);

a=na; for j=1:1:na

for i=1:1:hang

HL(i,j)=-z(i+a-1); end a=a-1; end b=nb;

for j=na+1:1:lie for i=1:1:hang

HL(i,j)=M(i+b-1); end b=b-1; end

%P是(na+nb)*(na+nb)维，K是(na+nb)*1维，h'是1*(na+nb)维 %递推求解

P=100*eye(na+nb);

Pstore=zeros(na+nb,hang);

Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; th=zeros(na+nb,hang);

% 预估参数初值 矩阵K初值

th(:,1)=[0.001;0.001;0.001;0.001]; K=[10;10;10;10];

for i=1:hang h=HL(i,:)';

K=P*h*inv(h'*P*h+1);

th(:,i+1)=th(:,i)+K*(z(i+na)-h'*th(:,i)); P=(eye(na+nb)-K*h')*P;

Pstore(:,i+1)=[P(1,1),P(2,2),P(3,3),P(4,4)]; end

re=[-1.5,0.7,2,0.5] th(:,hang+1)'

toc %计时结束

i=1:hang+1; figure(1)

plot(i,th(1,:),i,th(2,:),i,th(3,:),i,th(4,:))

title('待估参数过渡矩阵')

figure(2)

plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)) title('估计方差变化过程')