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
%--------%--------------------------历史全局最优的位置
ifFitnessleastgbest（:
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
升级会员 