# Java_OR

**Repository Path**: mizuki114/java_or

## Basic Information

- **Project Name**: Java_OR
- **Description**: Java学习运筹学代码
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 1
- **Forks**: 0
- **Created**: 2021-10-11
- **Last Updated**: 2023-03-05

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# Java运筹学算法

## 单纯形法（Simplex Algorithm）

整体来说比较简单，该代码目前存在以下几个问题：

1. 模型必须是标准型
2. 默认把最后几个变量当作初始基变量
3. 求解过程未考虑无解、无界解和无穷多解的情况
4. 浮点运算存在精度损失问题


## 修正单纯形法（Revised Simplex Algorithm）

相较于单纯形法，计算上更加复杂一些，需要理清楚逻辑关系，代码写的很烂，因为读文件的代码是复制单纯形法的，矩阵运算用了一个非常轻量化的Jama，但是Matrix好像不能根据ArrayList生成矩阵，所以很多数据都是ArrayList一套，Matrix一套。存在的问题和单纯形法基本上是一样的，所以不再写了。

## 列生成法（Column Generation）

### 一个没有列生成特点的列生成

这一次的代码花的时间比较长，一方面是中间有一个礼拜出差了，另一方面是在一个我意料之外的地方卡了一两天（下面细说）。

列生成法总体思路比较简单，就是把修正单纯形法（RSA）里计算进基变量的方法由逐个变量计算检验数改成了由求解器计算检验数，所以大部分代码还是直接复制RSA的，但是这一次改掉了RSA代码中冗余的部分，并调整了逻辑结构（我还发现RSA代码可能写错了。。。。），通过构造器重载的方法同时实现了RSA和CG，并进行了对比（艹感觉在写论文）。

下面简单讲讲写代码中的一些有意思的点：

1. 意外的问题

上课时讲列生成用的例子就是cut stock problem，所以代码里直接以这个问题作为输入，并生成待求解的模型。输入文件第一行表示木料的长度，第二行表示可切成的短木料的长度，第三行表示每个短木料的需求量。万万没想到遇到的第一个问题是如何根据这个输入生成对应的模型。理论上用穷举最简单，但问题是穷举需要用`for`嵌套，而有几个规格的短木料是输入的，也就是说有几层`for`在写代码的时候是不知道的。最后用了一个指针数组来解决这个问题：

> 首先用数组记录每一个规格切割根数的所有可能性，然后生成一个与规格个数相等的指针数组`int[] p`，每一个指针指向对应规格当前应切割的根数，控制第一个指针不断右移，并将每一个指针指向的数拼成一个数组（即对应一个方案）。当第一个指针指向最后一个数时，第二个指针右移一位，然后第一个指针重新指向首位。不断重复这个过程，当所有指针指向对应数组的最后一个数时，每一个情况就遍历完了。

这个过程有点像时钟，如果把每一个切割根数的数组首尾相接变成一个环，指针跑起来的过程就很像时针、分针、秒针。每当前一个指针走完一圈，下一个指针移动一位。我也不知道有没有更简单的方法解决这个问题，等后面如果刷leetcode的时候再看吧。

2. 关于列生成中确定进基变量方法的疑问

在使用求解器算出进基变量的系数列向量后，如何寻找对应的变量呢？我暂时能想到的就是通过循环，遍历每一个变量的列向量，看他们是否相等。但其实计算检验数也是一个循环，然后做一些简单的四则运算，所以两边的时间权衡就是`(求解器计算 + 通过循环找到对应变量) - 通过循环计算检验数`。所以在变量个数少的情况下，列生成的时间可能大于修正单纯形法。

此处还有一个问题，那就是在比较两个列向量是否相等的方法上，这个方法可能会影响找变量的时间。我目前想到了两种方法。一就是循环比较每一个分量，如果所有分量都相等那就相等。但是这个方法受基变量个数的影响，基变量个数等于约束个数，还等于规格个数，也就是说如果有很多规格，那这个计算方法也可能很慢。所以我又想到了第二个方法：首先将两个向量相减，然后求这个得到向量的1-范数，只要不为0那就是不相等。矩阵计算我用的是Jama，我没有看过源代码，不过我想这个方法虽然比较花里胡哨，但是可能会快一些。我还顺便做了一手对比：

