The purpose of this post is to develop a discrete-event simulation code using python. The objective of this code is to determine economic order quantity (EOQ) for inventory management. First the analytical solution of EOQ model will be used to validate python code. In the next step, the code will be applied to more realistic scenarios to maximize profit.

Regarding inventory control, economic order quantity (EOQ) is the optimum inventory order quantity that minimizes both holding cost and order cost.The model was initially developed by Ford W. Harris in 1913 and later popularized by R. H. Wilson and K. Andler.

Inventory holding costs originated from storing unsold or work-in-process (WIP) inventory. it includes the cost of damaged or spoiled inventory, storage space, handling, maintenance and insurance.

Ordering costs are fixed expenses to process an order. It includes cost to prepare requisition, purchase order, incoming inspection, received material handling, invoice processing and issuing payment. This is cost is fixed regardless of the order quantity.

Holding cost increases as order quantity increases as opposed to ordering cost that decreases with increased order quantity. The goal is to determine the optimum order quantity minimizing both costs. The total cost of a single inventory management consists of:

Purchasing Cost = unit price (P) X annual demand (D)

Order Setup Cost = demand (D) X setup (S)/ order quantity (Q)

Holding Cost = average stock (Q/2) X annual per unit holding cost (H)

So, Total Cost , TC = PD + DS/Q + HQ/2

An analytical solution can be derived for optimum order quantity, Q* where TC is minimum. This can be done in situation where the D, S and H are constants. The analytical solution will be used later in this article to validate the python code to comparing Q*. Taking first order derivative with respect to Q in cost equation:

0 = 0 + DS/Q2 + H/2

Rearranging the above equation:

Q* = 2DSH

Python code is provided below. First, define simulation class and necessary functions.

import numpy as np

class Simulation: 
    
    def __init__(self, reorder_point, order_quantity):
        self.inventory_level = order_quantity 
        self.order_size = 0 
        self.clock = 0.0  
        self.clock_order_rec = self.generate_interarrival()   
        self.clock_inventory_rec = float('inf')     
        self.sales = 0    
        self.purchase_cost = 0   
        self.order_cost = 0
        self.holding_cost = 0
        self.order_quantity = order_quantity
        self.reorder_point = reorder_point  
        
    def advance_time(self):
        clock_event = min(self.clock_order_rec, self.clock_inventory_rec)    
        self.holding_cost += self.inventory_level*1*(clock_event - self.clock)
        self.clock = clock_event 
        
        if self.clock_inventory_rec <= self.clock_order_rec:  
            self.handle_delivery_event()   
        else:
            self.handle_customer_event()
        
    def handle_customer_event(self):
        demand = self.generate_demand()   
        
        if self.inventory_level > demand:
            self.sales += 1000*demand   
            self.inventory_level -= demand   
        else:
            self.sales += 1000*self.inventory_level   
            self.inventory_level = 0    
            
        if self.inventory_level < self.reorder_point and self.order_size == 0: 
            self.order_size = self.order_quantity
            self.purchase_cost += 500*self.order_size  
            self.order_cost += 2000
            self.clock_inventory_rec = self.clock + 2 
            
        self.clock_order_rec = self.clock + self.generate_interarrival()
            
    def handle_delivery_event(self):
        self.inventory_level += self.order_size
        self.order_size = 0 
        self.clock_inventory_rec = float('inf')
    
    def generate_interarrival(self): 
        #return np.random.randint(1, 2) 
        return 1
    
    def generate_demand(self):
        #return np.random.randint(1, 2)  
        return 10

Next, build and run simulation steps.

np.random.seed(0)

obsTime = []
inventoryLevel = []

s = Simulation(20, 50)
        
while s.clock < 240:
    s.advance_time()  
    obsTime.append(s.clock)
    inventoryLevel.append(s.inventory_level)
            
    print (s.clock, s.inventory_level, s.order_cost, s.holding_cost)

Plot the results.

import matplotlib.pyplot as plt

plt.figure()
plt.step(obsTime, inventoryLevel, where = 'post')
plt.xlabel('Simulation time (days)')
plt.ylabel('Inventory level')