全国大学生数学建模竞赛2015B题“互联网+”时代的出租车资源配置MATLAB程序

目录

一、第1问

二、第2问

三、第3问

---

一、第1问

注：前五段代码为数据预处理代码，处理数据为当前广为流传的经过处理的数据。由于数据量太大以及本人水平有限，前五段代码运行时间较长，每段代码根据电脑性能不一需运行1~2小时。经过处理的数据为：地点南京，时间2016 年 8 月 6 日 - 12 日每日 0 点 - 12 点内的五种统计信息。 信息解释为 demand（打车需求量） distribute（出租车分布） money（车费） response（被抢单时间） satisfy（打车难易度） 五种信息。

2015B题经过处理的数据提取链接

链接：https://pan.baidu.com/s/1NSq8F8nscGLvL5achbp9EA 
提取码：6666

%% 定义变量（4292个demand）
data=xlsread('demand_2016.04.06_320100_.csv','demand_2016.04.06_320100_','B2:E4293');
hour=data(:,1);
lon=data(:,2);
lat=data(:,3);
value=data(:,4);
geobubble(lat,lon);
count=zeros(25,1);
count00=zeros(25,1);
%% 计数24小时
for t=0:23
    value000=0;
    for i0=1:4292
    if hour(i0)==t
       count(t+2,1)=count(t+2,1)+1;
    end
    end
end
for i=1:25
    count00(i)=sum(count(1:i,1));
end
count000=count00+1;
%% 分区域化，将小矩阵写入大矩阵（24小时）
for t=0:23
lon0=zeros(1000,20);%存储经度
lat0=zeros (1000,20);%存储纬度
m=zeros(1000,20);%存储标号
ll=cell(20,20);
mm=cell(20,20);
ll(:)={0};
mm(:)={0};
    for nn0=1:20 %lat遍历
    for n0=1:20 %lon遍历
        n=0;
    for i=count000(t+1):count00(t+2)
       if lon(i)<=118.25+0.05*n0 && lon(i)>=118.25+0.05*(n0-1) && lat(i)<=32.5-0.05*(nn0-1) && lat(i)>=32.5-0.05*nn0 
       n=n+1;
       lon0(n,n0)=lon(i);
       lat0(n,n0)=lat(i);
       m(n,n0)=i;
       
       ll{nn0,n0}=[lon0(:,n0),lat0(:,n0)];
       mm{nn0,n0}=m(:,n0);
       end 
    end
    end
    end
%% 分区域计算value （24小时）
value00=zeros(20,20);
for i=1:20
    for j=1:20
        count0=0;
        value0=zeros(1000,1);
            for jj=1:1000
                if mm{j,i}==0
                    continue
                end
                if mm{j,i}(jj,:)==0
                    continue
                end
                count0=count0+1;
                value0(count0,1)=value(mm{j,i}(jj,:),1);
            end
        value00(j,i)=sum(value0)/count0;
    end
end
ind=find(isnan(value00));%找出所有为NaN的位置标号
value00(ind)=0;
%% 依次填入表格内
ii=0+400*t;
for j=1:20
    for i=1:20
        ii=ii+1;
        value000=value00(j,i);
        xlswrite('didi2015B.xlsx',t,'Sheet1',['B',num2str(ii)]);
        xlswrite('didi2015B.xlsx',j,'Sheet1',['C',num2str(ii)]);
        xlswrite('didi2015B.xlsx',i,'Sheet1',['D',num2str(ii)]);
        xlswrite('didi2015B.xlsx',value000,'Sheet1',['E',num2str(ii)]);
    end
end
end

%% 定义变量（12000个）
data=xlsread('distribute_2016.04.06_320100_.csv','distribute_2016.04.06_320100_','B2:E12001');
hour=data(:,1);
lon=data(:,2);
lat=data(:,3);
value=data(:,4);
geobubble(lat,lon);
count=zeros(25,1);
count00=zeros(25,1);
%% 计数24小时
for t=0:23
    value000=0;
    for i0=1:12000
    if hour(i0)==t
       count(t+2,1)=count(t+2,1)+1;
    end
    end
