[PYTHON] Solve combinatorial optimization problems with Google's OR-Tools (5) Decide on a date course

Introduction

This article is the 5th article to solve the exercise of the reference text "Problem solving series by Python: How to make an optimization model using a data analysis library" about mathematical optimization.

The first article is below. Please see here first.

Solving mathematical optimization model exercises with Google's OR-Tools (1) The easiest mass filling problem https://qiita.com/ttlabo/private/7e8c93f387814f99931f

Let's decide a date course

This is an applied problem that appears in the reference text. Let's try the following problems at once.

: speech_balloon: Problem

ruby:8.4.problem

Dating at an amusement park with 8 attractions(8-Figure 4)。
Maximize total satisfaction within the 200 minute time limit.
You don't have to go all the date courses.
Optimize the date course that meets the conditions.

: question: ** Thinking **

■ Modeling policy Considering the time limit of 200 minutes, we will consider maximizing satisfaction as much as possible without going all around. To do this, prepare the following variables.

: point_right: Variable definition

\begin{eqnarray*}

S_i\;&\small{:=}&\;\small{Whether to choose the i-th attraction}\\
M_{ij}\;&\small{:=}&\;\small{Whether to move from the i-th attraction to the j-th attraction}\\
satisfy\;&\small{:=}&\;\small{Final satisfaction(Total)}\\
stayTime\;&\small{:=}&\;\small{Final stay time(Total)}\\
moveTime\;&\small{:=}&\;\small{Final travel time(Total)}\\
totalTime\;&\small{:=}&\;\small{Whole time(Sum of final satisfaction and staying time)}

\end{eqnarray*}

Solving the optimization gives the optimum value of these variables to determine the date course. In addition, values and information that are known in advance are defined as constants as follows.

: point_right: Constant definition

\begin{eqnarray*}

Satisfaction constant satisfactionDef\\
Stay time constant stayTimeDef\\
Move time constant moveTimeDef

\end{eqnarray*}

The order is reversed, but first set the values for the constants as follows: (1) Set the value of the satisfaction constant as follows.

ruby:8.4.renshu.py

satisfyDef = [0,50,36,45,79,55,63,71,42]

satisfyDef is an array, and the satisfaction of the i-th attraction is satisfyDef [i]. For example, the satisfaction level of the second attraction is 36. The index number of the attraction shall be counted from 0. (2) Set the value for the stay time constant as follows.

ruby:8.4.renshu.py

stayTimeDef = [0,20,28,15,35,17,18,14,22]

stayTimeDef is an array, and the staying time of the i-th attraction is stayTimeDef [i]. (3) Set the values for the movement time constant as follows.

ruby:8.4.renshu.py

moveTimeDef = {}
moveTimeDef[0,3] = 7
～

moveTimeDef is a Dictinary type, and the time to move from i to j is moveTimeDef [i, j]. For example, the travel time when moving from the 0th attraction (entrance) to the 3rd attraction (coffee cup) is moveTimeDef [0,3], which means 7 minutes.

(1) Satisfaction constant, (2) Stay time constant, and (3) Travel time constant are the values defined in the reference text.

Next, define the variables Si and Mij as follows. (1) Si is defined as follows.

ruby:8.4.renshu.py

#Variable of whether to choose an attraction
# 1:Choose/ 0:Not selected
S = {}
for i in range(9):
    S[i] = solver.BoolVar("S%i" % i)

Here, the 0th to 8th attractions are defined by either "1: Select" or "0: Do not select" variables (BoolVar type). For example, if S [4] = 1, it means to choose the 4th attraction.

(2) Mij is defined as follows.

ruby:8.4.renshu.py

#Variable of whether to move from i to j with respect to movement
# 1:Moving/ 0:Do not move
M = {}
for i in range(9):
    for j in range(9):
        M[i,j] = solver.BoolVar("M%i%i" % (i,j))

Here, each combination of attractions 0th to 8th is defined by a variable (BoolVar type) of either "1: move" or "0: do not move". There are 9 combinations x 9 combinations, 81 combinations in total. For example, if M [1,3] = 1, it means to move from the first attraction to the third attraction.

Using the above variables and constants, solve the optimization problem while setting various conditions (constraints) between the variables and constants.

: a: ** Answer **

