全国大学生数学建模竞赛2011B题交巡警服务平台的设置与调度MATLAB程序与LINGO程序

目录

一、第1问

二、第2问

三、第3问 

四、第4问 

五、第5问

代码中用到的数据

链接：https://pan.baidu.com/s/1YqLXLp6rhgl4_i8z0pYHzQ?pwd=6666
提取码：6666

---

一、第1问

%第1问MATLAB程序
%% 对应计算两点间距离（A区）
a=xlsread('2011B.xls','全市交通路口节点数据','B2:C93');
b=xlsread('2011B.xls','全市交通路口的路线','A2:B144');
c=zeros(143,2);
d=zeros(143,2);
e=zeros(143,1);
for i=1:143
c(i,:)=a(b(i,1),1:2);
end

for i=1:143
    if b(i,2)>92 
        continue
    end
d(i,:)=a(b(i,2),1:2);
end

for i=1:143
    if b(i,2)>92 
       e(i,1)=0;
        continue
    end
e(i,1)=sqrt((c(i,1)-d(i,1))^2+(c(i,2)-d(i,2))^2);
end

%% 画图（A区）
c1=xlsread('2011B.xls','全市交通路口的路线','L2:M141');
cc1=xlsread('2011B.xls','全市交通路口的路线','J2:J141');
d1=xlsread('2011B.xls','全市交通路口的路线','N2:O141'); 
dd1=xlsread('2011B.xls','全市交通路口的路线','K2:K141');
for i=1:140 
plot([c1(i,1) d1(i,1)],[c1(i,2) d1(i,2)],'linewidth',1.5)
hold on
plot(c1(i,1),c1(i,2),'ro','linewidth',1.5)
text(c1(i,1),c1(i,2),[num2str(cc1(i))],'FontSize',12);
hold on
plot(d1(i,1),d1(i,2),'ro','linewidth',1.5)
text(d1(i,1),d1(i,2),[num2str(dd1(i))],'FontSize',12);
xlabel('x','FontSize',26);
ylabel('y','FontSize',26);
set(gca,'FontSize',26,'linewidth',1);
end

%% 灰度处理（A区）
f=imread('D:\数学建模\2021暑期\2011B模拟3\第1问\1.jpg');
g=rgb2gray(f);
imshow(g);

%% floyd算法
function [D,path,min1,path1]=floyd(a,start,terminal)
D=a;n=size(D,1);path=zeros(n,n);
for i=1:n
    for j=1:n
        if D(i,j)~=inf
            path(i,j)=j;
        end
    end
end
for k=1:n
    for i=1:n
        for j=1:n
            if D(i,k)+D(k,j)<D(i,j)
                D(i,j)=D(i,k)+D(k,j);
                path(i,j)=path(i,k);
            end
        end
    end
end
if  nargin==92
    min1=D(start,terminal);
    m(1)=start;
    i=1;
    path1=[];
    while path(m(i),terminal)~=terminal
        k=i+1;
        m(k)=path(m(i),terminal);
        i=i+1;
    end
    m(i+1)=terminal;
    path1=m;
end

%% 对应构造邻接矩阵并计算最小值（A区）
a=xlsread('2011B.xls','处理过的全市交通路口的路线','A2:A141');
b=xlsread('2011B.xls','处理过的全市交通路口的路线','B2:B141');
d=xlsread('2011B.xls','处理过的全市交通路口的路线','G2:G141');
W=zeros(92);
for i=1:140
W(a(i),b(i))=d(i);
W(b(i),a(i))=d(i);
end
ind=find(W==0);%找出A中所有为0的位置标号
W(ind)=inf;
for i=1:92
W(i,i)=0;
end
[P,path]=floyd(W);
xlswrite('各点间最短距离2011B.xlsx',P)

%% 给交巡警平台分配最近点为管辖范围
d=xlsread('各点间最短距离2011B.xlsx','Sheet1','A1:T92');
d=0.1*d;
x=zeros(92,1);
y=zeros(92,1);
m=zeros(92,1);
for i=1:92
    m(i)=min(d(i,:));
    [x(i),y(i)]=find(d(i,:)==m(i));
end

二、第2问