| 项目  | 1    | 2    | 3    | 4    | 5    | 6    | 2-6均值 |
| ----- | ---- | ---- | ---- | ---- | ---- | ---- | ------- |
| RSA   | 21   | 14   | 15   | 13   | 14   | 13   | 13.8    |
| CG(n) | 26   | 34   | 30   | 22   | 25   | 23   | 26.8    |
| CG(i) | 27   | 33   | 32   | 34   | 27   | 22   | 29.6    |

> CG(n)表示使用范数对比，CG(i)表示使用循环对比，计算平均值时去掉了第一次的结果，因为感觉第一次的结果偏差有点大（指RSA）

结果非常的amazing啊，数据量小的时候，RSA果然还是快于CG，而使用范数对比的速度哪怕在只有3种规格的时候仍然略快于使用循环。

总的来说，列生成的代码并不是很难，这个代码有一个未解决的问题，那就是：虽然以cut stock problem为模板，但求出来的不一定是整数解。暂时我还不知道怎么求整数解，后面再看吧。

### 基于列生成求解下料问题（Cutting Stock Problem）

上次写完列生成的代码之后回寝室洗澡，明明应该很开心，但总觉得哪里不对劲。秦虎讲课的样子在我脑海里不断闪过：我不需要知道每一个系数是多少……最大的优势在于我只取我需要的列……列生成……生成……

我突然想起了什么：我现在的代码里哪里有生成需要的列的过程呢？反而是用花里胡哨的指针先把所有的情况列出来了，然后还需要逐一对比列向量选取进基变量，或许列生成的精髓我并没有领悟。

过去的几天一直在写论文，终于今天把论文初稿写完了，在经过了快乐的摸鱼下午，晚上又打开代码想这个问题。我之前的思路是先有下料问题的各个参数，然后通过参数生成数学模型，然后调用修正单纯形法，只是把计算出基变量系数的过程交给Gurobi。如果真的按照秦虎上课所说，不需要知道完整的模型也可以求解，那就需要绕过第二步，由下料问题的参数直接进入求解过程。真的可以吗？

真的可以。下料问题其实有一些特殊的性质。比如初始基可行解可以直接构造出来，每一个决策变量的c都是1。基于这两个性质，完全可以绕过模型，直接求解。

具体步骤如下：首先，根据下料问题的参数，可以直接构造出初始的B逆（每一个方案只切一种长度），然后可以直接丢给Gurobi计算出可以节约的方案对应的决策变量的系数（虽然此时根本没有决策变量这个概念了，因为连模型都没有）。根据这个算出来的系数和b'可以直接算出基变量中的哪一个该出基（只需要知道其在B逆中是第几列即可），然后更新B逆，进入下一次迭代。当Gurobi的最优解小于等于0时即达到最优解，此时将B逆再取逆即可得到B，也就是可以直接反应切割方案的系数矩阵。至此，我们并不需要列出所有方案，也不需要根据Gurobi的输出去判断哪个变量换入（直接把上一节的两个问题绕过去了，感觉上一节白写了），只需要维护一个spec个数阶的矩阵即可求解。

整个代码异常清爽，所有代码包括读入文件的加起来也才200行，要知道我上一次光生成所有情况都花了200行，所有代码加起来更是有将近500多行。我也统计了一下计算时间（纯计算的时间，不包括读入文件及生成所有情况），可以看到同样是列生成，这一次比上一次快得多，虽然在小规模下还是无法追上纯修正单纯形法，但在大规模数据下肯定是列生成有优势。Gilmore and Gomory 当年发明这个方法的时候所面对的下料问题，有40多种spec，超过1亿种切割方式。在这个规模下，维护一个40阶的矩阵可比维护一个40×1亿阶的系数矩阵要快得多了（更不谈先要把这个矩阵生成出来）。