Consider the program. The contents of the program basically follow Google's OR-Tools guide. (https://developers.google.com/optimization)

Write a spell at the beginning of the program.

ruby:8.4.renshu.py

from __future__ import print_function
from ortools.linear_solver import pywraplp

Since it is solved by the mixed integer programming solver, it is declared below.

ruby:8.4.renshu.py

# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

Next, we define constants.

ruby:8.4.renshu.py

#Constant definition=================================================================

#Attraction name
attrDef = ['','entrance','Cafe','boat','cup','Restaurant','Ferris wheel',
           'Haunted house','coaster','maze']

#Satisfaction constant
satisfyDef = [0,50,36,45,79,55,63,71,42]

#Staying time constant
stayTimeDef = [0,20,28,15,35,17,18,14,22]

The above first definition attrDef is not used in this program. This time I defined it for the time being. As I wrote in the explanation of the constant definition, satisfactionDef and stayTimeDef are arrays, and the index number and attraction number are linked. Index number 0 corresponds to attraction number 0 (entrance), satisfaction is 0, and staying time is 0. The travel time is defined as follows.

ruby:8.4.renshu.py

#Travel time constant
# moveTimeDef[i,j]← Time to move from i to j
moveTimeDef = {}
moveTimeDef[0,0] = 0
moveTimeDef[0,1] = 9
moveTimeDef[0,2] = 0
moveTimeDef[0,3] = 7
〜

Here, the travel time of the route that cannot be traveled is set to 0. In the above, only 4 routes are described, but all routes are defined in the entire source posted below. Next, define the variables.

ruby:8.4.renshu.py

#Variable definition=================================================================

#Variable of whether to choose an attraction
# 1:Choose/ 0:Not selected
S = {}
for i in range(9):
    S[i] = solver.BoolVar("S%i" % i)

#Variable of whether to move from i to j with respect to movement
# 1:Moving/ 0:Do not move
M = {}
for i in range(9):
    for j in range(9):
        M[i,j] = solver.BoolVar("M%i%i" % (i,j))

As explained earlier, S and M above are declared as variables that take 0/1. The OR-Tools solver calculates and decides which value to take.

ruby:8.4.renshu.py

#Declare satisfaction
satisfy = solver.IntVar(0.0, 1000, 'satisfy')
#Declare the length of stay
stayTime = solver.IntVar(0.0, 1000, 'stayTime')
#Travel time variable
moveTime = solver.IntVar(0.0, 1000, 'moveTime')
#Whole time variable
totalTime = solver.IntVar(0.0, 1000, 'totalTime')

Declare the above satisfaction, stay time, travel time, and total time (stay time + travel time) as variables that take integers from 0 to 100, respectively. Now that we have prepared constants and variables, let's set the conditions (set the constraint expression).

Regarding the setting of the constraint expression, the reference text has the following supplement, so we will assemble it including this condition.

ruby:8.4.problem

・ If you choose an attraction, visit it.
・ Exit when you enter the attraction.
-Connected from the entrance.(Do not end in the middle)

Below is the source of the constraints.

ruby:8.4.renshu.py

#Constraint expression===================================================================

#Be sure to pass the entrance
solver.Add(S[0] == 1)

Set the above constraint expression for the entrance. Since the entrance is always selected, set the value of S [0] to 1.

ruby:8.4.renshu.py

#Set restrictions on movement
#Set M to 0 for routes that do not move
solver.Add(M[0,0] == 0)
solver.Add(M[0,2] == 0)
solver.Add(M[0,5] == 0)
~

Set the route in the same way as the constant definition above. Since M is a variable, the solver may determine it, but here we set it explicitly. Set M [i, j] to 0 for routes that do not move. For example, you cannot move from 0th to 0th attraction, and so on from 0th to 5th attraction. Set a total of 41 routes to 0.

ruby:8.4.renshu.py

#Settings to avoid duplication(The attraction is visited only once. Cases where both are not visited are possible)
solver.Add(M[0,1] + M[1,0] <= 1)
solver.Add(M[0,3] + M[3,0] <= 1)
solver.Add(M[0,4] + M[4,0] <= 1)
~