end
for i=1:25
    count00(i)=sum(count(1:i,1));
end
count000=count00+1;
%% 分区域化，将小矩阵写入大矩阵（24小时）
for t=0:23
lon0=zeros(5000,20);%存储经度
lat0=zeros (5000,20);%存储纬度
m=zeros(5000,20);%存储标号
ll=cell(20,20);
mm=cell(20,20);
ll(:)={0};
mm(:)={0};
    for nn0=1:20 %lat遍历
    for n0=1:20 %lon遍历
        n=0;
    for i=count000(t+1):count00(t+2)
       if lon(i)<=118.25+0.05*n0 && lon(i)>=118.25+0.05*(n0-1) && lat(i)<=32.5-0.05*(nn0-1) && lat(i)>=32.5-0.05*nn0 
       n=n+1;
       lon0(n,n0)=lon(i);
       lat0(n,n0)=lat(i);
       m(n,n0)=i;
       
       ll{nn0,n0}=[lon0(:,n0),lat0(:,n0)];
       mm{nn0,n0}=m(:,n0);
       end 
    end
    end
    end
%% 分区域计算value （24小时）
value00=zeros(20,20);
for i=1:20
    for j=1:20
        count0=0;
        value0=zeros(5000,1);
            for jj=1:5000
                if mm{j,i}==0
                    continue
                end
                if mm{j,i}(jj,:)==0
                    continue
                end
                count0=count0+1;
                value0(count0,1)=value(mm{j,i}(jj,:),1);
            end
        value00(j,i)=sum(value0)/count0;
    end
end
ind=find(isnan(value00));%找出所有为NaN的位置标号
value00(ind)=0;
%% 依次填入表格内
ii=0+400*t;
for j=1:20
    for i=1:20
        ii=ii+1;
        value000=value00(j,i);
        xlswrite('didi2015B.xlsx',t,'Sheet1',['B',num2str(ii)]);
        xlswrite('didi2015B.xlsx',j,'Sheet1',['C',num2str(ii)]);
        xlswrite('didi2015B.xlsx',i,'Sheet1',['D',num2str(ii)]);
        xlswrite('didi2015B.xlsx',value000,'Sheet1',['E',num2str(ii)]);
    end
end
end

%% 定义变量（2825个）
data=xlsread('money_2016.04.06_320100_.csv','money_2016.04.06_320100_','B2:E2826');
hour=data(:,1);
lon=data(:,2);
lat=data(:,3);
value=data(:,4);
geobubble(lat,lon);
count=zeros(25,1);
count00=zeros(25,1);
%% 计数24小时
for t=0:23
    value000=0;
    for i0=1:2825
    if hour(i0)==t
       count(t+2,1)=count(t+2,1)+1;
    end
    end
end
for i=1:25
    count00(i)=sum(count(1:i,1));
end
count000=count00+1;
%% 分区域化，将小矩阵写入大矩阵（24小时）
for t=0:23
lon0=zeros(1000,20);%存储经度
lat0=zeros (1000,20);%存储纬度
m=zeros(1000,20);%存储标号
ll=cell(20,20);
mm=cell(20,20);
ll(:)={0};
mm(:)={0};
    for nn0=1:20 %lat遍历
    for n0=1:20 %lon遍历
        n=0;
    for i=count000(t+1):count00(t+2)
       if lon(i)<=118.25+0.05*n0 && lon(i)>=118.25+0.05*(n0-1) && lat(i)<=32.5-0.05*(nn0-1) && lat(i)>=32.5-0.05*nn0 
       n=n+1;
       lon0(n,n0)=lon(i);
       lat0(n,n0)=lat(i);
       m(n,n0)=i;
       
       ll{nn0,n0}=[lon0(:,n0),lat0(:,n0)];
       mm{nn0,n0}=m(:,n0);
       end 
    end
    end
    end
