Bi-objective evolutionary optimization of level of service in urban transportation based on traffic density

: This work investigates levels of service in urban transportation coupling a multi-objective evolutionary algorithm with the multi-agent traffic simulator MATSim. The evolutionary algorithm searches combinations of the number of private/public transportation users, capacity of buses, and time interval between bus departures minimizing traffic density and travel time simultaneously. MATSim simulates the movement of 27,000 agents according to the solutions of the evolutionary algorithm on a model of the traffic network of Quito city. We study the trade-off in objectives and analyze the solutions produced to gain knowledge about the conditions to achieve different levels of service. Also, we analyze fuel consumption and particulate matter emissions for the trade-off solutions. This work is useful for decision makers to suggest policies that can improve mobility combining private and public transportation.

ABOUT THE AUTHOR Rolando Armas has a PhD degree from the Interdisciplinary Graduate School of Science and Technology at Shinshu University, Nagano, Japan. His research interest includes sustainable systems, evolving design, evolutionary computation, traffic simulation, emission models, and data mining tools. Currently, he is a member of Frontiers in Massive Optimization and Computational Intelligence (MODO), an International Associated Laboratory between Japan and France which research includes massive optimization problems and massive optimization algorithms.

PUBLIC INTEREST STATEMENT
Approximately half of the world population currently lives in metropolitan areas, and for 2050 the projections show that will increase to twothirds. The people's mobilization demands traffic infrastructure and public transportation that does not grow at the same pace as the population, which cause problems like congestion and pollution. In this work, we simulate the movement of 27,000 agents on a model of the traffic network of Quito city to investigate levels of service in urban transportation coupling a multi-objective evolutionary algorithm with the multi-agent traffic simulator MATSim. The evolutionary algorithm searches combinations of the number of private/ public transportation users, capacity of buses, and time interval between bus departures minimizing traffic density and travel time simultaneously. We study the trade-off in objectives and analyze fuel consumption and particulate matter emissions of the trade-off solutions which assists to the decision makers for suggest policies that can improve mobility combining private and public transportation.

