matlab参数识别
最小二乘法：
%递推公式，更新 p0=p2;
for n=2:N-1%%递推最小二乘法
K0=p0*X(n,:)'*inv(1+X(n,:)*p0*X(n,:)');%计算Ｋ
Theta_abs=Theta_abs+K0*(Y(n)-X(n,:)*Theta_abs);%计算估计值Theta ； p3=p0-K0*X(n,:)*p0;%计算P p0=p3;
%误差平方和最小
Y1=X(n,:)*Theta_abs;%递推值 J=(Y(n,:)-Y1)*(Y(n,:)-Y1)'
if (J 对于()()()()()b n a n n k u b k u b n y a k y a k y b a -++=-+-+.......111 引进后移算子()()11-=-k y k y q 假定在初始条件0时z 变换得到 ()()()a b b n a n n n z a z a z b z b b z X z Y z H ----+++++==...1 (11110) ARX 模型有：()() ?????++=+++=------b b a n n n a n z b z b b q B z a z a q A (11) 1011 11 ()()() ()()k v k u q B q k y q A d +=---11；()k v 为均值为0的噪声项 上式可以改写为()()()()l k k v i k u b i k z a k z b a n i i n i i ,..,2,1,1 1 =+-+--=∑∑== 上式改写为最小二乘格式()())(k v k h k z T +=θ（3） 对于（3）式的l 次观测构成一个线性方程组 [][] ?????=------=T n n T b a n a b b b a a a n k u k u n k z k z k h ,...,,,,...,,)(),...,1(),(),...,1()(2121θ即 l l l V H Z +=θ. ()()()[]()()()[]T l T l l v v v V l z z z Z ,...2,1,,...,2,1= ()()()()()()()()()???? ? ?? ? ????-------------=????????????=)()1()(21)2()1(10)1(021b a b a b a l n l u l u n l z l z n u u n z z n u u n z z l h h h H 取极小化准则函数()()[]()()θθθL T L l l k T H z H z k h k z J --=-=∑=1 2，极小化()θJ ，求得参 数θ的估计值θ?，[] T n n n a b b b a a a ?,...,?,?,?,...,?,??2121=θ ()()[]()() θθθ???12l l T l l l k T H Z H Z k h k z J --=-=∑=表示为了确定使准则最小的条件，将该式对各参数求导，并令其结果等于零： ()() l T l l T l l l l Z H H H H Z H J 1 ?,0?2?-==--=??θθθ 即，只要矩阵l H 是满秩的，l T l H H 则是 正定的，使准则为极小的条件得到满足， 最小二乘估计的递推算法（RLS ） 最小二乘法，不仅占用大量内存，而且不适合于在线辨识，为了解决这个问题， 把它转化为递推算法：修正项+=+k k θθ??1 ()()()()()()()()()???? ? ?? ? ????-------------=????????????=)()1()(21)2()1(10)1(021b a b a b a l n l u l u n l z l z n u u n z z n u u n z z l h h h H 若令()1-=l T l l H H P ,则[][]l T l l T l l l l l T l l l P h Ph h I h P P h h P P 111111111+-+++-++++-=+= [][] l T l l l l l T l l l l l T l l l l l P h K I P h P h h P K h z K 1111 111 11111;1; ??++++++++++++-=+= -+=θθθ 加权递推最小二乘（RWLS ）： () ()() ()(),11k e k u z B k z z A +=--()()() (),11 k v z C z D k e --=e(k)为有色噪声，v(k)为白噪声。 () ()()111 ---=z A z B z G ()()() 11 1---z C z D z N 。取()()11--=z C z A ，() 11=-z D [][] l T l l l l l T l l l l l l T l l l l l P h K I P h P h w h P K h z K 1111111 1 11111;11 ;??+++++++++++++-=+= -+=θθθ 当噪声为有色噪声时，采用增广最小二乘法：其思路采用CARMA 模型。在实际应用中噪声v(k)有两种形式： ()()()()()()()()k k h k z k v k k h k z k v θθ??;1??-=--= 1 matlab 最小二乘法拟合 [a,Jm]=lsqcurvefit(fun,a0,x,y);fun 不支持句柄函数 a0为最优化的初始值，fun 为数据原型函数。 x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)； lb ≤ x ≤ ub [x,norm,res,ef,out,lam,jac]=lsqcurvefit(@F,x0,t,y,v1,v2,opt,P1,P2,...) 其中输出变量的含义为： 1） x ： 最优解 2） norm ： 误差的平方和 3）res: 误差向量 4） ef ： 程序结束时的状态指示： · >0：收敛 · 0：函数调用次数或迭代次数达到最大值（该值在options 中指定） · <0：不收敛 5) out: 包含以下数据的一个结构变量 · funcCount 函数调用次数 · iterations 实际迭代次数 · cgiterations 实际PCG 迭代次数（大规模计算用） · algorithm 实际使用的算法 · stepsize 最后迭代步长（中等规模计算用） · firstorderopt 一阶最优条件满足的情况（大规模计算用） 6) lam ：上下界所对应的Lagrange 乘子 7) jac ：结果（x 点）处的雅可比矩阵 输入参数 其中输入变量的含义为： · x0为初始解（缺省时程序自动取x0=0） · t,y: 拟合数据 · v1,v2: 参数待求x 的上下界 · options：包含算法控制参数的结构 LineSearchType 线搜索方法 (…cubicpoly?,?quadcubic?(缺省值)) opt=optimset(oldopts,newopts) 可以设定的参数比较多，对lsqnonlin 和lsqcurvefit ，常用的有以下一些参数： Diagnostics 是否显示诊断信息（ 'on' 或'off ） Display 显示信息的级别（'off' ， 'iter' ， 'final ，'notify ） LargeScale 是否采用大规模算法（ 'on' 或'off ）缺省值为on MaxIter 最大迭代次数 TolFun 函数计算的误差限 TolX 决策变量的误差限 Jacobian 目标函数是否采用分析Jacobi 矩阵（'on' ，'off ） MaxFunEvals 目标函数最大调用次数 LevenbergMarquardt 搜索方向选用LM 法(…on?), GN 法(…off?,缺省值) LineSearchType 线搜索方法(…cubicpoly?,?quadcubic?(缺省值)) LargeScale: [ on | off ] LevenbergMarquardt: [ {on} | off ] 例子1用matlab 实现对 ()4.0,4,sin 122sin sin 22===-+=k m m f λ? λ? λ??设 1 首先编写m 文件 function f =lsq(x,xdata) f=x(1)*sin(xdata)+0.5*x(2)*sin(2*xdata)./(1-x(2)^2*sin(xdata).^2).^0.5 2 利用lsqcurvefit 函数调用m 文件 m=4;k=0.4 o=[0:0.01*pi:2*pi]; xdata=o; ydata=m*sin(o)+0.5*k*sin(2*o)./(1-(k^2*sin(o).^2)).^0.5; x0 = [0; 0]; [x,resnorm] = lsqcurvefit(@lsq,x0,xdata,ydata) 结果得到：x(1)= 4.0000;x(2)=0.4000;resnorm = 6.3377e-016 2 nlinfit m=4;k=0.4; o=[0:0.005*pi:2*pi]; xdata=o; ydata=m*sin(o)+0.5*k*sin(2*o)./(1-(k^2*sin(o).^2)).^0.5; x0 = [0; 0]; beta = nlinfit(xdata,ydata,@lsq,x0) 例子1.1 用fminunc 函数； clear k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); F =@(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata); lb=[10;1;8]; ub=[20;10;15]; x0=[1 1 1]; [x,fr] = fminunc(F,x0)%fminunc [x,fr] = fminsearch(F,x0)%fminsearch options= saoptimset('Display','iter', 'TolFun',1e-30); [x,fval,exitFlag,output] = simulannealbnd(F,x0,lb,ub,options)%模拟退火[x,fval] = ga(F,x0,[],[],[],[],lb,ub)%遗传 例子1.2 用遗传算法的参数识别 k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); F =@(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata); [x,fval] = ga(F,3,[],[],[],[],[10;1;8],[20;10;15]) %ee=norm(E);%使用差平方和最小原则；或者使用sum(abs(E)); %ee=norm(E)/sqrt(n);%使用rms准则 例子1.3 利用multistart方法 k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); F =@(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata); ms=MultiStart; opts=optimset('Algorithm','interior-point', 'LargeScale','off'); problem=createOptimProblem('fmincon','x0',[10,1,8],'objective',F,'lb' ,[1,0,1],'ub',[20,10,15],'options',opts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,200) 例子1.4 利用globalsearch k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); F =@(x)norm(x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata); gs = GlobalSearch('Display','iter'); opts=optimset('Algorithm','interior-point'); problem=createOptimProblem('fmincon','x0',[10,1,8],'objective',F,'lb' ,[1,0,1],'ub',[20,10,15],'options',opts); [xming,fming,flagg,outptg,manyminsg] = run(gs,problem) 例子1.5 利用multistart和lsqcurvefit k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); ms=MultiStart; opts=optimset('Algorithm','trust-region-reflective'); problem=createOptimProblem('lsqcurvefit','x0',[10,1,8],'xdata',xdata, 'ydata',ydata,'objective',@myfun,'lb',[1,0,1],'ub',[20,10,15],'option s',opts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) function y=myfun(x,xdata) y=x(1).*exp(x(2).*xdata)+x(3).*sin(xdata); 例子1.6利用multistart和lsqnonlin(8s) k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); F =@(x)x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata; ms=MultiStart; opts=optimset('Algorithm','trust-region-reflective'); problem=createOptimProblem('lsqnonlin','x0',[10,1,8],'xdata',xdata, 'ydata',ydata,'objective',F,'lb',[1,0,1],'ub',[20,10,15],'options',op ts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) 1.7 利用matlabpool parallel加速 tic; k1=13;k2=1.3;k3=9.1; xdata=0:pi/100:pi; ydata= k1*exp(k2*xdata)+k3*sin(xdata); F =@(x)x(1).*exp(x(2).*xdata)+x(3).*sin(xdata)-ydata; ms=MultiStart('Display','iter','UseParallel','always'); opts=optimset('Algorithm','trust-region-reflective'); matlabpool open 2 problem=createOptimProblem('lsqnonlin','x0',[10,1,8],'xdata',xdata, 'ydata',ydata,'objective',F,'lb',[1,0,1],'ub',[20,10,15],'options',op ts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) time=toc;%47.3541s 函数模型 function y=myfun(x,xdata) y=x(1)*exp(x(2).*xdata)+x(3)*sin(xdata); 例子2（多输入变量） function F=myfun(k,xdata) F = k(1)*xdata(:,1)+k(2)*xdata(:,2)+k(3)*xdata(:,3); xdata = [3.6 7.7 9.3; 4.1 8.6 2.8; 1.3 7.9 10.0]' ydata = [16.5 150.6 263.1]' k0 = [0,0,0]' [x,resnorm]=lsqcurvefit(@parameter_identification,k0,xdata,ydata) 例子3 （带参数限制条件的） M函数： function F = myfun(k,xdata) F = k(1)*exp(k(2)*x)+k(3)*sin(x); Workspace 代码： k1=13;k2=1.3;k3=9.1; x=0:pi/100:pi; y = k1*exp(k2*x)+k3*sin(x); k0=[-1,-1,-1] [k,resnorm]=lsqcurvefit(@myfun,k0,x,y,[5,0.1,2],[20,8,15]) 结果k = 13.0000 1.3000 9.1000 例子3 lsqnonlin %编写M文件： function E = myfun(k,x,y) F = k(1)*exp(k(2)*x)+k(3)*sin(x); E=F-y; k1=13;k2=1.3;k3=9.1; x=0:pi/100:pi; y = k1*exp(k2*x)+k3*sin(x); k0=[-1,-1,-1] options=optimset('TolFun',1e-10,'TolX',1e-10,'MaxFunEvals ',500 ,'display','iter'); k=lsqnonlin(@myfun,k0,[5,0.1,2],[20,8,15],options,x,y) 1 变参数微分方程数值求解 例子2 求 function dydt=fun(t,y,u,v) r=u+2;s=v-2; dydt=[r+y(2); s*y(1)-2*s*y(2)]; u=[1;5;15;20;25]; v=[6;12;18;24;30]; tspan=0:1:4; y0=[0 2]; yy=y0; for i=1:length(tspan)-1 [t,y]=ode45(@fun,[tspan(i),tspan(i+1)],y0,[],u(i),v(i)); y0=y(end,: ); yy=[yy;y0]; end plot(tspan,yy,'-o') 2.1 匿名函数法 f=@(t,y,u,v) [u+2+y(2); (v-2)*y(1)-2*(v-2)*y(2)] u=[1;5;15;20;25]; v=[6;12;18;24;30]; tspan=0:1:4; y0=[0 2]; yy=y0; for i=1:length(tspan)-1 [t,y]=ode45(f,[tspan(i),tspan(i+1)],y0,[],u(i),v(i)); y0=y(end,: ); yy=[yy;y0]; end plot(tspan,yy,'-o') 2.2 修改加上时间tt(显示所有计算值) clear u=[1;5;15;20;25]; v=[6;12;18;24;30]; tspan=0:1:4; y0=[0 2]; tt =[];yy=[]; for i=1:length(tspan)-1 [t,y]=ode45(@fun,[tspan(i),tspan(i+1)],y0,[],u(i),v(i)); y0=y(end,: ); tt=[tt;t]; yy=[yy;y]; end plot(tt,yy);%所有的计算数值。 2.3 同过差值可以调节精度。如果u,v随着t是时刻变化的，但是通过测试手段 只能测得某一时刻的u,v. clear global yy t1=0:0.1:4;%如果u,v随着t是时刻变化的，可以通过此数值来调节精度 tspan=0:1:4; u=[1;5;15;20;25]; u1=spline(tspan,u,t1); v=[6;12;18;24;30];v1= spline(tspan,v,t1); y0=[0 2]; yy=y0; for i=1:length(t1)-1 [t,y]=ode45(@fun,[t1(i),t1(i+1)],y0,[],u1(i),v1(i)); y0=y(end,: ); yy=[yy;y0]; end plot(t1,yy) 2 适用matlab对一个常微分方程进行参数回归 问题如下： 已知实验数据x,y，并且x,y的关系满足以下常微分方程 Dy/dx=-k*(y-y0)*y^2 其中k是需要回归的参数，y0是一个常数，通常等于y向量中的最后一个数值。要求： 1.通过lsqcurvefit或者lsqnonlin回归出系数k 2.画出模型预测值和实验值的对比图，模型预测值可以通过得到常微分方程数值解后三次样条spline插值得到。我已经写好的程序如下，里面有错误，我自己找不出来，请高手帮帮忙，谢谢啊 可以加我的QQ交流：40231185 ======================================= function odetest clc;clear; global Je J0 data=xlsread('flux.xls'); xdata=data(:,1)';ydata=data(:,2)'; beta0=0.1;Je=ydata(end);J0=ydata(1); options=optimset('TolFun',1e-20,'TolX',1e-20,'MaxFunEvals',100,'Algorithm','trust-region-reflecti ve'); beta=lsqcurvefit(@cakefun,beta0,xdata,ydata,[],[],options); Jc=cakefun(beta,xdata); plot(xdata,ydata,'o',xdata,Jc); function y=cakefun(beta,x) global J0 tspan=[0 max(x)]; [m,n]=size(x); [tt yy] = ode23s(@modeleqs,tspan,J0,[],beta); yc=spline(tt',yy',x); y=yc; function dydt = modeleqs(t,y,beta) global Je dydt =-beta*(y-Je)*y.^2; =========================== 数据如下： X Y 0 1176.918115 1 637.4225727 2 326.1218103 3 265.7631528 4 220.666083 5 200.7988265 7 181.6298121 9 170.1634436 11 162.4684024 14 151.6322759 17 150.3811328 20 146.8242069 25 139.8735164 30 137.5861186 35 135.1093977 40 131.7195994 45 131.631951 53 126.4159865 60 125.0219926 70 123.041967 80 121.7344741 90 120.8522371 100 116.8166294 110 117.2211892 120 114.8271487 130 113.2291363 140 113.2365511 150 112.9655866 3变参数微分方程系数回归 您好斑竹，我的最终目的是要识别k参数，我的微分方程模型为 y'=k(1)(y-y0)*y^k(2)+v ; (我的参数识别方法是参考一个例子，后面给贴出来。但是这个例子是不带v输入变量的）1 我前面一段程序先假设k=3,然后利用变参量数值解法求解了（仿真）y值；(因为这个式子有个带外部输入的参数v，没法求得解析解，所以只好用数值微分方程方法求出ydata； 2 我后面一段程序就是利用前边得到的t值和y值去识别参数k； 3.1 生成数据xxdata为2 的x数据。 function dydt=modelsq(t,y,k,v) dydt =k(1)*(y-6).^k(2)+v; clear; k=[-3.548 3.236]; y0=6; x=[0 1 2 3 4 5 7 9 11 14 17 20 25 30 35 40 45 53 60 70 80 90 100 110 120 130 140 150]';%给出x值 fraction=500;%计算出500个数据； t1=min(x):(max(x)-min(x))/fraction:max(x); v1=1:length(x); v=spline(x,v1,t1); yy=y0; for i=1:length(t1)-1 [t y]=ode23s(@modelsq,[t1(i) t1(i+1)],y0,[],k,v(i)); y0=y(end,: ); yy=[yy;y0]; end xydata=[t1',yy,v']; save('xydata.txt','xydata','-ascii');%把tt和yy存为txt文件。 3.2 参数识别 方法1 lsqcurvefit function dydt=modelsq(t,y,k,v) dydt =k(1)*(y-6).^k(2)+v; function yi=num_value(k,x) global v; y0=6; yy=y0; for i=1:length(x)-1 [t y]=ode23s(@modelsq,[x(i) x(i+1)],y0,[] ,k,v(i)); y0=y(end,: ); yy=[yy;y0]; end yi=yy; clear; global v; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); v=xydata(:,3); k0=[4000,2000]'; %初值也可以改为k0=[-6,2]' options=optimset('TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',300, 'Algorithm','trust-region-reflective', 'display', 'iter'); lb=[-13 1];ub=[-1 5]; kp=lsqcurvefit(@num_value,k0,xdata,ydata,lb,ub,options); plot(xdata,ydata,xdata,num_value(kp,xdata)) %kp =[-3.5400,3.23000]; 方法2 multistart和lsqnonglin,全局算法（速度比较慢） global v; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); v=xydata(:,3); F =@(k)num_value(k,xdata)-ydata; ms=MultiStart('Display','iter', 'UseParallel','always');%开启并行运算matlabpool open 2; opts=optimset('Algorithm','trust-region-reflective'); problem=createOptimProblem('lsqnonlin','x0',[-4,2],'objective',F,'lb', [-6 1],'ub',[-1 5],'options',opts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) matlabpool close opts=optimset('Algorithm','trust-region-reflective','MaxTime',200，'Maxiter',400); 方法3 lsqcurvefit 和multistart global v; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); v=xydata(:,3); ms=MultiStart('Display','iter','UseParallel','always'); opts=optimset('Algorithm','trust-region-reflective'); matlabpool open problem=createOptimProblem('lsqcurvefit','x0',[-4,2], 'objective',@num_value,'xdata',xdata,'ydata',ydata,'lb',[-6,1],'ub',[ -1,5],'options',opts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) 方法4 globalsearch global v; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); v=xydata(:,3); F=@(k)norm(num_value(k,xdata)-ydata); gs=GlobalSearch('Display','iter'); matlabpool open opts=optimset('Algorithm','interior-point', 'UseParallel','always'); problem=createOptimProblem('fmincon','x0',[-4,2],'objective',F,'lb', [-6 1],'ub',[-1 5],'options',opts); [xming,fming,flagg,outptg,manyminsg] = run(gs,problem) 方法5 模拟退火 tic; global v; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); v=xydata(:,3); F=@(k)norm(num_value(k,xdata)-ydata); X0 =[-4,2];% Starting point lb=[-6 1];ub=[-1 5]; options= saoptimset('Display','iter'); [x,fval,exitFlag,output] = simulannealbnd(F,X0,lb,ub,options) time=toc; 二、关于求解变参数微分方程 看到有不少人问过如何处理"变参数微分方程", 所以抽了一点时间,写了一个例子, 希望能帮助到那些需要求解此类问题的人。 %%%================================%%% clear all fun=inline('[y(2);sin(w*t)-2*y(1)-3*y(2)]','t','y','flag','w'); tsp=[0 10]; y0=[1 1]; xlim(tsp) for w=1:10 [t,y]=ode45(fun,tsp,y0,[],w); plot(t,y) title(['w = ',num2str(w)]); pause(0.5); end 三、关于求解变上限积分问题 经常看到有人询问如何求解“变上限积分问题”，现举一个例子，说明该如何求解之：希望以后碰到类似问题时，能够自己动手解决． %%%===============================%%% clear all warning off x=linspace(0,150); for k=1:length(x) xx=x(k); fun=inline('exp(-lam.^2)'); w(k)=0.62*2/sqrt(pi)*quadl(fun,0,sqrt(pi)*xx/150); end plot(x,w) 5.1关于带参数的积分问题 例如以下问题: 函数为y=sin(k.*x).*x.^2，对x积分， 积分区域为【1，5】，目的是要画k 和y 的图形. 方法一（用循环）： 1.f=@(k) quad(@(x)sin(k.*x).*x.^2,0,5) 2.kk=linspace(0,5); 3.y=zeros(size(kk)); 4.for ii=1:length(kk) 5.y(ii)=f(kk(ii)); 6.end 7.plot(kk,y) 复制代码 方法二（不用循环）： 1.plot(linspace(0,5),arrayfun(@(k) quad(@(x) sin(k.*x).*x.^2,0,5),linspace(0,5))) 复制代码 三维参数拟合 function F=myfunf(k,x,y,z) %构造待拟合的函数模型，k为待回归系数 F=k(1)+k(2)*x+k(3)*y+k(4)*x.^2+k(5)*x.*y+k(6)*y.^2+k(7)*x.^3+k(8)*x.^2.*y-z; clear x=[1000:250:6000]'; y=[7.69,15.38,23.08,30.77,38.46,46.15,53.85,61.54,69.23,76.92,84.62,92.31,100]'; z= [9.1 10 11 11.4 11.4 11.1 12 12.3 12.7 12.6 13.3 13 13.6 13.6 13.7 13 12.7 11.9 11.6 11.2 10.8 13.9 14.8 15.7 15.2 16.1 17.1 17.1 17.3 16.74 17.4 17.6 18.4 20.9 20.4 21.7 21.9 21.3 20.5 19.7 18.9 18.3 19.6 20.1 21.7 24.1 25.5 25.6 25.9 25.7 25.3 26.7 26.7 26.6 27.3 27.9 29.3 28 27.5 26.8 25.8 24.9 23.9 26 26.4 28.6 28 28.5 28.7 28.7 28.9 29.3 29.4 29.4 32.4 34.3 35.2 36.4 35.2 34.2 33.4 32.3 31.4 30.2 30.4 31.4 32.2 33.2 33.6 34.9 35.8 35.5 35.6 36.4 36.3 40.9 41.6 41.1 44.3 43.6 42.3 41.5 40.6 39.3 38.7 36.5 38.4 41.1 41.9 43.1 43.1 42.6 42.1 42.1 43.8 44.3 48.3 50.4 47.7 51.2 51.3 50.5 49.6 48.8 47.7 47.1 42.8 46.3 47.8 47 49.1 50 52.8 53 53.3 55.3 55.2 55.7 57.5 57.8 58.9 60.1 59.4 58.7 57.8 56.8 55.6 48.1 49.5 54 53.8 55.3 55.8 55 54 54.2 56.8 57.3 63 63.6 62.4 64.4 66.1 65.2 64.3 63.4 62.5 61.3 54 56.8 59.4 59.4 60.8 61.2 62.9 62.5 63.9 64.4 64.5 68.9 69.6 69.4 70 72 72.7 71.2 69.8 68.1 66.9 57.8 61.4 65.6 65 67.6 70.7 69.8 68.2 68.9 68.3 69.8 75.9 76.7 76.4 76.2 76.8 76.5 75.8 73.5 71.6 70.4 63.1 65.6 68.6 69.4 73 73.4 74 74.4 74.2 75 75.7 82.7 83.2 83 82.5 85 85.7 83.8 81.6 80.1 79.4 67.8 70.9 73.2 74.8 79.2 80.9 81 81.7 81 80.8 83.7 89.9 90.5 88.9 90 90.9 90.3 89.1 87.3 85.5 84.1 74.6 75.6 83.5 83.4 87.6 90.7 91.2 92.6 94.4 96.8 99.8 100 101 100.5 100 100.4 99.8 98.2 96.74 94.5 92.307]; [x2,y2]=meshgrid(x,y); t11=interp2(x,y,z,x2,y2); k0=ones(1,8); lb=-1e5*ones(1,8); ub=1e5*ones(1,8); options=optimset('TolFun',1e-30,'TolX',1e-30,'MaxFunEvals',500 ,'display','iter'); [k,fminm]=lsqnonlin(@myfunf,k0,lb,ub,options,x2,y2,t11) %k=[0.1208 0.0097 0.0722 -0.0000 0.0003 0.0026 0.0000 -0.0000] %fminm=4.2986e+003 figure(1) surf(x2,y2,t11) % 画出原始数据三维图 figure(2) mesh(x2,y2,myfunf(k,x2,y2,t11)) %surf(x2,y2,myfunf(k,x2,y2,t11)) %画出误差曲面图 原始数据图 误差平面图 误差平面图 方法2全局搜索算法 function F=myfunf2(k,x,y) %另一种方法 F=k(1)+k(2)*x+k(3)*y+k(4)*x.^2+k(5)*x.*y+k(6)*y.^2+k(7)*x.^3+k(8)*x.^ 2.*y; clear x=[1000:250:6000]'; y=[7.69,15.38,23.08,30.77,38.46,46.15,53.85,61.54,69.23,76.92,84.62,92.31,100]'; z= [9.1 10 11 11.4 11.4 11.1 12 12.3 12.7 12.6 13.3 13 13.6 13.6 13.7 13 12.7 11.9 11.6 11.2 10.8 13.9 14.8 15.7 15.2 16.1 17.1 17.1 17.3 16.74 17.4 17.6 18.4 20.9 20.4 21.7 21.9 21.3 20.5 19.7 18.9 18.3 19.6 20.1 21.7 24.1 25.5 25.6 25.9 25.7 25.3 26.7 26.7 26.6 27.3 27.9 29.3 28 27.5 26.8 25.8 24.9 23.9 26 26.4 28.6 28 28.5 28.7 28.7 28.9 29.3 29.4 29.4 32.4 34.3 35.2 36.4 35.2 34.2 33.4 32.3 31.4 30.2 30.4 31.4 32.2 33.2 33.6 34.9 35.8 35.5 35.6 36.4 36.3 40.9 41.6 41.1 44.3 43.6 42.3 41.5 40.6 39.3 38.7 36.5 38.4 41.1 41.9 43.1 43.1 42.6 42.1 42.1 43.8 44.3 48.3 50.4 47.7 51.2 51.3 50.5 49.6 48.8 47.7 47.1 42.8 46.3 47.8 47 49.1 50 52.8 53 53.3 55.3 55.2 55.7 57.5 57.8 58.9 60.1 59.4 58.7 57.8 56.8 55.6 48.1 49.5 54 53.8 55.3 55.8 55 54 54.2 56.8 57.3 63 63.6 62.4 64.4 66.1 65.2 64.3 63.4 62.5 61.3 54 56.8 59.4 59.4 60.8 61.2 62.9 62.5 63.9 64.4 64.5 68.9 69.6 69.4 70 72 72.7 71.2 69.8 68.1 66.9 57.8 61.4 65.6 65 67.6 70.7 69.8 68.2 68.9 68.3 69.8 75.9 76.7 76.4 76.2 76.8 76.5 75.8 73.5 71.6 70.4 63.1 65.6 68.6 69.4 73 73.4 74 74.4 74.2 75 75.7 82.7 83.2 83 82.5 85 85.7 83.8 81.6 80.1 79.4 67.8 70.9 73.2 74.8 79.2 80.9 81 81.7 81 80.8 83.7 89.9 90.5 88.9 90 90.9 90.3 89.1 87.3 85.5 84.1 74.6 75.6 83.5 83.4 87.6 90.7 91.2 92.6 94.4 96.8 99.8 100 101 100.5 100 100.4 99.8 98.2 96.74 94.5 92.307]; [x2,y2]=meshgrid(x,y); t11=interp2(x,y,z,x2,y2); k0=ones(1,8); lb=-1e5*ones(1,8); ub=1e5*ones(1,8); F=@(k)myfunf2(k,x2,y2)-t11 ms=MultiStart; opts=optimset('Algorithm','trust-region-reflective','Display','iter'); problem=createOptimProblem('lsqnonlin','x0',k0,'objective',F,'lb',lb,'ub',ub,'options',o pts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) %xminm=[ 1.3571 -0.0006 0.6744 0.0000 0.0001 -0.0002 -0.0000 -0.0000] mesh(x2,y2,F(xminm)) 误差曲面图 fminm =1.0797e+003;%误差平方和 总结，方法1 没有搜索到全局最优值。方法2得到的最小平方误差和为方法1的1/3倍。 %fminm=4.2986e+003 % fminm =1.0797e+003; 最新加入，微分方程组参数识别 function dy=model1(t,y,k) dy=zeros(2,1); dy(1) =k(1)*y(2)+k(2); dy(2)=k(3)*y(2)-y(1)*k(4); function dy=model2(t,y) dy=zeros(2,1); dy(1) =1*y(2)+0.2.; dy(2)=0.5*y(2)-y(1)*k(4); clear; y0=[1,1]; t1=0:0.1:5; [t y]=ode45(@model2,t1,y0); clear; k=[ 1 0.2 0.5 0.7 ]; y0=[1,1]; t1=0:0.1:5; yy=y0; for i=1:length(t1)-1 [t y]=ode45(@model1,[t1(i) t1(i+1)],y0,[],k); y0=y(end,: ); yy=[yy;y0]; end xydata=[t1',yy]; save('xydata.txt','xydata','-ascii');%把tt和yy存为txt文件。 function yi=lsq33(k,x) y0=[1,1]; yy=y0; for i=1:length(x)-1 [t y]=ode45(@model1,[x(i) x(i+1)],y0,[] ,k); y0=y(end,: ); yy=[yy;y0]; end yi=yy(:,1); %只输出y的第一列 clear; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); x0=[1,1,1,1]; lb=[0, 0,0,0]; ub=[3,3,3,3]; options=optimset('TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',300,'Algorit hm','trust-region-reflective', 'display', 'iter'); kp=lsqcurvefit(@lsq33,x0,xdata,ydata,lb,ub,options) plot(xdata,ydata,xdata,lsq33(kp,xdata)) clear load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); F =@(k)lsq33(k,xdata)-ydata; x0=[1,1,1,1]; lb=[0, 0,0,0]; ub=[3,3,3,3]; opts=optimset('Algorithm','trust-region-reflective'); problem=createOptimProblem('lsqnonlin','x0',x0,'objective',F,'lb',lb, 'ub',ub,'options',opts); [xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,100) %参数识别，假设只能测量到x和y(1)。 clear; load xydata.txt; xdata=xydata(:,1); ydata=xydata(:,2); k0=[1,1,1,1]'; %初值也可以改为k0=[-6,2]' ms=MultiStart('Display','iter'); options=optimset('TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',300, 'Algorithm','trust-region-reflective', 'display', 'iter'); lb=[-13 1];ub=[-1 5]; kp=lsqcurvefit(@lsq33,k0,xdata,ydata,lb,ub,options); plot(xdata,ydata,xdata,num_value(kp,xdata)) %kp =[-3.5400,3.23000]; MatLab向量化编程:arrayfun函数的使用 自己的实验； 1 t=0:0.1:5; T=arrayfun(@(x) x.^2,t); 0-5的挨个平方 2 t=0:0.1:5; f=@(x)x.^2; f(t);同上 3 t=0:0.1:5; f=@(x)x.^2; y=arrayfun(f,t);同上 1 求f(k)=int(sin(k*x)*x^2,0,5,x),k取[0,5]上的不同值； 方法1 匿名函数 tic; f=@(k)quadl(@(x)sin(k.*x).*x.^2,0,5);%相当于多重匿名函数呵呵。 kk=linspace(0,50;%默认分为100等分 y2=zeros(size(kk)); for ii=1:length(kk) y2(ii)=f(kk(ii)); end time=toc; 方法2：向量化方法+匿名函数 tic ; y=arrayfun(@(k) quadl(@(x) sin(k*x).*x.^2,0,5),linspace(0,5));%相当于求一百组循环代入k。 time=toc; 例子1 匿名函数数值微分方程 f=@(t,x)[x(2); -t*x(1)+exp(t)*x(2)+3*sin(2*t)]; tspan=[3.9 4.0]; %求解区间 y0=[2 8]; %初值 [t,x]=ode45(f,tspan,y0); plot(t,x(:,1),'-o',t,x(:,2),'-*') ; legend('y1','y2'); title('y'' ''=-t*y + e^t*y'' +3sin2t') ; xlabel('t') ;ylabel('y'); 最小二乘参数估计 摘要： 最小二乘的一次性完成辨识算法（也称批处理算法），他的特点是直接利用已经获得的所有（一批）观测数据进行运算处理。这种算法在使用时，占用内存大，离线辨识，观测被辨识对象获得的新数据往往是逐次补充到观测数据集合中去的。在应用一次完成算法时，如果要求在每次新增观测数据后，接着就估计出系统模型的参数，则需要每次新增数据后要重新求解矩阵方程()Z l T l l T l ΦΦΦ-∧=1θ。 最小二乘辩识方法在系统辩识领域中先应用上已相当普及，方法上相当完善，可以有效的用于系统的状态估计，参数估计以及自适应控制及其他方面。 关键词： 最小二乘（Least-squares ），系统辨识（System Identification ） 目录： 1.目的 (1) 2.设备 (1) 3引言 (1) 3.1 课题背景 (1) 4数学模型的结构辨识 (2) 5 程序 (3) 5.1 M 序列子函数 ................................................................................. 错误！未定义书签。 5.2主程序............................................................................................... 错误！未定义书签。 6实验结果： ................................................................................................................................... 3 7参考文献： ................................................................................................. 错误！未定义书签。 1.目的 1.1掌握系统辨识的理论、方法及应用 1.2熟练Matlab 下最小二乘法编程 1.3掌握M 序列产生方法 2.设备 PC 机1台（含Matlab 软件） 3引言 3.1 课题背景 最小二乘理论是有高斯（K.F.Gauss ）在1795年提出：“未知量的最大可能值是这样一个数值，它使各次实际观测值和计算值之间的差值的平方乘以度量其精度的数值以后的和最小。”这就是最小二乘法的最早思想。 最小二乘辨识方法提供一个估算方法，使之能得到一个在最小方差意义上与实验数据最 %ITD法识别模态参数 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=');%指定计算模态阶数 Matlab系统辨识尝试之详细过程1 前面介绍了Matlab系统辨识工具箱的一些用法，这里拿一个直观的例子来尝试工具箱的具体用法。比较长，给个简单目录吧： 1.辨识的准备 2.辨识数据结构的构造 3.GUI辨识 4.辨识效果 5.对固有频率的辨识 6.结构化辨识 7.灰箱辨识 8.加入kalman滤波的灰箱辨识 1.辨识的准备 在辨识前，首先要根据自己辨识的情况，确定要辨识的状态空间模型的一些特点，如连续还是离散的；有无直通 分量（即从输入直通到输出的分量）；输入延迟；初始状态等。了解了这些情况就可以更快速的配置辨识时的一些设 置选项。 2.辨识数据结构的构造 使用原始数据构造iddata结构： data=iddata(y,u,Ts); 这里以一个弹簧质量系统的仿真为例 代码如下，其中用到了函数MDOFSolve，这在之前的博文介绍过（https://www.wendangku.net/doc/8d14329317.html,/?p=183），拿来用即可。如果发现运行有错误，可以将MDOFSolve函数开头的一句 omega2=real(eval(omega2)); 注释掉。 %弹簧质量系统建模 clc clear close all m=200; k=980*1000; c=1.5*1000; m1=1*m; m2=1.5*m; k1=1*k; k2=2*k; k3=k1; %%由振动力学知识求固有频率 M=[m10;0m2]; K=[k1+k2-k2;-k2k3+k2]; [omega,phi,phin]=MDOFSolve(M,K); fprintf('固有频率:%fHz\n',subs(omega/2/pi)); %%转化到状态空间 innum=2; outnum=2; statenum=4; A=[0100; -(k1+k2)/m10k2/m10; 0001; k2/m20-(k3+k2)/m20]; B=[00; 1/m10; 00; 01/m2]; C=[1000; 0010]; D=zeros(outnum,innum); K=zeros(statenum,innum); mcon=idss(A,B,C,D,K,'Ts',0);%连续时间模型 figure impulse(mcon) %%信号仿真，构造数据供辨识 n=511;%输入信号长度 Ts=0.001; t=0:Ts:(n-1)*Ts; u1=idinput(n,'prbs');%输入1为伪随机信号 u2=zeros(n,1);%输入2为空 u=[u1u2]; simdat=iddata([],u,Ts);%形成输入数据对象 e=randn(n,2)*1e-7; simopt=simOptions('AddNoise',true,'NoiseData',e);%添加噪声yn=sim(mcon,simdat,simopt);%加噪声仿真 y=sim(mcon,simdat);%无噪声仿真 模态参数识别方法的比较研究 发表时间：2017-09-07T14:07:39.937Z 来源：《防护工程》2017年第9期作者：安鹏强[导读] 本文将频域法、时域法和整体识别法识别模态参数的应用范围、存在的优缺点进行对比、分析和说明。 航天长征化学工程股份有限公司兰州分公司甘肃兰州 730050 摘要：本文将频域法、时域法和整体识别法识别模态参数的应用范围、存在的优缺点进行对比、分析和说明，对模态参数识别的研究方向具有指导意义。 关键词：模态参数识别；频域法；时域法；整体识别法 引言 多自由度线性振动系统的微分方程可以表达为[1]： [M]{x ?(t)}+[C]{x ?(t)}+[K]{x(t)}={f(t)} 通过将试验采集的系统输入与输出信号用于参数识别的方法中，进而对系统的模态质量、模态阻尼、模态刚度、模态固有频率及模态振型进行识别，这一过程称为结构的模态参数识别。本文将对模态参数识别的频域法、时域法及整体识别法三者的应用范围、存在的优缺点进行对比、分析和说明。 1频域法 模态参数识别的频域法是结合傅里叶变换理论[1]形成的，这种方法是从实测数据的频响函数曲线上对测试结构的模态参数进行估计。图解法[1]是最早的频域模态参数识别方法，随之，又陆续发展了导纳圆拟合法[2]、最小二乘迭代法[2]、有理式多项式法[2]等多种频域模态参数识别方法。 频域法的优点是直观、简便，噪声影响小，模态定阶问题易于解决。频域法识别模态参数的思路是首先借助实测频响函数曲线对模态参数进行粗略的估计，进而将初步观测的模态估计值作为一些频域识别法的最初输入值，通过反复的迭代获取最终的模态参数。频域识别方法对于实测频响函数的分布容易控制，其输人数据是主观人为的。频域中参数识别方法识别结果的精准度，取决于测试试验中获得的频响函数质量的好坏。判断实测频响函数的质量，就要看其曲线的光滑[2]和曲线的饱满程度[2]，曲线越光滑越饱满的实测频响函数，用其进行参数识别时，识别精度越高。 2时域法 模态参数识别的时域法的研究与应用比频域法晚，时域法可以克服频域法的一些缺陷。时域模态参数识别的技术优点在于无需获得激励力即可进行参数的识别[3-7]。对于一些大型的工程结构如大坝、桥梁等，获取激励荷载不太容易，但容易测得他们在风、地脉动等环境激励下的响应数据，把这些响应数据用于时域中一些参数识别的方法上，即可对测试结构的模态参数进行识别。 时域法的优点不仅在于其无需激励设备、减少测试费用而且可以避免由信号截断而造成对识别精度的影响，并且可实现对大型工程结构的在线参数识别，真实地反映结构的动力特性。但是由于响应信号中含有大量的噪声，这会使得所识别的模态中含有虚假模态。目前，对于如何剔除噪声模态、优化识别过程中的一些参数问题、以及怎样更稳定、可靠地进行模态定阶等成为时域法研究中的重要课题。目前常用的判定模态真假的方法是稳定图方法[8]，该方法的基本思想在于不同阶次的系统模型会对虚假模态的影响比较大，在稳定图中出现次数最多的模态可认为是系统的真实模态。 3整体识别法 结构模态参数识别的单输入单输出类型是针对单个响应点的数据进行相应的计算，从而得到该测点对应的模态频率、阻尼比和振型系数等动力参数，但是对于有多个测点的试验，若要用单输入单输出类型的识别方法对多自由度结构进行参数识别，则需要对各个测点单独计算来识别各个测点对应的模态参数，通过对各个测点分别计算处理，得到每一个测点数据所识别的模态参数，然后求取所有测点响应识别的算术平均值来作为整体结构最终的识别结果。理论上讲，用每个测点数据识别的结果应该是一样的，但实际测试实验中，因测试实验中测点布置位置的不同、测试中其他因素及识别方法上的不完善会使得各个测点的识别结果不同、识别精度不同及错误的识别结果等现象。因此，对于多测点的测试试验，用单输入单输出类型的识别方法进行参数识别不仅会因多次重复导致计算工作量复杂累赘而且识别结果的正确性及精度无法保证。 整体识别的方法避免了单输入单输出类型的一些不足之处。该方法通过将结构上的所有测点的实测数据同时进行识别计算，所识别得到的结果作为结构整体的模态参数，每阶模态的固有频率和阻尼比是唯一的，减小了随机误差，提高了识别进度，并且使得计算工作量大大减少。 4三种识别方法的比较分析 （1）频域内的模态参数识别方法方便、快捷，但在实际运用中人为的主观选择性对识别结果的影响较大； （2）基于环境激励的时域模态参数的识别方法具有测试试验的花费较少、测试相对安全，并且识别精度较高。因此，基于环境激励的时域模态参数的识别方法已成为科研工作者研究的热点问题。 （3）对于多测点的测试试验，用频域和时域的单输入单输出类型识别模态参数不仅会因多次重复导致计算工作量复杂累赘而且识别结果的正确性及精度无法保证。整体识别法将所有测点的数据同时进行处理计算，得到结构的整体识别结果。整体识别方法通过对所有测点数据同时进行识别计算，减小了随机误差，提高了识别进度，使得计算工作量大大减少。 （4）对比时域和频域识别方法对虚假模态的剔除，可以看出，频域中的剔除虚假模态主要依据模态频率在频幅曲线图上会出现峰值的原理，利用该峰值处的幅值角是否为0°或180°来剔除虚假模态；相对频域剔除虚假模态的方法来说，时域中的剔除虚假模态的方法有定量的精度判别指标。总体看来，时域识别方法无法判别是否已将系统的所有模态进行识别且对于阻尼比的确定还有待研究。参考文献 [1] 曹树谦，张德文，萧龙翔. 振动结构模态分析－理论、实验与应用[M]. 天津大学出版社，2001． [2] 王济，胡晓. Matlab在振动信号处理中的应用[M]. 水利水电出版社，2006. 例1、考虑仿真对象 )()2(5.0)1()2(7.0)1(5.1)(k v k u k u k z k z k z +-+-=-+-- 其中，)(k v 是服从正态分布的白噪声N )1,0(。输入信号采用4阶M 序列，幅度为1。选择如下形式的辨识模型 )()2()1()2()1()(2121k v k u b k u b k z a k z a k z +-+-=-+-+ 设输入信号的取值是从k =1到k =16的M 序列，则待辨识参数LS θ?为LS θ?=(T T -ΦΦΦ1)z 。其中，被辨识参数LS θ?、观测矩阵Φ的表达式为: ????? ???????=2121?b b a a LS θ (3)(4)(16)z z z ??????=???? ??z (2)(1)(2)((3)(2)(3)(2)(15)(14)(15)(14)z z u u z z u u z z u u --????--??Φ=????--?? 程序框图如图1所示。Matlab 仿真程序如下： %二阶系统的最小二乘一次完成算法辨识程序，文件名：LS.m u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画输入信号u的径线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线 subplot(3,1,3) %画三行一列图形窗口中的第三个图形 stem(z),grid on %画出输出观测值z的径线图形，并显示坐标网格u,z %显示输入信号和输出观测信号 %L=14 %数据长度 HL=[-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9) -z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14)] %给样本矩阵 赋值 模态参数辨识的频域方法 吕毅宁 目录 模态参数辨识的频域方法 (1) 单点输入单点输出(SISO) (1) 图解法............................................................................................................ 1 频域多参考点模态参数辨识（MIMO ） ............................................................ 2 频域模态测试和参数辨识的可控性和可观性. (5) 单点输入单点输出(SISO) 图解法 1) 峰值检测 半功率点 )(2 1 )()(21r j H j H j H ωωω= = （1） r r ωωωξ21 2-= （2） 2) 模态检测 () ir r jr r r r ij r jr ir r r r r jr ir r r ij Q A Q j j Q j H ψσψσσψψωσωψψω-= -= -= +-= ) ()( （3） 式中，r Q 是模态比例换算因子。 在上式中，() r ij A 是模态质量r m 和模态刚度r k 的函数，又由下面的关系 2r r r m k ω= （4） 联立即可求得模态质量和模态刚度。 3) 圆拟合法 固有频率 max ==ω ωωd ds r r （5） 振型 r er I ij g k H 1 -= （6） jr ir r er k k ??= （7） er k 是等效模态刚度，r r r k g η= 是等效结构阻尼。 ()r ij r I ij ir r r jr R g k )(2==-H ?? （8） 模态阻尼 r g ) 1(2tan 211 ωα-= （9） r g ) 1(2tan 222 -= ωα （10） 2 tan 2 tan 22 1 12ωωω+-= r r g （11） 模态刚度 由 r er r I ij g k H 1 )1(-= =ω （12） 可得 r r I ij er g H k )1(1 =-= ω （13） 模态质量 2 r r r k m ω= （14） 其他方法，如正交多项式曲线拟合法，非线性优化辨识方法。 频域多参考点模态参数辨识（MIMO ） 一个N 自由度粘性阻尼线性系统，对它施加P 个激励力，在N 个点上进行响应 一．模态参数识别的单模态法 常见的单模态识别有三种方法：直接读数法（分量分析法）、最小二乘圆拟合法和差分法。 所谓单模态识别法，是指一次只识别一阶模态的模态参数，所用数据为该阶模态共振频率附近的频响函数值。待识别的这阶模态称为主导模态，余模态称为剩余模态，剩余模态的影响可以全部忽略或简化处理。 1． 直接读数法（分量分析法） 1）基本公式 所谓分量分析法就是讲频响函数分成实部分量和虚部分量来进行分析。 N 自由度结构系统结构，p 点激励l 点响应的实模态频响函数可表示如下： 2222222111 ()(1)(1)N r r lp r er r r r r g H j K g g ωωωω=??--=+??-+-+?? ∑ (1.1) 其中r er lr pr K K φφ= ，为第二阶等效刚度 /r r ωωω= g 2r r r ζω= ，为第r 阶模态结构阻尼比 当ω趋近于某阶模态的固有频率时，该模态起主导作用，称为主导模态或者主模 态。 在主模态附近，其他模态影响较小。若模态密度不是很大，各阶模态比较远离，其余模态的频响函数值在该模态附近很小，且曲线比较平坦，即几乎不随频率而变化，因此其余模态的影响可以用一个复常数来表示，第r 阶模态附近可用剩余模态表示成： 222222211 ()()(1)(1)R I r r lp C C er r r r r g H j H H K g g ωωωω??-= -++??-+-+?? (1.2) ()lp H ω的实部和虚部可分别表示如下： 222211 ()(1)R R r lp C er r r H H K g ωωω??-= +??-+?? (1.3) 2221 ()(1)I I r lp C er r r g H H K g ωω??-= +??-+?? (1.4) 【1】随机序列产生程序 【2】白噪声产生程序 【3】M序列产生程序 【4】二阶系统一次性完成最小二乘辨识程序 【5】实际压力系统的最小二乘辨识程序 【6】递推的最小二乘辨识程序 【7】增广的最小二乘辨识程序 【8】梯度校正的最小二乘辨识程序 【9】递推的极大似然辨识程序 【10】Bayes辨识程序 【11】改进的神经网络MBP算法对噪声系统辨识程序【12】多维非线性函数辨识程序的Matlab程序【13】模糊神经网络解耦Matlab程序 【14】F-检验法部分程序 【1】随机序列产生程序 A=6; x0=1;M=255; for k=1:100 x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=v1; x0=x1; v0=v1; end v2=v k1=k; %grapher k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');title('(0,1)均匀分布的随机序列') 【2】白噪声产生程序 A=6; x0=1; M=255; f=2; N=100; for k=1:N x2=A*x0; x1=mod (x2,M); v1=x1/256; v(:,k)=(v1-0.5)*f; x0=x1; v0=v1; end v2=v k1=k; %grapher k=1:k1; plot(k,v,k,v,'r'); xlabel('k'), ylabel('v');title('(-1,+1)均匀分布的白噪声') 【3】M序列产生程序 X1=1;X2=0;X3=1;X4=0; %移位寄存器输入Xi初T态（0101），Yi为移位寄存器各级输出m=60; %置M序列总长度 for i=1:m %1# Y4=X4; Y3=X3; Y2=X2; Y1=X1; X4=Y3; X3=Y2; X2=Y1; X1=xor(Y3,Y4); %异或运算 if Y4==0 U(i)=-1; else U(i)=Y4; end end M=U %绘图 i1=i k=1:1:i1; plot(k,U,k,U,'rx') xlabel('k') ylabel('M序列') title('移位寄存器产生的M序列') 【4】二阶系统一次性完成最小二乘辨识程序 %FLch3LSeg1 u=[-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1]; %系统辨识的输入信号为一个周期的M序列 z=zeros(1,16); %定义输出观测值的长度 for k=3:16 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想输出值作为观测值 end subplot(3,1,1) %画三行一列图形窗口中的第一个图形 stem(u) %画出输入信号u的经线图形 subplot(3,1,2) %画三行一列图形窗口中的第二个图形 i=1:1:16; %横坐标范围是1到16，步长为1 plot(i,z) %图形的横坐标是采样时刻i, 纵坐标是输出观测值z, 图形格式为连续曲线 模态分析方法在发动机曲轴上的应用研究 xx （xx大学 xxxxxxxx学院 , 山西太原 030051） 摘要：综述模态分析在研究结构动力特性中的应用，介绍模态分析的两大方法：数值模态分析与试验模态分析。并着重介绍目前的研究热点一一工作模态分析。通过发动机曲轴的模态分析这一具体的实例，综述了运行模态分析国内外研究现状，指出了其关键技术、存在问题以及研究发展方向。 关键词：模态分析数值模态试验模态工作模态 Abstract ：Sums up methods of model analysis applied on the research of configuration dynamic;al characteristio. It introduces two methods of model analysis: numerical value model analysis and experimentation model analysis. Then it stresses the hotspot-working model analysis.Some key techniques, unsolved problems and research directions of OMA were also discussed. Key words:Model analysis Numerical value model analysis Experimentation model analysis Working model analysis 1、引言 1.1模态分析的基本概念 物体按照某一阶固有频率振动时，物体上各个点偏离平衡位置的位移是满足一定的比例关系的，可以用一个向量表示，这个就称之为模态。模态这个概念一般是在振动领域所用，你可以初步的理解为振动状态，我们都知道每个物体都具有自己的固有频率，在外力的激励作用下，物体会表现出不同的振动特性。 一阶模态是外力的激励频率与物体固有频率相等的时候出现的，此时物体的振动形态叫做一阶振型或主振型；二阶模态是外力的激励频率是物体固有频率的两倍时候出现，此时的振动外形叫做二阶振型，以依次类推。 各种模态分析方法总结与比较 一、模态分析 模态分析是计算或试验分析固有频率、阻尼比和模态振型这些模态参数的过程。 模态分析的理论经典定义：将线性定常系统振动微分方程组中的物理坐标变换为模态坐标，使方程组解耦，成为一组以模态坐标及模态参数描述的独立方程，以便求出系统的模态参数。坐标变换的变换矩阵为模态矩阵，其每列为模态振型。 模态分析是研究结构动力特性一种近代方法，是系统辨别方法在工程振动领域中的应用。模态是机械结构的固有振动特性，每一个模态具有特定的固有频率、阻尼比和模态振型。这些模态参数可以由计算或试验分析取得，这样一个计算或试验分析过程称为模态分析。这个分析过程如果是由有限元计算的方法取得的，则称为计算模记分析；如果通过试验将采集的系统输入与输出信号经过参数识别获得模态参数，称为试验模态分析。通常，模态分析都是指试验模态分析。振动模态是弹性结构的固有的、整体的特性。如果通过模态分析方法搞清楚了结构物在某一易受影响的频率围各阶主要模态的特性，就可能预言结构在此频段在外部或部各种振源作用下实际振动响应。因此，模态分析是结构动态设计及设备的故障诊断的重要方法。 模态分析最终目标是在识别出系统的模态参数，为结构系统的振动特性分析、振动故障诊断和预报以及结构动力特性的优化设计提供依据。 二、各模态分析方法的总结 （一）单自由度法 一般来说，一个系统的动态响应是它的若干阶模态振型的叠加。但是如果假定在给定的频带只有一个模态是重要的，那么该模态的参数可以单独确定。以这个假定为根据的模态参数识别方法叫做单自由度(SDOF)法n1。在给定的频带围，结构的动态特性的时域表达表示近似为： ()[]}{}{T R R t r Q e t h r ψψλ= 2-1 而频域表示则近似为： ()[]}}{ {()[]2ωλωψψωLR UR j Q j h r t r r r -+-= 2-2 单自由度系统是一种很快速的方法，几乎不需要什么计算时间和计算机存。 这种单自由度的假定只有当系统的各阶模态能够很好解耦时才是正确的。然而实际情况通常并不是这样的，所以就需要用包含若干模态的模型对测得的数据进行近似，同时识别这些参数的模态，就是所谓的多自由度(MDOF)法。 单自由度算法运算速度很快，几乎不需要什么计算和计算机存，因此在当前小型二通道或四通道傅立叶分析仪中，都把这种方法做成置选项。然而随着计算机的发展，存不断扩大，计算速度越来越快，在大多数实际应用中，单自由度方法已经让位给更加复杂的多自由度方法。 1、峰值检测 峰值检测是一种单自由度方法，它是频域中的模态模型为根据对系统极点进行局部估计(固有频率和阻尼)。峰值检测方法基于这样的事实：在固有频率附近，频响函数通过自己的极值，此时其实部为零(同相部分最 %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=');%指定计算模态阶数 设一非线性系统如下所示: ()()() ()()111y k y k e u k k βαε--=-+-+ 0.75α= 0.35β=0.25γ=,()k ε是零均值，方差为0.01的噪声序列（均匀白噪声）。 （1）试设计一种激励信号能持续激励系统的各工作点（平衡点） （2）用适当的方法辨识出系统的等价模型（用另一组数据来检验模型的泛化性） 说明：下面讨论的都是离散系统，所以时间坐标均采用离散时间节点k 。 解： (1) 线性化处理寻找系统的合理输入信号 可以求得系统的平衡点为： ()0.75y k α== (1.1) 按题意要求最后系统必须收敛于平衡点附近，即： ()lim 0.75k y k →∞ = (1.2) 为了找出系统的合理输入信号，使得系统最终工作在平衡点附近，这里可以将系统线性化处理，将上述非线性系统进行泰勒展开得： ()()()()()()()2323111111112!3!! n n y k y k y k y k y k u k k n αβαβαβαβγε--+---++-=-+ 因为 ，后面()1y k -的高阶项都可以扔掉（只作为寻找输入信号使用）， 所以系统可以化为下式： 此时不妨设系统输出()y k 的最后的极限为A , 从式1.2得0.75A =。 那么应该满足 ()()lim lim 1k k y k y k A →∞ →∞ =-= (1.5) 从而有 ()()()11A u k k αβγε-=-+ (1.6) 同时为了抵消系统的部分噪声，这里采用MA TLAB 软件编程产生另一服从同一分布的均匀噪声()1k ε，将式1.6变形得： ()()()0111 1u k u k k αβ εγ γ --= - (1.7) 式1.7中()0u k 是一个最后收敛于系统平衡点0.75的基本信号，这里可以采用一阶线性系统的阶跃响应曲线作为基本信号()0u k ，同时考虑系统的平衡点，即设计为： ()/00.75k T u k e -=- (1.8) T 是一阶线性系统对应的时间常数， 反应到输入基本信号()0u k 上就是过零点作()0u k 对应()()()()11y k y k u k k αβγε--=-+01αβ<（1.3） (1.4) 振动模态分析理论与应用 模态参数识别频域法 当系统阻尼为比例阻尼或小阻尼时，阻尼矩阵经模态坐标变换后可以对角化，模态参数为实数，频响函数可按实模态展开。若在p 点激励，在l 点测量，则频响函数可表示为对于粘性阻尼有 ∑ 1 2 ωω ξ2ωω1 )ω(N i i i i lp lp j D H =+= 对于结构阻尼有 ∑ 1 2ωω 1 )ω(N i i i lp lp jg D H =+= 以上两式即为实模态参数识别的基本公式 6.1 实模态识别图解法 6.1.1 共振法 这是一种经典的模态分析方法，其基本思想是：当激励频率在系统某阶固有频率r ω附近时， 该阶模态导纳便起主导作用，其余各阶模态导纳的影响可忽略不计。即 )ω(≈)ω(lpr lp H H 此时，整个系统等效于一个单自由度系统。利用幅频特性和相频特性，便可确定系统的模态参数（参看图6-1）。 在待测结构上选择l 个测试点，求其中某点P 对所有各点的位移导纳。点数l 一般应等于或大于拟选的模态数N （自由度数）。则p 点对任意点l 的位移导纳可作如下处理： 当激振频率在r 阶固有频率附近时有 () () 2 22 2∞ 1 2 ωωξ4ωω1≈ ωω ξ2ωω1 )ω(∑ ++==r r i r lp i i i i i lp lp j D j D H 因此，测得的幅频曲线)ω(lp H 的第r 个峰值位置（共振频率点），便可近似确定r 阶固有频率r ω。由r ω两侧半功率带宽，可以确定r 阶模态阻尼比)ω2/Δω(ξr r =。由r ω处位移 有 ()r r lp r lp D H ξ2)ω(= 所以 ()()r lp r r lp H D ωξ2= 由因为 ()r pr lr r lp k D φ φ= 故在令pr φ的值等于1（振型中各元素具有确定的比例，其绝对值可认为地指定，不妨取第r 阶振型第p 个元素pr φ的值等于1）时，由原点导纳曲线的峰值可得r 阶模态刚度为 ) ω(ξ21 r pp r r H k = 此外，当r ωω=时，l 个导纳的幅值分别为 r r pr r r p k H ξ2φφ|)ω(|11= r r pr r r p k H ξ2φφ|)ω(|22= r r pr lr r lp k H ξ2φφ|)ω(|= 写成矩阵形式 = lr r r r r pr r lp r p r p k H H H φφφξ2φ| )ω(|| )ω(||)ω(|2121 因此，第r 阶振型为 {}±±±==| )ω(||)ω(|| )ω(|φφ φφ2121r lp r p r p lr r r r H H H 为表示振型的几何形状，上试中各导纳幅值应考虑其相位，可用正负号表示同相或反相，对 于实模态，其振型向量的各分量都是实数，且只有大小和正负之差。因此，系统作固有振动时，各坐标点同时达到极值，同时通过平衡位置。用共振法确定模态参数，方法简单直观。但由于忽略了相邻模态的影响，识别出的模态精度不高，特别是识别振型和阻尼时，可能引起较大的误差。另外当各阶模态耦合较密时可能识别不出单个模态。因此这种方法一般只用于对模态的初步分析。 6.1.2分量分析法 分量分析法的思想是利用导纳的实频和虚频特性识别出系统的模态参数。其优点是能考虑其余模态的影响。 环境振动下模态参数识别方法综述 摘要：模态分析是研究结构动力特性的一种近代方法，是系统识别方法在工程振动领域中的应用。环境振动是一种天然的激励方式，环境振动下结构模态参数识别就是直接利用自然环境激励，仅根据系统的响应进行模态参数识别的方法。与传统模态识别方法相比，具有显著的优点。本文主要是做了环境振动下模态识别方法的一个综述报告。 关键词：环境振动模态识别综述 Abstract: The modal analysis is the study of structural dynamic characteristics of a modern method that is vibration system identification methods in engineering applications in the field. Ambient vibration is a natural way of incentives, under ambient vibration modal parameter identification is the direct use of the natural environment, incentives, based only on the response of the system for modal parameter identification method. With the traditional modal identification methods, has significant advantages. This paper is a summary report of the environmental vibration modal identification method. Keywords: Ambient vibration ；modal parameters ；Review 随着我国交通运输事业的发展，各种形式的大、中型桥梁不断涌现，由于大型桥梁结构具有结构尺大、造型复杂、不易人工激励、容易受到环境影响、自振频率较低等特点，传统模态参数识别技术在应用上的局限性越来越突出。传统的振动试验采用重振动器或落锤激励桥梁，需要投入大量人力和试验设备，激励成本增高，难度大，而且对于桥梁这样的大型复杂结构，激励(输入)往往很难测得，也不适合长期监测的实验模态分析。 环境振动是指振幅很小的环境地面运动。系由天然的和（或）人为的原因所造成，例如风、海浪、交通干扰或机械振动等，受激结构的振幅较小，但响应涵盖频率丰富。系统或者结构的模态参数包括：模态频率、模态阻尼、模态振型等。模态参数识别是系统识别的一部分，通过模态参数的识别可以了解系统或结构的动力学特性，这些动力特性可以作为结构有限元模型修正、故障诊断、结构实时监测的评定标准和基础。环境振动下的模态参数识别就是利用自然环境激励，根据结构的动力响应来进行模态参数识别的方法。 1 环境振动下模态参数识别的优点 传统的模态识别方法利用结构的输入和输出信号识别结构的模态参数。对于工作中的大型结构，无论是对其实施外部激励还是测试外部激励都十分困难。而环境振动方法仅仅利用被测试的输出数据识别结构的时间序列分析法模态参数。用环境振动对结构进行模态参数识别，具有明显的优点： M序列发生程序 function [out]=Mseries(order) %M-series function % by Chen Hao,2011/12/16 for i=1:order x(i)=1; end for n=1:252 y(1:2)=x(order-1:order); for i=order:-1:2 x(i)=x(i-1); end x(1)=xor(y(1),y(2)); out(n)=1-2*x(order); end stem(out); end 计算程序 function [out]=addfunction( A,B,order,n ) %addfunction % by Chen Hao,2011/12/9 order=n*order; for i=1:order out(i)=0; end for k=0:order-1 for i=0:order-1 out(k+1)=out(k+1)+A(i-k+order+1)*B(i+1); end end for k=1:order out(k)=(out(k)-out(order))/(order+1); end end 主程序 function main() %Main function % by Chen Hao,2011/12/9 out=Mseries(6); time=1:252; simin=[time' out']; assignin('base','simin',simin); sim('Mseries_distinguish'); z=simout.signals.values(2:253); out=[out out]; g=addfunction(out,z,63,4); g=g(1:63); plot(g); num=[1 2]; den=[3 1 1]; s=tf(num,den); t=0:62; h=impulse(s,t); 2015 年春季学期研究生课程考核 （读书报告、研究报告） 考核科目:实验模态分析 学生所在院（系）: 学生所在学科: 学生姓名: 学号: 学生类别: 考核结果阅卷人 实验模态分析与参数识别报告 模态分析可分为实验模态分析与工作模态分析等。模态分析的最终目标是识别出系统的模态参数，为结构系统的振动分析、振动故障诊断和预报、结构动力特性的优化设计提供依据。模态是机械结构的固有振动特性，每一个模态具有特定的固有频率、阻尼比和模态振型。这些模态参数可以由计算或试验分析取得，这样一个计算或试验分析过程称为模态分析。 1、模态分析原理 模态分析的过程是将线性时不变系统振动微分方程组中的物理坐标变换为模态坐标，使方程组解耦，成为一组以模态坐标及模态参数描述的独立方程，坐标变换的变换矩阵为振型矩阵，其每列即为各阶振型。 []{}[]{}[]{}{}()M X C X K X F t ++= （1） 其中：[]M —质量矩阵，[]K —刚度矩阵，[]C —粘性阻尼矩阵，{}()F t —激励力的列阵。 振动模态是弹性结构固有的、整体的特性。如果通过模态分析方法搞清楚了结构物在某一易受影响的频率范围内，各阶主要模态的特性，就可能预知结构在此频段内，在外部或内部各种振源作用下实际振动响应，而且一旦通过模态分析知道模态参数并给予验证，就可以把这些参数用于设计过程，优化系统动态特性，或者研究把该结构连接到其他结构上时所产生的影响。 方程（1）经傅氏变换，可得频域内的振动方程： [][][]{}{}2()()()M j C K X F w w w w -++= （2） 对应于固有频率1ω的固有振型或模态向量以幅值最大点为参考点的表达式为：{}{}11max 1()()X X w w w =。它们亦即简谐自由振动的主振型，满足以下关系式： [][]{}2()0i K M w j -= （3） 此代数方程组的系数行列式等于零，即为特征方程式；[]M ，[]K 为实数对称矩阵，[] M 正定，[]K 为非负定，其特征值20ω和对应的特征向量为实数。 主振型矩阵[]{}{}{}1 2,,,,n j j j j 轾=臌为实模态矩阵。根据振型的正交性： [][][][]1T M M j j =，[][][][]1T K K j j =；系统阻尼为比例阻尼时， [][][][]1T C C j j =。系统辨识最小二乘参数估计matlab
ITD模态参数识别matlab修改版
Matlab系统辨识尝试之详细过程1
模态参数识别方法的比较研究
Matlab_系统辨识_应用例子
模态参数辨识的频域方法
模态参数识别的单模态法,模态参数识别的导纳圆法
系统辨识及其matlab仿真(一些噪声和辨识算法)
模态分析与参数识别
各种模态分析方法总结与比较
ITD模态参数识别matlab修改版
系统辨识MATLAB仿真
模态参数识别频域法
环境振动下模态参数识别方法综述
M序列辨识 Matlab程序 系统辨识
实验模态分析与参数识别报告