CSP问题

CSP是蛮经典的np-hard问题。一维CSP就是说：已知原料木材长度是，我现在想加工出个长度为的木材，我应该如何加工能够使我所使用的原料木材的数量最少？我们想用列生成算法解这个问题，可以把这个问题建立成一个集覆盖模型。

Set

  • ：所需木材的集合，
  • ：切割方案的集合，

Parameter

  • ：切割方案中得到的木材的数量
  • ：木材的需求量

varible

  • ：binary，是否使用切割方案

因此，CSP问题的模型可以写成：

为什么说这个模型是集覆盖？因为可以改写成，这就是集覆盖约束。

但是MP这个模型是没办法直接解的，因为我们没有这个集合，通过枚举的方式得到所有和对于大规模问题是不现实的。这种情况使用列生成逐步把它们生成出来是最好的。我们只考虑的一个子集，将松弛成连续变量得到限制主问题RMP的数学模型：

令表示集覆盖约束的对偶值，则RMP的reduced cost是。定义表示原料木材的长度；定义表示表示木材的长度；定义表示所需原料木材数量的上界，。从而切割方案对应的可以通过下面这个模型生成，我们把这个模型叫做定价子问题：

SP实际上是一个背包问题，因为可以写成，最优解是一致的。把看成是物品的价值，看成物品的重量，看成是背包的容量，至此建模工作全部完毕。

Colossus框架使用方法

下面我们借助Colossus框架，通过一两百行代码实现一维CSP的近似SOTA的结果。

规模 最优值（根数） 求解时间
大规模（100种需求） 607 14.6 s
中等规模（50种需求） 244 7.0 s
小规模（15种需求） 24 1.4 s

Colossus框架

截止2026年7月23日，Colossus框架包含以下特性：

  1. 自定义工具箱管理、列池管理、简洁可拓展的列定义
  2. 列生成框架支持mini-batch策略、对偶平滑策略等，提供便捷的求解结果访问接口

Colossus框架目前还在持续开发中，后期考虑开源。 

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Author  : syuansheng (Dalian Maritime University)

"""
example：使用Colossus框架求解切割下料问题
"""

import sys
sys.path.append('../../')
from math import ceil,floor
from itertools import accumulate
from gurobipy import Model,GRB,quicksum
from matplotlib import pyplot as plt
from Colossus import Column,ColumnPool,ToolBox,Environment,Engine


# 小算例
I=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
l={1:80,2:95,3:125,4:160,5:210,6:245,7:280,8:320,9:375,10:410,11:460,12:520,13:610,14:705,15:850}
d={1:5,2:6,3:4,4:7,5:3,6:8,7:5,8:6,9:4,10:7,11:3,12:5,13:4,14:2,15:3}
L=1000
M=sum([ceil(d[i]/floor(L/l[i])) for i in I])

# 中等算例
# I = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50]
# l = {
#     1:85, 2:92, 3:110, 4:125, 5:135, 6:148, 7:155, 8:165, 9:175, 10:190,
#     11:205, 12:215, 13:225, 14:235, 15:245, 16:255, 17:265, 18:280, 19:295, 20:310,
#     21:325, 22:340, 23:355, 24:370, 25:385, 26:400, 27:420, 28:440, 29:460, 30:480,
#     31:505, 32:530, 33:550, 34:570, 35:590, 36:610, 37:635, 38:660, 39:685, 40:710,
#     41:740, 42:770, 43:800, 44:830, 45:855, 46:880, 47:905, 48:930, 49:955, 50:980
# }
# d = {
#     1:10, 2:9, 3:12, 4:8, 5:11, 6:14, 7:6, 8:13, 9:7, 10:10,
#     11:12, 12:9, 13:15, 14:8, 15:11, 16:10, 17:13, 18:9, 19:14, 20:12,
#     21:8, 22:11, 23:10, 24:15, 25:9, 26:12, 27:14, 28:8, 29:11, 30:10,
#     31:13, 32:9, 33:12, 34:10, 35:15, 36:8, 37:11, 38:14, 39:10, 40:9,
#     41:12, 42:11, 43:8, 44:10, 45:9, 46:13, 47:7, 48:11, 49:10, 50:12
# }
# L = 1000
# M=sum([ceil(d[i]/floor(L/l[i])) for i in I])

