最優(yōu)化原理及應(yīng)用

1、“最優(yōu)化原理及應(yīng)用”2008200388 姚遠1用C語言，或者Matlab, 或者Fortran等編寫一個完整的Simulated Annealing算法和Genetic 算法的優(yōu)化程序。解: 本題采用Matlab語言編寫一個完整的SA算法優(yōu)化程序。在該程序中選用的代價函數(shù)為：，初始的C0=1000，每一個階段的Lk選為20，接受概率設(shè)為0.6，迭代的終止條件為e0.00001) k=k+1; C=C0/k; for(i=1:Lk) w=2*rand(1,1)-1; x1=x0+w*delta_x; %產(chǎn)生一個x1 F1=x14-x13-15*x12+1; delta_f=F1-F; e=ab

2、s(delta_f); if (F1epsilon) x0=x1; F=F1; X(i)=i+1; end end endendx0F下面為采用遺傳算法的優(yōu)化程序。該程序中的代價函數(shù)與上面的SA算法所用的一致。設(shè)定變量的二進制碼鏈長度為10，基因庫中的二進制碼鏈個數(shù)為10。自變量的取值區(qū)間為，設(shè)定遺傳算法的迭代次數(shù)為500。在每次迭代保留基因的選擇中分別采用了均值法和輪盤賭的方法，在保留、交換和異化的過程中每次都將最好的基因保留在基因庫的最后一行，即精英操作。經(jīng)過GA算法得出的結(jié)果為：x1=5.6696 f(x)=-3.8851。從圖中我們可以看出SA和GA算法均找到了該代價函數(shù)的最小值點。在

3、運行的過程中GA的速度要明顯快于SA。程序如下：%遺傳算法% F=sin(x-1.5)2)+(x-6)2-3; 最優(yōu)化問題的代價函數(shù)% Q=1000-( sin(x-1.5)2)+(x-6)2-3); 遺傳算法中定義的fitnessclcclear allnum=10;length=10;total=2(length)-1;min=-5; %取值區(qū)間的長度max=5;G=100;genome=round(rand(num,length); %定義自變量的隨機二進制編碼，長度為lengthfor (k=1:G) X=zeros(1,num); %二進制轉(zhuǎn)化為十進制 for(i=1:num) fo

4、r(j=length:-1:1) X(i)=X(i)+genome(i,j)*2(length-j); end X(i)=X(i)/total*(max-min)+min; end for(i=1:num) %計算各變量的fitnessQ(i)=1000-(sin(X(i)-1.5)2)+(X(i)-6)2-3); end order,location=sort(Q); bestgene=location(num); %最好的基因位置 BEST(k,:)=genome(bestgene,:); %最好的基因 % %保留% q=sum(Q);% vector=(Q/q)*num;% vector

5、=floor(vector); %將品質(zhì)因數(shù)大于均值的gene選出來,保留到下次% j=1;% for (i=1:num)% if (vector(i)=1)% geneupdate(j,:)=genome(i,:);% j=j+1;% end% end% geneupdate(num,:)=BEST(k,:); %保留 %使用輪盤賭的篩選方法 q=sum(Q); selectvariable=ceil(q)*rand(1,1); j=1; while (1) a=0; i=1; a=a+Q(i); while (1) if(aselectvariable) i=i+1; a=a+Q(i);

6、else geneupdate(j,:)=genome(i,:); j=j+1; break; end end if(j=num) break; end end geneupdate(num,:)=BEST(k,:); %交換 kk=1; probc=0.5; for(i=1:2:num-1) variable1=rand(1,1); variable2=floor(length-1)*rand(1,1)+1); %找出截斷點 if (variable1probc) for (j=variable2:length) a(kk)=geneupdate(i,j); b(kk)=geneupdate

7、(i+1,j); kk=kk+1; end kk=1; for (j=variable2:length) geneupdate(i,j)=b(kk); geneupdate(i+1,j)=a(kk); kk=kk+1; end end end geneupdate(num,:)=BEST(k,:); genome=geneupdate; %異化 probt=0.1; for (i=1:num) for (n=1:length) variable3=rand(1,1); if (variable3probt) if (genome(i,n)=1) genome(i,n)=0; else geno

8、me(i,n)=1; end end end end geneupdate(num,:)=BEST(k,:); genome=geneupdate;end x=0; for(j=length:-1:1) x=x+genome(num,j)*2(length-j); end x=x/total*(max-min)+min; xF=sin(x-1.5)2)+(x-6)2-3; F2. 利用模擬退火和遺傳算法求函數(shù)的最小值點，x的精度控制在0.001范圍內(nèi)。解：在本題中所用的代價函數(shù)為：，初始的C0=1000，每一階段的Lk取為20，接受概率設(shè)為0.6，迭代的終止條件為e0.000001) k=k+

