[PYTHON] Optimization learned with OR-Tools [Linear programming: multi-stage model]

Formulation: determination variable

It seems to be a determinant of how much each oil is mixed each month, but this is not enough. When using oil, you can use what you bought or what you had in stock, so we distinguish between them. You also need to know how much you can buy for inventory, given the total limit of 1,000 tons.

Therefore, we need at least three determinants for the oil set $ O $ and the month set $ M $.

x_{i,j} \geq 0 \quad \forall i \in O, \forall j \in M \quad purchase\\
y_{i,j} \geq 0 \quad \forall i \in O, \forall j \in M \quad mixed\\
z_{i,j} \geq 0 \quad \forall i \in O, \forall j \in M \quad stock\\

$ x_ {i, j} $ is the amount of oil $ i $ purchased in the month $ j $, $ y_ {i, j} $ is the oil $ i $ mixed in the month $ j $ to make soap What is the amount of oil $ i $ that $ z_ {i, j} $ has at the beginning of the month $ j $? Represents. Also, although not strictly required, we also provide a variable for how many tons of soap were made each month. This makes it easier to describe constraints. Monthly $ j $ soap production

t_j \quad \forall j \in M

will do.

Formulation: Constraint

You need to specify how the inventory of each oil fluctuates each month. Therefore

z_{i,j} + x_{i,j} - y_{i,j}= z_{i,j+1} \quad \forall i \in O, \forall j \in M \ {last month}

Prepare the formula. This means that the inventory at the beginning of the month plus the purchase for the current month minus the usage for the current month will be the inventory for the next month.

Since there are lower and upper limits for the capacity of the warehouse each month,

C_{min} \leq \sum_{i} z_{i,j} \leq C_{max} \quad  \forall j \in M

Next is the mixing constraint. There was a regulation about the range of fatty acids in the finished product. We will use the monthly production $ t_j = \ sum_ {i} y_ {i, j} \ quad \ forall j \ in M $. There is a set of fatty acids $ A $, and the ratio of the specified range $ [l_k, u_k] $ for each fatty acid $ k \ in A $ to the fatty acid $ k $ contained in the oil $ i $ $ p_ {i, k} $ Think about. Then, two formulas, the lower limit and the upper limit, are created for the range of fatty acid content.

\sum_i y_{i,j} p_{i,k} \geq l_k t_j \quad \forall k \in A, \forall j \in M \\
\sum_i y_{i,j} p_{i,k} \leq u_k t_j \quad \forall k \in A, \forall j \in M

Finally, we have to meet the monthly demand (although this time it is constant at 5,000t). Given the monthly demand of $ j $ $ D_j $,

t_j \geq D_j \quad \forall j \in M

Formulation: Objective function

The purpose this time was to minimize costs. That is, to minimize the sum of oil purchase and storage costs. The purchase unit price varies depending on the oil and month, but the storage cost per ton is fixed. Therefore, the objective function is

\sum_i \sum_j x_{i,j} P_{i,j} + \sum_i \sum_j z_{i,j}p

$ P_ {i, j} $ is the monthly unit price of oil $ i $, and $ p $ is the storage cost per ton ($ 5 this time).

Implementation

Here is the code to solve this example. my_or_tools.py and tableutils.py Please see Author GitHub for the contents.

blend_multi.py

from my_or_tools import SolVal, ObjVal, newSolver

