Mathematical optimization that can be used for free work with Python + PuLP

Postscript (2020/05/25)

The mathematical optimization solver (COIN-CBC for 1 thread), which is included in PuLP and is used by default, does not support multithreading (at least for Windows), so it takes time to solve integer optimization problems. May take. There was an article that wrote how to install the multi-threaded version of COIN-CBC, so I will post the link. Of course, it is available for commercial use free of charge.

--Windows: make the calculation engine of Python mathematical optimization package PuLP faster (Windows version easy edition) & specify the calculation time --Mac: Use your own CBC solver with PuLP

Postscript (2020/05/11)

Another free solver called "MIPCL" seems to solve the problem faster than PuLP (the one-threaded version of COIN-CBC that comes with it).

(Reference entry) Finally a free mathematical optimization solver that can be used? MIPCL

However, MIPCL seems to have a little peculiarity in installation and grammar, so at the moment, I think that you can use MIPCL if you want speed and PuLP if you want ease. If it takes time to write and solve with PuLP, it is one way to rewrite it to MIPCL ~ ~, and considering the time and effort, it is also one way to introduce MIPCL from the beginning and get used to it ~ ~ I think.

Introduction

Such a mathematical problem called mathematical optimization or mathematical planning

Problem (1): Linear optimization problem (Linear programming problem)

\begin{alignat}{2}
 &\mbox{Maximize} 
 &\qquad x + y & \\
 &\mbox{Constraint} 
 & 2x + y &\leq 2 \\
 &
 & x + 2y &\leq 2 \\
 &
 & x &\geq 0 \\
 &
 & y &\geq 0 \\
\end{alignat}

a - コピー.png

And such a math problem

Problem (2): Integer optimization problem (integer programming problem)

\begin{alignat}{3}
 &\mbox{Minimize} 
 &\qquad \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & 
 &\qquad & \\
 &\mbox{subject to} 
 & \sum_{j \in J} x_{ij} &\leq 1 
 & &\forall i \in I \\
 &
 &\sum_{i \in I} x_{ij} &= 1
 & &\forall j \in J \\
 &
 & x_{ij} &\in \{0,1\}
 & &\forall i \in I, \quad \forall j \in J
\end{alignat}

Introducing Python's PuLP package that solves.

It is a target, but it is assumed that the objective function and the constraint expression can be described in linear form (linear expression), and the variable can be described as continuous value or discrete value or a mixture of them.

As an aside, problem (2) is sometimes called an integer linear optimization problem (integer linear programming problem). Also, problems with a mixture of continuous and discrete variables are called mixed integer optimization problems (mixed integer programming problems, mixed integer linear optimization problems, mixed integer linear programming problems).

What is PuLP

You should think of PuLP itself as a ** modeling API **, ** modeling language ** that makes it easy to code mathematical optimization problems on Python. It is different from the ** solver ** that actually solves the formula, but since the solver called COIN is included, you can solve the problem just by installing PuLP. PuLP itself supports several solvers other than COIN, so if those solvers are installed separately, you can (easily) switch the solver to call from PuLP.

reference

Qiita entry of ancestors

-Python in optimization -Introduction to solving linear programming problems with PuLP

Head family site

https://pypi.python.org/pypi/PuLP
http://pythonhosted.org/PuLP/ --https://www.coin-or.org/ (Solver COIN included)

Why Python + PuLP

It's just a ** personal opinion **, but ...

Price

Gurobi Optimizer, IBM ILOG CPLEX Optimization Studio, FICO Xpress Optimization Suite, Numerical Optimizer are expensive ... ** → PuLP is free! The same solver COIN is also free! !! ** **

Speed

GLPK and lp_solve are slow ... Microsoft Excel is slow and can handle a small number of variables and constraint expressions ... ** → The COIN included in PuLP is not bad for a free solver! You can call a commercial solver just by changing the code! !! ** **

Commercial use

It says that if you want to use SCIP for commercial purposes, please email me ... ** → PuLP & COIN is OK for commercial use! There is no obligation to publish the source code of the incorporated product! !! ** **

Model description language

