# 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 simpyimport numpy as npnp.random.seed(0)env = simpy.Environment()#let start with run of  5 iterationsenv.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_levelself.reorder_qty = reorder_qtyself.balance = 0self.num_ordered = 0self.inventory = self.reorder_qtyself.env = envself.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.inventoryself.num_ordered = ceil(self.num_ordered/self.reorder_qty)*self.reorder_qtyself.balance -= 50*self.num_orderedyield self.env.timeout(2.0)self.inventory += self.num_orderedself.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*interarrivalself.demand = self.generate_demand()if self.demand < self.inventory:self.balance += 100*self.demandself.inventory -= self.demand#print(f'customer comes I sold {self.demand} at time {round(self.env.now,2)}')else:self.balance += 100*self.inventoryself.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

`import numpy as npimport simpyimport matplotlib.pyplot as pltnp.random.seed(0)from math import ceilclass inventory_simulation:def __init__(self,env:simpy.Environment,reorder_level:int,reorder_qty:int) -> None:self.reorder_level = reorder_levelself.reorder_qty = reorder_qtyself.balance = 0self.num_ordered = 0self.inventory = self.reorder_qtyself.env = envself.obs_time = []self.inventory_level = []self.demand =0def 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.inventoryself.num_ordered = ceil(self.num_ordered/self.reorder_qty)*self.reorder_qtyself.balance -= 50*self.num_orderedyield self.env.timeout(2.0)self.inventory += self.num_orderedself.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*interarrivalself.demand = self.generate_demand()if self.demand < self.inventory:self.balance += 100*self.demandself.inventory -= self.demand#print(f'customer comes I sold {self.demand} at time {round(self.env.now,2)}')else:self.balance += 100*self.inventoryself.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()`