apis/examples_python/tsp_8py_source.html

from itertools import product

import os

import sys

from gams import GamsWorkspace, GamsModifier, SymbolUpdateType


GAMS_MODEL = """

$if not set tspdata $abort 'tspdata not set'


Set

   ii    'cities'

   i(ii) 'subset of cities';


Alias (ii,jj), (i,j,k);


Parameter c(ii,jj) 'distance matrix';


$gdxIn %tspdata%

$load ii c

$gdxIn


$if not set nrCities $set nrCities 20

i(ii)$(ord(ii) < %nrCities%) = yes;


Variable

   x(ii,jj) 'decision variables - leg of trip'

   z        'objective variable';


Binary Variable x;

x.fx(ii,ii) = 0;


Equation

   objective  'total cost'

   rowsum(ii) 'leave each city only once'

   colsum(jj) 'arrive at each city only once';


* the assignment problem is a relaxation of the TSP

objective.. z =e= sum((i,j), c(i,j)*x(i,j));


rowsum(i).. sum(j, x(i,j)) =e= 1;


colsum(j).. sum(i, x(i,j)) =e= 1;


$if not set cmax $set cmax 2

set cut /c0*c%cmax%/;


Parameter

   acut(cut,ii,jj) 'cut constraint matrix'

   rhscut(cut)     'cut constraint rhs';


Equation sscut(cut) 'sub_tour elimination cut';

sscut(cut).. sum((i,j), acut(cut,i,j)*x(i,j)) =l= rhscut(cut);


set cc(cut) 'previous cuts';

cc(cut) = no;


$if set cutdata execute_load '%cutdata%', cc, acut, rhscut;


acut(cut,i,j)$(not cc(cut)) = eps;


rhscut(cut)$(not cc(cut)) = card(ii);


Model assign /all/;


option optcr=0;

"""


if __name__ == "__main__":

    sys_dir = sys.argv[1] if len(sys.argv) > 1 else None

    work_dir = sys.argv[2] if len(sys.argv) > 2 else None

    ws = GamsWorkspace(system_directory=sys_dir, working_directory=work_dir)


    # number of cuts that can be added to a model instance

    cuts_per_round = 10


    # current cut

    cur_cut = 0


    # cut limit for current model instance (cmax = cur_cut + cuts_per_round)

    cmax = 0


    # database used to collect all generated cuts

    cut_data = ws.add_database()

    cc = cut_data.add_set("cc", 1)

    acut = cut_data.add_parameter("acut", 3)

    rhscut = cut_data.add_parameter("rhscut", 1)


    n = []  # list of cities (i1, i2, i3, ...)

    while True:

        # create a new model instance if the cut limit is reached

        if cur_cut >= cmax:

            cmax = cur_cut + cuts_per_round

            cut_data.export()


            # create the checkpoint

            job = ws.add_job_from_string(GAMS_MODEL)

            cp = ws.add_checkpoint()

            opt = ws.add_options()

            opt.defines["nrcities"] = "20"

            opt.defines["cmax"] = str(cmax - 1)

            opt.defines["cutdata"] = cut_data.name

            opt.defines["tspdata"] = (

                '"'

                + os.path.join(ws.system_directory, "apifiles", "Data", "tsp.gdx")

                + '"'

            )

            job.run(gams_options=opt, checkpoint=cp)


            # fill the cities list only once

            if not n:

                for i in job.out_db["i"]:

                    n += i.keys


            # instantiate the model instance with modifiers mi_acut and mi_rhscut

            mi = cp.add_modelinstance()

            mi_acut = mi.sync_db.add_parameter("acut", 3)

            mi_rhscut = mi.sync_db.add_parameter("rhscut", 1)

            mi.instantiate(

                "assign use mip min z", [GamsModifier(mi_acut), GamsModifier(mi_rhscut)]

            )


        # solve model instance using update type accumulate and clear acut and rhscut afterwards

        mi.solve(SymbolUpdateType.Accumulate)

        mi.sync_db["acut"].clear()

        mi.sync_db["rhscut"].clear()


        # collect graph information from the solution

        graph = {}

        for i, j in product(n, n):

            if mi.sync_db["x"].find_record([i, j]).level > 0.5:

                graph[i] = j


        # find all sub tours and add the required cuts by modifying acut and rhscut

        not_visited = list(n)  # make copy

        while not_visited:

            i = not_visited[0]

            sub_tour = [i]

            while graph[i] != not_visited[0]:

                i = graph[i]

                sub_tour.append(i)

            not_visited = list(filter(lambda x: x not in sub_tour, not_visited))


            # add the cuts to both databases (cut_data, mi.sync_db)

            for i, j in product(sub_tour, sub_tour):

                acut.add_record([f"c{cur_cut}", i, j]).value = 1

                mi_acut.add_record([f"c{cur_cut}", i, j]).value = 1

            rhscut.add_record(f"c{cur_cut}").value = len(sub_tour) - 0.5

            mi_rhscut.add_record([f"c{cur_cut}"]).value = len(sub_tour) - 0.5

            cc.add_record(f"c{cur_cut}")

            cur_cut += 1


        # quit the loop if a sub tour contains all cities and print the solution

        if len(sub_tour) == len(n):

            print(f"z={mi.sync_db['z'].first_record().level}")

            print("sub_tour: ")

            for i in sub_tour:

                print(f"{i} -> ", end="")

            print(sub_tour[0])

            break