第八講基于群體智能的參數(shù)辨識(shí)

基于群體智能算法的參數(shù)辨識(shí)遺傳算法(GeneticAlgorithm)優(yōu)勝劣汰遺傳過(guò)程常用術(shù)語(yǔ)

個(gè)體(染色體)

群體

復(fù)制

交叉

變異

選擇

適應(yīng)度

適應(yīng)值優(yōu)化目標(biāo)目標(biāo)函數(shù)適應(yīng)度函數(shù)優(yōu)化目標(biāo)函數(shù)適應(yīng)度函數(shù)適應(yīng)值通過(guò)適應(yīng)度函數(shù)來(lái)決定染色體的優(yōu)劣程度遺傳編碼二進(jìn)制編碼（編碼，譯碼）對(duì)每個(gè)參數(shù)確定它的變化范圍，并用一個(gè)二進(jìn)制數(shù)來(lái)表示，若參數(shù)a的變化范圍為(amin,amax)，用m位二進(jìn)制數(shù)b來(lái)表示，則兩者之間滿(mǎn)足實(shí)數(shù)編碼選擇方法輪盤(pán)賭交叉操作(二進(jìn)制)單點(diǎn)交叉：交叉點(diǎn)k的范圍為[1，N-1]，N為個(gè)體變量數(shù)目，在該點(diǎn)分界相互交換變量多點(diǎn)交叉，m個(gè)交叉位置可無(wú)重復(fù)隨機(jī)地選擇，在交叉點(diǎn)之間的變量間續(xù)地相互交換，產(chǎn)生兩個(gè)新的后代均勻交叉將每個(gè)點(diǎn)都作為潛在的交叉點(diǎn)，隨機(jī)地產(chǎn)生與個(gè)體等長(zhǎng)的0、1掩碼，掩碼的片斷表明了哪個(gè)父體向子體提供變量值交叉操作（實(shí)值）子個(gè)體1=a*父?jìng)€(gè)體1+(1-a)*父?jìng)€(gè)體2子個(gè)體2=(1-a)*父?jìng)€(gè)體1+a*父?jìng)€(gè)體2a為1到0之間的隨機(jī)數(shù)變異（二進(jìn)制）對(duì)于二進(jìn)制編碼個(gè)體而言，變異意味著變量的對(duì)某個(gè)二進(jìn)制位執(zhí)行求反運(yùn)算變異前：101100111001101變異后：101101111001101變異（實(shí)數(shù)）新個(gè)體=舊個(gè)體+(rand-0.5)*(amax-amin)

clc;clearall;closeall;Ge=100;%代數(shù)m=80;%一次選擇個(gè)體數(shù)目Pc=0.9;%交叉概率Pm=0.05;%變異概率ts=0.5;%采樣周期%用于記憶最優(yōu)個(gè)體對(duì)應(yīng)適應(yīng)值%對(duì)應(yīng)參數(shù)值bestQQ=10000;k_s=0;T1_s=0;T2_s=0;%被擬合對(duì)象參數(shù)st=70;lp=st/ts;yout0=0;yout=0;k1=4;t11=5;t22=9;timef=0;x00=0;x01=0;r=1.0;%階躍輸入a00=exp(-ts/t11);b00=1-a00;a01=exp(-ts/t22);b01=1-a01;forii=1:lpx00=a00*x00+k1*b00*r;x01=a01*x01+1*b01*x00;yout0(ii)=x01;timef(ii)=ii*ts;endforkk=1:Gex1l=3.8;x1r=4.2;%kx2l=4;x2r=6;%T1x3l=7;x3r=9;%T2fori=1:m%產(chǎn)生初始種群

pop1(i)=rand(1)*(x1r-x1l)+x1l;pop2(i)=rand(1)*(x2r-x2l)+x2l;pop3(i)=rand(1)*(x3r-x3l)+x3l;end

%按照誤差平方的最小值計(jì)算目標(biāo)函數(shù)

fori=1:mk=pop1(i);T1=pop2(i);T2=pop3(i);x=[0,0];yout=0;error=0;Q=0;a1=exp(-ts/T1);b1=1-a1;a2=exp(-ts/T2);b2=1-a2;forj=1:lpx(1)=a1*x(1)+k*b1*r;x(2)=a2*x(2)+1*b2*x(1);yout(j)=x(2);error(j)=yout(j)-yout0(j);Q=Q+(error(j)^2);endQQ(i)=Q;end%計(jì)算適應(yīng)度，由于是直接求最小值問(wèn)題，所以直接將Q取為個(gè)體的適應(yīng)度

