Optimal Transport Networks in Spatial Equilibrium

Pure Julia implementation of the model and algorithms presented in:

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

The library is based on the JuMP modeling framework for mathematical optimization in Julia, and roughly follows the MATLAB OptimalTransportNetworkToolbox (v1.0.4b) provided by the authors. Compared to MATLAB, the Julia library presents a simplified interface. Notably, the graph object contains both the graph structure and all data parameterizing the network, whereas the param object only contains (non-spatial) model parameters.

Exported Functions

Problem Setup

init_paramaters() - Create parameters dictionary

create_graph() - Create network graph dictionary

apply_geography() - (Optional) apply geographical features to alter the graph edge weights (network building and traversing costs)

Compute Optimal Network + Refine Solution in Non-Convex Cases

optimal_network() - Compute optimal network given parameters and graph. Can also be used to just solve the model on a given network.

annealing() - Refine solution using simulated annealing in non-convex cases (automatically called in optimal_network() if param[:annealing] == true))

Plot Graph (Incl. Network Solution)

plot_graph() - Plot network graph and optimal infrastructure levels

Helper Functions to Manipulate Graphs

find_node() - Find index of node that is closest to a given pair of coordinates

add_node() - Add new node to graph with given coordinates and connected neighbors

remove_node() - Remove node from graph

Examples

# Convex case
param = init_parameters(K = 10)
graph = create_graph(param)
graph[:Zjn][61] = 10.0
results = optimal_network(param, graph)
plot_graph(graph, results[:Ijk])

# Nonconvex case, disabling automatic annealing
param = init_parameters(K = 10, gamma = 2, annealing = false)
graph = create_graph(param)
graph[:Zjn][61] = 10.0
results = optimal_network(param, graph)

# Run annealing
results_annealing = annealing(param, graph, results[:Ijk])

# Comparison
plot_graph(graph, results[:Ijk])
plot_graph(graph, results_annealing[:Ijk])

init_parameters(; kwargs...) -> Dict

Returns a param dict with the model parameters. These are independent of the graph structure and dimensions.

Keyword Arguments

• alpha::Float64=0.5: Cobb-Douglas coefficient on final good c^alpha * h^(1-alpha)
• beta::Float64=1: Parameter governing congestion in transport cost
• gamma::Float64=1: Elasticity of transport cost relative to infrastructure
• K::Float64=1: Amount of concrete/asphalt (infrastructure budget)
• sigma::Float64=5: Elasticity of substitution across goods (CES)
• rho::Float64=2: Curvature in utility (c^alpha * h^(1-alpha))^(1-rho)/(1-rho)
• a::Float64=0.8: Curvature of the production function L^alpha
• N::Int64=1: Number of goods traded in the economy (used for checks in create_graph())
• labor_mobility::Any=false: Switch for labor mobility (true/false or 'partial')
• cross_good_congestion::Bool=false: Switch for cross-good congestion
• nu::Float64=2: Elasticity of substitution b/w goods in transport costs if cross-good congestion
• m::Vector{Float64}=ones(N): Vector of weights Nx1 in the cross congestion cost function
• annealing::Bool=true: Switch for the use of annealing at the end of iterations (only if gamma > beta)
• verbose::Bool=true: Switch to turn on/off text output (from Ipopt or other optimizers)
• duality::Bool=true: Switch to turn on/off duality whenever available (fixed labor and beta <= 1).
• warm_start::Bool=true: Use the previous solution as a warm start for the next iteration
• kappa_min::Float64=1e-5: Minimum value for road capacities K
• min_iter::Int64=20: Minimum number of iterations
• max_iter::Int64=200: Maximum number of iterations
• tol::Float64=1e-5: Tolerance for convergence of road capacities K

Optional Parameters

