%% Cutting Stock Problem: Problem-Based
% This example shows how to solve a cutting stock problem using linear programming with an integer linear programming subroutine. The example uses the problem-based approach. For the solver-based approach, see Cutting Stock Problem: Solver-Based.
%
%% Problem Overview
% A lumber mill starts with trees that have been trimmed to fixed-length logs. Specify the fixed log length.

logLength = 40;

% The mill then cuts the logs into fixed lengths suitable for further processing. The problem is how to make the cuts so that the mill satisfies a set of orders with the fewest logs.
% Specify these fixed lengths and the order quantities for the lengths.

lengthlist = [8; 12; 16; 20];
quantity = [90; 111; 55; 30];
nLengths = length(lengthlist);

% Assume that there is no material loss in making cuts, and no cost for cutting.
%
%% Linear Programming Formulation
% Several authors, including Ford and Fulkerson [1] and Gilmore and Gomory [2], suggest the following procedure, which you implement in the next section. A cutting pattern is a set of lengths to which a single log can be cut.
%
% Instead of generating every possible cutting pattern, it is more efficient to generate cutting patterns as the solution of a subproblem. Starting from a base set of cutting patterns, solve the linear programming problem of minimizing the number of logs used subject to the constraint that the cuts, using the existing patterns, satisfy the demands.
% After solving that problem, generate a new pattern by solving an integer linear programming subproblem. The subproblem is to find the best new pattern, meaning the number of cuts from each length in lengthlist that add up to no more than the total possible length logLength. The quantity to optimize is the reduced cost of the new pattern, which is one minus the sum of the Lagrange multipliers for the current solution times the new cutting pattern. If this quantity is negative, then bringing that pattern into the linear program will improve its objective. If not, then no better cutting pattern exists, and the patterns used so far give the optimal linear programming solution. The reason for this conclusion is exactly parallel to the reason for when to stop the primal simplex method: the method terminates when there is no variable with a negative reduced cost. The problem in this example terminates when there is no pattern with negative reduced cost. For details and an example, see Column generation algorithms and its references.
% After solving the linear programming problem in this way, you can have noninteger solutions. Therefore, solve the problem once more, using the generated patterns and constraining the variables to have integer values.
%
%% MATLAB Problem-Based Formulation
% A pattern, in this formulation, is a vector of integers containing the number of cuts of each length in lengthlist. Arrange a matrix named patterns to store the patterns, where each column in the matrix gives a pattern. For example,
% \textrm{patterns}=\left\lbrack \begin{array}{cc}
% 2 & 0\\
% 0 & 2\\
% 0 & 1\\
% 1 & 0
% \end{array}\right\rbrack \ldotp
% The first pattern (column) represents two cuts of length 8 and one cut of length 20. The second pattern represents two cuts of length 12 and one cut of length 16. Each is a feasible pattern because the total of the cuts is no more than logLength = 40.
% In this formulation, if x is a column vector of integers containing the number of times each pattern is used, then patterns*x is a column vector giving the number of cuts of each type. The constraint of meeting demand is patterns*x >= quantity. For example, using the previous patterns matrix, suppose that x=\left\lbrack \begin{array}{c}
% 45\\
% 56
% \end{array}\right\rbrack. (This x uses 101 logs.) Then
% \textrm{patterns}*x=\left\lbrack \begin{array}{c}
% 90\\
% 112\\
% 56\\
% 45
% \end{array}\right\rbrack ,
% which represents a feasible solution because the result exceeds the demand
% \textrm{quantity}=\left\lbrack \begin{array}{c}
% 90\\
% 111\\
% 55\\
% 30
% \end{array}\right\rbrack \ldotp
% To have an initial feasible cutting pattern, use the simplest patterns, which have just one length of cut. Use as many cuts of that length as feasible for the log.

patterns = diag(floor(logLength./lengthlist));
nPatterns = size(patterns,2);

