Luciano Gervasoni
email: luciano.gervasoni@inria.fr

Marti Bosch
email: marti.bosch@epfl.ch

Serge Fenet
email: ‡serge.fenet@liris.cnrs.fr

Peter Sturm
email: peter.sturm@inria.fr

A framework for evaluating urban land use mix from crowd-sourcing data

Population in urban areas has been increasing at an alarming rate in the last decades. This evidence, together with the rising availability of massive data from cities, has motivated research on sustainable urban development.

In this paper we present a GIS-based land use mix analysis framework to help urban planners to compute indices for mixed uses development, which may be helpful towards developing sustainable cities. Residential and activities land uses are extracted using OpenStreetMap crowd-sourcing data. Kernel density estimation is performed for these land uses, and then used to compute the mixed uses indices. The framework is applied to several cities, analyzing the land use mix output.

I. INTRODUCTION

The number of people living in cities has been increasing considerably since 1950, from 746 million to 3.9 billion in 2014 [START_REF]World urbanization prospects: The 2014 revision, highlights[END_REF]. More than 66% of the world's population are projected to live in urban areas by 2050, against 30% in 1950. Even though the global rural population is expected to only slightly decline by 2050, the continuing population growth and urbanization are thus projected to add 2.5 billion people to the world's urban population by 2050. This situation brings new challenges on how to conceive cities that host such amounts of population in a sustainable way. This sustainability question should address several aspects, ranging from economical to social and environmental matters among others.

In this paper, we focus on the formalization of a measure of mixed use development or land use mix (LUM) in a city, i.e. how different urban land uses are located close enough to each other. Such an urban land use mix has been largely proven to contain beneficial outcomes in terms of sustainability. Indeed, dense and compact cities have been largely discussed in the literature [START_REF] Gordon | Are compact cities a desirable planning goal?[END_REF], [START_REF] Breheny | The compact city and transport energy consumption[END_REF], [START_REF] Burgess | The compact city debate: A global perspective[END_REF], [START_REF] Neuman | The compact city fallacy[END_REF] in terms of sustainable development. While the desired degree of compactness remains an open question, a wide agreement exists on the positive impact of land use mix: it has been proven to be positive not only in terms of sustainable development, but also to contribute to several aspects such as societal, health and public transportation among others.

In order to aid urban planners, we propose a framework to capture land use mix by means of crowd-sourced data. The output is a Geographic Information System (GIS) containing the degree of land use mixture of cities. For this, land uses extraction is performed from the input data. Then, spatial statistics kernels are estimated, and used to define the final mix degree. Even though it is desirable to achieve a degree of land use mix as high as possible, in this study we do not perform numerical interpretations of its desired value in order to contribute to sustainable development, and we focus only on the computation of the spatialized land use mix map.

The manuscript is organized as follows. In Sect. II, an overview of the general context is presented, followed in Sect. III by an analysis of existing indicators. Sect. IV outlines the data used in our framework. The proposed framework is presented in Sect. V. Finally, we present an application to a few cities in Sect. VI, followed by the conclusions and a discussion of future work.

II. CONTEXT A. Motivation

A clear steady migratory pattern from rural to urban areas can be observed through history, becoming more and more important since the nineteenth century. This process, that can be observed in all countries, has lead to a constant increase of urban population: now more than half of the human population lives in urban areas, and cities are consequently facing an ever-increasing population and struggle to continue ensuring essential services and quality of life. This concentration leads to new demands on how to conceive cities in a way that promotes sustainability and efficient resource consumption. Still, urban processes take place as a consequence of different interacting factors, linked between them in such a way that the resulting process is complex.

Given the increasing number of people living in cities, understanding the underlying complexity of these urban patterns is becoming a pressing issue. As a result of the increasing availability of massive amounts of urban data, it is now possible to analyze the ways citizens interact within cities. More importantly, understanding the complex feedback loop linking citizens' and cities' development could allow us to better solve future sustainability questions. Consequently, considerable efforts have been emerging recently, aiming to understand cities' sustainability and increase it by using newly-available data sources.

The present work was motivated by three main observations:

• First, as sustainability is becoming the central common issue in several scientific fields, we advocate that computer science tools, and more specifically data mining and model building applied to real world data, can help to better understand the dynamics of complex socioeconomic systems. These tools will allow a better comprehension and efficient decision making towards improved sustainability and resilience. • Second, several measures of land use mix exist, see [START_REF] Song | Comparing measures of urban land use mix[END_REF] and references therein. Most can be categorized in divisional or integral measures. The former are more expressive than the latter, delivering a single value for each divisional unit. However, all of them are sensitive to the scale used to aggregate data to compute the measure.

In this sense, a desirable measure should be independent of the scale of analysis. • Third, even if crowd-sourcing data is growing with both increasing precision and frequency, missing data will always occur to some degree, impacting the above measures. Hence, statistical tools are vital in this context, in order to better capture the cities' underlying features.

B. Brief historical background

Before the 20th century, mixed use was a natural trend in city development, since scarcity of transportation possibilities imposed geographical proximity constraints on the location of every-day activities. Early in the 1900s in the United States, zoning practices started to assign unique land uses, inducing segregation between residential and activities uses. This has been occurring particularly from the 1910s to the 1950s, where mixed use development was quite infrequent. During this period, segregated development was the norm. In later years, starting from the 60s and 70s, and after having been neglected for several decades, mixed uses started appearing again, bringing altogether its advantages in various aspects to the society.