• optimizer = Ipopt.Optimizer: Optimizer to be used
• optimizer_attr::Dict: Dict of attributes passed to the optimizer (e.g. Dict(:tol => 1e-5))
• model_attr::Dict: Dict of tuples (length 2) passed to the model (e.g. Dict(:backend => (MOI.AutomaticDifferentiationBackend(), MathOptSymbolicAD.DefaultBackend())) to use Symbolic AD)
• model::Function: For custom models => a function that taks an optimizer and an 'auxdata' structure as created by create_auxdata() as input and returns a fully parameterized JuMP model
• recover_allocation::Function: For custom models => a function that takes a solved JuMP model and 'auxdata' structure as input and returns the allocation variables. In particular, it should return a dict with symbol keys returning at least objects :welfare => scalar welfare measure, :Pjn => prices, :PCj => aggregate consumption, and :Qjkn => flows.

Examples

param = init_parameters(K = 10, labor_mobility = true)

create_graph(param, w = 11, h = 11; type = "map", kwargs...) -> Dict

Initialize the underlying graph, population and productivity parameters.

Arguments

• param::Dict: Dict that contains the model parameters (only needed for checks)
• w::Int64=11: Number of nodes along the width of the underlying graph if type != "custom" (integer)
• h::Int64=11: Number of nodes along the height of the underlying graph if type != "custom" (integer, odd if triangle)

Keyword Arguments

• type::String="map": Either "map", "square", "triangle", or "custom"
• adjacency::BitMatrix: J x J Adjacency matrix (only used for custom network)
• x::Vector{Float64}: x coordinate (longitude) of each node (only used for custom network)
• y::Vector{Float64}: y coordinate (latitude) of each node (only used for custom network)
• omega::Vector{Float64}: Vector of Pareto weights for each node or region in partial mobility case (default ones(J or nregions))
• Zjn::Matrix{Float64}: J x N matrix of producties per node (j = 1:J) and good (n = 1:N) (default ones(J, N))
• Lj::Vector{Float64}: Vector of populations in each node (j = 1:J) (only for fixed labour case)
• Hj::Vector{Float64}: Vector of immobile good in each node (j = 1:J) (e.g. housing, default ones(J))
• Lr::Vector{Float64}: Vector of populations in each region (r = 1:nregions) (only for partial mobility)
• region::Vector{Int64}: Vector indicating region of each location (only for partial mobility)

Examples

graph = create_graph(init_parameters())

apply_geography(graph, geography; kwargs...) -> Dict

Update the network building costs of a graph based on geographical features and remove edges impeded by geographical barriers. Aversion to altitude changes rescales building infrastructure costs delta_i by (see also user manual to MATLAB toolbox):

1 + alpha_up * max(0, Z2-Z1)^beta_up + alpha_down * max(0, Z1-Z2)^beta_down

and similarly for graph traversal costs delta_tau.

Arguments

• graph: Dict or NamedTuple that contains the network graph to which the geographical features will be applied.

• geography: Dict or NamedTuple representing the geographical features, with the following fields:

  • z::Vector{Float64}: (Optional) J x 1 vector containing the z-coordinate (elevation) for each node, or nothing if no elevation data.

  • z_is_friction::Bool: (Optional) logical value indicate that z represents friction rather than elevation. In that case, the measure of building cost is the average friction of the two nodes mean(Z1,Z2) rather than the difference Z2-Z1.

  • obstacles::Matrix{Int64}: (Optional) Nobs x 2 matrix specifying (i, j) pairs of nodes that are connected by obstacles, where Nobs is the number of obstacles, or nothing if no obstacles.

Keyword Arguments

• across_obstacle_delta_i::Float64=Inf: Rescaling parameter for building cost that crosses an obstacle.
• along_obstacle_delta_i::Float64=Inf: Rescaling parameter for building cost that goes along an obstacle.
• across_obstacle_delta_tau::Float64=Inf: Rescaling parameter for transport cost that crosses an obstacle.
• along_obstacle_delta_tau::Float64=Inf: Rescaling parameter for transport cost that goes along an obstacle.
• alpha_up_i::Float64=0: Building cost scale parameter for roads that go up in elevation.
• beta_up_i::Float64=1: Building cost elasticity parameter for roads that go up in elevation.
• alpha_up_tau::Float64=0: Transport cost scale parameter for roads that go up in elevation.
• beta_up_tau::Float64=1: Transport cost elasticity parameter for roads that go up in elevation.
• alpha_down_i::Float64=0: Building cost scale parameter for roads that go down in elevation.
• beta_down_i::Float64=1: Building cost elasticity parameter for roads that go down in elevation.
• alpha_down_tau::Float64=0: Transport cost scale parameter for roads that go down in elevation.
• beta_down_tau::Float64=1: Transport cost elasticity parameter for roads that go down in elevation.

