# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020-2024, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using CSV
using DataFrames
using Gurobi
using JSON
using JuMP
using OrderedCollections
using Random

include("jumpext.jl")
include("dist.jl")

dict = OrderedDict
Time = Int
Random.seed!(42)

# Structs
# =========================================================================

Base.@kwdef struct Component
    name::String
end

Base.@kwdef struct Product
    name::String
    comp::Vector{Component}
end

Base.@kwdef struct Center
    name::String
    latitude::Float64
    longitude::Float64
    prod_out::Vector{Product}
end

Base.@kwdef struct Plant
    name::String
    latitude::Float64
    longitude::Float64
    prod_out::Vector{Product}
    prod_in::Product
end

Base.@kwdef struct Emission
    name::String
end

Base.show(
    io::IO,
    p::Union{Component,Product,Center,Plant,Emission},
) = print(io, p.name)

Base.@kwdef struct Instance
    T::UnitRange{Int}
    centers::Vector{Center}
    plants::Vector{Plant}
    products::Vector{Product}
    emissions::Vector{Emission}
    alpha_mix::dict{Tuple{Plant,Product,Component,Component},Float64}
    alpha_plant_emission::dict{Tuple{Plant,Emission,Time},Float64}
    alpha_tr_emission::dict{Tuple{Product,Emission,Time},Float64}
    c_acq::dict{Tuple{Center,Product,Time},Float64}
    c_center_disp::dict{Tuple{Center,Product,Time},Float64}
    c_emission::dict{Tuple{Emission,Time},Float64}
    c_fix::dict{Tuple{Plant,Time},Float64}
    c_open::dict{Tuple{Plant,Time},Float64}
    c_plant_disp::dict{Tuple{Plant,Product,Time},Float64}
    c_store::dict{Tuple{Center,Product,Time},Float64}
    c_tr::dict{Tuple{Product,Time},Float64}
    c_var::dict{Tuple{Plant,Time},Float64}
    m_cap::dict{Plant,Float64}
    m_center_disp::dict{Tuple{Product,Time},Float64}
    m_dist::dict{Tuple{Union{Center,Plant},Plant},Float64}
    m_emission::dict{Tuple{Emission,Time},Float64}
    m_init::dict{Tuple{Center,Product,Component,Time},Float64}
    m_plant_disp::dict{Tuple{Plant,Product,Time},Float64}
    m_store::dict{Tuple{Center,Product,Time},Float64}
end


# Generate
# ==============================================================================