More recently, since the late 90s, mixed use development has (re-)emerged as a major concept in the context of urban planning. The Congress for the New Urbanism ( [START_REF] Leccese | Charter of the new urbanism[END_REF], http://www.cnu.org) campaigns for "pedestrian-friendly, and mixed-use" neighborhoods, and the Smart Growth Network [START_REF]This is smart growth[END_REF] includes mixing uses as one of its ten principles. The concept of Transit-Oriented Development and several transit agencies also support the provision of a mix of land uses [START_REF] Cervero | Transit-oriented development and joint development in the United States: A literature review[END_REF]. In order to attain this state of highly mixed development, both intensity and diversity of land uses are promoted, together with the integration of the segregated uses.

C. Impacts of mixed use development

In this section, we briefly outline the principal issues and effects related to sustainability within the presence of high/low mixed use development. Both direct and indirect relations between sustainability and mixed uses has been intensively studied in the literature, as quickly described below:

• Urban sprawl: Urban sprawl refers to the process of spatial expansion of the population moving away from central urban areas into sparse, mono-functional and usually car-dependent communities. This suburbanization, i.e. the creation of suburban areas, has been growing with an ever-increasing rate in the last century. As a consequence, low urban land use mix is one of the distinctive characteristics related to this urban sprawl phenomenon (i.e. residential or industrial sprawl).

In the classical literature, the urban sprawl process has been largely linked to negative effects in terms of environment, health, society and economy [START_REF] Squires | Urban sprawl: Causes, consequences, & policy responses[END_REF]. More recently, and in conjunction with the appearance of the concept of sustainable development of cities, urban planners are increasingly taking into account the consequences of sprawling [START_REF] Johnson | Environmental impacts of urban sprawl: a survey of the literature and proposed research agenda[END_REF] -even if the negative consequences of sprawl have been already addressed early in the 70s.

The lack of coordinated land use planning and its negative consequences start to appear as one of the key components of urban sprawl [START_REF] Nelson | Growth management principles and practices[END_REF], [START_REF] Ewing | Characteristics, causes, and effects of sprawl: A literature review[END_REF], [START_REF]Sprawl: The dark side of the American dream[END_REF], [START_REF]State of the cities -1999[END_REF]. • Transportation: In terms of public transportation, sprawling areas with low mixed use development are not sustainable [START_REF] Cervero | Transit-oriented development and joint development in the United States: A literature review[END_REF]. In this direction, those sprawling areas with high dispersion involve the fact that large distances have to be driven for low demands, causing inefficiency in the transportation area. • Health: It was shown that land use diversity has an impact on different types and amounts of physical exercises [START_REF] Christian | How important is the land use mix measure in understanding walking behaviour? Results from the RESIDE study[END_REF]. This means that urban design may directly impact the physical activity of citizens, which is an important health issue nowadays. Thus, a higher mix of uses is associated to more walking trips [START_REF] Christian | How important is the land use mix measure in understanding walking behaviour? Results from the RESIDE study[END_REF] and less obesity [START_REF] Frank | Obesity relationships with community design, physical activity, and time spent in cars[END_REF]. For the latter, it was claimed to be effective as health intervention. • Car dependency: The choice of travel behavior, considering modal choice and distance traveled, has been strongly related to urban land use balance in [START_REF] Cervero | Mixed land-uses and commuting: evidence from the american housing survey[END_REF], [START_REF] Kockelman | Travel behavior as function of accessibility, land use mixing, and land use balance: evidence from San Francisco Bay Area[END_REF]. The presence of nearby commercial land uses was associated with both short commuting distances and low vehicle ownership rates [START_REF] Cervero | Mixed land-uses and commuting: evidence from the american housing survey[END_REF], which is positive in environmental terms since the sources of pollutant emissions are reduced. Studies on designing urban forms that reduce vehicle dependence [START_REF] Newman | Urban design to reduce automobile dependence[END_REF] determined that a lower automobile dependence requires a minimum value of "urban intensity", (i.e. residents and jobs per hectare), which is directly related to the urban land use mix concept.

III. LIMITS OF EXISTING INDICATORS

Efforts on measuring land use mix abound in the literature. A comprehensive review is provided in [START_REF] Song | Comparing measures of urban land use mix[END_REF], whereas each measure has its own strengths and weaknesses. The suitability of a certain measure of urban land use mix relies directly on two aspects. Firstly, it depends on the context of application, where the intended use is the major driver for the measure's behavior. Secondly, the input data characteristics condition the performance of the different measures, in aspects such as the presence of noise among others.

Urban land use mix measures are largely inspired by landscape ecology metrics, or spatial statistics analysis [START_REF] Raines | Description and comparison of geologic maps with FRAGSTATS -a spatial statistics program[END_REF]. Several indices have been proposed, such as the Atkinson index, the Clustering index, the Dissimilarity index, the Exposure index and the Gini index. In addition, Shannon entropy-based metrics were also proposed.

The input data-sets are rarely discussed. They are normally obtained to study a certain region without further questioning. On the one hand, some of these data-sets are private, obstructing the reproducibility of the research. On the other hand, open data-sets exist, but are usually limited to certain regions. It is well known that results are sensitive to the input data-sets, and the usage of heterogeneous data-sets across the world may have led to a trend in developing ad-hoc measures which model the desired behavior for the input data-set, while eventually limiting results to that region of analysis.

Several works in the literature have to cope with the modifiable areal unit problem (MAUP), naturally produced by the aggregation of data. Outputs of measures depend on the chosen geography of division. The sensitivity of the results to the chosen aggregation strategy is undesirable in any application and is a major concern when trying to compare metrics. The usage of geo-localized input data, together with the estimation of a continuous surface by means of a kernel density estimation method allows to cope with the MAUP problem [START_REF] Carlos | Density estimation and adaptive bandwidths: a primer for public health practitioners[END_REF].

IV. DATASET

The OpenStreetMap 1 (OSM) is a collaborative project to create a free editable map of the world; it is a prominent example of volunteered geographic information (VGI). It is a knowledge collective that provides user-generated street maps [START_REF] Haklay | OpenStreetMap: User-generated street maps[END_REF]. Volunteers across the world share geographic information on OSM in various ways and are sometimes considered as "intelligent sensors" [START_REF] Goodchild | Citizens as sensors: web 2.0 and the volunteering of geographic information[END_REF].

Since its creation, the project has been increasingly used across the world for a wide variety of purposes. Quality metrics have been proposed in [START_REF] Forghani | A quality study of the OpenStreetMap dataset for Tehran[END_REF], [START_REF] Barron | A comprehensive framework for intrinsic OpenStreetMap quality analysis[END_REF], [START_REF] Mooney | Towards quality metrics for OpenStreetMap[END_REF], [START_REF] Fan | Quality assessment for building footprints data on OpenStreetMap[END_REF], followed by different quality assessments, in particular for different countries. For instance, it has been concluded that the quality is "fairly accurate" for England [START_REF] Haklay | How good is volunteered geographical information? A comparative study of OpenStreetMap and Ordnance Survey datasets[END_REF], and it is even shown that OSM data is superior to the official data-set for Great Britain Meridian 2. Thus, the previous work has been extended for France [START_REF] Girres | Quality assessment of the French Open-StreetMap dataset[END_REF]. Studies focusing on the street network of Germany have been also been conducted in [START_REF] Neis | The street network evolution of crowdsourced maps: OpenStreetMap in Germany 2007-2011[END_REF], where it is concluded that the data-sets can be considered complete in relative comparison to a commercial data-set.

In addition, the OSM data-set for Hamburg already covers about 99.8% of the street network [START_REF] Over | Generating webbased 3D city models from OpenStreetMap: The current situation in Germany[END_REF] according to the surveying office of Hamburg. The latter study also remarks that "Besides the street network, the real advantage of the dataset is the availability of manifold points of interest". These points of interest allow for a deeper understanding of city dynamics, enriched with the provided location and embedded 1 http://www.openstreetmap.org/ information. In China, the volume of points of interest has been increasing substantially, e.g. nine-fold in the period 2007-2013 [START_REF] Liu | Automated identification and characterization of parcels with OpenStreetMap and points of interest[END_REF].

In this work, we use OSM data for analyzing cities' land use mix, which leads to numerous advantages. Firstly, worldwide coverage is a huge asset. It is the first step towards cities metrics for further comparisons using homogeneous data-sets. Secondly, it is being continuously updated by means of crowdsourcing, allowing to adapt its information to the rapid urban changes present nowadays. In addition, there exists a very active community, which iteratively improves the data precision and completeness. Thirdly, OSM is Open data. The fact that data is freely available for everyone allows for coherence and valid comparison between different contributions in the field of urban planning. As well, results are reproducible for anyone, an important aspect for the community, towards improval of these tools. Finally, geo-localized data provide a great advantage in terms of spatial location specificity, granting the possibility of a finer analysis in comparison to gridded data for instance.

Currently, the biggest limitation related to OSM is still missing data. Even though it can occur quite frequently, it is reduced drastically in big cities where lots of contributors exist, coming from a growing community of crowd-sourcing.

In contrast, Land Use Land Cover (LULC) data has several limitations compared to OSM. For LULC, availability is generally restricted to country level boundaries. Then, even though resolution has been improving considerably, intensities of activity and residential uses are not captured. Further, the classification by means of using a single allocation of land use, and the aggregated land use categories provided are insufficient to infer cities' land use mix [START_REF] Christian | How important is the land use mix measure in understanding walking behaviour? Results from the RESIDE study[END_REF]. Last but not least, accounting for the increasing speed of change in urban dynamics, LULC data-sets can get outdated in relatively short time.

V. METHOD

In this section, we describe the computing pipeline that allows us to compute the continuous spatial representation of the land use mix from raw OSM data.

A. Data extraction

We explain here how to retrieve the different land uses which will be analyzed later to define the land use mix degree. The OSM data for a given geographical area, defined by its geographical bounding box, is retrieved in the shapefile format, obtained from Mapzen Metro Extracts2 . These files are the result of the osm2pgsql process, where the OSM data is converted to postGIS-enabled PostgreSQL databases. From these files, points and polygons are used as input data for further processing.

We first perform a classification of both the points of interest (POIs) and the polygons based on residential or activity uses. It is done following the OSM Wiki 3 associated to a certain land use type according to their input information as denoted in Table I.

In OSM data-sets a great quantity of polygons are tagged as buildings, without any additional information. In such cases, the polygons' land uses are inferred as detailed in Table II. First, polygons containing a defined key for land use are sorted in residential, activity, and other uses (e.g. forest, water). Let P R , P A and P O denote the three sets of polygons associated with these land uses. All polygons tagged as buildings are processed to estimate their containing land use. For a polygon P , the land use is estimated as follows:

LU (P ) = LU (p) if ∃p with min p∈P R ∪P A ∪P O A(p) , P ⊆ p
where A(•) denotes a polygon's area, and LU (•) its land use. Note that polygons with different land uses might overlap. We consider as relevant the information contained in the smallest encompassing polygon whose land use is known. Sometimes, on the contrary, it might happen that a polygon which needs to be inferred is not contained in any polygon with a defined land use value. A residential purpose is then assumed for these polygons. This hypothesis is made in order to somewhat counter the fact that in OSM data-sets, still relatively few residential tags exist. In the future, this hypothesis may be dropped. Experimentally, it has been observed that polygons which are uniquely tagged as buildings (i.e. no other complementary information exists) correspond in great part to residential buildings. Firstly, it is important to note that otherwise, the building polygon would have been tagged with complementary information (e.g. any activity classification). Secondly, there exists no land use tag under the area of the building, which is related to a non-residential land use (i.e. natural or industrial land uses).

Later, all residential buildings with a computed squared footage smaller than 12 squared meters, are filtered out. This is done due to the high probability that such buildings are a false positive in the classification, rather than a true residential building.

Finally, the polygons are converted to POIs by computing their centroid. The objective is to obtain a full distribution of geo-referenced points which contain information on residential or activity land uses, throughout a city. The results of this procedure are depicted in Fig. 1 for London, England, where transparency is used to clearly distinguish highly concentrated zones.

B. Kernel Density Estimation

In this section we describe the process of performing a Kernel density estimation (KDE), a statistical process for density estimation, for both activity and residential uses. As stated in [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF], the KDE is a well-recognized technique for visualizing and analyzing complex and technical data in a clear and understandable way to non-mathematicians.

The KDE infers a density function f , given an i.i.d. sample x 1 , x 2 , ..., x n from the corresponding probability distribution. The kernel density estimator is defined as

f (x) = 1 nh n i=1 K x -x i h , (1)
with kernel K(•) and h a smoothing parameter denominated bandwidth. In our case, x i are two-dimensional vectors containing the latitude and longitude coordinates. Then, given a set of data points, the KDE interpolates a continuous surface using a given kernel (defined e.g. by a Gaussian/normal function). This procedure has two major advantages, namely spatial smoothing and spatial interpolation. Both aspects are relevant in our context, where missing data can occur. Thus, centers of activities and residential purposes can be smoothed, while interpolation is done locally to cope with missing data. As described before, the KDE is controlled by a weighting function K and the bandwidth h. For the weighting function, we assume that both residential and activity uses have a normally distributed spatial influence. For instance, the impact of an activity is maximal in its original location, and reduces with increasing distance, according to the normal distribution, which is the default choice in most works using KDE.

The bandwidth is associated with the distance of influence. It is important to mention that several techniques for an automatic computation of the bandwidth exist in the literature, but in several applications the desire is to reduce complexity during the estimation at the expense of precision. In our case, we are mostly interested in capturing information (i.e. residential and activities points) within local neighborhoods. It is thus natural to relate the bandwidth to the extent of walkable distances.

In [START_REF] Atash | Redesigning suburbia for walking and transit: Emerging concepts[END_REF], it is suggested that 400 meters correspond to the distance an average American will walk rather than drive. According to [START_REF] Aultman-Hall | Using GIS for evaluation of neighborhood pedestrian accessibility[END_REF], this distance is also considered as the greatest distance someone is likely to walk to a transit station. Similarly, in [START_REF]Liveable neighbourhoods: a Western Australian Government sustainable cities initiative[END_REF] the implementation of an average destination within a distance of between 400 and 450 meters is promoted. Nonetheless, some studies promote walkable distances slightly higher than the one adopted, as for example in [START_REF] Mccormack | The relationship between destination proximity, destination mix and physical activity behaviors[END_REF] that studied the association between physical activity and the mixture of destinations located within 400 and 1500 meters of residents' homes.

In our context, and for the purpose of favoring the neighborhoods with high uses mix, we decided to adopt the suggested value of 400 meters. This value defines the spatial bandwidth during the density estimation procedure.

As mentioned before, the KDE is a very efficient tool to smooth spatial data and interpolate locally-missing data. However, in the case of strict geographical borders like frontiers or coastal regions, it can lead to an undesirable over-smoothing and the estimation of non-existent information. Such a case will be studied in detail in further work.

As done in [START_REF] Song | Comparing measures of urban land use mix[END_REF], we chose to adopt two land use types: residential and activities (i.e. related to non-residential uses, such as shopping, leisure, etc.). Consequently, all geo-localized POIs from the procedure described above, are used to compute one KDE for each category: residential and activities KDE. The probability density function is then evaluated in a grid of points covering the region of analysis with any resolution that a user could want. In practice, we construct a grid with a step of 100 meters. This value was set in order to aid the visualization of the land uses densities on the different neighborhoods within a city. Finally, the computed grid points are normalized for each computed KDE.

The resulting KDE's are depicted in Fig. 2, applied to the city of London, England.

C. Computation of land use mix

In this section we focus on the computation of a land use mixture measure, given the input densities estimated for residential and activity uses respectively. The goal of this measure is to determine to what extent the spatial configuration of land uses is well distributed in a city, as well as to assess the co-occurrence of both residential and activity uses in neighborhoods. From a global point of view, it is desirable to achieve a high degree of land use mix in most neighborhoods.

As expressed in [START_REF] Song | Comparing measures of urban land use mix[END_REF], measures for land use mix must address two underlying concepts: distance and quantity. Distance describes the proximity between different land uses. We compute this closeness, or influence of a certain land use on another one, according to the walkability distance, as denoted in Sect. V-B. The quantity is modeled in the density estimation procedure, given the input of geo-localized activity and residential spots.

In accordance with the notion of distance and quantity, the land use mix measurement must depict the regions where a high use mix exists. By contrast, it will also point out the poorly distributed usages with low diversity of land uses.

As expressed before, this framework intends to be an aid for urban planners, and thus the visualization of geo-spatial land use mixture is of paramount importance to provide support during the decision-making procedure. In this context, this framework is also designed to expose the importance of the different sub-regions within a city. The grid meshes computed from KDE results allow for determining the different subregion's importance in terms of land use intensity. This provides a great aid for urban planners on the need for improving and planning mixed use development at particular locations.

In our first results, the land use intensity is illustrated using a bubble plot as shown in Fig. 4. The intuition behind this is to help visualize the land use mix in the core regions of a city, where land use intensity, either residential, activity, or both, is high. Other methods fail to highlight this important aspect, which is crucial in distinguishing important sub-regions relative to others.

In this work we chose to compute the Entropy Index (EI), as done in [START_REF] Song | Measuring the effects of mixed land uses on housing values[END_REF], to evaluate the degree of land use mix for every point along the two spatial dimensions of the KDE's. Still, our method allows to use any kind of index computation, given the local context of analysis built from OSM data.

Let P j be the percentage of each land use type, obtained from the respective KDE's. The amount of land uses k employed equals the number of KDE's computed (here, two, although with increasing completeness of OSM data, this may increase):

EI = - k j=1 P j ln(P j ) ln(k) (2)
The Entropy Index output will be between 0 and 1, where the higher the value, the higher the mixture of land uses. Fig. 3 presents the computed land use mix for London, England, while Fig. 4 shows a bubble plot which better conveys local land use intensities. A further analysis of these results is presented in Sect. VI.

D. Implementation

The framework's infrastructure is built over a Python stack and is composed of the following modules:

• extract_uses.py given the points and polygons shapefiles from the Mapzen Metro Extract, the extraction of their uses is done according to Section V-A. The geometry of the polygons is processed with the Shapely library, and the final output is a Pandas [START_REF] Mckinney | Data structures for statistical computing in python[END_REF] DataFrame with the geo-referenced categorized POIs pois_df. • kde.py given pois_df, it determines the KDE for both activity kde_act and residential kde_res types using the Statsmodels [START_REF] Seabold | Statsmodels: Econometric and statistical modeling with python[END_REF] library. • loaders.py manages a local Hadoop Distributed File Store (HDFS) through the [START_REF] Alted | Pytables: processing and analyzing extremely large amounts of data in python[END_REF] library. Such an HDFS stores for each studied city the results of the expensive computations of extract_uses.py and kde.py.

When for a given city such information is not stored in the local HDFS, it proceeds as follows: 1) A query for the corresponding city's shapefiles is issued to Mapzen 2) When the shapefiles are received, it calls extract_uses.py and kde.py in order to compute pois_df, kde_act and kde_res 3) The results of the computation are stored in the local HDFS under the corresponding city key.