9、1; C=C0/k; for(i=1:Lk) w=2*rand(1,1)-1; x1=x0+w*delta_x; F1=2*x12-x1-1; delta_f=F1-F; e=abs(delta_f); if (F1epsilon) x0=x1; F=F1; end end endend x0 F X=-1:0.1:1; y=2*X.2-X-1; plot(X,y); grid on;set(gcf,color,w);下面利用GA算法來計算f(x)的最小值，為了達到題目中給定的精度要求，這里設(shè)定二進制碼鏈的長度為11，基因庫中二進制碼鏈的個數(shù)為10，自變量的取值區(qū)間定為：。遺傳算法的迭代次數(shù)為

10、100。在每次迭代保留基因的選擇中分別采用了均值法和輪盤賭的方法，在保留、交換和異化的過程中每次都將最好的基因保留在基因庫的最后一行，即精英操作。經(jīng)過GA算法得出的結(jié)果為：x1=0.2496 f(x)=-1.1250。滿足題目中所給的精度要求。程序如下：%遺傳算法% F=2*x2-x-1; 最優(yōu)化問題的代價函數(shù)% Q=1000-(2*x2-x-1); 遺傳算法中定義的fitnessclcclear allnum=10;length=11;total=2(length)-1;min=-1; %取值區(qū)間的長度max=1;G=100;genome=round(rand(num,length); %定

11、義自變量的隨機二進制編碼，長度為lengthfor (k=1:G) X=zeros(1,num); %二進制轉(zhuǎn)化為十進制 for(i=1:num) for(j=length:-1:1) X(i)=X(i)+genome(i,j)*2(length-j); end X(i)=X(i)/total*(max-min)+min; end for(i=1:num) %計算各變量的fitness Q(i)=1000-(2*X(i)2-X(i)-1); end order,location=sort(Q); bestgene=location(num); %最好的基因位置 BEST(k,:)=genome

12、(bestgene,:); %最好的基因 % %保留% q=sum(Q);% vector=(Q/q)*num;% vector=floor(vector); %將品質(zhì)因數(shù)大于均值的gene選出來,保留到下次% j=1;% for (i=1:num)% if (vector(i)=1)% geneupdate(j,:)=genome(i,:);% j=j+1;% end% end% geneupdate(num,:)=BEST(k,:); %保留 %使用輪盤賭的篩選方法 q=sum(Q); selectvariable=ceil(q)*rand(1,1); j=1; while (1) a=0

13、; i=1; a=a+Q(i); while (1) if(aselectvariable) i=i+1; a=a+Q(i); else geneupdate(j,:)=genome(i,:); j=j+1; break; end end if(j=num) break; end end geneupdate(num,:)=BEST(k,:); %交換 kk=1; probc=0.5; for(i=1:2:num-1) variable1=rand(1,1); variable2=floor(length-1)*rand(1,1)+1); if (variable1probc) for (j=

14、variable2:length) a(kk)=geneupdate(i,j); b(kk)=geneupdate(i+1,j); kk=kk+1; end kk=1; for (j=variable2:length) geneupdate(i,j)=b(kk); geneupdate(i+1,j)=a(kk); kk=kk+1; end end end geneupdate(num,:)=BEST(k,:); genome=geneupdate; %異化 probt=0.1; for (i=1:num) for (n=1:length) variable3=rand(1,1);% n=flo

15、or(7*rand(1,1)+1); if (variable3probt) if (genome(i,n)=1) genome(i,n)=0; else genome(i,n)=1; end end end end geneupdate(num,:)=BEST(k,:); genome=geneupdate;end x=0; for(j=length:-1:1) x=x+genome(num,j)*2(length-j); end x=x/total*(max-min)+min; x F=2*x2-x-1; F3. 利用模擬退火和遺傳算法求函數(shù)的最小值點，x的精度控制在0.001范圍內(nèi)。解：

16、首先采用SA算法，取C0=1000，每階段的Lk=100，接受概率取為0.5，設(shè)定迭代次數(shù)K=100，為了加快收斂，令C0=C0/log(k)。x，y的擾動都定為0.1。從給定的代價函數(shù)我們?nèi)菀卓闯銎渥钚≈迭c是（0，0）。經(jīng)過SA運算可以得到優(yōu)化后的結(jié)果為：x1=2.2114e-004，y1=3.1916e-004。x，y的收斂趨勢如下圖所示：結(jié)果達到了題目中的精度要求。程序如下：% SA退火算法計算F=2*x2+y2的最小值clear allclcC0=1000;x0=0.5*(2*rand(1,1)-1);y0=0.5*(2*rand(1,1)-1);k=1;Lk=1000;F=2*x02

