MATLAB Examples

Create Multiperiod Inventory Model in Problem-Based Framework

This example shows how to create a multiperiod inventory model in the problem-based framework. The problem is to schedule production of fertilizer blends over a period of time using a variety of ingredients whose costs depend on time in a predictable way. Assume that you know in advance the demand for the fertilizers. The objective is to maximize profits while meeting demand, where the costs are for purchasing raw ingredients and for storing fertilizer over time. You can determine costs in advance by using futures or other contracts.

Contents

Fertilizers and Ingredients

Granular fertilizers have nutrients nitrogen (N), phosphorous (P), and potassium (K). You can blend the following raw materials to obtain fertilizer blends with the requisite nutrients.

load fertilizer
blends = blendDemand.Properties.VariableNames % Fertilizers to produce
nutrients = rawNutrients.Properties.RowNames
raws = rawNutrients.Properties.VariableNames % Raw materials
blends =

  1x2 cell array

    {'Balanced'}    {'HighN'}


nutrients =

  3x1 cell array

    {'N'}
    {'P'}
    {'K'}


raws =

  1x6 cell array

    {'MAP'}    {'Potash'}    {'AN'}    {'AS'}    {'TSP'}    {'Sand'}

The two fertilizer blends have the same nutrient requirements (10% N, 10% P, and 10% K by weight), except the "HighN" blend has an additional 10% N for a total of 20% N.

disp(blendNutrients) % Table is in percentage
         Balanced    HighN
         ________    _____

    N    10          20   
    P    10          10   
    K    10          10   

The raw materials have the following names and nutrient percentages by weight.

disp(rawNutrients) % Table is in percentage
         MAP    Potash    AN    AS    TSP    Sand
         ___    ______    __    __    ___    ____

    N    11      0        35    21     0     0   
    P    48      0         0     0    46     0   
    K     0     60         0     0     0     0   

The raw material Sand has no nutrient content. Sand dilutes other ingredients, if necessary, to obtain the requisite percentages of nutrients by weight.

Store the numbers of each of these quantities in variables.

nBlends = length(blends);
nRaws = length(raws);
nNutrients = length(nutrients);

Forecast Demand and Revenue

Assume that you know in advance the demand in weight (tons) for the two fertilizer blends for the time periods in the problem.

disp(blendDemand)
                 Balanced    HighN
                 ________    _____

    January      750         300  
    February     800         310  
    March        900         600  
    April        850         400  
    May          700         350  
    June         700         300  
    July         700         200  
    August       600         200  
    September    600         200  
    October      550         200  
    November     550         200  
    December     550         200  

You know the prices per ton at which you sell the fertilizer blends. These prices per ton do not depend on time.

disp(blendPrice)
    Balanced    HighN
    ________    _____

    400         550  

Prices of Raw Materials

Assume that you know in advance the prices in tons for the raw materials. These prices per ton depend on time according to the following table.

disp(rawCost)
                 MAP    Potash    AN     AS     TSP    Sand
                 ___    ______    ___    ___    ___    ____

    January      350    610       300    135    250    80  
    February     360    630       300    140    275    80  
    March        350    630       300    135    275    80  
    April        350    610       300    125    250    80  
    May          320    600       300    125    250    80  
    June         320    600       300    125    250    80  
    July         320    600       300    125    250    80  
    August       320    600       300    125    240    80  
    September    320    600       300    125    240    80  
    October      310    600       300    125    240    80  
    November     310    600       300    125    240    80  
    December     340    600       300    125    240    80  

Storage Cost

The cost for storing blended fertilizer applies per ton and per time period.

disp(inventoryCost)
    10

Capacity Constraints

You can store no more than inventoryCapacity tons of total fertilizer blends at any time period.

disp(inventoryCapacity)
        1000

You can produce a total of no more than productionCapacity tons in any time period.

disp(productionCapacity)
        1200

Connection Among Production, Sales, and Inventory

You start the schedule with a certain amount, or inventory, of fertilizer blends available. You have a certain target for this inventory at the final period. At each time period, the amount of fertilizer blend is the amount at the end of the previous time period, plus the amount produced, minus the amount sold. In other words, for times greater than 1:

inventory(time,product) = inventory(time-1,product) + production(time,product) - sales(time,product)

This equation implies that the inventory is counted at the end of the time period. The time periods in the problem are as follows.

months = blendDemand.Properties.RowNames;
nMonths = length(months);

The initial inventory affects the inventory at time 1 as follows.

inventory(1,product) = initialInventory(product) + production(1,product) - sales(1,product)

The initial inventory is in the data blendInventory{'Initial',:}. The final inventory is in the data blendInventory{'Final',:}.

Assume that unmet demand is lost. In other words, if you cannot fill all the orders in a time period, the excess orders do not carry over into the next time period.

Optimization Problem Formulation

The objective function for this problem is profit, which you want to maximize. Therefore, create a maximization problem in the problem-based framework.

inventoryProblem = optimproblem('ObjectiveSense','maximize');

The variables for the problem are the quantities of fertilizer blends that you make and sell each month, and the raw ingredients that you use to make those blends. The upper bound on sell is the demand, blendDemand, for each time period and each fertilizer blend.

make = optimvar('make',months,blends,'LowerBound',0);
sell = optimvar('sell',months,blends,'LowerBound',0,'UpperBound',blendDemand{months,blends});
use  = optimvar('use',months,raws,blends,'LowerBound',0);

