Model: Objective function and plant constraints

Project.toml

@@ -4,6 +4,7 @@ authors = ["Alinson S. Xavier <git@axavier.org>"]
 version = "0.1.0"
 
 [deps]
+Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"

src/RELOG.jl

@@ -3,6 +3,7 @@ module RELOG
 include("instance/structs.jl")
 include("instance/parse.jl")
 include("model/jumpext.jl")
+include("model/dist.jl")
 include("model/build.jl")
 
 end # module RELOG

src/model/build.jl

@@ -9,22 +9,35 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
 
     # Transportation edges
     # -------------------------------------------------------------------------
+
+    # Connectivity
     E = []
+    E_in = Dict(src => [] for src in plants ∪ centers)
+    E_out = Dict(src => [] for src in plants ∪ centers)
+
+    function push_edge!(src, dst, m)
+        @show src.name, dst.name, m.name
+        push!(E, (src, dst, m))
+        push!(E_out[src], (dst, m))
+        push!(E_in[dst], (src, m))
+    end
+
     for m in products
         for p1 in plants
-            m ∉ keys(p1.output) || continue
+            m ∈ keys(p1.output) || continue
 
             # Plant to plant
             for p2 in plants
                 p1 != p2 || continue
                 m ∉ keys(p2.input_mix) || continue
-                push!(E, (p1, p2, m))
+                push_edge!(p1, p2, m)
             end
 
             # Plant to center
             for c in centers
+                @show m.name, p1.name, c.name, m == c.input
                 m == c.input || continue
-                push!(E, (p1, c, m))
+                push_edge!(p1, c, m)
             end
         end
 
@@ -34,23 +47,33 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
             # Center to plant
             for p in plants
                 m ∈ keys(p.input_mix) || continue
-                push!(E, (c1, p, m))
+                push_edge!(c1, p, m)
             end
 
             # Center to center
             for c2 in centers
                 m == c2.input || continue
-                push!(E, (c1, c2, m))
+                push_edge!(c1, c2, m)
             end
         end
     end
 
+    # Distances
+    distances = Dict()
+    for (p1, p2, m) in E
+        d = _calculate_distance(p1.latitude, p1.longitude, p2.latitude, p2.longitude)
+        distances[p1, p2, m] = d
+        @show p1.name, p2.name, m.name, d
+    end
 
     # Decision variables
     # -------------------------------------------------------------------------
 
     # Plant p is operational at time t
     x = _init(model, :x)
+    for p in plants
+        x[p.name, 0] = p.initial_capacity > 0 ? 1 : 0
+    end
     for p in plants, t in T
         x[p.name, t] = @variable(model, binary = true)
     end
@@ -58,35 +81,155 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
     # Amount of product m sent from center/plant u to center/plant v at time T
     y = _init(model, :y)
     for (p1, p2, m) in E, t in T
-        y[p1.name, p2.name, m.name, t] = @variable(model, lower_bound=0)
+        y[p1.name, p2.name, m.name, t] = @variable(model, lower_bound = 0)
     end
 
     # Amount of product m produced by plant/center at time T
     z_prod = _init(model, :z_prod)
     for p in plants, m in keys(p.output), t in T
-        z_prod[p.name, m.name, t] = @variable(model, lower_bound=0)
+        z_prod[p.name, m.name, t] = @variable(model, lower_bound = 0)
     end
     for c in centers, m in c.outputs, t in T
-        z_prod[c.name, m.name, t] = @variable(model, lower_bound=0)
+        z_prod[c.name, m.name, t] = @variable(model, lower_bound = 0)
     end
 
     # Amount of product m disposed at plant/center p at time T
     z_disp = _init(model, :z_disp)
     for p in plants, m in keys(p.output), t in T
-        z_disp[p.name, m.name, t] = @variable(model, lower_bound=0)
+        z_disp[p.name, m.name, t] = @variable(model, lower_bound = 0)
     end
     for c in centers, m in c.outputs, t in T
-        z_disp[c.name, m.name, t] = @variable(model, lower_bound=0)
+        z_disp[c.name, m.name, t] = @variable(model, lower_bound = 0)
     end
-    
+
+    # Total plant input
+    z_input = _init(model, :z_input)
+    for p in plants, t in T
+        z_input[p.name, t] = @variable(model, lower_bound = 0)
+    end
+
 
     # Objective function
     # -------------------------------------------------------------------------
+    obj = AffExpr()
 
