Merge pull request #32 from OptimalTransportNetworks/development

Development

Author: SebKrantz

NEWS.md

@@ -1,4 +1,13 @@
+# 0.1.8
+
+* Fixed dual solution with `cross_good_congestion = true`, and set the default `duality = true` in `init_parameters()`. It is highly recommended to keep `beta <= 1` in fixed labor cases to harness the dual solutions, which yield a tremendous speedup.
+
+# 0.1.7
+
+* If `duality = true`, in `init_parameters()` (and `labor_mobility = false` - no dual solutions available for mobile labor cases), a new direct dual implementation of the model is used (with hard-coded sparse hessians passed directly to Ipopt) which is significantly faster than any other means of solving the model (10x speedup). However, with `cross_good_congestion = true`, the dual solution may be inaccurate. This may be fixed in the future.  
+
 # 0.1.6
+
 ## Breaking Changes
 * All spatial parameters, including `Lj`, `Lr`, `Hj`, `hj`, `Zjn`, `omegaj`, `omegar`, and `region` are now stored in the `graph` structure created by `create_graph()`. `create_graph()` therefore only returns the `graph` structure, instead of both the (updated) parameters and the graph. Converesely, `init_parameters()` only contains parameters that are independent of the particular geography defined by `create_graph()`.
 
@@ -7,7 +16,6 @@
 * Better spline options for plotting frictions surface (geography).
 * More faithful translation of `apply_geography()`. 
 
-
 # 0.1.5
 
 * Removed the MATLAB toolbox and corresponding documentation (PDF files) from the repo to decrease size. A new repo was created for the MATLAB toolbox at [SebKrantz/OptimalTransportNetworkToolbox](https://github.com/SebKrantz/OptimalTransportNetworkToolbox). This repo and especially the `docs` folder continue to be very useful for Julia users, but are no longer part of the Julia library. 

Project.toml

@@ -1,7 +1,7 @@
 name = "OptimalTransportNetworks"
 uuid = "e2b46e68-897f-4e4e-ba36-a93c9789fd96"
 authors = ["Sebastian Krantz <sebastian.krantz@graduateinstitute.ch>"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"

src/main/init_parameters.jl

@@ -20,7 +20,7 @@ Returns a `param` dict with the model parameters. These are independent of the g
 - `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=false`: Switch to turn on/off duality whenever available
+- `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
@@ -42,7 +42,7 @@ param = init_parameters(K = 10, labor_mobility = true)
 function init_parameters(; alpha = 0.5, beta = 1, gamma = 1, K = 1, sigma = 5, rho = 2, a = 0.8, N = 1, 
                          labor_mobility = false, cross_good_congestion=false, nu = 2, m = ones(N), 
                          annealing=true, 
-                         verbose = true, duality = false, warm_start = true, 
+                         verbose = true, duality = true, warm_start = true, 
                          kappa_min = 1e-5, min_iter = 20, max_iter = 200, tol = 1e-5, kwargs...)
     param = Dict()
 

src/models/solve_allocation_by_duality.jl

@@ -236,7 +236,7 @@ function hessian_duality(
             end
 
             # Assign result
-            values[ind] = -obj_factor * term + 1e-5 # Somehow need this increment to make it work
+            values[ind] = -obj_factor * term # + 1e-5 # Somehow need this increment to make it work
         end
     end
     return

src/models/solve_allocation_by_duality_cgc.jl

@@ -230,7 +230,7 @@ function hessian_duality_cgc(
             end
 
             # Assign result
-            values[ind] = -obj_factor * term + 1e-6
+            values[ind] = -obj_factor * term # + 1e-6
         end
     end
     return
@@ -280,7 +280,7 @@ function recover_allocation_duality_cgc(x, auxdata)
         LL = max.(Lambda' - Lambda, 0) # P^n_k - P^n_j (non-negative flows)
         Qjk += m[n] * (LL / m[n]) .^ (nu/(nu-1))
     end
-    Qjk = (Qjk .* temp) .^ (nu-1)/(nu*beta)
+    Qjk = (Qjk .* temp) .^ ((nu-1)/(nu*beta))
 
     # Calculate the flows Qjkn
     Qjkn = zeros(graph.J, graph.J, param.N)