Modeling urban taxi services with multiple user classes and vehicle modes

(1)

Modeling urban taxi services with multiple user classes and vehicle modes

K.I. Wong

a

, S.C. Wong

b,*

, Hai Yang

c

, J.H. Wu

d

a

Department of Transportation Technology and Management, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 30050, Taiwan b

Department of Civil Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong c

Department of Civil Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong d_{TJKM Transportation Consultants, 5960 Inglewood, Suite 100, Pleasanton, CA 94588, USA}

a r t i c l e

i n f o

Article history:

Received 5 September 2007

Received in revised form 13 March 2008 Accepted 13 March 2008

Keywords: Network model Urban taxi services Multiple user classes Multiple vehicle modes Hierarchical modal choice

a b s t r a c t

This paper extends the model of urban taxi services in congested networks to the case of multiple user classes, multiple taxi modes, and customer hierarchical modal choice. There are several classes of customers with different values of time and money, and several modes of taxi services with distinct combinations of service area restrictions and fare lev-els. The multi-class multi-mode formulation allows the modeling of both mileage-based and congestion-based taxi fare charging mechanisms in a uniﬁed framework, and can more realistically model most urban taxi services, which are charged on the basis of both time and distance. The introduction of multiple taxi modes can also be used to model the differ-entiation between luxury taxis and normal taxis by their respective service areas and cus-tomer waiting times. We propose a simultaneous mathematical formulation of two equilibrium sub-problems for the model. One sub-problem is a combined network equilib-rium model (CNEM) that describes the hierarchical logit mode choice model of occupied taxis and normal trafﬁc, together with the vacant taxi distributions in the network. The other sub-problem is a set of linear and nonlinear equations (SLNE), which ensures the sat-isfaction of the relation between taxi and customer waiting times, the relation between customer demand and taxi supply for each taxi mode, and taxi service time constraints. The CNEM can be formulated as a variational inequality program that is solvable by means of a block Gauss–Seidel decomposition approach coupled with the method of successive averages. The SLNE can be solved by a Newtonian algorithm with a line search. The CNEM is formulated as a special case of the general travel demand model so that it is possible to incorporate the taxi model into an existing package as an add-on module, in which the algorithm for the CNEM is built in practice. Most of the parameters are observable, given that such a calibrated transport planning model exists. A numerical example is used to demonstrate the effectiveness of the proposed methodology.

1. Introduction

In most large cities the taxi industry is subject to various types of regulation, such as entry restriction and price control,

and the economic consequences of regulatory restraints have been examined in different ways (Douglas, 1972; De Vany,

1975; Shreiber, 1975; Manski and Wright, 1976; Foerster and Gilbert, 1979; Beesley and Glaister, 1983; Schroeter, 1983;

Frankena and Pautler, 1986; Cairns and Liston-Heyes, 1996; Arnott, 1996; Yang et al., 2005a,b; Loo et al., 2007). The general

objective of these studies has been to understand the manner in which demand and supply are equilibrated in the presence

of such regulations, thereby providing information for government decision making (Beesley and Glaister, 1983). It is

*Corresponding author. Tel.: +852 2859 1964; fax: +852 2559 5337. E-mail address:hhecwsc@hkucc.hku.hk(S.C. Wong).

Contents lists available atScienceDirect

Transportation Research Part B

(2)

commonly realized that two principal characteristics distinguish the taxi market from the idealized market of conventional economic analyses: the role of customer waiting time and the complex relationship between the users (customers) and sup-pliers (ﬁrms) of the taxi services.

In the taxi market, the equilibrium quantity of the service supplied (total taxi-hours) is greater than the equilibrium quantity demanded (occupied taxi-hours) by a certain amount of slack (vacant taxi-hours). This amount of slack governs the average customer waiting time. The expected customer waiting time is generally considered to be an important value or quality of the services that are received by customers. This variable affects customer decisions about whether to take a taxi, and thus plays a crucial role in the determination of the price level and the resultant equilibrium of the market (De

Vany, 1975; Abe and Brush, 1976; Manski and Wright, 1976; Foerster and Gilbert, 1979). A reduction in the expected waiting

time increases the demand for taxi services. However, from the point of view of each taxi firm, the expected customer wait-ing time is different from the quality of the typical product. In most markets in which quality is a variable, each firm decides what quality to offer. In the taxi market, the expected waiting time is usually not amenable to differentiation, but depends on the total number of vacant taxi-hours. An individual firm cannot offer customers an expected waiting time that is different

from that offered by other ﬁrms, although a large ﬁrm may be able to reduce the expected waiting time (Frankena and

Pau-tler, 1986).

All of the above studies use an aggregate approach, and do not take into account the spatial structure of the taxi market. Although taxis constitute an important transportation mode that offers a speedy, comfortable, and direct transportation ser-vice, they make considerable demands on limited road space and significantly contribute to traffic congestion, including when empty and cruising for customers. In view of the necessity and importance of the network modeling of taxi traffic,

Yang and Wong (1998)made an initial attempt to characterize taxi movements in a road network for a given customer

ori-gin–destination (O–D) demand pattern. They proposed a simultaneous system of equations to describe the movements of both empty and occupied taxis, and solved the problem with a ﬁxed-point algorithm. The model explicitly deals with the effects of taxi ﬂeet size and the degree of taxi driver uncertainty about customer demand and service conditions through various system performance measures such as taxi utilization and taxi availability at equilibrium, and thus helped to provide

information for government decision making about taxi regulations.Wong and Yang (1998)reformulated the taxi service

problem in networks as an optimization problem, which led to a more efﬁcient and convergent iterative balancing algorithm for the case without trafﬁc congestion.

Wong et al. (2001)extended the simple network model of urban taxi services proposed byYang and Wong (1998)by

incorporating congestion effects and customer demand elasticity, reformulating the problem as a simultaneous optimization

of two sub-problems, and developing a solution algorithm.Wong et al. (2002)developed a sensitivity-based solution

algo-rithm for a taxi model with congested effects and elastic demand. The potential applications of the model have been dem-onstrated by several case studies of the urban area of Hong Kong.Yang et al. (2001)conducted a case study of the calibration and validation of the simple network model for the urban area of Hong Kong.Yang et al. (2002)then investigated the nature of the demand-supply equilibrium in a regulated market for taxi service in the urban area of Hong Kong. They concentrated on the effects of alternative regulatory restraints on market equilibrium by investigating the social surplus, ﬁrm proﬁt, and

customer demand at various levels of taxi fare and ﬂeet size in regulated, competitive, and monopoly markets.Wong et al.

(2005)studied the bilateral micro-searching behavior for urban taxi services using the absorbing Markov chain approach.

Although a network equilibrium approach for urban taxi services has been developed, some important problems remain unsolved. One of the main issues is the provision of several modes of taxi services in large cities. For instance, there has been concern about the provision of accessible taxis for the handicapped, whose travel characteristics are very different from

those of other customers, and demand from afﬂuent customers for luxury taxis (Hong Kong Transport Advisory Committee,

1992, 1998). Another difﬁculty is how to model the modal similarity of different kinds of taxi services for traveler

decision-making processes. In addition, operational considerations have led to the imposition of certain service area restrictions on taxi operations. For example, in response to the rapid urban expansion in Hong Kong, nine new towns were developed in

the New Territories (rural areas) in the past three decades (Loo and Chow, 2008). To cope with the new demand in these

new towns, rural taxis were introduced to ensure service quality and taxi availability. While rural taxis are only allowed to operate in the New Territories, urban taxis can serve the entire Hong Kong (see Fig. 1 ofLoo et al. (2007)for the geograph-ical operating areas). As urban and rural taxis make different contributions to trafﬁc congestion, different government reg-ulations such as license fees and unit charges are applied. Finally, the coexistence of mileage-based and congestion-based taxi fare charging mechanisms is not uncommon in most large cities.

To tackle these difficult problems, this paper extends the single-class network model of urban taxi services to incorporate multiple user classes, multiple taxi modes, and the hierarchical modal choice of customers for taxi services. We consider taxis and normal traffic to analyze the structure of the taxi model, in which normal traffic is assumed to capture all traffic other than taxis in the network. We assume that there are several classes of customers with different values of time and money, and several modes of taxi services with distinct combinations of service area restrictions and fare levels. For service area restrictions, the taxi services are geographically divided into several modes by restricting certain types of taxis to pick-ing up or settpick-ing down customers in certain areas (i.e., they operate only within specified zones). For taxi fares, the model allows different modes of taxis to charge at different levels, but the fares are independent of the types of customers who are using the service. Both mileage-based and congestion-based charging mechanisms are considered for the general taxi fare structure. This more realistically models most urban taxi services, which are charged on the basis of both time and distance. The introduction of multiple taxi modes can also be used to model the differentiation between luxury taxis and normal taxis