Introduction
Currently, around 54% of the world's population lives in urban areas. The world's population is expected to grow in the future and the proportion of the urban population as well. By 2050, the urban population is projected to increase to 66% (United Nations & Social Affairs, 2014). Population growth and urbanization trends have increased the demand for urban transportation and will continue to do so in the future. Often, existent traffic infrastructure and public transportation systems cannot cope with the increasing demand. Thus, dwellers tend to use their private means of transportation. Particularly, in emerging large cities, where the public transportation has not been able to satisfy the demand, the number of vehicles is increasing considerably and rendering mobility problems even worse (Gakenheimer, 1999). Congestion and delays add substantial costs that could be avoided, increases pollution and the risk of accidents.
The design optimization of urban transportation systems is a complex problem with multiple interrelated components that impose a continuous challenge for city planners. Evolutionary algorithms (EAs) and other heuristics have proved effective to optimize difficult problems and are being investigated in the mobility and transportation domain. Some applications of single and multi-objective EAs to road network design, bus network optimization, and public transportation routing scheduling are as follows.
In an early work (Bielli, Caramia, & Carotenuto, 2002), the authors use a single-objective genetic algorithm (GA) for bus network optimization. Hee, Kyung, and Hyuk (2012) solve a road network design problem (RNDP) using a bi-level optimization approach that reflects the different objectives between planners and network users. They focused on the design of a small network with six links and six nodes using a multi-objective GA and optimizing three objectives related to travel time, fuel consumption, and accessibility to network's nodes. Other works focus on public transportation routing and scheduling. Arbex and da Cunha (2015) use a bi-objective approach to minimize both passengers and operators costs on a benchmark network (15 nodes and 21 links). Ma et al. (2017) studied stop selection, line planning, and timetables for customized bus (CB). They used a single-objective GA to solve a weighted fitness function considering passenger, operator, and societal perspectives and test their model on real network areas that cover 3km 2 and 18km 2 . Other works focusing on frequency setting and vehicle scheduling problem employee GAs and a bi-level approach. Regarding frequency setting problem, in Yu, Yang, and Yao (2010), at the upper level, the frequency routes are optimized minimizing the passenger's travel time meanwhile in the lower level an assignation of transit trips to bus optimal route strategy is set. A single-objective GA is implemented to search the frequency of each route. The results obtained show that the optimal solution can decrease total travel time by about 6% compared with the current situation. They used a heuristic where the optimal strategy and the expected total travel time from each node to the destination node are computed, and the demand is assigned to the network according to the optimal strategy from all origins to the destination. In another study (Verbas & Mahmassani, 2015), the authors propose a bi-level model where at the upper level a frequency setting problem is addressed minimizing waiting times and use a heuristic model running a simulation at the lower level according to a demand and frequency. Both previous works used one pre-defined demand, a macro-model, and simulation schema to study the effects of a frequency setting in travel time and waiting time, respectively. Regarding scheduling problems, in Farhan Ahmad (2005), the authors state a bi-level approach to solve a bus scheduling problem for overlapping routes. In the first level, the model minimizes the number of buses required for each route individually and considering load feasibility. In the second level, given the fleet size computed in the first level, the fleet size is minimized again by considering all routes together by using a single-objective GA. They tested the model in a real-world network with a total of 60 nodes and 70 links. Link distance in time units and symmetric demand matrix gives the evaluation of each solution. The results showed that the reduction in fleet size is not significant as it was expected during model formulation stage due to little overlapping of routes in the test network.
In this work, we study the urban transportation system under various proportions of private and public transportation users, aiming to understand the conditions to achieve different levels of service and their implications on the optimal configuration of the public transportation, fuel consumption (energy), and the impact on the environment.
The level of service (LoS) is a quality measure describing operational conditions within a traffic stream, generally in terms of service measures such as speed and travel time, freedom of maneuver, traffic interruptions, and comfort and convenience (Board, 2000). In this work, we define LoS in terms of traffic density, which is determined by the total number of vehicles that circulate in the city given by the proportions of the population that use public and private transportation.
To make public transportation more attractive, city planners need to guarantee a reliable service with appropriate frequency and capacity of buses to satisfy the mobility demand. For a given proportion of public transportation users, we should find a proper configuration of bus capacities and departure times between buses. Different proportion and configurations will impact travel time, fuel consumption, costs, and emissions, in addition to traffic density. Most of these criteria are in conflict with each other. For example, the transportation service administration would like to maximize the use of the infrastructure minimizing costs and energy. Increasing the frequency of buses can increase the capacity of the system. However, this could carry additional costs. Users would like to minimize their travel times. This could be an incentive to use private transportation but could increase their costs. In addition, people living in zones where there are public transportation routes would like to minimize emissions or pollution. A way to improve the flow is reducing the volume of traffic or density. However, reducing density has impacts in the total travel time of the users because they need to use public transportation. These trade-offs suggest that compromised solutions can be found using a multi-objective optimization approach.
We use a multi-objective EA to search combinations of number of private/public transportation users, capacity of buses, and time interval between bus departures, minimizing traffic density and travel time simultaneously. In order to have a broad vision of the possible solutions to the problem and their LoS, at this time we assume that buses are always available in different capacities, do not constraint their total number, and do not consider cost.
We conduct our study using a scenario based on a real traffic network from Quito city (Ecuador). Quito's population has grown considerably in the last years. Currently, its urban transportation system is highly congested due to an increase of private vehicles combined with a poorly satisfied demand by public transportation. It is important to find alternative configurations of the urban transportation system to improve the LoS and quality of life in the city. In this work, we model the mobility of 27,000 agents that use private and public transportation on an area of 40km 2 , and simulate their mobility using MATSim (Horni, Nagel, & Axhausen, 2016), a multi-agent transport simulator.
The rest of this paper is organized as follows. In the next section, we describe the design optimization approach we follow, the problem formulation, and the simulation of mobility. In Section 3, we describe the EA, mainly its representation, operators, and fitness functions. Section 4 is devoted to simulation results and discussion. Finally, in Section 5, conclusions and future work are presented.

Components and overview
We follow the design optimization approach illustrated in Figure 1. We first formulate the biobjective optimization of LoS in urban transportation based on transit density, simulate the urban transportation in the city using a specialized transport simulator. We use a multi-objective EA to find optimal solutions to the formulated problem evaluating the quality of solutions using the outcome of the simulation, analyze the solutions produced by the EA, and extract knowledge about the system. The approach is iterative, where the solutions and obtained knowledge are made available to the expert for further analysis. The feedback from the expert could be used to reformulate the problem to study additional details. Once the expert is satisfied, the knowledge gained could be suggested to improve the real system.

