www.gusucode.com > 《LINGO软件及应用》程序和数据 > 《LINGO软件及应用》程序和数据/12第12章/Lanli12_9_9.m

    clc, clear
xy0=xlsread('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',1,'B2:C583'); %读入第1个表单中的路口坐标数据
nj=xlsread('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',1,'E2:E583');
uv=xlsread('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',2,'A2:B929'); %读入第2个表单中的道路顶点编号
a=zeros(582);
pt=xlsread('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls',3,'B2:B81'); %读入第3个表单中平台位置标号
n=length(pt); %警力平台个数
for i=1:length(uv)
a(uv(i,1),uv(i,2))=sqrt((xy0(uv(i,1),1)-xy0(uv(i,2),1))^2+(xy0(uv(i,1),2)-xy0(uv(i,2),2))^2)*0.1;
end
a=a+a'; %构造实对称的邻接矩阵
a=tril(a); a=sparse(a);
b=graphallshortestpaths(a,'directed',0);%求所有节点对之间的最短距离;
speed=60; sp=speed/60; fail=1; t=0; 
tic %计时开始
b32=b(32,:); %32号节点到其他所有节点的距离
while fail~=0;
ptD=[];%记录负责围堵的平台
B=[]; numD=0;
A=find(b32<(3+t)*sp), numA=length(A); %确定tmin无法封锁的节点集合A及个数
for i=1:numA
    indB1=find(uv(:,1)==A(i)); %求A中节点的相邻节点地址
    indB2=find(uv(:,2)==A(i)); %求A中节点的相邻节点地址
    B=[B,setdiff(uv(indB1,2)',A)]; %B集合为A集合的相邻节点构成的集合
    B=[B,setdiff(uv(indB2,1)',A)];
end
B=unique(B),numB=length(B); %去掉重复元素，并统计B中元素个数
for j=1:numB
    [tt,indT]=min(b(B(j),pt)/sp); %平台到该节点的最短时间
    if tt<t
        numD=numD+1;
        ptD(numD)=pt(indT); %B(j)由平台pt(indT)负责围堵
        %pt(indT)=[]; %删除该平台编号，以后不能使用
    end
end
if numD==numB
    fail=0;
else
    t=t+0.1;
end
end
B,t %显示围堵的节点和时间
c32=b(32,B); %提出嫌疑犯出发点到围堵点之间的最短距离 
d=b(B,pt); %提出围堵点与平台之间的最短距离
xlswrite('Ldata1295.xlsx',c32);
xlswrite('Ldata1295.xlsx',d,1,'A3') %把矩阵d写入表单Sheet1中A3开始的单元格中