OptimalTransportNetworks.jl

Optimal Transport Networks in Spatial Equilibrium - in Julia

Modern Julia (JuMP) translation of the MATLAB OptimalTransportNetworkToolbox (v1.0.4b) implementing the quantitative spatial economic model of:

Fajgelbaum, P. D., & Schaal, E. (2020). Optimal transport networks in spatial equilibrium. Econometrica, 88(4), 1411-1452.

The model/software uses duality principles to optimize over the space of networks, nesting an optimal flows problem and a neoclasical general-equilibrium trade model into a global network design problem to derive the optimal (welfare maximizing) transport network (extension) from any primitive set of economic fundamantals [population per location, productivity per location for each of N traded goods, endowment of a non-traded good, and (optionally) a pre-existing transport network].

For more information about the model see this folder and the MATLAB User Guide.

The model is the first of its kind and a pathbreaking contribution towards the welfare maximizing planning of transport infrastructure. Its creation has been funded by the European Union through an ERC Research Grant. The author of this Julia library has no personal connections to the authors, but has used their Matlab library for research purposes and belives that it deserves an accessible open-source implementation. Community efforts to further improve the code are welcome. In particular, there is a probabilistic extenstion to solving the model using MCMC methods which may be more suitable for large networks, implemented in:

Kreindler, G., Gaduh, A., Graff, T., Hanna, R., & Olken, B. A. (2023). Optimal Public Transportation Networks: Evidence from the World's Largest Bus Rapid Transit System in Jakarta (No. w31369). National Bureau of Economic Research.

Example

The code for this example is in example04.jl. See the examples folder for more examples.

This plot shows the endowments on a map-graph: circle size is population, circle colour is productivity (the central node is more productive), the black lines indicate geographic barriers, and the background is shaded according to the cost of network building (elevation), indicating a mountain in the upper right corner.

This plot shows the optimal network after 200 iterations, keeping population fixed and not allowing for cross-good congestion. The size of nodes indicates consumption in each node.

Performance Notes

The Julia implementation does not provide hard-coded Gradients, Jacobians, and Hessians as the MATLAB implementation does for some model cases, but relies solely on JuMP's automatic differentiation. This has proven ineffective for dual solutions to the model where the objective is quite complex. Thus, at present, duality does not help to speed up computations in Julia, and accordingly the default is duality = false. I expect this to change in when support for detecting nonlinear subexpressions will be added to JuMP.
Symbolic autodifferentiation via MathOptSymbolicAD.jl can provide significant performance improvements. The symbolic backend can be activated using:

import MathOptInterface as MOI
import MathOptSymbolicAD

param[:model_attr] = Dict(:backend => (MOI.AutomaticDifferentiationBackend(), 
                                        MathOptSymbolicAD.DefaultBackend())) 
                                # Or:  MathOptSymbolicAD.ThreadedBackend()

It is recommended to use Coin-HSL linear solvers for Ipopt to speed up computations. In my opinion the simplest way to use them is to get a (free for academics) license and download the binaries here, extract them somewhere, and then set the hsllib (place here the path to where you extracted libhsl.dylib, it may also be called libcoinhsl.dylib, in which case you may have to rename it to libhsl.dylib) and linear_solver options as follows:

param[:optimizer_attr] = Dict(:hsllib => "/usr/local/lib/libhsl.dylib", # Adjust path
                                :linear_solver => "ma57") # Use ma57, ma86 or ma97

The Ipopt.jl README suggests to use the larger LibHSL package for which there exists a Julia module and proceed similarly. In addition, users may try an optimized BLAS and see if it yields significant performance gains (and let me know if it does).