Beginning Inventory simulation using SimPy

Production simulation for the industrial process has been around since the 1950s and is a powerful technique for testing business conditions, logic, and changes in processes. Traditionally simulations were done using spreadsheets tools like excel or statistical-based packages such as SPSS or Minitab. In this blog, we will use the “Discrete event simulation” technique to model and test inventory simulations or scenarios. Discrete event simulation or DES creates events at a specific interval within the given time frame and these are then processed subsequently using blocks of processes based on some business conditions, this difference to the continuous event simulation like using differential equations to solve rocket trajectories, circuit analysis and other physical systems operating in real-time. The important point is the state of the program for DES is observed at fixed regular time points and is build on discrete stochastics mathematics such as queueing theory(e.g. works of Agner Erlang). The purpose of the simulation is to understand the interactions and impacts of events, processes, the sequence of actions to model physical processes.

There are several paid options when it comes to discrete event simulations such as Rockwell’s Arena, Matlab, etc., however, we will use the python simpy package, an open-source discrete-event simulation package based on Python.

Let's build a simple shop scenario where we are selling large cookie packs at 100 dollars a piece while the purchase and delivery cost is 50 dollars. We will use an (r, Q) model, a well-known policy of modeling over inventory strategy such that we reorder stuff only when the inventory reaches r unit level and order Q quantity. Obviously, in a real-world scenario, the overall model is much more complicated involving multiple strategies to optimize costs and lead-times.

The details of the simulation we are going to build are as follows

1. generate customer arrival based on exponential distribution with a rate factor
2. generate a random demand of widgets by the customer based on a uniform distribution
3. Reorder is generated based on the difference of reorder level and current inventory level. I will try to round off my reorder using the ceiling function

Let set up the environment first. Env create a simulation object and until variable sets the run cycle

import simpy
import numpy as np
np.random.seed(0)
env = simpy.Environment()
#let start with run of  5 iterations
env.run(until=5.0)

SimPy is based on generators or processes which are invoked based on logic provided in the program. We will create a shop_runner generator or process

We will create a simulation class with reorder level, reorder quantity, starting cost balance, num of units ordered, inventory level, observation time, and demand.

We will create a simple reorder level reorder Quantity model such that reorder is triggered when inventory reaches the reorder limit.

class inventory_simulation:
def __init__(self,env:simpy.Environment,reorder_level:int,reorder_qty:int) -> None:
self.reorder_level = reorder_level
self.reorder_qty = reorder_qty
self.balance = 0
self.num_ordered = 0
self.inventory = self.reorder_qty
self.env = env
self.obs_time = []
self.inventory_level = []
self.demand =0

Let's add demand and customer interarrival period modeled randomly. Customers enter the shop randomly based on rate parameter exp. distribution

def generate_interarrival(self) -> np.array:
return np.random.exponential(1./5)
def generate_demand(self) -> np.array:
return np.random.randint(1,5)

Let's create a function to handle reorders as generated based on reorder level and EOQ / reorder quantity

def handle_order(self) -> None:
#print(f'at {round(self.env.now)} placed order for {self.num_ordered}')
self.num_ordered = self.reorder_level + 1 -self.inventory
self.num_ordered = ceil(self.num_ordered/self.reorder_qty)*self.reorder_qty
self.balance -= 50*self.num_ordered
yield self.env.timeout(2.0)
self.inventory += self.num_ordered
self.num_ordered = 0

Let's create a runner function to orchestrate different sub-processes of customer arrivals, demands, and reorders

def runner_setup(self):
while True:
interarrival = self.generate_interarrival()
yield self.env.timeout(interarrival)
self.balance -= self.inventory*2*interarrival
self.demand = self.generate_demand()
if self.demand < self.inventory:
self.balance += 100*self.demand
self.inventory -= self.demand
#print(f'customer comes I sold {self.demand} at time {round(self.env.now,2)}')
else:
self.balance += 100*self.inventory
self.inventory = 0
#print(f'{self.inventory} is out of stock at {round(self.env.now,2)}')
if self.inventory < self.reorder_level and self.num_ordered ==0:
self.env.process(self.handle_order())

Finally, let's create an observation function to record inventory levels

def observe(self):
while True:
self.obs_time.append(self.env.now)
self.inventory_level.append(self.inventory)
yield self.env.timeout(0.1)

The plot shows a typical inventory behavior where items are sold and restocked based on lead-times and policy models

No. of cookie biscuits sold vs time

import numpy as np
import simpy
import matplotlib.pyplot as plt
np.random.seed(0)
from math import ceil
class inventory_simulation:
def __init__(self,env:simpy.Environment,reorder_level:int,reorder_qty:int) -> None:
self.reorder_level = reorder_level
self.reorder_qty = reorder_qty
self.balance = 0
self.num_ordered = 0
self.inventory = self.reorder_qty
self.env = env
self.obs_time = []
self.inventory_level = []
self.demand =0
def handle_order(self) -> None:
#print(f'at {round(self.env.now)} placed order for {self.num_ordered}')
self.num_ordered = self.reorder_level + 1 -self.inventory
self.num_ordered = ceil(self.num_ordered/self.reorder_qty)*self.reorder_qty
self.balance -= 50*self.num_ordered
yield self.env.timeout(2.0)
self.inventory += self.num_ordered
self.num_ordered = 0
#print(f'at {round(self.env.now)} recieved order for {self.num_ordered}')
def generate_interarrival(self) -> np.array:
return np.random.exponential(1./5)
def generate_demand(self) -> np.array:
return np.random.randint(1,5)
def observe(self):
while True:
self.obs_time.append(self.env.now)
self.inventory_level.append(self.inventory)
yield self.env.timeout(0.1)
def runner_setup(self):
while True:
interarrival = self.generate_interarrival()
yield self.env.timeout(interarrival)
self.balance -= self.inventory*2*interarrival
self.demand = self.generate_demand()
if self.demand < self.inventory:
self.balance += 100*self.demand
self.inventory -= self.demand
#print(f'customer comes I sold {self.demand} at time {round(self.env.now,2)}')
else:
self.balance += 100*self.inventory
self.inventory = 0
#print(f'{self.inventory} is out of stock at {round(self.env.now,2)}')
if self.inventory < self.reorder_level and self.num_ordered ==0:
self.env.process(self.handle_order())
def plot_inventory(self):
plt.figure()
plt.step(self.obs_time,self.inventory_level)
plt.xlabel('Time')
plt.ylabel('No of biscuits in inventory')
def run(simulation:inventory_simulation,until:float):
simulation.env.process(simulation.runner_setup())
simulation.env.process(simulation.observe())
simulation.env.run(until=until)
s = inventory_simulation(simpy.Environment(),10,30)
run(s,8)
#print(s.inventory_level)
s.plot_inventory()