17、+y02;Fdelta_x=0.1;delta_y=0.1;epsilon=0.8;K=100;I=1;for (j=1:K) k=k+1; C0=C0/log(k); for(i=1:Lk) x1=x0+(2*rand(1,1)-1)*delta_x; y1=y0+(2*rand(1,1)-1)*delta_y; F1=2*x12+y12; delta_f=F1-F; if (F1epsilon) x0=x1; y0=y1; F=F1; temp1(I)=x1; temp2(I)=y1; I=I+1; end end endendplot (temp1,r); grid on;set(gcf

18、,color,w);hold onplot (temp2,g); grid on;set(gcf,color,w); x0 y0 F采用GA算法，為了達到題中所給的精度要求，所以取二進制碼鏈的長度為22。其中前11位表示x，后11位表示y?；驇炖锒M制碼鏈的個數(shù)為25個。遺傳迭代200次。采用輪盤賭和均值篩選兩種方式，進行精英操作。得出的結(jié)果為： x1=4.8852e-004，y1=4.8852e-004。將每次迭代的精英代碼對應(yīng)的x，y值畫出來可以得到下圖：得到的結(jié)果達到了題目中所要求的精度要求。程序如下：%應(yīng)用遺傳算法計算F=2*x2+y2的最小值clcclear allnum=25;l

19、ength_x=11;length_y=11;length=length_x+length_y;total_x=2(length_x)-1;total_y=2(length_y)-1;min=-1; %取值區(qū)間的長度max=1;G=200;k=1;genome=round(rand(num,length); %自變量定義的隨機二進制編碼，長度為lengthfor (i=1:G) %二進制轉(zhuǎn)化為十進制 X=zeros(1,num); Y=zeros(1,num); for(i=1:num) for(j=length_x:-1:1) X(i)=X(i)+genome(i,j)*2(length_x

20、-j); end X(i)=X(i)/total_x*(max-min)+min; end for(i=1:num) for(j=length:-1:length-length_x+1) Y(i)=Y(i)+genome(i,j)*2(length-j); end Y(i)=Y(i)/total_y*(max-min)+min; end for(i=1:num) %計算各變量的fitness Q(i)=1000-(2*X(i)2+Y(i)2); end order,location=sort(Q); bestgene=location(num); %最好的基因位置 BEST(k,:)=geno

21、me(bestgene,:); %最好的基因 %保留 q=sum(Q); vector=Q/q*num; vector=floor(vector); %將品質(zhì)因數(shù)大于均值的gene選出來,保留到下次 j=1; for (i=1:num) if (vector(i)=1) geneupdate(j,:)=genome(i,:); j=j+1; end end geneupdate(num,:)=BEST(k,:); %交換 kk=1; probc=0.5; for(i=1:2:num-1) variable1=rand(1,1); variable2=floor(length-1)*rand(1

22、,1)+1); if (variable1probc) for (j=variable2:length) a(kk)=geneupdate(i,j); b(kk)=geneupdate(i+1,j); kk=kk+1; end kk=1; for (j=variable2:length) geneupdate(i,j)=b(kk); geneupdate(i+1,j)=a(kk); kk=kk+1; end end end geneupdate(num,:)=BEST(k,:); genome=geneupdate; %異化 probt=0.2; for (i=1:num) for (n=1:

23、length) variable3=rand(1,1);% n=floor(7*rand(1,1)+1); if (variable3probt) if (genome(i,n)=1) genome(i,n)=0; else genome(i,n)=1; end end end end geneupdate(num,:)=BEST(k,:); genome=geneupdate; k=k+1;end x=0; for(j=length_x:-1:1) x=x+genome(num,j)*2(length_x-j); end x=x/total_x*(max-min)+min; x y=0; f

24、or(j=length:-1:length-length_x+1) y=y+genome(num,j)*2(length-j); end y=y/total_x*(max-min)+min; y % X=zeros(1,G); Y=zeros(1,G); for(i=1:G) for(j=length_x:-1:1) X(i)=X(i)+BEST(i,j)*2(length_x-j); end X(i)=X(i)/total_x*(max-min)+min; end for(i=1:G) for(j=length:-1:length-length_x+1) Y(i)=Y(i)+BEST(i,j

25、)*2(length-j); end Y(i)=Y(i)/total_y*(max-min)+min; end plot (X,g); grid on;set(gcf,color,w); hold on plot (Y,r); grid on;set(gcf,color,w);4推銷員問題是很有名的最佳化問題。假定一個推銷員必須訪問某區(qū)域里的N個村莊，他從第一號村莊出發(fā)走訪其他N-1個村莊。如果村莊i和村莊j 之間的距離為dij ,那么他如何選擇路線才能使他走的距離最近？給出其最佳化解決方法。解：在該題中假設(shè)有10個村莊，代價函數(shù)定為走遍所有村莊的距離和，固定一個起始點，運用SA算法對其進行優(yōu)

