Hello, children's shoes. I'm Xiao Ming. A few days ago, I shared an SMT solver z3. See the link address:
https://xxmdmst.blog.csdn.net/article/details/120279521
Although the SMT solver is very powerful and can solve logic problems, Sudoku, equations and even reverse problems, there is a disadvantage that I can only find a feasible solution. If I want to find the maximum or minimum value of the feasible solution, I can't complete the planning solution function similar to Excel.
As mentioned earlier, scipy library can solve linear programming. Unfortunately, I found in the actual test this week that it does not support integer constraints and can only solve real numbers. I almost gave up writing this article, but later I found the library of PuLP, and I almost found the new world. There was such a library dedicated to planning and solving, so I spent a few days studying it.
Today, I will first demonstrate the usage of scipy library, and then demonstrate how PuLP library can solve linear programming.
Of course, PuLP library is not omnipotent. Although it can solve linear programming problems, it still has many difficult functions compared with Excel, such as nonlinear programming. If it is really necessary to carry out nonlinear programming, you can consider studying the CVXPY library. The official website address is: https://www.cvxpy.org/
Of course, for planning solutions, more than 95% of the scenarios are linear programming solutions, and PuLP is enough to deal with the scenarios we need to deal with.
Take a look at this article. Oh, this is the essence of nearly 50 hours of research. Basically, you can master 50 hours of effort in one hour. And at the end of this article, there are also technical eggs, which will help you master more simple and efficient skills.
Using scipy library for planning and solving
Download address of scipy and numpy offline documents: https://docs.scipy.org/doc/
Online document address: https://docs.scipy.org/doc/scipy/reference/
The function of SciPy library supporting linear programming is scipy.optimize.linprog
The function is declared as:
cipy.optimize.linprog(c, A_ub = None, b_ub = None, A_eq = None, b_eq = None, bounds = None, method = 'interior-point', callback = None, options = None, x0 = None)
This function is used to find the objective function under the following constraints(
c
T
x
c^{T} x
Minimum value of cTx):
A
u
b
x
≤
b
u
b
A
e
q
x
=
b
e
q
l
≤
x
≤
u
\begin{array}{r} A_{u b} x \leq b_{u b} \\ A_{e q} x=b_{e q} \\ l \leq x \leq u \end{array}
Aubx≤bubAeqx=beql≤x≤u
Description of main parameters:
- c: Unknown coefficient of objective function {one dimensional array}
- A_up: unknown coefficient of inequality constraint {two-dimensional array}
- b_up: upper bound of inequality constraint value {one dimensional array}
- A_eq: unknown coefficient of equality constraint {two dimensional array}
- B_eq: upper limit of equality constraint value {one dimensional array}
- bounds: minimum and maximum values of decision variables
Example 1:
Assume that the objective function is required z = 4 x 1 + 3 x 2 z=4x_1+3x_2 Z = maximum value of 4x1 + 3x2.
The constraints are:
{
2
x
1
+
x
2
≤
10
x
1
+
x
2
≤
8
x
2
≤
7
x
1
,
x
2
≥
0
\left\{\begin{array}{l} 2 x_{1}+x_{2} \leq 10 \\ x_{1}+x_{2} \leq 8 \\ x_{2} \leq 7 \\ x_{1}, x_{2} \geq 0 \end{array}\right.
⎩⎪⎪⎨⎪⎪⎧2x1+x2≤10x1+x2≤8x2≤7x1,x2≥0
The objective function can be understood as finding the minimum value of - z, and then converted to the parameters required by linprog. Pass in:
import numpy as np
from scipy import optimize as op
z = np.array([4, 3])
A_ub = np.array([[2, 1], [1, 1]])
b_ub = np.array([10, 8])
# The tuple distribution represents the boundary range of x1 and x2
bounds = [0, None], [0, 7]
# Find the minimum value after the inverse of the objective function
res = op.linprog(-z, A_ub, b_ub, bounds=bounds)
print(f"Maximum value of objective function z={-res.fun:.2f}，At this time, the decision variable of the objective function is{res.x.round(2)}")
Maximum value of objective function z=26.00，At this time, the decision variable of the objective function is[2. 6.]
Example 2:
Objective function: max z = 2 x 1 + 3 x 2 − 5 x 3 \max z=2 x_{1}+3 x_{2}-5 x_{3} maxz=2x1+3x2−5x3
Constraints:
{
x
1
+
x
2
+
x
3
=
7
2
x
1
−
5
x
2
+
x
3
≥
10
x
1
+
3
x
2
+
x
3
≤
12
x
1
,
x
2
,
x
3
≥
0
\left\{\begin{array}{l} x_{1}+x_{2}+x_{3}=7 \\ 2 x_{1}-5 x_{2}+x_{3} \geq 10 \\ x_{1}+3 x_{2}+x_{3} \leq 12 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right.
⎩⎪⎪⎨⎪⎪⎧x1+x2+x3=72x1−5x2+x3≥10x1+3x2+x3≤12x1,x2,x3≥0
First, convert the constraints to:
{ − 2 x 1 + 5 x 2 − x 3 ≤ − 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 + x 2 + x 3 = 7 x 1 , x 2 , x 3 ≥ 0 \left\{\begin{array}{l} -2 x_{1}+5 x_{2}-x_{3} \leq -10 \\ x_{1}+3 x_{2}+x_{3} \leq 12 \\ x_{1}+x_{2}+x_{3}=7 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧−2x1+5x2−x3≤−10x1+3x2+x3≤12x1+x2+x3=7x1,x2,x3≥0
Then write the following code:
z = np.array([2, 3, -5])
A_ub = np.array([[-2, 5, -1], [1, 3, 1]])
b_ub = np.array([-10, 12])
A_eq = np.array([[1, 1, 1]])
b_eq = np.array([7])
# Tuples represent the boundary ranges of x1, x2 and x3, respectively
bounds = [0, None], [0, None], [0, None]
# Find the minimum value after the inverse of the objective function, that is, the maximum value of the objective function
res = op.linprog(-z, A_ub, b_ub, A_eq, b_eq, bounds=bounds)
print(f"Maximum value of objective function z={-res.fun:.2f}，At this time, the decision variable of the objective function is{res.x.round(2)}")
result:
Maximum value of objective function z=14.57，At this time, the decision variable of the objective function is[6.43 0.57 0. ]
At this point, if you need to constrain variables to integers, there is nothing the scipy library can do.
Let's introduce the use of the library PuLP specially used for linear programming:
Using PuLP for programming
PULP is a library for solving linear programming (LP) problems. It describes optimization problems as mathematical models, generates MPS or LP files, and then calls LP solvers, such as CBC, GLPK, CPLEX, Gurobi, etc.
After installing the pump library, the CBC solver is available by default, and other solvers need additional installation to use.
To install the pump Library:
pip install pulp
See: https://pypi.org/project/PuLP/
Official documents: https://coin-or.github.io/pulp/
Use steps:
- Establish the optimization problem through pull.lpproblem:
prob = pulp.LpProblem(name='NoName', sense=LpMinimize)
Parameters:
- Name: the name of the problem written in the lp file;
- sense: maximum or minimum, which can be either LpMinimize\LpMaximize.
- Create an optimization variable through the LpVariable function:
pulp.LpVariable(name, lowBound=None, upBound=None, cat='Continuous', e=None)
Parameters:
- Name: variable name, which will be output to the. lp file;
- lowBound: lower bound of variable;
- upBound: upper bound of variable;
- cat: variable type, which can be one of LpInteger\LpBinary\LpContinuous (integer, binary or continuous);
- e: Indicates whether the variable exists in the objective function and constraint. It is mainly used to realize the column generation algorithm.
- Add objective functions and constraints to the optimization problem by + = for example:
prob += x1 + x2 == 100, "objective function " prob += 0.100*x1 + 0.200*x2 >= 8.0, "Constraint 1"
- Use solver to solve, and the default is CBC:
prob.solve()
If other LP solvers have been installed on this machine, you can specify solvers, such as GLPK solvers:
prob.solve(GLPK(msg = 0))
Download address of GLPK Solver: https://sourceforge.net/projects/winglpk/files/latest/download
If you need to use the default CBC to solve MIP problems, you can specify use_mps=False:
prob.solve(use_mps=False)
Let's look at how to use PuLP to solve the above linear programming problem requiring integer constraints:
Objective function:
max
z
=
2
x
1
+
3
x
2
−
5
x
3
\max z=2 x_{1}+3 x_{2}-5 x_{3}
maxz=2x1+3x2−5x3
Constraints:
{ x 1 + x 2 + x 3 = 7 2 x 1 − 5 x 2 + x 3 ≥ 10 x 1 + 3 x 2 + x 3 ≤ 12 x 1 , x 2 , x 3 ≥ 0 \left\{\begin{array}{l} x_{1}+x_{2}+x_{3}=7 \\ 2 x_{1}-5 x_{2}+x_{3} \geq 10 \\ x_{1}+3 x_{2}+x_{3} \leq 12 \\ x_{1}, x_{2}, x_{3} \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧x1+x2+x3=72x1−5x2+x3≥10x1+3x2+x3≤12x1,x2,x3≥0
The code is as follows:
from pulp import *
prob = LpProblem('max_z', sense=LpMaximize)
x1 = LpVariable('x1', 0, None, LpInteger)
x2 = LpVariable('x2', 0, None, LpInteger)
x3 = LpVariable('x3', 0, None, LpInteger)
# Set objective function
prob += 2*x1+3*x2-5*x3
# constraint condition
prob += x1+x2+x3 == 7
prob += 2*x1-5*x2+x3 >= 10
prob += x1+3*x2+x3 <= 12
status = prob.solve()
# print("solution status:", LpStatus[prob.status])
print(f"Maximum value of objective function z={value(prob.objective)}，At this time, the decision variable of the objective function is:",
{v.name: v.varValue for v in prob.variables()})
result:
Maximum value of objective function z=14.0，At this time, the decision variable of the objective function is: {'x1': 7.0, 'x2': 0.0, 'x3': 0.0}
Another example:
Production of two types of fans (fan A and fan B)
The production of fan A requires 3 hours of working hours, 4 kW of electricity and 9 tons of steel;
It takes 7 hours to produce fan B, 5 kW of electricity and 5 tons of steel.
The company can provide 300 hours of working hours, 250 kW of electricity and 420 tons of steel.
It is assumed that the unit profits of the two products are 2 million yuan and 2.1 million yuan respectively. How to arrange the production of the two products to maximize the profit?
Solution: if the output of fan A is x, the output of fan B is y and the maximum profit is Pmax, then
x,y>=0
3x+7y<=300
4x+5y<=250
9x+5y<=420
Pmax=200x+210y
Since the production quantity must be an integer, there must be an integer constraint. The code is as follows:
from pulp import *
prob = LpProblem('myProblem', sense=LpMaximize)
x = LpVariable('x', 0, 100, LpInteger)
y = LpVariable('y', 0, 100, LpInteger)
# Set objective function
prob += 200*x+210*y
# constraint condition
prob += 3*x+7*y <= 300
prob += 4*x+5*y <= 250
prob += 9*x+5*y <= 420
status = prob.solve()
# print("solution status:", LpStatus[prob.status])
print(f"Maximum value of objective function z={value(prob.objective)}，At this time, the decision variable of the objective function is:",
{v.name: v.varValue for v in prob.variables()})
Maximum value of objective function z=11460.0，At this time, the decision variable of the objective function is: {'x': 30.0, 'y': 26.0}
It can be seen that pump has solved the problem under integer constraints.
About constants involved in LpStatus:
LpStatus = {
LpStatusNotSolved:"Not Solved",
LpStatusOptimal:"Optimal",
LpStatusInfeasible:"Infeasible",
LpStatusUnbounded:"Unbounded",
LpStatusUndefined:"Undefined",
}
- Not Solved: the state before the solve() function has not been investigated.
- Optimal: the optimal solution is found.
- Impossible: there is no feasible solution to the problem (for example, constraints x < = 1 and x > = 2 are defined).
- Unbounded: the constraint condition is not bounded, and maximization will lead to infinity (for example, there is only one constraint such as x > = 3).
- Undefined: the optimal solution may exist but has not been solved.
Let's learn more cases to consolidate our knowledge:
PuLP solving linear programming case
Combination of drug ingredients to minimize total cost
A pharmaceutical factory produces three kinds of drugs a, B and C. The raw materials available are a, B, C and D, and the costs are 5 yuan, 6 yuan, 7 yuan and 8 yuan per kilogram respectively. The dosage of various drugs per kilogram of different raw materials is as follows. If the pharmaceutical factory produces exactly 100 grams of drug a, at least 530 grams of drug B and at most 160 grams of drug C every day. How to select various ingredients to meet the production and minimize the total cost?
First, we calculate the objective function and constraints.
Objective function:
min
z
=
5
x
1
+
6
x
2
+
7
x
3
+
8
x
4
\min z=5x_{1}+6x_{2}+7x_{3}+8x4
minz=5x1+6x2+7x3+8x4
Constraints:
{ x 1 + x 2 + x 3 + x 4 = 100 5 x 1 + 4 x 2 + 5 x 3 + 6 x 4 ≥ 530 2 x 1 + x 2 + x 3 + 2 x 4 ≤ 160 x 1 , x 2 , x 3 , x 4 ≥ 0 \left\{\begin{array}{l} x_{1}+x_{2}+x_{3}+x_4=100 \\ 5x_1+4x_2+5x_3+6x_4 \geq 530 \\ 2x_1+x_2+x_3+2x_4 \leq 160 \\ x_1, x_2, x_3, x_4 \geq 0 \end{array}\right. ⎩⎪⎪⎨⎪⎪⎧x1+x2+x3+x4=1005x1+4x2+5x3+6x4≥5302x1+x2+x3+2x4≤160x1,x2,x3,x4≥0
Let's use python to solve it:
from pulp import *
prob = LpProblem('Objective functions and constraints', sense=LpMinimize)
x1 = LpVariable('x1', 0, None, LpInteger)
x2 = LpVariable('x2', 0, None, LpInteger)
x3 = LpVariable('x3', 0, None, LpInteger)
x4 = LpVariable('x4', 0, None, LpInteger)
# Set objective function
prob += 5*x1+6*x2+7*x3+8*x4
# constraint condition
prob += x1+x2+x3+x4 == 100
prob += 5*x1+4*x2+5*x3+6*x4 >= 530
prob += 2*x1+x2+x3+2*x4 <= 160
print(prob)
status = prob.solve(use_mps=True)
print(f"Minimum value of objective function z={value(prob.objective)}，At this time, the decision variable of the objective function is:",
{v.name: v.varValue for v in prob.variables()})
result:
Objective functions and constraints:
MINIMIZE
5*x1 + 6*x2 + 7*x3 + 8*x4 + 0
SUBJECT TO
_C1: x1 + x2 + x3 + x4 = 100
_C2: 5 x1 + 4 x2 + 5 x3 + 6 x4 >= 530
_C3: 2 x1 + x2 + x3 + 2 x4 <= 160
VARIABLES
0 <= x1 Integer
0 <= x2 Integer
0 <= x3 Integer
0 <= x4 Integer
Minimum value of objective function z=670.0，At this time, the decision variable of the objective function is: {'x1': 30.0, 'x2': 0.0, 'x3': 40.0, 'x4': 30.0}
Distribution of nurses per day in the hospital
A new hospital in a county needs to be equipped with nurses according to the requirements of each department. The minimum number of nurses from Monday to Sunday is 34, 25, 36, 30, 28, 31 and 32. According to the regulations, a nurse should work continuously for 5 days a week. How many nurses does this hospital need at least?
answer:
A nurse has to work five consecutive days a week. We assume that the number of nurses working continuously from Monday is x1, Tuesday is x2,..., and so on. On Monday, except for the nurses who start working from Tuesday and Wednesday, other nurses will work on this day, and the same is true on other days. Then we can use linear programming to solve this simple problem:
from pulp import *
# Set the number of nurses working from Monday to Sunday
x = [LpVariable(f"x{i}", lowBound=0, upBound=18, cat=LpInteger)
for i in range(1, 8)]
min_nums = [34, 25, 36, 30, 28, 31, 32]
prob = LpProblem('Objective functions and constraints', sense=LpMinimize)
prob += lpSum(x)
for i, num in enumerate(x):
prob += lpSum(x)-x[(i+5) % 7]-x[(i+6) % 7] >= min_nums[i]
status = prob.solve()
print("Minimum number of nurses z=", value(prob.objective))
print("The number of nurses starting work from Monday to Sunday are:", [value(x[i]) for i in range(7)])
print("The number of people working from Monday to Sunday is:", [
value(lpSum(x)-x[(i+5) % 7]-x[(i+6) % 7]) for i in range(7)])
Minimum number of nurses z= 44.0 The number of nurses starting work from Monday to Sunday are: [8.0, 0.0, 14.0, 1.0, 12.0, 0.0, 9.0] The number of people working from Monday to Sunday is: [35.0, 27.0, 36.0, 30.0, 29.0, 31.0, 32.0]
Conclusion at least 44 nurses are needed. Of course, the distribution of nurses is not the only solution. There are many distribution methods. One of the feasible distribution methods is given above.
Data envelopment analysis (DEA)
Data envelopment analysis (DEA) is a new cross research field of operations research, management science and mathematical economics. It is a quantitative analysis method to evaluate the relative effectiveness of comparable units of the same type by using the method of linear programming according to multiple input indicators and multiple output indicators.
For example, there are six branches in a certain region of a bank, and their input and output are as follows:
Now it is necessary to analyze whether Banking Office 1 can achieve the effect of less input and more output through the combination of other banking offices.
This problem is actually quite complex. We can first see how to solve this problem with Excel.
Firstly, the Excel table is supplemented as follows:
(input) the formula of B11 cell is: = VLOOKUP($A11,$A:$F,COLUMN(),0)*$H
(output) the formula of D11 cell is: = VLOOKUP($A11,$A:$F,COLUMN(),0)
(input-output after combination) the formula of cell B12 is: = SUMPRODUCT(B:B,$H:$H)
The sum of x is set as the sum of the weights of six banking offices, and the formula of H2 cell is: = SUM(H3:H8)
Then set the solver parameters:
Final result:
So how do we use Python to solve this problem:
Read data first:
import pandas as pd
import numpy as np
df = pd.read_excel(r"C:\Users\ASUS\Desktop\Programming solution.xlsx",
index_col=0, usecols="A:F", sheet_name=1,
skiprows=1, nrows=6)
data = df.values
data
array([[ 20, 149, 1300, 636, 1570],
[ 18, 152, 1500, 737, 1730],
[ 23, 140, 1500, 659, 1320],
[ 22, 142, 1500, 635, 1420],
[ 22, 129, 1200, 626, 1660],
[ 25, 142, 1600, 775, 1590]], dtype=int64)
Then start solving:
from pulp import *
prob = LpProblem('Objective functions and constraints', sense=LpMinimize)
# Define 6 variables
x = [LpVariable(f"x{i}", lowBound=0, upBound=1) for i in range(1, 7)]
# Define expectation E
e = LpVariable("e", lowBound=0, upBound=1)
# Data for office 1
office1 = data[0]
# Weighted average of offices
office_wavg = np.sum(data*np.array(x)[:, None], axis=0)
# Define the target variable and expect it to be stored in the e variable
prob += e
# The sum of weights is 1
prob += lpSum(x) == 1
# Less investment
for i in range(2):
prob += office_wavg[i] <= office1[i]*e
# More output
for i in range(2, data.shape[1]):
prob += office_wavg[i] >= office1[i]
print(prob)
status = prob.solve()
print("Solution state:", LpStatus[prob.status])
print(f"Minimum value of objective function z={value(prob.objective)}，At this time, the decision variable of the objective function is:",
{v.name: v.varValue for v in prob.variables()})
print("Combined inputs and outputs:",[f"{value(office_wavg[i]):.2f}" for i in range(data.shape[1])])
result:
Objective functions and constraints:
MINIMIZE
1*e + 0
SUBJECT TO
_C1: x1 + x2 + x3 + x4 + x5 + x6 = 1
_C2: - 20 e + 20 x1 + 18 x2 + 23 x3 + 22 x4 + 22 x5 + 25 x6 <= 0
_C3: - 149 e + 149 x1 + 152 x2 + 140 x3 + 142 x4 + 129 x5 + 142 x6 <= 0
_C4: 1300 x1 + 1500 x2 + 1500 x3 + 1500 x4 + 1200 x5 + 1600 x6 >= 1300
_C5: 636 x1 + 737 x2 + 659 x3 + 635 x4 + 626 x5 + 775 x6 >= 636
_C6: 1570 x1 + 1730 x2 + 1320 x3 + 1420 x4 + 1660 x5 + 1590 x6 >= 1570
VARIABLES
e <= 1 Continuous
x1 <= 1 Continuous
x2 <= 1 Continuous
x3 <= 1 Continuous
x4 <= 1 Continuous
x5 <= 1 Continuous
x6 <= 1 Continuous
Solution state: Optimal
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
It can be seen that the solution result is consistent with Excel.
Of course, consistent results can be obtained for any banking office:
from pulp import *
# Define 6 variables
x = [LpVariable(f"x{i}", lowBound=0, upBound=1) for i in range(1, 7)]
# Define expectation E
e = LpVariable("e", lowBound=0, upBound=1)
# Data for office 1
for num, office in enumerate(data, 1):
prob = LpProblem('Objective functions and constraints', sense=LpMinimize)
# Weighted average of offices
office_wavg = np.sum(data*np.array(x)[:, None], axis=0)
# Define the target variable and expect it to be stored in the e variable
prob += e
# The sum of weights is 1
prob += lpSum(x) == 1
# Less investment
for i in range(2):
prob += office_wavg[i] <= office1[i]*e
# More output
for i in range(2, data.shape[1]):
prob += office_wavg[i] >= office1[i]
status = prob.solve()
print(f"---------Office?{num}----------")
print("Solution state:", LpStatus[prob.status])
print(f"Minimum value of objective function z={value(prob.objective)}，At this time, the decision variable of the objective function is:",
{v.name: v.varValue for v in prob.variables()})
print("Combined inputs and outputs:", [
f"{value(office_wavg[i]):.2f}" for i in range(data.shape[1])])
---------Office 1----------
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------Office 2----------
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------Office 3----------
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------Office 4----------
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------Office 5----------
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
---------Office 6----------
Minimum value of objective function z=0.96780303，At this time, the decision variable of the objective function is: {'e': 0.96780303, 'x1': 0.0, 'x2': 0.66098485, 'x3': 0.0, 'x4': 0.0, 'x5': 0.33901515, 'x6': 0.0}
Combined inputs and outputs: ['19.36', '144.20', '1398.30', '699.37', '1706.27']
From the results, all the results are consistent.
traveling salesman problem
A postman starts from the post office (located in community a) every day, needs to deliver letters to 7 communities in different locations, and then collects letters from the mailbox in these 7 places and returns to the post office. Through long-time observation and recording, the postman sorted out the average time required for cycling between the seven areas, as shown in the table below. How to plan a postman's route for a day to minimize the time spent on the road?
This problem can be solved by violence in Python, and the thinking is very direct and simple.
However, our main goal today is to explain in depth how Python solves linear programming, so let's demonstrate the use of binary variables in linear programming.
First, let's try using basic constraints:
- It can't be diagonal
- Origin uniqueness
- Goal uniqueness
from pulp import *
import numpy as np
import pandas as pd
distances = np.array([
[0, 19, 20, 16, 25, 20, 19],
[19, 0, 20, 17, 11, 12, 8],
[20, 20, 0, 25, 17, 18, 20],
[16, 17, 25, 0, 11, 8, 9],
[25, 11, 17, 11, 0, 13, 14],
[20, 12, 18, 8, 13, 0, 15],
[19, 8, 20, 9, 14, 15, 0]
])
# Create a binary variable. i represents the column (starting point), j represents the row (target), and xij represents whether the current path is selected
x = np.array([
[LpVariable(f"{i}{j}", cat=LpBinary)
for i in "ABCDEFG"]
for j in "ABCDEFG"
])
prob = LpProblem('Objective functions and constraints', sense=LpMinimize)
# Minimum time required
z = (distances*x).sum()
prob += z
# It can't be diagonal, that is, from yourself to yourself
for i in range(x.shape[0]):
prob += x[i, i] == 0
# Origin uniqueness
for i in x.sum(axis=0):
prob += i == 1
# Goal uniqueness
for i in x.sum(axis=1):
prob += i == 1
status = prob.solve()
print("z=", value(prob.objective))
print("Path selection:")
result = [[int(value(cell)) for cell in row]for row in x]
result = pd.DataFrame(result, index=list("ABCDEFG"), columns=list("ABCDEFG"))
print(result)
result:
z= 88.0 Path selection: A B C D E F G A 0 0 1 0 0 0 0 B 0 0 0 0 0 0 1 C 1 0 0 0 0 0 0 D 0 0 0 0 1 0 0 E 0 0 0 0 0 1 0 F 0 0 0 1 0 0 0 G 0 1 0 0 0 0 0
It can be seen that a < - > C, B < - > G, D - > F - > e - > d appear in the solution result path.
How to add constraints to avoid the problem of circular path?
hypothesis
U
i
U_i
Ui ＾ indicates that the i-th cell is the i-th one. Theoretical research (I don't understand its mathematical principle) shows that when the constraint is:
{
u
i
−
u
j
+
n
x
i
j
≤
n
−
1
1
≤
i
≠
j
≤
n
−
1
\left\{\begin{array}{l} u_i-u_j+nx_{ij} \leq n-1 \\ 1 \leq i \neq j \leq n-1 \end{array}\right.
{ui−uj+nxij≤n−11≤i=j≤n−1
Where n represents the number of destinations, i.e. 7,
x
i
j
x_{ij}
xij ， indicates whether the specified path is selected. There are two possibilities: 0 or 1.
To avoid the above loop path problem, we continue to add the above constraints to our program to test:
paths = [LpVariable(f"x{_}", lowBound=0, upBound=6,
cat=LpInteger) for _ in range(7)]
prob += paths[0] == 1
for j in range(1, 7):
for i in range(1, 7):
if i == j:
continue
prob += paths[i]-paths[j]+7*x[j, i] <= 6
status = prob.solve()
print("z=", value(prob.objective))
print("Arrival order:", [int(value(path)) for path in paths])
print("Path selection:")
result = [[int(value(cell)) for cell in row]for row in x]
result = pd.DataFrame(result, index=list("ABCDEFG"), columns=list("ABCDEFG"))
print(result)
result:
z= 93.0 Arrival order: [1, 2, 0, 5, 1, 6, 3] Path selection: A B C D E F G A 0 0 0 0 0 1 0 B 0 0 0 0 1 0 0 C 1 0 0 0 0 0 0 D 0 0 0 0 0 0 1 E 0 0 1 0 0 0 0 F 0 0 0 1 0 0 0 G 0 1 0 0 0 0 0
It can be seen from the results that the minimum time is 93 and the path is a - > C - > e - > b - > G - > D - > F - > a.
The complete code is:
from pulp import *
import numpy as np
import pandas as pd
distances = np.array([
[0, 19, 20, 16, 25, 20, 19],
[19, 0, 20, 17, 11, 12, 8],
[20, 20, 0, 25, 17, 18, 20],
[16, 17, 25, 0, 11, 8, 9],
[25, 11, 17, 11, 0, 13, 14],
[20, 12, 18, 8, 13, 0, 15],
[19, 8, 20, 9, 14, 15, 0]
])
# i represents the column (starting point) and j represents the row (target)
x = np.array([
[LpVariable(f"{i}{j}", cat=LpBinary)
for i in "ABCDEFG"]
for j in "ABCDEFG"
])
prob = LpProblem('Objective functions and constraints', sense=LpMinimize)
# Minimum time required
z = (distances*x).sum()
prob += z
# It can't be diagonal, that is, from yourself to yourself
for i in range(x.shape[0]):
prob += x[i, i] == 0
# Origin uniqueness
for i in x.sum(axis=0):
prob += i == 1
# Goal uniqueness
for i in x.sum(axis=1):
prob += i == 1
paths = [LpVariable(f"x{_}", lowBound=0, upBound=6,
cat=LpInteger) for _ in range(7)]
prob += paths[0] == 1
for j in range(1, 7):
for i in range(1, 7):
if i == j:
continue
prob += paths[i]-paths[j]+7*x[j, i] <= 6
status = prob.solve()
print("z=", value(prob.objective))
print("Arrival order:", [int(value(path)) for path in paths])
print("Path selection:")
result = [[int(value(cell)) for cell in row]for row in x]
result = pd.DataFrame(result, index=list("ABCDEFG"), columns=list("ABCDEFG"))
print(result)
So is the result of linear programming correct? We can obtain the most accurate minimum value through program traversal. Compare:
import itertools
import pandas as pd
import numpy as np
distances = np.array([
[0, 19, 20, 16, 25, 20, 19],
[19, 0, 20, 17, 11, 12, 8],
[20, 20, 0, 25, 17, 18, 20],
[16, 17, 25, 0, 11, 8, 9],
[25, 11, 17, 11, 0, 13, 14],
[20, 12, 18, 8, 13, 0, 15],
[19, 8, 20, 9, 14, 15, 0]
])
distances = pd.DataFrame(distances, index=list(
"ABCDEFG"), columns=list("ABCDEFG"))
min_time = distances.sum().sum()
min_paths = []
for path in itertools.permutations("BCDEFG", 6):
path = ("A",)+path+("A",)
distance_cur = 0
for src, dest in zip(path[:-1], path[1:]):
distance_cur += distances.at[dest, src]
if min_time > distance_cur:
min_time = distance_cur
min_paths.clear()
if min_time >= distance_cur:
min_paths.append(path)
print("The shortest paths are:")
for min_path in min_paths:
print("->".join(min_path))
print("The minimum time is:", min_time)
The shortest paths are: A->C->E->B->G->D->F->A A->F->D->G->B->E->C->A Minimum time: 93
From the results, we can see that the linear rule correctly solves one of the correct answers.
Egg: special skills enable the Z3 solver to find multiple solutions
Thank you for reading this article carefully. I'll give you a skill so that Z3 constraint solver can also obtain all feasible solutions.
When we introduced the z3 Solver (see the address at the beginning of the article for details), we mentioned that one disadvantage of z3 is that we can only find one feasible solution and can not find all feasible solutions. But in fact, as long as we use a little skill, we can let z3 find all the feasible solutions. Idea: every time a solution is found, it will be added to the constraint. The result is not allowed to meet this solution. Until no feasible solution can be found, all feasible solutions can be found.
Take the previous eight queens problem as an example. First, let's look at the results of the ordinary backtracking algorithm:
result = [0] * 8 # The subscript represents the number of rows where the queen is located, and the value stores the column where the queen is located
def print_eight_queen(result):
for column in result:
for i in range(8):
if i == column:
print(end="Q ")
else:
print(end="* ")
print()
n = 1 # The nth condition is satisfied
def cal8queen(row: int = 0):
global n
if row == 8:
print(f"----------{n}----------")
print_eight_queen(result)
n += 1
for column in range(8):
# Only when the row column placement does not conflict with the previous placement
if is_ok(row, column):
result[row] = column
cal8queen(row + 1)
def is_ok(row, column):
"""
Check on page row Place the queen on page column Whether the column meets the conditions
The condition for passing the verification is row All rows in front of the row do not have the same number of columns as the current column, nor are they on the diagonal of the current column
"""
leftup, rightup = column - 1, column + 1
for i in range(row - 1, -1, -1):
if column == result[i] or leftup == result[i] or rightup == result[i]:
return False
leftup -= 1
rightup += 1
return True
cal8queen(0)
The results of intercepting the beginning and end are printed as follows:
----------1---------- Q * * * * * * * * * * * Q * * * * * * * * * * Q * * * * * Q * * * * Q * * * * * * * * * * * Q * * Q * * * * * * * * * Q * * * * ... ----------92---------- * * * * * * * Q * * * Q * * * * Q * * * * * * * * * Q * * * * * * * * * * Q * * * Q * * * * * * * * * * * * Q * * * * * Q * * *
Then we let the z3 constraint solver help us find these 92 solutions:
from z3 import *
# Each queen must record the corresponding column position of the queen in different rows
Q = [Int(f'Q_{i}') for i in range(8)]
# Each queen is in columns 0,1,2,..., 7
val_c = [And(0 <= Q[i], Q[i] <= 7) for i in range(8)]
# Up to one queen per column
col_c = [Distinct(Q)]
# Diagonal constraint
diag_c = [If(i == j,
True,
And(Q[i] - Q[j] != i - j, Q[i] - Q[j] != j - i))
for i in range(8) for j in range(i)]
def print_eight_queen(result):
for column in result:
for i in range(8):
if i == column:
print(end="Q ")
else:
print(end="* ")
print()
s = Solver()
s.add(val_c + col_c + diag_c)
n = 1 # The nth condition is satisfied
while s.check() == sat:
result = s.model()
result = [result[Q[i]].as_long() for i in range(8)]
print(f"----------{n}----------")
print_eight_queen(result)
s.add(Not(And([qx == p for qx, p in zip(Q, result)])))
n += 1
----------1---------- * * * Q * * * * * Q * * * * * * * * * * * * * Q * * * * * Q * * Q * * * * * * * * * Q * * * * * * * * * Q * * * * * * * * * Q * ... ----------92---------- * * * * * * Q * * Q * * * * * * * * * * * Q * * * * Q * * * * * Q * * * * * * * * * * Q * * * * * * * * * * * Q * * * * Q * * *
You can see that it is also the solution of 92. Although the order is different, these 92 solutions have been found.