OR-Tools is fast and portable software for combinatorial optimization.
About OR-Tools
OR-Tools is an open source software suite for optimization, tuned for tackling the world's toughest problems in vehicle routing, flows, integer and linear programming, and constraint programming.
After modeling your problem in the programming language of your choice, you can use any of a half dozen solvers to solve it: commercial solvers such as Gurobi or CPLEX, or open-source solvers such as SCIP, GLPK, or Google's GLOP and award-winning CP-SAT.
A set of techniques for finding feasible solutions to a problem expressed as constraints (e.g., a room can't be used for two events simultaneously, or the distance to the crops must be less than the length of the hose, or no more than five TV shows can be recorded at once).
The Glop linear optimizer finds the optimal value of a linear objective function, given a set of linear inequalities as constraints (e.g., assigning people to jobs, or finding the best allocation of a set of resources while minimizing cost). Glop and the mixed-integer programming software SCIP are also available via the Google Apps Script Optimization Service.
Linear optimization (or linear programming) is the name given to computing the best solution to a problem modeled as a set of linear relationships. These problems arise in many scientific and engineering disciplines. (The word "programming" is a bit of a misnomer, similar to how "computer" once meant "a person who computes." Here, "programming" refers to the arrangement of a plan , rather than programming in a computer language.)
import math
from ortools.linear_solver.python import model_builder
# 这种方式不常用
model = model_builder.ModelBuilder()
x = model.new_int_var(0.0, math.inf, "x")
y = model.new_int_var(0.0, math.inf, "y")
print("Number of variables =", model.num_variables)
# x + 7 * y <= 17.5.
model.add(x + 7 * y <= 17.5)
# x <= 3.5.
model.add(x <= 3.5)
print("Number of constraints =", model.num_constraints)
# Maximize x + 10 * y.
model.maximize(x + 10 * y)
# Create the solver with the SCIP backend, and solve the model.
solver = model_builder.ModelSolver("scip")
status = solver.solve(model)
if status == model_builder.SolveStatus.OPTIMAL:
print("Solution:")
print("Objective value =", solver.objective_value)
print("x =", solver.value(x))
print("y =", solver.value(y))
else:
print("The problem does not have an optimal solution.")
# [END print_solution]
# [START advanced]
print("\nAdvanced usage:")
print("Problem solved in %f seconds" % solver.wall_time)
2.2.2 决策变量
Var(lb, ub, integer, name), 不确定类型变量的创建
NumVar(lb, ub, name), 连续型变量
IntVar(lb, ub, name), 整数型变量
BoolVar(name), 布尔型变量
def Var(self, lb: "double", ub: "double", integer: "bool", name: "std::string const &") -> "operations_research::MPVariable *":
r"""
Creates a variable with the given bounds, integrality requirement and
name. Bounds can be finite or +/- MPSolver::infinity(). The MPSolver owns
the variable (i.e. the returned pointer is borrowed). Variable names are
optional. If you give an empty name, name() will auto-generate one for you
upon request.
"""
return _pywraplp.Solver_Var(self, lb, ub, integer, name)
# 求解器的求解状态
status = solver.Solve()
if status == pywraplp.Solver.OPTIMAL
2.2.6 读取求解值
variables = {}
for i in range(N):
j = matrix[i][0]
y = 1 if matrix[i][3] > 0 else 0
variables[j, i] = solver.IntVar(y, 1000, f'x_{j}_{i}')
for key in variables.keys():
# 设置的决策变量有一个
print(variables[key].solution_value())
from ortools.linear_solver import pywraplp
N = 16
n = 6
matrix = [
[0,7,4,1],
[0,8,6,0],
[1,3,1,1],
[1,2,1,0],
[2,5,2,1],
[2,5,3,1],
[3,8,5,0],
[3,2,1,1],
[3,8,6,1],
[3,2,1,1],
[4,8,5,1],
[4,4,2,1],
[4,1,1,1],
[4,5,3,0],
[5,1,1,0],
[5,6,3,0]
]
solver = pywraplp.Solver.CreateSolver('SCIP')
variables = {}
for i in range(N):
j = matrix[i][0]
y = 1 if matrix[i][3] > 0 else 0
variables[j, i] = solver.IntVar(y, 1000, f'x_{j}_{i}')
for j in range(n):
data = []
for i in range(N):
if variables.get((j, i)):
data.append((j, i))
if len(data) > 0:
constrain = sum((variables[key] for key in data)) >= 1
solver.Add(constrain)
else: print('err')
constrain = sum((variables[j, i] * matrix[i][2] for (j, i) in variables.keys())) <= 1000
solver.Add(constrain)
solver.Maximize(sum((variables[j, i] * matrix[i][1] for (j, i) in variables.keys())))
status = solver.Solve()
0
solver.Objective().Value()
2965.0
for key in variables.keys():
print(variables[key].solution_value())
The CP-SAT solver[6] bundled with OR-Tools won a total of eleven gold medals between 2018 and 2020 in the MiniZinc Challenge,[7] an international constraint programming competition.
-- Wikipedia
MiniZinc Challenge, 约束规划求解竞赛的多次金牌玩家.
Constraint optimization, or constraint programming (CP), is the name given to identifying feasible solutions out of a very large set of candidates, where the problem can be modeled in terms of arbitrary constraints. CP problems arise in many scientific and engineering disciplines. (The word "programming" is a bit of a misnomer, similar to how "computer" once meant "a person who computes". Here, "programming" refers to the arrangement of a plan , rather than programming in a computer language.)
import pyscipopt as opt
# 创建上述数独
matrix = [
[8, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 3, 6, 0, 0, 0, 0, 0],
[0, 7, 0, 0, 9, 0, 2, 0, 0],
[0, 5, 0, 0, 0, 7, 0, 0, 0],
[0, 0, 0, 0, 4, 5, 7, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 3, 0],
[0, 0, 1, 0, 0, 0, 0, 6, 8],
[0, 0, 8, 5, 0, 0, 0, 1, 0],
[0, 9, 0, 0, 0, 0, 4, 0, 0]
]
# 普通参数
# matrix shape
N = 9
# litte matrix shape
n = 3
# 填入的数字
constants = [e for e in range(1, N + 1)]
model = opt.Model("sudoku")
# 创建变量, 需要 9 * 9 * 9 = 729个决策变量
variables = {}
for r in range(N):
for c in range(N):
for k in constants:
key = (r, c, k)
variables[key] = model.addVar(
name = '_'.join(str(e) for e in key),
# 数据类型为布尔型, 即0, 1
vtype = 'B'
)
# 同时这里实现上述模型中的决策变量 0-1约束
# 约束1, 对已经填写了数字变量坐标添加1策略约束
for r in range(N):
for c in range(N):
if matrix[r][c] != 0:
constrain = variables[(r, c, matrix[r][c])] == 1
model.addCons(constrain)
# 约束2, 每个网格内有且只有一个数字
for r in range(N):
for c in range(N):
constrain = opt.quicksum(variables[(r, c, k)] for k in constants) == 1
model.addCons(constrain)
# 约束3, 行列约束
for i in range(N):
for k in constants:
# 行约束
r_constrain = opt.quicksum(variables[i, c, k] for c in range(N)) == 1
model.addCons(r_constrain)
# 列约束
c_constrain = opt.quicksum(variables[r, i, k] for r in range(N)) == 1
model.addCons(c_constrain)
# 约束4, 子宫格约束
for row in range(n):
for col in range(n):
for k in constants:
constrain = opt.quicksum(variables[r + n * row, c + n * col, k] for r in range(n) for c in range(n)) == 1
model.addCons(constrain)
# 无目标函数
model.optimize()
if model.getStatus() == 'optimal':
print('sudoku solution:')
# 预先构造一个数独容器用于存放解析结果
solution = [[0 for _ in range(N)] for _ in range(N)]
for (r, c, k) in variables.keys():
if model.getVal(variables[r, c, k]) == 1:
solution[r][c] = k
# 绘制数独结果
sudo_data = (''.join(f'{"| " if c % 3 == 0 else " "}{solution[r][c]}{" " if matrix[r][c] == 0 else "*"}' for c in range(N)) + '|' + ('\n'+ '-' * 31 if (r + 1) % 3 == 0 and r != N - 1 else '') for r in range(N))
print('\n'.join(sudo_data))
else: print('warning: no solution')
from ortools.linear_solver import pywraplp
matrix = [
[8, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 3, 6, 0, 0, 0, 0, 0],
[0, 7, 0, 0, 9, 0, 2, 0, 0],
[0, 5, 0, 0, 0, 7, 0, 0, 0],
[0, 0, 0, 0, 4, 5, 7, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 3, 0],
[0, 0, 1, 0, 0, 0, 0, 6, 8],
[0, 0, 8, 5, 0, 0, 0, 1, 0],
[0, 9, 0, 0, 0, 0, 4, 0, 0]
]
# 普通参数
# matrix shape
N = 9
# litte matrix shape
n = 3
# 填入的数字
constants = [e for e in range(1, N + 1)]
# 使用scip求解器
solver = pywraplp.Solver.CreateSolver('SCIP')
variables = {}
for r in range(N):
for c in range(N):
for k in constants:
key = (r, c, k)
variables[key] = solver.BoolVar('_'.join(str(e) for e in key))
for r in range(N):
for c in range(N):
if matrix[r][c] != 0:
constrain = variables[(r, c, matrix[r][c])] == 1
solver.Add(constrain)
for r in range(N):
for c in range(N):
# 求和, 可以直接使用sum()函数, 也可以使用solver.Sum
constrain = sum(variables[(r, c, k)] for k in constants) == 1
solver.Add(constrain)
for i in range(N):
for k in constants:
# 行约束
r_constrain = sum(variables[i, c, k] for c in range(N)) == 1
solver.Add(r_constrain)
# 列约束
c_constrain = sum(variables[r, i, k] for r in range(N)) == 1
solver.Add(c_constrain)
for row in range(n):
for col in range(n):
for k in constants:
constrain = sum(variables[r + n * row, c + n * col, k] for r in range(n) for c in range(n)) == 1
solver.Add(constrain)
status = solver.Solve()
status
# 0
solution = [[0 for _ in range(N)] for _ in range(N)]
for (r, c, k) in variables.keys():
if variables[r, c, k].solution_value() == 1:
solution[r][c] = k
# 绘制数独结果
sudo_data = (''.join(f'{"| " if c % 3 == 0 else " "}{solution[r][c]}{" " if matrix[r][c] == 0 else "*"}' for c in range(N)) + '|' + ('\n'+ '-' * 31 if (r + 1) % 3 == 0 and r != N - 1 else '') for r in range(N))
print('\n'.join(sudo_data))
| 8* 1 2 | 7 5 3 | 6 4 9 |
| 9 4 3*| 6* 8 2 | 1 7 5 |
| 6 7* 5 | 4 9* 1 | 2* 8 3 |
-------------------------------
| 1 5* 4 | 2 3 7*| 8 9 6 |
| 3 6 9 | 8 4* 5*| 7* 2 1 |
| 2 8 7 | 1* 6 9 | 5 3* 4 |
-------------------------------
| 5 2 1*| 9 7 4 | 3 6* 8*|
| 4 3 8*| 5* 2 6 | 9 1* 7 |
| 7 9* 6 | 3 1 8 | 4* 5 2 |
4.2 约束求解
整体解决思路还是类似的, 但是由于cp_model模型提供了更为简单的方式来实现.
from ortools.sat.python import cp_model
# 在求解上稍微有点不一样, model并没有solver()
model = cp_model.CpModel()
# 需要创建一个求解器的实例
solver = cp_model.CpSolver()
variables = {}
for r in range(N):
for c in range(N):
variables[r, c] = cell if (cell := matrix[r][c]) != 0 else model.NewIntVar(1, N, f'v[{r}{c}]')
for r in range(N):
model.AddAllDifferent((variables[r, c] for c in range(N)))
for c in range(N):
model.AddAllDifferent((variables[r, c] for r in range(N)))
for row in range(n):
for col in range(n):
model.AddAllDifferent((variables[r + n * row, c + n * col] for r in range(n) for c in range(n)))
status = solver.Solve(model)
if status == cp_model.OPTIMAL:
print('sudoku solution:')
# 预先构造一个数独容器用于存放解析结果
solution = [[0 for _ in range(N)] for _ in range(N)]
for i in range(N):
for j in range(N):
solution[i][j] = int(solver.Value(variables[i, j]))
# 绘制数独结果
sudo_data = (''.join(f'{"| " if c % 3 == 0 else " "}{solution[r][c]}{" " if matrix[r][c] == 0 else "*"}' for c in range(N)) + '|' + ('\n'+ '-' * 31 if (r + 1) % 3 == 0 and r != N - 1 else '') for r in range(N))
print('\n'.join(sudo_data))
else: print('warning: no solution')
Option Explicit
'@name Eight_Queen_Algorithm
'@authot HLA
'@description
'8皇后算法
'棋盘 8x8, 将8个棋子皇后放置在特定的位置,使其无法互相攻击, 即位置, 相互不处于同一列, 行, 对角线, 反对角线
Private Const iBoundary As Long = 8
Private Const iFrontier As Long = 7
Dim arrQueen(iBoundary) As Long
Dim xPosition As Long
Dim yPosition As Long
Dim iCcount As Long '计算解法数量
Sub Eight_Queens()
Dim i As Long, sText As String
' 禁用事件加快运行的速度
DisEvents
For i = 0 To iFrontier
arrQueen(i) = -1
Next
iCcount = 0: yPosition = 1: xPosition = 1
placeQueen 0
Erase arrQueen
sText = "总共 " & iCcount & " 种解法!"
yPosition = 1: iCcount = 0
Cells(xPosition, yPosition).Value = sText
Cells(xPosition, yPosition).Font.Bold = True
Range("b2:i" & xPosition).Font.Italic = True
yPosition = 0: xPosition = 0
Cells.Columns.AutoFit
EnEvents
MsgBox sText, vbInformation + vbOKOnly, "Tips"
End Sub
Private Sub placeQueen(ByVal n As Long)
' 算法核心'
Dim i As Long, k As Long, j As Long, m As Long
Dim arrStatus(iBoundary) As Boolean
For i = 0 To iFrontier
arrStatus(i) = True
Next
i = 0
Do While i < n
arrStatus(arrQueen(i)) = False
k = n - i
m = arrQueen(i) + k
If m >= 0 And m < iBoundary Then arrStatus(m) = False
m = arrQueen(i) - k
If m >= 0 And m < iBoundary Then arrStatus(m) = False
i = i + 1
Loop
For i = 0 To iFrontier
If arrStatus(i) = True Then
If n < iFrontier Then
arrQueen(n) = i
placeQueen (n + 1)
Else
arrQueen(n) = i
iCcount = iCcount + 1
yPosition = 1
Cells(xPosition, yPosition).Value = "这是第" & iCcount & "个解法, 如下:"
Cells(xPosition, yPosition).Font.Bold = True
xPosition = xPosition + 1
Dim a As Long
a = xPosition
For m = 0 To iFrontier
If arrQueen(m) = 0 Then
yPosition = 2
Else
yPosition = arrQueen(m) + 2
End If
Cells(xPosition, yPosition).Value = "Queen"
Cells(xPosition, yPosition).Interior.Color = 65535
xPosition = xPosition + 1
Next
CreateBorders a, xPosition - 1
i = m
End If
End If
Next
End Sub
Private Sub CreateBorders(ByVal x As Long, ByVal y As Long)
' 绘制图形
Dim iRng As Range
Set iRng = Range("b" & x & ":" & "i" & y)
iRng.Borders(xlDiagonalDown).LineStyle = xlNone
iRng.Borders(xlDiagonalUp).LineStyle = xlNone
With iRng.Borders(xlEdgeLeft)
.LineStyle = xlContinuous
.ColorIndex = 0
.TintAndShade = 0
.Weight = xlThin
End With
With iRng.Borders(xlEdgeTop)
.LineStyle = xlContinuous
.ColorIndex = 0
.TintAndShade = 0
.Weight = xlThin
End With
With iRng.Borders(xlEdgeBottom)
.LineStyle = xlContinuous
.ColorIndex = 0
.TintAndShade = 0
.Weight = xlThin
End With
With iRng.Borders(xlEdgeRight)
.LineStyle = xlContinuous
.ColorIndex = 0
.TintAndShade = 0
.Weight = xlThin
End With
With iRng.Borders(xlInsideVertical)
.LineStyle = xlContinuous
.ColorIndex = 0
.TintAndShade = 0
.Weight = xlThin
End With
With iRng.Borders(xlInsideHorizontal)
.LineStyle = xlContinuous
.ColorIndex = 0
.TintAndShade = 0
.Weight = xlThin
End With
Set iRng = Nothing
End Sub
Private Sub DisEvents() '禁用干扰项
With ThisWorkbook.Application
.ScreenUpdating = False '禁止屏幕刷新
.EnableEvents = False '禁用事件
.Calculation = xlCalculationManual '禁用自动计算
.Interactive = False '禁止交互(在执行宏时,如果在表格输入内容会造成宏终止)
End With
End Sub
Private Sub EnEvents() '启用
With ThisWorkbook.Application
.ScreenUpdating = True
.EnableEvents = True
.Calculation = xlCalculationAutomatic
.Interactive = True
End With
End Sub
The diagonal constraint is a little trickier than the row and column constraints. First, if two queens lie on the same diagonal, one of the following conditions must be true:
The row number plus the column number for each of the two queens are equal. In other words, queens(j) + j has the same value for two different indices j.
The row number minus the column number for each of the two queens are equal. In this case, queens(j) - j has the same value for two different indices j.
One of these conditions means the queens lie on the same ascending diagonal ( going from left to right), while the other means they lie on the same descending diagonal. Which condition corresponds to ascending and which to descending depends on how you order the rows and columns. As mentioned in the previous section, the ordering has no effect on the set of solutions, just on how you visualize them.
So the diagonal constraint is that the values of queens(j) + j must all be different, and the values of queens(j) - j must all be different, for different j.
To apply the AddAllDifferent method to queens(j) + j, we put the N instances of the variable, for j from 0 to N-1, into an array, diag1, as follows:
The constraint for queens(j) - j is created similarly.
# Creates the variables.
# There are `board_size` number of variables, one for a queen in each column
# of the board. The value of each variable is the row that the queen is in.
queens = [model.NewIntVar(0, board_size - 1, f"x_{i}") for i in range(board_size)]
# 这里只创建了八个变量, 因为和数独不相同的地方在于, 八皇后不需要只需要建立对8个放置皇后的位置坐标变量进行约束
# Creates the constraints.
# All rows must be different.
model.AddAllDifferent(queens)
# 让后让八个变量各不相同, 这部分还勉强可以理解
# 为什么这里可以如此表示对角线的约束?
# No two queens can be on the same diagonal.
model.AddAllDifferent(queens[i] + i for i in range(board_size))
model.AddAllDifferent(queens[i] - i for i in range(board_size))
# 但是关于对角线的这部分, 不是很理解其中的逻辑
其解题逻辑和数独当中的求解很大的区别.
"""OR-Tools solution to the N-queens problem."""
import sys
import time
from ortools.sat.python import cp_model
class NQueenSolutionPrinter(cp_model.CpSolverSolutionCallback):
"""Print intermediate solutions."""
def __init__(self, queens):
cp_model.CpSolverSolutionCallback.__init__(self)
self.__queens = queens
self.__solution_count = 0
self.__start_time = time.time()
def solution_count(self):
return self.__solution_count
def on_solution_callback(self):
current_time = time.time()
print(
f"Solution {self.__solution_count}, "
f"time = {current_time - self.__start_time} s"
)
self.__solution_count += 1
all_queens = range(len(self.__queens))
for i in all_queens:
for j in all_queens:
if self.Value(self.__queens[j]) == i:
# There is a queen in column j, row i.
print("Q", end=" ")
else:
print("_", end=" ")
print()
print()
def main(board_size):
# Creates the solver.
model = cp_model.CpModel()
# Creates the variables.
# There are `board_size` number of variables, one for a queen in each column
# of the board. The value of each variable is the row that the queen is in.
queens = [model.NewIntVar(0, board_size - 1, f"x_{i}") for i in range(board_size)]
# Creates the constraints.
# All rows must be different.
model.AddAllDifferent(queens)
# No two queens can be on the same diagonal.
model.AddAllDifferent(queens[i] + i for i in range(board_size))
model.AddAllDifferent(queens[i] - i for i in range(board_size))
# Solve the model.
solver = cp_model.CpSolver()
solution_printer = NQueenSolutionPrinter(queens)
solver.parameters.enumerate_all_solutions = True
solver.Solve(model, solution_printer)
# Statistics.
print("\nStatistics")
print(f" conflicts : {solver.NumConflicts()}")
print(f" branches : {solver.NumBranches()}")
print(f" wall time : {solver.WallTime()} s")
print(f" solutions found: {solution_printer.solution_count()}")
if __name__ == "__main__":
# By default, solve the 8x8 problem.
size = 8
if len(sys.argv) > 1:
size = int(sys.argv[1])
main(size)
queens = [model.NewIntVar(0, board_size - 1, f"x_{i}") for i in range(board_size)]
# 列约束, 每一列都不相同
model.AddAllDifferent(queens)
# 这时候再看对角线约束就很容易理解
model.AddAllDifferent(queens[i] + i for i in range(board_size))
model.AddAllDifferent(queens[i] - i for i in range(board_size))
all_queens = range(len(self.__queens))
for i in all_queens:
for j in all_queens:
if self.Value(self.__queens[j]) == i:
# There is a queen in column j, row i.
print("Q", end=" ")
else:
print("_", end=" ")
print()
print()
from ortools.linear_solver import pywraplp
N = 8
solver = pywraplp.Solver('queens', pywraplp.Solver.SCIP_MIXED_INTEGER_PROGRAMMING)
variables = {}
for r in range(N):
for c in range(N):
key = (r, c)
variables[key] = solver.BoolVar('_'.join(str(e) for e in key))
for i in range(N):
# 行约束
r_constrain = sum(variables[i, c] for c in range(N)) == 1
solver.Add(r_constrain)
# 列约束
c_constrain = sum(variables[r, i] for r in range(N)) == 1
solver.Add(c_constrain)
# \ 向 对角线约束
for k in range(2 - N, N - 1):
c = sum(variables[i, i - k] for i in range(N) if 0 <= i - k < N) <= 1
solver.Add(c)
# / 向对角线约束
for k in range(1, 2 * N - 2):
c = sum(variables[i, k - i] for i in range(N) if 0 <= k - i < N) <= 1
solver.Add(c)
status = solver.Solve()
# 0
# ic = 0
# 无效特性
# while solver.NextSolution():
# ic += 1
# print(f"Solution #{ic}:")
# for (r, c) in variables.keys():
# if variables[r, c].solution_value() == 1:
# print(f" queen {r} {c}")
# print(ic)
# 由于存在多个可行解, 这里只能得到一个解, 每次的求解得到的可能是不同的解
for r in range(N):
for c in range(N):
if variables[r, c].solution_value() == 1:
print("Q", end=" ")
else:
print("_", end=" ")
print()
print()
As of 2020-02-10, only Gurobi and SCIP support NextSolution()
def NextSolution(self) -> bool
Some solvers (MIP only, not LP) can produce multiple solutions to the problem. Returns true when another solution is available, and updates the MPVariable* objects to make the new solution queryable. Call only after calling solve.
The optimality properties of the additional solutions found, and whether or not the solver computes them ahead of time or when NextSolution() is called is solver specific.
As of 2020-02-10, only Gurobi and SCIP support NextSolution(), see linear_solver_interfaces_test for an example of how to configure these solvers for multiple solutions. Other solvers return false unconditionally.
但实际上, scip求解器并不支持这个特性.
在Google ortools groups, 相关的问题.
六. 背包问题
The goal of packing problems is to find the best way to pack a set of items of given sizes into containers with fixed capacities. A typical application is loading boxes onto delivery trucks efficiently. Often, it's not possible to pack all the items, due to the capacity constraints. In that case, the problem is to find a subset of the items with maximum total size that will fit in the containers.
from ortools.sat.python import cp_model
workers = 4 # 工人
shifts = 3 # 每天排班次数
days = 3 # 持续天数
worker_ids = [_ for _ in range(workers)]
shift_ids = [_ for _ in range(shifts)]
day_ids = [_ for _ in range(days)]
model = cp_model.CpModel()
configs = {}
for w in worker_ids:
for d in day_ids:
for s in shift_ids:
configs[(w, d, s)] = model.NewBoolVar(f"shift_n{w}_d{d}_s{s}")
# 每天 - 每个班次 只安排一个工人
for d in day_ids:
for s in shift_ids:
model.AddExactlyOne(configs[(w, d, s)] for w in worker_ids)
# 每个工人, 每天, 每个班次只安排最多一次班
for w in worker_ids:
for d in day_ids:
model.AddAtMostOne(configs[(w, d, s)] for s in shift_ids)
# 均衡分配任务
# 每名工人的最少排班次数
min_shifts_per_worker = (shifts * days) // workers
# 每名工人的最大排班次数
# 假如刚刚好分配, 最大值就是最小值; 假如不能平均分配, 那么就需要有部分人分配多1次, 即平均值(取整数部分) + 1
max_shifts_per_worker = min_shifts_per_worker if shifts * days % workers == 0 else min_shifts_per_worker + 1
# 每个工人
for w in worker_ids:
shifts_worked = []
for d in day_ids:
# 每天 / 每班次
for s in shift_ids:
shifts_worked.append(configs[(w, d, s)])
# 每个工人这么多天当中的累积排班次数上限 & 下限的限制
constraint = sum(shifts_worked)
model.Add(min_shifts_per_worker <= constraint)
model.Add(constraint <= max_shifts_per_worker)
solver = cp_model.CpSolver()
# quick_restart_no_lp (kind of probing with linearization_level:0)
# 禁用线性规划求解器
# 需要注意单线程/多线程
# https://github.com/google/or-tools/discussions/2666
# 参数的作用
# solver.parameters.linearization_level = 0
# solver.SetSolverSpecificParametersAsString("linearization_level:0")
solver.parameters.enumerate_all_solutions = True
class WorkerShiftSolutionPrinter(cp_model.CpSolverSolutionCallback):
def __init__(self, configs, workers, days, shifts, limit):
super().__init__()
self._configs = configs
self._workers = workers
self._days = days
self._shifts = shifts
self._solution_limit = limit
self._solution_count = 0
def on_solution_callback(self):
self._solution_count += 1
print(f"solution {self._solution_count}")
for d in range(self._days):
print(f"day: {d}")
for w in range(self._workers):
is_working = False
for s in range(self._shifts):
if self.Value(self._configs[(w, d, s)]):
is_working = True
break
print(f" worker {w} works shift {s}" if is_working else f" worker {w} does not work")
if self._solution_count >= self._solution_limit:
print(f"stop search after {self._solution_limit} solutions")
self.StopSearch()
def solution_count(self):
return self._solution_count
solution_limit = 5
solution_printer = WorkerShiftSolutionPrinter(configs, workers, days, shifts, solution_limit)
solver.Solve(model, solution_printer)
print(f" - conflicts : {solver.NumConflicts()}") # Returns the number of conflicts since the creation of the solver.
print(f" - branches : {solver.NumBranches()}") # Returns the number of search branches explored by the solver.
print(f" - wall time : {solver.WallTime()} s")
solution 1
day: 0
worker 0 works shift 1
worker 1 does not work
worker 2 works shift 2
worker 3 works shift 0
day: 1
worker 0 works shift 2
worker 1 works shift 0
worker 2 does not work
worker 3 works shift 1
day: 2
worker 0 does not work
worker 1 works shift 2
worker 2 works shift 0
worker 3 works shift 1
...
model.NewIntervalVar?
Signature:
model.NewIntervalVar(
start: Union[ForwardRef('LinearExpr'), ForwardRef('IntVar'), numbers.Integral, numpy.integer, int],
size: Union[ForwardRef('LinearExpr'), ForwardRef('IntVar'), numbers.Integral, numpy.integer, int],
end: Union[ForwardRef('LinearExpr'), ForwardRef('IntVar'), numbers.Integral, numpy.integer, int],
name: str,
) -> ortools.sat.python.cp_model.IntervalVar
Docstring:
Creates an interval variable from start, size, and end.
An interval variable is a constraint, that is itself used in other
constraints like NoOverlap.
Internally, it ensures that `start + size == end`.
Args:
start: The start of the interval. It can be an affine or constant
expression.
size: The size of the interval. It can be an affine or constant
expression.
end: The end of the interval. It can be an affine or constant expression.
name: The name of the interval variable.
Returns:
An `IntervalVar` object.
这里引入一个专门的变量类型 - 间隔变量.
# 创建非重叠约束
for machine in all_machines:
model.AddNoOverlap(machine_to_intervals[machine])
# 确定执行的优先顺序约束
for job_id, job in enumerate(jobs_data):
for task_id in range(len(job) - 1):
model.Add(all_tasks[job_id, task_id + 1].start >= all_tasks[job_id, task_id].end)
model.AddNoOverlap?
Signature:
model.AddNoOverlap(
interval_vars: Iterable[ortools.sat.python.cp_model.IntervalVar],
) -> ortools.sat.python.cp_model.Constraint
Docstring:
Adds NoOverlap(interval_vars).
A NoOverlap constraint ensures that all present intervals do not overlap
in time.
Args:
interval_vars: The list of interval variables to constrain.
Returns:
An instance of the `Constraint` class.
obj_var = model.NewIntVar(0, horizon, "makespan")
model.AddMaxEquality(obj_var, [all_tasks[job_id, len(job) - 1].end for job_id, job in enumerate(jobs_data)])
model.Minimize(obj_var)
# Creates the solver and solve.
solver = cp_model.CpSolver()
status = solver.Solve(model)
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print("Solution:")
# Create one list of assigned tasks per machine.
assigned_jobs = collections.defaultdict(list)
for job_id, job in enumerate(jobs_data):
for task_id, task in enumerate(job):
machine = task[0]
assigned_jobs[machine].append(
assigned_task_type(
start=solver.Value(all_tasks[job_id, task_id].start),
job=job_id,
index=task_id,
duration=task[1],
)
)
# Create per machine output lines.
output = ""
for machine in all_machines:
# Sort by starting time.
assigned_jobs[machine].sort()
sol_line_tasks = "Machine " + str(machine) + ": "
sol_line = " "
for assigned_task in assigned_jobs[machine]:
name = f"job_{assigned_task.job}_task_{assigned_task.index}"
# Add spaces to output to align columns.
sol_line_tasks += f"{name:15}"
start = assigned_task.start
duration = assigned_task.duration
sol_tmp = f"[{start},{start + duration}]"
# Add spaces to output to align columns.
sol_line += f"{sol_tmp:15}"
sol_line += "\n"
sol_line_tasks += "\n"
output += sol_line_tasks
output += sol_line
# Finally print the solution found.
print(f"Optimal Schedule Length: {solver.ObjectiveValue()}")
print(output)
else:
print("No solution found.")
# Statistics.
print("\nStatistics")
print(f" - conflicts: {solver.NumConflicts()}")
print(f" - branches : {solver.NumBranches()}")
print(f" - wall time: {solver.WallTime()}s")
八. 分配问题
One of the most well-known combinatorial optimization problems is the assignment problem. Here's an example: suppose a group of workers needs to perform a set of tasks, and for each worker and task, there is a cost for assigning the worker to the task. The problem is to assign each worker to at most one task, with no two workers performing the same task, while minimizing the total cost.
from ortools.linear_solver import pywraplp
costs = [
[90, 80, 75, 70],
[35, 85, 55, 65],
[125, 95, 90, 95],
[45, 110, 95, 115],
[50, 100, 90, 100],
]
num_workers = 5
num_tasks = 4
solver = pywraplp.Solver.CreateSolver("SCIP")
configs = {}
for i in range(num_workers):
for j in range(num_tasks):
configs[i, j] = solver.BoolVar(f'{i}_{j}')
# 每个工作器最大分配一个任务
for i in range(num_workers):
solver.Add(solver.Sum([configs[i, j] for j in range(num_tasks)]) <= 1)
# 每个任务只能分配给一个工作器
for j in range(num_tasks):
solver.Add(solver.Sum([configs[i, j] for i in range(num_workers)]) == 1)
# 目标函数
cost_terms = []
for i in range(num_workers):
for j in range(num_tasks):
cost_terms.append(costs[i][j] * configs[i, j])
solver.Minimize(solver.Sum(cost_terms))
status = solver.Solve()
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
print(f"total cost = {solver.Objective().Value()}\n")
for i in range(num_workers):
for j in range(num_tasks):
if int(configs[i, j].solution_value()) == 1:
print(f"worker {i} assigned to task {j}." + f" cost: {costs[i][j]}")
else:
print("no solution found.")
total cost = 265.0
worker 0 assigned to task 3. cost: 70
worker 1 assigned to task 2. cost: 55
worker 2 assigned to task 1. cost: 95
worker 3 assigned to task 0. cost: 45
"""Solves an assignment problem for given group of workers."""
from ortools.sat.python import cp_model
def main() -> None:
# Data
costs = [
[90, 76, 75, 70, 50, 74],
[35, 85, 55, 65, 48, 101],
[125, 95, 90, 105, 59, 120],
[45, 110, 95, 115, 104, 83],
[60, 105, 80, 75, 59, 62],
[45, 65, 110, 95, 47, 31],
[38, 51, 107, 41, 69, 99],
[47, 85, 57, 71, 92, 77],
[39, 63, 97, 49, 118, 56],
[47, 101, 71, 60, 88, 109],
[17, 39, 103, 64, 61, 92],
[101, 45, 83, 59, 92, 27],
]
num_workers = len(costs)
num_tasks = len(costs[0])
# Allowed groups of workers:
group1 = [
[0, 0, 1, 1], # Workers 2, 3
[0, 1, 0, 1], # Workers 1, 3
[0, 1, 1, 0], # Workers 1, 2
[1, 1, 0, 0], # Workers 0, 1
[1, 0, 1, 0], # Workers 0, 2
]
group2 = [
[0, 0, 1, 1], # Workers 6, 7
[0, 1, 0, 1], # Workers 5, 7
[0, 1, 1, 0], # Workers 5, 6
[1, 1, 0, 0], # Workers 4, 5
[1, 0, 0, 1], # Workers 4, 7
]
group3 = [
[0, 0, 1, 1], # Workers 10, 11
[0, 1, 0, 1], # Workers 9, 11
[0, 1, 1, 0], # Workers 9, 10
[1, 0, 1, 0], # Workers 8, 10
[1, 0, 0, 1], # Workers 8, 11
]
# Model
model = cp_model.CpModel()
# Variables
x = {}
for worker in range(num_workers):
for task in range(num_tasks):
x[worker, task] = model.new_bool_var(f"x[{worker},{task}]")
# Constraints
# Each worker is assigned to at most one task.
for worker in range(num_workers):
model.add_at_most_one(x[worker, task] for task in range(num_tasks))
# Each task is assigned to exactly one worker.
for task in range(num_tasks):
model.add_exactly_one(x[worker, task] for worker in range(num_workers))
# Create variables for each worker, indicating whether they work on some task.
work = {}
for worker in range(num_workers):
work[worker] = model.new_bool_var(f"work[{worker}]")
for worker in range(num_workers):
for task in range(num_tasks):
model.add(work[worker] == sum(x[worker, task] for task in range(num_tasks)))
# Define the allowed groups of worders
model.add_allowed_assignments([work[0], work[1], work[2], work[3]], group1)
model.add_allowed_assignments([work[4], work[5], work[6], work[7]], group2)
model.add_allowed_assignments([work[8], work[9], work[10], work[11]], group3)
# Objective
objective_terms = []
for worker in range(num_workers):
for task in range(num_tasks):
objective_terms.append(costs[worker][task] * x[worker, task])
model.minimize(sum(objective_terms))
# Solve
solver = cp_model.CpSolver()
status = solver.solve(model)
# Print solution.
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
print(f"Total cost = {solver.objective_value}\n")
for worker in range(num_workers):
for task in range(num_tasks):
if solver.boolean_value(x[worker, task]):
print(
f"Worker {worker} assigned to task {task}."
+ f" Cost = {costs[worker][task]}"
)
else:
print("No solution found.")
if __name__ == "__main__":
main()
Total cost = 239.0
Worker 0 assigned to task 4. Cost = 50
Worker 1 assigned to task 2. Cost = 55
Worker 5 assigned to task 5. Cost = 31
Worker 6 assigned to task 3. Cost = 41
Worker 10 assigned to task 0. Cost = 17
Worker 11 assigned to task 1. Cost = 45
8.3 线性求和求解器
linear sum assignment solver, a specialized solver for the simple assignment problem, which can be faster than either the MIP or CP-SAT solver. However, the MIP and CP-SAT solvers can handle a much wider array of problems, so in most cases they are the best option.
Note: The linear sum assignment solver only accepts integer values for the weights and values. The section Using a solver with non-integer data shows how to use the solver if your data contains non-integer values.
np.meshgrid?
Return a list of coordinate matrices from coordinate vectors.
Make N-D coordinate arrays for vectorized evaluations of
N-D scalar/vector fields over N-D grids, given
one-dimensional coordinate arrays x1, x2,..., xn.
.. versionchanged:: 1.9
1-D and 0-D cases are allowed.
The array is the cost matrix, whose i, j entry is the cost for worker i to perform task j. There are four workers, corresponding to the rows of the matrix, and four tasks, corresponding to the columns.
the best path, 5634:
44->43->3->42->34->15->37->36->18->5->45->39->1->26->2->46->33->31->8->6->7->32->11->48->27->49->12->0->28->35->13->41->24->29->19->38->10->21->47->23->40->14->30->25->4->20->9->17->22->16
LKH is an effective implementation of the Lin-Kernighan heuristic for solving the traveling salesman problem.
Computational experiments have shown that LKH is highly effective. Even though the algorithm is approximate, optimal solutions are produced with an impressively high frequency. LKH has produced optimal solutions for all solved problems we have been able to obtain; including a 109399-city instance (at the time of writing, the largest nontrivial instance solved to optimality). Furthermore, the algorithm has improved the best known solutions for a series of large-scale instances with unknown optima, among these a 1,904,711-city instance (World TSP).
loc = np.array([[indexs[i], indexs[i + 1]] for i in range(50)])
distance_matrix[loc[:, 0], loc[:, 1]].sum()
# 5634
'->'.join(str(index) for index in indexs)
Many problems in computer science can be represented by a graph consisting of nodes and links between them. Examples are network flow problems, which involve transporting goods or material across a network, such as a railway system.
You can represent a network flow by a graph whose nodes are cities and whose arcs are rail lines between them. (They're called flows because their properties are similar to those of water flowing through a network of pipes.)
A key constraint in network flows is that each arc has a capacity - the maximum amount that can be transported across the arc in a fixed period of time.
The maximum flow problem is to determine the maximum total amount that can be transported across all arcs in the network, subject to the capacity constraints.
The first person to study this problem was the Russian mathematician A.N. Tolstoi, in the 1930s. The map below shows the actual railway network for which he wanted to find a maximum flow.
Return type: boolCompletes initialization after the problem is fully specified. Returns true if we successfully prove that arithmetic calculations are guaranteed not to overflow. ComputeAssignment() calls this method itself, so only clients that care about obtaining a warning about the possibility of arithmetic precision problems need to call this method explicitly. Separate from ComputeAssignment() for white-box testing and for clients that need to react to the possibility that arithmetic overflow is not ruled out. FinalizeSetup() is idempotent.
Return type: CostValueReturns the cost of the minimum-cost perfect matching. Precondition: success_ == true, signifying that we computed the optimum assignment for a feasible problem.
Return type: inline NodeIndexArguments: ArcIndex arcThese handy member functions make the code more compact, and we expose them to clients so that client code that doesn't have direct access to the graph can learn about the optimum assignment once it is computed.
Arguments: const GraphType& graph, NodeIndex num_left_nodesConstructor for the case in which we will build the graph incrementally as we discover arc costs, as might be done with any of the dynamic graph representations such as StarGraph or ForwardStarGraph.
Arguments: NodeIndex num_left_nodes, ArcIndex num_arcsConstructor for the case in which the underlying graph cannot be built until after all the arc costs are known, as is the case with ForwardStarStaticGraph. In this case, the graph is passed to us later via the SetGraph() method, below.
Return type: voidArguments: GraphType* graphOptimizes the layout of the graph for the access pattern our implementation will use. REQUIRES for LinearSumAssignment template instantiation if a call to the OptimizeGraphLayout() method is compiled: GraphType is a dynamic graph, i.e., one that implements the GroupForwardArcsByFunctor() member template method. If analogous optimization is needed for LinearSumAssignment instances based on static graphs, the graph layout should be constructed such that each node's outgoing arcs are sorted by head node index before the LinearSumAssignment::SetGraph() method is called.
Return type: voidArguments: CostValue factorSets the cost-scaling divisor, i.e., the amount by which we divide the scaling parameter on each iteration.
Return type: voidArguments: const GraphType* graphSets the graph used by the LinearSumAssignment instance, for use when the graph layout can be determined only after arc costs are set. This happens, for example, when we use a ForwardStarStaticGraph.
Return type: ArcIndexArguments: NodeIndex tail, NodeIndex head, FlowQuantity capacityAdds a directed arc with the given capacity from tail to head. * Node indices and capacity must be non-negative (>= 0). * Self-looping and duplicate arcs are supported. * After the method finishes, NumArcs() == the returned ArcIndex + 1.
Return type: FlowQuantityArguments: ArcIndex arcReturns the flow on the given arc in the last OPTIMAL Solve() context. Note: It is possible that there is more than one optimal solution. The algorithm is deterministic so it will always return the same solution for a given problem. However, there is no guarantee of this from one code version to the next (but the code does not change often).
Return type: voidArguments: std::vector<NodeIndex>* resultReturns the nodes that can reach the sink by non-saturated arcs, the outgoing arcs of this set form a minimum cut. Note that if this is the complement set of GetNodeReachableFromSource(), then the min-cut is unique. This works only if Solve() returned OPTIMAL.
Return type: voidArguments: std::vector<NodeIndex>* resultReturns the nodes reachable from the source by non-saturated arcs (.i.e. arc with Flow(arc) < Capacity(arc)), the outgoing arcs of this set form a minimum cut. This works only if Solve() returned OPTIMAL.
Return type: voidArguments: ArcIndex arc, FlowQuantity capacityChange the capacity of an arc. WARNING: This looks like it enables incremental solves, but as of 2018-02, the next Solve() will restart from scratch anyway. TODO(user): Support incrementality in the max flow implementation.
Return type: NodeIndexArguments: ArcIndex arcReturns user-provided data. The implementation will crash if "arc" is not in [0, NumArcs()).
10.1.3.1
示例
import numpy as np
from ortools.graph.python import max_flow
smf = max_flow.SimpleMaxFlow()
# Define three parallel arrays: start_nodes, end_nodes, and the capacities
# between each pair. For instance, the arc from node 0 to node 1 has a
# capacity of 20.
start_nodes = np.array([0, 0, 0, 1, 1, 2, 2, 3, 3])
end_nodes = np.array([1, 2, 3, 2, 4, 3, 4, 2, 4])
capacities = np.array([20, 30, 10, 40, 30, 10, 20, 5, 20])
# Add arcs in bulk.
# note: we could have used add_arc_with_capacity(start, end, capacity)
all_arcs = smf.add_arcs_with_capacity(start_nodes, end_nodes, capacities)
# Find the maximum flow between node 0 and node 4.
status = smf.solve(0, 4)
if status != smf.OPTIMAL:
print("There was an issue with the max flow input.")
print(f"Status: {status}")
exit(1)
print("Max flow:", smf.optimal_flow())
print("")
print(" Arc Flow / Capacity")
solution_flows = smf.flows(all_arcs)
for arc, flow, capacity in zip(all_arcs, solution_flows, capacities):
print(f"{smf.tail(arc)} / {smf.head(arc)} {flow:3} / {capacity:3}")
print("Source side min-cut:", smf.get_source_side_min_cut())
print("Sink side min-cut:", smf.get_sink_side_min_cut())
.png)
