В этом примере показано, как создать начальную точку для задачи оптимизации с именем index variables. Для именованных переменных индекса часто проще всего указать начальную точку, используя findindex функция.
Проблема заключается в многопериодном запасе, который включает смешивание сырых и рафинированных масел. Цель заключается в максимальном увеличении прибыли с учетом различных ограничений в отношении производственных и складских мощностей и «твердости» нефтяных смесей. Эта проблема взята у Уильямса [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 растительное масло в любой месяц. Задайте это и все последующие ограничения с помощью Выражения (Expressions) для ограничений и уравнений (Constraints and Equations).
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.