Examples

graph = create_graph(init_parameters())
geography = (z = 10*(rand(graph[:J]) .> 0.95), obstacles = [1 15; 70 72])
updated_graph = apply_geography(graph, geography)

plot_graph(updated_graph, geography = geography, obstacles = true)

optimal_network(param, graph; I0=nothing, Il=nothing, Iu=nothing, 
                verbose=false, return_model=0, solve_allocation = false) -> Dict

Solve for the optimal network by solving the inner problem and the outer problem by iterating over the FOCs.

Arguments

• param: Dict or NamedTuple that contains the model's parameters
• graph: Dict or NamedTuple that contains the underlying graph (created by create_graph() function)
• I0::Matrix{Float64}=nothing: (Optional) J x J matrix providing the initial guess for the iterations
• Il::Matrix{Float64}=nothing: (Optional) J x J matrix providing exogenous lower bound on infrastructure levels
• Iu::Matrix{Float64}=nothing: (Optional) J x J matrix providing exogenous upper bound on infrastructure levels
• verbose::Bool=false: (Optional) tell IPOPT to display results
• return_model::Int=0: (Optional) return the JuMP model and corresponding recover_allocation() function: 1 just returns these before solving the model, while 2 solves the model + optimal network and returns the two alongside the results.
• solve_allocation::Bool=false: (Optional) just solve the model with existing infrastructure I0 and return the results.

Examples

param = init_parameters(K = 10, duality = false)
graph = create_graph(param)
graph[:Zjn][61] = 10.0
results = optimal_network(param, graph)
plot_graph(graph, results[:Ijk])

annealing(param, graph, I0; kwargs...)

Runs the simulated annealing method starting from network I0. Only sensible if param.gamma > param.beta.

Arguments

• param: Dict or NamedTuple that contains the model's parameters
• graph: Dict or NamedTuple that contains the underlying graph
• I0: (optional) provides the initial guess for the iterations

Keyword Arguments

• verbose::Bool=false: (Optional) tell IPOPT to display results
• perturbation_method::String="random rebranching": Method to be used to perturbate the network ("random" is purely random, works horribly; "shake" applies a gaussian blur along a random direction, works alright; "rebranching" (deterministic) and "random rebranching" (default) is the algorithm described in Appendix A.4 in the paper, works nicely)
• preserve_central_symmetry::Bool=false: Only applies to "shake" method
• preserve_vertical_symmetry::Bool=false: Only applies to "shake" method
• preserve_horizontal_symmetry::Bool=false: Only applies to "shake" method
• smooth_network::Bool=true: Whether to smooth the network after each perturbation if perturbation_method is "random"
• smoothing_radius::Float64=0.25: Parameters of the Gaussian blur if perturbation_method is "random" or "shake"
• mu_perturbation::Float64=log(0.3): Parameters of the Gaussian blur if perturbation_method is "shake"
• sigma_perturbation::Float64=0.05: Parameters of the Gaussian blur if perturbation_method is "shake"
• num_random_perturbations::Int64=1: Number of links to be randomly affected ("random" and "random rebranching" only)
• display::Bool: Display the graph in each iteration as we go
• t_start::Float64=100: Initial temperature
• t_end::Float64=1: Final temperature
• t_step::Float64=0.9: Speed of cooling
• num_deepening::Int64=4: Number of FOC iterations between candidate draws
• Iu::Matrix{Float64}=Inf * ones(J, J): J x J matrix of upper bounds on network infrastructure Ijk
• Il::Matrix{Float64}=zeros(J, J): J x J matrix of lower bounds on network infrastructure Ijk
• final_model::JuMPModel: (Optionally) a readily parameterized JuMP model to be used (from optimal_network())
• recover_allocation::Function: The recover_allocation() function corresponding to final_model
• allocation::Dict: The result from recover_allocation() from a previous solution of the model: to skip an initial resolve without perturbations.

