Product Supply chain network design with uncertain demand

Logistics is the activity involving movement of a valuable transaction (product or empty!) from one location to another location. Need for such movement might arise from different factors depending upon the context. Movement can be used for various reasons such as meeting demand in distribution, movement of tenants, water for irrigation and drinking purposes among numerous other activities. There is a link that drives logistics and can be generically stated as movement of physical goods from point of supply to point of demand. Logistical activities involve movement of product and require resources that are capable of performing loading, transportation and unloading between consecutive supply and demand points, not necessarily in the same order always (empty movements). However, combination of successive supply and demand points forms a large chain of network commonly known as supply chain network from perspective of suppliers.

Supply chain management is the management of overall network for meeting customer demand of system for maximization of surplus. Supply chain surplus is the difference of revenue and costs of performing logistical, storage and production operations at various stages of the network. The management of supply chain can happen either at local or global level. Local supply chain management primarily deals with improving surplus of individual or subset of members while global supply chain management deals with improving overall network surplus. We put forth that logistical movement might comprise of loaded or empty movement in the network. The empty movement occurs when transport fleet returns back to supply location with vacant pallets. However, generally the fleets are managed by independent transportation providers and such empty movements are not part of network. The first and foremost factor is defining the system boundary for supply chain network. Such network would comprise of nodes and arcs in mathematical sense for analysis and decision making.

The article attempts to analyze product based supply chain comprising of bill of material (BOM) hierarchy in supply chain design for improving the supply chain surplus. Supply chain network comprises of regional cluster connected via ground transportation. Here, we have excluded other modes of transportation such as airfreight and sea freight to have a reasonable system boundary for decision making. To incorporate real world problems associated with inconsistent but predictive customer requirement the demand is stochastic.     

The problem deals with maximization of supply chain revenue and service level with novel fourfold dimensions focusing on BOM, uncertain demand and the movement of stock transfer order (STO) within the same echelon. Real life complexity is incorporated in the model with uncertain demand.