Problem definition
Given a total fixed number of travelers that can use either public or private transportation, a number of bus lines and a set of schedule periods, the problem consists in finding optimal configurations of the public transportation system in terms of bus capacities and headways (time between bus departures in the same line) for a range of proportions of public transportation users minimizing simultaneously travel time and traffic density.
The problem is formulated as follows: subject to x ¼ ðn Pc ; n Pt ; c; hÞ  to all bus lines and schedule periods, C is a finite set of bus capacities, and H is a finite set of headways expressed in minutes.

Area of study
We study the conditions to achieve different levels of service in Quito city (Ecuador). The geographical area of study includes the central zone of the city, where several universities, big malls, two parks, a stadium, and the main hub of the private and public transportation infrastructure are located. The area of study also includes residential zones and the business district. The total area represents approximately 40km 2 . Figure 2 shows the extension of the geographical area.
The population in Quito was around 2 million in 2010 (INEC, 2010), of which 80% move every day to accomplish various activities. Main mobility reasons according to DMQ (2012) are 31% for work, 32.5% for study (all educational levels), and 36.5% for other different activities. It is estimated that 47% of traffic flows go to the zone we study in this work (GrupoFaro, 2010). From this percentage, 74.5% use public transportation and 25.5% private cars. In our study, we simulate a sample of 27.000 agents.
The massive passenger transportation in Quito is performed by private and public companies. The Passenger Transportation Metropolitan Company is a public company in charge of the Bus Rapid Transit (BRT) service. In 2015, approximately 650,000 passengers were transported daily by the different BRT lines (DMQ, 2015). In our model, we consider five main BRT corridors (see Figure 2) which are the most demanded and congested routes.

Simulation
In this work, we use MATSim, an agent-based transport simulation framework (Horni et al., 2016). It performs a micro-simulation of agents that move on a transportation system producing information about agent routes and movements. MATSim requires as inputs the model of the transport network infrastructure, public transportation configuration, and the agents' mobility plans.
The network comes from Geofabrik and OpenStreetMap (Frederik, Topf, & Karch, 2007). The number of links of the network corresponding to the area of study is 8192. The links' attributes are length (m), flow capacity that defines how many cars can pass through the road during a unit of time, e.g., vehicles/hour, and free speed that represents the maximum flow velocity. For this work, we take into account all the main and secondary pathways with free speeds equal or above 30 km/h. Figure 3 shows the integration between the simulator and the EA. The EA searches optimal settings for capacity and headway bus line configuration for a given proportion of agents that will use public transportation.
To evaluate a solution, MATSim sets the bus line capacities and headways with the solution provided by the optimizer, computes the routes of the agents based on heuristics, and simulates the mobility of all agents (private and public transportation) following their mobility plans. MATSim simulates the public transportation by setting stop locations, lines, schedule, routes, departures, and vehicles. A line consists of two or more routes defined over the network links and a list of departure times with a reference to a particular bus specification (e.g., capacity, length, technology) (Horni et al., 2016). MATSim emission module (Hülsmann, Gerike, Benjamin, Nagel, & Luz, 2011) computes the fuel consumption and emissions for cars and buses by connecting MATSim output to comprehensive database Handbook on Emission Factors for Road Transportation (HBEFA) (Keller & Wuthrich, 2014). The output collected from the simulator is used to calculate travel time and traffic density, which are passed back to the optimizer as objective values to compute the fitness of the solution.

Evolutionary algorithm
This section explains the representation, algorithm, operators, and fitness functions that are used in this study.

Solution representation
In this work, we consider private cars and city buses as the only means of private and public transportation, respectively. The main components of our representation are the number of travelers that use public transportation n Pt and the bus lines configurations given by capacities c and headways h. In this work n Pt can take values in the range [0.2n A , 0.8n A ] in intervals of 0.1n A , where n A is the total number of travelers or agents. The number of private transportation travelers is n Pc ¼ n A À n Pt . The n B bus lines are identified by the set of indexes L ¼ L 1 ; Á Á Á ; L n B f gand the n S schedule periods by the set of indexes S ¼ M o ; M i ; A; N f g for morning, midday, afternoon, and night. A solution for the EA is represented by the following integer variables: Here, c j i is the capacity of bus (number of passengers) of line i for schedule period j, c j i 2 C ¼ 0; 1; 2 f gwhere 0, 1, and 2 denote small, medium, and large capacity, respectively. h j i is the headway departure or time between departure of buses on line i for schedule period j. h to represent headways from 5 to 60 min in steps of 5. Figure 4 illustrates the representation of solutions used by the EA.

