apis/examples_matlab_control/tsp_8m_source.html

function tsp(varargin)


    % check workspace info from arguments

    if nargin > 0

        wsInfo = gams.control.WorkspaceInfo();

        wsInfo.systemDirectory = varargin{1};

        ws = gams.control.Workspace(wsInfo);

    else

        ws = gams.control.Workspace();

    end


    model = {

    '$Title Traveling Salesman Problem                               '

    '$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;                                                 '};

    model = sprintf('%s\n', model{:});


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

    cutsPerRound = 10;

    % current cut

    curCut = 0;

    % cut limit for current model instance (cmax = curCut + cutsPerRound)

    cMax = 0;


    % database used to collect all generated cuts

    cutData =  ws.addDatabase();

    cc = cutData.addSet('cc', 1, '');

    acut = cutData.addParameter('acut', 3, '');

    rhscut = cutData.addParameter('rhscut', 1, '');


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

    n = {};

    subTour = {};


    while true

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

        if curCut >= cMax

            fprintf(',');

            cMax = curCut + cutsPerRound;

            cutData.export();


            % create the checkpoint

            tspJob = ws.addJobFromString(model);

            cp = ws.addCheckpoint();

            opt = ws.addOptions();

            opt.defines('nrcities', '20');

            opt.defines('cmax', int2str(cMax - 1));

            opt.defines('cutdata', cutData.name);


            % read input data from 'tsp.gdx'

            opt.defines('tspdata', ['"', fullfile(ws.systemDirectory, 'apifiles', 'Data', 'tsp.gdx'), '"']);

            tspJob.run(opt, cp);


            % fill the n list only once

            if numel(n) == 0

                for i = tspJob.outDB.getSet('i').records

                    n{end+1} = i{1}.key(1);

                end

            end


            % instantiate the model instance with modifiers miAcut and miRhscut

            mi = cp.addModelInstance();

            miAcut = mi.syncDB.addParameter('acut', 3, '');

            miRhscut = mi.syncDB.addParameter('rhscut', 1, '');

            modifiers = {gams.control.Modifier(miAcut), gams.control.Modifier(miRhscut)};

            mi.instantiate('assign use mip min z', modifiers);

        else

            fprintf('.');

        end


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

        mi.solve(gams.control.globals.SymbolUpdateType.ACCUMULATE);

        mi.syncDB.getParameter('acut').clear();

        mi.syncDB.getParameter('rhscut').clear();


        % collect graph information from the solution

        graph = containers.Map();

        notVisited = n;

        for i = n

            for j = n

                if mi.syncDB.getVariable('x').findRecord({i{1},j{1}}).level > 0.5

                    graph(i{1}) = j{1};

                end

            end

        end


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

        while numel(notVisited) ~= 0

            ii = notVisited{1};

            subTour = {ii};

            while ~strcmp(graph(ii), notVisited{1})

                ii = graph(ii);

                subTour{end+1} = ii;

            end


            for i = numel(notVisited):-1:1

                for j = 1:numel(subTour)

                    if strcmp(notVisited{i}, subTour{j})

                        notVisited(i) = [];

                        break

                    end

                end

            end


            % add the cuts to both databases (cutData, mi.SyncDB)

            for i = subTour

                for j = subTour

                    keys = {['c', int2str(curCut)], i{1}, j{1}};

                    rec = acut.addRecord(keys);

                    rec.value = 1;

                    rec = miAcut.addRecord(keys);

                    rec.value = 1;

                end

            end


            key = ['c', int2str(curCut)];

            rec = rhscut.addRecord(key);

            rec.value = numel(subTour) - 0.5;

            rec = miRhscut.addRecord(key);

            rec.value = numel(subTour) - 0.5;

            cc.addRecord(key);

            curCut = curCut + 1;

        end


        if numel(subTour) >= numel(n)

            break

        end

    end


    fprintf('\nz=%g\n', mi.syncDB.getVariable('z').record.level);

    fprintf('sub_tour: ');

    for i = subTour

        fprintf('%s -> ', i{1});

    end

    fprintf('%s\n', subTour{1});

end