| 项目  | 1    | 2    | 3    | 4    | 5    | 6    | 2-6均值 |
| ----- | ---- | ---- | ---- | ---- | ---- | ---- | ------- |
| RSA   | 21   | 14   | 15   | 13   | 14   | 13   | 13.8    |
| CG(n) | 26   | 34   | 30   | 22   | 25   | 23   | 26.8    |
| CG(i) | 27   | 33   | 32   | 34   | 27   | 22   | 29.6    |
| CSP   | 23   | 24   | 24   | 23   | 22   | 24   | 23.4    |

关于解不是整数的问题，在一本书中说，只需要把LP问题中非整数的决策变量向上取整即可得到一个近似最优解。或许直接求最优解太困难了吧，那我也就不再纠结这个问题了。

### 更加常见的列生成

原本我以为已经可以进入DW分解了，结果上网查资料发现网上的DW分解在列生成的步骤上跟我学的不一样，于是被迫又重新研究了一下列生成。

在秦虎讲的版本中，主问题需要我们自己写代码求解，要维护主问题的B_inverse，根据比值规则选择出基变量，只是把子问题给求解器计算。但是更加一般地，主问题也可以交给求解器计算，这样就相当于去掉了修正单纯形法的过程。

仍以下料问题为例，对于主问题的系数矩阵A来说，每一行对应一个成品的数量约束，每一列对应一种切割方案。如果我们把所有方案列出来，主问题的最优解就是这个问题的最优解。但是如果我们只给出一部分方案，此时主问题的最优解就是在这些方案中选择一个最优的，这种主问题成为受限制的主问题（Restricted Master Problem, RMP）。这时就该子问题登场了。子问题的解对应了一种切割方式，且只要子问题的最优值（Reduced Cost）大于0，这种切割方式一定比现在主问题中的那些更好。那么我们只要把子问题的解作为新的一列添加到主问题的A中，下一次求解主问题就可以选择这种方式切割。这样随着切割方式的增加，每一次主问题都会在现有的方案里选择最优的几个，然后求解子问题看是否存在更优的，然后把更优的方案添加到主问题的可选方案中（即在A中新增一列），不断迭代直到无法提升主问题的解，这时主问题的最优解即为全局最优解。

关于初始的A，其实可以随意指定一个可行解，比如只切一种规格，每种规格只切一根（例如17的木料，3、5、9的成品长度，每次只切3（3×1根）、5（5×1根）、9（9×1根），剩余的14、12、8的长度直接扔掉）。当然这个初始解的质量会影响求解速度，比如还是只切一种规格，但每种规格切到不能切为止（切15（3×5根）、15（5×3根）、9（9×1根），剩余的2、2、8的长度扔掉），这个解显然比之前那个好。

理解了上面的过程，算法本身就很简单了。

> 注意：子问题的目标函数中需要用到主问题的对偶变量，注意只有连续变量的线性规划才可以求对偶变量，所以主问题变量必须是连续的

## DW分解（Dantzig-Wolfe Decomposition）

在列生成熬了三版算法，终于到了DW分解。不过也得益于对列生成算法的一层一层深入理解，DW分解只花了半天就写完了。主体框架都是基于最终一版的列生成代码，只是在表示主问题和子问题的时候使用了DW分解的方法而已。

DW分解基于以下几个关于线性规划的定理：

1. 可行域内任意一个解可以表示为该可行域的极点（若可行域无界，还需加上极射线）的线性组合
2. 最优解一定在极点处取得

