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 monthbuy
, 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.1s 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.6s 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.06 (total) 0.00 (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.
Copyright 2012–2024 The MathWorks, Inc.