CN110457012A

CN110457012A - Multiple attribute decision making (MADM) software for sustainable transport network design

Info

Publication number: CN110457012A
Application number: CN201910627204.0A
Authority: CN
Inventors: 林宏志; 赵宇轩
Original assignee: Southeast University
Current assignee: Southeast University
Priority date: 2019-06-27
Filing date: 2019-06-27
Publication date: 2019-11-15

Abstract

The present invention realizes a kind of multiple attribute decision making (MADM) software for sustainable transport network design using R language.Upper layer policy maker generates a variety of feasible design schemes using Dirichlet distribution first.Then, lower layer's traveler makes a series of behavior reaction for every kind of feasible program, this is the four stage progressive die types for having feedback mechanism.The feedback procedure of the lower layer can converge to traffic system equilibrium, and this equilibrium result, which further feed back, arrives upper layer, with multiple attributes such as economy, environment, society, the safety of estimating transportation network.After calculating multiple attribute values finally, for each scheme, optimal feasible program is determined using multiple attributive decision making method.Due to using free R language of increasing income, this software is easily operated and propagates.The present invention for given Investment in Basic Facilities for Traffic budget in the case where to design sustainable transportation network particularly useful.

Description

Multiple attribute decision making (MADM) software for sustainable transport network design

Technical field:

The present invention has write a kind of multiple attribute decision making (MADM) software for sustainable transport network design using R language, belongs to Traffic engineering technical field.

Background technique:

Traffic system plays key effect in sustainable development, and therefore, sustainable city traffic has recently become one It is a important to grind the field of making internal disorder or usurp.Sustainable transport system is usually characterized by three aspects: economy, environment and society.However, existing Traffic study only focus on one or most two dimensions^[1].It is only least a portion of to study while considering all three dimensions. It should be noted that might have many measuring methods under each dimension of three dimensions.For example, economic dimension is in addition to system Total travel time except, can also include Consumer Surplus, network marginal capacity, overall cost of ownership, total charge income etc..This Outside, safety is also the important indicator of sustainability.However, it is in sustainable transport ignored a very long time, this is Very strange.In fact, traffic safety is influenced by traffic system each stage, including planning, design, construction and fortune Battalion.In view of Transportation Network Design Problem is to the profound impact for forming traffic system, it should sufficiently pay attention to traffic peace in network design phase Entirely.So-called active traffic safety planning just refers to the traffic safety examination in network design phase, to improve traffic system in advance Level of security.In view of traffic safety to the importance of sustainable development, this software is regarded as the 4th of sustainable transport Dimension is simultaneously brought into overall merit.

It making internal disorder or usurp although widely grind to one or two dimension, few people discuss three dimensions, needless to say Four dimensions.Due to not looked back them here there are many making internal disorder or usurp about two grinding for sustainable dimension, only summarize Research to three dimensions.Chen and Xu^[2]Three are considered in the double-deck goal programming frame of Traffic Micro-simulation Target: efficiency, environment and justice.Wherein, justice is a kind of main social factor Measurement Method.Yin etc.^[3]Propose one The double-deck Goal programming Model, with solving road charge collection pricing and traffic capacity investment problem, it considers following three dimensions: Economy, environment and justice.Feng and Timmermans^[4]It has been handled between mobility and fairness under environmental carrying capacity constraint It accepts or rejects.Li and Ge^[5]A kind of Bi-level Programming Models for Traffic Signal Timing design are proposed, being at the middle and upper levels one has The Multiobjective Programming of fairness constraints, i.e. maximization network marginal capacity simultaneously minimize total traffic emission.Szeto etc.^[1]Also it mentions A multiple target bi-level optimal model is gone out, to consider three dimensions in road network design, including economy, environment and society Meeting.Wang etc.^[6]The Bi-level Programming Models of a sustainable transport network design with stochastic demand are proposed, at the middle and upper levels It is maximization network marginal capacity problem, and with discharge and fair opportunity constraint.Sun etc.^[7]It is proposed a multiple target bilayer Plan model, to determine the best congestion price in terms of efficiency, environment and justice.

Compared with traditional three dimension of sustainable transport, safety is seldom evaluated.Yang etc.^[8]It proposes A kind of Bi-level Programming Models design effective, safe and environmentally friendly traffic system by section rate limitation.Haas and Bekhor^[9]The Bi-objective Bi-level Programming Models an of Transportation Network Design Problem are constructed, for minimizing system total travel time and maximization Road safety is horizontal.Wherein, security performance performance is a negative binomial model, lays particular emphasis on and submits interpreter to section therefore number Estimation.Possel etc.^[10]The multiple target bi-level optimal model for being used for Traffic Micro-simulation is defined, to minimize traffic Exhaust emissions, traffic accident sum and total travel time.There is certain accident rate in text for each type of section, by thing Therefore the traffic accident number in the section can be obtained multiplied by the vehicle total mileage (VMT) travelled on section in one day for rate.To this number Value is summed on all sections, that is, can determine the security level of transportation network.

Although the planning about sustainable transport system, design, management, operation and the various of control have been proposed to arrange It applies, but Traffic Micro-simulation (NDP) is most widely used one.It is intended to efficiently use limited resource (such as soil And investment) Lai Youhua traffic system performance, while by the behavior reaction of traveler (such as destination selection, mode of transportation select Select, route selection and departure time) it takes explicitly into account in design scheme.In existing software, Transportation Network Design Problem Main target is to mitigate traffic congestion, that is, minimizes total travel time or maximization network marginal capacity, and fails to fully consider Other dimensions of sustainability and conflicting with each other between them.

Bibliography:

[1] Szeto W Y, Jiang Y, Wang D Z W, Sumalee A.A Sustainable Road Network Design Problem with Land Use Transportation Interaction over Time[J].Netw Spat.Econ, 2015,15 (3): 791-822

[2] Chen A, Xu X.Goal programming approach to solving network design problem with multiple objectives and demand uncertainty[J].Expert Systems With Applications, 2012,39 (4): 4160-4170

[3] Yin Y, Li Z C, Lam W H K, Choi K.Sustainable Toll Pricing and Capacity Investment in a Congested Road Network:A Goal Programming Approach [J] .J.Transp.Eng, 2014,140 (12): 10

[4] Feng T, Timmermans H J P.Trade-offs between mobility and equity Maximization under environmental capacity constraints:A case study of an Integrated multi-objective model [J] .Transp.Res.Pt.C-Emerg.Technol, 2014,43: 267-279

[5] Li Z C, Ge X Y.Traffic signal timing problems with environmental And equity considerations [J] .J.Adv.Transp, 2014,48 (8): 1066-1086

[6] Wang H, Lam W H K, Zhang X N, Shno H.Sustainable Transportation Network Design with Stochastic Demands and Chance Constraints[J] .International Journal of Sustainable Transportation, 2015,9 (2): 126-144