(3)

by their respective service areas and customer waiting times. This extension has important implications for modeling taxi services with service area regulations such as taxi services in Hong Kong, where rural taxis are restricted to operating in rural areas, whereas urban taxis can provide service in both urban and rural areas. To model taxi traffic, it is assumed that a cus-tomer, having taken a taxi (similar to normal traffic), will try to minimize his or her individual travel cost from origin to des-tination; and a vacant taxi, having set down a customer, will try to minimize its individual expected search cost to meet the next customer. The probability that a vacant taxi meets a customer in a particular zone is specified by a logit model by assuming that the expected search time in each zone is an identically distributed random variable due to variations in the perceptions of taxi drivers and the random arrival of customers.

A simultaneous mathematical formulation of two sub-problems is proposed for the model. One sub-problem is a com-bined network equilibrium model (CNEM) that describes the hierarchical logit mode choice model of occupied taxis and nor-mal trafﬁc, together with the distribution of vacant taxis in the network that are searching for customers. Two common

approaches for this CNEM are the optimization approach (Shefﬁ, 1985; Lam and Huang, 1992; Wong et al., 2004) and the

variational inequality (VI) approach (Dafermos, 1982; Wu, 1997; Florian et al., 2002; Wu et al., 2006). When a

conges-tion-based taxi charge is neglected, the CNEM can be formulated as an optimization problem and solved by means of a

par-tial linearization approach (Wong et al., 2004). However, when both mileage-based and congestion-based fares are

considered, as in this paper, the problem can no longer be formulated as an optimization problem. Hence, we adopt the

var-iational inequality (VI) approach developed byFlorian et al. (2002)for a general multi-class multi-mode network

equilib-rium problem with a hierarchical (nested) logical structure. There are several multi-class combined models for travel

forecasting, notably those ofOppenheim (1995), Boyce and Bar-Gera (2001), and de Cea and Fernandez (2001). For a

com-prehensive review, readers may refer to a recent paper byBoyce and Bar-Gera (2004). The other sub-problem is a set of linear and nonlinear equations (SLNE), which ensures that the relation between taxi and customer waiting times, the relation be-tween customer demand and taxi supply for each taxi mode, and taxi service time constraints are satisﬁed. The VI program for the CNEM is solved by means of a block Gauss–Seidel decomposition approach coupled with the method of successive averages. The SLNE is solved by a Newtonian algorithm with a line search. On the solution algorithm, the block Gauss–Seidel decomposition algorithm is a widely used approach. We formulate the CNEM as a special case of the general travel demand model so that it is possible to incorporate our taxi model into an existing package as an add-on module, in which the algo-rithm for the CNEM is built in practice. The assumptions in the model are reasonable and can be determined in practice, and most of the parameters are observable, given that such a calibrated transport planning model exists. The additional effort in calibrating and validating the model with the inclusion of taxi ﬂows has been shown inYang et al. (2001, 2002), using the city of Hong Kong as an example. The model can be used for the assessment of the interrelationships between customer de-mand and taxi utilization for different modes of taxis that provide services in the network. A numerical example is used to demonstrate the effectiveness of the proposed methodology.

2. Model assumptions

2.1. Taxi movements in a road network

Consider a road network G(V, A) in which V is the set of nodes and A is the set of links in the network. Let I and J be the sets of customer origin and destination zones, respectively, and P and Q be the sets of user classes and taxi modes, respectively. In the following, the superscript ‘‘n” denotes normal trafﬁc, ‘‘o” denotes occupied taxis, and ‘‘v” denotes vacant taxis. We also use w 2 W = ((n, p), (o, pq), (v, q)) to represent the combination of user classes and transportation modes for p 2 P and q 2 Q.

The notation that is used throughout this paper is described as follows, and is also summarized inAppendix.

In any given hour, the total number of travelers of class p from origin zone i to destination zone j is Dp

ij, which is assumed

to be ﬁxed and known, and is expressed in vehicles per hour. The demand matrix can be split into two sub-matrices that are associated with non-taxi and taxi mode choices,

Dp ij¼ T n;p ij þ T o;p ij ; p 2 P; j 2 J; i 2 I; ð1Þ

where Tn;p_ij and To;p_ij are the numbers of non-taxi (normal trafﬁc) users and taxi users (customers), respectively. The matrix for taxi trafﬁc can be further split into several matrices for each taxi mode:

To;p ij ¼ X q2Q To;pq ij ; p 2 P; j 2 J; i 2 I; ð2Þ where To;pq

ij is the number of users who choose taxis of mode q 2 Q. A description of the mode split structure between normal

trafﬁc and taxis is shown inFig. 1, and the hierarchical mode choice is further described in Section2.6.

For each taxi mode, to meet customer demand, the occupied taxis carry customers from their origins to destinations. For each taxi mode q 2 Q, we have the following trip end constraints:

Oq i ¼ X j2J X p2P To;pq ij ; q 2 Q ; i 2 I V; ð3Þ Dq_j ¼X i2I X p2P To;pq_ij ; q 2 Q ; j 2 J V; ð4Þ

(4)

where Oq

i and D

q

j are the total customer demand from origin zone i 2 I and total customer demand to destination zone j 2 J,

respectively, for taxis of mode q 2 Q. In addition to occupied taxis, vacant taxis are also traveling in the network and search-ing for customers. The number of vacant taxis of mode q 2 Q that are travelsearch-ing from zone j to zone i and searchsearch-ing for cus-tomers is denoted as Tv;q

ji .

Trafﬁc in network Tw

