Add emission limits and penalties

docs/src/format.md

@@ -44,7 +44,7 @@ RELOG accepts as input a JSON file with four sections: `parameters`, `products`,
         "CO2": 0.052,
         "CH4": 0.003
       },
-      "disposal limit (tonne)": 100.0,
+      "disposal limit (tonne)": 100.0
     }
   }
 }
@@ -218,3 +218,31 @@ keys:
   }
 }
 ```
+
+## Emissions
+
+| Key                 | Description                                                                                                                            |
+| :------------------ | :------------------------------------------------------------------------------------------------------------------------------------- |
+| `limit (tonne)`     | Maximum amount of this greenhouse gas allowed to be emitted per year across the entire supply chain. Entry may be `null` if unlimited. |
+| `penalty ($/tonne)` | Penalty cost per tonne of this greenhouse gas emitted.                                                                                 |
+
+#### Example
+
+```json
+{
+  "emissions": {
+    "CO2": {
+      "limit (tonne)": 1000.0,
+      "penalty ($/tonne)": 50.0
+    },
+    "CH4": {
+      "limit (tonne)": null,
+      "penalty ($/tonne)": 1200.0
+    },
+    "N2O": {
+      "limit (tonne)": 10.0,
+      "penalty ($/tonne)": 15000.0
+    }
+  }
+}
+```

docs/src/problem.md

@@ -55,20 +55,22 @@ The mathematical model employed by RELOG is based on three main components:
 | $K^\text{cap}_{p}$          | Capacity of plant $p$, if the plant is open                                                                                                                                                                      | tonne          |
 | $K^\text{disp-limit}_{mt}$  | Maximum amount of material $m$ that can be disposed of (globally) at time $t$                                                                                                                                    | tonne          |
 | $K^\text{disp-limit}_{mut}$ | Maximum amount of material $m$ that can be disposed of at plant/center $u$ at time $t$                                                                                                                           | tonne          |
+| $K^\text{em-limit}_{gt}$    | Maximum amount of greenhouse gas $g$ allowed to be emitted (globally) at time $t$                                                                                                                                | tonne          |
+| $K^\text{em-plant}_{gpt}$   | Amount of greenhouse gas $g$ released by plant $p$ at time $t$ for each tonne of input material processed                                                                                                        | tonne/tonne    |
+| $K^\text{em-tr}_{gmt}$      | Amount of greenhouse gas $g$ released by transporting 1 tonne of material $m$ over one km at time $t$                                                                                                            | tonne/km-tonne |
 | $K^\text{mix}_{pmt}$        | If plant $p$ receives one tonne of input material at time $t$, then $K^\text{mix}_{pmt}$ is the amount of product $m$ in this mix. Must be between zero and one, and the sum of these amounts must equal to one. | tonne          |
+| $K^\text{out-fix}_{cmt}$    | Fixed amount of material $m$ collected at center $m$ at time $t$                                                                                                                                                 | \$/tonne       |
+| $K^\text{out-var-len}_{cm}$ | Length of the $K^\text{out-var}_{c,m,*}$ vector.                                                                                                                                                                 | --             |
+| $K^\text{out-var}_{c,m,i}$  | Factor used to calculate variable amount of material $m$ collected at center $m$. See `eq_z_collected` for more details.                                                                                         | --             |
 | $K^\text{output}_{pmt}$     | Amount of material $m$ produced by plant $p$ at time $t$ for each tonne of input material processed                                                                                                              | tonne          |
-| $K^\text{plant-em}_{gpt}$   | Amount of greenhouse gas $g$ released by plant $p$ at time $t$ for each tonne of input material processed                                                                                                        | tonne/tonne    |
-| $K^\text{tr-em}_{gmt}$      | Amount of greenhouse gas $g$ released by transporting 1 tonne of material $m$ over one km at time $t$                                                                                                            | tonne/km-tonne |
-| $R^\text{tr}_{mt}$          | Cost to send material $m$ at time $t$                                                                                                                                                                            | \$/km-tonne    |
 | $R^\text{collect}_{cmt}$    | Cost of collecting material $m$ at center $c$ at time $t$                                                                                                                                                        | \$/tonne       |
 | $R^\text{disp}_{umt}$       | Cost to dispose of material at plant/center $u$ at time $t$                                                                                                                                                      | \$/tonne       |
+| $R^\text{em}_{gt}$          | Penalty cost per tonne of greenhouse gas $g$ emitted at time $t$                                                                                                                                                 | \$/tonne       |
 | $R^\text{fix}_{ut}$         | Fixed operating cost for plant/center $u$ at time $t$                                                                                                                                                            | \$             |
 | $R^\text{open}_{pt}$        | Cost to open plant $p$ at time $t$                                                                                                                                                                               | \$             |
 | $R^\text{rev}_{ct}$         | Revenue for selling the input product of center $c$ at this center at time $t$                                                                                                                                   | \$/tonne       |
+| $R^\text{tr}_{mt}$          | Cost to send material $m$ at time $t$                                                                                                                                                                            | \$/km-tonne    |
 | $R^\text{var}_{pt}$         | Cost to process one tonne of input material at plant $p$ at time $t$                                                                                                                                             | \$/tonne       |
-| $K^\text{out-fix}_{cmt}$    | Fixed amount of material $m$ collected at center $m$ at time $t$                                                                                                                                                 | \$/tonne       |
-| $K^\text{out-var}_{c,m,i}$  | Factor used to calculate variable amount of material $m$ collected at center $m$. See `eq_z_collected` for more details.                                                                                         | --             |
-| $K^\text{out-var-len}_{cm}$ | Length of the $K^\text{out-var}_{c,m,*}$ vector.                                                                                                                                                                 | --             |
 
 ## Decision variables
 
@@ -80,8 +82,8 @@ The mathematical model employed by RELOG is based on three main components:
 | $z^{\text{disp}}_{umt}$      | `z_disp[u.name, m.name, t]`                  | Amount of product $m$ disposed of at plant/center $u$ at time $t$                                       | tonne  |
 | $z^{\text{input}}_{ut}$      | `z_input[u.name, t]`                         | Total plant/center input at time $t$                                                                    | tonne  |
 | $z^{\text{prod}}_{umt}$      | `z_prod[u.name, m.name, t]`                  | Amount of product $m$ produced by plant/center $u$ at time $t$                                          | tonne  |
-| $z^{\text{tr-em}}_{guvmt}$   | `z_tr_em[g.name, u.name, v.name, m.name, t]` | Amount of greenhouse gas $g$ released at time $t$ due to transportation of material $m$ from $u$ to $v$ | tonne  |
-| $z^{\text{plant-em}}_{gpt}$  | `z_plant_em[g.name, p.name, t]`              | Amount of greenhouse gas $g$ released by plant $p$ at time $t$                                          | tonne  |
+| $z^{\text{em-tr}}_{guvmt}$   | `z_em_tr[g.name, u.name, v.name, m.name, t]` | Amount of greenhouse gas $g$ released at time $t$ due to transportation of material $m$ from $u$ to $v$ | tonne  |
+| $z^{\text{em-plant}}_{gpt}$  | `z_em_plant[g.name, p.name, t]`              | Amount of greenhouse gas $g$ released by plant $p$ at time $t$                                          | tonne  |
 
 ## Objective function
 
@@ -151,6 +153,13 @@ The goals is to minimize a linear objective function with the following terms:
 \sum_{p \in P} \sum_{(u,m) \in E^-(p)} \sum_{t \in T} R^\text{var}_{pt} y_{upmt}
 ```
 