%第2问MATLAB程序
%% 根据距离矩阵对应分配点到交巡警平台的距离
d=xlsread('各点间最短距离2011B.xlsx','Sheet1','A1:CN92');
x=xlsread('2011B.xls','全市区出入口的位置','J2:J14');
y=xlsread('2011B.xls','全市区出入口的位置','K1:AD1');
D=zeros(13,20);
for i=1:13 %行
    for j=1:20 %列
        D(i,j)=d(x(i),y(j));
    end
end
D=0.1*D;

%第2问LINGO程序
model:  
sets:
Plat/1..20/;
Kou/1..13/:T;
Assign(Plat,Kou):dis,x,Time;
endsets
data:
v=1000;
dis=@file('D:\Downloads\Matlab2020a\bin\2.txt'); 
!打开平台和路口的距离矩阵20*13;
@text('D:\Downloads\Matlab2020a\bin\2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';'); 
Enddata
 
!min=aver;  
min=TT;
aver=@sum(Kou(j):T(j))/@size(Kou);!平均时间;
!min=aver;
TT<=8.0156;  ! 非线性化转化为线性化，时间约束转化为约束条件
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v); !到第j个路口的时间;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);
@for(kou(j):TT>=T(j));
@for(Plat(i):@sum(Kou(j):x(i,j))<=1);  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Assign(i,j):@bin(x(i,j))); 
End

三、第3问 

%第3问MATLAB程序
%% 找前2个点
d1=xlsread('各点间最短距离2011B.xlsx','Sheet1','A21:T92');
d2=xlsread('各点间最短距离2011B.xlsx','Sheet1','U21:CN92');
d1=0.1*d1;
d2=0.1*d2;
dq1=zeros(72,1);
dq2=zeros(72,1);
std0=zeros(5112,1);

x=zeros(5112,1);
y=zeros(5112,1);
m=zeros(72,5112);
mm=zeros(72,1);
d0=zeros(72,22);
i=0;
ii=0;

for zz1=1:72 %取列向量
    dq1(:,1)=d2(:,zz1);
    for zz2=1:72 %取列向量
        if zz1==zz2
            continue
        end
        dq2(:,1)=d2(:,zz2);
        d0=[d1,dq1,dq2];
        i=i+1;
        for z=1:72
        m(z,i)=min(d0(z,:));
        end
    end
end

for zz3=1:5112
    if max(m(:,zz3))<3
    [x(zz3,1),y(zz3,1)]=find(m(:,zz3)==max(m(:,zz3)));
    ii=ii+1;
    mm(:,ii)=m(:,zz3);
    end
end

%% 代入前4个点
d1=xlsread('各点间最短距离2011B.xlsx','Sheet1','A1:T92');
dq1=xlsread('各点间最短距离2011B.xlsx','Sheet1','A29:CN29');
dq2=xlsread('各点间最短距离2011B.xlsx','Sheet1','A39:CN39');
dq3=xlsread('各点间最短距离2011B.xlsx','Sheet1','A61:CN61');
dq4=xlsread('各点间最短距离2011B.xlsx','Sheet1','A92:CN92');
d1=0.1*d1;dq1=0.1*dq1';dq2=0.1*dq2';dq3=0.1*dq3';dq4=0.1*dq4';
d0=[d1,dq1,dq2,dq3,dq4];
x=zeros(92,1);
y=zeros(92,1);
m=zeros(92,1);

for z=1:92
        m(z,1)=min(d0(z,:));
        [x(z),y(z)]=find(d0(z,:)==m(z,1));
end

a=xlsread('2011B.xls','1.3','E2:E93');%92个路口最近的对应路口工作量
b=xlsread('2011B.xls','1.3','C2:C93');%92个路口最近的对应交巡警平台标号
c=zeros(24,1);
for i=1:24
    n=0;
    for j=1:92
        if b(j,1)==i 
        n=n+a(j,1);
        end
    end
    c(i,1)=n;
end
std4=std(c(:,1));

