系统辨识

最小二乘法的系统辨识

摘要：在研究一个控制系统过程中，建立系统的模型十分必要。因此，系统辨识在控制系统的研究中起到了至关重要的作用。本文主要介绍了系统辨识的最小二乘方法，最小二乘法的一次完成过程进行了推导，最小二乘法的一次完成的缺陷在于对于有色噪声并没有很好的辨识效果。其中系统辨识在工程中的应用非常广泛,系统辨识的方法有很多种,最小二乘法是一种应用极其广泛的系统辨识方法，阐述了动态系统模型的建立及其最小二乘法在系统辨识中的应用,并通过实例分析最小二乘法应用于直流调速系统的系统辨识。

关键词：系统辨识、最小二乘法

一、系统辨识的定义

系统辨识、状态估计和控制理论是现代控制理论三个相互渗透的环节。1962年，L．A．zadeh给出“辨识”的定义为：系统辨识是在对输入和输出观测的基础上，在指定的一类系统中，确定一个与被识别的系统等价的系统。[1]最先提出了系统辨识的定义。

随着科技的发展，数学建模对科学研究及指导及生产都有非常重要的意义。给一个系统建立数学模型是一个比较复杂的工作，其中关键的一个环节是系统辨识。系统辨识就是研究如何利用系统的输入、输出信号建立系统的数学模型。[7]系统数学模型是系统输入、输出及其相关变量间的数学关系式，它描述系统输入、输出及相关变量之间相互影响、变化的规律性。换句话说，系统辨识就是从系统的运算和实验数据建立系统的模型(模型结构和参数)。系统辨识的三要素：数据、模型类和准则。系统辨识的基本原理：在输入输出的基础上，从一类系统中确定一个与所测系统等价的系统。[2]

二、最小二乘法的引出

最小二乘法是1795年高斯在预测星体运行轨道最先提出的，它奠定了最小二乘估计理论的基础．到了20世纪60年代瑞典学者Austron把这个方法用于动态系统的辨识中，在这种辨识方法中，首先给出模型类型，在该类型下确定系统模型的最优参数。

我们可以将所研究的对象按照对其了解的程度分成白箱、灰箱和黑箱。于其内部结构、机制只了解一部分，对于其内部运行规律并不十分清楚，这样的研究对象通常称之为“灰箱”；如果我们对于研究对象的内部结构、内部机制及运行规律均一无所知的话，则把这样的研究对象称之为“黑箱”。研究灰箱和黑箱时，将研究的对象看作是一个系统，通过建立该系统的模型，对模型参数进行辨识来确定该系统的运行规律。对于动态系统辨识的方法有很多，但其中应用最广泛，辨识效果良好的就是最小二乘辨识方法，研究最小二乘法在系统辨识中的应用具有现实的、广泛的意义。[4]

应用最小二乘法对系统模型参数进行辨识的方法有离线辨识和在线辨识两种离线辨识是在采集到系统模型所需全部输入输出数据后，用最小二乘法对数据进行集中处理，从而获得模型参数的估计值；而在线辨识是一种在系统运行过程中进行的递推辨识方法，所应用的数据是实时采集的系统输入输出数据，应用递推算法对参数估计值进行不断修正，以取得更为准确的参数估计值。[8]

三、最小二乘法的原理

3.1 最小二乘法一次完成推导[1]

本文中以一个SISO 系统为例说明最小二乘法的原理。 假设一个SISO 系统如下图所示：

图1 SISO 系统结构图

其离散传递函数为：

3.1

输入输出的关系为：

)()()()(1

k y k e z G k u =+?- 3.2

进一步，我们可以得到：

)()()()()(1

1

k e z B k u z A k y +?=?-- 3.3

其中，扰动量)(k e 为均值为0，不相关的白噪声。 将式3.3写成差分方程的形式：

)