虽然任意一个解都可以表示为极点的线性组合，但对于一个大型的问题，我们可能根本无法知道其所有的极点。假设我们能列出所有极点，那么我们就能表示出整个可行域，但如果只知道部分极点，那我们就能表示出可行域的一部分。这就是主问题需要求解的内容：根据给出的极点，找出这个子区域内最好的解。找完之后根据流程，计算Reduced Cost，我们要在剩余所有未选择的极点内（虽然我们不知道这些极点的坐标，甚至都不知道它们的数量）选择一个能使Reduced Cost最大（若为最大化问题则为最小）的极点，这就是我们子问题需要求解的内容：在剩余极点中选择一个使Reduced Cost最大或最小的点。结合上面的定理2，意味着只要列出当前Constraint Set的所有约束，那么该问题的最优解一定是某个极点。这样我们再把这个新的极点加入主问题，重新求解主问题。当所有Constraint Set的子问题都无法改进主问题的解的时候，我们就得到了最优解，尽管这时候我们仍然不知道每一个极点的坐标，但已知的极点所围成的子区域已经覆盖了最优解，并且我们已经得到这个点了。

## Benders分解（Benders Decomposition）

想不到Benders分解也这么快就写完了，当初上课的时候这玩意儿可折磨了我好久，幸亏考试不考。这次的问题背景不再是一个实际问题了，因为秦虎的PPT上给的就是这个例子，所以整个代码就是围绕这个已知的模型写的。

Benders分解的具体原理就不多解释了，这里位置太小了写不下。大致就是通过求解子问题的对偶问题来给主问题添加约束，限制主问题中可行域的范围，相当于慢慢把最优解所在的区域给切出来。假设原问题是最小化问题，如果对偶子问题（最大化问题）有最优解，就根据极点添加一条Optimality Cut，只留下会使得对偶子问题最优解更大（即对应原子问题最优解更小）的区域；如果是无界解，说明存在一条极射线，沿着它的方向目标函数值会一直增大下去，那么就根据极射线添加一条Feasibility Cut，排除使对偶子问题无界解（即原子问题无解）的区域。不断进行下去直到子问题的原问题与对偶问题有同样的最优值，即得到最优解。

下个礼拜，大概就要开始学习分支定界了吧，如果论文的修改意见还没下来的话。

## 分支定界法（Branch and Bound）

老实讲我原本打算跳过Branch&Bound直接进入Branch&Price，结果发现步子迈大了扯着蛋了，Branch&Price里有的例子我是真不知道那些price模型是怎么来的，感觉要学的还有很多，所以先写一个B&B练练手吧。

因为一些意料之外的原因，我现在在校外的一家酒店隔离，十几天之后回学校了还要再隔离14天，估计就直接寒假了，所以接下来一个月我只能在隔离的酒店里用一台隔离前临时借来的笔记本写代码了。

总的来说，Branch&Bound还是比较简单的，无论是数学过程还是编码过程，我原本以为需要用到Callback，后来发现解这种简单的模型，如果手写B&B根本不用Callback。算法最核心的就是`solveModel`函数和`generateSubproblems`函数了，这两个函数互相调用，用递归的方式实现了对Branch形成的二叉树进行深度优先搜索的算法。

`solveModel`函数主要负责求解传给它的模型，记录最优解并选出最大的非整数变量，准备进行分支操作。`generateSubproblems`函数主要是针对`solveModel`传过去的变量下标获取其值进行分支，生成两个新的模型，并调用`solveModel`求解。当整个递归过程完成的时候，表示所有可能的情况都已经搜索完毕，输出记录的最优解即可。

其实还可以进行广度优先的搜索，但是先有广度优先的代码还是先有Branch&Price的代码，那就看我的心情了~

## 分支定价算法（Branch and Price）

### B&P for GAP

鸽了好几天（好像不止），终于在隔离点把第一版B&P算法写出来了（这个“第一版”好像有什么不好的预示 ^_^）

分支定价算法本质上就是把列生成和分支定界结合起来，先构造初始解，求解一个大型IP问题的LP松弛，并使用列生成产生可以改进当前解的方案（称作pricing）。求出LP的最优解后，进行分支，不断重复直到在LP问题中求出整数解，当所有分支都探索完毕后，最好的那个整数解就是原IP问题的最优解。