+    # Transportation cost
+    for (p1, p2, m) in E, t in T
+        obj += distances[p1, p2, m] * y[p1.name, p2.name, m.name, t]
+    end
+
+    # Center: Revenue
+    for c in centers, (p, m) in E_in[c], t in T
+        obj += c.revenue[t] * y[p.name, c.name, m.name, t]
+    end
+
+    # Center: Collection cost
+    for c in centers, (p, m) in E_out[c], t in T
+        obj += c.collection_cost[m][t] * y[c.name, p.name, m.name, t]
+    end
+
+    # Center: Disposal cost
+    for c in centers, m in c.outputs, t in T
+        obj += c.disposal_cost[m][t] * z_disp[c.name, m.name, t]
+    end
+
+    # Center: Operating cost
+    for c in centers, t in T
+        obj += c.operating_cost[t]
+    end
+
+    # Plants: Disposal cost
+    for p in plants, m in keys(p.output), t in T
+        obj += p.disposal_cost[m][t] * z_disp[p.name, m.name, t]
+    end
+
+    # Plants: Opening cost
+    for p in plants, t in T
+        obj += p.capacities[1].opening_cost[t] * (x[p.name, t] - x[p.name, t-1])
+    end
+
+    # Plants: Fixed operating cost
+    for p in plants, t in T
+        obj += p.capacities[1].fix_operating_cost[t] * x[p.name, t]
+    end
+
+    # Plants: Variable operating cost
+    for p in plants, (src, m) in E_in[p], t in T
+        obj += p.capacities[1].var_operating_cost[t] * y[src.name, p.name, m.name, t]
+    end
+
+    @objective(model, Min, obj)
 
     # Constraints
     # -------------------------------------------------------------------------
 
+    # Plants: Definition of total plant input
+    eq_z_input = _init(model, :eq_z_input)
+    for p in plants, t in T
+        eq_z_input[p.name, t] = @constraint(
+            model,
+            z_input[p.name, t] ==
+            sum(y[src.name, p.name, m.name, t] for (src, m) in E_in[p])
+        )
+    end
+
+    # Plants: Must meet input mix
+    eq_input_mix = _init(model, :eq_input_mix)
+    for p in plants, m in keys(p.input_mix), t in T
+        eq_input_mix[p.name, m.name, t] = @constraint(
+            model,
+            sum(y[src.name, p.name, m.name, t] for (src, m2) in E_in[p] if m == m2) ==
+            z_input[p.name, t] * p.input_mix[m][t]
+        )
+    end
+
+    # Plants: Calculate amount produced
+    eq_z_prod = _init(model, :eq_z_prod)
+    for p in plants, m in keys(p.output), t in T
+        eq_z_prod[p.name, m.name, t] = @constraint(
+            model,
+            z_prod[p.name, m.name, t] == z_input[p.name, t] * p.output[m][t]
+        )
+    end
+
+    # Plants: Produced material must be sent or disposed
+    eq_balance = _init(model, :eq_balance)
+    for p in plants, m in keys(p.output), t in T
+        eq_balance[p.name, m.name, t] = @constraint(
+            model,
+            z_prod[p.name, m.name, t] ==
+            sum(y[p.name, dst.name, m.name, t] for (dst, m2) in E_out[p] if m == m2) +
+            z_disp[p.name, m.name, t]
+        )
+    end
+
+    # Plants: Capacity limit
+    eq_capacity = _init(model, :eq_capacity)
+    for p in plants, t in T
+        eq_capacity[p.name, t] = @constraint(
+            model,
+            z_input[p.name, t] <= p.capacities[1].size * x[p.name, t]
+        )
+    end
+
+    # Plants: Disposal limit
+    eq_disposal_limit = _init(model, :eq_disposal_limit)
+    for p in plants, m in keys(p.output), t in T
+        isfinite(p.disposal_limit[m][t]) || continue
+        eq_disposal_limit[p.name, m.name, t] = @constraint(
+            model,
+            z_disp[p.name, m.name, t] <= p.disposal_limit[m][t]
+        )
+    end
+
+    # Plants: Plant remains open
+    eq_keep_open = _init(model, :eq_keep_open)
+    for p in plants, t in T
+        eq_keep_open[p.name, t] = @constraint(
+            model,
+            x[p.name, t] >= x[p.name, t-1]
+        )
+    end
 
     if variable_names
         _set_names!(model)

src/model/dist.jl

@@ -0,0 +1,11 @@
+# RELOG: Reverse Logistics Optimization
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+using Geodesy
+
+function _calculate_distance(source_lat, source_lon, dest_lat, dest_lon)::Float64
+    x = LLA(source_lat, source_lon, 0.0)
+    y = LLA(dest_lat, dest_lon, 0.0)
+    return round(euclidean_distance(x, y) / 1000.0, digits = 3)
+end

src/model/jumpext.jl

