function C = inventory(s,S,varargin)
% function C = inventory(s,S,varargin)
% The total cost of the (s,S)-inventory simulation
% - At s, order-up-to S
% - Inventory is controlled once a month
%
% ---------------------------------------------------------------
% INPUT
% s: Order level (1x1 int)
% S: Order-up-to level (1x1 in)
% varargin: Plot inventory level if 'plot' is given
% ---------------------------------------------------------------
% OUTPUT
% C: Total costs (1x1 double)
% ---------------------------------------------------------------
%
% Created: 2018-02-21 Heikki Puustinen
% Check if we wanted to plot
if ~isempty(varargin) && strcmp(varargin{1},'plot')
bPlot = true;
else
bPlot = false;
end
%% Initialization
y(1) = 60; % Initial inventory
Tmax = 120; % Total simulation time (months)
held = 0; % Number of held items during simulation
back = 0; % Number of backlogged items during simulation
c_f = 32; % Fixed order costs
c_v = 3; % Order cost per item
order_cost = 0; % Total order costs during simulation
hold_cost = 1; % Holding cost of items
short_cost = 5; % Shortage cost of items
mtd = 0.1; % Mean time between deliveries
dp = [1/6 3/6 5/6 1]; % Cumulative probabilities for size of demand
t(1) = 0; % Simulation time
tp = 0; % Time of previous event
D = 0; % Demand of items
td = 0; % Time of next delivery (items leave inventory)
tc = 1; % Time, when inventory is next controlled
ts = 0; % Time of next supply (items arrive at inventory)
oa = 0; % Order amount
% Time of first supply is first set to a large number (so that it doesn't
% happen before it's time is determined)
ts = Tmax+1;
% Time of first demand
td = exprnd(mtd);
%% Simulation
% Continue simulation while simulation time is less than the maximum time
while t(end) < Tmax
% The simulation has three kind of events:
% 1) Delivery is made out of the inventory
% 2) The inventory is controlled
% 3) The inventory is replenished
% Time of delivery is the next event, i.e., it has smallest time
if td <= ts && td <= tc
% Update simulation time to the time of current event time
t(end+1) = td;
% Update counters for held or backlogged items
% If there are items in inventory, add to held items the size of
% the inventory times the time interval between now and last event
if y(end) > 0
held = held + y(end) * (t(end) - tp);
elseif y(end) < 0
% If inventory is empty and items are backlogged, calculate the
% backlogged items similarly
back = back - y(end) * (t(end) - tp);
end
% Update time of previous event to now
tp = t(end);
% Generate random number from the uniform distribution in order to
% determine the demand.
u = rand;
% Loop through demand probabilities
for ii = 1:length(dp)
% If demand is less than ii:th element
if u < dp(ii)
% Set demand as the ii
D = ii;
% Break out from the loop as we found our demand
break;
end
end
% Update inventory level
y(end+1) = y(end) - D;
% Generate time of next demand
td = t(end) + exprnd(mtd);
elseif tc < td && tc < ts
% Inventory control is the next event, i.e., it has the smallest
% time
% Update simulation time to currenct event time
t(end+1) = tc;
% Update counters for held or backlogged items
if y(end) > 0
held = held + y(end) * (t(end) - tp);
elseif y(end) < 0
back = back - y(end) * (t(end) - tp);
end
% Update time of the previous event to now
tp = t(end);
% Check, if inventory level is under s
if y(end) < s
% Order amount (S-inventory)
oa = S - y(end);
% Order delay
ts = t(end) + unifrnd(0.5,1);
% Order costs (fixed + variable)
order_cost = order_cost + c_f + c_v * oa;
end
% Inventory level is not changed
y(end+1) = y(end);
% Update time of next control - at the start of next month
tc= t(end)+1;
elseif ts < td && ts < tc
% Order arrival (supply) is the next event
% Update Simulation time to current event time
t(end+1) = ts;
% Update counters for held or backlogged items
if y(end) > 0
held = held + y(end) * (t(end) - tp);
elseif y(end) < 0
back = back - y(end) * (t(end) - tp);
end
% Update time of previous event to now
tp = t(end);
% Update inventory level, i.e., add the order amount to inventory
y(end+1) = y(end) + oa;
% Reset time of supply
ts = Tmax+1;
% Reset order ammount
oa = 0;
end
% Update counters from the remaining simulation time
% in case all events are scheduled after Tmax
if min(td,min(ts,tc)) >= Tmax
t(end) = Tmax;
if y(end) > 0
held = held + y(end) * (t(end)-tp);
elseif y < 0
back = back - y(end) * (t(end)-tp);
end
end
end
% Calculate average monthly cost
C = held * hold_cost + back * short_cost + order_cost;
C = C / Tmax;
% Plot inventory level from the simulation
if bPlot
% Note that plotting (t,y) will connect the points with lines
% indicating that inventory would change linearly. How ever, the
% inventory changes in discrete time steps and in order to depict this,
% the following trick is made:
% Every data point is duplicated so that the size of the inventory is
% the same between consecutive events.
Y = zeros(2*length(y),1);
T = zeros(2*length(y)-2,1);
for ii = 1:length(y)
Y(2*ii-1) = y(ii);
Y(2*ii) = y(ii);
end
for ii=2:length(t)
T(2*(ii-2)+1) = t(ii);
T(2*(ii-2)+2) = t(ii);
end
T=[0; T; Tmax];
% Open figure etc.
h = figure('Name','MS-E2710, Exercise 2.2');
ax = axes('Parent',h);
xlabel(ax,'Time [months]')
ylabel(ax,'Inventory')
hold(ax,'on')
% Plot inventory level against time
plot(ax,T,Y,'b');
% Plot zero-level
plot(ax,[0 Tmax],[0 0],'r')
end