Examples

# Nonconvex case, disabling automatic annealing
param = init_parameters(K = 10, gamma = 2, annealing = false, duality = false)
graph = create_graph(param)
graph[:Zjn][61] = 10.0
results = optimal_network(param, graph)

# Run annealing
results_annealing = annealing(param, graph, results[:Ijk])

# Comparison
plot_graph(graph, results[:Ijk])
plot_graph(graph, results_annealing[:Ijk])

plot_graph(graph, edges = nothing; kwargs...) -> Plots.Plot

Plot a graph visualization with various styling options.

Arguments

• graph::Dict: The network graph (created with create_graph())
• edges::Matrix{Float64}=nothing: Matrix of edge weights (J x J)

Keyword Arguments

• grid::Bool=false: Show gridlines
• axis::Tuple=([], false): Axis ticks and labels (see Plots.jl docs, default disable axis)
• margin::Real=0mm: Margin around plot
• aspect_ratio::Symbol=:equal: Plot aspect ratio (set to a real number (h/w) if not :equal)
• height::Int=600: Plot height in pixels (width is proportional to height / aspect_ratio, but also depends on the relative ranges of the x and y coordinates of the graph)
• map::Vector=nothing: Values mapped to graph for background heatmap
• map_color::Symbol=:YlOrBr_4: Colorscale for background heatmap
• mesh::Bool=false: Show mesh lines between nodes
• mesh_color::Symbol=:grey90: Color for mesh lines
• mesh_style::Symbol=:dash: Linestyle for mesh lines
• mesh_transparency::Real=1: Opacity for mesh lines
• edges::Bool=true: Show edges between nodes
• edge_color::Symbol=:blue: Edge color or color gradient
• edge_scaling::Bool=false: Size edges based on raw values
• edge_transparency::Union{Bool,Real}=true: Transparency for edges
• edge_min::Real: Minimum edge value for scaling
• edge_max::Real: Maximum edge value for scaling
• edge_min_thickness::Real=0.1: Minimum thickness for edges
• edge_max_thickness::Real=2: Maximum thickness for edges
• arrows::Bool=false: Show arrowheads on edges
• arrow_scale::Real=1: Scaling factor for arrowheads
• arrow_style::String="long": Style of arrowheads ("long" or "thin")
• nodes::Bool=true: Show nodes
• node_sizes::Vector=ones(J): Sizes for nodes
• node_sizes_scale::Real=75: Overall scaling for node sizes
• node_shades::Vector=nothing: Shades mapped to nodes
• node_color::Symbol=:purple: Node color or color gradient
• node_stroke_width::Real=0: Stroke width for node outlines
• node_stroke_color::Symbol=nothing: Stroke color for node outlines
• geography=nothing: Dict or NamedTuple with geography data, see also apply_geography()
• obstacles::Bool=false: Show obstacles from geography
• obstacle_color::Symbol=:black: Color for obstacles
• obstacle_thickness::Real=3: Thickness for obstacles

Examples

param = init_parameters(K = 10, duality = false)
graph = create_graph(param)
graph[:Zjn][51] = 10.0
results = optimal_network(param, graph)
plot_graph(graph, results[:Ijk])

find_node(graph, x, y) -> Int64

Returns the index of the node closest to the coordinates (x,y) on the graph.

Arguments

• graph::Dict: Dict that contains the underlying graph (created by create_graph())
• x::Float64: x coordinate on the graph (between 1 and w)
• y::Float64: y coordinate on the graph (between 1 and h)

add_node(graph, x, y, neighbors) -> Dict

Add a node in position (x,y) and list of neighbors. The new node is given an index J+1. Returns an updated graph object.

Arguments

