matlab程序例子.docx
- 文档编号:7343645
- 上传时间:2023-05-11
- 格式:DOCX
- 页数:16
- 大小:18.20KB
matlab程序例子.docx
《matlab程序例子.docx》由会员分享,可在线阅读,更多相关《matlab程序例子.docx(16页珍藏版)》请在冰点文库上搜索。
matlab程序例子
%cou=q;
%thickness=0.016071;%cm
thickness=0.02;%cm
M=500;
step=5/(M-1);
wvl=[1547.5:
step:
1552.5];%%%4nm的波长范围分成1000个点
Neff=1.5;
wvlD2=1550*1e-7;
deviate2=2*pi*Neff*(1./(wvl*1e-7)-1/wvlD2);%%%失谐量
r=zeros(1,M);%r为反射率不是反射系数
form1=1:
M
ifm1>200&m1<300
r(1,m1)=1/1*step*(300-m1);
else
r(1,m1)=0;
end
end
%开始对折射率调制进行优化,采用的是量子粒子群算法
%------给定初始化条件----------------------------------------------
%N1---种群个数
%D---个体的维数
%------初始化种群的个体(可以在这里限定位置的范围)------------
N1=20;
D=60;
fork1=1:
N1
fork=1:
D
modu(k1,k)=-6+12*rand+i*(-6+12*rand);
end
end
%------先计算各个粒子的适应度,并初始化Pbest和gbest----------------------
%适应度f=[1/(M-1)*(所有点的(Rt-Ri)^2和)]^2
%-------传输矩阵法求各个粒子的适应值
%run=1;
%whilerun<=1000
absreflect1=zeros(1,M);
summ=0;
fork1=1:
N1
forj=1:
M
Matrix_g=[1,0;0,1];
fork=1:
D
delta1=deviate2(1,j);
kac=modu(k1,k);
alpha=sqrt(abs(kac)^2-delta1^2);
alphaL=alpha.*thickness;
T11=cosh(alphaL)+i*(delta1/alpha)*sinh(alphaL);
T12=(kac/alpha)*sinh(alphaL);
T21=(conj(kac)/alpha)*sinh(alphaL);
T22=cosh(alphaL)-i*(delta1/alpha)*sinh(alphaL);
Matrix_g=Matrix_g*[T11,T12;T21,T22];
end
TT11=Matrix_g
(1);TT12=Matrix_g(3);TT21=Matrix_g
(2);TT22=Matrix_g(4);
reflect1=-TT21/TT11;
absreflect1(1,j)=abs(reflect1);
Reflectivity1(j)=(abs(reflect1)).^2;
%%---------加入波长平衡因子
ifj>200&j<300
W(1,j)=[((300-j)*step)/1]^1;
else
W(1,j)=1;
end
Fit=[Reflectivity1(j)]-[r(1,j)];%反射
summ=abs(Fit^
(2))*W(1,j)+summ;
end
Fitness(1,k1)=summ/M;
summ=0;
end
%--------根据各个粒子的适应值判断粒子的初始化Pbest和gbest------------------------
Tiedaicishu=6001;
fitnessgbest=zeros(1,1);
gbest=zeros(D,1);
ggbest=zeros(D,Tiedaicishu);
Gbest=zeros(D,1);
mbest=zeros(1,D);
Fitnessbest=zeros(1,N1);
Fitnessleast=zeros(1,1);
fori3=1:
N1
pbest(i3,:
)=modu(i3,:
);%pbest为各个粒子当前最优
end
Fitnessbest(1,:
)=Fitness(1,:
);
SUM=0;
fork4=1:
D
fori3=1:
N1
SUM=pbest(i3,k4)+SUM;
end
mbest(1,k4)=SUM/N1;%mbest为D维空间中各个粒子自身最优位置的中心点
SUM=0;
end
%-----------------当前全局最优位置的适应值
Fitness1=Fitnessbest;
Fitnessleast=Fitness1(1,1);
foryyyy=2:
N1
ifFitnessleast>Fitness1(1,yyyy)
Fitnessleast=Fitness1(1,yyyy);%当前群体最优适应值
end
end
%--------------------------%-----------------当前全局最优粒子位置
foryyyy1=1:
N1
ifFitnessbest(1,yyyy1)==Fitnessleast
Gbest(:
1)=pbest(yyyy1,:
);
MOO=yyyy1;%Gbest为全局最优的初始值
end
end
gbest(:
1)=Gbest(:
1);%群体全局最优适应值
fitnessgbest=Fitnessleast;%第tiedaicishu次迭代后的全局最优值
ggbest(:
1)=gbest(:
1);
ff
(1)=fitnessgbest;
%不同的变异因子前保存初始值
gbest1(:
1)=Gbest(:
1);%群体全局最优适应值
fitnessgbest1=Fitnessleast;%第tiedaicishu次迭代后的全局最优值
ggbest1(:
1)=gbest1(:
1);
ff0
(1)=fitnessgbest;
Fitnessbest1=Fitnessbest;
mbest1=mbest;
%---------------------------
%改变以后重新初始
gbest=gbest1;
fitnessgbest=fitnessgbest1;
ggbest1=ggbest1;
Fitnessbest=Fitnessbest1;
ff
(1)=ff0;
mbest=mbest1;
%----------------------
%---------------------------
modu1=modu;
%-------------------------%------进入主要循环,按照公式依次迭代,直到满足精度要求------------
tiedaicishu=1;
TTiedaicishu=2500;
whiletiedaicishu<=TTiedaicishu
fori4=1:
N1
ifabs(Fitnessbest(1,i4))==abs(fitnessgbest)
xx(1,i4)=-230;
elseifabs(Fitnessbest(1,i4)) min_Fitness=abs(Fitnessbest(1,i4)); error_Fitness=(Fitnessbest(1,i4)-fitnessgbest)/min_Fitness; xx(1,i4)=log(error_Fitness); else min_Fitness=abs(fitnessgbest); error_Fitness=(Fitnessbest(1,i4)-fitnessgbest)/min_Fitness; xx(1,i4)=log(error_Fitness); end %%--------根据第tiedaicishu次迭代时xx的值判断a的值 ifxx(1,i4)>0.5 a(tiedaicishu,i4)=0.01; elseif(xx(1,i4)<=0.5)&(xx(1,i4)>0) a(tiedaicishu,i4)=0.05; elseif(xx(1,i4)<=0)&(xx(1,i4)>-0.5) a(tiedaicishu,i4)=0.1; elseif(xx(1,i4)<=-0.5)&(xx(1,i4)>-1) a(tiedaicishu,i4)=0.15; elseif(xx(1,i4)<=-1)&(xx(1,i4)>-1.5) a(tiedaicishu,i4)=0.2; elseif(xx(1,i4)<=-1.5)&(xx(1,i4)>-2) a(tiedaicishu,i4)=0.25; elseif(xx(1,i4)<=-2)&(xx(1,i4)>-2.5) a(tiedaicishu,i4)=0.3; elseif(xx(1,i4)<=-2.5)&(xx(1,i4)>-3) a(tiedaicishu,i4)=0.35; elseif(xx(1,i4)<=-3)&(xx(1,i4)>-3.5) a(tiedaicishu,i4)=0.4; elseif(xx(1,i4)<=-3.5)&(xx(1,i4)>-4) a(tiedaicishu,i4)=0.45; elseif(xx(1,i4)<=-4)&(xx(1,i4)>-4.5) a(tiedaicishu,i4)=0.5; elseif(xx(1,i4)<=-4.5)&(xx(1,i4)>-5) a(tiedaicishu,i4)=0.55; elseif(xx(1,i4)<=-5)&(xx(1,i4)>-5.5) a(tiedaicishu,i4)=0.6; elseif(xx(1,i4)<=-5.5)&(xx(1,i4)>-6) a(tiedaicishu,i4)=6.5; else a(tiedaicishu,i4)=0.7; end %a(tiedaicishu,i4)=0.1; fork1=1: D Pbest(i4,k1)=rand*pbest(i4,k1)+[1-rand]*gbest(k1,1);%Pbest迭代中粒子位置的变化的中间值 U=rand; ifU>0.5 modu1(i4,k1)=Pbest(i4,k1)-a(tiedaicishu,i4)*abs(mbest(1,k1)-modu1(i4,k1))*log(1/U); else modu1(i4,k1)=Pbest(i4,k1)+a(tiedaicishu,i4)*abs(mbest(1,k1)-modu1(i4,k1))*log(1/U); end end end %-----位置更替后判断新的Pbest和gbest------------ %s首先求适应值 Summ=0; fork1=1: N1 forj=1: M Matrix_g=[1,0;0,1]; fork=1: D delta1=deviate2(1,j); kac=modu(k1,k); alpha=sqrt(abs(kac)^2-delta1^2); alphaL=alpha.*thickness; T11=cosh(alphaL)+i*(delta1/alpha)*sinh(alphaL); T12=(kac/alpha)*sinh(alphaL); T21=(conj(kac)/alpha)*sinh(alphaL); T22=cosh(alphaL)-i*(delta1/alpha)*sinh(alphaL); Matrix_g=Matrix_g*[T11,T12;T21,T22]; end TT11=Matrix_g (1);TT12=Matrix_g(3);TT21=Matrix_g (2);TT22=Matrix_g(4); reflect1=-TT21/TT11; absreflect1(1,j)=abs(reflect1); Reflectivity1(j)=(abs(reflect1)).^2; %---------加入波长平衡因子 %%---------加入波长平衡因子 ifj>200&j<300 W(1,j)=[((300-j)*step)/1]^1; else W(1,j)=1; end Fit1=(Reflectivity1(j))-(r(1,j)); Summ=abs(Fit1^2)*W(1,j)+Summ;% end Fitness(tiedaicishu+1,k1)=Summ/M;%当前迭代个体适应值 Summ=0; end %-------------------------------- %更新的Pbest和gbest fori4=1: N1 ifFitness(tiedaicishu+1,i4) 1)历史取得的单个粒子最优目标函数值 pbest(i4,: )=modu1(i4,: );%更新的Pbest Fitnessbest(1,i4)=Fitness(tiedaicishu+1,i4);%个体最优适应值 end end %-------------------------------- %引入变异因子 fori4=1: N1 fork1=1: D %pbest1(i4,k1)=pbest(i4,k1)+0.1*(normrnd(0,1)+normrnd(0,1)*i);%rand*(1+0.5*i); %pbest1(i4,k1)=pbest(i4,k1)+normrnd(real(pbest(i4,k1)),0.1)+normrnd(imag(pbest(i4,k1)),0.1)*i; pbest1(i4,k1)=normrnd(real(pbest(i4,k1)),0.4)+normrnd(imag(pbest(i4,k1)),0.4)*i; end end Summm=0; fork1=1: N1 forj=1: M Matrix_g=[1,0;0,1]; fork=1: D delta1=deviate2(1,j); kac=pbest1(k1,k); alpha=sqrt(abs(kac)^2-delta1^2); alphaL=alpha.*thickness; T11=cosh(alphaL)+i*(delta1/alpha)*sinh(alphaL); T12=(kac/alpha)*sinh(alphaL); T21=(conj(kac)/alpha)*sinh(alphaL); T22=cosh(alphaL)-i*(delta1/alpha)*sinh(alphaL); Matrix_g=Matrix_g*[T11,T12;T21,T22]; end TT11=Matrix_g (1);TT12=Matrix_g(3);TT21=Matrix_g (2);TT22=Matrix_g(4); reflect1=-TT21/TT11; absreflect1(1,j)=abs(reflect1); Reflectivity1(j)=(abs(reflect1)).^2; %%---------加入波长平衡因子 %%---------加入波长平衡因子 ifj>200&j<300 W(1,j)=[((300-j)*step)/1]^1; else W(1,j)=1; end Fit11=Reflectivity1(j)-r(1,j); Summm=abs(Fit11^2)*W(1,j)+Summm; end Fitness1(tiedaicishu+1,k1)=Summm/M;%当前迭代个体适应值 Summm=0; end %----------------比较变异前后适应值更新位置 fori4=1: N1 ifFitness1(tiedaicishu+1,i4) 1)历史取得的单个粒子最优目标函数值 pbest(i4,: )=pbest1(i4,: );%更新的Pbest Fitnessbest(1,i4)=Fitness1(tiedaicishu+1,i4);%个体最优适应值 end end %--------------------------当前迭代次数的全局最优的位置的适应值 Fitness2=Fitnessbest; Fitnessleast=Fitness2(1,1); fori5=2: N1 ifFitnessleast>Fitness2(1,i5) Fitnessleast=Fitness2(1,i5); end end %--------%--------------------------当前迭代次数的全局最优的位置- fori51=1: N1 ifFitnessbest(1,i51)==Fitnessleast Gbest(: 1)=pbest(i51,: );%Gbest为第tiedaicishu次迭代时全局最优的 end end %--------%--------------------------历史全局最优的位置 ifFitnessleast gbest(: 1)=Gbest(: 1);%gbest为群体全局最优的 fitnessgbest=Fitnessleast; end ggbest(: tiedaicishu+1)=gbest(: 1); %--------------------------- SUMM=0; fork5=1: D fori6=1: N1 SUMM=pbest(i6,k5)+SUMM; end mbest(1,k5)=SUMM/N1;%mbest为D维空间中各个粒子自身最优位置的中心点 SUMM=0; end ff(tiedaicishu)=fitnessgbest; tiedaicishu=tiedaicishu+1; end %---------用传输矩阵法验证所叠加后的sneff的最优值gbest(tiedaicishu,: )给q,从而求反射谱. %figure,plot(wvl,10*log10(Reflectivity1),wvl,20*log10(r)) fff=0: 2/124: 2; absreflect1=zeros(1,M); forj=1: M Matrix_g=[1,0;0,1]; fork=1: D delta1=deviate2(1,j); kac=ggbest(k,1); %kac=coupling12mm422(k,1); alpha=sqrt(abs(kac)^2-delta1^2); alphaL=alpha.*thickness; T11=cosh(alphaL)+i*(delta1/alpha)*sinh(alphaL); T12=(kac/alpha)*sinh(alphaL); T21=(conj(kac)/alpha)*sinh(alphaL); T22=cosh(alphaL)-i*(delta1/alpha)*sinh(alphaL); Matrix_g=Matrix_g*[T11,T12;T21,T22]; end TT11=Matrix_g (1);TT12=Matrix_g(3);TT21=Matrix_g (2);TT22=Matrix_g(4); reflect1=-TT21/TT11; absreflect(1,j)=abs(reflect1); Reflectivity2(j)=(abs(reflect1)).^2; Q(j)=phase(reflect1); end figure,plot(wvl,Reflectivity2,wvl,r) figure,plot(wvl,Reflectivity) ggbest2(: 1)=ggbest(: 1)*(1+0.1); ggbest3(: 1)=ggbest(: 1)*(1-0.1); ggbest4(: 1)=ggbest(: 1)*(1+0.1*normrnd(0,1)); subplot(1,2,1);plot(fff,ggbest2(: 1),fff,ggbest3(: 1),fff,ggbest4(: 1),fff,ggbest(: 1)); subplot(1,2,2);plot(wvl,Reflectivity2,wvl,Reflectivity3,wvl,Reflectivity4,wvl,Reflectivity) fff=1: 1: 60; figure,plot(fff,ggbest(: 5410)); fffff=1: 1: 2500; figure,semilogy(fffff,ff(1,1: 2500)) sneffleft=zeros(1,D);%存q angle1eft=zeros(1,D);%存q forn1=1: D sneffleft(1,n1)=1550*abs(qqleft12mm430(n1,1)*10^(-7)/pi);%coupling12mm422 angle1eft(1,n1)=angle(qqleft12mm430(n1,1))/pi; end %%%%%%%%求时延 form1=1: M ifm1<=30 Reflectivity1(1,m1)=Reflec
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 程序 例子
![提示](https://static.bingdoc.com/images/bang_tan.gif)