Traveling Salesman Problem - Four (TSP4,SEQ=180)#

tsp4.py tsp4.gdx

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_tsp4.html
## LICENSETYPE: Demo
## MODELTYPE: MIP
## DATAFILES: tsp4.gdx
## KEYWORDS: mixed integer linear programming, traveling salesman problem, iterative subtour elimination


Traveling Salesman Problem - Four (TSP4,SEQ=180)

This is the fourth problem in a series of traveling salesman
problems. Here we revisit TSP1 and generate smarter cuts.
The first relaxation is the same as in TSP1.


Kalvelagen, E, Model Building with GAMS. forthcoming

de Wetering, A V, private communication.
"""

from __future__ import annotations

import os
from pathlib import Path

import gamspy.math as gams_math
from gamspy import Alias
from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Options
from gamspy import Ord
from gamspy import Parameter
from gamspy import Set
from gamspy import Smax
from gamspy import Sum
from gamspy import Variable
from gamspy.exceptions import GamspyException


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
        load_from=str(Path(__file__).parent.absolute()) + "/tsp4.gdx",
    )

    # SETS
    ii, i = m.getSymbols(["ii", "i"])

    # ALIASES
    jj, j = m.getSymbols(["jj", "j"])

    # PARAMETER
    c = m.getSymbols(["c"])[0]

    # VARIABLES
    x = Variable(
        m,
        name="x",
        type="binary",
        domain=[ii, jj],
        description="decision variables - leg of trip",
    )

    # EQUATIONS
    rowsum = Equation(
        m, name="rowsum", domain=ii, description="leave each city only once"
    )
    colsum = Equation(
        m,
        name="colsum",
        domain=jj,
        description="arrive at each city only once",
    )

    # the assignment problem is a relaxation of the TSP
    rowsum[i] = Sum(j, x[i, j]) == 1
    colsum[j] = Sum(i, x[i, j]) == 1

    # Objective Function
    objective = Sum([i, j], c[i, j] * x[i, j])

    # exclude diagonal
    x.fx[ii, ii] = 0

    # For this algorithm we can try a larger subset of 12 cities.
    # SETS
    i.setRecords([f"i{ir}" for ir in range(1, 13)])

    # options. Make sure MIP solver finds global optima.
    assign = Model(
        m,
        name="assign",
        equations=[rowsum, colsum],
        problem="mip",
        sense="min",
        objective=objective,
    )
    assign.solve(options=Options(relative_optimality_gap=0))

    # find and display tours
    t = Set(
        m,
        name="t",
        records=[f"t{t}" for t in range(1, 18)],
        description="tours",
    )

    if len(t) < len(i):
        raise GamspyException("Set t is possibly too small")

    # SETS
    tour = Set(m, name="tour", domain=[i, j, t], description="subtours")
    visited = Set(
        m,
        name="visited",
        domain=i,
        description="flag whether a city is already visited",
    )

    # SINGLETON SETS
    fromi = Set(
        m,
        name="fromi",
        domain=i,
        is_singleton=True,
        description="contains always one element: the from city",
    )
    nextj = Set(
        m,
        name="nextj",
        domain=j,
        is_singleton=True,
        description="contains always one element: the to city",
    )
    tt = Set(
        m,
        name="tt",
        domain=t,
        is_singleton=True,
        description="contains always one element: the current subtour",
    )

    # ALIASES ##
    ix = Alias(m, name="ix", alias_with=i)

    # initialize
    fromi[i].where[Ord(i) == 1] = True  # turn first element on
    tt[t].where[Ord(t) == 1] = True  # turn first element on

    # subtour elimination by adding cuts
    # Set
    cc = Set(m, name="cc", records=[f"c{c}" for c in range(1, 1001)])

    # Alias
    ccc = Alias(m, name="ccc", alias_with=cc)  # we allow up to 1000 cuts

    # Set
    curcut = Set(
        m,
        name="curcut",
        domain=cc,
        description="current cut always one element",
    )
    allcuts = Set(m, name="allcuts", domain=cc, description="total cuts")

    # Parameter
    cutcoeff = Parameter(m, name="cutcoeff", domain=[cc, i, j])
    rhs = Parameter(m, name="rhs", domain=cc)
    nosubtours = Parameter(
        m, name="nosubtours", description="number of subtours"
    )

    # Equation
    cut = Equation(m, name="cut", domain=cc, description="dynamic cuts")

    cut[allcuts] = (
        Sum([i, j], cutcoeff[allcuts, i, j] * x[i, j]) <= rhs[allcuts]
    )

    tspcut = Model(
        m,
        name="tspcut",
        equations=[rowsum, colsum, cut],
        problem="mip",
        sense="MIN",
        objective=objective,
    )

    curcut[cc].where[Ord(cc) == 1] = True

    for ccc_loop in ccc.toList():
        #  initialize
        fromi[i].where[Ord(i) == 1] = True  # turn first element on
        tt[t].where[Ord(t) == 1] = True  # turn first element on
        tour[i, j, t] = False
        visited[i] = False

        for i_loop in i.toList():
            nextj[j].where[
                x.l[fromi, j] > 0.5
            ] = True  # check x.l(fromi,j) = 1 would be dangerous
            tour[fromi, nextj, tt] = True  # store in table
            visited[fromi] = True  # mark city 'fromi' as visited
            fromi[j] = nextj[j]

            if nextj.toList()[0] in visited.toList():  # if already visited...
                tt[t] = tt[t.lag(1)]
                for ix_loop in ix.toList():
                    if (
                        ix_loop in visited.toList()
                    ):  # find starting point of new subtour
                        continue
                    fromi[ix_loop] = True

        nosubtours[...] = Sum(t, gams_math.Max(0, Smax(tour[i, j, t], 1)))

        if nosubtours.toValue() == 1:  # done: no subtours
            break

        # introduce cut
        for idx, t_loop in enumerate(t.toList()):
            if idx + 1 > nosubtours.toValue():
                continue
            rhs[curcut] = -1

            for i_loop, j_loop, t_loop2, _ in tour.records.itertuples(
                index=False
            ):
                if t_loop2 != t_loop:
                    continue

                cutcoeff[curcut, i_loop, j_loop].where[
                    x.l[i_loop, j_loop] > 0.5
                ] = 1
                # not needed due to nature of assignment constraints
                #        cutcoeff(curcut, i, j)$(x.l[i,j] < 0.5) = -1
                rhs[curcut] = rhs[curcut] + 1
            allcuts[curcut] = True  # include this cut in set
            curcut[cc] = curcut[cc.lag(1)]

        tspcut.solve()
        print(
            "Cut: ",
            ccc_loop,
            "\t\t # of subtours remaining: ",
            nosubtours.toValue() - 1,
        )

    if nosubtours.toValue() != 1:
        raise GamspyException("Too many cuts needed")

    print("No subtours remaining. Solution found!!\n")
    print("x: \n", x.pivot().round(1))

    import math

    assert math.isclose(tspcut.objective_value, 39, rel_tol=0.001)


if __name__ == "__main__":
    main()