【原创】二维光子晶体禁带的遗传优化算法

aiwa

GreenSim团队的博客开通了，更多原创程序请访问我们的博客
http://blog.sina.com.cn/greensim

1. function [Zp,Yp,Xp,LC1,LC2]=PBGGA(M,N,Pm)
   
2. %--------------------------------------------------------------------------
   
3. % PBGGA.m
   
4. % Photonic Band Gap Genetic Algorithm
   
5. % 二维光子晶体禁带的遗传优化算法
   
6. %--------------------------------------------------------------------------
   
7. % 函数功能
   
8. % 在平面波展开法的基础上，使用遗传算法寻找具有最大光子禁带的二维光子晶体结构
   
9. %--------------------------------------------------------------------------
   
10. % 主要参考文献
    
11. % 龚春娟,胡雄伟.快速遗传算法优化计算二维光子晶体[J].
    
12. % 半导体学报.2006,27(6):1098-1102
    
13. %--------------------------------------------------------------------------
    
14. % 输入参数列表
    
15. % M 遗传进化迭代次数
    
16. % N 种群规模(取偶数)
    
17. % Pm 变异概率
    
18. %--------------------------------------------------------------------------
    
19. % Zp 最大禁带
    
20. % Yp 最大禁带的区间上下限
    
21. % Xp 光子晶体结构决策变量（10×10的01矩阵）
    
22. % LC1 收敛曲线1，各代最优个体适应值的记录
    
23. % LC2 收敛曲线2，各代群体平均适应值的记录
    
24. %--------------------------------------------------------------------------
    
25. 
    
26. %第一步：变量初始化
    
27. epsa=11.56;%介质介电常数
    
28. epsb=1;%空气介电常数
    
29. Xp=zeros(10,10);%最优结构
    
30. Yp=zeros(1,2);%禁带上下限
    
31. Zp=0;%禁带宽度
    
32. LC1=zeros(1,M);%收敛曲线1
    
33. LC2=LC1;%收敛曲线2
    
34. 
    
35. %第二步：随机产生初始钟群
    
36. farm=cell(1,N);%采用细胞结构存储种群
    
37. for i=1:N
    
38. A=unidrnd(2,5,10)-ones(5,10);%由于基本单元结构是中心对称的，所以只需要对矩阵的上半部分进行初始化
    
39. B=flipud(fliplr(A));%剩余部分先做上下倒置，再做左右旋转
    
40. C=[A;B];%把两步分合并
    
41. farm{i}=reshape(C,100,1);%个体被存储在细胞结构farm中
    
42. end
    
43. 
    
44. counter=0;%设置迭代计数器
    
45. while counter<M%停止条件为达到最大迭代次数
    
46. 
    
47. %第三步：交叉
    
48. newfarm=cell(1,N);%交叉产生的新种群存在这里面
    
49. Ser=randperm(N);%用这个函数保证随机配对
    
50. for i=1:2:(N-1)
    
51. AA=reshape(farm{Ser(i)},10,10);% 待交叉的两个父代个体
    
52. BB=reshape(farm{Ser(i+1)},10,10);
    
53. A=AA(1:5,:);%为保证对称性，需先提取其中一半的基因
    
54. B=BB(1:5,:);
    
55. Manner=unidrnd(2);%随机选择交叉方式
    
56. if Manner==1%横交叉
    
57. cp=unidrnd(4);%随机选择交叉点
    
58. a=[A(1:cp,:);B((cp+1):5,:)];%双亲双子单点交叉
    
59. b=[B(1:cp,:);A((cp+1):5,:)];
    
60. else%列交叉
    
61. cp=unidrnd(9);%随机选择交叉点
    
62. a=[A(:,1:cp),B(:,(cp+1):10)];%双亲双子单点交叉
    
63. b=[B(:,1:cp),A(:,(cp+1):10)];
    
64. end
    
65. a=[a;flipud(fliplr(a))];%反演对称
    
66. b=[b;flipud(fliplr(b))];
    
67. newfarm{i}=reshape(a,100,1);%交叉后的子代存入newfarm
    
68. newfarm{i+1}=reshape(b,100,1);
    
69. end
    
70. %新旧种群合并
    
71. FARM=[farm,newfarm];
    
72. 
    
73. %第四步：选择复制
    
74. FITNESS=zeros(1,2*N);
    
75. fitness=zeros(1,N);
    
76. for i=1:(2*N)
    
77. X=FARM{i};
    
78. Band=pixelTM(X,10,epsa,epsb,0);%调用平面波展开法的程序，中间结果不画图
    
