代码手册之事件-活动网络模型(一)

vlambda
2021-06-06

代码手册之事件-活动网络模型(一)

该模型的求解程序主要由五个子程序组成。主程序文件main.m为程序的执行入口，子程序文件Event_Production.m、Activity_Production.m、Optimize.m、draw_timetable.m为主程序中所用到的自定义函数程序。

主程序

主程序main分为四部分：输入参数、生成集合、求解优化模型、绘制列车运行图。其中，输入参数有两种方式：一种采用直接对参数赋值的方式，如车辆数目、时刻表周期长度等等；另外一种采用导入数据的方式，使用内置的importdata函数导入文本数据，如车站的信息数据、线路信息数据等。在生成集合部分，通过自定义的两个子函数Event_Production.m、Activity_Production.m分别生成模型中需要的事件集合和活动集合。在求解优化模型部分，通过调用cplex求解原数学优化模型。最后一部分是绘制时刻表部分，通过自定义函数输入结果参数，即可绘制列车运行图。由于matlab在绘图时颜色较为单调，也可将计算结果输出，使用其他软件绘图。

1.1 main.m

主程序代码如下：

clc%----------输入参数-------------station=importdata('station.txt');station_name=importdata('station_name.txt');travel_time_min=importdata('travel_time_min.txt');travel_time_max=importdata('travel_time_max.txt');num_station=size(station,1);K=7; %列车数目t_plan=60; %规划时刻表周期t_min_w=1; %最小等待时间t_max_w=3; %最大等待时间t_max_headway=5; %最大车头时距t_min_headway=3; %最小车头时距%station中的车站属性车站0表示起始车站,1表示终点车站，2表示中间车站（需要停站），3表示中间车站（直接通过）。%在事件集合E中：0表示到达，1表示出发，2表示通过%------------生成集合---------------[E,E_arr,E_dep,E_pass] = Event_Production(K,station); %生成三类事件：出发事件、到达事件、通过事件[A_run,A_wait,A_headway]=Activity_Production(station,E); %生成三类活动：运行活动、等待活动、车头时距活动%----------求解优化模型--------------[result_v,result_z] = Optimize(E,t_plan,A_run,A_wait,A_headway,travel_time_min,travel_time_max,t_min_w, t_max_w,t_max_headway,t_min_headway);%----------绘制时刻表---------------[Time,Space] = draw_timetable(station_name,E,K,result_v);

子程序

子程序是根据主程序运行需要自定义的函数。在主程序中通过调用子程序完成整体计算。定义子程序的优势是可简化主程序的流程，对于某些重复性的运算操作可通过调用自定义的函数简化代码量。除此之外，通过封装程序也可提高程序可读性，更加便于读者阅读理解。本文中子程序分为四部分：事件产生子程序、活动产生子程序、求解优化模型子程序、绘制列车运行图子程序。

2.1 Event_Production.m

事件产生子程序（Event_Production）的输入参数为车辆数目K和车站信息station。该子程序输出三个集合：总事件集合E，到达事件集合E_arr和通过事件集合E_pass。该程序的编程思路为：按照列车的次序，重复填充在每辆地铁列车每个车站的事件。填充完第一辆列车的事件，接着填充第二辆，直到填充完所有列车的事件。需要注意的是：对于出发车站只有出发事件；对于到达车站，只有到达事件，对于通过车站，只有通过事件；而对于中间的非通过车站，需要先填充到达事件，再填充出发事件。最后集合E的元素通过将该三个事件的集合按照事件的次序编号汇总到该集合中。

事件产生的子程序代码如下：

function [E,E_arr,E_dep,E_pass] = Event_Production(K,station)num_station=size(station,1);temp=1; %事件计数temp_arr=1; %到达事件序号计数temp_dep=1; %出发事件序号计数temp_pass=1; %通过事件序号计数for k=1:K for i=1:num_station if station(i,2)==0 E_dep(1,temp_dep)=temp; E_dep(2,temp_dep)=i; temp_dep=temp_dep+1; temp=temp+1; elseif station(i,2)==1 E_arr(1,temp_arr)=temp; E_arr(2,temp_arr)=i; temp_arr=temp_arr+1; temp=temp+1; elseif station(i,2)==2 E_arr(1,temp_arr)=temp; E_arr(2,temp_arr)=i; temp_arr=temp_arr+1; E_dep(1,temp_dep)=temp+1; E_dep(2,temp_dep)=i; temp_dep=temp_dep+1; temp=temp+2; else  E_pass(1,temp_pass)=temp; E_pass(2,temp_pass)=i; temp_pass=temp_pass+1; temp=temp+1; end endend%填充事件集合EE=zeros(size(E_arr,2)+size(E_dep,2)+size(E_pass,2),3);for i=1:size(E_arr,2) E(E_arr(1,i),1)=E_arr(1,i); E(E_arr(1,i),2)=E_arr(2,i); E(E_arr(1,i),3)=0;endfor i=1:size(E_dep,2) E(E_dep(1,i),1)=E_dep(1,i); E(E_dep(1,i),2)=E_dep(2,i); E(E_dep(1,i),3)=1;endfor i=1:size(E_pass,2) E(E_pass(1,i),1)=E_pass(1,i); E(E_pass(1,i),2)=E_pass(2,i); E(E_pass(1,i),3)=2;endend

2.2 Activity_Production.m

活动产生子程序（Activity_Production）的输入参数为车站信息station和事件集合E。该子程序产生三个事件集合：列车运行活动集合A_run、列车等待活动集合A_wait和列车车头时距活动集合A_headway。该程序的编程思路为：对于E中的每一个事件，自该事件伊始，寻找下一个满足可组成活动条件的事件与之组成对应的运行活动、等待活动或通过活动，将该事件与找到的事件按照实数对(e1,e2)的形式存储到对应的活动集合中。重复该过程直到遍历所有集合E中的事件。需要注意的是：对于起点车站和通过的中间车站而言，没有等待活动；对于终点车站，没有列车运行活动和等待活动。