• measures.py determines the implemented indicators out of the pois_df, kde_act and kde_res (i.e. the entropy of Equation ( 2)) using Numpy arrays [START_REF] Van Der Walt | The numpy array: a structure for efficient numerical computation[END_REF]. • plots.py generates the desired plots from pois_df, kde_act and kde_res (i.e. longitude-latitude scatterplots, KDE plots...) using the Matplotlib library [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF]. • city_analysis.py is the Python class to interact with the Jupyter Server. To start an analysis for a given city, it will ask for the city's information to loaders.py. Once such information is available, the class manages the calls to the methods of measures.py and plots.py in order to get the outputs interactively. The implementation is represented in Figure 5. The source code of the framework is publicly available at 4 . A webinterface is also planned to be released to allow non computersavvy, urban planners among others, to evaluate the framework online.

VI. APPLICATION A. Grenoble, France

The city of Grenoble, France, is an ancient medium-size European city. It is located at the foot of the French Alps, and possesses a distinctive Y shape because of the surrounding mountain ranges of the Chartreuse, Belledonne and Vercors. The agglomeration of this city, containing a population size of 665, 000, was processed and its LUM was computed. A The different outputs are depicted in Fig. 7. In the activity's KDE, the high concentration of activities lying around the historical center can be observed. This pattern can be frequently seen in ancient European cities, where the historical centers contain the biggest concentration of activities due to a long historical process.

Regarding residential uses, it is noticeable that the highest concentration can also be found around the historical center, in agreement with the mentioned historical process. For this type of cities, urban development is the result of a long process dating back to long before the invention of the car. Consequently, one can find a good land use mix around the cities' historical center, as depicted visually in the different land use mix plots.

The results also exhibit clearly a particular observable structure: as the outskirts of the city are mainly composed of half-mountainous landscape, several residential collections can be found in areas remote enough from the existing activity centers, providing a good quality of life and access to nature, but still easily accessible by car or public transportation. The structures at 5.58E/45.13N and 5.84E/45.14N match respectively with the East side of the plateau of the Vercors and the valley of Vaulnaveys. Both are notoriously beautiful residential areas with few activities apart from arts and crafts shops and farms.

B. London, England

Another evaluation was performed on the city of London. As a result, a total number of 80, 702 and 238, 626 points were extracted for activity and residential uses respectively. The results of this procedure are depicted above in Fig. 1. The KDEs were then computed, as shown in Fig. 2. The LUM is shown in Fig. 3 and Fig. 4.

London city went across big changes across history, a process which molded the structure in which the current 8.674 million people live. In Fig. 2, the City of London and Soho are depicted as the areas containing the highest density of activities. Regarding residential land use, East Dulwich, the Higham Hill, and the area lying between the Wandsworth Common and the Clapham Common were found to be the most intense.

The Croydon Vision 2020 is an urban planning program seeking to develop Croydon as a hub of living, retailing, culture and business in the context of a local development framework. An important concentration of activities exists in this area, while an important residential area lies nearby, to its west. Our framework can detect that a relatively good mix of uses already exists for this region, as depicted in Fig. 3 at the longitude and latitude values of -0.100594 and 51.376495 respectively. Further, mixed-used development is acknowledged in this urban program. Thus, its mixed use development is expected to further improve in the near future with a higher integration of different land uses.

VII. CONCLUSIONS AND FUTURE WORK

In this paper, we presented a framework for capturing mixed uses development in cities. It uses crowd-sourcing data from OpenStreetMap to extract the geo-localized land uses. Due to the universality of this data source, we are able to process any geographical area in the world, as long as sufficient data are available in OSM. A Kernel Density Estimation is performed for each of the land uses, outputting the spatial distribution of its land uses. Based on this representation, a measure of land use mix is then calculated using the Entropy Index. The GIS output that results, shows enriched information for urban planners, supporting and aiding the decision-making procedure.

The framework was applied on the cities of London and Grenoble. Both activities and residential densities estimation were validated by means of inspecting local neighborhoods with a high concentration of the different land uses. Afterwards, their spatial LUM distribution was computed, allowing for an easy visualization concentrating on the mixed uses values, and showing the LUM relative importance in terms of land use intensity.

Future work includes integrating the LUM output for measuring the urban sprawl phenomenon and performing numerical interpretations of desirable mixed use values. We will also study the potential integration to transportation models, where land use mix correlation with the activities and residential uses can help improve demand estimation.

In addition, further investigation can be done by means of analyzing the different types of activities, which in this work are classified only as activities. Rich information provided from OSM data allows to divide the extracted activities into differing classifications such as shop, leisure, amenities, commercial and industrial among others. Finally, the estimation of LUM can be refined by taking into account, besides their location, the accessibility between different land uses, which is partly conditioned by the transportation infrastructure.

