<<MATLAB在數學建模中的應用>>

案例一.無約束目標函式最大值遺傳演算法求解策略

求解問題           

%主程式:用遺傳演算法求解y=200*exp(-0.05*x).*sin(x)在[-2,2]上的最大值
clc;
clear all;
close all;
global BitLength
global boundsbegin
global boundsend
bounds=[-2 2];   %一維自變數的取值範圍
precision=0.0001;%運算精度
boundsbegin=bounds(:,1);
boundsend=bounds(:,2);
%計算如果滿足求解精度至少需要多長的染色體
BitLength=ceil(log2((boundsend-boundsbegin)'./precision));
popsize=50;       %初始種群大小
Generationnmax=30;%最大代數
pcrossover=0.90;  %交配概率
pmutation=0.09;   %變異概率
%產生初始種群
population=round(rand(popsize,BitLength));
%計算適應度,返回適應度Fitvalue和累計概率cumsump
[Fitvalue,cumsump]=fitnessfun(population);
Generation=1;
while Generation<Generationnmax
   for j=1:2:popsize
       %選擇操作
       seln=selection(population,cumsump);
       %交叉操作
       scro=crossover(population,seln,pcrossover);
       scnew(j,:)=scro(1,:);
       scnew(j+1,:)=scro(2,:);
       %變異操作
       smnew(j,:)=mutation(scnew(j,:),pmutation);
       smnew(j+1,:)=mutation(scnew(j+1,:),pmutation);
   end
   population=smnew;%產生了新的種群
   %計算新種群的適應度
   [Fitvalue,cumsump]=fitnessfun(population);
   %記錄當前代最好的適應度和平均適應度
   [fmax,nmax]=max(Fitvalue);
   fmean=mean(Fitvalue);
   ymax(Generation)=fmax;
   ymean(Generation)=fmean;
   %記錄當前代最好的適應度和平均適應度
   x=transform2to10(population(nmax,:));
   %自變數取值範圍是[-2,2],需要把經過遺傳運算的最佳染色體整合到[-2,2]區間
   xx=boundsbegin+x*(boundsend-boundsbegin)/(power((boundsend),BitLength)-1);
   xmax(Generation)=xx;
   Generation=Generation+1;
end
Generation=Generation-1;
Bestpopulation=xx;
Besttargetfunvalue=targetfun(xx);
%繪製經過遺傳演算法後的適應度曲線.一般的,如果進化過程中種群的平均適應度與最大適應度
%在曲線上有相互趨同的形態,表示演算法收斂進行很順利,沒有出現震盪;
%在這種前提下,最大適應度個體連續若干代沒有發生進化表明種群已經成熟
figure(1);
hand1=plot(1:Generation,ymax);
set(hand1,'linestyle','-','linewidth',1.8,'marker','*','markersize',6);
hold on
hand2=plot(1:Generation,ymean);
set(hand2,'color','r','linestyle','-','linewidth',1.8,...
    'marker','h','markersize',6);
xlabel('進化代數');ylabel('最大/平均適應度');xlim([1 Generation]);
legend('最大適應度','平均適應度');
box off;hold off;
 

function scro=crossover(population,seln,pc)
BitLength=size(population,2);
pcc=IfCroIfMut(pc);%根據交叉概率決定是否進行交叉操作,1則是,0則否
if pcc == 1
    chb=round(rand*(BitLength-2))+1;%在[1,BitLength-1]範圍內隨機產生一個交叉位
    scro(1,:)=[population(seln(1),1:chb) population(seln(2),chb+1:BitLength)];
    scro(2,:)=[population(seln(2),1:chb) population(seln(1),chb+1:BitLength)];
else
    scro(1,:)=population(seln(1),:);
    scro(2,:)=population(seln(2),:);
end
end
 

function [Fitvalue,cumsump]=fitnessfun(population)
global BitLength
global boundsbegin
global boundsend
popsize=size(population,1);
for i=1:popsize
   x=transform2to10(population(i,:));%將二進位制轉化為10進位制
   %轉化為[-2,2]區間的實數
   xx=boundsbegin+x*(boundsend-boundsbegin)/(power((boundsend),BitLength)-1);
   Fitvalue(i)=targetfun(xx);       %計算函式值,即適應度
end
%給適應度加上一個大小合理的數以便保證種群適應值為正數
Fitvalue=Fitvalue'+230;
%計算選擇概率
fsum=sum(Fitvalue);
Pperpopulation=Fitvalue/fsum;
%計算累積概率
cumsump(1)=Pperpopulation(1);
for i=2:popsize
    cumsump(i)=cumsump(i-1)+Pperpopulation(i);
end
cumsump=cumsump';
end
 

function pcc=IfCroIfMut(mutORcro)
test(1:100)=0;
l=round(100*mutORcro);
test(1:l)=1;
n=round(rand*99)+1;
pcc=test(n);
end

function snnew=mutation(snew,pmutation)
    BitLength=size(snew,2);
    snnew=snew;
    pmm=IfCroIfMut(pmutation);%根據變異概率決定是否進行編譯操作,1則是,0則不是
    if pmm == 1
       chb=round(rand*(BitLength-1))+1;%在[1,BitLength]範圍內隨機產生一個變異位
       snnew(chb)=abs(snew(chb)-1);
    end
end

function seln=selection(population,cumsump)
%從種群中選擇兩個個體
for i=1:2
    r=rand;         %隨機產生一個隨機數
    prand=cumsump-r;
    j=1;
    while prand(j)<0
        j=j+1;
    end
    seln(i)=j;      %選中個體的序號
end
end

function x=transform2to10(Population)
    BitLength=size(Population,2);
    x=Population(BitLength);
    for i=1:BitLength-1
        x=x+Population(BitLength-i)*power(2,i);
    end
end

function y=targetfun(x)
y = 200*exp(-0.05*x).*sin(x);
end