()()2()1()()2()1()(2121k e n k u b k u b k u b n k y a k y a k y a k y n n +-?+-+-+--?-----= 3.4

令T n k u k u k u n k y k y k y k ])()2()1()()2()1([)(-?----?----=?

][2121n n b b b a a a ??=

θ

n

n n

n z a z a z a z b z b z b z A z B z G ---------+?++++?++==221122111111)

()()(

则式3.4可以写为：

)()()(k e k k y T +=θ? 3.5

将上述式子扩展到N 个输入、输出观测值{)(),(k y k u }，k=1,2,…，N+n 。将其代入到式3.5中，写成矩阵的形式为：

e Y +Φ=θ 3.6

其中，

取泛函)(θJ 为

)()()()()(21

2θθθθΦ-Φ-=?=+=Φ-=∑∑=Y Y e e i n e Y J T T N

i

最小二乘法原理既是使)(θJ 最小，对其求极值得：

由此可得系统的最小二乘法估计值为：

Y T T T ΦΦΦ=)(θ

这样，我们就得到了系统的最小二乘估计值。以上推导的最小二乘法存在一些缺点，比如：预先取得的观测值越多，系统参数估计的精度越高，但使得矩阵ΦΦT 的阶数越大，矩阵求逆计算量也越大，所需的存储空间也会越大；每增加一次观测值，必须重新计算Φ和1)(-ΦΦT ； 若T Φ列相关,即不满秩,则ΦΦT 为病态矩阵,无法求得最小二乘估计值。因此可以对该方法进行改进，从而可以得到递推最小二乘法，广义最小二乘法等改进后的方法。这里不做具体的介绍。 3.2最小二乘法的缺陷[ 5]

最小二乘是一种最基本的辨识方法，最小二乘法可以用于线性系统，也可以用于非线性系统；可用于离线估计和在线估计。在随机情况下，利用最小二乘法时，并不要求观测数据提供其概率统计方法的信息，而其估计结果，却有相当好的统计特性。但它具有两方面的缺陷：一是当模型噪声是有色噪声时，最小二乘估计不是无偏、一致估计；二是随着数据的增长，将出现所谓的“数据饱和”现象。针对这两个问题，出现了相应的辨识算法，如遗忘因子法、限定记忆法、偏差补偿法、增广最小二乘、广义最小二乘、

T

T

N n e n e n e e N n y n y n y Y ])()2()1([])()2()1([+?++=+?++=?

?

???

??

??

???+--+-+-??+-?---?-?---=Φ)()

()()

1()()2()1()2()()1()1()()1()1()

(N u N n u N y N n y N n y u n u y n y n y u n u y n y n y 0)]()[(=Φ-Φ-??=??θθθ

θY Y J T

辅助变量法、二步法及多级最小二乘法等。

四、最小二乘法应用在直流调速系统的系统辨识举例

闭环系统是把反映输出转速的电压信号反馈到系统输入端，与给定电压比较，形成一个闭环。由于反馈的作用，系统可以自行调整转速，这种方式也称为反馈控制。

引入转速负反馈信号以后，放大器的输入信号是给定信号s U 和反馈信号fn U 之差，即fn s n U U U -=?。当电动机负载增加时，电枢电阻压降必然增大。若是开环系统，电动机转速只能下降，但是在闭环系统中，转速稍有降落，反馈电压fn U 随之下降，接着是

↑↑→↓→↑→↑→?n U U U d ct n α

上述调节过程是当负载增加（或降低）时，相应地整流电压d U 就提高（或减小），从而补偿掉电动机电枢回路电阻电压的增量，维持电动机反电势E （或转速n ）几乎不变。

转速反馈闭环调速系统是一种基本的反馈控制系统，它具有三个基本特征。一，只用比例放大器的反馈控制系统，其被调量仍是有静差的。二，反馈控制系统的作用是：抵抗扰动，服从给定。三，系统的精度依赖于给定和反馈检测精度。图1.1是用集成电路运算放大器作为电压放大器的转速负反馈闭环控制有静差直流调速系统。检测的反馈信号fn U 与转速n 成正比，n U fn α=，α又称为转速反馈系数。由该调速系统的工作原理可以确定系统的输入量为电压S U ，输出量 为电动机转速n

电压比较环节：

fn s n U U U -=?

放大器： n P c K U U ?=

测速发电机：

n U fn α=

——放大器的电压放大系数；0

f

P R K R =-

α——测速反馈系数，单位为Vmin/r; 额定励磁下直流电动机

E dt

dI L

RI d

d d ++=0U （主电路，假定电流连续） n C E

e Φ=（额定励磁下的感应电动势）

dt

dn

GD T T L e ?=-3752（牛顿动力学定律，忽略粘性摩擦）

d m

e I C T Φ=（额定励磁下的电磁转矩）

式中 L T —包括电机空载转矩在内的负载转矩，单位为Nm ；

2GD —电力拖动系统运动部分折算到电机轴上的飞轮力矩，单位为2Nm ；

e m C C π

30

=

——电动机额定励磁下的转矩电流比，单位为Nm/A ；

定义下列时间常数：

R

L

T l =

——电枢回路电磁时间常数，单位为s ； m

e m C C GD T 3752

=

——电力拖动系统机电时间常数，单位为； 得电压与电流间的传递函数

1

/1)()()(I s 0+=

-s T R

s E s U s d d 电流与电动势间的传递函数为

s

T R

s I s I s E m dl d =

-)()()(

额定励磁下直流电动机的动态结构图如下：

图4.2 额定励磁下直流电动机的动态结构图

如果是理想空载，则0)(=s I dl ，直流电动机电枢回路的传递函数

()201()()1

e d m l m C N s U s T T s T s Φ=++

对以上各环节表达式整理，得到以电压为S U 输入，转速n 为输出的传递函数：

(*)

其中 该系统为典型的二阶系统。 在不影响数学模型的结构，性质的基础上简化数学模型，我们令 这样就有P K K α=，需要辨识的参数就有3个m T ,l T ,K .我们取观测的数据3N =. 以直流电压()S U s 为输入，以电动机转速n 为输出的系统模型为

2

()()1S m l m N s K

U s T T s T s K

=+++ (**) 其中P K K α=.

采用双极性变换法将式（**）离散化，即将1

1

211z s T z ---=+代入式（**）得

21222212()(12)()42(1)2(1)842(1)m l m m l m l m N z T K z z U z T T T T K T z K T T T z T T TT K T -----++=??????-++++-++++??????

将上式化为差分方程得到：

()()12345(1)2()(1)(2)N K a N K a N K a U K a U K a K =----++-+- 其中212

2(1)842(1)m l m l m K T T T a T T TT K T

-+-=

+++ T=1s 为采样周期

222

42(1)42(1)m l m m l m T T T T K T a T T TT K T

-++=

+++

()2()()1e S m l m K C N s U s T T s T s K Φ=

+++P e K K C α=Φ

1

1e C =Φ

2

32

42(1)m l m KT a T T TT K T

=+++ 2

42

242(1)m l m KT a T T TT K T

=+++

2

52

42(1)m l m KT a T T TT K T =+++ 现取电机调速的传递函数为 则采用双极性变换法将其离散化得：

现θ∧

的真实值已知,取.

则仿真对象为

其中e(K)为均值为0，方差为1的不相关随机序列。

辨识模型选择

()()12345(1)2()(1)(2)()N K a N K a N K a U K a U K a K e K =----++-+-+

观测数据长度400L =，则经最小二乘递推程序仿真后得到

估计方差变化：

估计参数过程：

=++2

()5

()25

S N s U s s s ()()=----++-+-25555

(1)2()(1)(2)

13

13

13

13

13

N K N K N K U K U K K =

====123452551051313131313

a a a a a ()()255105

(1)2()(1)(2)()

1313131313N K N K N K U K U K K e K =----++-+-+

五、结论

本文针对于系统辨识问题的某些方法进行了简要的介绍。详细推导了最小二乘法一次完成的算法，并举例说明了最小二乘法在模型参数的辨识中能得到良好的效果。

程序清单

clear

clc

%==========400 个产生M序列作为输入==

x=[0 1 0 1 1 0 1 1 1]; %initial valu

n=403; %n为脉冲数目

M=[]; %存放M 序列

for i=1:n

temp=xor(x(4),x(9));

M(i)=x(9);

for j=9:-1:2

x(j)=x(j-1);

end

x(1)=temp;

end

%===========产生均值为0，方差为1 的高斯白噪声=========

v=randn(1,400);

%==============产生观测序列z=================

z=zeros(402,1);

z(1)=-1;

z(2)=0;

for i=3:402

z(i)=(-0.1538)*z(i-1)-0.3846*z(i-2)+0.3846*M(i)+0.7692*M(i-1)+0.3846*M( i-2)+v(i-2);

end

%递推求解

P=100*eye(5); %估计方差

Pstore=zeros(5,401);

Pstore(:,1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5)];

Theta=zeros(5,401); %参数的估计值，存放中间过程估值

Theta(:,1)=[3;3;3;3;3];

% K=zeros(4,400); %增益矩阵

K=[10;10;10;10;10];

for i=3:402

h=[-z(i-1);-z(i-2);M(i);M(i-1);M(i-2)];

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

Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2));

P=(eye(5)-K*h')*P;

Pstore(:,i-1)=[P(1,1),P(2,2),P(3,3),P(4,4),P(5,5)];

end

i=1:401;

figure(1)

plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)) title('待估参数过渡过程')

figure(2)

plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:),i,Pstore(5 ,:))

title('估计方差变化过程')