%% 找第5个点
d2=xlsread('各点间最短距离2011B.xlsx','Sheet5','A1:BP92');%再设置寻找第5个平台
d2=0.1*d2;
dq5=zeros(92,1);
x=zeros(92,68);
y=zeros(92,68);
m=zeros(92,68);
i=0;
std5=zeros(68,1);
for zz1=1:68 %取列向量
        dq5(:,1)=d2(:,zz1);
        d0=[d1,dq1,dq2,dq3,dq4,dq5];
        i=i+1;
        for zz2=1:92
        m(zz2,i)=min(d0(zz2,:));
        [x(zz2,i),y(zz2,i)]=find(d0(zz2,:)==m(zz2,i));
        end
end

dd=xlsread('各点间最短距离2011B.xlsx','Sheet8','A1:BP92');%将序号25替换为对应标号
for nn1=1:68
    for nn2=1:92
        if (dd(nn2,nn1)==25)
            [x0,y0]=find((dd(:,nn1)==25));
            dd(x0,nn1)=(24+nn1);
        end
    end
end

a=xlsread('2011B.xls','1.3','E2:E93');%92个路口最近的对应路口工作量
b=xlsread('各点间最短距离2011B.xlsx','Sheet9','A1:BP92');%92个路口最近的对应交巡警平台标号
c=zeros(25,68);
for jj=1:68
   for i=1:24
    n=0;
        for j=1:92
        if b(j,jj)==i 
        n=n+a(j,1);
        end
        end
    c(i,jj)=n;
   end 
end
for zz4=1:68
    std5(zz4)=std(c(:,zz4),1);
end
std55=min(std5);
[x5]=find(std5==min(std5));
x5=x5+24;