After moving from a certain attraction i to j above, set the conditions so that the movement route will not be reversed and returned to i again. For example, if you move from the 0th attraction to the 1st attraction, M [0,1] = 1, so M [1,0] on the opposite route is set so that the value becomes 0 (does not move). I will. The same applies to the reverse movement route. Also, if M [0,1] + M [1,0] == 1, it will always move, so even if it does not move, M [0,1] + M [1,0] <= 1 I will. Set up a total of 20 routes.

ruby:8.4.renshu.py

#Conditions for M
#Do not move from one attraction to multiple attractions at the same time
for i in range(9):
    solver.Add(solver.Sum([M[i,j] for j in range(9)]) <= 1)

Set the above conditions for not moving from one attraction to multiple attractions at the same time. This is defined as that each attraction i can only take one route to each j (see the figure below).

Next, define the continuity condition of the route as follows. Route continuity means not selecting attractions that are not connected as shown in the figure below.

The constraint expression of the condition that the route is connected (continuous) is as follows.

ruby:8.4.renshu.py

#Route continuity condition
#0 starting point
step = [1,3,4]
solver.Add(solver.Sum([M[0,s] for s in step]) == solver.Sum([M[s,0] for s in step]))
#1 starting point
step = [0,2,3,4,5]
solver.Add(solver.Sum([M[1,s] for s in step]) == solver.Sum([M[s,1] for s in step]))
~

Let's look at the above starting from the first attraction. There are three routes to attraction 1: M [0,1], M [3,1], M [4,1]. As shown in the figure below (Pattern 1), let's say you enter attraction 1 by some route. Here, it is assumed that M [4,1] is entered among M [0,1], M [3,1], M [4,1](red arrow).

Only one of the three routes can be taken to enter attraction 1, so the value of any one of M [0,1], M [3,1], M [4,1] is 1, and the others The value of will be 0. Set the following as a constraint expression.

solver.Sum ([M [s, 1] for s in step]) = 1… Equation ①

Since the setting condition of the constraint expression is "Go out when you enter an attraction", consider the route to move from attraction 1 to another attraction. As shown in the above figure (Pattern 2), there are three attractions that can be moved from attraction 1, but you will have to select one from these. The constraint expression to select one from three is

solver.Sum ([M [1, s] for s in step]) = 1… Equation ②

Will be. However, since we assumed that we would take the route of M [4,1] earlier, the route that can actually be taken is either M [1,2] or M [1,5](above figure (pattern 3)). The solver calculates which route to take.

Also, it is possible that you will not be able to enter attraction 1, in which case you will not move to the next. In this case, the sum of Ms for both the incoming and outgoing routes is 0. Therefore, the conditions of the route to enter and exit attraction 1 should be as follows using equations (1) and (2).

solver.Sum([M[1,s] for s in step]) == solver.Sum([M[s,1] for s in step])

So far, we have explained the constraint expression of the continuity of the route for attraction 1. In the actual source, attractions 0 to 8 are defined in the same way. Constraint expressions for the remaining attractions are defined for the entire source listed below.

In addition, consider the relationship between the attractions you visit and the routes of entry and exit before and after. Similar to the path continuity condition, but defines the relationship between S and M as a constraint.

ruby:8.4.renshu.py

#Definition of the relationship between S and M
#0 starting point
step = [1,3,4]
solver.Add(S[0] - (solver.Sum([M[0,s] for s in step]) + solver.Sum([M[s,0] for s in step])) / 2 == 0)
#1 starting point
step = [0,2,3,4,5]
solver.Add(S[1] - (solver.Sum([M[1,s] for s in step]) + solver.Sum([M[s,1] for s in step])) / 2 == 0)
~

Regarding the above, the content of the constraint expression means "If you visit an attraction, you will come and go on the route. If you do not visit, you will not take the route." Let's look at this starting from attraction 1. First, check the following route map.

When you select attraction 1, you can consider the route to enter and the route to exit. Looking at pattern 1 for the route to enter, the following is established because it follows one of the routes of blue route M [0,1], green route M [3,1], and red route M [4,1]. To do.

solver.Sum ([M [s, 1] for s in step]) = 1… Equation ①

Similarly for the exit route, the condition of the sum of the routes of each color is

solver.Sum ([M [1, s] for s in step]) = 1… Equation ②

