Sailco Corporation needs to determine how many sailboats to produce during each of the next four quarters, to meet the seasonal demands (shown below) in time at a minimal cost.

First quarter: d1 = 40 sailboats; Second quarter: d2 = 60 sailboats; Third quarter: d3 = 75 sailboats; Fourth quarter: d4 = 25 sailboats.

At the beginning of the first quarter, Sailco has an inventory of 10 sailboats. During each quarter, Sailco can produce up to 40 sailboats with regular-time labor at a cost of 400 dollars per boat. By having employees work overtime during a quarter, Sailco can produce additional sailboats with overtime labor at a cost of 450 dollars per boat.

We assume that sailboats manufactured during a quarter can be used to meet the demand of this quarter. At the end of each quarter, after production has occurred and the current quarter's demand has been satisfied, a holding cost of 20 dollarsper boat is incurred.

In [1]:
import sys
import pandas as pd
import numpy as np

from gamspy import (
    Container,Set,Alias,Parameter,Variable,Equation,Model,Problem,Sense,Options,
    Domain,Number,Sum,Product,Smax,Smin,Ord,Card,SpecialValues,
    ModelStatus,SolveStatus,
)
from gamspy.exceptions import GamspyException

options = Options(variable_listing_limit=100, equation_listing_limit=0,seed=101)
m = Container(options=options)