ijas described above will choose paths R w ij R, w 2 W, where R ¼ [ pijR n;p ij [ pqijR o;pq ij [ qjiR v;q

ji is the set of all

routes. The corresponding path ﬂow variables are denoted as fw

r, r 2 Rwij. Let tabe the travel time on link a 2 A, and vwa be

the total vehicle ﬂow of the corresponding class-mode combination on link a 2 A. Further, let dw

ij;arbe the link route incidence

matrix, which is equal to 1 when route r between the O–D pair (i,j) of the corresponding users traverses link a and 0 other-wise. The total ﬂow on a link can then be obtained by

va¼ X p2P vn;p a þ X p2P X q2Q vo;pq a þ X q2Q vv;q a ; a 2 A; ð5Þ where vw a¼Pi2I P j2J P r2Rw ij dw_ij;arfw r.

2.2. Customer and taxi waiting times

Let Wq_i be the customer waiting time for taxis of mode q in zone i, and wq_i be the taxi waiting/searching time for taxis of mode q in zone i. The customer waiting time, which is an endogenous variable of the model, varies across zones. It is ex-pected that an additional customer waiting for a taxi of a particular mode will only increase the exex-pected waiting time for taxis of that mode, and will not be related to customers who are waiting for taxis of other modes. Thus, within each zone, the expected waiting time of customers who are waiting for taxis of mode q 2 Q can be described as a function of the density of vacant taxis of that mode in the zone Wq

i ¼ W q iðN q i;w q iÞ; i 2 I; q 2 Q , where N q

i is the number of vacant taxis of mode q per

hour that meet customers at zone i 2 I (note that Nq_i ¼ Oq_i at equilibrium). The function does not limit the value of customer waiting time. Our stationary model represents an equilibrium state over time; the modeling period has neither a beginning nor an end, and it is possible for the waiting time to be more than 1 h (the modeling period).

The speciﬁcation of the customer waiting time function depends on the distribution of taxi stands over individual zones. In the case with a continuous taxi stand distribution (taxis can pick up customers anywhere), we can assume that vacant taxis move randomly through the streets, and that for each taxi mode the expected average customer waiting time is pro-portional to the area of the zone and inversely propro-portional to the (cruising) vacant taxi hours:

Wq_i ¼ g Zi Nq iw q i ¼ g Zi Oq iw q i ; q 2 Q; i 2 I; ð6Þ

where Ziis the area of zone i 2 I and g is a model parameter that is common to all zones and all taxi modes. This approximate

distribution of waiting times can be derived theoretically (Douglas, 1972; Yang et al., 2002). Annual taxi service surveys (at sampled taxi stands and roadside observation points) have been conducted in Hong Kong since 1986 to gather information

on these customer and taxi waiting times (Transport Department, 1986–1998), and these surveys have also been used for

case studies of the validity of taxi models (Yang et al., 2001, 2002). 2.3. Cost structures

Let the generalized costs of travel on link a 2 A as perceived by different users be cw

a for the corresponding class-mode

combination, which is assumed to be a linear combination of the link travel time taand link length da. For simplicity, we

assume that the travel time ta(va) on link a 2 A is separable and an increasing function of total vehicular ﬂow vaon link Fig. 1. The hierarchical mode choice structure for normal trafﬁc and taxis.

(5)

a 2 A. Let bp₀be the value of time of users in class p (irrespective of whether they are taking taxis or normal trafﬁc), bn_{be the}

mileage costs that are charged to customers who are taking normal trafﬁc, bo;q₁ and bo;q₂ be the mileage and congestion-based costs that are charged to customers who are taking a taxi of mode q 2 Q, and bq0and bv,qbe the hourly and mileage (e.g., fuel)

operating costs of a taxi of mode q 2 Q. We have cn;p a ¼ b p 0taðvaÞ þ bnda; a 2 A; p 2 P; ð7Þ co;pq a ¼ b p 0taðvaÞ þ b o;q 1 daþ b o;q 2 taðvaÞ þ pqa; a 2 A; q 2 Q ; p 2 P; ð8Þ cv;q a ¼ b q 0taðvaÞ þ bv;qdaþ pqa; a 2 A; q 2 Q ; ð9Þ where pq

ais a cost equivalent speciﬁed charge that is applied to taxis of mode q that are entering link a 2 A. This

user-speciﬁed charge can serve several purposes. For example, in some cities, when a taxi enters a terminal (e.g., an airport) searching for customers, the taxi may have to pay a terminal charge. In the case of multiple taxi classes in which some of the classes are subject to area restrictions, this charge is a penalty to prevent the restricted mode of taxis from entering the restricted links. These charges are user speciﬁed according to the network characteristics and are zero when the taxi mode is not restricted.

We assume that each potential customer knows (from day-to-day experience) the expected waiting times for the differ-ent modes of taxi services that are available at the origin location, and can estimate the fares of differdiffer-ent taxi modes for his or her trip. The customer’s choice among various modes of taxis and normal traffic depends on the corresponding disutility val-ues, including the minimum cost of travel and other factors such as the waiting time and whether the user prefers a partic-ular mode of vehicle. The total generalized cost of travel (disutility) for users of class p who are taking normal traffic or taxis of mode q from zone i to zone j on path r can be defined as

Cn;p_r ¼X a2A dn;p_ij;arcn;p a ; r 2 R n;p ij ; p 2 P; j 2 J; i 2 I; ð10Þ Co;pqr ¼ X a2A

do;pq_ij;arco;pq a þ b p 1W q i q pq_; _{r 2 R}o;pq ij ; q 2 Q ; p 2 P; j 2 J; i 2 I; ð11Þ

respectively, where bp₁is the value of customer waiting time (for taxis) as perceived by user class p, and qpqis the inherent

inertia, which indicates that users of class p prefer taxis of mode q in relation to normal trafﬁc. The total searching cost for a vacant taxi of mode q that leaves zone j and goes to zone i on path r can be deﬁned as

Cv;q_r ¼X a2A dv;q_ji;arcv;q a ; r 2 R v;q ji ; q 2 Q; i 2 I; j 2 J: ð12Þ

With the above deﬁned path cost Cw

r, the corresponding minimum cost via the shortest route from origin i to destination j is

deﬁned as Cw ij, where C w ij¼ minðC w r;r 2 R w ij;w2 WÞ.

To facilitate our presentation, let the fare of a taxi ride be the monetary cost that is charged to a customer, which is a function of the fare level, time, and total distance that is traveled. The fare for taking a taxi of mode q 2 Q that is traveling from zone i to zone j along path r by a user of class p 2 P can be deﬁned as

Fpqr ¼

X

a2A

do;pq_ij;arFo;qa ¼

X

a2A

do;pq_ij;arðbo;q1 daþ bo;q2 taðvaÞÞ: ð13Þ

As the travel time on each path between each origin–destination (OD) pair may not necessarily be identical, we further de-ﬁne the average travel times of users who are taking various modes of transportation as hwij, where

hw ij¼ X r2Rw ij fw r X a2A dw ij;arta ! X r2Rw ij fw r; w2 W; j 2 J; i 2 I , : ð14Þ

These are the average travel times along all of the used paths between origin i and destination j, and they are weighted by the volume of ﬂows on the paths.

2.4. Taxi service time constraints

For each taxi mode q 2 Q, suppose that Nq_{cruising taxis are operating in the network, and consider one unit period (1 h) of}

taxi operations. The total occupied time of all taxis of mode q 2 Q is the taxi hours that are required to complete all To;pq ij ,

i 2 I, j 2 J, p 2 P trips, and is thus given byPp2P

P i2I P j2JT o;pq ij h o;pq

ij . In addition, the total unoccupied time consists of the moving

(searching) times of vacant taxis from zone to zone and the waiting (searching) times within zones, and is given by P j2J P i2IT v;q ji fh v;q ji þ w q

ig. The sum of total occupied taxi hours and total vacant taxi hours should be equal to the total taxi

ser-vice time. Therefore, for each taxi mode q 2 Q, the following taxi serser-vice time constraint must be satisﬁed in the (1 h) mod-eled period: X i2I X j2J X p2P To;pq ij h o;pq ij þ X j2J X i2I Tv;q ji fh v;q ji þ w q ig ¼ N q_; _{q 2 Q :} _ð15Þ

(6)

2.5. Behavior of vacant taxi drivers

The destination choice of empty taxis depends on the behavior of taxi drivers. We assume that all customers who are as-signed to a particular mode of taxi services are indistinguishable to the vacant taxis of that mode. That is, the taxi drivers of this mode will consider these potential customers to be identical when searching, and the customer search behavior and characteristics of a particular mode of vacant taxis are assumed to be independent of those for other modes of vacant taxis. We assume here that once a customer ride is completed, a taxi becomes vacant and cruises either in the same zone or goes to other zones in search of the next customer. During the search process, each taxi driver is assumed to attempt to minimize the individual expected search time/cost (including the operating cost) that is required to meet the next customer. The expected search time/cost in each zone is a random variable because of variations in the perceptions of taxi drivers and the random arrival of customers. We assume that this random variable is identically distributed, with a Gumbel density function. With these behavioral assumptions, the probability that a vacant taxi of mode q 2 Q that originated in zone j 2 J will eventually meet a customer in zone i 2 I is speciﬁed by the following logit model:

Pv;q i=j ¼ expfhq ðCv;qji þ b q 0w q iÞg P m2Iexpfh q ðCv;q_jm þ bq0w q mÞg ; q 2 Q; j 2 J; i 2 I; ð16Þ

where hq_{(1/h) is a non-negative parameter that reﬂects the degree of uncertainty in customer demand and taxi services of}

mode q in the whole market from the perspective of individual taxi drivers.

In a stationary equilibrium state, the movements of vacant taxis in the network should meet customer demand in all ori-gin zones, i.e., every customer will eventually receive taxi service after a certain waiting time. As there are Dq_j taxis of mode

q 2 Q to operate in destination zone j 2 J per hour, we have P_i2ITv;q

ji ¼ D q j; q 2 Q ; j 2 J and P j2JT v;q ji ¼ P j2JD q j Pv;q i=j ¼ O q i; q 2 Q; i 2 I.

2.6. Hierarchical logit mode choice

To determine the customer demand matrices for taxis, the modal split between normal traffic and taxis is determined in such a way that the selection of taxis is conditioned by a choice between taking a taxi and non-taxi traffic, based on the total generalized costs (disutility values) of the alternatives. We assume a hierarchical logit-type mode choice between taxis and normal traffic, where the selection of a taxi or normal traffic is in the upper choice level, and the choice of taxi mode is in the lower choice level. We also assume a logit-type mode choice function, which gives the proportion of trips taken by normal traffic Pn;p

ij and the proportion of trips taken by taxis P o;p ij : Pn;p ij ¼ expðbp 1C n;p ij Þ expðbp 1L o;p ij Þ þ expðb p 1C n;p ij Þ ; p 2 P; i 2 I; j 2 J; ð17Þ

where Po;p_ij ¼ 1 Pn;p_ij , and Lo;p_ij ¼ ð1=bp2Þ ln

P

q0_2Qexpðb

p 2C

o;pq0

ij Þ; p 2 P; j 2 J; i 2 I, is the logsum of the disutility of travel as

perceived by users of class p who are using taxis as their mode of transportation from zone i to zone j. bp 1and b

p

2are the

dis-persion coefﬁcients for the upper-level and lower-level logit mode choice model, respectively, for user-class p.

Given that taxis have been selected as the transportation mode, to determine the proportion of trips Po;pq_ij of taxis from origin zone i to destination zone j in taxis of mode q as chosen by users of class p, we introduce the following lower-level logit mode choice function:

Po;pq ij ¼ expðbp 2C o;pq ij Þ P q0_2Qexpðbp₂Co;pq 0 ij Þ ¼expðb p 2C o;pq ij Þ expðbp₂Lo;p_ij Þ; q 2 Q ; p 2 P; j 2 J; i 2 I: ð18Þ

The service area restrictions for different modes of taxis have been deﬁned implicitly in their (minimum) travel cost matrix, in which the travel cost will be very high (due to the penalty on the links based on restricted entry, as in Eqs.(8) and (9)) if the corresponding origin or destination is not accessible by that taxi mode, and the proportion of users who choose that mode of taxi will be very small (and negligible). Thus, for a given number of travelers of class p who are moving from origin zone i, to destination zone jDp

ij, the number of users of class p who are taking normal trafﬁc is T n;p ij ¼ D

p ij P

n;p

ij , and the number

of users of class p who are taking taxi services of mode q is To;pq_ij ¼ Dp_ij Po;p_ij Po;pq_ij . A description of the hierarchical mode choice

structure for normal trafﬁc and taxis is shown inFig. 1. Note that the dispersion coefﬁcient must be chosen such that

bp₂P bp₁8p 2 P for the consistency of the hierarchical logit model (Ortuzar and Willumsen, 1996). 3. Simultaneous mathematical formulation

In this section, we assume that two types of essential constraints can be temporarily relaxed: the relationship (Eq.(6))

between the customer and taxi waiting time at each origin zone (with Oq

i ¼ P j2J P p2PT o;pq ij Þ, wq iW q i X j2J X p2P To;pq ij gZi¼ 0; q 2 Q; i 2 I ð19Þ

(7)

and the conservation of flows between the zonal origin and destination subtotals of customer demands and occupied taxi movements (Eqs.(3) and (4)). Now, Wq_i, Oq_i, and Dq_j are temporarily treated as exogenous variables, i.e., Eqs.(3), (4) and (6)are not restricted to be satisfied at this stage. The internal inconsistency of these equations will be rectified later by a Newtonian algorithm for solving the SLNE.

3.1. Combined network equilibrium model (CNEM)

With the temporary relaxation of two types of essential constraints described in the previous section, we introduce the following mathematical program, which represents the combined network equilibrium model (CNEM) for multiple user clas-ses and multiple vehicle modes. The problem is formulated using the VI approach, which can handle asymmetric interactions of network ﬂows, especially for the problem of multiple user classes and vehicle modes. Note that the model presented here is a special case of that whichFlorian et al. (2002)presented. Our model considers only a road network in which customers (in the taxi or non-taxi mode) and vacant taxis are traveling. Whereas the general model ofFlorian et al. (2002)uses a ‘‘trip distribution-mode choice-assignment” hierarchy, our model uses a simultaneous ‘‘trip distribution-assignment” and ‘‘mode choice-assignment” hierarchy in which vacant taxis follow trip distribution and customers follow mode choice, and all of their movements are on the same road network. Therefore, it is a special case of the multiple user classes and vehicle modes equilibrium model in which the vacant taxi movements are not involved in mode choice and customer choices do not in-clude trip distribution.

The feasible region X of the demands and path ﬂows of the equilibrium model is as follows, where the dual variables are stated in brackets: X j2J Tv;q ji ¼ O q i; q 2 Q ; i 2 I; ða q iÞ ð20aÞ X i2I Tv;q ji ¼ D q j; q 2 Q ; j 2 J; ðb q jÞ ð20bÞ X q2Q To;pq ij ¼ T o;p ij ; p 2 P; j 2 J; i 2 I; ðL o;p ij Þ ð20cÞ To;p ij þ T n;p ij ¼ D p ij; p 2 P; j 2 J; i 2 I; ðL D;p ij Þ ð20dÞ X r2Rw ij fw r ¼ T w ij; w2 W; j 2 J; i 2 I; ðu w ijÞ ð20eÞ fw r P0; r 2 R w ij; w2 W; j 2 J; i 2 I; ðn w rÞ ð20fÞ Tw ij>0; w2 W; j 2 J; i 2 I: ð20gÞ

The constraints equations(20a) and (20b)are the conservation of ﬂows for vacant taxis; Eqs.(20c) and (20d)are the

con-servation of ﬂows for the total demand for taxis and the total demand for travel, respectively; Eq.(20e)deﬁnes that the

sum of all path ﬂows between each O–D pair should be equal to the OD matrix; and Eqs.(20f) and (20g)are the

non-neg-ativity constraints for the variables of path ﬂow and OD ﬂow, respectively. Let f ¼ ffw

r;r 2 Rwij;w2 W; j 2 J; i 2 Ig and T ¼ fT w

ij;w2 W; j 2 J; i 2 Ig. The VI program is given as follows. Find (f*, T*) 2 X

such that X i2I X j2J X w2W X r2Rw ij Cw rðf Þðfw r f w r Þ þ X p2P 1 bp₁ln T n;p ij ðT n;p ij T n;p ij Þ þ 1 bp₁ln T o;p ij ðT o;p ij T o;p ij Þ 0 B @ þ1 bp₂ X q2Q

ðlnðTo;pq_ij =To;p_ij ÞðTo;pq_ij To;pq_ij ÞÞ ! þX q2Q 1 hqln T v;q ji ðT v;q ji T v;q ji Þ ! P0 8ðf; TÞ 2 X ð21Þ where Cw r is deﬁned by Eqs.(10)–(12).

The Karush–Kuhn–Tucker (KKT) conditions of the VI program (Eq.(21)) are given in the following:

fw r : C w r u w ij n w r ¼ 0; r 2 R w ij; w2 W; j 2 J; i 2 I; ð22Þ Tn;p_ij : 1 bp₁ln T n;p ij þ u n;p ij L D;p ij ¼ 0; p 2 P; j 2 J; i 2 I; ð23Þ To;p_ij : 1 bp₁ln T o;p ij þ L o;p ij L D;p ij ¼ 0; p 2 P; j 2 J; i 2 I; ð24Þ To;pq_ij : 1 bp₂ln To;pq_ij To;p ij þ uo;pq_ij Lo;p_ij ¼ 0; q 2 Q ; p 2 P; j 2 J; i 2 I; ð25Þ Tv;q ji : 1 hq ln T v;q ji þ u v;q ji þ a q i þ b q j ¼ 0; q 2 Q ; i 2 I; j 2 J: ð26Þ

(8)

The complementarity conditions are fw rn w r ¼ 0; r 2 R w ij; w2 W; j 2 J; i 2 I; ð27Þ nw r P0; r 2 R w ij; w2 W; j 2 J; i 2 I ð28Þ

and Eqs.(20a)–(20g). At equilibrium, uw

ijis interpreted as the minimum generalized cost and C w ij¼ u

w

ij. From Eq.(25)and the fact that T o;pq ij >0, we have To;pq_ij To;p_ij ¼ expðb p 2u o;pq ij Þ expðb p 2L o;p ij Þ; q 2 Q ; p 2 P; j 2 J; i 2 I: ð29Þ

The deﬁnition of the conservation of ﬂows (Eq.(20c)) and the summation of Eq.(29)over q 2 Q implies the logsum of the

generalized cost of various taxi modes be Lo;p ij ¼ b1p 2 lnP_q0_2Qexpðbp₂uo;pq 0 ij Þ. As T o;p

ij >0, we can show that

To;pq ij To;p_ij ¼ expðbp 2u o;pq ij Þ P q0_2Qexp b p 2u o;pq0 ij ; q 2 Q ; p 2 P; j 2 J; i 2 I; ð30Þ

which satisﬁes the lower-level logit mode choice proportion of taxi trips by users of class p 2 P from origin zone i to desti-nation zone j who are taking taxis of mode q 2 Q. From Eq.(24), we have

To;p ij ¼ expðb p 1L o;p ij Þ expðb p 1L D;p ij Þ; p 2 P; j 2 J; i 2 I: ð31Þ

From Eq.(23)and the fact that Tn;p

ij >0, we have Tn;p ij ¼ expðb p 1u n;p ij Þ expðb p 1L D;p ij Þ; p 2 P; j 2 J; i 2 I: ð32Þ

Now, from the deﬁnition of the conservation of total ﬂows over the network (Eq.(20d)), we have

Tn;p ij Dp ij ¼ expðb p 1u n;p ij Þ

expðbp₁Lo;p_ij Þ þ expðbp₁un;p_ij Þ; p 2 P; j 2 J; i 2 I; ð33Þ

To;p ij Dp_ij ¼ expðbp 1L o;p ij Þ

expðbp₁Lo;p_ij Þ þ expðbp₁un;p_ij Þ; p 2 P; j 2 J; i 2 I; ð34Þ

which satisfies the upper-level logit mode choice proportion of normal traffic and taxi trips for each user class p 2 P. Com-bining Eqs.(31) and (34)gives the logsum of the generalized cost of travel in the network taking either normal traffic or taxis as LD;p ij ¼ b1p 1 lnðexpðbp 1L o;p ij Þ þ expðb p 1u n;p ij ÞÞ.

From Eq.(26)and the fact that Tv;q

ji >0, we have Tv;q ji ¼ expðh q_ðuv;q ji þ a q i þ b q jÞÞ; q 2 Q ; i 2 I; j 2 J: ð35Þ

Substituting Eq.(35)into Eq.(20b)gives rise to X i2I Tv;q ji ¼ X i2I expðhq ðuv;qji þ a q i þ b q jÞÞ ¼ D q j; q 2 Q; j 2 J: ð36Þ Thus, expðhq_bq jÞ ¼ Dq_j P i2Iexpðh q ðuv;qji þ a q iÞÞ ; q 2 Q; j 2 J: ð37Þ

Substituting Eq.(37)into Eq.(35)leads to Tv;q ji ¼ D q j expðhq ðuv;q_ji þ aq_iÞÞ P m2Iexpðh q ðuv;qjmþ a q mÞÞ ; q 2 Q ; i 2 I; j 2 J: ð38Þ

Comparing Eq.(38)with the logit-based searching behavioral assumption(16), aq

i can now be related to the taxi waiting time

in zone i for the taxis of mode q. As there is a set of dependent equations in the set of Eq.(20a) and (20b), the taxi waiting times are related to the taxi mode-speciﬁc variables by wq_i ¼ aq_i þ cq_{, where}_aq

i ¼ a q i=b q 0, and cq_¼N q Pi2I P j2J P p2PT o;pq ij h o;pq ij P j2J P i2IT v;q ji h v;q ji P i2IO q ia q i P i2IO q i ; q 2 Q; ð39Þ

which can be derived from the taxi service time constraint (Eq.(15)), are the taxi mode-speciﬁc variables (for details, see

Wong and Yang, 1998).

We showed that the solution of the VI program (Eq.(21)) satisﬁes all of the functional relationships that are required by the taxi model as deﬁned in Section2, except the set of temporarily relaxed constraints. It is noteworthy that the VI program

(9)

which case more efficient algorithms are available. The existence of the solution for the VI program is demonstrated by the fact that the VI problem admits at least one solution if its feasible region is a compact convex set and the function of the VI problem is continuous in the feasible region. In our case, the feasible region X as defined in Eq.(20)is a set of linear con-straints, positive and non-negative concon-straints, and the function defined in Eq.(21)is continuous with the region of non-neg-ative flows and positive OD matrices. Thus, we can show that there is at least one solution for the VI program. Other analytical properties of this VI program can be found inWu (1997) and Florian et al. (2002).

The ﬁnal output of the CNEM includes a complete set of link ﬂow patterns, the minimum travel cost and the correspond-ing average travel time between each OD pair, the customer demand for taxis and the distribution of vacant taxis within the

network, and the corresponding Lagrange multipliers ðaq

iÞ in the gravity model for the distribution of vacant taxis.

3.2. A set of linear and nonlinear equations (SLNE)

Solving the CNEM does not, in general, ensure the satisfaction of the relaxed constraints, Eqs.(3), (4) and (6), and the taxi service time constraint, Eq.(15). This creates internal inconsistency in the modeling and solution approach. To ensure that the desirable characteristics of the taxi model are satisfied, this inconsistency can be rectified by formulating the constraints as a set of equations, and the rectification of the set of equations is designated as the SLNE. Note that there is one dependent

equation for each taxi mode q 2 Q in the system of Eqs. (3) and (4) in view of the following fact.

P i2I P j2J P p2PT o;pq ij ¼ P i2IO q i ¼ P j2J P i2I P p2PT o;pq ij ¼ P j2JD q

j;q 2 Q. Therefore, one equation for each taxi mode can be dropped

in this set of ﬂow conservation equations, which is associated with an origin zone z 2 I or a destination zone z 2 J, where z is an arbitrarily selected zone. Assume that Dqzis an arbitrarily chosen dependent variable in the deleted constraint that

satis-ﬁes Dq z¼ P i2IO q i P j2fJzgD q j;q 2 Q .

The control variables of the SLNE are W ¼ ColðWq

i;i 2 I; q 2 QÞ, O ¼ ColðO q

i;i 2 I; q 2 Q Þ, D ¼ ColðD q

j;j 2 fJ zg; q 2 Q Þ, and

c = Col(cq, q 2 Q). The results that are obtained from the CNEM include a¼ ðaq_iðW; O; DÞ; i 2 I; q 2 Q Þ, ho¼

ðho;pqij ðW; O; DÞ; i 2 I; j 2 J; p 2 P; q 2 Q Þ, h v ¼ ðhv;qji ðW; O; DÞ; j 2 J; i 2 I; q 2 Q Þ, T o ¼ ðTo;pqij ðW; O; DÞ; i 2 I; j 2 J; p 2 P; q 2 Q Þ, and Tv_{¼ ðT}v;q ji ðW; O; DÞ; j 2 J; i 2 I; q 2 Q Þ.

Considering the relaxed constraints, Eqs.(3), (4) and (6), and the taxi service time constraint, Eq.(15), and dropping one equation for each taxi mode q 2 Q in the set of ﬂow conservation equations, the SLNE can be formulated as follows:

R1i;qðW; O; D; cÞ ¼ ðaqiðW; O; DÞ þ cqÞW q i X j2J X p2P To;pq ij ðW; O; DÞ gZi¼ 0; q 2 Q; i 2 I; ð40aÞ R2i;qðW; O; D; cÞ ¼ Oqi X j2J X p2P To;pq ij ðW; O; DÞ ¼ 0; q 2 Q; i 2 I; ð40bÞ R3j;qðW; O; D; cÞ ¼ Dqj X i2I X p2P To;pq ij ðW; O; DÞ ¼ 0; q 2 Q; j 2 fJ zg; ð40cÞ R4;qðW; O; D; cÞ ¼ X i2I X j2J X p2P To;pq ij ðW; O; DÞh o;pq ij ðW; O; DÞ þ X j2J X i2I Tv;q ji ðW; O; DÞfh v;q ji ðW; O; DÞ þ a q iðW; O; DÞ þ c q g Nq_{¼ 0; q 2 Q :} _ð40dÞ

Denote X ¼ ColðW; O; D; cÞ as the solution vector and R = Col(R1, R2, R3,R4) as the column residual vector that contains all

of the residues, where R1= Col(R1i,q,i 2 I,q 2 Q), R2= Col(R2i,q,i 2 I,q 2 Q), R3= Col(R3j,q,j 2 {J z},q 2 Q), and R4= Col(R4,q,q 2 Q).

The problem becomes one of ﬁnding a solution vector X such that the corresponding residual vector R vanishes. This is a ﬁxed-point problem within the solution space of X.

4. Solution algorithm 4.1. Algorithm for the CNEM

An extensive review of the solution algorithms for solving trafﬁc equilibrium models is given byPatriksson (2004). To

solve the VI program (Eq.(21)) for the CNEM, the block Gauss–Seidel decomposition approach coupled with the method

of successive averages (MSA) is adopted (seeFlorian et al. (2002)for details). The conditions for convergence are that the demand and link travel cost functions are strongly monotonic and that the Jacobian matrix of the cost functions with respect to the link flow of each mode is only mildly asymmetric. This is true in this multi-class multi-mode model because the link cost functions are separable for the total flow on each link. In reality, even if it is not clear whether or not this condition is satisfied, this method can be used to find an approximate solution. Readers are referred toWu et al. (2006)for a discussion of the conditions for convergence in a multi-class network equilibrium model.

Step 1. Initialization

Set the iteration number r = 0. Select an initial feasible solution T(r)and the corresponding f(r). Generally, we can select the initial solution that corresponds to the free-ﬂow network.

Step 2. Computation of generalized costs Set r = r + 1.

Compute the minimum generalized costs uw

ijand then the logsums L o;p ij and L

D;p ij .

(10)

Step 3. Computation of demand ﬂows in the combined model

Compute the demand ﬂows Tn;p

ij , T o;p ij , and T

o;pq

ij based on the generalized costs that were obtained in Step 2.

Compute the vacant taxi ﬂow Tv;q_ji by solving the doubly constrained gravity sub-model, i.e., Eq.(35)subject to Eqs.

(20a), (20b) and (20g).

Step 4. Equilibrium assignment of demand ﬂows on the network

With the demand matrices T(r)_{obtained in Step 3 above, solve the ﬁxed demand network equilibrium problem}

below to obtain the auxiliary ﬂow pattern f, such that X i2I X j2J X w2W X r2Rw_ij Cw rðfÞðfrw fwrÞ P 0 8f 2 XðTÞ

subject to Eqs.(20e) and (20f). Step 5. Convergence test

If kf fðrÞk 6 ef_{, then go to Step 7.}

Step 6. Method of successive averages

Obtain the ﬂow pattern of the next iteration by fðrþ1Þ¼ r rþ1f

ðrÞ

þ 1

rþ1f; return to Step 2.

Step 7. Computation of ﬁnal results

Compute the ﬁnal f and C, where C is a collection of the path costs.

4.2. Sensitivity analysis of the CNEM

Based on the implicit function theorem and the work ofTobin and Friesz (1988), the sensitivity analysis of the equilibrium results in the CNEM with respect to the perturbation vector (small changes in the SLNE control variables) can be formulated as a system of linear equations, in which the number of equations equals the number of variables, and the perturbed CNEM results can be obtained by solving this system of equations. The variables in the system of equations include the equilibrium path flows and the OD flows for normal traffic and (occupied and vacant) taxis, the Lagrange multipliers that are associated with the hierarchical logit modal split and vacant taxi distribution, and the minimum travel cost between each OD pair.

However, in practical applications, the number of paths is enormous for large networks. This leads to several computa-tional difﬁculties in the sensitivity analysis. First, the memory requirement for storing the arc/path incidence matrix for each OD pair of a large network is huge. Second, the sensitivity analysis for the network equilibrium problems that are discussed above requires inverting a very large matrix. The computational effort is extremely demanding when the number of paths in the network is large.

To overcome these difficulties we adopt a diagonalization approach, in which the travel time of each link in the network is assumed to be less sensitive to the change in the flow on that link when the network is perturbed from the equilibrium solu-tion, i.e., the derivative of the travel time with respect to link flow is neglected from one SLNE iteration to the next. In other words, we conduct sensitivity analysis of the VI program (Eq.(21)) in the region "T 2 X(f) (with Eq.(20e)relaxed) rather than sensitivity analysis of the whole network equilibrium problem. The path flow variables, and thus the travel cost matrix, are fixed when sensitivity analysis is conducted to determine the descent direction for the next SLNE iteration. Although without proof, this sensitivity analysis at the restricted region is very useful and sufficiently accurate for the determination of the Jacobian matrix in the SLNE. This is because in the SLNE, from computational experience, the feasible region that is governed by the set of equations is largely dominated by the OD demand T.

Assuming negligible changes in travel times on all links, the derivatives of the minimum travel costs for normal traf-fic, occupied taxi flow, and vacant taxi flow between each OD pair vanish. Thus, the derivatives of the link flow are unconstrained in the system of equations, and the derivatives of the path flow variables are no longer needed in the analysis. The derivatives of the hierarchical modal split between normal traffic and taxis and the derivatives of the va-cant taxi distribution become independent of the path flow variables, and can be solved separately. In contrast, the aver-age travel times of users in the SLNE depend on the path flow variables and the travel times along the paths. Although we make the assumption of negligible marginal travel times for all links, the average travel times between the OD pairs are not necessarily constant. It is tedious and difficult to determine the actual changes of the average travel times with respect to the perturbation parameters. Hence, we make a further approximation, that the average travel time is fixed in the iteration. This approximation is equivalent to the requirement that the change of path flows due to the change of total flows between the OD pairs is distributed among all paths on a pro-rata basis, so that the average travel times do not change.

Given the above assumptions in the diagonalization approach, the following Jacobian matrix for the SLNE can be derived by using an approach that is similar to that ofWong et al. (2002):

B ¼ B11 B12 B13 B14 B21 B22 B23 B24 B31 B32 B33 B34 B41 B42 B43 B44 2 6 6 6 4 3 7 7 7 5¼ rWR1 rOR1 rDR1 rcR1 rWR2 rOR2 rDR2 rcR2 rWR3 rOR3 rDR3 rcR3 rWR4 rOR4 rDR4 rcR4 2 6 6 6 4 3 7 7 7 5; ð41Þ

(11)

where B11;B22 are jI Qj jI Qj matrices, B33is a j{J z} Qj j{J z} Qj matrix, B44is a jQj jQj matrix, and B is a

j{I + I + J} Qj j{I + I + J} Qj matrix. The elements of the Jacobian matrix can be obtained by substituting the results of the partial sensitivity analysis of the CNEM into the direct derivative of the SLNE, with respect to the perturbation param-eters in the SLNE. Explicit expressions for the elements in the Jacobian matrix are shown inAppendix.

4.3. Solution procedure

This section presents a Newtonian solution procedure to solve the problem. At the ﬁrst iteration and starting from an ini-tial guess solution for the SLNE X(0), we compute the residual vector R(0)using the CNEM results. By applying the results of sensitivity analysis, we determine the Jacobian matrix B(0)_{. Let the current iteration be k. The improved solution can be}

ob-tained by

Xðkþ1Þ_{¼ X}ðkÞ_kðkÞ_lðkÞ_AðkÞ_RðkÞ_; _ð42Þ

where A(k)= [B(k)]1_{, l}(k)

is a control factor between 0 and 1 so that X(k) k(k)l(k)A(k)R(k)> 0, which guarantees that the up-dated control variables will not fall into the negative side, and k(k)_{is the optimal step size between 0 and 1 so that the}

resid-ual error is minimized. We adopt the golden section method for the line search. Here, the residresid-ual error is deﬁned as the weighted norm of the residual vector R(k + 1)(X(k) k(k)l(k)A(k)R(k)), where a diagonal matrix of weighting factors is adopted to normalize the different dimensions that appear in the residual vector. Denote x = Diag(x1,x2,x3,x4) as the weighting

vec-tor, where x1= Diag(1/gXi,i 2 I,q 2 Q), x2¼ Diagð1=Pj2J

P p2PD p ij;i 2 I; q 2 QÞ, x3¼ Diagð1=Pi2I P p2PD p ij;j 2 fJ zg; q 2 Q Þ, and

x4= Diag(1/Nq,q 2 Q), and E(R(k+1)) = kx R(k+1)(X(k) k(k)l(k)A(k)R(k))k as the total residual error. With the new residual

vec-tor, the improved solution can be determined using Eq. (42)for the next iteration. The solution process continues until

E(R(k + 1)_{) < e}R_{, which is an acceptable error. The procedure can be summarized as follows:}

Step 1: Set k = 0. Select an initial guess solution, X(k)_{, and compute the residual vector R}(k)_{using the CNEM results.}

Step 2: Determine the Jacobian matrix B(k)_{at the current solution X}(k)_.

Step 3: Solve for the descent direction Y(k)_{from the set of simultaneous equations B}(k)_Y(k)_{= R}(k)_{. This is more}

computation-ally efﬁcient than directly inverting the matrix B(k).

Step 4: Compute the maximum control factor l(k)_{in the range (0, 1) so that X}(k)

l(k)Y(k)_{> 0.}

Step 5: Perform a line search to determine the optimal step size k(k) _{in the range (0, 1) so that the total residual error}

E(R(k + 1)) is minimized, where R(k+1)is evaluated at X(k+1)= X(k) k(k)l(k)Y(k). Step 6: If E(R(k+1)_{) < e}R_{, then stop; otherwise, set k = k + 1 and go to Step 2.}

5. Numerical example

Consider an 8 8 square grid network with bidirectional links between each adjacent node, as shown inFig. 2. We have designed a scenario with two classes of customers and three modes of taxis. The customers are divided into a high-income group and a low-income group, where the high-income customers place a higher value on time and prefer to take luxury taxis, and the low-income customers place a lower value on time and prefer to take normal taxis. In the taxi market, there are normal taxis, luxury taxis, and restricted area taxis. Normal taxis provide the main means of transportation in the net-work and have a reasonable fare level; luxury taxis provide a better transportation service through less waiting time and more comfortable seats or extra facilities and have a higher fare level; and restricted area taxis provide services with a lower fare level, but pick up or set down customers only in certain areas. The normal taxis and luxury taxis can operate anywhere in the network, whereas the restricted area taxis can provide services only in the upper triangular portion of the network, as shown inFig. 2c. In the one-hour study period, the total number of customers that is generated in each zone is 1000 veh/h, of which 20% is in the high-income group and 80% is in the low-income group, with their destinations being evenly distributed among all of the zones in the network. To compare the differences between the multiple-class multiple-mode model and models in previous studies, we compute another scenario with a single user class and single taxi mode. In this scenario, the total number of customers generated in each zone is 1000 veh/h, and the taxis can travel anywhere in the network with-out area restrictions.

The travel impedance functions for the links are given as ta¼ t0að1 þ 0:5ðva=saÞ2Þ, where the free-ﬂow travel times of all

links, t0

a, are 0.04 h, and the capacities, sa, of all links are 3000 veh/h. The lengths of all links, da, are 3 km. The generalized link

cost functions and disutility functions take the forms that are shown in Eqs.(7)–(12), and the relationship between the cus-tomer and taxi waiting times takes the form that is shown in Eq.(6). The other input data in the example are given inTable 1. The parameters for the single-class single-mode scenario are simply taken to be the weighted average of the parameters in the multiple-class multiple-mode scenario.

To illustrate the numerical convergence characteristics of the solution algorithm at different congestion levels, all O–D demand matrices are multiplied by a ‘‘scaling factor” to generate various levels of demand, and hence congestion, in the net-work. The initial solution for the SLNE is taken as follows. For all taxi classes, the customer waiting time W is taken as 0.05 h for all zones; the taxi customer demands at all origins and destinations, O and D, respectively, are taken as the corresponding

(12)

demands that are evaluated at the free-ﬂow travel time in the network; and the slack variable c is evaluated at the corre-sponding demand ﬂows. The problem is solved by the sensitivity-based solution algorithm, and the solution is said to have converged when the residual error E(R) is less than 1%.Fig. 3shows the convergence characteristics of the problem, in which all cases at different demand levels converge quickly in a few iterations. Nevertheless, the required number of iterations in-creases with the level of congestion due to the approximation in the evaluation of the Jacobian matrix for the SLNE. Similar convergence characteristics are observed in all of the analyses that are discussed in the following paragraphs.

To compare the multiple-class multiple-mode scenario with the single-class single-mode scenario, the service level and performance of each taxi mode of the two scenarios are investigated with the total customer demand matrix being uniformly scaled up or down.Fig. 4displays the average taxi utilization versus the value of the scaling factor on demand, where there

Fig. 2. The example network: (a) an 8 8 grid network; (b) service area for normal taxis and luxury taxis; and (c) service area for restricted area taxis.

Table 1

The input parameters for the example problems

Input parameters Scenario 1. Two user classes and three taxi modes Scenario 2. Single class single mode

Values Values

Users’ value of time ðbp0;p 2 PÞ ¼ ð100; 50Þ ð$=hÞ ðb p

0;p 2 PÞ ¼ ð60Þ ð$=hÞ Value of customer waiting time for taxis ðbp1;p 2 PÞ ¼ ð200; 100Þ ð$=hÞ ðb

p

1;p 2 PÞ ¼ ð120Þ ð$=hÞ Mileage cost to a user in normal trafﬁc (bn_{) = 3 ($/km)} _(bn_{) = 3 ($/km)} Mileage charge to a taxi customer ðbo;q1 ;q 2 Q Þ ¼ ð3; 4; 2Þ ð$=kmÞ ðb

o;q

1 ;q 2 Q Þ ¼ ð3Þ ð$=kmÞ Congestion-based charge to a taxi customer ðbo;q2 ;q 2 Q Þ ¼ ð60; 80; 40Þ ð$=hÞ ðb

o;q

2 ;q 2 Q Þ ¼ ð60Þ ð$=hÞ Bias coefﬁcients for the preference to take a taxi (qpq,p = 1,q 2 Q) = (0, 40, 0) ($) and

(qpq,p = 2,q 2 Q) = (20, 0, 20) ($)

(qpq,p = P,q 2 Q) = (0) ($) [Not applicable in this scenario]

Hourly operating costs of taxis ðbq0;q 2 Q Þ ¼ ð80; 100; 80Þð$=hÞ ðb q

0;q 2 Q Þ ¼ ð85Þð$=hÞ Mileage operating costs of taxis (bv,q_{, q 2 Q) = (0.5, 0.5, 0.5) ($/km)} _(bv,q_{, q 2 Q) = (0.5) ($/km)} Dispersion coefﬁcient for the upper-level logit mode

choice

ðbp1;p 2 PÞ ¼ ð0:01; 0:03Þð1=$Þ ðb p

1;p 2 PÞ ¼ ð0:026Þð1=$Þ Dispersion coefﬁcient for the lower-level logit mode

choice

ðbp

2;p 2 PÞ ¼ ð0:02; 0:06Þð1=$Þ ðbp2;p 2 PÞ ¼ ð0:052Þð1=$Þ [Not applicable in this scenario]

Dispersion coefﬁcient for vacant taxi search behavior

(hq

,q 2 Q) = (0.2, 0.2, 0.2) (1/$) (hq

,q 2 Q) = (0.2) (1/$) Taxi ﬂeet size (Nq_{,q 2 Q) = (10,000, 5000, 5000) (veh)} _(Nq_{,q 2 Q) = (20,000) (veh)} Parameter for the relationship of customer and taxi

waiting times

(13)

are three taxi modes for scenario 1 (S1) and one taxi mode for scenario 2 (S2). The average taxi utilization is the ratio of the total occupied taxi time to the total taxi service time. For all cases, the average utilization increases with the scaling factor, and the rate of increase may be greater or less than the rate of the scaling factor because of two effects. As the scaling factor increases, the total demand for taxis increases, but traffic congestion will suppress the demand for taxis because of the con-gestion-based charges, hence the demand for taxis may increase or decrease, and the occupied taxi time for each customer ride should increase because of traffic congestion. For the first scenario, the utilization of luxury taxis and restricted area

0 1 2 3 4 5 6 Number of iterations 0 1 2 Residual error Scaling factor = 0.50 Scaling factor = 0.75 Scaling factor = 1.00 Scaling factor = 1.25 Scaling factor = 1.50

Fig. 3. The convergence characteristics of the SLNE with different scaling factors on the demand matrix.

0 20 40 60 80 100 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

Scaling factor on demand matrix

Taxi utilization (%

)

S1: Normal taxi S1: Luxury taxi S1: Restricted area taxi S2: Single Class Single Mode

(14)

taxis increases sharply with the scaling factor, but for normal taxis the utilization increases at a lower rate than the scaling factor when the scaling factor exceeds 1. This can be explained by the shift in the increased demand from normal taxis to luxury taxis, restricted area taxis, and normal traffic. Although the utilization of normal taxis increases at a lower rate than the scaling factor, it does not decrease (and is expected to approach the maximum utilization rate, or the total taxi service time), because the increasing effect on the total trip time due to increased total demand and traffic congestion should always dominate the decreasing effect because of decreased demand through traffic congestion and thus a longer customer waiting time. As discussed inWong et al. (2001), taxi utilization should always increase with traffic congestion for a wide range of parameters. For the second scenario, the case of a single taxi mode also shows increasing utilization with the scaling on cus-tomer demand, but the first scenario can provide more insightful analysis into the multi-class nature of the problem.

Figs. 5 and 6plot, respectively, the average waiting times of taxis and customers for each taxi mode versus the scaling

factor on the total customer demand. For both scenarios, the average taxi waiting times decrease sharply with the scaling factor with an increase in the total customer demand for travel. This is because when the taxi utilization and total occupied taxi hours increase, the vacant taxi hours decrease, and thus the taxi waiting times decrease. However, the average customer waiting times increase with the scaling factor, as shown inFig. 6. For the first scenario, the average customer waiting times of normal taxis and restricted area taxis increase more sharply with the scaling factor than does the average customer waiting time of luxury taxis. This is because most taxi service hours for normal taxis and restricted area taxis have already been occu-pied by customers because of their high utilization rates, and any additional occuoccu-pied taxi hours will significantly decrease taxi availability and hence increase customer waiting times. The second scenario shows the customer waiting time of the single taxi model operating as that of the normal taxis of S1 but with a fleet size of 20,000, the total number of taxis in S1. As expected, the customer waiting time is relatively short compared to the customer waiting times of the taxi modes in the first scenario, as there is no area restriction and the taxi utilization is lower than that of the normal or restricted area taxis in S1.

The rest of the section focuses on the multiple-class multiple-mode scenario. To show the effect of congestion-based charges in the fare structure of taxi services, the congestion-based charges for all of the taxi modes are uniformly scaled and the mileage costs are ﬁxed at a level, as shown inTable 1.Fig. 7portrays the average taxi utilizations for the normal, luxury, and restricted area taxis versus the ratio of the congestion-based charge to the mileage cost. The taxi utilizations for all taxi modes generally decrease as the ratio increases, and those of the normal taxis and luxury taxis decrease more rapidly than the taxi utilization of the restricted area taxis. This is because as the ratio increases, the fare levels of the taxi services increase, and thus the demand for a taxi mode generally decreases and may shift to other modes of transport (i.e., normal trafﬁc here) or other modes of taxis. In this example, because the congestion level of the restricted area (which is geographically remote from the city centre) is lower, the increase in the congestion-based charge has a lesser effect on the restricted area taxis. In addition, because the fare level of the restricted area taxis is generally lower, the increase in

0.0 0.5 1.0 1.5 2.0 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

Average taxi waiting time (hr)

(15)

the congestion-based charge has a lower impact on this taxi mode in accordance with the hierarchical mode split function. At the point at which the ratio equals zero, the congestion in the network has no effect on the taxi charges and does not affect the split in demand between taxi modes. However, when the ratio increases beyond 50 (km/h) here (although this value is too large and unreasonable in a real-world situation), the taxi utilization, which also represents the taxi market share, of re-stricted area taxis becomes larger than that of normal taxis. This shows that the level of congestion-based charges in a multi-class taxi market can also play an important role in the operational performance of taxi services.

0.00 0.05 0.10 0.15 0.20 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

Average customer waiting time (hr)

Fig. 6. Average customer waiting times against the scaling factor on the demand matrix.

0 20 40 60

Ratio of congestion-based charge to mileage-based charge (km/h) 0 20 40 60 80 100 Taxi utilization (%) Normal taxis Luxury taxis Restricted area taxis

(16)

To investigate the effects of taxi ﬂeet size and fare level on the demand for and supply of taxi services, the performance measures of the different taxi services are calculated for various numbers of taxis in service and fare levels.Figs. 8 and 9

depict the level of demands for normal traffic and taxis versus the ratio of taxi fleet sizes for luxury taxis, N2, and restricted area taxis, N3, to the taxi fleet size of normal taxis, N1.Fig. 8shows the demands for normal traffic and taxi services as the number of luxury taxis varies. The level of demand for luxury taxis initially grows with the ratio N2/N1 because luxury taxis provide better transportation service through a larger fleet size and decrease in the customer waiting time. Therefore, luxury

0.0 0.5 1.0 1.5

Ratio of luxury taxi fleet size to that of normal taxis 0 10000 20000 30000 40000 Demand (veh/h ) Normal traffic Normal taxis Luxury taxis Restricted area taxis

Fig. 8. Demands for normal traffic and taxis against the ratio of the luxury taxi fleet size to the normal taxi fleet size.

0.0 0.5 1.0 1.5

Ratio of restricted area taxi fleet size to that of normal taxis 0 10000 20000 30000 40000 Demand (veh/h ) Normal traffic Normal taxis Luxury taxis Restricted area taxis

(17)

taxis can attract more customers from other taxis or non-taxi trafﬁc. The demand for luxury taxis becomes steady as the ratio exceeds a certain value (approximately 1 here) because for this example problem, most potential customers of luxury taxis have already been exploited through the provision of better service by a large ﬂeet size, and any additional taxis that are put into service may not substantially increase the demand. A similar trend is also observed when the number of restricted area

0.0 0.5 1.0 1.5 2.0 2.5

Ratio of luxury taxi fare level to that of normal taxis 0 500 1000 1500 2000 Revenue (' 000 $) Normal taxi Luxury taxi Restricted area taxi

Fig. 10. Taxi revenues against the ratio of the luxury taxi fare level to the normal taxi fare level.

0.0 0.5 1.0 1.5 2.0 2.5

Ratio of restricted area taxi fare level to that of normal taxis

0 500 1000 1500 2000

Revenue

('

000 $

)

Normal taxi Luxury taxi Restricted area taxi

(18)

taxis varies, as shown inFig. 9. When the ratio of the restricted area taxi fleet size to the fleet size of normal taxis, N3/N1 exceeds 1, although the restricted area taxis provide better transportation service through a large fleet size with a lower fare level, the total level of demand for restricted area taxis becomes steady, and is much lower than that for normal taxis. This occurs because in the example network, restricted area taxis cannot pick up or set down customers in certain areas.

Figs. 10 and 11show the revenues of taxis of various modes against the ratio of the fare levels of luxury taxis, F2, and

restricted area taxis, F3, to the fare level of normal taxis, F1. The number of customers, the distance of each customer trip, and the fare level affect the revenues of the taxis. The revenue of luxury taxis initially increases with the corresponding fare level with the ratio F2/F1, up to approximately 0.75, and decreases thereafter, as shown inFig. 10. The reason is that for this example problem, when the fare level is at the low end, the gain in revenue that is generated by increasing the fare level is greater than the loss of revenue due to the loss of customer demand. However, when the fare increases beyond a certain level, the demand decreases at a higher rate and hence the revenue drops. This observation is characterized by the logit func-tional form that we adopt for the hierarchical mode choice model, in which the demand initially remains at a high level when the fare is low and decays exponentially when the fare exceeds a certain level (an inverse ‘‘S” shape). As restricted area taxis cannot capture the long trips outside of their service areas that have shifted from the luxury taxis, normal taxis beneﬁt the most and show a sharp increase in revenue when the revenue of luxury taxis drops. A similar trend is observed for the rev-enues of taxis against the ratio of the fare level of restricted area taxis to that of normal taxis, F3/F1, as shown inFig. 11. In this case, normal taxis still derive the greatest beneﬁt and capture the main demand shift from restricted area taxis, because of their fare competitiveness in relation to luxury taxis.

6. Conclusions

We have proposed a taxi model with multiple user classes, multiple taxi modes, and customer hierarchical modal choice in a congested road network. The contributions of this paper include the explicit consideration of multiple user classes that model user heterogeneity and multiple taxi modes that model both time-based and distance-based fare charges. The simul-taneous consideration of time-based and distance-based taxi fares at the network level is original, and can more realistically model most urban taxi services, which are charged on the basis of both time and distance. This extension has important implications for modeling taxi services with service area regulations such as taxi services in Hong Kong, where rural taxis are restricted to operating in rural areas, but urban taxis can provide service anywhere.

In the simultaneous mathematical formulation of two equilibrium sub-problems, one sub-problem is a combined net-work equilibrium model (CNEM) that is formulated as a VI program and solved by the widely used block Gauss-Seidel decomposition approach coupled with the method of successive averages, and the other sub-problem is a set of linear and nonlinear equations (SLNE) and is solved by a Newtonian solution algorithm. In the solution algorithm, the Jacobian ma-trix in the SLNE is obtained as a function of the solution from the CNEM results, and is renewed at each SLNE iteration to maintain the best quality, most likely descent direction. The proposed Newtonian solution algorithm makes the taxi model applicable to large-scale networks. We formulated the CNEM as a special case of the general travel demand model so that it is possible to incorporate our taxi model into an existing package as an add-on module, in which the algorithm for the CNEM is built in practice. We have presented a numerical example to illustrate the proposed model and algorithm. Interesting re-sults have been obtained for the example problem, and phenomena that may not hold for the abstract aggregate demand-supply model, in which all customers and taxis are assumed to have homogenous values of time and money, can be observed from our model with multiple user classes and taxi modes.

The taxi model that has been presented in this paper provides a very useful tool for the planning and evaluation of dif-ferent policy options and scenarios for urban taxi services. Its ability to model multiple user classes and multiple taxi modes makes the model applicable to a wide range of taxi problems, such as the modeling of accessible taxis for providing special services to handicapped passengers, and luxury taxis with better services and facilities for affluent customers. For example, accessible taxis, which are equipped with special facilities and equipment such as wheelchairs and/or are operated by spe-cially trained taxi drivers, can generally serve both regular customers and customers with disabilities, and disabled custom-ers will have a high preference for them (although they may still use normal taxis). It is possible that those taxis will charge at a higher fare level to ensure availability, and that handicapped passengers, in some cities, may receive various forms of transportation subsidization from the government or charity organizations. In the case of luxury taxis, more affluent custom-ers (such as those on business trips or tourists) will prefer to take a luxury taxi, which is equipped with better seats or enter-tainment facilities, rather than a normal taxi, despite the higher fare. This can be reflected in the model by adopting a higher value in the corresponding inertia parameter. The values of the inertia parameters can be obtained by the calibration of the cost functions, or user-specified according to the situation.

Acknowledgements

The work that is described in this paper was jointly supported by grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project Nos: HKU7019/99E and HKUST6212/07E), the William Mong Young Researcher Award in Engineering 2002–2003 from the University of Hong Kong, and the National Basic Research Program of China (2006CB705503).