ANL-CEEESA
/
RELOG
mirror of https://github.com/ANL-CEEESA/RELOG.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Implement stochastic model

Browse Source
feature/stochastic
Alinson S. Xavier 3 years ago committed by Alinson S. Xavier
parent e915a57e58
commit 07ca3abb4f
Signed by: isoron
GPG Key ID: 0DA8E4B9E1109DCA
6 changed files with 427 additions and 332 deletions
Show Stats Download Patch File Download Diff File

8
Project.toml
Unescape View File

@ -6,13 +6,12 @@ version = "0.5.1"
[deps] [deps]
CRC = "44b605c4-b955-5f2b-9b6d-d2bd01d3d205" CRC = "44b605c4-b955-5f2b-9b6d-d2bd01d3d205"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63" GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d" Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d"
HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692" JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@ -24,16 +23,13 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568" ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StochasticPrograms = "8b8459f2-c380-502b-8633-9aed2d6c2b35"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
[compat] [compat]
CRC = "4" CRC = "4"
CSV = "0.7" CSV = "0.7"
Cbc = "0.6"
Clp = "0.8"
DataFrames = "0.21"
DataStructures = "0.17"
GZip = "0.5" GZip = "0.5"
Geodesy = "0.5" Geodesy = "0.5"
JSON = "0.21" JSON = "0.21"

4
src/graph/build.jl
Unescape View File

@ -59,7 +59,7 @@ function build_graph(instance::Instance)::Graph
dest.location.longitude, dest.location.longitude,
) )
values = Dict("distance" => distance) values = Dict("distance" => distance)
arc = Arc(source, dest, values) arc = Arc(length(arcs) + 1, source, dest, values)
push!(source.outgoing_arcs, arc) push!(source.outgoing_arcs, arc)
push!(dest.incoming_arcs, arc) push!(dest.incoming_arcs, arc)
push!(arcs, arc) push!(arcs, arc)
@ -72,7 +72,7 @@ function build_graph(instance::Instance)::Graph
for dest in shipping_nodes_by_plant[plant] for dest in shipping_nodes_by_plant[plant]
weight = plant.output[dest.product] weight = plant.output[dest.product]
values = Dict("weight" => weight) values = Dict("weight" => weight)
arc = Arc(source, dest, values) arc = Arc(length(arcs) + 1, source, dest, values)
push!(source.outgoing_arcs, arc) push!(source.outgoing_arcs, arc)
push!(dest.incoming_arcs, arc) push!(dest.incoming_arcs, arc)
push!(arcs, arc) push!(arcs, arc)

1
src/graph/structs.jl
Unescape View File

@ -7,6 +7,7 @@ using Geodesy
abstract type Node end abstract type Node end
mutable struct Arc mutable struct Arc
index::Int
source::Node source::Node
dest::Node dest::Node
values::Dict{String,Float64} values::Dict{String,Float64}

430
src/model/build.jl
Unescape View File

