3.4 Putting it all together with the lot sizing example: (Hart et al., 2017)

We will now write a complete model from scratch using a well-known multi-period optimization problem for optimal lot-sizing adapted from Hagen et al. (2001) shown below.

\[\begin{split}\begin{align} min \sum _{t \in T} c_{t} y_{t} + h_{t}^{+}I_{t}^{+} + h_{t}^{-}I_{t}^{-} \\ s.t. \;\;\;I_{t} &= I_{t-1} + X_{t} - d_{t} \\ I_{t} &= I_{t}^{+} - I_{t}^{-} \\ X_{t} &\leq Py_{t} \\ X_{t}, I_{t}^{+}, I_{t}^{-} &\geq 0 \\ y_{t} &\in \{0,1\} \end{align} \end{split}\]

Our goal is to find the optimal production \(X_{t}\) given known demands \(d_{t}\), fixed cost \(c_{t}\) associated with active production in a particular time period, an inventory holding cost \(h_{t}^{+}\) and a shortage cost \(h_{t}^{-}\) (cost of keeping a backlog) of orders. The variable \(y_{t}\) (binary) determines if we produce in time \(t\) or not, and \(I_{t}^{+}\) represents inventory that we are storing across time period \(t\), while \(h_{t}^{-}\) represents the magnitude of the backlog. Note that equation (4) is a constraint that only allows production in time period t if the indicator variable \(y_{t}=1\). Below we write a Pyomo model for this problem and solve it using glpk using the data provided below.

Parameter	Description	Value
\(c\)	fixed cost of production	4.6
\(I_{0}^{+}\)	initial value of positive inventory	5.0
\(I_{0}^{-}\)	initial value of backlogged orders	0.0
\(h^{+}\)	cost (per unit) of holding inventory	0.7
\(h^{-}\)	shortage cost (per unit)	1.2
\(P\)	maximum production amount (big-M value)	5
\(d\)	demand	[5,7,6.2,3.1,1.7]

import pyomo.environ as pyo

model = pyo.ConcreteModel()
model.T = pyo.RangeSet(5)    # time periods

i0 = 5.0           # initial inventory
c = 4.6            # setup cost
h_pos = 0.7        # inventory holding cost
h_neg = 1.2        # shortage cost
P = 5.0            # maximum production amount

# demand during period t
d = {1: 5.0, 2:7.0, 3:6.2, 4:3.1, 5:1.7}

# define the variables
model.y = pyo.Var(model.T, domain=pyo.Binary)
model.x = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.i = pyo.Var(model.T)
model.i_pos = pyo.Var(model.T, domain=pyo.NonNegativeReals)
model.i_neg = pyo.Var(model.T, domain=pyo.NonNegativeReals)

# define the inventory relationships
def inventory_rule(m, t):
    if t == m.T.first():
        return m.i[t] == i0 + m.x[t] - d[t]
    return m.i[t] == m.i[t-1] + m.x[t] - d[t]
model.inventory = pyo.Constraint(model.T, rule=inventory_rule)

def pos_neg_rule(m, t):
    return m.i[t] == m.i_pos[t] - m.i_neg[t]
model.pos_neg = pyo.Constraint(model.T, rule=pos_neg_rule)

# create the big-M constraint for the production indicator variable
def prod_indicator_rule(m,t):
    return m.x[t] <= P*m.y[t]
model.prod_indicator = pyo.Constraint(model.T, rule=prod_indicator_rule)

# define the cost function
def obj_rule(m):
    return sum(c*m.y[t] + h_pos*m.i_pos[t] + h_neg*m.i_neg[t] for t in m.T)
model.obj = pyo.Objective(rule=obj_rule)

# solve the problem
solver = pyo.SolverFactory('glpk')
solver.solve(model)

# print the results
for t in model.T:
    print('Period: {0}, Prod. Amount: {1}'.format(t, pyo.value(model.x[t]))) 

Period: 1, Prod. Amount: 3.0
Period: 2, Prod. Amount: 5.0
Period: 3, Prod. Amount: 5.0
Period: 4, Prod. Amount: 5.0
Period: 5, Prod. Amount: 0.0