%第3问LINGO程序
model:
sets:
Plat/1..92/:Y;
Kou/1..92/:T;
Assign(Plat,Kou):dis,x;
endsets
data:
v=1000;
dis=@ole('D:\Downloads\Matlab2020a\bin\各点间最短距离2011B1.xlsx','A1:CN92'); !打开92*92个路口间的距离矩阵;
@text('D:\Downloads\Matlab2020a\bin\A1.txt')=@writefor(Plat(i)|Y(i)#GT#0:'Y(',i,')=',Y(i),';');
@text('D:\Downloads\Matlab2020a\bin\A2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';');
k=5;  !增加平台数;
enddata
min=k;
TT<3;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);!到第j个路口的时间;
@for(kou(j):TT>=T(j));
@for(Plat(i):@for(Kou(j):x(i,j)<=Y(i)));  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Plat(i)|i#LE#20:Y(i)=1);  !前20个路口已经为平台;
@sum(Plat(i):Y(i))=20+k;
@for(Assign(i,j):@bin(x(i,j))); 
@for(Plat(i):@bin(Y(i))); 
End

四、第4问 

%第4问MATLAB程序
%% 对应计算两点间距离（全市）
a=xlsread('2011B.xls','全市交通路口节点数据','B2:C583');
b=xlsread('2011B.xls','全市交通路口的路线','A2:B929');
c=zeros(928,2);%起点坐标
d=zeros(928,2);%终点坐标
e=zeros(928,1);
for i=1:928
c(i,:)=a(b(i,1),1:2);
end

for i=1:928
d(i,:)=a(b(i,2),1:2);
end

for i=1:928
e(i,1)=sqrt((c(i,1)-d(i,1))^2+(c(i,2)-d(i,2))^2);
end

%% 对应构造邻接矩阵并计算最小值（全市）
a1=xlsread('2011B.xls','全市交通路口的路线','A2:A929');
b1=xlsread('2011B.xls','全市交通路口的路线','B2:B929');
d1=xlsread('2011B.xls','全市交通路口的路线','G2:G929');
W=zeros(582);
for i=1:928
W(a1(i),b1(i))=d1(i);
W(b1(i),a1(i))=d1(i);
end
ind=find(W==0);%找出所有为0的位置标号
W(ind)=inf;
for i=1:582
W(i,i)=0;
end
[P,path]=floyd(W);
xlswrite('各点间最短距离2011B2.xlsx',W,'邻接矩阵')
xlswrite('各点间最短距离2011B2.xlsx',P,'最短路径')

%% 数据处理后可代入LINGO（单位处理）
a=xlsread('各点间最短距离2011B.xlsx','Sheet1','A1:VJ582');
a=a*100;
xlswrite('各点间最短距离2011B2.xlsx',a)

%% 对应构造邻接矩阵并计算最小值（B区）
a2=xlsread('2011B.xls','处理过的全市交通路口的路线B','A2:A114');
b2=xlsread('2011B.xls','处理过的全市交通路口的路线B','B2:B114');
d2=xlsread('2011B.xls','处理过的全市交通路口的路线B','G2:G114');
WB=zeros(73);
for i=1:113
WB(a2(i)-92,b2(i)-92)=d2(i);
WB(b2(i)-92,a2(i)-92)=d2(i);
end
indB=find(WB==0);%找出所有为0的位置标号
WB(indB)=inf;
for i=1:73
WB(i,i)=0;
end
[PB,pathB]=floyd(WB);
xlswrite('各点间最短距离2011B3.xlsx',PB,'B')

%% 对应构造邻接矩阵并计算最小值（C区）
a3=xlsread('2011B.xls','处理过的全市交通路口的路线C','A2:A248');
b3=xlsread('2011B.xls','处理过的全市交通路口的路线C','B2:B248');
d3=xlsread('2011B.xls','处理过的全市交通路口的路线C','G2:G248');
WC=zeros(154);
for i=1:247
WC(a3(i)-165,b3(i)-165)=d3(i);
WC(b3(i)-165,a3(i)-165)=d3(i);
end
indC=find(WC==0);%找出所有为0的位置标号
WC(indC)=inf;
for i=1:154
WC(i,i)=0;
end
[PC,pathC]=floyd(WC);
xlswrite('各点间最短距离2011B3.xlsx',PC,'C')

%% 对应构造邻接矩阵并计算最小值（D区）
a4=xlsread('2011B.xls','处理过的全市交通路口的路线D','A2:A77');
b4=xlsread('2011B.xls','处理过的全市交通路口的路线D','B2:B77');
d4=xlsread('2011B.xls','处理过的全市交通路口的路线D','G2:G77');
WD=zeros(52);
for i=1:76
WD(a4(i)-319,b4(i)-319)=d4(i);
WD(b4(i)-319,a4(i)-319)=d4(i);
end
indD=find(WD==0);%找出所有为0的位置标号
WD(indD)=inf;
for i=1:52
WD(i,i)=0;
end
[PD,pathD]=floyd(WD);
xlswrite('各点间最短距离2011B3.xlsx',PD,'D')

%% 对应构造邻接矩阵并计算最小值（E区）
a5=xlsread('2011B.xls','处理过的全市交通路口的路线E','A2:A152');
b5=xlsread('2011B.xls','处理过的全市交通路口的路线E','B2:B152');
d5=xlsread('2011B.xls','处理过的全市交通路口的路线E','G2:G152');
WE=zeros(103);
for i=1:151
WE(a5(i)-371,b5(i)-371)=d5(i);
WE(b5(i)-371,a5(i)-371)=d5(i);
end
indE=find(WE==0);%找出所有为0的位置标号
WE(indE)=inf;
for i=1:103
WE(i,i)=0;
end
[PE,pathE]=floyd(WE);
xlswrite('各点间最短距离2011B3.xlsx',PE,'E')

%% 对应构造邻接矩阵并计算最小值（F区）
a6=xlsread('2011B.xls','处理过的全市交通路口的路线F','A2:A169');
b6=xlsread('2011B.xls','处理过的全市交通路口的路线F','B2:B169');
d6=xlsread('2011B.xls','处理过的全市交通路口的路线F','G2:G169');
WF=zeros(108);
for i=1:168
WF(a6(i)-474,b6(i)-474)=d6(i);
WF(b6(i)-474,a6(i)-474)=d6(i);
end
indF=find(WF==0);%找出所有为0的位置标号
WF(indF)=inf;
for i=1:108
WF(i,i)=0;
end
[PF,pathF]=floyd(WF);
xlswrite('各点间最短距离2011B3.xlsx',PF,'F')

%第4问LINGO程序
1.寻找B区需要加入的工作点LINGO程序：
model:
sets:
Plat/1..73/:Y;
Kou/1..73/:T;
Assign(Plat,Kou):dis,x;
endsets
data:
v=10;
dis=@ole('D:\Downloads\Matlab2020a\bin\各点间最短距离2011B3.xlsx','A1:BU73'); !打开73*73个路口间的距离矩阵;
@text('D:\Downloads\Matlab2020a\bin\B1.txt')=@writefor(Plat(i)|Y(i)#GT#0:'Y(',i,')=',Y(i),';');
@text('D:\Downloads\Matlab2020a\bin\B2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';');
k=2;  !增加平台数;
enddata
min=k;
TT<3;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);!到第j个路口的时间;
@for(kou(j):TT>=T(j));
@for(Plat(i):@for(Kou(j):x(i,j)<=Y(i)));  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Plat(i)|i#LE#8:Y(i)=1);  !前8个路口已经为平台;
@sum(Plat(i):Y(i))=8+k;
@for(Assign(i,j):@bin(x(i,j))); 
@for(Plat(i):@bin(Y(i))); 
End

2.寻找C区需要加入的工作点LINGO程序：
model:
sets:
Plat/1..154/:Y;
Kou/1..154/:T;
Assign(Plat,Kou):dis,x;
endsets
data:
v=10;
dis=@ole('D:\Downloads\Matlab2020a\bin\各点间最短距离2011B3.xlsx','A1:EX154'); !打开154*154个路口间的距离矩阵;
@text('D:\Downloads\Matlab2020a\bin\C1.txt')=@writefor(Plat(i)|Y(i)#GT#0:'Y(',i,')=',Y(i),';');
@text('D:\Downloads\Matlab2020a\bin\C2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';');
k=15;  !增加平台数;
enddata
min=k;
TT<3;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);!到第j个路口的时间;
@for(kou(j):TT>=T(j));
@for(Plat(i):@for(Kou(j):x(i,j)<=Y(i)));  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Plat(i)|i#LE#17:Y(i)=1);  !前17个路口已经为平台;
@sum(Plat(i):Y(i))=17+k;
@for(Assign(i,j):@bin(x(i,j))); 
@for(Plat(i):@bin(Y(i))); 
End

3.寻找D区需要加入的工作点LINGO程序：
model:
sets:
Plat/1..52/:Y;
Kou/1..52/:T;
Assign(Plat,Kou):dis,x;
endsets
data:
v=10;
dis=@ole('D:\Downloads\Matlab2020a\bin\各点间最短距离2011B3.xlsx','A1:AZ52'); !打开52*52个路口间的距离矩阵;
@text('D:\Downloads\Matlab2020a\bin\D1.txt')=@writefor(Plat(i)|Y(i)#GT#0:'Y(',i,')=',Y(i),';');
@text('D:\Downloads\Matlab2020a\bin\D2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';');
k=8;  !增加平台数;
enddata
min=k;
TT<3;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);!到第j个路口的时间;
@for(kou(j):TT>=T(j));
@for(Plat(i):@for(Kou(j):x(i,j)<=Y(i)));  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Plat(i)|i#LE#9:Y(i)=1);  !前9个路口已经为平台;
@sum(Plat(i):Y(i))=9+k;
@for(Assign(i,j):@bin(x(i,j))); 
@for(Plat(i):@bin(Y(i))); 
End

4.寻找E区需要加入的工作点LINGO程序：
model:
sets:
Plat/1..103/:Y;
Kou/1..103/:T;
Assign(Plat,Kou):dis,x;
endsets
data:
v=10;
dis=@ole('D:\Downloads\Matlab2020a\bin\各点间最短距离2011B3.xlsx','A1:CY103'); !打开103*103个路口间的距离矩阵;
@text('D:\Downloads\Matlab2020a\bin\E1.txt')=@writefor(Plat(i)|Y(i)#GT#0:'Y(',i,')=',Y(i),';');
@text('D:\Downloads\Matlab2020a\bin\E2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';');
k=17;  !增加平台数;
enddata
min=k;
TT<3;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);!到第j个路口的时间;
@for(kou(j):TT>=T(j));
@for(Plat(i):@for(Kou(j):x(i,j)<=Y(i)));  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Plat(i)|i#LE#9:Y(i)=1);  !前9个路口已经为平台;
@sum(Plat(i):Y(i))=9+k;
@for(Assign(i,j):@bin(x(i,j))); 
@for(Plat(i):@bin(Y(i))); 
End

5.寻找F区需要加入的工作点LINGO程序：
model:
sets:
Plat/1..108/:Y;
Kou/1..108/:T;
Assign(Plat,Kou):dis,x;
endsets
data:
v=10;
dis=@ole('D:\Downloads\Matlab2020a\bin\各点间最短距离2011B3.xlsx','A1:DD108'); !打开108*108个路口间的距离矩阵;
@text('D:\Downloads\Matlab2020a\bin\F1.txt')=@writefor(Plat(i)|Y(i)#GT#0:'Y(',i,')=',Y(i),';');
@text('D:\Downloads\Matlab2020a\bin\F2.txt')=@writefor(Assign(i,j)|x(i,j)#GT#0:'x(',i,',',j,')=',x(i,j),';');
k=13;  !增加平台数;
enddata
min=k;
TT<3;
@for(kou(j):T(j)=@sum(Plat(i):x(i,j)*dis(i,j))/v);!到第j个路口的时间;
@for(kou(j):TT>=T(j));
@for(Plat(i):@for(Kou(j):x(i,j)<=Y(i)));  !第i个平台最多到一个路口;
@for(Kou(j):@sum(Plat(i):x(i,j))=1); !第j个路口恰好有一个平台到达;
@for(Plat(i)|i#LE#11:Y(i)=1);  !前11个路口已经为平台;
@sum(Plat(i):Y(i))=11+k;
@for(Assign(i,j):@bin(x(i,j))); 
@for(Plat(i):@bin(Y(i))); 
End

五、第5问

%第5问MATLAB程序
%% 定义变量
%求给定t时的集合At,Bt,Ct,Pt(可选平台集合)
%给定时间t
t=9.05;  %给定时间
v=10; %警察速度
v1=10;%罪犯速度
pos=32; %罪犯位置

%% 对应构造邻接矩阵并计算最小值（全市）
a1=xlsread('2011B.xls','全市交通路口的路线','A2:A929');
b1=xlsread('2011B.xls','全市交通路口的路线','B2:B929');
d1=xlsread('2011B.xls','全市交通路口的路线','G2:G929');
W=zeros(582);
for i=1:928
W(a1(i),b1(i))=d1(i);
W(b1(i),a1(i))=d1(i);
end
for i=1:582
W(i,i)=0;
end

%% 定义变量
A=xlsread('各点间最短距离2011B2.xlsx','最短路径','A1:VJ582');
P=W;
At=[];
DAt=[];
num1=0;

%% step1 求t分钟内形成的范围At
for i=1:582    
    if A(i,pos)/v1<t        %pos=32; 罪犯逃跑点 A(582, 582)任意两点之间的最短路径
        num1=num1+1; 
        At=[At,i];
        DAt=[DAt,A(i,pos)];
    end      
end
 
hold on 
for i=1:num1
    plot(At(i),At(i),'bo');   %画At点
end

%% step2 求与At相邻的点Bt
Bt=[];
DBt=[];
U=zeros(582,1);
for i=1:num1
    i1=At(i);
    for j=1:582     %对列进行遍历
        if (P(i1,j)~=0)   % 对邻接矩阵进行判断，以i1为起点，j为终点是否有路即是否相邻
            U(j)=1;     % U的第j行为1，代表第j个路口是At的相邻路口
        end  %找出与At(i)相连的点
    end
end
   
for i=1:582   % 对U进行行遍历
    if U(i)==1
        Bt=[Bt,i];
        DBt=[DBt,A(i,pos)];
    end
end
num2=length(Bt);

% for i=1:num2
%     plot(x(Bt(i)),y(Bt(i)),'r+');   %画At点
% end

%% step3 求C(t)=B(t)-A(t)
U=zeros(582,1);     % U(i)==0  指嫌犯不能到达该路口       U(i)==1 指嫌犯能够到达该路口
  
for i=1:num1    % 以At为范围        At 1x71矩阵   num1=71
    k=At(i);        %k 嫌犯所能到达的第i个路口
    U(k)=1;         %将其置为1，代表嫌犯能够到达该路口
end
Ct=[];
DCt=[];
for i=1:num2    % 以Bt为范围    Bt 1x95矩阵       num2=95
    i1=Bt(i);           % 提取Bt中第i个与At相邻的路口即为i1
    if U(i1)==0     % 如果i1路口嫌犯不能到达，则为最小外包围圈的路口路口
        Ct=[Ct,i1];     % 将i1路口存入外包围圈Ct中
    end
end
num3=length(Ct);
% for i=1:num3
%     DCt=[DCt,
DB=xlsread('2011B.xls','全市交通路口的路线','A2:B929');
DC=xlsread('2011B.xls','全市交巡警平台','B2:B81');
kp=length(DB);  % 获得连接线总数，附件sheet2中可得     kp=928
px=zeros(2,1); py=zeros(2,1);   % 建立矩阵存储两点的x与y坐标
hold on
 
for i=1:kp  % kp=928
    px(1)=DB(i,1); py(1)=DB(i,1)  ;
    px(2)=DB(i,2); py(2)=DB(i,2)  ;
    plot(px,py,'k');   %画原始点及线路
end
 
pos=32;  %罪犯逃跑点
plot(pos,pos,'bo');    
%求到包围点最近的平台
kp=length(DC);      % 获得平台个数    length(DC)=80
U=zeros(kp,1);        % 判断该平台是否需要参与围捕活动  U(s)==1 参加 U(s)==0 不参加
Dp=zeros(length(Ct),1); % 对应Ct记录包围点获得指定的平台    length(Ct)=24个需要封堵的路口
TT=zeros(length(Ct),1);  % 对应Ct记录存储平台到包围点的最短时间
% 遍历所有交巡警平台到每个路口的时间，取最小值
for i=1:length(Ct);    % 遍历封堵路口
    tmin=1000;      % 设置t初值，优化最长时间
    for j=1:length(DC);     % 遍历交巡警服务平台
        tt=A(Ct(i),DC(j,1))/v;      % 从最短路矩阵计算时间
        % 如果第j个平台到第i路口的时间tt小于tmin，并且第j个平台目前没有任务，tmin 迭代更新，记录平台号
       if tt<tmin &&U(j)==0     
           s=j; tmin=tt;      
       end
    end
    % 对应记录第i个需要封堵路口的最短时间，Dp对应记录处理路口的交巡警平台号
    TT(i)=tmin; U(s)=1; Dp(i)=DC(s,1); 
end
 
for i=1:length(Ct);     % 再遍历需围堵路口
    T1(i)=A(Dp(i),Ct(i))/v;     % 计算第i个围堵平台到第i个需要维度的路口路口所需时间
    T2(i)=A(pos,Ct(i))/v1;       % 计算罪犯从A区起点32到第i个路口所需时间
end
 
  for i=1:length(Ct)
      plot(Ct(i),Ct(i),'ro');   %画包围圈点
  end
  
% for i=1:length(Dp)
%     plot(x(Dp(i)),y(Dp(i)),'bs');   %画平台点
% end
xlabel('x坐标');
ylabel('y坐标');
fprintf('序号 新包围点,平台，封锁时间,罪犯到达时间\n');
 
for i=1:length(Ct)
    fprintf('%4d %4d  %4d   %6.4f    %6.4f\n',i,Ct(i),Dp(i),TT(i)+3,T2(i));
end 
Info=zeros(length(Ct),4);
 
for i=1:length(Ct)  
   Info(i,1)=i;Info(i,2)=Ct(i); Info(i,3)=Dp(i); Info(i,4)=TT(i)+3; Info(i,5)=T2(i);
end     
xlswrite('第五问答案2011B.xlsx',Info); %输出到XLSX文件中

%% 灵敏度分析
x=xlsread('第五问答案2011B（总）.xlsx','Sheet1','G4:G2');
y=xlsread('第五问答案2011B（总）.xlsx','Sheet1','H4:H2');
p=polyfit(x,y,3);
y2=polyval(p,x);
p=vpa(poly2sym(p),2);

x0=40:1:60;
y0=- 5.5e-5.*x0.^3 + 6.9e-4.*x0.^2 + 0.31.*x0;
plot(x0,y0,'k','linewidth',1.5)
xlabel('罪犯速度(km/h)','FontSize',16);
ylabel('给定时间(min)','FontSize',16);
set(gca,'FontSize',16,'linewidth',1);