Continuous transportation as a material distribution topology optimization problem

The problem of moving a commodity with a given initial mass distribution to a pre-specified target mass distribution so that the total work is minimized can be traced back at least to Monge’s work from 1781. Here, we consider a version of this problem aiming to minimize a combination of road construction and transportation cost by determining, at each point, the local direction of transportation. This paper covers the modeling of the problem, highlights how it can be formulated as a material distribution topology optimization problem, and shows some results.

play also on a smaller scale, for instance in agriculture or forestry, when temporary or otherwise designated roads must be paved. For general consideration of distribution of matter, it is natural to employ a fluid dynamics formalism. An early attempt to do so was made by Beckmann (1952), who minimizes transportation cost under a conservation of matter constraint. Beckmann seems to have been unfamiliar with the Monge-Kantorovich problem at the time, but there is a close link between the formulations, noted by Igbida (2013), among others. For a survey of continuous transportation modeling (see for instance Puu 2009).
The literature on the Monge problem and its derivation is extensive (albeit small compared to network counterparts), but considers mainly the theoretical properties of the problem formulations and its solutions. Much less seems to have been written regarding solution methods applicable for practical transportation problems. The intention of this paper is to bridge this gap so that practically useful algorithms can be formulated ultimately.
In a general setting, the objective is to find a transportation map T that optimally rearranges the measure μ + into μ − as detailed below. Let T be the set of all transportation maps T pushing μ + onto μ − , that is, μ + (T −1 (B)) = μ − (B) for each Borel set B ⊂ . Given a cost function c that represents the unit transportation cost, the so-called (generalized) Monge problem is then to find inf T ∈T c(x, T (x)) dμ + (x).
( 1 ) This problem is in general ill-posed. In 1942, Kantorovich suggested a relaxed version of the transport problem that always admits a solution (Kantorovitch 1958). More precisely, the Monge-Kantorovich problem searches for a so-called transportation plan γ ∈ M that minimizes × c(x, y) dγ (x, y).
( 2 ) Here, M is the set of all non-negative Radon measures γ on × that have μ + and μ − as marginals. That is, γ (B × ) = μ + (B) and γ ( × B) = μ − (B) for each Borel set B ⊂ . In essence, for each pair of subsets B i ⊂ and B j ⊂ , γ (B i × B j ) can be viewed as the amount to be transported from B i to B j . During the last two decades, there has been a renewed activity in studying these transportation problems (Ambrosio and Pratelli 2003;Braso and Petrahe 2014;Evans and Gangbo 1999;Igbida 2009). Evans and Gangbo (1999) made a breakthrough, showing existence of a transportation map provided that c is the Euclidean distance, μ + = f + dx, and μ − = f − dx, where f + and f − are given functions that specify the initial and target mass distribution and satisfy supp{f + }∩supp{f − } = ∅. However, the numerical treatment of the problem has this far received little, if any, attention.
Here, we aim to solve numerically a continuous transportation problem by using material distribution based topology optimization. The rationale for this is that road design is effectively a material distribution problem, and transportation is nothing but flow of matter. The material distribution method was originally introduced by Bendsøe and Kikuchi (1988), who considered optimization problems in solid mechanics. The monograph by Bendsøe and Sigmund (2003) gives a comprehensive presentation of material distribution based topology optimization techniques and highlights many applications. The review by Sigmund and Maute (2013) gives an overview of different approaches to topology optimization and compares their strengths and weaknesses. The "flagship problem" for topology optimization is to minimize the compliance of an elastic body subject to a constraint on the amount of available material. There are many well-established techniques to solve this problem and some of these are currently used in the design process of advanced components in the automotive and aeronautical industries. Today, large-scale material distribution problems have many millions or even billions of design variables Schmidt and Schulz 2011;Wadbro and Berggren 2009). For other problems, such as those with state constraints or those concerning wave propagation, the methodologies are still maturing. During the last decades, much work has been focused on developing the methodologies (Le et al. 2010) and extending the ideas to other fields, such as fluid flow (Borrvall and Petersson 2003), fluid-structure interaction (Andreasen and Sigmund 2013;Yoon 2010), and acoustic (Dühring et al. 2008;Wadbro 2014;Wadbro and Berggren 2006) as well as electromagnetic (Aage and Johansen 2017;Andkjaer and Sigmund 2011;Erentok and Sigmund 2011;Nomura et al. 2007;Hassan et al. 2014) wave propagation. Although various kinds of systems including transportation processes (thermal, electric, convective, etc.) have been targeted for topology optimization in the past, the application of topology optimization to continuous logistic transportation models is, to the best of our knowledge, novel. The contribution closest in style to our work is that of Ryu et al. (2012), who used topology optimization techniques to attack a path planning problem for a robot that moves in an environment with obstacles.

