Main Content

カッティング ストック問題: ソルバーベース

ライブ スクリプトを開く

この例は、整数線形計画法のサブルーチンと共に線形計画法を使用して、カッティング ストック問題を解く方法を説明します。この例では、ソルバーベースの最適化問題の設定のアプローチを使用します。問題ベースのアプローチについては、カッティング ストック問題: 問題ベースを参照してください。

問題概要

製材所は、木から固定長の木材を切り出すことから始めます。木材の固定長を指定します。

logLength = 40;

その後、以降の処理に適した固定長に木材を切ります。問題は、製材所が最小限の木材で一連の注文に対応できるように、どのように木材を切り出すかということです。

これらの固定長と長さごとの注文の数量を指定します。

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

切り出す際に素材の無駄は生じず、カッティングのコストもかからないものとします。

線形計画法の定式化

Ford と Fulkerson [1] および Gilmore と Gomory [2] を含む複数の著者が以下の手順を提案しています。これを次のセクションで実装します。カッティング パターンとは 1 本の木材から切り出すことができる長さのセットです。

すべての考えられるカッティング パターンを生成するより、部分問題の解としてカッティング パターンを生成する方が効率的です。カッティング パターンの基本セットから始め、既存のパターンを使用したカッティングで要求を満たすという制約の下で、使用する木材数を最小限にする線形計画問題を解きます。

その問題を解いた後に、整数線形計画法の部分問題を解いて新しいパターンを生成します。この部分問題では最適な新しいパターン、つまり、考えられる長さ合計の logLength を総和が上回らないように、lengthlist の各長さで切り出せる木材の数を求めます。最適化するための数量は、新しいパターンの被約費用であり、現在の解のラグランジュ乗数と新しいカッティング パターンの積の和を、1 から引いたものです。この数量が負の場合、そのパターンを線形計画問題に持ち込むと、その目的関数が改善されます。それ以外の場合は、より良いカッティング パターンは存在せず、今までに使用されたパターンが最適な線形計画問題の解を与えます。この結論の理由は、主シンプレックス法 (被約費用が負になる変数がないときに停止する方法) を停止する場合の理由と正確に対応します。この例の問題は、負の被約費用のパターンがないときに終了します。詳細と例については、列生成アルゴリズムとその参考文献を参照してください。

この方法で線形計画問題を解くと、非整数解が得られることがあります。したがって、生成されたパターンを使用し、変数が整数値を持つという制約の下で、もう一度問題を解きます。

MATLAB ソルバーベースでの定式化

この定式化では、パターンは lengthlist の各長さのカット (切り出した木材) 数を含む整数のベクトルです。パターンを保存するための行列 patterns を準備します。行列の各列が 1 つのパターンを表します。以下に例を示します。

patterns=[20020110].

最初のパターン (列) は長さ 8 の 2 本のカットと長さ 20 の 1 本のカットを表します。2 番目のパターンは、長さ 12 の 2 本のカットと長さ 16 の 1 本のカットを表します。それぞれは、カットの長さの合計が logLength = 40 を上回らないため、実行可能なパターンです。

この定式化では、x が、各パターンの使用回数を含む整数の列ベクトルである場合、patterns*x は各タイプのカット数を示す列ベクトルとなります。要求を満たすための制約は、patterns*x >= quantity です。たとえば、前出の行列 patterns を使用して、x=[4556] であると仮定します (この x は 101 本の木材を使用します)。この場合

patterns*x=[901125645],

これは要求を上回っているため、実行可能な解を表しています。

quantity=[901115530].

初期の実行可能なカッティング パターンを得るため、カットの長さが 1 つだけの最も簡単なパターンを使用します。その長さで、木材に対して実行可能なだけ多くのカット数を使用します。

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

現在のラグランジュ乗数に基づいて既存のパターンから新しいパターンを生成するには、部分問題を解きます。それ以上の改善が見つからなくなるまで、ループ内で部分問題を呼び出してパターンを生成します。部分問題の目的関数は、現在のラグランジュ乗数のみに依存します。変数は、各長さのカット数を表す非負の整数です。唯一の制約は、1 つのパターン内のカットの長さの合計が木材の長さを超えないことです。下限ベクトル lb2 と行列 A2、およびこれらの範囲と線形制約を表す b2 を作成します。

lb2 = zeros(nLengths,1);
A2 = lengthlist';
b2 = logLength;

ソルバーからの不要なフィードバックを避けるために、外側のループと内側の部分問題の解の両方に対して Display オプションを 'off' に設定します。

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

ループの変数を初期化します。

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

ループを呼び出します。

while reducedCost < reducedCostTolerance && exitflag > 0
    lb = zeros(nPatterns,1);
    f = lb + 1;
    A = -patterns;
    b = -quantity;
    
    [values,nLogs,exitflag,~,lambda] = linprog(f,A,b,[],[],lb,[],lpopts); 
    if exitflag > 0
        fprintf('Using %g logs\n',nLogs);
        % Now generate a new pattern, if possible
        f2 = -lambda.ineqlin;
        [values,reducedCost,pexitflag] = intlinprog(f2,1:nLengths,A2,b2,[],[],lb2,[],ipopts);
        reducedCost = 1 + reducedCost; % continue if this reducedCost is negative
        newpattern = round(values);
        if pexitflag > 0 && reducedCost < reducedCostTolerance
            patterns = [patterns newpattern];
            nPatterns = nPatterns + 1;
        end
    end
end
Using 97.5 logs
Using 92 logs
Using 89.9167 logs
Using 88.3 logs

これで、線形計画問題の解が得られます。解を完全なものにするために、すべての変数が整数である intlinprog を使用して、最終パターンでもう一度問題を解きます。また、各パターンと問題全体について、使用されない木材の数量 (フィート単位) として無駄を計算します。

if exitflag <= 0 
    disp('Error in column generation phase')
else
    [values,logsUsed,exitflag] = intlinprog(f,1:length(lb),A,b,[],[],lb,[],[],ipopts);
    if exitflag > 0
        values = round(values);
        logsUsed = round(logsUsed);
        fprintf('Optimal solution uses %g logs\n', logsUsed);
        totalwaste = sum((patterns*values - quantity).*lengthlist); % waste due to overproduction
        for j = 1:numel(values)
            if values(j) > 0
                fprintf('Cut %g logs with pattern\n',values(j));
                for w = 1:size(patterns,1)
                    if patterns(w,j) > 0
                        fprintf('    %d 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
Optimal solution uses 89 logs

Cut 18 logs with pattern

    5 cut(s) of length 8

    Waste of this pattern is 0

Cut 15 logs with pattern

    2 cut(s) of length 20

    Waste of this pattern is 0

Cut 1 logs with pattern

    2 cut(s) of length 8
    2 cut(s) of length 12

    Waste of this pattern is 0

Cut 55 logs with pattern

    2 cut(s) of length 12
    1 cut(s) of length 16

    Waste of this pattern is 0

Total waste in this problem is 28.

長さ、各長さの必要数、各長さの生成数を表すテーブルを作成します。

table(lengthlist,quantity,patterns*values,'VariableNames',["Length" "Number Required" "Number Produced"])
ans=4×3 table
    Length    Number Required    Number Produced
    ______    _______________    _______________

       8             90                 92      
      12            111                112      
      16             55                 55      
      20             30                 30

解では、長さ 8 が 2 つ、長さ 12 が 1 つ余分に生成され、合計で 28 フィートの過剰生成となります。

参考文献

[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–2022 The MathWorks, Inc.

関連するトピック