ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

名前付きインデックス変数による最適化の初期点の作成

この例では、名前付きインデックス変数がある最適化問題の初期点の作成方法を説明します。名前付きインデックス変数の場合、通常、初期点を指定する最も簡単な方法は、関数 findindex を使用することです。

問題は、原料となる油と精製された油のブレンドを含む、複数期間の在庫問題です。目的は、生産能力と在庫能力、およびオイル ブレンドの "硬度" への制約下において、収益を最大化することです。この問題は、Williams [1] からの引用です。

問題の説明

問題には、メーカーの精製食用油の原料となる植物油 2 種類と非植物油 3 種類が含まれます。メーカーは、毎月、最大 200 トンの植物油と最大 250 トンの非植物油を精製できます。メーカーは、各原料油を 1000 トンずつ保管できます。これは、原料油の購入コストがその月と油の種類によって異なるため、有益です。"硬度" と呼ばれる品質が、各油に関連付けられています。ブレンドされた油の硬度は、成分となる油の硬度を線形に重み付けしたものになります。

処理の限界のため、メーカーは 1 か月に精製する油を 3 種類以下に制限しています。また、1 種類の油を精製した月は、その油を少なくとも 20 トンは精製しなければなりません。最後に、1 種類の植物油を精製した月は、3 番目の非植物油も精製しなければなりません。

収益は、販売された油の各トンに対する定数です。コストは油の購入コストであり、油と月に応じて異なります。また、月ごとに各油を保管するためのトンあたりの固定コストがあります。油の精製に関するコストはありませんが、メーカーは精製後の油を保管できません (販売しなければなりません)。

問題のデータの入力

計画期間と油の名前付きインデックス変数を作成します。

months = {'January','February','March','April','May','June'};
oils = {'veg1','veg2','non1','non2','non3'};
vegoils = {'veg1','veg2'};
nonveg = {'non1','non2','non3'};

ストレージと使用量のパラメーターで変数を作成します。

maxstore = 1000; % Maximum storage of each type of oil
maxuseveg = 200; % Maximum vegetable oil use
maxusenon = 250; % Maximum nonvegetable oil use
minuseraw = 20; % Minimum raw oil use
maxnraw = 3; % Maximum number of raw oils in a blend
saleprice = 150; % Sale price of refined and blended oil
storecost = 5; % Storage cost per period per oil quantity
stockend = 500; % Stock at beginning and end of problem
hmin = 3; % Minimum hardness of refined oil
hmax = 6; % Maximum hardness of refined oil

原料油の硬度を次のベクトルとして指定します。

h = [8.8,6.1,2,4.2,5.0];

原料油のコストを次の配列として指定します。各配列の各行は、月あたりの原料油のコストを表します。最初の行は 1 月のコストを表し、最終行は 6 月のコストを表しています。

costdata = [...
110 120 130 110 115
130 130 110  90 115
110 140 130 100  95
120 110 120 120 125
100 120 150 110 105
 90 100 140  80 135];

変数の作成

これらの問題の変数を作成します。

  • sell、各月に販売された各油の数量

  • store、各月の末に保管されていた各油の数量

  • buy、各月に購入された各油の数量

さらに、各月に精製および販売される油の種類の数と、最小生産量に対する制約を説明するため、1 種類の油が 1 か月に販売される場合に必ず 1 となる補助 2 値変数 induse を作成します。

sell = optimvar('sell', months, oils, 'LowerBound', 0);
buy = optimvar('buy', months, oils, 'LowerBound', 0);
store = optimvar('store', months, oils, 'LowerBound', 0, 'UpperBound', maxstore);

induse = optimvar('induse', months, oils, 'Type', 'integer', ...
    'LowerBound', 0, 'UpperBound', 1);

各月に販売された油の合計量を、produce と名付けます。

produce = sum(sell,2);

目的関数の作成

問題の目的関数を作成するには、収益を計算し、油の購入と保管のコストを減算します。

最大化を目的とする最適化問題を作成し、目的関数を Objective プロパティとして含めます。

prob = optimproblem('ObjectiveSense', 'maximize');
prob.Objective = sum(saleprice*produce) - sum(sum(costdata.*buy)) - sum(sum(storecost*store));

目的関数の式はかなり長くなります。必要に応じて、showexpr(prob.Objective) コマンドを使用して、それを確認することができます。

制約の作成

問題には、設定が必要な制約がいくつかあります。

6 月に保管された各油の量は 500 です。下限と上限を使用して、この制約を設定します。

store('June', :).LowerBound = 500;
store('June', :).UpperBound = 500;

