Merge pull request #33 from OptimalTransportNetworks/development

Development

Author: SebKrantz

examples/example01.jl

@@ -8,16 +8,14 @@ using OptimalTransportNetworks
 # Set parameters: try with labor mobility: true/false, convex: beta>=gamma,
 # nonconvex: gamma>beta, cross good congestion: true/false
 
-param = init_parameters(labor_mobility = true, K = 10, gamma = 1, beta = 1, verbose = true,
-                        N = 1, cross_good_congestion = false, nu = 1, duality = false)
-
-# tol is the tolerance in distance b/w iterations for road capacity
-# kappa=I^gamma/delta_i, default is 1e-7 but we relax it a bit here
+param = init_parameters(labor_mobility = true, K = 10, gamma = 1, beta = 1, N = 1, 
+                        cross_good_congestion = true)
 
 # ------------
 # Init network
 
-graph = create_graph(param, 11, 11, type = "map") # create a map network of 11x11 nodes located in [0,10]x[0,10]
+# create a map network of 11x11 nodes located in [0,10]x[0,10]
+graph = create_graph(param, 11, 11, type = "map") 
 # note: by default, productivity and population are equalized everywhere
 
 # Customize graph
@@ -33,7 +31,6 @@ graph[:Zjn][Ni, :] .= 1 # central node more productive
 
 
 # Plot the network with the optimal transportation plan
-
 plot_graph(graph, results[:Ijk], node_sizes = results[:Cj])
 
 # The size/shade of the nodes reflects the total consumption at each node

examples/example04.jl

@@ -9,7 +9,7 @@ import Random
 # Init and Solve network
 
 param = init_parameters(K = 100, labor_mobility = false,
-                        N = 1, gamma = 1, beta = 1, duality = false)
+                        N = 1, gamma = 1, beta = 1, duality = true)
 w = 13; h = 13
 graph = create_graph(param, w, h, type = "map")
 

src/main/annealing.jl

@@ -37,7 +37,7 @@ Runs the simulated annealing method starting from network `I0`. Only sensible if
 # Examples
 ```julia
 # Nonconvex case, disabling automatic annealing
-param = init_parameters(K = 10, gamma = 2, annealing = false)
+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)

src/main/optimal_network.jl

@@ -21,7 +21,7 @@ Solve for the optimal network by solving the inner problem and the outer problem
 
 # Examples
 ```julia
-param = init_parameters(K = 10)
+param = init_parameters(K = 10, duality = false)
 graph = create_graph(param)
 graph[:Zjn][61] = 10.0
 results = optimal_network(param, graph)

src/main/plot_graph.jl

@@ -49,7 +49,7 @@ Plot a graph visualization with various styling options.
 
 # Examples
 ```julia
-param = init_parameters(K = 10)
+param = init_parameters(K = 10, duality = false)
 graph = create_graph(param)
 graph[:Zjn][51] = 10.0
 results = optimal_network(param, graph)

src/models/solve_allocation_by_duality_cgc.jl

@@ -146,8 +146,10 @@ function hessian_duality_cgc(
 
         # Constant term: first part of Bprime
         cons = m1dbeta * n1dnum1 * numbetam1
-        # Constant in T'
-        cons2 = numbetam1 / (numbetam1+nu*beta)
+        ## Constant in T'
+        # cons2 = numbetam1 / (numbetam1+nu*beta)
+        # New: constant in T' -> More accurate!
+        cons3 = m1dbeta * n1dnum1 * (numbetam1+nu*beta) / (1+beta)
 
         # Precompute elements
         res = recover_allocation_duality_cgc(x, auxdata)
@@ -201,10 +203,13 @@ function hessian_duality_cgc(
                     if jd == j
                         KPprimeAB1 = m1dbeta * Pjn[j, nd]^(-sigma) * PCj[j]^(sigma-1)
                         term += Qjkn[j, k, n] * (KPprimeAB1 - KPAprimeB1 + KPABprime1) # Derivative of Qjkn
-                        term += T * (((1+beta) * KPprimeAB1 + cons2 * KPABprime1) * G + Gprime) # T'G + TG'
+                        # term += T * (((1+beta) * KPprimeAB1 + cons2 * KPABprime1) * G + Gprime) # T'G + TG'
+                        term += T * ((1+beta) * KPprimeAB1 * G + Gprime) # T'G (first part) + TG'
+                        term += cons3 * Qjkn[j, k, nd] / PCj[j] * G # Second part of T'G
                     elseif jd == k
                         term += Qjkn[j, k, n] * (KPAprimeB1 - KPABprime1) # Derivative of Qjkn
-                        term -= T * cons2 * KPABprime1 * G # T'G: second part [B'(k) has opposite sign]
+                        # term -= T * cons2 * KPABprime1 * G # T'G: second part [B'(k) has opposite sign]
+                        term -= cons3 * Qjkn[j, k, nd] / PCj[j] * G # Second part of T'G 
                     end
                 end
                 if Qjkn[k, j, n] > 0 # Flows in the direction of j