Fig. 1 .

 1 Fig. 1. (a) OpenStreetMap of London, England. (b) and (c) depict, respectively, the extracted Activity and Residential uses.

Fig. 2 .

 2 Fig. 2. Estimated KDE's for London, England, for (a), (c) activity and (b), (d) residential uses.

Fig. 3 .

 3 Fig. 3. Land use mixture for London.

Fig. 4 .

 4 Fig. 4. Bubble plot of LUM for London. Each point's size denotes the local land use intensity.

Fig. 5 .

 5 Fig. 5. Implementation of the framework

Fig. 7 .

 7 Fig. 7. The framework's results for Grenoble, France. (a), (b) stand for Activity and Residential KDE's, while (c) and (d) stand for LUM.

  . Points and polygons are , restaurant, biergarten, cafe, fast food, food court, ice cream, pub, restaurant, college, kindergarten, library, public bookcase, school, music school, driving school, language school, university, fuel, bicycle rental, bus station, car rental, taxi, car wash, ferry terminal, atm, bank, bureau de change, baby hatch, clinic, dentist, doctors, hospital, nursing home, pharmacy, social facility, veterinary, arts centre, brothel, casino, cinema, community centre, foun-

	Key	Value
		Activities classification
	amenity	bar, pubtain, gambling, nightclub, planetarium, social centre, strip-
		club, studio, swingerclub, theatre, animal boarding, an-
		imal shelter, courthouse, coworking space, crematorium,
		dive centre, dojo, embassy, fire station, gym, internet cafe,
		marketplace, police, post office, townhall
	shop	*
	building commercial, office, industrial, retail, warehouse, cathedral,
		chapel, church, mosque, temple, synagogue, shrine, civic,
		hospital, school, stadium, train station, transportation, univer-
		sity, public, kiosk, garage, garages, hangar, stable, cowshed,
		digester
	leisure	adult gaming centre,	amusement arcade,	beach resort,
		dance, hackerspace, ice rink, pitch, sports centre, stadium,
		summer camp, swimming area, water park
	landuse	commercial, industrial, retail, port, quarry, salt pond, con-
		struction, military, garages
		Residential classification
	building hotel, farm, apartment, apartments, dormitory, house, res-
		idential, retirement home, terrace, houseboat, bungalow,
		static caravan, detached
		Other land use classification
	landuse	cemetery, landfill, railway, water, reservoir, basin, allot-
		ments, conservation, farmland, farmyard, forest, grass, green-
		field, greenhouse horticulture, meadow, orchard, pasture,
		peat cutting, plant nursery, recreation ground, village green,
		vineyard
		TABLE I
		LAND USES CLASSIFICATION.

https://mapzen.com/data/metro-extracts/

https://wiki.openstreetmap.org/

https://github.com/martibosch/landusemix


Thomas Capelle

Peter Sturm

Arthur Vidard

Brian Morton

Optimisation-Based Calibration and Model Selection for the Tranus Land Use Module

Instantiating land use and transport integrated models (LUTI modelling) is a complicated task, requiring substantial data collection, parameter estimation and expert analysis. In this work, we present a partial effort towards the automation of the calibration of Tranus, one of the most popular LUTI models. First, we give a detailed mathematical description of the activity module and the usual calibration approach. Secondly, we reformulate the estimation of the endogenous parameters called shadow prices as an optimisation problem. We also propose an optimisation algorithm for the calibration of the substitution submodel, setting a base for future fully integrated calibration. We analyse the case of transportable and nontransportable economic sectors and propose a detailed mathematical scheme for each case. We also discuss how to validate calibration results and propose to use synthetic data generated from real world problems in order to assess convergence properties and accuracy of calibration methods. Results of this methodology are presented for realistic scenarios. Finally, we propose a model selection scheme to reduce the number of shadow prices that need to be calibrated.

Introduction

LUTI (land-use and transportation integrated) models aim at representing the complex interactions between land use and transportation offer and demand within a territory. They are mainly used to evaluate different alternative planning scenarios, by simulating their tendential impacts on patterns of land use and travel behaviour. Since the early 60's LUTI modelling has attracted researchers that aimed to model the complex economical relations in urban areas; a good overview of the evolution and history of LUTI modelling can be found in [START_REF] Wegener | Overview of land-use transport models[END_REF]. Setting up a LUTI model requires the estimation of several types of parameters to reproduce as closely as possible, observations gathered on the studied area (socio-economic data, transport surveys, etc.). The vast majority of available calibration approaches are semi-automatic, estimating one subset of parameters at a time, without a global integrated estimation. Automatic calibration of LUTI model is not a common practice; an exception has been proposed for the Meplan model [START_REF] Abraham | Parameter Estimation in Urban Models: Theory and Application to a Land Use Transport Interaction Model of the Sacramento, California Region[END_REF].

We consider Tranus [START_REF] De La Barra | Mathematical description of TRANUS[END_REF], an open source LUTI model that is widely used. Tranus is a classical LUTI model, with two separated modules: the activity module and the transport module. The activity module, is an equilibrium based model based in micro-economic principles that balance the offer and demand of the different economical sectors that interact at each level. Economical sectors are considered in the broad sense, amongst them we have: land, goods, salaries, housing, transportation demand, etc. Also, the price paid for each economical sector has to be balanced with offer and demand, thus there are two equilibria that have to be achieved, offer versus demand and (production) cost versus prices. The transportation module, computes de costs of transportation and assigns the demand to the network. Both modules interact back and forth until a general equilibrium is achieved.

The calibration process is usually done by an expert modeller who iteratively tunes a group of parameters to reproduce as close as possible the observations gathered in the area of study. This process is usually done manually, with little to no automation, adjusting the different economical parameters (for example, the demand curves for different goods in a specific geographical zone). At the same time, Tranus computes internally a set of adjustment coefficients (called shadow prices in Tranus) that correct the utilities and account for the un-modelled effects. These endogenous variables help the model achieve a better response and fit more precisely to the observed data.

In this paper we will address several shortcomings of the classical approach of calibration used in Tranus. In our previous work [START_REF] Capelle | Formulating LUTI calibration as an optimisation problem: Example of tranus shadow price estimation[END_REF] we proposed a first approach for the reformulation of the heuristic calibration algorithm used in the land use and activity module as an optimisation problem. In this paper we extend this approach by having a closer look at the inner loop that computes the shadow prices and propose an efficient methodology for their estimation by decoupling the calibration in smaller problems. To be able to do this, we have to carefully investigate the system of equations that are computed in the activity module. We also introduce auxiliary variables, which enables the a closed form computation instead of an iterative one. This in turn makes it possible to use sophisticated numerical optimisation methods and opens the door to the simultaneous estimation of different parameter types of the model. The ultimate goal of this approach is to simultaneously calibrate the various parameters of Tranus' inner and outer loops. In this direction, this scheme gives valuable results for usually hard-to-calibrate parameters (so-called penalising factors).

We also present a detailed methodology for the construction of synthetic scenarios based on real calibrated study areas. These synthetic scenarios have a perfect fit without the need of the so-called shadow prices (usually we set their value to zero), enabling us to validate our optimisation algorithms knowing the ground truth values of the shadow prices. A simple example is presented to understand the problematic of synthetic scenario generation and the corresponding equilibrium prices problem.

Finally, we question the rationale of usual calibration approaches for Tranus (and other LUTI models), which consist in estimating parameters for which the model reproduces observations exactly. In Tranus, this is achieved by enriching the underlying macro-economic model with the already mentioned auxiliary variables, the shadow prices. While this allows to correct for unavoidable un-modeled effects, it also bears the risk of overparameterisation/overfitting. We propose a model selection scheme, aiming at a compromise between model complexity (here, number of shadow prices) and goodness of fit to observations, reducing the risk of overfitting and increasing the likelihood of achieving good predictions with a model.

Description of Tranus land use module

Tranus is a land use and transportation integrated model (LUTI), providing a framework for modelling land use and transportations in an integrated manner. It can be used at urban, regional or even national scale. The area of study is divided in spatial zones and economical sectors; the basic concepts of the original input-output model (see [START_REF] Leontief | Multi-Regional Input-Output Analysis[END_REF]) have been generalised and given a spatial dimension. The concept of sectors is more general than in the traditional definition. It may include the classical sectors in which the economy is divided (agriculture, manufacturing, mining, etc.), factors of production (capital, land and labour), population groups, employment, floorspace, land, energy, or any other that is relevant to the spatial system being represented. Tranus combines two main modules: the land use and activity and the transportation modules.

The land use module

The land use and activity module simulates a spatial economic system by modelling the locations of activities and the interactions between economic sectors for a specific time period. The transportation module, on the other hand, dispatches the travel demand induced by the activity model and assigns it to the transport supply. Both modules are linked together, serving both as input and output for each other. In this way the movements of people or freight are explained as the results of the economic and spatial interaction between activities, the transport system and the real estate market. In turn, the accessibility that results from the transport system influences the location and interaction between activities, also affecting land rent. The two modules use discrete choice logit models [START_REF] Mcfadden | Conditional logit analysis of qualitative choice behaviour[END_REF][START_REF] Mcfadden | Mixed MNL models for discrete response[END_REF], linked together in a consistent way. This includes activitylocation, land-choice, and multi-modal path choice and trip assignment.

To attain a convergence status, Tranus runs both modules iteratively until an equilibrium is found. The land use and transportation modules need to reach their own respective equilibrium status. First, the land use module needs to achieve equilibrium between offer and demand, and equilibrium between the price paid and the cost of producing each economic sector. This is done at current transportation costs and disutilities. Secondly, the transportation module takes as input the transport demand and equilibrates the transportation network to satisfy the given demand. Both modules are ran iteratively until a general equilibrium status is found. This is achieved when neither land use nor transportation, evolve anymore, as illustrated in Figure 1.