[7] Sun X, Liu Z Y, Thompson R G, Bie Y M, Weng J X, Chen S Y.A multi- objective model for cordon-based congestion pricing schemes with nonlinear Distance tolls [J] .Journal of Central South University, 2016,23 (5), 1273-1282

[8] Yang Y N, Lu H P, Yin Y F, Yang H.Optimization of Variable Speed Limits for Efficient, Safe, and Sustainable Mobility [J] .Transportation Research Record, 2013 (2333): 37-45

[9] Haas I, Bekhor S.Network design problem considering system time Minimization and road safety maximization:formulation and solution approaches [J] .Transportmetrica A, 2017,13 (9): 829-851

[10] Possel B, Wismans L J J, Van Berkum E C, Bliemer M C J.The multi- objective network design problem using minimizing externalities as Objectives:comparison of a genetic algorithm and simulated annealing Framework [J] .Transportation, 2018,45 (2): 545-572

Summary of the invention:

Technical problem: current Transportation Network Design Problem software at least has three.Firstly, pacifying in sustainable transport Full dimension is seldom taken into account.In this software, measured using the security performance function of highway safety handbook (HSM) section and The security performance of road network.Secondly, the model of lower layer is usually the user equilibrium model for being used for traffic flow distribution, this is quite to limit to , because the decision behavior of traveler includes traffic generation, traffic distribution, traffic modal splitting and traffic flow distribution.Therefore, originally Software uses a kind of four-phase model of band feedback, and the traffic system to realize lower layer is balanced, to sufficiently show traveler Decision behavior.Finally, the complexity of multiple target bilayer model is well-known.In order to solve proposed bilayer model, This software is distributed in Dirichlet, the straight average method (MSA) that weight is successively decreased, Frank-Wolfe algorithm and dijkstra's algorithm On the basis of, devise multiple attribute decision making (MADM) (MADM) method.

Technical solution: the present invention has write a kind of multiple attribute decision making (MADM) for sustainable transport network design using R language Software, the software specifically includes the following steps:

(1) overall algorithm designs

Step 1: upper layer policy maker generates a random transportation network traffic capacity using Dirichlet distribution Lifting scheme Δ c, Protocol Numbers m=0, lower layer's traveler make a series of behavior reactions；

Step 2: underlying model is the progressive die that traffic generation, traffic distribution, traffic modal splitting and traffic flow distribute Type reaches traffic system balance by feedback iteration, can calculate the road section traffic volume flow and transit time when equilibrium state；

Step 3: four attribute values of sustainable transport when calculating traffic system balance:

Economic attribution

Environment attribute

Social property

Security attribute；

Step 4: return step one is divided into another random transportation network traffic capacity using Dirichlet and is promoted Scheme Δ c, Protocol Numbers m=m+1 when m=M, stop circulation when scheme number reaches (M >=200) i.e. M predetermined, Go to step 5；

Step 5: optimal Transportation Network Design Problem scheme is found out using multiple attributive decision making method.

(2) the specific calculating process of step 2

Step 2 is the underlying model in traffic system balance and bilayer model.For each design scheme on upper layer, There is corresponding lower layer's equilibrium state.Algorithm flow chart is as shown in Fig. 2, detailed calculating process is as follows:

Step 1: obtaining traffic capacity Lifting scheme Δ c from Dirichlet distribution Dir (α)；

Step 2: by being uniformly distributed initialization traffic distribution matrixN=0 is set, indicates the number of iterations；

Step 3: user equilibrium is based on by Frank-Wolfe algorithm, traffic distribution matrix is distributed into transportation network, with The magnitude of traffic flow on each section a and travel time are calculated, later, the most short travel time between starting point i and destination j, i.e.,It can be acquired by dijkstra's algorithm；

Step 4: being based onTraffic distribution matrix is updated using destination preference pattern

Step 5: the MSA to be successively decreased using weight is to traffic distribution matrixWithIt is averaging

Step 6: the convergence of traffic distribution matrix is checked using opposite root square error (RRSE)

If meeting the condition of convergence, 8 are gone to step, otherwise goes to step 7；

Step 7: enabling n=n+1, user equilibrium is then based on for traffic distribution matrix by Frank-Wolfe algorithmPoint Dispensing transportation network, to calculate the magnitude of traffic flow and the travel time on each section a, later, between starting point i and destination j The most short travel time, i.e.,It can be calculated by dijkstra's algorithm, feedback arrives step 4；

Step 8: output traffic distribution matrixAnd the magnitude of traffic flow v on a of section_aBetween starting point i and destination j Travel time

(3) the specific calculating process of step 3

Step 1: economic attribution, with total travel time (TTT) calculating

Wherein A is the section set on network；t_aIt is the travel time on a of section, as unit of minute；v_aIt is on a of section Flow, as unit of pcu/h.

Step 2；Environment attribute is calculated with the total CO discharge (TCO) of vehicle

Wherein e_a(v_a) indicate section a Plays vehicle equivalent (pcu) CO contaminant capacity.Following macroscopical function is with generation One function of table is usually used in estimating the vehicle CO discharge based on section:

Wherein l_a, t_aAnd e_aRespectively with km, minute and gram be per hour unit.

Step 3: fair property, the maximum rate of O-D travel time is determined before and after being enhanced with road passage capability Adopted space is fair:

Wherein π_w(Δ c) and π_wIt (0) is when O-D goes on a journey to the minimum of w with before after implementing network design scheme Δ c Between；W is O-D pairs of set.

Step 4: security attribute obtains road section traffic volume according to the security performance function suggested in highway safety handbook (HSM) Accident number:

SPF_a=[β₁×(v_a/c_a)²+β₂×(v_a/c_a)+β₃]×v_a×l_a×365×10^-8

Wherein, SPF_aIndicate average accident number (acc) per hour of annual section a；v_aIt is that the average of section a is handed over per hour Flux (pcu/h)；l_aIt is the length (km) of section a；c_aIt is the practical road section capacity (pcu/h) on a of section；And β₁, β₂ And β₃It is the parameter that traffic survey data calibration can be used.Therefore, the summation of all road section traffic volume accidents can be used as network peace Full measurement is as follows:

(4) the specific calculating process of step 5

Step 1: calculate four attributes of traffic system equilibrium under each scheme, i.e., economic, environment, social and safety；

Step 2: since the dimension of attribute is different, they must carry out standard before integration, be the bigger the better for numerical value Attribute, standardization formula be

Wherein y_ijIt is the value of the attribute j of mode i；It is the worst case of attribute j；It is the optimal cases of attribute j； z_ijIt is the standardized value of the attribute j of mode i, the attribute the smaller the better for numerical value, standardizing formula is

Step 3: the weight w of attribute is determined using feature vector method_j；

Step 4: the comprehensive score of network design scheme i is