def solve_model(Part,Target,Cost,Inventory,D,SC,SL):
    #Part:Table of oils and fatty acids,Target:Regulation of fatty acid content,Cost:Oil purchase cost,Inventory:Initial inventory,
    #D:Soap demand,SC:Oil storage costs,SL:Range of oils in stock

  s = newSolver('Multi-period soap blending problem')   #Declare solver
  Oils= range(len(Part))    #Number of oil types
  Periods, Acids = range(len(Cost[0])), range(len(Part[0]))     #Number of periods,Types of fatty acids
  Buy = [[s.NumVar(0,D,'') for _ in Periods] for _ in Oils]     #Determinant variable x. How many oil i to buy in month j
  Blnd = [[s.NumVar(0,D,'') for _ in Periods] for _ in Oils]    #The coefficient of determination y. How many oils i to use for month j
  Hold = [[s.NumVar(0,D,'') for _ in Periods] for _ in Oils]    #The coefficient of determination z. How many oils i have in stock per month j
  Prod = [s.NumVar(0,D,'') for _ in Periods]    #Auxiliary variables. Soap produced on month j
  CostP= [s.NumVar(0,D*1000,'') for _ in Periods]   #Auxiliary variables. Month j purchase cost
  CostS= [s.NumVar(0,D*1000,'') for _ in Periods]   #Auxiliary variables. Storage cost for month j
  Acid = [[s.NumVar(0,D*D,'') for _ in Periods] for _ in Acids] #Auxiliary variables. Amount of fatty acid k produced in month j
  for i in Oils:                             
    s.Add(Hold[i][0] == Inventory[i][0])    #Constraints. Inventory at the time of decision is given
  for j in Periods:                                      
    s.Add(Prod[j] == sum(Blnd[i][j] for i in Oils)) #Constraints. Month j production is the sum of the amount of oil used
    s.Add(Prod[j] >= D) #Constraints. Month j production needs to exceed that month's demand
    if j < Periods[-1]: #Excluding the last month
      for i in Oils:
        s.Add(Hold[i][j]+Buy[i][j]-Blnd[i][j] == Hold[i][j+1])  #Constraints.(Inventory at the beginning of the month)+(Purchase amount)-(amount to use)=(Inventory at the beginning of next month)
    s.Add(sum(Hold[i][j] for i in Oils) >= SL[0])   #Constraints. Lower limit of oil that can be stored
    s.Add(sum(Hold[i][j] for i in Oils) <= SL[1])   #Constraints. Upper limit of oil that can be stored
    for k in Acids: 
      s.Add(Acid[k][j]==sum(Blnd[i][j]*Part[i][k] for i in Oils))   #Constraints. The amount of fatty acid k produced in month j is the sum of the products of the amount of oil used and the content of fatty acid k in that oil.
      s.Add(Acid[k][j] >= Target[0][k] * Prod[j])   #Constraints. Lower limit of fatty acid k production
      s.Add(Acid[k][j] <= Target[1][k] * Prod[j])   #Constraints. Upper limit of fatty acid k production
    s.Add(CostP[j] == sum(Buy[i][j] * Cost[i][j] for i in Oils))    #Constraints. The purchase cost for month j is the sum of the product of the purchase amount and the price.
    s.Add(CostS[j] == sum(Hold[i][j] * SC for i in Oils))           #Constraints. The storage cost for month j is the sum of the product of inventory and storage cost.
  Cost_product = s.Sum(CostP[j] for j in Periods)   #Objective function. Total purchase cost
  Cost_storage = s.Sum(CostS[j] for j in Periods)   #Objective function. Total storage cost
  s.Minimize(Cost_product+Cost_storage) #Minimize the sum of total purchase costs and total storage costs
  rc = s.Solve()
  B,L,H,A = SolVal(Buy),SolVal(Blnd),SolVal(Hold),SolVal(Acid)
  CP,CS,P = SolVal(CostP),SolVal(CostS),SolVal(Prod)
  return rc,ObjVal(s),B,L,H,P,A,CP,CS

my_blend_multi.py

from blend_multi import solve_model

Content = [#Oil and fatty acids contained
    [36,20,33,6,4,0,1],
    [0,68,13,0,0,8,11],
    [0,6,0,66,16,5,7],
    [0,32,0,0,0,14,54],
    [0,0,49,3,39,7,2],
    [45,0,40,0,15,0,0,0],
    [0,0,0,0,0,28,72],
    [36,55,0,0,0,0,9],
    [12,48,34,0,4,2,0]
]

Target = [#Regulation of fatty acid content
    [13.3,23.2,17.8,3.7,4.6,8.8,23.6],
    [14.6,25.7,19.7,4.1,5.0,9.7,26.1]
]