In the activity model, there are two general types of economic sectors; transportable and non-transportable sectors. The main difference between these two types of sectors, is that transportable sectors can be consumed in a different place from where they were produced. As an example, the demand for coal from a metal industry can be satisfied by a mining industry located in another region. On the other hand, a typical non-transportable sector is floorspace: land is consumed where it is "produced".

Transportable sectors generate flux, that induces transport demand, which ultimately influence transportation costs. Non-transportable sectors, on the other hand, neither require transportation nor generate fluxes. Usually, three types of economic sectors are classified: land or floorspace, households and businesses. Land is usually composed of two or three types of residential floorspace (e.g. detached houses, apartments, mobile homes), and commercial floorspace of offices and stores. Households are usually classified by socio-economic level, based on income or the household composition. Business sectors comprise industries (whose output is mainly destined for exportation), services (schools, universities, recreational) and commerce. The standard approach for the consumption chain is as follows: Industry has a demand for labour (households) and service businesses. Households also consume services, and services also require labour, thus "consume households". Finally, all businesses and households consume land. For instance, households will locate in residential zones, and the feedback of household and business "consumption" will induce home-to-work trips (see [START_REF] Lowry | A model of metropolis[END_REF]). This process results in economic exchanges, sometimes inducing flux (transportable sectors) and sometimes in-place consumption (land). The offer and demand is equilibrated and a set of equilibrium prices for each economic sector is attained.

The land use module's objective is to find an equilibrium between the production and demand of all economic sectors and zones of the modelled region. To attain the equilibrium, various parameters and functions are used to represent the behaviour of the different economic agents. Among these parameters are demand elasticities, attractiveness of geographical zones, technical coefficients, etc. In the following, we introduce the parts of the terminology, parameters and equations used in Tranus that are relevant to this paper. See [START_REF] De La Barra | Mathematical description of TRANUS[END_REF] for a complete description. • Productions: X n i expresses how many "items" of an economic sector n are present/produced in a zone i.

• Demands: D mn i expresses how many items of a sector n are demanded by the part of sector m located in zone i.

• Prices: p n i defines the price of (one item of) sector n located in zone i.

It is important to realise that "price" in the case of land, is the actual rent, whereas the price of a household is derived from the salary.

Productions, demands and prices form part of a dynamic system of equations. These equations depend on one another, and are linked by a list of equations that need to be computed one after another. This is detailed in [START_REF] De La Barra | Mathematical description of TRANUS[END_REF]. A graphical representation of this feedback is represented in Figure 2. For instance, demand induces production and vice-versa. The iteration scheme is as follows: prices of a current iteration translate into intermediate variables (that will not be detailed here) which enables the computation of demand and consumption costs (noted as c in Figure 2). This is done based on the current transportation costs and disutilities. Once demand and costs are known, the current production is computed and fed back to compute a new set of prices, for a next iteration.

The process is bottom-up, starting with land use prices and exogenous production and demand up to the production and prices of transportable sectors. All the above computations are repeated until convergence is attained in productions X and prices p at the same time (convergence in these two sets of variables implies convergence in all others).

We are only going to deal with a subset of model equations relevant to this paper. Demand is computed for all combinations of zone i, demanding (consuming) sector m and demanded sector n:

D mn i = (X * m i + X m i ) a mn i S mn i (2.1)
D n i = D * n i + m D mn i (2.2)
where X * m i is the given exogenous production (for exports), X m i the induced endogenous production obtained in the previous iteration (or initial values), and D * n i exogenous demand. D n i in (2.2) then gives the total demand for sector n in zone i. a mn i is a technical demand coefficient and S mn i is the substitution proportion of sector n when consumed by sector m on zone i (explained later in more detail).

In parallel to demand, one computes the utility of all pairs of production and consumption zones, j and i:

U

n ij = λ n (p n j + h n j ) + t n ij . (2.3)
Here, λ n is the marginal utility of income for sector n and t n ij represents transport disutility. Since utilities and disutilities are difficult to model mathematically (they include subjective factors such as the value of time spent in transportation), Tranus incorporates adjustment parameters h n j , so-called shadow prices, amongst the model parameters to be estimated.

From utility, we compute the probability that the production of sector n demanded in zone i, is located in zone j. Every combination of n, i and j is computed:

P r n ij = A n j exp -β n U n ij h A n h exp (-β n U n ih )
.

(2.4)

Here, h ranges over all zones, A n j represents attractiveness of zone j for sector n and β n is the dispersion parameter for the multinomial logit model expressed by the above equation.

From these probabilities, new productions are then computed for every combination of sector n, production zone j and consumption zone i:

X n ij = D n i P r n ij . (2.5)
Total production of sector n in zone j, is then:

X n j = i X n ij (2.6) = i D n i P r n ij .
(2.7)

Given the computed demand and production, consumption costs are computed as

cn i = j X n ij p n j + tm n ij D n i (2.8)
where tm n ij is the monetary cost of transporting one item of sector n from zone j to zone i.

These finally determine the new prices:

p m i = V A m i + n a mn i S mn i cn i (2.9)
where V A m i is value added by the production of an item of sector m in zone i, to the sum of values of the input items.

Calibration

The calibration process consists in adjusting the model parameters to be able to reproduce a base year's data in the study area. Obtaining a good calibration is a long process, that is usually performed by experts and can take months. A mix of tools are used to estimate the various parameters of the model. Econometrical, ad-hoc procedures and interactive trial-and-error can be counted among the tools used by experts to obtain a good fitting model.

For the calibration phase, parameters are separated in three sets: i. Parameters that are computed externally using the appropriate data and econometrical techniques.

ii. The adjustment parameters h n j of the utilities (2.3), known as shadow prices.

iii. The remaining parameters (for example the penalisation factors and logit dispersions parameters).

After computing the external parameters (set i), and giving initial values to set iii, the model iterates until convergence. The iteration process is constructed in such a way, that the shadow prices will be adjusted to force the productions to reproduce the observed productions X 0 in the study area. These variables will "try" to compensate for the other parameters to have a perfect fit; they act as correction terms to compensate for parts of the utility that are not represented by the model. One wants to make the values of the shadow prices as small as possible. This process of parameter calibration is done repeatedly until the expert modeller is satisfied with the parameters and the values of the shadow prices.

The computation of the shadow prices is automatically done as follows at the end of each iteration (cf. figure 2 and the above equations):

h n,t+1 i = (h n,t i + p n,t i ) X n,t i X n 0,i -p n,t+1 i .
(2.10)

The shadow prices for the next iteration t + 1 increase proportionally to the excess of computed, as compared to observed, productions.

Proposed Calibration Approaches

Our main motivations are to replace the sequential calibration process outlined above by a process that rigorously estimates as many parameters as possible, taking into account all available constraints and assumptions in a systematic manner, to automise as much as possible the calibration process, and to make it more reproducible. We believe that a natural way of achieving these goals is to explicitly formulate the calibration process in terms of a cost function (or possibly, as a multi-criteria decision problem) that is to be minimised or maximised, with respect to a set of constraints, when given. This is for example not directly the case in the existing approach, where the estimation of shadow prices and other parameters is done without a definition of a clearly defined cost function. Formulating calibration via explicit cost functions enables to use the rich variety of optimisation algorithms existing in the literature and in numerical libraries.

A first step in this direction concerns the estimation of shadow prices, a second step deals with the automatic estimation of both shadow prices and other parameters; these two steps are described in the following.

Reformulating as an optimisation problem

It is important to notice that a calibration of the land use module involves the estimation of all the parameters of the model to make productions as close as possible to the base year data. To reformulate the calibration as an optimisation problem, we must compute the shadow prices that make the productions as similar as possible to the observed productions. This can be achieved by the following optimisation problem:

min h X(h) -X 0 2 . (3.1)
Here, h is a vector containing all shadow prices, X 0 the vector of observed productions, and X(h) the vector of productions computed by the model, after convergence of the iterative process shown in figure 2. The dependency of these on the shadow prices is visible from equations (2.3) to (2.7). Each evaluation of the productions X(h) involves the convergence of the dynamic system exposed in Figure 2. Each evaluation of the cost function involves the convergence of the dynamic system in productions as well as prices.

The dual convergence can be avoided by taking the prices as input variables. Moreover, one can compute directly productions that are in equilibrium for a given set of shadow prices and prices. To do so, we observe that the computation of demand and production involves a set of linear equations (2.1), (2.2), (2.5), and (2.7). If we re-organise these equations, knowing that only productions are needed in our cost function, we may only need to compute these. To do so, we substitute D n i in equation (2.5) using equations (2.1) and (2.2), giving:

X n ij = D * n i + m (X * m i + X m i ) a mn i S mn i P r n ij . (3.2)
Upon substituting this into (2.7), we obtain the following linear system in X n j :

X n j = i D * n i + m X * m i a mn i S mn i P r n ij + i m a mn i S mn i P r n ij X m i . (3.3)
By construction, the solution of this linear system represents an equilibrium of production and demand.

Nevertheless, computing the gradient of cost function (3.1) is still difficult. Each evaluation of the productions involves solving a linear system of the type (3.3). Estimating the gradient numerically via finite differences, is possible but rather costly. Moreover, even if productions computed this way are in equilibrium, the prices p still need to iterate until convergence is obtained.

Another observation can be made: we already have a base year production data set. For a successful calibration, we want to have the computed productions equal to the base year productions. Hence, we can simply impose this condition by replacing productions in the right hand side of (3.3), with the observed base year productions. This approach enables us to compute the productions directly, without the need to solve a linear system. To address the second problem (equilibrium of prices), we add the prices explicitly to the set of parameters to be optimised. We use the current values of prices, and compare them against the prices computed by the model in the next iteration, cf. (2.9). The difference between the current prices and the ones computed by the model trough equations (2.3) to (2.8), is added to (3.1), in order to form a new cost function:

min h,p X(h, p, X 0 ) -X 0 2 + p(h, p, X 0 ) -p 2 . (3.4)
Here, p is the vector of prices computed by the model using (2.9) and the notation X(h, p, X 0 ) shows that modelled productions are computed as explained above by substituting observed productions X 0 into the right-hand side of (3.3).

The above cost function has a closed-form that permits us to compute the derivatives directly. No more iterations or waiting for convergence is required in this approach. The cost function (3.4) is of (non-linear) least squares type, meaning that any least squares optimisation approach can be used; in our work we apply the Levenberg-Marquardt method [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF].

Let us also note that other choices than the L 2 norm would of course be possible to define the cost function of (3.4). We may also weight the two terms differently, in order to favour equilibrium in production over that in prices or vice-versa in cases where a global equilibrium cannot be reached.

Land use sectors

Land is a very particular economical sector, it must be consumed where it is produced. Moreover, land does not consume other economical sectors and the amount of available land is fixed. For the calibration purpose, the prices for the land use sectors are known, this means the p n i variables for the calibration year are considered as input and do not need to be computed. This translates into a simplified set of equations for the computation of production of land. We have to detail two extra equations to understand how this enters our optimisation scheme. First, as land is non-transportable, the location probability (2.4) vanishes, so equation (3.3) can be re-written as:

X n i = D * n i + m (X * m 0i + X m 0i ) a mn i S mn i . (3.5)
We now detail the two variables a mn i (technical demand coefficients) and S mn i (substitution probabilities):

a mn i = min mn + gap mn exp(-δ mn (p n i + h n i )) . (3.6)
Where min mn , gap mn and δ mn are parameters of the demand function estimated externally using data from the rents. The substitution probability S mn i is only a function of shadow prices associated to the same zone i, see the next section for more details. As land prices are known, we can clearly see that the production of land X n i is only a function of the shadow prices of the same geographical zone i (there is no sum or dependence involving other zones j). Of course we have interactions between other economical sectors in the same zone, but in practice, the number of economical sectors is much smaller than the number of zones, which leads to optimisation problems (one per zone) that are very small, with the number of variables equal to that of land use sectors. We can re-write the optimisation problem (3.4) as one optimisation problem for each geographical zone i:

∀i min h i X i (h i , X 0 ) -X i0 2 (3.7)
Just as an example, for the North-Carolina-1 model (see later), there are only 3 land sectors: apartments, mobile-homes, detached houses.

Once the optimisation is done for each geographical zone and the shadow prices for land use are computed, we can proceed to computing the optimal shadow prices of the transportable sectors. We will further exploit this feature of the model to obtain an automated calibration of the substitution parameters. The details of the derivative estimation for the optimisation algorithm for the land use sectors can be found in the appendix.

New approach for computing the induced production of transportable sectors

Transportable sectors, are economical sectors that consume (and can be consumed) in a different location from where they are produced. Housing and commerce are examples of such sectors. For this type of sector, the technical coefficients a mn i are considered constant and in all practical usages of Tranus, the substitution probabilities S mn i = 1, i.e., there is no substitution considered in the transportable economic sectors (this is not a limitation of the model, rather it is common practice). Realising this implies that the total demand D n i is not a function of the prices or the shadow prices, and enables the computation of the total demand D n i (2.2) for each transportable sector n and geographical zone i only as a function of the base year data X 0 . In fact, this is done just after computing the land use sectors' shadow prices.

Thus, in the computation of the induced productions of a transportable sector n, X n j (cf. (2.7): X n j = i D n i P r n ij ),we only need to determine the values of the location probabilities P r n ij . If we go back to the definition of the location probabilities (2.4) and the underlying utilities (2.3), we realise that the utility makes no distinction between the price and shadow price part, so if we set:

φ n j = p n j + h n j (3.8)
the location probability can be computed as a function of φ. Instead of posing the induced production as a function of (h, p, X 0 ), we can look at the induced production X(φ, X 0 ) as a function of φ. Obtaining the optimal values of φ that minimise the difference between computed and observed productions, is the solution to the following problem:

min φ X(φ) -X 0 2 . (3.9)
Since the location probability P r n ij is a function only of the φ n variables for the same sector n, we get one optimisation problem for each economical sector n. Each of these optimisation problems is relatively simple and small in size, there are as many variables as geographical zones. The gradient of the cost function can be computed analytically using the well known derivatives of the logit probability P r. We use the Levenberg-Marquardt method to solve each problem. Once all the optimal values of P r have been computed, we can compute the prices by solving the linear system (2.9) for prices. Doing so allows us to recover the shadow prices from φ, subtracting the prices from the respective optimal φ values.

One consideration that one has to deal with, is that the location probabilities P r follow a logit formulation, so they can only be identified up to a constant per economical sector. This is a known property of logit models. As the prices are obtained from equation (2.9), this approach is considerably simpler and more stable than solving the double-objective optimisation approach proposed in (3.4), moreover, it exploits every little detail of the formulation of each function of the model. It also permits to calibrate incrementally, starting by the land use sectors and then obtaining the calibration in the transportable sectors. From the mathematical point of view it is also simpler, because the large optimisation problem in (3.4) is now decoupled into smaller optimisation problems, with fewer variables, allowing the modeller to finish the calibration of one set of variables before moving to the next stage.

To summarise, in the case of land use (non-transportable) sectors, there is one small optimisation problem to be solved for each geographical zone, whereas for the transportable sectors, we have one optimisation problem per economic sector.

We encountered some numerical issues relied to the fact the the location probability would vanish for large values of the utility function; this behaviour is explained in D.

Simultaneous estimation of shadow prices and land use substitution parameters

As stated in the introduction, one would like to have a simultaneous estimation of the whole set of parameters. In this section we present one step in this direction, we have constructed a two-phase algorithm that permits the estimation of the shadow prices and the substitution parameters at the same time. We have chosen the penalising factors in the substitution sub model because these are very hard to calibrate parameters, as relevant data are not readily available. The scheme we propose exploits the fact that the substitution sub model is used for land sectors, where we have already a simplified computation of the productions, as explained in Section 3.2. Tranus models include a discrete choice sub-model that represents the households' ability to choose among different types of residential buildings (i.e., floor space). That choice is captured in the "substitution model". The model is driven by the substitution probabilities:

S mn i = W n i exp(-ω mn a mn i • (p n i + h n i )) l∈K m W l i exp (-ω ml a ml i • (p l i + h l i )) . (4.1)
Here, K m represents the set of substitutes that sector m has access to, for example, for "rich" households m, this could be K m = {condos, detached houses}. Using Tranus terminology, W n i is an "attractor", a parameter that represents attributes of floor space sector n other than cost (utility); it is specified (and potentially calibrated) for each zone in which sector n is present. From equation (3.6) we can see that the demand coefficient a mn i is also a function of the prices and shadow prices. It is important to remember that prices are known for land sectors.

Simultaneous calibration of shadow prices and penalising factors

1. Phase 1: estimating parameters' initial values with multinomial logistic regression. The substitution model's parameters are estimated with multinomial logistic regression [START_REF] Mcfadden | Conditional logit analysis of qualitative choice behaviour[END_REF]. The data that are essential for estimation are household level observations on floor space consumption, housing expenditure, and the Tranus sector to which the household belongs. The dependent variable in a regression will be the choice of floor space sector, and the independent variable is the housing expenditure. The regressions are conducted separately for each household sector, and they yield estimates of -ω mn for each combination of floor space sector and household sector.

2. Phase 2: fine tuning the penalising factors. The penalising factors estimated in Phase 1 probably still need to be fine tuned to reduce the differences between the predicted production of floor space and the observed production of floor space. Fine tuning probably would also be necessary to achieve reasonable values of the floor space sectors' shadow prices.

If we consider all of Tranus' parameters fixed except the parameters ω, and include these parameters in the optimisation problem presented in (3.4), we obtain the following cost function:

f (h, ω) = X(X 0 , h, ω) -X 0 2 . (4.2)
We would like to find the values of ω that reduce the corresponding shadow prices. We propose to solve the following equation:

min ω∈Ω f (h = 0, ω) (4.3)
where Ω is a set of bounds over the penalising factors ω. We used a conjugategradient algorithm to solve this problem, and the starting points are the values obtained from the Multinomial Logistic regression of Phase 1. If we call ω * the solution of (4.3), then the final values for the shadow prices for the land use sectors are:

h * = arg min h f (h, w * ) .
The details and derivative estimation are presented in the appendix.

Generation of Synthetic Scenarios

The evaluation of a LUTI calibration is a difficult process, mainly due to the noise in the data and the fact that obtaining ground truth information is almost impossible.

Our optimisation scheme needs as input the base years' productions and parameters (X 0 , parameters). Then, the calibration is done against this information. We could think of a model that does not need the shadow prices to attain a perfect fit, hence, create a synthetic scenario where a "perfect" fit is achieved with shadow prices set to zero. To generate this "perfect fit" scenario, we have to solve a subproblem of the original calibration optimisation problem exposed in (3.4), where we do not consider the observed productions. We only need to obtain equilibrium in prices, and the values of the induced productions afterwards. To do so, we replace the consumption cost equation (2.8) in the prices (2.9), and by identifying the location probability as P r n ij = X n ij /D n i , we obtain the following system:

p m i = V A m i + n a mn i S mn i j P r n ij (h, p) • p n j + tm n ij p(h,p)
.

(5.1)

The dependence of the location-choice probability makes this system hard to solve even for small models. Our approach to solve this fixed-point equation is to solve the following optimisation problem:

min p p(h, p) -p 2 .
(5.2)