Step 5: after calculating M all scores, network design scheme can be ranked up, and can selected The network design of highest scoring is selected, here it is optimal designs.

(5) programming

# is used for the multiple attribute decision making (MADM) software of sustainable transport network design: bilayer model and algorithm

# step 1: initialization, by format input data and necessary packet；

#1.1 load calculates the packet of shortest path, prepares to call dijkstra shortest path first, notices that igraph packet is first Secondary use needs to install, and is then just capable of calling；

#install.packages (" igraph ") # installs igraph packet

library(igraph)

Options (digits=3)

#1.2 creates the distance matrix of figure, and comprising all candidate road sections, first is classified as section label (Road), secondary series Piont mark (Road origin) is played for section, third is classified as road segment end label (Road destination), and the 4th is classified as this The section free flow time (free flow time), the 5th is classified as road passage capability (capacity), and the 6th to be classified as road long It spends (length), this is sentenced in traffic flow distribution for common Nguyen-Dupuis network, and detailed parameter setting can refer to journey Preface shelves；

# can also be replicated in Excel, then be executed

#e=read.delim (" clipboard ", header=F)

E=matrix (c (1,1,5,7.0,900,4.00,2,1,12,9.0,700,4.00,3,4,5,9.0,700, 4.00,4,4,9,12.0,900,7.00,5,5,6,3.0,800,2.00,6,5,9,9.0,600,4.00,7,6,7,5.0,900, 4.00,8,6,10,13.0,500,8.00,9,7,8,5.0,300,4.00,10,7,11,9.0,400,5.00,11,8,2,9.0, 700,5.00,12,9,10,10.0,700,6.00,13,9,13,9.0,600,5.00,14,10,11,6.0,700,4.00,15, 11,2,9.0,700,5.00,16,11,3,8.0,700,4.00,17,12,6,7.0,300,4.00,18,12,8,14.0,700, 9.00,19,13,3,11.0,700,6.00), ncol=6, byrow=T)

Colnames (e)=c (" Road ", " Road origin ", " Road destination ", " Free Time ", " Road capacity ", " Road length ")

#e# is used to check the breakpoint of program

#1.3 inputs initial transport need matrix d0, and first is classified as the label (OD pair) of origin and destination pair, and second has been classified as Piont mark (origin), third are classified as terminal label (destination), and the 4th is classified as transport need (demand)；

Tge=3000# total status transport need

D0=matrix (c (1,1,2,0.2*tge, 2,1,3,0.4*tge, 3,4,2,0.3*tge, 4,4,3,0.1* Tge), ncol=4, byrow=T) # original allocation scheme

Colnames (d0)=c (" OD pair ", " Origin ", " Destination ", " Demand ")

#d0# is used to check the breakpoint of program

The customized Frank-Wolfe algorithmic function of # pays attention to the requirement matrix d form such as d0 of input, transportation network e's Form e for example above, relative error 0.001；

Fw=function (e, d)

