diff --git a/Project.toml b/Project.toml
index f4b5119..b03521d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,13 +6,12 @@ version = "0.5.1"
 [deps]
 CRC = "44b605c4-b955-5f2b-9b6d-d2bd01d3d205"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
-Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GZip = "92fee26a-97fe-5a0c-ad85-20a5f3185b63"
 Geodesy = "0ef565a4-170c-5f04-8de2-149903a85f3d"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -24,16 +23,13 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StochasticPrograms = "8b8459f2-c380-502b-8633-9aed2d6c2b35"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [compat]
 CRC = "4"
 CSV = "0.7"
-Cbc = "0.6"
-Clp = "0.8"
-DataFrames = "0.21"
-DataStructures = "0.17"
 GZip = "0.5"
 Geodesy = "0.5"
 JSON = "0.21"
diff --git a/src/graph/build.jl b/src/graph/build.jl
index 27a323e..ac9e1a9 100644
--- a/src/graph/build.jl
+++ b/src/graph/build.jl
@@ -59,7 +59,7 @@ function build_graph(instance::Instance)::Graph
                 dest.location.longitude,
             )
             values = Dict("distance" => distance)
-            arc = Arc(source, dest, values)
+            arc = Arc(length(arcs) + 1, source, dest, values)
             push!(source.outgoing_arcs, arc)
             push!(dest.incoming_arcs, arc)
             push!(arcs, arc)
@@ -72,7 +72,7 @@ function build_graph(instance::Instance)::Graph
         for dest in shipping_nodes_by_plant[plant]
             weight = plant.output[dest.product]
             values = Dict("weight" => weight)
-            arc = Arc(source, dest, values)
+            arc = Arc(length(arcs) + 1, source, dest, values)
             push!(source.outgoing_arcs, arc)
             push!(dest.incoming_arcs, arc)
             push!(arcs, arc)
diff --git a/src/graph/structs.jl b/src/graph/structs.jl
index e4f1149..840fe33 100644
--- a/src/graph/structs.jl
+++ b/src/graph/structs.jl
@@ -7,6 +7,7 @@ using Geodesy
 abstract type Node end
 
 mutable struct Arc
+    index::Int
     source::Node
     dest::Node
     values::Dict{String,Float64}
diff --git a/src/model/build.jl b/src/model/build.jl
index 1acda5f..9d8f3db 100644
--- a/src/model/build.jl
+++ b/src/model/build.jl
@@ -2,278 +2,326 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # Released under the modified BSD license. See COPYING.md for more details.
 