Adaptive ε-Sampling and v-Hood (AεSvH) algorithm
In this section, we explain the proposed AεSvH algorithm, an extension of the multi-and manyobjective (AεSεH) algorithm (Aguirre, Oyama, & Tanaka, 2013;Aguirre, Yazawa, Oyama, & Tanaka, 2014). The proposed AεSvH creates neighborhoods in variable space to bias mating for recombination and uses ε-dominance principles to truncate the population when the number of non-dominated solutions is larger than the population size.
Procedure 1 AεSvH Require: Population size P size Ensure: F 1 , set of Pareto non-dominated solutions 1: ε s 0, Δ s Δ 0 //set ε s -dominance factor and its step of adaptation 2: P random//initial population P, jPj ¼ P size 3: evaluation(P) 4: repeat 5://Parent selection 6: v-hood mating(H, Psize) 8://Offspring creation 9: Q recombination and mutation (P 0 )//jQj ¼ jPj ¼ P size 10://Evaluation and front sorting 11: evaluation(Q) 12: F non-dominated sorting ( The algorithm's flow is illustrated in Procedure 1. It starts by setting initial default values for the parameter ε s used in and its step of adaptation Δ s . The population P is randomly initialized and evaluated. The main loop starts by creating neighborhoods of solutions, grouping them by proximity in variable space using a v-hood creation procedure. Next, v-hood mating creates a pool of mates P 0 selecting them from the same neighborhood. The selected mates are recombined and mutated to create the offspring population Q. Since mates come from the same neighborhood, a solution would recombine only with a solution that is close by in decision space. After offspring Q is evaluated, nondominated sorting is performed on the population that results from joining the current population P and its offspring Q. The population of size 2P size sorted in non-dominated fronts F is then truncated to obtain the surviving population P of size P size using a ε-sampling truncation procedure set with parameter ε s . The parameter ε s and its step of adaptation Δ s are adapted every generation so that the number of ε-sampled solutions approach the size of the population. In the following we explain with more detail v-hood creation, v-hood mating and ε-sampling truncation.
The v-hood creation procedure splits the surviving population P into neighborhoods as illustrated in Procedure 2. First, this procedure randomly selects an individual z from the population P. Then, solutions in P with the same proportion of public transportation users as z N Pt are grouped in Y. Next, z and Y are assigned to neighborhood H i and removed from the population P. Neighborhood creation is repeated until all solutions in P have been assigned to a neighborhood.
Mating for recombination is implemented by the procedure v-hood mating illustrated in Procedure 3. Neighborhoods are considered to be elements of a list. To select two mates, first a neighborhood from the list is specified deterministically in a round-robin schedule. Then, two individuals are selected randomly within the specified neighborhood, so that an individual will recombine with other individual that is located close by in variable n Pt space. Due to the round-robin schedule, the next two mates will be selected from the next neighborhood in the list. When the end of the neighborhood lists is reached, mating continues with the first neighborhood in the list. Thus, all individuals have the same probability of being selected within a specified neighborhood, but due to the round-robin scheduling individuals belonging to neighborhoods with fewer members have more recombination opportunities that those belonging to neighborhoods with more members. Once the pool of all mates P 0 has been established, they are recombined and mutated according to the order they were selected during mating.
The ε-sampling truncation (Aguirre et al., 2013) receives the sets of solutions F created by non-dominated sorting and selects exactly P size surviving solutions from them. If the number of non-dominated solutions F 1 j j< P size , the fronts F i are copied iteratively to P until the population is filled. If F i overfills P, the required number of solutions are selected randomly from F i and copied to P. On the other hand, if the number of non-dominated solutions F 1 j j> P size , it calls ε-sampling with parameter ε s to get from F 1 a sample of ε s non-dominated solutions F ε 1 F 1 and the set of ε s -dominated solutions D εs . F ε 1 is copied to P. If jPj > P size , solutions are randomly eliminated from P until its size is P size . If jPj < P size , solutions are randomly copied from the ε s -dominated solutions D εs to P until its size is P size . The ε sampling procedure is done in way that the surviving population P always includes the extreme solutions present in F 1 .
Procedure 2 υ-hood creation (P) ; 2: i 1 3: j < P size 4: while j < P size do 5: fx r 1 ; x r 2 2 H i r 1^r2 ¼ r and ð 1; j j H i j Þ; r 1 Þr 2 g 6: y tournamentðx r1 ; x r 2 Þ // decide based on dominance rank 7: fx r 3 ; x r 4 2 H i r 3^r4 ¼ r and ð 1; j j H i j Þ; r 3 Þr 4 g 8: z tournamentðx r 3 ; x r 4 Þ // decide based on dominance rank 9: P 0 P 0 [ fy; zg 10: i 1 þ ði mod N H Þ 11: j j þ 1 12: end while 13: return P 0 The difference between the proposed AεSvH algorithm and the AεSεH is the space where neighborhoods for parent selection is performed. In AεSεH neighborhoods are created in objective space, whereas in AεSvH neighborhoods are created in decision space. Specifically, we use the variable n Pt to create the neighborhoods so that all solutions with the same proportion of public transportation users belong to the same neighborhood. In this way, the method allows a less-disruptive recombination compared to AεSεH for the problem dealt with in this work.

