Policy Research Working Paper 9526 Roads Development Optimization for All-Season Service Accessibility Improvement in Rural Nepal Using a Novel Cost-Time Model and Evolutionary Algorithm Andries Heyns Robert Banick Suraj Regmi Poverty and Equity Global Practice January 2021 Policy Research Working Paper 9526 Abstract Existing methods of prioritizing rural roads for construc- monsoon seasons—within Karnali’s infrastructure budget tion in hilly and mountainous settings require expensive constraints. Road-specific improvements in accessibility to data collection or major simplifications of ground con- services are measured by estimating accessibility changes ditions. Traditional social surplus based-methods favor resulting from each proposed road within a multimodal economic and political decision criteria over social criteria, accessibility model. In this paper, walking across Karnali’s despite evidence of the latter’s importance, and struggle mountainous, high-elevation terrain is incorporated as a to scale beyond major roads to feeder roads, forcing local primary modality—a rarity in related accessibility literature. governments with limited capacity to adopt ad-hoc alter- These improvements are implemented within heuristic and native criteria. Using roads proposed for construction in integer-linear programming optimization models. Optimi- Nepal’s remote Karnali province, this paper develops a zation-determined solutions were calculated within a day, scalable method to prioritize these roads for inclusion in and substantially outperformed the actual roads selected construction plans with the aim of optimizing potential by Karnali’s provincial government in terms of accessibility, accessibility improvements to specified services in dry and efficiency, and political economy. This paper is a product of the Poverty and Equity Global Practice. It is part of a larger effort by the World Bank to provide open access to its research and make a contribution to development policy discussions around the world. Policy Research Working Papers are also posted on the Web at http://www.worldbank.org/prwp. The authors may be contacted at rbanick@worldbank.org, sregmi@worldbank.org, and andriesheyns@gmail.com. The Policy Research Working Paper Series disseminates the findings of work in progress to encourage the exchange of ideas about development issues. An objective of the series is to get the findings out quickly, even if the presentations are less than fully polished. The papers carry the names of the authors and should be cited accordingly. The findings, interpretations, and conclusions expressed in this paper are entirely those of the authors. They do not necessarily represent the views of the International Bank for Reconstruction and Development/World Bank and its affiliated organizations, or those of the Executive Directors of the World Bank or the governments they represent. Produced by the Research Support Team Roads Development Optimization for All-Season Service Accessibility Improvement in Rural Nepal Using a Novel Cost-Time Model and Evolutionary Algorithm Andries Heynsa,b,∗, Robert Banickc , Suraj Regmic a Humanitarian Logistics and Supply Chain Research Institute, Hanken School of Economics, Helsinki, Finland b Laboratoryfor Location Science, University of Alabama, Tuscaloosa, AL, USA c World Bank, Poverty and Equity Global Practice 1. Background Due to Nepal’s lack of infrastructure and accessibility challenges brought upon by severe terrain and isolated rural populations, roadworks construction is an important element of development plans made by the Nepali government. Improved transportation links are seen as key to reducing food insecurity and unlocking the economic potential of rural areas (Thapa and Shively, 2017). As a result, Nepal’s total road length has proliferated from a base of 2,700 km in the 1970s to over 80,000km in 2015; over 7,500km of new roads were built in 2017 and 2018 alone (World Bank, 2017; Ministry of Finance, 2015). Further propelling this spree are local political economies wherein every incentive aligns towards more road building: constituents routinely cite roads as the highest priority in meetings or surveys (Sapkota, 2018; World Bank, 2016). A significant portion of recent road construction, particularly at the local level, is characterized by a lack of policy coordination, little to no study of anticipated benefits, haphazard construction standards, and dis- regard for environmental safeguards (World Bank, 2018). Accessibility improvements are unstructured and fragmented since roads are often built based on political expediency, without reference to studies or master plans (Coburn, 2020; World Bank, 2018). Furthermore, the government’s capacity is so constrained that it usually only spends 70%-80% of its capital budget and the average lifespan of roadworks projects is 12 years (World Bank, 2017). Besides perpetuating an infrastructure gap, these delays degrade the accuracy and resulting informational value of standard measures of project worth – e.g. Net Present Value (NPV) (Belli et al., 2001) or Adjusted Present Value (APV) (Myers, 1974) – and hinder the synchronization of infrastruc- ture investments with other development activities. This is unfortunate, given that Nepal’s scarce resources and harsh terrain merit meticulous, targeted approaches to maximize gains. Mechanisms for integrated plan- ning at the provincial and municipal levels are under development but their adoption is inconsistent and the methods are relatively basic (Rural Access Programme (RAP3), 2019) when compared to established ones previously implemented by the national government in collaboration with donors (World Bank, 2019). Nepal’s northwestern Karnali province (also known as Province 6) exhibits extreme levels of inaccessi- bility even by Nepal’s standards. With a predominantly rural population, access to critical services is often measured in hours, walking is a primary transport modality, and it contains the only two districts (of 77) whose headquarters are not connected by vehicular roads of any sort, Humla and Dolpa (Devkota et al., 2012). An illustration of the magnitude of Karnali’s accessibility challenges is provided in Figure 1, in which the provincial population share with long average travel times to financial institutions is markedly higher in Karnali than all other provinces (Banick and Kawasoe, 2019). Like much of Nepal, damages from rains and landslides further constrain access to services during Nepal’s June-September monsoon season, when heavy rains intensely damage almost all of Karnali’s poorly main- tained roads (Huber, 2015). In Figure 2, examples of typical monsoon road conditions in Karnali are provided, showing their vulnerability to bad weather. Significant reductions in travel speeds (40%-80%) for affected road segments are common during this season, particularly for poorly engineered roads with inadequate drainage, and/or poor surfaces (Pantha et al., 2010). Since terrain and available budget effectively prohibit redundancy in roads, damage to a single road can severely impact entire regional populations. ∗ Corresponding author, contactable at andries.heyns@hanken.fi Figure 1: Comparing access to financial institutions in all provinces of Nepal. Each dot represents the share of one municipality’s population in that travel time category (Banick and Kawasoe, 2019). (a) (b) Figure 2: An example of typical monsoon road conditions in Nepal’s Karnali Province, which are vulnerable to rain and landslides. Images from Chapagain and Limbu (2018) and provided by Nepal’s Rural Access Programme 3. 2 Because Karnali’s severe terrain impedes accessibility, many remote areas are poorly covered by existing roads and services (see Figure 3(a)), while branch roads usually only reach small populations and are costlier to construct in remote, mountainous areas in its sparsely populated north and northwest (such as the Dolpo region of Karnali in Figure 3(b)). Quantifying this remoteness has always challenged analysts; linear distances and routes are inaccurate and misleading in a mountainous environment with substantial ups, downs, and detours around impassable features. Karnali’s Multidimensional Poverty Index (MPI) headcount ratio is the worst of any province in Nepal and development is a policy imperative (National Planning Commission and Central Bureau of Statistics, 2020).1 To deliver rapid economic growth and alleviate poverty government officials have expressed interest in adopting more data driven approaches to planning development projects – especially road building projects, given the aforementioned infrastructure gap and high citizen demand. To support their efforts we investigated a more data driven approach for prioritizing rural roads construction and strategically aligning it with other key components of development. This paper describes that approach, which accommodates the basic realities of budget, data availability, and institutional capacity in Nepal. The novel contributions resulting from the work produced in this paper are primarily a) a workflow to generate and compare multi-modal accessibility models for planned infrastructure in rural, mountainous contexts, and b) an optimization approach that takes the resulting accessibility data as inputs and returns numerous possible roads construction alternatives with optimal characteristics across a budget range, as opposed to traditional optimization methods which typically return a single solution with minimal comparative insight. We present the work as a complement or easily scalable alternative to traditional roads investment decision tools in rural, mountainous settings where data collection may be prohibitively costly. By applying our methods, we are able to determine numerous candidate roads development plans within a large budget range, selected from 92 candidate roads for inclusion in Karnali’s 5-year Provincial Transport Master Plan, a final draft of which was shared by Nepal’s Rural Access Programme 3 (RAP3). When compared to the actual roads selected by Karnali’s government, our solutions outperformed their plan by 1.5 to 3 times when evaluated with respect to accessibility improvements to health care, financial services, and government/market centers at equivalent cost levels, or yielded 1 to 2 times the improvements at half the cost. The rest of this paper is structured as follows. In Section 2 we present an overview of relevant literature from several disciplines considered in the development of our approach, including development and roads in rural areas, geospatial accessibility models, optimization approaches investigated for roads planning, and rural roads planning specifically in Nepal. In Section 3 we reveal the input data that were used, describe our geospatial analysis methodology, and elucidate our optimization methodology. In Section 4 we present our optimization results and compare the quality of the solutions to the roads development plan recently selected for Karnali’s Provincial Master Transport Plan. In Section 5 we elaborate further on the significance of the results, discuss the applicability of our methods within Nepal, identify limitations, and suggest future improvements and research directions. We close with a brief conclusion in Section 6. 2. Literature 2.1. Development, roads and accessibility in rural areas Close links between rural development and access to markets, services, and employment centers are a common theme of development literature (Edmonds, 1998; Booth et al., 2000; Porter, 1997, 2002a,b; Bird et al., 2002; Christiaensen et al., 2003; Stifel et al., 2003). As physical access is primarily dependent on transportation infrastructure, traditionally focus has fallen on quantifying the economic impacts of rural roads investment (Jalan and Ravallion, 2002). Decades of academic research indicate that the realized benefits of rural roads – like improvements in household consumption, school attendance, food and non- food item market prices, and non-farm employment – vary considerably based on contextual characteristics. These include agronomic potential (Binswanger et al., 1993), initial level of development (Khandker et al., 1 Multidimensional poverty is a composite measure of poverty based on multiple potential deprivations, including non- economic measures like access to health care, education, or overall standard of living. More information can be found at http://hdr.undp.org/en/content/what-multidimensional-poverty-index. 3 (a) (b) Figure 3: (a) Karnali’s population, roads, and services cluster in the lower altitude south, and (b) roads vs. population density in the remote Dolpo region of northeastern Karnali. 4 2009; Bryceson et al., 2008), local market characteristics (Casaburi et al., 2013), and lack of agglomeration economies and/or human capital (Asher et al., 2018). Given this contextual uncertainty, evidence exists that the most consistent living standards improvements from rural roads are not economic in nature but instead linked to improved access to services such as health (Van de Walle and Cratty, 2002), finance (Binswanger et al., 1993), education (Adukia et al., 2020; Kristjanson et al., 2005), or administration. Comparing Asian Development Bank projects in the Philippines, Indonesia, and Sri Lanka, Hettige (2006) finds that rural roads most consistently benefit the poor through improved access to services, whereas poverty reduction is conditional on the contextual factors listed above or complementary investments. The relationship between rural accessibility and human development is pronounced in South Asia com- pared to global norms. Whereas 29% and 11% of rural and urban populations respectively live in multidimen- sional poverty globally, in South Asia 64% and 25% do because of the excessive concentration of economic opportunity in urban areas (Asher et al., 2018; United Nations Development Programme, 2017). In India, increasing distance from towns and their services correlates with reduced household earnings (Asher et al., 2018) and improving rural road connectivity increases rural household consumption and reduces rural poverty (Rao and Jayasree, 2003; Khandker and Koolwal, 2011). In Nepal, Jacoby (2000) and Fafchamps and Shilpi (2003) demonstrated that income, cropping intensity, access to non-farm employment, and education levels in rural areas progressively fall as travel times to markets increase. In a small-scale study of 3 Nepali village development committees, Sapkota (2018) found that households with longer travel times to services have higher poverty levels and lower levels of health, education, and reported happiness. Shrestha (2012) found that a 1% increase in travel times to paved roads decreased rural farmland values and implied profits by 0.25%. 2.2. Measures of accessibility Measuring accessibility baselines and improvements is critical to establishing rural roads’ benefits. Simple metrics like linear distance, densities of roads and services, or reported travel times in household surveys are easy to compute or parse from existing data sources but are often difficult to scale, only tangentially related to access, influenced by idiosyncratic survey responses (Ahlstr¨ om et al., 2011), or inaccurate in mountainous areas (Huber, 2015). Nevertheless, these are commonly used in economists’ evaluations of rural roads: for instance, Banerjee et al. (2012) and Thapa and Shively (2017) measure accessibility in terms of road length density per administrative area in China and Nepal, respectively. Planners and researchers looking for greater accuracy and scalability will typically employ geospatial models of accessibility. The foundational models employed are network and cost-time models, which may be further elaborated in gravity models or measures of geographic and potential accessibility (Rodrigue, 2016). Network analysis calculates the anticipated travel time between one or many origins and destinations over an integrated transportation network. Transportation networks can take the form of roads, paths, railways, waterways, flights, or combinations thereof. Speeds are assigned to segments of this network which are used to calculate travel times between origin-destination pairs, with adjustments for impedances like traffic or border crossings. Implementations of network analysis range from the study of P´ aez (2004) on network accessibility’s relationship to economic activity in East Asia, the study of intra-city accessibility in Colombia by Sabogal et al. (2018), World Bank studies of accessibility to jobs in Cox’s Bazar Bangladesh (World Bank, 2020), and even Google Maps’s Directions tool. Geospatial cost-time models represent travel in terms of the minimum time to cross each cell in a grid (raster). These times are calculated based on some combination of the underlying roads and road speeds, landcover, walking speeds, terrain, and other factors (Zhu et al., 2016). The advantage of this approach is that grid cells can model travel as on-network, off-network (e.g. walking), or as a combination of the two (walk to the road and then take a bus). Relevant examples leveraging these off-road capabilities to analyze access in difficult contexts include the analysis of Kosmidou-Bradley and Blankespoor (2019) on accessibility to services in rural Afghanistan, the study by Ahlstr¨ om et al. (2011) on poverty’s relationship to accessibility in southern Sri Lanka, models of accessibility to markets in Africa by Harvest Choice and International Food Policy Research Institute (IFPRI) (2016), and global-scale analysis of wealth stratification by accessibility to urban centers by Weiss et al. (2018). The accuracy with which models reflect seasonal access changes depends on the model and metrics em- ployed. Transport geographers (Jafino et al., 2020; Rodrigue et al., 2016; Jenelius, 2009) and civil engineers 5 (Shrestha, 2018) have developed concepts of criticality and redundancy to explain the impact of seasonal disruptions on road network performance within a network analysis. Because these network-scale analyses approach seasonal impacts in terms of binary open-closed changes in roads’ status they obscure quantita- tive changes in network and road-scale performance from reduced speeds over damaged but still-functioning roads. Many are also data intensive and difficult to replicate in rural, developing world contexts. Where seasonal access is considered by social sciences researchers, they tend to rely on survey data, as in the model of seasonal market access in Ethiopia by Hirvonen et al. (2017). In Nepal, Huber (2015) adjusts speeds for major roads and airports for the monsoon based on expert guidance. Walking travel must be modeled in rural areas of the developing world where lack of roads or simple poverty compel many to travel long distances on foot to reach services or transport networks. They can therefore represent the entire travel time to a service or be integrated with road transport for a multi-modal accessibility model. Representations of walking speeds vary in sophistication: while the hiking functions of Tobler (1993) and Van Wagtendonk and Benedict (1980) provide foundations for modeling walking speeds, they translate better to cost-time modeling, which treat walking cells like any other, than network analysis, which ignores off-network travel. In practice planners and researchers tend to ignore or simplify walking times in their analysis of rural roads given roads’ primacy in most economic and transport systems. In Nepal, Devkota et al. (2012) modeled walking as a primary transport modality when recommending rural trail bridge investments. To the authors’ knowledge however, the model we adapt here (from Banick and Kawasoe (2019)) remains the only published, scalable model of travel times which incorporates monsoon impacts and walking speeds in a rugged and sometimes high-altitude context. 2.3. Traditional rural transport valuation Since available budgets for rural road building rarely match the observed needs, transportation planners and policy makers must choose the most productive subset of roads from a collection of potential investments. Valuation and prioritization exercises for rural roads investments traditionally define such productivity ac- cording to economic cost-benefit analysis (CBA). Van der Tak and Anandarup (1971) framed the classic method of evaluating rural roads construction benefits in terms of welfare gains within a CBA, as measured by society’s net savings from reductions in the cost of supplying goods, increased employment, and increases in commercial traffic linked to new roads. Improvements in travel times are principally used to measure changes in vehicle operating costs and opportunity costs of lost labor time. Lebo and Schelling (2001) up- dated this method to better emphasize year-round access to rural areas by incorporating seasonal access cost variations, the monetary value of improved access to services, and improved valuation of time savings into traditional cost-benefit analysis. As noted by Belli et al. (2001), the large budgets and timeframes required for these data-intensive methods make them difficult to scale to large sets of roads or road alternatives. Accordingly, Lebo and Schelling argued for analyzing roads’ cost-effectiveness – relative to stated develop- ment objectives, instead of economic outcomes – when evaluating lightly trafficked but locally critical rural transport infrastructure where comprehensive CBA may be prohibitively costly. Accessibility indicators based on travel times may complement or outperform CBAs in rural areas as they present easily parsable decision criteria to lightly resourced rural governments, scale easily across large sets of roads, and readily translate to monetary opportunity costs of travel. Additionally, the travel time indicators generated by cost-time models have demonstrated correlations with changes in fields as varied as conservation (Nelson and Chomitz, 2011), poverty reduction (Porter, 2007; Ahlstr¨ om et al., 2011), trade (Schmitz et al., 2012), and disease prevention (Gilbert et al., 2014). For these reasons and the aforementioned difficulty of appropriately resourcing CBAs for many roads over large spatial scales, we evaluate roads according to their cost-effectiveness versus accessibility indicators in this paper. 2.4. Optimization approaches to (rural) transport planning Previous sections detailed examples of criteria considered when planning roads development, and the manners in which the benefits of proposed projects may be quantified – specifically in a rural developing- world context, with a focus on accessibility. It is only natural that more than one potential approach (e.g. candidate development plan) for the upgrading of a roads network exist (the work of Scaparra and Church (2005) provides excellent examples of this), from which it is necessary to identify at least one which provides superior trade-offs (in cost, accessibility improvement, etc.) compared to the others. This section discusses 6 the optimization approaches which exist to identify these superior alternatives and how they have been employed in works concerning rural roads Optimizing across rural networks presents a unique and understudied challenge since traditional trans- portation optimization models are designed for urban, relatively developed settings and hence are likely to underperform in providing solutions for rural, sparsely populated settings (Shrestha, 2018). However, the optimization of rural roads development has not been entirely overlooked and has seen recent advances, building upon earlier foundational work in this context by Kumar and Tillotson (1991), Kumar and Kumar (1999), and Kumar and Jain (2000) – noteworthy is that these studies are all focused in South Asia. The demand for transportation arises from the spatial mismatch between where people are and the locations of services they want to reach, so it is common for transportation network design studies to include the service locations in their evaluation process. Kumar and Tillotson (1991) investigated the use of a mathematical model that is capable of determining numerous high-quality road network layouts that provide trade-offs in cost and travel times to fixed service/facility locations – market centers in their case. They investigated a rural area in India, and the travel time to market centers (which serves as a proxy for other essential services) was investigated as the trade-off criterion against cost. Their objective was to provide all- weather access to every village in their study area, achieved by minimizing the cumulative costs of upgrades to existing roads and the construction of new ones, while minimizing the sum of travel costs required from each village to reach a selected service/facility location. This work was extended by Kumar and Kumar (1999) and Kumar and Jain (2000), in which the same model was extended to be capable of determining trade-offs for three accessibility criteria (educational institutions in this instance, but possible to consider others) measured against cost. Scaparra and Church (2005) introduced a model that is capable of determining candidate road network development approaches in a rural setting, with consideration given to the improvement of all-season con- nectivity between villages and the efficiency of the routes between them – given a fixed budget. Their work delves deeper into heuristic approaches when compared to previous works, which is a necessity as roads development projects increase in complexity (discussed in more detail later). Their model did not, however, consider access to specific service locations. Murawski and Church (2009) proposed an approach which deter- mines solutions that simultaneously upgrade the road network’s accessibility to health care locations, while minimizing the associated financial costs. Their model was solved using commercial optimization software and was used to search for optimal solutions at different maximum budget levels. They went on to show that treating existing facility locations as fixed and upgrading transport network links to all-weather roads can substantially increase all-season access to health services – even with modest investments. Unfortunately, their model only considered one service type, which is not practical in a real-world environment where the building of roads should consider accessibility to more than just one service. Specifically in the context of Nepal, Shrestha et al. (2014) explored a bi-objective model for rural road network construction and upgrading in hilly regions of Nepal. Their objectives were to minimize the overall user operating costs (which can be thought of as a roads usage difficulty measure) and to maximize the population covered (serviced) by the network. The model then approximates solutions which provide optimal trade-offs in these objectives, determined for different maximum budget values. Similarly, Shrestha (2018) presented a multi-objective decision support model for rural road networks that suggests road network links for improvement or new construction. Using a similar approach to Shrestha et al. (2014) and with the objectives minimizing four population- and distance-based criteria, optimal trade-off solutions were approximated at different budget levels. In this paper we implement two solution approaches that are popular for roads planning, and which were used in the aforementioned literature. Heuristics are considered because they are easily compatible with rural roads planning (Ma et al., 2014; Scaparra and Church, 2005; Kumar and Kumar, 1999) and can obtain numerous well-performing solutions within reasonable computation times when the complexity of the problem is too overwhelming for traditional methods. An integer-linear programming (ILP) approach is the second method that we implement, because it provides a simpler way to solve optimization problems compared to our heuristic approach and has been used in rural roads planning before (Murawski and Church, 2009; Shrestha et al., 2014; Shrestha, 2018). An ILP approach is, however, more sensitive to computational complexity and the size of problems which it can solve, and returns fewer solutions when compared to a heuristic approach – resulting in less data and less insightful solution-analyses. 7 2.5. Existing techniques of rural road planning in Nepal Nepal’s forbidding terrain, low level of development, and underdeveloped data ecosystem complicate effective measurement of economic gains or accessibility changes from planned roads investments, let alone the implementation of advanced optimization models. Moreover, models built for planners must balance technical sophistication with communicative efficiency to policy makers and the public. As such, policy makers in Nepal and India have favored rough, easily calculable measures of accessibility to set policy and measure progress. Nepali policy makers traditionally have measured accessibility policy goals in travel times, e.g. “85% of children should live within 30 minutes walking distance of a school” (Ministry of Education, Nepal and UNESCO Office Kathmandu, 2015). These figures are typically collected by interviews with local administrators or in one-off survey instruments, making them prone to manipulation and impossible to scale. Other accessibility metrics include reaching target dates for all villages of certain population sizes to be connected by all-weather roads (Kumar and Tillotson, 1991) or connecting a given percentage of key population groups (caste, ethnic, or age-based) to roads (LAHURNIP and IWGIA, 2014). Accessibility models are rarely applied in the planning of Nepalese roadworks, where the level of method- ological sophistication falls sharply from major to minor roads. Central authorities and collaborating donors typically evaluate major proposed “strategic” roads using welfare-based economic valuation methods. For instance, the World Bank’s forthcoming Strategic Road Connectivity and Trade Improvement Project (SR- TIC) prepared a cost-benefit analysis of 12 road segments evaluating with- and without-project scenarios on the basis of combined vehicle operating cost and travel time savings (World Bank, 2019). But in practice this rigor is unusual and at a local level little to no proper evaluation is done (Coburn, 2020); in interviews in August 2018 and March 2019, staff at the Rural Access Programme spoke of how municipal authorities keen to demonstrate progress but lacking planning capacity frequently plowed roads into hillsides with little regard for environmental outcomes, investment efficiency, or even basic engineering principles. Some efforts are underway to equip local authorities with access-based planning tools that accommodate their limited data collection and analysis capacity. Over the last 15 years the Rural Access Programme has worked with local governments in Nepal to design and implement Transport Master Plans that provide a measurable, data-driven basis for prioritizing proposed roads. The most recent iteration is the Provincial Transport Master Plan, a standardized tool for provincial governments to prioritize roads based on a spread- sheet matrix of cost projections, one day traffic density studies, policy priority, and the population totals of administrative centers touched by roads (Rural Access Programme (RAP3), 2019). An extract of the prioritization tool within Karnali’s draft Provincial Transport Master Plan is shown in Figure 4. Figure 4: Ranking road building options in Karnali’s draft Provincial Transport Master Plan (Rural Access Programme (RAP3), 2019). 8 3. Data and methods 3.1. Inputs The conditions of data production in Nepal mirror the physical conditions of the roads themselves: frag- mented, incomplete, and least comprehensive for the most remote areas (Bikas Udhyami, 2018). In this environment, gathering sufficient primary data to accurately assess travel times and accessibility is difficult. Our model manages these issues by using global and freely available datasets – Shuttle Radar Topology Mission (SRTM) 30 m terrain elevation model, High Resolution Settlement Layer (HRSL) raster population data sets, OpenStreetMap, national land cover products (ICIMOD 2013) – and core governmental datasets like roads, health care, government offices, and financial institutions that a partner government agency can be reasonably expected to possess. We employ the HRSL data set instead of the WorldPop data set used in Banick and Kawasoe (2019) because it represents uninhabited areas as “Null” instead of fractional values as in WorldPop, and thus better capture Karnali’s unpopulated expanses, illustrated in Figure 3(b). In all cases OpenStreetMap data can provide a suitable “backstop” data set for missing or inaccurate governmental datasets. The data sources that were used to model accessibility to services are summarized in Table 1. Table 1: Data sources used in this research. Resolution Data Source Date or Type Official road network of Nepal Survey Department of Nepal 2019 Vector (Lines) Road and path network of Nepal OpenStreetMap Dec 2019 Vector (Lines) Trail Bridge Infrastructure Trail bridges Aug 2019 Vector (Points) Support Program (TBIS) Karnali provincial government Proposed roads in Karnali province Aug 2019 Vector (Lines) (via RAP) NASA Shuttle Radar 1 arc-second Terrain Global Topology Mission (≈ 30m × 30m) Integrated Center for Mountain Landcover of Nepal 2020 1km × 1km Development (ICIMOD) Rivers of Nepal OpenStreetMap Dec 2019 Vector (Lines) Health facilities Nepal Ministry of Health 2016 Vector (Points) Karnali municipal headquarters World Bank 2019 Vector (Points) Formal financial institutions Nepal Rastra Bank 2017 Vector (Points) (classes A - D) (central bank) High Resolution Settlement 1 arc-second CIESIN and Facebook 2018 Layer Maps (HRSL) (≈ 30m × 30m) Geospatial data for health services and formal financial institutions were sourced from the Ministry of Health (2016) and Nepal Rastra Bank (2017), respectively, while municipal headquarters were compiled from assorted Nepali-language policy documents and then geolocated with reference to OpenStreetMap. This basket of services was selected as a representative of a cross-section of developmental, economic, and political needs; Figure 3(a) shows their current distribution versus existing roads. Effective health care is an essential input to human development, particularly in the region of Nepal most afflicted by food insecurity and consequent malnutrition. Formal financial institutions (hereafter referred to as financial institutions) contribute to economic development by enabling the growth of local small to medium enterprises (SMEs) – unfortunately, no suitably comprehensive data set of informal financial institutions exists. Finally, municipal headquarters not only provide essential administrative services but are acceptable proxies for principal local markets. Education was considered for inclusion but a suitably comprehensive geospatial data set could not be located. The Rural Access Programme (RAP), a donor-funded infrastructure support program active in Karnali, shared shapefiles for the 92 individual proposed road segments specified in Karnali’s draft Provincial Trans- port Master Plan as of early 2019 and the average per-kilometer budget charge associated with the different 9 physio-geographic regions in Karnali: Low Hills, Middle Hills, High Hills, and Mountains. These charges ag- gregate materials, labor, management, and administrative expenses. We followed the government’s practice as described by RAP in an August 2018 interview, allocating per-kilometer budget charges to road segments based on their length and the majority geography of the district they primarily fell within. Thus Mountain roads were 33% more expensive than Low Hill roads. The values that were used are provided in Table 2. Table 2: Average per-kilometer budget charge associated with different physio-geographic regions. Values were provided in NPR and converted to USD using an exchange rate of 120 NPR to 1 USD (2 May 2020). Road cost Road cost Geography (NPR / km) (USD / km) Low Hills 15 million 126 thousand Middle Hills 17 million 143 thousand High Hills 18 million 151 thousand Mountains 20 million 168 thousand 3.2. Roads vs. road sequences Roads do not exist in isolation, but as part of networks. As such, proposed roads may implicitly depend on the construction of other roads linking them to the existing road network. We reviewed the proposed roads and grouped them into combination sets with internal hierarchies (tiers) reflecting these spatial dependencies. Thus a tier 2 road enables a tier 3 road to be built but depends in turn on a tier 1 road being built. Applying these hierarchical concepts to generate all plausible combinations of road tiers expands the number of candidates from 92 roads to 261 combinations – called sequences – and corresponding GIS files were created from visually observable hierarchies. In a real-world scenario, such hierarchies would be identified by experts in road construction and regional accessibility; the priority here was to illustrate our process and its accommodation of such hierarchies. The accessibility improvement modeling process was then run for these sequences as well. The notion of such tiers that form part of the sequences identification process is illustrated in Figure 5 for a number of connecting roads ranked into three tiers. For narrative simplicity, in the remainder of this paper all potential additions to a candidate layout are referred to as a road sequence – whether a standalone road or a combination of roads. 3.3. Geospatial cost-time analyses A recently released geospatial model of Nepal representing distances in units of time offers a scalable, replicable solution to the problem of modeling travel in Karnali (Banick and Kawasoe, 2019).2 The model calculates the minimum travel time to cross every 30 m2 of Nepal’s two dimensional surface based on the underlying terrain, landcover, altitude, roads, and road conditions; an accompanying monsoon model adjusts these factors and resulting travel times for monsoon conditions. The two dimensional output is packaged into a 30 m2 gridded (“raster”) friction surface that can be used to calculate accessibility to various services from any place on the grid. These accessibility rasters can be aggregated into descriptive statistics at almost any useful level of analysis – administrative areas, catchments for specific markets or hospitals, or catchments for individual roads. Applying the friction surface to services of interest allowed us to specify the precise accessibility gains associated with each planned road or sequence of roads in Karnali province. We measured the increase in accessibility associated with specific proposed roads by generating a friction surface for every proposed road incorporating its projected travel time improvements, then using said model to calculate accessibility to selected nearby services. Comparing these to similar measures from a reference data set based on the existing road network yielded the accessibility improvement associated with an individual road. 2 Readers with a deeper interest in the process of cost-time modeling and subsequent accessibility analysis in Nepal are encouraged to refer to Banick and Kawasoe (2019), where the modeling process and parameters are discussed at length. 10 Figure 5: Roads do not exist in isolation, but as part of networks, and some proposed roads may depend on the construction of other roads – resulting in internal hierarchies (tiers). Numerous possible combinations exist for such roads when selected with consideration given to their tier rank, and the resulting combinations of roads are called sequences. 3.4. Cost-time modeling The first component of our methodology involves calculating the accessibility to services by means of cost-time models. The backbone of our cost-time models comes from a layer of on-road/on-path speeds. The principal data input is a transportation layer produced by merging roads and paths from the Nepal Survey Department (2019) and OpenStreetMap with roads surfaces and conditions data compiled in 2015 for the Rural Access Index. Road segments were assigned speeds based on their surface type and Department of Roads classification (Department of Roads, 2013). These speeds were then modified according to the condition and gradient of the road. Using modifiers provided by a focus group of World Food Programme civil engineers, we then calculated a separate monsoon speed for each segment based on the average deterioration in speeds over given road surface/condition combinations under monsoon conditions. A separate interim model of off-road, off-path walking speeds was then prepared. To do so, we first calcu- lated the slope for each grid cell in the SRTM 1 arc-second (≈ 30 m2 ) digital elevation model. These slopes were input into Tobler’s hiking formula to derive a base walking speed for each grid cell. We multiplied this speed by a fractional landcover walking speed modifier using a 1 km2 gridded landcover model of Nepal from the International Centre for Integrated Mountain Development (ICIMOD, 2013). A downward adjustment in walking speeds was applied to extreme altitudes (greater than 3500 m) where lower oxygen counts reduce aerobic capacity. Rivers and lakes from OpenStreetMap were rasterized and “burned” into the walking data set as impassable features that could only be crossed where roads or trail bridges existed. The off-road and on-road interim models were merged into final friction surfaces showing the average speed to cross each 30 m2 of Karnali province using the fastest available means of transport (Figure 6(a)). The models were merged by overlapping them and taking the lowest positive value (quickest speed) in each cell. This ensures road or path speeds are used where they exist and walking speeds where not. Two final mesh speed models were produced, one each for the dry season and monsoon season. The decrease in speeds during monsoon season varies significantly based on the condition and types of underlying roads, ranging from 15% to 90% of dry season speeds. Walking speeds are uniformly reduced by 25%. These mesh speeds are multiplied against a cost-distance surface representing the combined linear-vertical distance for each cell (Figure 6(b)). The purpose of cost-time modeling is to prepare the necessary tools for analyzing accessibility to chosen 11 (a) (b) (c) Figure 6: (a) Combining on-road and off-road speeds to find a fastest-possible mesh speed, (b) modifying the mesh speed by vertical distances to arrive at a final friction surface, and (c) using the friction surface to generate a service accessibility model. services. By referencing the final, friction surface against a geospatial point data set of service locations the average travel time from every cell in the raster to the nearest service can be calculated (Figure 6(c)) and then visualized. An example of the results of the discussed methods are presented in Figure 7, in which accessibility to health facilities in the Dolpo region of Karnali is provided in Figure 7(a), and the updated accessibility if road sequence 44 is built in Figure 7(b). The net gain in accessibility can be visualized as the difference between accessibility before and after construction of the sequence, as is illustrated in Figure 8 for sequence 44 that was introduced in Figure 7. We did so for each possible sequence with respect to health care facilities, financial institutions, and municipal headquarters, for both dry and monsoon seasons. In our analysis we scale this to a batch process, repeating the cost-time and accessibility modeling for each proposed road and sequence (Figure 9(a)). To improve performance, small scale friction surfaces confined to the extent of the proposed road/sequence are quickly generated applying the road speed parameters to the new, improved roads, then merged with the existing roads’ friction surface. This proposed road-specific friction surface is used to calculate one proposed road-specific accessibility models for each service, covering the whole of Karnali province (Figure 9(b)). A total of 261 proposed roads and sequences resulted in 783 12 (a) (b) Figure 7: Comparing models of accessibility to health facilities in Dolpo using (a) current roads, and (b) if sequence 44 is constructed. 13 Figure 8: Visualizing the improvement in accessibility to health facilities after building sequence 44 in Dolpo. The net effect of the sequence’s construction is determined as the difference in accessibility before construction (Figure 7(a)) and after construction (Figure 7(b)). proposed road service layers, and performed with respect to each season this resulted in 1,566 layers in total. 3.5. Calculating accessibility improvements Where friction surface cells overlap with populated places (such as the black cells in Figures 7(a) and (b)) they provide a precise measure of accessibility to each service for the underlying population. Multiplying the travel times in an accessibility model cell by the population density underneath yields a composite measure of person-hours at a given location (illustrated in Figure 10(a)) – for brevity, person-hours is referred to as “p-h” in all figures from here on. This measure is advantageous as it balances population density against overall remoteness: high travel times will offset low population densities, ensuring remote and sparsely settled populations’ needs are not drowned out by aggregate numbers from densely settled populations nearer to towns. Using the HRSL improves the accuracy of this analysis as HRSL imputes Null values into unpopulated areas, ensuring that Karnali’s remote expanses do not unduly weight results. We aggregate the cumulative person-hours for each road across the whole of Karnali. Then we take the difference in aggregate person-hours between existing and proposed roads / season / service accessibility combinations (Figure 9). This yields the improvement in accessibility to a service resulting from a proposed road or sequence in a given season. By repeating this process for every proposed sequence (Figure 11) we generate six objective measures of accessibility improvements associated with each – that are tabulated in a spreadsheet for entry into optimization models (Figure 12). 3.6. Optimization The selection of which road sequences to include in a candidate plan is a process that is combinatorially complex. Budgetary constraints make it infeasible to construct all proposed roads so a subset needs to be selected in a manner that optimizes the resulting accessibility improvements. The consideration of multiple accessibility criteria, considered as multiple maximization objectives, results in the existence of numerous candidate plans comprising various sequence combinations – each providing trade-offs in optimality with respect to the different objectives. The approach that we followed to obtain such candidate solutions is described in this section. 14 (a) (b) Figure 9: (a) Efficient batch generation of season-specific friction surfaces for every proposed road, and (b) batch generation of accessibility models for 3 services for each of 261 potential road sequences in 2 seasons. 15 (a) (b) Figure 10: (a) Computing person-hours saved by multiplying accessibility models against population models, and (b) Calculating the relative improvement in accessibility (measured in person-hours) per road/service/season combination compared to the existing roads configuration. Figure 11: Batch generation of person hour rasters for each road sequence/season/service combination Figure 12: Final outputs of the accessibility improvement modeling process are packaged into tabular forms for use in optimization models 16 Our approach considers a maximum and minimum budget specification and candidate plans therefore provide alternatives within this range. A benefit of following such an approach is that it allows for policy makers to examine trade-offs in construction costs and avoid unnecessary spending. For example, if a plan that costs 5% less than a plan at the maximum budget, yet achieves accessibility improvements that are similar to those which may be achieved at the maximum budget, then the 5% in savings could be allocated to other development projects. Our approach does not allow for the same road to be selected for inclusion in a candidate plan more than once – when a candidate plan is configured during the optimization process, sequences that share roads with those already included in the plan are excluded from consideration. For example, if multiple sequences include road DR5007, then a candidate development plan can only include one of these road sequences and all the others (and the standalone use of DR5007) are excluded from consideration. As more sequences are added to a potential plan, the list of exclusions grows. The constraint on duplicate roads selection as discussed here is an important factor to consider during the mathematical problem formulation and optimization processes so that practical solutions are obtained, and is further elucidated in the remainder of this section, where relevant. Finally, our approach considers 6 accessibility maximization objectives (3 services per season) in this study alone and potentially more in future work. Depending on the solution approach followed, the criterion of cost may be considered either as a budgetary constraint (enforcing that a candidate plan’s costs are between maximum and minimum limits) or as an additional cost minimization objective. 3.6.1. Candidate solutions and evaluation In multi-objective optimization, a candidate solution (roads plan) is evaluated by multiple objective func- tions – mathematical functions which calculate the plan’s performance with respect to each of the objectives. In our approach the objective functions determine the sum of the accessibility improvements (in person- hours) of all the sequences in a plan, for each service per season. Similarly, based on feedback from key RAP staff during a December 2019 interview with key staff, the cost of a plan is the sum of the estimated cost of all its constituent sequences. The objective function values of a candidate plan correspond to a single point in objective function space – this is illustrated in Figure 13 by (hypothetical) example candidate plans represented as points and evaluated with respect to two accessibility maximization objectives (plotted on the axes). This example figure shows the solutions evaluated with respect to only two objectives, but the same principles apply for 3 or more objectives. Minimizing cost is often considered as a maximization objective, by maximizing 1/cost. In multi-objective optimization, solutions such as those in Figure 13 are classified either as non-dominated (superior) or dominated (inferior) (Zitzler et al., 2004). The set of non-dominated solutions exhibit superior trade-off alternatives when compared to the dominated solutions, and form what is commonly known as the Pareto-optimal front, or simply the Pareto front, as may be observed in the figure (Zitzler et al., 2003). Only the solutions in the Pareto front need to be presented to decision makers because of their superior quality; that is, dominated solutions are avoided. The Pareto front is determined using the same principles of domination for problems involving more than two objectives, and the notion of Pareto optimality remains applicable, albeit harder to visualize. 3.6.2. Heuristic solution approach The search for the exact (true) Pareto front for optimization problems that include multiple objectives involves a significant computational challenge and may become prohibitively large to solve within realistic computation times (Jia et al., 2007; Owen and Daskin, 1998; ReVelle and Eiselt, 2005; Tong et al., 2009; Xiao et al., 2002). In such instances, powerful heuristics are often employed in order to approximate the set of solutions on the Pareto front (Xiao et al., 2002; Zitzler et al., 2004; Heyns and van Vuuren, 2018). These heuristics target and explore promising regions of the solution space3 in order to arrive at solutions that are approximated to be Pareto-optimal, and in the process avoid the computationally expensive consideration of solutions in inferior regions of the solution space. It has been demonstrated that heuristics are well capable 3 The solution space in the context of this paper is the set of all feasible plans, e.g. all the combinations of sequences that avoid road repetitions (the exclusions) and with costs that fall within the monetary budget. 17 Figure 13: An illustration of the notion of Pareto optimality and the translation of candidate roads plans to objective function space. Numerous candidate solutions are evaluated with respect to the objective functions, and those that are determined to be Pareto optimal are sought for decision making purposes. of determining the true Pareto front within reduced computation times when compared to exact approaches (Kim et al., 2008; Heyns and van Vuuren, 2016, 2018). The class of Multi-Objective Evolutionary Algorithms (MOEAs) is one heuristic approach that may be considered to approximate a diverse set of trade-off solutions on the Pareto front (Fonseca and Fleming, 1993; Purshouse and Fleming, 2003). These algorithms iteratively evolve an initial, randomly generated population of candidate solutions through carefully controlled evolution (based on natural principles) over multiple generations, finally arriving at a set of solutions that approximate the Pareto front (Deb et al., 2002; Cheshmehgaz et al., 2015). Due to their population-based nature (numerous solutions simultaneously being discovered and evolved), numerous solutions along the Pareto front can be approximated in a single run. Example applications of Multi-Objective Evolutionary Algorithms include the optimization of route planning (Lau et al., 2009; Maji and Jha, 2017), emergency response (Georgiadou et al., 2010), financial engineering (Branke et al., 2009), and the deployment and power assignment of wireless sensor networks (Konstantinidis et al., 2010). The non-dominated sorting genetic algorithm-II (NSGA-II) (Deb et al., 2002) is a popular MOEA that was selected for implementation here (albeit a modified version described below), building on the previous research and applications of Heyns (2016), Heyns and van Vuuren (2018), and Heyns et al. (2019) in which the NSGA-II has been used for the placement of numerous facility types, such as radar, defense systems, and fire-detection cameras. The NSGA-II falls within the class of genetic algorithms (Deb et al., 2002). Genetic algorithms have been implemented for the optimization of rural roads maintenance strategies (Mathew and Isaac, 2014) and rural highway layout planning (Ma et al., 2014) because their candidate solution representation scheme is easily applied in the context of rural roads planning. Genetic algorithms represent solutions as chromosomes, and a candidate roads plan is represented as a chromosome string of sequence numbers. For example, a chromosome [23, 53, 99, 104, 156, 201] represents a plan with sequences 23, 53, 99, 104, 156 and 201 selected for construction. Evolution-inspired selection processes and modification operators are iteratively performed on a randomly generated population of such chromosomes and terminates when the algorithm converges – e.g. when it ceases to show significant improvement in the solution quality of successive generations (Deb et al., 2002; Heyns, 2016). 18 Figure 14: Examples of typical crossover procedures that exist for VLCs. In one-point crossover, parent chromosomes are “cut” and combined, resulting in two new offspring chromosomes. In uniform crossover, parent chromosomes of different lengths are selected and entries between parents are randomly paired and swapped, also resulting in two new offspring chromosomes. Two methods are employed to generate new solutions during the evolution of populations: crossover and mutation. Crossovers (swaps) performed between sub-strings of selected “parent” chromosomes create new “offspring” solutions by mixing the parents’ sequences with each other (Deb et al., 2002). Mutation promotes the diversity of sequences in the chromosomes – and the population in general – by stochastically introducing new, unexplored sequences into the chromosomes (Deb et al., 2002). Our problem required an unconventional approach, so that traditional crossover and mutation operators were not implemented. We desire diverse solution alternatives across minimum and maximum budget values – a requirement that is made possible to achieve in a single run because of MOEA solution processes, together with variable- length chromosomes (VLCs) (Cramer, 1985; Srikanth et al., 1995; Ting et al., 2009). Naturally, crossovers and mutations are subject to the constraint that a single road may only be included once in a candidate plan, which is not a sophisticated requirement and simply requires the implementation of duplicate-checking procedures within constituent sequences, discussed in more detail later. Road sequence costs vary, which means that feasible candidate plans with varying sequence lengths exist (whether considering a specific fixed budget or a budget range) making a VLC application particularly essential. Crossover approaches that have been investigated for VLCs include one-point crossover (Ting et al., 2009; Ripon et al., 2006; Srikanth et al., 1995) and uniform crossover (Ting et al., 2009). In one- point crossover, illustrated in Figure 14, one point along each of two chromosomes is selected, where the chromosome is “cut” and each of the two parts from one chromosome is combined with one of the two parts from the other, resulting in two new offspring chromosomes (this is also a common approach in fixed-length chromosome applications). In uniform crossover, also illustrated in Figure 14, two chromosomes of different lengths are selected and each sequence in the shorter chromosome is swapped with a randomly chosen sequence in the longer chromosome, also resulting in two new offspring chromosomes. Mutation processes for VLC applications, illustrated in Figure 15, follow the standard approach of randomly selecting an entry (in our case, road sequence) in a chromosome and replacing it with an entry that is not included in the chromosome (Deb et al., 2002).4 [h] In order to explore our unique, budget-spanning solution space, our crossover approach was inspired by an implementation of Jia et al. (2007). In the first step of this approach, two entire chromosomes are combined and any duplicate entries are removed. Jia et al. then follow a “best-drop” entry-removal process from the Global/Regional Interchange Algorithm of Densham and Rushton (1992). In this approach, each entry along the chromosome is separately removed, reducing cost and objective function values. The removal of a strong entry will typically result in significant weakening of the objective function values, while the removal of a less significant entry will typically have a small impact. In other words, a weaker entry does not contribute 4 In binary chromosome approaches the selected entry is flipped, e.g. 0 to 1, or 1 to 0 (Ting et al., 2009; Srikanth et al., 1995). 19 Figure 15: Mutation processes for VLC applications follow the standard approach of randomly selecting an entry in a chromosome and replacing it with another randomly selected entry that is not currently in the chromosome. much to the solution’s objective function values and the weakest (least significant) entry may be identified and removed from the chromosome. In our VLC crossover approach, illustrated in Figure 16, the combination of two chromosomes may result in a new chromosome that exceeds the maximum budget or which is near the maximum limit. Such over- or high-budget solutions are, in fact, desired, since a method of repeated entry removal is introduced here and serves two purposes: 1) to return over-budget solutions back within feasible range, and 2) to explore solutions at various cost levels across the budget range. VLC crossover comprises three steps. These are 1) selecting two parents and removing any possible duplicate sequence entries, 2) removing sequences from one chromosome that are in the exclusion list of the other (the choice of which chromosome to remove entries from is randomly selected), and the chromosomes are then combined, resulting in the first potential offspring solution (Offspring 1 in Figure 16), and 3) repeated removal of the weakest sequence until the minimum budget is reached. As opposed to the implementations of Jia et al. (2007) and Densham and Rushton (1992) where only one entry is removed, in the third step of our implementation we continue determining and removing the weakest entries – thereby discovering a new candidate plan with a lower cost after each removal – until the chromosome cost goes below the minimum budget. In this manner, from the initial combined chromosome, multiple solutions that explore various cost levels within the budget range are discovered by the repeated removal of entries, while the removal of only the weakest entries promotes the discovery of stronger solutions. All the offspring solutions in the third step that are within the budgetary limits progress to an offspring population, which is then subjected to our novel mutation approach. The mutation process followed in our approach is inspired by a multi-objective modification to the classic Teitz-Bart algorithm (TBA) (Teitz and Bart, 1968), as investigated in the Hybrid Geospatial Algorithm (HGA) of Heyns (2016). The TBA is in essence a single-solution, single-objective algorithm that starts by randomly generating a single starting solution and repetitively improving it with multiple swaps. The typical approach is to identify all the candidate entries that are not included in the starting solution. One entry in this set is randomly selected and consecutively swapped with each entry in the starting solution (Teitz and Bart, 1968). If any objective function value improvements are achieved as a result of candidate swaps, the swap that results in the best improvement in the objective function value is accepted and results in the new current solution, and the entry that has been selected is removed from the set of remaining candidate entries. This procedure is then repeated, relative to a continually updated current solution and until all entries in the candidate set have been swapped. The current solution, after all candidate swaps have been completed, is accepted as the final, best solution. In our implementation the HGA’s mutation process is performed on each of the multiple offspring solutions from crossover, and utilizes a similar swap procedure to that implemented by the TBA (Teitz and Bart, 1968; Heyns, 2016). Our mutation approach is illustrated in Figure 17. Each entry in the chromosome is swapped with an arbitrarily selected sequence that is not in the duplicate-road sequence exclusion list of the other entries, and each swap that is performed results in a new offspring. Compared to the TBA – in which candidate solutions are evaluated with respect to a single objective – new chromosomes in our problem are evaluated with respect to multiple objectives. Whereas the TBA selects the single solution resulting from the swaps that results in the best improvement (if any improvement is found), our implementation accepts all offspring strings emanating from mutation that exhibit improvements over the parent with respect to any of the objective functions. The approach serves to rapidly explore very large solution spaces and returned 20 Figure 16: The VLC crossover process that was followed in this paper. In the first step, two parents are selected and any possible duplicate entries are removed. Then, sequences from one chromosome that are in the exclusion list of the other are removed and the chromosomes combined, resulting in the first potential offspring solution (Offspring 1). In the final step, weakest sequences are repeatedly removed until the minimum budget is reached. impressive results for large geospatial optimization problems when implemented within the HGA and when compared to the NSGA-II in the test problems of Heyns (2016). The HGA is, in fact, mostly the same as the NSGA-II, with the most significant difference being this mutation approach. Our heuristic’s processes from start to convergence are now described, and are illustrated in Figure 18. A population of N feasible solutions are stochastically generated, forming the current population of VLCs which cost within the budget range. Crossovers and drops are then performed on parents selected from the current population – parents which perform well with respect to the objective functions are favored for crossover selection – meaning that the offspring solutions typically exhibit some of the strong properties (strong sequences) of their parents (Deb et al., 2002). The crossover process first repeatedly combines randomly selected chromosomes until N feasible offspring are discovered (steps 1 and 2 from Figure 16), after which weakest drops are performed on each of these and those that are within the budget range proceed to the offspring population and to mutation. During mutation, all new offspring that show improvements and which are within budget progress to the larger offspring population together with the solutions from crossover. The current and offspring populations are combined (size larger than 2N ), and the N best solutions from this combined population are identified – performed so that the population size does not spiral out of control (resulting in intractable complexity) following repeated crossover and mutation processes of successive generations. This is achieved by the Fast Non-dominating Sorting Algorithm (Deb et al., 2002), which applies non-dominated principles and preferences for solution diversity to filter populations down to a specified size. If convergence is observed (when no significant improvement in the solution quality of successive generations is observed), the algorithm terminates, and the Pareto front approximation of the final population is returned 21 Figure 17: The VLC mutation process followed in this paper. Each entry in the chromosome is swapped with an arbitrarily selected sequence that is not in the exclusions of the other entries, and each swap that is performed results in a new offspring. for analysis. If the algorithm does not yet converge, the N solutions become the current population of the next generation, and the process is repeated until convergence is observed. Weighted-sum approach ILP solvers take as input mathematical formulations of an objective function and constraints. They return a single solution that optimizes the selection of decision variables (the candidate sequences) so that the objective function is maximized or minimized according to specification, while ensuring that the constraints are enforced. A strong characteristic of an ILP approach is that the corner points of the Pareto front, which effectively optimize with respect to only a single objective while ignoring the others, may be determined exactly (Arora, 2011). This is because ILP solvers are inherently designed for single-objective optimization and determining the corner points only considers one objective, while determining points with respect to multiple objectives (between the corner points) limits the solvers’ Pareto-approximation abilities (Stewart, 2007). These corner points provide an indication of the extent of the true Pareto front – and avoids a known weakness of heuristics which struggle to reach corner point regions (Kim et al., 2008). An illustration of such corner points is provided in Figure 19 for a problem with three maximization objectives. The shaded surface is the general shape of where other Pareto optimal solutions may be expected (the corners are also Pareto optimal), such as illustrated in the figure for five additional Pareto optimal solutions. In practical applications, the Pareto optimal solutions do not necessarily lie on the surface, but tend to follow the general shape. The points are projected onto the obj-1/obj-2 plane to help with viewing perspective and to provide an indication of the solutions’ objective function values. In order to determine the solutions between the corner points with these packages requires transforming the multiple objective functions into a single objective function using a weighted sum (Cohon, 1978; Murray et al., 2007). In its simplest form, a weighted-sum objective function for a problem with Ni objectives is given by Ni maximize Os = w i Vi , (1) i=1 22 Figure 18: An overview of the custom NSGA-II-inspired VLC heuristic process followed in this paper. The approach borrows solution search methods from Jia et al. (2007), Teitz and Bart (1968), Densham and Rushton (1992), and Heyns (2016), and is implemented within the NSGA-II framework of Deb et al. (2002). where the objective function values achieved with respect to each objective, Vi , are combined using weights wi that reflect the relative importance of each objective. By varying the objective weights in multiple runs, the solver may be “directed” to a specific region between the corner points, and Pareto optimal solutions along the shaded area (such as those in Figure 19) may be traced out. This weighted-sum method provides a simple way to solve MO optimization problems, since the mathe- matical formulations are relatively simple to provide as input compared to the sophisticated code and multiple parameters and sensitivity tests that are required for heuristics. Due to its accessibility and historical promi- nence the weighted-sum approach remains a popular choice to solve MO optimization problems from various recent applications, such as building design (Machairas et al., 2014), carpool matching (Xia et al., 2019), and land-use allocation (Yao et al., 2018), and has also been used for rural roads planning research (Shrestha et al., 2014; Shrestha, 2018). A popular ILP alternative to the weighted-sum approach is the ǫ-constraint method (Mavrotas, 2009; Mavrotas and Florios, 2013), which has been shown to be capable of obtaining superior solution spread along the Pareto front approximation in certain problem instances Amin and Zhang (2013). However, the ǫ-constraint method requires constant fine-tuning and interaction during the solution process – a burden which is compounded as the number of objectives increase. The method is therefore not followed here because of the large number of objectives considered – but remains an option to investigate in future work. The weighted-sum approach does, nevertheless, hold disadvantages. Evenly distributed weights may nevertheless result in an unevenly distributed Pareto front approximation (Das and Dennis, 1997), and while a truly optimal solution may be found for a specific weight combination, there is no guarantee that the optimal solution for a specific weight combination lies near the true Pareto front (Marler and Arora, 2004; Arora, 2011; Khan and Rehman, 2013). Furthermore, determining points on the Pareto front in this manner may require a large number of weight combinations when many objectives are considered (ReVelle and Eiselt, 2005; Tong et al., 2009) – as is the case here for six accessibility objectives – and determining suitable weight combinations is a sensitive and laborious process (Marler and Arora, 2004). When different cost levels are to be considered (such as in our problem), the multiple weight combination runs also need to be repeated at each cost level, with cost either considered as a maximum constraint or as an additional minimization objective. The mathematical formulation to our problem is now presented. 23 Figure 19: An example of a three-dimensional Pareto front determined with respect to three maximization objectives. The corner points are determined with consideration given to a single objective only, while other solutions along the front may be determined using a weighted-sum objective function. Ns denotes the number of sequences available for consideration. Na denotes the number of accessibility criteria considered. i, j denote indices of sequences, where i, j ∈ {1, . . . , Ns }. a denotes the index of accessibility criteria, where a ∈ {1, . . . , Na }. vi,a is the accessibility improvement value that sequence i achieves with respect to criterion a. xi is 1 if sequence i is selected, and 0 otherwise. xj is 1 if sequence j is selected, and 0 otherwise. Xi,j is pre-determined for all i, j pairs, and equals 1 if sequence i has an overlap (common road) with sequence j , and 0 otherwise. ci denotes the cost of sequence i. C denotes the maximum plan cost (maximum budget). The objective can then be written as Ns maximize Vc = xi vi,a ∀ a ∈ {1, . . . , Na } (2) i=1 subject to the constraints Ns Ns xi xj Xi,j ≤ 0 (3) i=1 j =i+1 Ns xi c i ≤ C (4) i=1 xi , xj ∈ {0, 1}. (5) 24 Figure 20: Objective function values can vary significantly in minimum and maximum values of magnitude, and it is required to transform these values so that their magnitudes are similar (e.g. to normalize the values). This allows the objective weights in the weighted-sum approach to more effectively direct the search towards desired solution regions. The objective in (2) is to maximize the total improvement with respect to each accessibility criterion a ∈ {1, . . . , Na }. Constraint (3) ensures that no sequence overlap (duplicate road selection) occurs – if sequences at i and j are selected, i.e. xi = 1 and xj = 1, the constraint requires that Xi,j = 0 (no overlap between i and j ). Constraint (4) ensures that the combined cost of all the sequences selected in the solution does not exceed the maximum cost, while constraint (5) specifies binary requirements on the decision variables (a sequence is either selected or not, so roads are not “half built”). In the second stage of ILP optimization, the objectives in (2) – which were determined separately with respect to each individual criterion – are reduced to a single objective function using various combinations Na of weights, wa , assigned to each accessibility criterion, where a=1 wa = 1. When following a weighted- sum approach, the values of the weights are, however, not only significant relative to each other but also relative to their corresponding objective functions (Marler and Arora, 2010). This is because the relative magnitudes of the separate objective functions may differ significantly. Consider a scenario such as that illustrated in Figure 20, where the maximum objective function values that may be achieved for objective 1 is an accessibility improvement of 100k person-hours, and for objective 2 is an accessibility improvement of 250k person-hours. If we want to determine a solution that lies between the two corner points of objectives 1 and 2 – e.g. a solution that reflects equal importance between these corner points and does not consider objective 3 – weights of w1 = 0.5, w2 = 0.5 and w3 = 0 would be assigned. However, the objective function effectively aims to maximize (w1 V1 + w2 V2 ), which is not dictated solely by weights or relative magnitudes alone. In this case, the weighted sum will be dominated by criterion 2 because of the larger magnitude of which it is capable and the returned solution will be located towards (or even on) the corner point of objective 2. When using weights to represent the relative importance of the objectives, it is therefore required to transform the objective function values so that their magnitudes are similar and so that some do not naturally dominate the aggregate objective function (Marler and Arora, 2010). In order to achieve this, the objective functions are normalized so that they range between zero and one, so that the importance of weights is reflected more accurately on the resulting solutions. The maximum objective function values with respect to 25 each criterion are already determined in stage 1 of the ILP approach, in (2), and are rescaled to a value of one. The minimum value with respect to each objective function is rescaled to zero – these minimum values remain to be determined. However, this process is not overly sophisticated since the values may be found at the corner points not associated with the objective. Consider Figure 20 once more, where we aim to find a solution between the corners of objective 1 and objective 2. The minimum value that may be achieved with respect to objective 1 is simply found at the corner of objective 2, and the minimum of objective 2 is found at the corner of objective 1 – when considering objectives 1 and 2 only. Using the minimum and maximum values for each objective function, values determined between these values may be scaled to between zero (at the minimum) and one (at the maximum). In this manner, the effect of objective function magnitudes on a weighted-sum may be eliminated and the search for solutions is directed by the weights in the weighted-sum. If we were to search for solutions between all three corners, the minimum of each objective still lies at one of the other corners, and the value at each corner needs to be compared to the others before the minimum may be identified. For example, when only considering objectives 1 and 2, the minimum for objective 2 lies at the corner of objective 1. Now, if objective 3 is also considered, then the minimum of objective 2 may be found at either the corner of objective 1 or objective 3. If we look at the figure and the projections of the corner points onto the obj1/obj2-plane, then the minimum for objective 2 is now found at the corner of objective 3 and this minimum value becomes zero relative to other values. Let A ⊆ {1, . . . , Na } denote the set of criteria for which a weighted objective function is sought, with the weights of the selected criteria adding to 1. Furthermore, let va,max denote the maximum value of accessibility criterion a ∈ A, determined in ILP stage 1. Similarly, let va,min denote the minimum value of accessibility criterion a ∈ A, which is pre-determined as the lowest value found with respect to the objective at the Pareto corner points of the criteria included in A. The normalized objective function is then to maximize Ns ( i=1 xi vi,a − va,min ) maximize V = wa ∀a∈A (6) (va,max − va,min ) a∈ A The objective in (6) is subject to the same constraints (3) to (5). 3.7. Optimization process setup The accessibility improvement in person-hours was determined for each of the 261 sequences, for each season and with respect to health facilities, financial institutions, and municipal headquarters – therefore 6 accessibility criteria in total. Minimum and maximum budget limits of 6 and 13 billion NPR were specified. We used 13 billion NPR as the upper bound because this was the budget for the roads that were eventually selected by the provincial government. The much lower limit of 6 billion NPR provided a large budget span, which allows us to compare rates of accessibility improvement that may be achieved for different services at different budget ranges, and to compare how improvements between services relate to each other across the budget range (e.g. linear or exponential). The objectives in both the ILP and heuristic approaches were to find solutions so that person-hour accessibility improvements with respect to all 6 criteria are maximized. 3.7.1. Heuristic Cost was implemented as an additional, seventh objective by the heuristic, with the aim of minimizing total cost – implemented as a maximization objective, maximizing 1/cost. By adding the cost minimization objective, the heuristic is forced to not only find solutions that maximize the accessibility objectives, but also to explore the maximization of these objectives at lower cost levels, from 13 down to 6 billion NPR. Because of the stochastic nature of the heuristic’s Pareto-front approximation process, the solutions returned by different optimization runs will generally vary in quality. It is therefore standard practice to repeat the process multiple times (Knowles et al., 2006; Kim et al., 2008; Heyns et al., 2019). The results of all the runs are then combined and a final attainment front – essentially the Pareto front of all the combined Pareto front approximations – is identified. We repeated our process 40 times, from which a final attainment front of 5,000 solutions that are approximately Pareto optimal with respect to all 7 objectives was identified. The authors’ personal optimization code was used for the heuristic runs, executed in MATLAB. 26 3.7.2. ILP weighted-sum In the weighted-sum ILP approach, cost was considered as a maximum constraint and solutions were sought for various weight combinations at maximum cost levels beginning at 6 billion NPR and increasing in increments of 1 billion NPR up to 13 billion NPR, resulting in eight cost levels. In the first ILP stage, no weights are assigned because the objective in (2) is repeated for each individual accessibility criterion, at each cost level in order to determine the corner points on the Pareto front. In the second stage, equal weights were assigned with respect to all pairs of two criteria that can be compared from the six, e.g. w1 = 0.5 and w2 = 0.5, then w1 = 0.5 and w3 = 0.5, and so forth for each pair, resulting in twenty such weight combinations at each cost level. Then, all combinations of three criteria were equally weighted together, e.g. w1 = 0.33, w2 = 0.33, and w3 = 0.33, and so forth, resulting in 15 such weight combinations at each cost level. No more weight combinations were investigated, as these were primarily used for visual comparison to the heuristic results and visualizing solutions with respect to more than three objectives is challenging. The six runs from the first stage together with the thirty five weighted runs from the second stage totaled 41 runs at each of the eight cost levels, resulting in a global number of 328 weighted-sum runs. The runs were all executed with CPLEX and called from within MATLAB. This is because data preparation was done in MATLAB, and for each CPLEX run a different set of variables (depending on the criteria for which weights are assigned) are selected and provided as input – this was managed by MATLAB. 4. Results 4.1. Pareto fronts The solutions obtained by the heuristic and weighted-sum runs are presented in this section. It was possible to determine each set of heuristic and ILP results (for all heuristic runs and weight combinations) within one day.5 The 40 heuristic runs each terminated in around 10 minutes or less, while the ILP run times showed more variation because of different weight sets and budget limits in each run, but averaged at 13 minutes. Figure 21 visualizes the results of our heuristic runs in the form of the 5,000 solutions plotted simultaneously in 3D objective space for both seasons. It should be noted that these 5,000 solutions are the same solutions in both these 3D plots – although all solutions were determined with consideration given to 7 objectives, they are plotted against different axes because plotting them with respect to more than 4 objectives (3 accessibility axes and cost in color as in Figure 21) is visually challenging. The same 5,000 solutions are examined in the remainder of this paper. The cost of the solutions in Figure 21 are indicated as increasing in darkness from light gray at 6 billion NPR, to dark gray at 13 billion NPR. Solutions number 4,845 and 3,449 are explicitly linked between their solution “dots” and their roads layouts show representative examples of optimal solutions at various budget ranges. These layouts are shown alongside the layout and corresponding objective function values of the actual roads selected by Karnali’s provincial government. Figure 21 exemplifies how solutions may be presented to decision-makers to assist their analysis of options. A more evolved presentation might provide this data in interactive form. These results are analyzed in more detail in the following section. In Figure 22(a)–(c), the 5,000 solutions are now presented in 2D, along axes which are pairs of accessibility improvements in person-hours with respect to different services, for the dry season. A legend to describe the colors and markers used in Figure 22(a)–(c) is provided in (d). Selected solutions are indicated in the colors orange, magenta, blue, and red, which represent solutions at budget levels of 7, 9, 11 and 13 billion NPR, respectively. When examining the color solutions at each level, it may be observed that three are larger circles – these are the ILP solutions that were determined with the budget set at the indicated maximum value. Of these three ILP solutions at each budget level, the two outer points at each budget level are the values determined by single-objective optimization, while the middle solutions were determined by the weighted-sum with respect to the two objectives in each plot (i.e. determined with two equal weights for the two objectives). The smaller color circles are heuristic solutions that are determined at each budget level and collectively form a Pareto front extracted with respect to the two objectives in each plot. To aid in understanding these 5 All optimization processes in this paper were run on a Dell 7820 Precision desktop PC, running Windows 10 Pro with an Intel Xeon Silver 4110 processor and 64 GB memory. 27 28 Figure 21: A summary of the heuristic optimization outputs. Solutions 4845 and 3449 are plotted against all 6 accessibility maximization and the cost minization objectives, with objective function values displayed. The actual roads selected by Karnali’s provincial government in the draft Provincial Transport Master Plan are visualized for comparison. (a) (b) (c) (d) Figure 22: (a)–(c) Objective function values (in person-hours) of the 5000 heuristic solutions and selected ILP solutions, plotted against values achieved with respect to pairs of different services types, during the dry season. A legend is provided in (d) for the images in (a)–(c). 29 additional “sub”-fronts, one can consider the 5,000 heuristic solutions as a subset of the entire solution space of all possible solutions, and which have been identified to return optimal trade-offs with respect to the 6 accessibility criteria, and cost (i.e. a 7-dimensional Pareto front of 5,000 solutions). Then, similar to the process of determining a Pareto front from the entire solution space with respect to some objectives (see Section 3.6.1), we are now determining sub-fronts from the 5,000-solution space with respect to selected objective pairs. In Figure 23, the results are presented in the same manner as in Figure 22, for accessibility improvements during monsoon season. The same 5,000 heuristic solutions and selected weighted-sum solutions are now presented in Figure 24, with the axes in each graph representing accessibility improvements for each service with respect to the two seasons – (a) for health facilities, (b) for financial institutions, and (c) for municipal headquarters. The same legend in Figures 22 and 23 is again used here, for budget levels and ILP and heuristic solutions. 4.2. Analysis When analyzing the Pareto fronts in Figures 22 and 23, the trade-offs between improving accessibility to services become clear. At each budget level, moving along solutions that improve in quality with respect to one objective causes a degradation in solution quality with respect to another objective. This provides insight into the decision-making process, as it is clear that the objectives are conflicting and that the trade-offs should be examined. The solutions exhibit constant improvement in the objective function values when moving from one budget level to the next, indicating a relatively linear relationship between cost and accessibility improvement. Particularly encouraging is that the heuristic solutions at each budget level are similar in solution quality to those yielded by the ILP approach. This is significant, since the heuristic solutions are determined thousands at a time, while each of the ILP solutions are determined separately (and our results were determined within similar computation times for heuristic and ILP solutions). The fact that our heuristic returns solutions of similar quality and which follow the shape of the sub-fronts returned by the ILP solutions indicates that the heuristic solution process is powerful and returns solutions of very strong quality. When analyzing the Pareto fronts in Figure 24, the trade-offs between improving accessibility to the same services during different seasons reveals interesting characteristics. The trade-offs in the objective function values are small – indicating that gains in accessibility in one season scale favorably to gains in another season. This is particularly noteworthy for the comparison of health facilities in Figure 24(a), in which the improvements in health access is highly linear, possibly due to the wider spatial distribution and higher absolute number of health facilities, which reduces the occurrence of relatively extreme accessibility gains for any specific road sequence. Such observations may lead to improved decision-making processes. In this case, it may be possible to eliminate the consideration of seasons in health care access. In Figures 25 and 26, the heuristic solutions are plotted with respect to each individual accessibility objective during dry and monsoon season, respectively. The horizontal axis represents each solution number from 1 to 5,000, which were sorted to increase in cost. Therefore, as may be seen in the figures, the budget levels intervals along the solution numbers is non-uniform and displayed by vertical lines where the solutions crossed a budget interval value. The reason for the non-uniformity in cost distribution lies with the heuristic solution methodology, which appears to struggle to find more solutions between 9 and 12 billion NPR – yet manages to obtain solutions of high quality despite the lower distribution. Such restrictions may, nevertheless, be improved in future versions of the heuristic. Also displayed in the figures are the ILP solutions determined at each budget level – these are the single-objective solutions determined with respect to the single accessibility objective in each plot. Finally, the accessibility improvement achieved by the actual roads selected by the Karnali provincial government are also displayed in each plot in red – since their roads plan was selected at a budget of 13 billion NPR, it is located along the corresponding budget line at the right of each plot. Comparing the heuristic solution quality to that of the ILP once more reveals that the heuristic solutions return high-quality solutions that nearly match the quality of the ILP solutions. An example of interesting observations that may be made when analyzing these plots for decision-making purposes is when the solutions for health facilities and municipal headquarters are compared in the dry season. For health facilities, between 9 and 12 billion NPR, a number of solutions are concentrated “lower” along the vertical axis. This indicates that a number of solutions do not perform well with respect to health facilities in this budget range, and this is clarified when examining the municipal headquarters solutions. Here, these solutions are clustered to perform well with respect to municipal access – solutions between 9 and 12 billion NPR are therefore generally more beneficial to municipal access than health facilities, while the improvements to financial institutions in 30 (a) (b) (c) (d) Figure 23: (a)–(c) Similar to Figure 22, objective function values (in person-hours) of the optimization-determined solutions, plotted against service accessibility improvements during the monsoon season. A legend is provided in (d) for the images in (a)–(c). 31 (a) (b) (c) (d) Figure 24: Objective function values (in person-hours) of the 5000 heuristic solutions and selected ILP solutions, plotted against values achieved with respect to different seasons, per service type. 32 this budget range is more universal during the dry season. Performing a similar analysis in the same budget range for the monsoon season in Figure 26, it is interesting to observe that financial institutions are more sensitive to the effects of the monsoon season, as these solutions now exhibit similar clustering in the “lower” objective function value range which was not observed for the dry season. Besides the implications for road building solutions, these findings suggest to decision-makers that seasonal accessibility may impact financial inclusion or access to finance policy objectives. Continuing the analysis of Figures 25 and 26, conspicuous improvement in solution quality is observed when comparing the heuristic solutions to the actual roads plan selected by the Karnali provincial govern- ment. The maximum accessibility improvements for solutions at or near the maximum 13 billion NPR budget are approximately 3, 3, and 1.5 times those yielded by the actual selected roads for health facilities, finan- cial institutions, and municipal headquarters, respectively. Even at half the budget the mean accessibility improvement is double that of the actual roads selected for health services and financial institutions and roughly equal for municipal headquarters. Assessing solution accessibility results in the aggregate reveals important differences in the potential accessibility gains to each type of service. Figure 27(a) plots the mean per NPR gain in accessibility for each service in each season across all 5,000 approximately Pareto-optimal heuristic solutions at different budget levels. Across all solutions, it is observed that mean accessibility gain per NPR is highest for financial institutions and weakest for health facilities. However, this likely reflects the over-concentration of financial institutions in Karnali’s major towns where most profitable businesses reside and the comparative dispersion of public services like health facilities and municipal headquarters. To account for these spatial patterns we compute each service’s nearest neighbor index, a spatial statistic that describes the distribution of a point feature set relative to a completely spatially random process as clustered (approaching 0), random (1), or regular (exceeding 1) (Cover and Hart, 1967), using a common geographic extent. We multiply each service’s index by its mean accessibility figure to establish a dispersion-adjusted accessibility index, with the results presented in Figure 27(b). So adjusted, health facilities retain the lowest accessibility gain per NPR while municipal headquarters gain two times more efficiently than other services. However, the continued upward trend in accessibility per NPR at the maximum budget levels in Figures 25 and 26 suggests that additional budgetary investments would continue significantly improving accessibility to the studied services. The charts in Figures 28 reveal a considerable disparity in investment efficiency amongst proposed roads. These charts visualize the person-hours improvements per 100,000 NPR associated with each road sequence, as generated by the cost-time modeling and analysis routine, with sequences ascending in efficiency from left to right. Sequences are colored from light to dark depending on the frequency of their selection within the 5,000 heuristic solution sets to illustrate the relationship between efficiency and solution selection. Setting aside extreme values outside the 5th and 95th percentiles, the most efficient road during the dry season is approximately 6.6, 18.3, and 9.3 times more efficient than the least efficient road for health facilities, financial institutions, and municipal headquarters respectively. The same pattern holds for the monsoon season, where the most efficient road or sequence is 9.0, 28.1, and 9.0 times more efficient than the least efficient road. Although the most efficient roads are generally the most frequently selected, certain efficient combinations are not selected often and certain less efficient roads are frequently selected. Efficient combinations may be superseded by even more efficient combinations using the same road sequences, and cheap but less efficient roads may be useful complements to more efficient roads and combinations. The last point underlines the utility of our method: judging roads in isolation ignores the significant complementarities that exist between choices. 5. Discussion 5.1. Summary The marked gains in cost-effectiveness from our methods justify their additional complexity, time invest- ment, and expense. Our methods substantially outperform current planning approaches in simultaneously improving accessibility to all services studied. Geospatial cost-time models accurately quantify mean ac- cessibility improvements in a challenging environment and heuristic optimization models efficiently process and highlight the most promising combinations of rural roads based on the provided objectives (accessibility improvements) and constraints (budget). The system provides flexibility to policy makers and analysts by proposing solutions at various budget levels and accommodating additional objectives, accessibility-related 33 Figure 25: Comparing projected accessibility improvements (dry season) from heuristic and weighted sum solutions to the actual road plan endorsed by Karnali’s provincial government in the draft Provincial Transport Master Plan. 34 Figure 26: Comparing projected accessibility improvements (monsoon season) from heuristic and weighted sum solutions to the actual road plan endorsed by Karnali’s provincial government in the draft Provincial Transport Master Plan. 35 (a) (b) Figure 27: Adjusting per NPR gains in accessibility by the relative dispersion of different services (Nearest Neighbor Index (Cover and Hart, 1967)) shows a more complex picture. Municipal headquarters gain the most per NPR relative to their high level of dispersion. or otherwise, with low marginal effort. The end result is a workflow better equipped than existing tools to maximize the human development improvements from rural roads construction in mountainous environments. Table 3 provides an overview of the effort (measured in days) required to arrive at the results presented in the previous section – for a data set of 1,566 raster layers covering Karnali’s 27,984 km2 extent. Larger areas and/or additional service accessibility criteria (added as additional optimization objectives) would increase the computation times for geospatial and optimization analysis, and the interpretation and visualization of results in such instances may similarly require more time. 5.2. Expert feedback Expert feedback was sought from World Bank economists, geographers, and transportation specialists in a series of presentations and one-on-one meetings at the World Bank’s headquarters in February-March 2020. The reception was positive, with strong appreciation for the marked improvement in optimization-based results over the actual roads selected by the Karnali provincial government. Critical feedback focused on the need to incorporate political economy considerations into our analysis to better reflect decision-makers’ real-world constraints. The actual roads selected likely incorporated unobserved policy or social objectives, such as connecting important towns like Simikot, better integrating specific ethnic groups, boosting tourist access, or facilitating military movements to border regions. To accommodate this feedback we re-evaluated our results in terms of a seventh political objective. Some political considerations are context-specific and cannot be predicted ahead of time but the political imperative to ensure equitable investment across a province, “spreading the money out”, is universal. Therefore we re-evaluated our optimization results to identify solutions maximizing political dispersion. Political dispersion was measured in terms of the number of municipalities holding at least five kilometers of one road sequence in a solution set. For example, two example heuristic solutions with large numbers of municipalities touched (solutions 3,449 and 4,845 of 5,000) are visualized in Figure 29. These reach many more municipalities than the actual roads selected in Karnali’s draft Provincial Transport Master Plan (which reaches 20, as illustrated in Figure 30), while solution 4,845 reaches 46 of a maximum 52 possible municipalities (from all the potential roads) at roughly two-thirds the maximum cost of 18.9 billion NPR (to construct all proposed roads). 5.3. Advances in geospatial methods The methodology described offers several improvements on existing practices for prioritizing rural road construction in mountainous areas with heavy seasonal changes in accessibility. Measuring roadwork plans in terms of the cost-effectiveness of improving access to specific services links infrastructure spending more 36 Figure 28: Comparing the efficiency of individual roads sequences to their frequency of selection. 37 (a) (b) Figure 29: Visualizing accessibility improvements from solutions at (a) the maximum budget level, and (b) 8.7 billion NPR budget range. 38 Table 3: Level of effort required for the different processes that were followed to arrive at the results presented in this paper. Step Level of Effort Accessibility modeling 8 days Collecting, preparing data layers 2 days 1 day computation time Building friction surface (with model already prepared) Building accessibility layers 3 days computation time Person-hours analysis 2 days computation time Optimization 2 days 1 day Heuristics (with heuristic already fully developed) 1 day ILP (running separate weighted-sums in multiple CPLEX instances simultaneously) Analysis 13-20 days Interpreting results 3-5 days Creating charts, maps, and other visuals 5-10 days Modifying analysis with additional or 5 days updated parameters Figure 30: Visualizing accessibility improvements from the actual roads plan endorsed by Karnali’s provincial government in the draft Provincial Transport Master Plan 39 explicitly to development objectives than prevailing Nepali practitioner models of infrastructure project eval- uation, which typically simplify around administrative aggregates (i.e. x schools per ward), focus exclusively on one specific service type, or disregard services entirely. More rigorous cost-benefit evaluation methods may similarly disregard services, despite their developmental importance, while presupposing data inputs (e.g. international roughness index, traffic flows per segment) that are prohibitively expensive to collect ac- curately at scale for large numbers of rural roads. Finally, both approaches either ignore, simplify, or require painstaking effort to measure the effects of seasonal changes on transport speeds and overlaps in proposed road accessibility improvements. Vector-based geospatial models offer comparable scalability to cost-time models but fail to manage the inherent complexity of rural, mountainous contexts. Where network models lack geospatial data on walking paths (as is common in Nepal) they must simplify considerably; in Nepal such simplifications are likely to route travelers over difficult or impassable terrain. By contrast, raster-based models explicitly record travel times across every surface, allowing accurate modeling of walking via least-cost paths to the nearest destination or road. Modeling walking also enables walking-specific speed reductions at high altitudes. Furthermore, raster-based travel models can measure improvements in accessibility more accurately and specifically by pinpointing the precise improvement for households in a grid cell via interaction with raster- based population models. This is advantageous in Nepal’s sparsely populated hills where administrative aggregates can mask extremes. Beyond making the analysis more accurate, this specificity allows policy makers to compare roads that best reduce travel times for the most remote populations to roads that improve travel times for the greatest overall population, weighting them according to policy objectives if so desired. This research’s final contribution to geospatial analysis methods is to prove the utility and feasibility of producing and applying cost-time models at scale. Advances in parallel processing and computing power increasingly enable batch production of complex, multi-modal, multi-scenario (e.g. seasonal road speeds) geospatial accessibility models in reasonable timeframes. The advent of affordable cloud computing solutions places such methods within the reach of analysts in Nepal and similar contexts where relatively expensive computing power and internet speeds have traditionally constrained analysis possibilities. Producing multiple accessibility models for each proposed roadway has several advantages: a) it eases the modeling of interactions (overlapping improvements) between roads; b) the analysis can scale to exponentially greater numbers of roads, services, and seasonal scenarios given sufficient geospatial data and computing power; and c) because individual roads are discrete analytical elements they easily fit within the solution sets of optimization models. 5.4. Advances in optimization methods Several novel optimization approaches were followed in this research. Most notable is the consideration of a budget span – instead of a fixed maximum value – which was achieved by novel variable-length chromosome representations of candidate roads plans in the heuristic. Furthermore, the heuristic is capable of determining thousands of solutions which exhibit optimality with respect to multiple conflicting objectives, in reasonable computation times. Seven objectives were considered by the heuristic here, which is an unusually large number even when compared to traditional methods in multi-objective optimization. A significant benefit to obtaining such a large number of solutions that are optimal with respect to many objectives is that it provides decision-makers with rapid analysis capabilities. By determining a large number of solutions once-off with respect to many objectives, the behavior of solutions may later be examined with respect to only selected objectives and criteria. From the set of solutions determined by the heuristic, sub- fronts may simply be extracted (as was performed in Section 4), instead of having to re-perform heuristic or ILP optimization with respect to the selected objectives. Such capabilities may provide decision makers with rapid insight during the analysis process, and can aid in finalizing a budget and most significant objectives – which may then be followed by a final round of optimization to “fine-tune” solutions to these specifications. Other benefits to having a large number of solutions (which were discussed in Section 4), include added insight into a) the relative linearity or concavity of a specific service’s accessibility improvements versus cost increases, b) the trade-offs that exist between accessibility services and seasons, c) the roads’ investment efficiency and selection frequency, and d) the relative improvements to accessibility to different services in different seasons yield at various investment budget levels. Such analyses of large numbers of solutions have not been performed in the literature – most likely because methods capable of obtaining such solutions have not existed until now. 40 The analysis of road sequences’ gains, efficiency, and importance across solutions may help policy makers identify the key roadworks in potential plans. These strategies may further help compensate for the delays and inefficiencies inherent to Nepal’s capital infrastructure planning bureaucracy, by focusing policy makers’ attention on the roadworks most critical to success within an overall optimized portfolio. 5.5. Limitations We must note certain limitations to our analysis. The most significant is that due to resource and data constraints, we present the accessibility gains as a simultaneous phenomenon instead of staggered over the provincial government’s 5-year road construction window. A deeper analysis would specify a yearly budgetary max within the overall budget max and stagger the solution roads across the time period, with gains appropriately discounted. Additionally, we employed flat per-kilometer road costs per geographic area due to lack of data, and to better align and compare with Nepali practitioners and academics (Shrestha, 2018); more specific road cost data would improve our method’s accuracy. Finally, this analysis treats population as spatially fixed, while in reality population is dynamic and new road construction can induce population movements; for example, in Nepal new markets often spring up at newly built roads heads. This could create additional unobserved objectives if the government speculatively selected roads to create new market hubs or connect underpopulated remote villages to reduce population outflows. The complexity of the analysis routines described here places them effectively out of reach of even national government planning agencies in Nepal, let alone provincial or municipal ones. Thus government will need to employ universities and private sector companies with specialized expertise to replicate this analysis in new areas. Given that the studied investment planning exercise happens on a 5 year basis this seems a manageable extra expense. As stated above, care was taken to use free global or common national data sets to ease adoption and replication by other parties. A further limitation is the difficulty of communicating our results to policy makers compared to traditional instruments. Welfare-based evaluation methods ease the cognitive burden on decision makers by measuring investment returns in terms of familiar NPR values or easily parsable rates of return. To compensate we have produced thorough explanatory visuals for our methods and results. Nevertheless, correctly interpreting a 3D bubble chart – itself a simplification of 6D results – of person-hours is inherently more challenging than a table of NPV or Internal Rate of Return (IRR) values. This erodes decision makers’ ability and willingness to use our results. Therefore real-world implementations will need to invest additional time and effort in correctly explaining results and guiding policy makers to use them appropriately. 5.6. Potential improvements and future research opportunities Given the complexity of the multi-disciplinary analysis methods undertaken we employed a relatively straightforward series of inputs, metrics, and decision criteria for roads investment. Having proven the value of this work, we have several suggestions for enlarging its focus, improving its specificity, and aligning it better with traditional planning methods. Additional georeferenced services would nuance the analysis, each one representing two (dry and monsoon season) additional objectives. For instance, the analysis could incorporate accessibility to schools, water points, agricultural extension services, or other important services, even at a village or household level. A variation on this theme is to differentiate qualitatively between services, for instance by separating hospitals from health posts or major from minor financial institutions and then re-running the analysis with each as an objective. A more involved modification would involve calculating cumulative person-hours within service-specific constrained catchments that reflect travel time limits on the utility of services, both generally and for specific demographic groups. For instance, in an August 2018 interview RAP staff described how individuals might travel 8 hours for important health problems but never to bring produce to weekly markets and Nepalese women are typically warier of overnight travel (in excess of 4 hours one way). Accessibility indicators could also be calculated solely for target demographic groups (i.e. pregnant women to maternal health centers). Particularly ambitious would be replacing cumulative person-hours with more data intensive changes in gravity models of service access based on the quality (attractiveness) of services vs. the distance and size of reached populations. Alternatively, the geospatial models could be adapted to provide inputs to a comprehensive cost-benefit analysis. While this paper has explored our methods as an alternative to CBA for planning remote rural roads, 41 given its difficulty and expense in this context, with additional investment in data production and analysis it could instead improve CBA. With information on annual trips per service, time elasticities of service usage, and/or economic benefits of service usage the economic benefits of time saved, increased health, and/or better access to markets and financing could be quantified and discounted within a time-bound cost-benefit analysis. Sensitivity tests could be conducted on optimization results to verify the stability of preferred solutions to cost overruns or accessibility miscalculations. Recommending optimal solutions with the greatest degree of stability across reasonable cost overruns would help guard against the construction industry’s notoriously poor cost control. Similarly, the optimization routines could offset Nepal’s budgetary and planning inefficiencies by recommending which core road sequences to finish first. 6. Conclusion This research has demonstrated the utility of a multidisciplinary approach to arrive at meaningful and substantially improved roadworks plans in a mountainous and data-poor environment characterized by walk- ing as a primary transport modality. In doing so we provide an alternative to traditional planning tools for policy makers in Nepal and other similar contexts. This research uses a cost-time model to measure travel time distances as it credibly manages walking speeds and off-road travel in Nepal’s rugged terrain, is computable with free or Nepali government-standard data sets, and returns intuitive travel times. While cost-time models are standard geospatial analysis tools, we were not able to find prior examples of deriving accessibility metrics for proposed construction via the batch production and comparison of road-specific cost-time models, let alone applications of this model in mountainous areas with walking as a primary travel modality. Facility optimization models are similarly standard tools within the field of operations research, but to our knowledge none have previously incorporated the outputs of cost-time models, much less to analyze access improvements to rural roads. The optimization methods followed in this paper introduced a number of novel approaches that improve on previous roads optimization work and provide more intricate and insightful decision-making possibilities. Thus this research represents a methodological advance within both disciplines and the Nepali country context. To enable replication and adoption within Nepal we have employed open data and developed open source tools where possible; where not we have employed governmental data sources and industry standard software packages common in Nepal. As a result our models provide more precise estimates of populations’ accessibility to specific services in Nepal at a lower cost and greater scalability than available academic or practitioner alternatives. Acknowledgments We are grateful to the UK Department for International Development in Nepal for funding this work via its Evidence for Development program. Hiroki Uematsu’s belief in a more data-driven approach to roads building in Nepal inspired this research and his guidance and review helped ensure its quality. Nethra Palaniswamy and Nandini Krishnan arranged crucial feedback sessions in Washington DC; the participants to these sessions are too many to name here, but their inputs greatly improved the “real world” applicability of this research. Benjamin Stewart, Tom Gertin, and Andres Chamorro of the World Bank’s Geospatial Operations Support Team and Walker Kosmidou-Bradley of the Poverty GP helped us validate and improve the technical modeling. Jonas Parby and Drs. Jagat Shresthra and Seth Appiah-Opoku generously reviewed the initial draft of this paper. The counsel, data, and documents provided by Diva Malla, Arjun Poudel, Michael Green, Kirsteen Merrilees, and Philippa Jeffries from IMC Worldwide’s Rural Access Programme 3 (RAP3) were crucial for enabling this research. References Adukia, A., Asher, S., Novosad, P., 2020. Educational Investment Responses to Economic Opportunity: Evidence from Indian Road Construction. American Economic Journal: Applied Economics 12, 348–376. doi:10.1257/app.20180036. 42 om, A., Pilesj¨ Ahlstr¨ o, P., Lindberg, J., 2011. Improved accessibility modeling and its relation to poverty – A case study in Southern Sri Lanka. Habitat International 35, 316–326. doi:10.1016/j.habitatint.2010. 11.002. Amin, S.H., Zhang, G., 2013. A multi-objective facility location model for closed-loop supply chain network under uncertain demand and return. Applied Mathematical Modelling 37, 4165–4176. doi:10.1016/j. apm.2012.09.039. Arora, J.S., 2011. Introduction to optimum design. 3rd ed ed., Academic Press, Boston, MA. Asher, S., Nagpal, K., Novosad, P., 2018. The Cost of Distance: Geography and Gov- ernance in Rural India. Working Paper. World Bank. Washington, DC. Available from https://caligari.dartmouth.edu/˜knagpal/KaranNagpal JMP.pdf. Banerjee, A., Duflo, E., Qian, N., 2012. On the Road: Access to Transportation Infrastructure and Eco- nomic Growth in China. Working Paper 17897. National Bureau of Economic Research. Available from http://www.nber.org/papers/w17897. Banick, R.S., Kawasoe, Y., 2019. Measuring Inequality of Access : Modeling Phys- ical Remoteness in Nepal. Report 140353. The World Bank. Available from http://documents.worldbank.org/curated/en/605991565195559324/Measuring-Inequality-of-Access- Modeling-Physical-Remoteness-in-Nepal. Belli, P., Anderson, J.R., Tan, J.P., Barnum, H.N., Dixon, J.A., 2001. Economic Analysis of Investment Operations: Analytical Tools and Practical Applications. The World Bank, Washington, DC. doi:10. 1596/0-8213-4850-7. Bikas Udhyami, 2018. A Study Into Development Data in Nepal. Report study. Bikas Ud- hyami. Lalitpur, Nepal. Available from https://www.nepalindata.com/media/resources/items/19/ bA Study Into Development Data in Nepal.pdf. Binswanger, H.P., Khandker, S.R., Rosenzweig, M.R., 1993. How infrastructure and financial institutions affect agricultural output and investment in India. Journal of Development Economics 41, 337–366. doi:10. 1016/0304-3878(93)90062-R. Bird, K., Hulme, D., Shepherd, A., Moore, K., 2002. Chronic Poverty and Remote Rural Areas. SSRN Scholarly Paper ID 1754490. Social Science Research Network. Rochester, NY. Available from https://papers.ssrn.com/abstract=1754490. Booth, D., Hanmer, L., Lovell, E., 2000. Poverty and Transport. Report. Overseas Development Institute. Available from https://www.odi.org/sites/odi.org.uk/files/odi-assets/publications-opinion-files/3554.pdf. Branke, J., Scheckenbach, B., Stein, M., Deb, K., Schmeck, H., 2009. Portfolio optimization with an envelope- based multi-objective evolutionary algorithm. European Journal of Operational Research 199, 684–693. doi:10.1016/j.ejor.2008.01.054. Bryceson, D.F., Bradbury, A., Bradbury, T., 2008. Roads to Poverty Reduction? Exploring Rural Roads’ Im- pact on Mobility in Africa and Asia. Development Policy Review 26, 459–482. doi:10.1111/j.1467-7679. 2008.00418.x. Casaburi, L., Glennerster, R., Suri, T., 2013. Rural Roads and Intermediated Trade: Regression Disconti- nuity Evidence from Sierra Leone. SSRN Scholarly Paper ID 2161643. Social Science Research Network. Rochester, NY. Available from https://papers.ssrn.com/abstract=2161643. Chapagain, B., Limbu, B., 2018. Rural roads: Repaired in winter but unfit for monsoon. Available from https://myrepublica.nagariknetwork.com/news/47052/. Cheshmehgaz, H.R., Haron, H., Sharifi, A., 2015. The review of multiple evolutionary searches and multi-objective evolutionary algorithms. Artificial Intelligence Review 43, 311–343. doi:10.1007/ s10462-012-9378-3. 43 Christiaensen, L., Demery, L., Paternostro, S., 2003. Reforms, remoteness and risk in Africa: Under- standing inequality and poverty during the 1990s. Working Paper 2003/70. The United Nations Uni- versity World Institute for Development Economics Research (UNU-WIDER). Helsinki. Available from https://www.wider.unu.edu/sites/default/files/dp2003-070.pdf. Coburn, B., 2020. Nepal’s Road-Building Spree Pushes into the Heart of the Himalayas. Yale Environment 360 Available from https://e360.yale.edu/features/paving-the-himalayas-a-road-building-spree-rolls-over- nepal. Cohon, J.L., 1978. Multiobjective Programming and Planning. Academic Press, New York (NY). Cover, T., Hart, P., 1967. Nearest neighbor pattern classification. IEEE Transactions on Information Theory 13, 21–27. doi:10.1109/TIT.1967.1053964. Cramer, N.L., 1985. A Representation for the Adaptive Generation of Simple Sequential Programs, in: Proceedings of the 1st International Conference on Genetic Algorithms, L. Erlbaum Associates Inc., USA. pp. 183–187. Das, I., Dennis, J.E., 1997. A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems. Structural optimization 14, 63–69. doi:10.1007/ BF01197559. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6, 182–197. doi:10.1109/4235.996017. Densham, P.J., Rushton, G., 1992. A More Efficient Heuristic for Solving Large P-Median Problems. Papers in Regional Science 71, 307–329. doi:10.1111/j.1435-5597.1992.tb01849.x. Department of Roads, 2013. Nepal Road Standard 2070. Report. Government of Nepal, Ministry of Physical Infrastructure & Transport, Department of Roads. Kathmandu, Nepal. Available from https://tid.p3.gov.np/wp-content/uploads/2019/06/Nepal-Road-Standard-2070.pdf. Devkota, B., Dudycha, D., Andrey, J., 2012. Planning for non-motorized travel in rural Nepal: a role for geographic information systems. Journal of Transport Geography 24, 282–291. doi:10.1016/j.jtrangeo. 2012.03.007. Edmonds, G., 1998. Wasted time: the price of poor access. Rural Access Technical Papers (RATP) No. 3. Development Policies Department, International Labour Office Geneva. Geneva, Switzerland. Available from https://www.ilo.org/wcmsp5/groups/public/—ed emp/—emp policy/— invest/documents/publication/wcms 142663.pdf. Fafchamps, M., Shilpi, F., 2003. The spatial division of labour in Nepal. Journal of Development Studies 39, 23–66. doi:10.1080/00220380312331293577. Fonseca, C.M., Fleming, P.J., 1993. Genetic algorithms for multiobjective optimization: Formulation, dis- cussion and generalization, in: Proceedings of the Fifth International Conference on Genetic Algorithms, pp. 416–423. Georgiadou, P.S., Papazoglou, I.A., Kiranoudis, C.T., Markatos, N.C., 2010. Multi-objective evolutionary emergency response optimization for major accidents. Journal of Hazardous Materials 178, 792–803. doi:10. 1016/j.jhazmat.2010.02.010. Gilbert, M., Golding, N., Zhou, H., Wint, G.R.W., Robinson, T.P., Tatem, A.J., Lai, S., Zhou, S., Jiang, H., Guo, D., Huang, Z., Messina, J.P., Xiao, X., Linard, C., Van Boeckel, T.P., Martin, V., Bhatt, S., Gething, P.W., Farrar, J.J., Hay, S.I., Yu, H., 2014. Predicting the risk of avian influenza A H7N9 infection in live-poultry markets across Asia. Nature Communications 5, 1–7. doi:10.1038/ncomms5116. Harvest Choice and International Food Policy Research Institute (IFPRI), 2016. Travel Time to Markets in Africa South of the Sahara. Dataset. Harvard Dataverse. Available from https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/YKDWJD. 44 Hettige, H., 2006. When do rural roads benefit the poor and how? An in-depth analysis based on case studies. Operations Evaluation Dept., Asian Development Bank, Metro Manila, Philippines. Available from https://www.adb.org/sites/default/files/publication/29406/when-rural-roads-benefit-poor.pdf. Heyns, A.M., 2016. A multi-objective approach towards geospatial facility location. PhD Thesis. Stellenbosch University. Heyns, A.M., du Plessis, W., Kosch, M., Hough, G., 2019. Optimisation of tower site locations for camera- based wildfire detection systems. International Journal of Wildland Fire 28, 651–665. Heyns, A.M., van Vuuren, J.H., 2016. A multi-resolution approach towards point-based multi-objective geospatial facility location. Computers, Environment and Urban Systems 57, 80–92. doi:10.1016/j. compenvurbsys.2016.01.007. Heyns, A.M., van Vuuren, J.H., 2018. Multi-type, multi-zone facility location. Geographical Analysis 50, 3–31. doi:10.1111/gean.12131. Hirvonen, K., Hoddinott, J., Minten, B., Stifel, D., 2017. Children’s Diets, Nutrition Knowledge, and Access to Markets. World Development 95, 303–315. doi:10.1016/j.worlddev.2017.02.031. Huber, S., 2015. Accessibility of peripheral areas in Nepal: the role of infrastructure development and environmental constraints as limiting factors. Research Paper. TU Dresden. Berlin. Available from http://rgdoi.net/10.13140/RG.2.1.3241.0962. ICIMOD, 2013. Land cover of Nepal 2010 [Data set]. Available from https://doi.org/10.26066/rds.9224. Jacoby, H.G., 2000. Access to Markets and the Benefits of Rural Roads. The Economic Journal 110, 713–737. doi:10.1111/1468-0297.00562. Jafino, B.A., Kwakkel, J., Verbraeck, A., 2020. Transport network criticality metrics: a comparative analysis and a guideline for selection. Transport Reviews 40, 241–264. doi:10.1080/01441647.2019.1703843. Jalan, J., Ravallion, M., 2002. Geographic poverty traps? A micro model of consumption growth in rural China. Journal of Applied Econometrics 17, 329–346. doi:10.1002/jae.645. Jenelius, E., 2009. Network structure and travel patterns: explaining the geographical disparities of road network vulnerability. Journal of Transport Geography 17, 234–244. doi:10.1016/j.jtrangeo.2008.06. 002. on Jia, H., Ord´ ˜ez, F., Dessouky, M.M., 2007. Solution approaches for facility location of medical supplies for large-scale emergencies. Computers & Industrial Engineering 52, 257–276. doi:10.1016/j.cie.2006.12. 007. Khan, S.A., Rehman, S., 2013. Iterative non-deterministic algorithms in on-shore wind farm design: A brief survey. Renewable and Sustainable Energy Reviews 19, 370–384. doi:10.1016/j.rser.2012.11.040. Khandker, S., Bakht, Z., Koolwal, G., 2009. The Poverty Impact of Rural Roads: Evidence from Bangladesh. Economic Development and Cultural Change 57, 685–722. doi:10.1086/598765. Khandker, S.R., Koolwal, B.G., 2011. Estimating the Long-Term Impacts of Rural Roads: A Dynamic Panel Approach. Policy Research Working Paper. World Bank. Available from https://elibrary.worldbank.org/doi/pdf/10.1596/1813-9450-5867. Kim, K., Murray, A.T., Xiao, N., 2008. A multiobjective evolutionary algorithm for surveillance sensor placement. Environment and Planning B: Planning and Design 35, 935–948. doi:10.1068/b33139. Knowles, J.D., Thiele, L., Zitzler, E., 2006. A tutorial on the performance assessment of stochastic multiob- jective optimizers. Technical Report 214. Computer Engineering and Networks Laboratory (TIK). ETH Zurich. 45 Konstantinidis, A., Yang, K., Zhang, Q., Zeinalipour-Yazti, D., 2010. A multi-objective evolutionary algo- rithm for the deployment and power assignment problem in wireless sensor networks. Computer Networks 54, 960–976. doi:10.1016/j.comnet.2009.08.010. Kosmidou-Bradley, W., Blankespoor, B., 2019. Measuring Mobility in Afghanistan using Time- Cost Raster Models. Technical Report Methodology Note. World Bank Group. Avail- able from http://documents.worldbank.org/curated/en/854001546846115537/pdf/Measuring-Mobility-in- Afghanistan-Using-Time-Cost-Raster-Models-Methodology-Note.pdf. Kristjanson, P., Radeny, M., Baltenweck, I., Ogutu, J., Notenbaert, A., 2005. Livelihood mapping and poverty correlates at a meso-level in Kenya. Food Policy 30, 568–583. doi:10.1016/j.foodpol.2005.10.002. Kumar, A., Kumar, P., 1999. User-Friendly Model for Planning Rural Roads. Transportation Research Record: Journal of the Transportation Research Board 1652, 31–39. doi:10.3141/1652-05. Kumar, A., Tillotson, H.T., 1991. Planning Model for Rural Roads. Transportation Research Record: Journal of the Transportation Research Board 1291, 171–181. Kumar, P., Jain, S.S., 2000. Optimal rural road network planning for developing countries. Road and Transport Research 9, 51–66. LAHURNIP, IWGIA, 2014. A Study on the Socio-Economic Status of Indigenous Peoples in Nepal. Lawyers’ Association for Human Rights of Nepalese Indigenous Peoples (LAHURNIP), The International Work Group for Indigenous Affairs (IWGIA), Kathmandu, Nepal. Lau, H.C.W., Chan, T.M., Tsui, W.T., Chan, F.T.S., Ho, G.T.S., Choy, K.L., 2009. A fuzzy guided multi- objective evolutionary algorithm model for solving transportation problem. Expert Systems with Applica- tions 36, 8255–8268. doi:10.1016/j.eswa.2008.10.031. Lebo, J., Schelling, D., 2001. Design and appraisal of rural transport infrastructure : ensuring basic access for rural communities. Technical Report WTP496. The World Bank. Washington, DC. Available from http://documents.worldbank.org/curated/en/227731468184131693/pdf/multi0page.pdf. Ma, C., Ma, C., Ye, Q., He, R., Song, J., 2014. An Improved Genetic Algorithm for the Large-Scale Rural Highway Network Layout. Mathematical Problems in Engineering 2014, 1–6. doi:10.1155/2014/267851. Machairas, V., Tsangrassoulis, A., Axarli, K., 2014. Algorithms for optimization of building design: A review. Renewable and Sustainable Energy Reviews 31, 101–112. doi:10.1016/j.rser.2013.11.036. Maji, A., Jha, M.K., 2017. Multi-Objective Evolutionary Algorithm framework for highway route planning with case study. Advances in Transportation Studies 41, 51–72. Marler, R., Arora, J., 2004. Survey of multi-objective optimization methods for engineering. Structural and Multidisciplinary Optimization 26, 369–395. doi:10.1007/s00158-003-0368-6. Marler, R.T., Arora, J.S., 2010. The weighted sum method for multi-objective optimization: new insights. Structural and Multidisciplinary Optimization 41, 853–862. doi:10.1007/s00158-009-0460-7. Mathew, B.S., Isaac, K.P., 2014. Optimisation of maintenance strategy for rural road network using genetic algorithm. International Journal of Pavement Engineering 15, 352–360. doi:10.1080/10298436.2013. 806807. Mavrotas, G., 2009. Effective implementation of the epsilon-constraint method in Multi-Objective Mathe- matical Programming problems. Applied Mathematics and Computation 213, 455–465. Mavrotas, G., Florios, K., 2013. An improved version of the augmented epsilon-constraint method (AUG- MECON2) for finding the exact pareto set in multi-objective integer programming problems. Applied Mathematics and Computation 219, 9652–9669. doi:10.1016/j.amc.2013.03.002. 46 Ministry of Education, Nepal and UNESCO Office Kathmandu, 2015. Education for All Na- tional Review Report: 2001-2015 (Nepal). National Review Report. Nepal Ministry of Education and United Nations Children’s Education Fund (UNICEF). Kathmandu, Nepal. Available from https://unesdoc.unesco.org/ark:/48223/pf0000232769. Ministry of Finance, 2015. Fiscal Year 2014/2015. Economic Survey. Government of Nepal Ministry of Fi- nance. Available from https://mof.gov.np/uploads/document/file/Final Economic Survey 2071-72 English (Final) 20150716082638.pdf. Murawski, L., Church, R.L., 2009. Improving accessibility to rural health services: The maximal covering network improvement problem. Socio-Economic Planning Sciences 43, 102–110. doi:10.1016/j.seps. 2008.02.012. Murray, A.T., Kim, K., Davis, J.W., Machiraju, R., Parent, R., 2007. Coverage optimization to sup- port security monitoring. Computers, Environment and Urban Systems 31, 133–147. doi:10.1016/j. compenvurbsys.2006.06.002. Myers, S.C., 1974. Interactions of Corporate Financing and Investment Decisions—Implications for Capital Budgeting. The Journal of Finance 29, 1–25. doi:10.1111/j.1540-6261.1974.tb00021.x. National Planning Commission and Central Bureau of Statistics, 2020. Provincial SDG Dashboard. Available from https://dataviz.worldbank.org/views/ProvincialSDGDashboardv8/ Province?iframeSizedToWindow=true&:embed=y&:showAppBanner=false&:display count=no &:showVizHome=no. Nelson, A., Chomitz, K.M., 2011. Effectiveness of Strict vs. Multiple Use Protected Areas in Reducing Tropical Forest Fires: A Global Analysis Using Matching Methods. PLoS ONE 6, e22722. doi:10.1371/ journal.pone.0022722. Owen, S.H., Daskin, M.S., 1998. Strategic facility location: A review. European Journal of Operational Research 111, 423–447. aez, A., 2004. Network Accessibility and the Spatial Distribution of Economic Activity in Eastern Asia:. P´ Urban Studies doi:10.1080/0042098042000268429. Pantha, B.R., Yatabe, R., Bhandary, N.P., 2010. GIS-based highway maintenance prioritization model: an integrated approach for highway maintenance in Nepal mountains. Journal of Transport Geography 18, 426–433. doi:10.1016/j.jtrangeo.2009.06.016. Porter, G., 1997. Mobility and Inequality in Rural Nigeria: The Case of Off-Road Communities. Tijdschrift voor economische en sociale geografie 88, 65–76. doi:10.1111/j.1467-9663.1997.tb01584.x. Porter, G., 2002a. Improving mobility and access for the off-road rural poor through Intermediate Means of Transport. World transport policy and practice. 8, 6–19. Porter, G., 2002b. Living in a Walking World: Rural Mobility and Social Equity Issues in Sub-Saharan Africa. World Development 30, 285–300. doi:10.1016/S0305-750X(01)00106-1. Porter, G., 2007. Transport planning in sub-Saharan Africa. Progress in Development Studies 7, 251–257. doi:10.1177/146499340700700305. Purshouse, R.C., Fleming, P.J., 2003. Evolutionary many-objective optimisation: an exploratory analysis, in: The 2003 Congress on Evolutionary Computation, 2003, pp. 2066–2073. doi:10.1109/CEC.2003.1299927. Rao, K.M.L., Jayasree, K., 2003. Rural Infrastructure Planning with emphasis on road network connectivity by Coplanar Concurrent Theory, in: Proceedings of the Map India Conference 2003, Map India Conference 2003, New Delhi, India. p. 12. ReVelle, C., Eiselt, H., 2005. Location analysis: A synthesis and survey. European Journal of Operational Research 165, 1–19. doi:10.1016/j.ejor.2003.11.032. 47 Ripon, K.N., Tsang, C.H., Kwong, S., Ip, M.K., 2006. Multi-objective evolutionary clustering using variable- length real jumping genes genetic algorithm, in: Pattern Recognition, 2006. ICPR 2006. 18th International Conference on, IEEE. pp. 1200–1203. Rodrigue, J.P., 2016. The Role of Transport and Communication Infrastructure in Realising Development Outcomes, in: Grugel, J., Hammett, D. (Eds.), The Palgrave Handbook of International Development. Palgrave Macmillan UK, London, pp. 595–614. doi:10.1057/978-1-137-42724-3_33. Rodrigue, J.P., Comtois, C., Slack, B., 2016. The Geography of Transport Systems. 4 ed., Routledge. doi:10.4324/9781315618159. Rural Access Programme (RAP3), 2019. Provincial Transport Master Plan Guidelines. Internal document. Rural Access Programme (RAP3). Kathmandu, Nepal. Sabogal, O., Escobar, D., Oviedo, D., 2018. Making Accessibility Visible: Visualizing Spatial Accessibil- ity Through Multi-Dimensional Scaling Model. Modern Applied Science 12, 70–85. doi:10.5539/mas. v12n6p70. Sapkota, J.B., 2018. Access to infrastructure and human well-being: evidence from rural Nepal. Development in Practice 28, 182–194. doi:10.1080/09614524.2018.1424802. Scaparra, M.P., Church, R.L., 2005. A GRASP and Path Relinking Heuristic for Rural Road Network Development. Journal of Heuristics 11, 89–108. doi:10.1007/s10732-005-7000-4. Schmitz, C., Biewald, A., Lotze-Campen, H., Popp, A., Dietrich, J.P., Bodirsky, B., Krause, M., Weindl, I., 2012. Trading more food: Implications for land use, greenhouse gas emissions, and the food system. Global Environmental Change 22, 189–209. doi:10.1016/j.gloenvcha.2011.09.013. Shrestha, J.K., 2018. Rural Road Network Decision Model for Hilly Regions of Nepal. Journal of Advanced College of Engineering and Management 4, 51–64. doi:10.3126/jacem.v4i0.23178. Shrestha, J.K., Benta, A., Lopes, R.B., Lopes, N., 2014. A multi-objective analysis of a rural road network problem in the hilly regions of Nepal. Transportation Research Part A: Policy and Practice 64, 43–53. doi:10.1016/j.tra.2014.03.005. Shrestha, S.A., 2012. Access to the North-South Roads and Farm Profits in Rural Nepal. Unpublished manuscript. University of Michigan. Ann Arbor, Michigan. Available from https://www.dartmouth.edu/neudc2012/docs/paper 45.pdf. Srikanth, R., George, R., Warsi, N., Prabhu, D., Petry, F., Buckles, B., 1995. A variable-length ge- netic algorithm for clustering and classification. Pattern Recognition Letters 16, 789–800. doi:10.1016/ 0167-8655(95)00043-G. Stewart, T., 2007. The essential multiobjectivity of linear programming. ORiON 23, 1–15. doi:10.5784/ 23-1-43. Stifel, D.C., Minten, B., Dorosh, P., 2003. Transactions Costs and Agricultural Productivity: Implications of Isolation for Rural Poverty in Madagascar. MSSD Discussion Paper No. 56. International Food Policy Research Institute. Available from http://www.ssrn.com/abstract=449220. Van der Tak, H.G., Anandarup, R., 1971. The economic benefits of road transport projects. Technical Report OCP13. The World Bank. Washington, DC. Available from http://documents.worldbank.org/curated/en/983941468766182594/pdf/multi0page.pdf. Teitz, M.B., Bart, P., 1968. Heuristic Methods for Estimating the Generalized Vertex Median of a Weighted Graph. Operations Research 16, 955–961. doi:10.1287/opre.16.5.955. Thapa, G., Shively, G., 2017. Road and Market Access, and Household Food Security in Nepal. Report. World Food Programme. Lalitpur, Nepal. Available from https://www.wfp.org/publications/road-and- market-access-and-household-food-security-nepal. 48 Ting, C.K., Lee, C.N., Chang, H.C., Wu, J.S., 2009. Wireless Heterogeneous Transmitter Placement Using Multiobjective Variable-Length Genetic Algorithm. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 39, 945–958. doi:10.1109/TSMCB.2008.2010951. Tobler, W., 1993. Three Presentations on Geographical Analysis and Modeling: Non- Isotropic Geographic Modeling; Speculations on the Geometry of Geography; and Global Spatial Analysis. Technical Report 93-1. National Center for Geographic Information and Analysis. University of California Santa Barbara. Available from http://ncgia.ucsb.edu/technical-reports/PDF/93-1.pdf. Tong, D., Murray, A., Xiao, N., 2009. Heuristics in spatial analysis: A genetic algorithm for coverage maximization. Annals of the Association of American Geographers 99, 698–711. United Nations Development Programme, 2017. Human Development Report 2016: Human Development for Everyone. Human Development Report, United Nations, New York, NY, USA. Available from https://www.un-ilibrary.org/economic-and-social-development/human-development- report-2016 b6186701-en. Van Wagtendonk, J.W., Benedict, J.M., 1980. Travel Time Variation on Backcountry Trails. Journal of Leisure Research 12, 99–115. doi:10.1177/004728758001900266. Van de Walle, D., Cratty, D., 2002. Impact Evaluation of a Rural Road Re- habilitation Project. Working Paper 44472. World Bank. Washington, DC. Available from http://documents.worldbank.org/curated/en/860291468129595286/pdf/ 444720WP0Dwall10Box334044B01PUBLIC1.pdf. Weiss, D.J., Nelson, A., Gibson, H.S., Temperley, W., Peedell, S., Lieber, A., Hancher, M., Poyart, E., Belchior, S., Fullman, N., Mappin, B., Dalrymple, U., Rozier, J., Lucas, T.C.D., Howes, R.E., Tusting, L.S., Kang, S.Y., Cameron, E., Bisanzio, D., Battle, K.E., Bhatt, S., Gething, P.W., 2018. A global map of travel time to cities to assess inequalities in accessibility in 2015. Nature 553, 333–336. doi:10.1038/nature25181. World Bank, 2016. Strengthening Connectivity in Nepal. Available from https://www.worldbank.org/en/news/feature/2016/02/17/strengthening-rural-connectivity-in-nepal. World Bank, 2017. Climbing Higher: Towards a Middle-Income Nepal. Eco- nomic Memorandum. The World Bank Group. Washington, DC. Available from http://documents.worldbank.org/curated/en/358501495199225866/pdf/115156-CEM-PUBLIC-SAREC- 70p-Country-Economic-Memorandum-19-May-2017.pdf. World Bank, 2018. Implementation and Completion and Results Report on a Credit in the Amount of SDR 38.7 million to the Government of Nepal for a Bridges Improve- ment and Maintenance Program. Report 121710-NP. The World Bank Group. Available from http://documents.worldbank.org/curated/en/190621516371769416/pdf/Nepal-BIMP-P125495-ICR- Main-Document-Approved-1-15-2018-01162018.pdf. World Bank, 2019. NEPAL: Strategic Road Connectivity and Trade Improvement Project. Internal document. World Bank Group. Washington, DC. World Bank, 2020. Cox’s Bazar Growth Diagnostic. Growth Diagnostic. The World Bank. Washington, DC. Xia, J., Curtin, K.M., Huang, J., Wu, D., Xiu, W., Huang, Z., 2019. A carpool matching model with both social and route networks. Computers, Environment and Urban Systems 75, 90–102. doi:10.1016/j. compenvurbsys.2019.01.008. Xiao, N., Bennett, D.A., Armstrong, M.P., 2002. Using evolutionary algorithms to generate alternatives for multiobjective site-search problems. Environment and Planning A 34, 639–656. doi:10.1068/a34109. Yao, J., Zhang, X., Murray, A.T., 2018. Spatial Optimization for Land-use Allocation: Accounting for Sus- tainability Concerns. International Regional Science Review 41, 579–600. doi:10.1177/0160017617728551. 49 Zhu, G., Wang, L., Zhang, P., Song, K., 2016. A Kind of Urban Road Travel Time Forecasting Model with Loop Detectors:. International Journal of Distributed Sensor Networks 2016, 1–11. doi:10.1155/2016/ 9043835. Zitzler, E., Laumanns, M., Bleuler, S., 2004. A Tutorial on Evolutionary Multiobjective Optimization, in: Gandibleux, X., Sevaux, M., S¨ orensen, K., T’kindt, V. (Eds.), Metaheuristics for Multiobjective Optimi- sation, Springer Berlin, Heidelberg. pp. 3–37. Zitzler, E., Thiele, L., Laumanns, M., Fonseca, C., da Fonseca, V., 2003. Performance assessment of multiob- jective optimizers: an analysis and review. IEEE Transactions on Evolutionary Computation 7, 117–132. doi:10.1109/TEVC.2003.810758. 50