%% 分区域计算value （24小时）
value00=zeros(20,20);
for i=1:20
    for j=1:20
        count0=0;
        value0=zeros(1000,1);
            for jj=1:1000
                if mm{j,i}==0
                    continue
                end
                if mm{j,i}(jj,:)==0
                    continue
                end
                count0=count0+1;
                value0(count0,1)=value(mm{j,i}(jj,:),1);
            end
        value00(j,i)=sum(value0)/count0;
    end
end
ind=find(isnan(value00));%找出所有为NaN的位置标号
value00(ind)=0;
%% 依次填入表格内
ii=0+400*t;
for j=1:20
    for i=1:20
        ii=ii+1;
        value000=value00(j,i);
        xlswrite('didi2015B.xlsx',t,'Sheet1',['B',num2str(ii)]);
        xlswrite('didi2015B.xlsx',j,'Sheet1',['C',num2str(ii)]);
        xlswrite('didi2015B.xlsx',i,'Sheet1',['D',num2str(ii)]);
        xlswrite('didi2015B.xlsx',value000,'Sheet1',['E',num2str(ii)]);
    end
end
end

%% 定义变量（3832个）
data=xlsread('response_2016.04.06_320100_.csv','response_2016.04.06_320100_','B2:E3833');
hour=data(:,1);
lon=data(:,2);
lat=data(:,3);
value=data(:,4);
geobubble(lat,lon);
count=zeros(25,1);
count00=zeros(25,1);
%% 计数24小时
for t=0:23
    value000=0;
    for i0=1:3832
    if hour(i0)==t
       count(t+2,1)=count(t+2,1)+1;
    end
    end
end
for i=1:25
    count00(i)=sum(count(1:i,1));
end
count000=count00+1;
%% 分区域化，将小矩阵写入大矩阵（24小时）
for t=0:23
lon0=zeros(1000,20);%存储经度
lat0=zeros (1000,20);%存储纬度
m=zeros(1000,20);%存储标号
ll=cell(20,20);
mm=cell(20,20);
ll(:)={0};
mm(:)={0};
    for nn0=1:20 %lat遍历
    for n0=1:20 %lon遍历
        n=0;
    for i=count000(t+1):count00(t+2)
       if lon(i)<=118.25+0.05*n0 && lon(i)>=118.25+0.05*(n0-1) && lat(i)<=32.5-0.05*(nn0-1) && lat(i)>=32.5-0.05*nn0 
       n=n+1;
       lon0(n,n0)=lon(i);
       lat0(n,n0)=lat(i);
       m(n,n0)=i;
       
       ll{nn0,n0}=[lon0(:,n0),lat0(:,n0)];
       mm{nn0,n0}=m(:,n0);
       end 
    end
    end
    end
%% 分区域计算value （24小时）
value00=zeros(20,20);
for i=1:20
    for j=1:20
        count0=0;
        value0=zeros(1000,1);
            for jj=1:1000
                if mm{j,i}==0
                    continue
                end
                if mm{j,i}(jj,:)==0
                    continue
                end
                count0=count0+1;
                value0(count0,1)=value(mm{j,i}(jj,:),1);
            end
        value00(j,i)=sum(value0)/count0;
    end
end
ind=find(isnan(value00));%找出所有为NaN的位置标号
value00(ind)=0;
%% 依次填入表格内
ii=0+400*t;
for j=1:20
    for i=1:20
        ii=ii+1;
        value000=value00(j,i);
        xlswrite('didi2015B.xlsx',t,'Sheet1',['B',num2str(ii)]);
        xlswrite('didi2015B.xlsx',j,'Sheet1',['C',num2str(ii)]);
        xlswrite('didi2015B.xlsx',i,'Sheet1',['D',num2str(ii)]);
        xlswrite('didi2015B.xlsx',value000,'Sheet1',['E',num2str(ii)]);
    end
end
end

