<jupyter_start><jupyter_text>Traditional Model Solution<jupyter_code>import math
import numpy
from gurobipy import *
import os

p = 0.8
sigma=math.sqrt((1-p)/29/29+p/300-((1-p)/29+p/30)*((1-p)/29+p/30))
mu_model = (1-p)/29+p/30
c = 0
h_List = [182,660,1155,465,46,218,270,345,1292,2235,2770,1350,195,2129,8697,2775,172,27,900,17201,1279,21780]
b_List = [1000,1000,1000,1000,1000,1000,1000,1000,5000,1000,5000,1000,1000,5000,10000,5000,1000,1000,1000,20000,5000,40000]
L_List = [19,15,8,9,9,6,9,8,16,15,6,7,2,7,5,10,11,35,10,8,7,4]

S_MAX = 101
N = 22
Structure = [[-1],[-1],[-1],[-1],[-1],[-1],[-1],[-1],[0,1],[-1],[2,3,4],[-1],[-1],[5,6,7],[8,9,10],[11,12],[-1],[-1],[-1],[13,14,15],[16,17],[18,19,20]]
S_in_supplier = [10,10,10,10,10,10,10,10,1e20,10,1e20,10,10,1e20,1e20,1e20,10,10,10,1e20,1e20,1e20]
S_out_customer = [1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,1e20,5]

Non_Start = [] #find the set of stages with upstream stages
for j in range (N):
    if (Structure[j][0] > -1):
        Non_Start.append(j)

for i in range (N): #re-scaling
    b_List[i] = b_List[i] *100 / 100.0
    h_List[i] = h_List[i] / 100.0
            
def optimal_solution (k):
    F_S_S = []
    solution = []
    for S_in in range (0,S_MAX):
        temp = []
        for S_out in range (0,S_MAX):
            SubCost = []
            for j in range (N):
                if ((S_in+L_List[j]-S_out) >= 0):
                    SubCost.append(h_List[j]*k*sigma*math.sqrt(S_in+L_List[j]-S_out))
                else: 
                    SubCost.append(5000)
            temp.append(SubCost)
        F_S_S.append(temp)

    F_S_S = numpy.array(F_S_S)
    print("F_S_S generated!")

    #Solve MILP
    try:
        m = Model("Algo-Trad")
        # Create variables
        chi = m.addVars(N, S_MAX, S_MAX, vtype=GRB.BINARY, name="chi")
        # Set objective
        obj = 0
        for j in range (N):
            for mu in range (S_MAX):
                for nu in range (S_MAX):
                    obj += F_S_S[mu,nu,j]*chi[j,mu,nu]
        m.setObjective(obj, GRB.MINIMIZE)
        # Add constraint
        for j in range (N):
            m.addConstr(quicksum(chi[j,mu,nu] for mu in range (S_MAX) for nu in range (S_MAX)) == 1)
            if (S_out_customer[j] < 1e20):
                m.addConstr(quicksum(nu*chi[j,mu,nu] for mu in range (S_MAX) for nu in range (S_MAX)) == S_out_customer[j])
            if (S_in_supplier[j] < 1e20):
                m.addConstr(quicksum(mu*chi[j,mu,nu] for mu in range (S_MAX) for nu in range (S_MAX)) == S_in_supplier[j])
        for i in range (len(Non_Start)):
            for j in Structure[Non_Start[i]]:
                m.addConstr(quicksum(mu*chi[Non_Start[i],mu,nu] for mu in range (S_MAX) for nu in range (S_MAX)) >= quicksum(nu*chi[j,mu,nu] for mu in range (S_MAX) for nu in range (S_MAX)))

        m.optimize()

    except GurobiError as e:
        print('Error code ' + str(e.errno) + ': ' + str(e))

    except AttributeError:
        print('Encountered an attribute error')

    for j in range (N):
        for mu in range (S_MAX):
            for nu in range (S_MAX):
                if(chi[j,mu,nu].X == 1):
                    solution.append([F_S_S[mu,nu,j],mu,nu]) #obj,S_in,S_out
    
    return solution.copy()

