Blending Problem I (BLEND)

blend.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_blend.html
## LICENSETYPE: Demo
## MODELTYPE: LP
## KEYWORDS: linear programming, blending problem, manufacturing, alloy blending


Blending Problem I (BLEND)

A company wishes to produce a lead-zinc-tin alloy at minimal cost.
The problem is to blend a new alloy from other purchased alloys.


Dantzig, G B, Chapter 3.4. In Linear Programming and Extensions.
Princeton University Press, Princeton, New Jersey, 1963.
"""

from __future__ import annotations

import os

import numpy as np

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():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
    )

    # Set
    alloy = Set(
        m,
        name="alloy",
        records=["a", "b", "c", "d", "e", "f", "g", "h", "i"],
        description="products on the market",
    )
    elem = Set(
        m,
        name="elem",
        records=["lead", "zinc", "tin"],
        description="required elements",
    )

    # Data
    compdat = Parameter(
        m,
        name="compdat",
        domain=[elem, alloy],
        records=np.array([
            [10, 10, 40, 60, 30, 30, 30, 50, 20],
            [10, 30, 50, 30, 30, 40, 20, 40, 30],
            [80, 60, 10, 10, 40, 30, 50, 10, 50],
        ]),
        description="composition data (pct)",
    )
    price = Parameter(
        m,
        name="price",
        domain=alloy,
        records=np.array([4.1, 4.3, 5.8, 6.0, 7.6, 7.5, 7.3, 6.9, 7.3]),
        description="composition data (price)",
    )
    rb = Parameter(
        m,
        name="rb",
        domain=elem,
        records=np.array([30, 30, 40]),
        description="required blend",
    )

    # Variable
    v = Variable(
        m,
        name="v",
        domain=alloy,
        type="Positive",
        description="purchase of alloy (pounds)",
    )

    # Equation
    objective = Sum(alloy, price[alloy] * v[alloy])
    pc = Equation(m, name="pc", domain=elem, description="purchase constraint")
    mb = Equation(m, name="mb", description="material balance")

    pc[elem] = Sum(alloy, compdat[elem, alloy] * v[alloy]) == rb[elem]
    mb[...] = Sum(alloy, v[alloy]) == 1

    b1 = Model(
        m,
        name="b1",
        equations=[pc],
        problem="LP",
        sense=Sense.MIN,
        objective=objective,
    )
    b2 = Model(
        m,
        name="b2",
        equations=[pc, mb],
        problem="LP",
        sense=Sense.MIN,
        objective=objective,
    )

    report = Parameter(
        m,
        name="report",
        domain=[alloy, "*"],
        description="comparison of model 1 and 2",
    )

    b1.solve()

    report[alloy, "blend-1"] = v.l[alloy]
    b2.solve()
    report[alloy, "blend-2"] = v.l[alloy]

    print(report.pivot())

    import math

    assert math.isclose(b2.objective_value, 4.980000, rel_tol=0.001)


if __name__ == "__main__":
    main()