Problem description
In this section, we present a model that accounts for the road construction cost as well as the transport cost of a single commodity in a steady state setting. Since, the model is new and the modeling includes many different symbols, we start by providing the following nomenclature list: x Position (x = (x 1 , x 2 )) ρ Density distribution of commodity q Rate of production/consumption of the commodity u Transport velocity v Transport speed (|u| = v) Potential κ Conductivity α Road design s A function specifying the relation between α and v F Set of feasible conductivities A Set of feasible road designs J T Transport cost J T -relaxed, differentiable version of transport cost J R Road cost β Tradeoff parameter.

Transportation of a product
On a microscopic level, transportation of goods is a discrete process where point charges (representing trucks, backhaulers etc.) move along certain paths (on roads or offroad) between sources and destinations. On a macroscopic level, and seen over long periods of time, it is relevant to model transportation as a continuous process in the same way as continuous fields are used to describe electric current or gas flow that is in reality movement of microscopic particles. This view on transportation enables the use of a mathematically tractable formalism from continuum mechanics. The model we propose is similar in spirit to that of Beckmann (1952), who also considers flow under a conservation of matter condition. While (Beckmann 1952) set the theoretical foundation for continuous transportation modeling with the aim of minimizing transportation, our aim is to consider also road construction cost, and to employ topology optimization in order to find solutions for general scenarios.
As our model problem, we consider a commodity that is produced or consumed at the space dependent rate q [kg m −2 s −1 ] and transported with velocity u [m s −1 ], where |u| = v, and v : R 2 → R + is a space dependent transportation speed. Denote by ρ [kg m −2 ] the instantaneous density field of the commodity. Moreover, we assume that the production, transport, and consumption of the product are all confined to be inside a region ⊂ R 2 . That is, we assume that is large enough so that ρ| ∂ ≡ 0. In the numerical experiments, is selected to be the unit square. Figure 1 shows an image of the transportation scenario according to the continuous model. By conservation of mass, we have that ∂ρ ∂t + ∇ · (ρu) = q in for all times t.
Here, we assume that the system is at steady state, so Furthermore, since we have neither inflow nor outflow of the product over the boundary ∂ , we also have that of the product in . This amount can be viewed as the total amount of the product in storage, en route from producer to consumer. Here, the total transportation (or storage) cost is proportional to the total mass of product in transportation at steady state. Hence, the total transportation cost is where c T is the unit transportation cost, equivalent to the sum of the capital cost of the commodity, and the cost of keeping the commodity in transportation. It is assumed that transportation can be viewed on a macroscopic level, so that the unit transportation cost is independent of ρ.
To make the transportation as efficient as possible, we are thus interested in finding the transportation velocity u that minimizes J T . To find the optimal transportation velocity u ∈ U = u ∈ L ∞ ( ) 2 | |u| = v a.e. in given the transportation speed v and the production-consumption rate q satisfying condition (5), we may solve min u∈U J T subject to steady state mass conservation (4) and that the density ρ ≥ 0 in . As a consequence of the conditions above, not all u ∈ U are feasible; after all, the overall velocity field u must be directed from the production area to the consumption area and not in the reverse direction. Formally, the unknown in optimization (8) is the transportation velocity u, but since the transportation speed v is given, we are effectively solving for the local transportation direction (u/v).

