东南大学自动化学院
实 验 报 告
课程名称： 系统辨识
实验名称： 系统辨识与仿真
院（系）： 自动化 专 业： 自动化
姓 名： 学 号：
同组人员： 实验时间： 2013 年 月 日
评定成绩： 审阅教师：
实验一：
给定系统G(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)
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
双线性变化为
，代入原传递函数求得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
B1(q-1)=0.2q-2 B2(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);
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)<max(m-1),m(1)=max(m-1);end %限制A(z)的阶次不低于任何B(z)的阶次
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)<ll(2)
%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 =
0
taoc =
0
c =
1
即(1-1.3348q-1+0.6985q-2)y = 0.5008u + W;
误差曲线如下：