AMPL, a modeling language that Mr. K of K & R, who is famous for C language, was involved in the development, and original modeling languages such as IBM ILOG CPLEX Optimization Studio, FICO Xpress Optimization Suite, and Numerical Optimizer. It's easy to understand because you can describe it. But it can't be used anywhere else, and if you want to do something complicated, you have to look at the manual, and even if you google, there is little information ... ** → PuLP can also write models like mathematical formulas! Since it is Python, it is easy to link with other algorithms, and if you don't understand it, you can get a lot of information by google! !! Since the description of the model part can be changed as it is and only the solver to be called can be changed, PuLP can actually become a common language! !! !! ** **

Development environment

The integrated development environment that comes with IBM ILOG CPLEX Optimization Studio and FICO Xpress Optimization Suite can display various information, and it's convenient to be able to click with the mouse. But I can't use it anywhere else ... ** → Because it is Python, it can be used on Visual Studio Code, Visual Studio (even for free Community edition and Express edition), Jupyter Notebook and Spyder! Works on various machines and environments! !! Anyone can co-develop and test anywhere! !! !! ** **

To researchers and students

Researchers and students may get a commercial solver cheaply or for free for research and learning purposes.

If the research content seems to be for commercial use in the future, researchers should consider using Python + PuLP from now on, considering the cost involved. As I've written many times, the solver you call from PuLP can be easily switched, so if you don't need to touch the depths of commercial solvers (such as tuning integer optimization parameters), PuLP is a good modeling API. think. Also, the amount of coding required is not so large for porting and compatibility from commercial solvers that have APIs for Python.

Students should consider using Python + PuLP (COIN) from now on, considering the use of mathematical optimization at their current part-time job or future employment or their customers. Of course, it depends on the scale of the project, but it is quite a hurdle to explain the cost-effectiveness of raising the cost of purchasing and maintaining a commercial solver. The same applies even if you are conducting research commissioned by a company or joint research. Think of options that corporate people can use.

That was my personal opinion.

I'm not sure about the recent free solver MIPCL, so I haven't touched it. I'm sorry.

Assumed environment

The following is written assuming the case of ** Anaconda ** that supports ** Python 3.7 ** on ** Windows ** (supported). Replace it appropriately on Mac or Linux.

Python 2 series (2.7 etc.) is old, so don't use it.

Installation

In advance!

If you have installed Gurobi Optimizer or simulation software in the past in addition to the Python you installed yourself, there is a possibility that the Python processing system is stuck. In that case, start the command prompt (Windows start button → Windows system tools → command prompt) and type where python or where pip to check which processing system is prioritized. ..

If the ʻAnaconda3` folder is displayed first, as shown below, it means that Anaconda has priority, as expected in this article.

C:\Users\(username)>where python
C:\Users\(username)\Anaconda3\python.exe
C:\Users\(username)\AppData\Local\Microsoft\WindowsApps\python.exe

C:\Users\(username)>where pip
C:\Users\(username)\Anaconda3\Scripts\pip.exe

C:\Users\(username)>

If something that is not the ʻAnaconda3` folder is displayed first, it means that the processing system that is not Anaconda is prioritized. Please decide for yourself whether it is okay to keep it as it is, and change the order of the items listed in "Path" of "Environment Variables" if necessary.