活动产生的子程序代码如下：

function [A_run,A_wait,A_headway] = Activity_Production(station,E)num_station=size(station,1);temp_run=1; %run活动计数temp_wait=1; %wait活动计数temp_headway=1; %headway活动计数for i=1:size(E,1) for j=i+1:size(E,1) if E(i,3)==1&&E(j,3)==0||E(i,3)==1&&E(j,3)==2||E(i,3)==2&&E(j,3)==0||E(i,3)==2&&E(j,3)==2 A_run(1,temp_run)=i; A_run(2,temp_run)=j; temp_run=temp_run+1; break end If E(i,3)==2&&E(j,3)==2||E(i,3)==2&&E(j,3)==1||E(i,3)==0&&E(j,3)==2&&E(i,2)~=num_station||E(i,3)==0&&E(j,3)==1&&E(i,2)~=num_station A_wait(1,temp_wait)=i; A_wait(2,temp_wait)=j; temp_wait=temp_wait+1; break end endendfor i=1:size(E,1) for j=i+1:size(E,1) If E(i,3)==1&&E(j,3)==1&&E(i,2)==1&&E(j,2)==1||E(i,3)==0&&E(j,3)==0&&E(i,2)==num_station&&E(j,2)==num_station||E(i,3)==1&&E(j,3)==0&&E(i,2)==E(j,2)&&E(i,2)~=1&&E(i,2)~=num_station||E(i,3)==2&&E(j,3)==2&&E(i,2)==E(j,2)&&E(i,2)~=1&&E(i,2)~=num_station A_headway(1,temp_headway)=i; A_headway(2,temp_headway)=j; temp_headway=temp_headway+1; break end endendEnd

2.3 Optimize.m

yalmip是由Lofberg开发的一种免费的优化求解工具，其最大特色在于集成许多外部的最优化求解器，形成一种统一的建模求解语言，并提供了Matlab的调用API。本文使用matlab的yalmip工具箱调用cplex，基本语法、环境及配置可参考相关的资料。

求解优化模型子程序（Optimize）的输入参数为所有的事件及活动集合、车头时距、列车区间运行时间、车站等待时间，程序的输出为事件的发生时刻、目标函数值。函数的内部的程序按照优化模型的语法要求书写即可。

求解优化模型的子程序代码如下：

function [result_v,result_z] = Optimize(E,t_plan,A_run,A_wait,A_headway,travel_time_min,travel_time_max, t_min_w,t_max_w,t_max_headway,t_min_headway)num_E=size(E,1);%定义决策变量v=sdpvar(1,size(E,1));%添加约束条件C=[]; %约束条件集合%单个事件约束for i=1:size(E,1) s=v(1,i); C=[C;s>=0];endfor i=1:size(E,1) s=v(1,i);    C=[C;s<=t_plan];           %规划1个小时的时刻表end%同一车次的事件约束for i=1:size(A_run,2) s=v(1,A_run(2,i))-v(1,A_run(1,i)); C=[C;s<=travel_time_max(E(A_run(1,i),2),E(A_run(2,i),2))];endfor i=1:size(A_run,2) s=v(1,A_run(2,i))-v(1,A_run(1,i)); C=[C;s>=travel_time_min(E(A_run(1,i),2),E(A_run(2,i),2))];end
for i=1:size(A_wait,2) s=v(1,A_wait(2,i))-v(1,A_wait(1,i)); C=[C;s<=t_max_w];endfor i=1:size(A_wait,2) s=v(1,A_wait(2,i))-v(1,A_wait(1,i)); C=[C;s>=t_min_w];end%不同车次的事件约束for i=1:size(A_headway,2) s=v(1,A_headway(2,i))-v(1,A_headway(1,i)); C=[C;s<=t_max_headway];endfor i=1:size(A_headway,2) s=v(1,A_headway(2,i))-v(1,A_headway(1,i)); C=[C;s>=t_min_headway];end
% 配置ops = sdpsettings('verbose',0,'solver','cplex');%z=0;%for i=1:size(E,1)% z=z+v(1,i);%endz=v(1,size(E,1));% 目标函数reuslt = optimize(C,z); %默认求最小值if reuslt.problem == 0 % problem =0 代表求解成功; vresult=value(v); zresult=value(z);else disp('求解出错');endresult_v=vresult;result_z=zresult;end

2.4 draw_timetable.m

绘制列车运行图子程序（draw_timetable）的输入参数为车站名称station _na me, E、列车运行线数量K、事件的发生时刻vresult，该程序无数值结果输出，只是在函数内部运行绘制图像。绘制列车时刻表的基本原理与绘制折线图的原理相同。值得注意的是，对于不同列车的运行线要区别绘制。

绘制列车运行图子程序（draw_timetable）代码如下：

function [ Time,Space ] = draw_timetable(station_name,E,K,vresult)num_E=size(E,1);num_e_train=num_E/K;Time=zeros(K,num_e_train);Space=zeros(K,num_e_train);temp=1;k=1;for i=1:size(E,1) if temp<=num_e_train Time(k,temp)=vresult(i); Space(k,temp)=E(i,2); temp=temp+1; end if temp>num_e_train temp=1; k=k+1; endendfor k=1:K plot(Time(k,:),Space(k,:)); hold onendfor i=1:size(station_name,1) L(i)=i;endset(gca,'ytick',L);set(gca,'yticklabel',station_name');end