Potential approach
Below, we assume that we are given the productionconsumption rate q and are interested in finding the best transportation strategy to minimize J T . Given q that satisfies condition (5), we still need to find a density ρ and a velocity u to satisfy (4). Using an analogy from electromagnetism, we restrict the set of feasible transportation velocities by assuming that there exists a potential and an associated conductivity κ ≥ 0 so that Here, can be interpreted as an economic potential that drives the commodity flow. The flow rate ρ|u| = ρv will then be proportional to the conductivity κ and ∇ . Through the layout of κ, it is possible to guide the flow. The steady state mass balance (4) can thus be written as Remark 1 The formulation above is similar to the one used by Evans and Gangbo (1999), who proved the existence of a transportation map provided that transportation cost between two points is equal to the Euclidean distance between these points (essentially |u| ≡ 1 in our setting).
In particular, they proved that this map can be obtained by building a flow by solving an ODE involving λ, ∇z, and q, related through the Euler-Lagrange formulation of Kantorovich's dual problem where z is the dual variable and λ is a Lagrange multiplier for a constraint of the type |∇z| ≤ 1.
Note that, by construction, the conductivity will at best be determined up to a constant by state equation (10). For any reasonable transportation map, the value of the potential far away from the area of interest will neither affect the transportation strategy nor the total amount of the product in the system. To avoid the ambiguity of only having the potential determined up to a constant, we choose to add the boundary conditions where D ⊂ ∂ is a small boundary portion with strictly positive measure (| D | > 0). D can be interpreted as a "shunt to ground" that annihilates the effect of numerical deviation from condition (5). A variational form of (10), (11), and (12) is: Having solved variational problem (13), we can identify The restriction of the feasible transportation velocities provided by using potential and conductivity κ ensures that we have a velocity-density pair (u, ρ) that satisfies mass conservation (4) and the condition ρ ≥ 0 in . The problem of determining the best conductivity distribution given a road design can hence be written where F is the set of feasible conductivities and solves variational problem (13). Note that by scaling κ by a positive factor and scaling the potential by the inverse of the same factor, the value of the objective function remains unchanged. Hence, we can without loss of generality define the set of feasible conductivities as Variational problem (13) may not have a unique solution for any κ ∈ F 0 . A standard procedure in topology optimization for structures, which we also will employ here, is to modify the lower bound for the conductivities. More precisely, we define the set of feasible conductivities as where κ is a small, strictly positive number. If the transportation speed v satisfies v min ≤ v ≤ v max , for some strictly positive constants v min and v max , then the bilinear form associated with variational problem (13) is coercive.
To avoid the non-differentiability of the absolute value at 0, we approximate the objective function in problem (15) by where is a small positive parameter ( = 10 −8 for the numerical experiments in this paper) and solves variational problem (13). To conclude, the optimization problem that we solve with the aim of finding the best transportation direction at each point in by determining the conductivity distribution given a road design is

Adding road design
The main aim of this study is to determine the best road layout and transportation direction to minimize a combination of transportation cost and road construction cost. To describe the road layout, we use a so-called material indicator function α defined so that α = 1 where a road is present and α = 0 else. The set of admissible road layouts is Throughout this paper, we assume that the transportation speed v is directly determined by α; the transportation speed is higher when transporting goods on roads than when transporting goods off-road. That is, the transportation speed only depends on whether the ground is covered by roads or not; spatial properties and, for example, slopes, their directions, and which material that covers the ground does not influence the transportation speed. Hence, the transportation speed is given by the composite function where the speed function s : [0, 1] → R + . Moreover, we assume that the road construction cost J R is linear in the road design, that is where c R is the unit road cost. Henceforth, we assume that the unit transportation cost and the unit road cost both are 1; that is, c T = 1 [kg −1 ] and c R = 1 [m −2 ]. To study the tradeoff between transportation cost and road construction cost, we minimize a weighted sum of the road construction cost and the transportation cost. That is, we formulate our main optimization problem as

