matlab使用相关分析法和hankel矩阵法进行系统辨识

最新推荐文章于 2023-10-07 15:42:30 发布
最新推荐文章于 2023-10-07 15:42:30 发布
阅读量4.7k 收藏 51
点赞数 11
分类专栏： 系统辨识 文章标签： matlab 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_39647439/article/details/104124972
版权

待辨识模型：
%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 : n
        X( i1 ) = A( i1 -1 ) ;
    end
    OUT( 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 : 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矩阵进行模型的辨识
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');

参考书籍：
系统辨识理论及应用

知识网站：
https://wenku.baidu.com/view/e4e37083581b6bd97f19ea91.html

我名字里有只鸟
关注
  • 11
    点赞
  • 踩
  • 51
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
专栏目录
线性代数笔记：汉克尔矩阵（Hankel matrix）
qq_40206371的博客
02-28 5771
汉克尔矩阵 (Hankel Matrix) 是指每一条副对角线上的元素都相等的矩阵
用相关法辨识系统的脉冲响应 matlab,利用相关分析法辨识脉冲响应
weixin_39870664的博客
03-18 524
利用相关分析法辨识脉冲响应自1205 刘彬 412511411 实验方案设计1.1 生成输入数据和噪声用M 序列作为辨识的输入信号，噪声采用标准正态分布的白噪声。 生成白噪声时，首先利用乘同余法生成U[0,1]均匀分布的随机数，再利用U[0,1]均匀分布的随机数生成标准正态分布的白噪声。 1.2 过程仿真模拟过程传递函数)(s G ，获得输出数据y(k)。)(s G 采取串联传递函数仿真，2121...
系统建模与仿真之利用相关分析法辨识脉冲响应
最新发布
programminginc的博客
10-07 827
在这里分享我们《系统建模与仿真》课程的实验课的matlab编程题目，虽然我写的可能不太好，但是代码也可以给初学者提供参考。，但是分成了多个函数，看着不方便，而且有些可以修改，供同校以及其他有需要的人参考。
Hankel矩阵法进行系统辨识
weixin_43789096的博客
08-09 5630
Hankel矩阵
论文笔记：Hankel Matrix Factorization for Tagged Time Series to Recover Missing Values during Blackouts
qq_40206371的博客
02-28 1380
ICDE 2019 0 摘要 在执行时间序列分析时，恢复时间序列中的缺失值至关重要。而本文研究的blackouts问题，即在一定时期内丢失所有数据，是最紧迫和最具挑战性的问题之一。现有的时间序列缺失值恢复方法无法正确处理这个问题，但在这项工作中，我们提出了一种基于 Hankel 矩阵分解的标记时间序列方法，称为 HKMF-T，遵循将数据序列分解为平稳趋势和外部影响成分。通过将数据序列转换为其 Hankel 矩阵形式，HKMF-T 将高阶时间相关性隐含的平滑趋势建模为两个低秩矩阵的乘积，并...
相关分析法和hankel矩阵法进行系统辨识
01-31
基于matlab的采用相关分析法和hankel矩阵法进行系统辨识，使用的是《系统辨识理论及应用》这本书的一个传递函数，通过自己的理解，写出来的进行系统辨识的算法。有兴趣可以自己去买一下这本书回来看，很好的一本教材
Hankel.rar_Hankel矩阵_hankel matrix_传递函数辨识_汉克尔矩阵_系统辨识 Hankel
07-14
《Hankel矩阵及其在系统辨识中的应用》 Hankel矩阵，又称汉克尔矩阵，是线性系统理论中的一个...对于工程师和研究人员来说，掌握Hankel矩阵的相关理论和应用技巧，无疑将有助于提升他们在系统分析和设计中的专业水平。
SystemParamter.rar_levy_levy辨识_传递函数辨识_系统参数辨识_面积法
07-15
总的来说，系统参数辨识是控制系统设计和分析的基础，通过Hankel矩阵法、Levy法、Levy辨识和面积法等技术，我们可以从实验数据中提取出系统的动态特性，为控制系统的建模、仿真和优化提供有力支持。这些方法的熟练...
用Levy或hankel矩阵等方法系统辨识程序
02-23
总的来说，系统辨识是一个涉及多种数学工具和算法的复杂过程，包括但不限于面积法、Levy法、Hankel矩阵的使用和脉冲响应的分析。掌握这些方法可以帮助我们更好地理解和预测系统的动态行为，从而在控制设计、故障诊断...
Hankel 变换：分别是 Hankel 变换和逆 Hankel 变换的有效实现。-matlab开发
06-01
"hat" 和 "ihat" 使用矩阵乘积直接积分执行 n 阶汉克尔变换。 “ht”和“iht”通过先验积分贝塞尔核来执行0阶汉克尔变换。 "dht" 和 "idht" 实现了整数阶 n 的准离散 Hankel 变换。 而且，最后但并非最不重要的是，...
经典的系统辨识方法
05-16
通过给予脉冲激励信号可以求得系统的时域输出信号 。 对传递函数求反拉普拉斯变换
相关分析法辨识
04-12
系统辨识中的相关分析法 大家互相学习 仅供参考 共同进步
基于MATLAB利用相关分析法辨识脉冲响应
10-17
学习系统辨识课程时候的实验作业，仅供大家参考~~~
故障诊断——奇异值分解法&hankel矩阵
yezi12123的博客
03-23 6516
基本思路 将离散的时序信号通过矩阵的变换，转换成二维的矩阵，利用奇异值分解，提取主要特征，根据奇异值大小判断信号是否有故障。 信号转换矩阵 hanker矩阵 对于离散信号X=[x(1)x(2)...x(N)]X=\begin {bmatrix}x(1)&x(2)&...&x(N)\end{bmatrix}X=[x(1)​x(2)​...​x(N)​] hanker矩阵可以表...
Toplitz矩阵 Hankel矩阵 Hilbert矩阵
qykshr的专栏
03-20 6038
 Toplitz矩阵是指 矩阵的对角线上的元素相等 以及平行于对角线的从左上到右下的一斜列上的所有元素相等. 即 T = ( t_{i-j} ), i,j=1,2,...,n 例 T = [ 1 2 3 4 5 6 7 1 2 3 4 5 8 7 1 2 3 4 9 8 7 1 2 3 0 9 8 7 1 2 x 0 9 8 7 1
MATLAB经典建模分享一hankel矩阵
临风暖阳的博客
07-10 1437
%% hankel矩阵 clc,clear,close all c=[1:5], r=[3:9], H=hankel(c,r) c = 1 2 3 4 5 r = 3 4 5 6 7 8 9 警告: 输入列的最后一个元素与输入行的第一个元素不匹配。 在反对角线冲突中，列具有更高优先级。 In hankel (line 27) In BJMan (line 5) H = 1 2 3 4
相关分析法辨识系统脉冲响应
weixin_43789096的博客
08-11 4140
背景介绍 在实际工程中，由于系统的测量都是载噪的，而且噪声对观测数据的影响常常达到不可忽略的地步，因此当噪声影响足以使得要求的精度不足时，就必须考虑噪声的影响。实际中，系统噪声存在各种难以精确描述的因素：数学模型中未加以考虑的各种干扰作用、模型线性化及其他类似的假设所引起的误差等（系统误差）、输入输出的量测误差等（量测误差）等等，这些误差都难以用精准的数学模型来描述，且难以使用物理工具来测量。这些难以用确定性模型来精确描述的随机因素对过程的影响不容忽视，因此处理它们的方式应尽量简化。一般我们在确定性模型上以
利用相关分析法辨识脉冲响应
假电工的真的博客
06-05 6030
主程序: %% 参数初始化 times=1; Np=63;%2^6-1,输入序列循环周期 N=252*times; a=1;%输入序列幅值 T0=1;%采样时间 delta_g=zeros(200,1); ratio_vy=zeros(200,1); %% 人机对话 sigma=input('请输入噪声标准差: '); r=input('请输入循环周期数(本程序中，输入2或3): ...
使用python实现的相关分析法系统辨识
weixin_42439763的博客
06-16 1054
最近系统辨识课要求做实验来进行系统辨识，第一个方法是用相关分析法。目前有的案例基本上都是matlab程序，所以写了一个python程序来实现。首先导入需要的第三方库，其中最重要的是scipy.signal库，里面提供了构建传递函数、阶跃响应等等功能，如果没有这个库，python很难完成matlab中的简单功能。 接着构建了异或函数，后来发现其实python有第三库里面有内嵌，不过这个函数简单，自己写的也能用 生成M序列的函数在我另一个帖子里有介绍，这里直接拿来用了。M序列的python实现 总的来说这个函数

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值