+- Emissions penalty cost, incurred for each tonne of greenhouse gas emitted:
+
+```math
+\sum_{g \in G} \sum_{t \in T} R^\text{em}_{gt} \left(
+  \sum_{p \in P} z^{\text{em-plant}}_{gpt} + \sum_{(u,v,m) \in E} z^{\text{em-tr}}_{guvmt}
+\right)
+```
 ## Constraints
 
 - Definition of plant input (`eq_z_input[p.name, t]`):
@@ -271,20 +280,29 @@ The goals is to minimize a linear objective function with the following terms:
 \end{align*}
 ```
 
-- Computation of transportation emissions
-  (`eq_tr_em[g.name, u.name, v.name, m.name, t`):
+- Computation of transportation emissions (`eq_emission_tr[g.name, u.name, v.name, m.name, t`):
 
 ```math
 \begin{align*}
-& z^{\text{tr-em}}_{guvmt} = K^{\text{dist}}_{uv} K^\text{tr-em}_{gmt} y_{uvmt}
+& z^{\text{em-tr}}_{guvmt} = K^{\text{dist}}_{uv} K^\text{em-tr}_{gmt} y_{uvmt}
 & \forall g \in G, (u, v, m) \in E, t \in T
 \end{align*}
 ```
 
-- Computation of plant emissions (`eq_plant_em[g.name, p.name, t]`):
+- Computation of plant emissions (`eq_emission_plant[g.name, p.name, t]`):
+
 ```math
 \begin{align*}
-& z^{\text{plant-em}}_{gpt} = \sum_{(u,m) \in E^-(p)} K^\text{plant-em}_{gpt} y_{upmt}
+& z^{\text{em-plant}}_{gpt} = \sum_{(u,m) \in E^-(p)} K^\text{em-plant}_{gpt} y_{upmt}
 & \forall g \in G, p \in P, t \in T
+\end{align*}
 ```
 
+- Global emissions limit (`eq_emission_limit[g.name, t]`):
+
+```math
+\begin{align*}
+& \sum_{p \in P} z^{\text{em-plant}}_{gpt} + \sum_{(u,v,m) \in E} z^{\text{em-tr}}_{guvmt} \leq K^\text{em-limit}_{gt}
+& \forall g \in G, t \in T
+\end{align*}
+```

src/instance/parse.jl

@@ -128,6 +128,19 @@ function parse(json)::Instance
         plants_by_name[name] = plant
     end
 
+    # Read emissions
+    emissions = Emissions[]
+    emissions_by_name = OrderedDict{String,Emissions}()
+    if haskey(json, "emissions")
+        for (name, edict) in json["emissions"]
+            limit = timeseries(edict["limit (tonne)"], null_val = Inf)
+            penalty = timeseries(edict["penalty (\$/tonne)"])
+            emission = Emissions(; name, limit, penalty)
+            push!(emissions, emission)
+            emissions_by_name[name] = emission
+        end
+    end
+
     return Instance(;
         time_horizon,
         building_period,
@@ -138,5 +151,7 @@ function parse(json)::Instance
         centers_by_name,
         plants,
         plants_by_name,
+        emissions,
+        emissions_by_name,
     )
 end

src/instance/structs.jl

@@ -54,6 +54,12 @@ Base.@kwdef struct Plant
     initial_capacity::Float64
 end
 
+Base.@kwdef struct Emissions
+    name::String
+    limit::Vector{Float64}
+    penalty::Vector{Float64}
+end
+
 Base.@kwdef struct Instance
     building_period::Vector{Int}
     centers_by_name::OrderedDict{String,Center}
@@ -64,4 +70,6 @@ Base.@kwdef struct Instance
     time_horizon::Int
     plants::Vector{Plant}
     plants_by_name::OrderedDict{String,Plant}
+    emissions_by_name::OrderedDict{String,Emissions}
+    emissions::Vector{Emissions}
 end

src/model/build.jl

@@ -119,15 +119,15 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
     end
 
     # Transportation emissions by greenhouse gas
-    z_tr_em = _init(model, :z_tr_em)
+    z_em_tr = _init(model, :z_em_tr)
     for (p1, p2, m) in E, t in T, g in keys(m.tr_emissions)
-        z_tr_em[g, p1.name, p2.name, m.name, t] = @variable(model, lower_bound = 0)
+        z_em_tr[g, p1.name, p2.name, m.name, t] = @variable(model, lower_bound = 0)
     end
 
     # Plant emissions by greenhouse gas
-    z_plant_em = _init(model, :z_plant_em)
+    z_em_plant = _init(model, :z_em_plant)
     for p in plants, t in T, g in keys(p.emissions)
-        z_plant_em[g, p.name, t] = @variable(model, lower_bound = 0)
+        z_em_plant[g, p.name, t] = @variable(model, lower_bound = 0)
     end
 
 
@@ -192,6 +192,30 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
         )
     end
 
+    # Emissions penalty cost
+    for emission in instance.emissions, t in T
+        # Plant emissions penalty
+        for p in plants
+            if emission.name in keys(p.emissions)
+                add_to_expression!(
+                    obj,
+                    emission.penalty[t],
+                    z_em_plant[emission.name, p.name, t],
+                )
+            end
+        end
+        # Transportation emissions penalty
+        for (p1, p2, m) in E
+            if emission.name in keys(m.tr_emissions)
+                add_to_expression!(
+                    obj,
+                    emission.penalty[t],
+                    z_em_tr[emission.name, p1.name, p2.name, m.name, t],
+                )
+            end
+        end
+    end
+
     @objective(model, Min, obj)
 
     # Constraints
@@ -323,25 +347,41 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
     end
 
     # Transportation emissions
-    eq_tr_em = _init(model, :eq_tr_em)
+    eq_emission_tr = _init(model, :eq_emission_tr)
     for (p1, p2, m) in E, t in T, g in keys(m.tr_emissions)
-        eq_tr_em[g, p1.name, p2.name, m.name, t] = @constraint(
+        eq_emission_tr[g, p1.name, p2.name, m.name, t] = @constraint(
             model,
-            z_tr_em[g, p1.name, p2.name, m.name, t] ==
+            z_em_tr[g, p1.name, p2.name, m.name, t] ==
             distances[p1, p2, m] * m.tr_emissions[g][t] * y[p1.name, p2.name, m.name, t]
         )
     end
 
     # Plant emissions
-    eq_plant_em = _init(model, :eq_plant_em)
+    eq_emission_plant = _init(model, :eq_emission_plant)
     for p in plants, t in T, g in keys(p.emissions)
-        eq_plant_em[g, p.name, t] = @constraint(
+        eq_emission_plant[g, p.name, t] = @constraint(
             model,
-            z_plant_em[g, p.name, t] ==
+            z_em_plant[g, p.name, t] ==
             p.emissions[g][t] * sum(y[src.name, p.name, m.name, t] for (src, m) in E_in[p])
         )
     end
 
+    # Global emissions limit
+    eq_emission_limit = _init(model, :eq_emission_limit)
+    for emission in instance.emissions, t in T
+        isfinite(emission.limit[t]) || continue
+        eq_emission_limit[emission.name, t] = @constraint(
+            model,
+            sum(
+                z_em_plant[emission.name, p.name, t] for
+                p in plants if emission.name in keys(p.emissions)
+            ) + sum(
+                z_em_tr[emission.name, p1.name, p2.name, m.name, t] for
+                (p1, p2, m) in E if emission.name in keys(m.tr_emissions)
+            ) <= emission.limit[t]
+        )
+    end
+
     if variable_names
         _set_names!(model)
     end

test/fixtures/simple.json

@@ -155,5 +155,15 @@
       ],
       "initial capacity (tonne)": 0
     }
+  },
+  "emissions": {
+    "CO2": {
+      "limit (tonne)": [1000.0, 1100.0, 1200.0, 1300.0],
+      "penalty ($/tonne)": [50.0, 55.0, 60.0, 65.0]
+    },
+    "CH4": {
+      "limit (tonne)": null,
+      "penalty ($/tonne)": 1200.0
+    }
   }
 }

test/src/instance/parse_test.jl

@@ -68,6 +68,19 @@ function instance_parse_test_1()
     @test c2.opening_cost == [1000, 1000, 1000, 1000]
     @test c2.fix_operating_cost == [400, 400, 400, 400]
     @test c2.var_operating_cost == [5, 5, 5, 5]
+
+    # Emissions
+    @test length(instance.emissions) == 2
+    co2 = instance.emissions[1]
+    @test co2.name == "CO2"
+    @test co2.limit == [1000.0, 1100.0, 1200.0, 1300.0]
+    @test co2.penalty == [50.0, 55.0, 60.0, 65.0]
+    @test instance.emissions_by_name["CO2"] === co2
+    ch4 = instance.emissions[2]
+    @test ch4.name == "CH4"
+    @test ch4.limit == [Inf, Inf, Inf, Inf]
+    @test ch4.penalty == [1200.0, 1200.0, 1200.0, 1200.0]
+    @test instance.emissions_by_name["CH4"] === ch4
 end
 
 

test/src/model/build_test.jl

@@ -9,8 +9,8 @@ function model_build_test()
     y = model[:y]
     z_disp = model[:z_disp]
     z_input = model[:z_input]
-    z_tr_em = model[:z_tr_em]
-    z_plant_em = model[:z_plant_em]
+    z_em_tr = model[:z_em_tr]
+    z_em_plant = model[:z_em_plant]
     x = model[:x]
     obj = objective_function(model)
     # print(model)
@@ -47,16 +47,16 @@ function model_build_test()
     )
 
     # Variables: Transportation emissions
-    @test haskey(z_tr_em, ("CO2", "L1", "C3", "P4", 1))
-    @test haskey(z_tr_em, ("CH4", "L1", "C3", "P4", 1))
-    @test haskey(z_tr_em, ("CO2", "C2", "L1", "P1", 1))
-    @test haskey(z_tr_em, ("CH4", "C2", "L1", "P1", 1))
+    @test haskey(z_em_tr, ("CO2", "L1", "C3", "P4", 1))
+    @test haskey(z_em_tr, ("CH4", "L1", "C3", "P4", 1))
+    @test haskey(z_em_tr, ("CO2", "C2", "L1", "P1", 1))
+    @test haskey(z_em_tr, ("CH4", "C2", "L1", "P1", 1))
 
     # Variables: Plant emissions
-    @test haskey(z_plant_em, ("CO2", "L1", 1))
-    @test haskey(z_plant_em, ("CO2", "L1", 2))
-    @test haskey(z_plant_em, ("CO2", "L1", 3))
-    @test haskey(z_plant_em, ("CO2", "L1", 4))
+    @test haskey(z_em_plant, ("CO2", "L1", 1))
+    @test haskey(z_em_plant, ("CO2", "L1", 2))
+    @test haskey(z_em_plant, ("CO2", "L1", 3))
+    @test haskey(z_em_plant, ("CO2", "L1", 4))
 
     # Plants: Definition of total plant input
     @test repr(model[:eq_z_input]["L1", 1]) ==
@@ -134,10 +134,23 @@ function model_build_test()
     @test ("P4", 1) ∉ keys(model[:eq_disposal_limit])
 
     # Products: Transportation emissions
-    @test repr(model[:eq_tr_em]["CH4", "L1", "C3", "P4", 1]) ==
-          "eq_tr_em[CH4,L1,C3,P4,1] : -0.333354 y[L1,C3,P4,1] + z_tr_em[CH4,L1,C3,P4,1] = 0"
+    @test repr(model[:eq_emission_tr]["CH4", "L1", "C3", "P4", 1]) ==
+          "eq_emission_tr[CH4,L1,C3,P4,1] : -0.333354 y[L1,C3,P4,1] + z_em_tr[CH4,L1,C3,P4,1] = 0"
 
     # Plants: Plant emissions
-    @test repr(model[:eq_plant_em]["CO2", "L1", 1]) ==
-          "eq_plant_em[CO2,L1,1] : -0.1 y[C2,L1,P1,1] - 0.1 y[C1,L1,P2,1] + z_plant_em[CO2,L1,1] = 0"
+    @test repr(model[:eq_emission_plant]["CO2", "L1", 1]) ==
+          "eq_emission_plant[CO2,L1,1] : -0.1 y[C2,L1,P1,1] - 0.1 y[C1,L1,P2,1] + z_em_plant[CO2,L1,1] = 0"
+
+    # Objective function: Emissions penalty costs
+    @test obj.terms[z_em_plant["CO2", "L1", 1]] == 50.0  # CO2 penalty at time 1
+    @test obj.terms[z_em_plant["CO2", "L1", 2]] == 55.0  # CO2 penalty at time 2
+    @test obj.terms[z_em_plant["CO2", "L1", 3]] == 60.0  # CO2 penalty at time 3
+    @test obj.terms[z_em_plant["CO2", "L1", 4]] == 65.0  # CO2 penalty at time 4
+    @test obj.terms[z_em_tr["CO2", "L1", "C3", "P4", 1]] == 50.0  # CO2 transportation penalty at time 1
+    @test obj.terms[z_em_tr["CH4", "L1", "C3", "P4", 1]] == 1200.0  # CH4 transportation penalty at time 1
+
+    # Global emissions limit constraints
+    @test repr(model[:eq_emission_limit]["CO2", 1]) ==
+          "eq_emission_limit[CO2,1] : z_em_tr[CO2,C2,L1,P1,1] + z_em_tr[CO2,C2,C1,P1,1] + z_em_tr[CO2,C1,L1,P2,1] + z_em_tr[CO2,L1,C3,P4,1] + z_em_plant[CO2,L1,1] ≤ 1000"
+    @test ("CH4", 1) ∉ keys(model[:eq_emission_limit])
 end