apis/examples_matlab_control/cutstock_8m_source.html

function cutstock(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


    % instantiate Options and define parameters

    opt = ws.addOptions();

    cutstockData = ws.addDatabase('csdata');

    opt.setAllModelTypes('Cplex');

    opt.optCR = 0; % Solve to optimality

    maxpattern = 35;

    opt.defines('pmax', sprintf('%d', maxpattern));

    opt.defines('solveMasterAs', 'RMIP');


    % define input data

    d = containers.Map({'i1', 'i2', 'i3', 'i4'}, {97, 610, 395, 211});

    w = containers.Map({'i1', 'i2', 'i3', 'i4'}, {45, 36, 31, 14});

    r = 100;


    widths = cutstockData.addSet('i', 1, 'widths');

    rawWidth = cutstockData.addParameter('r', 0, 'raw width');

    demand = cutstockData.addParameter('d', 1, 'demand');

    width = cutstockData.addParameter('w', 1, 'width');


    for key = keys(d)

        widths.addRecord(key{1});

        rec = demand.addRecord(key{1});

        rec.value = d(key{1});

    end

    for key = keys(w)

        rec = width.addRecord(key{1});

        rec.value = w(key{1});

    end

    rec = rawWidth.addRecord();

    rec.value = 100;


    masterModel = {

    '$Title Cutting Stock - Master problem               '

    '                                                    '

    'Set  i    widths                                    '

    'Parameter                                           '

    ' w(i) width                                         '

    ' d(i) demand                                        '

    'Scalar                                              '

    ' r    raw width;                                    '

    '$gdxin csdata                                       '

    '$load i w d r                                       '

    '                                                    '

    '$if not set pmax $set pmax 1000                     '

    'Set  p        possible patterns  /1*%pmax%/         '

    '     pp(p)    dynamic subset of p                   '

    'Parameter                                           '

    ' aip(i,p) number of width i in pattern growing in p;'

    '                                                    '

    '* Master model                                      '

    'Variable xp(p)     patterns used                    '

    '         z         objective variable               '

    'Integer variable xp; xp.up(p) = sum(i, d(i));       '

    '                                                    '

    'Equation numpat    number of patterns used          '

    '     demand(i) meet demand;                         '

    '                                                    '

    'numpat..     z =e= sum(pp, xp(pp));                 '

    'demand(i)..  sum(pp, aip(i,pp)*xp(pp)) =g= d(i);    '

    '                                                    '

    'model master /numpat, demand/;                      '

    '                                                    '};

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


    subModel = {

    '$Title Cutting Stock - Pricing problem is a knapsack model '

    '                                                           '

    'Set  i    widths                                           '

    'Parameter                                                  '

    ' w(i) width;                                               '

    'Scalar                                                     '

    ' r    raw width;                                           '

    '                                                           '

    '$gdxin csdata                                              '

    '$load i w r                                                '

    '                                                           '

    'Parameter                                                  '

    ' demdual(i) duals of master demand constraint /#i eps/;    '

    '                                                           '

    'Variable  z, y(i) new pattern;                             '

    'Integer variable y; y.up(i) = ceil(r/w(i));                '

    '                                                           '

    'Equation defobj                                            '

    '         knapsack;                                         '

    '                                                           '

    'defobj..     z =e= 1 - sum(i, demdual(i)*y(i));            '

    'knapsack..   sum(i, w(i)*y(i)) =l= r;                      '

    'model pricing /defobj, knapsack/;                          '

    '                                                           '};

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


    % create initial checkpoint

    masterCP = ws.addCheckpoint();

    masterInitJob = ws.addJobFromString(masterModel);

    masterInitJob.run(opt, masterCP, cutstockData);


    masterJob = ws.addJobFromString('execute_load ''csdata'', aip, pp; solve master min z using %solveMasterAs%;', masterCP);


    pattern = cutstockData.addSet('pp', 1, 'pattern index');

    patternData = cutstockData.addParameter('aip', 2, 'pattern data');


    % Initial pattern: pattern i hold width i

    patternCount = 0;

    for key = keys(w)

        patternCount = patternCount + 1;

        patternCount_str = sprintf('%d', patternCount);

        pattern.addRecord(patternCount_str);

        rec = patternData.addRecord(key{1}, patternCount_str);

        rec.value = floor(r / w(key{1}));

    end


    % create model instance for sub job

    subCP = ws.addCheckpoint();

    ws.addJobFromString(subModel).run(opt, subCP, cutstockData);


    subMI = subCP.addModelInstance();


    % define modifier demdual

    demandDual = subMI.syncDB.addParameter('demdual', 1, 'dual of demand from master');

    subMI.instantiate('pricing min z using mip', opt, gams.control.Modifier(demandDual));


    % find new pattern

    patternAdded = true;

    while patternAdded

        masterJob.run(opt, masterCP, cutstockData);


        % Copy duals into gmssubMI.syncDB DB

        demandDual.clear();

        for d = masterJob.outDB.getEquation('demand').records

            rec = demandDual.addRecord(d{1}.key(1));

            rec.value = d{1}.marginal;

        end


        subMI.solve();


        z_level = subMI.syncDB.getVariable('z').record.level;

        if z_level < -0.00001

            if patternCount == maxpattern

                fprintf('Out of pattern. Increase maxpattern (currently [%d])\n', maxpattern);

                patternAdded = false;

            else

                fprintf('New pattern Value: %g\n', z_level);

                patternCount = patternCount + 1;

                s = pattern.addRecord(sprintf('%d', patternCount));

                for y = subMI.syncDB.getVariable('y').records

                    if y{1}.level > 0.5

                        rec = patternData.addRecord(y{1}.key(1), s.key(1));

                        rec.value = round(y{1}.level);

                    end

                end

            end

        else

            patternAdded = false;

        end

    end


    % solve final MIP

    opt.defines("solveMasterAs", "MIP");

    masterJob.run(opt, cutstockData);


    fprintf("Optimal Solution: %g\n", masterJob.outDB.getVariable("z").record.level);

    for xp = masterJob.outDB.getVariable("xp").records

        if xp{1}.level > 0.5

            fprintf("  pattern [%s] [%g] times: ", xp{1}.key(1), xp{1}.level);

            for aip = masterJob.outDB.getParameter("aip").records

                fprintf(" [%s] : [%g]\n", aip{1}.key(1), aip{1}.value);

            end

            fprintf("\n");

        end

    end


    % clean up of unmanaged resources

    cutstockData.dispose();

    subMI.dispose();

    opt.dispose();

end