www.pudn.com > ppso_train.zip > ppso_train.m, change:2012-02-21,size:5212b


%此算法为捕食搜索粒子群算法在列车运行调整中的应用 
clear all 
clc 
skb=xlsread('D:\开题与论文\论文1\列车运行调整时刻表');%读入列车时刻表 
[hang,lie]=size(skb); 
change_time=skb(:,1:2:lie)*60+skb(:,2:2:lie);%将列车时刻换算为分钟 
detal=zeros(hang,lie/2); 
detal(1,1)=detal(1,1)+10;%晚点时间长度 
ST0=change_time(1:2:hang,:);AR0=change_time(2:2:hang,:);%读取原定时刻表的到、发时间 
re_max=5;%初始化最大限制(限制总数),每30分钟为一限制。 
x=[]; 
x=[x;change_time+detal];%初始化,晚点时间+原计划时间[384,814+150] 
re=[3,10,15,20,25];%限制间隔 
x_ori=x; 
xc=x;%中心点  
restrict=1;%初始化最小限制区域 
Gbestmat=[]; 
w=0.8;c1=2.05;c2=2.05;%粒子群算法参数 
n=0;%迭代终止条件 
%alpha=[0.8,0.8,0.8,0.6,0.6,0.2,0.4,0.6,0.8,0.4,0.6,0.6];%列车等级权值 
alpha=[0.8,0.8,0.6,0.6,0.2,0.4]; 
q=0.8;%列车晚点时间与晚点数权值 
tau_b=2;tau_s=1;%起、停附加时间 
beta=4;%车站追踪间隔时间 
L=[3,2,2,3,2,3];%各车站股道数 
A=1440;%一天的时间,也用作罚函数 
T=[17,5,8,6,22;17,5,8,6,22;21,6,10,8,28;21,6,10,8,28;35,10,17,13,47;26,8,12,10,35]';%各列车在各区间最小运行时间矩阵 
B=[0,0,3,0;0,0,3,0;2,0,3,0;2,0,3,0;2,1,3,1;2,0,3,0]';%各列车在各站的最小作业时间矩阵 
 
aa=0;bb=0;cc=0;dd=0;ff=0;gg=0; 
 
while (restrict<=re_max) 
    for coun=1:5 
        %以下在最小限制区域内初始化粒子 
        ps=29;%随后生成粒子数 
        for i=2:ps+1 %随机生成粒子 
            st=xc(1:2:hang,:); 
            st(1,2:lie/2)=st(1,2:lie/2)+round(rand(1,lie/2-1).*10); 
            for j=1:4 
                ar(j,:)=st(j,:)+T(j,:)+round(rand(1,lie/2).*re(restrict));%按行生成列车到、发时刻 
                st(j+1,:)=ar(j,:)+B(j,:)+tau_b+tau_s+round(rand(1,lie/2).*re(restrict)); 
            end 
            ar(hang/2,:)=st(hang/2,:)+T(hang/2,:)+round(rand(1,lie/2).*re(restrict)); 
            x_a(1:2:hang,:)=st; 
            x_a(2:2:hang,:)=ar; 
            x=[x;x_a]; %生成随机点 
        end 
        ST=x(1:2:hang*(ps+1),:);%提取出发矩阵 
        AR=x(2:2:hang*(ps+1),:);%提取到达矩阵 
         
         
        %开始粒子群算法 
        pl=x;%初始赋值各个粒子历史最优 
        v=2*0.3*re(restrict)*rand(hang*(ps+1),lie/2)-0.3*re(restrict);%随机初始化粒子速度,范围为当前限制*0.3 
        for particle=1:ps+1 
            ar=AR(((hang/2)*particle-4):(hang/2)*particle,:); 
            st=ST(((hang/2)*particle-4):(hang/2)*particle,:); 
            [aa,bb,cc,dd,ff,gg,hh,fitness(particle)]=shiyingzhi(aa,bb,cc,dd,ff,gg,hang,lie,ar,st,AR0,ST0,tau_s,tau_b,beta,B,T,L,alpha,q,A); 
        end 
        fit_c=fitness(1);%中心点适应值 
        L_best=fitness';%初始赋值各个粒子历史最优适应值 
        [G_best,G_num]=min(L_best);%求取全局最优适应值及其序号 
        pg=x((hang*G_num-9):hang*G_num,:);%初始赋值全局最优 
        Gbestmat=[Gbestmat,G_best];%防止每次循环时Gbestmat矩阵都被更新 
 
 
        for it=1:500 
            for i=1:ps+1 
                v((hang*i-9):hang*i,:)=w*v((hang*i-9):hang*i,:)+c1.*rand().*(pl((hang*i-9):hang*i,:)-x((hang*i-9):hang*i,:))+c2.*rand()*(pg-x((hang*i-9):hang*i,:)); 
                  
                v_com=v((hang*i-9):hang*i,:);%中间变量 
                v_xiaoyu=find(v_com<-0.3*re(restrict));%速度越界设定为上下限*0.3 
                v_com(v_xiaoyu)=-0.3*re(restrict); 
                v_dayu=find(v_com>0.3*re(restrict)); 
                v_com(v_dayu)=0.3*re(restrict); 
                v((hang*i-9):hang*i,:)=v_com; 
                 
                x((hang*i-9):hang*i,:)=x((hang*i-9):hang*i,:)+round(v((hang*i-9):hang*i,:)); 
                 
                x_com=x((hang*i-9):hang*i,:); 
                x_bijiao=find(x_com>x_ori+re(restrict)|x_com<x_ori);%位置越界判定 
                x_com(x_bijiao)=x_ori(x_bijiao)+round(rand(size(x_bijiao,1),1).*re(restrict)); 
                x((hang*i-9):hang*i,:)=x_com; 
                 
                sk_com=x((hang*i-9):hang*i,:);%使用此中间变量目的为提取ar,st矩阵 
                [ha,li]=size(sk_com);%求取行数 
                st=sk_com(1:2:ha,:); 
                ar=sk_com(2:2:ha,:); 
                [aa,bb,cc,dd,ff,gg,hh,fit(i)]=shiyingzhi(aa,bb,cc,dd,ff,gg,hang,lie,ar,st,AR0,ST0,tau_s,tau_b,beta,B,T,L,alpha,q,A);%再次计算出每个粒子的适应值 
                if fit(i)<L_best(i) 
                    L_best(i)=fit(i);%更新各个粒子历史最优适应值 
                    pl((hang*i-9):hang*i,:)=x((hang*i-9):hang*i,:);%更新各个粒子历史最优值 
                end 
                if fit(i)<G_best 
                   G_best=fit(i);%每个粒子的适应值都与全局最优比较,更新全局最优适应值 
                   Gbestmat=[Gbestmat,G_best]; 
                   G_best=min(Gbestmat); 
                   pg=x((hang*i-9):hang*i,:);%更新全局最优 
                end 
            end 
        end 
        n=n+1;%设置迭代终止条件 
        if fit_c-G_best>0  
             xc=pg;%更新中心点 
             x=xc;%必须将xc赋值为x 
             restrict=0; 
             break 
        else 
            x=xc; 
            continue 
        end 
    end 
    if n>=15%终止条件为迭代了800次以上 
        break 
    end 
    restrict=restrict+1; 
end 
plot(Gbestmat); 
xiaoshi=fix(pg./60);%求取小时 
fenzhong=pg-fix(pg./60).*60;%求取分钟 
tiaozhenghou(:,1:2:12)=xiaoshi; 
tiaozhenghou(:,2:2:12)=fenzhong; 
success=xlswrite('D:\开题与论文\论文1\调整后',tiaozhenghou);