メーカーは、どの月であっても、maxuseveg を超える植物油は精製できません。これ以降のすべての制約を、制約の式を使用して設定します。

vegoiluse = sell(:, vegoils);
vegused = sum(vegoiluse, 2) <= maxuseveg;

メーカーは、どの月であっても、maxusenon を超える非植物油は精製できません。

nonvegoiluse = sell(:,nonveg);
nonvegused = sum(nonvegoiluse,2) <= maxusenon;

ブレンドされた油の硬度は、hmin through hmax から取得しなければなりません。

hardmin = sum(repmat(h, 6, 1).*sell, 2) >= hmin*produce;
hardmax = sum(repmat(h, 6, 1).*sell, 2) <= hmax*produce;

月末に保管されていた各油の量は、月初の量に購入した油の量を加え、販売した量を引いた量に等しくなります。

initstockbal = 500 + buy(1, :) == sell(1, :) + store(1, :);
stockbal = store(1:5, :) + buy(2:6, :) == sell(2:6, :) + store(2:6, :);

1 か月の間にいずれかの油を精製する場合は、その油を少なくとも minuseraw 精製し、販売しなければなりません。

minuse = sell >= minuseraw*induse;

対応する油が精製される場合は、変数 induse が必ず 1 になるようにします。

maxusev = sell(:, vegoils) <= maxuseveg*induse(:, vegoils);
maxusenv = sell(:, nonveg) <= maxusenon*induse(:, nonveg);

メーカーが各月に販売できる油の量は、maxnraw 以下です。

maxnuse = sum(induse, 2) <= maxnraw;

1 種類の植物油を精製する場合は、油 non3 も精製および販売しなければなりません。

deflogic1 = sum(induse(:,vegoils), 2) <= induse(:,'non3')*numel(vegoils);

問題に制約式を含めます。

prob.Constraints.vegused = vegused;
prob.Constraints.nonvegused = nonvegused;
prob.Constraints.hardmin = hardmin;
prob.Constraints.hardmax = hardmax;
prob.Constraints.initstockbal = initstockbal;
prob.Constraints.stockbal = stockbal;
prob.Constraints.minuse = minuse;
prob.Constraints.maxusev = maxusev;
prob.Constraints.maxusenv = maxusenv;
prob.Constraints.maxnuse = maxnuse;
prob.Constraints.deflogic1 = deflogic1;

問題を解く

初期点を使用した場合と使用しなかった場合の最終的な相違を確認するために、ヒューリスティックな方法を使用しないようにオプションを設定します。その後、問題を解きます。

opts = optimoptions('intlinprog','Heuristics','none');
[sol1,fval1,exitstatus1,output1] = solve(prob,'options',opts)
LP:                Optimal objective value is -1.075130e+05.                                        

Cut Generation:    Applied 26 Gomory cuts, 4 clique cuts,                                           
                   and 1 cover cut.                                                                 
                   Lower bound is -1.056500e+05.                                                    

Branch and Bound:

   nodes     total   num int        integer       relative                                          
explored  time (s)  solution           fval        gap (%)                                         
      17      0.10         1  -9.755000e+04   8.170810e+00                                          
      20      0.11         2  -9.982870e+04   5.701713e+00                                          
     180      0.17         3  -9.996667e+04   5.214782e+00                                          
     206      0.19         4  -1.001917e+05   4.978504e+00                                          
     387      0.26         5  -1.002787e+05   2.821848e+00                                          
    1325      0.57         5  -1.002787e+05   0.000000e+00                                          

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the
optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon
variables are integer within tolerance, options.IntegerTolerance = 1e-05 (the
default value).
sol1 = struct with fields:
       buy: [6x5 double]
    induse: [6x5 double]
      sell: [6x5 double]
     store: [6x5 double]

fval1 = 1.0028e+05
exitstatus1 = 
    OptimalSolution

output1 = struct with fields:
        relativegap: 0
        absolutegap: 0
      numfeaspoints: 5
           numnodes: 1325
    constrviolation: 1.2223e-12
            message: 'Optimal solution found....'
             solver: 'intlinprog'

初期点の使用

この問題の場合、初期点を使用すると、分枝限定の反復を節減できます。正しい次元の初期点を作成します。

x0.buy = zeros(size(buy));
x0.induse = zeros(size(induse));
x0.store = zeros(size(store));
x0.sell = zeros(size(sell));

植物油 veg2 と非植物油 non3 のみを販売するための初期点を設定します。この初期点を適切に設定するには、関数 findindex を使用します。

numMonths = size(induse,1);
[idxMonths,idxOils] = findindex(induse,1:numMonths,{'veg2','non3'});
x0.induse(idxMonths,idxOils) = 1;

