Gurobi的Python接口学习笔记和TSP问题求解案例

SACUKI

墨衍会员 · AI 创作全网分发码龄5年

关注

原创 2025-01-20 22:42:31 1.2k 阅读 AI 写同款· GEO 优化 ›

Gurobi求解器基础入门官方教程 Gurobi是一种数学规划（线性和凸二次规划）优化器。支持多种语言接口，本文以python+gurobi为主。阅读详情

文章目录

Basic Gurobi Concepts
Python API Fundamentals
Model Generator Best Practices
Iterating in Python and Aggregate Sums[^1]
TSP Example[^2]

Basic Gurobi Concepts

在这里插入图片描述

Parameters：参数用于控制求解器的行为，相关参数参考： https://docs.gurobi.com/projects/optimizer/en/current/concepts/parameters/groups.html
Attributes：属性是查询和修改相关对象特征的机制，相关属性参考：https://docs.gurobi.com/projects/optimizer/en/current/concepts/attributes/types.html
Environment：环境是Model的容器，gurobi的python接口会自动创建环境。设置环境的Parameters会让环境中的每个Model都默认这些Parameters，也可以单独为每个模型设置Parameters。

Python API Fundamentals

Important global functions
- setParam()：设置环境的参数
- read()：从一个model文件中读取，创建
Primary objects
- Model
- Var
- Constr
Other important methods
- Model.setObjective(expression)：传入Expr对象，创建目标函数
- Model.optimize()：开启优化
- Model.reset()：重置模型的状态，discarding any previously computed solution information.
Model summary information methods
- Model.printAtrr()：Print the value of one or more attributes.
- Model.printStatus()：Print statistics about the model (number of constraints and variables, number of non-zeros in constraint matrix, smallest and largest coefficients, etc.).
- Model.printQuality()：Print statistics about the quality of the computed solution (constraint violations, integrality violations, etc.).
- Model.getVars()：Retrieve a list of all variables in the model.
- Model.getConstrs()：Retrieve a list of all linear constraints in the model.
- Model.getRow()：Retrieve the list of variables that participate in a constraint, and the associated coefficients. The result is returned as a LinExprobject.
- Model.getCol()：Retrieve the list of constraints in which a variable participates, and the associated coefficients. The result is returned as a Column object.

Model Generator Best Practices

在这里插入图片描述

使用Gurobi求解优化模型比较有条理的步骤如下：

先建模，建模可按以下步骤一步一步分析，避免脑袋绕晕
1. 定义模型中要使用的集合
2. 定义模型中要使用的参数
3. 定义模型中要使用的决策变量
4. 定义模型的目标函数
5. 分析处模型设计的所有约束条件
然后使用Gurobi求解，具体又可按以下步骤进行
1. 使用list或tuplelist对象保存集合（附：集合在模型中一般用作索引）
2. 使用dict或tupledict对象保存参数
3. 创建模型对象
  1. 添加变量
  2. 设置目标函数
  3. 添加约束条件
  4. 进行优化

tuplelist和tupledict相比list和dict的关系是什么？前者时专门为gurobi设计的，索引和运算能力更强大。

Iterating in Python and Aggregate Sums¹

当遇到for-all这种建模结构时，使用Python的for循环语句很有帮助，考虑下面这样一个约束条件:
$\sum_{j\in J}x_{ij}\le5,\quad \forall i\in I$

for i in I:
	model.addConstr(gp.quicksum(x[i,j] for j in J)<=5)

先版本支持一种使用生成器的更简洁写法，只不过没有第一种好理解，

model.addConstrs((x.sum(i,'*')<=5 for i in I))

全局函数quicksum和tupledict对象的sum/prod方法的关系：（1）quicksum与sum方法相比更灵活，因为它可以接受任何可迭代的表达式，可以处理不同类型的变量和系数组合，甚至是复杂的表达式，只要它们最终可以相加。tupledict的 sum方法则局限于对 tupledict对象中的变量进行求和。（2）对于点积用tupledict对象的prod方法和全局函数quicksum进行计算没啥区别。

TSP Example²

