单纯形算法（部分枚举方法）

可以证明（1）线性规划问题的最优解如果存在，那么一定可以在可行域的极点上面找到；（2）方程的特解如果非负，则称为线性规划问题的一个基本可行解，任何一个基本可行解都对应一个极点，且若两个基本可行解只有一个基变量不一样，那么它们对应的一定是相邻的极点。所以，一种朴素的求解线性规划问题的方法就是枚举出所有的基本可行解，然后比较其对应目标函数值的大小，这其实就是枚举算法。但对于一个的系数矩阵A，最多能找到个基本可行解……它几乎是一个几何数级复杂的的算法，效率是非常低下的。

Fig1. 单纯性算法的几何直观

虽然很笨，但枚举算法的思想是很棒的，聪明的数学家对枚举算法进行优化，得到了部分枚举方法——单纯性算法。要实现部分枚举，最重要的一点就是判断当前的基本可行解是否最优，如果不是最优，怎么找到一个更优（至少不是更差）的解。

Fig2. 单纯性算法
Fig3. 对偶单纯性算法

单纯性算法的核心

假设我们有一个基本可行解,那么对于任意一个基本可行解,不妨将其对应分成。那么对应的目标函数值和对应的目标函数值之间的关系是： 习惯上一般把第二部分叫做检验数，从而，只要检验数非负，那么对于最小化问题就可以确定是最优解了，否则，就一定可以通过换基找到一个更好（至少不更差）的基本可行解。那换基又该怎么换呢？一般是取最小检验数对应的非基变量作为进基变量（目标函数往往降低的更快，但也不一定），然后最小比值规则确定出基变量和移动的步长（步长无穷大，说明问题无界解），做完这个操作，我们就移动到可行域的另一个与先前极点相邻的极点上了，如Fig1。

分支定界算法（部分枚举方法）

Fig4. 分支定界算法的步骤
Fig5. 分支定界树

整数规划问题也一样可以用枚举算法求解，但显然这种算法的效率非常低下。和单纯性算法一样，分支定界算法也是从枚举算法得到启发而实现的一种部分枚举方法。在单纯形法中实现部分枚举的思想是让解越来越好，而在分支定界算法中实现部分枚举的基本思想则是剪枝，终止对某些节点的分支。

分支定界算法的核心——如何剪枝、选择那个节点以及哪个变量分支

  1. 节点是整数解，剪掉，不需要再分支了。
  2. 节点是不可行解，剪掉，不需要再分支了，因为子问题可行域更小，也一定是不可行解。
  3. 通过bound来实现剪枝，对于最小化问题，当前最好的整数解对应目标函数值作为上界，所有松弛问题最优解对应的最小目标函数值作为下界，如果某个节点的松弛解目标函数值甚至比上界还大就可以不用划分了，因为可行域越来越小，子问题不可能又更好的解。
  4. 节点选择上推荐混合策略——没有找到可行解前先深度优先，快速找到可行解，方便使用bound剪枝。有整数行解之后使用Bset优先策略，即选择最有希望找到最优解的节点（比如最小化问题，直观上就是目标函数值更小的节点更有希望）。
  5. 常用的分支是对离整数最远的变量进行的（一般效果更好），分别添加如下约束，就可以构造两个子问题：

计算复杂性理论

计算复杂度包括空间复杂度和时间复杂度，一般我们更关心时间复杂度。时间复杂度的作用是定性分析算法运行时间和问题规模的关系，它是问题规模的函数。一般情况下我们分析算法的时间复杂度都是在最坏情况下考虑的：（1）首先计算算法的语句频度，即算法中的语句所需执行的总次数；（2）然后去找当趋于无穷时和同一个量级函数，再使用运算符去掉低阶项、和高阶项前的系数；（3）最终得到的就是算法的时间复杂度。

Fig6. 学习计算复杂性理论有什么用
Fig7. 几类问题之间的关系

P、NP、NPC和NP-Hard问题的定义

  • P问题是指存在多项式时间算法求解的问题。
  • NP问题是指不确定是否存在多项式时间算法求解，但是可以使用多项式时间算法验证一个解到底是不是该问题的解的问题。
  • 一个问题时NPC问题需要满足两个条件：（1）这个问题时NP问题（2）所有NP问题都可以用多项式时间的约化算法约化为该问题。
  • 如果一个问题满足“所有NP问题都可以用多项式时间的约化算法约化为该问题”，那么它就是一个NP-Hard问题。

