apis/examples_python/tsp_8py_source.html

from gams import *

from itertools import product

import os

import sys


def get_model_text():

    return '''

$Title Traveling Salesman Problem Instance with Python


$Ontext


The sub_tour elimination constraints are generated by a Python

script. The MIP is solved over and over, but GAMS have to

generate the model only after n cuts have been added.


$Offtext


$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


$if not set nrCities $set nrCities 20

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


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

          z         objective variable;

binary variable x; x.fx(ii,ii) = 0;


equations 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__":

    if len(sys.argv) > 1:

        ws = GamsWorkspace(system_directory = sys.argv[1])

    else:

        ws = GamsWorkspace()


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

    cuts_per_round = 10

    # current cut

    curcut = 0

    # cut limit for current model instance (cmax = curcut + 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, "")

    # list of cities (i1, i2, i3, ...)

    n = []


    # do-while loop

    while True:

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

        if curcut >= cmax:

            cmax = curcut + cuts_per_round

            cut_data.export();


            # create the checkpoint

            tsp_job = ws.add_job_from_string(get_model_text())

            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(os.path.dirname(os.path.realpath(__file__)), "..", "Data", "tsp.gdx") + '"'

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


            # fill the n list only once

            if not n:

                for i in tsp_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 = {}

        not_visited = list(n)

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

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

                graph[i] = j


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

        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(["c"+str(curcut), i, j]).value = 1

                mi_acut.add_record(["c"+str(curcut), i, j]).value = 1

            rhscut.add_record("c"+str(curcut)).value = len(sub_tour)-0.5

            mi_rhscut.add_record(["c"+str(curcut)]).value = len(sub_tour)-0.5

            cc.add_record("c"+str(curcut))

            curcut += 1


        # quit the do/while loop if a subtour contains all cities and print the solution

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

            print("z=" + str(mi.sync_db["z"].first_record().level))

            print("sub_tour: ")

            for i in sub_tour:

                print(i + " -> ", end="")

            print(sub_tour[0])

            break;