Numerical treatment
We solve variational problem (13) The discretized variational problem (24) can be written in matrix form as where the components of the N × N stiffness matrix K and the N × 1 right hand side vector f are respectively. Similarly as for the conductivity, we use a piecewise constant function α h to define the road layout in , and consequently we approximate the transportation speed by a piecewise constant function v h . For future notation, we let κ = (κ 1 , κ 2 (18) numerically, we use the composite midpoint quadrature rule to obtain the approximation where x c m is the midpoint (or center) of element E m and h solves variational problem (24) with the conductivity distribution associated with κ. Note that the discretized transportation cost J h T is an approximation of the differentiable transportation cost J T and depends on parameter . Further, since we have approximated the road design by a piecewise constant function, the discretized road construction cost is In our conceptual optimization (23), the road design α is constrained to only attain the values 0 and 1 almost everywhere. Design optimization problems with this type of binary constraints are often associated with various issues. From a mathematical viewpoint, these problems are typically ill-posed in that the problem lacks solutions within the set of feasible designs; and from a computational point of view the problem is a large-scale non-linear mixed integer program and hence computationally intractable. We remark that for real-life problems there are often requirements, such as a minimal width of structural members of the design, which assures that the problem has solutions within the space of admissible designs.
Here, we choose to attack the problem by following a standard procedure applied for material distribution topology optimization. First, we relax the constraints on the road network design to the continuum [0, 1]. That is, we replace the set of feasible road network designs A by Since the conductivity already was allowed to take values in the continuum [κ, 1], we simply replace the set of feasible conductivities F by its elementwise constant counterpart This far, we have not discussed how the speed function s is selected. In principle, the replacement of A by A h likely gives rise to non-binary optimal road network designs. To promote binary designs, we will use the socalled SIMP (Solid Isotropic Material with Penalisation) approach (Bendsøe 1989), which is typically employed in material distribution topology optimization for linear elastic structures to suppress intermediate densities. Accordingly, we define the transportation speed function as s( where v min and v max are the transportation speeds outside and on roads, respectively, and p ≥ 1 is a so-called SIMP-penalty parameter. Hence, for m = 1, 2, . . . , M, the transportation speed in element E m is The idea of SIMP in this context is to let intermediate values of α contribute more to the road construction cost compared to the relative increase that they provide in transportation speed. When adding a penalty term to suppress intermediate values in material distribution topology, the numerical solutions often exhibit a strong mesh-dependence in that, as the mesh is refined, the optimized designs exhibit finer and finer structures. Inspired by its successful use for many applications, starting with linear elasticity problems about two decades ago, we will use a filtering method. More precisely, we use a linear filter as suggested by Bruns and Tortorelli (2001) and define the physical road design and conductivity through respectively, where the M × M filter matrices F R and F C for the road design and conductivity, respectively. Element ij of these filter matrices are given by respectively, where τ R and τ C are the filter radii for the road design and conductivity, respectively, and c R and c C are constant M ×1 normalization vectors selected to ensure that F R 1 M = F C 1 M = 1 M , where 1 M is the M × 1 vector with all ones. To conclude, the problem that we want to solve is That is, the problem is to determine both α and κ to minimize a weighted sum of road construction cost J h R and transportation cost J h T . In the numerical experiments, the method of moving asymptotes (MMA) (Svanberg 1987) solves problem (34). The required gradients are computed by using the expressions provided in the Appendix and the chain rule to take the filtering into account.

Numerical experiments
We consider three test cases denoted TC1, TC2, and TC3, respectively. These test cases are depicted in Fig. 2 and described below. In all test cases, is selected to be the unit square whose lower left and upper right corner are located at (0, 0) and (1, 1), respectively. That is = {x ∈ R 2 | 0 ≤ x 1 ≤ 1, 0 ≤ x 2 ≤ 1}. The boundary portion D ⊂ ∂ on which we stipulate that | D ≡ 0 is selected to be a small portion of ∂ around (0, 0).
For all test cases, q, the total production and consumption of the product, is selected so that the road construction cost J R and the transportation cost J C are of the same order. For both the conductivity and the road design variable, we use the initial guess if |x 1 − 1/2| < 7/16 and |x 2 − 1/2| < 7/16, 0 otherwise. Put simply, the initial guess for both α h and κ h is 1/2 over the whole unit square, except from a border of width 1/16, where it is 0. The reason that we set both these variables to 0 around the edge of the domain is that we expect that there will be no benefit from placing neither roads nor adding any conductivity to these parts. However, since we do not want to influence the design between the supply and demand positions, we let the initial road and conductivity distributions all be equal to 1/2 in a region including these locations and their nearest surrounding for all test cases. For all experiments, we set the parameters controlling the transportation speed, as given by expression (31) as: v min = 1, v max = 5, and p = 3. That is, the transportation speed is five times as large on roads as outside them and it is possible to transport goods even if no roads are placed in . In all experiments, we set the filter radii for the road design and conductivity to τ R = τ C = 1/128.
As a first experiment, we solve optimization (34) for parameter β = 0.5 for all test cases and using a resolution of 256×256 elements. Figure 3 shows the resulting road network designα p (left column) as well as the mass flux ρu (arrows in the right column) for TC1 (top row), TC2 (middle row), and TC3 (bottom row). In all figures showing the road network design and the mass flux in this section, the outline of the computational region is illustrated by a thin line, and the supply and demand positions for the different cases are shown by plus and minus signs, respectively. The optimized road design for TC1 (top row, left image in Fig. 3) is rather intuitive, with essentially three straight lines from the supply positions to the demand position. The optimized road design for TC2 (middle row, left image in Fig. 3) is similar to a binary tree structure with the root located at the single demand position and branch points located between the demand and the supply region. The optimized road design for TC3 (bottom row, left image in Fig. 3) consists of two disjoint parts: the upper main part has a main transportation streak connected to multiple smaller roads that branch further toward the supply as well as the demand region; the lower part is a smaller road that primarily appears to serve the parts of the supply and demand region that are located closest to each other. We remark that an intuitively better design could be obtained by replacing the smaller lower road by a straight one. The obtained design is however a local optimum. Further tests suggest that such non-straight roads are less pronounced at higher resolutions. Figure 4 shows the evolution of the optimization for TC2 with parameter β = 0.5 during the first 30 iterations using a resolution of 256×256 elements. The solid line shows the evolution of the objective function (J h R + J h T )/2, the dotted line shows the evolution of the road construction cost J h R , and the dash-dotted line shows the evolution of the transport cost J h T . The value axis uses a logarithmic scale, so the horizontal dashed lines correspond to the values 0.1, 0.2, . . . , 1.4 and the lines for 0.1 are thicker than the others. All values are normalized with respect to the initial objective function value. The images in Fig. 4 show the road network designα p as well as the mass flux ρu at iterations 1, 5, 10, 15, and 30. In the images for iterations 5, 10, 15, and 30, the road design is on the left and the mass flux is on the right; while, for iteration 1: the road design is on top and the mass flux is at the bottom. The final road network design and mass flux, illustrated in the middle column of Fig. 3 were obtained after 988 iterations and has an objective function value of 0.1926 times the initial objective function value. The objective function value after 30 iterations was 0.2528 times the initial objective function value. After the first 30 iterations, the main features of the road design and mass flux are in place and the optimization progresses slowly in terms of improvement of the objective function: already after 145 iterations, the objective function value is less than 0.2 times the initial objective function value. Fig. 4 Snapshots of the road design and mass flux at iterations 1, 5, 10, 15, and 30 together with the evolution of the transport cost J h T , the road construction cost J h R and the total objective function during the first 30 iterations of the optimization that produced the results in the middle row of Fig. 3 To study the tradeoff between road construction cost and the transportation cost, we solve problem (34) for all three test cases for sequence values of parameter β on two resolutions: 256 × 256 elements and 1024 × 1024 elements discretizing . In this tradeoff study, we solve L = 101 problems and let the parameter values we consider for β be linearly placed between 0 and 1. More precisely, we solve optimization (34) for all β ∈ U L = {0, 1/(L − 1), . . . , 1}.
(36) Figure 5 shows the relative road construction cost and transportation cost for the optimized designs from these experiments. Here, the relative road construction cost and transportations cost for a particular test case and resolution are computed as respectively, whereα * (η) = F R α * (η) andκ * (η) = F C κ * (η), in which κ * (η) and α * (η) solves optimization (23) with β = η. The spacing between the points representing different tradeoffs in Fig. 5 is not uniform with respect to β. To obtain evenly distributed tradeoff points (more precisely, points on the so-called Pareto set), the method proposed by Das and Dennis (1998) may be used. Tracking the Pareto set is a matter deserving careful treatment for an industrial code to minimize the need for user intervention. The focus of the present paper is on the fundamental of the transportation problem, and the authors consider the kinks of the presented Pareto sets to be sufficiently well resolved, albeit using a perhaps unnecessarily high number of points. Selected optimized road designs and mass fluxes, two for each test case, obtained during this tradeoff study on a resolution of 1024×1024 elements are shown in Figs. 6, 7, and 8. The corresponding tradeoff values for road construction cost and transportation cost are marked in Fig. 5. Figure 6 shows results from TC1 obtained with parameter β = 0.18 (bottom row) and β = 0.68 (top row).  Fig. 6; the road's width decreases as β increases until it is comparable with the filter size τ R (cf. the top row in Fig. 3; the results on both resolutions share the same main characteristics). Finally, as β becomes larger, the road construction cost gets increasingly costly compared to the transportation cost and the optimized road design no longer connects the supply and demand, so an increasing share of the transportation needs to be carried out off-road. This is illustrated by the optimized road design obtained for β = 0.68, the bottom left image in Fig. 6, where the high-speed road material only covers a small area close to the demand location. Figure 7 shows results from TC2 obtained with parameter β = 0.16 (bottom row) and β = 0.54 (top row). Just as for TC1, when β is large (β ≥ 0.87 for this experiment) in problem (34), the road construction cost dominates completely and no roads are placed inside , and when β = 0 essentially the full domain is filled with roads. For very small values of β the road design covers the convex envelope of the supply and demand regions, but as β increases, roads start to form. The bottom row shows the road design and mass flux for β = 0.16. In this case, the road design reminds of a root structure with wide roots near the demand positions that branches out to more and finer roots close and over the supply region. As β increases further, the overall root structure remains. However, the individual roots become thinner and the distance between subsequent branch points increases. This procedure continues until the road tree only has few main roots or roads from the demand region to the edge of the supply region. For this test case, this occurs at about β = 0.54, corresponding to the results in the top row of Fig. 7. Increasing β further results in that the roads become shorter from the end closest to the supply region. Figure 8 shows results from TC3 obtained with parameter β = 0.24 (bottom row) and β = 0.73 (top row). Just as for TC1 and TC2, when β is large (β ≥ 0.74 for this experiment) in problem (34), the road construction cost dominates completely and no roads are placed inside , and when β = 0 essentially the full domain is filled with roads. Moreover, just as for TC2, for very small values of β the road design covers the convex envelope of the supply and demand regions but as β increases, roads start to form. For this experiment, the road design for small to intermediate values of parameter β (0.16 ≤ β ≤ 0.23 for this experiment), the road design does not fully cover the convex envelope of the supply and demand regions, and smaller roads start to form near the remote edges of these regions. For intermediate values of β (0.24 ≤ β ≤ 0.73 for Fig. 8 The optimized road design (left column) and the mass flux (right column) for TC3 with β = 0.24 (bottom row) and β = 0.73 (top row) computed on a resolution of 1024×1024 elements this experiment), the road design consists of one connected component that branches out to come close to all parts of the supply and demand regions. As β increases, the width of the individual roads, as well as the number of road segments in each branch, decreases. The bottom and top row in Fig. 8 show the road design and mass flux for β = 0.24 and β = 0.73. In contrast to the two other test cases, the road design connects the supply and demand regions up until the point when β becomes large enough not to motivate the construction of any roads. This is likely due to the fact that in this case both supply and demand are distributed, so that there is no location where one can place a small piece of road that would decrease the transportation cost for most of the goods; for the optimized road designs at high road cost TC1 and TC2 all goods was moved through a small region just in front of the demand region. Perhaps one could find a particular range of values of parameter β so that the optimized road design only consists of a small part of the region between supply and demand, but the numerical experiments presented in this section suggest that if such a parameter region exist it would be very narrow.

Discussion
This paper models and solves continuous transportation problems as material distribution topology optimization problems. The end goal is to find the optimal placement of roads (material) in a region as well as the local transportation direction for a commodity. Here, the road layout determines the local transport speed (we have a high transport speed on the roads and a lower transport speed off the roads). In contrast to typical material distribution topology optimization approaches, the original method, proposed in this paper, uses two design fields that represent the road design and the local transportation direction, respectively.
Given a road design, the suggested potential approach yields problem (19)-a problem that is related to the minimum heat compliance problem, which is a standard test problem for material distribution topology optimization. Both problems are governed by Poisson's equation, but with different source terms and boundary conditions. However, there are two major differences. The first is that the objective in our setting is to minimize a weighted L 1 -norm of the gradient of the solution field of the governing equation (13), while the aim in the minimum heat compliance problem is to minimize a weighted L 2 -norm of the temperature field. The second difference is that in the minimum heat compliance problem, one sets out to place two materials with conductivities κ and 1 so the original problem is to determine for each point which material type it should contain; here, on the contrary, the conductivities may attain any value in the continuous range from κ to 1. (Recall that for the transportation problem, the road design is determined by α-which at each point holds the value 1 or 0 to signify whether the point is occupied by a road or not, respectively-controls the transportation speed, while the conductivities determine the local transportation direction.) The focus on this paper is to simultaneously find the best road design and local transport direction to minimize a total cost including the transport cost and the road construction cost.
By using the strategy proposed in this paper, we have successfully optimized road designs for a few test cases, in which the supply and demand positions are concentrated around given points as well as distributed over given regions. By solving problem (34) for a sequence of values for parameter β, we obtain a sequence of road network designs corresponding to different tradeoffs between the road construction cost and the transportation cost. These solutions and the corresponding relative road construction and transportation cost may provide guidance when making decisions regarding if and where to construct roads to increase transportation efficiency.
An interesting tradeoff problem occurring in many practical situations is when, given a supply and demand function q as well as an existing road network covering the region R ⊂ , one faces the decision if new roads should be constructed and if so where. To solve this problem within the current framework, the only change that needs to be done is to modify the road construction cost, so that the integral in definition (22) is taken over \ R instead of over and modify the set of admissible road designs to ensure that α ≡ 1 in R .
To take for example the presence of regions of water (lakes and rivers) covering an area W ⊂ into account one can either modify the set of feasible road designs to ensure that no roads are placed in W that is, α| W ≡ 0 or modify the road construction cost so that the unit road cost c R becomes spatially dependent and thus the road construction cost becomes Other terrain features, such as slopes and different soil compositions that may affect the suitability and cost of road construction as well as the transportation speed either on-or off-road can be taken into account by using road construction cost (38) and also making the unit transportation cost c T spatially dependent, so that the total transportation cost becomes and also to modify the transportation speed so that is does not only depend on the speed function s but also on the spatial properties. One possible such definition would be to let the off-road transportation speed v min depend on spatial properties, but keeping the on-road transportation speed v max . Such a change can easily be accomplished, as the road design to transportation speed mapping would still remain a point-wise operation in this case.