Operators
To create offspring we follow the representation described above and apply uniform crossover (UX) with probability P c , using a mixing ratio pCv. For mutation, we first decide whether to mutate n Pt with probability P n Pt m or capacities and headways (c, h) with probability 1 À P n Pt m . When capacities and headways are mutated, we traverse capacities and headways for all lines and schedule periods and mutate them with probability P m per variable in (c, h). The mutation operators increase or decrease randomly with equal probability the code of the variable in which operates.

Evaluation functions
The fitness functions f 1 ðxÞ and f 2 ðxÞ for travel time and traffic density referred in the problem definition of Section 2.2 are computed from the outcome of the simulation. In the following, we express them in terms of the data generated from the simulator.

Travel time f 1
The travel time tt is given by the total travel time of private and public transportation and expressed by The average travel time for private (cars) is expressed by where t il represents the travel time on link l for vehicle i and L is the total links of the route.
The average travel time for public transit users is computed by where t a i represents the arrival time to its destination and t d i is the departure time from the origin of the trip for an agent i. This time includes walking time to the bus stop, waiting time, and bus riding time.

Density f 2
One measure of the effectiveness to define levels of service is density. Density describes the proximity of vehicles to each other, which is the principal influence on freedom to maneuver. Thus, it is an appropriate descriptor of service quality (Roess, Prassas, & McShane, 2011). Density is a relationship given by D ¼ v/S where v is the volume or flow rate, S is the average passenger-car speed (km/h) and D is the density (pc/km/ln) (Roess et al., 2011). For our scenario, density is computed by where N TripsPc and N TripsPt are the total number of trips of private vehicles and public transportation, respectively, computed on the routes taken by agents. AvgSpeed is the total average speed of vehicles during the scenario duration time. AvgNLinks is the average of links traversed for each agent during its trip.