%% 定义变量（12000个）
data=xlsread('satisfy_2016.04.06_320100_.csv','satisfy_2016.04.06_320100_','B2:E12001');
hour=data(:,1);
lon=data(:,2);
lat=data(:,3);
value=data(:,4);
geobubble(lat,lon);
count=zeros(25,1);
count00=zeros(25,1);
%% 计数24小时
for t=0:23
    value000=0;
    for i0=1:12000
    if hour(i0)==t
       count(t+2,1)=count(t+2,1)+1;
    end
    end
end
for i=1:25
    count00(i)=sum(count(1:i,1));
end
count000=count00+1;
%% 分区域化，将小矩阵写入大矩阵（24小时）
for t=0:23
lon0=zeros(5000,20);%存储经度
lat0=zeros (5000,20);%存储纬度
m=zeros(5000,20);%存储标号
ll=cell(20,20);
mm=cell(20,20);
ll(:)={0};
mm(:)={0};
    for nn0=1:20 %lat遍历
    for n0=1:20 %lon遍历
        n=0;
    for i=count000(t+1):count00(t+2)
       if lon(i)<=118.25+0.05*n0 && lon(i)>=118.25+0.05*(n0-1) && lat(i)<=32.5-0.05*(nn0-1) && lat(i)>=32.5-0.05*nn0 
       n=n+1;
       lon0(n,n0)=lon(i);
       lat0(n,n0)=lat(i);
       m(n,n0)=i;
       
       ll{nn0,n0}=[lon0(:,n0),lat0(:,n0)];
       mm{nn0,n0}=m(:,n0);
       end 
    end
    end
    end
%% 分区域计算value （24小时）
value00=zeros(20,20);
for i=1:20
    for j=1:20
        count0=0;
        value0=zeros(5000,1);
            for jj=1:5000
                if mm{j,i}==0
                    continue
                end
                if mm{j,i}(jj,:)==0
                    continue
                end
                count0=count0+1;
                value0(count0,1)=value(mm{j,i}(jj,:),1);
            end
        value00(j,i)=sum(value0)/count0;
    end
end
ind=find(isnan(value00));%找出所有为NaN的位置标号
value00(ind)=0;
%% 依次填入表格内
ii=0+400*t;
for j=1:20
    for i=1:20
        ii=ii+1;
        value000=value00(j,i);
        xlswrite('didi2015B.xlsx',t,'Sheet1',['B',num2str(ii)]);
        xlswrite('didi2015B.xlsx',j,'Sheet1',['C',num2str(ii)]);
        xlswrite('didi2015B.xlsx',i,'Sheet1',['D',num2str(ii)]);
        xlswrite('didi2015B.xlsx',value000,'Sheet1',['E',num2str(ii)]);
    end
end
end

%% 归一化处理
function y=guiyi(x,type,ymin,ymax)
%实现正向或负向指标归一化，返回归一化后的数据矩阵
%x为原始数据矩阵, 一行代表一个样本, 每列对应一个指标
%type设定正向指标1,负向指标2
%ymin,ymax为归一化的区间端点
[n,m]=size(x);
y=zeros(n,m);
xmin=min(x);
xmax=max(x);
switch type
    case 1
        for j=1:m
            y(:,j)=(ymax-ymin)*(x(:,j)-xmin(j))/(xmax(j)-xmin(j))+ymin;
        end
    case 2
        for j=1:m
            y(:,j)=(ymax-ymin)*(xmax(j)-x(:,j))/(xmax(j)-xmin(j))+ymin;
        end
end

