系统辨识与仿真实验大报告汇总

东南大学自动化学院

实验报告

课程名称：系统辨识

实验名称：系统辨识与仿真

院（系）：自动化专业：自动化

姓名：学号：

同组人员：实验时间：2013 年月日评定成绩：审阅教师：

实验一：

给定系统G(s)= 4102+++s s s ， 1.

用连续系统仿真方法，求零时刻的单位脉冲响应g(k)（采样步长为0.1s ，持续5s ）； 2.

用双线性变换求离散传递函数G(z)； 3.

编程产生7级逆重复M 序列u ，幅值：±1，仿真时间：0：30.5s ，仿真步长：0.1s ； 4.

以u 为输入，求G(z)的输出y ；（保存y,u ，作为实验二的数据） 5. 编程求自相关函数Ruu(k)和互相关函数Ruy(k)，k=1, (51)

(化成单位脉冲及其响应的形式)

6. 将g(k)、Ruu(k)、Ruy(k)在一个图中绘出，并比较g(k)和Ruy(k)的差别。

实验步骤及实验结果：

（1）建立simulink 仿真系统图如下：

设置脉冲周期为10s ，幅值为10，脉宽为1%，运行后观察示波器结果如下，并输出到文件：

（2）输入如下代码进行双线性变换：

>> n1=[0 1 10];

>> d1=[1 1 4];

>> f=tf(n1,d1)

84.392.724

.41

.02.03.022+--+z z z z

Transfer function:

s + 10

-----------

s^2 + s + 4

>> [n2,d2]=c2dm(n1,d1,0.1,'tustin');

>> printsys(n2,d2,'z')

num/den =

0.070755 z^2 + 0.04717 z - 0.023585

-----------------------------------

z^2 - 1.8679 z + 0.90566