When you visit attraction 1, the value of S [1] is 1, so you can connect this value with equations (1) and (2) as follows.

Value of S [1]-(Formula ① + Formula ②) ÷ 2 = 0 ... Formula ③

In addition, this formula ③ holds even if you do not visit attraction 1. Since the value of S [1] is 0, and the values of equations ① and ② are also 0, respectively.

0　-　(0　+　0)　÷　2　=　0

Therefore, the relationship between S and M can be defined using equation (3), and is as follows.

S[1] - (solver.Sum([M[1,s] for s in step]) + solver.Sum([M[s,1] for s in step])) / 2 == 0

In the actual source, attractions 0 to 8 are defined in the same way. Constraint expressions for the remaining attractions are defined for the entire source listed below.

That concludes the explanation of complicated constraint expressions. The remaining constraints are defined below.

ruby:8.4.renshu.py

#Calculation of satisfaction
satisfy = solver.Sum([S[i] * satisfyDef[i] for i in range(9)])

#Calculation of staying time
stayTime = solver.Sum([S[i] * stayTimeDef[i] for i in range(9)])

#Calculation of travel time
moveTime = solver.Sum([M[i,j] * moveTimeDef[i,j] for i in range(9) for j in range(9)])

#Calculation of total time
totalTime = stayTime + moveTime

The above satisfaction and staying time are multiplied by the selected attraction S [i] value and each constant value, and the sum of them is taken. If you select an attraction, the S [i] value is 1, otherwise S [i] is 0, so if you take the sum, only the value of the selected attraction will be calculated.

For the travel time, multiply each i, j combination by the M [i, j] value and the constant value of the travel path, and take the sum of them. Also, the total time is "stay time + travel time".

Finally, the condition is "Maximize total satisfaction within the 200-minute time limit."

ruby:8.4.renshu.py

#Overall time constraints
solver.Add(totalTime <= 200)

#Maximize satisfaction
solver.Maximize(satisfy)

The above conditions for the total time of 200 minutes or less and the conditions for maximizing satisfaction are defined separately.

That's all for setting the constraint expression. Hereafter, the solution is executed.

ruby:8.4.renshu.py

status = solver.Solve()

In the following sources, if a solution is found, the S [i] value at each i and the M [i, j] value at i, j when moving will be displayed.

ruby:8.4.renshu.py

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    print("culculate Time = ", solver.WallTime(), " milliseconds")

    print('satisfy', satisfy.solution_value())
    print('stayTime', stayTime.solution_value())
    print('moveTime', moveTime.solution_value())
    print('totalTime', totalTime.solution_value())
    for i in range(9):
        print(i,S[i].solution_value())
    for i in range(9):
        for j in range(9):
            if M[i,j].solution_value() >= 0.5:
                print('M',i,j,'=',M[i,j].solution_value())

else:
    print('The problem does not have an optimal solution.')

This is the end of programming. The actual execution result is as follows.

Solution: ok Objective value = 405.0 culculate Time = 427 milliseconds satisfy 405.0 stayTime 141.0 moveTime 59.0 totalTime 200.0 0 1.0 1 1.0 2 0.0 3 1.0 4 1.0 5 1.0 6 1.0 7 1.0 8 1.0 M 0 3 = 1.0 M 1 0 = 1.0 M 3 6 = 1.0 M 4 1 = 1.0 M 5 4 = 1.0 M 6 7 = 1.0 M 7 8 = 1.0 M 8 5 = 1.0

As a result of the optimization calculation, the date course was decided. The route takes just 200 minutes and the satisfaction level is 405.

In addition, the breakdown of staying time and traveling time was found to be 141 minutes for staying time and 51 minutes for traveling time.

Looking at the attractions to select, the values other than S [2] are 1, so It is a route to visit attractions other than boats.

In addition, the route to go around is based on the value of M [i, j]. 　0 → 3 → 6 → 7 → 8 → 5 → 4 →1→0 I got the answer.

The entire program is shown below.

ruby:8.4.renshu.py

from __future__ import print_function
from ortools.linear_solver import pywraplp

# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)


#Constant definition=================================================================

#Attraction name
attrDef = ['','entrance','Cafe','boat','cup','Restaurant','Ferris wheel',
           'Haunted house','coaster','maze']

