问题描述¶
在 Revenue Management Problem 中，一家公司希望在一个分为若干时段的时间范围内，通过销售某种产品来最大化其收入。它必须决定在时间范围开始时购买的产品总量。然后，在每个时段，它必须决定该时段内销售的产品数量。整个时间范围内销售的产品总数不得超过最初购买的数量。产品价格随时间上涨。为了获得更多利润，公司应当为后续客户保留一些产品，而不是在早期全部售出。为了在每个时段做出最明智的决策，它必须考虑后续时段的需求。由于需求是随机的，公司运行大量仿真以获得对给定单位分配下收入的稳健估计。
学习要点¶
- 使用整数变量表示采购量与各时段的预留量。
- 使用
create_double_external_function将 Python 仿真函数接入目标表达式。 - 使用
add_evaluation_point提供已知评价点，并用external_evaluation_limit限制新评价次数。
数据¶
时间范围由 3 个时段组成，产品的初始成本为 $80。
每个时段 t 的需求由方程 Dₜ=μₜXYₜ 定义，其中：
- Yₜ 服从速率参数 λ=1 的指数分布。
- X 服从形状参数 k=1、尺度参数 θ=1 的 Gamma 分布，等价于标准指数分布。
- μₜ 是该时段的平均需求。
每个时段的价格和平均需求见下表： **时段12****3**价格100300400平均需求502030每个时段的价格和平均需求 为了获得收入的稳健估计，仿真需要使用 Monte Carlo method 运行大量次数（1,000,000 次）。因此每次仿真需要数秒钟。由于该评估函数开销极大，我们不能负担大量运行次数，必须谨慎选择每个评估点。
建模思路¶
OptAgent 模型定义三个取值在 0 到 100 之间的整数变量，分别表示初始采购量、为后两个时段预留的数量、为最后一个时段预留的数量。后一个变量不得大于前一个变量。
Python 仿真函数根据这三个变量计算各时段的销售收入，减去采购成本，返回平均净收益。模型通过 create_double_external_function 注册回调，并最大化调用表达式的值。仿真固定随机种子，以便比较候选决策。
代码使用 add_evaluation_point 注册一个已知评价点。main 的 evaluation_limit 参数传给 solve 的 external_evaluation_limit，time_limit 则传给 time_limit_s。单次仿真本身有计算开销；运行时应结合这两个预算检查终止状态。
Python 实现¶
# Copyright (c) Hexaly. Permission is hereby granted to use, copy,
# and modify this code for applications developed with Hexaly.
import math
import random
from pathlib import Path
from optagent import OptModel, solve
class RevenueManagementFunction:
def __init__(self, seed):
self.nb_periods = 3
self.prices = [100, 300, 400]
self.mean_demands = [50, 20, 30]
self.purchase_price = 80
self.evaluated_points = [{
"point": [100, 50, 30],
"value": 4740.99
}]
self.nb_simulations = int(1e6)
self.seed = seed
# External function
def evaluate(self, argument_values):
variables = [argument_values.get(i) for i in range(argument_values.count())]
# Initial quantity purchased
nb_units_purchased = variables[0]
# Number of units that should be left for future periods
nb_units_reserved = variables[1:] + [0]
# Set seed for reproducibility
random.seed(self.seed)
# Create distribution
X = [gamma_sample() for i in range(self.nb_simulations)]
Y = [[exponential_sample() for i in range(self.nb_periods)]
for j in range(self.nb_simulations)]
# Run simulations
sum_profit = 0.0
for i in range(self.nb_simulations):
remaining_capacity = nb_units_purchased
for j in range(self.nb_periods):
# Generate demand for period j
demand_j = int(self.mean_demands[j] * X[i] * Y[i][j])
nb_units_sold = min(
max(remaining_capacity - nb_units_reserved[j], 0),
demand_j)
remaining_capacity = remaining_capacity - nb_units_sold
sum_profit += self.prices[j] * nb_units_sold
# Calculate mean revenue
mean_profit = sum_profit / self.nb_simulations
mean_revenue = mean_profit - self.purchase_price * nb_units_purchased
return mean_revenue
def exponential_sample(rate_param=1.0):
u = random.random()
return math.log(1 - u) / (-rate_param)
def gamma_sample(scale_param=1.0):
return exponential_sample(scale_param)
def main(output_file=None, time_limit=30, evaluation_limit=None):
model = OptModel()
# Generate data
revenue_management = RevenueManagementFunction(1)
nb_periods = revenue_management.nb_periods
# Declare decision variables
variables = [model.int(0, 100) for i in range(nb_periods)]
# Create the function
func_expr = model.create_double_external_function(revenue_management.evaluate)
for evaluation in revenue_management.evaluated_points:
func_expr.add_evaluation_point(evaluation["point"], evaluation["value"])
func_call = func_expr(*variables)
# Declare constraints
for i in range(1, nb_periods):
model.constraint(variables[i] <= variables[i - 1])
model.maximize(func_call)
try:
solution = solve(
model,
time_limit_s=float(time_limit),
external_evaluation_limit=int(evaluation_limit) if evaluation_limit is not None else None,
)
except RuntimeError as error:
if "evaluation cancelled" not in str(error).lower():
raise
print(f"External simulation cancelled by time limit: {error}")
return None
if not solution.feasible:
print(f"No feasible solution found; Status = {solution.feasible}")
return solution
print(f"Mean revenue = {func_call.value}; Status = {solution.feasible}")
# Write the solution in a file
if output_file is not None:
lines = [f"obj={func_call.value}", f"b={variables[0].value}"]
lines.extend(f"r{i + 1}={variables[i].value}" for i in range(1, nb_periods))
Path(output_file).write_text("\n".join(lines) + "\n", encoding="utf-8")
return solution
实例调用¶
该示例不依赖外部实例文件，直接运行仿真收益模型。
solution_revenue = main(time_limit=1)
getattr(solution_revenue, "feasible", None)