This file shows that the dynamic programming problem from Question 5 on the Final Exam (Deterministic Periodic Review Inventory Model) can be turned into a linear integer programming problem that can be solved with Pyomo/GLPK.
from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("glpk")
We start with setting up the index set and the parameters, and initialize a concrete Pyomo model.
I = [1, 2, 3, 4]
I0 = [0, 1, 2, 3, 4]
r = {1:1, 2:6, 3:2, 4:3}
K = 1
h = 0.1
max_production = 6
model = ConcreteModel()
In order to obtain a linear model, we treat $x_i$, the number of items produced in period $i$, $s_i$, the number of items in stock at the end of period $i$, and $k_i$, a boolean variable which indicates whether or not to set up production in period $i$ as independent decision variables. For convenience of notation, we allow index $i=0$ and define $s_0=0$, i.e., we are starting with zero items in the inventory.
model.x = Var(I, within=NonNegativeIntegers)
model.s = Var(I0, within=NonNegativeIntegers)
model.k = Var(I, within=Boolean)
model.init = Constraint(expr = model.s[0]==0)
The objective function is simply the sum of all setup and holding costs.
model.z = Objective(expr=sum(model.k[i]*K + model.s[i]*h for i in I), sense=minimize)
The supply rule constrains the decision variables such that inventory at the beginning of each period plus number items produced within the period balance the number of items requested within the period plus the inventory that remains at the end of the period.
def supply_rule (model,i):
return model.s[i-1] + model.x[i] == r[i] + model.s[i]
model.supply = Constraint(I, rule=supply_rule)
The number of items produced $x_i$ and the setup decision variable $k_i$ satisfy a surprisingly simple linear constraint. (If the problem does not have an explicit maximal production, max_production can be set to the total number of items to be produced!)
def limit_rule (model,i):
return model.x[i] <= max_production*model.k[i]
model.limit = Constraint(I, rule=limit_rule)
Now we solve and inspect the result. GLPK returns one of the solutions we had hand-derived via dynamic programming.
results = opt.solve(model)
model.x.get_values()
model.z.expr()