TSP 问题最早的描述是1759年欧拉研究的骑士周游问题, 即对于国际象棋棋盘中的64 个方格, 走访64个方格一次且仅一次, 并最终返回起点。其一般性描述为: 有一旅行商要访问 𝑛个城市, 每个城市必须访问且只能访问一次, 需要寻求到一条包含所有 𝑛个城市的最短访问路线。TSP问题的模型可以写为：
$\begin{align} &min\quad z=\sum_{(i,j)\in A,i\ne j}d_{ij}x_{ij}\\s.t. &\left\{\begin{matrix} \sum_{j\ne i}x_{ij}&=1,&\forall i \in V\\ \sum_{i\ne j}x_{ij}&=1,&\forall j \in V\\ u_{i}-u_{j}+|V|x_{ij}&\le|V|-1,&\forall (i,j)\in A,i\ne j,i\ne o,j\ne 0\\ x_{ij}&=0或1,&(i,j)\in A\\ u_{i}&\ge0,&i\in V \end{matrix}\right. \end{align}$
建模说明：

集合
- V：图中节点集合
- A：途中所有弧的集合
参数
- $d_{ij}$ ：节点i和j之间的距离
决策变量
- $x_{ij}$ 是01变量，如果为1则表示从i走到j，否则不从i走到j
- $u_{i}$ 是为了使用MTZ约束引入的决策变量，可以理解为节点的访问次序？³
目标函数：最小化总的旅行距离
约束条件：从起点经过每一个点最后回到起点，即哈密顿回路约束——前两个约束表示图中任何一个点都离开1次到达1次，第三个约束叫做MTZ约束，用于消除子回路。此外，消除子回路还可以使用subtour-elimination,这个约束原理比较简单，通过枚举消除所有子回路，只留下最大的那个回路，但是在在gurobi里面得用回调函数写²，不方便。

求解说明，仍然以TSPLIB提供的brilin52算例进行实验，算例获取方式参考。

首先，导入需要使用的包，

import math
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import gurobipy as gp
from gurobipy import GRB

再把berlin52的数据读进来，

# 读取berlin52算例
def read_locations(loc):
    matrix = []
    with open(loc, 'r') as file:
        lines = file.readlines()
        start_reading = False
        for line in lines:
            if line.startswith('NODE_COORD_SECTION'):
                start_reading = True
                continue
            if line.startswith('EOF'):
                break
            if start_reading:
                parts = line.strip().split()
                if len(parts) == 3:
                    _, x, y = parts
                    matrix.append([float(x), float(y)])
    return np.array(matrix)

# 计算任意两地的距离
def calc_distances(city_locations):
    return np.array(([[math.sqrt(np.sum((loc1-loc2)**2)) for loc2 in city_locations]
             for loc1 in city_locations]))

city_locations = read_locations(loc=r"D:/待处理的文件/CodeSpace_for_spyder/遗传算法案例/berlin52/berlin52.tsp") # 读取到berlin52处地点的欧几里得坐标
distances = calc_distances(city_locations) # 计算出任意两地之间的距离矩阵

然后，按照我们上文提到过的Gurobi求解步骤，一步一步的写求解代码，避免脑袋绕晕，

# 定义模型中要用到的集合
V = [i for i in range(distances.shape[0])]
A = [(i,j) for i in range(distances.shape[0]) for j in range(distances.shape[0])]
# 定义模型中要用到的参数
d = gp.tupledict({(i,j):distances[i,j] for i in range(distances.shape[0]) for j in range(distances.shape[0])})
# 创建模型
model = gp.Model(name='TSP')
# 添加决策变量
x = model.addVars(A,vtype=GRB.BINARY,name='x') # x_{ij}
u = model.addVars(V,vtype=GRB.CONTINUOUS,name='u') # u_{i}
# 设置目标函数
model.setObjective(x.prod(d),sense=GRB.MINIMIZE)
# 添加约束条件
## 约束（1）
for i in V:
    model.addConstr(
        gp.quicksum(x[i,j] for j in V if j!=i) == 1
        )
## 约束（2）
for j in V:
    model.addConstr(
        gp.quicksum(x[i,j] for i in V if i!=j) == 1
        )
## 约束(3)
V_abs = len(V)
for i in V[1:]:
    for j in V[1:]:
        if i != j :
            model.addConstr(
                u[i]-u[j]+V_abs*x[i,j] <= V_abs-1
                )
# 进行优化
model.optimize()

最后，分析并可视化以下求解结果，

# 结果分析与可视化
if model.Status == GRB.OPTIMAL:
    print('最小巡游距离:{:.2f}'.format(model.ObjVal))
best_x = np.zeros((V_abs,V_abs)) # 用于存放最优解
for i in range(V_abs):
    for j in range(V_abs):
        best_x[i,j] = int(x[i,j].X)
pos1,pos2 = np.where(best_x) # 找到所有非0值对应索引,pos1是行索引，pos2是列索引
order = [] # 将最优解转换为访问次序，存储在order中
previous_point = pos1[0]
order.append(previous_point)
for _ in range(V_abs):
    order.append(pos2[previous_point])
    previous_point = pos2[previous_point]
print('巡游次序:');print(0,end='')
for i in range(1,V_abs):
    print('->{}'.format(order[i]),end='')
print('->0',end='')
mpl.rcParams['font.sans-serif']=['Times New Roman']
plt.close('all')
fig,ax = plt.subplots(dpi=150)
ax.scatter(city_locations[:,0],city_locations[:,1],c='r',alpha=0.5)
for i in range(V_abs):
    ax.plot([city_locations[order[i]][0],city_locations[order[i+1]][0]],
            [city_locations[order[i]][1],city_locations[order[i+1]][1]],c='b')
ax.set_title('TSP Tour');ax.set_xlabel('x');ax.set_ylabel('y')
plt.show()

在这里插入图片描述

求解结果，最优解与TSPLIB公布的最优值7542基本相近，比咱之前用遗传算法求解的结果好不少，这么看来Gurobi确实厉害，实验结束。

Optimal solution found (tolerance 1.00e-04)
Best objective 7.544365901904e+03, best bound 7.544365901904e+03, gap 0.0000%
最小巡游距离:7544.37
巡游次序:
0->21->30->17->2->16->20->41->6->1->29->22->19->49->28->15->45->43->33->34->35->38->39->36->37->47->23->4->14->5->3->24->11->27->26->25->46->12->13->51->10->50->32->42->9->8->7->40->18->44->31->48->0