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.

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

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

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

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

disp(inventoryCost)

10

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

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.

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;

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;

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: 162
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=`*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