function generate_data()
    # Time window
    T = 1:2

    # Cities
    cities_a = dict(
        "Chicago" => [41.881832, -87.623177],
        "New York City" => [40.712776, -74.005974],
        "Los Angeles" => [34.052235, -118.243683],
        "Houston" => [29.760427, -95.369804],
        "Phoenix" => [33.448376, -112.074036],
        "Philadelphia" => [39.952583, -75.165222],
        "San Antonio" => [29.424122, -98.493629],
        "San Diego" => [32.715736, -117.161087],
        "Dallas" => [32.776664, -96.796988],
        "San Jose" => [37.338208, -121.886329],
    )

    cities_b = dict(
        "Chicago" => [41.881832, -87.623177],
        "Phoenix" => [33.448376, -112.074036],
        "Dallas" => [32.776664, -96.796988],
    )

    # Components
    film = Component("Film")
    paper = Component("Paper")
    cardboard = Component("Cardboard")

    # Products
    waste = Product(
        name="Waste",
        comp=[film, paper, cardboard],
    )
    film_bale = Product(
        name="Film bale",
        comp=[film, paper, cardboard],
    )
    cardboard_bale = Product(
        name="Cardboard bale",
        comp=[paper, cardboard],
    )
    cardboard_sheets = Product(
        name="Cardboard sheets",
        comp=[cardboard],
    )
    products = [waste, film_bale, cardboard_bale, cardboard_sheets]

    # Centers
    centers = [
        Center(
            name="Collection ($city_name)",
            latitude=city_lat,
            longitude=city_lon,
            prod_out=[waste],
        )
        for (city_name, (city_lat, city_lon)) in cities_a
    ]

    # Plants
    plants_a = [
        Plant(
            name="MRF ($city_name)",
            latitude=city_lat,
            longitude=city_lon,
            prod_in=waste,
            prod_out=[film_bale, cardboard_bale],
        )
        for (city_name, (city_lat, city_lon)) in cities_b
    ]
    plants_b = [
        Plant(
            name="Paper Mill ($city_name)",
            latitude=city_lat,
            longitude=city_lon,
            prod_in=cardboard_bale,
            prod_out=[cardboard_sheets],
        )
        for (city_name, (city_lat, city_lon)) in cities_b
    ]
    plants = [plants_a; plants_b]

    # Emissions
    emissions = [
        Emission("CO2"),
    ]

    alpha_mix = dict(
        (p, r, cin, cout) => 0.0
        for p in plants
        for cin in p.prod_in.comp
        for r in p.prod_out
        for cout in r.comp
    )
    for p in plants_a
        alpha_mix[p, film_bale, film, film] = 0.98
        alpha_mix[p, film_bale, paper, paper] = 0.02
        alpha_mix[p, film_bale, cardboard, cardboard] = 0.02
        alpha_mix[p, cardboard_bale, paper, paper] = 0.02
        alpha_mix[p, cardboard_bale, cardboard, cardboard] = 0.75
    end
    for p in plants_b
        alpha_mix[p, cardboard_sheets, cardboard, cardboard] = 0.95
    end
    alpha_plant_emission = dict(
        (p, s, t) => 0.01
        for p in plants, s in emissions, t in T
    )
    alpha_tr_emission = dict(
        (r, s, t) => 0.01
        for r in products, s in emissions, t in T
    )
    c_acq = dict(
        (q, r, t) => 1.0
        for q in centers
        for r in q.prod_out, t in T
    )
    c_center_disp = dict(
        (q, r, t) => 0
        for q in centers
        for r in q.prod_out
        for t in T
    )
    c_emission = dict(
        (s, t) => 0.01
        for s in emissions, t in T
    )
    c_fix = dict(
        (p, t) => 1_000.0
        for p in plants, t in T
    )
    c_open = dict(
        (p, t) => 10_000.0
        for p in plants, t in T
    )
    c_plant_disp = dict(
        (p, r, t) => (r == cardboard_sheets ? -100.0 : -10.0)
        for p in plants
        for r in p.prod_out, t in T
    )
    c_store = dict(
        (q, r, t) => 1.0
        for q in centers
        for r in q.prod_out, t in T
    )
    c_tr = dict(
        (r, t) => 0.05
        for r in products, t in T
    )
    c_var = dict(
        (p, t) => 1.0
        for p in plants, t in T
    )
    m_cap = dict(
        p => 50000.0
        for p in plants
    )
    m_center_disp = dict(
        (r, t) => 0.0
        for r in products, t in T
    )
    metric = KnnDrivingDistance()
    m_dist = dict(
        (p, q) => _calculate_distance(
            p.latitude,
            p.longitude,
            q.latitude,
            q.longitude,
            metric,
        )
        for p in [plants; centers], q in plants
    )
    m_emission = dict(
        (s, t) => 1_000_000
        for s in emissions, t in T
    )
    m_init = dict()
    for q in centers, r in q.prod_out
        ratio = dict(
            c => rand(1:10)
            for c in r.comp
        )
        total = dict(
            t => rand(1:1000)
            for t in T
        )
        for c in r.comp, t in T
            m_init[q, r, c, t] = ratio[c] * total[t]
        end
    end
    m_plant_disp = dict(
        (p, r, t) => 1_000_000
        for p in plants
        for r in p.prod_out
        for t in T
    )
    m_store = dict(
        (q, r, t) => 1_000
        for q in centers
        for r in q.prod_out, t in T
    )

    return Instance(;
        T,
        centers,
        plants,
        products,
        emissions,
        alpha_mix,
        alpha_plant_emission,
        alpha_tr_emission,
        c_acq,
        c_center_disp,
        c_emission,
        c_fix,
        c_open,
        c_plant_disp,
        c_store,
        c_tr,
        c_var,
        m_cap,
        m_center_disp,
        m_dist,
        m_emission,
        m_init,
        m_plant_disp,
        m_store
    )
end

# Write
# ==============================================================================

