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