名前付きインデックス変数による最適化の初期点の作成
この例では、名前付きインデックス変数がある最適化問題の初期点の作成方法を説明します。名前付きインデックス変数の場合、通常、初期点を指定する最も簡単な方法は、関数 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;
問題を解く
初期点を使用した場合と使用しなかった場合の最終的な相違を確認するために、ヒューリスティックな方法は使用せずに legacy アルゴリズムを使用するようにオプションを設定します。その後、問題を解きます。
opts = optimoptions('intlinprog','Heuristics','none','Algorithm','legacy'); [sol1,fval1,exitstatus1,output1] = solve(prob,'options',opts)
Solving problem using intlinprog.
LP: Optimal objective value is 107512.962963.
Cut Generation: Applied 40 Gomory cuts, 2 clique cuts,
and 1 mir cut.
Upper bound is 103812.342726.
Cut Generation: Applied 19 Gomory cuts.
Upper bound is 102950.537141.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
9 0.25 1 7.920556e+04 2.997325e+01
14 0.25 2 9.068704e+04 1.351810e+01
34 0.27 3 9.982870e+04 3.075075e+00
64 0.31 4 1.002139e+05 2.678897e+00
64 0.31 5 1.002139e+05 2.678897e+00
94 0.32 6 1.002139e+05 2.652285e+00
136 0.34 7 1.002787e+05 1.757033e+00
651 0.47 8 1.002787e+05 0.000000e+00
651 0.47 8 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 intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.
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: 8
numnodes: 651
constrviolation: 4.8004e-11
algorithm: 'legacy'
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 の初期保管量、最初の月の購入なし、および veg2 と non3 の一定の使用量に基づいて、最初の月の制約 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 UniformOutput を false に設定します。
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)Solving problem using intlinprog.
LP: Optimal objective value is 107512.962963.
Cut Generation: Applied 40 Gomory cuts, 2 clique cuts,
and 1 mir cut.
Upper bound is 103812.342726.
Relative gap is 164.49%.
Cut Generation: Applied 19 Gomory cuts.
Upper bound is 102950.537141.
Relative gap is 162.29%.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
9 0.16 2 7.920556e+04 2.997325e+01
14 0.16 3 9.068704e+04 1.351810e+01
30 0.18 4 9.982870e+04 3.075075e+00
60 0.20 5 1.002139e+05 2.678897e+00
60 0.20 6 1.002139e+05 2.678897e+00
90 0.21 7 1.002139e+05 2.652285e+00
134 0.22 8 1.002787e+05 2.087195e+00
627 0.32 9 1.002787e+05 1.451133e-14
627 0.32 9 1.002787e+05 1.451133e-14
Optimal solution found.
Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance, options.IntegerTolerance = 1e-05.
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: 9
numnodes: 627
constrviolation: 3.4788e-11
algorithm: 'legacy'
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 627 branch-and-bound steps, but using no initial point took 651 steps.
参考文献
[1] Williams, H. Paul.Model Building in Mathematical Programming.Fourth edition.J. Wiley, Chichester, England.Problem 12.1, "Food Manufacture1."1999.
参考
関連するトピック
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)