# # 大规模算例
# I = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,
# 32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,
# 63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100]
# l = {
#     1:35, 2:42, 3:55, 4:68, 5:72, 6:88, 7:95, 8:105, 9:118, 10:125,
#     11:132, 12:145, 13:155, 14:162, 15:175, 16:185, 17:192, 18:205, 19:215, 20:225,
#     21:235, 22:245, 23:258, 24:265, 25:275, 26:285, 27:295, 28:305, 29:318, 30:325,
#     31:335, 32:345, 33:358, 34:365, 35:375, 36:385, 37:395, 38:408, 39:415, 40:425,
#     41:435, 42:445, 43:458, 44:465, 45:475, 46:485, 47:495, 48:505, 49:515, 50:525,
#     51:535, 52:545, 53:558, 54:565, 55:575, 56:585, 57:595, 58:608, 59:615, 60:625,
#     61:635, 62:645, 63:658, 64:665, 65:675, 66:685, 67:695, 68:705, 69:715, 70:725,
#     71:735, 72:745, 73:758, 74:765, 75:775, 76:785, 77:795, 78:808, 79:815, 80:825,
#     81:835, 82:845, 83:858, 84:865, 85:875, 86:885, 87:895, 88:905, 89:915, 90:925,
#     91:935, 92:945, 93:955, 94:965, 95:975, 96:982, 97:988, 98:992, 99:996, 100:998
# }
# d = {
#     1:15, 2:14, 3:16, 4:13, 5:17, 6:12, 7:18, 8:14, 9:15, 10:13,
#     11:16, 12:14, 13:15, 14:13, 15:17, 16:12, 17:18, 18:14, 19:16, 20:13,
#     21:15, 22:14, 23:16, 24:12, 25:17, 26:13, 27:15, 28:14, 29:16, 30:12,
#     31:14, 32:13, 33:15, 34:12, 35:16, 36:14, 37:13, 38:15, 39:12, 40:14,
#     41:13, 42:12, 43:15, 44:14, 45:13, 46:16, 47:12, 48:14, 49:13, 50:15,
#     51:12, 52:11, 53:13, 54:12, 55:14, 56:11, 57:12, 58:13, 59:11, 60:12,
#     61:11, 62:12, 63:13, 64:10, 65:12, 66:11, 67:13, 68:10, 69:12, 70:11,
#     71:10, 72:11, 73:12, 74:10, 75:11, 76:12, 77:10, 78:11, 79:12, 80:10,
#     81:9, 82:10, 83:11, 84:9, 85:10, 86:11, 87:9, 88:10, 89:11, 90:9,
#     91:8, 92:9, 93:10, 94:8, 95:9, 96:10, 97:8, 98:9, 99:10, 100:8
# }
# L = 1000
# M=sum([ceil(d[i]/floor(L/l[i])) for i in I])

# 定义如何生成初始可行方案
def generateInitialCols():
	columnpool = ColumnPool()
	check_dict = _get_omega_use_for() # 方案omega是针对哪个木材切的呢
	for omega in range(1,M+1): 
		column = Column()
		column.omega = omega
		column.tau=1
		for (lb,ub),idx in check_dict.items():
			if lb<=omega<=ub:
				# omega是为需求i=idx切的
				column.coefficients = [0 if i!=idx else floor(L/l[i]) for i in I]
		columnpool.addCol(column)
	return columnpool

def _get_omega_use_for():
	lst = [ceil(d[i]/floor(L/l[i])) for i in I] # 原料木材只切一种木材时各种木材所需原料木材的数量构成的列表
	check_dict = {}
	accumulate_lst = accumulate(lst)
	for idx,ub in enumerate(accumulate_lst,start=1):
		check_dict[(ub-lst[idx-1]+1,ub)] = idx
	return check_dict

# 定义限制主问题的数学模型
def buildRestrictedMP(columnpool):
	rmp_model = Model()
	labda = rmp_model.addVars([omega for omega in range(1,columnpool.col_num+1)],vtype=GRB.CONTINUOUS,name='labda')
	rmp_model.setObjective(
		labda.sum('*')
		,sense = GRB.MINIMIZE
		)
	a = {}
	check_dict = _get_omega_use_for()
	for omega in columnpool.keys():
		for i in I:
			a[omega,i] = columnpool[omega].coefficients[i-1]
	rmp_model.addConstrs(
		(quicksum(a[omega,i]*labda[omega] for omega in range(1,columnpool.col_num+1))>=d[i] for i in I)
		,name = 'cons_rmp'
		)
	return rmp_model

# 定义如何通过精确解算法求解fingerprint对应的定价子问题获得改善列
def updateAndSolveSP_Exact(y,fingerprint,omega_next_col):
	# 因为这个定价子问题只有1个，所以fingerprint参数实际上用不上
	# 下面用动态规划算法做精确求解，寻找改善列
	improve_column_lst = []
	item_price,item_weight = {},{}
	for i in I:
		item_price[i],item_weight[i] = y[i-1],l[i] # 物品i的价格和重量
	knapsack_solution,knapsack_obj = _solve_knapsack(L,item_set=I,item_price=item_price,item_weight=item_weight)
	if (1-knapsack_obj) > -1e-6:
		# 认为找不到改善列
		return improve_column_lst
	else:
		# 这个算法只能单列生成
		column = Column()
		column.omega = omega_next_col
		column.reduced_cost = 1-knapsack_obj
		column.tau = 1
		column.coefficients = list(knapsack_solution.values())
		improve_column_lst.append(column)
		return	improve_column_lst