%% 熵权法（内嵌函数）
function [s,w]=shang(x,ind)
%实现用熵值法求各指标(列）的权重及各数据行的得分
%x为原始数据矩阵, 一行代表一个样本, 每列对应一个指标
%ind指示向量，指示各列正向指标还是负向指标，1表示正向指标，2表示负向指标
%s返回各行（样本）得分，w返回各列权重
[n,m]=size(x); % n个样本, m个指标
%%数据的归一化处理
for i=1:m
    if ind(i)==1 %正向指标归一化
        X(:,i)=guiyi(x(:,i),1,0.002,0.996);    %若归一化到[0,1], 0会出问题
    else %负向指标归一化
        X(:,i)=guiyi(x(:,i),2,0.002,0.996);
    end
end
%%计算第j个指标下，第i个样本占该指标的比重p(i,j)
for i=1:n
    for j=1:m
        p(i,j)=X(i,j)/sum(X(:,j));
    end
end
%%计算第j个指标的熵值e(j)
k=1/log(n);
for j=1:m
    e(j)=-k*sum(p(:,j).*log(p(:,j)));
end
d=ones(1,m)-e; %计算信息熵冗余度
w=d./sum(d); %求权值w
s=100*w*X'; %求综合得分

%% 熵权法
data=xlsread('didi2015B总.xlsx','Sheet1','E2:I9601');
Ind=[2 1 1 2 2]; %指定各指标的正向or负向
[S,W]=shang(data,Ind);
S=S';
xlswrite('didi2015B总.xlsx',S,'Sheet1','J2');

%% 绘图（hour=8）
x8=xlsread('didi2015B总.xlsx','Sheet1','C3202:C3601');
y8=xlsread('didi2015B总.xlsx','Sheet1','D3202:D3601');
shang8=xlsread('didi2015B总.xlsx','Sheet1','J3202:J3601');
huitu8=zeros(20);
for i=1:400
    huitu8(x8(i),y8(i))=shang8(i);
end
x=(1:20)';
y=(1:20)';
[x1,y1]=meshgrid(1:0.01:20,1:0.01:20);
z1=interp2(x,y,huitu8,x1,y1,'spline');
%subplot(1,2,1);%可分区域画图
surf(x1,y1,z1)
shading interp
colormap gray;
axis([0 20 0 20 -20 60])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
zlabel('F(x,y,8)','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);

%% 绘图（hour=23）
x23=xlsread('didi2015B总.xlsx','Sheet1','C9202:C9601');
y23=xlsread('didi2015B总.xlsx','Sheet1','D9202:D9601');
shang23=xlsread('didi2015B总.xlsx','Sheet1','J9202:J9601');
huitu23=zeros(20);
for i=1:400
    huitu23(x23(i),y23(i))=shang23(i);
end
x=(1:20)';
y=(1:20)';
[x1,y1]=meshgrid(1:0.01:20,1:0.01:20);
z1=interp2(x,y,huitu23,x1,y1,'spline');
%subplot(1,2,2);%可分区域画图
figure %若分区域画图，删去本行
surf(x1,y1,z1)
shading interp
colormap gray;
axis([0 20 0 20 -20 60])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
zlabel('F(x,y,23)','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);

%% 熵权法1.1(dis+money)
data1=xlsread('didi2015B总.xlsx','Sheet1','F2:G9601');
Ind1=[1 1]; %指定各指标的正向or负向
[S1,W1]=shang(data1,Ind1);
S1=S1';
xlswrite('didi2015B总.xlsx',S1,'Sheet2','A2');
%% 熵权法1.2(demand+res+satisfy)
data21=xlsread('didi2015B总.xlsx','Sheet1','E2:E9601');
data22=xlsread('didi2015B总.xlsx','Sheet1','H2:I9601');
data2=[data21 data22];
Ind2=[1 1 1]; %指定各指标的正向or负向
[S2,W2]=shang(data2,Ind2);
S2=S2';
xlswrite('didi2015B总.xlsx',S2,'Sheet2','D2');
%% 回归分析1.1
Y1=S1;X1=[ones(9600,1) data1];
[b1,bint1,r1,rint1,stats1]=regress(Y1,X1);
%% 回归分析1.2
Y2=S2;X2=[ones(9600,1) data2];
[b2,bint2,r2,rint2,stats2]=regress(Y2,X2);

二、第2问

%% 数据预处理（标准化并计算补贴）
x=xlsread('didi2015B总.xlsx','Sheet1','E2:F9601');
ymin=0.002;
ymax=0.996;
k=6;
%实现正向或负向指标归一化，返回归一化后的数据矩阵
%x为原始数据矩阵, 一行代表一个样本, 每列对应一个指标
%type设定正向指标1,负向指标2
%ymin,ymax为归一化的区间端点
[n,m]=size(x);
y=zeros(n,m);
xmin=min(x);
xmax=max(x);
for j=1:m
    y(:,j)=(ymax-ymin)*(x(:,j)-xmin(j))/(xmax(j)-xmin(j))+ymin;
end
%计算补贴
xx=y(:,1)./y(:,2);
ind=find(xx==1);%找出所有为0的位置标号
xx(ind)=0;
xx=xx*k;
xlswrite('didi2015B总.xlsx',xx,'Sheet1','K2:K9601');

%% 熵权法
data1=xlsread('didi2015B总.xlsx','Sheet1','E2:I9601');
data2=xlsread('didi2015B总.xlsx','Sheet1','K2:K9601');
data=[data1 data2];
Ind=[2 1 1 2 2 1]; %指定各指标的正向or负向
[S,W]=shang(data,Ind);
S=S';
xlswrite('didi2015B总.xlsx',S,'Sheet1','L2');

%% 绘图（hour=8）
x8=xlsread('didi2015B总.xlsx','Sheet1','C3202:C3601');
y8=xlsread('didi2015B总.xlsx','Sheet1','D3202:D3601');
shang8=xlsread('didi2015B总.xlsx','Sheet1','L3202:L3601');
huitu8=zeros(20);
for i=1:400
    huitu8(x8(i),y8(i))=shang8(i);
end
x=(1:20)';
y=(1:20)';
[x1,y1]=meshgrid(1:0.01:20,1:0.01:20);
z1=interp2(x,y,huitu8,x1,y1,'spline');
%subplot(1,2,1);%可分区域画图
surf(x1,y1,z1)
shading interp
colormap gray;
axis([0 20 0 20 -20 60])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
zlabel('F(x,y,8)','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);

%% 绘图（hour=23）
x23=xlsread('didi2015B总.xlsx','Sheet1','C9202:C9601');
y23=xlsread('didi2015B总.xlsx','Sheet1','D9202:D9601');
shang23=xlsread('didi2015B总.xlsx','Sheet1','L9202:L9601');
huitu23=zeros(20);
for i=1:400
    huitu23(x23(i),y23(i))=shang23(i);
end
x=(1:20)';
y=(1:20)';
[x1,y1]=meshgrid(1:0.01:20,1:0.01:20);
z1=interp2(x,y,huitu23,x1,y1,'spline');
%subplot(1,2,2);%可分区域画图
figure %若分区域画图，删去本行
surf(x1,y1,z1)
shading interp
colormap gray;
axis([0 20 0 20 -20 60])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
zlabel('F(x,y,23)','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);

%% 熵权法2.1(dis+money+butie2)
data0=xlsread('didi2015B总.xlsx','Sheet1','K2:K9601');
data3=[data1 data0];
Ind3=[1 1 1]; %指定各指标的正向or负向
[S3,W3]=shang(data3,Ind3);
S3=S3';
xlswrite('didi2015B总.xlsx',S3,'Sheet2','H2');
%% 熵权法2.2(demand+res+satisfy)
data4=data2;
Ind4=[1 1 1]; %指定各指标的正向or负向
[S4,W4]=shang(data4,Ind4);
S4=S4';
xlswrite('didi2015B总.xlsx',S4,'Sheet2','L2');
%% 回归分析2.1
Y3=S3;X3=[ones(9600,1) data3];
[b3,bint3,r3,rint3,stats3]=regress(Y3,X3);
%% 回归分析2.2
Y4=S4;X4=[ones(9600,1) data4];
[b4,bint4,r4,rint4,stats4]=regress(Y4,X4);

三、第3问

%% 数据预处理（代入补贴计算）
x=xlsread('didi2015B总.xlsx','Sheet1','E2:F9601');
butie=xlsread('didi2015B总.xlsx','Sheet1','K2:K9601');
butie=butie/6;
ymin=0.002;
ymax=0.996;
nn=1170;%zhibiao1小于1的量=zhibiao3小于1的量
xx=zeros(9600,1);
for i=1:9600
%实现正向或负向指标归一化，返回归一化后的数据矩阵
%x为原始数据矩阵, 一行代表一个样本, 每列对应一个指标
%type设定正向指标1,负向指标2
%ymin,ymax为归一化的区间端点
[n,m]=size(x);
y=zeros(n,m);
xmin=min(x);
xmax=max(x);
for j=1:m
    y(:,j)=(ymax-ymin)*(x(:,j)-xmin(j))/(xmax(j)-xmin(j))+ymin;
end
%计算补贴
xx(i,1)=y(i,1)./y(i,2);
if xx(i,1)==1%找出所有为0的位置标号
xx(i,1)=0;
end
if butie(i)<1
xx(i)=xx(i)*10;
end
if butie(i)>=1 && butie(i)<50
xx(i)=xx(i)*6;
end
if butie(i)>=50 && butie(i)<100
xx(i)=xx(i)*4;
end
if butie(i)>=20
xx(i)=xx(i)*2;
end
end
xlswrite('didi2015B总.xlsx',xx,'Sheet1','M2');

%% 熵权法
data1=xlsread('didi2015B总.xlsx','Sheet1','E2:I9601');
data2=xx;
data=[data1 data2];
Ind=[1 1 1 1 1 1]; %指定各指标的正向or负向
[S,W]=shang(data,Ind);
S=S';
xlswrite('didi2015B总.xlsx',S,'Sheet1','N2');
zhibiao3=S;
w=zeros(9600,1);
for ii=1:9600
if zhibiao3(ii)<1
    continue
end
if 20-zhibiao3(ii)>0 %不满足大于20
    w(ii)=20-zhibiao3(ii);
end
if 20-zhibiao3(ii)<0 %满足大于20
    w(ii)=0;
end
end
min1=sum(w)/(20*nn);%最终指标（供求关系匹配度）

%% 绘图（hour=8）
x8=xlsread('didi2015B总.xlsx','Sheet1','C3202:C3601');
y8=xlsread('didi2015B总.xlsx','Sheet1','D3202:D3601');
shang8=xlsread('didi2015B总.xlsx','Sheet1','N3202:N3601');
huitu8=zeros(20);
for i=1:400
    huitu8(x8(i),y8(i))=shang8(i);
end
x=(1:20)';
y=(1:20)';
[x1,y1]=meshgrid(1:0.01:20,1:0.01:20);
z1=interp2(x,y,huitu8,x1,y1,'spline');
%subplot(1,2,1);%可分区域画图
surf(x1,y1,z1)
shading interp
colormap gray;
axis([0 20 0 20 -10 40])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
zlabel('F(x,y,8)','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);

%% 绘图（hour=23）
x23=xlsread('didi2015B总.xlsx','Sheet1','C9202:C9601');
y23=xlsread('didi2015B总.xlsx','Sheet1','D9202:D9601');
shang23=xlsread('didi2015B总.xlsx','Sheet1','N9202:N9601');
huitu23=zeros(20);
for i=1:400
    huitu23(x23(i),y23(i))=shang23(i);
end
x=(1:20)';
y=(1:20)';
[x1,y1]=meshgrid(1:0.01:20,1:0.01:20);
z1=interp2(x,y,huitu23,x1,y1,'spline');
%subplot(1,2,2);%可分区域画图
figure %若分区域画图，删去本行
surf(x1,y1,z1)
shading interp
colormap gray;
axis([0 20 0 20 -10 40])
xlabel('x','FontSize',28);
ylabel('y','FontSize',28);
zlabel('F(x,y,23)','FontSize',28);
set(gca,'FontSize',28,'linewidth',1);