Create Initial Point for Optimization with Named Index Variables

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].

Problem Description

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).

Enter Problem Data

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 Variables

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);

Create Objective

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.

Create Constraints

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;

Solve Problem

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.046990e+05.                                                    

Branch and Bound:

   nodes     total   num int        integer       relative                                          
explored  time (s)  solution           fval        gap (%)                                         
      28      0.26         1  -9.945833e+04   4.818992e+00                                          
      70      0.29         2  -9.998889e+04   4.262813e+00                                          
     106      0.32         3  -1.002787e+05   3.891142e+00                                          
    1111      0.64         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: 1111
    constrviolation: 1.5176e-11
            message: 'Optimal solution found....'
             solver: 'intlinprog'

Use Initial Point

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.046990e+05.                                                    
                   Relative gap is 166.74%.                                                        

Branch and Bound:

   nodes     total   num int        integer       relative                                          
explored  time (s)  solution           fval        gap (%)                                         
      28      0.27         2  -9.945833e+04   4.818992e+00                                          
      74      0.30         3  -9.987222e+04   4.384608e+00                                          
     120      0.33         4  -1.002787e+05   3.891142e+00                                          
    1124      0.69         4  -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: 4
           numnodes: 1124
    constrviolation: 4.3485e-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 1124 branch-and-bound steps,
but using no initial point took 1111 steps.

Reference

[1] Williams, H. Paul. Model Building in Mathematical Programming. Fourth edition. J. Wiley, Chichester, England. Problem 12.1, "Food Manufacture1." 1999.

See Also

|

Related Topics