kittyzc的博客

MTZ是Miller-Tucker-Zemlin inequalities的缩写。除了定义是否用到边xij​外，还需要定义一个ui​用来表示此时车辆的当前载货量。注意这里x变量需要定义为有向。MTZ的求解速度不快，10个点3辆车都需要3秒左右时间。

详见：https://github.com/VROOM-Project/vroom/blob/master/docs/API.md需要定义如果没有指定经纬度和地图server的话，则需要定义matrices。

PyVRP使用HGS（hybrid genetic search）算法求解VRP类问题。

我们不妨假设仓库有3种物料需要搬运到生产线，搬运时间都是3min。三种物料假设都需要搬运3件，有2个叉车可以进行搬运，每个叉车每次智能搬运一件。

这里概述一下help文档中Supply chain GIS model例子的要点：触发事件的agent和执行任务的agent。在这个案例中，触发事件的agent是retailer，不断有订单生成；而执行任务的agent是vehicle，不断将订单运从distributor运送到retailer。模型的最终效果如下图：首先需要设计agent。在供应链的例子中，设计了4种agent。其中dis和order最简单，order就是个classret需要包含一个，用于生成order，寻找veh。

在constraint programming中，objective是不起作用的，调用optimize!(model)后，会给出满足约束的可行解。使用JuMP的约束规划API，可以极大简化模型。

看一个简单的数值优化的例子：我们将其建立为一个N阶段的问题：初始值为M。结果为非常接近理论最优值。

添加辅助变量y。比如Either3x1​2x2​≤18orx1​4x2​≤16可以用3x1​2x2​≤18Myx1​4x2​≤16M1−y来代替。

关于列生成的列子，可以参考《运筹系列8：Set partitioning问题的列生成方法》。我们这里使用priority queue存储分枝节点，按照最简单的下标顺序，对所有非整数变量进行分枝。求解主问题，然后找到检验数

我们得到reduce cost，可以作为tsp的lb = (10 + 10 + 15 + 20 + 5 + 10) = 70。我们用行表示从每个点出去的距离，用列表示从每个点进来的距离。

Node算子用来存储搜索树的状态。其中level等于path的长度，path是当前节点已经访问过的vertex清单，bound则是当前的lb。这里的bound函数是一种启发式方法，等于当前路径的总长度，再加上往后走两步的最小值。