Experimental setup
We optimize five BRT lines, splitting the schedule of buses in four time periods from 06h00 to 09h00, 09h00 to 16h00, 16h00 to 20h00, and 20h00 to 22h00. These correspond to the schedule periods set S ¼ M o ; M i ; A; N f gfor morning, midday, afternoon, and night. For each schedule period, the algorithm searches appropriate capacities and headways (frequencies) of buses. Thus, the number of variables is 41, one for the proportion of public transportation users n Pt , 5lines Â 4scheduleperiods for capacities c and 5lines Â 4scheduleperiods for headways h. Capacity categories for BRT vehicles are small, medium and large with 25, 70, and 115 passengers, respectively. Headway is split into 12 categories from 5 to 60 min in steps of five minutes.
For the EA, we set population size to 100 and initialize it randomly. For the operators, we set crossover rate to Pc = 1.0 with mixing ratio probability pCv=0.3. This means that all individuals undergo uniform crossover and offspring keeps in average 70% of the variables of one parent and 30% of the other. Mutation rate P nPt m is set to 0.1 to mutate n Pt and 1 À P nPt m ¼ 0:9 to mutate variables in (c, h). Mutation rate per variable in (c, h) is Pm = 1/n, where n ¼ 40 is the number of variables in (c, h). Thus, in average, mutation changes the variable n Pt corresponding to the number of public transportation users for 10% of the population and the variables corresponding to capacity or headways of one of the bus lines for 90% of the population. We conduct 10 runs of the EA fixing the number of generations to 100. We perform the experiments on a 64 bits Rocks cluster with 12 nodes, of which nine are 6-core 2 GHz and 3 are 16-core 2.8 GHz. The population is evaluated in parallel and the evaluation of one individual takes in average 14 min.
To evaluate each solution, we run MATSim simulating the mobility of n A ¼ 27:000 agents. The mobility plan for each agent consists of two main trips or legs: (1) from home to work, study, or others and (2) from work, study, or others to home. The plans are designed so that all agents move first from each home location to different points along the area of study. Those points are facility locations such as universities, workplaces, and others like malls, and parks. In their second trip, the agents move back home. The distribution of home locations, workplaces, and education facilities for the mobility plan were chosen taking into account census data and a previous mobility study (Demoraes, 2005). We prepare in advance the settings of public transportation configuration: lines, routes, and stop locations. The public transportation configuration consists of five BRT lines, 112 stop locations, and two routes per each bus line. When MATSim receives the solution proposed by the EA, i.e., a vector of values corresponding to variables (n Pt , c, h), the agents that will use public transportation are chosen randomly with probability n Pt =n A . The other agents are set to use private cars. Also, the capacities and headways for the five lines are set accordingly to the values suggested by the EA.
To compute emissions, we have selected seven categories of vehicles for private and one for public transportation. Table 1 shows the distribution of vehicle categories for private transportation chosen for the scenario, which is in accordance with census transportation data (INEC, 2010) for Quito city. We assign a category of vehicle to each agent randomly following this distribution. For public transportation, we have chosen a heavy goods vehicle HGV type which weight goes from 14 to 20 tons and diesel EuroII as fuel/technology type. Vehicle characteristics enable to identify the emission factors from HBEFA database. To compute warm emissions, MATSim combines the vehicle information with the kinematic features such as speed and stop duration obtained from the simulation.

Comparison between AεSεH and AεSvH
In this section, we compare the proposed algorithm AεSvH that creates neighborhoods in decision space to perform recombination with the AεSεH that creates neighborhoods in objective space. First, we compute the hypervolume of the non-dominated set of solutions at each generation to verify their evolution. The reference point to calculate the hypervolume is fixed to the worst values in each objective, obtained from the Pareto set of non-dominated solutions found by the algorithms in all runs. Figure 5 shows box plots of the hypervolume over the generations for the 10 runs of both algorithms. In blue we show results by AεSεH and in red results by AεSvH. Note that the hypervolume rapidly increases with the generations, showing that evolution is proceeding as expected. However, the hypervolume by AεSvH is significantly larger than AεSεH. This means that AεSvH can find sets of non-dominated solutions with better properties of convergence and diversity than AεSεH. Figure 6 shows the hypervolume of a typical run by both algorithms starting them with the same random seed. Figure 7 shows the sets of non-dominated solutions found by the algorithms for the same runs. Solutions by AεSεH are shown in black, whereas solutions by AεSvH are colored according to n Pt values. Note that AεSvH can find solutions for all proportions of public transportation users n Pt , whereas AεSεH loses diversity in variable space and is not able to find solutions for large values of n Pt , particularly for 0:7 and 0:8. In addition, note that for all values of n Pt solutions by AεSvH converge to better values in objective space than solutions by AεSεH.
The two set coverage index is a metric used to compare non-dominated sets found by mutiobjective EAs and expressed in Equation (8). Given two Pareto non-dominated sets, A and B, setcoverage C(A,B) counts the number of solutions in the set B that are dominated by the solutions in the set A upon the number of solutions in the set B.
CðA; BÞ ¼ jfb 2 Bj 9a 2 A : a0 À bgj B j j (8) Table 2 shows the average set-coverage values of the Pareto non-dominated solutions found by both algorithms. Note that in average solutions by AεSvH dominate many more solutions by AεSεH, which is in accordance with the hypervolume values reported above.
In the following we focus our analysis on solutions found by AεSvH.

