matlab使用相关分析法和hankel矩阵法进行系统辨识
待辨识模型：
%H = 180.60 / ( 48.96s^2 + 14.00s + 1.0 )
H = tf( [180.60],[ 48.96 14 1.0 ] ) ;
b = [ 180.60 ] ;
a = [ 48.96 14 1.0 ];
待辨识模型转化为离散模型,可以适用脉冲响应不变法，也可以使用双线性变换法方法，采样时间为1S
采样时间的标准，采样时间选择标准如下，fmax为系统截止频率
[bz,az] = impinvar(b,a,1 );%使用脉冲响应不变法将S域转到Z域
[bz1,az1] = bilinear(b,a,1) %使用双线性变换的方法，将S域转到Z域
通过系统的阶跃响应的图形，观察系统的过渡时间
取过渡时间为50s
选择M序列的寄存器个数
n = 6 ;
N = 2^n - 1 ;
即M序列的一个周期为63个点
%生成M序列用于辨识系统的模型
n = 6 ;%M序列取4个移位寄存器
N = 2^n - 1 ; %总的M序列的点数
%Tuse * （ Np-1）> T95 ( 在范围(1.2 ~ 1.5 )T95 )
a = 1;
delta = T0 ;
%M序列出来的值会跟初值有关，所以出来的值根据初值不同有多种情况
A = zeros( n , 1 ) ; %A数组
X = zeros( n , 1 ) ; %X数组
XOR = zeros( n , 1 ) ; %XOR数组
A(1) = 1 ;
A(2) = 1 ;
A(3) = 1 ;
A(4) = 1 ;
A(5) = 1 ;
A(6) = 0 ;
XOR( 6 ) = 1 ;
XOR( 5 ) = 1 ;
r = 4 ; %周期数
for i = 1 : rN
X(1) = xor( A(6), A(5) ) ;
for i1 = 2 : n
X( i1 ) = A( i1 -1 ) ;
end
OUT( i ) = A(n);
t( i ) = deltai ;
if OUT( i ) > 0.5
u( i ) = a ;
else u( i ) = -a;
end
for j = 1 : n
A( j ) = X( j ) ;
end
end
通过生成的M序列，输入系统的离散模型中，生成离散模型的输出数据Y
%Ny 表示方程的阶次，进行输入输出数组的填充
Ny = 2 ;
Y = zeros( r * N + Ny , 1 );
X = zeros( r * N + Ny , 1 );
YPulse = zeros( r * N + Ny , 1 );
XPulse = zeros( r * N + Ny , 1 );
XPulse(3) = a ;
% YPulse1 = zeros( r * N + Ny , 1 );
% XPulse1 = zeros( r * N + Ny , 1 );
% XPulse1(3) = 1/(T0) ;
for i = 1 : r * N
%X( i + Ny ) = u( i ) ;
X( i + Ny ) = u( i ) ;
end
for i = 1 : r * N
Y( i + Ny ) = 0.8033 * X( i + Ny ) + …
1.6065 X( i + Ny - 1 ) + …
0.8033 * X( i + Ny - 2 ) + …
1.7332 * Y( i + Ny - 1 ) - …
0.7510 * Y( i + Ny - 2 );
end
接下来采用相关分析法，仅仅根据输入数据和输出数据，得到系统的脉冲响应曲线
Y1 = Y( Ny + ( r - 1 ) * N + 1 : r * N + Ny ) ;
X1 = X( Ny + ( r - 1 ) * N + 1 : r * N + Ny ) ;
YSum = 0 ;
for i = 1 : N
YSum = YSum + Y1( i ) ;
end
YAvr = YSum / N ;
Y1 = Y1 - YAvr ;%减去直流分量
Rmz = zeros( N , 1) ;
%计算互相关矩阵，输入与输出的互相关矩阵，进行卷积
for i = 1 : N
Rmz( i )= 0 ;
sum = 0 ;
for j = 1 : N
if( j >= i )
sum = sum + X1( j - i + 1 ) * Y1( j ) ;
else
sum = sum + X1( N - ( i - j ) + 1 )*Y1( j ) ;
end
end
Rmz( i ) = sum / N ;
end
c = -Rmz( N ) ;
gK = ( N * ( Rmz + c ))/( ( N + 1 ) * ( a ^ 2 ) )
gY = impz( bz1,az1 ,50,1) ;
figure(1)
plot( 1 : N , gK , ‘r’)
接下来使用Hankel矩阵的方法，获取系统的模型的参数
%通过Hankel矩阵进行模型的辨识
Nsys = 2 ;%待辨识的系统的系统的阶次
Gsys1 = zeros( Nsys , Nsys ) ;
Asys1 = zeros( Nsys , 1 ) ;
Gsys2 = zeros( Nsys , 1 ) ;
%注意Gsys是从1开始的
for i = 1 : Nsys
for j = 1 : Nsys
Gsys1( i , j ) = gK( i + j ) ;
end
Gsys2( i ) = -gK( Nsys + i + 1 ) ;
end
Asys1 = inv( Gsys1 ) * Gsys2
Bsys1 = zeros( Nsys + 1 , 1 ) ;
Gsys2 = [ gK( 1 : Nsys + 1 ) ] ;
Asys2 = eye(Nsys + 1 ) ;
for i = 2 : Nsys + 1
Asys2( i , : ) = [ Asys1( Nsys - ( i - 2 ):Nsys )’ , Asys2( i , i : Nsys + 1 ) ] ;
end
Asys2;
Gsys2;
Bsys1 = Asys2 * Gsys2 ;
Aresult = [ 1 , flip( Asys1’ ) ];
Bsys1
Asys1
GY1 = impz( Bsys1 ,[ 1 , flip( Asys1’ ) ] , 50 , 1);
figure(2);
% plot( 1 : size(gY) , gY , ‘r’ , 1 : size(YPulse) , YPulse , ‘b’, 1:size(YPulse1),YPulse1,‘g’);
plot( 1 : N , gK , ‘r’ , 1 : size(GY1) , GY1 , ‘b’);
拟合出的曲线与使用相关分析法分析出来的曲线能很好的重合
辨识出来的模型参数，asys1为传递函数的分母的参数
待辨识的系统的离散模型的参数
附上代码段：
%H = 180.60 / ( 48.96s^2 + 14.00s + 1.0 )
H = tf( [180.60],[ 48.96 14 1.0 ] ) ;
figure(3) ;
step( H ) ;
T95 = 35 ; %上升时间是35S
T0 = 1; %采样时间取上升时间的( 1/5 ~ 1/15 ) ;
%如果采样时间取整形，则出来的阶跃响应的波形会存在误差
%改为使用浮点形，则可以减小误差
b = [ 180.60 ] ;
a = [ 48.96 14 1.0 ];
f = 0 : 0.01 : 2 ;
w = 2*pi*f ;
% figure(6)
% freqs(b,a,w)
[bz,az] = impinvar(b,a,1 );%使用脉冲响应不变法将S域转到Z域
[bz1,az1] = bilinear(b,a,1) %使用双线性变换的方法，将S域转到Z域
% [bz2,az2] = bilinear(b,a,1/T0);
%使用脉冲响应不变法得到的模型
G = tf( bz ,az , T0 ) ;
%figure(3) ;
%step( G ) ;
%使用双线性变换后得到的模型
% G1 = tf( bz1 , az1 , 1 ) ;
% G2 = tf( bz2 , az2 , T0 ) ;
% figure(3) ;
% step(G1) ;
% figure(4) ;
% step(G2) ;
% figure(5);
% step(G)
%生成M序列用于辨识系统的模型
n = 6 ;%M序列取4个移位寄存器
N = 2^n - 1 ; %总的M序列的点数
%Tuse * （ Np-1）> T95 ( 在范围(1.2 ~ 1.5 )*T95 )
a = 1;
delta = T0 ;
%M序列出来的值会跟初值有关，所以出来的值根据初值不同有多种情况
A = zeros( n , 1 ) ; %A数组
X = zeros( n , 1 ) ; %X数组
XOR = zeros( n , 1 ) ; %XOR数组
A(1) = 1 ;
A(2) = 1 ;
A(3) = 1 ;
A(4) = 1 ;
A(5) = 1 ;
A(6) = 0 ;
XOR( 6 ) = 1 ;
XOR( 5 ) = 1 ;
r = 4 ; %周期数
for i = 1 : r*N X(1) = xor( A(6), A(5) ) ;for i1 = 2 : nX( i1 ) = A( i1 -1 ) ;endOUT( i ) = A(n); t( i ) = delta*i ; if OUT( i ) > 0.5 u( i ) = a ; else u( i ) = -a; end for j = 1 : n A( j ) = X( j ) ; end
end
%Ny 表示方程的阶次，进行输入输出数组的填充
Ny = 2 ;
Y = zeros( r * N + Ny , 1 );
X = zeros( r * N + Ny , 1 );
YPulse = zeros( r * N + Ny , 1 );
XPulse = zeros( r * N + Ny , 1 );
XPulse(3) = a ;
% YPulse1 = zeros( r * N + Ny , 1 );
% XPulse1 = zeros( r * N + Ny , 1 );
% XPulse1(3) = 1/(T0) ;
for i = 1 : r * N %X( i + Ny ) = u( i ) ;X( i + Ny ) = u( i ) ;
end
%进行差分方程的计算
% for i = 1 : r * N
% Y( i + Ny ) = 3.688 * X( i + Ny ) + ...
% 7.376 * X( i + Ny - 1 ) + ...
% 3.688 * X( i + Ny - 2 ) + ...
% 1.488 * Y( i + Ny - 1 ) - ...
% 0.5099 * Y( i + Ny - 2 );
% end
for i = 1 : r * N Y( i + Ny ) = 0.8033 * X( i + Ny ) + ...1.6065* X( i + Ny - 1 ) + ...0.8033 * X( i + Ny - 2 ) + ...1.7332 * Y( i + Ny - 1 ) - ...0.7510 * Y( i + Ny - 2 );
end
% for i = 1 : N
% YPulse1( i + Ny ) = 3.688 * XPulse1( i + Ny ) + ...
% 7.376 * XPulse1( i + Ny - 1 ) + ...
% 3.688 * XPulse1( i + Ny - 2 ) + ...
% 1.488 * YPulse1( i + Ny - 1 ) - ...
% 0.5099 * YPulse1( i + Ny - 2 );
% end
for i = 1 : N YPulse( i + Ny ) =0.8033 * XPulse( i + Ny ) + ...1.6065 * XPulse( i + Ny - 1 ) + ...0.8033 * XPulse( i + Ny - 2 ) + ...1.7332 * YPulse( i + Ny - 1 ) - ...0.7510 * YPulse( i + Ny - 2 );
end
figure(4)
plot( 1:size( Y( Ny + 1 : r * N + Ny ) ) , Y( Ny + 1 : r * N + Ny ), 'r ' , t , u , 'b' ) ;
Y1 = Y( Ny + ( r - 1 ) * N + 1 : r * N + Ny ) ;
X1 = X( Ny + ( r - 1 ) * N + 1 : r * N + Ny ) ;
YSum = 0 ;
for i = 1 : N YSum = YSum + Y1( i ) ;
end
YAvr = YSum / N ;
Y1 = Y1 - YAvr ;%减去直流分量
Rmz = zeros( N , 1) ;
%计算互相关矩阵，输入与输出的互相关矩阵，进行卷积
for i = 1 : N Rmz( i )= 0 ; sum = 0 ;for j = 1 : Nif( j >= i )sum = sum + X1( j - i + 1 ) * Y1( j ) ;elsesum = sum + X1( N - ( i - j ) + 1 )*Y1( j ) ; end endRmz( i ) = sum / N ;
end
c = -Rmz( N ) ;
gK = ( N * ( Rmz + c ))/( ( N + 1 ) * ( a ^ 2 ) )
gY = impz( bz1,az1 ,50,1) ;
figure(1)
plot( 1 : N , gK , 'r')%通过Hankel矩阵进行模型的辨识
Nsys = 2 ;%待辨识的系统的系统的阶次
Gsys1 = zeros( Nsys , Nsys ) ;
Asys1 = zeros( Nsys , 1 ) ;
Gsys2 = zeros( Nsys , 1 ) ;
%注意Gsys是从1开始的
for i = 1 : Nsysfor j = 1 : NsysGsys1( i , j ) = gK( i + j ) ;endGsys2( i ) = -gK( Nsys + i + 1 ) ;
end
Asys1 = inv( Gsys1 ) * Gsys2
Bsys1 = zeros( Nsys + 1 , 1 ) ;
Gsys2 = [ gK( 1 : Nsys + 1 ) ] ;
Asys2 = eye(Nsys + 1 ) ;
for i = 2 : Nsys + 1Asys2( i , : ) = [ Asys1( Nsys - ( i - 2 ):Nsys )' , Asys2( i , i : Nsys + 1 ) ] ;
end
Asys2;
Gsys2;
Bsys1 = Asys2 * Gsys2 ;
Aresult = [ 1 , flip( Asys1' ) ];
Bsys1
Asys1
GY1 = impz( Bsys1 ,[ 1 , flip( Asys1' ) ] , 50 , 1);
figure(2);
% plot( 1 : size(gY) , gY , 'r' , 1 : size(YPulse) , YPulse , 'b', 1:size(YPulse1),YPulse1,'g');
plot( 1 : N , gK , 'r' , 1 : size(GY1) , GY1 , 'b');
参考书籍：
系统辨识理论及应用
知识网站：
https://wenku.baidu.com/view/e4e37083581b6bd97f19ea91.html
matlab使用相关分析法和hankel矩阵法进行系统辨识相关推荐
- 用matlab的编程法和游动鼠标法求二阶传递函数的上升时间、峰值时间、超调量和调节时间 - Gavin_Hall的博客 - CSDN博客
1. 准备 终值:c(∞) 上升时间 tr:响应从峰值的10%上升到峰值的90%所需要的时间:而阶跃响应则是从终值的10%上升到终值的90%所需要的时间:对有振荡的系统,也可以定义为从0到第一次到达终 ...
- 【光学】基于matlab GUI矩阵法和等效界面法光学薄膜对反射率影响【含Matlab源码 2102期】
⛄一.获取代码方式 获取代码方式1: 完整代码已上传我的资源:[光学]基于matlab GUI矩阵法和等效界面法光学薄膜对反射率影响[含Matlab源码 2102期] 点击上面蓝色字体,直接付费下载, ...
- matlab riccati法 临界转速,利用传递矩阵法和Riccati传递矩阵法分析转子临界转速...
利用传递矩阵法和Riccati传递矩阵法分析转子临界转速 利用传递矩阵法和Riccati传递矩阵法分析转子临界转速 一. 所需求解转子参数 将转子简化为如下所示: 三个盘的参数为: 另,阶梯轴的三段轴 ...
- matlab层次分析法代码_基于主成分分析法和层次分析法的工程项目经理胜任力评价研究...
摘 要:根据工程项目经理胜任力评价指标,运用主成分分析法和层次分析法相结合的数学方法对工程项目经理的胜任力进行合理公正的评价.首先运用主成分分析法筛选重要指标,再运用层次分析法对工程项目经理进行定量与 ...
- matlab可视化大学物理学_传输矩阵法在大学物理波动光学教学中的应用
1 提出问题 在现代光学技术中,从基本的光学元件增反膜和增透膜[1],到超快光路中用来补偿飞秒激光色散的啁啾镜,以及半导体微腔领域中广泛使用的分布式布拉格反射器(DBR)[2],这些光学元件基本的特征 ...
- 用相关法辨识系统的脉冲响应 matlab,基于相关分析法的系统辨识算法对比及仿真...
计算机工程应用技术 ComputerKnowledgeand Technology 电脑知识第12卷第9期 (2016年3月) 基于相关分析法的系统辨识算法对比及仿真 冀征难 (国防科技大学 机电工程 ...
- 线性代数 --- Gauss消元的部分主元法和完全主元法
Gauss消元的部分主元法和完全主元法 心怀二意的人,在他一切所行的路上都没有定见.----雅各书1章8节 笔者的一些话:刚开始写这篇文章的时候,我觉得高斯消元很简单.因为,这时的我已经完成了我一直想 ...
- matlab max与min获取矩阵最大最小值函数
1.matlab中Max的用法1(Min类似) Matlab中max函数在矩阵中求函数大小的实例如下: C = max(A) 1)返回一个数组各不同维中的最大元素. 2)如果A是一个向量,max(A) ...
- matlab2c使用c++实现matlab函数系列教程-hankel函数
全栈工程师开发手册 (作者:栾鹏) matlab2c动态链接库下载 matlab库函数大全 matlab2c基础教程 matlab2c开发全解教程 matlab2c调用方法: 1.下载动态链接库 2. ...
最新文章
- android 保存 用户名和密码 设置等应用信息优化
- 第十六届全国大学生智能车提问与回复 |7月10日
- //yield return用于无缝实现迭代模式。
- maven中snapshot快照库与maven-metadata.xml
- 字符串暴力匹配算法+思路分析
- 将python代码编译成.so文件
- 第 8 天 多线程与多进程
- 蓝桥杯历届试题-六角填数（12）
- sitemesh2.4
- 2019最新猎豹网校JAVA语言数据结构与算法教程(Java语言 )
- 步进电机驱动器细分原理_步进驱动器细分设置表说明
- YOLOV5 + 双目测距（python）
- 含有一个量词的命题的否命题_第三节：简单的逻辑联结词、全称量词与存在量词...
- 计算机功能自定义,设计大师学教学：自定义鼠标右键功能提升CAD绘图效率-鼠标右键菜单设置...
- java 读取word模板文件路径_Java 读取Word模板替换内容并另存
- 七大江河水系－－黄河(二)
- 懂车帝与蛋蛋订车两大平台对比
- airtest获取设备号和获取设备宽度、高度、绝对坐标 相对坐标、滑动屏幕
- C#基础 uint,long,ulong,float,decimal 定义并初始化
- 扫盲贴 家用无线路由的安全设置