#Satisfaction constant
satisfyDef = [0,50,36,45,79,55,63,71,42]

#Staying time constant
stayTimeDef = [0,20,28,15,35,17,18,14,22]

#Travel time constant
# moveTimeDef[i,j]← Time to move from i to j
moveTimeDef = {}
moveTimeDef[0,0] = 0
moveTimeDef[0,1] = 9
moveTimeDef[0,2] = 0
moveTimeDef[0,3] = 7
moveTimeDef[0,4] = 12
moveTimeDef[0,5] = 0
moveTimeDef[0,6] = 0
moveTimeDef[0,7] = 0
moveTimeDef[0,8] = 0
moveTimeDef[1,0] = 9
moveTimeDef[1,1] = 0
moveTimeDef[1,2] = 11
moveTimeDef[1,3] = 12
moveTimeDef[1,4] = 7
moveTimeDef[1,5] = 13
moveTimeDef[1,6] = 0
moveTimeDef[1,7] = 0
moveTimeDef[1,8] = 0
moveTimeDef[2,0] = 0
moveTimeDef[2,1] = 11
moveTimeDef[2,2] = 0
moveTimeDef[2,3] = 0
moveTimeDef[2,4] = 14
moveTimeDef[2,5] = 8
moveTimeDef[2,6] = 0
moveTimeDef[2,7] = 0
moveTimeDef[2,8] = 0
moveTimeDef[3,0] = 7
moveTimeDef[3,1] = 12
moveTimeDef[3,2] = 0
moveTimeDef[3,3] = 0
moveTimeDef[3,4] = 11
moveTimeDef[3,5] = 0
moveTimeDef[3,6] = 7
moveTimeDef[3,7] = 12
moveTimeDef[3,8] = 0
moveTimeDef[4,0] = 12
moveTimeDef[4,1] = 7
moveTimeDef[4,2] = 14
moveTimeDef[4,3] = 11
moveTimeDef[4,4] = 0
moveTimeDef[4,5] = 9
moveTimeDef[4,6] = 13
moveTimeDef[4,7] = 9
moveTimeDef[4,8] = 13
moveTimeDef[5,0] = 0
moveTimeDef[5,1] = 13
moveTimeDef[5,2] = 8
moveTimeDef[5,3] = 0
moveTimeDef[5,4] = 9
moveTimeDef[5,5] = 0
moveTimeDef[5,6] = 0
moveTimeDef[5,7] = 13
moveTimeDef[5,8] = 7
moveTimeDef[6,0] = 0
moveTimeDef[6,1] = 0
moveTimeDef[6,2] = 0
moveTimeDef[6,3] = 7
moveTimeDef[6,4] = 13
moveTimeDef[6,5] = 0
moveTimeDef[6,6] = 0
moveTimeDef[6,7] = 7
moveTimeDef[6,8] = 0
moveTimeDef[7,0] = 0
moveTimeDef[7,1] = 0
moveTimeDef[7,2] = 0
moveTimeDef[7,3] = 12
moveTimeDef[7,4] = 9
moveTimeDef[7,5] = 13
moveTimeDef[7,6] = 7
moveTimeDef[7,7] = 0
moveTimeDef[7,8] = 6
moveTimeDef[8,0] = 0
moveTimeDef[8,1] = 0
moveTimeDef[8,2] = 0
moveTimeDef[8,3] = 0
moveTimeDef[8,4] = 13
moveTimeDef[8,5] = 7
moveTimeDef[8,6] = 0
moveTimeDef[8,7] = 6
moveTimeDef[8,8] = 0


#Variable definition=================================================================

#Variable of whether to choose an attraction
# 1:Choose/ 0:Not selected
S = {}
for i in range(9):
    S[i] = solver.BoolVar("S%i" % i)

#Variable of whether to move from i to j with respect to movement
# 1:Moving/ 0:Do not move
M = {}
for i in range(9):
    for j in range(9):
        M[i,j] = solver.BoolVar("M%i%i" % (i,j))

#Declare satisfaction
satisfy = solver.IntVar(0.0, 1000, 'satisfy')
#Declare the length of stay
stayTime = solver.IntVar(0.0, 1000, 'stayTime')
#Travel time variable
moveTime = solver.IntVar(0.0, 1000, 'moveTime')
#Whole time variable
totalTime = solver.IntVar(0.0, 1000, 'totalTime')