植物油と非植物油の最大の制約をそれぞれ満たします。

x0.sell(:,idxOils) = repmat([200,250],numMonths,1)
x0 = struct with fields:
       buy: [6x5 double]
    induse: [6x5 double]
     store: [6x5 double]
      sell: [6x5 double]

最初の月に油を購入しないように初期点を設定します。

x0.buy(1,:) = 0;

油の種類ごとの 500 の初期保管量、最初の月の購入なし、および veg2non3 の一定の使用量に基づいて、最初の月の制約 initstockbal を満たします。

x0.store(1,:) = [500 300 500 500 250];

関数 findindex を使用して、残りの在庫残量の制約 stockbal を満たします。

[idxMonths,idxOils] = findindex(store,2:6,{'veg2'});
x0.store(idxMonths,idxOils) = [100;0;0;0;500];

[idxMonths,idxOils] = findindex(store,2:6,{'veg1','non1','non2'});
x0.store(idxMonths,idxOils) =  500;

[idxMonths,idxOils] = findindex(store,2:6,{'non3'});
x0.store(idxMonths,idxOils) = [0;0;0;0;500];

[idxMonths,idxOils] = findindex(buy,2:6,{'veg2'});
x0.buy(idxMonths,idxOils) = [0;100;200;200;700];

[idxMonths,idxOils] = findindex(buy,2:6,{'non3'});
x0.buy(idxMonths,idxOils) = [0;250;250;250;750];

初期点が実行可能であることを確認します。制約にはさまざまな次元があるため、実行不可能性を確認するときには、名前と値のペア cellfun UniformOutputfalse に設定します。

inf{1} = infeasibility(vegused,x0);
inf{2} = infeasibility(nonvegused,x0);
inf{3} = infeasibility(hardmin,x0);
inf{4} = infeasibility(hardmax,x0);
inf{5} = infeasibility(initstockbal,x0);
inf{6} = infeasibility(stockbal,x0);
inf{7} = infeasibility(minuse,x0);
inf{8} = infeasibility(maxusev,x0);
inf{9} = infeasibility(maxusenv,x0);
inf{10} = infeasibility(maxnuse,x0);
inf{11} = infeasibility(deflogic1,x0);
allinfeas = cellfun(@max,inf,'UniformOutput',false);
anyinfeas = cellfun(@max,allinfeas);
disp(anyinfeas)
     0     0     0     0     0     0     0     0     0     0     0

すべての実行不可能性が 0 であり、初期点が実行可能であることを示しています。

初期点を使用して問題を再実行します。

[sol2,fval2,exitstatus2,output2] = solve(prob,x0,'options',opts)
LP:                Optimal objective value is -1.075130e+05.                                        

Cut Generation:    Applied 26 Gomory cuts, 4 clique cuts,                                           
                   and 1 cover cut.                                                                 
                   Lower bound is -1.056500e+05.                                                    
                   Relative gap is 169.17%.                                                        

Cut Generation:    Applied 11 Gomory cuts.                                                          
                   Lower bound is -1.052476e+05.                                                    
                   Relative gap is 168.14%.                                                        

Branch and Bound:

   nodes     total   num int        integer       relative                                          
explored  time (s)  solution           fval        gap (%)                                         
      30      0.17         2  -1.002787e+05   4.933132e+00                                          
    1490      0.58         2  -1.002787e+05   0.000000e+00                                          

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the
optimal value, options.AbsoluteGapTolerance = 0 (the default value). The intcon
variables are integer within tolerance, options.IntegerTolerance = 1e-05 (the
default value).
sol2 = struct with fields:
       buy: [6x5 double]
    induse: [6x5 double]
      sell: [6x5 double]
     store: [6x5 double]

fval2 = 1.0028e+05
exitstatus2 = 
    OptimalSolution

output2 = struct with fields:
        relativegap: 0
        absolutegap: 0
      numfeaspoints: 2
           numnodes: 1490
    constrviolation: 9.3773e-12
            message: 'Optimal solution found....'
             solver: 'intlinprog'

今度は、solve はより少ない分枝限定ステップで解を求めました。

fprintf(['Using the initial point took %d branch-and-bound steps,\nbut ',...
    'using no initial point took %d steps.'],output2.numnodes,output1.numnodes)
Using the initial point took 1490 branch-and-bound steps,
but using no initial point took 1325 steps.

参考文献

[1] Williams, H. Paul.Model Building in Mathematical Programming.Fourth edition.J. Wiley, Chichester, England.Problem 12.1, "Food Manufacture1."1999.

参考

|

関連するトピック