We have to make sure that the solution of (5.2) is a set of prices that are in equilibrium, that is for which p = p. After obtaining convergence in the prices, we compute the productions and then use them as observed base year productions in our synthetic scenario. This methodology produces a scenario where the optimal value of the shadow prices is zero (by construction) and that reproduces the base years' productions perfectly. We could also set the shadow prices value to any other value than zero here. The generation of such a synthetic scenario enables us to test our calibration methodology and optimisation algorithms against a known optimal value (shadow prices equal to zero). The details of the implementation are in C.

Equilibrium prices in synthetic scenario: 1 economical sector, 2 zones

In this section we present a simple example that shows that the equilibrium pricing problem can be complicated, and that uniqueness of the solution is not guaranteed. Let us consider only one economical sector m = 1 (we will just drop the exponent m in the following) and two geographical zones i, j ∈ {1, 2}. Let us consider no substitution, i.e. S 11 1 = S 11 2 = 1. The equilibrium condition (5.1) can be re-written by two equations:

p 1 = V A 1 + a 1 • P r 11 • p 1 + a 1 • P r 12 • (p 2 + tm 12 ) p 2 = V A 2 + a 2 • P r 21 • (p 1 + tm 21 ) + a 2 • P r 22 • p 2 (5.3)
It is important to notice that tm 11 = tm 22 = 0. This simple case is very sensitive to the values of the different parameters. We managed to find combinations of the different parameters (V A, a i and tm ij ) that give rise to multiple solutions, one solution or no solution at all, for the prices, see Figure 3. The curves where very sensitive to the demand coefficient a i . This example shows that modifying a certain parameter can shift the whole set of prices to a different equilibrium, and that the modeller has to be aware of this behaviour.

We were curious to know if this problem had multiple solutions, and even in this simple case it had proven to be complex. This gives us a starting point to further investigate the problem of potential existence of multiple fixed points for our calibration problems, even though in practice we have not observed problems of convergence to wrong local minima.

6 Results

Results of Shadow Price Calibration

We consider the example exposed in [START_REF] Capelle | Formulating LUTI calibration as an optimisation problem: Example of tranus shadow price estimation[END_REF], a simple model that allows to illustrate the methodology for generating synthetic scenarios with perfect fit (with "ground truth" shadow prices equal to zero). We applied our approach to the Example C model from the Tranus website1 , a small model with 3 zones and 5 sectors. First, we generated synthetic data from that model as described just above, with shadow prices h n i = 0 for each sector n and zone i. As expected, the cost function is zero at h = 0, and increases its value when we get away from the optimum. The cost function appears to be locally convex near the optimal value, cf. figure 4. If we consider for example sector 1 and zone 1, we can plot a "slice" of the cost function (3.4) near the optimal value h 1 1 = 0, p 1 1 = 2.676 as shown in figure 4. Here we can observe that as the shadow price gets larger the cost increases up to a plateau state (X 1 1 (h) → 0). In the case of the price p, if we move away from the optimal value p = 2.676, the cost increases quadratically.

We tested the robustness of the optimisation scheme with 1,000 random initial sets of shadow price values; the optimisation procedure outlined in this paper always converged to the ground truth solution. The initial values of shadow prices in these random trials were generated from a uniform distribution in [-10, 10], which is highly representative (prices are in the interval [0, 4] and nearly all shadow prices of a model are in practice smaller than the corresponding prices). ) near the optimal value (0, 2.676).

We also applied the same procedure to a Tranus model for North-Carolina (North-Carolina-1 model, consisting of102 zones and 12 economical sectors) modified by our synthetic data methodology. After setting the desired value for the shadow prices to h = 0, we tried 10,000 random initial sets of shadow prices values; and the algorithm proved to converge to the correct shadow prices for every single starting point. As for the calibration procedure implemented in Tranus' release, it failed to converge when starting values were too far away from the zero vector. We considered initial shadow prices uniformly distributed in the interval [

• -p max , • p max ],
where is a parameter in [0, 1] and p max is the maximum observed price. As increases, the initial shadow prices can take values further away from the optimal solution h * = 0. These initial values are representative of the expected values of shadow prices as one would like that shadow prices do not exceed prices. Table 1 presents the convergence status for each value of (1,000 random values where taken for each ). We observe that the iterative approach of Tranus fails to converge as the initial values get further away from the true solution.

Reducing the number of shadow prices, early results

We propose a model selection scheme to reduce the number of shadow prices needed to have a model reproduce the observed data. The classical approach is to iteratively modify the parameters until a perfect fit is achieved (with near zero cost function), when this is achieved the modeller will look at the values of the shadow prices as a quality test for the calibration. If the shadow prices are large, it means that the model has to compensate the various effects of the other economical parameters to attain a perfect fit. Thus, the modeller will tweak economical parameters (such as dispersion parameters) to maintain the perfect fit but with smaller shadow prices. A calibration will be completed when the model reproduces the observed data perfectly and the values of the shadow prices are small (for some economical sectors we will ask their variance to be small instead).

As there are as many shadow prices as observations we are trying to fit (there is one shadow price per observed production) the risk of overfitting is possible, which will in turn undermine the predictive capabilities of the model. What we propose, is sacrificing the perfect fit of the cost function, in order to lower the number of parameters to calibrate, particularly of shadow prices. To find the optimal trade-off between how many shadow prices we keep in the model and the desired value of the cost function is something that will have to be discussed with the community of modellers. What we propose here, is a simple model selection scheme that instead of having one shadow price per economical sector and zone, we keep only one shadow price per geographical zone. In this case, for the North-Carolina-1 model, as there are only 3 economical sectors for the land, we reduce by two thirds the complexity of the model. The calibration of the remaining third of shadow prices gave rise to a residual fit of the cost function of only 3% (ratio of residuals over observed productions).

Selecting the shadow prices to be kept in the model is easily done for the land sectors (nontransportable), because the prices (rents) are known. We achieved this, by computing the optimal shadow prices from (3.1) and setting to zero the "small" shadow prices. One can adjust the threshold used to declare shadow prices as being small to a desired level in the cost function, thus keeping more or less shadow prices. We also see in Figure 5 that the values of the shadow prices relatively to the prices have not seen a large increase.

Figure 5: The graphs show ratios of calibrated shadow prices over prices. Left: when all shadow prices are estimated (in this case, the fit of computed to observed productions is perfect). Right: here only one third of shadow prices are estimated, the others are set to zero. The fit is not perfect but good (3%, see text). Note that the scales of the two graphs are different. One can observe that on the right-hand side, shadow price to price ratios are not much larger than on the left-hand side, another indicator that it is plausible to exclude many shadow prices from the calibration.

We are currently working on testing this methodology for two time periods, to determine if noise in the base year productions could really propagate to the shadow prices and undermine the predictive capabilities of the model.

Results of Estimation of Shadow Prices and Substitution Parameters

We applied this procedure to a second real-scale LUTI model for North-Carolina (North-Carolina-2 model), with 38 zones, 3 floorspace and 9 other economic sectors. Figure 6 shows the shadow prices for all zones and floorspace sectors, after the two phases of our process, cf. section 4. After each phase, a global equilibrium of demand, production and prices, is achieved, however after the novel second phase, shadow prices are much smaller, meaning that the model represents reality much better (small ratios of shadow prices over prices is a crucial criterion used by practitioners to assess the quality of a Tranus LUTI model).

Conclusions and final remarks

The Tranus LUTI framework is a very powerful tool and the modelling possibilities are endless. However, the complexity of these large scale models is something that can not be underestimated, making the calibration and utilisation of these tools very expensive.

We have contributed with a reformulation of the land use module that simplifies the calibration process, exploiting the very basics of the mathematics that are behind the microeconomic models used, permitting the expert to incrementally calibrate the variables (from land use sectors to transportable sectors). The optimisation approach is more stable and clear than the classical approach, and enables the use of powerful optimisation algorithms currently available, solving the non convergence issues of the previous approach.

The procedure exposed for generating synthetic data is simple and straightforward, enabling us to try and benchmark our methodologies. We are currently preparing a set of benchmarks of calibrated models.

The proposed methodology for reducing the number of shadow prices needs additional fine tuning, but is a first step in what we consider a promising direction. We believe that the model "as it is" with one shadow price per observation bears a risk of overfitting.

Determining which shadow prices have to be removed may not be completely automisable, and the expert eye of the modeller has to have the last call.

Finally, the simultaneous calibration of different parameter types is a potentially very powerful tool. The results that we have for the North Carolina models have proven to be useful and saved many trial and error sessions. We would like to apply this idea of simultaneous optimisation to other "hard" to calibrate parameters, and we are working with modellers to identify them. A fully integrated and automatic calibration is our dream.

A Detailed mathematical description A.1 Land use economical sectors

If we consider land-use sectors, many simplifications are done to the equations of the model. Production is estimated using equation (3.5). The technical coefficient a mn i represents the demanded amount of sector n by sector m in zone i and is described by the equation (3.6). For the land use sectors, the consumption costs cn i are equal to the price paid: p n i + h n i , where the prices p n i are known. Then equation (4.1) is as follows:

S mn i = W n i exp(-ω mn a mn i (p n i + h n i )) p∈K m W p i exp (-ω mp a mp i (p n i + h n i )) . (A.1)
In the following, we will note U mn i = -ω mn a mn i (p n i + h n i ).

Derivative estimation:

Let us consider m and m as consuming sectors, n as land use sector and q ∈ K m ∂a mn i ∂h q i = -δ mn gap mn e -δ mn (p n i +h n i )

q = n 0 q = n (A.2)
The well known logit derivatives for S mn i are:

∂S mn i ∂h q i =            ∂U mn i ∂h n i S mn i -S mn i 2 q = n - ∂U mq i ∂h q i S mn i S mq i q = n (A.3)
if q = n, then ∂U mn i ∂h q i = 0, then for q = n:

∂U mn i ∂h n i = -ω mn [ ∂a mn i ∂h n i (p n i + h n i ) + a mn i ] = -ω mn [a mn i -δ mn gap mn (p n i + h n i )e -δ mn (p n i +h n i ) ] (A.4)
With these results, we can compute the gradient of the production function exposed in (3.5):

∂X n i ∂h q k = ∂ m (X m 0i + X * m 0i ) a mn i S mn i ∂h q k = m (X m 0i + X * m 0i ) ∂ ∂h q k [a mn i S mn i ] = m (X m 0i + X * m 0i ) [ ∂a mn i ∂h q k S mn i + a mn i ∂S mn i ∂h q k ] (A.5)
We can do the same computations for ω and get:

∂X n i ∂ω m q = m∈K n (X m 0i + X * m 0i )a mn i ∂S mn i ∂ω m q = (X m 0i + X * m 0i )a m n i ∂S m n i ∂ω m q (A.6)
where:

∂S mn i ∂ω mq = -a mn i (p n i + h n i ) S mn i -S mn i 2 q = n a mq i (p q i + h q i )S mn i S mq i q = n replacing (A.7) in Equation (A.6): ∂X n i ∂ω mq = -(X m 0i + X * m 0i )(a mn i ) 2 (p n i + h n i ) S mn i -S mn i 2 q = n (X m 0i + X * m 0i )a mn i a mq i (p q i + h q i )S mn i S mq i q = n A.

Transportable Sectors

The location probability derivatives follow the well-known logit derivatives:

∂P r n ij ∂h n k = -λ n β n P r n ij -P r n ij 2 k = j λ n β n P r n ij P r n ik k = j (A.7)
If we consider the optimisation function exposed in (3.4), we need to compute derivatives for p. From equation (2.9) the derivatives are as follows:

∂ pm i ∂h q k = n ∂a mn i ∂h q k S mn i cn i + a mn i ∂S mn i ∂h q k cn i + a mn i S mn i ∂c n i ∂h q k = ∂a mq i ∂h q k S mq i cq i + n a mn i ∂S mn i ∂h q k cn i + a mn i S mn i ∂c n i ∂h q k
and equation (2.8):

∂c n i ∂h q k = j ∂P r n ij ∂h n k (p n j + tm n ij ) if q = n is transportable 0 else as ∂P r n ij ∂h q k = 0 for q = n; ∂c n i ∂h n k = -λ n β n (P r n ik -P r n ik 2 )(p n k + tm n ik ) + λ n β n j =k P r n ij P r n ik (p n j + tm n ij )
The utility function is linear in the prices and the shadow prices, so the derivative estimation where the utility is involved are the same for h q k and p q k . In this case, we have to compute differently for the prices:

∂c n i ∂p n k = P r n ik + j ∂P r n ij ∂p n k (p n j + tm n ij ) if n is transportable 1 else
Finally, the partial derivatives for the complete cost function, are as follows:

∂f ∂h n k = ∂ ∂h n k [||X -X 0 || 2 + ||p -p|| 2 ] = ∂ ∂h n k [ n j (X n j -X n j0 ) 2 + n j (p n j -p n j ) 2 ] = 2[ n j (X n j -X n j0 ) ∂X n j ∂h n k + n j (p n j -p n j ) ∂p n j ∂h n k
They are all that is needed to apply standard non-linear least squares optimisers.

B New approach for computing the induced production of transportable sectors

Let us consider a transportable sector n. The demand for sector n in zone i follows the equation:

D i n = D * n i + m a mn i S mn i (X m i + X * m i )
Tranus is a very general model with a very large potential, but in practice many simplifications are made. For example, in the models we have considered, there is no substitution for transportable sectors. This condition simplifies the equation above, due to leading to S mn i = 1 for n. Also, the technical coefficient a mn i does not depend nor on the prices or shadow prices and is estimated externally. As we want to compute the demand for the base year production, we replace the X m i from above by the base year production X m 0i , thus transforming the demand into a constant that does not change between iterations (it depends neither on prices nor on shadow prices).

If we come back to the initial computations of induced production X from equation (2.7),

X n j = i D n i P r n ij
if we consider that the probability P r n ij is a function of φ n j = p n j + h n j , we can rewrite the calibration as an optimisation problem independently for each economic sector n as exposed in (3.9), more precisely:

min φ n 1 ,...,φ n J j |X n j (φ n 1 , . . . , φ n J ) -X n 0j | 2 (B.1)
Since this is an independent problem for each economical sector n, we will drop the index. If we name Φ the optimisation function of (B.1):

Φ(φ 1 , φ 2 , . . . , φ J ) = j |[ i D m i P r ij (φ j )] -X 0j | 2 .
Knowing that the demand D n i is constant, the derivatives are as follows:

∂Φ ∂φ k = 2 j ( i D m i P r ij (φ j ) -X 0j )D m i ∂P r n ij ∂φ k
where ∂P r n ij ∂φ k is analogous to (A.7). Once the optimal value φ * of the optimisation problem exposed in (B.1) is obtained, the prices can be extracted from the linear system (2.9). Finally the shadow prices can be computed subtracting the corresponding p n j from φ * n j . Due to the nature of the logit model, the value of the shadow prices is defined up to a constant. One could desire to fix the shadow price of one particular zone to zero for example, or what we actually do, just subtract the average value thus obtaining centered shadow prices.

C Synthetic Scenario Construction

For testing purpose, we propose a way of generating data sets that correspond to productions X equal to the inputed productions X 0 for a given set of shadow prices h 0 .

If we consider a fixed set of shadow prices h 0 , we need to obtain prices that are in equilibrium. To do so, we solve (2.9) until convergence is attained for the give set of shadow prices h 0 . We iterate until the value of X and p from one iteration to the next one remains constant, obtaining a pair ( X, p) that has attained convergence. It is important to notice that given h 0 is fixed, when p attains an equilibrium status, the production X is in equilibrium too. Doing so, is equivalent to solving the optimisation problem exposed in (5.2)

The first approach was for a given scenario, with in the X 0 case, the output X attained on the convergence of the system. This does not work, as the output value of X is calculated using the input value X 0 . More specifically, the attractor A n i of the logit formulation for the location probability (2.4) follows the equation:

A n i = W n i k b n k X k 0i α n
where W n i is the inputed attractor variable for zone i and sector n, b n k and α n are technical coefficients. As the computed production depends on the location probabilities, the computed X is a function of X 0 , the base year production. We need to modify the values of W n i as follows:

Ŵ n i = A n i k b n k Xk i α n
So, to set X as the equilibrium production, we rewrite both X 0 and W :

X 0 ← X W ← Ŵ
Doing this, we can have a scenario that perfectly reproduces the base year's production for a given set of shadow prices h 0 . In practice we use h 0 = 0.

D Numerical aspects

Local optimisation may converge to local minima. We observed this in practice, depending on the setting of the parameter β and on the starting point for the φ. An observation that seemed strange at first sight, was as follows. When estimating the φ for one sector, after convergence, the residuals of the cost function were all non-zero (besides for two zones for which observed production was zero). Further, all these residuals, besides for one zone, were exactly equal to one another and the residuals summed up to approximately zero. This seemingly strange behaviour has an explanation, as follows.

First, it must be noted that the sum of computed productions, does not depend on the values of the φ:

j X n j = j i D n i P r n ij = i D n i j P r n ij 1 = i D n i
If the data are coherent, then the sum of computed productions must equal that of observed ones:

j X n j = j X n 0j
Hence, the sum of residuals must be equal to zero, as was observed in practice.

The other issue concerned the fact that all non-zero residuals but one, were exactly equal to one another. This can be explained as follows. For one zone j, the value of φ j was sufficiently large at some stage of the estimation, so that the computed probabilities P r ij effectively became equal to zero, for all i: the absolute value of the argument of the exponential exp(-β(λφ j + t ij )) became so large that the exponential effectively got evaluated to zero. This in turn means that the computed production for that zone, also was computed as zero since j X n j = j i

D n i P r n ij 0 = 0
Hence, the residual for zone j is non-zero, and actually equal the (negative of the) observed production X 0j . Since the sum of residuals over all zones must equal zero, as shown above, we must have:

k =j (X k -X 0k ) = -X 0j
Remember that the cost function to be minimised is the sum of squared residuals; as for the zones other than j, this means:

min k =j (X k -X 0k ) 2
It can be shown that given the constraint that the sum of residuals must equal a known value, the cost function is a minimum if that known value is equally apportioned to the residuals, i.e. if all the residuals are equal to that value, divided by the number of residuals:

∀k = j : X k -X 0k = - X 0j k =j 1
This explains the observation made in practice, described above.

Figure 1 :

 1 Figure 1: Schematic overview of Tranus.

Figure 2 :

 2 Figure 2: Sketch of computations in the land use and activity module.

Figure 3 :

 3 Figure 3: Contour plot of equations (5.3). From left to right: no solution, one solution and multiple solutions.

Figure 4 :

 4 Figure 4: Plot of cost function for a given pair (h 1 1 , p 1 1) near the optimal value (0, 2.676).

Figure 6 :

 6 Figure 6: Ratios of shadow prices and prices after phase 1 (left) and 2 (right). Note the different scales of the graphs.

Table 1 :

 1 Comparison of calibration algorithms for the North-Carolina-1 model. Shown are the percentage of random trials for which the algorithms converged to the correct solution.

	value:	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
	Tranus software	100%	100%	100%	63%	3%	0%	0%	0%	0%	0%	0%
	Our algorithm	100%	100%	100%	100%	100%	100%	100%	100%	100%	100%	100%

http://www.tranus.com/tranus-english

Acknowledgment: This work is supported by the CITiES project (ANR-12-MONU-0020).