sumQQ=sum(QQ);cfitness=QQ/sumQQ;fori=1:mcfitness(i)=(1-cfitness(i))/(m-1);end%根據(jù)適應(yīng)度大小選擇下一代，（比例選擇）

fori=2:mcfitness(i)=cfitness(i-1)+cfitness(i);end%選擇運(yùn)算fori=1:mp=rand(1);index=1;while(p>cfitness(index))index=index+1;endnewpop1(i)=pop1(index);newpop2(i)=pop2(index);newpop3(i)=pop3(index);end%放回原來(lái)數(shù)組里fori=1:m

pop1(i)=newpop1(i);pop2(i)=newpop2(i);pop3(i)=newpop3(i);endfori=1:m%交叉運(yùn)算

p=rand(1);index1(i)=floor(rand(1)*(m-1))+1;index2(i)=floor(rand(1)*(m-1))+1;index3(i)=floor(rand(1)*(m-1))+1;if(p<Pc)A=rand(1);pop1(i)=pop1(i)*A+(1-A)*pop1(index1(i));B=rand(1);pop2(i)=pop2(i)*B+(1-B)*pop2(index2(i));C=rand(1);pop3(i)=pop3(i)*C+(1-C)*pop3(index3(i));endend

%變異運(yùn)算fori=1:mp=rand(1);ifp<Pmpop1(i)=pop1(i)+Leng1*(rand(1)-0.5)*(x1r-x1l);pop2(i)=pop2(i)+Leng2*(rand(1)-0.5)*(x2r-x2l);pop3(i)=pop3(i)+Leng3*(rand(1)-0.5)*(x3r-x3l);if(pop1(i)>x1r)pop1(i)=x1r;endif(pop1(i)<x1l)pop1(i)=x1l;endif(pop2(i)>x2r)pop2(i)=x2r;endif(pop2(i)<x2l)pop2(i)=x2l;endif(pop3(i)>x3r)pop3(i)=x3r;endif(pop3(i)<x3l)pop3(i)=x3l;endendend

