系统辨识作业及答案
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
一. 问答题
1. 介绍系统辨识的步骤。
答:(1)先验知识和建模目的的依据;(2)实验设计;(3)结构辨识;(4)参数估计;(5)模型适用性检验。
2. 考虑单输入单输出随机系统,状态空间模型
[])
()(11)()(11)(0201)1(k v k x k y k u k x k x +=⎥⎦⎤
⎢⎣⎡+⎥⎦⎤⎢⎣⎡=+ 转换成ARMA 模型。
答:ARMA 模型的特点是u(k)=0,
[])
()(11)()(0201)1(k v k x k y k x k x +=⎥⎦
⎤
⎢⎣⎡=+
3. 设有一个五级移位寄存器,反馈取自第2级和第3级输出的模2加法和。
试说明:
(1) 其输出序列是什么? (2) 是否是M 序列?
(3) 它与反馈取自第4级与第3级输出模2加法和所得的序列有何不同? (4) 其逆M 序列是什么? 答:(1)设设输入序列1 1 1 1 1
111018110107101006010015100114001113011112111111)()()()()()()()(()()()()()()()01110161110115110101410100)13(010011210011110011110011109()()()()()()()001112401110)23(111012211010211010020010011910011180011117()()()()()()()()10011
3200111310111030001112911010281010027010012610011
25 其输出序列为:1 1 1 1 1 0 0 1 0 1
⑵不是M 序列
⑶第4级与第3级模2相加结果
100108001007010006100015000114001113011112111111)()()()()()()()(()()()()()()()11110161110115110101410101)13(010111210110110110010110019()()()()()()()110012410010)23(001002201000211000120000111900111180111117()()()()()()()()01111
3211110311110130110102910101280101127101102601100
25 不同点:第2级和第3级模二相加产生的序列,是从第4时刻开始,每隔7个时刻重复一次;第4级与第3级模2相加产生的,序列,是从第2时刻开始每隔15个时刻重复一次。
⑷第5级与第4级模2相加结果如下:已知其为M 序列。
001008010007100006000015000114001113011112111111)()()()()()()()(()()()()()()()10100
160100115100111400110)13(011001211000111000110000109()()()()()()11111221111021111012011010
1910101180101017
M 序列: 1 1 1 1 1 0 0 0 0 1 0 0 0 1 1 0 0 1 0 1 0 1 方波信号: 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 逆重复M: 0 1 0 1 0 0 1 0 1 1 1 0 1 1 0 0 1 1 1 1 1 1
4. 画出广义最小二乘法的离线迭代算法的简单计算框图。
答:广义最小二乘法的离线迭代算法的简单计算框图如下:
5. 考虑如下数学模型x b x a y cos sin +=,试用N k k y k x ,...,3,2,1),(),(=输入输出数据
估计系统参数b a ,。
答:()()()
1cos 1sin 1x b x a y +=
()()()2cos 2sin 2x b x a y +=
()()()N x b N x a N y cos sin +=
Y E Y T T φφφθφθ1-)(=⇒+=⇒∧
6. 利用最小二乘算法辨识如下模型参数
z(k)-1.5z(k-1)+0.7z(k-2)=u(k-1)+0.5u(k-2)+v(k)
其中,v(k)是零均值白噪声。
当模型阶次为2时,可以获得准确的辨识结果,而模型阶次取3时,只能得到如下一组模型参数辨识结果(括号内为模型参数真值): a 1=-1.08884(-1.5) a 2=0.08326(0.7) a 3=0.28781(0.0) b 1=1.00000(1.0) b 2=0.91116(0.5) b 3=0.20558(0.0) 显然,辨识结果已经远远偏离了模型参数真值,试从理论上解释为什么会出现这种现象。
答:对于n 阶系统与n+1阶系统参数估计之间有如下的关系:
对于n+1阶系统
()()()11()()A z y k B z u k e k --=+
设其待估参数为
()011111...(1)(2)T T T
n n n n n b a b a b a b θθθ++⎡⎤⎡⎤+==⎣
⎦⎣⎦ 则(1)()[()]T
n A Y n θθθ=-Φ-Φ
由题目知n=2时系统参数为准确值,则n=3时按照上式去计算,估算出的系数必远
远偏离系统模型参数值。
7. 请说明闭环系统不可辨识的原因。
答:闭环系统不可辨识的原因:反馈使得一个闭环系统对不同的输入常产生差不多相同的输出,观测的输入输出数据所包含的信息比开环辨识少的多;输入信号与噪声因反馈而相关:有偏估计,非一致性估计;在闭环条件下,用开环辨识方法系统的参数有时也是不可唯一辨识的。
8. 设闭环系统前向通道模型为
)()2(7.0)1()2(45.0)1(4.1)(k k u k u k y k y k y ε+-+-+----=
反馈调节器为)1(2.0)()(-+=k y k y k u
试画出其闭环系统框图,并判断系统是否可辨识?
答:系统是可以辨识的,由于为非奇异,故在)1(2.0)()(-+=k y k y k u 条件下,参数是可以辨识的。
闭环系统框图如下图所示:
9. 对系统模型阶次进行辨识,得到1阶-4阶的参数估计,性能指标与系统模型阶次的关系
解:由F 检验法原理知
()(1)(,1)(1)J n J n t n n J n -++=
+
若
(,1) 3.09 t n n+≤
则可以接受系统阶数。
由计算得,t(1,2)=4.13 , t(2,3)=0.49 , t(3,4)=0.0034, t(4,5)=0
所以系统的阶数为3。
二.编程题
1.(1)编程产生一组正态分布的白噪声信号,它的均值和方差以及长度可随意调整,将
产生的白噪声信号存入数据文件data1.txt
(2)编程产生一组M序列信号,它的幅值和长度可随意调整,将产生的M序列存入数据文件data2.txt
(3)编程产生一组逆重复M序列信号,它的幅值和长度可随意调整,将产生的逆重复M序列存入数据文件data3.txt
解:(1)function y=WNoise(N,E,V AR)
% N为长度E为均值V AR为方差
y=randn(1,N);
y=y-mean(y);
y=y/std(y);
y=E+sqrt(V AR)*y;
plot(y)
title('严晓龙实验:产生一组正态分布的白噪声信号')
save data1.txt y -ascii
调用函数实验:WNoise(400,0,1),得到数据见data1.txt,如图所示:
(2)function seq=mseq(a,L,N)
% a为M序列幅值N为长度L为移位单位数
register=randint(1,L) %寄存器初始化
p=zeros(1,L); %特征向量
p(L-1:L)=1; %默认最后两个寄存器相加
temp=0;
for i=1:N
seq(i)=a*register(L);
temp=sum(register.*p);
register(2:L)=register(1:L-1); %移位
register(1)=mod(temp,2);
end
x=0:1/5:1.2;
stairs(seq);grid;
set(gca,'ylim',[-0.2,1.2]);
ylabel('M序列')
title('严晓龙实验:移位寄存器产生的M序列')
save data2.txt seq -ascii
调用函数实验:mseq(2,40,15),得到数据见data2.txt,和下图:
(3)function seq=invM(a,L,N)
% a为M序列幅值N为长度L为移位单位数register=randint(1,L) %寄存器初始化
p=zeros(1,L); %特征向量
p(L-1:L)=1; %默认最后两个寄存器相加temp=0;
for i=1:2^L-1
seq(i)=register(L);
temp=sum(register.*p);
register(2:L)=register(1:L-1); %移位
register(1)=mod(temp,2);
end
seq=[seq seq];
for i=1:2*(2^L-1)
if mod(i,2)==1
invm(i)=1;
else
invm(i)=0;
end
seq(i)=a*xor(seq(i),invm(i));
end
for i=1:N
if mod(i,2*(2^L-1))==0
mseq(i)=seq(2*(2^L-1));
else
mseq(i)=seq(mod(i,2*(2^L-1)));
end
end
seq=mseq;
stairs(seq);grid;
set(gca,'ylim',[-0.2,1.2]);
title('严晓龙实验:产生一组逆重复M序列信号')
save data3.txt seq -ascii
调用函数实验:invM (1,10,40),得到数据见data3.txt,和下图
2.12.mat中的数据是单输入单输出系统进行采样后100对输入输出数据,其中input表示
系统的输入数据,output表示受到噪声污染后的系统的输出数据。
在辨识过程中,可以认为噪声具有正态分布,其均值为0。
(1)判断该系统的阶次(方法不限)
(2)利用递推最小二乘法进行参数估计。
解:
模型阶数的辨识,一般说来低阶模型描述粗糙,高阶模型精度高。
残差平方总和J(n)是
模型阶数的函数∑
=-
=
N
k
T
K k
y
n
J
1
2
)
)
(
(
)
(θ
ϕ
在不同的模型阶数的假设下,参数估计得到的J(n)值亦不同。
讨论如下(1)当n=1时程序如下:
启动matlab,打开12.mat;运行下面程序
u=zeros(100,1);%构造输入矩阵
z=zeros(100,1);%构造输出矩阵
i=1:1:100;
u(i,1)=input(i);
z=zeros(100,1);%构造输出矩阵
i=1:1:100;
z(i,1)=output(i);
r=100;
for p=1:(r-2) %利用循环生成观测矩阵
h(p,:)=[-z(p+1) u(p+1)]; %
end
hl=h;
for b=1:(r-2) %生成输出矩阵
zl(b,:)=[z(b+2)];
zl'
end
zl'
%根据最小二乘法公式进行参数辩识
c1=hl'*hl;
c2=inv(c1);
c3=hl'*zl;
c=c2*c3;
a1=c(1)
a2=c(2)
j=0;
for k=4:100;
hl=[-z(k-1);u(k-1)]';
x=hl*c;
y=z(k)-x;
s=y*y;
j=j+s;
end
j
仿真结果如下
a1 = -0.2576 a2 = 0.6985 j = 0.8556
(2)当 n=2时程序如下(输入输出数据同上,只给出不同于一阶系统的程序不同之处)
其中U、Z分别是作业要求给出得的输入输出,数据输入同上。
启动matlab,打开12.mat;运行下面程序
u=zeros(100,1);%构造输入矩阵
z=zeros(100,1);%构造输出矩阵
i=1:1:100;
u(i,1)=input(i);
z=zeros(100,1);%构造输出矩阵
i=1:1:100;
z(i,1)=output(i);
r=100;%利用循环生成观测矩阵。
for p=1:(r-2)
h(p,:)=[-z(p+1) -z(p) u(p+1) u(p)];
end
hl=h;
%生成输出矩阵。
for b=1:(r-2)
zl(b,:)=[z(b+2)];
zl'
end
zl'
%根据最小二乘法公式进行参数辩识
c1=hl'*hl;
c2=inv(c1);
c3=hl'*zl;
c=c2*c3;
%输出辩识参数
a2=c(2)
b1=c(3)
b2=c(4)
j=0;
%求J(n)
for k=4:100; %开始求K
h1=[-z(k-1),-z(k-2),u(k-1),u(k-2)]';
x=h1'*c;
y=z(k)-x;
s=y*y;
j=j+s;
end
j
仿真结果如下:
a1 = -0.4362 a2 = 0.2407 b1 = 1.8844 b2 = -1.1313 j = 0.5977
(3)当n=3时程序如下
启动matlab,打开12.mat;运行下面程序
u=zeros(100,1);%构造输入矩阵
z=zeros(100,1);%构造输出矩阵
i=1:1:100;
u(i,1)=input(i);
z=zeros(100,1);%构造输出矩阵
i=1:1:100;
z(i,1)=output(i);
r=100;
for p=2:(r-1)
h(p,:)=[-z(p+1) -z(p) -z(p-1) u(p+1) u(p) u(p-1)];
end
hl=h;
for b=2:(r-1)
zl(b,:)=[z(b+1)];
zl';
end
zl';
c1=hl'*hl;
c2=inv(c1);
c3=hl'*zl;
c=c2*c3;
a1=c(1)
a2=c(2)
a3=c(3)
b2=c(5)
b3=c(6)
j=0;
for k=4:100;
hl=[-z(k-1);-z(k-2);-z(k-3);u(k-1);u(k-2);u(k-3)]';
x=hl*c;
y=z(k)-x;
s=y*y;
j=j+s;
end
j
仿真结果如下
a1 = -1.0000 a2 = 0 a3 = 2.1316e-14 b1 = -2.2737e-12 b2 = 1.3642e-12 b3 = 1.8190e-12 j = 2.6845
数据分析如下利用LS法先对系统参数进行初步辩识并根据其确定残差平方
总和
我们可以计算出相应的J(N)变化率:
由此我可以看出,在由二阶到三阶的T(n)很大,说明J(N)有一个很明
显的变化,而T(2)很小,说明随着阶数的增加J(N)明显变大,据此我们有
理由相信该系统应该是二阶的。
综上所述:第12组的数据的输入输出模型是2阶的,且a1 = -0.4362 a2 = 0.2407 b1 = 1.8844 b2 = -1.1313 j = 0.5977。