定义csk为当前在k，待访问点的集合scskmini∈s​cs−iidik​2kk1​...kn​S2k1​...2kn​我们有初始状态c({0kd0k​0。逐步扩大set的尺寸，遍历所有可能的subset，使用bellman方程迭代计算。

基础的node定义在LKH.h中用于2-level tree的segment定义如下：LKH可以使用3种数据结构，默认是2-level tree：2-level tree的flip操作（即2-opt算子），在Flip_SL.c中，特殊的地方在于flip过程中涉及到重新分配segment的操作，其逻辑如下，其中P开头的是node对应的segment编号。

使用flow约束构建基础模型。

我们对cutting stock问题进行建模。

model的方法和solver的方法很多是一致的，比如model.getNumCols() 以及model.solver()->getNumCols()，但是同步性有时会有差别，比如getColSolution() ，CbcModel可能不是最新的结果，这时可以用CbcModel::bestSolution()来获取结果。MaxS : 默认Infinity : 做多允许存储的可行解数量。MaxN: 默认Infinity : 最多允许搜索的节点数。需要首先定义一个线性规划的solver，传入model中。

从本期开始阅读开源求解器coin_or的源码。整体项目介绍见https://blog.youkuaiyun.com/kittyzc/article/details/100182189第一篇从clp开始。

通过一个例子来认识：max x1 + 2 x2 + 3 x3 + x4s.t. − x1 + x2 + x3 + 10 x4 ≤ 20x1 − 3 x2 + x3 ≤ 30x2 - 3.5 x4 = 00 ≤ x1 ≤ 40，x2 ≥ 0，x3 ≥ 0，2 ≤ x4 ≤3对应的mps文件为：2. 文件解析使用pysmps解析文件，例如上面的例子的结果是：前四行是命名行，分别是：项目名、目标函数名、行名（约束条件名）、列名（变量名）接下来两行是变量类型和约束条件类型接

这个模块是输入距离矩阵c和pi值，给出mst1tree。其中minimum_spanning_tree模块，可以使用prim算法自己实现，cpu上性能和scipy差不多。

首先定义MIP问题如下：1. 圆整启发式算法圆整启发式算法可以看做是广度优先的搜索方法，尽量一次把变量都变为整数，有两个事情要做：选择圆整变量、选择上圆整还是下圆整。1.1 简单圆整若A.j≥0A_{.j}\ge0A.j​≥0，称变量xjx_jxj​平凡向下可圆整（即向下取整仍保持约束可行）；反之叫平凡向上可圆整。为方便描述，定义A.jA_{.j}A.j​中负系数的个数称为下锁数，正系数的个数称为上锁数，两者的最小值称为变量xjx_jxj​的锁数。依次搜索变量，遇到平凡向下可圆整，则向下进行圆整

19、26预处理文献1. 基本思路常见的松弛策略包括对整型变量松弛到连续空间 (NLP 子问题) 和非线性约束的线性化 (MILP 子问题)，非线性约束的线性化称为外逼近割平面：整数变量松弛：求解凸 MINLP 的算法通常采用迭代法和分支定界法 (branch-and-bound, BB) 两种思想。迭代法是不断更新子问题和迭代点列, 直到算法收敛。2. BB算法1960 年, Land 和 Doig首次将 BB 算法应用到 MILP 问题，1965年由 Dakin首次应用到求解凸 MI

1. 数据以及当前最优结果介绍网址：http://www.math.uwaterloo.ca/tsp/data/ml/monalisa.htmltsp地址：http://www.math.uwaterloo.ca/tsp/data/ml/mona-lisa100K.tsp目前最优解是2009年3月17日,由Yuichi Nagata找到的，地址为：http://www.math.uwaterloo.ca/tsp/data/ml/tour/monalisa_5757191.tour，使用的是自己设计的遗

1. 示例来个简单的例子（读取mps文件并求解）>>> from cylp.cy.CyClpSimplex import CyClpSimplex, getMpsExample>>> s = CyClpSimplex()>>> f = getMpsExample()>>> s.readMps(f)0>>> s.initialSolve()'optimal'再来个复杂的例子（从头开始构建）impor

1. LU分解LU分解(LU Factorization)可以将一个矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积。LU分解非常简单，下图可以直观的看出LU分解的求解方法。LU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式，起源于高斯消元法，也即所谓的杜尔里特算法（Doolittle algorithm）：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下

应用：源码地址：https://github.com/exanauts/Simplex.jl，使用julia实现。git clone https://github.com/exanauts/Simplex.jl下载源码，然后cd netlib,运行gcc emps.c -o emps运行./get.sh下载例子，如果下载不下来，按照下图修改代码：测试报告使用netlib.jl进行打印：...

另一种方式是将所有的点拿出来做v0，重新计算mst，速度会慢不少，代码如下。使用简单的2-opt算子进行优化，以一定的概率接受退化解。第一种下界非常简单，每个点取最短的两条边加起来即可。我们可以用类似训练深度学习的方法，来逐步调节所有的。

main.cpp- Edge assembly crossover，核心程序***.tsp还可以做一些自定义算法配置，主要在env。fStage = 1;fStage = 2;fStage = 1;TerminationCondition() 里可以修改停止条件。

纯C语言编译，用于对称TSP问题求解。目前最大的应用案例为85900个城市。Concorde支持使用QSopt线性规划求解器获得bound。

这是可以给出上界的一个方法，可以证明构造出的路线不超过最优路线的1.5倍。步骤为：构造MST（最小生成树）将里面的奇点连接起来构成欧拉回路称为完美匹配。Edmonds给出了多项式时间内构造最小代价完美匹配的方法，其长度不超过最优解的1.5倍。证明方法也很直观，奇点最短路径可以拆分成两条完美匹配，其中总有一条的长度≤最优奇点路径长度/2≤最优完整路径长度/2按欧拉回路顺序逐个扫描点，跳过重复经过的点，即可构造一条完整的路径。

1. 工具调研总结一下，faiss用起来最方便，使用的人也多；集团的proxima还待完善，如果要上生产需谨慎；几个月前开源的Milvus已经看到使用案例了，看起来效果不错；如果用blink，也可以使用alink的内置向量检索，使用LSH实现。2. 部署服务faiss是不支持java的，如果不是对java有什么执念的，建议用python搭个服务，gevent+gunicorn实测下来很稳定。非要用java的话，github上有个faiss4java，采用c版本rpc通信的形式。需要下载faiss源码

参见这篇访谈：https://www.oberlin.edu/bob-bosch文章见：https://www2.oberlin.edu/math/faculty/bosch/tspart-page.html1. TSP ART概述将灰度图片先转为离散图，然后将离散图用一笔画画出来转为离散点的方法，使用的是weighted Voronoi stippling方法，参考这篇论文：https://www.cs.ubc.ca/labs/imager/tr/2002/secord2002b/secord.

1. Rapids1.1 CuPy：用GPU加速Numpyimport cupy as cp N = 10 ** 7data = cp.random.uniform(-1, 1, size=(N, 2))inside = (cp.sqrt((data ** 2).sum(axis=1)) < 1).sum()pi = 4 * inside / Nprint('pi: %.5f' % pi)在GPU上仅用640us，在cpu上需要用757ms，提升效果异常显著。1.2 cuDF：

1. 基础模块1.1 读取数据标准的tsp格式可以用tsplib95库进行读取：import numpy as npimport tsplib95import matplotlib.pyplot as pltdata = tsplib95.load('ALNS/examples/xqf131.tsp')coords = np.array([(coord[0],coord[1]) for city, coord in data.node_coords.items()])solution = t

1. 2-opt方法2-opt本质是一个边交换启发式算法，非常简单。首先用myopic的方式生成一条路线，然后逐点使用2-opt算法进行优化。我们使用kaggle上的一个数据，https://www.kaggle.com/jsaguiar/fast-2-opt-with-cython/notebook?select=cities.csv，5000个点。1.1 基础版本import numpy as npimport pandas as pdcities = pd.read_csv("/Users

1. GPU硬件架构简单理解，GPU就是很多很多非常弱的cpu在做并行计算。个人桌面电脑CPU只有2到8个CPU核心，GPU却有上千个核心。在英伟达的设计理念里，CPU和主存被称为Host，GPU被称为Device。Host和Device概念会贯穿整个英伟达GPU编程。GPU核心在做计算时，只能直接从显存中读写数据，程序员需要在代码中指明哪些数据需要从内存和显存之间相互拷贝。这些数据传输都是在总线上，因此总线的传输速度和带宽成了部分计算任务的瓶颈。Intel的CPU目前不支持NVLink，只能使用PC

jitclassimport numpy as npfrom numba import int32, float32 # import the typesfrom numba.experimental import jitclassspec = [ ('value', int32), # a simple scalar field ('array', float32[:]), # an array field]@jitcla

1. mars原理2. mars基础mars可以作为替代numpy的一个库，安装使用都很简单：pip install pymars2.1 使用mars.tensor模拟numpy2.2 使用mars.dataframe模拟pandas2.3 eager模式用户可以通过选项启用 Eager 模式，在程序开头或控制台会话中设置选项。>>> from mars.config import options>>> options.eager_mode = T

1. 单纯形从方程组的方阵衍生开去，brute方式是遍历所有的方阵，单纯形法则是使用了逐步搜索的方法。一般的问题都是变量比约束多，因此竖着切一刀，非基变量=0；相对应的对偶问题则是横着切一刀，基变量对应的约束可以取到等号。当出现等式约束时，需要添加人工变量（相当于升维），使用两阶段法求解。对偶单纯形，将变量≥约束去除，很容易求得初始解。2. 内点法3. 整数规划分支定界法，对小数变量进行分支，不断减少求解空间。4. 分解法添加约束叫行生成（benders分解）添加变量叫列生成（DW分解）

1. 介绍性能上不必苛求了，毕竟用的都是python原生的包。基于NetworkX包（操作图）和cspy包（启发式算法）开发，NetworkX的性能其实不是很好。可以用来学习column generation的写法，以及一系列VRP的heuristic算法。支持如下类型的问题：the Capacitated VRP (CVRP),the CVRP with resource constraints,the CVRP with time windows (CVRPTW),the CVRP with

1. 介绍cspy是好几个列生成开源算法包中用于列生成算法的基础算法包，主要是用于生成新的列。最开始用于公交司机调度问题，后面被广泛用于交通相关的问题中。cspy是基于networkx的，性能上堪忧啊，可以考虑进行改造。2. 算法清单BiDirectional: Bidirectional and monodirectional algorithmsTabu Heuristic Tabu SearchGreedyElim Greedy Elimination ProcedureGRASP GR