两级供应链一订多送的补货策略

两级供应链一订多送的补货策略

#main.py
from config import *
config  = None

def main(config):
    ETC_star = 10000000000
    Q_new = 0
    k_new = 0
    while ETC_star >= config.ETC:
        ETC_star = config.ETC
        config.m = config.m + 1
        k_star = config.k
        config.k = 0.5
        config.Q = Q(config)#Eq.(6)
        Q_new = 2*config.Q
        k_new = 2*config.k
        kn = 2*config.k

        # 绝对值函数
        def absf(x):
            if x >= 0:
                return x
            if x < 0:
                return -1 * x

        while absf(Q_new-config.Q)/Q_new>config.epsilon or absf(config.k-kn)/config.k>config.epsilon:
            config.Q = Q_new
            kn = k_new
            k_new = k_mean(config)#Eq.(7)
            config.k = k_new
            Q_new = Q(config) #Eq.(6)
        config.Q = Q_new
        config.ETC = ETC(config)#Eq.(3)
    m_star = config.m-1

    print('beta:%s,Q:%s,m:%s, k:%s,ETC:%s'%(config.beta,config.Q,m_star,config.k,config.ETC))

if __name__=='__main__':
    #load parameters
    config,unparsed = get_config()
    main(config)

#config.py
import argparse
import math
from scipy import integrate

parser = argparse.ArgumentParser()
###############################'const'##########################
arg = parser.add_argument_group('const')
arg.add_argument('--m',type=float,default=1,help='the number of batches in which the product is delivered from the supplier to the retailer in one replenishment cycle')
arg.add_argument('--k',type=float,default=0,help='safety factor')
arg.add_argument('--D',type=float,default=1000,help='the average demand per unit time of the retailer')
arg.add_argument('--P',type=float,default=2500,help='constant production rate of the supplier')
arg.add_argument('--r',type=float,default=1000,help='retailer’s reorder point')
arg.add_argument('--s',type=float,default=1000,help='the shipping point')
arg.add_argument('--A',type=float,default=50,help='the constant ordering cost per order')
arg.add_argument('--q',type=float,default=1000,help='the average batch size, q=Q/m')
arg.add_argument('--hs',type=float,default=3,help='the holding cost per unit product per unit time for the supplier')
arg.add_argument('--hr',type=float,default=5.5,help='the holding cost per unit product per unit time for the retailer, where hrXhs40 (Hill, 1999)')
arg.add_argument('--Cs',type=float,default=210,help='the supplier’s set-up cost per set-up')
arg.add_argument('--beta',type=float,default=1,help='the backorder ratio')
arg.add_argument('--pi0',type=float,default=100,help='gross marginal profit per unit')
arg.add_argument('--pi',type=float,default=50,help='penalty cost per unit short')
arg.add_argument('--F',type=float,default=25,help='fixed transportation cost per batch')
arg.add_argument('--Tt',type=float,default=0.04,help='constant transportation time')
arg.add_argument('--L',type=float,default=0.12,help='length of lead time')
arg.add_argument('--Q',type=float,default=100000000,help='the ordering number')
arg.add_argument('--ETC',type=float,default=10000000000,help='total cost')
arg.add_argument('--sita',type=float,default=30)
arg.add_argument('--epsilon',type=float,default=0.0001)

def get_config():
    config, unparsed = parser.parse_known_args()
    # config = update_config(config)

    return config, unparsed

def G(config):
    t = config.A + config.Cs + config.m*config.F
    return t
config, unparsed=get_config()
def H(config):
    a1 = config.hs
    a2 = config.hs*config.D/config.P
    a3 = config.hr/config.m
    a4 = config.hs/config.m
    a5 = 2*config.hs*config.D/(config.m*config.P)
    hm = a1-a2 + a3 - a4 +a5
    #hm = config.hs-config.hs*config.D/config.P+config.hr/config.m-config.hs(1-(2*config.D/config.P))/config.m
    #print(hm)
    return hm
print(H(config))
def fai(config):
    def Normal(x):
        result = 1 / math.sqrt(2 * math.pi) * math.exp(-x * x / 2)
        return result
    v,err = integrate.quad(Normal,-1000,config.k)
    return v
config, unparsed=get_config()
print(fai(config))

def cai(config):
    def Normal(x):
        result = 1 / math.sqrt(2 * math.pi) * math.exp(-x * x / 2)
        return result
    rl = Normal(config.k) - config.k + config.k*fai(config)
    return rl

def ETC(config):
    a1 = config.D/config.Q*(G(config)+config.m*(config.pi+config.pi0*(1-config.beta))*config.sita*math.sqrt(config.Tt)*cai(config))
    a2 = config.hr*(config.k*config.sita*math.sqrt(config.Tt)+(1-config.beta)*config.sita*math.sqrt(config.Tt)*cai(config))
    t = a1 + a2 +config.Q/2*H(config)
    return t

#Eq.(6)
def Q(config):
    b1 = 2*config.D*(G(config)+config.m*(config.pi+config.pi0*(1-config.beta))*config.sita*math.sqrt(config.Tt)*cai(config))
    result = math.sqrt(b1/H(config))

    return result

#根据累积概率值得到正态分布的分位数Eq.(7)
def k_mean(config):
    #b = config.hr(1-config.beta)+(config.m*config.D/config.Q)*(config.pi+config.pi0*(1-config.beta))
    b1 = config.hr*(1 - config.beta)
    b2 = config.m*config.D/config.Q
    b3 = config.pi+config.pi0*(1-config.beta)
    b = b1 + b2*b3

    quantile = 1-config.hr/b

    def Nl(x):
        def Normal(x):
            result = 1 / math.sqrt(2 * math.pi) * math.exp(-x * x / 2)
            return result

        v, err = integrate.quad(Normal, -1000, x)
        return v
    x1 = 0
    x2 = 10
    x_star = (x1+x2)/2
    while Nl(x_star)-quantile>0.0001 or Nl(x_star)-quantile<-0.0001:
        if Nl(x_star) > quantile:
            x2 = x_star
        if Nl(x_star)<quantile:
            x1 = x_star
        x_star = (x1+x2)/2
    return x_star