79. B1=Band(:,1);
    
80. B2=Band(:,2);
    
81. for j=1:length(B1)
    
82. if B1(j)>1||B2(j)>1
    
83. B1(j)=0;
    
84. B2(j)=0;
    
85. end
    
86. end
    
87. BB=abs(B1-B2);
    
88. pos=find(BB==max(BB));
    
89. Z=BB(pos(1));%最宽的禁带作为适应值
    
90. FITNESS(i)=Z;
    
91. end
    
92. Ser=randperm(2*N);%选择复制采取两两随机配对竞争的方式，具有保留最优个体的能力
    
93. for i=1:N
    
94. f1=FITNESS(Ser(2*i-1));
    
95. f2=FITNESS(Ser(2*i));
    
96. if f1>=f2
    
97. farm{i}=FARM{Ser(2*i-1)};
    
98. fitness(i)=FITNESS(Ser(2*i-1));
    
99. else
    
100. farm{i}=FARM{Ser(2*i)};
     
101. fitness(i)=FITNESS(Ser(2*i));
     
102. end
     
103. end
     
104. %记录最佳个体和收敛曲线
     
105. maxfitness=max(fitness)
     
106. meanfitness=mean(fitness)
     
107. LC1(counter+1)=maxfitness;%收敛曲线1，各代最优个体适应值的记录
     
108. LC2(counter+1)=meanfitness;%收敛曲线2，各代群体平均适应值的记录
     
109. pos=find(fitness==maxfitness);
     
110. xp=farm{pos(1)};
     
111. Zp=maxfitness;
     
112. Xp=reshape(xp,10,10);
     
113. 
     
114. %第五步：变异
     
115. for i=1:N
     
116. if Pm>rand;%变异概率为Pm
     
117. X=farm{i};
     
118. XX=reshape(X,10,10);%为保证对称性，需先将个体基因变形
     
119. XA=XX(1:5,:);
     
120. posX=unidrnd(5,1,1);%随机选取变异位
     
121. posY=unidrnd(10,1,1);
     
122. XA(posX(1),posY(1))=mod(XA(posX(1),posY(1))+1,2);%变异，0变成1，1编程0
     
123. %XA(posX(2),posY(2))=mod(XA(posX(2),posY(2))+1,2);
     
124. XB=flipud(fliplr(XA));
     
125. XX=[XA;XB];
     
126. farm{i}=reshape(XX,100,1);%还原个体形态
     
127. end
     
128. end
     
129. farm{pos(1)}=reshape(Xp,100,1);
     
130. 
     
131. counter=counter+1
     
132. end
     
133. 
     
134. %输出结果并绘图
     
135. figure(1);
     
136. xp=reshape(Xp,100,1);
     
137. Band=pixelTM(xp,10,epsa,epsb,1);%画最优个体的禁带图
     
138. B1=Band(:,1);
     
139. B2=Band(:,2);
     
140. BB=abs(B1-B2);
     
141. pos=find(BB==max(BB));
     
142. Zp=BB(pos(1));
     
143. Yp=[B1(pos(1)),B2(pos(1))];
     
144. figure(2);
     
145. plot(LC1);%画收敛曲线1
     
146. figure(3);
     
147. plot(LC2);%画收敛曲线2

复制代码

--------------------------------------------------------------------------------

【GreenSim团队简介】

——更多原创程序，请访问GreenSim团队的主页→http://blog.sina.com.cn/greensim

    GreenSim团队长期从事建模仿真、算法设计、代写程序等业务，为您提供实验、课题、论文、项目等方面的仿真服务。
    GreenSim团队擅长的学科领域主要有
【智能算法】遗传算法、模拟退火、神经网络、蚁群算法、免疫优化、支持向量机、拟物拟人
【运筹优化】数学建模、数值优化算法、大规模/非线性问题、组合优化、NP完全问题
【信息与通信】无线通信、信号处理、计算机网络
【其它学科】任何问题，只要能归结为数学问题或者算法仿真，我们都有实力和信心解决
    GreenSim团队的联系方式：
QQ: 330264704，535115369  Email：greensim@163.com（来 信 必 复）

[ 本帖最后由 aiwa 于 2007-3-18 13:34 编辑 ]

junziyang

好东西。建议版主加分！

zhangbit

运行不了啊  说M没有定义

[ 本帖最后由 zhangbit 于 2008-1-10 01:44 编辑 ]

筆調零亂

pixelTM这个函数能分享一下吗，谢谢您！

handb

楼主太牛了，赞