Aircraft Allocation under uncertain Demand (AIRCRAF)#

aircraft.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_aircraft.html
## LICENSETYPE: Demo
## MODELTYPE: LP
## KEYWORDS: linear programming, aircraft managing, allocation problem


Aircraft Allocation under uncertain Demand (AIRCRAF)

The objective of this model is to allocate aircrafts to routes to maximize
the expected profit when traffic demand is uncertain. Two different
formulations are used, the delta and the lambda formulation.


Dantzig, G B, Chapter 28. 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 Alias
from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Ord
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)),
    )

    # Sets
    i = Set(
        m,
        name="i",
        records=["a", "b", "c", "d"],
        description="aircraft types and unassigned passengers",
    )
    j = Set(
        m,
        name="j",
        records=[f"route-{i}" for i in range(1, 6)],
        description="assigned and unassigned routes",
    )
    h = Set(
        m, name="h", records=list(range(1, 6)), description="demand states"
    )

    # Alias
    hp = Alias(m, name="hp", alias_with=h)

    # Parameters
    dd = Parameter(
        m,
        name="dd",
        domain=[j, h],
        records=np.array([
            [200, 220, 250, 270, 300],
            [50, 150, 0, 0, 0],
            [140, 160, 180, 200, 220],
            [10, 50, 80, 100, 340],
            [580, 600, 620, 0, 0],
        ]),
        description="demand distribution on route j",
    )

    lamda = Parameter(
        m,
        name="lamda",
        domain=[j, h],
        records=np.array([
            [0.2, 0.05, 0.35, 0.2, 0.2],
            [0.3, 0.7, 0, 0, 0],
            [0.1, 0.2, 0.4, 0.2, 0.1],
            [0.2, 0.2, 0.3, 0.2, 0.1],
            [0.1, 0.8, 0.1, 0, 0],
        ]),
        description="probability of demand state h on route j",
    )

    c = Parameter(
        m,
        name="c",
        domain=[i, j],
        records=np.array([
            [18, 21, 18, 16, 10],
            [0, 15, 16, 14, 9],
            [0, 10, 0, 9, 6],
            [17, 16, 17, 15, 10],
        ]),
        description="costs per aircraft (1000s)",
    )

    p = Parameter(
        m,
        name="p",
        domain=[i, j],
        records=np.array([
            [16, 15, 28, 23, 81],
            [0, 10, 14, 15, 57],
            [0, 5, 0, 7, 29],
            [9, 11, 22, 17, 55],
        ]),
        description="passenger capacity of aircraft i on route j",
    )

    aa = Parameter(
        m,
        name="aa",
        domain=i,
        records=np.array([10, 19, 25, 15]),
        description="aircraft availability",
    )
    k = Parameter(
        m,
        name="k",
        domain=j,
        records=np.array([13, 13, 7, 7, 1]),
        description="revenue lost (1000 per 100 bumped)",
    )
    ed = Parameter(m, name="ed", domain=j, description="expected demand")
    gamma = Parameter(
        m,
        name="gamma",
        domain=[j, h],
        description="probability of exceeding demand increment h on route j",
    )
    deltb = Parameter(
        m,
        name="deltb",
        domain=[j, h],
        description="incremental passenger load in demand states",
    )

    ed[j][...] = Sum(h, lamda[j, h] * dd[j, h])
    gamma[j, h] = Sum(hp.where[(Ord(hp) >= Ord(h))], lamda[j, hp])
    deltb[j, h] = dd[j, h] - dd[j, h.lag(1, "linear")].where[dd[j, h]]

    # Variables
    x = Variable(
        m,
        name="x",
        domain=[i, j],
        type="positive",
        description="number of aircraft type i assigned to route j",
    )
    y = Variable(
        m,
        name="y",
        domain=[j, h],
        type="positive",
        description="passengers actually carried",
    )
    b = Variable(
        m,
        name="b",
        domain=[j, h],
        type="positive",
        description="passengers bumped",
    )
    oc = Variable(m, name="oc", type="positive", description="operating cost")
    bc = Variable(m, name="bc", type="positive", description="bumping cost")

    # Equations
    ab = Equation(m, name="ab", domain=i, description="aircraft balance")
    db = Equation(m, name="db", domain=j, description="demand balance")
    yd = Equation(
        m,
        name="yd",
        domain=[j, h],
        description="definition of boarded passengers",
    )
    bd = Equation(
        m,
        name="bd",
        domain=[j, h],
        description="definition of bumped passengers",
    )
    ocd = Equation(m, name="ocd", description="operating cost definition")
    bcd1 = Equation(
        m, name="bcd1", description="bumping cost definition: version 1"
    )
    bcd2 = Equation(
        m, name="bcd2", description="bumping cost definition: version 2"
    )

    phi = oc + bc
    ab[i] = Sum(j, x[i, j]) <= aa[i]
    db[j] = Sum(i, p[i, j] * x[i, j]) >= Sum(h.where[deltb[j, h]], y[j, h])
    yd[j, h] = y[j, h] <= Sum(i, p[i, j] * x[i, j])
    bd[j, h] = b[j, h] == dd[j, h] - y[j, h]
    ocd[...] = oc == Sum([i, j], c[i, j] * x[i, j])
    bcd1[...] = bc == Sum(j, k[j] * (ed[j] - Sum(h, gamma[j, h] * y[j, h])))
    bcd2[...] = bc == Sum([j, h], k[j] * lamda[j, h] * b[j, h])

    alloc1 = Model(
        m,
        name="alloc1",
        equations=[ab, db, ocd, bcd1],
        problem="LP",
        sense=Sense.MIN,
        objective=phi,
    )
    alloc2 = Model(
        m,
        name="alloc2",
        equations=[ab, yd, bd, ocd, bcd2],
        problem="LP",
        sense=Sense.MIN,
        objective=phi,
    )

    y.up[j, h] = deltb[j, h]
    alloc1.solve()

    print("Number of passengers carried with limit:")
    print(y.pivot())

    y.up[j, h] = np.inf

    alloc2.solve()
    print()
    print("Number of passengers carried without limit:")
    print(y.pivot())

    import math

    assert math.isclose(alloc2.objective_value, 1566.042189, rel_tol=0.001)


if __name__ == "__main__":
    main()