apis/examples_matlab_control/benders2stage_8m_source.html

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


    data = {

    'Sets                                                                        '

    'i factories                                   /f1*f3/                       '

    'j distribution centers                        /d1*d5/                       '

    '                                                                            '

    'Parameter                                                                   '

    'capacity(i) unit capacity at factories                                      '

    '             /f1 500, f2 450, f3 650/                                       '

    'demand(j)   unit demand at distribution centers                             '

    '             /d1 160, d2 120, d3 270, d4 325, d5 700 /                      '

    'prodcost    unit production cost                   /14/                     '

    'price       sales price                            /24/                     '

    'wastecost   cost of removal of overstocked products /4/                     '

    '                                                                            '

    'Table transcost(i,j) unit transportation cost                               '

    '      d1    d2    d3    d4    d5                                            '

    ' f1   2.49  5.21  3.76  4.85  2.07                                          '

    ' f2   1.46  2.54  1.83  1.86  4.76                                          '

    ' f3   3.26  3.08  2.60  3.76  4.45;                                         '

    '                                                                            '

    '$ifthen not set useBig                                                      '

    'Set                                                                         '

    ' s scenarios /lo,mid,hi/                                                    '

    '                                                                            '

    'Table ScenarioData(s,*) possible outcomes for demand plus probabilities     '

    '     d1  d2  d3  d4  d5  prob                                               '

    ' lo  150 100 250 300 600 0.25                                               '

    ' mid 160 120 270 325 700 0.50                                               '

    ' hi  170 135 300 350 800 0.25;                                              '

    '$else                                                                       '

    '$if not set nrScen $set nrScen 10                                           '

    'Set       s                 scenarios                        /s1*s%nrScen%/;'

    'parameter ScenarioData(s,*) possible outcomes for demand plus probabilities;'

    'option seed=1234;                                                           '

    'ScenarioData(s,''prob'') = 1/card(s);                                       '

    'ScenarioData(s,j)      = demand(j)*uniform(0.6,1.4);                        '

    '$endif                                                                      '};

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


    masterModel = {

    'Sets                                                            '

    'i factories                                                     '

    'j distribution centers                                          '

    '                                                                '

    'Parameter                                                       '

    '  capacity(i)    unit capacity at factories                     '

    '  prodcost       unit production cost                           '

    '  transcost(i,j) unit transportation cost                       '

    '                                                                '

    '$if not set datain $abort ''datain not set''                    '

    '$gdxin %datain%                                                 '

    '$load i j capacity prodcost transcost                           '

    '                                                                '

    '* Benders master problem                                        '

    '$if not set maxiter $set maxiter 25                             '

    'Set                                                             '

    '  iter             max Benders iterations /1*%maxiter%/         '

    '                                                                '

    'Parameter                                                       '

    '  cutconst(iter)   constants in optimality cuts                 '

    '  cutcoeff(iter,j) coefficients in optimality cuts              '

    '                                                                '

    'Variables                                                       '

    '  ship(i,j)        shipments                                    '

    '  product(i)       production                                   '

    '  received(j)      quantity sent to market                      '

    '  zmaster          objective variable of master problem         '

    '  theta            future profit                                '

    'Positive Variables ship;                                        '

    '                                                                '

    'Equations                                                       '

    '  masterobj        master objective function                    '

    '  production(i)    calculate production in each factory         '

    '  receive(j)       calculate quantity to be send to markets     '

    '  optcut(iter)     Benders optimality cuts;                     '

    '                                                                '

    'masterobj..                                                     '

    '  zmaster =e=  theta -sum((i,j), transcost(i,j)*ship(i,j))      '

    '                     - sum(i,prodcost*product(i));              '

    '                                                                '

    'receive(j)..       received(j) =e= sum(i, ship(i,j));           '

    '                                                                '

    'production(i)..    product(i) =e= sum(j, ship(i,j));            '

    'product.up(i) = capacity(i);                                    '

    '                                                                '

    'optcut(iter)..  theta =l= cutconst(iter) +                      '

    '                       sum(j, cutcoeff(iter,j)*received(j));    '

    '                                                                '

    'model masterproblem /all/;                                      '

    '                                                                '

    '* Initialize cut to be non-binding                              '

    'cutconst(iter) = 1e15;                                          '

    'cutcoeff(iter,j) = eps;                                         '};

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


    subModel = {

    'Sets                                                            '

    '  i factories                                                   '

    '  j distribution centers                                        '

    '                                                                '

    'Parameter                                                       '

    '  demand(j)   unit demand at distribution centers               '

    '  price       sales price                                       '

    '  wastecost   cost of removal of overstocked products           '

    '  received(j) first stage decision units received               '

    '                                                                '

    '$if not set datain $abort ''datain not set''                    '

    '$gdxin %datain%                                                 '

    '$load i j demand price wastecost                                '

    '                                                                '

    '* Benders'' subproblem                                          '

    '                                                                '

    'Variables                                                       '

    '  sales(j)         sales (actually sold)                        '

    '  waste(j)         overstocked products                         '

    '  zsub             objective variable of sub problem            '

    'Positive variables sales, waste                                 '

    '                                                                '

    'Equations                                                       '

    '  subobj           subproblem objective function                '

    '  selling(j)       part of received is sold                     '

    '  market(j)        upperbound on sales                          '

    ';                                                               '

    '                                                                '

    'subobj..                                                        '

    '  zsub =e= sum(j, price*sales(j)) - sum(j, wastecost*waste(j)); '

    '                                                                '

    'selling(j)..  sales(j) + waste(j) =e= received(j);              '

    '                                                                '

    'market(j)..   sales(j) =l= demand(j);                           '

    '                                                                '

    'model subproblem /subobj,selling,market/;                       '

    '                                                                '

    '* Initialize received                                           '

    'received(j) = demand(j);                                        '};

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


    dataJob = ws.addJobFromString(data);


    optData = ws.addOptions();

    optData.defines('useBig', '1');

    optData.defines('nrScen', '100');

    dataJob.run(optData);

    optData.dispose();


    scenarioData = dataJob.outDB.getParameter('ScenarioData');

    opt = ws.addOptions();

    opt.defines('datain', dataJob.outDB.name);

    maxiter = 40;

    opt.defines('maxiter', int2str(maxiter));

    opt.setAllModelTypes('cplex');


    cpMaster = ws.addCheckpoint();

    cpSub = ws.addCheckpoint();


    ws.addJobFromString(masterModel).run(opt, cpMaster, dataJob.outDB);


    masteri = cpMaster.addModelInstance();

    cutconst = masteri.syncDB.addParameter('cutconst', 1, 'Benders optimality cut constant');

    cutcoeff = masteri.syncDB.addParameter('cutcoeff', 2, 'Benders optimality coefficients');

    theta = masteri.syncDB.addVariable('theta', 0, gams.control.globals.VarType.FREE, ...

        'Future profit function variable');

    thetaFix = masteri.syncDB.addParameter('thetaFix', 0, '');

    masteri.instantiate('masterproblem max zmaster using lp', opt, {...

        gams.control.Modifier(cutconst), gams.control.Modifier(cutcoeff), ...

        gams.control.Modifier(theta, gams.control.globals.UpdateAction.FIXED, thetaFix)});


    ws.addJobFromString(subModel).run(opt, cpSub, dataJob.outDB);


    subi = cpSub.addModelInstance();

    received = subi.syncDB.addParameter('received', 1, 'units received from master');

    demand = subi.syncDB.addParameter('demand', 1, 'stochastic demand');

    subi.instantiate('subproblem max zsub using lp', opt, {...

        gams.control.Modifier(received), gams.control.Modifier(demand)});


    opt.dispose();


    lowerbound = -Inf;

    upperbound = Inf;

    objmaster = Inf;

    iter = 1;


    while true

        fprintf('Iteration: %d\n', iter);


        % Solve master

        if iter == 1 % fix theta for first iteration

            thetaFix.addRecord().value = 0;

        else

            thetaFix.clear();

        end

        masteri.solve(gams.control.globals.SymbolUpdateType.BASECASE);


        fprintf(' Master %s : obj=%g\n', masteri.modelStatus.select, masteri.syncDB.getVariable('zmaster').record.level);

        if iter > 1

            upperbound = masteri.syncDB.getVariable('zmaster').record.level;

        end

        objmaster = masteri.syncDB.getVariable('zmaster').record.level - theta.record.level;


        % Set received from master

        received.clear();

        for r = masteri.syncDB.getVariable('received').records

            received.addRecord(r{1}.keys).value = r{1}.level;

            cutcoeff.addRecord(int2str(iter), r{1}.key(1));

        end


        cutconst.addRecord(int2str(iter));

        objsub = 0.0;

        for s = dataJob.outDB.getSet('s').records

            demand.clear();

            for j = dataJob.outDB.getSet('j').records

                demand.addRecord(j{1}.keys).value = scenarioData.findRecord({s{1}.key(1), j{1}.key(1)}).value;

            end


            subi.solve(gams.control.globals.SymbolUpdateType.BASECASE);

            fprintf(' Sub %s : obj=%g\n', subi.modelStatus.select, subi.syncDB.getVariable('zsub').record.level);


            probability = scenarioData.findRecord({s{1}.key(1), 'prob'}).value;

            objsub = objsub + probability * subi.syncDB.getVariable('zsub').record.level;

            for j = dataJob.outDB.getSet('j').records

                keys = j{1}.keys;

                record = cutconst.findRecord(int2str(iter));

                record.value = record.value + probability * subi.syncDB.getEquation('market').findRecord(keys).marginal * demand.findRecord(keys).value;

                record = cutcoeff.findRecord({int2str(iter), j{1}.key(1)});

                record.value = record.value + probability * subi.syncDB.getEquation('selling').findRecord(keys).marginal;

            end

        end


        lowerbound = max(lowerbound, objmaster + objsub);

        iter = iter + 1;

        if iter > maxiter

            error('Benders out of iterations');

        end


        fprintf(' lowerbound: %g upperbound: %g objmaster: %g\n', lowerbound, upperbound, objmaster);


        if upperbound - lowerbound < 0.001 * (1 + abs(upperbound))

            break;

        end

    end


    masteri.dispose();

    subi.dispose();

end