数据¶
所提供的实例来自 Espinoza 等人,其格式如下:
一个 .cpit 文件存储与带约束的露天矿坑极限问题 (CPIT) 相关的数据:
块的数量
时间段的数量
资源的数量
折现率(用于计算块开采的利润)
对每种资源和每个时间段,该时间段内该资源消耗的上下界
对每个块,其利润
对每个块和每种资源,该块的开采对该资源的消耗(仅当非零时写入文件)
一个 .prec 文件存储每个块的前驱块。
建模思路¶
带约束的露天矿坑极限问题 (CPIT) 的 OptAgent 模型使用 集合决策变量。每个集合表示在某个时间段内开采的块。我们添加一个虚拟时间段来表示保持未开采的块。借助 partition 算子,我们确保每个块最多被开采一次。
对于给定类型的资源,一个 lambda 函数对所有被开采的块应用 sum 算子来计算该时间段内该资源的总消耗量。注意,该求和中项的数量在搜索过程中是变化的,集合的大小也随之变化。然后我们可以根据为问题定义的上下界对该数量施加约束。
然后,使用 find 算子,我们可以获取每个块被选定的开采时间段的索引。这使我们能够写出块之间的优先关系约束:每个块必须在与其前驱块的同一时间段或更晚的时间段被开采。
一个时间段的利润通过一个 lambda 函数对该时间段内被开采的块应用 sum 算子来计算。同样地,该求和中项的数量在搜索过程中是变化的,集合的大小也随之变化。我们最大化所有时间段的总利润。
Python 实现¶
from pathlib import Path
from optagent import OptModel, solve
def read_instance(instance_file):
lines = Path(instance_file).read_text().splitlines()
index = 2
nb_blocks = int(lines[index].split()[1])
index += 1
nb_periods = int(lines[index].split()[1])
index += 1
nb_resources = int(lines[index].split()[1])
index += 1
discount_rate = float(lines[index].split()[1])
index += 2
lower_bounds = [[None] * nb_periods for _ in range(nb_resources)]
upper_bounds = [[None] * nb_periods for _ in range(nb_resources)]
for _ in range(nb_resources * nb_periods):
resource, period, relation, value = (
lines[index].split()[0],
lines[index].split()[1],
lines[index].split()[2],
lines[index].split()[-1],
)
index += 1
target = lower_bounds if relation == "G" else upper_bounds
target[int(resource)][int(period)] = float(value)
index += 1
profit = [0.0] * nb_blocks
for _ in range(nb_blocks):
block, value = lines[index].split()
index += 1
profit[int(block)] = float(value)
index += 2
resource_use = [[0.0] * nb_blocks for _ in range(nb_resources)]
while lines[index] != "EOF":
block, resource, value = lines[index].split()
resource_use[int(resource)][int(block)] = float(value)
index += 1
precedence = {}
prec_file = Path(instance_file).parent / "prec" / f"{Path(instance_file).stem}.prec"
for line in prec_file.read_text().splitlines():
values = [int(value) for value in line.split()]
if values[1]:
precedence[values[0]] = values[2:]
discounted_profit = [[value / (1 + discount_rate) ** period for value in profit] for period in range(nb_periods)]
return nb_blocks, nb_periods, lower_bounds, upper_bounds, discounted_profit, resource_use, precedence
def main(instance_file, time_limit=60):
nb_blocks, nb_periods, lower_bounds, upper_bounds, discounted_profit, resource_use_data, precedence = read_instance(
instance_file
)
model = OptModel()
# 最后一个集合( period_nb_periods )是 虚拟时间段 ,表示"该矿块选择不开采"
period_sets = [
model.set(nb_blocks)
for period in range(nb_periods + 1)
]
periods = model.array(period_sets)
model.constraint(model.partition(period_sets))
resource_use = [model.array(values) for values in resource_use_data]
for resource, (lower, upper) in enumerate(zip(lower_bounds, upper_bounds)):
for period in range(nb_periods):
consumption = model.sum(
period_sets[period], model.lambda_function(lambda block: resource_use[resource][block])
)
if lower[period] is not None:
model.constraint(consumption >= lower[period])
if upper[period] is not None:
model.constraint(consumption <= upper[period])
extraction_period = [model.find(periods, block) for block in range(nb_blocks)]
for block, predecessors in precedence.items():
for predecessor in predecessors:
model.constraint(extraction_period[block] >= extraction_period[predecessor])
profit_per_period = []
for period in range(nb_periods):
profit = model.array(discounted_profit[period])
profit_per_period.append(model.sum(period_sets[period], model.lambda_function(lambda block: profit[block])))
objective = model.sum(profit_per_period)
model.maximize(objective)
solution = solve(model, time_limit_s=float(time_limit))
values = {'discounted_profit': objective.value}
print(f"Discounted profit = {values['discounted_profit']}; Status = {solution.feasible}")
return solution
INSTANCE_DIR = Path.cwd() / "instances"
solution = main(INSTANCE_DIR / "zuck_small.cpit", time_limit=5)