From b2d49b3d03c26c0f151928692ed05b09aabe1a89 Mon Sep 17 00:00:00 2001 From: webartifex Date: Fri, 31 Mar 2017 01:53:16 +0200 Subject: [PATCH] Create first version for final assignment --- inventory.py | 285 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 285 insertions(+) create mode 100644 inventory.py diff --git a/inventory.py b/inventory.py new file mode 100644 index 0000000..53f63ab --- /dev/null +++ b/inventory.py @@ -0,0 +1,285 @@ +""" +inventory.py +~~~~~~~~~~~~ + +This module solves an inventory decision problem. + +Run with "python inventory.py" + +Problem description: + + "Suppose I sell iPads. I plan on selling iPads for 10 weeks. Every week, + I place an order with Apple and the iPads arrive instantly. If I run out + of inventory, I agree to give customers a “rain check.” This means + items are backordered." +""" + +import random + +import numpy as np + +# Use the tqdm package to show a nice status bar. +try: + from tqdm import tqdm +except ImportError: + tqdm = None + + +class InventoryDP: + """A dynamic programming solution to an inventory decision problem. + + Parameters + ---------- + time_horizon : integer + Number of time periods to be simulated in the decision tree + holding_cost : float + Cost in Dollars to keep one item on stock for one period + price : float + Sales price per item used to model "rain check" costs + max_inventory : integer + Maximum number of units that can be kept on stock at any time + max_demand : integer + Maximum number of units that can be demanded in a time period + + Returns + ------- + self : InventoryDP + The instance is returned to allow method chaining. + """ + + def __init__(self, time_horizon=10, holding_cost=5, price=100, + max_inventory=100, max_demand=10): + self.time_horizon = time_horizon + self.holding_cost = holding_cost + self.price = price + self.max_inventory = max_inventory + self.max_demand = max_demand + self._tree_height = None + + def _initialize(self): + """Initialize the data structures that represent the decision tree""" + + # The decision tree is modeled implicitly by storing the nodes' data + # in two two-dimensional arrays of appropriate sizes. The first + # dimension represents the inventory level. Since the problem + # formulation allows for negative inventories (= backorders), the + # "center" index in the first dimension equals a zero inventory level. + self._tree_height = 2 * self.max_inventory + 1 + self._shape = (self._tree_height, self.time_horizon + 1) + self._values = np.zeros(self._shape, dtype=np.float64) + self._decisions = np.zeros(self._shape, dtype=np.float64) + + def recurse(self): + """Built the decision tree and solve it by backwards induction""" + + self._initialize() + + # The values in the terminal nodes are assumed to be 0. + # Since the values array is initialized to zeroes, + # nothing needs to be computed for the terminal nodes. + + # Initialize the progress bar. + if tqdm: + total = self.time_horizon * self._tree_height + pbar = tqdm(total=total) + + # Iterate chronologically backwards over time and over all states. + # This implicitly visits all of the tree's nodes. + for period in reversed(range(self.time_horizon)): + for state in range(self._tree_height): + + # Update the progress bar. + if tqdm: + pbar.update() + + # Since the state dimension of the decision tree includes + # negative values for the inventory (= backorders), the + # index variable needs to be adjusted. + current_inventory = state - self.max_inventory + + # Create a two-dimensional array to store the temporary values + # that go into the expexted value calculations further below. + # The first dimension is of a size equal to the potential + # demands whereas the second dimension is equal to the number + # of maximum possible order size (i.e., the maximum number such + # that the inventory does not exceed max_inventory in the edge + # case of zero demand). For example, when max_inventory = 100, + # current_inventory = 95, and max_demand = 10, the decision + # maker could order 15 items at most (if there is no demand). + # Since it is assumed that the demand realization occurs + # after the decision maker places and receives the order, the + # potential excess of inventory over max_inventory will be + # forfeited. + possible_orders = (self.max_inventory - current_inventory + + self.max_demand) + shape = (self.max_demand + 1, possible_orders) + local_values = np.zeros(shape, dtype=np.float64) + + # Iterate over the cartesian product of possible orders and + # demand realizations and update the state variables. + for order in range(possible_orders): + for demand in range(self.max_demand + 1): + + # Calculate the inventory for the consequent period. + new_inventory = current_inventory + order - demand + + # If the calculated inventory is above or below + # the max_inventory, the value is clipped and a + # penalty for either loosing a customer entirely + # or not being able to keep the excess on stock + # is introduced. + if new_inventory < -1 * self.max_inventory: + new_inventory = -1 * self.max_inventory + throw_away_penalty = self.price * 2 + elif new_inventory > self.max_inventory: + new_inventory = self.max_inventory + throw_away_penalty = self.price * 2 + else: + throw_away_penalty = 0 + + # Calculate the realized cost. + local_values[demand, order] = ( + self.holding_cost * max(0, new_inventory) + + self.price * max(0, -1 * new_inventory) + + self._values[new_inventory + self.max_inventory, + period + 1] + + throw_away_penalty + ) + + # Calculate the expected value for each possible decision. + expected_values = (local_values.sum(axis=0) / (self.max_demand + 1)) + + # Save the minimum expected value in the tree matrices. + self._values[state, period] = expected_values.min() + self._decisions[state, period] = expected_values.argmin() + + if tqdm: + pbar.close() + + return self + + def simulate(self, initial_inventory): + """Simulate the decision problem and calculate the actual costs. + + Parameters + ---------- + initial_inventory : integer + Level of inventory at the beginning of the time horizon. + Must be between -max_inventory and +max_inventory. + If negative, the decision problem starts with a backorder. + + Returns + ------- + (costs, orders, demands) : float, ndarray [time_horizon], + ndarray [time_horizon] + The return value is a tuple that consists of the total + costs and two arrays with the orders and realized + demands for each period. + + Raises + ------ + ValueError + If initial_inventory is not between -max_inventory and + +max_inventory. + RuntimeError + If InventoryDP.recurse() was not run before. + """ + + # Check the initial_inventory parameter. + if not (0 <= initial_inventory <= abs(self.max_inventory)): + raise ValueError('initial_inventory out of range') + + # Ensure that the dynamic program was run before. + if self._tree_height is None: + raise RuntimeError('Call InventoryDP.recurse() first') + + # Initialize the data structures for the final output. + shape = (self.time_horizon, ) + orders = np.zeros(shape, dtype=np.int32) + demands = np.zeros(shape, dtype=np.int32) + costs = 0 + + # Iterate over the time horizon in chronological order. + for period in range(self.time_horizon): + + # Determine the optimal order quantity. + state = initial_inventory + self.max_inventory + orders[period] = self._decisions[state, period] + + # Get a random demand realization. + demands[period] = random.randint(0, self.max_demand) + + # Calculate the inventory for the consequent period. + initial_inventory += orders[period] - demands[period] + + # If the calculated inventory is above or below + # the max_inventory, the value is clipped and a + # penalty for either loosing a customer entirely + # or not being able to keep the excess on stock + # is introduced. + if initial_inventory < -1 * self.max_inventory: + initial_inventory = -1 * self.max_inventory + throw_away_penalty = self.price * 2 + elif initial_inventory > self.max_inventory: + initial_inventory = self.max_inventory + throw_away_penalty = self.price * 2 + else: + throw_away_penalty = 0 + + # Accumulate the total incurred costs. + costs += ( + self.holding_cost * max(0, initial_inventory) + + self.price * max(0, -1 * initial_inventory) + + throw_away_penalty + ) + + return costs, orders, demands + + +# Create a list of experiments with varied parameter settings. +experiments = [ + { + 'time_horizon': 10, + 'holding_cost': 5, + 'price': 5, + 'max_inventory': 100, + 'max_demand': 10, + }, + { + 'time_horizon': 10, + 'holding_cost': 10, + 'price': 5, + 'max_inventory': 100, + 'max_demand': 10, + }, + { + 'time_horizon': 10, + 'holding_cost': 5, + 'price': 10, + 'max_inventory': 100, + 'max_demand': 10, + }, +] + + +# Run the experiments. +for i, parameters in enumerate(experiments): + print() + print('Experiment #{:d}'.format(i+1)) + print('==============') + print() + print(' Parameters: ', parameters) + print() + + dp = InventoryDP(**parameters).recurse() + + # Simulate two independent runs for each experiment. + for j in range(2): + costs, orders, demands = dp.simulate(initial_inventory=10) + + print() + print(' Simulation #{:d}'.format(j+1)) + for k in range(len(orders)): + print(' Period #{:d}: {:d} ordered and {:d} sold' + .format(k+1, orders[k], demands[k])) + print(' Total costs: {:.2f}'.format(costs))