Price = [#Monthly oil price
    [118,128,182,182,192],
    [161,152,149,156,174],
    [129,191,118,198,147],
    [103,110,167,191,108],
    [102,133,179,119,140],
    [127,100,110,135,163],
    [171,166,191,159,164],
    [171,131,200,113,191],
    [147,123,135,156,116]
]

Inventory = [#Initial inventory
    [15],[52],[193],[152],[70],[141],[43],[25],[89]
]

def main():
    import sys
    import tableutils
    m=len(Content)
    n=len(Content[0])
    k=len(Price[0])
    if len(sys.argv)<=1:
        print('Usage is main [content|target|cost|inventory|run] [seed]')
        return
    """ elif len(sys.argv)>2:
        random.seed(int(sys.argv[2])) """
    C = Content
    T = Target
    K = Price
    I = Inventory
    if sys.argv[1]=='content':  #Output table of oils and fatty acids contained
        for j in range(m):
            C[j].insert(0,'O'+str(j))
        C.insert(0,['']+['A'+str(i) for i in range(n)])
        tableutils.printmat(C,False,1)
    elif sys.argv[1]=='target': #Outputs fatty acid content specifications
        T.insert(0,['']+['A'+str(i) for i in range(n)])
        T[1].insert(0,'Min')
        T[2].insert(0,'Max') 
        tableutils.printmat(T,True,1)
    elif sys.argv[1]=='cost':   #Output monthly oil price
        for j in range(m):
            K[j].insert(0,'O'+str(j))
        K.insert(0,['']+['Month '+str(i) for i in range(k)])
        tableutils.printmat(K)
    elif sys.argv[1]=='inventory':  #Output inventory held at the time of decision making
        for j in range(m):
            I[j].insert(0,'O'+str(j))
        I.insert(0,['Oil','Held'])
        tableutils.printmat(I)
    elif sys.argv[1]=='run':    #Run solver
        Demand=5000 #Monthly soap demand
        Limits=[0,1000]   #Lower and upper limits you can have in stock
        Cost=5
        rc,Value,B,L,H,P,A,CP,CS=solve_model(C,T,K,I,Demand,Cost,Limits)
        if len(B):
            A.append([0 for l in range(len(A[0]))])
            for j in range(len(A)-1):
                for l in range(len(A[0])):
                    A[j][l] = A[j][l]/P[l]
                    A[-1][l] += A[j][l] 
            for j in range(m):
                B[j].insert(0,'O'+str(j))
                L[j].insert(0,'O'+str(j))
                H[j].insert(0,'O'+str(j))
            for l in range(n):
                A[l].insert(0,'A'+str(l))
            A[-1].insert(0,'Total')
            B.insert(0,['Buy qty']+['Month '+str(i) for i in range(k)])
            L.insert(0,['Blend qty']+['Month '+str(i) for i in range(k)])
            H.insert(0,['Hold qty']+['Month '+str(i) for i in range(k)])
            A.insert(0,['Acid %']+['Month '+str(i) for i in range(k)])
            P=[P]
            P[0].insert(0,'Prod qty')
            CP=[CP]
            CP[0].insert(0,'P. Cost')
            CS=[CS]
            CS[0].insert(0,'S. Cost')

            tableutils.printmat(B,True,1)
            tableutils.printmat(L,True,1)
            tableutils.printmat(H,True,1)
            tableutils.printmat(P,True,1)
            tableutils.printmat(CP,True,2)
            tableutils.printmat(CS,True,2)
            tableutils.printmat(A,True,1)

main()

result

Running the above code gives the following results:

Takashi was able to operate the soap factory in a good way.

development

Here are some of the possible developments of this example. ・ Monthly demand fluctuates • Prioritize maximizing profits rather than meeting demand. In this case, the selling price of the product is required. ・ Inventory level is managed for each oil, not for total ・ There may be uncertainty about the amount of fatty acids contained in certain oils. Etc.

Summary

I'm really tired this time