下载文档原格式
Matlab中的系统辨识与参数估计技术Matlab（Matrix Laboratory）是一款强大的数学软件，被广泛应用于科学计算、数据处理和工程设计等领域。
在实际工程项目中，经常需要通过已有的数据来推断系统的行为模型，这就涉及到系统辨识与参数估计技术。
本文将介绍在Matlab中使用系统辨识与参数估计技术的方法和步骤。
一、系统辨识与参数估计的概念系统辨识和参数估计是在给定输入输出数据的前提下，通过数学或统计方法来推断系统的动态模型和参数值的过程。
系统辨识旨在从实验数据中提取出模型的结构信息，而参数估计则是为了获得模型的具体参数值。
二、离散时间系统的辨识与参数估计对于离散时间系统，常用的辨识方法有ARX、ARMA和ARMAX等。
以ARX 模型为例，其数学表达式为：y(t) = -a(1)y(t-1) - a(2)y(t-2) - … - a(na)y(t-na) + b(1)u(t-1) + b(2)u(t-2) + … +b(nb)u(t-nb)其中，y(t)表示系统的输出，u(t)表示系统的输入，a和b分别是系统的参数。
在Matlab中，可以使用System Identification Toolbox来进行辨识和参数估计。
首先，需要将实验数据导入到Matlab中，然后根据数据的性质选择合适的辨识方法和模型结构。
接下来，使用辨识工具箱提供的函数，通过最小二乘法或最大似然估计等算法来得到系统的参数估计值。
三、连续时间系统的辨识与参数估计对于连续时间系统，常用的辨识方法有传递函数模型、状态空间模型和灰色系统模型等。
以传递函数模型为例，其数学表达式为：G(s) = num(s)/den(s)其中，num(s)和den(s)分别是系统的分子和分母多项式。
在Matlab中，可以使用System Identification Toolbox或Control System Toolbox 来进行连续时间系统的辨识和参数估计。
ITD（Intervaled Time Domain）法是一种在时域内识别模态参数的方法，它通过对信号的时域特征进行分析，得到系统的模态参数。
在MATLAB环境下，可以使用以下步骤实现ITD法识别模态参数的编程：
1. 读入信号数据
首先需要读入需要进行模态参数识别的信号数据。
可以使用MATLAB中的`load`函数或直接从文件中读取数据。
2. 对信号进行预处理
在进行模态参数识别之前，需要对信号进行预处理，以消除噪声等干扰因素。
可以使用MATLAB中的各种信号处理函数进行滤波、去噪等操作。
3. 计算信号的时域特征
ITD法需要计算信号的时域特征，包括信号的均值、方差、峰值等。
可以使用MATLAB中的`mean`、`var`、`peak`等函数计算这些特征值。
4. 确定模态参数
根据ITD法的原理，通过对信号的时域特征进行分析，可以得到系统的模态参数。
具体来说，可以根据信号的峰值分布情况确定系统的模态频率和阻尼比等参数。
可以使用MATLAB中的`histogram`函数绘制信号的峰值分布图，并手动或自动确定模态参数。
5. 实现ITD法识别程序
将以上步骤编写成程序，即可实现ITD法识别模态参数的功能。
可以使用MATLAB中的脚本文件或函数文件进行编写。
需要注意的是，ITD法是一种比较简单的模态参数识别方法，适用于某些特定的应用场景。
在实际应用中，需要根据具体情况选择合适的模态参数识别方法。
基于模态参数识别的ITD算法改进李玉刚;叶庆卫;周宇;王晓东【摘要】固有时间尺度分解(ITD)算法在前处理和系统定阶方面存在一定的人为因素,对模态参数的提取会造成误差,且对噪声较为敏感.针对上述问题,提出一种改进的ITD算法.利用基于数据驱动的随机子空间算法对原始数据进行处理,将正交三角分解得到的数据作为ITD法的输入数据,采用稀疏优化正交匹配追踪算法求出特征矩阵,并通过特征矩阵计算特征值、模态频率和阻尼比.通过统计的方法,从众多模态参数中选取真实模态,有效避免虚假模态的产生.实验结果表明,与ITD算法相比,改进ITD算法可降低噪声的影响,解决系统模型阶次必须准确定阶的要求,使模态参数的提取更加精确.%The pre-processing and system order determination of Intrinsic Time-scale Decomposition(ITD) algorithm involve some human factors,which causes error on extraction of modal parameters,and ITD algorithm is also sensitive to noise.To solve these problems,an improved ITD algorithm is proposed.Firstly,the stochastic subspace algorithm based on data driving is used to deal with the original data.The data obtained by orthogonal triangular decomposition is used as the input data of ITD method.The sparse optimization orthogonal matching pursuit algorithm is used to find the feature matrix characteristic value,modal frequency and damping ratio can be calculated.The real modal is selected from many modal parameters through the method of statistics,which effectively avoids the false pared with ITD algorithm,the improved ITD algorithm reduces the influence of noise,solves the problem of systemorder determination and makes the extraction of modal parameters more precise.【期刊名称】《计算机工程》【年(卷),期】2017(043)004【总页数】6页(P298-303)【关键词】固有时间尺度分解算法;模态参数;模型阶次;稀疏优化;相对误差【作者】李玉刚;叶庆卫;周宇;王晓东【作者单位】宁波大学信息科学与工程学院,浙江宁波 315211;宁波大学信息科学与工程学院,浙江宁波 315211;宁波大学信息科学与工程学院,浙江宁波315211;宁波大学信息科学与工程学院,浙江宁波 315211【正文语种】中文【中图分类】TP391.4模态参数[1]识别是桥梁健康监测的重要组成部分,现阶段发展比较成熟的时域方法主要有固有时间尺度分解(Intrinsic Time-SCAle Decomposition,ITD)法[2-3]和随机子空间(Stochastic Subspace Identification,SSI)[4-5]算法等,其优点在于无需贵重的激励识别,使得激励未知的系统辨识问题成为可能。
最近在做一个项目的方案设计，应各位老总的要求，只有系统框图和器件选型可不行，为了凸显方案设计的高大上，必须上理论分析，炫一下“技术富”，至于具体有多大实际指导意义，那就不得而知了！本人也是网上一顿百度，再加几日探索，现在对用matlab 实现系统辨识有了一些初步的浅薄的经验，在此略做一小节。
必须要指出的是，本文研究对象是经典控制论理最简单最常用的线性时不变的siso 系统，而且是2阶的哦，也就是具有如下形式的传递函数：121)(22++=Ts s T s G ξ 本文要做的就是，对于有这样传递函数的一个系统，要辨识得到其中的未知数T ， ξ!!这可是控制系统设计分析的基础哦，没有系统模型，啥理论、算法都是白扯，在实际工程中非常重要哦！ 经过总结研究，在得到系统阶跃响应实验数据之后（当然如果是其他响应，也有办法可以辨识，在此还是只讨论最简单的阶跃响应实验曲线，谁让你我是菜鸟呢），利用matlab 至少可以有两种方法实现实现（目前我只会两种，呵呵）！一、函数法二、GUI 系统辨识工具箱下面分别作详细介绍！一、 函数法看官别着急，先来做一段分析（请看下面两排红*之间部分），这段分析是网上找来的，看看活跃一下脑细胞吧，如果不研读一下，对于下面matlab 程序，恐怕真的就是一头雾水咯！*******************************************************************************G(s)可以分解为：))((1)(212ωω++=s s T s G其中， [][]11112221--=-+=ξξωξξωTT1ω、2ω都是实数且均大于零。
则有：211ωω=T ，21212ωωωωξ+=传递函数进一步化为：))(()(2121ωωωω++=s s s G 因此，辨识传递函数就转化为求解1ω、2ω。
当输入为单位阶跃函数时，对上式进行拉普拉斯反变换，得系统时域下的单位阶跃响应为：tteet y 212111221)(ωωωωωωωω---+--=即 tteet y 21211122)(1ωωωωωωωω-----=-令1ω=2ωk )1(>k，得tk t ek e k k t y 22111)(1ωω-----=-⎥⎦⎤⎢⎣⎡--=---t k t e k e k k 2)1(2111ωω 对上式两边取以e 为底的对数得[]⎥⎦⎤⎢⎣⎡-+--=---t k e k t k k t y 2)1(211ln 1ln )(1ln ωω 当∞→t 时，⎥⎦⎤⎢⎣⎡---t k e k 2)1(11ln ω0→，则上式化简为[]t k k t y 21ln )(1ln ω--=-该式的形式满足直线方程b at t y +=)(*其中，)(*t y =[])(1ln t y -，1ln ,2-=-=k kb a ω)1(>k通过最小二乘算法实现直线的拟合，得到a ，b 的值，即可得到1ω、2ω的值，进而可得系统的传递函数。
基于速度响应信号的模态参数识别ITD算法研究
武晓飞;郑海起;栾军英;赵亚星
【期刊名称】《机械制造》
【年(卷),期】2009(047)009
【摘要】以环境随机激励作用下的大型塔架结构为研究背景,提出了一种基于速度响应信号的ITD方法.首先应用随机减量技术从速度随机响应信号中提取速度自由响应信号,然后推导出基于速度信号的ITD法并应用于识别系统的模态参数.该方法不用将速度信号积分成位移信号,而是直接对速度信号进行模态参数识别,避免了积分过程中引入误差.将试验结果与传统的基于位移响应的ITD法和峰值谱法进行对比,证明该方法具有较高的计算精度.
【总页数】4页(P29-32)
【作者】武晓飞;郑海起;栾军英;赵亚星
【作者单位】军械工程学院火炮工程系,石家庄,050003;军械工程学院火炮工程系,石家庄,050003;军械工程学院火炮工程系,石家庄,050003;军械工程学院火炮工程系,石家庄,050003
【正文语种】中文
【中图分类】TH113
【相关文献】
1.基于响应协方差小波变换和SVD的结构工作模态参数识别 [J], 徐晓霞;任伟新;韩建刚
2.基于传递函数有理分式的模态参数识别算法研究 [J], 于亚婷;李敏;姜敏;邹宇
3.基于响应功率谱传递比的应变模态参数识别方法研究 [J], 曹林波; 颜王吉; 任伟新
4.基于脉冲响应数据的ARMA法建模以及模态参数识别 [J], 郭永刚;许亮华;水小平
5.基于速度信号的大型结构模态参数识别ITD方法 [J], 栾军英;张晓涛;武晓飞因版权原因，仅展示原文概要，查看原文内容请购买。
模态参数是指结构动力特性的基本参数，是描述结构动力特性的基本概念，包括固有频率、阻尼比、振型等。
结构模态参数的准确识别，是进行结构健康监测及故障诊断的重要基础，直接关系到结构安全，因此，开展结构模态参数识别技术研究具有重要的理论意义与工程实用价值。
近年来，利用环境激励已大量应用于土木工程的结构动力特性测试中。
环境激励测试能够在结构的实际工作状态下进行，更真实地了解结构的动力特性和结构性能。
本文将对各种模态识别方法进行分类汇总、论述，并对环境激励下模态参数识别算法有待进一步研究的问题进行了展望。
1频域识别算法1.1峰值拾取法基于结构的频响函数在其固有频率位置处会出现峰值的特征，可以实现对结构的模态参数识别。
由于环境激励下无法得到结构的频响函数，用功率谱密度函数代替结构的频响函数实现模态参数的识别，功率谱由实测的随机振动信号快速傅立叶变化转化得到。
姜蕾蕾[1]将幂指数窗应用于多种结构中，并与其他五种窗函数对比研究，确定能够有效改善傅立叶变换后频谱的质量,从而提高峰值拾取法的频率和阻尼比识别精度,拓宽峰值拾取法对阻尼比的适用范围。
陈涛[2]将测点传递率函数矩阵的第2阶奇异值倒数的均值为模态指示函数，建立基于多参考测点平均的峰值拾取法，准确识别系统的模态频率及振型。
在实际应用中，该方法只需计算少量的局部极值点，识别速度快，适用性广泛，被大量使用在实测实验中。
但由于峰值拾取法对峰值的选择较为敏感，对于峰值存在干扰或者峰值较小的信号，可能导致参数提取不准确，并且输出结果可能受到峰值选择的主观性影响，存在一定的不确定性。
因此，在使用时需要综合考虑实际需求和信号特征，选择合适的峰值。
1.2频域分解法频域分解法是峰值拾取法的优化算法，基本原理是根据振动响应构建谱函数矩阵，通过奇异值分解，将多自由度系统转换为单自由度体系，依靠峰值法选取特征频率，进而对系统进行识别。
频域分解法在20世纪80年代由Prevosto[3]所提出。
实用标准文案4. 设某物理量Y 与X 满足关系式Y=aX 2+bX+c ，实验获得一批数据如下表，试辨识模型参数a ，b 和c 。
（50分）报告要求：要有问题描述、参数估计原理、程序流程图、程序清单，最后给出结果及分析。
(1)问题描述：由题意知，这是一个已知模型为Y=aX 2+bX+c ，给出了10组实验输入输出数据，要求对模型参数a ，b ，c 进行辨识。
这里对该模型参数辨识采用递推最小二乘法。
(2)参数估计原理对该模型参数辨识采用递推最小二乘法，即RLS （ recurisive least square ），它是一种能够对模型参数进行在线实时估计的辨识方法。
其基本思想可以概括为：新的估计值)(ˆk θ=旧的估计值)1(ˆ-k θ+修正项 下面将批处理最小二乘法改写为递推形式即递推最小二乘参数估计的计算方法。
批处理最小二乘估计θˆ为Y T TΦΦΦ=-1)(ˆθ，设k 时刻的批处理最小二乘估计为：k T k k T k Y ΦΦΦ=-1)(ˆθ令111)]1()()1([)()(----+-=ΦΦ=k k k P k P T kT k ϕϕ K 时刻的最小二乘估计可以表示为k T k Y k P k Φ=)()(ˆθ=)]()()[(11k y k Y k P k T k ϕ+Φ-- =)]1(ˆ)()()[()1(ˆ--+-k k k y k K k Tθϕθ；式中)()()(k k P k K ϕ=，因为要推导出P(k)和K(k)的递推方程，因此这里介绍一下矩阵求逆引理：设A 、（A+BC ）和（I +B CA 1-）均为非奇异方阵，则111111)()(------+-=+CA B CA I B A A BC A 通过运用矩（3)程序流程图(如右图1所示）递推最小二乘法（RLS）步骤如下：已知：n、b n和d。
aStep 1 ：设置初值)0(ˆθ和P(0)，输入初始数据；Step2 ：采样当前输出y(k)、和输入u(k)Step3 ：利用上面式①②③计算)(k K、)(ˆkθ和)(k P；Step4 ：k→k+1，返回step2，继续循环。
使用Matlab技术进行系统辨识的基本方法概述：系统辨识是指通过对已知输入输出数据的分析和处理，推断出系统的动态性质和数学模型的过程。
在科学研究、工程设计和控制应用中，系统辨识扮演着重要的角色。
而Matlab作为一种强大的数值计算和数据分析软件，为系统辨识提供了便利且高效的工具。
本文将介绍使用Matlab进行系统辨识的基本方法，并结合实例进行讲解。
一、数据采集与准备在进行系统辨识之前，首先需要采集到对应的输入输出数据。
一般来说，输入信号是已知的，可以通过外部激励或者系统自身的变动来获取；而输出信号则是根据输入信号通过系统响应得到的。
在采集数据时，需要注意数据的质量和采样频率的选择。
二、数据预处理在进行系统辨识之前，数据通常需要进行一些预处理，以去除噪声、平滑数据和调整时间序列等。
这可以通过Matlab中的数据处理函数和滤波器实现。
例如，可以使用高斯滤波器对数据进行平滑处理，或者使用降噪算法去除不必要的噪声。
三、参数估计参数估计是系统辨识的核心步骤之一，它通过对已知数据进行分析和处理，推断出系统的数学模型和参数。
在Matlab中，有多种方法和工具可供选择，如最小二乘法、最大似然法、系统辨识工具箱等。
这些工具可以根据不同的模型和数据类型灵活选择，并提供相应的算法和函数。
四、模型验证与优化根据估计得到的系统模型和参数，可以使用Matlab进行模型验证和优化。
模型验证是指将估计得到的模型与真实系统进行对比，检验其拟合程度和预测能力。
如果模型的拟合程度较差，则需要对参数进行调整和优化，以提高模型的准确性和稳定性。
五、模型预测与应用在系统辨识完成之后，可以使用得到的模型进行系统预测和应用。
通过对未知输入信号进行预测，可以得到相应的输出响应，进而实现对系统动态性质的分析和控制。
Matlab提供了丰富的预测和应用函数，例如时域模拟、频域分析、控制系统设计等，可以满足不同应用场景的需求。
六、案例分析为了更好地理解和掌握使用Matlab进行系统辨识的基本方法，下面通过一个简单的案例进行分析。
%ITDxx识别模态参数
clear
clc
close all hidden
format long
%% txt文件下输入
fni=input('ITD法模态参数识别-输入数据文件名：','s');
fid=fopen(fni,'r');
mn=fscanf(fid,'%d',1);%模态阶数
%定义输入实测数据类型
%ig=1时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果%ig=2频域数据如频响函数实部和虚部数据
ig=fscanf(fid,'%f',1);
%ig=1时，f为采样频率sf，ig=2时，f为频率间隔df
f=fscanf(fid,'%f',1);
fno=fscanf(fid,'%s',1);%输出数据文件名
b=fscanf(fid,'%f',[ig,inf]);%实测时域或频域数据
status=fclose(fid);
%%
clc;
clear all;
format long
[FileName,PathName] = uigetfile('*.mat', 'Select the Mat-files of time signal'); %窗口读文件,并获取包含路径的文件名
if isequal(FileName,0)
disp('User cancel the selection');%如果取消选择则显示提示
return;
else
FULLFILE=fullfile(PathName,FileName);
Signal_str= sprintf('User selected signal file:%s',FULLFILE);
disp(Signal_str);
Struct=load(FULLFILE);
end
c=fieldnames(Struct);%得到一个元胞数组，包含Struct中各个域名（倘若有多个的话）
b=getfield(Struct,c{1}); %获取c{1}对应的域中的内容
b=b(3601:9600);
%%
%ig=1时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果
%ig=2频域数据如频响函数实部和虚部数据
ig=input('数据类型ig=');
f=input('采样频率f=');%指定采样频率
mn=input('计算模态阶数mn=');%指定计算模态阶数
%建立特征方程矩阵的阶数（为模态阶数的2倍）
nm=2*mn;
%组织识别计算多用的时域数据及参数
if ig==1%实测时域数据
sf=f;%采样频率
n=fix(length(b)/2);%向0靠拢取整，取时域数据的长度
h=b(1,1:2*n)';%将输入时域数据赋值给列向量h
dt=1/sf;%时间间隔
t=0:dt:(2*n-1)*dt;%建立离散时间向量
else %实测频域数据
df=f;%取频率间隔
n=length(b(1,:));
f=0:df:(n-1)*df;%建立离散频率向量
H=b(1,:)'+b(2,:)'*i;%建立对应正负频率的实测频响函数向量
H(n+1)=real(H(n));
H(n+2:2*n)=conj(H(n:-1:2));%conj求负数的共轭值
h=real(ifft(H));%频响函数经IFFT并取实部变换成脉冲响应函数t=linspace(0,1/df,2*n);%建立离散时间向量
dt=t(2)-t(1);%计算时间间隔
end
%计算自由振动响应矩阵
L=length(h);
M=L/2;
for k=1:nm
x1(k,:)=h(k:L-(nm-k+1))';
x2(k,:)=h(k+1:L-(nm-k))';
end
%用最小二乘法求解特征方程矩阵
B=x1\x2;%B=x2*x1'*inv(x1*x1');
[A,V]=eig(B);%计算特征值及特征向量（特征值V，特征向量A）%变换特征值对角阵为一向量
for k=1:nm
U(k)=V(k,k);
end
F1=abs(log(U'))./(2*pi*dt);%计算模态频率向量
D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));%计算阻尼比向量%计算振型系数向量
l=1;
for k=0:(2*n-1)
Va(k+1,:)=[conj(U).^k];
end
S1=(inv(conj(Va')*Va)*conj(Va')*h);%inv矩阵求逆
h1=real(Va*S1);%计算生成的脉冲响应函数%绘制脉冲响应函数拟合曲线图
figure(1);
plot(t,h,':',t,h1);
xlabel('时间（s）');
ylabel('幅值');
legend('实测','拟合');
grid on;
if ig>1
H1=fft(Va*S1);%计算生成的频响函数
%绘制频响函数实部拟合曲线图
figure(2);
nn=1:n;
subplot(2,1,1);
plot(f,real(H(nn)),':',f,real(H1(nn)),'r'); xlabel('频率(Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,b(2,:),':',f,imag(H1(nn)),'r');
xlabel('频率(Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
end
[F2,I]=sort(F1);%将自振频率从小到大排列
%剔除方程中的非模态项（非共轭根）和共轭项（重复项）m=0;
for k=1:1:nm-1
if F2(k)~=F2(k+1)
continue;
end
m=m+1;
l=I(k);
F(m)=F1(l);%自振频率
D(m)=D1(l);%阻尼比
S(m)=S1(l);%振型系数
end
%打开文件输出识别的模态参数数据
fno='out.txt';
fid=fopen(fno,'w');
fprintf(fid,'频率（Hz）阻尼比（%%）振型系数\n');
for k=1:m
fprintf(fid,'%10.4f %10.4f%10.6f\n',F(k),D(k)*100.0,imag(S(k))); end
status=fclose(fid);。
