A Dichotomic Algorithm for Transportation Network and Land Use Problem

Redeveloping sites to accommodate housing or new economic activities is a major urban policy challenge. This problem belongs to the class of NP-hard problems. In this paper, we present an attractive mixed integer nonlinear programming formulation for the transportation network and land use problem. We first introduce a new useful nonlinear formulation of this challenging combinatorial optimization problem. Then, an alternative to considering a linearity of the constraints is to reformulate the problem as a new exact, compact discrete linear model. The problem is solved by our new algorithm and numerical results are presented for a number of test problems in academics instances.


Introduction
The Transportation Network and Land-Use (TNLU) problem is a combination of the transportation network optimization problem and land use planning or Quadratic Assignment Problem (QAP).This model appears, for example, if we decided to find locations or relocations of some a city amenities (home, shop, work or leisure places), that may reduce travel time or travel distance.By solving this problem we want to have the best layout of the city in the direction of the location of activities and roads linking these activities.A better configuration of activity locations and transportation for a city would reduce congestion and pollution.For example, we can consider the Koopman and Beckmann (1957) version of the problem, which can informally be stated with reference to the following practical situation: a municipality must decide to assign n activities to an equal number of locations and want to minimize the total cost of location and transportation.For each pair of activities (i, j) a flow of communication f i j is known, and for each pair of locations (l, k) the corresponding distance L kl is known.The transportation cost between activities i and j is f i j L kl , where i and j are assigned to locations k and l, respectively.The objective of the municipality is to find (i) an assignment that minimizes the sum of the total linear cost for installing an activity to a location, (ii) the total quadratic interaction cost, (iii) the flow between the activities and (iv) total cost that link roads between activities.
Since the 60s many authors are interested in these two problems separately: for example Scott (1969), Boyce, Farhi andWeischedel (1973)andHoang Hai Hoc (1973) for network design problems; Lawler (1963), Maniezzo (1997), Xia and Yuan (2006), Huizhen Zhang (2010) and many others for the activities's location problems.Marc Los and more recently Lin and Feng are interested in the combined problem.They have considered a more complex model in the sense that their formulation allows the location of several activities in the same area.Marc Los gave an exact solving method and heuristic methods.The Los's model take into account to the OD flow, once the activities have been fixed, via the shortest paths in the sense of a fixed cost independent of the flow on the arcs (see Gueye (2012)).For more details or a review of the literature, we refer the reader to Baldé and Ndiaye (2016), Duranton and Puga (2015), and Gueye (2012).The TNLU can be formulated as a Mixed Integer NonLinear Programming (MINLP) problem.Many algorithms and softwares (as Cplex (2016), Gurobi (2016), Couenne (2018), etc.) are proposed in the literature to solve non linear (quadratic) programmming problem, but present some limitations for large number of variables.The quadratic assignment problem constitutes the natural formulation of a large number of concrete problems, belonging to various fields.Its complexity is such that the proposed methods so far to solve it were not of sufficient effectiveness.
In this paper, we propose finite linearization and an algorithm to find good solutions with fewer iterations and running times.It is a new resolution method based on a dichotomic translation of a hyperplane and can give better results compared to those obtained by other methods of solvers.We refer the reader to see first Lavallée, Ndiaye and Seck (2011); and Ndiaye, Lavallée and Seck (2013) for an application of the algorithm to the Klee and Minty problem, for example.As the algorithm is designed to solve a linear programming problem, we first linearize the TNLU problem and introduce (arbitrary constraint) support hyperplane guaranteed by the objective function.In other words, we add in the constraint a hyperplane parameterized by real α 0 , which serves as a barrier.Secondly, to start with the resolution, we set an initial condition and a stopping criterion.To verify the obtained solution of the senegaulois algorithm, we compare it with solutions given by Cplex and Gurobi (we give to Cplex and Gurobi the quadratic problem).
The paper is organized as follows.In section 2 we present definition and assumptions.In section 3, we give the formulation of the TNLU.In section 4 we introduce a new linearization for the TNLU.In section 5 we develop a brief presentation of senegaulois algorithm.The section 6 provides numerical simulations using test problems in academics instances.Finally, in section 7, we provide conclusions and suggestions for future works.