Trade-offs
Next, we analyze trade-offs between objectives. Unless stated otherwise, we use in our analysis the set of non-dominated solutions found in the 10 runs of the algorithm. The correlation values between objectives travel time (tt) and density (dens) are shown in Table 3 correlations suggest trade-offs between objectives. This table also includes the correlation between the objectives and fuel consumption (fc), which is measured but not used for optimization.    Figure 8 shows the trade-off between travel time and density, which is in agreement with the negative correlation shown in Table 3. This figure is colored according to the proportion of agents n Pt that use public transportation. Also, it includes vertical lines to mark ranges of the LoS computed from the traffic density. To compute LoS, we use the six levels of service defined by The Highway Capacity Manual (Board, 2000), designated by the letters A through F, with A being the highest LoS and F the lowest. Table 4 shows the different levels according to density. In the figure, the lines correspond to ranges of LoS A-D.
From this figure, it can be seen that high densities of traffic correspond to smaller proportions of public transportation users n Pt and vice versa, as expected. Regarding levels of service, note that the best LoS A requires that 70% and 80% of the agents use public transportation, LoS B can be achieved with 50% and 60% of these agents, level C with 30% and 40%, and the worst level D requires 20% when more people uses private cars that leads to the highest density of traffic.
From the same figure, we note that solutions for a given value of n Pt are clustered in density and travel time. An analysis of these clusters and their trade-off ranges provide valuable information. To illustrate, note that within each cluster there is no much difference in density, but travel time can change significantly, especially for larger fractions of public transportation users. The difference in travel time for each value of n Pt is due to the various configurations of bus headways and bus capacities. Within each cluster, solutions with lower travel time are achieved by scheduling large capacity buses with high frequency. As mentioned above, larger trade-off ranges in travel time are seen for larger n Pt , where public transportation is used more. For example, the trade-off range in average travel time is around 500 s for cluster n Pt ¼ 70% and 150 s for cluster n Pt ¼ 30%. This clearly shows ranges of improvement in travel time that can be achieved by optimizing public transportation with almost no impact to traffic density for a given value of n Pt . On the other hand, notable changes in density are observed between clusters, for different values of n Pt . This suggests that to reduce density a significant number of users should be encouraged to use public transportation instead of private cars.
In addition to travel time and density, we also estimate fuel consumption though at this time we do not use it as an objective value for optimization. Figure 9 shows travel time and fuel consumption, coloring solutions by the LoS. Note that there is a perceptible trade-off between fuel consumption and travel time, in agreement with the negative correlation shown in Table 3. Looking at the levels of service, note that in general solutions with the best LoS, such as A and B, have lower values of fuel consumption, but higher travel times. This is because high levels of service imply low traffic densities, and vice versa. High levels of service are achieved under scenarios with larger fractions of public transportation users, where, on average takes the user longer to get to their destinations, but less private vehicles circulate the city. On the other hand, solutions with worst levels of service, C and D, in general show better travel time but worse fuel consumption. Thus, the increase in fuel consumption is because more people use their cars instead of public transportation.  Armas et al., Cogent Engineering (2018), 5: 1466671 https://doi.org/10.1080/23311916.2018 Due to the high positive correlation between density and fuel consumption, shown in Table 3, a clear correspondence between Fig. 8 and 9 emerges. Looking at these two figures, note that for a given fraction n Pt of public transportation users the difference in fuel consumption can become significant although density does not change much. This can be seen clearer in Figure 10 that shows fuel consumption and density, coloring solutions according to n Pt . The strong positive correlation between density and fuel consumption suggests that the latter could be used instead of the former as objective value for optimization.

