Blending Problem

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.

Importing Modules

The gamspy Python package is loaded in order to access GAMSPy syntax.

from gamspy import Container, Set, Parameter, Variable, Equation, Model, Sum, Sense
import numpy as np

Data

The company can purchase nine different alloys, labeled ‘a’ through ‘i’. Each alloy has a specific composition of lead, zinc, and tin, which is detailed in the compdat parameter. The cost of each alloy is listed in the price parameter, and the required amounts of lead, zinc, and tin are given by the rb parameter.

m = Container()

# Set
alloy = Set(
    container=m, name="alloy", records=["a", "b", "c", "d", "e", "f", "g", "h", "i"]
)
elem = Set(container=m, name="elem", records=["lead", "zinc", "tin"])

# Data
compdat = Parameter(
    container=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],
        ]
    ),
)
price = Parameter(
    container=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]),
)
rb = Parameter(container=m, name="rb", domain=elem, records=np.array([30, 30, 40]))

# Variable
v = Variable(container=m, name="v", domain=alloy, type="Positive")

Model

pc = Equation(container=m, name="pc", domain=elem)
mb = Equation(container=m, name="mb")

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

b1 = Model(
    container=m,
    name="b1",
    equations=[pc],
    problem="LP",
    sense=Sense.MIN,
    objective=Sum(alloy, price[alloy] * v[alloy]),
)

b2 = Model(
    container=m,
    name="b2",
    equations=[pc, mb],
    problem="LP",
    sense=Sense.MIN,
    objective=Sum(alloy, price[alloy] * v[alloy]),
)

Solve

report = Parameter(container=m, name="report", domain=[alloy, "*"])

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

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

Reporting

report.pivot()

	blend-1	blend-2
b	0.6	0.6
d	0.4	0.4