名前付きインデックス変数による最適化の初期点の作成
この例では、名前付きインデックス変数がある最適化問題の初期点の作成方法を説明します。名前付きインデックス変数の場合、通常、初期点を指定する最も簡単な方法は、関数 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;
実行可能な初期点の作成と使用
正しい次元の初期点を作成します。
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: [6×5 double]
induse: [6×5 double]
store: [6×5 double]
sell: [6×5 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];初期点が実行可能であることを確認します。
issatisfied(prob,x0)
ans = logical
1
すべての問題の制約が x0 で満たされていて、初期点が実行可能であることを示しています。
初期点を使用して問題を解きます。
[sol,fval,exitstatus,output] = solve(prob,x0)
Solving problem using intlinprog.
Running HiGHS 1.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
Matrix [1e-01, 2e+02]
Cost [5e+00, 2e+02]
Bound [1e+00, 1e+03]
RHS [3e+00, 5e+02]
Assessing feasibility of MIP using primal feasibility and integrality tolerance of 1e-06
Solution has num max sum
Col infeasibilities 0 0 0
Integer infeasibilities 0 0 0
Row infeasibilities 0 0 0
Row residuals 0 0 0
Presolving model
126 rows, 115 cols, 368 nonzeros 0s
96 rows, 85 cols, 508 nonzeros 0s
96 rows, 85 cols, 508 nonzeros 0s
MIP start solution is feasible, objective value is 39250
Solving MIP model with:
96 rows
85 cols (30 binary, 0 integer, 0 implied int., 55 continuous)
508 nonzeros
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 1074750 39250 2638.22% 0 0 0 0 0.0s
R 0 0 0 0.00% 107512.962963 57078.571429 88.36% 0 0 0 80 0.0s
C 0 0 0 0.00% 105440.917048 94225 11.90% 142 27 0 150 0.0s
L 0 0 0 0.00% 104903.646216 99741.666667 5.18% 303 41 0 178 0.1s
10.0% inactive integer columns, restarting
Model after restart has 90 rows, 82 cols (27 bin., 0 int., 0 impl., 55 cont.), and 481 nonzeros
0 0 0 0.00% 104744.514099 99741.666667 5.02% 20 0 0 633 0.2s
0 0 0 0.00% 104744.514099 99741.666667 5.02% 20 20 0 685 0.2s
L 0 0 0 0.00% 104643.726006 100278.703704 4.35% 113 22 0 986 0.5s
B 18 0 6 3.52% 104643.726006 100278.703704 4.35% 132 22 27 3282 0.6s
Solving report
Status Optimal
Primal bound 100278.703704
Dual bound 100279.004373
Gap 0.0003% (tolerance: 0.01%)
Solution status feasible
100278.703704 (objective)
0 (bound viol.)
2.16004991671e-14 (int. viol.)
0 (row viol.)
Timing 1.04 (total)
0.01 (presolve)
0.00 (postsolve)
Nodes 405
LP iterations 7677 (total)
1240 (strong br.)
445 (separation)
2491 (heuristics)
Optimal solution found.
Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.
sol = struct with fields:
buy: [6×5 double]
induse: [6×5 double]
sell: [6×5 double]
store: [6×5 double]
fval = 1.0028e+05
exitstatus =
OptimalSolution
output = struct with fields:
relativegap: 2.9983e-06
absolutegap: 0.3007
numfeaspoints: 6
numnodes: 405
constrviolation: 1.8929e-11
algorithm: 'highs'
message: 'Optimal solution found.↵↵Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.'
solver: 'intlinprog'
参考文献
[1] Williams, H. Paul.Model Building in Mathematical Programming.Fourth edition.J. Wiley, Chichester, England.Problem 12.1, "Food Manufacture1."1999.
Copyright 2012–2024 The MathWorks, Inc.