Emissions
Another important aspect of sustainable public transportation is its environmental impact. As stated above, larger fractions of public transportation users n Pt reduce overall fuel consumption in the city as shown in Figure 10. Also, best levels of service and low levels of fuel consumption can be achieved simultaneously, as shown in Figure 9.   Figure 13. Public transportation PM emission by n Pt .
In addition to fuel consumption, we collect data from MATSim's emissions module. In this work, we analyze particulate matter (PM) emissions. The concentration of PM is composed of engine exhaust gas emissions as well as abrasion and suspension. Regarding health effects, the smaller the particle, the worse the impact on humans. The effects range from damages of the respiratory system to carcinogenic effects. Figure 11 shows travel time against PM emissions of the non-dominated solutions found by the EA, coloring by n Pt . Note that there is a trade-off between travel time and PM emissions within each cluster n Pt . For a given fraction of n Pt the number of private cars is the same and the range of PM emissions observed is determined by the different configurations of headways (frequency) of buses. Thus, for a given n Pt , increasing the frequency of buses reduces travel time but increases PM. The difference in travel time is more notorious for large n Pt , where most people use public transportation. In this work we model buses with diesel as fuel, which reflects the current technology for most buses in Quito. This result suggests that a change of technology for buses, for example, from diesel to electric, is essential to achieve reasonable travel time and low traffic densities with minimum PM emissions. Figure 11 shows that the upper bound of PM increases with smaller fractions of n Pt . This is due to configurations with high frequency of buses combined with a large number of private cars. In fact, the solution with highest PM emissions for n Pt ¼ 80% includes a configuration where buses are run with slightly higher frequency than the solution with highest PM emissions for n Pt ¼ 30%. Therefore, the difference in PM emissions between these two extreme solution is mostly due to the difference in the number of private cars. High-frequency configurations of buses for low demand scenarios slightly reduce travel time and are likely to have a lower occupancy rate. These solutions appear as non-dominated because at this time we do not consider cost nor take into account occupancy rate of the buses. Figures 12 and 13 show box plots of the number of public transportation trips and PM emissions, respectively, grouped by n Pt . From these figures, we note that the median of number of trips and PM emissions increase with n Pt from 20% to 50%, and are similar for fractions 60% and above.
We analyze with more detail two solutions found by the algorithm for n Pt ¼ 80%, denoted as A8b and A8a as shown on Figure 11. Note that there is a considerable difference in PM emissions between them. Figure 14 shows separately the configurations of bus capacities (top) and headways (bottom) of all five lines and schedules for these solutions. Capacities and headways are colored according to their code. Together to the identifier of the solution we also include the value of PM emissions computed for public transportation, which are lower than those reported in Figure 11 where the value of PM emissions include private cars as well. Note that overall the frequency of buses is higher in A8b than in A8a, which is reflected in the values of PM emissions. In addition, note that to satisfy the high demand of public transportation these solutions require low capacities and relatively low frequencies for some lines and schedules, see, for example, Line3 in the afternoon schedule. This illustrates the relevance of assigning buses of different capacities and frequencies to different schedules.
The geo-located differences in PM emissions between A8b and A8a are illustrated in Figure 15, where red color indicates an increase of PM emissions in A8b and blue a reduction. Note, for example, that the levels of PM increase considerably on the network segments where bus Lines 4 and 3 are located (see Figure 2 for the location of bus lines), where A8b has higher frequency than A8a as shown in Figure 14.

Conclusions and future work
This work studied levels of service in urban transportation on a scenario of a real-world traffic network. We coupled a multi-objective EA with the multi-agent transport simulator MATSim to perform a bi-level search of the proportions of private/public transportation users and optimal  configurations of buses capacities and headways associated to each proportion. The multi-objective EA minimized travel time and density simultaneously. We defined levels of service in terms of density and analyzed it in conjunction with the trade-offs between objectives on the set of nondominated solutions found in 10 runs of the algorithm. We also analyzed fuel consumption and the effects of PM emission over the area. The results of this study provide valuable insights to understand better the conditions that can lead to different levels of service in the city. It also suggests ways to improve mobility combining private and public transportation. Our results clearly show that to improve LoS, more travelers must use public transportation, but better LoS do not necessary implies a reduction in PM because it depends on the BRT headways and the diesel engines of the buses. Results also show that the trade-off between travel time and fuel consumption agree with several configurations in capacities and headways according to the proportion of users. Also, PM emission shows that some of that settings can be environmentally friendly and user convenient. That information is worthy for decision makers due to the compromise between that objectives. In addition, results suggest that a change of technology for buses is essential to achieve reasonable travel time and low traffic densities with minimum PM emissions.
In the future, we would like to refine the model to consider constraints in resources. Also, we include other objectives such as cost, waiting time, bus occupancy rates, and evaluate the impact of a transition towards hybrid and electrical cars for private and public transportation. We would also like to extend the definition of service to cover other criteria besides density.

Funding
This work was supported by Shinshu University (Japan) and the National Secretariat of Higher Education, Science, Technology and Innovation (Senescyt), Ecuador.