%循環(huán)結(jié)束，求出最優(yōu)個(gè)體輸出a1=exp(-ts/T1);b1=1-a1;a2=exp(-ts/T2);b2=1-a2;x=[0,0];yout=0;forj=1:lpx(1)=a1*x(1)+k*b1*r;x(2)=a2*x(2)+1*b2*x(1);yout(j)=x(2);endfigure;plot(timef,yout,‘b--’,timef,yout0,’k’);figure;Plot(errorhistor);%一次迭代結(jié)束，%找出最優(yōu)個(gè)體存儲(chǔ)下來(lái)[f,I]=min(QQ);iff<bestQQbestQQ=f;k_s=pop1(I);T1_s=pop2(I);T2_s=pop3(I);enderrorhistory=[errorhistorybestQQ];end%TPSO(ParticleSwarmOptimization)ParticleSwarmOptimizationTom:WhereshouldIflyto?Tom:Thisisthebestplaceinmysearchingfoodhistory.John:Thisisbestplaceforourswarmnow.PSO算法采用速度一位置搜索模型。每個(gè)粒子代表解空間的一個(gè)候選解，解的優(yōu)劣程度由適應(yīng)度函數(shù)決定，而適應(yīng)度函數(shù)根據(jù)優(yōu)化目標(biāo)定義。PSO算法隨機(jī)初始化一群粒子，然后通過(guò)迭代找到最優(yōu)解。在每一次迭代中，粒子通過(guò)跟蹤兩個(gè)“極值”來(lái)更新自己。一個(gè)是粒子本身所找到的最優(yōu)解，即個(gè)體極值pbest；，另一個(gè)是整個(gè)種群目前找到的最優(yōu)解，稱(chēng)之為全局極值gbest。PSO的初始化（1）設(shè)定群體規(guī)模m，確定搜索空間[xmin,xmax]和單步前進(jìn)的最大值vmax（一般應(yīng)小于最大搜索空間的20%）；（2）對(duì)任意粒子，在[xmin,xmax]內(nèi)服從均勻分布產(chǎn)生個(gè)體xi；（3）對(duì)任意粒子，在[-vmax,vmax]內(nèi)服從均勻分布產(chǎn)生個(gè)體vin；（4）對(duì)任意粒子，設(shè)Pi=Xi。%對(duì)一個(gè)系統(tǒng)進(jìn)行仿真，并且計(jì)算絕對(duì)誤差的矩積分，評(píng)價(jià)函數(shù)functionQ=jianmoobj(k,T1,T2);globaltimefyoutyout0rints=0.5;st=100;lp=st/ts;rin=1.0;x0=[0,0];x=[0,0];yout0=0;yout=0;error=0;Q=0;k0=1.5;t01=6;t02=16;a01=exp(-ts/t01);b01=1-a01;a02=exp(-ts/t02);b02=1-a02;a11=exp(-ts/T1);b11=1-a11;a12=exp(-ts/T2);b12=1-a12;fori=1:lptimef(i)=i*ts;r(i)=rin;x(1)=x(1)*a11+k*b11*r(i);x(2)=x(2)*a12+1*b12*x(1);yout(i)=x(2);x0(1)=x0(1)*a01+k0*b01*r(i);x0(2)=x0(2)*a02+1*b02*x0(1);yout0(i)=x0(2);error(i)=yout(i)-yout0(i);Q=Q+(error(i)^2);end%繪制最終的輸出曲線(xiàn)functionQ=jianmoout(k,T1,T2);globaltimefyoutyout0rints=0.5;rin=1.0;x0=[0,0];x=[0,0];yout0=0;yout=0;error=0;st=100;lp=st/ts;k0=1.5;t01=6;t02=16;a01=exp(-ts/t01);b01=1-a01;a02=exp(-ts/t02);b02=1-a02;a11=exp(-ts/T1);b11=1-a11;a12=exp(-ts/T2);b12=1-a12;fori=1:lptimef(i)=i*ts;r(i)=rin;x(1)=x(1)*a11+k*b11*r(i);x(2)=x(2)*a12+1*b12*x(1);yout(i)=x(2);x0(1)=x0(1)*a01+k0*b01*r(i);x0(2)=x0(2)*a02+1*b02*x0(1);yout0(i)=x0(2);endfigure(2);plot(timef,rin,'r',timef,yout,'b',timef,yout0,'g');xlabel('Time(s)');ylabel('rin,yout,yout0');%PSO算法closeall;clearall;globaltimefyoutyout0rin%隨機(jī)初始化Gp=50;%循環(huán)代數(shù)errorhistory=[];%初始化各個(gè)粒子x=rand(30,3);[n,m]=size(x);%初始化各個(gè)粒子的第一步前進(jìn)的速度v=rand(size(x));dd=0.2;xmin=[1514];xmax=[2718];vmax=dd.*(xmax-xmin);fori=1:nforj=1:mx(i,j)=rand(1)*(xmax(1,j)-xmin(1,j))+xmin(1,j);v(i,j)=rand(1)*vmax(1,j);endendfori=1:ny(i)=jianmoobj(x(i,1),x(i,2),x(i,3));%計(jì)算原始粒子的目標(biāo)值

py=zeros(1,n);endp=x;%原始粒子的歷史最優(yōu)位置既是本身ymin=inf;%原始粒子中的最優(yōu)值k=1;%找出未迭代前所有粒子中的最優(yōu)值fori=1:nif((ymin)>y(i))ymin=y(i);k=i;endendg=x(k,:);%未迭代前原始粒子中的最優(yōu)值c1=1.5;%常數(shù)c2=1.5;y0=y;%保留各個(gè)粒子目標(biāo)值的前一時(shí)刻值gy=inf;%全局目標(biāo)最優(yōu)值fork=1:gp%慣性因子

%w=1-0.8*k/lp;w=0.73;%個(gè)體歷史最優(yōu)值

fori=1:nif(y0(i)>=y(i))forj=1:mp(i,j)=x(i,j);%p為個(gè)體最優(yōu)值

endpy(i)=y(i);%py為個(gè)體最優(yōu)目標(biāo)值

endendy0=y;%全局最優(yōu)值

fori=1:nif((gy)>(py(i)))g=p(i,:);gy=py(i);endend%計(jì)算個(gè)體前進(jìn)速度，并限制最大值

fori=1:nforj=1:mv(i,j)=w*v(i,j)+c1*rand(1)*(p(i,j)-x(i,j))+c2*rand(1)*(g(j)-x(i,j));end%速度最大值限制

forj=1:m