Stochastic Benders - Sequential GAMS Loop (SPBENDERS1)#

spbenders1.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_spbenders1.html
## LICENSETYPE: Demo
## MODELTYPE: LP
## KEYWORDS: linear programming, stochastic Benders algorithm, transportation problem


Stochastic Benders - Sequential GAMS Loop (SPBENDERS1)

This example demonstrates a stochastic Benders implementation for the
simple transport example.

This is the first example of a sequence of stochastic Benders
implementations using various methods to solve the master and
subproblem.

This first example implements the stochastic Benders algorithm using
sequential solves of the master and subproblems in a GAMS loop.
"""

from __future__ import annotations

import os

import pandas as pd

import gamspy.math as gams_math
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)),
    )

    # Prepare data
    cost = pd.DataFrame([
        ["f1", "d1", 2.49],
        ["f1", "d2", 5.21],
        ["f1", "d3", 3.76],
        ["f1", "d4", 4.85],
        ["f1", "d5", 2.07],
        ["f2", "d1", 1.46],
        ["f2", "d2", 2.54],
        ["f2", "d3", 1.83],
        ["f2", "d4", 1.86],
        ["f2", "d5", 4.76],
        ["f3", "d1", 3.26],
        ["f3", "d2", 3.08],
        ["f3", "d3", 2.60],
        ["f3", "d4", 3.76],
        ["f3", "d5", 4.45],
    ])

    scenarios = pd.DataFrame([
        ["lo", "d1", 150],
        ["lo", "d2", 100],
        ["lo", "d3", 250],
        ["lo", "d4", 300],
        ["lo", "d5", 600],
        ["lo", "prob", 0.25],
        ["mid", "d1", 160],
        ["mid", "d2", 120],
        ["mid", "d3", 270],
        ["mid", "d4", 325],
        ["mid", "d5", 700],
        ["mid", "prob", 0.50],
        ["hi", "d1", 170],
        ["hi", "d2", 135],
        ["hi", "d3", 300],
        ["hi", "d4", 350],
        ["hi", "d5", 800],
        ["hi", "prob", 0.25],
    ])

    cut_coefficients = pd.DataFrame([
        [idx, f"d{center}", 0]
        for idx in range(1, 26)
        for center in range(1, 6)
    ])

    # Set
    i = Set(m, name="i", records=["f1", "f2", "f3"], description="factories")
    j = Set(
        m,
        name="j",
        records=["d1", "d2", "d3", "d4", "d5"],
        description="distribution centers",
    )
    s = Set(m, name="s", records=["lo", "mid", "hi"], description="scenarios")

    # Data
    capacity = Parameter(
        m,
        name="capacity",
        domain=i,
        records=pd.DataFrame([["f1", 500], ["f2", 450], ["f3", 650]]),
        description="unit capacity at factories",
    )
    demand = Parameter(
        m,
        name="demand",
        domain=j,
        records=pd.DataFrame(
            [["d1", 160], ["d2", 120], ["d3", 270], ["d4", 325], ["d5", 700]]
        ),
        description="unit demand at distribution centers",
    )
    prodcost = Parameter(
        m, name="prodcost", records=14, description="unit production cost"
    )
    price = Parameter(m, name="price", records=24, description="sales price")
    wastecost = Parameter(
        m,
        name="wastecost",
        records=4,
        description="cost of removal of overstocked products",
    )
    transcost = Parameter(
        m,
        name="transcost",
        domain=[i, j],
        records=cost,
        description="unit transportation cost",
    )
    ScenarioData = Parameter(
        m,
        name="scenariodata",
        domain=[s, "*"],
        records=scenarios,
        description="possible outcomes for demand plus probabilities",
    )

    # Set
    iter = Set(
        m,
        name="iter",
        records=[f"{idx}" for idx in range(1, 26)],
        description="max Benders iterations",
    )
    itActive = Set(
        m, name="itActive", domain=iter, description="active Benders cuts"
    )

    # Parameter
    cutconst = Parameter(
        m,
        name="cutconst",
        domain=iter,
        records=pd.DataFrame([[f"{idx}", 0] for idx in range(1, 26)]),
        description="constants in optimality cuts",
    )
    cutcoeff = Parameter(
        m,
        name="cutcoeff",
        domain=[iter, j],
        records=cut_coefficients,
        description="coefficients in optimality cuts",
    )

    # Variable
    ship = Variable(
        m, name="ship", domain=[i, j], type="Positive", description="shipments"
    )
    product = Variable(m, name="product", domain=i, description="production")
    received = Variable(
        m, name="received", domain=j, description="quantity sent to market"
    )
    zmaster = Variable(
        m, name="zmaster", description="objective variable of master problem"
    )
    theta = Variable(m, name="theta", description="future profit")

    # Equation
    masterobj = Equation(
        m, name="masterobj", description="master objective function"
    )
    production = Equation(
        m,
        name="production",
        domain=i,
        description="calculate production in each factory",
    )
    receive = Equation(
        m,
        name="receive",
        domain=j,
        description="calculate quantity to be send to markets",
    )
    optcut = Equation(
        m, name="optcut", domain=iter, description="Benders optimality cuts"
    )

    masterobj[...] = zmaster == theta - Sum(
        (i, j), transcost[i, j] * ship[i, j]
    ) - Sum(i, prodcost * product[i])
    receive[j] = received[j] == Sum(i, ship[i, j])
    production[i] = product[i] == Sum(j, ship[i, j])
    optcut[itActive] = theta <= cutconst[itActive] + Sum(
        j, cutcoeff[itActive, j] * received[j]
    )
    product.up[i] = capacity[i]

    masterproblem = Model(
        m,
        name="masterproblem",
        equations=m.getEquations(),
        problem="LP",
        sense=Sense.MAX,
        objective=zmaster,
    )

    # Variable
    sales = Variable(
        m,
        name="sales",
        domain=j,
        type="Positive",
        description="sales (actually sold)",
    )
    waste = Variable(
        m,
        name="waste",
        domain=j,
        type="Positive",
        description="overstocked products",
    )
    zsub = Variable(
        m, name="zsub", description="objective variable of sub problem"
    )

    # Equation
    subobj = Equation(
        m, name="subobj", description="subproblem objective function"
    )
    selling = Equation(
        m, name="selling", domain=j, description="part of received is sold"
    )
    market = Equation(
        m, name="market", domain=j, description="upperbound on sales"
    )

    subobj[...] = zsub == Sum(j, price * sales[j]) - Sum(
        j, wastecost * waste[j]
    )
    selling[j] = sales[j] + waste[j] == received.l[j]
    market[j] = sales[j] <= demand[j]

    subproblem = Model(
        m,
        name="subproblem",
        equations=[subobj, selling, market],
        problem="LP",
        sense=Sense.MAX,
        objective=zsub,
    )

    # Scalar
    rgap = Parameter(m, name="rgap", records=0, description="relative gap")
    lowerBound = Parameter(
        m,
        name="lowerBound",
        records=float("-inf"),
        description="global lower bound",
    )
    upperBound = Parameter(
        m,
        name="upperBound",
        records=float("inf"),
        description="global upper bound",
    )
    objMaster = Parameter(m, name="objMaster", records=0)
    objSub = Parameter(m, name="objSub", records=0)
    rtol = 0.001

    received.l[j] = 0

    for iteration, _ in iter.records.itertuples(index=False):
        objSub[...] = 0

        for scenario, _ in s.records.itertuples(index=False):
            demand[j] = ScenarioData[scenario, j]
            subproblem.solve()
            objSub[...] = objSub + ScenarioData[scenario, "prob"] * zsub.l
            cutconst[iteration] = cutconst[iteration] + ScenarioData[
                scenario, "prob"
            ] * Sum(j, market.m[j] * demand[j])
            cutcoeff[iteration, j] = (
                cutcoeff[iteration, j]
                + ScenarioData[scenario, "prob"] * selling.m[j]
            )

        itActive[iteration] = True

        if (
            lowerBound.records.values[0][0]
            < objMaster.records.values[0][0] + objSub.records.values[0][0]
        ):
            lowerBound[...] = objMaster + objSub

        rgap[...] = (upperBound - lowerBound) / (1 + gams_math.abs(upperBound))
        if rgap.records.values[0][0] < rtol:
            break

        masterproblem.solve()
        upperBound.setRecords(zmaster.records["level"])
        objMaster.setRecords(zmaster.records["level"] - theta.records["level"])

    print(
        "Master problem Objective Function Value: ",
        masterproblem.objective_value,
    )
    print("Subproblem Objective Function Value: ", subproblem.objective_value)


if __name__ == "__main__":
    main()