Merge pull request #44 from OptimalTransportNetworks/development

Development

Author: SebKrantz

NEWS.md

@@ -1,3 +1,7 @@
+# 0.1.9
+
+* The exact dual solution for the Hessian with `cross_good_congestion = true` does not have good numerical properties in some cases. Therefore, by default now an approximate solution is used which works better for most problems. Users can set `duality = 2` to use the exact solution in the CGC case. 
+
 # 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.

Project.toml

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

README.md

@@ -31,7 +31,7 @@ This plot shows the optimal network after 200 iterations, keeping population fix
 
 ## 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. One exception is dual solutions which are implemented directly using Ipopt and hard-coded sparse hessians - they are super fast. By default `duality = true`, and if labor is fixed and `beta <= 1`, users will experience very fast dual solves. I also expect non-dual solutions to speed up when [support for detecting nonlinear subexpressions](https://github.com/jump-dev/JuMP.jl/issues/3738) will be added to JuMP.  
+* The Julia implementation for the most part 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. One exception is dual solutions which are implemented directly using Ipopt and hard-coded sparse Hessians - they are super fast. By default `duality = true`, and if labor is fixed and `beta <= 1`, users will experience very fast dual solves. I also expect non-dual solutions to speed up when [support for detecting nonlinear subexpressions](https://github.com/jump-dev/JuMP.jl/issues/3738) will be added to JuMP.  
 
 * Symbolic autodifferentiation via [MathOptSymbolicAD.jl](https://github.com/lanl-ansi/MathOptSymbolicAD.jl) can also provide significant performance improvements for non-dual cases. The symbolic backend can be activated using:
 

docs/src/index.md

@@ -26,7 +26,7 @@ This plot shows the optimal network after 200 iterations, keeping population fix
 
 ## 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. One exception is dual solutions which are implemented directly using Ipopt and hard-coded sparse hessians - they are super fast. By default `duality = true`, and if labor is fixed and `beta <= 1`, users will experience very fast dual solves. I also expect non-dual solutions to speed up when [support for detecting nonlinear subexpressions](https://github.com/jump-dev/JuMP.jl/issues/3738) will be added to JuMP.  
+* The Julia implementation for the most part 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. One exception is dual solutions which are implemented directly using Ipopt and hard-coded sparse Hessians - they are super fast. By default `duality = true`, and if labor is fixed and `beta <= 1`, users will experience very fast dual solves. I also expect non-dual solutions to speed up when [support for detecting nonlinear subexpressions](https://github.com/jump-dev/JuMP.jl/issues/3738) will be added to JuMP.  
 
 * Symbolic autodifferentiation via [MathOptSymbolicAD.jl](https://github.com/lanl-ansi/MathOptSymbolicAD.jl) can also provide significant performance improvements for non-dual cases. The symbolic backend can be activated using:
 

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=true`: Switch to turn on/off duality whenever available (fixed labor and beta <= 1)
+- `duality::Bool=true`: Switch to turn on/off duality whenever available (fixed labor and beta <= 1). Note that if `cross_good_congestion == true`, setting `duality = 2` uses an exact algorithm to compute the hessian, which has however not shown good numerical properties in some cases. 
 - `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

src/models/solve_allocation_by_duality_cgc.jl

@@ -129,7 +129,7 @@ function hessian_duality_cgc(
         graph = auxdata.graph
         nodes = graph.nodes
         kappa = auxdata.kappa
-        inexact_algo = param.duality == 2
+        inexact_algo = param.duality == true # param.duality == 2 for exact algorithm
         beta = param.beta
         m1dbeta = -1 / beta
         sigma = param.sigma