双线性变化为)1()1(2+-=z T z s ，代入原传递函数求得G(z)= ，

将其分子分母同除以4.24即为程序结果。

（3）逆重复M 序列的生成代码：

clear;

t=[0:0.1:30.5];

lenth=size(t,2);

u(1:7)=1;

for i=8:lenth

u(i)=(mod(u(i-1)+u(i-7),2));

end

u=u*2-1;

a=1;

for i=1:lenth %与方波相乘

u(i)=u(i)*a;

a=-a;

end

uu(1,:)=t;

uu(2,:)=u;

save u11 uu

plot(u(1:140))

ylim([-2 2])

显示结果如下图所示：

（4）采用simulink仿真，建立仿真模型如下：

运行后示波器输出如下波形：

（5）输入如下代码：

load u11;

uy(1,:)=uu(2,:);

a=[0.3 0.2 -0.1];

b=[4.24 -7.92 3.84];

uy(2,:)=filter(a,b,uy(1,:));

length=size(uy,2);

t=[0.1:0.1:5.1];

lt=size(t,2);

for tao=1:lt

for i=1:2

mu(i,tao)=0;

for j=1:127

mu(i,tao)=mu(i,tao)+uy(i,j+tao)*uy(1,j+2); %加2是为了把脉冲画完整end;

mu(i,tao)=mu(i,tao)*10/(127); %乘10是为了化为单位脉冲end;

end;

plot(t,mu(1,:),'r-',t,mu(2,:),'b-')

显示结果如下图所示：

（6）运行结果如下图：

图中Ryu（t）和y（k）两者波形相近，但是前者对于后者有滞后。

实验二：

运行sj11.p ，按任务一的要求产生仿真数据y 、u ，并编程求出模型的参数估计；其中：

A(q -1)=1-1.5q -1+0.7q -2 B(q -1)=0.2q -2

B 1(q -1)=0.2q -2 B 2(q -1)=-0.5q -1 C(q -1)=1-0.5q -1

实验步骤及实验结果：

任务一：请仿真一个Ay=Bu+Aw....模型,并用递推随机逼近法——及AIC 准则法来辨识这个模型。

设广义误差e(k)是参数估计值θ的函数，参数辨识问题可通过极小化e(k)的方差来实现。即求参

数θ使下列准则函数最小：J(θ)=1/2E{e2(k)}。J(θ)的负梯度为： =E{-e(k) }。如果可求 解 =0，则可求得参数的估计。但当e(k)的分布未知时，实际上是不可求解的。 在计算数学中，求二次函数的极小值常采用迭代法。首先给出参数的一个估计值，以二次函数在该参数估计值处的负梯度为修正方向，选取适当的步长后，修正参数估计值，直到收敛。

仿此，我们有：θ(k+1)= θ(k)+ρ(k) 。 如果在求 时不求期望，则得到一个随机的迭代算法，称之为随机逼近法。 load u11;

>> u=uu(2,:);

a=[1 -1.5 0.7];

b=[0 0 0.2];

y=filter(b,a,u);

lenth=size(uu,2);

τ

θ??

??????-J τθ????????-J τθ????????-J τθ????????-J

a1=sqrt(cov(y)); %求得信号的标准差

r=normrnd(0,1,1,lenth);

>> a=[1];

>> b=[1];

x=filter(b,a,r);

b1=sqrt(cov(x));

y1=y+x*0.2*a1/b1; %噪信比设为0.1

plot(y1)

uyr(2,:)=u;

uyr(1,:)=y1;

>> save y6 uyr

>> load y6;

>> y=uyr(1,:)';

u=uyr(2,:)';

clear uyr;

% 识别输入维数和数据总数

[langth,r]=size(u);

% 设定模型结构

tao(1)=0;

m(1)=input('输入A(z)的阶次');

mn=m(1);

for i=1:r

i

m(i+1)=1+input('输入B(z)的阶次');

mn=mn+m(i+1); % 计算参数维数

tao(i+1)=input('输入B(z)的纯时延');

end

mm=max(m)+max(tao)+1; % 求得方程式的起始位置m(r+2)=input('输入C(z)的阶次');

mn=mn+m(r+2);

% 初始化：THETA(0) P(0) L(0) FI(0) PSI(0)

P=1000000*eye(mn);

THETAa=zeros(mn,1);

L=THETAa;

FI=L;

PSI=L;

%---------------------

cv=1; % 迭代标志(数据复用）

while cv==1

v=0; % 损失函数初值

lamda=0.95; %遗忘因子

THETA=THETAa;

THETAa=zeros(mn,1); %参数平均值保存

%构造初始FI

for i=1:m(1)

FI(i)=-y(mm-i);

end

mk=m(1);

for j=1:r

for i=1:m(j+1)

FI(mk+i)=u(mm-tao(j)-i+1,j);

end

mk=mk+m(j+1);

end

PSI=FI;

for k=mm:langth % 递推估计参数

EPS=y(k)-FI'*THETA; % 计算预报误差

AMY=stable(L,THETA,EPS,m);%判断A(z)的稳定性

THETA=THETA+L*EPS*AMY; % 计算新参数估计

if k>langth-15

THETAa=THETA+THETAa;

end

FII=FI;

for i=mk+1:mk+m(r+2)

FII(i)=0;

end

y1=FI'*THETA; %带噪声模型的输出估计

y2=FII'*THETA; %不带噪声模型的输出估计

EPS1=y(k)-y1; % d--白色噪声估计

EPS2=y(k)-y2; % e--有色噪声估计

%==================================================== %

yy(1)=-y(k); % a--方程误差法

%yy(1)=-y1; % b--输出误差法

%yy(1)=-y2; % c--输出误差法2

%=================================================

%yy(r+2)=-EPS2; %可有的组合：（c,e)，（b,e)

%

yy(r+2)=EPS1; %可有的组合：（a,d)增广最小二乘，（b,d) %yy(r+2)=EPS2; %可有的组合：（c,e)，（b,e)

%yy(r+2)=-EPS1; %可有的组合：（a,d)增广最小二乘，（b,d) %====================================================

for i=1:r

if k-tao(i+1)+1>0

yy(i+1)=u(k-tao(i+1)+1,i);

else

yy(i+1)=0;

end

end

ii=0; %% 更新FI PSI

for i=1:r+2

if m(i)>1

for j=2:m(i)

ij=m(i)+2-j;

FI(ii+ij)=FI(ii+ij-1);

% PSI(ii+ij)=PSI(ii+ij-1);

end

end

if m(i)>0,FI(ii+1)=yy(i);end

ii=ii+m(i);

end

[L,P,v]=UD(FI,L,P,v,mn,lamda,EPS); %% 计算增益L，修正P和V end

THETAa=THETAa/15

