diff --git a/Manifest.toml b/Manifest.toml index eba600e..f1fd005 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.1" +julia_version = "1.10.2" manifest_format = "2.0" -project_hash = "4170463ca454c6921d0f3b304e50e00eea5e30ba" +project_hash = "c30276d4b440bf43a36bc652205ac642814efc55" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" @@ -31,6 +31,12 @@ git-tree-sha1 = "b7ba7f0d727433c961909b329c4d2263268da4c9" uuid = "44b605c4-b955-5f2b-9b6d-d2bd01d3d205" version = "4.0.0" +[[deps.CSTParser]] +deps = ["Tokenize"] +git-tree-sha1 = "b544d62417a99d091c569b95109bc9d8c223e9e3" +uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f" +version = "3.4.2" + [[deps.CSV]] deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "PrecompileTools", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings", "WorkerUtilities"] git-tree-sha1 = "679e69c611fff422038e9e21e270c4197d49d918" @@ -55,6 +61,12 @@ git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" version = "0.7.4" +[[deps.CommonMark]] +deps = ["Crayons", "JSON", "PrecompileTools", "URIs"] +git-tree-sha1 = "532c4185d3c9037c0237546d817858b23cf9e071" +uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6" +version = "0.8.12" + [[deps.CommonSubexpressions]] deps = ["MacroTools", "Test"] git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7" @@ -154,6 +166,11 @@ deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" version = "1.6.0" +[[deps.ExprTools]] +git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec" +uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04" +version = "0.1.10" + [[deps.FilePathsBase]] deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"] git-tree-sha1 = "9f00e42f8d99fdde64d40c8ea5d14269a2e2c1aa" @@ -183,6 +200,11 @@ git-tree-sha1 = "ed98a4429bf0a033ccc5e036120181dd52f06d31" uuid = "0ef565a4-170c-5f04-8de2-149903a85f3d" version = "1.1.0" +[[deps.Glob]] +git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496" +uuid = "c27321d9-0574-5035-807b-f59d2c89b15c" +version = "1.3.1" + [[deps.Gurobi]] deps = ["LazyArtifacts", "Libdl", "MathOptInterface"] git-tree-sha1 = "5995b72d385235f3fe55f8f0c4ad61049f867814" @@ -250,6 +272,12 @@ version = "1.20.0" [deps.JuMP.weakdeps] DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" +[[deps.JuliaFormatter]] +deps = ["CSTParser", "CommonMark", "DataStructures", "Glob", "Pkg", "PrecompileTools", "Tokenize"] +git-tree-sha1 = "1c4880cb70a5c6c87ea36deccc3d7f9e7969c18c" +uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899" +version = "1.0.56" + [[deps.JuliaInterpreter]] deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] git-tree-sha1 = "7b762d81887160169ddfc93a47e5fd7a6a3e78ef" @@ -572,6 +600,17 @@ version = "1.10.0" deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[deps.TimerOutputs]] +deps = ["ExprTools", "Printf"] +git-tree-sha1 = "f548a9e9c490030e545f72074a41edfd0e5bcdd7" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.5.23" + +[[deps.Tokenize]] +git-tree-sha1 = "5b5a892ba7704c0977013bd0f9c30f5d962181e0" +uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624" +version = "0.5.28" + [[deps.TranscodingStreams]] git-tree-sha1 = "54194d92959d8ebaa8e26227dbe3cdefcdcd594f" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" @@ -581,6 +620,11 @@ weakdeps = ["Random", "Test"] [deps.TranscodingStreams.extensions] TestExt = ["Test", "Random"] +[[deps.URIs]] +git-tree-sha1 = "67db6cc7b3821e19ebe75791a9dd19c9b1188f2b" +uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" +version = "1.5.1" + [[deps.UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" diff --git a/Project.toml b/Project.toml index 0b343e9..5dbeee9 100644 --- a/Project.toml +++ b/Project.toml @@ -7,8 +7,10 @@ Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" diff --git a/dist.jl b/dist.jl index 1008650..5e8b016 100644 --- a/dist.jl +++ b/dist.jl @@ -8,6 +8,7 @@ using DataFrames using CRC using ZipFile using Statistics +using TimerOutputs crc32 = crc(CRC_32) @@ -77,31 +78,37 @@ function _calculate_distance( csv_filename = joinpath(basedir, "dist_driving.csv") # Download pre-computed driving data - if !isfile(csv_filename) - _download_zip( - "https://axavier.org/RELOG/0.6/data/dist_driving_0b9a6ad6.zip", - basedir, - csv_filename, - 0x0b9a6ad6, - ) + @timeit "Download data" begin + if !isfile(csv_filename) + _download_zip( + "https://axavier.org/RELOG/0.6/data/dist_driving_0b9a6ad6.zip", + basedir, + csv_filename, + 0x0b9a6ad6, + ) + end end - # Fit kNN model - df = DataFrame(CSV.File(csv_filename, missingstring="NaN")) - dropmissing!(df) - coords = Matrix(df[!, [:source_lat, :source_lon, :dest_lat, :dest_lon]])' - metric.ratios = Matrix(df[!, [:ratio]]) - metric.tree = KDTree(coords) + @timeit "Fit KNN model" begin + df = DataFrame(CSV.File(csv_filename, missingstring="NaN")) + dropmissing!(df) + coords = Matrix(df[!, [:source_lat, :source_lon, :dest_lat, :dest_lon]])' + metric.ratios = Matrix(df[!, [:ratio]]) + metric.tree = KDTree(coords) + end end - # Compute Euclidean distance - dist_euclidean = - _calculate_distance(source_lat, source_lon, dest_lat, dest_lon, EuclideanDistance()) + @timeit "Compute Euclidean distance" begin + dist_euclidean = + _calculate_distance(source_lat, source_lon, dest_lat, dest_lon, EuclideanDistance()) + end + + @timeit "Predict driving distance" begin + idxs, _ = knn(metric.tree, [source_lat, source_lon, dest_lat, dest_lon], 5) + ratio_pred = mean(metric.ratios[idxs]) + dist_pred = round(dist_euclidean * ratio_pred, digits=3) + isfinite(dist_pred) || error("non-finite distance detected: $dist_pred") + end - # Predict ratio - idxs, _ = knn(metric.tree, [source_lat, source_lon, dest_lat, dest_lon], 5) - ratio_pred = mean(metric.ratios[idxs]) - dist_pred = round(dist_euclidean * ratio_pred, digits=3) - isfinite(dist_pred) || error("non-finite distance detected: $dist_pred") return dist_pred end diff --git a/jumpext.jl b/jumpext.jl index 5ac3ac4..9247393 100644 --- a/jumpext.jl +++ b/jumpext.jl @@ -27,11 +27,9 @@ function _init(model::JuMP.Model, key::Symbol)::OrderedDict end function _set_names!(model::JuMP.Model) - @info "Setting variable and constraint names..." time_varnames = @elapsed begin _set_names!(object_dictionary(model)) end - @info @sprintf("Set names in %.2f seconds", time_varnames) end function _set_names!(dict::Dict) diff --git a/model-3.jl b/model-3.jl index 15d8106..93dc044 100644 --- a/model-3.jl +++ b/model-3.jl @@ -8,6 +8,8 @@ using JSON using JuMP using OrderedCollections using Random +using TimerOutputs +using Printf include("jumpext.jl") include("dist.jl") @@ -47,12 +49,9 @@ Base.@kwdef struct Emission name::String end -Base.show( - io::IO, - p::Union{Component,Product,Center,Plant,Emission}, -) = print(io, p.name) +Base.show(io::IO, p::Union{Component,Product,Center,Plant,Emission}) = print(io, p.name) -Base.@kwdef struct Instance +Base.@kwdef mutable struct Instance T::UnitRange{Int} centers::Vector{Center} plants::Vector{Plant} @@ -113,69 +112,49 @@ function generate_data() 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], - ) + 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 + 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 + 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 + 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"), - ] + 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 + (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 @@ -187,63 +166,22 @@ function generate_data() 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 - ) + 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 + (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( @@ -252,38 +190,20 @@ function generate_data() 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 + ) 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 - ) + 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 - ) + 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, @@ -309,7 +229,7 @@ function generate_data() m_emission, m_init, m_plant_disp, - m_store + m_store, ) end @@ -318,29 +238,17 @@ end function write_json(data, filename) json = dict() - json["parameters"] = dict( - "time horizon (years)" => data.T.stop, - ) + 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 - ], + "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 + 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( @@ -349,30 +257,18 @@ function write_json(data, filename) "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 + data.m_init[q, r, c, t] for c in r.comp, 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 + "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( @@ -382,55 +278,30 @@ function write_json(data, filename) "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 + data.alpha_mix[p, r, c_in, c_out] for + c_in in p.prod_in.comp, c_out in r.comp ], - ) - for r in p.prod_out + "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 ), - "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 - ], + "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 + 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 + "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) @@ -439,7 +310,7 @@ end # Read # ============================================================================== -function read_json(filename) +function read_json(filename, max_centers = 10, max_plants = 10) json = JSON.parsefile(filename) T = 1:json["parameters"]["time horizon (years)"] centers = [] @@ -468,99 +339,122 @@ function read_json(filename) m_plant_disp = dict() m_store = dict() - 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] + @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 - for (name, prod_data) in json["products"] - comp = [] - for (comp_name) in prod_data["components"] - if comp_name ∉ keys(components_by_name) - components_by_name[comp_name] = Component(comp_name) + @timeit "Read: Products" begin + for (name, prod_data) in json["products"] + comp = [] + for (comp_name) in prod_data["components"] + if comp_name ∉ keys(components_by_name) + components_by_name[comp_name] = Component(comp_name) + end + push!(comp, components_by_name[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] + r = Product(name, comp) + products_by_name[name] = r + push!(products, r) for t in T - alpha_tr_emission[r, s, t] = s_data[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 - for (name, center_data) in json["centers"] - 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_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] + @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_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 - for (plant_name, plant_data) in json["plants"] - 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)"] - 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] + @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 + 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)"] for t in T - alpha_plant_emission[p, s, t] = s_data[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 - metric = KnnDrivingDistance() - m_dist = dict( - (p, q) => _calculate_distance( - p.latitude, - p.longitude, - q.latitude, - q.longitude, - metric, + @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 ) - for p in [plants; centers], q in plants - ) + end return Instance(; T, @@ -586,7 +480,7 @@ function read_json(filename) m_emission, m_init, m_plant_disp, - m_store + m_store, ) end @@ -603,8 +497,11 @@ function generate_json() end function solve(filename, optimizer) - @info "Reading JSON file" - data = read_json(filename) + reset_timer!() + + @timeit "Read JSON" begin + data = read_json(filename) + end T = data.T centers = data.centers @@ -616,397 +513,415 @@ function solve(filename, optimizer) # Graph # ------------------------------------------------------------------------- - E = [] - E_in = dict(src => [] for src in plants) - E_out = dict(src => [] for src in plants ∪ centers) - function push_edge!(src, dst, r) - push!(E, (src, dst, r)) - 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 + @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) + push!(E, (src, dst, r)) + push!(E_out[src], (dst, r)) + push!(E_in[dst], (src, r)) 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) + 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 + @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)) + # Decision variables # ------------------------------------------------------------------------- - 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 - - 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 - - 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 - - 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 - - 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) + @timeit "Model: Add variables" begin + @timeit "y" begin + y = _init(model, :y) + for (q, p, r) in E + y[q, p, r] = @variable(model, [r.comp, T], lower_bound = 0) + end end - end - - 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) + @timeit "y_total" begin + y_total = _init(model, :y_total) + for (q, p, r) in E + y_total[q, p, r] = @variable(model, [T], lower_bound = 0) end end - end - - 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) + @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_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 end - 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 - - 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 - - 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 - - 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 - - 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 - - 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 # Objective function # ------------------------------------------------------------------------- - obj = AffExpr() - - # Center disposal - 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 + @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 - 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] - ) - 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]) + end + end - # Center storage - 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 + # 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 - for (q, p, r) in E, c in r.comp, t in T - add_to_expression!( - obj, - data.m_dist[q, p] * data.c_tr[r, t], - y[q, p, r, c, t] - ) - 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_qpr[c, t]) + add_to_expression!(obj, data.c_var[p, t], y_qpr[c, t]) + end + end + end - # Transportation emissions - 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 + # 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 - for p in plants, t in T - add_to_expression!( - obj, - data.c_fix[p, t], - x_open[p, t] - ) - 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 - 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 + # 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 - # Variable operating cost - for (q, p, r) in E, c in r.comp, t in T - add_to_expression!( - obj, - data.c_var[p, t], - y[q, p, r, c, t] - ) - 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 disposal - 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 + # 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 - # Plant emissions - 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] - ) + @objective(model, Min, obj) end - @objective(model, Min, obj) - # Constraints # ------------------------------------------------------------------------- - 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(data.m_init[q, r, c, t] for c in r.comp) + - z_store_total[q, r, t-1] - ) - 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) - ) - - 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 - - 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 + @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(data.m_init[q, r, c, t] for c in r.comp) + z_store_total[q, r, t-1] + ) + end 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 - - 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 - - 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 - - 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 + 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 - 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 - - 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 - - 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 - - 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 - - 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 - ) - ) - end - - 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 - - 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 - - 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 - - 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 - - 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 - - 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 - - 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 # Optimize # ------------------------------------------------------------------------- - _set_names!(model) optimize!(model) # Report: Transportation @@ -1037,8 +952,8 @@ function solve(filename, optimizer) data.m_dist[q, p], value(y[q, p, r, c, t]), data.m_dist[q, p] * data.c_tr[r, t] * value(y[q, p, r, c, t]), - data.c_var[p, t] * value(y[q, p, r, c, t]) - ] + data.c_var[p, t] * value(y[q, p, r, c, t]), + ], ) end CSV.write("$output_dir/transp.csv", df) @@ -1066,17 +981,13 @@ function solve(filename, optimizer) c.name, t, data.m_init[q, r, c, t], - sum( - value(y[q, p, r, c, t]) - for (p, r2) in E_out[q] - if r == r2 - ), + sum(value(y[q, p, r, c, t]) for (p, r2) in E_out[q] if r == r2), value(z_store[q, r, c, t]), value(z_center_disp[q, r, c, t]), data.m_init[q, r, c, t] * data.c_acq[q, r, t], data.c_store[q, r, t] * value(z_store[q, r, c, t]), data.c_center_disp[q, r, t] * value(z_center_disp[q, r, c, t]), - ] + ], ) end CSV.write("$output_dir/centers.csv", df) @@ -1098,7 +1009,7 @@ function solve(filename, optimizer) value(x_open[p, t]), data.c_open[p, t] * (value(x_open[p, t]) - value(x_open[p, t-1])), data.c_fix[p, t] * value(x_open[p, t]), - ] + ], ) end CSV.write("$output_dir/plants.csv", df) @@ -1127,14 +1038,9 @@ function solve(filename, optimizer) t, value(z_prod[p, r, c, t]), value(z_plant_disp[p, r, c, t]), - sum( - value(y[p, q, r, c, t]) - for (q, r2) in E_out[p] - if r == r2; - init=0.0 - ), + sum(value(y[p, q, r, c, t]) for (q, r2) in E_out[p] if r == r2; init = 0.0), data.c_plant_disp[p, r, t] * value(z_plant_disp[p, r, c, t]), - ] + ], ) end CSV.write("$output_dir/plant-outputs.csv", df) @@ -1156,7 +1062,7 @@ function solve(filename, optimizer) t, value(z_plant_emissions[p, s, t]), data.c_emission[s, t] * value(z_plant_emissions[p, s, t]), - ] + ], ) end CSV.write("$output_dir/plant-emissions.csv", df) @@ -1189,11 +1095,10 @@ function solve(filename, optimizer) value(y_total[q, p, r, t]), value(z_tr_emissions[q, p, r, s, t]), data.c_emission[s, t] * value(z_tr_emissions[q, p, r, s, t]), - ] + ], ) end CSV.write("$output_dir/transp-emissions.csv", df) return end -