系统辨识的几种方法实现MATLAB代码
系统辨识，即在已知系统阶次的情况下辨识出系统的参数或辨识系统阶次，一般用到的方法有最小二乘法、辅助变量法等，对于离线与在线辨识，只是加入一个递推的问题。下面本文给出系统辨识方法的实现代码。
首先是处理数据：
U1=importdata('uy1.txt');
W1=U1.data;%第一列u,第二列y
s=length(W1(:,1));n=2;N=s-n;
A1=zeros(N,1);B1=zeros(N,4);
for i=1:N
A1(i)=W1(i+2,2);%y3,y4,y5,....
B1(i,:)=[-W1(i+1,2),-W1(i,2),W1(i+1,1),W1(i,1)];%[-y2,-y1,u2,u1]...
end
%解析法
for i=1:N
X1(:,i)=A1(i)*pinv(B1(i,:));
end
a1=sum(X1(1,:));a2=sum(X1(2,:));
b1=sum(X1(3,:));b2=sum(X1(4,:));
X=[a1/N,a2/N,b1/N,b2/N]';
%最小二乘法
% fi=B;y=A;
Xls=inv((B1'*B1))*B1'*A1;
%递推最小二乘法
c=100;
Xdls=[0;0;0;0];P=c^2*eye(4);
for i=1:N
ps=B1(i,:)';
K=P*ps*1/(1+ps'*P*ps);
P=P-P*ps*1/(1+ps'*P*ps)*ps'*P;
Xdls=Xdls+K*(A1(i)-ps'*Xdls);
end
%辅助变量法
Xiv=Xls;
u1=W1(:,1);
Z1=B1;
for j=1:1000
y1=Z1*Xiv;
y1=cat(1,W1(1:2,2),y1);
for i=1:N
Z1(i,:)=[-y1(i+1),-y1(i),u1(i+1),u1(i)];
end
Xiv0=Xiv;
Xiv=inv((Z1'*B1))*Z1'*A1;
if sum(abs(Xiv0-Xiv))<0.1%收敛
break
end
end
%递推辅助变量法
c=100;
Xdiv=Xls;P=c^2*eye(4);
u1=W1(:,1);
Z1=B1;
for i=1:N
ps=B1(i,:)';
y1=Z1*Xdiv;
y1=cat(1,W1(1:2,2),y1);
Z1(i,:)=[-y1(i+1),-y1(i),u1(i+1),u1(i)];
K=P*Z1(i,:)'*1/(1+ps'*P*Z1(i,:)');
P=P-K*ps'*P;
Xdiv=Xdiv+K*(A1(i)-ps'*Xdiv);
end
%广义最小二乘法
m=2;
Xols=Xls;
e1(1)=0;e1(2)=0;
for p=1:100
for k=n+1:n+N
e1(k)=A1(k-2)-B1(k-2,:)*Xols;
end
for i=1:N
for j=1:m
O1(i,j)=-e1(n+i-j);
end
end
f1=inv((O1'*O1))*O1'*e1(n+1:n+N)';
for k=n+1:n+N
y1(k)=W1(k,2)+f1(1)*W1(k-1,2)+f1(2)*W1(k-2,2);
u1(k)=W1(k,1)+f1(1)*W1(k-1,1)+f1(2)*W1(k-2,1);
end
y1(1)=W1(1,2);y1(2)=W1(2,2);
for i=1:N
fi(i,:)=[-y1(i+1),-y1(i),u1(i+1),u1(i)];
end
Xols0=Xols;
Xols=inv((fi'*fi))*fi'*y1(n+1:n+N)';
if sum(abs(Xols-Xols0))<0.1
break
end
end
%递推广义最小二乘法
c=100;m=2;
Xdols=Xls;
P1=c^2*eye(4);
f=[0;0];
y1(1)=W1(1,2);y1(2)=W1(2,2);
e1(1)=0;e1(2)=0;
for k=n+1:n+N
y1(k)=W1(k,2)+f1(1)*W1(k-1,2)+f1(2)*W1(k-2,2);
u1(k)=W1(k,1)+f1(1)*W1(k-1,1)+f1(2)*W1(k-2,1);
ps=[-y1(k-1),-y1(k-2),u1(k-1),u1(k-2)]';
e1(k)=A1(k-2)-B1(k-2,:)*Xdols;
w=[-e1(k-1),-e1(k-2)]';
for i=1:k-n
for j=1:m
O(i,j)=-e1(n+i-j);
end
end
for i=1:k-n
fi(i,:)=[-y1(i+1),-y1(i),u1(i+1),u1(i)];
end
P1=inv((fi'*fi));
K1=P1*ps*1/(1+ps'*P1*ps);
P1=P1-K1*ps'*P1;
Xdols=Xdols+K1*(y1(k)-ps'*Xdols);
if det(O'*O)~=0
P2=inv((O'*O));
K2=P2*w*1/(1+w'*P2*w);
P2=P2-K2*w'*P2;
f=f+K2*(e1(k)-w'*f);
end
end
%夏氏偏差修正法
Xls=inv((B1'*B1))*B1'*A1;
L=inv((B1'*B1))*B1';
M=eye(N)-B1*inv((B1'*B1))*B1';
Xx=Xls;Xb=0;
for p=1:1000
e=A1-B1*Xx;
for i=1:N-2
for j=1:2
O(i,j)=-e(n+i-j);
end
end
D=O'*M*O;
f=inv(D)*O'*M*A1;
Xb0=Xb;
Xb=L*O*f;
Xx=Xx-Xb;
if abs(sum(Xb-Xb0))<0.1
break
end
end
%夏氏改良法
L=inv((B1'*B1))*B1';
Xxg=Xls;Xb=0;
for p=1:1000
e=A1-B1*Xxg;
for i=1:N-2
for j=1:2
O(i,j)=-e(n+i-j);
end
end
f=inv((O'*O))*O'*e;
Xx0=Xxg;
Xxg=Xxg-L*O*f;
if abs(sum(Xxg-Xx0))<0.1
break
end
end
%递推的夏氏法
e(1)=0;e(2)=0;
P=eye(6);beta=zeros(6,1);
for k=1:N
ps=[-W1(k+1,2),-W1(k,2),W1(k+1,1),W1(k,1),-e(k+1),-e(k)]';
e(k+2)=A1(k)-ps'*beta;
r=P*ps*1/(1+ps'*P*ps);
P=P-r*ps'*P;
beta=beta+r*(W1(k+2,2)-ps'*beta);
end
Xdx=beta;
%增广矩阵法
Xm=[0;0;0;0;0;0];
P=eye(6);
for k=1:N
ps=[-W1(k+1,2),-W1(k,2),W1(k+1,1),W1(k,1),-e(k+1),-e(k)]';
e(k+2)=A1(k)-ps'*Xm;
K=P*ps*1/(1+ps'*P*ps);
P=P-K*ps'*P;
Xm=Xm+K*(W1(k+2,2)-ps'*Xm);
End
%残差序列相关的极大似然法
A2=A3;B2=B3;W2=W3;
Xls2=inv((B2'*B2))*B2'*A2;
c1=0.2;c2=0.1;
Xml=[Xls2;c1;c2];
sgm=0;
e(1)=0.1;e(2)=0.1;
de_a1=zeros(N+n,1);de_a2=zeros(N+n,1);
de_b1=zeros(N+n,1);de_b2=zeros(N+n,1);
de_c1=zeros(N+n,1);de_c2=zeros(N+n,1);
dJ_X=zeros(6,1);dJ_X2=zeros(6,6);
for p=1:1000
for k=n+1:n+N
ps=[-W2(k-1,2),-W2(k-2,2),W2(k-1,1),W2(k-2,1),-e(k-1),-e(k-2)]';
e(k)=A2(k-2)-ps'*Xml;
de_a1(k)=W2(k-1,2)-c1*de_a1(k-1)-c2*de_a1(k-2);
de_a2(k)=W2(k-2,2)-c1*de_a2(k-1)-c2*de_a2(k-2);
de_b1(k)=-W2(k-1,1)-c1*de_b1(k-1)-c2*de_b2(k-2);
de_b2(k)=-W2(k-2,1)-c1*de_b2(k-1)-c2*de_b2(k-2);
de_c1(k)=-e(k-1)-c1*de_c1(k-1)-c2*de_c2(k-2);
de_c2(k)=-e(k-2)-c1*de_c1(k-1)-c2*de_c2(k-2);
de_X=[de_a1(k),de_a2(k),de_b1(k),de_b2(k),de_c1(k),de_c2(k)]';
dJ_X=dJ_X+e(k)*de_X;
dJ_X2=dJ_X2+de_X*(de_X)';
end
sgm0=sgm;
sgm=S/N;
Xml0=Xml;
Xml=Xml-inv(dJ_X2)*dJ_X;
if (sgm-sgm0)/sgm0<1e-3
break
end
end
下面我们进行系统阶数的确定
%按残差方差定阶
clear;
U1=importdata('uy1.txt');
W1=U1.data;%第一列u,第二列y
U2=importdata('uy2.txt');
W2=U2.data;%第一列u,第二列y
U3=importdata('uy3.txt');
W3=U3.data;%第一列u,第二列y
V=U1.textdata;
%n=1;
n1=1;N=length(W1)-n1;
A1=zeros(N,1);B1=zeros(N,3);
A2=zeros(N,1);B2=zeros(N,3);
A3=zeros(N,1);B3=zeros(N,3);e=zeros(500,1);
J1=0;
A1=A2;B1=B2;W1=W2;
% A1=A3;B1=B3;W1=W3;
for i=1:N
A1(i)=W1(i+1,2);%y3,y4,y5,....
B1(i,:)=[-W1(i,2),W1(i+1,1),W1(i,1)];%[-y2,-y1,u2,u1]...
A2(i)=W2(i+1,2);
B2(i,:)=[-W2(i,2),W2(i+1,1),W2(i,1)];
A3(i)=W3(i+1,2);
B3(i,:)=[-W3(i,2),W3(i+1,1),W3(i,1)];
end
Xls1=inv((B1'*B1))*B1'*A1;
for k=2:N+1
e(k)=A1(k-1)-B1(k-1,:)*Xls1;
J1=J1+e(k)^2;
end
sgm1=e'*e/500
%n=2
n2=2;N=length(W1)-n2;
A1=zeros(N,1);B1=zeros(N,5);
A2=zeros(N,1);B2=zeros(N,5);
A3=zeros(N,1);B3=zeros(N,5);e(1)=0;e=zeros(500,1);
J2=0;
A1=A2;B1=B2;W1=W2;
% A1=A3;B1=B3;W1=W3;
for i=1:N
A1(i)=W1(i+2,2);%y3,y4,y5,....
B1(i,:)=[-W1(i+1,2),-W1(i,2),W1(i+2,1),W1(i+1,1),W1(i,1)];%[-y2,-y1,u2,u1]...
A2(i)=W2(i+2,2);
B2(i,:)=[-W2(i+1,2),-W2(i,2),W2(i+2,1),W2(i+1,1),W2(i,1)];
A3(i)=W3(i+2,2);
B3(i,:)=[-W3(i+1,2),-W3(i,2),W3(i+2,1),W3(i+1,1),W3(i,1)];
end
Xls2=inv((B1'*B1))*B1'*A1;
for k=3:N+2
e(k)=A1(k-2)-B1(k-2,:)*Xls2;
J2=J2+e(k)^2;
end
sgm2=e'*e/500
%n=3
n3=3;N=length(W1)-n3;
A1=zeros(N,1);B1=zeros(N,7);
A2=zeros(N,1);B2=zeros(N,7);
A3=zeros(N,1);B3=zeros(N,7);e=zeros(N+n3,1);
J3=0;
A1=A2;B1=B2;W1=W2;
% A1=A3;B1=B3;W1=W3;
for i=1:N
A1(i)=W1(i+3,2);%y3,y4,y5,....
B1(i,:)=[-W1(i+2,2),-W1(i+1,2),-W1(i,2),W1(i+3,1),W1(i+2,1),W1(i+1,1),W1(i,1)];%[-y2,-y1,u2,u1]...
A2(i)=W2(i+3,2);
B2(i,:)=[-W2(i+2,2),-W2(i+1,2),-W2(i,2),W2(i+3,1),W2(i+2,1),W2(i+1,1),W2(i,1)];
A3(i)=W3(i+3,2);
B3(i,:)=[-W3(i+2,2),-W3(i+1,2),-W3(i,2),W3(i+3,1),W3(i+2,1),W3(i+1,1),W3(i,1)];
end
Xls3=inv((B1'*B1))*B1'*A1;
for k=4:N+3
e(k)=A1(k-3)-B1(k-3,:)*Xls3;
J3=J3+e(k)^2;
end
sgm3=e'*e/500
%n=3
n4=4;N=length(W1)-n4;
A1=zeros(N,1);B1=zeros(N,9);
A2=zeros(N,1);B2=zeros(N,9);
A3=zeros(N,1);B3=zeros(N,9);e=zeros(N+n4,1);
J4=0;
A1=A2;B1=B2;W1=W2;
% A1=A3;B1=B3;W1=W3;
for i=1:N
A1(i)=W1(i+4,2);%y3,y4,y5,....
B1(i,:)=[-W1(i+3,2),-W1(i+2,2),-W1(i+1,2),-W1(i,2),W1(i+4,1),W1(i+3,1),W1(i+2,1),W1(i+1,1),W1(i,1)];%[-y2,-y1,u2,u1]...
A2(i)=W2(i+4,2);
B2(i,:)=[-W2(i+3,2),-W2(i+2,2),-W2(i+1,2),-W2(i,2),W2(i+4,1),W2(i+3,1),W2(i+2,1),W2(i+1,1),W2(i,1)];
A3(i)=W3(i+4,2);
B3(i,:)=[-W3(i+3,2),-W3(i+2,2),-W3(i+1,2),-W3(i,2),W3(i+4,1),W3(i+3,1),W3(i+2,1),W3(i+1,1),W3(i,1)];
end
Xls4=inv((B1'*B1))*B1'*A1;
for k=5:N+4
e(k)=A1(k-4)-B1(k-4,:)*Xls4;
J4=J4+e(k)^2;
end
sgm4=e'*e/500
plot([n1,n2,n3,n4],[J1,J2,J3,J4])
hold on
plot(n1,J1,'*',n2,J2,'*',n3,J3,'*',n4,J4,'*')
xlabel('n')
ylabel('J')
title('按估计误差方差最小法')
%F检验法
t2=(J1-J2)/J2*(500-2*2)/2*(2-1);
t3=(J2-J3)/J3*(500-2*3)/2*(3-2);
t4=(J3-J4)/J4*(500-2*4)/2*(4-3)
%白噪声的AIC法
%na=1,nb=2;
n1=1;N=length(W1)-n1;
A1=zeros(N,1);B1=zeros(N,4);
e=zeros(500,1);
for i=1:N-1
A1(i)=W1(i+1,2);
B1(i,:)=[-W1(i,2),W1(i+2,1),W1(i+1,1),W1(i,1)];
end
Xls1=inv((B1'*B1))*B1'*A1;
for k=2:N+1
e(k)=A1(k-1)-B1(k-1,:)*Xls1;
end
sgm12=e'*e/N
%na=1;nb=3;
A1=zeros(N,1);B1=zeros(N,5);
e=zeros(500,1);
for i=1:N-2
A1(i)=W1(i+1,2);
B1(i,:)=[-W1(i,2),W1(i+3,1),W1(i+2,1),W1(i+1,1),W1(i,1)];
end
Xls1=inv((B1'*B1))*B1'*A1;
for k=2:N+1
e(k)=A1(k-1)-B1(k-1,:)*Xls1;
end
sgm13=e'*e/N
%na=1;nb=4;
A1=zeros(N,1);B1=zeros(N,6);
e=zeros(500,1);
for i=1:N-3
A1(i)=W1(i+1,2);
B1(i,:)=[-W1(i,2),W1(i+4,1),W1(i+3,1),W1(i+2,1),W1(i+1,1),W1(i,1)];
end
Xls1=inv((B1'*B1))*B1'*A1;
for k=2:N+1
e(k)=A1(k-1)-B1(k-1,:)*Xls1;
end
sgm14=e'*e/N
系统辨识的几种方法实现MATLAB代码相关推荐
- python演示验证图像叠加过程_Python叠加矩形框图层2种方法及效果代码实例
本篇文章小编给大家分享一下Python叠加矩形框图层2种方法及效果代码实例,代码介绍的很详细,小编觉得挺不错的,现在分享给大家供大家参考,有需要的小伙伴们可以来看看. 两种方式以及效果: 方式一,使用 ...
- 由浅入深CIL系列:5.抛砖引玉：判断string是否为空的四种方法的CIL代码看看效率如何？...
本节将接触几个新的CIL操作码如下 ldc.i4.0 将整数值 0 作为 int32 推送到计算堆栈上 Ceq 比较两个值.如果这两个值相等,则将整数值 1 (int32) ...
- matlab 仿真光学实验报告,光学实验数值仿真的三种方法及MATLAB实现
光学实验数值仿真的三种方法及 MATLAB实现 5 结 论 (1)数值模拟结果表明三种方法都能对光学 实验现象进行正确地仿 真,因此在课 堂教学 中适 当应用这种仿真模拟 ,将光学实验 中复杂的数学 ...
- 图像灰度化的三种方法（matlab、C++、Python实现）
灰度化处理就是将一幅色彩图像转化为灰度图像的过程.彩色图像分为R,G,B三个分量,分别显示出红绿蓝等各种颜色,灰度化就是使彩色的R,G,B分量相等的过程.灰度值大的像素点比较亮(像素值最大为255,为 ...
- 随机信号功率谱密度函数理论、估计方法及MATLAB代码
文章的内容整理自网络,仅Matlab代码部分进行了部分修正,具体而言: 理论部分来自:现代通信原理2.5:确定信号的能量谱密度.功率谱密度与自相关函数 估计和代码部分来自: 随机信号功率谱密度估计 P ...
- 独家 | 使用Python实现机器学习特征选择的4种方法（附代码）
作者:Sugandha Lahoti 翻译:李洁 校对:杨光 本文约3500字,建议阅读13分钟. 本文中,我们将研究从数据集中选择特征的不同方法;同时通过使用Python中Scikit-learn ...
- java 隐藏标题栏_两种方法一句代码隐藏Activity的标题栏
把Activity的标题栏隐藏有两种方法.一种是在在Activity里面设置javacode.还有一种是在项目的清单文件AndroidManifest.xml中设置模版样式. 一.在Activity中 ...
- matlab 求矩阵秩,求矩阵秩的两种方法及MATLAB的应用
摘 要: 高等代数是一门逻辑思维比较强和理论知识比较深的学科, 它具有丰富的数学知识, 涉及许多重要的数学思想, 其在数学领域的应用很广泛, 如行列式.矩阵的相关计算和求解线性方程组的解方面的应 ...
- Caputo 分数阶一维问题基于 L1 逼近的空间二阶方法(附Matlab代码)
Caputo 分数阶一维问题基于 L1 逼近的空间二阶方法 Caputo 分数阶一维问题基于 L1 逼近的快速差分方法(附Matlab程序) 文章目录 Caputo 分数阶一维问题基于 L1 逼近的空 ...
最新文章
- Google发布三大新品，Pixel手机价格直逼苹果
- 带有支付功能的产品如何进行测试
- POJ2455 Secret Milking Machine【二分，最大流】
- 题目1188：约瑟夫环
- Matplotlib实例教程（十一）堆栈图
- bzoj 1968: [Ahoi2005]COMMON 约数研究【枚举】
- linux g++ 关闭 ‘typedef’ 警告_Linux学习13CentOS安装mysql5.6环境
- how is my real odata request hijacked by Mock server
- Python 开发者在迁移到 Go(lang) 时需要知道哪些事？
- hdu 4280 最大流sap
- java treemap_Java TreeMap lastEntry（）方法与示例
- Python 输入一些数，统计最大值及其出现的频率，求一个数的全部质因数
- QT设置画笔/画刷颜色
- svn —— 版本回退
- 计算机复制功能快捷键,电脑复制快捷键是什么（全部复制粘贴的快捷键是什么）...
- 数字化住宅小区对计算机网络有需求,浅谈智能小区宽带接入及其技术发展趋势...
- conda创建管理虚拟环境
- 天下武功唯快不破------实验吧
- jdk1.8的安装教程
- HTML转义字符、Javascript转义字符、HTML特殊字符对照表
热门文章
- 决战秋招 -- 经典面试题集锦
- android开发常用的ADB命令
- linuxoracle图形界面无法跳出_Linux 7图形化安装Oracle或者其他软件，打不开图形界面的问题 | 信春哥，系统稳，闭眼上线不回滚！...
- 【python】Dpark源码分析
- 【python】批量实现modis数据的辐射定标，大气校正及地形校正
- vscode在html看到图片的插件_利用花瓣插件 下载高清大图
- GPU 编程与CG 语言之阳春白雪下里巴人——CG学习读书笔记之数学函数（之一）。
- atoi()函数的实现
- 高数(下) 第八章：空间解析集合与向量代数
- Excel合并多列增加指定字符指定字符替换为换行符调整行高步骤