遗传算法工具箱识别(GA)Bouc-Wen模型参数辨识_识别

Bouc-Wen模型因数字处理方便简单而得到较为广泛的应用，力可以表示为：

利用遗传算法工具箱（GA）对Bouc-Wen模型进行参数识别。

实验数据来源于对磁流变阻尼器（MR damper）进行性能测试，试验获得的数据包括力F，位移x，采用频率已知，速度和加速度可以由位移求导得出。

参数识别出现程序如下：（文件名：Copy_0_of_BoucWen）

function j=myfung(x)

y0=[0];

yy=y0;

tspan=[]';

s=[]';

v=[]';

Ft=[]';

rr=max(size(s));%计算数据个数

i=1;

while

(i1e5))%%判断是否出现奇异点，具体忘了。。

[t

y]=ode45(@uubird,[tspan(i),tspan(i+1)],y0,[],v(i),x);%参考论坛的

y0=y(end,:);

yy=[yy;y0];

i=i+1;

kk=max(size(y));

if kk>150 %微分方程计算，停止是有条件的（具体没去研究），这边设置150次，不管有没有收敛，都停止，不然整个程序运行的实际太久，你也可以改成其他的，慢慢研究

break;

end

end

if (i==rr)&(~isnan(yy(1,1)))==1%判断是否出现奇异点（就是NAN），如果没有出现，就是正常的

F=x(:,4)*yy(:,1)+x(:,5)*(s-ones(size(s)) *x(:,6))+x(:,7)*v;%x(:,4)代表alpha 5代表k0,6代表s0 7代表c0 位移s就是公式中的x

j=sum((F-Ft).*(F-Ft));

i=i+1;

else i<(rr-1)%出现奇异点（ＮＡＮ）

j=1e10;%因为出现奇异值，所以随便给一个目标函数值（这个要足够大），目的是排除这组优化值（也就是这个种群不要）

i=rr;

end

function dy=uubird(t,y,v,x)

dy=-x(:,1)*abs(v)*y*abs(y)-x(:,2)*(v)*y^2+x(:,3)*v;% 1r 2beita 3A 取n=2

运行结果

拟合效果对比程序（文件名BoucWenjianyan247）：

clc

clear

datasorce=load('247-0.txt');

II=1;

a=2000;

b=3000; %为了拟合。。

datasorce=datasorce(a:b,:);%去一个周期进行对比

s=datasorce(:,1)*0.7857-0.027860;%位移换算（你不用管）

Ft=datasorce(:,4)*1000/3.44; %力换算

sa=(max(s)-min(s))/2-max(s);

s=s+sa;

sFt=(max(Ft)-min(Ft))/2-max(Ft);

Ft=Ft+sFt;

for i=1:(b-a+1)

A(i)=i*0.001;

end

T=A';

p=polyfit(T,s(:),10);%拟合，用于求解

y10=polyval(p,T);

syms x

y= p(1)*x^10 + p(2)*x^9 +p(3)*x^8 + p(4)*x^7 + p(5)*x^6 + p(6)*x^5 + p(7)*x^4 + p(8)*x^3 +p(9)*x^2 + p(10)*x +p(11);

f=diff(y);

ac=diff(y,2);

for i=1:(b-a+1)

x=i*0.001;

v(i)=eval(f);

acc(i)=eval(ac);

vm(i)=(v(i)*v(i)-s(i)*acc(i))^0.5;

I(i)=0;

end

v=v';

T=T';

y0=[0];

yy=y0;

rr=max(size(s));

tspan=0:0.001:(rr-1)*0.001;

i=1;

x=[0.764 -0.7651 2.291 3.108

4.994 2.163925 3.002];%%参数识别出出来的结果

while (i

[t y]=ode45(@uubird,[tspan(i),tspan(i+1)],y0,[],v(i),x);

y0=y(end,:);

yy=[yy;y0];

i=i+1;

end

F=x(:,4)*yy(:,1)+x(:,5)*(s-ones(size(s))*x(:,6))+x(:,7)*v;%4alpha 5k0 6x0 7c0

s=s(100:2:800,:);

v=v(100:2:800,:);

F=F(100:2:800,:);

Ft=Ft(100:2:800,:);

figure(1),plot(v,Ft,'k',v,F,'r--'),hold on

figure(2),plot(s,Ft,'k',s,F,'r--'),hold on

实线是试验测得的结果。

[attach]188480[/attach][attach]188480[/attach][attach]188480[/attach]

如果觉得效果不够好，可以增加迭代次数。