% To generate new patterns from the existing ones based on the current Lagrange multipliers, solve a subproblem. Call the subproblem in a loop to generate patterns until no further improvement is found. The subproblem objective depends only on the current Lagrange multipliers. The variables are nonnegative integers representing the number of cuts of each length. The only constraint is that the sum of the lengths of the cuts in a pattern is no more than the log length.

subproblem = optimproblem();
cuts = optimvar('cuts', nLengths, 1, 'Type','integer','LowerBound',zeros(nLengths,1));
subproblem.Constraints = dot(lengthlist,cuts) <= logLength;

% To avoid unnecessary feedback from the solvers, set the Display option to 'off' for both the outer loop and the inner subproblem solution.

lpopts = optimoptions('linprog','Display','off');
ipopts = optimoptions('intlinprog',lpopts);

% Initialize the variables for the loop.

reducedCost = -inf;
reducedCostTolerance = -0.0001;
exitflag = 1;

% Call the loop.

while reducedCost < reducedCostTolerance && exitflag > 0
    logprob = optimproblem('Description','Cut Logs');     
    % Create variables representing the number of each pattern used
    x = optimvar('x', nPatterns, 1, 'LowerBound', 0);
    % The objective is the number of logs used
    logprob.Objective.logsUsed = sum(x);
    % The constraint is that the cuts satisfy the demand
    logprob.Constraints.Demand = patterns*x >= quantity;
    
    [values,nLogs,exitflag,~,lambda] = solve(logprob,'options',lpopts);
    
    if exitflag > 0
        fprintf('Using %g logs\n',nLogs);
        % Now generate a new pattern, if possible
        subproblem.Objective = 1.0 - dot(lambda.Constraints.Demand,cuts);
        [values,reducedCost,pexitflag] = solve(subproblem,'options',ipopts);
        newpattern = round(values.cuts);
        if double(pexitflag) > 0 && reducedCost < reducedCostTolerance
            patterns = [patterns newpattern];
            nPatterns = nPatterns + 1;
        end
    end
end

%
% You now have the solution of the linear programming problem. To complete the solution, solve the problem again with the final patterns, changing the solution variable x to the integer type. Also, compute the waste, which is the quantity of unused logs (in feet) for each pattern and for the problem as a whole.

if exitflag <= 0 
    disp('Error in column generation phase')
else
    x.Type = 'integer';
    [values,logsUsed,exitflag] = solve(logprob,'options',ipopts);
    if double(exitflag) > 0
        values.x = round(values.x); % in case some values were not exactly integers
        logsUsed = sum(values.x);
        fprintf('Optimal solution uses %g logs\n', logsUsed);
        totalwaste = sum((patterns*values.x - quantity).*lengthlist); % waste due to overproduction
        for j = 1:size(values.x)
            if values.x(j) > 0
                fprintf('Cut %g logs with pattern\n',values.x(j));
                for w = 1:size(patterns,1)
                    if patterns(w,j) > 0
                        fprintf('    %g cut(s) of length %d\n', patterns(w,j),lengthlist(w));
                    end
                end
                wastej = logLength - dot(patterns(:,j),lengthlist); % waste due to pattern inefficiency
                totalwaste = totalwaste + wastej;
            fprintf('    Waste of this pattern is %g\n',wastej);
            end
        end
        fprintf('Total waste in this problem is %g.\n',totalwaste);
    else 
        disp('Error in final optimization')
    end
end

% Part of the waste is due to overproduction, because the mill cuts one log into three 12-foot pieces, but uses only one. Part of the waste is due to pattern inefficiency, because the three 12-foot pieces are 4 feet short of the total length of 40 feet.
%
%% References
% [1] Ford, L. R., Jr. and D. R. Fulkerson. A Suggested Computation for Maximal Multi-Commodity Network Flows. Management Science 5, 1958, pp. 97-101.
% [2] Gilmore, P. C., and R. E. Gomory. A Linear Programming Approach to the Cutting Stock Problem--Part II. Operations Research 11, No. 6, 1963, pp. 863-888.
% Copyright 2012 The MathWorks, Inc.