已知的NPC和NP-Hard问题

求解P、NP、NPC、NP-Hard问题的难度是递增的，如果遇到了NPC、NP-Hard问题可以大胆尝试各种启发式算法，如果遇到P、NP问题则应该考虑是否存在好的精确算法。

Fig8. 已知NPC问题

割平面方法（以Gomory的割平面方法为例）

Fig9. 凸包示意图
Fig10. Gomory割平面方法

根据凸包的相关理论，如果能找到整数规划问题的凸包，那么求解对应的整数规划问题就等价于求解对应的线性规划问题，可是找到整数规划问题可行域(一群整数点)的凸包往往本身就是一个NP-Hard问题。我们往往只能找到一个近似的凸包，这也就是割平面算法要做的事情。一般加了割平面之后问题会变成对偶可行而原始不可行，这时非常适合用对偶单纯性算法进行求解，比用单纯性法重新算要快很多

凸包和有效不等式

  • 凸包：凸包conv(X)是只可以包含X的最小凸集。
  • 有效不等式：满足两个条件的约束——（1）能够排除在松弛问题可行域但是不在原问题可行域中的解（2）不能排除原问题的可行解
  • 强有效不等式：贴着整数点构成的facet（刻面）切的有效不等式（Gomory割不是强有效不等式，割的效果并不一定很好），这个方法还在研究，很多构造强有效不等式的方法不是公开的，这是商业求解器公司的商业秘密。

在实际中，单纯使用割平面算法求解整数规划问题还是很慢的，目前所有主流求解器的核心框架都是Branch and Cut（分支切割算法），也就是将分支定界算法和割平面算法结合起来。分支切割算法的基本思想简单来说就是将分支定界算法和割平面算法结合起来，在分支定界树的每个节点处插入若干个割平面，从而更快地去逼近整数规划问题的凸包，显著减少分支定界树中节点的数量。以Gurobi求解混合整数规划的算法框架为例，首先，Gurobi会选择一个节点，对这个节点进行预求解操作（所谓预求解就是删除模型中的冗余变量和约束，精简模型）。然后，Gurobi会对节点做松弛处理（如果是非根节点，就不需要再松弛了，因为本身就是松弛的）。根据节点松弛问题解的信息，Gurobi会调用内置的19种割平面算法，试图逼近可行域（指原整数规划问题的可行域）的凸包，然后再对节点进行求解。如果松弛问题求解结果是整数可行解那么就剪枝，如果是小数解，以Max问题为例，我们就得到了原整数规划问题最优目标函数值的一个上界，这时Gurobi为了算法能有效剪枝还会调用内部几十种的启发式算法快速找到一个较好的整数可行解（在Max问题中就是找下界，也叫做incumbent）。最后，执行分支操作，如此循环，直到所有节点都不再活跃（被剪掉）。

Fig11. Gurobi求解MIP的算法框架

列生成算法——以下料问题为例

有长度为9，14，16米且单价分别为5、9、10元的三种木材，市场有需求——30个4米的木材、20个5米的木材、40个7米的木材。要求在满足市场需求的情况下，最小化总成本，请问应该使用什么切割方案？一种直观的思路是枚举出所有可能的切割方案，将每种切割方案的使用次数作为决策变量，然后就可以将问题建立成一个简单的线性整数规划模型。但是这样建模存在两个问题：

  1. 当问题规模增大时，可能的切割方案会指数级增长，难以全部枚举。
  2. 并不是所有的切割方案都是好的方案，某些切割方案的使用次数在最优解一定是0。

这种在建模前期难以事先确定所有决策变量及其系数列向量的问题就非常适合使用列生成算法求解，列生成算法的基本思想是——首先建立一个主模型（MP，有时也成为限制主模型，记作RMP），这个主模型是原模型的线性松弛，并且只包含了部分决策变量，并保证可行。将主模型求解后，判断是否存在变量及其系数列向量加入主模型后能够使得主模型的解更优，如果能就加入，并重复该过程，否则就说明已经找到原模型线性松弛的最优解了。

如何判断是否存在变量及其系数列向量加入主模型后能够使得主模型的解更优