Trad_Binomial = open("Trad-Binomial-Simulation.txt", mode='w')
Trad_Poisson = open("Trad-Poisson-Simulation.txt", mode='w')
Trad_Solution = open("Trad-Solution.txt", mode='w')
for k in [1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
         2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,
          3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.0,
         4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5.0]:
    opt_temp = optimal_solution(k)
    opt_temp = numpy.array(opt_temp)
    S_in_List = opt_temp[:,1]
    S_out_List = opt_temp[:,2]
    S_in_List = S_in_List.tolist()
    S_out_List = S_out_List.tolist()
    W_List = S_in_List
    B_List = []
    for i in range (N):
        B_List.append(round(mu_model*(S_in_List[i]+L_List[i]-S_out_List[i])+k*sigma*math.sqrt(S_in_List[i]+L_List[i]-S_out_List[i])))        
        
    Trad_Solution.write(str(k)+'\t')
    for i in range (N):
        Trad_Solution.write(str(S_in_List[i])+'\t')
    for i in range (N):
        Trad_Solution.write(str(S_out_List[i])+'\t')
    for i in range (N):
        Trad_Solution.write(str(B_List[i])+'\t')
    Trad_Solution.write('\n') 
    #simulation-binomial
    for w in range (100):
        numpy.random.seed(w+10)
        Z = numpy.random.binomial(1,p,500).tolist() #Z-29
        demand = [] #real-time demand 0 or 1
        for z in Z:
            for tempintdemand in range (28):
                demand.append(0)
            if (z == 0):
                demand.append(1)
            else:
                demand.append(0)
                demand.append(1)
        cost = numpy.zeros(N).tolist()
        backordertimes = numpy.zeros(N).tolist()
        for i in range (N):
            stock = []
            stock.append(B_List[i])
            for j in range (1, len(demand)):
                stock.append(stock[j-1])
                if ((j-S_out_List[i]) >= 0):
                    stock[j] -= demand[int(j-S_out_List[i])]
                if ((j-W_List[i]-L_List[i]) >= 0):
                    stock[j] += demand[int(j-W_List[i]-L_List[i])]
                if (stock[j] < 0):
                    cost[i] -= stock[j]*b_List[i]
                    backordertimes[i] -= stock[j]
                    stock[j] = 0
                else:
                    cost[i] += h_List[i]*stock[j]
        Trad_Binomial.write(str(k)+'\t'+str(w)+'\t')
        for i in range (N):
            Trad_Binomial.write(str(cost[i])+'\t')
        for i in range (N):
            Trad_Binomial.write(str(backordertimes[i])+'\t')
        Trad_Binomial.write('\n')
    print(k, "Simulation-Binomial Ended!")
    #simulation-poisson
    for w in range (100):
        numpy.random.seed(w+10)
        demand = numpy.random.poisson(mu_model,int(((1-p)*29+p*30)*500)).tolist()
        cost = numpy.zeros(N).tolist()
        backordertimes = numpy.zeros(N).tolist()
        for i in range (N):
            stock = []
            stock.append(B_List[i])
            for j in range (1, len(demand)):
                stock.append(stock[j-1])
                if ((j-S_out_List[i]) >= 0):
                    stock[j] -= demand[int(j-S_out_List[i])]
                if ((j-W_List[i]-L_List[i]) >= 0):
                    stock[j] += demand[int(j-W_List[i]-L_List[i])]
                if (stock[j] < 0):
                    cost[i] -= stock[j]*b_List[i]
                    backordertimes[i] -= stock[j]
                    stock[j] = 0
                else:
                    cost[i] += h_List[i]*stock[j]
        Trad_Poisson.write(str(k)+'\t'+str(w)+'\t')
        for i in range (N):
            Trad_Poisson.write(str(cost[i])+'\t')
        for i in range (N):
            Trad_Poisson.write(str(backordertimes[i])+'\t')
        Trad_Poisson.write('\n')
    print(k, "Simulation-Poisson Ended!")
Trad_Binomial.close() #every row: k,w,cost of stages 0~21, backorder times of stages 0~21
Trad_Poisson.close()
Trad_Solution.close() #every row: k, inbound service time of stages 0~21, outbound service time of stages 0~21, base stock of stages 0~21
import os
import numpy

N = 22

File = open("Trad-Solution.txt", mode='r')
temp = []
while True:
    Line = File.readline()
    if Line:
        Line = Line[:-2]
        Line = list(map(float,Line.split('\t')))
        Line = Line[1:]
        temp.append(Line)
    else:
        break
File.close()
print("inbound service time")
print(numpy.average(temp, axis=0)[0:N])
print(numpy.max(temp, axis=0)[0:N])
print(numpy.min(temp, axis=0)[0:N])
print("outbound service time")
print(numpy.average(temp, axis=0)[N:(2*N)])
print(numpy.max(temp, axis=0)[N:(2*N)])
print(numpy.min(temp, axis=0)[N:(2*N)])
print("base stock")
print(numpy.average(temp, axis=0)[(2*N):(3*N)])
print(numpy.max(temp, axis=0)[(2*N):(3*N)])
print(numpy.min(temp, axis=0)[(2*N):(3*N)])<jupyter_output>inbound service time
[10. 10. 10. 10. 10. 10. 10. 10. 25. 10.  0. 10. 10.  4.  6.  1. 10. 10.
 10. 11.  0.  1.]
[10. 10. 10. 10. 10. 10. 10. 10. 25. 10.  0. 10. 10.  4.  6.  1. 10. 10.
 10. 11.  0.  1.]
[10. 10. 10. 10. 10. 10. 10. 10. 25. 10.  0. 10. 10.  4.  6.  1. 10. 10.
 10. 11.  0.  1.]
outbound service time
[25. 25.  0.  0.  0.  4.  4.  4.  6.  6.  6.  1.  1. 11. 11. 11.  0.  0.
  1.  1.  1.  5.]
[25. 25.  0.  0.  0.  4.  4.  4.  6.  6.  6.  1.  1. 11. 11. 11.  0.  0.
  1.  1.  1.  5.]
[25. 25.  0.  0.  0.  4.  4.  4.  6.  6.  6.  1.  1. 11. 11. 11.  0.  0.
  1.  1.  1.  5.]
base stock
[0.17073171 0.         1.         1.09756098 1.09756098 1.
 1.         1.         1.90243902 1.09756098 0.         1.
 1.         0.         0.         0.         1.2195122  2.3902439
 1.09756098 1.         0.53658537 0.        ]
[1. 0. 1. 2. 2. 1. 1. 1. 2. 2. 0. 1. 1. 0. 0. 0. 2. 3. 2. 1. 1. 0.]
[0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 0. 0. 0. 1. 2. 1. 1. 0. 0.]
{'fraction_non_alphanumeric': 0.18596141051230872, 'fraction_numerical': 0.10201818585052118, 'mean_word_length': 2.049019607843137, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '249005', 'license': 'MIT', 'n_tokens_mistral': 4347, 'n_tokens_neox': 3809, 'n_words': 755, 'path': 'Exp-Bio29,30/Traditional-Model.ipynb', 'repo_name': 'ZhiyanFANG/inventory-placement-experiment', 'size': '9018'}