v

cv=input('0--停止，1--再迭代');

for i=1:mn

P(i,i)=P(i,i)*10;

end

end

输出结果为：

THETAa =

-1.4577

0.6005

0.1759

-1.2732

-0.0204

v =

59.7182

实验三：

以实验二中的数据y,u，编程求模型的结构和参数估计；

clear all; %最小二乘法for MISO

load y6;

tt(:,1)=uyr(1,:)'; %读入数据，并赋给变量tt。

tt(:,2)=uyr(2,:)'; %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uyr;

%plot(tt(:,1))

ll=size(tt); %得到数据维数

r=ll(2)-1;

nn=2;

et=0;

while(et==0)

clear m tao ff fa zz zta fb yy;

m(1)=input('输入A(z)的阶次'); %指定模型结构

tao(1)=0;

for i=1:r

i %给出输入编号

m(i+1)=input('输入B(z)的阶次')+1; %阶次加一，表示参数个数

tao(i+1)=input('输入B(z)的时延');

end

if m(1)

n=m(1)+max(tao)+1; %算出一个方程最多使用的数据

lll=ll(1)-n; %算出可列出的方程数

%lll=200;

in1=r+1; %构造观测数据矩阵ff

[zta,m,tao]=ls0(tt,m,tao,ll);

taomax=max(tao+m);

a=[1 zta(1:m(1))']; %输出仿真

b=[zeros(1,tao(2)) zta(m(1)+1:m(1)+m(2))'];

tf=filter(b,a,tt(:,2));

e=tt(:,1)-tf;

ee=0;

for k=1+taomax:taomax+190

ee=ee+e(k)*e(k);

end

ee=ee/190

AIC=2*size(zta,1)+190*log(ee) %AIC准则plot(e)

et=input('改变结构？YES=0,NO=1 ');

end;

输出结果为

zta =

-1.2233

0.4389

0.2047

m =

2 1

tao =

0 2

ee =

0.0696

AIC =

-500.2515

误差曲线如下：

实验四：

运行sj11.p，按任务二得到数据文件，并进行结构和参数估计。

任务二：未知系统辨识。您的数据文件名是:335.mat。其中，数组uu的第一行是输出，第二行是输入。

clear all;%逐步回归法

load 335;

ll=size(uu);

if ll(1)

%tt=uyr';

tt = uu';

else

%tt=ury;

tt = uu;

end

clear uu;

ll=size(tt);

in1=ll(2);

[zta,m,tao]=ls10(tt); %辨识高阶模型

%----------------------------------------------

a=[1 zta(1:m(1))']; %输出仿真

ttr=tt;

ttr(:,1)=0;

k=m(1);

for i=1:in1-1

b=[zta(k+1:k+m(i+1))'];

k=k+m(i+1);

ttr(:,1)=ttr(:,1)+filter(b,a,tt(:,i+1));

end

%---------------------------------------------

%---------------------------------------------

[zta,m,tao]=stepreg1(ttr); %逐步回归法

%计算输出误差

taomax=max(tao+m);

for k=1:taomax

y(k)=tt(k,1); %用原观测数据计算输出误差

end

for k=1+taomax:ll(1)

mm=0;

for i=1:in1

if i==1

for j=1:m(i)

% fb(mm+j)=-y(k-j);

fb(mm+j)=-tt(k-j,1);

end;

else

for j=1:m(i)

if k-tao(i)-j+1>0

fb(mm+j)=tt(k-tao(i)-j+1,i);

else

fb(mm+j)=0;

end

end;end;

mm=mm+m(i);

end;

y(k)=fb*zta';

end

e=tt(:,1)-y';

%计算损失函数值

ee=0;

for k=1:ll(1)

ee=ee+e(k)*e(k);

end

xci=ee/ll(1);

%=ee/ll(1)

%估计噪声模型

taoc=0;

lc=size(e);

c(1) = 1;

[cc,mc,taoc]=stepreg1(e);

%显示参数

m

tao

zta

mc

taoc

c(2:mc+1)=cc

plot (e)

输出结果为：

m =

2 1

tao =

0 1

zta =

-1.3348 0.6985 0.5008

mc =

taoc =

c =

1

即(1-1.3348q-1+0.6985q-2)y = 0.5008u + W; 误差曲线如下：