Python Code:
(Link: https://gist.github.com/adarsh70/6ce9bf9c21e3c9c7514ad8db5bb486ef)

#Import pulp and pandas packages

from pulp import*

import pandas as pd

#Raw Input Data

# Input Data

scenarios = [1,2,3]

probs = [0.1,0.6,0.3]

nodes = ['S1','S2','M1','M2','D1','D2','R1','R2']

S = [['S1','RM1'],['S2','RM2']]

M = [['M1','P1'],['M1','P2'],['M2','P2']]

D = [['D1','P1'],['D1','P2'],['D2','P2']]

R = [['R1','P1'],['R1','P2'],['R2','P2']]

T = [1,2,3,4,5]

RM = ['RM1','RM2']

P = ['P1','P2']

BOM = [['P1','RM1',8],['P1','RM2',5],['P2','RM2',4]]

SL = [['R1', 'P1', 1, 0.95],

 ['R1', 'P1', 2, 0.95],

 ['R1', 'P1', 3, 0.95],

 ['R1', 'P1', 4, 0.95],

 ['R1', 'P1', 5, 0.95],

 ['R1', 'P2', 1, 0.90],

 ['R1', 'P2', 2, 0.90],

 ['R1', 'P2', 3, 0.90],

 ['R1', 'P2', 4, 0.90],

 ['R1', 'P2', 5, 0.90],

 ['R2', 'P2', 1, 0.85],

 ['R2', 'P2', 2, 0.85],

 ['R2', 'P2', 3, 0.85],

 ['R2', 'P2', 4, 0.85],

 ['R2', 'P2', 5, 0.85]]

# Selling Price

price = [['R1', 'P1', 1, 21000],

 ['R1', 'P1', 2, 21000],

 ['R1', 'P1', 3, 11000],

 ['R1', 'P1', 4, 22000],

 ['R1', 'P1', 5, 26000],

 ['R1', 'P2', 1, 19000],

 ['R1', 'P2', 2, 10000],

 ['R1', 'P2', 3, 11000],

 ['R1', 'P2', 4, 24000],

 ['R1', 'P2', 5, 13000],

 ['R2', 'P2', 1, 12000],

 ['R2', 'P2', 2, 28000],

 ['R2', 'P2', 3, 24000],

 ['R2', 'P2', 4, 14000],

 ['R2', 'P2', 5, 19000]]

# Transportation Cost

trans = [['S1', 'S2', 550],

 ['S1', 'M1', 610],

 ['S1', 'M2', 280],

 ['S1', 'D1', 520],

 ['S1', 'D2', 640],

 ['S1', 'R1', 80],

 ['S1', 'R2', 360],

 ['S2', 'S1', 210],

 ['S2', 'M1', 70],

 ['S2', 'M2', 540],

 ['S2', 'D1', 760],

 ['S2', 'D2', 100],

 ['S2', 'R1', 500],

 ['S2', 'R2', 530],

 ['M1', 'S1', 610],

 ['M1', 'S2', 660],

 ['M1', 'M2', 60],

 ['M1', 'D1', 360],

 ['M1', 'D2', 520],

 ['M1', 'R1', 510],

 ['M1', 'R2', 430],

 ['M2', 'S1', 580],

 ['M2', 'S2', 670],

 ['M2', 'M1', 380],

 ['M2', 'D1', 200],

 ['M2', 'D2', 740],

 ['M2', 'R1', 520],

 ['M2', 'R2', 60],

 ['D1', 'S1', 530],

 ['D1', 'S2', 120],

 ['D1', 'M1', 370],

 ['D1', 'M2', 490],

 ['D1', 'D2', 290],

 ['D1', 'R1', 430],

 ['D1', 'R2', 660],

 ['D2', 'S1', 200],

 ['D2', 'S2', 50],

 ['D2', 'M1', 130],

 ['D2', 'M2', 620],

 ['D2', 'D1', 260],

 ['D2', 'R1', 350],

 ['D2', 'R2', 140],

 ['R1', 'S1', 460],

 ['R1', 'S2', 590],

 ['R1', 'M1', 670],

 ['R1', 'M2', 770],

 ['R1', 'D1', 370],

 ['R1', 'D2', 280],

 ['R1', 'R2', 210],

 ['R2', 'S1', 480],

 ['R2', 'S2', 240],

 ['R2', 'M1', 400],

 ['R2', 'M2', 480],

 ['R2', 'D1', 380],

 ['R2', 'D2', 510],

 ['R2', 'R1', 730]]

# Raw material inventory holding cost at manufacturer

hcr = [['M1', 'RM1', 1, 20],

 ['M1', 'RM1', 2, 20],

 ['M1', 'RM1', 3, 20],

 ['M1', 'RM1', 4, 20],

 ['M1', 'RM1', 5, 20],

 ['M1', 'RM2', 1, 10],

 ['M1', 'RM2', 2, 10],

 ['M1', 'RM2', 3, 10],

 ['M1', 'RM2', 4, 10],

 ['M1', 'RM2', 5, 10],

 ['M2', 'RM2', 1, 10],

 ['M2', 'RM2', 2, 10],

 ['M2', 'RM2', 3, 10],

 ['M2', 'RM2', 4, 10],

 ['M2', 'RM2', 5, 10]]

# Product inventory holding cost at manufacturer

hcp = [['M1', 'P1', 1, 70],

 ['M1', 'P1', 2, 70],

 ['M1', 'P1', 3, 70],

 ['M1', 'P1', 4, 70],

 ['M1', 'P1', 5, 70],

 ['M1', 'P2', 1, 50],

 ['M1', 'P2', 2, 50],

 ['M1', 'P2', 3, 50],

 ['M1', 'P2', 4, 50],

 ['M1', 'P2', 5, 50],

 ['M2', 'P2', 1, 50],

 ['M2', 'P2', 2, 50],

 ['M2', 'P2', 3, 50],

 ['M2', 'P2', 4, 50],

 ['M2', 'P2', 5, 50]]   

# Product inventory holding cost at distributor

hcd = [['D1', 'P1', 1, 80],

 ['D1', 'P1', 2, 80],

 ['D1', 'P1', 3, 80],

 ['D1', 'P1', 4, 80],

 ['D1', 'P1', 5, 80],

 ['D1', 'P2', 1, 60],

 ['D1', 'P2', 2, 60],

 ['D1', 'P2', 3, 60],

 ['D1', 'P2', 4, 60],

 ['D1', 'P2', 5, 60],

 ['D2', 'P2', 1, 60],

 ['D2', 'P2', 2, 60],

 ['D2', 'P2', 3, 60],

 ['D2', 'P2', 4, 60],

 ['D2', 'P2', 5, 60]]

# Supplier Capacity

scap = [['S1', 'RM1', 20000], ['S2', 'RM2', 20000]]

# Manufacturer production capacity

    

mprcap = [['M1', 'P1', 3000], ['M1', 'P2', 7000], ['M2', 'P2', 4000]]

# Manufacturer product storage capacity

mscp = [['M1', 'P1', 1800], ['M1', 'P2', 1800], ['M2', 'P2', 1800]]

# Manufacturer raw material storage capacity

mscr = [['M1', 'RM1', 800], ['M1', 'RM2', 750], ['M2', 'RM2', 900]]

# Distributor product storage capacity

mscd = [['D1', 'P1', 2500], ['D1', 'P2', 2000], ['D2', 'P2', 2000]]

# Demand 

dem = [['R1','P1',400, 0.1,600, 0.6,800, 0.3],['R1','P2',300, 0.1,400, 0.6,500, 0.3],['R2','P2',200, 0.1,300, 0.6,400, 0.3]]

# Production Cost

prodcost = [['M1','P1',600],['M1','P2',480],['M2','P2',680]]

# Initial raw material inventory at manufacturer 

iimr = [['M1','RM1',40],['M1','RM2',100],['M2','RM2',80]]

# Initial product inventory at manufacturer

iimp = [['M1','P1',25],['M1','P2',50],['M2','P2',60]]

# Initial product inventory at distributor

iid = [['D1','P1',180],['D1','P2',350],['D2','P2',80]]

#Pre-processing of data

# Enumeration of Supplier to Manufacturer Flow List

smarc1 = []

for i,j in enumerate(S):

    for k,l in enumerate(M):

        for m,n in enumerate(BOM):

            if n[0] == l[1] and n[1] == j[1]:

                for t in T:

                    for s in scenarios:

                        smarc1.append([j[0],l[0],j[1],t,s])

                    

smarc = []

[smarc.append(item) for item in smarc1 if item not in smarc]

# Enumeration of Manufacturer to Distributor Flow List

mdarc1 = []

for i,j in enumerate(M):

    for k,l in enumerate(D):

        if j[1] == l[1]:

            for t in T:

                for s in scenarios:

                    mdarc1.append([j[0],l[0],j[1],t,s])

                

mdarc = []

[mdarc.append(item) for item in mdarc1 if item not in mdarc]

# Enumeration of Distributor to Retailer Flow List

drarc1 = []

for i,j in enumerate(D):

    for k,l in enumerate(R):

        if j[1] == l[1]:

            for t in T:

                for s in scenarios:

                    drarc1.append([j[0],l[0],j[1],t,s])

                

drarc = []

[drarc.append(item) for item in drarc1 if item not in drarc]

# Enumeration of Manufacturer Raw Material Stock Transfer Flow List

stmrarc1 = []

for i,j in enumerate(M):

    for k,l in enumerate(M):

        for m,n in enumerate(BOM):

            if j[0] != l[0] and n[0] == l[1] and n[0] == j[1]:

                for t in T:

                    for s in scenarios:

                        stmrarc1.append([j[0],l[0],n[1],t,s])

                    

stmrarc = []

[stmrarc.append(item) for item in stmrarc1 if item not in stmrarc]

# Enumeration of Manufacturer Product Stock Transfer Flow List

stmparc1 = []

for i,j in enumerate(M):

    for k,l in enumerate(M):

        if j[0] != l[0] and j[1] == l[1]:

            for t in T:

                for s in scenarios:

                    stmparc1.append([j[0],l[0],j[1],t,s])

                    

stmparc = []

[stmparc.append(item) for item in stmparc1 if item not in stmparc]

# Enumeration of Distributor Product Stock Transfer Flow List

stdarc1 = []

for i,j in enumerate(D):

    for k,l in enumerate(D):

        if j[0] != l[0] and j[1] == l[1]:

            for t in T:

                for s in scenarios:

                    stdarc1.append([j[0],l[0],j[1],t,s])

                    

stdarc = []

[stdarc.append(item) for item in stdarc1 if item not in stdarc]

# Enumeration of Raw Material Inventory at Manufacturer

imr1 = []

for i,j in enumerate(M):

    for k,l in enumerate(BOM):

        if j[1] == l[0]:

            for t in T:

                for s in scenarios:

                    imr1.append([j[0],l[1],t,s])

                

imr = []

[imr.append(item) for item in imr1 if item not in imr]

# Enumeration of Product Inventory at Manufacturer

imp1 = []

for i,j in enumerate(M):

    for t in T:

        for s in scenarios:

            imp1.append([j[0],j[1],t,s])

                

imp = []

[imp.append(item) for item in imp1 if item not in imp]

# Enumeration of Product Inventory at Distributor

ind1 = []

for i,j in enumerate(D):

    for t in T:

        for s in scenarios:

            ind1.append([j[0],j[1],t,s])

                

ind = []

[ind.append(item) for item in ind1 if item not in ind]

# Enumeration of supplier, raw material and time list

srt1 = []

for i,j in enumerate(S):

    for t in T:

        for s in scenarios:

            srt1.append([j[0],j[1],t,s])

        

srt = []

[srt.append(item) for item in srt1 if item not in srt]

# Enumeration of manufacturer, raw material, product and time list

mrpt1 = []

for i,j in enumerate(M):

    for k,l in enumerate(BOM):

        if j[1] == l[0]:

            for t in T:

                for s in scenarios:

                    mrpt1.append([j[0],l[1],j[1],t,s])

        

mrpt = []

[mrpt.append(item) for item in mrpt1 if item not in mrpt]

# Enumeration of retailer, product and time list

rpt1 = []

for i,j in enumerate(R):

        for t in T:

            for s in scenarios:

                rpt1.append([j[0],j[1],t,s])

        

rpt = []

[rpt.append(item) for item in rpt1 if item not in rpt]

# Converting the data list to dataframe

SL = pd.DataFrame(SL,columns=['Retailer','Product','Time','SL'])

SL = SL.set_index(['Retailer','Product','Time'])

price = pd.DataFrame(price,columns=['Retailer','Product','Time','Price'])

price = price.set_index(['Retailer','Product','Time'])

trans = pd.DataFrame(trans,columns=['From','To','TransCost'])

trans = trans.set_index(['From','To'])

hcr = pd.DataFrame(hcr,columns=['Manufacturer','RawMaterial','Time','InvCost'])

hcr = hcr.set_index(['Manufacturer','RawMaterial','Time'])

hcp = pd.DataFrame(hcp,columns=['Manufacturer','Product','Time','InvCost'])

hcp = hcp.set_index(['Manufacturer','Product','Time'])

hcd = pd.DataFrame(hcd,columns=['Distributor','Product','Time','InvCost'])

hcd = hcd.set_index(['Distributor','Product','Time'])

qty = pd.DataFrame(BOM,columns=['Product','RawMaterial','Quantity'])

qty = qty.set_index(['Product','RawMaterial'])

scap = pd.DataFrame(scap,columns=['Supplier','RawMaterial','Capacity'])

scap = scap.set_index(['Supplier','RawMaterial'])

mprcap = pd.DataFrame(mprcap,columns=['Manufacturer','Product','ProdCap'])

mprcap = mprcap.set_index(['Manufacturer','Product'])

mscp = pd.DataFrame(mscp,columns=['Manufacturer','Product','ProdStorCap'])

mscp = mscp.set_index(['Manufacturer','Product'])

mscr = pd.DataFrame(mscr,columns=['Manufacturer','RawMaterial','ProdRawCap'])

mscr = mscr.set_index(['Manufacturer','RawMaterial'])

mscd = pd.DataFrame(mscd,columns=['Distributor','Product','ProdStorCap'])

mscd = mscd.set_index(['Distributor','Product'])

dem = pd.DataFrame(dem,columns=['Retailer','Product','LowDemand','LowProb','MedDemand','MedProb','HighDemand','HighProb'])

dem = dem.set_index(['Retailer','Product'])

prodcost = pd.DataFrame(prodcost,columns=['Manufacturer','Product','ProdCost'])

prodcost = prodcost.set_index(['Manufacturer','Product'])

iimr = pd.DataFrame(iimr,columns=['Manufacturer','RawMaterial','Inventory'])

iimr = iimr.set_index(['Manufacturer','RawMaterial'])

iimp = pd.DataFrame(iimp,columns=['Manufacturer','Product','Inventory'])

iimp = iimp.set_index(['Manufacturer','Product'])

iid = pd.DataFrame(iid,columns=['Distributor','Product','Inventory'])

iid = iid.set_index(['Distributor','Product'])

#Product Supply Chain Network Design

x = LpVariable.dicts("Supp_to_Mfg_Flow_variable", ((i,j,r,t,s) for i,j,r,t,s in smarc),lowBound=0,cat='Integer')

y = LpVariable.dicts("Mfg_to_Distr_Flow_variable", ((j,k,p,t,s) for j,k,p,t,s in mdarc),lowBound=0,cat='Integer')

z = LpVariable.dicts("Distr_to_Retlr_Flow_variable", ((k,l,p,t,s) for k,l,p,t,s in drarc),lowBound=0,cat='Integer')

stmr = LpVariable.dicts("Stock_Transfer_Mfg_Rawm_Flow_variable", ((i,j,r,t,s) for i,j,r,t,s in stmrarc),lowBound=0,cat='Integer')

stmp = LpVariable.dicts("Stock_Transfer_Mfg_Prod_Flow_variable", ((i,j,p,t,s) for i,j,p,t,s in stmparc),lowBound=0,cat='Integer')

std = LpVariable.dicts("Stock_Transfer_Distr_Flow_variable", ((j,k,p,t,s) for j,k,p,t,s in stdarc),lowBound=0,cat='Integer')

imfgr = LpVariable.dicts("Inv_Mfg_Rawm_variable", ((j,r,t,s) for j,r,t,s in imr),lowBound=0,cat='Integer')

imfgp = LpVariable.dicts("Inv_Mfg_Prod_variable", ((j,p,t,s) for j,p,t,s in imp),lowBound=0,cat='Integer')

indr = LpVariable.dicts("Inv_Distr_Prod_variable", ((k,p,t,s) for k,p,t,s in ind),lowBound=0,cat='Integer')

rq = LpVariable.dicts("RawMtrl_Used_Mfg_variable", ((j,r,p,t,s) for j,r,p,t,s in mrpt),lowBound=0,cat='Integer')

pq = LpVariable.dicts("Product_Produced_Mfg_variable", ((j,p,t,s) for j,p,t,s in imp),lowBound=0,cat='Integer')

dsl = LpVariable.dicts("Demand_Slack_variable", ((l,p,t,s) for l,p,t,s in rpt),lowBound=0,cat='Integer')

model = LpProblem("Product Supply Chain Network Design",LpMaximize)

# Objective Function

model += lpSum([probs[int(s-1)] * price.loc[(l,p,t),'Price'] * z[k,l,p,t,s] for k,l,p,t,s in drarc]) - lpSum([probs[int(s-1)] * trans.loc[(i,j),'TransCost'] * x[i,j,r,t,s] for i,j,r,t,s in smarc]

                                                                                      + [probs[int(s-1)] * trans.loc[(j,k),'TransCost'] * y[j,k,p,t,s] for j,k,p,t,s in mdarc]

                                                                                      + [probs[int(s-1)] * trans.loc[(k,l),'TransCost'] * z[k,l,p,t,s] for k,l,p,t,s in drarc]

                                                                                      + [probs[int(s-1)] * trans.loc[(i,j),'TransCost'] * stmr[i,j,r,t,s] for i,j,r,t,s in stmrarc]

                                                                                      + [probs[int(s-1)] * trans.loc[(i,j),'TransCost'] * stmp[i,j,p,t,s] for i,j,p,t,s in stmparc]

                                                                                      + [probs[int(s-1)] * trans.loc[(j,k),'TransCost'] * std[j,k,p,t,s] for j,k,p,t,s in stdarc]

                                                                                      + [probs[int(s-1)] * hcr.loc[(j,r,t),'InvCost'] * imfgr[j,r,t,s] for j,r,t,s in imr]

                                                                                      + [probs[int(s-1)] * hcp.loc[(j,p,t),'InvCost'] * imfgp[j,p,t,s] for j,p,t,s in imp]

                                                                                      + [probs[int(s-1)] * hcd.loc[(k,p,t),'InvCost'] * indr[k,p,t,s] for k,p,t,s in ind]

                                                                                      + [probs[int(s-1)] * prodcost.loc[(j,p),'ProdCost'] * pq[j,p,t,s] for j,p,t,s in imp])

# supplier capacity constraint

for i,r,t,s in srt:

    model += lpSum([x[a,j,b,c,d] for a,j,b,c,d in x if a == i and b == r and c == t and d == s]) <= scap.loc[(i,r),'Capacity']

# BOM conversion constraint

for j,r,p,t,s in mrpt:

    model += rq[j,r,p,t,s] == qty.loc[(p,r),'Quantity'] * pq[j,p,t,s]

    

# manufacturer raw material balance constraint 

for j,r,t,s in imr:

    if t == 1:

        model += lpSum([stmr[i,a,b,c,d] for i,a,b,c,d in stmr if a == j and b == r and c == t and d == s] + [x[i,a,b,c,d] for i,a,b,c,d in x if a == j and b == r and c == t and d == s]) + iimr.loc[(j,r),'Inventory'] == imfgr[j,r,t,s] + lpSum([stmr[a,i,b,c,d] for a,i,b,c,d in stmr if a == j and b == r and c == t and d == s] + [rq[a,b,p,c,d] for a,b,p,c,d in rq if a == j and b == r and c == t and d == s]) 

for j,r,t,s in imr:

    if t > 1:

        model += lpSum([stmr[i,a,b,c,d] for i,a,b,c,d in stmr if a == j and b == r and c == t and d == s] + [x[i,a,b,c,d] for i,a,b,c,d in x if a == j and b == r and c == t and d == s]) + imfgr[j,r,int(t-1),s] == imfgr[j,r,t,s] + lpSum([stmr[a,i,b,c,d] for a,i,b,c,d in stmr if a == j and b == r and c == t and d == s] + [rq[a,b,p,c,d] for a,b,p,c,d in rq if a == j and b == r and c == t and d == s]) 

# manufacturer product balance constraint 

for j,p,t,s in imp:

    if t == 1:

        model += lpSum([stmp[i,a,b,c,d] for i,a,b,c,d in stmp if a == j and b == p and c == t and d == s]) + pq[j,p,t,s] + iimp.loc[(j,p),'Inventory'] == imfgp[j,p,t,s] + lpSum([stmp[a,i,b,c,d] for a,i,b,c,d in stmp if a == j and b == p and c == t and d == s] + [y[a,k,b,c,d] for a,k,b,c,d in y if a == j and b == p and c == t and d == s])

for j,p,t,s in imp:

    if t > 1:

        model += lpSum([stmp[i,a,b,c,d] for i,a,b,c,d in stmp if a == j and b == p and c == t and d == s]) + pq[j,p,t,s] + imfgp[j,p,int(t-1),s] == imfgp[j,p,t,s] + lpSum([stmp[a,i,b,c,d] for a,i,b,c,d in stmp if a == j and b == p and c == t and d == s] + [y[a,k,b,c,d] for a,k,b,c,d in y if a == j and b == p and c == t and d == s])

# manufacturer production capacity constraint

for j,p,t,s in imp:

    model += pq[j,p,t,s] <= mprcap.loc[(j,p),'ProdCap']

    

# manufacturer product storage capacity constraint

for j,p,t,s in imp:

    model += imfgp[j,p,t,s] <= mscp.loc[(j,p),'ProdStorCap']

    

# manufacturer raw material storage capacity constraint

for j,r,t,s in imr:

    model += imfgr[j,r,t,s] <= mscr.loc[(j,r),'ProdRawCap']

    

# distributor product balance constraint 

for k,p,t,s in ind:

    if t == 1:

        model += lpSum([std[i,a,b,c,d] for i,a,b,c,d in std if a == k and b == p and c == t and d == s]) + lpSum([y[j,a,b,c,d] for j,a,b,c,d in y if a == k and b == p and c == t and d == s]) + iid.loc[(k,p),'Inventory'] == indr[k,p,t,s] + lpSum([std[a,i,b,c,d] for a,i,b,c,d in std if a == k and b == p and c == t and d == s] + [z[a,l,b,c,d] for a,l,b,c,d in z if a == k and b == p and c == t and d == s])

for k,p,t,s in ind:

    if t > 1:

        model += lpSum([std[i,a,b,c,d] for i,a,b,c,d in std if a == k and b == p and c == t and d == s]) + lpSum([y[j,a,b,c,d] for j,a,b,c,d in y if a == k and b == p and c == t and d == s]) + indr[k,p,int(t-1),s] == indr[k,p,t,s] + lpSum([std[a,i,b,c,d] for a,i,b,c,d in std if a == k and b == p and c == t and d == s] + [z[a,l,b,c,d] for a,l,b,c,d in z if a == k and b == p and c == t and d == s])

# distributor product storage capacity constraint

for k,p,t,s in ind:

    model += indr[k,p,t,s] <= mscd.loc[(k,p),'ProdStorCap']

    

# demand balance constraint

for l,p,t,s in rpt:

    if s == 1:

        model += lpSum([z[k,a,b,c,d] for k,a,b,c,d in z if a == l and b == p and c == t and d == s]) + dsl[l,p,t,s] == dem.loc[(l,p),'LowDemand']  

for l,p,t,s in rpt:

    if s == 2:

        model += lpSum([z[k,a,b,c,d] for k,a,b,c,d in z if a == l and b == p and c == t and d == s]) + dsl[l,p,t,s] == dem.loc[(l,p),'MedDemand']  

for l,p,t,s in rpt:

    if s == 3:

        model += lpSum([z[k,a,b,c,d] for k,a,b,c,d in z if a == l and b == p and c == t and d == s]) + dsl[l,p,t,s] == dem.loc[(l,p),'HighDemand'] 

        

# demand service level constraint

for l,p,t,s in rpt:

    if s == 1:

        model += lpSum([z[k,a,b,c,d] for k,a,b,c,d in z if a == l and b == p and c == t and d == s]) >= dem.loc[(l,p),'LowDemand'] * SL.loc[(l,p,t),'SL'] 

        

for l,p,t,s in rpt:

    if s == 2:

        model += lpSum([z[k,a,b,c,d] for k,a,b,c,d in z if a == l and b == p and c == t and d == s]) >= dem.loc[(l,p),'MedDemand'] * SL.loc[(l,p,t),'SL']    

        

for l,p,t,s in rpt:

    if s == 3:

        model += lpSum([z[k,a,b,c,d] for k,a,b,c,d in z if a == l and b == p and c == t and d == s]) >= dem.loc[(l,p),'HighDemand'] * SL.loc[(l,p,t),'SL']    

#Solve the integer linear programming problem and get the solution

model.writeLP("ProductSCND.lp")

model.solve()

LpStatus[model.status]

# Print our objective function value (Total Revenue)

print (value(model.objective))

for i,j,r,t,s in x:

    if x[(i,j,r,t,s)].varValue > 0:

        v = (x[(i,j,r,t,s)].varValue)

        print ("Supplier %s to Manufacturer %s of raw material %s in period %s for scenario %s with quantity %s is moved!" % (i,j,r,t,s,v))

for j,k,p,t,s in y:

    if y[(j,k,p,t,s)].varValue > 0:

        v = (y[(j,k,p,t,s)].varValue)

        print ("Manufacturer %s to Distributor %s of end product %s in period %s for scenario %s with quantity %s is moved!" % (j,k,p,t,s,v))

for k,l,p,t,s in z:

    if z[(k,l,p,t,s)].varValue > 0:

        v = (z[(k,l,p,t,s)].varValue)

        print ("Distributor %s to Retailer %s of end product %s in period %s for scenario %s with quantity %s is moved!" % (k,l,p,t,s,v))

for j,r,p,t,s in rq:

    if rq[(j,r,p,t,s)].varValue > 0:

        v = (rq[(j,r,p,t,s)].varValue)

        print ("Manufacturer %s use Raw material %s to produce end product %s in period %s for scenario %s with quantity %s for production!" % (j,r,p,t,s,v))

for j,p,t,s in pq:

    if pq[(j,p,t,s)].varValue > 0:

        v = (pq[(j,p,t,s)].varValue)

        print ("Manufacturer %s produces end product %s in period %s for scenario %s with quantity %s" % (j,p,t,s,v))

for l,p,t,s in dsl:

    if dsl[(l,p,t,s)].varValue > 0:

        v = (dsl[(l,p,t,s)].varValue)

        print ("Demand of retailer %s for end product %s in period %s for scenario %s is not met %s quantity!" % (l,p,t,s,v))

for j,r,t,s in imfgr:

    if imfgr[(j,r,t,s)].varValue > 0:

        v = (imfgr[(j,r,t,s)].varValue)

        print ("Inventory at manufacturer %s for raw material %s in period %s for scenario %s is stored with %s quantity!" % (j,r,t,s,v))

for j,p,t,s in imfgp:

    if imfgp[(j,p,t,s)].varValue > 0:

        v = (imfgp[(j,p,t,s)].varValue)

        print ("Inventory at manufacturer %s for end product %s in period %s for scenario %s is stored with %s quantity!" % (j,p,t,s,v))

for k,p,t,s in indr:

    if indr[(k,p,t,s)].varValue > 0:

        v = (indr[(k,p,t,s)].varValue)

        print ("Inventory at distributor %s for end product %s in period %s for scenario %s is stored with %s quantity!" % (k,p,t,s,v))

for i,j,r,t,s in stmr:

    if stmr[(i,j,r,t,s)].varValue > 0:

        v = (stmr[(i,j,r,t,s)].varValue)

        print ("Stock Transfer from manufacturer %s to manufacturer %s of raw material %s in period %s for scenario %s with quantity %s is moved!" % (i,j,r,t,s,v))

for i,j,p,t,s in stmp:

    if stmp[(i,j,p,t,s)].varValue > 0:

        v = (stmp[(i,j,p,t,s)].varValue)

        print ("Stock Transfer from manufacturer %s to manufacturer %s of end product %s in period %s for scenario %s with quantity %s is moved!" % (i,j,p,t,s,v))

for j,k,p,t,s in std:

    if std[(j,k,p,t,s)].varValue > 0:

        v = (std[(j,k,p,t,s)].varValue)

        print ("Stock Transfer from distributor %s to distributor %s of end product %s in period %s for scenario %s with quantity %s is moved!" % (j,k,p,t,s,v))

Dijkastra Algorithm for Shortest Path Problem (SPP) in R (Code to be tested and validated)

R Code for SPP algorithm

ShortestPathProblemDijkastra <- function(SPData) {

ArcData <- SPData

colnames(ArcData) <- c("fromNode","toNode","arcDist")

RawNodeList <- c(ArcData[,1],ArcData[,2])

RawNodeList <- unique(RawNodeList)

RawNodeList <- data.frame(RawNodeList)

PreNodes <- unique(c(ArcData[,2]))

PreNodes <- data.frame(PreNodes)

PostNodes <- unique(c(ArcData[,1]))

PostNodes <- data.frame(PostNodes)

library(dplyr)

colnames(PreNodes) <- c("Node")

colnames(PostNodes) <- c("Node")

colnames(RawNodeList) <- c("Node")

NoPreNodes <- anti_join(RawNodeList,PreNodes,by=c("Node"))

NoPostNodes <- anti_join(RawNodeList, PostNodes,by=c("Node"))

Startnode <- NoPreNodes[1]

Endnode <- NoPostNodes[1]

NodeLabel <- data.frame(matrix(NA,ncol = 4))

colnames(NodeLabel) <- c("fromNode","Dist","PreNode","Type")

NodeLabel[1,1] <- NoPreNodes[1]

NodeLabel[1,2] <- 0

NodeLabel[1,4] <- ("Permanent")

i <- 2

dist <- 0

NodeLabelMod <- data.frame(matrix(NodeLabel[1,1],nrow = 1,ncol = 1))

colnames(NodeLabelMod) <- ("fromNode") 

while (nrow(NodeLabel) < nrow(RawNodeList)) {

ds <- merge(NodeLabelMod,ArcData)

aa <- data.frame(ds$fromNode[ds$arcDist == min(ds$arcDist)],

ds$toNode[ds$arcDist == min(ds$arcDist)], ds$arcDist[ds$arcDist == 

min(ds$arcDist)])

NodeLabel[i,1] <- aa[1,2]

NodeLabel[i,2] <- NodeLabel[i-1,2] + aa[1,3]

NodeLabel[i,3] <- NoPreNodes[1]

NodeLabel[i,4] <- ("Permanent")

NodeLabelMod[1] <- NodeLabel[i,1]

NoPreNodes[1] <- NodeLabel[i,1]

dist <- dist + aa[1,3] 

i <- i + 1

}

rm(list=c("ArcData","RawNodeList","PreNodes","PostNodes",

"NoPreNodes","NoPostNodes","Startnode","Endnode","i","dist",

"ds","aa"))

NodeLabel <- NodeLabel[,-4] 

colnames(NodeLabel) <- c("Node","Dist","PreNode")

return(NodeLabel)

}

R Code for finding shortest path to a node based on preceding result (res) and node

NodeNetwork <- function(res,node) { 

seq <- NULL

n <- which(res$Node == as.character(node))

for (i in 1:n) {

aa <- res[i,1]

seq <- paste(seq,aa,sep="---->")

}

seq <- substring(seq,6)

MinDist <- res[n,2]

return(paste("The sequence for flow is: ",seq, " with distance of: 

",MinDist,sep=""))

}

Illustration:

ReadData <- read.csv("SPP_Data.csv")

SPP_Data.csv:

 
 


  
fromnode
  
tonode
  
arcdistance
 

  
1
  
2
  
5
 

  
1
  
3
  
1
 

  
2
  
4
  
7
 

  
2
  
5
  
1
 

  
2
  
6
  
6
 

  
3
  
2
  
2
 

  
3
  
4
  
6
 

  
3
  
5
  
7
 

  
4
  
3
  
7
 

  
4
  
6
  
4
 

  
4
  
7
  
6
 

  
5
  
4
  
3
 

  
5
  
6
  
5
 

  
5
  
7
  
9
 

  
6
  
2
  
7
 

  
6
  
7
  
2

res <- ShortestPathProblemDijkastra(ReadData)

NodeNetwork(res,7)


Result:



[1] "The sequence for flow is: 1---->3---->2---->5---->4---->6---->7 with 
distance of: 13"

Manpower Planning using Shiny in R

Manpower  planning is an operational problem that requires determining number of manpower for performing standardised work. The post focuses on implementation of such system with use of functional programming in R. Shiny is a powerful library available in R that helps in intertwining code with user interface (ui) for building a GUI (graphical user interface). We will describe R scripts for server and ui as we go along with the implementation. 

Script for ui in R

First step is to load the shiny package.

library(shiny)

Next define the UI for random distribution application

shinyUI(fluidPage(

The title of application in GUI window is created as:

titlePanel("Workstation Manpower Planning"),

The sidebar layout is created with controls to select the random distribution type and number of observations to generate. Note that the use of the br() element is to introduce extra vertical spacing. In side bar panel, input files are loaded as csv. Product changeover data sets account for change over time in switching from one product to another product. ProductChangeover.csv file consists of workstation, from product, to product and associated changeover time. Demand file has associated demand for the three products. Work station available time depicts the time available at workstation in a shift, by removing the downtime. Work content provides statistics for operation time of different products at various workstations. The statistics are mean and standard deviation of set up and manual time, while fixed processing times on machines, if any. Data sets are shown at end of posts for your understanding.

  sidebarLayout(

    sidebarPanel(

      fileInput('WCfile', 'Select the Workcontent.csv file',

                accept=c('text/csv','text/comma-separated-values,text/plain','.csv')),

      tags$hr(),

      fileInput('COfile', 'Select the ProductChangeover.csv file',

                accept=c('text/csv','text/comma-separated-values,text/plain','.csv')),

      tags$hr(),

      fileInput('DEMfile', 'Select the Demand.csv file',

                accept=c('text/csv','text/comma-separated-values,text/plain','.csv')),

      tags$hr(),

      fileInput('WSATfile', 'Select the Workstationavailtimes.csv file',

                accept=c('text/csv','text/comma-separated-values,text/plain','.csv')),

      tags$hr(),

      helpText("All .csv files must be uploaded in required format under headers")

    ),

Next, we have the main panel that is declare in tab form with three categories: first one is output table on planned manpower allocation at the workstations, second tab is balancing efficiency plot and third tab is plotting of product mix losses.

    mainPanel(

      tabsetPanel(type="tab",

                  tabPanel("Summary - Workstation Planned Manpower",dataTableOutput("mytable")),

                  tabPanel("Plot- Workstation Manpower Efficiency", plotOutput("myplot")),

                  tabPanel("Plot - Workstation Product Mix Loss", plotOutput("myplot1"))

                  )

   

    )

  ))

)

Script for server in R:

First we load the shiny library:

library(shiny)

Next, we create a function under shinyServer that declares input from ui script.

shinyServer(function(input, output) {

  observe({
    WCfile = input$WCfile
    COfile = input$COfile
    DEMfile = input$DEMfile
    WSATfile = input$WSATfile
    if (is.null(WCfile) || is.null(COfile) || is.null(DEMfile) || is.null(WSATfile)) {
      return(NULL)
    }

The input data is stored in different data frames from csv files for work content, product changeover times, demand and workstation available times.
 
    workcontent.df <- data.frame(read.csv(WCfile$datapath))
    productco.df <- data.frame(read.csv(COfile$datapath))
    dem.df <- data.frame(read.csv(DEMfile$datapath))
    wsavailtim.df <- data.frame(read.csv(WSATfile$datapath))
 
Step 1 calculations require determining maximum manual work content for each workstation among all products which is MAX(Manual time + Set-up time) + MAX(Change-over time)

 
    workcontent.df$totalmanualtime <- workcontent.df$Mean_Set.up + workcontent.df$SD_Set.up + workcontent.df$Mean_Manual + workcontent.df$SD_Manual
    co.df <- data.frame(Workstation=productco.df$Workstation, time=productco.df$Changeover_time)
    maxcotime.df <- aggregate(.~Workstation,co.df,FUN = max)
    wc.df <- data.frame(Workstation=workcontent.df$Workstation,time=workcontent.df$totalmanualtime)
    maxwc.df <- aggregate(.~Workstation,wc.df,FUN = max)
    aggregwc.df <- rbind(maxwc.df,maxcotime.df)
    manualwc.df <- aggregate(.~Workstation,aggregwc.df,FUN = sum)
 

Step 2 Calculations is to determine maximum production units for each workstation i.e. SUM(Workstation production units for all products)
 
    wswc.df <-data.frame(Workstation=workcontent.df$Workstation,Product=workcontent.df$Product)
    demwspr.df <- merge(wswc.df,dem.df,by="Product")
    aggrdemws.df <- data.frame(Workstation=demwspr.df$Workstation,Demand=demwspr.df$Demand)
    aggrdemws.df <- aggregate(.~Workstation,aggrdemws.df,FUN = sum)
 
 
Step 3 calculations determine cycle time for each workstation i.e. Max production units / Available time. Step 4 determines number of workers for each workstation as ROUNDUP(Max work content/Cycle time) and Step 5 determine Manpower Efficiency for each workstation as (# Operators/(Max work content/Cycle time))
 
    cycletime.df <- merge(wsavailtim.df,aggrdemws.df,by="Workstation")
    cycletime.df$ct <- cycletime.df$Demand/cycletime.df$Availabletime
    NoofOper.df <- merge(manualwc.df,cycletime.df,by="Workstation")
    NoofOper.df$Operators <- ceiling(NoofOper.df$time/NoofOper.df$ct)
    NoofOper.df$Efficiency <- round(((NoofOper.df$time/NoofOper.df$ct)/NoofOper.df$Operators)*100,0)
    NoofOpersop.df <- data.frame(Workstation = NoofOper.df$Workstation, Number_Of_Manpower = NoofOper.df$Operators,Manpower_Efficiency = paste(NoofOper.df$Efficiency,"%"))
    NoofOpersef.df <- data.frame(Workstation = NoofOper.df$Workstation, Manpower_Efficiency = NoofOper.df$Efficiency)
 
Step 6 calculates product mix loss percentage for the workstations. First, calculate the manual time by adding mean and standard deviations of manual time and set up time. This will aid in design of workstation for maximum time. Next, multiply the total work content with demand of the products for the workstations. After that, get the total work content on the workstations, obtained by summing it up for all products. The design of workstation requires work content that is multiplication of Step 1 and Step 2 time and quantities. The product mix loss percentage is [1 - (Actual Work content / Design Work content)].
 
    wcwspr <- data.frame(Workstation = workcontent.df$Workstation,Product = workcontent.df$Product, time_des = workcontent.df$totalmanualtime)
    wcwspr <- merge(wcwspr,dem.df,by="Product")
    wcwspr$ActualWC <- wcwspr$time_des * wcwspr$Demand
    aggrwcws <- aggregate(.~Workstation,data = wcwspr,FUN = sum)
    aggrwcws <- merge(aggrwcws,manualwc.df,by="Workstation")
    aggrwcws$DesgnWC <- aggrwcws$time * aggrwcws$Demand
    aggrwcws$ProductMixLoss <- (1 - (aggrwcws$ActualWC/aggrwcws$DesgnWC)) * 100
    aggrwcws$ProductMixLoss <- round(aggrwcws$ProductMixLoss)
    aggrwcwsref <- data.frame(Workstation = aggrwcws$Workstation, Product_Mix_Loss = paste(aggrwcws$ProductMixLoss,"%"))
    aggrwcwsplot <- data.frame(Workstation = aggrwcws$Workstation, Product_Mix_Loss = aggrwcws$ProductMixLoss)

Merge two data frames to get number of operators by workstation.
 
    NoofOpersop.df <- merge(NoofOpersop.df,aggrwcwsref,by="Workstation")
 
 
Finally, delete unwanted data frames
 
    rm(list=c("aggrdemws.df","aggregwc.df","co.df","cycletime.df","demwspr.df","manualwc.df","maxcotime.df","maxwc.df","wc.df","wswc.df"))
 
 
The last step is to display and plot the output for the three tabs. First, render data table to display number of operators, load lattice package and plot bar chart for efficiency and product mix loss for the workstations.
 
    output$mytable <- renderDataTable({NoofOpersop.df})
 
    library(lattice)
    output$myplot <- renderPlot({
     barchart(Manpower_Efficiency~Workstation,data=NoofOpersef.df,label = NoofOpersef.df$Manpower_Efficiency,col="grey",main="Workstation Manpower Efficiency",xlab = "Workstation",ylab = "Efficiency (%)",ylim=c(0,100))
    })
 
    output$myplot1 <- renderPlot({
      barchart(Product_Mix_Loss ~ Workstation,data=aggrwcwsplot,xlab="Workstation",ylab="Product Loss Percentage",col="blue",ylim=c(0,100))
    })
 
  })

})

Datasets (.csv)

1. Workstationavailtimes.csv

Workstation	Availabletime
WS1	420
WS2	420
WS3	420
WS4	420

2. Workcontent.csv

Workstation	Product	Mean_Set-up	SD_Set-up	Mean_Manual	SD_Manual	Machine Time
WS1	A	11	3	28	11	5
WS1	B	4	1	12	3	6
WS1	C	3	1	8	2	4
WS2	A	5	2	12	10	12
WS2	B	7	3	15	5	8
WS3	A	5	3	37	26	0
WS4	A	14	5	34	6	24

3. Demand.csv

Product	Demand
A	4000
B	2500
C	1000

4. ProductChangeover.csv

Workstation	From_Product	To_Product	Changeover_time
WS1	A	B	4
WS1	B	A	2
WS1	B	C	1
WS1	C	B	3
WS1	A	C	2
WS1	C	A	3
WS2	A	B	5
WS2	B	A	8

Output

The code will run the app to give three outputs on Summary, Plot of workstation manpower efficiency and Plot of workstation product mix loss.

Efficiency versus Productivity

For years, efficiency has been a curious parameter in manufacturing systems. But, there are soft factors in the shop floor that are not considered for analyzing the system in a more comprehensible manner. While, efficiency has been replaced with productivity it is still not able to come to terms with the way in which it has to be implemented in the system. Before, going ahead it is important that we understand the fundamental difference that differentiates the two manufacturing jargons with focus on their specific and generic interpretation.

Efficiency

Efficiency is measured by actual output production to aggregate capacity of resources in the manufacturing systems. Resources can be anything that is involved in the manufacturing operation, be it human workforce (or labor) and machines. Thus, we calculate efficiency as:

Efficiency = Actual output / Capacity

For illustration purpose, say capacity of a machine is 200 units per day and the output produced is 150 units per day on day 1 and 180 units per day on day 2. The efficiency of the machine is 75% and 90% for day 1 and day 2 respectively. The machine is a part of an assembly line with stated capacity of 120 units per day with assembled components into a finished product and the production for same day 1 and day 2 are 110 units and 100 units. The efficiency of the assembly line is 92% and 83% on day 1 and day 2. This doesn’t make more sense here as 75% efficiency of machine resulted in 92% efficiency of the line, while 90% efficiency resulted in 83% efficiency of the line. The actual reason for contradictory results will be explained later for better understanding. From, above formula we calculate the efficiency in a generic sense by aggregation of resources for a system and with single resource for a sub-system. The efficiency of a sub-system may not always provide the right answer, but would certainly make sense by determining efficiency of a system.

Determining the system boundaries is essential to understand what is within and beyond the capabilities of the system. The system has to operate within the stated external constraints and other resources are aggregated to work towards achieving the output by exploiting the internal constraints.  

Productivity                                                                                                               

Productivity is a rather better measure of the system performance as compared to efficiency, which is more like a performance tracker used for hiding the actual performance. Productivity of a manufacturing system is calculated as:

Productivity = Output produced / Inputs used

The productivity provides a rather better measure than efficiency which doesn’t consider the impact on output produced by the system. But, is there a consideration for difference between the system and sub-system for productivity calculation. To calculate the productivity of a sub-system, one has to build the link between sub-systems and output. If, one identifies that there is no direct value addition by a resource and it doesn’t contribute to output produced, than it is a clear indication of deficiency of the system. This where the concept of work content (or labor content) comes into the picture and provides a better understanding of direct labor, indirect labor, overheads, materials and depreciation among others. The different factors involved in the manufacturing must be tied together for a better understanding of the system capabilities. But, how does one calculate aggregate production cost and tie it down with labor, capital and materials in the system. It is wrong to compare a unit of measurement with different unit of measurement. To, avoid this it is important that everything is translated in the same scale for comparison of productivity and it is done first at aggregate level and then move down the ladder for identifying gaps in sub-systems for improvement opportunities.