#Constraint expression===================================================================

#Be sure to pass the entrance
solver.Add(S[0] == 1)

#Set restrictions on movement
#Set M to 0 for routes that do not move
solver.Add(M[0,0] == 0)
solver.Add(M[0,2] == 0)
solver.Add(M[0,5] == 0)
solver.Add(M[0,6] == 0)
solver.Add(M[0,7] == 0)
solver.Add(M[0,8] == 0)
solver.Add(M[1,1] == 0)
solver.Add(M[1,6] == 0)
solver.Add(M[1,7] == 0)
solver.Add(M[1,8] == 0)
solver.Add(M[2,0] == 0)
solver.Add(M[2,2] == 0)
solver.Add(M[2,3] == 0)
solver.Add(M[2,6] == 0)
solver.Add(M[2,7] == 0)
solver.Add(M[2,8] == 0)
solver.Add(M[3,2] == 0)
solver.Add(M[3,3] == 0)
solver.Add(M[3,5] == 0)
solver.Add(M[3,8] == 0)
solver.Add(M[4,4] == 0)
solver.Add(M[5,0] == 0)
solver.Add(M[5,3] == 0)
solver.Add(M[5,5] == 0)
solver.Add(M[5,6] == 0)
solver.Add(M[6,0] == 0)
solver.Add(M[6,1] == 0)
solver.Add(M[6,2] == 0)
solver.Add(M[6,5] == 0)
solver.Add(M[6,6] == 0)
solver.Add(M[6,8] == 0)
solver.Add(M[7,0] == 0)
solver.Add(M[7,1] == 0)
solver.Add(M[7,2] == 0)
solver.Add(M[7,7] == 0)
solver.Add(M[8,0] == 0)
solver.Add(M[8,1] == 0)
solver.Add(M[8,2] == 0)
solver.Add(M[8,3] == 0)
solver.Add(M[8,6] == 0)
solver.Add(M[8,8] == 0)

#Settings to avoid duplication(The attraction is visited only once. Cases where both are not visited are possible)
solver.Add(M[0,1] + M[1,0] <= 1)
solver.Add(M[0,3] + M[3,0] <= 1)
solver.Add(M[0,4] + M[4,0] <= 1)
solver.Add(M[1,2] + M[2,1] <= 1)
solver.Add(M[1,3] + M[3,1] <= 1)
solver.Add(M[1,4] + M[4,1] <= 1)
solver.Add(M[1,5] + M[5,1] <= 1)
solver.Add(M[2,4] + M[4,2] <= 1)
solver.Add(M[2,5] + M[5,2] <= 1)
solver.Add(M[3,4] + M[4,3] <= 1)
solver.Add(M[3,6] + M[6,3] <= 1)
solver.Add(M[3,7] + M[7,3] <= 1)
solver.Add(M[4,5] + M[5,4] <= 1)
solver.Add(M[4,6] + M[6,4] <= 1)
solver.Add(M[4,7] + M[7,4] <= 1)
solver.Add(M[4,8] + M[8,4] <= 1)
solver.Add(M[5,7] + M[7,5] <= 1)
solver.Add(M[5,8] + M[8,5] <= 1)
solver.Add(M[6,7] + M[7,6] <= 1)
solver.Add(M[7,8] + M[8,7] <= 1)


#Conditions for M
#Do not move from one attraction to multiple attractions at the same time
for i in range(9):
    solver.Add(solver.Sum([M[i,j] for j in range(9)]) <= 1)

