数据¶
混流问题实例来自此 GitHub 仓库，并采用 JSON 格式：
“components”：表征原材料的对象数组，包含：
“name”：源节点的名称
“lower”：最小流出量（未使用，等于 0）
“upper”：最大流出量（供应量）
“price”：一个单位原材料的价格
“quality”：包含原材料中每种属性浓度的对象
“products”：表征最终产品的对象数组，包含：
“name”：目标节点的名称
“lower”：最小流入量（需求量）
“upper”：最大流入量（容量）
“price”：一个单位最终产品的销售价格
“quality_lower”：包含最终产品中每种属性最小浓度的对象
“quality_upper”：包含最终产品中每种属性最大浓度的对象
“pool_size”：包含池容量的对象
“component_to_product_bound”：表征源节点和目标节点之间边的对象数组，包含：
“component”：源节点 c 的名称
“product”：目标节点 p 的名称
“bound”：边 (c, p) 上的最大流量
“cost”：边 (c, p) 上每单位流量的成本
“component_to_pool_fraction”：表征源节点和池节点之间边的对象数组，包含：
“component”：源节点 c 的名称
“pool”：池节点 o 的名称
“fraction”：来自源节点 c 的池 p 流入量的最大比例
“cost”：边 (c, o) 上每单位流量的成本
“pool_to_product_bound”：表征池节点和目标节点之间边的对象数组，包含：
“pool”：池节点 o 的名称
“product”：目标节点 p 的名称
“bound”：边 (o, p) 上的最大流量
“cost”：边 (o, p) 上每单位流量的成本
原材料和最终产品的价格可以通过边上的 “price” 或 “cost” 字段表示。“cost” 表示用于 “randstd” 实例，而 “price” 表示用于其他实例。
Python 实现使用标准库 json 读取实例。
建模思路¶
混流问题的 OptAgent 模型对应于问题的 Q 形式，同时使用流量比例和流量量。我们首先声明三个浮点决策变量数组，分别表示：
- 从每个源节点 c 到每个目标节点 p 的流量；
- 从每个池 o 到每个目标节点 p 的流量；
- 每个池 o 中来自每个源节点 c 的流入量比例。
通过将后两个量相乘，我们计算出从每个源节点 c 经过每个池 o 到每个目标节点 p 的流量。
然后我们可以开始编写约束。首先，我们确保每个池中来自源节点的流入比例之和等于 1。然后，每个目标节点的总流入量必须满足需求同时不超出产品容量。类似地，每个源节点的总流出量不得超过原材料的供应量。最后，每个最终产品的每种属性浓度必须在一定的容差区间内。为了强制执行该约束，我们计算直接来自源节点和来自池的属性量，以及目标节点的总流入量。
最后，我们计算目标函数的值。总利润等于销售产品获得的总收入减去初始组件的成本再减去除流量成本。
Python 实现¶
from pathlib import Path
import json
from optagent import OptModel, solve
class PoolingInstance:
#
# Read instance data
#
def __init__(self, instance_file):
with Path(instance_file).open(encoding="utf-8") as problem:
problem = json.load(problem)
self.nbComponents = len(problem["components"])
self.nbProducts = len(problem["products"])
self.nbAttributes = len(problem["components"][0]["quality"])
self.nbPools = len(problem["pool_size"])
# Components
self.componentPrices = [problem["components"][c]["price"]
for c in range(self.nbComponents)]
self.componentSupplies = [problem["components"][c]["upper"]
for c in range(self.nbComponents)]
self.componentQuality = [list(problem["components"][c]["quality"].values())
for c in range(self.nbComponents)]
self.componentNames = [problem["components"][c]["name"]
for c in range(self.nbComponents)]
componentsIdx = {}
for c in range(self.nbComponents):
componentsIdx[problem["components"][c]["name"]] = c
# Final products (blendings)
self.productPrices = [problem["products"][p]["price"]
for p in range(self.nbProducts)]
self.productCapacities = [problem["products"][p]["upper"]
for p in range(self.nbProducts)]
self.demand = [problem["products"][p]["lower"]
for p in range(self.nbProducts)]
self.productNames = [problem["products"][p]["name"]
for p in range(self.nbProducts)]
productIdx = {}
for p in range(self.nbProducts):
productIdx[problem["products"][p]["name"]] = p
self.minTolerance = [[0 for _ in range(self.nbAttributes)]
if (problem["products"][p]["quality_lower"] == None)
else list(problem["products"][p]["quality_lower"].values())
for p in range(self.nbProducts)]
self.maxTolerance = [list(problem["products"][p]["quality_upper"].values())
for p in range(self.nbProducts)]
# Intermediate pools
self.poolNames = list(problem["pool_size"].keys())
self.poolCapacities = [problem["pool_size"][o] for o in self.poolNames]
poolIdx = {}
for o in range(self.nbPools):
poolIdx[self.poolNames[o]] = o
# Flow graph
# Edges from the components to the products
self.upperBoundComponentToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbComponents)]
self.costComponentToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbComponents)]
# Edges from the components to the pools
self.upperBoundFractionComponentToPool = [[0 for _ in range(self.nbPools)]
for _ in range(self.nbComponents)]
self.costComponentToPool = [[0 for _ in range(self.nbPools)]
for _ in range(self.nbComponents)]
# Edges from the pools to the products
self.upperBoundPoolToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbPools)]
self.costPoolToProduct = [[0 for _ in range(self.nbProducts)]
for _ in range(self.nbPools)]
# Bound and cost on the edges
for edge in problem["component_to_product_bound"]:
self.upperBoundComponentToProduct[componentsIdx[edge["component"]]] \
[productIdx[edge["product"]]] = edge["bound"]
if len(edge) > 3:
self.costComponentToProduct[componentsIdx[edge["component"]]] \
[productIdx[edge["product"]]] = edge["cost"]
for edge in problem["component_to_pool_fraction"]:
self.upperBoundFractionComponentToPool[componentsIdx[edge["component"]]] \
[poolIdx[edge["pool"]]] = edge["fraction"]
if len(edge) > 3:
self.costComponentToPool[componentsIdx[edge["component"]]] \
[poolIdx[edge["pool"]]] = edge["cost"]
for edge in problem["pool_to_product_bound"]:
self.upperBoundPoolToProduct[poolIdx[edge["pool"]]] \
[productIdx[edge["product"]]] = edge["bound"]
if len(edge) > 3:
self.costPoolToProduct[poolIdx[edge["pool"]]] \
[productIdx[edge["product"]]] = edge["cost"]
class _OptAgentRunner:
def __init__(self, time_limit):
self.model = OptModel()
self.time_limit = time_limit
def __enter__(self):
return self
def __exit__(self, exc_type, exc_value, traceback):
return False
def solve(self):
return solve(self.model, time_limit_s=float(self.time_limit))
def main(instance_file, output_file=None, time_limit=20):
data = PoolingInstance(instance_file)
with _OptAgentRunner(time_limit) as optimizer:
#
# Declare the optimization model
#
model = optimizer.model
# Zero-capacity arcs are constants, shared by every derived expression.
zero = model.constant(0.0)
def bounded_flow(upper):
return model.float(0, upper) if upper > 0 else zero
# Decision variables
# Flow from the components to the products
flowComponentToProduct = [[bounded_flow(data.upperBoundComponentToProduct[c][p])
for p in range(data.nbProducts)] for c in range(data.nbComponents)]
# Fraction of the total flow in pool o coming from the component c
fractionComponentToPool = [[bounded_flow(data.upperBoundFractionComponentToPool[c][o])
for o in range(data.nbPools)]
for c in range(data.nbComponents)]
# Flow from the pools to the products
flowPoolToProduct = [
[bounded_flow(data.upperBoundPoolToProduct[o][p])
for p in range(data.nbProducts)] for o in range(data.nbPools)]
# Flow from the components to the products and passing by the pools
flowComponentToProductByPool = [
[[(fractionComponentToPool[c][o] * flowPoolToProduct[o][p]
if data.upperBoundFractionComponentToPool[c][o] > 0
and data.upperBoundPoolToProduct[o][p] > 0 else zero)
for p in range(data.nbProducts)] for o in range(data.nbPools)]
for c in range(data.nbComponents)]
# Constraints
# Proportion
for o in range(data.nbPools):
proportion = model.sum(fractionComponentToPool[c][o]
for c in range(data.nbComponents))
model.constraint(proportion == 1)
# Component supply
for c in range(data.nbComponents):
flowToProducts = model.sum(flowComponentToProduct[c][p]
for p in range(data.nbProducts))
flowToPools = model.sum(flowComponentToProductByPool[c][o][p]
for p in range(data.nbProducts) for o in range(data.nbPools))
totalOutFlow = model.sum(flowToPools, flowToProducts)
model.constraint(totalOutFlow <= data.componentSupplies[c])
# Pool capacity (bounds on edges)
for c in range(data.nbComponents):
for o in range(data.nbPools):
flowComponentToPool = model.sum(flowComponentToProductByPool[c][o][p]
for p in range(data.nbProducts))
edgeCapacity = model.prod(data.poolCapacities[o],
fractionComponentToPool[c][o])
model.constraint(flowComponentToPool <= edgeCapacity)
# Product capacity
for p in range(data.nbProducts):
flowFromPools = model.sum(flowPoolToProduct[o][p] for o in range(data.nbPools))
flowFromComponents = model.sum(flowComponentToProduct[c][p]
for c in range(data.nbComponents))
totalInFlow = model.sum(flowFromComponents, flowFromPools)
model.constraint(totalInFlow <= data.productCapacities[p])
model.constraint(totalInFlow >= data.demand[p])
# Product tolerance
for p in range(data.nbProducts):
for k in range(data.nbAttributes):
# Attribute from the components
attributeFromComponents = model.sum(
data.componentQuality[c][k] * flowComponentToProduct[c][p]
for c in range(data.nbComponents)
if data.upperBoundComponentToProduct[c][p] > 0 and data.componentQuality[c][k] != 0)
# Attribute from the pools
attributeFromPools = model.sum(
data.componentQuality[c][k] * flowComponentToProductByPool[c][o][p]
for o in range(data.nbPools) for c in range(data.nbComponents)
if data.upperBoundFractionComponentToPool[c][o] > 0
and data.upperBoundPoolToProduct[o][p] > 0 and data.componentQuality[c][k] != 0)
# Total flow in the blending
totalFlowIn = model.sum(flowComponentToProduct[c][p]
for c in range(data.nbComponents)) \
+ model.sum(flowPoolToProduct[o][p] for o in range(data.nbPools))
totalAttributeIn = model.sum(attributeFromComponents, attributeFromPools)
model.constraint(totalAttributeIn >= data.minTolerance[p][k] * totalFlowIn)
model.constraint(totalAttributeIn <= data.maxTolerance[p][k] * totalFlowIn)
# Objective function
# Cost of the flows from the components directly to the products
directFlowCost = model.sum(
data.costComponentToProduct[c][p] * flowComponentToProduct[c][p]
for c in range(data.nbComponents) for p in range(data.nbProducts)
if data.upperBoundComponentToProduct[c][p] > 0 and data.costComponentToProduct[c][p] != 0)
# Cost of the flows from the components to the products passing by the pools
undirectFlowCost = model.sum(
(data.costComponentToPool[c][o] + data.costPoolToProduct[o][p]) *
flowComponentToProductByPool[c][o][p]
for c in range(data.nbComponents)
for o in range(data.nbPools) for p in range(data.nbProducts)
if data.upperBoundFractionComponentToPool[c][o] > 0 and data.upperBoundPoolToProduct[o][p] > 0
and data.costComponentToPool[c][o] + data.costPoolToProduct[o][p] != 0)
# Gain of selling the final products
productsGain = model.sum((model.sum(flowComponentToProduct[c][p]
for c in range(data.nbComponents))
+ model.sum(flowPoolToProduct[o][p]
for o in range(data.nbPools)))
* data.productPrices[p] for p in range(data.nbProducts))
# Cost of buying the components
componentsCost = model.sum(
(model.sum(flowComponentToProduct[c][p]
for p in range(data.nbProducts))
+ model.sum(flowComponentToProductByPool[c][o][p]
for p in range(data.nbProducts)
for o in range(data.nbPools)))
* data.componentPrices[c] for c in range(data.nbComponents))
profit = productsGain - componentsCost - (directFlowCost + undirectFlowCost)
# Maximize the total profit
model.maximize(profit)
# Solve with OptAgent
solution = optimizer.solve()
if not solution.feasible:
print(f"No feasible solution found; Status = {solution.feasible}")
return solution
#
# Write the solution
#
if output_file is not None:
with open(output_file, 'w', encoding="utf-8") as f:
component_to_poduct = []
component_to_pool_fraction = []
pool_to_product = []
# Solution flows from the components to the products
for c in range(data.nbComponents):
for p in range(data.nbProducts):
component_to_poduct.append(
{"component": data.componentNames[c],
"product": data.productNames[p],
"flow": flowComponentToProduct[c][p].value})
# Solution fraction of the inflow at pool o coming from the component c
for c in range(data.nbComponents):
for o in range(data.nbPools):
component_to_pool_fraction.append(
{"component": data.componentNames[c],
"pool": data.poolNames[o],
"flow": fractionComponentToPool[c][o].value})
# Solution flows from the pools to the products
for o in range(data.nbPools):
for p in range(data.nbProducts):
pool_to_product.append(
{"pool": data.poolNames[o],
"product": data.productNames[p],
"flow": flowPoolToProduct[o][p].value})
json.dump({"objective": profit.value, "solution":
{"component_to_pool_fraction": component_to_pool_fraction,
"component_to_product": component_to_poduct,
"pool_to_product": pool_to_product}}, f)
print(f"Profit = {profit.value}; Status = {solution.feasible}")
return solution
本地运行¶
在 notebook 所在目录执行以下 cell，即可调用一个混流实例。
INSTANCE_DIR = Path.cwd() / "instances"
print("Instances:", INSTANCE_DIR)
solution_pooling = main(INSTANCE_DIR / "haverly1.json", time_limit=1)