Этот пример показывает, как создать начальную точку для задачи оптимизации, которая назвала индексные переменные. Для именованных индексных переменных часто самый легкий способ задать начальную точку состоит в том, чтобы использовать функцию findindex.
Проблемой является проблема материально-технических ресурсов мультипериода, которая включает смешивающиеся сырые данные и усовершенствованные масла. Цель состоит в том, чтобы максимизировать прибыль, подвергающуюся различным ограничениям на производство и мощности материально-технических ресурсов и на "твердость" нефтяных смешений. Эта проблема взята от Уильямса [1].
Проблема включает два типа нефти сырых овощей и три типа необработанной неовощной нефти, которую производитель может совершенствовать в пищевое масло. Производитель может совершенствовать до 200 тонн растительных масел и до 250 тонн неовощных масел в месяц. Производитель может сохранить 1 000 тонн каждой необработанной нефти, которая выгодна, потому что стоимость покупательных необработанных масел зависит от месяца, а также типа нефти. Названная "твердость" качества сопоставлена с каждой нефтью. Твердость смешанной нефти является линейно взвешенной твердостью составляющих масел.
Из-за обработки ограничений производитель ограничивает количество масел, усовершенствованных в любом месяце к не больше, чем три. Кроме того, если нефть усовершенствована за месяц, по крайней мере 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)
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.046990e+05.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
28 0.11 1 -9.945833e+04 4.818992e+00
70 0.14 2 -9.998889e+04 4.262813e+00
114 0.17 3 -1.002787e+05 3.891142e+00
1086 0.43 3 -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: 3
numnodes: 1086
constrviolation: 1.5111e-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];
Удовлетворите остающиеся ограничения баланса запаса stockbal при помощи функции 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];Проверяйте, что начальная точка выполнима. Поскольку ограничения имеют различные размерности, устанавливают пару "имя-значение" UniformOutput cellfun на false при проверке infeasibilities.
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
Все infeasibilities являются нулем, который показывает, что начальная точка выполнима.
Повторно выполните проблему с помощью начальной точки.
[sol2,fval2,exitstatus2,output2] = solve(prob,x0,'options',opts)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.046990e+05.
Relative gap is 166.74%.
Branch and Bound:
nodes total num int integer relative
explored time (s) solution fval gap (%)
28 0.11 2 -9.945833e+04 4.818992e+00
66 0.14 3 -9.987222e+04 4.384608e+00
159 0.17 4 -9.996667e+04 4.120028e+00
166 0.18 5 -1.001917e+05 3.886209e+00
249 0.20 6 -1.002139e+05 2.927655e+00
408 0.26 7 -1.002787e+05 2.538777e+00
1107 0.46 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: 1107
constrviolation: 3.6096e-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 1107 branch-and-bound steps, but using no initial point took 1086 steps.
[1] Уильямс, Х. Пол. Построение моделей в Математическом программировании. Четвертый выпуск. Дж. Вайли, Чичестер, Англия. Проблема 12.1, "Еда Manufacture1". 1999.