26、化。得到以下的結(jié)果：沒有優(yōu)化之前的路徑：一號村莊的坐標為（46，26）。總長度為：554.6568優(yōu)化之后的路徑：總長度為：237.0437程序如下：% 利用SA算法實現(xiàn)推銷員問題clear allclcC0=500000000;Lk=2000;N=10;% position=round(100*rand(2,N); %隨機給出10個村莊的位置坐標position=46,4,61,92,11,32,13,96,42,71;26,27,77,44,99,93,20,32,92,79;%給出10個村莊的固定位置F0=0;F1=0;k=1;epsilon=0.5; posorder0=1:N; po

27、sorder1=1:N; for (i=1:N-1) F0=F0+sqrt(position(1,posorder0(i+1)-position(1,posorder0(i)2 +(position(2,posorder0(i+1)-position(2,posorder0(i)2); end F0 %代價函數(shù) %循環(huán)開始的地方while (1) k=k+1; C0=C0/log(k); if (C0eps) break; end for (i=1:Lk) posorder1=posorder0; %隨機交換除第一個村莊之外的其余村莊的位置 p=ceil(N-1)*rand(1,1)+1);

28、while (p=1) %避免p=1，即避免交換第一個村莊 p=ceil(N-1)*rand(1,1)+1); end q=ceil(N-1)*rand(1,1)+1); while (q=1) q=ceil(N-1)*rand(1,1)+1); end temp=posorder1(p); posorder1(p)=posorder1(q); posorder1(q)=temp; F1=0; for (i=1:N-1) F1=F1+sqrt(position(1,posorder1(i+1)-position(1,posorder1(i)2 +(position(2,posorder1(i+

29、1)-position(2,posorder1(i)2); end if (F1epsilon F0=F1; posorder0=posorder1; end end end end F0 for (i=1:N) position1(:,i)=position(:,posorder0(i); end figure (1) set(gcf,color,w); plot(position1(1,:),position1(2,:); hold on; plot(position1(1,:),position1(2,:),*); title=(優(yōu)化后的行走路線); figure (2) set(gcf

30、,color,w); plot(position(1,:),position(2,:); hold on plot(position(1,:),position(2,:),*); title=(沒有經(jīng)過優(yōu)化的路線);5如下遞推公式：是有名的Widrow LMS算法遞推公式。其中Wn 是權(quán)向量，Xn 是輸入向量，是常數(shù)，n 是誤差函數(shù)。該公式的推導是一個典型的二次方程式加隨機編程處理最佳化問題，請給出該公式的完整推導。（參閱B.Widrow的Adaptive signal processing 教材和龔耀寰的“自適應(yīng)原理及應(yīng)用”或其他有關(guān)自適應(yīng)信號處理書籍）解：由各變量的定義，最小均方誤差MSE

31、=E，根據(jù)遞推公式：，所以，將代入上式，并對兩邊求數(shù)學期望，可以得到：，又因為，所以可以求得。6基陣方向圖設(shè)計是陣設(shè)計的基本問題，對于等間隔線列陣，已經(jīng)證明切比雪夫加權(quán)是最佳權(quán)系數(shù)，其最佳化意義在于對于相同次瓣級，主瓣寬度最窄；而對于固定主瓣寬度次瓣級最低（參見Urick“工程水聲原理”）。假設(shè)有一6元等間隔線列陣，試利用模擬退火或遺傳算法求其最佳加權(quán)系數(shù)。并畫出有關(guān)波束圖* 由切比雪夫加權(quán)法求得的最佳加權(quán)系數(shù)為0.3, 0.69, 1, 1, 0.69, 0.3,可利用該系數(shù)求出方向圖函數(shù)G0(),為空間角，范圍-90,90,構(gòu)成代價函數(shù)f=|G()G0()|和品質(zhì)函數(shù)F=C-f。解：應(yīng)用SA退火算法求一個6元等間隔線列陣的最佳加權(quán)系數(shù)，假設(shè)6元的等間隔線列陣相鄰陣元的間距為半波長，信號入射角度為theta，聲程差計算以線列陣的中心為參考點。所以可以得到該時延對應(yīng)的相位phi為,由于相位為余弦表達式，且在推導過程中以6元直線陣的中點為參考點，所以可以得到6元陣的輸出為：余弦為對稱函數(shù)，且權(quán)值也是左右對稱，所以可以簡化上式得到：上式即為所要的代價函數(shù)。采用SA算法優(yōu)化權(quán)值，設(shè)定C0=1000，迭代過程使C0=C0*log(k)。每個過程Lk=100。運算得到的結(jié)果如下：0.2927 0.6977