目录
1.2 计算模型中等待时间t2(灵敏度分析)MATLAB程序
2.2 代入t2到模型进行计算(灵敏度分析)MATLAB程序
前言
使用的数据为:上海虹桥机场T2航站楼航班数据2021.8.26(上海机场集团官网可以找到)
链接:https://pan.baidu.com/s/18VhpM8wssX2rDgSqEMwqwg
提取码:6666
一、第1问
1.1 计算模型中等待时间t2MATLAB程序
%% 计算模型中等待时间t2 %% 每1小时划分到达航班 t1=xlsread('上海虹桥机场T2航站楼航班数据2021.8.26.xlsx','到达','A2:A211'); t1=datevec(t1);%处理时间函数 zz1=zeros(24,1);%每1小时飞机航班 f1=zeros(24,1);%每1小时每架飞机的人数 R1=zeros(24,1);%每1小时乘客乘坐出租车的比例 z1=zeros(24,1);%每1小时需要车的个数 for i=1:24 %航班划分 tt1=find(t1(:,4)==i); size1=size(tt1,1); zz1(i,1)=size1; end for i=1:24 if i>=8 && i<=18 %乘坐飞机人数划分 f1(i)=180; else f1(i)=120; end if i>=8 && i<=21 %乘坐出租车比例划分 R1(i)=0.15; else R1(i)=0.3; end end for i=1:24 %最终乘坐出租车的人数/1.8=多少辆出租车 z1(i)=zz1(i)*f1(i)*R1(i)/1.8; end z1=round(z1);%每1小时需要车的个数 %% 每1小时划分出发航班 t2=xlsread('上海虹桥机场T2航站楼航班数据2021.8.26.xlsx','出发','A2:A215'); t2=datevec(t2);%处理时间函数 zz2=zeros(24,1);%每1小时飞机航班 f2=zeros(24,1);%每1小时每架飞机的人数 R2=zeros(24,1);%每1小时乘客乘坐出租车的比例 z2=zeros(24,1);%每1小时到达车的个数 for i=1:24 %航班划分 tt2=find(t2(:,4)==i); size2=size(tt2,1); zz2(i,1)=size2; end for i=1:24 if i>=8 && i<=18 %乘坐飞机人数划分 f2(i)=180; else f2(i)=120; end if i>=8+2 && i<=21+2 %乘坐出租车比例划分(需提前2小时出发,2小时后到达机场) R2(i)=0.15; else R2(i)=0.3; end end for i=1:24 %最终乘坐出租车的人数/1.8=多少辆出租车 z2(i)=zz2(i)*f2(i)*R2(i)/1.8; end z2=round(z2*0.5);%每1小时到达车的个数 plot(zz1,'linewidth',2); xlabel('t (h)','FontSize',26); ylabel('航班数 (架)','FontSize',26); set(gca,'FontSize',26,'linewidth',1); set(gca,'XTick',0:2:24); figure plot(zz2,'linewidth',2); xlabel('t (h)','FontSize',26); ylabel('航班数 (架)','FontSize',26); set(gca,'FontSize',26,'linewidth',1); set(gca,'XTick',0:2:24); %% 代入模型计算 lamda=zeros(24,1); C=zeros(24,1); u=zeros(24,1); F=zeros(24,1); p=zeros(24,1); w=zeros(24,1); z=zeros(24,1); p0=zeros(24,1); c=3;%上车点个数 for t=1:24 z(t)=z1(t); lamda(t)=z2(t); t1=2; t2=0.5; t3=0.5; u0=t1+t2+t3;%一辆车服务一位乘客 n0=round(60/u0); if z(t)>n0 C(t)=0; else C(t)=n0/z(t); end u(t)=u0+C(t); F(t)=1-exp(-u(t)); p(t)=lamda(t)/(c*u(t)); kk=zeros(3,1); for k=0:c-1 kk(k+1,1)=(1/factorial(k))*(lamda(t)/u(t))^k; end kkk=sum(kk); p0(t)=(kkk+(1/factorial(c))*(1/(1-p(t)))*(lamda(t)/u(t))^c)^(-1); w(t)=abs(((c*p(t))^c*p(t))/(factorial(c)*(1-p(t))^2*lamda(t))*p0(t)); end ind=find(isnan(w));%找出A中所有为NaN的位置标号 w(ind)=0; w=60*w; save('w.mat','w');%保存等待时间t2
1.2 计算模型中等待时间t2(灵敏度分析)MATLAB程序
%% 计算模型中等待时间t2(灵敏度分析) lmd1=1.0; lmd2=1.0; lmd3=1.0; lmd4=1.0; imd5=1.0; %% 每1小时划分到达航班 t1=xlsread('上海虹桥机场T2航站楼航班数据2021.8.26.xlsx','到达','A2:A211'); t1=datevec(t1);%处理时间函数 zz1=zeros(24,1);%每1小时飞机航班 f1=zeros(24,1);%每1小时每架飞机的人数 R1=zeros(24,1);%每1小时乘客乘坐出租车的比例 z1=zeros(24,1);%每1小时需要车的个数 for i=1:24 %航班划分 tt1=find(t1(:,4)==i); size1=size(tt1,1); zz1(i,1)=size1; end for i=1:24 if i>=8 && i<=18 %乘坐飞机人数划分 f1(i)=180*lmd1; else f1(i)=120*lmd1; end if i>=8 && i<=21 %乘坐出租车比例划分 R1(i)=0.15; else R1(i)=0.3; end end for i=1:24 %最终乘坐出租车的人数/1.8=多少辆出租车 z1(i)=zz1(i)*f1(i)*R1(i)/1.8; end z1=round(z1);%每1小时需要车的个数 %% 每1小时划分出发航班 t2=xlsread('上海虹桥机场T2航站楼航班数据2021.8.26.xlsx','出发','A2:A215'); t2=datevec(t2);%处理时间函数 zz2=zeros(24,1);%每1小时飞机航班 f2=zeros(24,1);%每1小时每架飞机的人数 R2=zeros(24,1);%每1小时乘客乘坐出租车的比例 z2=zeros(24,1);%每1小时到达车的个数 for i=1:24 %航班划分 tt2=find(t2(:,4)==i); size2=size(tt2,1); zz2(i,1)=size2; end for i=1:24 if i>=8 && i<=18 %乘坐飞机人数划分 f2(i)=180*lmd1; else f2(i)=120*lmd1; end if i>=8+2 && i<=21+2 %乘坐出租车比例划分(需提前2小时出发,2小时后到达机场) R2(i)=0.15; else R2(i)=0.3; end end for i=1:24 %最终乘坐出租车的人数/1.8=多少辆出租车 z2(i)=zz2(i)*f2(i)*R2(i)/1.8; end z2=round(z2*0.5);%每1小时到达车的个数 %% 代入模型计算 lamda=zeros(24,1); C=zeros(24,1); u=zeros(24,1); F=zeros(24,1); p=zeros(24,1); w=zeros(24,1); z=zeros(24,1); p0=zeros(24,1); c=3;%上车点个数 for t=1:24 z(t)=z1(t); lamda(t)=z2(t); t1=2; t2=0.5; t3=0.5; u0=t1+t2+t3;%一辆车服务一位乘客 n0=round(60/u0); if z(t)>n0 C(t)=0; else C(t)=n0/z(t); end u(t)=u0+C(t); F(t)=1-exp(-u(t)); p(t)=lamda(t)/(c*u(t)); kk=zeros(3,1); for k=0:c-1 kk(k+1,1)=(1/factorial(k))*(lamda(t)/u(t))^k; end kkk=sum(kk); p0(t)=(kkk+(1/factorial(c))*(1/(1-p(t)))*(lamda(t)/u(t))^c)^(-1); w(t)=abs(((c*p(t))^c*p(t))/(factorial(c)*(1-p(t))^2*lamda(t))*p0(t)); end ind=find(isnan(w));%找出A中所有为NaN的位置标号 w(ind)=0; w=60*w; save('wlmd.mat','w');%保存等待时间t2
二、第2问
2.1 代入t2到模型进行计算MATLAB程序
%第2问MATLAB程序 %% 代入t2到模型进行计算 load('w');%载入等待时间t2 t1=5;%送客区至排队区时间 t2=w;%排队区等待时间 t3=30;%向市中心运客时间 t4=t3;%向市中心运客时间 t5=t1+t2+t3-t4;%在市中心运客时间 s=0.5;%每公里0.5元的燃油费 x1=8;%虹桥机场到市区的距离 v=40/60;%出租车的均速40KM/h m=4;%出租车平均每公里收费 ph1=zeros(24,1);%里程利用率 for tt=1:4 ph1(tt)=0.30; end for tt=5:6 ph1(tt)=0.45; end for tt=7:10 ph1(tt)=0.55; end for tt=11:12 ph1(tt)=0.75; end for tt=13:14 ph1(tt)=0.55; end for tt=15:17 ph1(tt)=0.45; end for tt=18:19 ph1(tt)=0.55; end for tt=20:21 ph1(tt)=0.35; end for tt=22:24 ph1(tt)=0.30; end tA=zeros(24,1); DQA=zeros(24,1); tB=zeros(24,1); PB=zeros(24,1); CB=zeros(24,1); QB=zeros(24,1); DQB=zeros(24,1); for t=1:24 %A方案 tA(t)=t1+t2(t)+t3;%A方案的总时间 PA=x1*m;%A方案载客至市区赚的钱 CA=s*x1;%A方案油耗 QA=PA-CA;%A方案收益 DQA(t)=QA/tA(t);%A方案单位时间收益 %B方案 tB(t)=t4+t5(t);%B方案的总时间 PB(t)=ph1(t)*(t5(t)*v)*m;%B方案载客至市区赚的钱 CB(t)=s*(x1+t5(t)*v);%B方案油耗 QB(t)=PB(t)-CB(t);%B方案收益 DQB(t)=QB(t)/tB(t);%B方案单位时间收益 end DQ=DQA-DQB; for i=1:7 %无飞机时间段赋0 DQ(i)=0; end %% 绘图(位于红线下的为可选择返回市区B方案) x=8:24; v=DQ(8:24); xq=8:0.01:24; vq=interp1(x,v,xq,'spline');%一维三次样条插值 figure plot(xq,vq,'linewidth',2) hold on plot([0 24],[0 0],'linewidth',2);%%10:22~12:38可直接去市区接单(边界值) xlabel('t (h)','FontSize',26); ylabel('DQ (¥/min)','FontSize',26); set(gca,'FontSize',26,'linewidth',1); set(gca,'XTick',0:2:24); axis([0 24 -0.5 1])
2.2 代入t2到模型进行计算(灵敏度分析)MATLAB程序
%% 代入t2到模型进行计算(灵敏度分析) lmd1=1.0; lmd2=1.0; lmd3=1.0; lmd4=1.0; lmd5=1.0; load('wlmd');%载入等待时间t2 t1=5;%送客区至排队区时间 t2=w;%排队区等待时间 t3=30;%向市中心运客时间 t4=t3;%向市中心运客时间 t5=t1+t2+t3-t4;%在市中心运客时间 s=0.5*lmd5;%每公里0.5元的燃油费 x1=8*lmd4;%虹桥机场到市区的距离 v=40/60;%出租车的均速40KM/h m=4*lmd3;%出租车平均每公里收费 ph1=zeros(24,1);%里程利用率 for tt=1:4 ph1(tt)=0.30; end for tt=5:6 ph1(tt)=0.45; end for tt=7:10 ph1(tt)=0.55; end for tt=11:12 ph1(tt)=0.75; end for tt=13:14 ph1(tt)=0.55; end for tt=15:17 ph1(tt)=0.45; end for tt=18:19 ph1(tt)=0.55; end for tt=20:21 ph1(tt)=0.35; end for tt=22:24 ph1(tt)=0.30; end ph1=ph1*lmd2; tA=zeros(24,1); DQA=zeros(24,1); tB=zeros(24,1); PB=zeros(24,1); CB=zeros(24,1); QB=zeros(24,1); DQB=zeros(24,1); for t=1:24 %A方案 tA(t)=t1+t2(t)+t3;%A方案的总时间 PA=x1*m;%A方案载客至市区赚的钱 CA=s*x1;%A方案油耗 QA=PA-CA;%A方案收益 DQA(t)=QA/tA(t);%A方案单位时间收益 %B方案 tB(t)=t4+t5(t);%B方案的总时间 PB(t)=ph1(t)*(t5(t)*v)*m;%B方案载客至市区赚的钱 CB(t)=s*(x1+t5(t)*v);%B方案油耗 QB(t)=PB(t)-CB(t);%B方案收益 DQB(t)=QB(t)/tB(t);%B方案单位时间收益 end DQ=DQA-DQB; for i=1:7 %无飞机时间段赋0 DQ(i)=0; end %% 绘图(位于红线下的为可选择返回市区B方案) x=8:24; v=DQ(8:24); xq=8:0.01:24; vq=interp1(x,v,xq,'spline');%一维三次样条插值 ind0=find(vq<0);%找出A中小于0的位置标号 count=size(ind0,2); t0=count*0.6;%前往市区时间段长度(min) plot(xq,vq,'linewidth',2) hold on plot([0 24],[0 0],'linewidth',2); xlabel('t (h)','FontSize',26); ylabel('DQ (¥/min)','FontSize',26); set(gca,'FontSize',26,'linewidth',1); set(gca,'XTick',0:2:24); axis([0 24 -0.5 1])
三、第3问
3.1 第3问求解MATLAB程序
%第3问MATLAB程序 b1=0.1; b2=0.2; b3=0.3; a1=@(k) 1-exp(-b1*k); a2=@(k) 1-exp(-b2*k); a3=@(k) 1-exp(-b3*k); Q=zeros(1,1); for k=1:5 t1=120;%出租车从蓄车池到上车点时间 t2=30;%乘客从上车点到上车时间 t3=60;%出租车从上车点启动离开时间 T1=t1+(k-1)*a1(k)*t1;%T1表达式 T2=t2+(k-1)*a2(k)*t2;%T2表达式 T3=t3+(k-1)*a3(k)*t3;%T3表达式 Q(k,1)=(T1+T2+T3)/(2*k); end
3.2 第3问仿真MATLAB程序
%% 灵敏度分析 k=66.1; s=1; %服务台的个数 mu=3600/k; %单位时间内能服务的顾客数 lambda=32; %单位时间内到达的顾客数 ro=lambda/mu; ros=ro/s; sum1=0; for i=0:(s-1) sum1=sum1+ro.^i/factorial(i); end sum2=ro.^s/factorial(s)/(1-ros); p0=1/(sum1+sum2); p=ro.^s.*p0/factorial(s)/(1-ros); Lq=p.*ros/(1-ros); L=Lq+ro; W=L/lambda; Wq=Lq/lambda; fprintf('排队等待的平均人数为%5.2f人\n',Lq) fprintf('系统内的平均人数为%5.2f人\n',L) fprintf('平均逗留时间为%5.2f分钟\n',W*60) fprintf('平均等待时间为%5.2f分种\n',Wq*60)
四、第4问
4.1 第4问仿真MATLAB程序
%第4问MATLAB程序 P=zeros(100,1); lamda0=zeros(24,1); for t=1:24 for lamda=0.01:0.01:1 T1=50;%送长途时间 m=5;%出租车平均每公里收费 v=40/60;%出租车的均速40KM/h m1=0.5;%%每公里0.5元的燃油费 ph1=zeros(24,1);%里程利用率 for tt=1:4 ph1(tt)=0.30; end for tt=5:6 ph1(tt)=0.45; end for tt=7:10 ph1(tt)=0.55; end for tt=11:12 ph1(tt)=0.75; end for tt=13:14 ph1(tt)=0.55; end for tt=15:17 ph1(tt)=0.45; end for tt=18:19 ph1(tt)=0.55; end for tt=20:21 ph1(tt)=0.35; end for tt=22:24 ph1(tt)=0.30; end T0=15;%排队时间 T2=lamda*T1;%送短途时间 P1=v*T1*(m-m1)+2*T2*v*(ph1(t)*m-m1); P2=(T1+T2)*v*m-(2*T2+T1)*v*m1; % P(round(lamda*100),1)=abs(P1-P2); P(round(lamda*100),1)=abs((P1+P2)/(T0+2*T2+T1)-P1/(T0+T1)); end [lamda00]=find(P==min(P)); lamda0(t,1)=lamda00/100; end