/Users/mas/anaconda3/lib/python3.11/site-packages/pandas/core/arrays/masked.py:60: UserWarning: Pandas requires version '1.3.6' or newer of 'bottleneck' (version '1.3.5' currently installed).
  from pandas.core import (

Model as inventory flow:¶

  • One node per month, with given demand at each month.

  • Inflow to each month node from "regular labor", with capacity 40 and cost 400

  • Inflow to each month node from a "overtime labor", with unlimited capacity and cost 450

  • month-to-month arcs representing inventory, with cost 20

  • balance constraints applied at month nodes.

In [3]:
months = Set(m,'months') 
demand = Parameter(m,'demand',months,domain_forwarding=True,
    records=[('month1', 40), ('month2', 60), ('month3', 75), ('month4', 25)])
supply = Parameter(m,'supply',months,description='External inventory supply',
    records=[('month1', 10)])

maxRegularLabor = Parameter(m, records=40, description='maximum boats produced by regular labor per month')
costRegular = Parameter(m, records=400, description='cost to build a boat with regular labor')
costOvertime = Parameter(m, records=450, description='cost to build a boat with overtime labor')
costInv = Parameter(m, records=20, description='cost to store a boat for one month')

# VARIABLES #
Regular = Variable(m,"Regular","positive",domain=months,description="boats produced by regular labor")
Inventory = Variable(m,"Inventory","positive",domain=months,description="inventory passed to next month")
Overtime = Variable(m,"Overtime","positive",domain=months,description="boats produced by overtime labor")

# EQUATIONS #
balance = Equation(m,'balance',domain=months,description='balance constraints at month nodes')
# demand = inventory  in and out, plus labor s in
balance[months]= ( Inventory[months] == Inventory[months.lag(1)] + 
    supply[months] - demand[months] + Regular[months] + Overtime[months] )

sailco = Model(m,
    name="sailco",
    equations=m.getEquations(),
    problem=Problem.LP,
    sense=Sense.MIN,
    # cost = inventory costs for first 3 months, plus regular and overtime labor costs
    objective=Sum(months, costInv*Inventory[months].where[~months.last] +
        costRegular*Regular[months]+ costOvertime*Overtime[months]),
)

# capacity constraint on regular labor arcs
Regular.up[months] = maxRegularLabor

sailco.solve(solver='cplex',solver_options={'lpmethod': 3, 'netfind': 2, 'preind': 0, 'names': 'no'},output=None)

print("Objective Function Value:  ", round(sailco.objective_value, 4), "\n")
display("Inventory:", Inventory.toList())
display("Regular:", Regular.toList())
display("Overtime:", Overtime.toList())

Objective Function Value:   78450.0 


'Inventory:'
[('month1', 10.0), ('month2', 0.0), ('month3', 0.0), ('month4', 0.0)]
'Regular:'
[('month1', 40.0), ('month2', 40.0), ('month3', 40.0), ('month4', 25.0)]
'Overtime:'
[('month1', 0.0), ('month2', 10.0), ('month3', 35.0), ('month4', 0.0)]

Model as a network:¶

  • One node per month, with given demand at each month.

  • Inflow to each month node from a "regular labor" node, with capacity 40 and cost 400

  • Inflow to each month node from a "overtime labor" node, with unlimited capacity and cost 450

  • Additional node for sailboats needed (to split into ones produced regular and overtime)

  • Month-to-month arcs representing inventory, (cost 20)

  • Balance constraints applied at all nodes.

In [4]:
# Now build the network (independently)
node = Set(m,'node',records=['Sailboat', 'Regular', 'Overtime'] + [f"m{i}" for i in range(5)])
BuildType = Set(m,'BuildType',domain=node,records=['Regular', 'Overtime'])
T = Set(m,'T',domain=node,records=[f"m{i}" for i in range(5)])

arc = Set(m,'arc',domain=[node,node]) 
arc['Sailboat', BuildType] = True #Connects Sailboat node to the buildtype nodes (regular and overtime)
arc[BuildType, T].where[~T.first] = True
arc[T,T.lead(1)] = True

display(arc.pivot(fill_value=''))

i = Alias(m,'i',node)
j = Alias(m,'j',node)
k = Alias(m,'k',node)

b = Parameter(m,'b',domain=node)
b[T].where[T.first] = supply['month1']
# transform data from one set into an other (months into T)
b[T].where[~T.first] = -Sum(months.where[T.ord-1==months.ord], demand[months])
# Need to figure out the b for the master supply node
b['Sailboat'] = -Sum(T, b[T]) #Generate the 190 boats required because 10 are already in supply
display(b.toList())

c = Parameter(m,'c',domain=[node,node])
c['Regular',T].where[~T.first] = costRegular
c['Overtime',T].where[~T.first] = costOvertime
c[T,T.lead(1)].where[~T.first] = costInv

u = Parameter(m,'u',domain=[node,node])
u[arc] = SpecialValues.POSINF 
u['Regular',T].where[~T.first] = maxRegularLabor

# Now it is just a regular MCNF (using A,b,c,u)

# VARIABLES #
x = Variable(m,"x","positive",domain=[node,node],description="flow")

# EQUATIONS #
# ensure flow out - flow in = b, use dynamic set 
f_balance = Equation(m,'f_balance',domain=node)
f_balance[i]= Sum(arc[i,k], x[i,k]) - Sum(arc[j,i], x[j,i]) == b[i]

mcf = Model(m,
    name="mcf",
    equations=[f_balance],
    problem=Problem.LP,
    sense=Sense.MIN,
    objective=Sum(arc[i,j], c[i,j]*x[i,j]),
)

x.up[arc] = u[arc] 

mcf.solve(solver='cplex',solver_options={'lpmethod': 3, 'netfind': 2, 'preind': 0, 'names': 'no'},output=None)

print("Objective Function Value:  ", round(mcf.objective_value, 4), "\n\nf:")
display(x.pivot())
Regular Overtime m1 m2 m3 m4
Sailboat True True
Regular True True True True
Overtime True True True True
m0 True
m1 True
m2 True
m3 True 
[('Sailboat', 190.0),
 ('m0', 10.0),
 ('m1', -40.0),
 ('m2', -60.0),
 ('m3', -75.0),
 ('m4', -25.0)]
Objective Function Value:   78450.0 

f:
Regular Overtime m1 m2 m3 m4
Sailboat 145.0 45.0 0.0 0.0 0.0 0.0
Regular 0.0 0.0 40.0 40.0 40.0 25.0
Overtime 0.0 0.0 0.0 10.0 35.0 0.0
m0 0.0 0.0 10.0 0.0 0.0 0.0
m1 0.0 0.0 0.0 10.0 0.0 0.0
m2 0.0 0.0 0.0 0.0 0.0 0.0
m3 0.0 0.0 0.0 0.0 0.0 0.0

Update initial stock and rerun¶

In [5]:
T1 = Set(m,'T1',domain=node,is_singleton=True,records=['m0']) #Only holds one element, throws out the last one every time there is a new one
b[T1] = 30 #saying there are 30 boats to start now
b[T].where[~T.first] = -Sum(months.where[T.ord-1==months.ord], demand[months])
b['Sailboat'] = -Sum(T, b[T])

mcf.solve(solver='cplex',solver_options={'lpmethod': 3, 'netfind': 2, 'preind': 0, 'names': 'no'},output=None)

print("Objective Function Value:  ", round(mcf.objective_value, 4), "\n\nf:")
display(x.pivot())

Objective Function Value:   70050.0 

f:
Regular Overtime m1 m2 m3 m4
Sailboat 145.0 25.0 0.0 0.0 0.0 0.0
Regular 0.0 0.0 40.0 40.0 40.0 25.0
Overtime 0.0 0.0 0.0 0.0 25.0 0.0
m0 0.0 0.0 30.0 0.0 0.0 0.0
m1 0.0 0.0 0.0 30.0 0.0 0.0
m2 0.0 0.0 0.0 0.0 10.0 0.0
m3 0.0 0.0 0.0 0.0 0.0 0.0

Do for random initial stock¶

In [6]:
def mysample(size):
    return rng.normal(loc=10.0, scale=4.0, size=size)

rng = np.random.default_rng(seed=101)
sample = mysample(10)
result = pd.DataFrame(columns=['iter','b-T1','cost'])

for cnt in range(len(sample)):
    if sample[cnt] < 0:
        raise Exception("Must have positive initial inventory", b)
    b[T1] = sample[cnt]
    b[T].where[~T.first] = -Sum(months.where[T.ord-1==months.ord], demand[months])
    b['Sailboat'] = -Sum(T, b[T])

    mcf.solve(solver='cplex',solver_options={'lpmethod': 3, 'netfind': 2, 'preind': 0, 'names': 'no'},output=None)
    if mcf.solve_status != SolveStatus.NormalCompletion:
        raise GamspyException("bad data given", mcf.solve_status)
    if mcf.status == ModelStatus.OptimalGlobal:
        result.loc[cnt] = ['iter'+str(cnt+1),sample[cnt],mcf.objective_value]
    else:
        result.loc[cnt] =  ['iter'+str(cnt+1),-1,0]

print(result)

     iter       b-T1          cost
0   iter1   6.839390  79809.062300
1   iter2   1.861498  81949.555829
2   iter3  12.413207  77412.320995
3   iter4  12.977178  77169.813409
4   iter5   8.761253  78982.661296
5   iter6  11.469285  77818.207239
6   iter7  16.841577  75508.121819
7   iter8  14.243191  76625.427715
8   iter9  12.830556  77232.860884
9  iter10  12.750998  77267.071052
In [7]:
import plotly.express as px

fig = px.line(result,x='b-T1',y='cost',markers=True)
fig.show()