#Route continuity condition
#0 starting point
step = [1,3,4]
solver.Add(solver.Sum([M[0,s] for s in step]) == solver.Sum([M[s,0] for s in step]))
#1 starting point
step = [0,2,3,4,5]
solver.Add(solver.Sum([M[1,s] for s in step]) == solver.Sum([M[s,1] for s in step]))
#2 starting point
step = [1,4,5]
solver.Add(solver.Sum([M[2,s] for s in step]) == solver.Sum([M[s,2] for s in step]))
#3 starting points
step = [0,1,4,6,7]
solver.Add(solver.Sum([M[3,s] for s in step]) == solver.Sum([M[s,3] for s in step]))
#4 starting point
step = [0,1,2,3,5,6,7,8]
solver.Add(solver.Sum([M[4,s] for s in step]) == solver.Sum([M[s,4] for s in step]))
#5 starting points
step = [1,2,4,7,8]
solver.Add(solver.Sum([M[5,s] for s in step]) == solver.Sum([M[s,5] for s in step]))
#6 starting point
step = [3,4,7]
solver.Add(solver.Sum([M[6,s] for s in step]) == solver.Sum([M[s,6] for s in step]))
#7 starting point
step = [3,4,5,6,8]
solver.Add(solver.Sum([M[7,s] for s in step]) == solver.Sum([M[s,7] for s in step]))
#8 starting points
step = [4,5,7]
solver.Add(solver.Sum([M[8,s] for s in step]) == solver.Sum([M[s,8] for s in step]))

#Definition of the relationship between S and M
#0 starting point
step = [1,3,4]
solver.Add(S[0] - (solver.Sum([M[0,s] for s in step]) + solver.Sum([M[s,0] for s in step])) / 2 == 0)
#1 starting point
step = [0,2,3,4,5]
solver.Add(S[1] - (solver.Sum([M[1,s] for s in step]) + solver.Sum([M[s,1] for s in step])) / 2 == 0)
#2 starting point
step = [1,4,5]
solver.Add(S[2] - (solver.Sum([M[2,s] for s in step]) + solver.Sum([M[s,2] for s in step])) / 2 == 0)
#3 starting points
step = [0,1,4,6,7]
solver.Add(S[3] - (solver.Sum([M[3,s] for s in step]) + solver.Sum([M[s,3] for s in step])) / 2 == 0)
#4 starting point
step = [0,1,2,3,5,6,7,8]
solver.Add(S[4] - (solver.Sum([M[4,s] for s in step]) + solver.Sum([M[s,4] for s in step])) / 2 == 0)
#5 starting points
step = [1,2,4,7,8]
solver.Add(S[5] - (solver.Sum([M[5,s] for s in step]) + solver.Sum([M[s,5] for s in step])) / 2 == 0)
#6 starting point
step = [3,4,7]
solver.Add(S[6] - (solver.Sum([M[6,s] for s in step]) + solver.Sum([M[s,6] for s in step])) / 2 == 0)
#7 starting point
step = [3,4,5,6,8]
solver.Add(S[7] - (solver.Sum([M[7,s] for s in step]) + solver.Sum([M[s,7] for s in step])) / 2 == 0)
#8 starting points
step = [4,5,7]
solver.Add(S[8] - (solver.Sum([M[8,s] for s in step]) + solver.Sum([M[s,8] for s in step])) / 2 == 0)

#Calculation of satisfaction
satisfy = solver.Sum([S[i] * satisfyDef[i] for i in range(9)])

#Calculation of staying time
stayTime = solver.Sum([S[i] * stayTimeDef[i] for i in range(9)])

#Calculation of travel time
moveTime = solver.Sum([M[i,j] * moveTimeDef[i,j] for i in range(9) for j in range(9)])

#Calculation of total time
totalTime = stayTime + moveTime

#Overall time constraints
solver.Add(totalTime <= 200)

#Maximize satisfaction
solver.Maximize(satisfy)


#Solution execution
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    print("culculate Time = ", solver.WallTime(), " milliseconds")

    print('satisfy', satisfy.solution_value())
    print('stayTime', stayTime.solution_value())
    print('moveTime', moveTime.solution_value())
    print('totalTime', totalTime.solution_value())
    for i in range(9):
        print(i,S[i].solution_value())
    for i in range(9):
        for j in range(9):
            if M[i,j].solution_value() >= 0.5:
                print('M',i,j,'=',M[i,j].solution_value())

else:
    print('The problem does not have an optimal solution.')

Exhibit

This article is based on the exercises described in the reference text "Problem Solving Series with Python: How to Create an Optimization Model Using a Data Analysis Library" about mathematical optimization.

■ Reference text "Problem-solving series using Python: How to create an optimization model using a data analysis library" Tsutomu Saito [Author] Modern Science Company [Publishing]

Opinions etc.

If you have any opinions or corrections, please let us know.