这块和咱本科学的灵敏度分析——判断一个新变量和系数列向量的加入对原模型的影响思路是一样的，就是去算这个新变量对应的检验数是否小于0（对于Min问题）。前面提到过检验数的计算公式是，其中是对偶列向量。对于最小化问题，咱们就可以以为决策列向量，以为目标函数和相关的约束构造子问题，判断是否存在这样的新变量及其系数列向量，如果存在则求出它们。

话不多说，咱们直接用上述下料问题说明列生成算法的步骤。

(1)构造主模型
解：设、、分别代表三种切割模式的使用次数，这三种切割模式分别是,,。

(2)求解主模型
计算出对偶列向量。

(3)构造所有子模型
解：设，，，。则第个子模型如下,k=1,2,3。

(4)求解子模型
判断是否存在子模型的目标函数值小于0，如果存在取最小目标函数值对应的加入主模型重复上述过程计算，否则主模型已经找到原模型线性松弛的最优解了，列生成算法终止。

(5)将最终的主模型转换为IP问题继续求解
这个过程可能会遇到一些问题，不能保证得到的就一定是整数最优解，只是一个较好的解，得用分支定价的方式解？？？，这块儿还是不太理解。

上述过程用gurobi实现的方式如下：

列生成算法求解下料问题
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Date    : 2025-07-25 14:34:23
# @Author  : syuansheng (Dalian Maritime University)

from gurobipy import *

class Data:
	def __init__(self):
		self.l = [4,5,7] # 需求类型
		self.L = [9,14,16] # 可用原木的长度
		self.c = [5,9,10] # 各种长度原木的价格

class CGModel:
	def __init__(self):
		self.data = Data()
		# 初始化主模型
		y = self.initiate_main_model()
		# 初始化所有子模型
		self.initiate_sub_model(y)

	def initiate_main_model(self):
		self.main_model = Model()
		self.x = self.main_model.addVars(3,vtype=GRB.CONTINUOUS,name='x')
		self.main_model.setObjective(
			5*self.x[0]+5*self.x[1]+5*self.x[2],
			sense = GRB.MINIMIZE
			)
		self.main_model.addConstr(2*self.x[0]>=30)
		self.main_model.addConstr(self.x[1]>=20)
		self.main_model.addConstr(self.x[2]>=40)
		self.main_model.optimize()
		# 获取对偶变量的值
		y = self.main_model.getAttr('Pi',self.main_model.getConstrs())
		return y

	def initiate_sub_model(self,y):
		self.submodels = []
		for k in range(len(self.data.L)):
			# k子模型
			submodel = Model()
			p = submodel.addVars(3,vtype=GRB.INTEGER,name='p')
			submodel.setObjective(
				self.data.c[k]-quicksum(y[i]*p[i] for i in range(3)),
				sense=GRB.MINIMIZE
				)
			submodel.addConstr(quicksum(p[i]*self.data.l[i] for i in range(3))<=self.data.L[k])
			self.submodels.append(submodel)

	def update(self):
		"""
		不断更新主模型和子模型并求解，直到得到原模型的最优解
		"""
		while True:
			best_ObjVal,best_P,submodel_idx = 0,None,0
			# 求解子模型
			for k,submodel in enumerate(self.submodels):
				submodel.optimize()
				if submodel.ObjVal < best_ObjVal:
					# 说明存在列加入主模型后使得主模型更好
					best_ObjVal = submodel.ObjVal # 子模型最优解目标函数值
					best_P = submodel.X # 子模型最优解
					submodel_idx = k
			# 退出
			if best_ObjVal >= 0 :
				break
			# 将新的列添加到主模型中，求解主模型并更新子模型
			else:
				col = Column(best_P,self.main_model.getConstrs()) # 创建要加入到主模型的列
				self.main_model.addVar(obj=self.data.c[submodel_idx],column=col,vtype=GRB.CONTINUOUS,name='new_x') # 将新的列加入主模型
				self.main_model.optimize() # 求解主模型
				self.main_model.write('ttest.lp')
				y = self.main_model.getAttr('Pi',self.main_model.getConstrs()) # 新的对偶变量值
				print('----------------------------',y,best_P,submodel_idx)
				self.initiate_sub_model(y) # 更新子模型
		# 将模型转换为对应的IP模型继续求解
		self.main_model.setAttr('VType',self.main_model.getVars(),GRB.INTEGER)
		self.main_model.optimize()


if __name__ == '__main__':
	model = CGModel()
	model.update()
	print(model.main_model.ObjVal)
	print(model.main_model.X)