{

#1.4 calculates each OD pairs of shortest path and path flow according to the path free flow time

G=add.edges (graph.empty (13), t (e [, 2:3]), and weight=e [, 4]) # creation figure, 13 be section Point number, using the time as weight the length of non-path

B12=get.shortest.paths (g, from=" 1 ", to=" 2 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 1 to the shortest path of terminal 2

B13=get.shortest.paths (g, from="] ", to=" 3 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 1 to the shortest path of terminal 3

B42=get.shortest.paths (g, from=" 4 ", to=" 2 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 4 to the shortest path of terminal 2

B43=get.shortest.paths (g, from=" 4 ", to=" 3 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 4 to the shortest path of terminal 3

# creates a matrix, for saving each OD pairs of shortest path and flow

V=cbind (e [, 1])

The initial each OD of st0=numeric (4) # storage is to most short running time

Colnames (V)=" Road "

V

#OD to 12 shortest path and flow

Sp12=as.vector (b12) # is converted into section label (Road)

St0 [1]=sum (e [sp12,4]) each section temporal summation of #

X12=cbind (e [sp12,1], rep (d [Isosorbide-5-Nitrae], length (sp12))) section # label and flow, in algorithm Iteration starting point

Colnames (x12)=c (" Road ", " V12 ")

x12

It is the special matrix for saving iteration starting point that V=merge (V, x12, by=" Road ", all=TRUE) #, which defines V,

V [is.na (V)]=0

V

#OD to 13 shortest path and flow

Sp13=as.vector (b13) # is converted into section label (Road)

St0 [2]=sum (e [sp13,4]) each section temporal summation of #

X13=cbind (e [sp13,1], rep (d [2,4], length (sp13))) section # label and flow, in algorithm Iteration starting point

Colnames (x13)=c (" Road ", " V13 ")

x13

It is the special matrix for saving iteration starting point that V=merge (V, x13, by=" Road ", all=TRUE) #, which defines V,

V [is.na (V)]=0

V

#OD to 42 shortest path and flow

Sp42=as.vector (b42) # is converted into section label (Road)

St0 [3]=sum (e [sp42,4]) each section temporal summation of #

X42=cbind (e [sp42,1], rep (d [3,4], length (sp42))) section # label and flow, in algorithm Iteration starting point

Colnames (x42)=c (" Road ", " V42 ")

x42

It is the special matrix for saving iteration starting point that V=merge (V, x42, by=" Road ", all=TRUE) #, which defines V,

V [is.na (V)]=0

V

#OD to 43 shortest path and flow

Sp43=as.vectot (043) # is converted into section label (Road)

St0 [4]=sum (e [sp43,4]) each section temporal summation of #

X43=cbind (e [sp43,1], rep (d [4,4], length (sp43))) section # label and flow, in algorithm Iteration starting point

Colnames (x43)=c (" Road ", " V43 ")

x43

It is the special matrix for saving iteration starting point that V=merge (V, x43, by=" Road ", all=TRUE) #, which defines V,

V [is.na (V)]=0

V

# obtains initial flow when the flow summation on all shortest paths

VS=rowSums (V [, seq (ncol (V) -3, ncol (V))])

VS

# step 2: the impedance in each section is updated

T0=e [, 4] and the # free flow time

C=e [, 5] and # road passage capability

A=0.15

B=4

Tp=function (v)

t0*(1+a*(v/c)^b)

}

repeat{

# step 3: next iteration direction is found

G2=add.edges (graph.empty (13), t (e [, 2:3]), weight=tp (VS)) # structural map, 13 are The number of node updates link proportion

B12=get.shortest.paths (g2, from=" 1 ", to=" 2 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 1 to the shortest path of terminal 2

B13=get.shortest.paths (g2, from=" 1 ", to=" 3 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 1 to the shortest path of terminal 3

B42=get.shortest.paths (g2, from=" 4 ", to=" 2 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 4 to the shortest path of terminal 2

B43=get.shortest.paths (g2, from=" 4 ", to=" 3 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 4 to the shortest path of terminal 3

# creates a provisional matrix, for saving each OD pairs of shortest path and flow

V=cbind (e [, 1])

Colnames (V)=" Road "

V

#OD to 12 shortest path and flow

Sp12=as.vector (b12) # is converted into section label (Road)

Colnames (x12)=c (" Road ", " V12 ")

x12

V [is.na (V)]=0

V

#OD to 13 shortest path and flow

Sp13=as.vector (b13) # is converted into section label (Road)

Colnames (x13)=c (" Road ", " V13 ")

x13

V [is.na (V)]=0

V

#OD to 42 shortest path and flow

Sp42=as.vector (b42) # is converted into section label (Road)

Colnames (x42)=c (" Road ", " V42 ")

x42

V [is.na (V)]=0

V

#OD to 43 shortest path and flow

Sp43=as.vector (b43) # is converted into section label (Road)

Colnames (x43)=c (" Road ", " V43 ")

x43

V [is.na (V)]=0

V

# obtains iteration direction when the flow summation on all shortest paths

VS2=rowSums (V [, seq (ncol (V) -3, ncol (V))])

VS2

# step 4: iteration step length is calculated

Step=function (lamda)

X2=VS2

X1=VS

Q=x1+lamda* (x2-x1)

sum((x2-x1)*tp(q))

}

#lamda=uniroot (step, c (0,1)) $ root# pays attention to the value range of lamda, and step-length cannot be too long, Uniroot requires the functional value symbol at both ends on the contrary, some functions not necessarily meet, and may insure to look for using optimize function To the optimal value of function of a single variable；

G=function (lamda) { step (lamda) ^2 }

Lamda=optimize (g, c (0,1)) $ minimum

lamda

# step 5: new iteration starting point is determined

VS3=VS+lamda* (VS2-VS)

VS3

# step 6: test for convergence

If ((sqrt (sum ((VS3-VS) ^2))/sum (VS)) < 0.001) break

VS=VS3# substitutes origin VS with new point VS3 if being unsatisfactory for the condition of convergence, and so circulation is until convergence

}

# step 7: the eigenmatrix result and OD running time matrix u of equilibrium state is exported；

# step 7.1: flow, transit time and the speed in output each path of equilibrium state；

Result=cbind (e [, 1], round (VS, 0), tp (VS), e [, 6]/(tp (VS)/60), and e [, 5], round (VS, 0)/e [, 5])

Colnames (result)=c (" Road ", " Volume ", " Time ", " Speed ", " Road Capacity ", " Level of Service″)

# step 7.2: each OD running time matrix u is exported；

G=add.edges (graph.empty (13), t (e [, 2:3]), and weight=result [, 3]) # creation figure, 13 For the number of node, result is the matrix that step 7 generates

B13=get.shortest.paths (g, from=" 1 ", to=" 3 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 1 to the shortest path of terminal 3

B43=get.shortest.paths (g, from=" 4 ", t0=" 3 ", mode=" out ", output=" Epath ") $ epath [[1]] # is from starting point 4 to the shortest path of terminal 3

# creates a running time matrix, for saving each OD pairs of journey time, each OD journey time of original hypothesis It is 0

U=matrix (c (1,1,2,0,2,1,3,0,3,4,2,0,4,4,3,0), ncol=4, byrow=T)

#OD to 12 journey time

Sp12=as.vector (b12) # is converted into section label (Road)

U [Isosorbide-5-Nitrae]=sum (result [sp12,3]) each section temporal summation of #

#OD to 13 journey time

Sp13=as.vector (b13) # is converted into section label (Road)

U [2,4]=sum (result [sp13,3]) each section temporal summation of #

#OD to 42 journey time

Sp42=as.vector (b42) # is converted into section label (Road)

U [3,4]=sum (result [sp42,3]) each section temporal summation of #

#OD to 43 journey time

Sp43=as.vector (b43) # is converted into section label (Road)

U [4,4]=sum (result [sp43,3]) each section temporal summation of #

Accessible running time st0 between #OD couples of u=cbind (u, st0)

# exports result matrix and OD running time matrix in the form of a list

List (result, u)

}

#fw (e, d0) # is used to check the breakpoint of program

# step 8: the multinomial Logit function mlogit of destination selection is defined, is inputted as each OD running time time square Battle array u and various regions transport need sg, exports as new traffic distribution matrix；

Mlogit=function (u, sg)

{

D=nuneric (4)

D [1]=sg [1] * exp (- 0.1*u [Isosorbide-5-Nitrae])/(exp (- 0.1*u [Isosorbide-5-Nitrae])+exp (1-0.1*u [2,4]))

D [2]=sg [1] * exp (1-0.1*u [2,4])/(exp (- 0.1*u [Isosorbide-5-Nitrae])+exp (1-0.1*u [2,4]))

D [3]=sg [2] * exp (- 0.1*u [3,4])/(exp (- 0.1*u [3,4])+exp (1-0.1*u [4,4]))

D [4]=sg [2] * exp (1-0.1*u [4,4])/(exp (- 0.1*u [3,4])+exp (1-0.1*u [4,4]))

Cbind (u [, 1:3], d)

}

# step 9: traffic distribution and traffic point comprehensive under a given transport need d0 and sg and transportation network e are defined Function is balanced with alternating iteration, the running time matrix of each OD when can be in the hope of user equilibrium state for initial traffic distribution, User reselects destination according to the matrix, and new running time matrix, the mistake can be generated again for new traffic distribution Cheng Yizhi circulation carries out, until traffic distribution matrix no longer changes；

Cda=function (e, d0, sg)

K=3

repeat{

D1=mlogit (fw (e, d0) [[2]], sg)

K=k+1

If (sqrt (sum ((d1 [, 4]-d0 [, 4]) ^2))/sum (d0 [, 4]) < 0.01) break# meets certain essence Degree requires just to stop

D0 [, 4]=d0 [, 4]+(1/k) * (d1 [, 4]-d0 [, 4]) # uses iteration weighting method (Method of here Successive Average, MSA), subtracted herein as weight with the increase of cycle-index using the inverse of cycle-index It is few；

#print (d1) # is used to check the breakpoint of program

#print (k) # is for checking program

If (k==100) also jumps out circulation but without required precision is met if break# cycle-index reaches 100 times

}

#print(d1)

D2=fw (e, d1)

#print(d2)

List (d1, d2)

}

The performance of # status traffic system, the Before-After for policy compare

Ge1=c (sum (d0 [1:2,4]), sum (d0 [3:4,4])) # is used to check the breakpoint of program

Before=cda (e, d0, ge1) [[2]] [[2]] # is used to check the breakpoint of program

# step 10:Dirichlet distribution method main program；

#install.packages (" DIRECT ") # installation generates the packet of Dirichlet, needs to install for the first time, again not It needs

Library (DIRECT) # load packet

Set.seed (100) # sets random seed, the same random number generated every time in this way

Nsd=500# indicates the number of Dirichilet sample, is compared for each, finds out optimal one

The number of candicate=7# candidate road section

RD=rDirichlet (nsd, rep (1, candicate)) # generates Dirichlet random distribution

Inv=2000# investment budgey

The value added of len=inv*rD/ (0.3*c (5,8,2,3,5,5,7)) alternative road section capacity of #

St=matrix (numeric (nsd*4), nc=4) # is used to store 4 attributes of each Dirichlet scheme

It is the parameter for measuring safety below #

B1=358.6

B2=-407.7

B3=175.3

The operation starting time of # calculation procedure！！Front be all in defined function, and be not take up calculate the time；

Timestart <-Sys.time ()

For (i in (1:nsd))

{

#i=2# is for checking program

It is that a new network e is constructed to each distribution model below #

# one update scheme of Ren=cbind (c (3,4,5,9,13,16,19), len [i ,])

Colnames (Ren)=c (" Road ", " Enhancement ") # and section number combine

Ren2=merge (e, Ren, by=" Road ", all=TRUE) # merges

Ren2 [, 7] [and is.na (Ren2 [, 7])]=0# replaces NA with 0, because NA participates in the result calculated still NA

Ren2 [, 5]=Ren2 [, 5]+Ren2 [, 7] # is added with the original traffic capacity, becomes the improved traffic capacity

E1=Ren2 [, 1:6] # one new network e1

D3=cda (e1, d0, ge1) # calls fixed transport need (1200,800) and traffic distribution d0 previously defined Cda function

The value of 4 attributes is sought below #, and is stored in st matrix

#1. economic attribution: total travel time (minute)

St [i, 1]=d3 [[2]] [[1]] [, 2] %*%d3 [[2]] [[1]] [, 3]

#2. environment attribute: the discharge amount of CO

St [i, 2]=(0.2038*d3 [[2]] [[1]] [, 3] * exp (0.7962*e [, 6]/d3 [[2]] [[1]] [, 3])) %*%d3 [[2]] [[1]] [, 2]

#3. social property: the variation before and after policy

St [i, 3]=max (d3 [[2]] [[2]] [, 4]/before [, 4])

#4. security attribute: accident number

St [i, 4]=sum ((b1*d3 [[2]] [[1]] [, 6] and ^2+b2*d3 [[2]] [[1]] [, 6] and+b3) * d3 [[2]] [[1]] [, 2] * e1 [, 6] and * 365*10^ (- 8))

Print (paste (" the ", i, " 4 attributes of a scheme are respectively ", round (st [i ,], 2)))

}

# step 11: by each attribute, this 4 attributes are all the smaller the better

Zst=st

Zst [, 1]=(max (st [, 1])-st [, 1])/(max (st [, 1])-min (st [, 1]))

Zst [, 2]=(max (st [, 2])-st [, 2])/(max (st [, 2])-min (st [, 2]))

Zst [, 3]=(max (st [, 3])-st [, 3])/(max (st [, 3])-min (st [, 3]))

Zst [, 4]=(max (st [, 4])-st [, 4])/(max (st [, 4])-min (st [, 4]))

# step 12: the weight of each attribute is calculated

A=c (1,1/3,1/2,1/4,3,1,2,1,2,1/2,1,1/2,4,1,2,1)

A=matrix (A, nc=4, byrow=T)

A

A=apply (A, 1, prod)

Ai=a^ (1/4)

W=ai/sum (ai)

W# weight vectors

Lamda=sum ((A%*%w)/w)/4

CI=(lamda-4)/3

CR=CI/0.89

CR# consistency check

# step 13: the final score of each scheme is exported

Zst%*%w

Si=zst%*%w

Which.max (zst%*%w)

Si [which.max (zst%*%w) ,]

# step 14: the result that output finally needs

The optimal update side # Ren=cbind (c (3,4,5,9,13,16,19), len [which.max (zst%*%w) ,]) Case

Colnames (Ren)=c (" Road ", " Enhancement ") # and section number combine

Ren2=merge (e, Ren, by=" Road ", all=TRUE) # merges

E1=Ren2 [, 1:6] # optimal network e

D3=cda (e1, d0, ge1) # calls previously defined cda function to transport need ge1 and traffic distribution d01

#print (st [which.max (zst%*%w) ,], digits=7) # exports 4 attributes of network at this time, protects Stay 7 decimals；

FormatC (st [which.max (zst%*%w) ,], format=' f ', digits=3)

# can specify C Format Type in formatC () function with format=parameter, such as " d " (integer), " f " ' it is (fixed Point real number), " e " (scientific notation), " E ", " g " (the less output format of selection digit), " G ", " fg " (fixed point real number but use Digits specifies number of significant digit), " s " (character string) can use the specified output width of width, specify number of significant digit with digits Digit (format f) after (when format e, E, g, G, fg) or decimal point

Write.csv (d3 [[1]], file=" traffic distribution matrix .csv ") # remains to work at present mesh with csv format Record

Write.csv (d3 [[2]] [[1]], file=" section balances result .csv ")

Write.csv (d3 [[2]] [[2]], file=" OD travel time .csv ")

Write.csv (Ren2, file=" Transportation Network Design Problem scheme .csv ")

Getwd () # checks the preservation address of output file

The runing time of ### calculation procedure

Timeend <-Sys.time ()

Runningtime <-timeend-timestart

Print (runningtime) # exports runing time.

The utility model has the advantages that current Transportation Network Design Problem software only includes the dimension of one or two sustainable development, fail Sustainable development is estimated comprehensively.This software combines four important dimensions of the sustainable transports such as society, economy, environment, safety Get up, carries out Transportation Network Design Problem.This software has preferable validity and operability, for making full use of limited traffic Infrastructure investment designs sustainable transportation network, has significant effect.In addition, this software has following technical characterstic: 1) open source is free: R language is open source software, and traffic programme and management software are mostly external software, expensive；2) it is easy to grasp Make: the universal interface of external business software is complicated, using cumbersome, it is difficult to it is grasped for general engineer and planner, and this software It is easily operated；3) matrix algorithm is quickly and efficiently: external business software mostly uses C language to develop, and compares when carrying out algorithm realization It is inefficient, calculating process is greatly accelerated using matrix language.

Detailed description of the invention:

Fig. 1 is Figure of abstract.

Fig. 2 is the algorithm flow chart of underlying model

Fig. 3 is Nguyen-Dupuis test network.

Specific embodiment:

The hardware environment that this software uses is Inter 2.60GHz CPU/4G memory, software environment is Windows 8.1/ R 3.3.3/Rstudio, software are programmed using R 3.3.3.Use Nguyen-Dupuis transportation network as shown in Figure 3 Illustrate specific embodiment.Section parameter includes flowing freely travel time, road section capacity and road section length, such as table 1 It is shown

The section parameter of table 1.Nguyen-Dupuis network

There are two starting area 1 and 4 and two destinations regions 2 and 3 in Nguyen-Dupuis network.Assuming that at peak The traffic for being worth starting point 1 and 4 in the time occurs to be respectively 1800pcu/h and 1200pcu/h.That is, O₁=1800pcu/h, O₄ =1200pcu/h.The candidate road section for promoting the traffic capacity is 3,4,5,9,13,16,19.Problem is determined most according to four attributes Good road section capacity Lifting scheme, best to utilize investment budgey.Using from α=(1,1,1,1,1,1,1) Dirichlet is distributed Dir (α) and generates M=500 Investment Allocation mode at random, and wherein its dimension is the quantity of candidate road section. It the use of predetermined equation and parameter is that each allocation model calculates four sustainable attributes.

Assuming that section investment function is

g_a(Δc_a)=0.3 × Δ c_a×l_a, a ∈ A

Wherein l_aFor the road section length as unit of km.Therefore, for each realization of Dirichlet distribution, it will Have

Wherein, investment budgey is set as B=2000.

In the traffic distribution of underlying model, it is well known that do not include many in traditional Gravity Models with significant solution Release the key variables of power.Wherein, most influential is traveler destination preference.For example, traveler generally prefer that it is traditional Destination rather than area newly developed.Therefore, destination selection herein is using the multinomial logit with traveler preference Model.For simplicity, the destination preference pattern in feedback procedure is reduced to

Wherein β_jIt is preference of the traveler to destination j, β_tIt is O-D to the Path Travel Time coefficient between ij.It can make With real example data calibration β_jAnd β_tValue.β is arranged in we herein₂=0, β₃=1 and β_t=-0.1.That is, purpose is regional Traveler preference on domain 2 is 0, and the traveler preference in destination region 3 is 1, it means that traveler usually prefers mesh Ground 3.Travel time coefficient is -0.1, it means that the travel time is disutility.

In the traffic flow distribution of underlying model, traditional user equilibrium method is used herein, this method is by section performance Function is applied in equilibrium state.The common section performance function developed by Bureau of Public Roads (BPR) is as follows:

Wherein t_a(v_a, c_a+Δc_a) it is with magnitude of traffic flow v_a, existing traffic capacity c_aΔ c is promoted with the traffic capacity_aGive Determine the impedance function of section a；It is the free flow impedance of section a；α and β is can be with the traffic delay system of real example data calibration Number.For α and β, traditional BPR value is respectively 0.15 and 4.0, this is also used in our simulation study.Therefore, we can be with Road section traffic volume flow and section transit time are acquired by Frank-Wolfe algorithm.

Convergence is set as 0.01, i.e. ε=0.01 with respect to root square error (RRSE).By to each Dirichlet Allocation model use above-mentioned parameter, can restrain to obtain a stable solution, this solution have consistent travel time/cost and Traffic distribution matrix.In addition it is possible to use dijkstra's algorithm calculates the shortest path travel time between i and j.Finally, may be used To calculate the value of four sustainable attributes.Take care attribute calculating in used parameter beta₁=358.6, β₂=-407.7 And β₃=175.3.

After the performance for listing all M=500 modes, it is standardized.In order to determine weight vectors, Policymaker is invited to complete pairwise comparison matrix according to 9 grades of scales shown in table 1.Assuming that comparator matrix is

Wherein row title is followed successively by " economic attribution ", " environment attribute ", " social property " and " security attribute ".Comparator matrix Consistency ratio (CR) be 0.004, be less than threshold value 0.1, therefore comparator matrix meet consistency check and can be used for really Determine weight.Weight vectors are (0.10,0.34,0.19,0.37) w=, it means that policy maker more focus on environment attribute and Security attribute.The comprehensive score of each mode can be calculated.Top score is exactly that optimal road section capacity is promoted Mode, as shown in table 2.In this mode, TTT=116056min, TCO=34 873g, Equity=0.94, Safety= 28.16acc。

The best Transportation Network Design Problem of table 2. and section performance

Claims

1. the present invention is using a kind of multiple attribute decision making (MADM) software for sustainable transport network design of R language development, the software Algorithm frame the following steps are included:

Step 1: upper layer policy maker generates a random transportation network traffic capacity using Dirichlet distribution and is promoted Scheme Δ c, Protocol Numbers m=0, lower layer's traveler make a series of behavior reactions；

Step 2: underlying model is the sequence model that traffic generation, traffic distribution, traffic modal splitting and traffic flow distribute, and is led to It crosses feedback iteration and reaches traffic system balance, the road section traffic volume flow and transit time when equilibrium state can be calculated；

Economic attribution

Environment attribute

Social property

Security attribute；

Step 4: return step one is divided into another random transportation network traffic capacity lifting scheme using Dirichlet Δ c, Protocol Numbers m=m+1 when m=M, stop circulation, go to when scheme number reaches (M >=200) i.e. M predetermined Step 5；

2. for the step two in right 1, using following calculating process:

Step 3: user equilibrium being based on by Frank-Wolfe algorithm, traffic distribution matrix is distributed into transportation network, to calculate The magnitude of traffic flow and travel time on each section a, later, the most short travel time between starting point i and destination j, i.e.,It can To be acquired by dijkstra's algorithm；

Step 7: enabling n=n+1, user equilibrium is then based on for traffic distribution matrix by Frank-Wolfe algorithmIt distributes to Transportation network, it is later, most short between starting point i and destination j to calculate the magnitude of traffic flow and the travel time on each section a Travel time, i.e.,It can be calculated by dijkstra's algorithm, feedback arrives step 4；

Step 8: output traffic distribution matrixAnd the magnitude of traffic flow v on a of section_aTrip between starting point i and destination j Time

3. for the step five in right 1, using following calculating process:

Step 2: since the dimension of attribute is different, they must carry out standard, the category being the bigger the better for numerical value before integration Property, standardization formula is

Wherein y_ijIt is the value of the attribute j of mode i；It is the worst case of attribute j；It is the optimal cases of attribute j；z_ijIt is The standardized value of the attribute j of mode i, the attribute the smaller the better for numerical value, standardization formula are

Step 4: the comprehensive score of network design scheme i is

Step 5: after calculating M all scores, network design scheme can be ranked up, and can choose Divide highest network design, here it is optimal designs.

4. the algorithm in right 1,2,3 is write with R language, specifically it is summarized as follows:

The Multiple Attribute Decision Problems of # sustainable transport network design: bilayer model and algorithm

# step 1: initialization, by format input data and necessary packet；

#1.1 load calculates the packet of shortest path, prepares to call diikstra shortest path first, notices that igraph packet makes for the first time With needing to install, then just it is capable of calling；

#install.packages (" igraph ") # installs igraph packet

library(igraph)

Options (digits=3)

#1.2 creates the distance matrix of figure, and comprising all candidate road sections, first is classified as section label (Road), and second is classified as road Section plays piont mark (Road origin), and third is classified as road segment end label (Road destination), and the 4th is classified as the section The free flow time (free flow time), the 5th is classified as road passage capability (capacity), and the 6th is classified as link length (length), this is sentenced in traffic flow distribution for common Nguyen-Dupuis network, and detailed parameter setting can refer to program Document；

# can also be replicated in Excel, then be executed

#e=read.delim (" clipboard ", header=F)

E=matrix (c (1,1,5,7.0,900,4.00,2,1,12,9.0,700,4.00,3,4,5,9.0,700,4.00,4, 4,9,12.0,900,7.00,5,5,6,3.0,800,2.00,6,5,9,9.0,600,4.00,7,6,7,5.0,900,4.00,8, 6,10,13.0,500,8.00,9,7,8,5.0,300,4.00,10,7,11,9.0,400,5.00,11,8,2,9.0,700, 5.00,12,9,10,10.0,700,6.00,13,9,13,9.0,600,5,00,14,10,11,6.0,700,4.00,15,11, 2,9.0,700,5.00,16,11,3,8.0,700,4.00,17,12,6,7.0,300,4.00,18,12,8,14.0,700, 9.00,19,13,3,11.0,700,6.00), ncol=6, byrow=T)

#e# is used to check the breakpoint of program

#1.3 inputs initial transport need matrix d0, and first is classified as the label (OD pair) of origin and destination pair, and second is classified as starting point mark Number (origin), third are classified as terminal label (destination), and the 4th is classified as transport need (demand)；

Tge=3000# total status transport need

D0=matrix (c (1,1,2,0.2*tge, 2,1,3,0.4*tge, 3,4,2,0.3*tge, 4,4,3,0.1*tge), Ncol=4, byrow=T) # original allocation scheme

Colnames (d0)=c (" OD pair ", " Origin ", " Destination ", " Demand ")

#d0 # is used to check the breakpoint of program

The customized Frank-Wolfe algorithmic function of # pays attention to the requirement matrix d form such as d0 of input, the form of transportation network e Such as e above, relative error 0.001；

Fw=function (e, d)

{

G=add.edges (graph.empty (13), t (e [, 2:3]), and weight=e [, 4]) # creation figure, 13 be node Number, using the time as weight the length of non-path

# creates a matrix, for saving each OD pairs of shortest path and flow

V=cbind (e [, 1])

The initial each OD of st0=numeric (4) # storage is to most short running time

Colnames (V)=" Road "

V

#OD to 12 shortest path and flow

Sp12=as.vector (b12) # is converted into section label (Road)

St0 [1]=sum (e [sp12,4]) each section temporal summation of #

X12=cbind (e [sp12,1], rep (d [Isosorbide-5-Nitrae], length (sp12))) section # label and flow, in algorithm repeatedly For starting point

Colnames (x12)=c (" Road ", " V12 ")

x12

V [is.na (V)]=0

V

#OD to 13 shortest path and flow

Sp13=as.vector (b13) # is converted into section label (Road)

St0 [2]=sum (e [sp13,4]) each section temporal summation of #

X13=cbind (e [sp13,1], rep (d [2,4], length (sp13))) section # label and flow, in algorithm repeatedly For starting point

Colnames (x13)=c (" Road ", " V13 ")

x13

V [is.na (V)]=0

V

#OD to 42 shortest path and flow

Sp42=as.vector (b42) # is converted into section label (Road)

St0 [3]=sum (e [sp42,4]) each section temporal summation of #

X42=cbind (e [sp42,1], rep (d [3,4], length (sp42))) section # label and flow, in algorithm repeatedly For starting point

Colnames (x42)=c (" Road ", " V42 ")

x42

V [is.na (V)]=0

V

#OD to 43 shortest path and flow

Sp43=as.vector (b43) # is converted into section label (Road)

St0 [4]=sum (e [sp43,4]) each section temporal summation of #

X43=cbind (e [sp43,1], rep (d [4,4], length (sp43))) section # label and flow, in algorithm repeatedly For starting point

Colnames (x43)=c (" Road ", " V43 ")

x43

V [is.na (V)]=0

V

# obtains initial flow when the flow summation on all shortest paths

VS=rowSums (V [, seq (ncol (V) -3, ncol (V))])

VS

# step 2: the impedance in each section is updated

T0=e [, 4] and the # free flow time

C=e [, 5] and # road passage capability

A=0.15

B=4

Tp=function (v)

t0*(1+a*(v/c)^b)

}

repeat{

# step 3: next iteration direction is found

G2=add.edges (graph.empty (13), t (e [, 2:3]), weight=tp (VS)) # structural map, 13 be node Number, update link proportion

V=cbind (e [, 1])

Colnames (V)=" Road "

V

#OD to 12 shortest path and flow

Sp12=as.vector (b12) # is converted into section label (Road)

Colnames (x12)=c (" Road ", " V12 ")

x12

V [is.na (V)]=0

V

#OD to 13 shortest path and flow

Sp13=as.vector (b13) # is converted into section label (Road)

Colnames (x13)=c (" Road ", " V13 ")

x13

V [is.na (V)]=0

V

#OD to 42 shortest path and flow

Sp42=ad.vector (b42) # is converted into section label (Road)

Colnames (x42)=c (" Road ", " V42 ")

x42

V [is.ns (V)]=0

V

#OD to 43 shortest path and flow

Sp43=as.vector (b43) # is converted into section label (Road)

Colnames (x43)=c (" Road ", " V43 ")

x43

V [is.na (V)]=0

V

# obtains iteration direction when the flow summation on all shortest paths

VS2=rowSums (V [, seq (ncol (V) -3, ncol (V))])

VS2

# step 4: iteration step length is calculated

Step=function (lamda)

X2=VS2

X1=VS

Q=x1+lamda* (x2-x1)

sum((x2-x1)*tp(q))

}

G=function (lamda) { step (lamda) ^2 }

Lamda=optimize (g, c (0,1)) $ minimum

lamda

# step 5: new iteration starting point is determined

VS3=VS+lamda* (VS2-VS)

VS3

# step 6: test for convergence

If ((sqrt (sum ((VS3-VS) ^2))/sum (VS)) < 0.001) break

}

Result=cbind (e [, 1], round (VS, 0), tp (VS), e [, 6]/(tp (VS)/60), e [, 5], round (VS, 0)/e [, 5])

# step 7.2: each OD running time matrix u is exported；

G=add.edges (graph.empty (13), t (e [, 2:3]), and weight=result [, 3]) # creation figure, 13 be section The number of point, result are the matrix that step 7 generates

# creates a running time matrix, and for saving each OD pairs of journey time, each OD journey time of original hypothesis is 0

U=matrix (c (1,1,2,0,2,1,3,0,3,4,2,0,4,4,3,0), ncol=4, byrow=T)

#OD to 12 journey time

Sp12=as.vector (b12) # is converted into section label (Road)

#OD to 13 journey time

Sp13=as.vector (b13) # is converted into section label (Road)

U [2,4]=sum (result [sp13,3]) each section temporal summation of #

#OD to 42 journey time

Sp42=as.vector (b42) # is converted into section label (Road)

U [3,4]=sum (result [sp42,3]) each section temporal summation of #

#OD to 43 journey time

Sp43=as.vector (b43) # is converted into section label (Road)

U [4,4]=sum (result [sp43,3]) each section temporal summation of #

Accessible running time st0 between #OD couples of u=cbind (u, st0)

# exports result matrix and OD running time matrix in the form of a list

List (result, u)

}

#fw (e, d0) # is used to check the breakpoint of program

# step 8: the multinomial Logit function mlogit of destination selection is defined, is inputted as each OD running time time matrix u With various regions transport need sg, export as new traffic distribution matrix；

Mlogit=function (u, sg)

{

D=numeric (4)

D [3]=sg [2] * exp (- 0.1*u [3,4])/(exp (- 0.1*u [3,4])+exp (1-0.1*u [4,4]))

D [4]=sg [2] * exp (1-0.1*u [4,4])/(exp (- 0.1*u [3,4])+exp (1-0.1*u [4,4]))

Cbind (u [, 1:3], d)

}

# step 9: traffic distribution and traffic assignation comprehensive under a given transport need d0 and sg and transportation network e are defined and is handed over Function is balanced for iteration, the running time matrix of each OD, user when can be in the hope of user equilibrium state for initial traffic distribution Destination is reselected according to the matrix, new running time matrix can be generated again for new traffic distribution, the process one Straight circulation carries out, until traffic distribution matrix no longer changes；

Cda=function (e, d0, sg)

K=3

repeat{

D1=mlogit (fw (e, d0) [[2]], sg)

K=k+1

If (sqrt (sum ((d1 [4]-d0 [, 4]) ^2))/sum (d0 [, 4]) < 0.01) break# meets certain required precision With regard to stopping

#print (d1) # is used to check the breakpoint of program

#print (k) # is for checking program

}

#print(d1)

D2=fw (e, d1)

#print(d2)

List (d1, d2)

}

The performance of # status traffic system, the Before-After for policy compare

Ge1=c (sum (d0 [1；2,4]), sum (d0 [3:4,4])) # is used to check the breakpoint of program

# step 10:Dirichlet distribution method main program；

#install.packages (" DIRECT ") # installation generates the packet of Dirichlet, needs to install for the first time, not need again

Library (DIRECT) # load packet

Hsd=500 # indicates the number of Dirichilet sample, is compared for each, finds out optimal one

The number of candicate=7 # candidate road section

Inv=2000# investment budgey

It is the parameter for measuring safety below #

B1=358.6

B2=-407.7

B3=175.3

timestart<-Sys.time()

For (i in (1:nsd))

{

#i=2# is for checking program

It is that a new network e is constructed to each distribution model below #

# one update scheme of Ren=cbind (c (3,4,5,9,13,16,19), len [i ,])

Colnames (Ren)=c (" Road ", " Enhancement ") # and section number combine

Ren2=merge (e, Ren, by=" Road ", all=TRUE) # merges

Ren2 [, 7] [and is.na (Ren2 [, 7])]=0# replaces NA with 0, because NA participates in the result calculated or NA

E1=Ren2 [, 1:6] # one new network e1

D3=cda (e1, d0, ge1) # calls previously defined cda to fixed transport need (1200,800) and traffic distribution d0 Function

The value of 4 attributes is sought below #, and is stored in st matrix

#1. economic attribution: total travel time (minute)

St [i, 1]=d3 [[2]] [[1]] [, 2] %*%d3 [[2]] [[1]] [, 3]

#2. environment attribute: the discharge amount of CO

#3. social property: the variation before and after policy

St [i, 3]=max (d3 [[2]] [[2]] [, 4]/before [, 4])

#4. security attribute: accident number

}

# step 11: by each attribute, this 4 attributes are all the smaller the better

Zst=st

Zst [, 1]=(max (st [, 1])-st [, 1])/(max (st [, 1])-min (st [, 1]))

Zst [, 2]=(max (st [, 2])-st [, 2])/(max (st [, 2])-min (st [, 2]))

Zst [, 3]=(max (st [, 3])-st [, 3])/(max (st [, 3])-min (st [, 3]))

Zst [, 4]=(max (st [, 4])-st [, 4])/(max (st [, 4])-min (st [, 4]))

# step 12: the weight of each attribute is calculated

A=c (1,1/3,1/2,1/4,3,1,2,1,2,1/2,1,1/2,4,1,2,1)

A=matrix (A, nc=4, byrow=T)

A

A=apply (A, 1, prod)

Ai=a^ (1/4)

W=ai/sum (ai)

W# weight vectors

Lamda=sum ((A%*%w)/w)/4

CI=(lamda-4)/3

CR=CI/0.89

CR# consistency check

# step 13: the final score of each scheme is exported

Zst%*%w

Si=zst%*%w

Which.max (zst%*%w)

Si [which.max (zst%*%w) ,]

# step 14: the result that output finally needs

Ren=cbind (c (3,4,5,9,13,16,19), len [which.max (zst%*%w) ,]) optimal update scheme of #

Colnames (Ren)=c (" Road ", " Enhancement ") # and section number combine

Ren2=merge (e, Ren, by=" Road ", all=TRUE) # merges

E1=Ren2 [, 1:6] # optimal network e

#print (st [which.max (zst%*%w) ,], digits=7) # exports 4 attributes of network at this time, retains 7 Decimal；

FormatC (st [which.max (zst%*%w) ,], format=' f ', digits=3)

# can specify C Format Type in formatC () function with format=parameter, such as " d " (integer), " f " ' (fixed point is real Number), " e " (scientific notation), " E ", " g " (the less output format of selection digit), " G ", " fg " (fixed point real number but use Digits specifies number of significant digit), " s " (character string) can use the specified output width of width, specify number of significant digit with digits Digit (format f) write.csv (d3 [[1]], file=" traffic distribution after (when format e, E, g, G, fg) or decimal point Matrix .csv ") # with csv format remains to work at present catalogue

Write.csv (d3 [[2]] [[1]], file=" section balances result .csv ")

Write.csv (d3 [[2]] [[2]], file=" OD travel time .csv ")

Write.csv (Ren2, file=" Transportation Network Design Problem scheme .csv ")

Getwd () # checks the preservation address of output file

The runing time of ### calculation procedure

timeend<-Sys.time()

runningtime<-timeend-timestart

Print (runningtime) # exports runing time.