Definitions, Assumptions and Theorems
Let us formulate the main concepts used in the paper, mainly for sections 3 and 5.
respectively the set of vertices of the graph G 1 and G 2 , and E 1 and E 2 represent respectively all edges (or arcs) of G 1 and G 2 .
The graph matching plays a central role in solving correspondence problems in computer vision.Graph matching problems that incorporate pair-wise constraints can be cast as a quadratic assignment problem (QAP) (see Koopmans and Beckmann(1957)).Unfortunately, (QAP) is known to be NP-hard (see Sahni and Gonzalez (1976)) and many algorithms have been proposed to solve different relaxations.
Definition 2.2 (Exact Graph Matching).The graph matching problem can be stated as follows: Given two graphs When such a mapping π exists, this is called an isomorphism.
A bijection of a finite set of cardinal n on itself is a permutations of n elements.We denote by S n the set of permutations de {1, ..., n}.
S n = {π : Given an activity i ∈ E(G 1 ), then π(i) = k implies that activity i (i − th component in the ordered set N 2 ) is at location k.In the optimization problem, we are looking for the best assignment, i.e., we have to optimize some objective function which depends on the assignment π.Assignments can be represented in different ways.The bijective mapping between two finite sets N 1 and N 2 can be represented in a straight forward way by a perfect matching in a bipartite graph G = (N 1 , N 2 ; E), where the vertex sets N 1 and N 2 have n vertices.Edge (i, k) ∈ E is an edge of perfect matching iff k = π(i).
Definition 2.4 (Final graph).Let G = (N, E) be the general graph.G is determined by the data of two sets: • a non-empty finite set N whose elements are called vertices The arcs ((i, k), ( j, l)) for all i j, k l of G are connected if starting from the nodes i to the zone k = π(i) it is possible to reach the node j at the zone l = π( j).We can see, in the Figure 1 the location of the set of activities on the zones as a bijective mapping π(.) of all activities to all zones.Now, let's define two theorems dealing with the question of existence of a hyperplane that separates two given disjoint convex subsets.
(where C ∞ 1 and C ∞ 2 are asymptotic cones of C 1 and C 2 , respectively).Then, it is possible to separate C 1 and C 2 .There exists a vector ξ ∈ E such that: and C 2 are strictly separated by a hyperplane.

Localization area and network design
Theorem 2 (Hahn-Banach II).Let E be a finite-dimensional vector space with inner product ⟨., .⟩.Let C 1 , C 2 ⊂ E be two non-empty disjoint convex subsets.Then, there exists a nonzero vector ξ ∈ E\{0} such that : • Let C be a non-empty closed convex set.∀u ∈ ∂C (the border of C), C is supported at u.If C is a convex set of a topological vector space, any point of the border of C is in a supporting hyperplane.
where p is the set of the half-space containing C.

Indices, Sets and Parameters
• N 1 : set of activities (= {1, ..., n}), • N 2 : set of zones (= {1, ..., n}), • N: set of assignment, set of arcs between each activity pairs, • E 2 = {(k, l) : k, l ∈ N 2 , k l}: set of arcs between each pairs of the area, • F = ( f i j ) n×n : the flow matrix volume from activity i and j, • D = (d i j ) n×n : the distance matrix between location i and j, • C = (c ik ) n×n : location cost matrix of the activity i on the area k, • B = (b i j ) n×n : construction cost matrix of "build" the link or road (i, j).

Variables
The variables of the model can be divided into three groups.The first group includes the binary variables x ik .The second group y i j represents the network configuration.The last group contains the continuous variables L kl that determine the length of the path k to l.
• Let us start with the observation that every permutation π of the set N = {1, ..., n} can be represented by an n × n matrix X π = (x i j ), such that Matrix X π = (x i j ) is called a permutation matrix and caracterized by following assignment contraints: x ik ∈ {0, 1}, for all (i, k) ∈ E.
The 3 and 4 equalities of assigning each element of N 1 to an element of N 2 .
• To determine the existence of a path in the network (whether or not there is a direct road link between i and j), we define binary variables (boolean) y i j here below: • Then it is necessary to introduce variables of shortest path L kl between the location k and location l and variables z kl i j which take value 1 if the shortest path between k and l goes by link (i, j) and 0 otherwise.We have: If we consider the distances of the arcs (i.e. between each pair of activity) as the costs associated with the flow on the arcs, then the problem of the shortest path from k to l can be formulated as a minimum-cost flow problem (FCM) to move a flow unit from the activity i to k to the activity j to l.
The out-degree of i: is the number of arcs that have i as an initial vertex (e.g. the number of arcs coming out of i).
The in-degree of i: is the number of arcs that have i as a terminal vertex (e.g. the number of arcs coming into i).
The minimum-cost network flow problem (FCM) is a fundamental problem in network analysis that involves sending flow over a network at minimal cost.Let G 1 = (N 1 , E 1 ) be a directed graph.For each link (i, j) ∈ E 1 , associate a cost per unit of flow, designated by d i j .
(FCM) : min where, o(k): is origin (area of origin) and d(l): is destination (area of destination) .
The shortest-path problem (SPP) (Kelly and O'Neill (1991))is to find the directed paths of shortest length from a given root node to all other nodes.Moreover, the (SPP) is also a special case of the (FCM), where the objective is to find the minimum distance or length between two given nodes (e.g., nodes k and l).Here, the values of the n × 1 vector b are now restricted only to 0, 1 or −1.The mathematical formulation of the (SPP) is as follows: The optimum solution to this problem sends unit flow from the root to every other node along a shortest path.

Transportation Network and Land Use (TNLU) Problem
This problem is a combined model composed of the following subproblems: an activity localization problem for the optimal location of activities; a problem of configuration and dimensioning of the transport network to define a transport network or its evolutions, considering a set of potential road, having a cost of construction and use, according to the observed flows.The objective is to assign each activity to a location such that the total cost is minimized: The constraints 3.14 and 3.15 are the so-called assignment polytope or assignment constraints (see Burkard, Cela, Pardalos and Pitsoulis (1999)), they make it possible to ensure that each activity is assigned in a zone and each zone accommodates one and only one activity.Equalities 3.16 and 3.17 define the origin and destination of the path.The constraint 3.18 is the conservation constraint at each node.The constraints 3.19 ensure that the sum of the distances of the activities assigned to the different zones does not exceed the maximum distance (or length) L kl between the zones.The constraints 3.20 make it possible to define the roads to construires or not.The 3.21 and 3.22 constraints ensure that the decision variable is binary.The 3.23 constraints ensure that the decision variables are continuous.

Linearizations of the QAP
To linearize the objective function we introduce new variables and linear contraints for the QAP to eliminate the quadratic term of 3.13.This methods transform the problem into a mixed linear program (see Burkard, Cela, Pardalos and Pitsoulis (1999)).The first linearization, proposed independently by several authors Fortet (1960), Watters (1967) , involved the addition of one binary variable and two linear constraints for each product of two variables.In 1974, Glover and Woolsey (1974) made a major breakthrough by proposing a linearization where the additional variable was not required to be explicitly defined as an integer.This linearization will be referred to, in this paper, as the standard linearization.Consider the problem of assigning a set of activities to a set of locations, with the cost being a function of the distance and flow between the activities, plus costs associated with a facility being placed at a certain location.Hence the name Quadratic Assignment Problem by Koopmans and Beckmann (1957): x ik , y i j ∈ {0, 1}, Here the feasible region X is the assignment polytope.

Standard Linearization (Bin × Bin)
Standard Linearization (S L) consists of linearizing the product of two binary variables (see Glover (1975)).The principle of "standard" (or "classic") linearization is to replace remplace x ik x jl par X ik jl .Then add the three inequalities below: The problem (QAP) is then rewritten thus adding to the assignment constraint the system 17: (QAP) : subject to 15, − 16. (18)

Product Linearization (Bin × Bin)
We always consider the formulation of Koopmans and Beckmann (1957).We multiply both sets of assignment constraints 4.3 by each problem variable x ik , ∀ (i, k) ∈ E. We then linearize the resulting products by requiring X ik jl = x ik x jl , ∀ ((i, k), ( j, l)) ∈ E. We obtain the following linearization constraints (see Blanchard, Elloumi, Faye and Wicker (2002): The reformulation of the problem (QAP) is given as follows: (20)

Linearizing the Product of a Binary and a Continuous Variable (Bin × Cont)
In the objective-function, the product L kl X ik jl is replaced by a non-negative variable Z ik jl which will take the same value as L kl when X ik jl = 1, and will be zero when X ik jl = 0.If L kl is bounded below by zero and above by L max (or any BigM), i.e. 0 ≤ L kl ≤ L max .Suppose now that we have Z ik jl = L kl X ik jl , where L kl is a continuous variable and X ik jl is your binary variable.Add the three inequalities below (Coelho (2018)) Note that if X ik jl is zero, than the first inequality ensures that Z ik jl will be zero as well (note that the third inequality only states that Z ik jl has to be greater than a negative number) .On the other hand, if X ik jl is 1, the first inequality ensures that Z ik jl is less than our BigM = L max , which is further tightened by our second inequality.The last inequality states then that Z ik jl has to be greater than or equal to L max , exactly as we wanted.Finally, if L max is not bounded below by 0, but has bounds [0, L max ] (with L max assumed to be positive due to the last inequality), then the form of the linearized inequalities are the following:

Binary to Continuous Variables
Consider a problem (IP) 0−1 involving a binary variable x ∈ {0, 1}.This can be reformulated as follows (see Billionnet, Elloumi and Lambert (2013)): (IP) 0−1 : Since a binary variable x i ∈ X can only take values in {0, 1}, any univariate equation in x that has exactly x = 0 and x = 1 as solutions can replace the binary constraint x i ∈ {0, 1} n .The most commonly used is the quadratic constraint The linear reformulation of (IP) 0−1 consists of asking X ii = x i x i , ∀ i, and then using standard linearization 4.1.1.
Additionally, for all Constraints (a) and (b) imply X ii = x i .Consequently, contraints (c) and (d) become 0 ≤ X ii ≤ 1.
This process would reduce all binary problems to nonconvex quadratically constrained problems.We refer the reader to Liberti, Cafieri and Tarissan (2009) for more details.
Finally, we obtain the continue linearization problem:

The Senegaulois Algorithm
We consider a linear programming problem below : x ≥ 0 where x represents the vector of variables (to be determined), c and b are vectors of (known) coefficients, A is a (known) matrix of coefficients, and (.) T is the matrix transpose.Recall that convex polyhedron on E is a set P as: Now we briefly recall some concepts, such as the Hahn-Banach theorems.The proofs can be obtained in Achmanov (1984), and in Brezis, Ciarlet and Lions (1999) for polyhedron.
Using the separation theorem, every such hyperplane H divides the whole space in two convex sets: known as halfspaces.One part contains the polytope K defined from the constraints in 1.
Then, the following cases hold: 1. Hyperplane H does not divide the whole polytope K of constraints, then H is outside the feasible solutions, i.e.H ∩ K = ∅ (H ∩ K does not contain any solution).In such case, it is necessary to choose another hyperplane H ′ // H (parallel) such that H ′ lies between the origin 0 and H.
2. Hyperplane H divides the whole polytope K of constraints, then there is a non-empty intersection between H and K, i.e.H ∩ K ∅.In this situation, the problem contains some feasible solutions and K is separated by H.If the stopping test fails, then a new separation hyperplane H * // H is chosen such that H * lies in the half-space which does not contain the origin 0.
3. Polytope K is empty, the problem has no solution.

The Algorithm
Suppose H : f (x) = α and K the initial polytope defined by the set of constraints of problem 5.1.If H ∩ K = ∅, then we state α := α 2 .That is to say that we divide the gap between 0 and α into two parts.It is the main idea of the algorithm.
Now suppose that we have two values of α denoted by α 0 and α 1 such that: MILP1: Mixed Integer Linear Problem with continuous variables Conservation of flows for each pair (k, l) The variables • H : In addition, two cases are possible : 1. H : f (x) = α and K ∩ H = ∅, then put α 1 := α 2. H : f (x) = α and K ∩ H ∅, then put α 0 := α And so on.
In this algorithm the hyperplane is defined by: where n ∑ i=1 c i x i is the objective function (the one to maximize in our purpose) of the LP.Remark: From a geometrical point of view, at each step of the algorithm, the new hyperplane is obtained from the previous one by a translation.Thus, we search vectors that characterizes polytope H defined by ( 27); then we search an unitary vector − → v as vectorial step.

The Stopping Test
As stated above, our algorithm works with successive approximations but the result is obtained with the desired accuracy1 .That is to say that it is the operator who gives the acceptable tolerance.Suppose that margin tolerance is ϵ ∈ Q, with equation notations in 26, then the stopping test is :

To Provide a Solution
When both H ∩ K ∅ and | α 0 − α 1 |≤ ϵ are true, a feasible solution is obtained in the part of the polytope as we can see it in Figure 2.
Figure 2. The dichotomic procedure Geometrical explanation: Unless we have obtained the optimum, the hyperplane H 0 : f (x) = α separate the other constraints.Basis are vertices of the residual polytope.There are n such vertices associated with the separating hyperplane satisfying the constraints.But all interior points of the residual polytope satisfy the constraints and are in the permitted interval.All convex combination of these points is a valid solution2 .
More details about the algorithm can be found in Lavallée, Ndiaye and Seck (2011); and Ndiaye, Lavallée and Seck (2013)

Implementation and Numerical Results
In this section we present computational studies on the linearization approaches and the senegaulois algorithm for the TNLU.

Implementation
Aiming to obtain the best possible bounds, we solve the linearization model globally with the senegaulois algorithm and the mixed integer model by Cplex (2016) andGurobi (2016).Simulations are done on a CPU@ 3.40 GHz Intel(R) Core(TM)i7-370 HP-Pro3500 with 16 GB of main memory running Linux ubuntu.The termination criteria, i.e. the difference between the lower and upper bounds in our algorithm, is set to 10 −5 .
To model the TNLU problem, we use a modeling tool known A Mathematical Programming Language (AMPL) (2018), version 20130609.Developed at Bell Laboratories, AMPL is a high level, an algebraic modeling language for linear and nonlinear optimization problems.Besides its ability to model optimization problems, AMPL also helps in comprehension and completeness of the model logic, due to the similarity of its syntax to algebraic notation.It provides for an easy way to express common linear programming structures (e.g., flow constraints).Although AMPL allows us to mathematically model the optimization problem, we use two external solvers such as Cplex (version 12.5.)and Gurobi to solve the problem.
The entry in AMPL includes a template file (*.mod), a data file (*.dat), and a runtime program (*.run).The AMPL groups them in a format that the solver understands.Our algorithm based on the dichotomous search as an AMPL script, is implemented in the (*.run) file.The AMPL execution program calls the script file (*.run) which provides the mathematical description of the model and the data into AMPL.Then, it passes this instance of the problem to the solver, which in turn, resolves the instance, and makes the solution to AMPL.

Numerical Results
In order to access the convergence of the senegaulois algorithm presented in section 5, we have applied it to eight examples.The obtained solutions, with our algorithm, are compared with those obtained with Cplex and Gurobi.The Table 1, with 8 test problems (tnlu1 to tnlu8), shows the result of the computational comparison and the performance of each resolution technics, where: • Inst: is the name of the test problem, • NbVar: the number of variables, • OptVal: the obtained optimal value, • Times(s): the execution time in seconds, • Iter: the number of iterations, • α 0 : the initial value for the senegaulois algorithm.
We see that the senegaulois algorithm provides good optimal value for the TNLU problem.For some test problems on academic instances, the number of iterations is quite small.For instances tnlu1, tnlu4 and tnlu8 the iterations are quite small compared to Cplex and Gorubi.However, for all instances, the obtained Optimal value remains nearly the same for the 3 methods.The Figures 3 and 4 show the evolution of the number of iterations and computational time of each algorithm, respectively (Cplex, Senegaulois and Gurobi).From these figures can see that iteration and time frequencies are maintained in senegaulois algorithm.
In addition, to evaluated the results we use the free software R project (2018), and display the 8 observations and the 11 variables:  We firt compute only the average of 2 obervations (times and iterations), then the number of variables and the firt value of the senegaulois algorithm.
On average, the same upper bound is appreciably found, but the running time and the number of iterations of Cplex and Gurobi remain large (see Table 2).Our numerical experiments show convergence with a fewer number of iterations and computational times.

Conclusion and Future Works
This study aims at describing and simulating a new type of optimization problem for the TNLU problem.Our algorithm for TNLU was presented and applied to problems from the literature designed to highlight some intrinsic difficulties of TNLU.The Senegaulois algorithm, based on a dichotomic search is compared to two programming solvers such as Cplex and Gurobi.We investigated the average of iterations and computational times and established interesting results characterizing some instances.In future work, we plan to study the combination between the senegaulois algorithm and genetic or greedy algorithm.

Table 2 .
Average of the values.