|
|
# RELOG: Reverse Logistics Optimization
|
|
|
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
|
|
|
# Released under the modified BSD license. See COPYING.md for more details.
|
|
|
|
|
|
using CSV
|
|
|
using DataFrames
|
|
|
using JSON
|
|
|
using JuMP
|
|
|
using OrderedCollections
|
|
|
using Random
|
|
|
using TimerOutputs
|
|
|
using Printf
|
|
|
using Gurobi
|
|
|
|
|
|
include("jumpext.jl")
|
|
|
include("dist.jl")
|
|
|
|
|
|
dict = OrderedDict
|
|
|
Time = Int
|
|
|
Random.seed!(42)
|
|
|
|
|
|
# Structs
|
|
|
# =========================================================================
|
|
|
|
|
|
Base.@kwdef mutable struct Component
|
|
|
name::String
|
|
|
end
|
|
|
|
|
|
Base.@kwdef mutable struct Product
|
|
|
name::String
|
|
|
comp::Vector{Component}
|
|
|
end
|
|
|
|
|
|
Base.@kwdef mutable struct Center
|
|
|
name::String
|
|
|
latitude::Float64
|
|
|
longitude::Float64
|
|
|
prod_out::Vector{Product}
|
|
|
end
|
|
|
|
|
|
Base.@kwdef mutable struct Plant
|
|
|
name::String
|
|
|
#county::String
|
|
|
latitude::Float64
|
|
|
longitude::Float64
|
|
|
prod_out::Vector{Product}
|
|
|
prod_in::Product
|
|
|
end
|
|
|
|
|
|
Base.@kwdef mutable struct Emission
|
|
|
name::String
|
|
|
end
|
|
|
|
|
|
Base.show(io::IO, p::Union{Component,Product,Center,Plant,Emission}) = print(io, p.name)
|
|
|
|
|
|
Base.@kwdef mutable struct Instance
|
|
|
T::UnitRange{Int}
|
|
|
centers::Vector{Center}
|
|
|
plants::Vector{Plant}
|
|
|
products::Vector{Product}
|
|
|
emissions::Vector{Emission}
|
|
|
alpha_mix::dict{Tuple{Plant,Product,Component,Component},Float64}
|
|
|
alpha_plant_emission::dict{Tuple{Plant,Emission,Time},Float64}
|
|
|
alpha_tr_emission::dict{Tuple{Product,Emission,Time},Float64}
|
|
|
c_acq::dict{Tuple{Center,Product,Time},Float64}
|
|
|
#c_interv::dict{Tuple{Center,Product,Time},Float64}
|
|
|
c_center_disp::dict{Tuple{Center,Product,Time},Float64}
|
|
|
c_emission::dict{Tuple{Emission,Time},Float64}
|
|
|
c_fix::dict{Tuple{Plant,Time},Float64}
|
|
|
c_open::dict{Tuple{Plant,Time},Float64}
|
|
|
c_plant_disp::dict{Tuple{Plant,Product,Time},Float64}
|
|
|
c_store::dict{Tuple{Center,Product,Time},Float64}
|
|
|
c_tr::dict{Tuple{Product,Time},Float64}
|
|
|
c_var::dict{Tuple{Plant,Time},Float64}
|
|
|
m_cap::dict{Plant,Float64}
|
|
|
m_center_disp::dict{Tuple{Product,Time},Float64}
|
|
|
m_dist::dict{Tuple{Union{Center,Plant},Plant},Float64}
|
|
|
m_emission::dict{Tuple{Emission,Time},Float64}
|
|
|
m_init::dict{Tuple{Center,Product,Component,Time},Float64}
|
|
|
m_plant_disp::dict{Tuple{Plant,Product,Time},Float64}
|
|
|
m_store::dict{Tuple{Center,Product,Time},Float64}
|
|
|
selected_edges::Union{Nothing,Set} = nothing
|
|
|
max_cap::dict{Plant,Float64}
|
|
|
max_OpC::dict{Plant,Float64}
|
|
|
max_OpenC::dict{Plant,Float64}
|
|
|
min_cap::dict{Plant,Float64}
|
|
|
min_OpC::dict{Plant,Float64}
|
|
|
min_OpenC::dict{Plant,Float64}
|
|
|
unit_open::dict{Plant,Float64}
|
|
|
unit_OpC::dict{Plant,Float64}
|
|
|
end
|
|
|
|
|
|
|
|
|
# Generate
|
|
|
# ==============================================================================
|
|
|
|
|
|
function generate_data()
|
|
|
# Time window
|
|
|
T = 1:2
|
|
|
|
|
|
# Cities
|
|
|
cities_a = dict(
|
|
|
"Chicago" => [41.881832, -87.623177],
|
|
|
"New York City" => [40.712776, -74.005974],
|
|
|
"Los Angeles" => [34.052235, -118.243683],
|
|
|
"Houston" => [29.760427, -95.369804],
|
|
|
"Phoenix" => [33.448376, -112.074036],
|
|
|
"Philadelphia" => [39.952583, -75.165222],
|
|
|
"San Antonio" => [29.424122, -98.493629],
|
|
|
"San Diego" => [32.715736, -117.161087],
|
|
|
"Dallas" => [32.776664, -96.796988],
|
|
|
"San Jose" => [37.338208, -121.886329],
|
|
|
)
|
|
|
|
|
|
cities_b = dict(
|
|
|
"Chicago" => [41.881832, -87.623177],
|
|
|
"Phoenix" => [33.448376, -112.074036],
|
|
|
"Dallas" => [32.776664, -96.796988],
|
|
|
)
|
|
|
|
|
|
# Components
|
|
|
film = Component("Film")
|
|
|
paper = Component("Paper")
|
|
|
cardboard = Component("Cardboard")
|
|
|
|
|
|
# Products
|
|
|
waste = Product(name="Waste", comp=[film, paper, cardboard])
|
|
|
film_bale = Product(name="Film bale", comp=[film, paper, cardboard])
|
|
|
cardboard_bale = Product(name="Cardboard bale", comp=[paper, cardboard])
|
|
|
cardboard_sheets = Product(name="Cardboard sheets", comp=[cardboard])
|
|
|
products = [waste, film_bale, cardboard_bale, cardboard_sheets]
|
|
|
|
|
|
# Centers
|
|
|
centers = [
|
|
|
Center(
|
|
|
name="Collection ($city_name)",
|
|
|
latitude=city_lat,
|
|
|
longitude=city_lon,
|
|
|
prod_out=[waste],
|
|
|
) for (city_name, (city_lat, city_lon)) in cities_a
|
|
|
]
|
|
|
|
|
|
# Plants
|
|
|
plants_a = [
|
|
|
Plant(
|
|
|
name="MRF ($city_name)",
|
|
|
latitude=city_lat,
|
|
|
longitude=city_lon,
|
|
|
prod_in=waste,
|
|
|
prod_out=[film_bale, cardboard_bale],
|
|
|
) for (city_name, (city_lat, city_lon)) in cities_b
|
|
|
]
|
|
|
plants_b = [
|
|
|
Plant(
|
|
|
name="Paper Mill ($city_name)",
|
|
|
latitude=city_lat,
|
|
|
longitude=city_lon,
|
|
|
prod_in=cardboard_bale,
|
|
|
prod_out=[cardboard_sheets],
|
|
|
) for (city_name, (city_lat, city_lon)) in cities_b
|
|
|
]
|
|
|
plants = [plants_a; plants_b]
|
|
|
|
|
|
# Emissions
|
|
|
emissions = [Emission("CO2")]
|
|
|
|
|
|
alpha_mix = dict(
|
|
|
(p, r, cin, cout) => 0.0 for p in plants for cin in p.prod_in.comp for
|
|
|
r in p.prod_out for cout in r.comp
|
|
|
)
|
|
|
for p in plants_a
|
|
|
alpha_mix[p, film_bale, film, film] = 0.98
|
|
|
alpha_mix[p, film_bale, paper, paper] = 0.02
|
|
|
alpha_mix[p, film_bale, cardboard, cardboard] = 0.02
|
|
|
alpha_mix[p, cardboard_bale, paper, paper] = 0.02
|
|
|
alpha_mix[p, cardboard_bale, cardboard, cardboard] = 0.75
|
|
|
end
|
|
|
for p in plants_b
|
|
|
alpha_mix[p, cardboard_sheets, cardboard, cardboard] = 0.95
|
|
|
end
|
|
|
alpha_plant_emission = dict((p, s, t) => 0.01 for p in plants, s in emissions, t in T)
|
|
|
alpha_tr_emission = dict((r, s, t) => 0.01 for r in products, s in emissions, t in T)
|
|
|
c_acq = dict((q, r, t) => 1.0 for q in centers for r in q.prod_out, t in T)
|
|
|
c_center_disp = dict((q, r, t) => 0 for q in centers for r in q.prod_out for t in T)
|
|
|
c_emission = dict((s, t) => 0.01 for s in emissions, t in T)
|
|
|
c_fix = dict((p, t) => 1_000.0 for p in plants, t in T)
|
|
|
c_open = dict((p, t) => 10_000.0 for p in plants, t in T)
|
|
|
c_plant_disp = dict(
|
|
|
(p, r, t) => (r == cardboard_sheets ? -100.0 : -10.0) for p in plants for
|
|
|
r in p.prod_out, t in T
|
|
|
)
|
|
|
c_store = dict((q, r, t) => 1.0 for q in centers for r in q.prod_out, t in T)
|
|
|
c_tr = dict((r, t) => 0.05 for r in products, t in T)
|
|
|
c_var = dict((p, t) => 1.0 for p in plants, t in T)
|
|
|
m_cap = dict(p => 50000.0 for p in plants)
|
|
|
m_center_disp = dict((r, t) => 0.0 for r in products, t in T)
|
|
|
metric = KnnDrivingDistance()
|
|
|
m_dist = dict(
|
|
|
(p, q) => _calculate_distance(
|
|
|
p.latitude,
|
|
|
p.longitude,
|
|
|
q.latitude,
|
|
|
q.longitude,
|
|
|
metric,
|
|
|
) for p in [plants; centers], q in plants
|
|
|
)
|
|
|
m_emission = dict((s, t) => 1_000_000 for s in emissions, t in T)
|
|
|
m_init = dict()
|
|
|
for q in centers, r in q.prod_out
|
|
|
ratio = dict(c => rand(1:10) for c in r.comp)
|
|
|
total = dict(t => rand(1:1000) for t in T)
|
|
|
for c in r.comp, t in T
|
|
|
m_init[q, r, c, t] = ratio[c] * total[t]
|
|
|
end
|
|
|
end
|
|
|
m_plant_disp =
|
|
|
dict((p, r, t) => 1_000_000 for p in plants for r in p.prod_out for t in T)
|
|
|
m_store = dict((q, r, t) => 1_000 for q in centers for r in q.prod_out, t in T)
|
|
|
|
|
|
return Instance(;
|
|
|
T,
|
|
|
centers,
|
|
|
plants,
|
|
|
products,
|
|
|
emissions,
|
|
|
alpha_mix,
|
|
|
alpha_plant_emission,
|
|
|
alpha_tr_emission,
|
|
|
c_acq,
|
|
|
c_center_disp,
|
|
|
c_emission,
|
|
|
c_fix,
|
|
|
c_open,
|
|
|
c_plant_disp,
|
|
|
c_store,
|
|
|
c_tr,
|
|
|
c_var,
|
|
|
m_cap,
|
|
|
m_center_disp,
|
|
|
m_dist,
|
|
|
m_emission,
|
|
|
m_init,
|
|
|
m_plant_disp,
|
|
|
m_store,
|
|
|
)
|
|
|
end
|
|
|
|
|
|
# Write
|
|
|
# ==============================================================================
|
|
|
|
|
|
function write_json(data, filename)
|
|
|
json = dict()
|
|
|
json["parameters"] = dict("time horizon (years)" => data.T.stop)
|
|
|
json["products"] = dict(
|
|
|
r.name => dict(
|
|
|
"components" => [c.name for c in r.comp],
|
|
|
"disposal limit (tonne)" => [data.m_center_disp[r, t] for t in data.T],
|
|
|
"transportation cost (\$/km/tonne)" => [data.c_tr[r, t] for t in data.T],
|
|
|
"transportation emissions (tonne/km/tonne)" => dict(
|
|
|
s => [data.alpha_tr_emission[r, s, t] for t in data.T] for
|
|
|
s in data.emissions
|
|
|
),
|
|
|
) for r in data.products
|
|
|
)
|
|
|
json["centers"] = dict(
|
|
|
q.name => dict(
|
|
|
"latitude" => q.latitude,
|
|
|
"longitude" => q.longitude,
|
|
|
"output" => dict(
|
|
|
r.name => dict(
|
|
|
"initial amount (tonne)" => [
|
|
|
data.m_init[q, r, c, t] for c in r.comp, t in data.T
|
|
|
],
|
|
|
"disposal cost (\$/tonne)" =>
|
|
|
[data.c_center_disp[q, r, t] for t in data.T],
|
|
|
"storage cost (\$/tonne)" =>
|
|
|
[data.c_store[q, r, t] for t in data.T],
|
|
|
"storage limit (tonne)" => [data.m_store[q, r, t] for t in data.T],
|
|
|
"acquisition cost (\$/tonne)" =>
|
|
|
[data.c_acq[q, r, t] for t in data.T],
|
|
|
) for r in q.prod_out
|
|
|
),
|
|
|
) for q in data.centers
|
|
|
)
|
|
|
json["plants"] = dict(
|
|
|
p.name => dict(
|
|
|
"latitude" => p.latitude,
|
|
|
"longitude" => p.longitude,
|
|
|
"input" => p.prod_in.name,
|
|
|
"output" => dict(
|
|
|
r.name => dict(
|
|
|
"output matrix" => [
|
|
|
data.alpha_mix[p, r, c_in, c_out] for
|
|
|
c_in in p.prod_in.comp, c_out in r.comp
|
|
|
],
|
|
|
"disposal limit (tonne)" =>
|
|
|
[data.m_plant_disp[p, r, t] for t in data.T],
|
|
|
"disposal cost (\$/tonne)" =>
|
|
|
[data.c_plant_disp[p, r, t] for t in data.T],
|
|
|
) for r in p.prod_out
|
|
|
),
|
|
|
|
|
|
"Min_capacity" => data.min_cap[p],
|
|
|
"Min_FixOperatingCost" => data.min_OpC[p],
|
|
|
"Min_OpeningCost" => data.min_OpenC[p],
|
|
|
"Max_capacity" => data.max_cap[p],
|
|
|
"Max_FixOperatingCost" => data.max_OpC[p],
|
|
|
"Max_OpeningCost" => data.max_OpenC[p],
|
|
|
"fixed operating cost (\$)" => [data.c_fix[p, t] for t in data.T],
|
|
|
"variable operating cost (\$/tonne)" => [data.c_var[p, t] for t in data.T],
|
|
|
"opening cost (\$)" => [data.c_open[p, t] for t in data.T],
|
|
|
"capacity (tonne)" => data.m_cap[p],
|
|
|
"emissions (tonne/tonne)" => dict(
|
|
|
s => [data.alpha_plant_emission[p, s, t] for t in data.T] for
|
|
|
s in data.emissions
|
|
|
),
|
|
|
) for p in data.plants
|
|
|
)
|
|
|
json["emissions"] = dict(
|
|
|
s.name => dict(
|
|
|
"penalty (\$/tonne)" => [data.c_emission[s, t] for t in data.T],
|
|
|
"limit (tonne)" => [data.m_emission[s, t] for t in data.T],
|
|
|
) for s in data.emissions
|
|
|
)
|
|
|
open(filename, "w") do io
|
|
|
JSON.print(io, json, 2)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Read
|
|
|
# ==============================================================================
|
|
|
function read_json(filename; max_centers=Inf, max_plants=Inf)::Instance
|
|
|
json = JSON.parsefile(filename)
|
|
|
T = 1:json["parameters"]["time horizon (years)"]
|
|
|
centers = []
|
|
|
components_by_name = dict()
|
|
|
emissions = []
|
|
|
emissions_by_name = dict()
|
|
|
plants = []
|
|
|
products = []
|
|
|
products_by_name = dict()
|
|
|
alpha_mix = dict()
|
|
|
alpha_plant_emission = dict()
|
|
|
alpha_tr_emission = dict()
|
|
|
c_acq = dict()
|
|
|
#c_interv = dict()
|
|
|
c_center_disp = dict()
|
|
|
c_emission = dict()
|
|
|
c_fix = dict()
|
|
|
c_open = dict()
|
|
|
c_plant_disp = dict()
|
|
|
c_store = dict()
|
|
|
c_tr = dict()
|
|
|
c_var = dict()
|
|
|
m_cap = dict()
|
|
|
m_center_disp = dict()
|
|
|
m_emission = dict()
|
|
|
m_init = dict()
|
|
|
m_plant_disp = dict()
|
|
|
m_store = dict()
|
|
|
max_cap = dict()
|
|
|
max_OpC = dict()
|
|
|
max_OpenC = dict()
|
|
|
min_cap = dict()
|
|
|
min_OpC = dict()
|
|
|
min_OpenC = dict()
|
|
|
unit_open = dict()
|
|
|
unit_OpC = dict()
|
|
|
|
|
|
@timeit "Read: Emissions" begin
|
|
|
for (emission_name, emission_data) in json["emissions"]
|
|
|
s = Emission(emission_name)
|
|
|
emissions_by_name[emission_name] = s
|
|
|
push!(emissions, s)
|
|
|
for t in T
|
|
|
c_emission[s, t] = emission_data["penalty (\$/tonne)"][t]
|
|
|
m_emission[s, t] = emission_data["limit (tonne)"][t]
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "Read: Products" begin
|
|
|
for (name, prod_data) in json["products"]
|
|
|
print(name)
|
|
|
comp = []
|
|
|
for (comp_name) in prod_data["components"]
|
|
|
if comp_name ∉ keys(components_by_name)
|
|
|
print(comp_name)
|
|
|
components_by_name[comp_name] = Component(comp_name)
|
|
|
end
|
|
|
push!(comp, components_by_name[comp_name])
|
|
|
end
|
|
|
r = Product(name, comp)
|
|
|
products_by_name[name] = r
|
|
|
push!(products, r)
|
|
|
for t in T
|
|
|
m_center_disp[r, t] = prod_data["disposal limit (tonne)"][t]
|
|
|
c_tr[r, t] = prod_data["transportation cost (\$/km/tonne)"][t]
|
|
|
end
|
|
|
for (s_name, s_data) in prod_data["transportation emissions (tonne/km/tonne)"]
|
|
|
s = emissions_by_name[s_name]
|
|
|
for t in T
|
|
|
alpha_tr_emission[r, s, t] = s_data[t]
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "Read: Centers" begin
|
|
|
for (name, center_data) in json["centers"]
|
|
|
if length(centers) >= max_centers
|
|
|
@warn "Maximum number of centers reached. Skipping remaining ones."
|
|
|
break
|
|
|
end
|
|
|
latitude = center_data["latitude"]
|
|
|
longitude = center_data["longitude"]
|
|
|
prod_out = [products_by_name[r] for r in keys(center_data["output"])]
|
|
|
q = Center(; name, latitude, longitude, prod_out)
|
|
|
push!(centers, q)
|
|
|
for r in prod_out, t in T
|
|
|
c_acq[q, r, t] =
|
|
|
center_data["output"][r.name]["acquisition cost (\$/tonne)"][t]
|
|
|
#c_interv[q, r, t] =
|
|
|
# center_data["output"][r.name]["intervention costs ($)"][t]
|
|
|
c_center_disp[q, r, t] =
|
|
|
center_data["output"][r.name]["disposal cost (\$/tonne)"][t]
|
|
|
c_store[q, r, t] =
|
|
|
center_data["output"][r.name]["storage cost (\$/tonne)"][t]
|
|
|
m_store[q, r, t] = center_data["output"][r.name]["storage limit (tonne)"][t]
|
|
|
for (c_idx, c) in enumerate(r.comp)
|
|
|
m_init[q, r, c, t] =
|
|
|
center_data["output"][r.name]["initial amount (tonne)"][t][c_idx]
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "Read: Plants" begin
|
|
|
for (plant_name, plant_data) in json["plants"]
|
|
|
if length(plants) >= max_plants
|
|
|
@warn "Maximum number of plants reached. Skipping remaining ones."
|
|
|
break
|
|
|
end
|
|
|
#county = plant_data["county"]
|
|
|
latitude = plant_data["latitude"]
|
|
|
longitude = plant_data["longitude"]
|
|
|
prod_in = products_by_name[plant_data["input"]]
|
|
|
prod_out = [products_by_name[r] for r in keys(plant_data["output"])]
|
|
|
p = Plant(plant_name, latitude, longitude, prod_out, prod_in)
|
|
|
push!(plants, p)
|
|
|
m_cap[p] = plant_data["capacity (tonne)"]
|
|
|
min_cap[p] = plant_data["Min_capacity"]
|
|
|
min_OpC[p] = plant_data["Min_FixOperatingCost"]
|
|
|
min_OpenC[p] = plant_data["Min_OpeningCost"]
|
|
|
max_cap[p] = plant_data["Max_capacity"]
|
|
|
max_OpC[p] = plant_data["Max_FixOperatingCost"]
|
|
|
max_OpenC[p] = plant_data["Max_OpeningCost"]
|
|
|
|
|
|
#Compute unit cost for expansion calculations----------------------
|
|
|
#unit opening cost
|
|
|
if max_cap[p] > min_cap[p]
|
|
|
unit_open[p] = (max_OpenC[p] - min_OpenC[p])/(max_cap[p] - min_cap[p])
|
|
|
else
|
|
|
unit_open[p] = (max_OpenC[p]/max_cap[p])
|
|
|
end
|
|
|
#unit operating cost
|
|
|
if max_cap[p] > min_cap[p]
|
|
|
unit_OpC[p] = (max_OpC[p] - min_OpC[p])/(max_cap[p] - min_cap[p])
|
|
|
else
|
|
|
unit_OpC[p] = (max_OpC[p]/max_cap[p])
|
|
|
end
|
|
|
#------------------------------------------------------------------
|
|
|
|
|
|
for t in T
|
|
|
c_fix[p, t] = plant_data["fixed operating cost (\$)"][t]
|
|
|
c_var[p, t] = plant_data["variable operating cost (\$/tonne)"][t]
|
|
|
c_open[p, t] = plant_data["opening cost (\$)"][t]
|
|
|
end
|
|
|
for r in prod_out,
|
|
|
(cin_idx, c_in) in enumerate(prod_in.comp),
|
|
|
(cout_idx, c_out) in enumerate(r.comp)
|
|
|
|
|
|
alpha_mix[p, r, c_in, c_out] =
|
|
|
plant_data["output"][r.name]["output matrix"][cout_idx][cin_idx]
|
|
|
end
|
|
|
for r in prod_out, t in T
|
|
|
c_plant_disp[p, r, t] =
|
|
|
plant_data["output"][r.name]["disposal cost (\$/tonne)"][t]
|
|
|
m_plant_disp[p, r, t] =
|
|
|
plant_data["output"][r.name]["disposal limit (tonne)"][t]
|
|
|
end
|
|
|
for (s_name, s_data) in plant_data["emissions (tonne/tonne)"]
|
|
|
s = emissions_by_name[s_name]
|
|
|
for t in T
|
|
|
alpha_plant_emission[p, s, t] = s_data[t]
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "Calculate distances" begin
|
|
|
metric = KnnDrivingDistance()
|
|
|
m_dist = dict(
|
|
|
(p, q) => _calculate_distance(
|
|
|
p.latitude,
|
|
|
p.longitude,
|
|
|
q.latitude,
|
|
|
q.longitude,
|
|
|
metric,
|
|
|
) for p in [plants; centers], q in plants
|
|
|
)
|
|
|
end
|
|
|
|
|
|
return Instance(;
|
|
|
T,
|
|
|
centers,
|
|
|
plants,
|
|
|
products,
|
|
|
emissions,
|
|
|
alpha_mix,
|
|
|
alpha_plant_emission,
|
|
|
alpha_tr_emission,
|
|
|
c_acq,
|
|
|
#c_interv,
|
|
|
c_center_disp,
|
|
|
c_emission,
|
|
|
c_fix,
|
|
|
c_open,
|
|
|
c_plant_disp,
|
|
|
c_store,
|
|
|
c_tr,
|
|
|
c_var,
|
|
|
m_cap,
|
|
|
m_center_disp,
|
|
|
m_dist,
|
|
|
m_emission,
|
|
|
m_init,
|
|
|
m_plant_disp,
|
|
|
m_store,
|
|
|
max_cap,
|
|
|
max_OpC,
|
|
|
max_OpenC,
|
|
|
min_cap,
|
|
|
min_OpC,
|
|
|
min_OpenC,
|
|
|
unit_open,
|
|
|
unit_OpC
|
|
|
)
|
|
|
end
|
|
|
|
|
|
# Multi-period Heuristic
|
|
|
# ==============================================================================
|
|
|
|
|
|
function compress(original::Instance)::Instance
|
|
|
T = original.T
|
|
|
alpha_plant_emission = Dict(
|
|
|
(p, s, 1) => mean([original.alpha_plant_emission[p, s, t] for t in T]) for
|
|
|
p in original.plants, s in original.emissions
|
|
|
)
|
|
|
alpha_tr_emission = Dict(
|
|
|
(r, s, 1) => mean([original.alpha_tr_emission[r, s, t] for t in T]) for
|
|
|
r in original.products, s in original.emissions
|
|
|
)
|
|
|
c_acq = Dict(
|
|
|
(q, r, 1) => mean([original.c_acq[q, r, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out
|
|
|
)
|
|
|
#c_interv = Dict(
|
|
|
# (q, r, 1) => mean([original.c_interv[q, r, t] for t in T]) for
|
|
|
# q in original.centers for r in q.prod_out
|
|
|
#)
|
|
|
c_center_disp = Dict(
|
|
|
(q, r, 1) => mean([original.c_center_disp[q, r, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out
|
|
|
)
|
|
|
c_emission = Dict(
|
|
|
(s, 1) => mean([original.c_emission[s, t] for t in T]) for s in original.emissions
|
|
|
)
|
|
|
c_fix = Dict((p, 1) => sum([original.c_fix[p, t] for t in T]) for p in original.plants)
|
|
|
c_open =
|
|
|
Dict((p, 1) => mean([original.c_open[p, t] for t in T]) for p in original.plants)
|
|
|
c_plant_disp = Dict(
|
|
|
(p, r, 1) => mean([original.c_plant_disp[p, r, t] for t in T]) for
|
|
|
p in original.plants for r in p.prod_out
|
|
|
)
|
|
|
c_store = Dict(
|
|
|
(q, r, 1) => mean([original.c_store[q, r, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out
|
|
|
)
|
|
|
c_tr = Dict((r, 1) => mean([original.c_tr[r, t] for t in T]) for r in original.products)
|
|
|
c_var = Dict((p, 1) => mean([original.c_var[p, t] for t in T]) for p in original.plants)
|
|
|
m_cap = Dict(p => original.m_cap[p] * length(T) for p in original.plants)
|
|
|
max_cap = Dict(p => original.max_cap[p] * length(T) for p in original.plants)
|
|
|
min_cap = Dict(p => original.min_cap[p] * length(T) for p in original.plants)
|
|
|
|
|
|
max_OpC = Dict(p => original.max_OpC[p] * length(T) for p in original.plants)
|
|
|
min_OpC = Dict(p => original.min_OpC[p] * length(T) for p in original.plants)
|
|
|
|
|
|
max_OpenC = Dict(p => original.max_OpenC[p] * length(T) for p in original.plants)
|
|
|
min_OpenC = Dict(p => original.min_OpenC[p] * length(T) for p in original.plants)
|
|
|
|
|
|
unit_open = Dict(p => original.unit_open[p] * length(T) for p in original.plants)
|
|
|
#unit_open = Dict(p => ((max_OpenC[p] - min_OpenC[p])/(max_cap[p] - min_cap[p])) for p in original.plants)
|
|
|
#unit_open = Dict(p => original.unit_open[p] for p in original.plants)
|
|
|
|
|
|
unit_OpC = Dict(p => original.unit_OpC[p] * length(T) for p in original.plants)
|
|
|
#unit_OpC = Dict(p => ((max_OpC[p] - min_OpC[p])/(max_cap[p] - min_cap[p])) for p in original.plants)
|
|
|
#unit_OpC = Dict(p => original.unit_OpC[p] for p in original.plants)
|
|
|
|
|
|
m_center_disp = Dict(
|
|
|
(r, 1) => sum([original.m_center_disp[r, t] for t in T]) for r in original.products
|
|
|
)
|
|
|
m_emission = Dict(
|
|
|
(s, 1) => sum([original.m_emission[s, t] for t in T]) for s in original.emissions
|
|
|
)
|
|
|
m_init = Dict(
|
|
|
(q, r, c, 1) => length(T) * maximum([original.m_init[q, r, c, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out for c in r.comp
|
|
|
)
|
|
|
m_plant_disp = Dict(
|
|
|
(p, r, 1) => length(T) * minimum([original.m_plant_disp[p, r, t] for t in T]) for
|
|
|
p in original.plants for r in p.prod_out
|
|
|
)
|
|
|
m_store = Dict(
|
|
|
(q, r, 1) => length(T) * minimum([original.m_store[q, r, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out
|
|
|
)
|
|
|
#=
|
|
|
m_init = Dict(
|
|
|
(q, r, c, 1) => sum([original.m_init[q, r, c, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out for c in r.comp
|
|
|
)
|
|
|
m_plant_disp = Dict(
|
|
|
(p, r, 1) => sum([original.m_plant_disp[p, r, t] for t in T]) for
|
|
|
p in original.plants for r in p.prod_out
|
|
|
)
|
|
|
m_store = Dict(
|
|
|
(q, r, 1) => sum([original.m_store[q, r, t] for t in T]) for
|
|
|
q in original.centers for r in q.prod_out
|
|
|
)
|
|
|
=#
|
|
|
return Instance(;
|
|
|
T=1:1,
|
|
|
centers=original.centers,
|
|
|
plants=original.plants,
|
|
|
products=original.products,
|
|
|
emissions=original.emissions,
|
|
|
alpha_mix=original.alpha_mix,
|
|
|
alpha_plant_emission,
|
|
|
alpha_tr_emission,
|
|
|
c_acq,
|
|
|
#c_interv,
|
|
|
c_center_disp,
|
|
|
c_emission,
|
|
|
c_fix,
|
|
|
c_open,
|
|
|
c_plant_disp,
|
|
|
c_store,
|
|
|
c_tr,
|
|
|
c_var,
|
|
|
m_cap,
|
|
|
m_center_disp,
|
|
|
m_dist=original.m_dist,
|
|
|
m_emission,
|
|
|
m_init,
|
|
|
m_plant_disp,
|
|
|
m_store,
|
|
|
max_cap,
|
|
|
max_OpC,
|
|
|
max_OpenC,
|
|
|
min_cap,
|
|
|
min_OpC,
|
|
|
min_OpenC,
|
|
|
unit_open,
|
|
|
unit_OpC
|
|
|
)
|
|
|
end
|
|
|
|
|
|
function benchmark_compress(filename, optimizer; max_centers=[Inf], max_plants=[Inf])
|
|
|
stats = []
|
|
|
for mc in max_centers, mp in max_plants
|
|
|
# Solve original
|
|
|
orig = read_json(filename; max_centers=mc, max_plants=mp)
|
|
|
reset_timer!()
|
|
|
_, stats_orig = solve(orig; optimizer)
|
|
|
stats_orig["Filename"] = filename
|
|
|
stats_orig["Method"] = "Original"
|
|
|
push!(stats, stats_orig)
|
|
|
|
|
|
# Solve compressed
|
|
|
compressed = compress(orig)
|
|
|
reset_timer!()
|
|
|
model_comp, stats_comp = solve(compressed; optimizer)
|
|
|
stats_comp["Filename"] = filename
|
|
|
stats_comp["Method"] = "Compressed"
|
|
|
push!(stats, stats_comp)
|
|
|
|
|
|
end
|
|
|
|
|
|
return DataFrame(stats)
|
|
|
end
|
|
|
|
|
|
|
|
|
# Run
|
|
|
# ==============================================================================
|
|
|
|
|
|
function generate_json()
|
|
|
@info "Generating data"
|
|
|
data = generate_data()
|
|
|
|
|
|
@info "Writing JSON file"
|
|
|
write_json(data, "output-3/case.json")
|
|
|
end
|
|
|
|
|
|
function solve(filename, optimizer; max_centers=Inf, max_plants=Inf, heuristic=false)
|
|
|
reset_timer!()
|
|
|
@timeit "Read JSON" begin
|
|
|
data = read_json(filename; max_centers, max_plants)
|
|
|
write_json(data, "output-3/BaselinePlastics_capExTestCase10C7P.json")
|
|
|
end
|
|
|
return solve(data; optimizer=optimizer, output_dir=dirname(filename), heuristic)
|
|
|
print_timer()
|
|
|
end
|
|
|
|
|
|
function solve(data::Instance; optimizer, output_dir=nothing, heuristic=false)
|
|
|
|
|
|
if heuristic
|
|
|
@info "Solving compressed instance..."
|
|
|
comp_data = compress(data)
|
|
|
#print("data compressed")
|
|
|
comp_model, comp_stats = solve(comp_data; optimizer, heuristic=false)
|
|
|
#print("compressed model solved")
|
|
|
|
|
|
# Filter plants
|
|
|
selected_plants = Plant[]
|
|
|
for p in comp_data.plants
|
|
|
if value(comp_model[:x_open][p, 1]) > 0.5
|
|
|
push!(selected_plants, p)
|
|
|
end
|
|
|
end
|
|
|
@info "Selected $(length(selected_plants)) out of $(length(comp_data.plants)) plants"
|
|
|
|
|
|
|
|
|
# Filter edges
|
|
|
selected_edges = Set()
|
|
|
for (src, dst, r) in comp_model[:E]
|
|
|
if value(comp_model[:y_total][src, dst, r, 1]) > 0.5
|
|
|
push!(selected_edges, (src, dst, r))
|
|
|
end
|
|
|
end
|
|
|
@info "Selected $(length(selected_edges)) out of $(length(comp_model[:E])) transportation edges"
|
|
|
|
|
|
|
|
|
data = Instance(;
|
|
|
data.T,
|
|
|
data.centers,
|
|
|
plants=selected_plants,
|
|
|
data.products,
|
|
|
data.emissions,
|
|
|
data.alpha_mix,
|
|
|
data.alpha_plant_emission,
|
|
|
data.alpha_tr_emission,
|
|
|
data.c_acq,
|
|
|
#data.c_interv,
|
|
|
data.c_center_disp,
|
|
|
data.c_emission,
|
|
|
data.c_fix,
|
|
|
data.c_open,
|
|
|
data.c_plant_disp,
|
|
|
data.c_store,
|
|
|
data.c_tr,
|
|
|
data.c_var,
|
|
|
data.m_cap,
|
|
|
data.m_center_disp,
|
|
|
data.m_dist,
|
|
|
data.m_emission,
|
|
|
data.m_init,
|
|
|
data.m_plant_disp,
|
|
|
data.m_store,
|
|
|
selected_edges,
|
|
|
data.max_cap,
|
|
|
data.max_OpC,
|
|
|
data.max_OpenC,
|
|
|
data.min_cap,
|
|
|
data.min_OpC,
|
|
|
data.min_OpenC,
|
|
|
data.unit_open,
|
|
|
data.unit_OpC
|
|
|
|
|
|
)
|
|
|
end
|
|
|
|
|
|
|
|
|
T = data.T
|
|
|
centers = data.centers
|
|
|
plants = data.plants
|
|
|
products = data.products
|
|
|
emissions = data.emissions
|
|
|
|
|
|
#model = Model(optimizer)
|
|
|
model = Model(Gurobi.Optimizer)
|
|
|
set_optimizer_attribute(model, "MIPGap", 0.005)
|
|
|
init_time = time()
|
|
|
|
|
|
# Graph
|
|
|
# -------------------------------------------------------------------------
|
|
|
@timeit "Build graph" begin
|
|
|
E = []
|
|
|
E_in = dict(src => [] for src in plants)
|
|
|
E_out = dict(src => [] for src in plants ∪ centers)
|
|
|
function push_edge!(src, dst, r)
|
|
|
e = (src, dst, r)
|
|
|
if data.selected_edges !== nothing && e ∉ data.selected_edges
|
|
|
return
|
|
|
end
|
|
|
push!(E, e)
|
|
|
push!(E_out[src], (dst, r))
|
|
|
push!(E_in[dst], (src, r))
|
|
|
end
|
|
|
for r in products
|
|
|
# Plant to plant
|
|
|
for p1 in plants
|
|
|
r ∈ p1.prod_out || continue
|
|
|
for p2 in plants
|
|
|
p1 != p2 || continue
|
|
|
r == p2.prod_in || continue
|
|
|
push_edge!(p1, p2, r)
|
|
|
end
|
|
|
end
|
|
|
# Center to plant
|
|
|
for q in centers
|
|
|
r ∈ q.prod_out || continue
|
|
|
for p in plants
|
|
|
r == p.prod_in || continue
|
|
|
push_edge!(q, p, r)
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
model[:E] = E
|
|
|
model[:E_in] = E_in
|
|
|
model[:E_out] = E_out
|
|
|
|
|
|
@printf("Building optimization problem with:\n")
|
|
|
@printf(" %8d plants\n", length(plants))
|
|
|
@printf(" %8d centers\n", length(centers))
|
|
|
@printf(" %8d products\n", length(products))
|
|
|
@printf(" %8d time periods\n", length(T))
|
|
|
@printf(" %8d transportation edges\n", length(E))
|
|
|
|
|
|
available = Dict((r, t) => 0.0 for r in data.products, t in T)
|
|
|
for q in centers, t in T
|
|
|
for r in q.prod_out, s in r.comp
|
|
|
available[r, t] += data.m_init[q, r, s, t]
|
|
|
end
|
|
|
end
|
|
|
|
|
|
capacity = Dict(r => 0.0 for r in data.products)
|
|
|
for p in plants
|
|
|
capacity[p.prod_in] += data.max_cap[p]
|
|
|
end
|
|
|
|
|
|
for r in products, t in T
|
|
|
if available[r, t] > capacity[r]
|
|
|
@warn "Not enough capacity to process $(r.name) at time $t: $(available[r,t]) > $(capacity[r])"
|
|
|
end
|
|
|
end
|
|
|
|
|
|
|
|
|
# Decision variables
|
|
|
# -------------------------------------------------------------------------
|
|
|
@timeit "Model: Add variables" begin
|
|
|
@timeit "y" begin
|
|
|
y = _init(model, :y)
|
|
|
for (q, p, r) in E, c in r.comp, t in T
|
|
|
y[q, p, r, c ,t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "y_total" begin
|
|
|
y_total = _init(model, :y_total)
|
|
|
for (q, p, r) in E, t in T
|
|
|
y_total[q, p, r, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_center_disp" begin
|
|
|
z_center_disp = _init(model, :z_center_disp)
|
|
|
for q in centers, r in q.prod_out, c in r.comp, t in T
|
|
|
z_center_disp[q, r, c, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_center_disp_total" begin
|
|
|
z_center_disp_total = _init(model, :z_center_disp_total)
|
|
|
for q in centers, r in q.prod_out, t in T
|
|
|
z_center_disp_total[q, r, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_store" begin
|
|
|
z_store = _init(model, :z_store)
|
|
|
for q in centers, r in q.prod_out, c in r.comp, t in T
|
|
|
if t == T.stop
|
|
|
z_store[q, r, c, t] = 0.0
|
|
|
else
|
|
|
z_store[q, r, c, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_store_total" begin
|
|
|
z_store_total = _init(model, :z_store_total)
|
|
|
for q in centers, r in q.prod_out
|
|
|
z_store_total[q, r, 0] = 0.0
|
|
|
for t in T
|
|
|
if t == T.stop
|
|
|
z_store_total[q, r, t] = 0.0
|
|
|
else
|
|
|
z_store_total[q, r, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
#=
|
|
|
@timeit "x_open" begin
|
|
|
x_open = _init(model, :x_open)
|
|
|
for p in plants
|
|
|
x_open[p, 0] = 0
|
|
|
for t in T
|
|
|
x_open[p, t] = @variable(model, binary = true)
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
=#
|
|
|
|
|
|
@timeit "x_open" begin
|
|
|
x_open = _init(model, :x_open)
|
|
|
for p in plants
|
|
|
if data.m_cap[p]>0
|
|
|
x_open[p,0] = 1
|
|
|
else
|
|
|
x_open[p,0] = 0
|
|
|
end
|
|
|
for t in T
|
|
|
x_open[p, t] = @variable(model, binary = true)
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "x_send" begin
|
|
|
x_send = _init(model, :x_send)
|
|
|
for p in plants, (q, r) in E_out[p], t in T
|
|
|
x_send[p, q, r, t] = @variable(model, binary = true)
|
|
|
end
|
|
|
end
|
|
|
@timeit "x_disp" begin
|
|
|
x_disp = _init(model, :x_disp)
|
|
|
for p in plants, r in p.prod_out, t in T
|
|
|
x_disp[p, r, t] = @variable(model, binary = true)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_prod" begin
|
|
|
z_prod = _init(model, :z_prod)
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
z_prod[p, r, c, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_plant_disp" begin
|
|
|
z_plant_disp = _init(model, :z_plant_disp)
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
z_plant_disp[p, r, c, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_tr_emissions" begin
|
|
|
z_tr_emissions = _init(model, :z_tr_emissions)
|
|
|
for (q, p, r) in E, s in emissions, t in T
|
|
|
z_tr_emissions[q, p, r, s, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
@timeit "z_plant_emissions" begin
|
|
|
z_plant_emissions = _init(model, :z_plant_emissions)
|
|
|
for p in plants, s in emissions, t in T
|
|
|
z_plant_emissions[p, s, t] = @variable(model, lower_bound = 0)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
########################## FOR CAPACITY EXPANSION ##############################################################
|
|
|
#Binary variable - if a new plant is built
|
|
|
@timeit "x_build" begin
|
|
|
x_build = _init(model, :x_build)
|
|
|
for p in plants, t in T
|
|
|
x_build[p, t] = @variable(model, binary = true)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Plant capacity expansion Variable
|
|
|
@timeit "z_capExp" begin
|
|
|
z_capExp = _init(model, :z_capExp)
|
|
|
for p in plants, t in T
|
|
|
z_capExp[p, t] = @variable(model, lower_bound = 0, upper_bound = (data.max_cap[p] - data.m_cap[p]))
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Total plant capacity Variable
|
|
|
@timeit "z_capTot" begin
|
|
|
z_capTot = _init(model, :z_capTot)
|
|
|
for p in plants
|
|
|
z_capTot[p, 0] = data.m_cap[p]
|
|
|
for t in T
|
|
|
z_capTot[p, t] = @variable(model, lower_bound = 0, upper_bound = data.max_cap[p])
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
##############################################################################################################
|
|
|
end
|
|
|
|
|
|
# Objective function
|
|
|
# -------------------------------------------------------------------------
|
|
|
@timeit "Model: Objective function" begin
|
|
|
obj = AffExpr()
|
|
|
|
|
|
# Center disposal
|
|
|
@timeit "c_center_disp" begin
|
|
|
for q in centers, r in q.prod_out, t in T
|
|
|
add_to_expression!(
|
|
|
obj,
|
|
|
data.c_center_disp[q, r, t],
|
|
|
z_center_disp_total[q, r, t],
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Center acquisition
|
|
|
@timeit "c_acq" begin
|
|
|
for q in centers, r in q.prod_out, c in r.comp, t in T
|
|
|
add_to_expression!(obj, data.m_init[q, r, c, t] * data.c_acq[q, r, t]) #+data.c_interv[q,r,t])
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Center storage
|
|
|
@timeit "c_store" begin
|
|
|
for q in centers, r in q.prod_out, t in T
|
|
|
add_to_expression!(obj, data.c_store[q, r, t], z_store_total[q, r, t])
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Transportation, variable operating cost
|
|
|
@timeit "c_tr, c_var" begin
|
|
|
for (q, p, r) in E
|
|
|
dist = data.m_dist[q, p]
|
|
|
#y_qpr = y[q, p, r]
|
|
|
for c in r.comp, t in T
|
|
|
add_to_expression!(obj, dist * data.c_tr[r, t], y[q, p, r,c, t])
|
|
|
add_to_expression!(obj, data.c_var[p, t], y[q, p, r, c, t])
|
|
|
end
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Transportation emissions
|
|
|
@timeit "c_emission" begin
|
|
|
for (q, p, r) in E, s in emissions, t in T
|
|
|
add_to_expression!(
|
|
|
obj,
|
|
|
data.c_emission[s, t],
|
|
|
z_tr_emissions[q, p, r, s, t],
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Fixed cost
|
|
|
@timeit "c_fix" begin
|
|
|
for p in plants, t in T
|
|
|
add_to_expression!(obj, data.c_fix[p, t], x_open[p, t])
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#=
|
|
|
# Opening cost
|
|
|
@timeit "c_open" begin
|
|
|
for p in plants, t in T
|
|
|
add_to_expression!(obj, data.c_open[p, t], x_open[p, t] - x_open[p, t-1])
|
|
|
end
|
|
|
end
|
|
|
=#
|
|
|
|
|
|
# Plant disposal
|
|
|
@timeit "c_plant_disp" begin
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
add_to_expression!(
|
|
|
obj,
|
|
|
data.c_plant_disp[p, r, t],
|
|
|
z_plant_disp[p, r, c, t],
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# Plant emissions
|
|
|
@timeit "c_emission" begin
|
|
|
for p in plants, s in emissions, t in T
|
|
|
add_to_expression!(obj, data.c_emission[s, t], z_plant_emissions[p, s, t])
|
|
|
end
|
|
|
end
|
|
|
|
|
|
######################## FOR CAPACITY EXPANSION ####################################################
|
|
|
|
|
|
#Plant opening cost
|
|
|
@timeit "c_open" begin
|
|
|
for p in plants, t in T
|
|
|
add_to_expression!(
|
|
|
obj,
|
|
|
data.min_OpenC[p],
|
|
|
x_build[p, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Plant Expansion cost
|
|
|
@timeit "c_exp" begin
|
|
|
for p in plants, t in T
|
|
|
add_to_expression!(
|
|
|
obj,
|
|
|
data.unit_open[p],
|
|
|
(z_capExp[p, t])
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Plant additional fixed operating cost due to expansion
|
|
|
@timeit "c_fix_exp" begin
|
|
|
for p in plants, t in T
|
|
|
add_to_expression!(
|
|
|
obj,
|
|
|
data.unit_OpC[p],
|
|
|
sum(z_capExp[p, i] for i = 1:t)
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
####################################################################################################
|
|
|
|
|
|
@objective(model, Min, obj)
|
|
|
end
|
|
|
|
|
|
# Constraints
|
|
|
# -------------------------------------------------------------------------
|
|
|
@timeit "Model: Constraints" begin
|
|
|
@timeit "eq_balance" begin
|
|
|
eq_balance = _init(model, :eq_balance)
|
|
|
for q in centers, r in q.prod_out, t in T
|
|
|
eq_balance[q, r, t] = @constraint(
|
|
|
model,
|
|
|
sum(y_total[q, p, r2, t] for (p, r2) in E_out[q] if r == r2) +
|
|
|
z_store_total[q, r, t] + sum(z_center_disp[q, r, c, t] for c in r.comp) ==
|
|
|
sum(data.m_init[q, r, c, t] for c in r.comp) + z_store_total[q, r, t-1]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
ratio(q, r, c, t) =
|
|
|
(data.m_init[q, r, c, t] / sum(data.m_init[q, r, d, t] for d in r.comp))
|
|
|
|
|
|
@timeit "eq_y_total" begin
|
|
|
eq_y_total = _init(model, :eq_y_total)
|
|
|
for (q, p, r) in E, t in T
|
|
|
eq_y_total[q, p, r, t] = @constraint(
|
|
|
model,
|
|
|
y_total[q, p, r, t] == sum(y[q, p, r, c, t] for c in r.comp)
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_split_y" begin
|
|
|
eq_split_y = _init(model, :eq_split_y)
|
|
|
for (q, p, r) in E, c in r.comp, t in T
|
|
|
if q ∉ centers
|
|
|
continue
|
|
|
end
|
|
|
eq_split_y[q, p, r, c, t] = @constraint(
|
|
|
model,
|
|
|
y[q, p, r, c, t] == ratio(q, r, c, t) * y_total[q, p, r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_split_center_disp" begin
|
|
|
eq_split_center_disp = _init(model, :eq_split_center_disp)
|
|
|
for q in centers, r in q.prod_out, c in r.comp, t in T
|
|
|
eq_split_center_disp[q, r, c, t] = @constraint(
|
|
|
model,
|
|
|
z_center_disp[q, r, c, t] ==
|
|
|
ratio(q, r, c, t) * z_center_disp_total[q, r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_split_store" begin
|
|
|
eq_split_store = _init(model, :eq_split_store)
|
|
|
for q in centers, r in q.prod_out, c in r.comp, t in T
|
|
|
eq_split_store[q, r, c, t] = @constraint(
|
|
|
model,
|
|
|
z_store[q, r, c, t] == ratio(q, r, c, t) * z_store_total[q, r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "eq_center_disposal" begin
|
|
|
eq_center_disposal = _init(model, :eq_center_disposal)
|
|
|
for r in products, t in T
|
|
|
centers_r = [q for q in centers if r ∈ q.prod_out]
|
|
|
if isempty(centers_r)
|
|
|
continue
|
|
|
end
|
|
|
eq_center_disposal[r, t] = @constraint(
|
|
|
model,
|
|
|
sum(z_center_disp_total[q, r, t] for q in centers_r) <=
|
|
|
data.m_center_disp[r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "eq_center_storage" begin
|
|
|
eq_center_storage = _init(model, :eq_center_storage)
|
|
|
for q in centers, r in q.prod_out, t in T
|
|
|
eq_center_storage[q, r, t] =
|
|
|
@constraint(model, z_store_total[q, r, t] <= data.m_store[q, r, t])
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_tr_emissions" begin
|
|
|
eq_tr_emissions = _init(model, :eq_tr_emissions)
|
|
|
for (q, p, r) in E, s in emissions, t in T
|
|
|
eq_tr_emissions[q, p, r, s, t] = @constraint(
|
|
|
model,
|
|
|
z_tr_emissions[q, p, r, s, t] ==
|
|
|
data.m_dist[q, p] *
|
|
|
data.alpha_tr_emission[r, s, t] *
|
|
|
y_total[q, p, r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#=
|
|
|
@timeit "eq_plant_capacity" begin
|
|
|
eq_plant_capacity = _init(model, :eq_plant_capacity)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_capacity[p, t] = @constraint(
|
|
|
model,
|
|
|
sum(y_total[q, p, r][t] for (q, r) in E_in[p]) <=
|
|
|
data.m_cap[p] * x_open[p, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
=#
|
|
|
|
|
|
@timeit "eq_plant_prod" begin
|
|
|
eq_plant_prod = _init(model, :eq_plant_prod)
|
|
|
for p in plants, r_out in p.prod_out, c_out in r_out.comp, t in T
|
|
|
eq_plant_prod[p, r_out, c_out, t] = @constraint(
|
|
|
model,
|
|
|
z_prod[p, r_out, c_out, t] == sum(
|
|
|
data.alpha_mix[p, r_out, c_in, c_out] * y[q, p, r_in, c_in, t]
|
|
|
for (q, r_in) in E_in[p] for
|
|
|
c_in in r_in.comp if data.alpha_mix[p, r_out, c_in, c_out] > 0
|
|
|
)
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_plant_disp" begin
|
|
|
eq_plant_disp = _init(model, :eq_plant_disp)
|
|
|
for p in plants, r in p.prod_out, t in T
|
|
|
eq_plant_disp[p, r, t] = @constraint(
|
|
|
model,
|
|
|
sum(z_plant_disp[p, r, c, t] for c in r.comp) <=
|
|
|
data.m_plant_disp[p, r, t] * x_disp[p, r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_plant_emissions" begin
|
|
|
eq_plant_emissions = _init(model, :eq_plant_emissions)
|
|
|
for p in plants, s in emissions, t in T
|
|
|
eq_plant_emissions[p, s, t] = @constraint(
|
|
|
model,
|
|
|
z_plant_emissions[p, s, t] ==
|
|
|
sum(y_total[q, p, r, t] for (q, r) in E_in[p]) *
|
|
|
data.alpha_plant_emission[p, s, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
@timeit "eq_emissions_limit" begin
|
|
|
eq_emissions_limit = _init(model, :eq_emissions_limit)
|
|
|
for s in emissions, t in T
|
|
|
eq_emissions_limit[s, t] = @constraint(
|
|
|
model,
|
|
|
sum(z_plant_emissions[p, s, t] for p in plants) +
|
|
|
sum(z_tr_emissions[q, p, r, s, t] for (q, p, r) in E) <=
|
|
|
data.m_emission[s, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "eq_plant_remains_open" begin
|
|
|
eq_plant_remains_open = _init(model, :eq_plant_remains_open)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_remains_open[p, t] =
|
|
|
@constraint(model, x_open[p, t] >= x_open[p, t-1])
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#=
|
|
|
@timeit "eq_plant_single_dest" begin
|
|
|
eq_plant_single_dest = _init(model, :eq_plant_single_dest)
|
|
|
for p in plants, r in p.prod_out, t in T
|
|
|
eq_plant_single_dest[p, r, t] = @constraint(
|
|
|
model,
|
|
|
sum(x_send[p, q, r, t] for (q, r2) in E_out[p] if r == r2) +
|
|
|
x_disp[p, r, t] <= 1
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
@timeit "eq_plant_send_limit" begin
|
|
|
eq_plant_send_limit = _init(model, :eq_plant_send_limit)
|
|
|
for p in plants, (q, r) in E_out[p], t in T
|
|
|
eq_plant_send_limit[p, q, r, t] = @constraint(
|
|
|
model,
|
|
|
y_total[p, q, r][t] <= data.m_cap[q] * x_send[p, q, r, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
=#
|
|
|
|
|
|
@timeit "eq_plant_balance" begin
|
|
|
eq_plant_balance = _init(model, :eq_plant_balance)
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
eq_plant_balance[p, r, c, t] = @constraint(
|
|
|
model,
|
|
|
z_prod[p, r, c, t] ==
|
|
|
z_plant_disp[p, r, c, t] +
|
|
|
sum(y[p, q, r, c, t] for (q, r2) in E_out[p] if r == r2)
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
############################################### FOR CAPACITY EXPANSION ###########################################
|
|
|
#Total amount send to a plant should satisfy the capacity constraints of the plant
|
|
|
|
|
|
@timeit "eq_plant_capacity" begin
|
|
|
eq_plant_capacity = _init(model, :eq_plant_capacity)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_capacity[p, t] = @constraint(
|
|
|
model,
|
|
|
sum(y_total[q, p, r, t] for (q, r) in E_in[p]) <=
|
|
|
z_capTot[p, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
# total plant capacity
|
|
|
@timeit "eq_plant_totalcapacity" begin
|
|
|
eq_plant_totalcapacity = _init(model, :eq_plant_totalcapacity)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_totalcapacity[p, t] = @constraint(
|
|
|
model,
|
|
|
z_capTot[p, t] == z_capTot[p, t-1] + z_capExp[p, t] + (data.min_cap[p] * x_build[p, t])
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Expansion capacity limitation constraints for plants
|
|
|
@timeit "eq_plant_Expcapacity" begin
|
|
|
eq_plant_ExpCapacity = _init(model, :eq_plant_Expcapacity)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_ExpCapacity[p, t] = @constraint(
|
|
|
model,
|
|
|
z_capExp[p, t] <= (data.max_cap[p] - data.m_cap[p]) * x_open[p, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Plant Feasibility constraint - A plant is open iff it was open at t-1or new plant is built at t
|
|
|
@timeit "eq_plant_feasibility" begin
|
|
|
eq_plant_feasibility = _init(model, :eq_plant_feasibility)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_feasibility[p, t] = @constraint(
|
|
|
model,
|
|
|
x_open[p, t] == x_open[p, t-1] + x_build[p, t]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Plant Initial opening condition Constraints
|
|
|
@timeit "eq_plant_opencondition" begin
|
|
|
eq_plant_opencondition = _init(model, :eq_plant_opencondition)
|
|
|
for p in plants, t in T
|
|
|
eq_plant_opencondition[p, t] = @constraint(
|
|
|
model,
|
|
|
x_build[p, t] <= 1 - x_open[p, t-1]
|
|
|
)
|
|
|
end
|
|
|
end
|
|
|
|
|
|
#Total capacity can only increase or stay the same
|
|
|
@timeit "eq_plant_totalcap_continuity" begin
|
|
|
eq_plant_totalcap_continuity = _init(model, :eq_plant_totalcap_continuity)
|
|
|
for p in plants, t in T
|
|
|
#if t>1
|
|
|
eq_plant_totalcap_continuity[p, t] = @constraint(
|
|
|
model,
|
|
|
z_capTot[p, t] >= z_capTot[p, t-1]
|
|
|
)
|
|
|
#end
|
|
|
end
|
|
|
end
|
|
|
################################################################################################################
|
|
|
|
|
|
end
|
|
|
|
|
|
# Optimize
|
|
|
# -------------------------------------------------------------------------
|
|
|
_set_names!(model)
|
|
|
|
|
|
model_build_time = time() - init_time
|
|
|
#write_to_file(model, "model.lp")
|
|
|
#return
|
|
|
#error("STOP")
|
|
|
|
|
|
optimize!(model)
|
|
|
|
|
|
stats = OrderedDict{String,Any}(
|
|
|
"Plants" => length(plants),
|
|
|
"Centers" => length(centers),
|
|
|
"Products" => length(products),
|
|
|
"Time periods" => length(T),
|
|
|
"Model build time (s)" => model_build_time,
|
|
|
"Variables" => num_variables(model),
|
|
|
"Constraints" =>
|
|
|
num_constraints(model, count_variable_in_set_constraints=false),
|
|
|
"Objective Value" => objective_value(model),
|
|
|
"Solve time (s)" => solve_time(model),
|
|
|
)
|
|
|
|
|
|
write_to_file(model, "model.lp")
|
|
|
|
|
|
#open("model.lp", "w") do f
|
|
|
#print(f, model)
|
|
|
#end
|
|
|
#write_to_file(model, "model.ilp")
|
|
|
|
|
|
#=
|
|
|
#Add check for composition consistency
|
|
|
#-----------------------------------------------------------------------------
|
|
|
pdt_ratio = dict((p, r, c, t) => 0.0 for p in plants, r in p.prod_out, c in r.comp, t in data.T)
|
|
|
transp_ratio = dict((p, q, r, c, t) => 0.0 for p in plants, (q,r) in E_out[p], c in r.comp, t in T)
|
|
|
|
|
|
#compute the proportion of components produced in the plants
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
pdt_ratio[p,r,c,t] = value(z_prod[p,r,c,t])/sum(value(z_prod[p,r,d,t]) for d in r.comp)
|
|
|
end
|
|
|
#compute the proportion of components transported from a plant to another plant
|
|
|
for p in plants, (q,r) in E_out[p], c in r.comp, t in T
|
|
|
transp_ratio[p,q,r,c,t] = y[q,p,r][c,t]/sum(value(y[q,p,r][d,t]) for d in r.comp)
|
|
|
end
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
if (pdt_ratio[p,r,c,t] - transp_ratio[q,p,r,c,t] for (q, r2) in E_out[p] if r == r2) >0.1
|
|
|
@warn "composition inconsistency detected for $r.name transported from $p.name to $q.name at time $t"
|
|
|
end
|
|
|
end
|
|
|
=#
|
|
|
if output_dir === nothing
|
|
|
return model, stats
|
|
|
end
|
|
|
|
|
|
# Report: Transportation
|
|
|
# -------------------------------------------------------------------------
|
|
|
df = DataFrame()
|
|
|
df."source" = String[]
|
|
|
df."source longitude (deg)" = Float64[]
|
|
|
df."source latitude (deg)" = Float64[]
|
|
|
df."destination" = String[]
|
|
|
df."destination longitude (deg)" = Float64[]
|
|
|
df."destination latitude (deg)" = Float64[]
|
|
|
df."product" = String[]
|
|
|
df."component" = String[]
|
|
|
df."time" = Int[]
|
|
|
df."distance (km)" = Float64[]
|
|
|
df."amount sent (tonne)" = Float64[]
|
|
|
df."transportation cost (\$)" = Float64[]
|
|
|
df."variable operating cost (\$)" = Float64[]
|
|
|
for (q, p, r) in E, c in r.comp, t in T
|
|
|
if value(y[q, p, r, c, t]) ≈ 0
|
|
|
continue
|
|
|
end
|
|
|
push!(
|
|
|
df,
|
|
|
[
|
|
|
q.name,
|
|
|
q.longitude,
|
|
|
q.latitude,
|
|
|
p.name,
|
|
|
p.longitude,
|
|
|
p.latitude,
|
|
|
r.name,
|
|
|
c.name,
|
|
|
t,
|
|
|
round(data.m_dist[q, p], digits=2),
|
|
|
round(value(y[q, p, r, c, t]), digits=2),
|
|
|
round(
|
|
|
data.m_dist[q, p] * data.c_tr[r, t] * value(y[q, p, r, c, t]),
|
|
|
digits=2,
|
|
|
),
|
|
|
round(data.c_var[p, t] * value(y[q, p, r, c, t]), digits=2),
|
|
|
],
|
|
|
)
|
|
|
end
|
|
|
CSV.write("$output_dir/transp.csv", df)
|
|
|
|
|
|
# Report: Centers
|
|
|
# -------------------------------------------------------------------------
|
|
|
df = DataFrame()
|
|
|
df."center" = String[]
|
|
|
df."product" = String[]
|
|
|
df."component" = String[]
|
|
|
df."time" = Int[]
|
|
|
df."amount available (tonne)" = Float64[]
|
|
|
df."amount sent (tonne)" = Float64[]
|
|
|
df."amount stored (tonne)" = Float64[]
|
|
|
df."amount disposed (tonne)" = Float64[]
|
|
|
df."acquisition cost (\$)" = Float64[]
|
|
|
df."storage cost (\$)" = Float64[]
|
|
|
df."disposal cost (\$)" = Float64[]
|
|
|
for q in centers, r in q.prod_out, c in r.comp, t in T
|
|
|
push!(
|
|
|
df,
|
|
|
[
|
|
|
q.name,
|
|
|
r.name,
|
|
|
c.name,
|
|
|
t,
|
|
|
round(data.m_init[q, r, c, t], digits=2),
|
|
|
round(
|
|
|
sum(value(y[q, p, r, c, t]) for (p, r2) in E_out[q] if r == r2; init=0),
|
|
|
digits=2,
|
|
|
),
|
|
|
round(value(z_store[q, r, c, t]), digits=2),
|
|
|
round(value(z_center_disp[q, r, c, t]), digits=2),
|
|
|
round(data.m_init[q, r, c, t] * data.c_acq[q, r, t], digits=2),
|
|
|
round(data.c_store[q, r, t] * value(z_store[q, r, c, t]), digits=2),
|
|
|
round(
|
|
|
data.c_center_disp[q, r, t] * value(z_center_disp[q, r, c, t]),
|
|
|
digits=2,
|
|
|
),
|
|
|
],
|
|
|
)
|
|
|
end
|
|
|
CSV.write("$output_dir/centers.csv", df)
|
|
|
|
|
|
# Report: Plants
|
|
|
# -------------------------------------------------------------------------
|
|
|
df = DataFrame()
|
|
|
df."plant" = String[]
|
|
|
df."Plant latitude (deg)" = Float64[]
|
|
|
df."Plant longitude (deg)" = Float64[]
|
|
|
df."Maximum Capacity" = Float64[]
|
|
|
df."time(year)" = Int[]
|
|
|
df."is open?" = Float64[]
|
|
|
df."is newly built" = Float64[]
|
|
|
df."Initial Capacity" = Float64[]
|
|
|
df."Expanded Capacity" = Float64[]
|
|
|
df."Total Capacity" = Float64[]
|
|
|
df."opening cost (\$)" = Float64[]
|
|
|
df."Expansion Cost" = Float64[]
|
|
|
df."fixed operating cost (\$)" = Float64[]
|
|
|
df."Exp fixed op cost" = Float64[]
|
|
|
df."Total fixed operating cost" = Float64[]
|
|
|
for p in plants, t in T
|
|
|
push!(
|
|
|
df,
|
|
|
[
|
|
|
p.name,
|
|
|
p.latitude,
|
|
|
p.longitude,
|
|
|
data.max_cap[p],
|
|
|
t,
|
|
|
value(x_open[p, t]),
|
|
|
value(x_build[p, t]),
|
|
|
data.m_cap[p],
|
|
|
value(z_capExp[p, t]),
|
|
|
value(z_capTot[p, t]),
|
|
|
#data.c_open[p, t] * (value(x_open[p, t]) - value(x_open[p, t-1])),
|
|
|
data.c_open[p, t] * value(x_build[p, t]),
|
|
|
data.unit_open[p] * value(z_capExp[p, t]),
|
|
|
data.c_fix[p, t] * value(x_open[p, t]),
|
|
|
data.unit_OpC[p] * sum(value(z_capExp[p, i]) for i in 1:t),
|
|
|
((data.c_fix[p, t] * value(x_open[p, t])) + (data.unit_OpC[p] * sum(value(z_capExp[p, i]) for i in 1:t))),
|
|
|
]
|
|
|
)
|
|
|
end
|
|
|
CSV.write("$output_dir/plants.csv", df)
|
|
|
|
|
|
# Report: Plant Outputs
|
|
|
# -------------------------------------------------------------------------
|
|
|
df = DataFrame()
|
|
|
df."plant" = String[]
|
|
|
#df."county" = String[]
|
|
|
df."Plant Longitude (deg)" = Float64[]
|
|
|
df."Plant Latitude (deg)" = Float64[]
|
|
|
df."product" = String[]
|
|
|
df."component" = String[]
|
|
|
df."time" = Int[]
|
|
|
df."amount produced (tonne)" = Float64[]
|
|
|
df."amount disposed (tonne)" = Float64[]
|
|
|
df."amount sent (tonne)" = Float64[]
|
|
|
df."disposal cost (\$)" = Float64[]
|
|
|
for p in plants, r in p.prod_out, c in r.comp, t in T
|
|
|
if value(z_prod[p, r, c, t]) ≈ 0
|
|
|
continue
|
|
|
end
|
|
|
push!(
|
|
|
df,
|
|
|
[
|
|
|
p.name,
|
|
|
#p.county,
|
|
|
p.longitude,
|
|
|
p.latitude,
|
|
|
r.name,
|
|
|
c.name,
|
|
|
t,
|
|
|
round(value(z_prod[p, r, c, t]), digits=2),
|
|
|
round(value(z_plant_disp[p, r, c, t]), digits=2),
|
|
|
round(
|
|
|
sum(
|
|
|
value(y[p, q, r, c, t]) for (q, r2) in E_out[p] if r == r2;
|
|
|
init=0.0,
|
|
|
),
|
|
|
digits=2,
|
|
|
),
|
|
|
round(
|
|
|
data.c_plant_disp[p, r, t] * value(z_plant_disp[p, r, c, t]),
|
|
|
digits=2,
|
|
|
),
|
|
|
],
|
|
|
)
|
|
|
end
|
|
|
CSV.write("$output_dir/plant-outputs.csv", df)
|
|
|
|
|
|
# Report: Plant Emissions
|
|
|
# -------------------------------------------------------------------------
|
|
|
df = DataFrame()
|
|
|
df."plant" = String[]
|
|
|
df."emission" = String[]
|
|
|
df."time" = Int[]
|
|
|
df."amount emitted (tonne)" = Float64[]
|
|
|
df."emission cost (\$)" = Float64[]
|
|
|
for p in plants, s in emissions, t in T
|
|
|
push!(
|
|
|
df,
|
|
|
[
|
|
|
p.name,
|
|
|
s.name,
|
|
|
t,
|
|
|
round(value(z_plant_emissions[p, s, t]), digits=2),
|
|
|
round(
|
|
|
data.c_emission[s, t] * value(z_plant_emissions[p, s, t]),
|
|
|
digits=2,
|
|
|
),
|
|
|
],
|
|
|
)
|
|
|
end
|
|
|
CSV.write("$output_dir/plant-emissions.csv", df)
|
|
|
|
|
|
# Report: Transportation Emissions
|
|
|
# -------------------------------------------------------------------------
|
|
|
df = DataFrame()
|
|
|
df."source" = String[]
|
|
|
df."destination" = String[]
|
|
|
df."product" = String[]
|
|
|
df."emission" = String[]
|
|
|
df."time" = Int[]
|
|
|
df."distance (km)" = Float64[]
|
|
|
df."amount sent (tonne)" = Float64[]
|
|
|
df."amount emitted (tonne)" = Float64[]
|
|
|
df."emission cost (\$)" = Float64[]
|
|
|
for (q, p, r) in E, s in emissions, t in T
|
|
|
if value(y_total[q, p, r, t]) ≈ 0
|
|
|
continue
|
|
|
end
|
|
|
push!(
|
|
|
df,
|
|
|
[
|
|
|
q.name,
|
|
|
p.name,
|
|
|
r.name,
|
|
|
s.name,
|
|
|
t,
|
|
|
round(data.m_dist[q, p], digits=2),
|
|
|
round(value(y_total[q, p, r, t]), digits=2),
|
|
|
round(value(z_tr_emissions[q, p, r, s, t]), digits=2),
|
|
|
round(
|
|
|
data.c_emission[s, t] * value(z_tr_emissions[q, p, r, s, t]),
|
|
|
digits=2,
|
|
|
),
|
|
|
],
|
|
|
)
|
|
|
end
|
|
|
CSV.write("$output_dir/transp-emissions.csv", df)
|
|
|
|
|
|
return model, stats
|
|
|
end
|
