|
楼主 |
发表于 2007-5-30 10:55
|
显示全部楼层
好的,我现在把重要的代码都帖上来。
%这是主程序
clear all;
clc;
global X K
popsize=2; %种群规模
c1=2; %个体最优导向系数
c2=2; %全局最优导向系数
%gbest_x; %全局最优解
%best_fitness; %最优解
%best_in_history; %最优解变化轨迹
w=inf;
iw1=0.9;
iw2=0.4;
X_min=5*10^(-5); %x的下限
X_max=5*10^(-5)+2.0*10^(-4); %x的上限
gen=2; %迭代次数
%exetime; t; %当前迭代次数
max_velocity=4.1*10^(-5); %最大速度
%eps=10^(-10);%设置精度(在已知最小值时候用)
%最大速度
best_in_history(gen)=inf; %初始化全局历史最优解
best_in_history(:)=inf; %初始化全局历史最优解
best_fitness=inf;
pop(popsize,5)=0;
for K=1:popsize
pop(K,1)=1.5*10^(-4)*rand(1);%初始化种群中的粒子位置,值为-1—1,步长为其速度
% X(K)=pop(K,1)
pop(K,3)=pop(K,1); %初始状态下个体最优值等于初始位置
pop(K,2)=1.0*10^(-6)+4*10^(-5)*rand(1); %初始化种群微粒速度,值为-0.01—0.01,间隔为0.0001
pop(K,4)=inf;
pop(K,5)=inf;
end
gbest_X=pop(1,1); %全局最优初始值为种群第一个粒子的位置
for t=1:gen
%计算权值
if t<=3
iwt(t)=((iw2-iw1)/(gen-1))*(t-1)+iw1
else
iwt(t)=iw2
end
adapting1 %计算适应值
updatepop1 %更新粒子位置
pause(0.01);
end
gbest_X
disp('最后得到的优化极值为:')
best_fitness
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
这是adaping1文件
%适值计算
global X K
%计算适应值并赋值
for K=1:popsize;
X(K)=pop(K,1)
w=PHASE1(pop(K,1),K)
pop(K,5)=w;
% pop(K,5)=100*pop(K,1)^2+(1-pop(K,1))^2;
if pop(K,4)>pop(K,5) %若当前适应值优于个体最优值,则进行个体最优信息的更新
pop(K,4)=pop(K,5); %适值更新
pop(K,3)=pop(K,1); %位置坐标更新
end
end
%计算完适应值后寻找当前全局最优位置并记录其坐标
if best_fitness>min(pop(:,4))
best_fitness=min(pop(:,4)); %全局最优值
gbest_X=pop(find(pop(:,4)==min(pop(:,4))),1) %全局最优粒子的位置
end
best_in_history(t)=best_fitness %记录当前全局最优
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
这是updatepop1文件
%粒子群速度与位置更新
global X K
%更新粒子速度
for K=1:popsize
pop(K,2)=iwt(t)*pop(K,2)+c1*rand(1)*(pop(K,3)-pop(K,1))+c2*rand(1)*(gbest_X-pop(K,1)); %更新速度
if abs(pop(K,2))>max_velocity
if pop(K,2)>0
pop(K,2)=max_velocity;
else
pop(K,2)=-max_velocity;
end
end
end
%更新粒子位置
for K=1:popsize
pop(K,1)=pop(K,1)+pop(K,2);
q=pop(K,1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z=PHASE1(k1,K)
r=0*10^(-6):1*10^(-6):24*10^(-6);
global X K
k1=X(K);
for m=1:length(r);
An(m)=quadl(@FUN1,0,8,1e-7,[],r(m));
double(An(m));
AM(m)=abs(An(m));
PH2es(m)=angle(An(m))*180/3.14;
end
s=0;
k1=1*10^(-4);k2=5*10^(-7);Rth=1*10^(-7);
PH1=PHASE(k1,k2,Rth);
PH1ex=PH1;
A=[ 0.0077 0.0195 0.0281 0.0159 -0.0216
-0.0063 0.0290 0.0278 -0.0060 0.0262
-0.0295 -0.0098 0.0209 -0.0327 -0.0241
-0.0198 -0.0201 0.0069 -0.0152 -0.0201
-0.0323 0.0165 -0.0037 0.0288 -0.0023];
for m=1:length(PH1ex);
PH2ex(m)=PH1ex(m).*(1+A([m]));
diff(m)=PH2ex(m)-PH2es(m);
s=s+(abs(diff(m)))^2;
end
z=s;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=FUN1(deta,r, k1)
global X K
f=1*10^5;
w=2*pi*f;
L=0.1*10^(-6);
k2=5*10^(-7);Rth=1*10^(-7);
ra=3*10^(-6);rb=0*10^(-6);
reff=sqrt(ra^2+rb^2);
K0=2.61*10^(-2);k0=2.0*10^(-5);
K1=2.5*10^6*k1;
K2=3*10^6*k2;
alpha1=5.31*10^7;Rp=0.817;P=0.1;
sigma0=sqrt(deta.*deta./(reff^2)+j*w/k0);
sigma1=sqrt(deta.*deta./(reff^2)+j*w/k1);
sigma2=sqrt(deta.*deta./(reff^2)+j*w/k2);
g1=(K0*sigma0)./(K1*sigma1);
g2=(K2*sigma2)./(K1*sigma1);
b=Rth*K1*sigma1;
s=alpha1./sigma1;
E=alpha1*(1-Rp)*P*exp(-deta.*deta/8)./[2*pi*K1*(sigma1.*sigma1-alpha1^2)];
H=(1+g1).*[1+g2.*(1+b)].*exp(sigma1*L)-(1-g1).*[1-g2.*(1-b)].*exp(-sigma1*L);
A1=[(1-g1).*[s-g2.*(1-b.*s)].*exp(-alpha1*L)-(g1+s).*[1+g2.*(1+b)].*exp(sigma1*L)].*E./H;
B1=[(1+g1).*[s-g2.*(1-b.*s)].*exp(-alpha1*L)-(g1+s).*[1-g2.*(1-b)].*exp(-sigma1*L)].*E./H;
F=[A1+B1+E].*deta.*besselj(0,deta.*r/reff)/(reff^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 由上到下依次为内部需要调用的函数文件。注这里面的PHASE又是一个函数,运行时调用一直以来没显示错误的,就不贴了。
[ 本帖最后由 xiaoshuizhu 于 2007-5-30 10:59 编辑 ] |
|