В этом примере показано, как создать начальную точку для задачи оптимизации, которая имеет именованные индексные переменные. Для именованных индексных переменных часто самый легкий способ задать начальную точку - это использовать findindex функция.
Проблема заключается в многопериодной инвентаризационной проблеме, которая включает смешивание необработанных и рафинированных масел. Цель заключается в максимальном увеличении прибыли с учетом различных ограничений на производственные и товарно-материальные мощности и на «твердость» нефтяных смесей. Эта задача взята из Williams [1].
Проблема связана с двумя видами необработанного растительного масла и тремя видами необработанного неведомого масла, которое производитель может перерабатывать в пищевое масло. Производитель может дорабатывать до 200 тонн растительных масел, и до 250 тонн нежестких масел в месяц. Производитель может хранить 1000 тонн каждого сырого масла, что выгодно, потому что стоимость приобретения сырого масла зависит как от месяца, так и от типа масла. Качество, называемое «твердость», связано с каждым маслом. Твердость смешанного масла является линейно взвешенной твердостью составных масел.
Из-за ограничений по переработке производитель ограничивает количество масел, очищенных за любой один месяц, не более чем тремя. Кроме того, если нефть очищается за месяц, то необходимо очистить не менее 20 тонн этой нефти. Наконец, если растительное масло рафинировано через месяц, то необходимо также рафинировать не вегетируемое масло 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];
Укажите стоимость необработанных масел в качестве этого массива. Каждая строка массива представляет стоимость сырого масла за месяц. Первая строка представляет затраты в январе, а последняя - затраты в июне.
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];Создайте следующие переменные задачи:
продать, количество каждой нефти продается каждый месяц
store, количество каждого масла, хранящегося в конце каждого месяца
buy, количество каждого масла, закупаемого каждый месяц
Кроме того, чтобы учесть ограничения на количество нефтей, очищенных и проданных каждый месяц, и минимальное количество произведенного, создайте вспомогательную двоичную переменную induse то есть 1 именно тогда, когда масло продается за месяц.
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) команда.
Задача имеет несколько ограничений, которые вы должны задать.
Количество каждой нефти, хранившейся в июне, составляет 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, :);
Если масло вообще очищается за месяц, по крайней мере 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;
Если растительное масло рафинировано, масло 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)
Solving problem using intlinprog.
LP: Optimal objective value is -1.075130e+05.
Cut Generation: Applied 41 Gomory cuts, 2 cover cuts,
1 mir cut, and 1 clique cut.
Lower bound is -1.047522e+05.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
19 0.12 1 -7.190185e+04 4.535003e+01
25 0.14 2 -7.205000e+04 4.505117e+01
25 0.14 3 -8.290000e+04 2.606702e+01
29 0.15 4 -8.775370e+04 1.909426e+01
40 0.15 5 -9.937870e+04 5.163142e+00
91 0.19 6 -9.987222e+04 4.643483e+00
105 0.19 7 -1.002139e+05 4.286718e+00
630 0.38 8 -1.002787e+05 2.283810e+00
1112 0.51 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 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: 8
numnodes: 1112
constrviolation: 1.1867e-11
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;
Удовлетворите initstockbal ограничение для первого месяца, основанное на начальном хранении 500 для каждого типа масла, и отсутствие покупки в первый месяц, и постоянное использование veg2 и non3.
x0.store(1,:) = [500 300 500 500 250];
Удовлетворите оставшимся ограничениям запаса запаса при помощи findindex функция.
[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
Все недопустимости равны нулю, что показывает, что начальная точка допустима.
Повторите задачу с помощью начальной точки.
[sol2,fval2,exitstatus2,output2] = solve(prob,x0,'options',opts)Solving problem using intlinprog.
LP: Optimal objective value is -1.075130e+05.
Cut Generation: Applied 41 Gomory cuts, 2 cover cuts,
1 mir cut, and 1 clique cut.
Lower bound is -1.047522e+05.
Relative gap is 166.88%.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
15 0.12 2 -7.190185e+04 4.535003e+01
21 0.14 3 -8.290000e+04 2.606702e+01
25 0.14 4 -8.775370e+04 1.909426e+01
30 0.14 5 -9.953889e+04 4.993907e+00
56 0.16 6 -9.982870e+04 4.689100e+00
138 0.20 7 -1.002787e+05 3.851839e+00
1114 0.49 7 -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: 7
numnodes: 1114
constrviolation: 1.7580e-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 1114 branch-and-bound steps, but using no initial point took 1112 steps.
[1] Уильямс, Х. Пол. Построение моделей в математическом программировании. Четвертое издание. Дж. Уайли, Чичестер, Англия. Проблема 12.1, «Продовольственная Manufacture1.» 1999 года.