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

This blog is what

We will learn mathematical optimization using OR-Tools developed by Google. Most of the content is based on this book. Practical Python AI Projects: Mathematical Models of Optimization Problems with Google OR-Tools (English Edition)

Click here for a list of codes by author The code described in this blog is almost the same as the original one of the author, and it is only partially corrected to Japanese. Mathematical rigor and theoretical commentary are of low priority. Please forgive me.

This time: Linear programming-multi-stage model-

Incorporate "when" into decision making. Taking a warehouse as an example, the order quantity for next month will change depending on how many stocks are in stock at the end of this month. Similarly, the order quantity for the next month will change depending on the order quantity for the next month. In this way, consider the optimal decision-making in situations where the decision-making at one stage affects the decision-making at the subsequent stages.

Scene setting

Takashi runs a soap factory. Since the performance of the factory affects Takashi's pocket money, I decided to look for an efficient operation method.

Soap is manufactured by scouring and mixing various oils. Oils have various scents such as apricot and avocado, and each oil contains various fatty acids (lauric acid, linoleic acid, oleic acid, etc.). Below is a table of fatty acids $ A_j $ contained in oil $ i . ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/308043/4b374123-89b4-be7d-38f9-1dbbb14c7c1e.png) Also, considering the properties of the soap (cleansing power, foaming, dryness of the skin, etc.), mix properly so that the final fatty acid content is within a certain range. Here, it is assumed that it falls within the range of the table below. ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/308043/f9d26fdd-5e95-050e-1ccb-eec94d4d285f.png) Here we consider a decision over several months. Each oil can be purchased for immediate use or for the next month. The table below shows the price of each oil ( / t) for each month. image.png You can store up to 1,000 tons of oil for later use. However, there is a monthly storage cost of $ 5 per ton. In addition, Takashi's factory must meet the monthly demand for 5,000 tons of soap. At this time, Takashi's factory has the oils listed in the table below. image.png

It's finally the main subject. ** How much oil should be smelted and mixed in which month to minimize costs? </ font> ** Ah, it was long.

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: image.png

image.png

image.png

image.png

image.png

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

Recommended Posts

Optimization learned with OR-Tools [Linear programming: multi-stage model]
Optimization learned with OR-Tools [Linear programming: project management]
Optimization learned with OR-Tools [Linear programming: Let's refine oil]
Optimization learned with OR-Tools Part0 [Introduction]
Regression with linear model
Linear Programming with PuLP
Solving Mathematical Optimization Model Exercises with Google's OR-Tools (3) Production Optimization Problems
Solve Mathematical Optimization Model Exercises with Google's OR-Tools (4) Solve Number Places
[Python] Object-oriented programming learned with Pokemon
Algorithm learned with Python 9th: Linear search
[Mathematical optimization problem] Linear programming using PuLP
Production planning optimization using linear programming (Python + PuLP)
Predict hot summers with a linear regression model
Portfolio optimization with Python (Markowitz's mean variance model)
Learn Wasserstein GAN with Keras model and TensorFlow optimization