Create Initial Point for Optimization with Named Index Variables

Open Live Script

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;

Create and Use Feasible Initial Point

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

issatisfied(prob,x0)

ans = logical
   1

All of the problem constraints are satisfied at x0, meaning the initial point is feasible.

Solve the problem using the initial point.

[sol,fval,exitstatus,output] = solve(prob,x0)

Solving problem using intlinprog.
Running HiGHS 1.7.0: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-01, 2e+02]
  Cost   [5e+00, 2e+02]
  Bound  [1e+00, 1e+03]
  RHS    [3e+00, 5e+02]
Assessing feasibility of MIP using primal feasibility and integrality tolerance of       1e-06
Solution has               num          max          sum
Col     infeasibilities      0            0            0
Integer infeasibilities      0            0            0
Row     infeasibilities      0            0            0
Row     residuals            0            0            0
Presolving model
126 rows, 115 cols, 368 nonzeros  0s
96 rows, 85 cols, 508 nonzeros  0s
96 rows, 85 cols, 508 nonzeros  0s

MIP start solution is feasible, objective value is 39250

Solving MIP model with:
   96 rows
   85 cols (30 binary, 0 integer, 0 implied int., 55 continuous)
   508 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   1074750         39250           2638.22%        0      0      0         0     0.0s
 R       0       0         0   0.00%   107512.962963   57078.571429      88.36%        0      0      0        80     0.0s
 C       0       0         0   0.00%   105440.917048   94225             11.90%      142     27      0       150     0.0s
 L       0       0         0   0.00%   104903.646216   99741.666667       5.18%      303     41      0       178     0.2s

10.0% inactive integer columns, restarting
Model after restart has 90 rows, 82 cols (27 bin., 0 int., 0 impl., 55 cont.), and 481 nonzeros

         0       0         0   0.00%   104744.514099   99741.666667       5.02%       20      0      0       633     0.2s
         0       0         0   0.00%   104744.514099   99741.666667       5.02%       20     20      0       685     0.2s
 L       0       0         0   0.00%   104643.726006   100278.703704      4.35%      113     22      0       986     0.6s
 B      18       0         6   3.52%   104643.726006   100278.703704      4.35%      132     22     27      3282     0.7s

Solving report
  Status            Optimal
  Primal bound      100278.703704
  Dual bound        100279.004373
  Gap               0.0003% (tolerance: 0.01%)
  Solution status   feasible
                    100278.703704 (objective)
                    0 (bound viol.)
                    2.16004991671e-14 (int. viol.)
                    0 (row viol.)
  Timing            1.16 (total)
                    0.01 (presolve)
                    0.00 (postsolve)
  Nodes             405
  LP iterations     7677 (total)
                    1240 (strong br.)
                    445 (separation)
                    2491 (heuristics)

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.RelativeGapTolerance = 0.0001. The intcon variables are integer within tolerance, options.ConstraintTolerance = 1e-06.

sol = struct with fields:
       buy: [6x5 double]
    induse: [6x5 double]
      sell: [6x5 double]
     store: [6x5 double]

fval = 
1.0028e+05

exitstatus = 
    OptimalSolution

output = struct with fields:
        relativegap: 2.9983e-06
        absolutegap: 0.3007
      numfeaspoints: 6
           numnodes: 405
    constrviolation: 1.8929e-11
          algorithm: 'highs'
            message: 'Optimal solution found....'
             solver: 'intlinprog'

Reference

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