@ -2,278 +2,326 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures using JuMP, LinearAlgebra, Geodesy, ProgressBars, Printf, DataStructures, StochasticPrograms
function build_model(instance::Instance, graph::Graph, optimizer)::JuMP.Model function build_model(
model = Model(optimizer) instance::Instance,
model[:instance] = instance graphs::Vector{Graph},
model[:graph] = graph probs::Vector{Float64};
create_vars!(model) optimizer,
create_objective_function!(model) )
create_shipping_node_constraints!(model) T = instance.time
create_process_node_constraints!(model)
return model
end
@stochastic_model model begin
# Stage 1: Build plants
# =====================================================================
@stage 1 begin
pn = graphs[1].process_nodes
PN = length(pn)
function create_vars!(model::JuMP.Model) # Var: open_plant
graph, T = model[:graph], model[:instance].time @decision(
model[:flow] =
Dict((a, t) => @variable(model, lower_bound = 0) for a in graph.arcs, t = 1:T)
model[:plant_dispose] = Dict(
(n, t) => @variable(
model, model,
lower_bound = 0, open_plant[n in 1:PN, t in 1:T],
upper_bound = n.location.disposal_limit[n.product][t], binary = true,
) for n in values(graph.plant_shipping_nodes), t = 1:T
) )
model[:collection_dispose] = Dict(
(n, t) => @variable( # Var: is_open
@decision(
model, model,
lower_bound = 0, is_open[n in 1:PN, t in 1:T],
upper_bound = n.location.amount[t], binary = true,
) for
n in values(graph.collection_shipping_nodes), t = 1:T
)
model[:store] = Dict(
(n, t) =>
@variable(model, lower_bound = 0, upper_bound = n.location.storage_limit)
for n in values(graph.process_nodes), t = 1:T
)
model[:process] = Dict(
(n, t) => @variable(model, lower_bound = 0) for
n in values(graph.process_nodes), t = 1:T
) )
model[:open_plant] = Dict(
(n, t) => @variable(model, binary = true) for n in values(graph.process_nodes), # Objective function
t = 1:T @objective(
)
model[:is_open] = Dict(
(n, t) => @variable(model, binary = true) for n in values(graph.process_nodes),
t = 1:T
)
model[:capacity] = Dict(
(n, t) => @variable(
model, model,
lower_bound = 0, Min,
upper_bound = n.location.sizes[2].capacity
) for n in values(graph.process_nodes), t = 1:T # Opening, fixed operating costs
sum(
pn[n].location.sizes[1].opening_cost[t] * open_plant[n, t] +
pn[n].location.sizes[1].fixed_operating_cost[t] * is_open[n, t]
for n in 1:PN
for t in 1:T
),
) )
model[:expansion] = Dict(
(n, t) => @variable( for t = 1:T, n in 1:PN
# Plant is currently open if it was already open in the previous time period or
# if it was built just now
if t > 1
@constraint(
model, model,
lower_bound = 0, is_open[n, t] == is_open[n, t-1] + open_plant[n, t]
upper_bound = n.location.sizes[2].capacity - n.location.sizes[1].capacity
) for n in values(graph.process_nodes), t = 1:T
) )
end
function slope_open(plant, t)
if plant.sizes[2].capacity <= plant.sizes[1].capacity
0.0
else else
(plant.sizes[2].opening_cost[t] - plant.sizes[1].opening_cost[t]) / @constraint(model, is_open[n, t] == open_plant[n, t])
(plant.sizes[2].capacity - plant.sizes[1].capacity)
end
end end
function slope_fix_oper_cost(plant, t) # Plant can only be opened during building period
if plant.sizes[2].capacity <= plant.sizes[1].capacity if t ∉ instance.building_period
0.0 @constraint(model, open_plant[n, t] == 0)
else
(plant.sizes[2].fixed_operating_cost[t] - plant.sizes[1].fixed_operating_cost[t]) /
(plant.sizes[2].capacity - plant.sizes[1].capacity)
end end
end end
function create_objective_function!(model::JuMP.Model)
graph, T = model[:graph], model[:instance].time
obj = AffExpr(0.0)
# Process node costs
for n in values(graph.process_nodes), t = 1:T
# Transportation costs
for a in n.incoming_arcs
c = n.location.input.transportation_cost[t] * a.values["distance"]
add_to_expression!(obj, c, model[:flow][a, t])
end end
# Opening costs # Stage 2: Flows, disposal, capacity & storage
add_to_expression!( # =====================================================================
obj, @stage 2 begin
n.location.sizes[1].opening_cost[t], @uncertain graph
model[:open_plant][n, t], pn = graph.process_nodes
psn = graph.plant_shipping_nodes
csn = graph.collection_shipping_nodes
arcs = graph.arcs
A = length(arcs)
PN = length(pn)
CSN = length(csn)
PSN = length(psn)
# Var: flow
@recourse(
model,
flow[a in 1:A, t in 1:T],
lower_bound = 0,
) )
# Fixed operating costs (base) # Var: plant_dispose
add_to_expression!( @recourse(
obj, model,
n.location.sizes[1].fixed_operating_cost[t], plant_dispose[n in 1:PSN, t in 1:T],
model[:is_open][n, t], lower_bound = 0,
upper_bound = psn[n].location.disposal_limit[psn[n].product][t],
) )
# Fixed operating costs (expansion) # Var: collection_dispose/
add_to_expression!(obj, slope_fix_oper_cost(n.location, t), model[:expansion][n, t]) @recourse(
model,
# Processing costs collection_dispose[n in 1:CSN, t in 1:T],
add_to_expression!( lower_bound = 0,
obj, upper_bound = graph.collection_shipping_nodes[n].location.amount[t],
n.location.sizes[1].variable_operating_cost[t],
model[:process][n, t],
) )
# Storage costs # Var: store
add_to_expression!(obj, n.location.storage_cost[t], model[:store][n, t]) @recourse(
model,
# Expansion costs store[
if t < T n in 1:PN,
add_to_expression!( t in 1:T,
obj, ],
slope_open(n.location, t) - slope_open(n.location, t + 1), lower_bound = 0,
model[:expansion][n, t], upper_bound = pn[n].location.storage_limit,
) )
else
add_to_expression!(obj, slope_open(n.location, t), model[:expansion][n, t])
end
end
# Plant shipping node costs # Var: process
for n in values(graph.plant_shipping_nodes), t = 1:T @recourse(
# Disposal costs model,
add_to_expression!( process[
obj, n in 1:PN,
n.location.disposal_cost[n.product][t], t in 1:T,
model[:plant_dispose][n, t], ],
lower_bound = 0,
) )
end
# Collection shipping node costs # Var: capacity
for n in values(graph.collection_shipping_nodes), t = 1:T @recourse(
# Disposal costs model,
add_to_expression!( capacity[
obj, n in 1:PN,
n.location.product.disposal_cost[t], t in 1:T,
model[:collection_dispose][n, t], ],
lower_bound = 0,
upper_bound = pn[n].location.sizes[2].capacity,
) )
end
@objective(model, Min, obj)
end
function create_shipping_node_constraints!(model::JuMP.Model) # Var: expansion
graph, T = model[:graph], model[:instance].time @recourse(
model[:eq_balance] = OrderedDict()
for t = 1:T
# Collection centers
for n in graph.collection_shipping_nodes
model[:eq_balance][n, t] = @constraint(
model, model,
sum(model[:flow][a, t] for a in n.outgoing_arcs) == expansion[
n.location.amount[t] - model[:collection_dispose][n, t], n in 1:PN,
t in 1:T,
],
lower_bound = 0,
upper_bound = (
pn[n].location.sizes[2].capacity -
pn[n].location.sizes[1].capacity
),
) )
end
for prod in model[:instance].products
if isempty(prod.collection_centers)
continue
end
expr = AffExpr()
for center in prod.collection_centers
n = graph.collection_center_to_node[center]
add_to_expression!(expr, model[:collection_dispose][n, t])
end
@constraint(model, expr <= prod.disposal_limit[t])
end
# Plants # Objective function
for n in graph.plant_shipping_nodes @objective(
@constraint(
model, model,
sum(model[:flow][a, t] for a in n.incoming_arcs) == Min,
sum(model[:flow][a, t] for a in n.outgoing_arcs) + sum(
model[:plant_dispose][n, t] # Transportation costs
) pn[n].location.input.transportation_cost[t] *
end a.values["distance"] *
end flow[a.index,t]
for n in 1:PN
for a in pn[n].incoming_arcs
for t in 1:T
) + sum(
# Fixed operating costs (expansion)
slope_fix_oper_cost(pn[n].location, t) * expansion[n, t] +
end # Processing costs
pn[n].location.sizes[1].variable_operating_cost[t] * process[n, t] +
# Storage costs
pn[n].location.storage_cost[t] * store[n, t] +
# Expansion costs
(
t < T ? (
(
slope_open(pn[n].location, t) -
slope_open(pn[n].location, t + 1)
) * expansion[n, t]
) : slope_open(pn[n].location, t) * expansion[n, t]
)
function create_process_node_constraints!(model::JuMP.Model) for n in 1:PN
graph, T = model[:graph], model[:instance].time for t in 1:T
) + sum(
# Disposal costs (plants)
psn[n].location.disposal_cost[psn[n].product][t] * plant_dispose[n, t]
for n in 1:PSN
for t in 1:T
) + sum(
# Disposal costs (collection centers)
csn[n].location.product.disposal_cost[t] * collection_dispose[n, t]
for n in 1:CSN
for t in 1:T
)
)
for t = 1:T, n in graph.process_nodes # Process node constraints
input_sum = AffExpr(0.0) for t = 1:T, n in 1:PN
for a in n.incoming_arcs node = pn[n]
add_to_expression!(input_sum, 1.0, model[:flow][a, t])
end
# Output amount is implied by amount processed # Output amount is implied by amount processed
for a in n.outgoing_arcs for arc in node.outgoing_arcs
@constraint( @constraint(
model, model,
model[:flow][a, t] == a.values["weight"] * model[:process][n, t] flow[arc.index, t] == arc.values["weight"] * process[n, t]
) )
end end
# If plant is closed, capacity is zero # If plant is closed, capacity is zero
@constraint( @constraint(
model, model,
model[:capacity][n, t] <= n.location.sizes[2].capacity * model[:is_open][n, t] capacity[n, t] <= node.location.sizes[2].capacity * is_open[n, t]
) )
# If plant is open, capacity is greater than base # If plant is open, capacity is greater than base
@constraint( @constraint(
model, model,
model[:capacity][n, t] >= n.location.sizes[1].capacity * model[:is_open][n, t] capacity[n, t] >= node.location.sizes[1].capacity * is_open[n, t]
) )
# Capacity is linked to expansion # Capacity is linked to expansion
@constraint( @constraint(
model, model,
model[:capacity][n, t] <= capacity[n, t] <=
n.location.sizes[1].capacity + model[:expansion][n, t] node.location.sizes[1].capacity + expansion[n, t]
) )
# Can only process up to capacity # Can only process up to capacity
@constraint(model, model[:process][n, t] <= model[:capacity][n, t]) @constraint(model, process[n, t] <= capacity[n, t])
if t > 1 if t > 1
# Plant capacity can only increase over time # Plant capacity can only increase over time
@constraint(model, model[:capacity][n, t] >= model[:capacity][n, t-1]) @constraint(model, capacity[n, t] >= capacity[n, t-1])
@constraint(model, model[:expansion][n, t] >= model[:expansion][n, t-1]) @constraint(model, expansion[n, t] >= expansion[n, t-1])
end end
# Amount received equals amount processed plus stored # Amount received equals amount processed plus stored
store_in = 0 store_in = 0
if t > 1 if t > 1
store_in = model[:store][n, t-1] store_in = store[n, t-1]
end end
if t == T if t == T
@constraint(model, model[:store][n, t] == 0) @constraint(model, store[n, t] == 0)
end end
@constraint( @constraint(
model, model,
input_sum + store_in == model[:store][n, t] + model[:process][n, t] sum(
flow[arc.index, t]
for arc in node.incoming_arcs
) + store_in == store[n, t] + process[n, t]
) )
end
# Plant is currently open if it was already open in the previous time period or # Material flow at collection shipping nodes
# if it was built just now
if t > 1
@constraint( @constraint(
model, model,
model[:is_open][n, t] == model[:is_open][n, t-1] + model[:open_plant][n, t] eq_balance_centers[
n in 1:CSN,
t in 1:T,
],
sum(
flow[arc.index, t]
for arc in csn[n].outgoing_arcs
) == csn[n].location.amount[t] - collection_dispose[n, t]
) )
else
@constraint(model, model[:is_open][n, t] == model[:open_plant][n, t]) # Material flow at plant shipping nodes
@constraint(
model,
eq_balance_plant[
n in 1:PSN,
t in 1:T,
],
sum(flow[a.index, t] for a in psn[n].incoming_arcs) ==
sum(flow[a.index, t] for a in psn[n].outgoing_arcs) +
plant_dispose[n, t]
)
# Enforce product disposal limit at collection centers
for t in 1:T, prod in instance.products
if isempty(prod.collection_centers)
continue
end
@constraint(
model,
sum(
collection_dispose[n, t]
for n in 1:CSN
if csn[n].product == prod
) <= prod.disposal_limit[t]
)
end
end
end end
# Plant can only be opened during building period ξ = [
if t ∉ model[:instance].building_period @scenario graph = graphs[i] probability = probs[i]
@constraint(model, model[:open_plant][n, t] == 0) for i in 1:length(graphs)
]
sp = instantiate(model, ξ; optimizer)
return sp
end end
function slope_open(plant, t)
if plant.sizes[2].capacity <= plant.sizes[1].capacity
0.0
else
(plant.sizes[2].opening_cost[t] - plant.sizes[1].opening_cost[t]) /
(plant.sizes[2].capacity - plant.sizes[1].capacity)
end
end
function slope_fix_oper_cost(plant, t)
if plant.sizes[2].capacity <= plant.sizes[1].capacity
0.0
else
(plant.sizes[2].fixed_operating_cost[t] - plant.sizes[1].fixed_operating_cost[t]) /
(plant.sizes[2].capacity - plant.sizes[1].capacity)
end end
end end