If the following is displayed, or if other processing systems are displayed but the ʻAnaconda3` folder is not displayed,

C:\Users\(username)>where python
information:The file with the given pattern was not found.

C:\Users\(username)>

This is the so-called Anaconda "pass does not pass" state. Probably because I didn't check "Add Anaconda to my PATH environment variable" when I installed Anaconda. (Reference) https://weblabo.oscasierra.net/python-anaconda-install-windows/ In that case, pass the path (either by yourself or by uninstalling and reinstalling Anaconda), or do the following work with Anaconda Prompt (Windows start button → Anaconda3 (64bit) → Anaconda Prompt) instead of the command prompt. Please give me.

Was good. Actual work from

Just in case, close the command prompt and reopen it (because it is necessary when the environment variables are rearranged). If Python is installed in a folder that requires writing with administrator privileges (for example, if you select "All users (requires admin privileges)" in the installation of Anaconda), select "Run as administrator". start up.

If you try to install from ** inside the company ** and get in the way of ** proxy **, you should temporarily set the proxy in the command prompt as follows (IP address and port number). Please replace with the value in your environment).

set HTTP_PROXY=111.222.111.222:3333
set HTTPS_PROXY=111.222.111.222:3333

Reference: http://d.hatena.ne.jp/showhey810/20140905/1409892787

Type pip install pulp to download and install the pulp. If ** security software ** asks for communication permission, please do so. If all goes well, you should see something like this:

C:\Users\(username)>pip install pulp
Collecting pulp
  Downloading PuLP-2.1-py3-none-any.whl (40.6 MB)
    |■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■| 40.6 MB 6.4 MB/s
Requirement already satisfied: pyparsing>=2.0.1 in c:\users\(username)\anaconda3\lib\site-packages (from pulp) (2.4.6)
Installing collected packages: pulp
Successfully installed pulp-2.1
C:\Users\(username)>

If it is already installed, you should see something like this: It's already installed, so there's no problem.

C:\Users\(username)>pip install pulp
Requirement already satisfied: pulp in c:\users\(username)\anaconda3\lib\site-packages (2.1)
Requirement already satisfied: pyparsing>=2.0.1 in c:\users\(username)\anaconda3\lib\site-packages (from pulp) (2.4.6)
C:\Users\(username)>

(Memo)

--You can also install with python -m pip install pulp. --You can update pip with pip install -U pulp or python -m pip install -U pulp.

Example 1)

Solve the opening problem (1).

\begin{alignat}{2}
 &\mbox{Maximize} 
 &\qquad x + y & \\
 &\mbox{Constraint} 
 & 2x + y &\leq 2 \\
 &
 & x + 2y &\leq 2 \\
 &
 & x &\geq 0 \\
 &
 & y &\geq 0. \\
\end{alignat}

code

`pulp_problem_1.py`


# coding: UTF-8

#linear/Import PuLP to solve integer linear optimization problems
import pulp
# sys.Import to get the maximum integer (int) that python can handle with maxsize
import sys

#Declare a mathematical optimization problem (maximization)
problem = pulp.LpProblem("Problem-1", pulp.LpMaximize)

#Declare variables (continuous)
#   *The variable object itself (x,y) and
#A string representation of a variable ("xx", "yy") To distinguish
#I dare to use the character string expression"xx", "yy"Is written
x = pulp.LpVariable("xx", 0, sys.maxsize, pulp.LpContinuous)
y = pulp.LpVariable("yy", 0, sys.maxsize, pulp.LpContinuous)

#Declare the objective function
#   *Parentheses are not required, but are listed for clarity
problem += ( x + y, "Objective" )

#Declare constraints
problem += ( 2 * x + y <= 2 , "Constraint_1" )
problem += ( x + 2 * y <= 2 , "Constraint_2" )

#Show all problem expressions
#   *The string representation of the variable is used when the problem is printed in a print statement.
print("Problem formula")
print("-" * 8)
print(problem)
print("-" * 8)

#Calculation
result_status = problem.solve()

#Display objective function value and solution (if solution is available)
print("")
print("Calculation result")
print("*" * 8)
print(f"Optimality= {pulp.LpStatus[result_status]}")
print(f"Objective function value= {pulp.value(problem.objective)}")
print(f"Solution x= {pulp.value(x)}")
print(f"　 y = {pulp.value(y)}")
print("*" * 8)

output

Problem formula
--------
Problem-1:
MAXIMIZE
1*xx + 1*yy + 0
SUBJECT TO
Constraint_1: 2 xx + yy <= 2

Constraint_2: xx + 2 yy <= 2

VARIABLES
xx <= 9.22337203685e+18 Continuous
yy <= 9.22337203685e+18 Continuous

--------

Calculation result
********
Optimality= Optimal
Objective function value= 1.33333334
Solution x= 0.66666667
　 y = 0.66666667
********

Digression (1)

Languages that can use operator overloading can easily write mathematical optimization expressions. Languages that aren't ... Java ...

Example (2)

Solve the opening problem (2).

\begin{alignat}{3}
 &\mbox{Minimize} 
 &\qquad \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & 
 &\qquad & \\
 &\mbox{subject to} 
 & \sum_{j \in J} x_{ij} &\leq 1 
 & &\forall i \in I \\
 &
 &\sum_{i \in I} x_{ij} &= 1
 & &\forall j \in J \\
 &
 & x_{ij} &\in \{0,1\}
 & &\forall i \in I, \quad \forall j \in J.
\end{alignat}

Digression (2)

You often use $ \ sum $ in mathematical optimization formulas, don't you? Mathematical optimization modeling languages are based on how you can write $ \ sum $. In that respect, PuLP is safe. By the way, in the introductory text of mathematical optimization, the set $ I $ and $ J $ in this problem is $ I: = \ {1, \ ldots, m \} $ and $ J: = \ Numbering from 1 by force, such as {1, \ ldots, n \} $, for example, constraint conditions

\sum_{j = 1}^{n} x_{ij} \leq 1 \quad \mbox{for} \;\; i = 1, \ldots, m

There is something that says, but I am confused as to what the elements of $ I $ and $ J $ represent, and which element of the number $ 2 $ is $ I $ or $ J $. It is recommended that you do not replace the set with a column of positive integers, and use the original names of the elements of the set as well. PuLP allows you to model naturally by using language specifications such as Python dictionaries and tuples.

Also, like the formula written with the $ \ LaTeX $ command in this problem, it would be nice to align the height of each formula to make it easier to see.

Also, this problem is mitigated by changing the 0-1 variable $ x_ {ij} \ in \ {0,1 \} $ to the continuous variable $ 0 \ leq x_ {ij} \ leq 1 $ (changed to a looser condition). It is well known that even if you solve a problem, there is an optimal solution for integers because of the total unimodularity of the problem, and you can get it. However, when you actually model it at work, the constraints are increasing steadily, and in most cases all unimodularity does not hold.

code

`pulp_problem_2.py`


# coding: UTF-8

#linear/Import PuLP to solve integer linear optimization problems
import pulp
#Import time to measure calculation time
import time



#Set of workers (use a list for convenience)
I = ["Mr. A", "Mr. B", "Mr. C"]

print(f"Assembly of workers I= {I}")


#Set of tasks (use a list for convenience)
J = ["Work i", "Work b", "Work c"]

print(f"Set of tasks J= {J}")


#Set of costs (temporary list) when worker i is assigned to task j
cc = [
      [ 1,  2,  3],
      [ 4,  6,  8],
      [10, 13, 16],
     ]

#Because cc is a list and the subscript is a number
#Define dictionary c, for example cc[0][0]Is c["Mr. A","Work i"]Make it accessible with
c = {} #Empty dictionary
for i in I:
    for j in J:
        c[i,j] = cc[I.index(i)][J.index(j)]

print("Cost c[i,j]: ")
for i in I:
    for j in J:
        print(f"c[{i},{j}] = {c[i,j]:2d},  ", end = "")
    print("")
print("")



#Declare a mathematical optimization problem (minimization)
problem = pulp.LpProblem("Problem-2", pulp.LpMinimize)
# pulp.LpMinimize :Minimize
# pulp.LpMaximize :Maximize


#Dictionary representing variable set
x = {} #Empty dictionary
       # x[i,j]Or x[(i,j)]so,(i,j)Read and write values using the tuple as a key

# 0-Declare one variable
for i in I:
    for j in J:
        x[i,j] = pulp.LpVariable(f"x({i},{j})", 0, 1, pulp.LpInteger)
        #On the variable label'['Or']'Or'-'For some reason'_'It changes to ...?
# lowBound,If you do not specify upBound, each-Infinity, +Infinity になる

#Inclusive notation can also be used
# x_suffixes = [(i,j) for i in I for j in J]
# x = pulp.LpVariable.dicts("x", x_suffixes, cat = pulp.LpBinary) 

# pulp.LpContinuous :Continuous variable
# pulp.LpInteger    :Integer variable
# pulp.LpBinary     : 0-1 variable


#Declare the objective function
problem += pulp.lpSum(c[i,j] * x[i,j] for i in I for j in J), "TotalCost"
# problem += sum(c[i,j] * x[i,j] for i in I for j in J)
#OK


#Declare constraints
#For each worker i, no more than one task can be assigned
for i in I:
    problem += sum(x[i,j] for j in J) <= 1, f"Constraint_leq_{i}"
    #On the constraint label'['Or']'Or'-'For some reason'_'It changes to ...?

#Exactly one worker is assigned for each task j
for j in J:
    problem += sum(x[i,j] for i in I) == 1, f"Constraint_eq_{j}"


#Show all problem expressions
print("Problem formula")
print(f"-" * 8)
print(problem)
print(f"-" * 8)
print("")



#Calculation
#Solver designation
solver = pulp.PULP_CBC_CMD()
# pulp.PULP_CBC_CMD() :Coin attached to PuLP-CBC
# pulp.GUROBI_CMD()   :Launch Gurobi from the command line(.Temporarily generate lp file)
# pulp.GUROBI()       :Launch Gurobi from the library(Library location required)
#Supports several other solvers
# (Example of use)
# if pulp.GUROBI_CMD().available():
#     solver = pulp.GUROBI_CMD()

#Start time measurement
time_start = time.perf_counter()

result_status = problem.solve(solver)
# solve()of()You can specify the solver within
#Pulp if nothing is specified.PULP_CBC_CMD()

#Time measurement finished
time_stop = time.perf_counter()



#Display objective function value and solution (if solution is available)
print("Calculation result")
print(f"*" * 8)
print(f"Optimality= {pulp.LpStatus[result_status]}, ", end="")
print(f"Objective function value= {pulp.value(problem.objective)}, ", end="")
print(f"Calculation time= {time_stop - time_start:.3f} (Seconds)")
print("Solution x[i,j]: ")
for i in I:
    for j in J:
        print(f"{x[i,j].name} = {x[i,j].value()},  ", end="")
    print("")
print(f"*" * 8)

output

Assembly of workers I= ['Mr. A', 'Mr. B', 'Mr. C']
Set of tasks J= ['Work i', 'Work b', 'Work c']
Cost c[i,j]:
c[Mr. A,Work i] =  1,  c[Mr. A,Work b] =  2,  c[Mr. A,Work c] =  3,  
c[Mr. B,Work i] =  4,  c[Mr. B,Work b] =  6,  c[Mr. B,Work c] =  8,  
c[Mr. C,Work i] = 10,  c[Mr. C,Work b] = 13,  c[Mr. C,Work c] = 16,  

Problem formula
--------
Problem-2:
MINIMIZE
1*x(Mr. A,Work i) + 3*x(Mr. A,Work c) + 2*x(Mr. A,Work b) + 4*x(Mr. B,Work i) + 8*x(Mr. B,Work c) + 6*x(Mr. B,Work b) + 10*x(Mr. C,Work i) + 16*x(Mr. C,Work c) + 13*x(Mr. C,Work b) + 0
SUBJECT TO
Constraint_leq_Mr. A: x(Mr. A,Work i) + x(Mr. A,Work c) + x(Mr. A,Work b) <= 1

Constraint_leq_Mr. B: x(Mr. B,Work i) + x(Mr. B,Work c) + x(Mr. B,Work b) <= 1

Constraint_leq_Mr. C: x(Mr. C,Work i) + x(Mr. C,Work c) + x(Mr. C,Work b) <= 1

Constraint_eq_Work i: x(Mr. A,Work i) + x(Mr. B,Work i) + x(Mr. C,Work i) = 1

Constraint_eq_Work b: x(Mr. A,Work b) + x(Mr. B,Work b) + x(Mr. C,Work b) = 1

Constraint_eq_Work c: x(Mr. A,Work c) + x(Mr. B,Work c) + x(Mr. C,Work c) = 1

VARIABLES
0 <= x(Mr. A,Work i) <= 1 Integer
0 <= x(Mr. A,Work c) <= 1 Integer
0 <= x(Mr. A,Work b) <= 1 Integer
0 <= x(Mr. B,Work c) <= 1 Integer
0 <= x(Mr. B,Work b) <= 1 Integer
0 <= x(Mr. C,Work i) <= 1 Integer
0 <= x(Mr. C,Work c) <= 1 Integer
0 <= x(Mr. C,Work b) <= 1 Integer

--------

Calculation result
********
Optimality= Optimal,Objective function value= 19.0,Calculation time= 0.040 (Seconds)
Solution x[i,j]:
x(Mr. A,Work i) = 0.0,  x(Mr. A,Work b) = 0.0,  x(Mr. A,Work c) = 1.0,
x(Mr. B,Work i) = 0.0,  x(Mr. B,Work b) = 1.0,  x(Mr. B,Work c) = 0.0,
x(Mr. C,Work i) = 1.0,  x(Mr. C,Work b) = 0.0,  x(Mr. C,Work c) = 0.0,
********

in conclusion

I think it's easier than Bob's painting class, so please try Python + PuLP. If the problem you want to solve ** takes a long time to calculate, see the multithreaded solver article ** with the link at the top and give it a try.