说起来好像挺简单，但是中文互联网上的资料确实比较少，不得已只能求助于英文资料。找了一篇1997年发在OR上的论文，用B&P求解广义分配问题（Generalized Assignment Problem，GAP）。学到一半感觉有一些问题没有搞明白，又去看了同作者1998年发在OR上的另一篇像教学一样的论文，勉强看懂到能写代码的程度，就动手了。当然对于B&P的数学原理还是有一些一知半解，看后面再慢慢补吧（如果不影响使用，不知道就不知道吧）。

其实这次代码写的很不好看，有一些内容并不是很清晰，有点像当初第一次写列生成的时候，这也是为什么我上面说“第一版”的一个原因。算法整体框架还是`solveModel()`和`branching()`两个函数递归调用完成分支的过程，整体流程在`solveModel()`里还是比较直观的。但是仍然有以下几个地方感觉写的比较复杂：

1. 在B&B的代码里我是将Gurobi模型作为参数传递，遇到分支就要复制整个模型，感觉这个过程比较麻烦而且有点消耗资源，所以这里我采用的是复制构建模型所需的内容（比如分配方案，已有分支条件等），但是这样就会有很多的代码行重复，其实我感觉可以把一个节点封装成一个对象（比如节点的方案、已有分支条件、当前解等），可能写起来会简洁很多
2. 在生成RMP的过程中，决策变量我按照习惯还是把索引一维化了（因为Java下Gurobi好像不原生支持多维下标），但是这就有一个问题，因为对于分配方案的存储，我采用是一个三维数组，第一维是机器索引，第二维存着一个机器下所有分配方案，第三维表示一个方案里每一个任务的分配与否，但是一维编号与这个方案的存放结构有冲突的地方，因为生成的方案是哪一个机器取决于列生成的结果，这个方案放入这个三维数组里是放在对应机器的所有方案的最后，但由于1号机器在2号前面诸如此类，新生成的1号机器的方案一定在2号机器所有方案的前面（这个跟遍历方案的顺序也有关，我直接就是一层一层遍历的），也就是说上一次2号机器的第一个方案的一维索引是k，那当1号机器生成了一个方案以后，2号机器第一个方案的一维索引就成了k+1，但是这个修改不能同步到对于分支变量的记录上（因为我记录分支的方法是分开记的，几号机器的几号方案），这个从二维到一维的对应关系其实是随着方案的生成一直在变的，我还没仔细想这个问题到底是一维索引的问题还是方案记录形式的问题，或许用链表记录方案更好？

其实这些只是“第一版”的一个原因，真正让我觉得这是第一版的，还是下面这个原因。

我按照论文里的方法，以class D的标准生成了大约5组不同规模的算例，用这个算法求解，出现了一个很奇怪的现象：每一组算例都可以求解，但是求解过程都是不断生成列，然后没有经过一次分支就直接得到一组整数解，且无法再通过pricing改进，而且这个解就是全局最优。按照我现在的理解，B&P应该是无法通过pricing改进后，得到的可能是一个小数解，再去分支，重新pricing。而现在的情况却是从根节点的初始解开始，不断pricing直接把最优的方案生成出来了，所以没有分支，这个方案是LP松弛问题的最优解，解又正好是整数，所以就是IP问题的最优解。

我不知道是问题特点就是如此（这也太扯淡了，5组算例虽然不多，但是都这样也太诡异了）还是我对算法流程的理解有问题。因此我下一步打算再研究研究那两篇论文，同时我想起来了我当时用CG求解下料问题的时候，最后得到的就是一个小数解，如果用B&P重新求解，应该就可以得到整数解了，正好可以拿来当我的“第二版”。

不过距离我刑满释放还有5天，这期间估计还得改论文，放出去后估计还有别的事要做，看能不能在寒假前把第二版写出来吧（毕竟不指望寒假放得早 ^_^）。