Mexico Steel - Small Static (MEXSS)#

mexss.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_mexss.html
## LICENSETYPE: Demo
## MODELTYPE: LP
## KEYWORDS: linear programming, production problem, distribution problem, scheduling, micro economics, steel industry


Mexico Steel - Small Static (MEXSS)

A simplified representation of the Mexican steel sector is used
to introduce a process industry production and distribution
scheduling problem.


Kendrick, D, Meeraus, A, and Alatorre, J, The Planning of Investment
Programs in the Steel Industry. The Johns Hopkins University Press,
Baltimore and London, 1984.

A scanned version of this out-of-print book is accessible at
http://www.gams.com/docs/pdf/steel_investment.pdf
"""

from __future__ import annotations

import os

import pandas

from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Parameter
from gamspy import Sense
from gamspy import Set
from gamspy import Sum
from gamspy import Variable


def main():
    cont = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
    )

    # Prepare data
    steel_plants = ["ahmsa", "fundidora", "sicartsa", "hylsa", "hylsap"]
    markets = ["mexico-df", "monterrey", "guadalaja"]
    commodities = [
        "pellets",
        "coke",
        "nat-gas",
        "electric",
        "scrap",
        "pig-iron",
        "sponge",
        "steel",
    ]
    final_products = ["steel"]
    intermediate_products = ["sponge", "pig-iron"]
    raw_materials = ["pellets", "coke", "nat-gas", "electric", "scrap"]
    processes = ["pig-iron", "sponge", "steel-oh", "steel-el", "steel-bof"]
    productive_units = [
        "blast-furn",
        "openhearth",
        "bof",
        "direct-red",
        "elec-arc",
    ]

    io_coefficients = pandas.DataFrame([
        ["pellets", "pig-iron", -1.58],
        ["pellets", "sponge", -1.38],
        ["coke", "pig-iron", -0.63],
        ["nat-gas", "sponge", -0.57],
        ["electric", "steel-el", -0.58],
        ["scrap", "steel-oh", -0.33],
        ["scrap", "steel-bof", -0.12],
        ["pig-iron", "pig-iron", 1.00],
        ["pig-iron", "steel-oh", -0.77],
        ["pig-iron", "steel-bof", -0.95],
        ["sponge", "sponge", 1.00],
        ["sponge", "steel-el", -1.09],
        ["steel", "steel-oh", 1.00],
        ["steel", "steel-el", 1.00],
        ["steel", "steel-bof", 1.00],
    ])

    capacity_utilization = pandas.DataFrame([
        ["blast-furn", "pig-iron", 1.0],
        ["openhearth", "steel-oh", 1.0],
        ["bof", "steel-bof", 1.0],
        ["direct-red", "sponge", 1.0],
        ["elec-arc", "steel-el", 1.0],
    ])

    capacities_of_units = pandas.DataFrame([
        ["blast-furn", "ahmsa", 3.25],
        ["blast-furn", "fundidora", 1.40],
        ["blast-furn", "sicartsa", 1.10],
        ["openhearth", "ahmsa", 1.50],
        ["openhearth", "fundidora", 0.85],
        ["bof", "ahmsa", 2.07],
        ["bof", "fundidora", 1.50],
        ["bof", "sicartsa", 1.30],
        ["direct-red", "hylsa", 0.98],
        ["direct-red", "hylsap", 1.00],
        ["elec-arc", "hylsa", 1.13],
        ["elec-arc", "hylsap", 0.56],
    ])

    rail_distances = pandas.DataFrame([
        ["ahmsa", "mexico-df", 1204],
        ["ahmsa", "monterrey", 218],
        ["ahmsa", "guadalaja", 1125],
        ["ahmsa", "export", 739],
        ["fundidora", "mexico-df", 1017],
        ["fundidora", "guadalaja", 1030],
        ["fundidora", "export", 521],
        ["sicartsa", "mexico-df", 819],
        ["sicartsa", "monterrey", 1305],
        ["sicartsa", "guadalaja", 704],
        ["hylsa", "mexico-df", 1017],
        ["hylsa", "guadalaja", 1030],
        ["hylsa", "export", 521],
        ["hylsap", "mexico-df", 185],
        ["hylsap", "monterrey", 1085],
        ["hylsap", "guadalaja", 760],
        ["hylsap", "export", 315],
        ["import", "mexico-df", 428],
        ["import", "monterrey", 521],
        ["import", "guadalaja", 300],
    ])

    product_prices = pandas.DataFrame([
        ["pellets", "domestic", 18.7],
        ["coke", "domestic", 52.17],
        ["nat-gas", "domestic", 14.0],
        ["electric", "domestic", 24.0],
        ["scrap", "domestic", 105.0],
        ["steel", "import", 150],
        ["steel", "export", 140],
    ])

    demand_distribution = pandas.DataFrame(
        [["mexico-df", 55], ["monterrey", 30], ["guadalaja", 15]]
    )

    dt = 5.209  # total demand for final goods in 1979
    rse = 40  # raw steel equivalence
    eb = 1.0  # export bound

    # Set
    i = Set(
        cont,
        name="i",
        records=pandas.DataFrame(steel_plants),
        description="steel plants",
    )
    j = Set(
        cont,
        name="j",
        records=pandas.DataFrame(markets),
        description="markets",
    )
    c = Set(
        cont,
        name="c",
        records=pandas.DataFrame(commodities),
        description="commidities",
    )
    cf = Set(
        cont,
        name="cf",
        records=pandas.DataFrame(final_products),
        domain=c,
        description="final products",
    )
    ci = Set(
        cont,
        name="ci",
        records=pandas.DataFrame(intermediate_products),
        domain=c,
        description="intermediate products",
    )
    cr = Set(
        cont,
        name="cr",
        records=pandas.DataFrame(raw_materials),
        domain=c,
        description="raw materials",
    )
    p = Set(
        cont,
        name="p",
        records=pandas.DataFrame(processes),
        description="processes",
    )
    m = Set(
        cont,
        name="m",
        records=pandas.DataFrame(productive_units),
        description="productive units",
    )

    # Data
    a = Parameter(
        cont,
        name="a",
        domain=[c, p],
        records=io_coefficients,
        description="input-output coefficients",
    )
    b = Parameter(
        cont,
        name="b",
        domain=[m, p],
        records=capacity_utilization,
        description="capacity utilization",
    )
    k = Parameter(
        cont,
        name="k",
        domain=[m, i],
        records=capacities_of_units,
        description="capacities of productive units",
    )
    dd = Parameter(
        cont,
        name="dd",
        domain=j,
        records=demand_distribution,
        description="distribution of demand",
    )
    d = Parameter(
        cont, name="d", domain=[c, j], description="demand for steel in 1979"
    )

    d["steel", j] = dt * (1 + rse / 100) * dd[j] / 100

    rd = Parameter(
        cont,
        name="rd",
        domain=["*", "*"],
        records=rail_distances,
        description="rail distances from plants to markets",
    )

    muf = Parameter(
        cont,
        name="muf",
        domain=[i, j],
        description="transport rate: final products",
    )
    muv = Parameter(
        cont, name="muv", domain=j, description="transport rate: imports"
    )
    mue = Parameter(
        cont, name="mue", domain=i, description="transport rate: exports"
    )

    muf[i, j] = (2.48 + 0.0084 * rd[i, j]).where[rd[i, j]]
    muv[j] = (2.48 + 0.0084 * rd["import", j]).where[rd["import", j]]
    mue[i] = (2.48 + 0.0084 * rd[i, "export"]).where[rd[i, "export"]]

    prices = Parameter(
        cont,
        name="prices",
        domain=[c, "*"],
        records=product_prices,
        description="product prices (us$ per unit)",
    )

    pd = Parameter(cont, name="pd", domain=c, description="domestic prices")
    pv = Parameter(cont, name="pv", domain=c, description="import prices")
    pe = Parameter(cont, name="pe", domain=c, description="export prices")

    pd[c] = prices[c, "domestic"]
    pv[c] = prices[c, "import"]
    pe[c] = prices[c, "export"]

    # Variable
    z = Variable(
        cont,
        name="z",
        domain=[p, i],
        type="Positive",
        description="process level",
    )
    x = Variable(
        cont,
        name="x",
        domain=[c, i, j],
        type="Positive",
        description="shipment of final products",
    )
    u = Variable(
        cont,
        name="u",
        domain=[c, i],
        type="Positive",
        description="purchase of domestic materials",
    )
    v = Variable(
        cont, name="v", domain=[c, j], type="Positive", description="imports"
    )
    e = Variable(
        cont, name="e", domain=[c, i], type="Positive", description="exports"
    )
    phipsi = Variable(cont, name="phipsi", description="raw material cost")
    philam = Variable(cont, name="philam", description="transport cost")
    phipi = Variable(cont, name="phipi", description="import cost")
    phieps = Variable(cont, name="phieps", description="export revenue")

    # Equation declaration
    mbf = Equation(
        cont,
        name="mbf",
        domain=[c, i],
        description="material balances: final products",
    )
    mbi = Equation(
        cont,
        name="mbi",
        domain=[c, i],
        description="material balances: intermediates",
    )
    mbr = Equation(
        cont,
        name="mbr",
        domain=[c, i],
        description="material balances: raw materials",
    )
    cc = Equation(
        cont,
        name="cc",
        domain=[m, i],
        description="capacity constraint",
    )
    mr = Equation(
        cont,
        name="mr",
        domain=[c, j],
        description="market requirements",
    )
    me = Equation(
        cont,
        name="me",
        domain=c,
        description="maximum export",
    )
    apsi = Equation(
        cont,
        name="apsi",
        description="accounting: raw material cost",
    )
    alam = Equation(
        cont,
        name="alam",
        description="accounting: transport cost",
    )
    api = Equation(cont, name="api", description="accounting: import cost")
    aeps = Equation(
        cont,
        name="aeps",
        description="accounting: export cost",
    )

    # Equation definition
    obj = phipsi + philam + phipi - phieps  # Total Cost

    mbf[cf, i] = Sum(p, a[cf, p] * z[p, i]) >= Sum(j, x[cf, i, j]) + e[cf, i]
    mbi[ci, i] = Sum(p, a[ci, p] * z[p, i]) >= 0
    mbr[cr, i] = Sum(p, a[cr, p] * z[p, i]) + u[cr, i] >= 0
    cc[m, i] = Sum(p, b[m, p] * z[p, i]) <= k[m, i]
    mr[cf, j] = Sum(i, x[cf, i, j]) + v[cf, j] >= d[cf, j]
    me[cf] = Sum(i, e[cf, i]) <= eb
    apsi[...] = phipsi == Sum((cr, i), pd[cr] * u[cr, i])
    alam[...] = philam == Sum((cf, i, j), muf[i, j] * x[cf, i, j]) + Sum(
        (cf, j), muv[j] * v[cf, j]
    ) + Sum((cf, i), mue[i] * e[cf, i])
    api[...] = phipi == Sum((cf, j), pv[cf] * v[cf, j])
    aeps[...] = phieps == Sum((cf, i), pe[cf] * e[cf, i])

    mexss = Model(
        cont,
        name="mexss",
        equations=cont.getEquations(),
        problem="LP",
        sense=Sense.MIN,
        objective=obj,
    )

    mexss.solve()

    import math

    assert math.isclose(mexss.objective_value, 538.8112, rel_tol=0.001)
    print(mexss.objective_value)


if __name__ == "__main__":
    main()