@@ -14,7 +14,7 @@ function fix(x::Float64, v::Float64; force)
     return abs(x - v) < 1e-6 || error("Value mismatch: $x != $v")
 end
 
-function set_name(x::Float64, n::String)
+function set_name(x::Number, n::String)
     # nop
 end
 
@@ -44,4 +44,4 @@ function _set_names!(dict::Dict)
             set_name(dict[name][idx], "$name[$idx_str]")
         end
     end
-end
+end

test/fixtures/simple.json

@@ -68,32 +68,32 @@
       }
     },
     "C2": {
-      "latitude (deg)": 41.881,
+      "latitude (deg)": 42.881,
       "longitude (deg)": -87.623,
       "input": null,
-      "outputs": ["P4"],
+      "outputs": ["P1"],
       "variable output (tonne/tonne)": {
-        "P4": 0
+        "P1": 0
       },
       "fixed output (tonne)": {
-        "P4": [50, 60, 70, 80]
+        "P1": [50, 60, 70, 80]
       },
       "revenue ($/tonne)": null,
       "collection cost ($/tonne)": {
-        "P4": 0.25
+        "P1": 0.25
       },
       "operating cost ($)": [150.0, 150.0, 150.0, 150.0],
       "disposal limit (tonne)": {
-        "P4": null
+        "P1": null
       },
       "disposal cost ($/tonne)": {
-        "P4": 0
+        "P1": 0
       }
     },
     "C3": {
-      "latitude (deg)": 41.881,
+      "latitude (deg)": 43.881,
       "longitude (deg)": -87.623,
-      "input": "P1",
+      "input": "P4",
       "outputs": [],
       "variable output (tonne/tonne)": {},
       "constant output (tonne)": {},
@@ -106,7 +106,7 @@
   },
   "plants": {
     "L1": {
-      "latitude (deg)": 41.881,
+      "latitude (deg)": 44.881,
       "longitude (deg)": -87.623,
       "input mix (%)": {
         "P1": 95.3,
@@ -138,7 +138,7 @@
       "capacities": [
         {
           "size (tonne)": 100,
-          "opening cost ($)": 500,
+          "opening cost ($)": [300, 400, 450, 475],
           "fixed operating cost ($)": 300,
           "variable operating cost ($/tonne)": 5.0
         },
@@ -149,7 +149,7 @@
           "variable operating cost ($/tonne)": 5.0
         }
       ],
-      "initial capacity (tonne)": 150
+      "initial capacity (tonne)": 0
     }
   }
 }

test/src/RELOGT.jl

@@ -6,6 +6,7 @@ using JuliaFormatter
 
 include("instance/parse_test.jl")
 include("model/build_test.jl")
+include("model/dist_test.jl")
 
 basedir = dirname(@__FILE__)
 
@@ -18,6 +19,7 @@ function runtests()
         instance_parse_test_1()
         instance_parse_test_2()
         model_build_test()
+        model_dist_test()
     end
 end
 

test/src/instance/parse_test.jl

@@ -45,7 +45,7 @@ function instance_parse_test_1()
     # Plants
     @test length(instance.plants) == 1
     l1 = instance.plants[1]
-    @test l1.latitude == 41.881
+    @test l1.latitude == 44.881
     @test l1.longitude == -87.623
     @test l1.input_mix ==
           Dict(p1 => [0.953, 0.953, 0.953, 0.953], p2 => [0.047, 0.047, 0.047, 0.047])
@@ -56,11 +56,11 @@ function instance_parse_test_1()
     @test l1.disposal_cost == Dict(p3 => [0, 0, 0, 0], p4 => [0.86, 0.86, 0.86, 0.86])
     @test l1.disposal_limit ==
           Dict(p3 => [Inf, Inf, Inf, Inf], p4 => [1000.0, 1000.0, 1000.0, 1000.0])
-    @test l1.initial_capacity == 150
+    @test l1.initial_capacity == 0
     @test length(l1.capacities) == 2
     c1 = l1.capacities[1]
     @test c1.size == 100
-    @test c1.opening_cost == [500, 500, 500, 500]
+    @test c1.opening_cost == [300, 400, 450, 475]
     @test c1.fix_operating_cost == [300, 300, 300, 300]
     @test c1.var_operating_cost == [5, 5, 5, 5]
     c2 = l1.capacities[2]
@@ -74,4 +74,4 @@ end
 function instance_parse_test_2()
     # Should not crash
     RELOG.parsefile(fixture("boat_example.json"))
-end
+end

test/src/model/build_test.jl

@@ -5,6 +5,79 @@ using JuMP
 
 function model_build_test()
     instance = RELOG.parsefile(fixture("simple.json"))
