function [FullEntry,TNUM,PM,SPN,STFRM,STFRP,STNRP,STPM,STRM,STSDM,STSM,STSVM,STTRP]=MAIN_Procedure_new
%This is the main procedure,trying to resolve the 100 pedes in 360s in 30*10 matrix.
tic;
clear;clc;
SM=zeros(30,10);PM=[];SVM=[];SDM=[];PP=[];%the initialization of [] matrix.
STSVM=zeros(100,360);STSDM=zeros(100,360);STPM=zeros(100,360);
STSM=zeros(30,10,360);STFRP=cell(100,360);STTRP=cell(100,360);STNRP=cell(100,360);STFRM=cell(100,360);% the initialization of {} matrix.
SFI=Static_field_density;%call the function Static_field_density.m to get the Static field density(SFI)
SPN=Pedes_Create_Exponential_Distribution;%call this function to create the pedes in exponential distribution
TNUM=sum(sum(SPN));STRM=cell(TNUM,9,360);
%r means the seconds,total is 360.There is the beginning of the timing cycling.
for r=1:360
%confirm the position of initial pedes.
[FullEntry,PM,SM,PP,SVM,SDM]=Comfirm_Initial_Pedes(SM,SPN,r,PM,PP,SVM,SDM);
% to store all the pedestrians'SVM.SDM.PM into STSVM.STSDM.STPM
for k=1:TNUM
if k<=length(PP)
STSVM(k,r)=SVM(k);STSDM(k,r)=SDM(k);STPM(k,r)=PM(k);
else
STSVM(k,r)=NaN;STSDM(k,r)=NaN;STPM(k,r)=NaN;
end
end
%the valuation of The Dynamic Field Intensity(DFI)
CFI=Get_CFI(SM,SFI);
%the NaSch model to update Speed
SVM=NaSch_model(SVM,PP,PM);
%Each existing peds would get a row in TRP. TRP means the total route matrix.
TRP=Get_TRP(PP,PM,SVM,SDM);
%use the SM to get FRP and NRP and sort them by CFI
[FRP,NRP]=Get_FRP_NRP(TRP,SM,CFI);
%now we will define the RM for NRP
RM=Get_RM(FRP,PM,TNUM,SVM,SDM);
%calculate the probability of choosing NRP or FRP
PFRP=Get_Probability(FRP,NRP,CFI);
%determine the final set
SFRP=randperm(length(PP),round(PFRP*length(PP)))';%SFRP means the subscript of the people who belongs to the FRP final.
%find the first man in RM.
[fsub,~]=find((cellfun(@isempty,RM)==1)==0,1);
if isempty(fsub)
fsub=0;
end
%now we should resolve the conflict.
FRM=Solve(PP,fsub,RM,SVM);
%store the FRM into STFRM,the FRP into STFRP,and so on.
for k=1:TNUM
if k<=length(PP)
STFRM(k,r)=FRM(k);STFRP(k,r)=FRP(k);STTRP(k,r)=TRP(k);STNRP(k,r)=NRP(k);
else
STFRM(k,r)={NaN};STFRP(k,r)={NaN};STTRP(k,r)={NaN};STNRP(k,r)={NaN};
end
end
%update the direction of pedes.
SDM=Update_Direction(PP,FRM,PM,SVM,SDM);
%now we can update the existing pedestrians'status.
for k=1:length(PP)
if ismember(PM(k,1),[1 31 61 91])%that means the pedes arrive to the destination.
SDM(k)= 0;SVM(k)=0;PM(k)= PM(k);
elseif ismember(k,SFRP)&&~isempty(FRM{k})
SM(PM(k))=0;%as the people will depart in this condition,so the SM in this place should be 0.
PM(k)=FRM{k}(1,end);%the latest PM is followed with the end of FRM.
else %that means the pede now has no place to go or choose the NRP.
PM(k)= PM(k);
end
end
SM(PM)=1;%update the SM using PM.
SM([1,31,61,91])=0;%the entry should always keep empty.
STSM(:,:,r)=SM;%the STSM is used to storing the SM at the end of each second.
STRM(:,:,r)=RM; %the STRM is used to storing the RM at the end of each second.
end
toc;
end
function SFI=Static_field_density
clear
SFI=200*ones(32,12);
SFI(2,2:5)=0;
SFI([1,32],:)=NaN;
SFI(:,[1,12])=NaN;
for k=2:11
for s=2:31
SFI=Assign(SFI,s,k);
end
end
SFI=SFI(2:31,2:11);
end
function SFI=Assign(SFI,s,k)
if SFI(s,k+1)>SFI(s,k)+1
SFI(s,k+1)=SFI(s,k)+1;
end
if SFI(s,k-1)>SFI(s,k)+1
SFI(s,k-1)=SFI(s,k)+1;
end
if SFI(s-1,k)>SFI(s,k)+1
SFI(s-1,k)=SFI(s,k)+1;
end
if SFI(s+1,k)>SFI(s,k)+1
SFI(s+1,k)=SFI(s,k)+1;
end
if SFI(s-1,k-1)>SFI(s,k)+1.5
SFI(s-1,k-1)=SFI(s,k)+1.5;
end
if SFI(s-1,k+1)>SFI(s,k)+1.5
SFI(s-1,k+1)=SFI(s,k)+1.5;
end
if SFI(s+1,k-1)>SFI(s,k)+1.5
SFI(s+1,k-1)=SFI(s,k)+1.5;
end
if SFI(s+1,k+1)>SFI(s,k)+1.5
SFI(s+1,k+1)=SFI(s,k)+1.5;
end
end
function [SPN]=Pedes_Create_Exponential_Distribution
%this function can create the 1000*1 matrix indicate the time interval for each pedestrian which submit to negative exponential distribution
SPN=zeros(360,1);
CSTime=cumsum(exprnd(360/100,100,1));
for k=1:360
SPN(k)=length(find(CSTime>=k-1&CSTime<k));%SPN means the number of pedestrians for every second.
end
end
function [FullEntry,PM,SM,PP,SVM,SDM]=Comfirm_Initial_Pedes(SM,SPN,r,PM,PP,SVM,SDM)
%confirm the position of initial pedes.
AllNoPed=find(SM==0);FullEntry=zeros(360,1);
EnNoPed=AllNoPed(AllNoPed>210);%210=30*(10-3),if 30 and 10 is changed,210 should also be changed.EnNoPed means the entrys which don't have pedestrians.
if EnNoPed~=0
IPos=randsrc(SPN(r),1,EnNoPed');%IPos means the initial pedes' position.
PM=[PM;IPos];%PM means the position matrix,the people previous in time would be in front at PM.
SM(IPos)=1;
PP=(1:(length(PP)+SPN(r)))';%PP is the Pedestrians Pool.
%confirm the speed of initial pedes.
ISV=randsrc(SPN(r),1,[1 2 3]);
SVM=[SVM;ISV];
ISD=3*ones(SPN(r),1);%ISD=randsrc(SPN(i),1,[1 2 3]);
SDM=[SDM;ISD];
else
FullEntry(r)=SPN(r);
end
end
function CFI=Get_CFI(SM,SFI)
DFI=zeros(30,10);
for k=1:10
for s=1:30
DFI(s,k)=(max(s-1,k-1)*size(find(SM(1:s,1:k)==1),1)*size(find(SM(1:s,1:k)==1),2))/(s*k);
end
end
CFI=DFI+SFI;%the valuation of The Comprehensive Field Intensity(CFI)
end
function [SVM] = NaSch_model(SVM,PP,PM)
%the NaSch model to update Speed
log=SVM<3;
SVM=SVM+log;
Rdde=randperm(length(PP),round(0.1*length(SVM)))';%0.1 means the p in NaSch model,Rdde means the Random decrease subscript of pedes.
SVM(Rdde)=SVM(Rdde)-1;
%if the pede has arrived to the exit,then we should put his speed to 0.
SVM(ismember(PM,[1 31 61 91]))=0;
end
function [TRP]=Get_TRP(PP,PM,SVM,SDM)
%j means the existing pedes in the field now.the following is to find out the ITRP.
%Each existing peds would get a row in ITRP.
%ITRP means the inital total route matrix.
ITRP=cell(length(PP),1);
if ~isempty(PP)
for k=1:length(PP)
if SVM(k)==0
ITRP(k,1)={PM(k,1)};
elseif ismember(PM(k,1),[1 31 61 91])
ITRP(k,1)={PM(k,1)};
elseif SVM(k,1)==3&&SDM(k,1)==1
ITRP(k,1)={[PM(k,1) PM(k,1)-1 PM(k,1)-2 PM(k,1)-3 PM(k,1)-30 PM(k,1)-31 PM(k,1)-32 PM(k,1)-60]};
elseif SVM(k,1)==2&&SDM(k,1)==1
ITRP(k,1)={[PM(k,1) PM(k,1)-1 PM(k,1)-2 PM(k,1)-30 PM(k,1)-31]};
elseif SVM(k,1)==1&&SDM(k,1)==1
ITRP(k,1)={[PM(k,1) PM(k,1)-1]};
elseif SVM(k,1)==3&&SDM(k,1)==2
ITRP(k,1)={[PM(k,1) PM(k,1)-1 PM(k,1)-2 PM(k,1)-30 PM(k,1)-31 PM(k,1)-32 PM(k,1)-60 PM(k,1)-61 PM(k,1)-62]};
elseif SVM(k,1)==2&&SDM(k,1)==2
ITRP(k,1)={[PM(k,1) PM(k,1)-1 PM(k,1)-30 PM(k,1)-31]};
elseif SVM(k,1)==1&&SDM(k,1)==2
ITRP(k,1)={PM(k,1)};
elseif SVM(k,1)==3&&SDM(k,1)==3
ITRP(k,1)={[PM(k,1) PM(k,1)-30 PM(k,1)-60 PM(k,1)-90 PM(k,1)-1 PM(k,1)-31 PM(k,1)-61 PM(k,1)-2]};
elseif SVM(k,1)==2&&SDM(k,1)==3
ITRP(k,1)={[PM(k,1) PM(k,1)-30 PM(k,1)-60 PM(k,1)-1 PM(k,1)-31]};
elseif SVM(k,1)==1&&SDM(k,1)==3
ITRP(k,1)={[PM(k,1) PM(k,1)-30]};
end
end
else
ITRP={[]};
end
%we should choose from ITRP to get the final matrix TRP
TRP=cellfun(@(x) x(x>0&x<=300),ITRP, 'UniformOutput',false);
end
function [FRP,NRP]=Get_FRP_NRP(TRP,SM,CFI)
%use the SM to get FRP and NRP and sort them by CFI
if ~isempty(TRP)
FRP=cell(length(TRP),1);MNRP=cell(length(TRP),1);NRP=cell(length(TRP),1);
for k=1:length(TRP)
FRP(k,1)={TRP{k}(SM(TRP{k})==0)};
MNRP(k,1)={TRP{k}(SM(TRP{k})==1)};
[~,RA0]=sort(CFI(FRP{k}));%sort FRP by CFI (ascend),RA0 is the origin rank.
FRP(k,1)={FRP{k}(RA0)};
if ~isempty(FRP{k})
NRP(k,1)={MNRP{k}(CFI(MNRP{k})<min(CFI(FRP{k})))};%the CFI of NRP should less than min(CFI(FRP))
else
NRP(k,1)={[]};
end
[~,RA1]=sort(CFI(NRP{k}));%sort NRP by CFI (ascend),RA1 is the origin rank.
NRP(k,1)={NRP{k}(RA1)};
end
else
FRP={[]};
end
end
function [RM]=Get_RM(FRP,PM,TNUM,SVM,SDM)
%now we will define the RM for NRP
RM=cell(length(FRP),9);
for k=1:length(FRP)
for s=1:length(FRP{k})
if (FRP{k}(s)-PM(k))==0%keep static
RM(k,s)={[]};
elseif ismember(PM(k),[1 31 61 91])
RM(k,s)={[]};
elseif SVM(k)==1&&SDM(k)==3&& (FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==1&&SDM(k)==1&& (FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==2&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==2&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-2
RM(k,s)={[PM(k) PM(k)-1 PM(k)-2]};
elseif SVM(k)==2&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==2&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-31
RM(k,s)={[PM(k) PM(k)-31]};
elseif SVM(k)==2&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-31
RM(k,s)={[PM(k) PM(k)-31]};
elseif SVM(k)==2&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==2&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==2&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==2&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==2&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-31
RM(k,s)={[PM(k) PM(k)-31]};
elseif SVM(k)==2&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-60
RM(k,s)={[PM(k) PM(k)-30 PM(k)-60]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-2
RM(k,s)={[PM(k) PM(k)-1 PM(k)-2]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-3
RM(k,s)={[PM(k) PM(k)-1 PM(k)-2 PM(k)-3]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-60
RM(k,s)={[PM(k) PM(k)-30 PM(k)-60]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-31
RM(k,s)={[PM(k) PM(k)-31]};
elseif SVM(k)==3&&SDM(k)==1&&(FRP{k}(s)-PM(k))==-32
RM(k,s)={[PM(k) PM(k)-1 PM(k)-32]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-31
RM(k,s)={[PM(k) PM(k)-31]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-62
RM(k,s)={[PM(k) PM(k)-31 PM(k)-62]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-2
RM(k,s)={[PM(k) PM(k)-1 PM(k)-2]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-60
RM(k,s)={[PM(k) PM(k)-30 PM(k)-60]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-32
RM(k,s)={[PM(k) PM(k)-31 PM(k)-32]};
elseif SVM(k)==3&&SDM(k)==2&&(FRP{k}(s)-PM(k))==-61
RM(k,s)={[PM(k) PM(k)-31 PM(k)-61]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-30
RM(k,s)={[PM(k) PM(k)-30]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-60
RM(k,s)={[PM(k) PM(k)-30 PM(k)-60]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-90
RM(k,s)={[PM(k) PM(k)-30 PM(k)-60 PM(k)-90]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-31
RM(k,s)={[PM(k) PM(k)-31]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-61
RM(k,s)={[PM(k) PM(k)-30 PM(k)-61]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-1
RM(k,s)={[PM(k) PM(k)-1]};
elseif SVM(k)==3&&SDM(k)==3&&(FRP{k}(s)-PM(k))==-2
RM(k,s)={[PM(k) PM(k)-1 PM(k)-2]};
else
RM(k,s)={10000};
end
end
RM(k,(length(FRP{k})+1):9)={[]};
end
RM((length(FRP)+1):TNUM,9)= {[]};
end
function [PFRP]=Get_Probability(FRP,NRP,CFI)
%calculate the probability of choosing NRP or FRP
SFRP=0;LFRP=0;SNRP=0;LNRP=0;
if ~isempty(FRP)
for k=1:length(FRP)
SFRP=SFRP+sum(CFI(FRP{k}));
LFRP=LFRP+length(FRP{k});
end
AFRP=SFRP/LFRP;
else
LFRP=0;
end
if ~isempty(FRP)
for k=1:length(NRP)
SNRP=SNRP+sum(CFI(NRP{k}));
LNRP=LNRP+length(NRP{k});
end
ANRP=SNRP/LNRP;
else
LNRP=0;
end
if LFRP~=0&&LNRP~=0
SURP=1/((2*(1/AFRP))+(1/ANRP));
PFRP=SURP*2*(1/AFRP);
elseif LFRP~=0&&LNRP==0
PFRP=1;
else
PFRP=0;
end
end
function [FRM]=Solve(PP,fsub,RM,SVM)
TCR=ones(length(PP),1);FRM=cell(length(PP),1);
if fsub~=0&&length(PP)>=fsub+1%for the case that has at least 2 PP.
for k=fsub:length(PP)%k means the main man now.
for s=1:9%s can indicate the different destination of the main man.
for t=(fsub+1):length(PP) %t can indicate the different pedes behind the main man and the first man.
CR=Conflict_Resolve_Main(RM{k,s},RM{t,1},SVM(k),SVM(t));
TCR(t-fsub,1)=CR;
end
if k~=fsub%for the case that it's not the first man.
for v=fsub:(k-1) %the main man should compare with the man before him.
CR=isempty(intersect(RM{k,s},FRM{v,1}));
TCR(length(PP)+v-2*fsub+1,1)=CR;
end
end
if ~ismember(0,TCR)
FRM{k,1}=RM{k,s};
break;
else
FRM{k,1}=[];
end
end
end
elseif fsub==length(PP)&&~isempty(PP) %this means only 1 PP.
FRM{fsub,1}=RM{fsub,1};
elseif fsub==0&&~isempty(PP) %for the case that fsub is 0 but PP is not 0.
FRM{length(PP),1}=[];
end
end
function CR=Conflict_Resolve_Main(A,B,VA,VB)
%A and B is the both side of conflict waiting to be resolved.VA and VB is
%the speed of each side.
[C,a1,b1]=intersect(A,B);
if isempty(C)
CR=1;%CR=1 means A get the passport
elseif a1(1)~=1&&b1(1)~=1
A3=A(2:a1(1))-A(1:a1(1)-1);
DISTa(abs(A3)==1|abs(A3)==30)=1;
DISTa(abs(A3)==31)=1.5;
B3=B(2:b1(1))-B(1:b1(1)-1);
DISTb(abs(B3)==1|abs(B3)==30)=1;
DISTb(abs(B3)==31)=1.5;
CR=(sum(DISTa)/VA)<=(sum(DISTb)/VB);%when A arrive at the same time with B,this function would transfer the passport to A.That means the A is more important.
elseif a1(1)==1
CR=1;
else
CR=0;
end
end
function [SDM]=Update_Direction(PP,FRM,PM,SVM,SDM)
%update the direction of pedes.
for k=1:length(PP)
if isempty(FRM{k})
SDM(k,1)=SDM(k,1);
elseif SVM(k)==2&&SDM(k)==1&&(FRM{k}(end)-PM(k))==-30
SDM(k,1)=3;
elseif SVM(k)==2&&SDM(k)==1&&(FRM{k}(end)-PM(k))==-31
SDM(k,1)=2;
elseif SVM(k)==2&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-1
SDM(k,1)=1;
elseif SVM(k)==2&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-30
SDM(k,1)=3;
elseif SVM(k)==2&&SDM(k)==3&&(FRM{k}(end)-PM(k))==-31
SDM(k,1)=2;
elseif SVM(k)==2&&SDM(k)==3&&(FRM{k}(end)-PM(k))==-1
SDM(k,1)=1;
elseif SVM(k)==3&&SDM(k)==1&&(FRM{k}(end)-PM(k))==-30
SDM(k,1)=3;
elseif SVM(k)==3&&SDM(k)==1&&(FRM{k}(end)-PM(k))==-60
SDM(k,1)=3;
elseif SVM(k)==3&&SDM(k)==1&&(FRM{k}(end)-PM(k))==-31
SDM(k,1)=2;
elseif SVM(k)==3&&SDM(k)==1&&(FRM{k}(end)-PM(k))==-32
SDM(k,1)=2;
elseif SVM(k)==3&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-1
SDM(k,1)=1;
elseif SVM(k)==3&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-2
SDM(k,1)=1;
elseif SVM(k)==3&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-30
SDM(k,1)=3;
elseif SVM(k)==3&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-60
SDM(k,1)=3;
elseif SVM(k)==3&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-32
SDM(k,1)=1;
elseif SVM(k)==3&&SDM(k)==2&&(FRM{k}(end)-PM(k))==-61
SDM(k,1)=3;
elseif SVM(k)==3&&SDM(k)==3&&(FRM{k}(end)-PM(k))==-31
SDM(k,1)=2;
elseif SVM(k)==3&&SDM(k)==3&&(FRM{k}(end)-PM(k))==-61
SDM(k,1)=2;
elseif SVM(k)==3&&SDM(k)==3&&(FRM{k}(end)-PM(k))==-1
SDM(k,1)=1;
elseif SVM(k)==3&&SDM(k)==3&&(FRM{k}(end)-PM(k))==-2
SDM(k,1)=1;
end
end
end
评论