110
src/model/getsol.jl
Unescape View File

@ -2,12 +2,31 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures using JuMP, LinearAlgebra, Geodesy, ProgressBars, Printf, DataStructures
function get_solution(
instance,
graph,
model,
scenario_index::Int;
marginal_costs=false,
)
value(x) = StochasticPrograms.value(x, scenario_index)
shadow_price(x) = StochasticPrograms.shadow_price(x, scenario_index)
function get_solution(model::JuMP.Model; marginal_costs = true)
graph, instance = model[:graph], model[:instance]
T = instance.time T = instance.time
pn = graph.process_nodes
psn = graph.plant_shipping_nodes
csn = graph.collection_shipping_nodes
arcs = graph.arcs
A = length(arcs)
PN = length(pn)
CSN = length(csn)
PSN = length(psn)
output = OrderedDict( output = OrderedDict(
"Plants" => OrderedDict(), "Plants" => OrderedDict(),
"Products" => OrderedDict(), "Products" => OrderedDict(),
@ -29,42 +48,52 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
), ),
) )
plant_to_process_node = OrderedDict(n.location => n for n in graph.process_nodes) pn = graph.process_nodes
plant_to_shipping_nodes = OrderedDict() psn = graph.plant_shipping_nodes
for p in instance.plants
plant_to_shipping_nodes[p] = [] plant_to_process_node_index = OrderedDict(
for a in plant_to_process_node[p].outgoing_arcs pn[n].location => n
push!(plant_to_shipping_nodes[p], a.dest) for n in 1:length(pn)
end )
plant_to_shipping_node_indices = OrderedDict(p => [] for p in instance.plants)
for n in 1:length(psn)
push!(plant_to_shipping_node_indices[psn[n].location], n)
end end
# Products # Products
for n in graph.collection_shipping_nodes for n in 1:CSN
node = csn[n]
location_dict = OrderedDict{Any,Any}( location_dict = OrderedDict{Any,Any}(
"Latitude (deg)" => n.location.latitude, "Latitude (deg)" => node.location.latitude,
"Longitude (deg)" => n.location.longitude, "Longitude (deg)" => node.location.longitude,
"Amount (tonne)" => n.location.amount, "Amount (tonne)" => node.location.amount,
"Dispose (tonne)" => "Dispose (tonne)" => [
[JuMP.value(model[:collection_dispose][n, t]) for t = 1:T], value(model[2, :collection_dispose][n, t])
"Disposal cost (\$)" => for t = 1:T
[JuMP.value(model[:collection_dispose][n, t]) * n.location.product.disposal_cost[t] for t = 1:T] ],
"Disposal cost (\$)" => [
value(model[2, :collection_dispose][n, t]) *
node.location.product.disposal_cost[t]
for t = 1:T
]
) )
if marginal_costs if marginal_costs
location_dict["Marginal cost (\$/tonne)"] = [ location_dict["Marginal cost (\$/tonne)"] = [
round(abs(JuMP.shadow_price(model[:eq_balance][n, t])), digits = 2) for round(abs(shadow_price(model[2, :eq_balance_centers][n, t])), digits=2) for t = 1:T
t = 1:T
] ]
end end
if n.product.name ∉ keys(output["Products"]) if node.product.name ∉ keys(output["Products"])
output["Products"][n.product.name] = OrderedDict() output["Products"][node.product.name] = OrderedDict()
end end
output["Products"][n.product.name][n.location.name] = location_dict output["Products"][node.product.name][node.location.name] = location_dict
end end
# Plants # Plants
for plant in instance.plants for plant in instance.plants
skip_plant = true skip_plant = true
process_node = plant_to_process_node[plant] n = plant_to_process_node_index[plant]
process_node = pn[n]
plant_dict = OrderedDict{Any,Any}( plant_dict = OrderedDict{Any,Any}(
"Input" => OrderedDict(), "Input" => OrderedDict(),
"Output" => "Output" =>
@ -75,39 +104,39 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
"Latitude (deg)" => plant.latitude, "Latitude (deg)" => plant.latitude,
"Longitude (deg)" => plant.longitude, "Longitude (deg)" => plant.longitude,
"Capacity (tonne)" => "Capacity (tonne)" =>
[JuMP.value(model[:capacity][process_node, t]) for t = 1:T], [value(model[2, :capacity][n, t]) for t = 1:T],
"Opening cost (\$)" => [ "Opening cost (\$)" => [
JuMP.value(model[:open_plant][process_node, t]) * value(model[2, :open_plant][n, t]) *
plant.sizes[1].opening_cost[t] for t = 1:T plant.sizes[1].opening_cost[t] for t = 1:T
], ],
"Fixed operating cost (\$)" => [ "Fixed operating cost (\$)" => [
JuMP.value(model[:is_open][process_node, t]) * value(model[2, :is_open][n, t]) *
plant.sizes[1].fixed_operating_cost[t] + plant.sizes[1].fixed_operating_cost[t] +
JuMP.value(model[:expansion][process_node, t]) * value(model[2, :expansion][n, t]) *
slope_fix_oper_cost(plant, t) for t = 1:T slope_fix_oper_cost(plant, t) for t = 1:T
], ],
"Expansion cost (\$)" => [ "Expansion cost (\$)" => [
( (
if t == 1 if t == 1
slope_open(plant, t) * JuMP.value(model[:expansion][process_node, t]) slope_open(plant, t) * value(model[2, :expansion][n, t])
else else
slope_open(plant, t) * ( slope_open(plant, t) * (
JuMP.value(model[:expansion][process_node, t]) - value(model[2, :expansion][n, t]) -
JuMP.value(model[:expansion][process_node, t-1]) value(model[2, :expansion][n, t-1])
) )
end end
) for t = 1:T ) for t = 1:T
], ],
"Process (tonne)" => "Process (tonne)" =>
[JuMP.value(model[:process][process_node, t]) for t = 1:T], [value(model[2, :process][n, t]) for t = 1:T],
"Variable operating cost (\$)" => [ "Variable operating cost (\$)" => [
JuMP.value(model[:process][process_node, t]) * value(model[2, :process][n, t]) *
plant.sizes[1].variable_operating_cost[t] for t = 1:T plant.sizes[1].variable_operating_cost[t] for t = 1:T
], ],
"Storage (tonne)" => "Storage (tonne)" =>
[JuMP.value(model[:store][process_node, t]) for t = 1:T], [value(model[2, :store][n, t]) for t = 1:T],
"Storage cost (\$)" => [ "Storage cost (\$)" => [
JuMP.value(model[:store][process_node, t]) * plant.storage_cost[t] value(model[2, :store][n, t]) * plant.storage_cost[t]
for t = 1:T for t = 1:T
], ],
) )
@ -120,7 +149,7 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
# Inputs # Inputs
for a in process_node.incoming_arcs for a in process_node.incoming_arcs
vals = [JuMP.value(model[:flow][a, t]) for t = 1:T] vals = [value(model[2, :flow][a.index, t]) for t = 1:T]
if sum(vals) <= 1e-3 if sum(vals) <= 1e-3
continue continue
end end
@ -178,19 +207,20 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
end end
# Outputs # Outputs
for shipping_node in plant_to_shipping_nodes[plant] for n2 in plant_to_shipping_node_indices[plant]
shipping_node = psn[n2]
product_name = shipping_node.product.name product_name = shipping_node.product.name
plant_dict["Total output"][product_name] = zeros(T) plant_dict["Total output"][product_name] = zeros(T)
plant_dict["Output"]["Send"][product_name] = product_dict = OrderedDict() plant_dict["Output"]["Send"][product_name] = product_dict = OrderedDict()
disposal_amount = disposal_amount =
[JuMP.value(model[:plant_dispose][shipping_node, t]) for t = 1:T] [value(model[2, :plant_dispose][n2, t]) for t = 1:T]
if sum(disposal_amount) > 1e-5 if sum(disposal_amount) > 1e-5
skip_plant = false skip_plant = false
plant_dict["Output"]["Dispose"][product_name] = plant_dict["Output"]["Dispose"][product_name] =
disposal_dict = OrderedDict() disposal_dict = OrderedDict()
disposal_dict["Amount (tonne)"] = disposal_dict["Amount (tonne)"] =
[JuMP.value(model[:plant_dispose][shipping_node, t]) for t = 1:T] [value(model[2, :plant_dispose][n2, t]) for t = 1:T]
disposal_dict["Cost (\$)"] = [ disposal_dict["Cost (\$)"] = [
disposal_dict["Amount (tonne)"][t] * disposal_dict["Amount (tonne)"][t] *
plant.disposal_cost[shipping_node.product][t] for t = 1:T plant.disposal_cost[shipping_node.product][t] for t = 1:T
@ -200,7 +230,7 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
end end
for a in shipping_node.outgoing_arcs for a in shipping_node.outgoing_arcs
vals = [JuMP.value(model[:flow][a, t]) for t = 1:T] vals = [value(model[2, :flow][a.index, t]) for t = 1:T]
if sum(vals) <= 1e-3 if sum(vals) <= 1e-3
continue continue
end end

90
src/model/solve.jl
Unescape View File

@ -2,14 +2,14 @@
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved. # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details. # Released under the modified BSD license. See COPYING.md for more details.
using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures using JuMP, LinearAlgebra, Geodesy, HiGHS, ProgressBars, Printf, DataStructures
function _get_default_milp_optimizer() function _get_default_milp_optimizer()
return optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0) return optimizer_with_attributes(HiGHS.Optimizer)
end end
function _get_default_lp_optimizer() function _get_default_lp_optimizer()
return optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0) return optimizer_with_attributes(HiGHS.Optimizer)
end end
@ -25,61 +25,81 @@ function _print_graph_stats(instance::Instance, graph::Graph)::Nothing
return return
end end
function solve_stochastic(;
scenarios::Vector{String},
probs::Vector{Float64},
optimizer,
)
@info "Reading instance files..."
instances = [parsefile(sc) for sc in scenarios]
@info "Building graphs..."
graphs = [build_graph(inst) for inst in instances]
@info "Building stochastic model..."
sp = RELOG.build_model(instances[1], graphs, probs; optimizer)
@info "Optimizing stochastic model..."
optimize!(sp)
@info "Extracting solution..."
solutions = [
get_solution(instances[i], graphs[i], sp, i)
for i in 1:length(instances)
]
return solutions
end
function solve( function solve(
instance::Instance; instance::Instance;
optimizer=nothing, optimizer=nothing,
output = nothing,
graph = nothing,
marginal_costs=true, marginal_costs=true,
return_model = false, return_model=false
) )
milp_optimizer = lp_optimizer = optimizer
if optimizer === nothing
milp_optimizer = _get_default_milp_optimizer()
lp_optimizer = _get_default_lp_optimizer()
end
if graph === nothing
@info "Building graph..." @info "Building graph..."
graph = RELOG.build_graph(instance) graph = RELOG.build_graph(instance)
_print_graph_stats(instance, graph) _print_graph_stats(instance, graph)
end
@info "Building optimization model..." @info "Building model..."
model = RELOG.build_model(instance, graph, milp_optimizer) model = RELOG.build_model(instance, [graph], [1.0]; optimizer)
@info "Optimizing MILP..." @info "Optimizing model..."
JuMP.optimize!(model) optimize!(model)
if !has_values(model) if !has_values(model)
error("No solution available") error("No solution available")
end end
solution = get_solution(model, marginal_costs = false)
@info "Extracting solution..."
solution = get_solution(instance, graph, model, 1)
if marginal_costs if marginal_costs
@info "Re-optimizing with integer variables fixed..." @info "Re-optimizing with integer variables fixed..."
all_vars = JuMP.all_variables(model) open_plant_vals = value.(model[1, :open_plant])
vals = OrderedDict(var => JuMP.value(var) for var in all_vars) is_open_vals = value.(model[1, :is_open])
JuMP.set_optimizer(model, lp_optimizer)
for var in all_vars for n in 1:length(graph.process_nodes), t in 1:instance.time
if JuMP.is_binary(var) unset_binary(model[1, :open_plant][n, t])
JuMP.unset_binary(var) unset_binary(model[1, :is_open][n, t])
JuMP.fix(var, vals[var]) fix(
end model[1, :open_plant][n, t],
open_plant_vals[n, t]
)
fix(
model[1, :is_open][n, t],
is_open_vals[n, t]
)
end end
JuMP.optimize!(model) optimize!(model)
if has_values(model) if has_values(model)
@info "Extracting solution..." @info "Extracting solution..."
solution = get_solution(model, marginal_costs = true) solution = get_solution(instance, graph, model, 1, marginal_costs=true)
else else
@warn "Error computing marginal costs. Ignoring." @warn "Error computing marginal costs. Ignoring."
end end
end end
if output !== nothing
write(solution, output)
end
if return_model if return_model
return solution, model return solution, model
else else
@ -93,7 +113,7 @@ function solve(filename::AbstractString; heuristic = false, kwargs...)
if heuristic && instance.time > 1 if heuristic && instance.time > 1
@info "Solving single-period version..." @info "Solving single-period version..."
compressed = _compress(instance) compressed = _compress(instance)
csol, model = solve(compressed; output = nothing, marginal_costs = false, return_model = true, kwargs...) csol, model = solve(compressed; marginal_costs=false, return_model=true, kwargs...)
@info "Filtering candidate locations..." @info "Filtering candidate locations..."
selected_pairs = [] selected_pairs = []
for (plant_name, plant_dict) in csol["Plants"] for (plant_name, plant_dict) in csol["Plants"]

Write Preview
Loading…
Cancel
Save