This example shows how to create an initial point for an optimization problem that has named index variables. For named index variables, often the easiest way to specify an initial point is to use the findindex function.
The problem is a multiperiod inventory problem that involves blending raw and refined oils. The objective is to maximize profit subject to various constraints on production and inventory capacities and on the "hardness" of oil blends. This problem is taken from Williams [1].
The problem involves two types of raw vegetable oil and three types of raw nonvegetable oil that a manufacturer can refine into edible oil. The manufacturer can refine up to 200 tons of vegetable oils, and up to 250 tons of nonvegetable oils per month. The manufacturer can store 1000 tons of each raw oil, which is beneficial because the cost of purchasing raw oils depends on the month as well as the type of oil. A quality called "hardness" is associated with each oil. The hardness of blended oil is the linearly weighted hardness of the constituent oils.
Because of processing limitations, the manufacturer restricts the number of oils refined in any one month to no more than three. Also, if an oil is refined in a month, at least 20 tons of that oil must be refined. Finally, if a vegetable oil is refined in a month, then nonvegetable oil 3 must also be refined.
The revenue is a constant for each ton of oil sold. The costs are the cost of purchasing the oils, which varies by oil and month, and there is a fixed cost per ton of storing each oil for each month. There is no cost for refining an oil, but the manufacturer cannot store refined oil (it must be sold).
Create named index variables for the planning periods and oils.
months = {'January','February','March','April','May','June'};
oils = {'veg1','veg2','non1','non2','non3'};
vegoils = {'veg1','veg2'};
nonveg = {'non1','non2','non3'};Create variables with storage and usage parameters.
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
Specify the hardness of the raw oils as this vector.
h = [8.8,6.1,2,4.2,5.0];
Specify the costs of the raw oils as this array. Each row of the array represents the cost of the raw oils in a month. The first row represents the costs in January, and the last row represents the costs in June.
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];Create these problem variables:
sell, the quantity of each oil sold each month
store, the quantity of each oil stored at the end of each month
buy, the quantity of each oil purchased each month
Additionally, to account for constraints on the number of oils refined and sold each month and the minimum quantity produced, create an auxiliary binary variable induse that is 1 exactly when an oil is sold in a month.
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);
Name the total quantity of oil sold each month produce.
produce = sum(sell,2);
To create the objective function for the problem, calculate the revenue, and subtract the costs of purchasing and storing oils.
Create an optimization problem for maximization, and include the objective function as the Objective property.
prob = optimproblem('ObjectiveSense', 'maximize'); prob.Objective = sum(saleprice*produce) - sum(sum(costdata.*buy)) - sum(sum(storecost*store));
The objective expression is quite long. If you like, you can see it using the showexpr(prob.Objective) command.
The problem has several constraints that you need to set.
The quantity of each oil stored in June is 500. Set this constraint by using lower and upper bounds.
store('June', :).LowerBound = 500; store('June', :).UpperBound = 500;
The manufacturer cannot refine more than maxuseveg vegetable oil in any month. Set this and all subsequent constraints by using Expressions for Constraints and Equations.
vegoiluse = sell(:, vegoils); vegused = sum(vegoiluse, 2) <= maxuseveg;
The manufacturer cannot refine more than maxusenon nonvegetable oil any month.
nonvegoiluse = sell(:,nonveg); nonvegused = sum(nonvegoiluse,2) <= maxusenon;
The hardness of the blended oil must be from hmin through hmax.
hardmin = sum(repmat(h, 6, 1).*sell, 2) >= hmin*produce; hardmax = sum(repmat(h, 6, 1).*sell, 2) <= hmax*produce;
The amount of each oil stored at the end of the month is equal to the amount at the beginning of the month, plus the amount bought, minus the amount sold.
initstockbal = 500 + buy(1, :) == sell(1, :) + store(1, :); stockbal = store(1:5, :) + buy(2:6, :) == sell(2:6, :) + store(2:6, :);
If an oil is refined at all in a month, at least minuseraw of the oil must be refined and sold.
minuse = sell >= minuseraw*induse;
Ensure that the induse variable is 1 exactly when the corresponding oil is refined.
maxusev = sell(:, vegoils) <= maxuseveg*induse(:, vegoils); maxusenv = sell(:, nonveg) <= maxusenon*induse(:, nonveg);
The manufacturer can sell no more than maxnraw oils each month.
maxnuse = sum(induse, 2) <= maxnraw;
If a vegetable oil is refined, oil non3 must also be refined and sold.
deflogic1 = sum(induse(:,vegoils), 2) <= induse(:,'non3')*numel(vegoils);Include the constraint expressions in the problem.
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;
To show the eventual difference between using an initial point and not using one, set options to use no heuristics. Then solve the problem.
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.10 1 -7.190185e+04 4.535003e+01
25 0.11 2 -7.205000e+04 4.505117e+01
25 0.11 3 -8.290000e+04 2.606702e+01
29 0.11 4 -8.775370e+04 1.909426e+01
40 0.12 5 -9.937870e+04 5.163142e+00
91 0.14 6 -9.987222e+04 4.643483e+00
105 0.15 7 -1.002139e+05 4.286718e+00
630 0.30 8 -1.002787e+05 2.283810e+00
1112 0.41 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'
For this problem, using an initial point can save branch-and-bound iterations. Create an initial point of the correct dimensions.
x0.buy = zeros(size(buy)); x0.induse = zeros(size(induse)); x0.store = zeros(size(store)); x0.sell = zeros(size(sell));
Set the initial point to sell only vegetable oil veg2 and nonvegetable oil non3. To set this initial point appropriately, use the findindex function.
numMonths = size(induse,1);
[idxMonths,idxOils] = findindex(induse,1:numMonths,{'veg2','non3'});
x0.induse(idxMonths,idxOils) = 1;Satisfy the maximum vegetable and nonvegetable oil constraints.
x0.sell(:,idxOils) = repmat([200,250],numMonths,1)
x0 = struct with fields:
buy: [6x5 double]
induse: [6x5 double]
store: [6x5 double]
sell: [6x5 double]
Set the initial point to buy no oil the first month.
x0.buy(1,:) = 0;
Satisfy the initstockbal constraint for the first month based on the initial store of 500 for each oil type, and no purchase the first month, and constant usage of veg2 and non3.
x0.store(1,:) = [500 300 500 500 250];
Satisfy the remaining stock balance constraints stockbal by using the findindex function.
[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];Check that the initial point is feasible. Because the constraints have different dimensions, set the cellfun UniformOutput name-value pair to false when checking the 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
All of the infeasibilities are zero, which shows that the initial point is feasible.
Rerun the problem using the initial point.
[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.11 2 -7.190185e+04 4.535003e+01
21 0.12 3 -8.290000e+04 2.606702e+01
25 0.12 4 -8.775370e+04 1.909426e+01
30 0.12 5 -9.953889e+04 4.993907e+00
56 0.14 6 -9.982870e+04 4.689100e+00
138 0.17 7 -1.002787e+05 3.851839e+00
1114 0.41 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'
This time, solve took fewer branch-and-bound steps to find the solution.
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] Williams, H. Paul. Model Building in Mathematical Programming. Fourth edition. J. Wiley, Chichester, England. Problem 12.1, "Food Manufacture1." 1999.