function write_json(data, filename)
    json = dict()
    json["parameters"] = dict(
        "time horizon (years)" => data.T.stop,
    )
    json["products"] = dict(
        r.name => dict(
            "components" => [c.name for c in r.comp],
            "disposal limit (tonne)" => [
                data.m_center_disp[r, t]
                for t in data.T
            ],
            "transportation cost (\$/km/tonne)" => [
                data.c_tr[r, t]
                for t in data.T
            ],
            "transportation emissions (tonne/km/tonne)" => dict(
                s => [
                    data.alpha_tr_emission[r, s, t]
                    for t in data.T
                ]
                for s in data.emissions
            )
        )
        for r in data.products
    )
    json["centers"] = dict(
        q.name => dict(
            "latitude" => q.latitude,
            "longitude" => q.longitude,
            "output" => dict(
                r.name => dict(
                    "initial amount (tonne)" => [
                        data.m_init[q, r, c, t]
                        for c in r.comp, t in data.T
                    ],
                    "disposal cost (\$/tonne)" => [
                        data.c_center_disp[q, r, t]
                        for t in data.T
                    ],
                    "storage cost (\$/tonne)" => [
                        data.c_store[q, r, t]
                        for t in data.T
                    ],
                    "storage limit (tonne)" => [
                        data.m_store[q, r, t]
                        for t in data.T
                    ],
                    "acquisition cost (\$/tonne)" => [
                        data.c_acq[q, r, t]
                        for t in data.T
                    ],
                )
                for r in q.prod_out
            )
        )
        for q in data.centers
    )
    json["plants"] = dict(
        p.name => dict(
            "latitude" => p.latitude,
            "longitude" => p.longitude,
            "input" => p.prod_in.name,
            "output" => dict(
                r.name => dict(
                    "output matrix" => [
                        data.alpha_mix[p, r, c_in, c_out]
                        for c_in in p.prod_in.comp, c_out in r.comp
                    ],
                    "disposal limit (tonne)" => [
                        data.m_plant_disp[p, r, t]
                        for t in data.T
                    ],
                    "disposal cost (\$/tonne)" => [
                        data.c_plant_disp[p, r, t]
                        for t in data.T
                    ],
                )
                for r in p.prod_out
            ),
            "fixed operating cost (\$)" => [
                data.c_fix[p, t]
                for t in data.T
            ],
            "variable operating cost (\$/tonne)" => [
                data.c_var[p, t]
                for t in data.T
            ],
            "opening cost (\$)" => [
                data.c_open[p, t]
                for t in data.T
            ],
            "capacity (tonne)" => data.m_cap[p],
            "emissions (tonne/tonne)" => dict(
                s => [
                    data.alpha_plant_emission[p, s, t]
                    for t in data.T
                ]
                for s in data.emissions
            )
        )
        for p in data.plants
    )
    json["emissions"] = dict(
        s.name => dict(
            "penalty (\$/tonne)" => [
                data.c_emission[s, t]
                for t in data.T
            ],
            "limit (tonne)" => [
                data.m_emission[s, t]
                for t in data.T
            ]
        )
        for s in data.emissions
    )
    open(filename, "w") do io
        JSON.print(io, json, 2)
    end
end

# Read
# ==============================================================================
function read_json(filename)
    json = JSON.parsefile(filename)
    T = 1:json["parameters"]["time horizon (years)"]
    centers = []
    components_by_name = dict()
    emissions = []
    emissions_by_name = dict()
    plants = []
    products = []
    products_by_name = dict()
    alpha_mix = dict()
    alpha_plant_emission = dict()
    alpha_tr_emission = dict()
    c_acq = dict()
    c_center_disp = dict()
    c_emission = dict()
    c_fix = dict()
    c_open = dict()
    c_plant_disp = dict()
    c_store = dict()
    c_tr = dict()
    c_var = dict()
    m_cap = dict()
    m_center_disp = dict()
    m_emission = dict()
    m_init = dict()
    m_plant_disp = dict()
    m_store = dict()

    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

    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
        r = Product(name, comp)
        products_by_name[name] = r
        push!(products, r)
        for t in T
            m_center_disp[r, t] = prod_data["disposal limit (tonne)"][t]
            c_tr[r, t] = prod_data["transportation cost (\$/km/tonne)"][t]
        end
        for (s_name, s_data) in prod_data["transportation emissions (tonne/km/tonne)"]
            s = emissions_by_name[s_name]
            for t in T
                alpha_tr_emission[r, s, t] = s_data[t]
            end
        end
    end

    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]
            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]
            for t in T
                alpha_plant_emission[p, s, t] = s_data[t]
            end
        end
    end

    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
    )

    return Instance(;
        T,
        centers,
        plants,
        products,
        emissions,
        alpha_mix,
        alpha_plant_emission,
        alpha_tr_emission,
        c_acq,
        c_center_disp,
        c_emission,
        c_fix,
        c_open,
        c_plant_disp,
        c_store,
        c_tr,
        c_var,
        m_cap,
        m_center_disp,
        m_dist,
        m_emission,
        m_init,
        m_plant_disp,
        m_store
    )