def _solve_knapsack01(capacity,item_set,item_price,item_weight):
	"""
	一维01背包问题的动态规划算法
	"""
	knapsack_solution = {item:0 for item in item_set} # 各个物品选多少
	knapsack_obj = -float('inf') # 最大价值
	f = {item:{} for item in item_set} # 指标函数的最优值
	# 设置边界条件
	f[item_set[-1]+1]={} 
	for w in range(0,capacity+1):
		f[item_set[-1]+1][w]=0
	# 根据递推方程做逆推计算
	for item in reversed(item_set):
		for w in range(0,capacity+1):
			part1 = f[item+1][w]
			part2 = -float('inf') if (w-item_weight[item])<0 else item_price[item]+f[item+1][(w-item_weight[item])]
			f[item][w]=max([part1,part2])
	# 记录最大价值
	knapsack_obj = f[1][capacity]
	# 回溯得到01背包最优解
	left_capacity = capacity
	for item in item_set:
		if f[item+1][left_capacity]<item_price[item]+(-float('inf') if left_capacity-item_weight[item]<0 else f[item+1][left_capacity-item_weight[item]]):
			knapsack_solution[item]=1
			left_capacity = left_capacity-item_weight[item]
		else:
			knapsack_solution[item]=0
	return knapsack_solution,knapsack_obj

def _solve_knapsack(capacity,item_set,item_price,item_weight):
	"""
	完全背包求解思路：把完全背包转换为01背包再调用01背包的动态规划方法解，最后将结果转换成完全背包的结果
	"""
	max_item_num = {} # 背包装单种物品最多装几个
	for item in item_set:
		max_item_num[item]=floor(capacity/item_weight[item])
	# print(max_item_num)
	# 完全背包转换成01背包
	len_dummy_item_set = sum([num for num in max_item_num.values()])
	dummy_item_set = [dummy_item for dummy_item in range(1,len_dummy_item_set+1)]
	# print(dummy_item_set)
	dummy_item_price,dummy_item_weight = {},{}
	counter = 1
	for item in item_set:
		dummy_item_num = max_item_num[item]
		for dummy_item in range(counter,counter+dummy_item_num):
			dummy_item_price[dummy_item]=item_price[item]
			dummy_item_weight[dummy_item]=item_weight[item]
		counter+=dummy_item_num
	# print(dummy_item_price)
	# print(dummy_item_weight)
	# 调用01背包方法求解
	knapsack_solution,knapsack_obj = _solve_knapsack01(capacity,dummy_item_set,dummy_item_price,dummy_item_weight)
	# print(knapsack_solution,knapsack_obj)
	# 将knapsack_solution调整回完全背包的解
	new_knapsack_solution = {}
	counter = 1
	for item in item_set:
		new_knapsack_solution[item]=0
		dummy_item_num = max_item_num[item]
		for dummy_item in range(counter,counter+dummy_item_num):
			if knapsack_solution[dummy_item]==1:
				new_knapsack_solution[item]+=1
		counter+=dummy_item_num
	return new_knapsack_solution,knapsack_obj

# 定义如何通过启发式法求解fingerprint对应的定价子问题获得改善列
def updateAndSolveSP_Heuristic(y,fingerprint,omega_next_col):
	# 实现思路：贪婪策略，优先装单位价值高的物品直至装不下再换下一物品
	improve_column_lst = []
	solution,obj = {},0
	item_price,item_weight,unit_price = {},{},{}
	for i in I:
		item_price[i],item_weight[i],unit_price[i] = y[i-1],l[i],y[i-1]/l[i] # 物品i的价格、重量以及单位重量的价值
	ordered_item_lst = sorted([item for item in unit_price.keys()],key=lambda item:unit_price[item],reverse=True)
	left_capacity = L
	for item in ordered_item_lst:
		solution[item]=floor(left_capacity/l[item])
		left_capacity = left_capacity-solution[item]*item_weight[item]
		obj+=solution[item]*item_price[item]
	if (1-obj)>-1e-6:
		return improve_column_lst
	else:
		column = Column()
		column.omega = omega_next_col
		column.reduced_cost = 1-obj
		column.tau = 1
		column.coefficients=[]
		for i in I:
			column.coefficients.append(solution[i])
		improve_column_lst.append(column)
		return improve_column_lst

if __name__ == '__main__':
	toolbox = ToolBox()
	toolbox.addTool('generateInitialCols',generateInitialCols)
	toolbox.addTool('buildRestrictedMP',buildRestrictedMP)
	toolbox.addTool('updateAndSolveSP_Exact',updateAndSolveSP_Exact)
	toolbox.addTool('updateAndSolveSP_Heuristic',updateAndSolveSP_Heuristic)
	environment = Environment(['labda'],['INTEGER'],[1])
	environment.col_file_output_on=False
	engine = Engine(environment,toolbox)
	engine.runEngine()
	# plt.plot(engine.result.history_incumbent)
	# plt.show()