Widgetco¶

Widgetco is about to introduce a new product. One unit of this product is produced by assembling subassembly 1 and subassembly 2. Before production begins on either subassembly, raw materials must be purchased and workers must be trained. Before the subassemblies can be assembled into the final product, the finished subassembly 2 must be inspected. The GAMS code that sets up this data is given below.

In [1]:
import numpy as np
import pandas as pd
import math
from gamspy import (
    Container,Set,Alias,Parameter,Variable,Equation,Model,Problem,Sense,Options,
    Domain,Number,Sum,Product,Smax,Smin,Ord,Card,SpecialValues,
)

options = Options(variable_listing_limit=0, equation_listing_limit=20)

m = Container(options=options)

activity = Set(m,'activity',records=[('A',"Train Workers"),('B',"Purchase Raw Materials"),
    ('C',"Make Subassembly 1"),('D',"Make Subassembly 2"),
    ('E',"Inspect Subassembly 2"),('F',"Assemble Subassemblies")])
I = Alias(m,'I',activity)
J = Alias(m,'J',activity)
pred = Set(m,'pred',domain=[I,J],description="I preceeds J", records=[('A','C'), ('B','C'), 
    ('A','D'), ('B','D'), ('D','E'), ('C','F'), ('E','F')])
duration = Parameter(m,'duration',domain=I,records=np.array([6, 9, 8, 7, 10, 12]))

/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 (
In [2]:
display(activity.records,pred.toList(),duration.records)
uni element_text
0 A Train Workers
1 B Purchase Raw Materials
2 C Make Subassembly 1
3 D Make Subassembly 2
4 E Inspect Subassembly 2
5 F Assemble Subassemblies 
[('A', 'C'),
 ('B', 'C'),
 ('A', 'D'),
 ('B', 'D'),
 ('D', 'E'),
 ('C', 'F'),
 ('E', 'F')]
I value
0 A 6.0
1 B 9.0
2 C 8.0
3 D 7.0
4 E 10.0
5 F 12.0
In [3]:
t = Variable(m,'t','positive',domain=I,description='time activity starts')
projDur = Variable(m,'projDur','free')

incidence = Equation(m,'incidence',domain=[I,J])
incidence[pred[I,J]]= t[J] >= t[I] + duration[I]

endTime = Equation(m,'endTime',domain=I)
endTime[I] = projDur >= t[I] + duration[I]

cpm = Model(m,"cpm",
    equations=m.getEquations(),
    problem=Problem.LP,
    sense=Sense.MIN,
    objective=projDur,
)

cpm.solve()
Out[3]:
Solver Status Model Status Objective Num of Equations Num of Variables Model Type Solver Solver Time
0 Normal OptimalGlobal 38 13 7 LP CPLEX 0.005
In [4]:
critical = Set(m,'critical',domain=I,description="critical activities")
critical[I] = Number(1).where[
    (Smax(J.where[pred[J,I]],incidence.m[J,I]) >= 1) |
    (Smax(J.where[pred[I,J]],incidence.m[I,J]) >= 1)]

display(critical.toList(),incidence.pivot(fill_value=0,value='marginal'))

['B', 'D', 'E', 'F']
C D E F
A 0.0 0.0 0.0 0.0
B -0.0 1.0 0.0 0.0
C 0.0 0.0 0.0 0.0
D 0.0 0.0 1.0 0.0
E 0.0 0.0 0.0 1.0

Now let's set up the dual

In [5]:
lamda = Variable(m,'lamda','positive',domain=[I,J],description='Dual var for incidence (critical arc)')
pi = Variable(m,'pi','positive',domain=I)

teq = Equation(m,'teq',domain=I,description='Equations for primal t vars')
teq[I]= ( -pi[I] - Sum(pred[I,J], lamda[I,J]) 
    + Sum(pred[J,I], lamda[J,I]) <= 0 )

projDur_dual_eq = Equation(m,'projDur_dual_eq',description='Dual equation for projDur')
projDur_dual_eq[:]= Sum(I, pi[I]) == 1 
                             
dual_cpm = Model(m,"dual_cpm",
    equations=[teq, projDur_dual_eq],
    problem=Problem.LP,
    sense=Sense.MAX,
    objective= (Sum(pred[I,J], lamda[I,J] * duration[I])
        + Sum(I,duration[I]*pi[I])),
)

dual_cpm.solve()
Out[5]:
Solver Status Model Status Objective Num of Equations Num of Variables Model Type Solver Solver Time
0 Normal OptimalGlobal 38 8 14 LP CPLEX 0.005
In [6]:
critical[I] = Number(1).where[
    (Smax(J.where[pred[J,I]],lamda.l[J,I]) >= 1) |
    (Smax(J.where[pred[I,J]],lamda.l[I,J]) >= 1)]

display(lamda.pivot(fill_value=0),critical.toList())
C D E F
A 0.0 0.0 0.0 0.0
B 0.0 1.0 0.0 0.0
C 0.0 0.0 0.0 0.0
D 0.0 0.0 1.0 0.0
E 0.0 0.0 0.0 1.0 
['B', 'D', 'E', 'F']

Now try to find early and late event times

In [7]:
projDur.fx[:] = projDur.l

eventtimes = m.addModel("eventtimes",
    equations=[incidence,endTime],
    problem=Problem.LP,
    sense=Sense.MAX,
    objective=Sum(I,t[I])
)

eventtimes.solve()
results = pd.DataFrame(
    t.toDense(),
    columns=['leTime'],
    index=I.records.uni)

eventtimes = Model(m,"eventtimes",
    equations=[incidence,endTime],
    problem=Problem.LP,
    sense=Sense.MIN,
    objective=Sum(I,t[I])
)
eventtimes.solve()
results['eeTime'] = t.toDense()

results['critical'] = results['eeTime'] >= results['leTime']
display(results)
leTime eeTime critical
uni
A 3.0 0.0 False
B 0.0 0.0 True
C 18.0 9.0 False
D 9.0 9.0 True
E 16.0 16.0 True
F 26.0 26.0 True