二维光子晶体禁带的遗传优化算法

1、二维光子晶体禁带的遗传优化算法(2007-03-10 10:21:26)function Zp,Yp,Xp,LC1,LC2=PBGGA(M,N,Pm)%-% PBGGA.m% Photonic Band Gap Genetic Algorithm% 二维光子晶体禁带的遗传优化算法%-% 函数功能% 在平面波展开法的基础上，使用遗传算法寻找具有最大光子禁带的二维光子晶体结构 %-% 主要参考文献% 龚春娟,胡雄伟.快速遗传算法优化计算二维光子晶体J.% 半导体学报.2006,27(6):1098-1102%-% 输入参数列表% M 遗传进化迭代次数% N 种群规模(取偶数)% Pm 变异概率%-

2、% Zp 最大禁带% Yp 最大禁带的区间上下限% Xp 光子晶体结构决策变量（10×10的01矩阵）% LC1 收敛曲线1，各代最优个体适应值的记录% LC2 收敛曲线2，各代群体平均适应值的记录%-%第一步：变量初始化epsa=11.56;%介质介电常数epsb=1;%空气介电常数Xp=zeros(10,10);%最优结构Yp=zeros(1,2);%禁带上下限Zp=0;%禁带宽度LC1=zeros(1,M);%收敛曲线1LC2=LC1;%收敛曲线2%第二步：随机产生初始钟群farm=cell(1,N);%采用细胞结构存储种群for i=1:NA=unidrnd(2,5,10)-

3、ones(5,10);%由于基本单元结构是中心对称的，所以只需要对矩阵的上半部分进行初始化B=flipud(fliplr(A);%剩余部分先做上下倒置，再做左右旋转C=A;B;%把两步分合并farmi=reshape(C,100,1);%个体被存储在细胞结构farm中endcounter=0;%设置迭代计数器while counter<M%停止条件为达到最大迭代次数%第三步：交叉newfarm=cell(1,N);%交叉产生的新种群存在这里面Ser=randperm(N);%用这个函数保证随机配对for i=1:2:(N-1)AA=reshape(farmSer(i),10,10);%

4、待交叉的两个父代个体 BB=reshape(farmSer(i+1),10,10);A=AA(1:5,:);%为保证对称性，需先提取其中一半的基因B=BB(1:5,:);Manner=unidrnd(2);%随机选择交叉方式if Manner=1%横交叉cp=unidrnd(4);%随机选择交叉点a=A(1:cp,:);B(cp+1):5,:);%双亲双子单点交叉b=B(1:cp,:);A(cp+1):5,:);else%列交叉cp=unidrnd(9);%随机选择交叉点a=A(:,1:cp),B(:,(cp+1):10);%双亲双子单点交叉b=B(:,1:cp),A(:,(cp+1):10)

5、;enda=a;flipud(fliplr(a);%反演对称b=b;flipud(fliplr(b);newfarmi=reshape(a,100,1);%交叉后的子代存入newfarm newfarmi+1=reshape(b,100,1);end%新旧种群合并FARM=farm,newfarm;%第四步：选择复制FITNESS=zeros(1,2*N);fitness=zeros(1,N);for i=1:(2*N)X=FARMi;Band=pixelTM(X,10,epsa,epsb,0);%调用平面波展开法的程序，中间结果不画图B1=Band(:,1);B2=Band(:,2);for

6、 j=1:length(B1)if B1(j)>1|B2(j)>1B1(j)=0;B2(j)=0;endendBB=abs(B1-B2);pos=find(BB=max(BB);Z=BB(pos(1);%最宽的禁带作为适应值FITNESS(i)=Z;endSer=randperm(2*N);%选择复制采取两两随机配对竞争的方式，具有保留最优个体的能力for i=1:Nf1=FITNESS(Ser(2*i-1);f2=FITNESS(Ser(2*i);if f1>=f2farmi=FARMSer(2*i-1);fitness(i)=FITNESS(Ser(2*i-1);else

7、farmi=FARMSer(2*i);fitness(i)=FITNESS(Ser(2*i);endend%记录最佳个体和收敛曲线maxfitness=max(fitness)meanfitness=mean(fitness)LC1(counter+1)=maxfitness;%收敛曲线1，各代最优个体适应值的记录 LC2(counter+1)=meanfitness;%收敛曲线2，各代群体平均适应值的记录 pos=find(fitness=maxfitness);xp=farmpos(1);Zp=maxfitness;Xp=reshape(xp,10,10);%第五步：变异for i=1:N

8、if Pm>rand;%变异概率为PmX=farmi;XX=reshape(X,10,10);%为保证对称性，需先将个体基因变形 XA=XX(1:5,:);posX=unidrnd(5,1,1);%随机选取变异位posY=unidrnd(10,1,1);XA(posX(1),posY(1)=mod(XA(posX(1),posY(1)+1,2);%变异，0变成1，1编程0%XA(posX(2),posY(2)=mod(XA(posX(2),posY(2)+1,2); XB=flipud(fliplr(XA);XX=XA;XB;farmi=reshape(XX,100,1);%还原个体形态endendfarmpos(1)=reshape(Xp,100,1);counter=counter+1end%输出结果并绘图figure(1);xp=reshape(Xp,100,1);Band=pixelTM(xp,10,epsa,epsb,1);%画最优个体的禁带图B1=Band(:,1);B2=Ba