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.

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

blends = 1x2 cell
    {'Balanced'}    {'HighN'}

nutrients = rawNutrients.Properties.RowNames

nutrients = 3x1 cell
    {'N'}
    {'P'}
    {'K'}

raws = rawNutrients.Properties.VariableNames % Raw materials

raws = 1x6 cell
    {'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)

Solving problem using linprog.

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: 110
          algorithm: 'dual-simplex-highs'
    constrviolation: 9.0949e-13
            message: 'Optimal solution found.'
      firstorderopt: 1.1369e-13
             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=12×6 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    

Related Topics

• Problem-Based Optimization Workflow