• graph::Dict: Dict that contains the underlying graph (created by create_graph())
• x::Float64: x coordinate of the new node (any real number)
• y::Float64: y coordinate of the new node (any real number)
• neighbors::Vector{Int64}: Vector of nodes to which it is connected (1 x n list of node indices between 1 and J, where n is an arbitrary # of neighbors)

Notes

The cost matrices delta_tau and delta_i are parametrized as a function of Euclidean distance between nodes. The new node is given population 1e-6 and productivity equal to the minimum productivity in the graph.

remove_node(graph, i) -> Dict

Removes node i from the graph, returning an updated graph object.

Arguments

• graph::Dict: Dict that contains the underlying graph (created by create_graph())
• i::Int64: index of the mode to be removed (integer between 1 and graph[:J])

dict_to_namedtuple(dict)

Convert a dictionary to a NamedTuple.

If the input is already a NamedTuple, it is returned unchanged. Otherwise, it creates a new NamedTuple from the dictionary's keys and values.

Arguments

• dict: A dictionary or NamedTuple to be converted.

Returns

• A NamedTuple equivalent to the input dictionary.

namedtuple_to_dict(namedtuple)

Convert a NamedTuple to a dictionary.

If the input is already a dictionary, it is returned unchanged. Otherwise, it creates a new dictionary from the NamedTuple's pairs.

Arguments

• namedtuple: A NamedTuple or dictionary to be converted.

Returns

• A dictionary equivalent to the input NamedTuple.

create_map(w, h) -> Dict

Creates a square graph structure with width w and height h (nodes have 8 neighbors in total, along horizontal and vertical dimensions and diagonals)

Arguments

• w: Width of graph (i.e. the number of nodes along horizontal dimension), must be an integer
• h: Height of graph (i.e. the number of nodes along vertical dimension), must be an integer

create_square(w, h) -> Dict

Creates a square graph structure with width w and height h (nodes have 4 neighbors in total, along horizontal and vertical dimensions, NOT diagonals)

Arguments

• w: width of graph (ie. the number of nodes along horizontal dimension), must be an integer
• h: height of graph (ie. the number of nodes along vertical dimension), must be an integer

create_triangle(w, h) -> Dict

Creates a triangular graph structure with width w and height h (each node is the center of a hexagon and each node has 6 neighbors, horizontal and along the two diagonals)

Arguments

• w: Width of graph (i.e. the max number of nodes along horizontal dimension), must be an integer
• h: Height of graph (i.e. the max number of nodes along vertical dimension), must be an odd integer

create_custom(adjacency, x, y) -> Dict

Creates a custom graph structure with given adjacency matrix, x and y vectors of coordinates.

Arguments

• adjacency: Adjacency matrix
• x: Vector of x coordinates of locations
• y: Vector of y coordinates of locations

represent_edges(graph)

Creates a NamedTuple providing detailed representation of the graph edges.

Arguments

• graph: NamedTuple that contains the underlying graph (created by dict_to_namedtuple(create_graph()))

Returns

• A NamedTuple with the following fields:
  • A: J x ndeg matrix where each column represents an edge. The value is 1 if the edge starts at node J, -1 if it ends at node J, and 0 otherwise.
  • Apos: J x ndeg matrix where each column represents an edge. The value is the positive part of the edge flow.
  • Aneg: J x ndeg matrix where each column represents an edge. The value is the negative part of the edge flow.
  • edge_start: J x ndeg matrix where each column represents an edge. The value is the starting node of the edge.
  • edge_end: J x ndeg matrix where each column represents an edge. The value is the ending node of the edge.

create_auxdata(param, graph, edges, I)

Creates the auxdata structure that contains all the auxiliary parameters for estimation

Arguments

• param: NamedTuple that contains the model's parameters (created by dict_to_namedtuple(init_parameters()))
• graph: NamedTuple that contains the underlying graph (created by dict_to_namedtuple(create_graph()))
• edges: NamedTuple that contains the edges of the graph (created by represent_edges())
• I: J x J symmetric matrix of current infrastructure (investments)

Returns

• A NamedTuple with the following fields:
  • param: The input parameter NamedTuple.
  • graph: The input graph NamedTuple.
  • edges: The edges of the graph.
  • kappa: The kappa matrix: I^gamma / delta_tau
  • kappa_ex: The extracted kappa values (ndeg x 1)

get_model(auxdata)

Construct the appropriate JuMP model based on the auxiliary data.

Arguments

• auxdata: Auxiliary data required for constructing the model (created by create_auxdata()).

Returns

• model: The constructed JuMP model.
• recover_allocation: A function to recover the allocation from the model solution.