Additionally, create a variable that represents the inventory at each time.

inventory = optimvar('inventory',months,blends,'LowerBound',0,'UpperBound',inventoryCapacity);

To calculate the objective function in terms of the problem variables, calculate the revenue and costs. The revenue is the amount you sell of each fertilizer blend times the price, added over all time periods and blends.

revenue = sum(blendPrice{1,:}.*sum(sell(months,blends),1));

The cost of ingredients is the cost for each ingredient used at each time, added over all time periods. Because the amount used at each time is separated into the amount used for each blend, also add over the blends.

blendsUsed = sum(use(months,raws,blends),3);
ingredientCost = sum(sum(rawCost{months,raws}.*blendsUsed));

The storage cost is the cost for storing the inventory over each time period, added over time and blends.

storageCost = inventoryCost*sum(inventory(:));

Now place the objective function into the Objective property of the problem by using dot notation.

inventoryProblem.Objective = revenue - ingredientCost - storageCost;

Problem Constraints

The problem has several constraints. First, express the inventory equation as a set of constraints on the problem variables.

materialBalance = optimconstr(months,blends);
timeAbove1 = months(2:end);
previousTime = months(1:end-1);
materialBalance(timeAbove1,:) = inventory(timeAbove1,:) == inventory(previousTime,:) +...
    make(timeAbove1,:) - sell(timeAbove1,:);
materialBalance(1,:) = inventory(1,:) == blendInventory{'Initial',:} +...
    make(1,:) - sell(1,:);

Express the constraint that the final inventory is fixed as well.

finalC = inventory(end,:) == blendInventory{'Final',:};

The total inventory at each time is bounded.

boundedInv = sum(inventory,2) <= inventoryCapacity;

You can produce a limited amount in each time period.

processLimit = sum(make,2) <= productionCapacity;

The amount that you produce each month of each blend is the amount of raw materials that you use. The squeeze function converts the sum from a nmonths-by-1-by- nblends array to a nmonths-by- nblends array.

rawMaterialUse = squeeze(sum(use(months,raws,blends),2)) == make(months,blends);

The nutrients in each blend must have the requisite values. In the following inner statement, the multiplication rawNutrients{n,raws}*use(m,raws,b)' adds the nutrient values at each time over the raw materials used.

blendNutrientsQuality = optimconstr(months,nutrients,blends);
for m = 1:nMonths
    for b = 1:nBlends
        for n = 1:nNutrients
            blendNutrientsQuality(m,n,b) = rawNutrients{n,raws}*use(m,raws,b)' == blendNutrients{n,b}*make(m,b);
        end
    end
end

Place the constraints into the problem.

inventoryProblem.Constraints.materialBalance = materialBalance;
inventoryProblem.Constraints.finalC = finalC;
inventoryProblem.Constraints.boundedInv = boundedInv;
inventoryProblem.Constraints.processLimit = processLimit;
inventoryProblem.Constraints.rawMaterialUse = rawMaterialUse;
inventoryProblem.Constraints.blendNutrientsQuality = blendNutrientsQuality;

Solve Problem

The problem formulation is complete. Solve the problem.

[sol,fval,exitflag,output] = solve(inventoryProblem)
Optimal solution found.


sol = 

  struct with fields:

    inventory: [12x2 double]
         make: [12x2 double]
         sell: [12x2 double]
          use: [12x6x2 double]


fval =

   2.2474e+06


exitflag = 

    OptimalSolution


output = 

  struct with fields:

         iterations: 164
    constrviolation: 5.4570e-12
            message: 'Optimal solution found.'
          algorithm: 'dual-simplex'
      firstorderopt: 6.5235e-12
             solver: 'linprog'

Display the results in tabular and graphical form.

if exitflag > 0
    fprintf('Profit: %g\n',fval);
    makeT = array2table(sol.make,'RowNames',months,'VariableNames',strcat('make',blends));
    sellT = array2table(sol.sell,'RowNames',months,'VariableNames',strcat('sell',blends));
    storeT = array2table(sol.inventory,'RowNames',months,'VariableNames',strcat('store',blends));
    productionPlanT = [makeT sellT storeT]
    figure
    subplot(3,1,1)
    bar(sol.make)
    legend('Balanced','HighN','Location','eastoutside')
    title('Amount Made')
    subplot(3,1,2)
    bar(sol.sell)
    legend('Balanced','HighN','Location','eastoutside')
    title('Amount Sold')
    subplot(3,1,3)
    bar(sol.inventory)
    legend('Balanced','HighN','Location','eastoutside')
    title('Amount Stored')
    xlabel('Time')
end
Profit: 2.24739e+06

productionPlanT =

  12x6 table

                 makeBalanced    makeHighN    sellBalanced    sellHighN    storeBalanced    storeHighN
                 ____________    _________    ____________    _________    _____________    __________

    January      1100            100          750             300          550                0       
    February      600            310          800             310          350                0       
    March         550            650          900             600            0               50       
    April         850            350          850             400            0                0       
    May           700            350          700             350            0                0       
    June          700            300          700             300            0                0       
    July          700            200          700             200            0                0       
    August        600            200          600             200            0                0       
    September     600            200          600             200            0                0       
    October       550            200          550             200            0                0       
    November      550            200          550             200            0                0       
    December      750            400          550             200          200              200