Update models to reflect new param and graph division.

Author: SebKrantz

src/models/model_fixed.jl

@@ -10,7 +10,7 @@ function model_fixed(optimizer, auxdata)
     Apos = auxdata.edges.Apos
     Aneg = auxdata.edges.Aneg
     psigma = (param.sigma - 1) / param.sigma
-    Lj = param.Lj
+    Lj = graph.Lj
 
     # Model
     model = Model(optimizer)
@@ -27,13 +27,13 @@ function model_fixed(optimizer, auxdata)
 
     # Defining Utility Funcion: from Cjn + parameters (by operator overloading)
     @expression(model, Cj, sum(Cjn .^ psigma, dims=2) .^ (1 / psigma))
-    @expression(model, cj, ifelse.(param.Lj .== 0, 0.0, Cj ./ param.Lj))
-    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (param.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
-    @expression(model, U, sum(param.omegaj .* param.Lj .* uj))
+    @expression(model, cj, ifelse.(graph.Lj .== 0, 0.0, Cj ./ graph.Lj))
+    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (graph.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
+    @expression(model, U, sum(graph.omegaj .* graph.Lj .* uj))
     @objective(model, Max, U)
 
     # Define Yjn (production) as expression
-    @expression(model, Yjn[j=1:graph.J, n=1:param.N], param.Zjn[j, n] * Ljn[j, n]^param.a)
+    @expression(model, Yjn[j=1:graph.J, n=1:param.N], graph.Zjn[j, n] * Ljn[j, n]^param.a)
     # Balanced flow constraints
     @constraint(model, Pjn[j in 1:param.J, n in 1:param.N],
         Cjn[j, n] + sum(A[j, i] * Qin[i, n] for i in 1:graph.ndeg) -
@@ -62,9 +62,9 @@ function recover_allocation_fixed(model, auxdata)
     results[:Cjn] = value.(model_dict[:Cjn])
     results[:Cj] = dropdims(value.(model_dict[:Cj]), dims = 2)
     results[:Ljn] = value.(model_dict[:Ljn])
-    results[:Lj] = param.Lj
+    results[:Lj] = graph.Lj
     results[:cj] = dropdims(value.(model_dict[:cj]), dims = 2)
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = dropdims(value.(model_dict[:uj]), dims = 2)
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_fixed_armington.jl

@@ -10,9 +10,9 @@ function model_fixed_armington(optimizer, auxdata)
     Apos = auxdata.edges.Apos
     Aneg = auxdata.edges.Aneg
     psigma = (param.sigma - 1) / param.sigma
-    Lj = param.Lj
+    Lj = graph.Lj
     # Production: Fixed because just one good is produced by each location
-    Yjn = param.Zjn .* Lj .^ param.a
+    Yjn = graph.Zjn .* Lj .^ param.a
 
     # Model
     model = Model(optimizer)
@@ -29,8 +29,8 @@ function model_fixed_armington(optimizer, auxdata)
     # Defining Utility Funcion: from Cjn + parameters (by operator overloading)
     @expression(model, Cj, sum(Cjn .^ psigma, dims=2) .^ (1 / psigma))
     @expression(model, cj, ifelse.(Lj .== 0, 0.0, Cj ./ Lj))
-    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (param.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
-    @expression(model, U, sum(param.omegaj .* Lj .* uj))
+    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (graph.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
+    @expression(model, U, sum(graph.omegaj .* Lj .* uj))
     @objective(model, Max, U)
 
     # Balanced flow constraints
@@ -51,18 +51,18 @@ function recover_allocation_fixed_armington(model, auxdata)
     model_dict = model.obj_dict
     results = Dict()
     # Production: Fixed because just one good is produced by each location
-    Yj = dropdims(sum(param.Zjn, dims = 2) .* param.Lj .^ param.a, dims = 2)
-    WY = param.Zjn .> 0
+    Yj = dropdims(sum(graph.Zjn, dims = 2) .* graph.Lj .^ param.a, dims = 2)
+    WY = graph.Zjn .> 0
 
     results[:welfare] = value(model_dict[:U])
     results[:Yjn] = WY .* Yj
     results[:Yj] = Yj
     results[:Cjn] = value.(model_dict[:Cjn])
     results[:Cj] = dropdims(value.(model_dict[:Cj]), dims = 2)
-    results[:Ljn] = WY .* param.Lj
-    results[:Lj] = param.Lj
+    results[:Ljn] = WY .* graph.Lj
+    results[:Lj] = graph.Lj
     results[:cj] = dropdims(value.(model_dict[:cj]), dims = 2)
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = dropdims(value.(model_dict[:uj]), dims = 2)
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_fixed_cgc.jl

@@ -9,7 +9,7 @@ function model_fixed_cgc(optimizer, auxdata)
     A = auxdata.edges.A
     Apos = auxdata.edges.Apos
     Aneg = auxdata.edges.Aneg
-    Lj = param.Lj
+    Lj = graph.Lj
     m = param.m # 1:param.N: Vector of weights on each goods flow for aggregate congestion term
     psigma = (param.sigma - 1) / param.sigma
     beta_nu = (param.beta + 1) / param.nu
@@ -30,19 +30,19 @@ function model_fixed_cgc(optimizer, auxdata)
     set_parameter_value.(kappa_ex, kappa_ex_init)
 
     # Defining Utility Funcion: from cj + parameters (by operator overloading)
-    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (param.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
-    @expression(model, U, sum(param.omegaj .* param.Lj .* uj))      # Overall Utility
+    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (graph.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
+    @expression(model, U, sum(graph.omegaj .* graph.Lj .* uj))      # Overall Utility
     @objective(model, Max, U)
 
     # Create the matrix B_direct (resp. B_indirect) of transport cost along the direction of the edge (resp. in edge opposite direction)
     B_direct = @expression(model, ((Qin_direct .^ param.nu) * m) .^ beta_nu ./ kappa_ex)
     B_indirect = @expression(model, ((Qin_indirect .^ param.nu) * m) .^ beta_nu ./ kappa_ex)
     # Final good constraints
     @expression(model, Dj, sum(Djn .^ psigma, dims=2) .^ (1 / psigma))
-    @constraint(model, cj .* param.Lj + Apos * B_direct + Aneg * B_indirect - Dj .<= -1e-8)
+    @constraint(model, cj .* graph.Lj + Apos * B_direct + Aneg * B_indirect - Dj .<= -1e-8)
 
     # Balanced flow constraints
-    @expression(model, Yjn, param.Zjn .* (Ljn .^ param.a))
+    @expression(model, Yjn, graph.Zjn .* (Ljn .^ param.a))
     @constraint(model, Pjn, Djn + A * Qin_direct - A * Qin_indirect - Yjn .<= -1e-8)
 
     # Local labor availability constraints ( sum Ljn <= Lj )
@@ -61,12 +61,12 @@ function recover_allocation_fixed_cgc(model, auxdata)
     results[:Yjn] = value.(model_dict[:Yjn])
     results[:Yj] = dropdims(sum(results[:Yjn], dims=2), dims = 2) 
     results[:Ljn] = value.(model_dict[:Ljn])
-    results[:Lj] = param.Lj
+    results[:Lj] = graph.Lj
     results[:Djn] = value.(model_dict[:Djn]) # Consumption per good pre-transport cost
     results[:Dj] = dropdims(value.(model_dict[:Dj]), dims = 2)
     results[:cj] = value.(model_dict[:cj])
-    results[:Cj] = results[:cj] .* param.Lj
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:Cj] = results[:cj] .* graph.Lj
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = value.(model_dict[:uj])
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_fixed_cgc_armington.jl

@@ -9,7 +9,7 @@ function model_fixed_cgc_armington(optimizer, auxdata)
     A = auxdata.edges.A
     Apos = auxdata.edges.Apos
     Aneg = auxdata.edges.Aneg
-    Lj = param.Lj
+    Lj = graph.Lj
     m = param.m # 1:param.N: Vector of weights on each goods flow for aggregate congestion term
     psigma = (param.sigma - 1) / param.sigma
     beta_nu = (param.beta + 1) / param.nu
@@ -29,8 +29,8 @@ function model_fixed_cgc_armington(optimizer, auxdata)
     set_parameter_value.(kappa_ex, kappa_ex_init)
 
     # Defining Utility Funcion: from cj + parameters (by operator overloading)
-    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (param.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
-    @expression(model, U, sum(param.omegaj .* Lj .* uj))      # Overall Utility
+    @expression(model, uj, ((cj / param.alpha) .^ param.alpha .* (graph.hj / (1-param.alpha)) .^ (1-param.alpha)) .^ (1-param.rho) / (1-param.rho))
+    @expression(model, U, sum(graph.omegaj .* Lj .* uj))      # Overall Utility
     @objective(model, Max, U)
 
     # Create the matrix B_direct (resp. B_indirect) of transport cost along the direction of the edge (resp. in edge opposite direction)
@@ -41,7 +41,7 @@ function model_fixed_cgc_armington(optimizer, auxdata)
     @constraint(model, cj .* Lj + Apos * B_direct + Aneg * B_indirect - Dj .<= -1e-8)
 
     # Balanced flow constraints
-    Yjn = param.Zjn .* Lj .^ param.a
+    Yjn = graph.Zjn .* Lj .^ param.a
     @constraint(model, Pjn, Djn + A * Qin_direct - A * Qin_indirect - Yjn .<= -1e-8)
 
     return model
@@ -53,19 +53,19 @@ function recover_allocation_fixed_cgc_armington(model, auxdata)
     model_dict = model.obj_dict
     results = Dict()
     # Production: Fixed because just one good is produced by each location
-    Yj = dropdims(sum(param.Zjn, dims = 2) .* param.Lj .^ param.a, dims = 2)
-    WY = param.Zjn .> 0
+    Yj = dropdims(sum(graph.Zjn, dims = 2) .* graph.Lj .^ param.a, dims = 2)
+    WY = graph.Zjn .> 0
 
     results[:welfare] = value(model_dict[:U])
     results[:Yjn] = WY .* Yj
     results[:Yj] = Yj
-    results[:Ljn] = WY .* param.Lj
-    results[:Lj] = param.Lj
+    results[:Ljn] = WY .* graph.Lj
+    results[:Lj] = graph.Lj
     results[:Djn] = value.(model_dict[:Djn]) # Consumption per good pre-transport cost
     results[:Dj] = dropdims(value.(model_dict[:Dj]), dims = 2)
     results[:cj] = value.(model_dict[:cj])
-    results[:Cj] = results[:cj] .* param.Lj
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:Cj] = results[:cj] .* graph.Lj
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = value.(model_dict[:uj])
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_fixed_duality.jl

@@ -10,14 +10,14 @@ function model_fixed_duality(optimizer, auxdata)
     Aneg = auxdata.edges.Aneg
     es = auxdata.edges.edge_start
     ee = auxdata.edges.edge_end
-    omegaj = param.omegaj
-    Lj = param.Lj
+    omegaj = graph.omegaj
+    Lj = graph.Lj
     alpha = param.alpha
     sigma = param.sigma
     rho = param.rho
-    hj = param.hj
+    hj = graph.hj
     hj1malpha = (hj / (1-alpha)) .^ (1-alpha)
-    Zjn = param.Zjn
+    Zjn = graph.Zjn
     beta = param.beta
     a = param.a
 
@@ -92,12 +92,12 @@ function recover_allocation_fixed_duality(model, auxdata)
     results[:welfare] = value(model_dict[:U])
     results[:Yjn] = value.(model_dict[:Yjn])
     results[:Yj] = dropdims(sum(results[:Yjn], dims=2), dims = 2)
-    results[:Cjn] = value.(model_dict[:cjn]) .* repeat(param.Lj, 1, param.N)
+    results[:Cjn] = value.(model_dict[:cjn]) .* repeat(graph.Lj, 1, param.N)
     results[:cj] = value.(model_dict[:cj])
-    results[:Cj] = results[:cj] .* param.Lj
+    results[:Cj] = results[:cj] .* graph.Lj
     results[:Ljn] = value.(model_dict[:Ljn])
-    results[:Lj] = param.Lj
-    results[:hj] = param.hj
+    results[:Lj] = graph.Lj
+    results[:hj] = graph.hj
     results[:uj] = value.(model_dict[:uj]) # param.u.(results[:cj], results[:hj])
     # Prices
     results[:Pjn] = value.(model_dict[:Pjn])

src/models/model_mobility.jl

@@ -10,7 +10,7 @@ function model_mobility(optimizer, auxdata)
     Apos = auxdata.edges.Apos
     Aneg = auxdata.edges.Aneg
     psigma = (param.sigma - 1) / param.sigma
-    Hj = param.Hj
+    Hj = graph.Hj
 
     # Model
     model = Model(optimizer)
@@ -37,8 +37,8 @@ function model_mobility(optimizer, auxdata)
     end
 
     # balanced flow constraints
-    # Yjn = @expression(model, param.Zjn .* Ljn .^ param.a) # Same thing
-    @expression(model, Yjn[j=1:graph.J, n=1:param.N], param.Zjn[j, n] * Ljn[j, n]^param.a)
+    # Yjn = @expression(model, graph.Zjn .* Ljn .^ param.a) # Same thing
+    @expression(model, Yjn[j=1:graph.J, n=1:param.N], graph.Zjn[j, n] * Ljn[j, n]^param.a)
     @constraint(model, Pjn[j in 1:param.J, n in 1:param.N],
         Cjn[j, n] + sum(A[j, i] * Qin[i, n] for i in 1:graph.ndeg) -
         Yjn[j, n] + sum(
@@ -71,7 +71,7 @@ function recover_allocation_mobility(model, auxdata)
     results[:Ljn] = value.(model_dict[:Ljn])
     results[:Lj] = value.(model_dict[:Lj])
     results[:cj] = ifelse.(results[:Lj] .== 0, 0.0, results[:Cj] ./ results[:Lj])
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = param.u.(results[:cj], results[:hj])
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_mobility_armington.jl

@@ -10,7 +10,7 @@ function model_mobility_armington(optimizer, auxdata)
     Apos = auxdata.edges.Apos
     Aneg = auxdata.edges.Aneg
     psigma = (param.sigma - 1) / param.sigma
-    Hj = param.Hj
+    Hj = graph.Hj
 
     # Model
     model = Model(optimizer)
@@ -36,7 +36,7 @@ function model_mobility_armington(optimizer, auxdata)
     end
 
     # balanced flow constraints
-     @expression(model, Yjn[j=1:graph.J, n=1:param.N], param.Zjn[j, n] * Lj[j]^param.a)
+     @expression(model, Yjn[j=1:graph.J, n=1:param.N], graph.Zjn[j, n] * Lj[j]^param.a)
     @constraint(model, Pjn[j in 1:param.J, n in 1:param.N],
         Cjn[j, n] + sum(A[j, i] * Qin[i, n] for i in 1:graph.ndeg) -
         Yjn[j, n] + sum(
@@ -64,9 +64,9 @@ function recover_allocation_mobility_armington(model, auxdata)
     results[:Cjn] = value.(model_dict[:Cjn])
     results[:Cj] = dropdims(sum(results[:Cjn] .^ ((param.sigma-1)/param.sigma), dims=2), dims = 2) .^ (param.sigma/(param.sigma-1))
     results[:Lj] = value.(model_dict[:Lj])
-    results[:Ljn] = (param.Zjn .> 0) .* results[:Lj]
+    results[:Ljn] = (graph.Zjn .> 0) .* results[:Lj]
     results[:cj] = ifelse.(results[:Lj] .== 0, 0.0, results[:Cj] ./ results[:Lj])
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = param.u.(results[:cj], results[:hj])
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_mobility_cgc.jl

@@ -34,7 +34,7 @@ function model_mobility_cgc(optimizer, auxdata)
     @objective(model, Max, u)
 
     # Utility constraint (Lj*u <= ... )
-    @constraint(model, Lj .* u - (cj .* Lj / param.alpha) .^ param.alpha .* (param.Hj / (1 - param.alpha)) .^ (1 - param.alpha) .<= -1e-8)
+    @constraint(model, Lj .* u - (cj .* Lj / param.alpha) .^ param.alpha .* (graph.Hj / (1 - param.alpha)) .^ (1 - param.alpha) .<= -1e-8)
 
     # Create the matrix B_direct (resp. B_indirect) of transport cost along the direction of the edge (resp. in edge opposite direction)
     B_direct = @expression(model, ((Qin_direct .^ param.nu) * m) .^ beta_nu ./ kappa_ex)
@@ -44,7 +44,7 @@ function model_mobility_cgc(optimizer, auxdata)
     @constraint(model, cj .* Lj + Apos * B_direct + Aneg * B_indirect - Dj .<= -1e-8)
 
     # Balanced flow constraints
-    @expression(model, Yjn, param.Zjn .* (Ljn .^ param.a))
+    @expression(model, Yjn, graph.Zjn .* (Ljn .^ param.a))
     @constraint(model, Pjn, Djn + A * Qin_direct - A * Qin_indirect - Yjn .<= -1e-8)
 
     # Labor resource constraint
@@ -71,7 +71,7 @@ function recover_allocation_mobility_cgc(model, auxdata)
     results[:Dj] = dropdims(value.(model_dict[:Dj]), dims = 2)
     results[:cj] = value.(model_dict[:cj])
     results[:Cj] = results[:cj] .* results[:Lj]
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = param.u.(results[:cj], results[:hj])
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])

src/models/model_mobility_cgc_armington.jl

@@ -33,7 +33,7 @@ function model_mobility_cgc_armington(optimizer, auxdata)
     @objective(model, Max, u)
 
     # Utility constraint (Lj*u <= ... )
-    @constraint(model, Lj .* u - (cj .* Lj / param.alpha) .^ param.alpha .* (param.Hj / (1 - param.alpha)) .^ (1 - param.alpha) .<= -1e-8)
+    @constraint(model, Lj .* u - (cj .* Lj / param.alpha) .^ param.alpha .* (graph.Hj / (1 - param.alpha)) .^ (1 - param.alpha) .<= -1e-8)
 
     # Create the matrix B_direct (resp. B_indirect) of transport cost along the direction of the edge (resp. in edge opposite direction)
     B_direct = @expression(model, ((Qin_direct .^ param.nu) * m) .^ beta_nu ./ kappa_ex)
@@ -43,7 +43,7 @@ function model_mobility_cgc_armington(optimizer, auxdata)
     @constraint(model, cj .* Lj + Apos * B_direct + Aneg * B_indirect - Dj .<= -1e-8)
 
     # Balanced flow constraints
-    @expression(model, Yjn, param.Zjn .* (Lj .^ param.a))
+    @expression(model, Yjn, graph.Zjn .* (Lj .^ param.a))
     @constraint(model, Pjn, Djn + A * Qin_direct - A * Qin_indirect - Yjn .<= -1e-8)
 
     # Labor resource constraint
@@ -62,12 +62,12 @@ function recover_allocation_mobility_cgc_armington(model, auxdata)
     results[:Yjn] = value.(model_dict[:Yjn])
     results[:Yj] = dropdims(sum(results[:Yjn], dims=2), dims = 2) 
     results[:Lj] = value.(model_dict[:Lj])
-    results[:Ljn] = (param.Zjn .> 0) .* results[:Lj]
+    results[:Ljn] = (graph.Zjn .> 0) .* results[:Lj]
     results[:Djn] = value.(model_dict[:Djn]) # Consumption per good pre-transport cost
     results[:Dj] = dropdims(value.(model_dict[:Dj]), dims = 2)
     results[:cj] = value.(model_dict[:cj])
     results[:Cj] = results[:cj] .* results[:Lj]
-    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, param.Hj ./ results[:Lj])
+    results[:hj] = ifelse.(results[:Lj] .== 0, 0.0, graph.Hj ./ results[:Lj])
     results[:uj] = param.u.(results[:cj], results[:hj])
     # Prices
     results[:Pjn] = shadow_price.(model_dict[:Pjn])