-    model = RELOG.build_model(instance, optimizer=HiGHS.Optimizer, variable_names=true)
+    model = RELOG.build_model(instance, optimizer = HiGHS.Optimizer, variable_names = true)
-    print(model)
+    y = model[:y]
-end
+    z_disp = model[:z_disp]
+    z_input = model[:z_input]
+    x = model[:x]
+    obj = objective_function(model)
+    # print(model)
+
+    @test obj.terms[y["L1", "C3", "P4", 1]] == (
+        111.118 + # transportation
+        12.0 # revenue
+    )
+    @test obj.terms[y["C1", "L1", "P2", 4]] == (
+        333.262 +  # transportation
+        0.25 + # center collection cost
+        5.0 # plant operating cost
+    )
+    @test obj.terms[z_disp["C1", "P2", 1]] == 0.23
+    @test obj.constant == (
+        150 * 4 * 3 # center operating cost
+    )
+    @test obj.terms[z_disp["L1", "P4", 2]] == 0.86
+    @test obj.terms[x["L1", 1]] == (
+        -100.0 + # opening cost
+        300 # fixed operating cost
+    )
+    @test obj.terms[x["L1", 2]] == (
+        -50.0 + # opening cost
+        300 # fixed operating cost
+    )
+    @test obj.terms[x["L1", 3]] == (
+        -25.0 + # opening cost
+        300 # fixed operating cost
+    )
+    @test obj.terms[x["L1", 4]] == (
+        475.0 + # opening cost
+        300 # fixed operating cost
+    )
+
+    # Plants: Definition of total plant input
+    @test repr(model[:eq_z_input]["L1", 1]) ==
+          "eq_z_input[L1,1] : -y[C2,L1,P1,1] - y[C1,L1,P2,1] + z_input[L1,1] = 0"
+
+    # Plants: Must meet input mix
+    @test repr(model[:eq_input_mix]["L1", "P1", 1]) ==
+        "eq_input_mix[L1,P1,1] : y[C2,L1,P1,1] - 0.953 z_input[L1,1] = 0"
+    @test repr(model[:eq_input_mix]["L1", "P2", 1]) ==
+        "eq_input_mix[L1,P2,1] : y[C1,L1,P2,1] - 0.047 z_input[L1,1] = 0"
+
+    # Plants: Calculate amount produced
+    @test repr(model[:eq_z_prod]["L1", "P3", 1]) ==
+        "eq_z_prod[L1,P3,1] : z_prod[L1,P3,1] - 0.25 z_input[L1,1] = 0"
+    @test repr(model[:eq_z_prod]["L1", "P4", 1]) ==
+        "eq_z_prod[L1,P4,1] : z_prod[L1,P4,1] - 0.12 z_input[L1,1] = 0"
+
+    # Plants: Produced material must be sent or disposed
+    @test repr(model[:eq_balance]["L1", "P3", 1]) ==
+        "eq_balance[L1,P3,1] : z_prod[L1,P3,1] - z_disp[L1,P3,1] = 0"
+    @test repr(model[:eq_balance]["L1", "P4", 1]) ==
+        "eq_balance[L1,P4,1] : -y[L1,C3,P4,1] + z_prod[L1,P4,1] - z_disp[L1,P4,1] = 0"
+
+    # Plants: Capacity limit
+    @test repr(model[:eq_capacity]["L1", 1]) ==
+        "eq_capacity[L1,1] : -100 x[L1,1] + z_input[L1,1] ≤ 0"
+
+    # Plants: Disposal limit
+    @test repr(model[:eq_disposal_limit]["L1", "P4", 1]) ==
+        "eq_disposal_limit[L1,P4,1] : z_disp[L1,P4,1] ≤ 1000"
+    @test ("L1", "P3", 1) ∉ keys(model[:eq_disposal_limit])
+
+    # Plants: Plant remains open
+    @test repr(model[:eq_keep_open]["L1", 4]) ==
+        "eq_keep_open[L1,4] : -x[L1,3] + x[L1,4] ≥ 0"
+    @test repr(model[:eq_keep_open]["L1", 1]) == "eq_keep_open[L1,1] : x[L1,1] ≥ 0"
+
+end

test/src/model/dist_test.jl

@@ -0,0 +1,10 @@
+# RELOG: Reverse Logistics Optimization
+# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+# Released under the modified BSD license. See COPYING.md for more details.
+
+using RELOG
+
+function model_dist_test()
+    # Euclidean distance between Chicago and Indianapolis
+    @test RELOG._calculate_distance(41.866, -87.656, 39.764, -86.148) == 265.818
+end