end


# Run
# ==============================================================================

function generate_json()
    @info "Generating data"
    data = generate_data()

    @info "Writing JSON file"
    write_json(data, "output-3/case.json")
end

function solve()
    @info "Reading JSON file"
    data = read_json("output-3/case.json")

    T = data.T
    centers = data.centers
    plants = data.plants
    products = data.products
    emissions = data.emissions

    model = Model(Gurobi.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
        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

    # 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)
        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)
            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)
        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

    # 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 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

    # 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 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

    # Fixed cost
    for p in plants, t in T
        add_to_expression!(
            obj,
            data.c_fix[p, t],
            x_open[p, t]
        )
    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

    # 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
    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
    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

    @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
        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
        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
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."source" = String[]
    df."destination" = String[]
    df."product" = String[]
    df."component" = String[]
    df."time" = Int[]
    df."distance (km)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."transportation cost (\$)" = Float64[]
    df."variable operating cost (\$)" = Float64[]
    for (q, p, r) in E, c in r.comp, t in T
        if value(y[q, p, r, c, t]) ≈ 0
            continue
        end
        push!(
            df,
            [
                q.name,
                p.name,
                r.name,
                c.name,
                t,
                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])
            ]
        )
    end
    CSV.write("output-3/transp.csv", df)

    # Report: Centers
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."center" = String[]
    df."product" = String[]
    df."component" = String[]
    df."time" = Int[]
    df."amount available (tonne)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."amount stored (tonne)" = Float64[]
    df."amount disposed (tonne)" = Float64[]
    df."acquisition cost (\$)" = Float64[]
    df."storage cost (\$)" = Float64[]
    df."disposal cost (\$)" = Float64[]
    for q in centers, r in q.prod_out, c in r.comp, t in T
        push!(
            df,
            [
                q.name,
                r.name,
                c.name,
                t,
                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
                ),
                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-3/centers.csv", df)

    # Report: Plants
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."plant" = String[]
    df."time" = Int[]
    df."is open?" = Float64[]
    df."opening cost (\$)" = Float64[]
    df."fixed operating cost (\$)" = Float64[]
    for p in plants, t in T
        push!(
            df,
            [
                p.name,
                t,
                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-3/plants.csv", df)

    # Report: Plant Outputs
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."plant" = String[]
    df."product" = String[]
    df."component" = String[]
    df."time" = Int[]
    df."amount produced (tonne)" = Float64[]
    df."amount disposed (tonne)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."disposal cost (\$)" = Float64[]
    for p in plants, r in p.prod_out, c in r.comp, t in T
        if value(z_prod[p, r, c, t]) ≈ 0
            continue
        end
        push!(
            df,
            [
                p.name,
                r.name,
                c.name,
                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
                ),
                data.c_plant_disp[p, r, t] * value(z_plant_disp[p, r, c, t]),
            ]
        )
    end
    CSV.write("output-3/plant-outputs.csv", df)

    # Report: Plant Emissions
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."plant" = String[]
    df."emission" = String[]
    df."time" = Int[]
    df."amount emitted (tonne)" = Float64[]
    df."emission cost (\$)" = Float64[]
    for p in plants, s in emissions, t in T
        push!(
            df,
            [
                p.name,
                s.name,
                t,
                value(z_plant_emissions[p, s, t]),
                data.c_emission[s, t] * value(z_plant_emissions[p, s, t]),
            ]
        )
    end
    CSV.write("output-3/plant-emissions.csv", df)

    # Report: Transportation Emissions
    # -------------------------------------------------------------------------
    df = DataFrame()
    df."source" = String[]
    df."destination" = String[]
    df."product" = String[]
    df."emission" = String[]
    df."time" = Int[]
    df."distance (km)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."amount emitted (tonne)" = Float64[]
    df."emission cost (\$)" = Float64[]
    for (q, p, r) in E, s in emissions, t in T
        if value(y_total[q, p, r, t]) ≈ 0
            continue
        end
        push!(
            df,
            [
                q.name,
                p.name,
                r.name,
                s.name,
                t,
                data.m_dist[q, p],
                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-3/transp-emissions.csv", df)
end