-using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures
-
-function build_model(instance::Instance, graph::Graph, optimizer)::JuMP.Model
-    model = Model(optimizer)
-    model[:instance] = instance
-    model[:graph] = graph
-    create_vars!(model)
-    create_objective_function!(model)
-    create_shipping_node_constraints!(model)
-    create_process_node_constraints!(model)
-    return model
-end
-
-
-function create_vars!(model::JuMP.Model)
-    graph, T = model[:graph], model[:instance].time
-    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,
-            lower_bound = 0,
-            upper_bound = n.location.disposal_limit[n.product][t],
-        ) for n in values(graph.plant_shipping_nodes), t = 1:T
-    )
-    model[:collection_dispose] = Dict(
-        (n, t) => @variable(
-            model,
-            lower_bound = 0,
-            upper_bound = n.location.amount[t],
-        ) 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),
-        t = 1:T
-    )
-    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,
-            lower_bound = 0,
-            upper_bound = n.location.sizes[2].capacity
-        ) for n in values(graph.process_nodes), t = 1:T
-    )
-    model[:expansion] = Dict(
-        (n, t) => @variable(
-            model,
-            lower_bound = 0,
-            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
-        (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
+using JuMP, LinearAlgebra, Geodesy, ProgressBars, Printf, DataStructures, StochasticPrograms
+
+function build_model(
+    instance::Instance,
+    graphs::Vector{Graph},
+    probs::Vector{Float64};
+    optimizer,
+)
+    T = instance.time
+
+    @stochastic_model model begin
+        # Stage 1: Build plants
+        # =====================================================================
+        @stage 1 begin
+            pn = graphs[1].process_nodes
+            PN = length(pn)
+
+            # Var: open_plant
+            @decision(
+                model,
+                open_plant[n in 1:PN, t in 1:T],
+                binary = true,
+            )
 
-function create_objective_function!(model::JuMP.Model)
-    graph, T = model[:graph], model[:instance].time
-    obj = AffExpr(0.0)
+            # Var: is_open
+            @decision(
+                model,
+                is_open[n in 1:PN, t in 1:T],
+                binary = true,
+            )
 
-    # Process node costs
-    for n in values(graph.process_nodes), t = 1:T
+            # Objective function
+            @objective(
+                model,
+                Min,
+
+                # 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
+                ),
+            )
 
-        # 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])
+            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,
+                        is_open[n, t] == is_open[n, t-1] + open_plant[n, t]
+                    )
+                else
+                    @constraint(model, is_open[n, t] == open_plant[n, t])
+                end
+        
+                # Plant can only be opened during building period
+                if t ∉ instance.building_period
+                    @constraint(model, open_plant[n, t] == 0)
+                end
+            end
         end
 
-        # Opening costs
-        add_to_expression!(
-            obj,
-            n.location.sizes[1].opening_cost[t],
-            model[:open_plant][n, t],
-        )
-
-        # Fixed operating costs (base)
-        add_to_expression!(
-            obj,
-            n.location.sizes[1].fixed_operating_cost[t],
-            model[:is_open][n, t],
-        )
-
-        # Fixed operating costs (expansion)
-        add_to_expression!(obj, slope_fix_oper_cost(n.location, t), model[:expansion][n, t])
-
-        # Processing costs
-        add_to_expression!(
-            obj,
-            n.location.sizes[1].variable_operating_cost[t],
-            model[:process][n, t],
-        )
-
-        # Storage costs
-        add_to_expression!(obj, n.location.storage_cost[t], model[:store][n, t])
-
-        # Expansion costs
-        if t < T
-            add_to_expression!(
-                obj,
-                slope_open(n.location, t) - slope_open(n.location, t + 1),
-                model[:expansion][n, t],
+        # Stage 2: Flows, disposal, capacity & storage
+        # =====================================================================
+        @stage 2 begin
+            @uncertain graph
+            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,
             )
-        else
-            add_to_expression!(obj, slope_open(n.location, t), model[:expansion][n, t])
-        end
-    end
 
-    # Plant shipping node costs
-    for n in values(graph.plant_shipping_nodes), t = 1:T
-        # Disposal costs
-        add_to_expression!(
-            obj,
-            n.location.disposal_cost[n.product][t],
-            model[:plant_dispose][n, t],
-        )
-    end
+            # Var: plant_dispose
+            @recourse(
+                model,
+                plant_dispose[n in 1:PSN, t in 1:T],
+                lower_bound = 0,
+                upper_bound = psn[n].location.disposal_limit[psn[n].product][t],
+            )
 
-    # Collection shipping node costs
-    for n in values(graph.collection_shipping_nodes), t = 1:T
-        # Disposal costs
-        add_to_expression!(
-            obj,
-            n.location.product.disposal_cost[t],
-            model[:collection_dispose][n, t],
-        )
-    end
+            # Var: collection_dispose/
+            @recourse(
+                model,
+                collection_dispose[n in 1:CSN, t in 1:T],
+                lower_bound = 0,
+                upper_bound = graph.collection_shipping_nodes[n].location.amount[t],
+            )
 
-    @objective(model, Min, obj)
-end
+            # Var: store
+            @recourse(
+                model,
+                store[
+                    n in 1:PN,
+                    t in 1:T,
+                ],
+                lower_bound = 0,
+                upper_bound = pn[n].location.storage_limit,
+            )
 
+            # Var: process
+            @recourse(
+                model,
+                process[
+                    n in 1:PN,
+                    t in 1:T,
+                ],
+                lower_bound = 0,
+            )
 
-function create_shipping_node_constraints!(model::JuMP.Model)
-    graph, T = model[:graph], model[:instance].time
-    model[:eq_balance] = OrderedDict()
-    for t = 1:T
-        # Collection centers
-        for n in graph.collection_shipping_nodes
-            model[:eq_balance][n, t] = @constraint(
+            # Var: capacity
+            @recourse(
                 model,
-                sum(model[:flow][a, t] for a in n.outgoing_arcs) ==
-                n.location.amount[t] - model[:collection_dispose][n, t],
+                capacity[
+                    n in 1:PN,
+                    t in 1:T,
+                ],
+                lower_bound = 0,
+                upper_bound = pn[n].location.sizes[2].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
-        for n in graph.plant_shipping_nodes
-            @constraint(
+            # Var: expansion
+            @recourse(
                 model,
-                sum(model[:flow][a, t] for a in n.incoming_arcs) ==
-                sum(model[:flow][a, t] for a in n.outgoing_arcs) +
-                model[:plant_dispose][n, t]
+                expansion[
+                    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
-    end
 
-end
+            # Objective function
+            @objective(
+                model,
+                Min,
+                sum(
+                    # Transportation costs
+                    pn[n].location.input.transportation_cost[t] * 
+                        a.values["distance"] *
+                        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] +
+
+                    # 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]
+                    )
+
+                    for n in 1:PN
+                    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
+                )
+            )
 
+            # Process node constraints
+            for t = 1:T, n in 1:PN
+                node = pn[n]
+
+                # Output amount is implied by amount processed
+                for arc in node.outgoing_arcs
+                    @constraint(
+                        model,
+                        flow[arc.index, t] == arc.values["weight"] * process[n, t]
+                    )
+                end
+        
+                # If plant is closed, capacity is zero
+                @constraint(
+                    model,
+                    capacity[n, t] <= node.location.sizes[2].capacity * is_open[n, t]
+                )
+        
+                # If plant is open, capacity is greater than base
+                @constraint(
+                    model,
+                    capacity[n, t] >= node.location.sizes[1].capacity * is_open[n, t]
+                )
+        
+                # Capacity is linked to expansion
+                @constraint(
+                    model,
+                    capacity[n, t] <=
+                        node.location.sizes[1].capacity + expansion[n, t]
+                )
+        
+                # Can only process up to capacity
+                @constraint(model, process[n, t] <= capacity[n, t])
+        
+                if t > 1
+                    # Plant capacity can only increase over time
+                    @constraint(model, capacity[n, t] >= capacity[n, t-1])
+                    @constraint(model, expansion[n, t] >= expansion[n, t-1])
+                end
+        
+                # Amount received equals amount processed plus stored
+                store_in = 0
+                if t > 1
+                    store_in = store[n, t-1]
+                end
+                if t == T
+                    @constraint(model, store[n, t] == 0)
+                end
+                @constraint(
+                    model,
+                    sum(
+                        flow[arc.index, t]
+                        for arc in node.incoming_arcs
+                    ) + store_in == store[n, t] + process[n, t]
+                )
 
-function create_process_node_constraints!(model::JuMP.Model)
-    graph, T = model[:graph], model[:instance].time
+            end
 
-    for t = 1:T, n in graph.process_nodes
-        input_sum = AffExpr(0.0)
-        for a in n.incoming_arcs
-            add_to_expression!(input_sum, 1.0, model[:flow][a, t])
-        end
+            # Material flow at collection shipping nodes
+            @constraint(
+                model,
+                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]
+            )
 
-        # Output amount is implied by amount processed
-        for a in n.outgoing_arcs
+            # Material flow at plant shipping nodes
             @constraint(
                 model,
-                model[:flow][a, t] == a.values["weight"] * model[:process][n, t]
+                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]
             )
-        end
 
-        # If plant is closed, capacity is zero
-        @constraint(
-            model,
-            model[:capacity][n, t] <= n.location.sizes[2].capacity * model[:is_open][n, t]
-        )
-
-        # If plant is open, capacity is greater than base
-        @constraint(
-            model,
-            model[:capacity][n, t] >= n.location.sizes[1].capacity * model[:is_open][n, t]
-        )
-
-        # Capacity is linked to expansion
-        @constraint(
-            model,
-            model[:capacity][n, t] <=
-            n.location.sizes[1].capacity + model[:expansion][n, t]
-        )
-
-        # Can only process up to capacity
-        @constraint(model, model[:process][n, t] <= model[:capacity][n, t])
-
-        if t > 1
-            # Plant capacity can only increase over time
-            @constraint(model, model[:capacity][n, t] >= model[:capacity][n, t-1])
-            @constraint(model, model[:expansion][n, t] >= model[:expansion][n, t-1])
+            # 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
 
-        # Amount received equals amount processed plus stored
-        store_in = 0
-        if t > 1
-            store_in = model[:store][n, t-1]
-        end
-        if t == T
-            @constraint(model, model[:store][n, t] == 0)
-        end
-        @constraint(
-            model,
-            input_sum + store_in == model[:store][n, t] + model[:process][n, t]
-        )
+    ξ = [
+        @scenario graph = graphs[i] probability = probs[i]
+        for i in 1:length(graphs)
+    ]
 
+    sp = instantiate(model, ξ; optimizer)
 
-        # 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[:is_open][n, t] == model[:is_open][n, t-1] + model[:open_plant][n, t]
-            )
-        else
-            @constraint(model, model[:is_open][n, t] == model[:open_plant][n, t])
-        end
+    return sp
+end
 
-        # Plant can only be opened during building period
-        if t ∉ model[:instance].building_period
-            @constraint(model, model[:open_plant][n, t] == 0)
-        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
diff --git a/src/model/getsol.jl b/src/model/getsol.jl
index 4ed0097..095e35b 100644
--- a/src/model/getsol.jl
+++ b/src/model/getsol.jl
@@ -2,12 +2,31 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # 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
 
+    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(
         "Plants" => 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)
-    plant_to_shipping_nodes = OrderedDict()
-    for p in instance.plants
-        plant_to_shipping_nodes[p] = []
-        for a in plant_to_process_node[p].outgoing_arcs
-            push!(plant_to_shipping_nodes[p], a.dest)
-        end
+    pn = graph.process_nodes
+    psn = graph.plant_shipping_nodes
+
+    plant_to_process_node_index = OrderedDict(
+        pn[n].location => n
+        for n in 1:length(pn)
+    )
+
+    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
 
     # Products
-    for n in graph.collection_shipping_nodes
+    for n in 1:CSN
+        node = csn[n]
         location_dict = OrderedDict{Any,Any}(
-            "Latitude (deg)" => n.location.latitude,
-            "Longitude (deg)" => n.location.longitude,
-            "Amount (tonne)" => n.location.amount,
-            "Dispose (tonne)" =>
-                [JuMP.value(model[:collection_dispose][n, t]) for t = 1:T],
-            "Disposal cost (\$)" =>
-                [JuMP.value(model[:collection_dispose][n, t]) * n.location.product.disposal_cost[t] for t = 1:T]
+            "Latitude (deg)" => node.location.latitude,
+            "Longitude (deg)" => node.location.longitude,
+            "Amount (tonne)" => node.location.amount,
+            "Dispose (tonne)" => [
+                value(model[2, :collection_dispose][n, 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
             location_dict["Marginal cost (\$/tonne)"] = [
-                round(abs(JuMP.shadow_price(model[:eq_balance][n, t])), digits = 2) for
-                t = 1:T
+                round(abs(shadow_price(model[2, :eq_balance_centers][n, t])), digits=2) for t = 1:T
             ]
         end
-        if n.product.name ∉ keys(output["Products"])
-            output["Products"][n.product.name] = OrderedDict()
+        if node.product.name ∉ keys(output["Products"])
+            output["Products"][node.product.name] = OrderedDict()
         end
-        output["Products"][n.product.name][n.location.name] = location_dict
+        output["Products"][node.product.name][node.location.name] = location_dict
     end
 
     # Plants
     for plant in instance.plants
         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}(
             "Input" => OrderedDict(),
             "Output" =>
@@ -75,39 +104,39 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
             "Latitude (deg)" => plant.latitude,
             "Longitude (deg)" => plant.longitude,
             "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 (\$)" => [
-                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
             ],
             "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] +
-                JuMP.value(model[:expansion][process_node, t]) *
+                value(model[2, :expansion][n, t]) *
                 slope_fix_oper_cost(plant, t) for t = 1:T
             ],
             "Expansion cost (\$)" => [
                 (
                     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
                         slope_open(plant, t) * (
-                            JuMP.value(model[:expansion][process_node, t]) -
-                            JuMP.value(model[:expansion][process_node, t-1])
+                            value(model[2, :expansion][n, t]) -
+                            value(model[2, :expansion][n, t-1])
                         )
                     end
                 ) for t = 1:T
             ],
             "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 (\$)" => [
-                JuMP.value(model[:process][process_node, t]) *
+                value(model[2, :process][n, t]) *
                 plant.sizes[1].variable_operating_cost[t] for t = 1:T
             ],
             "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 (\$)" => [
-                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
             ],
         )
@@ -120,7 +149,7 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
 
         # Inputs
         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
                 continue
             end
@@ -178,19 +207,20 @@ function get_solution(model::JuMP.Model; marginal_costs = true)
         end
 
         # 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
             plant_dict["Total output"][product_name] = zeros(T)
             plant_dict["Output"]["Send"][product_name] = product_dict = OrderedDict()
 
             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
                 skip_plant = false
                 plant_dict["Output"]["Dispose"][product_name] =
                     disposal_dict = OrderedDict()
                 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["Amount (tonne)"][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
 
             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
                     continue
                 end
diff --git a/src/model/solve.jl b/src/model/solve.jl
index 12fb30d..81af576 100644
--- a/src/model/solve.jl
+++ b/src/model/solve.jl
@@ -2,14 +2,14 @@
 # Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 # 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()
-    return optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
+    return optimizer_with_attributes(HiGHS.Optimizer)
 end
 
 function _get_default_lp_optimizer()
-    return optimizer_with_attributes(Clp.Optimizer, "LogLevel" => 0)
+    return optimizer_with_attributes(HiGHS.Optimizer)
 end
 
 
@@ -25,61 +25,81 @@ function _print_graph_stats(instance::Instance, graph::Graph)::Nothing
     return
 end
 
-function solve(
-    instance::Instance;
-    optimizer = nothing,
-    output = nothing,
-    graph = nothing,
-    marginal_costs = true,
-    return_model = false,
+function solve_stochastic(;
+    scenarios::Vector{String},
+    probs::Vector{Float64},
+    optimizer,
 )
+    @info "Reading instance files..."
+    instances = [parsefile(sc) for sc in scenarios]
 
-    milp_optimizer = lp_optimizer = optimizer
-    if optimizer === nothing
-        milp_optimizer = _get_default_milp_optimizer()
-        lp_optimizer = _get_default_lp_optimizer()
-    end
+    @info "Building graphs..."
+    graphs = [build_graph(inst) for inst in instances]
 
-    if graph === nothing
-        @info "Building graph..."
-        graph = RELOG.build_graph(instance)
-        _print_graph_stats(instance, graph)
-    end
+    @info "Building stochastic model..."
+    sp = RELOG.build_model(instances[1], graphs, probs; optimizer)
+
+    @info "Optimizing stochastic model..."
+    optimize!(sp)
 
-    @info "Building optimization model..."
-    model = RELOG.build_model(instance, graph, milp_optimizer)
+    @info "Extracting solution..."
+    solutions = [
+        get_solution(instances[i], graphs[i], sp, i)
+        for i in 1:length(instances)
+    ]
 
-    @info "Optimizing MILP..."
-    JuMP.optimize!(model)
+    return solutions
+end
+
+function solve(
+    instance::Instance;
+    optimizer=nothing,
+    marginal_costs=true,
+    return_model=false
+)
+    @info "Building graph..."
+    graph = RELOG.build_graph(instance)
+    _print_graph_stats(instance, graph)
+
+    @info "Building model..."
+    model = RELOG.build_model(instance, [graph], [1.0]; optimizer)
+
+    @info "Optimizing model..."
+    optimize!(model)
     if !has_values(model)
         error("No solution available")
     end
-    solution = get_solution(model, marginal_costs = false)
+
+    @info "Extracting solution..."
+    solution = get_solution(instance, graph, model, 1)
 
     if marginal_costs
         @info "Re-optimizing with integer variables fixed..."
-        all_vars = JuMP.all_variables(model)
-        vals = OrderedDict(var => JuMP.value(var) for var in all_vars)
-        JuMP.set_optimizer(model, lp_optimizer)
-        for var in all_vars
-            if JuMP.is_binary(var)
-                JuMP.unset_binary(var)
-                JuMP.fix(var, vals[var])
-            end
+        open_plant_vals = value.(model[1, :open_plant])
+        is_open_vals = value.(model[1, :is_open])
+
+        for n in 1:length(graph.process_nodes), t in 1:instance.time
+            unset_binary(model[1, :open_plant][n, t])
+            unset_binary(model[1, :is_open][n, t])
+            fix(
+                model[1, :open_plant][n, t],
+                open_plant_vals[n, t]
+            )
+            fix(
+                model[1, :is_open][n, t],
+                is_open_vals[n, t]
+            )
+            
         end
-        JuMP.optimize!(model)
+        optimize!(model)
         if has_values(model)
             @info "Extracting solution..."
-            solution = get_solution(model, marginal_costs = true)
+            solution = get_solution(instance, graph, model, 1, marginal_costs=true)
         else
             @warn "Error computing marginal costs. Ignoring."
         end
     end
 
-    if output !== nothing
-        write(solution, output)
-    end
-
     if return_model
         return solution, model
     else
@@ -87,13 +107,13 @@ function solve(
     end
 end
 
-function solve(filename::AbstractString; heuristic = false, kwargs...)
+function solve(filename::AbstractString; heuristic=false, kwargs...)
     @info "Reading $filename..."
     instance = RELOG.parsefile(filename)
     if heuristic && instance.time > 1
         @info "Solving single-period version..."
         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..."
         selected_pairs = []
         for (plant_name, plant_dict) in csol["Plants"]