From ca054292032e4528f6a77a34f23ebc72e03cf00b Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Tue, 23 Sep 2025 11:39:39 -0500 Subject: [PATCH] uc: Add quadratic terms --- docs/guide/problems.ipynb | 116 ++++++++++++++++++++++---------------- miplearn/problems/uc.py | 7 +++ 2 files changed, 73 insertions(+), 50 deletions(-) diff --git a/docs/guide/problems.ipynb b/docs/guide/problems.ipynb index 10ef906..1fe751c 100644 --- a/docs/guide/problems.ipynb +++ b/docs/guide/problems.ipynb @@ -1221,7 +1221,6 @@ "id": "7048d771", "metadata": {}, "source": [ - "\n", "
\n", "Note\n", "\n", @@ -1230,7 +1229,7 @@ "\n", "### Formulation\n", "\n", - "Let $T$ be the number of time steps, $G$ be the number of generation units, and let $D_t$ be the power demand (in MW) at time $t$. For each generating unit $g$, let $P^\\max_g$ and $P^\\min_g$ be the maximum and minimum amount of power the unit is able to produce when switched on; let $L_g$ and $l_g$ be the minimum up- and down-time for unit $g$; let $C^\\text{fixed}$ be the cost to keep unit $g$ on for one time step, regardless of its power output level; let $C^\\text{start}$ be the cost to switch unit $g$ on; and let $C^\\text{var}$ be the cost for generator $g$ to produce 1 MW of power. In this formulation, we assume linear production costs. For each generator $g$ and time $t$, let $x_{gt}$ be a binary variable which equals one if unit $g$ is on at time $t$, let $w_{gt}$ be a binary variable which equals one if unit $g$ switches from being off at time $t-1$ to being on at time $t$, and let $p_{gt}$ be a continuous variable which indicates the amount of power generated. The formulation is given by:" + "Let $T$ be the number of time steps, $G$ be the number of generation units, and let $D_t$ be the power demand (in MW) at time $t$. For each generating unit $g$, let $P^\\max_g$ and $P^\\min_g$ be the maximum and minimum amount of power the unit is able to produce when switched on; let $L_g$ and $l_g$ be the minimum up- and down-time for unit $g$; let $C^\\text{fixed}$ be the cost to keep unit $g$ on for one time step, regardless of its power output level; let $C^\\text{start}$ be the cost to switch unit $g$ on; let $C^\\text{prod-lin}$ be the linear cost coefficient for generator $g$ to produce 1 MW of power; and let $C^\\text{prod-quad}$ be the quadratic cost coefficient for unit $g$. For each generator $g$ and time $t$, let $x_{gt}$ be a binary variable which equals one if unit $g$ is on at time $t$, let $w_{gt}$ be a binary variable which equals one if unit $g$ switches from being off at time $t-1$ to being on at time $t$, and let $p_{gt}$ be a continuous variable which indicates the amount of power generated. The formulation is given by:" ] }, { @@ -1238,14 +1237,14 @@ "id": "bec5ee1c", "metadata": {}, "source": [ - "\n", "$$\n", "\\begin{align*}\n", "\\text{minimize} \\;\\;\\;\n", " & \\sum_{t=1}^T \\sum_{g=1}^G \\left(\n", " x_{gt} C^\\text{fixed}_g\n", " + w_{gt} C^\\text{start}_g\n", - " + p_{gt} C^\\text{var}_g\n", + " + p_{gt} C^\\text{prod-lin}_g\n", + " + p_{gt}^2 C^\\text{prod-quad}_g\n", " \\right)\n", " \\\\\n", "\\text{such that} \\;\\;\\;\n", @@ -1287,12 +1286,11 @@ "id": "01bed9fc", "metadata": {}, "source": [ - "\n", "### Random instance generator\n", "\n", "The class `UnitCommitmentGenerator` can be used to generate random instances of this problem.\n", "\n", - "First, the user-provided probability distributions `n_units` and `n_periods` are sampled to determine the number of generating units and the number of time steps, respectively. Then, for each unit, the probabilities `max_power` and `min_power` are sampled to determine the unit's maximum and minimum power output. To make it easier to generate valid ranges, `min_power` is not specified as the absolute power level in MW, but rather as a multiplier of `max_power`; for example, if `max_power` samples to 100 and `min_power` samples to 0.5, then the unit's power range is set to `[50,100]`. Then, the distributions `cost_startup`, `cost_prod` and `cost_fixed` are sampled to determine the unit's startup, variable and fixed costs, while the distributions `min_uptime` and `min_downtime` are sampled to determine its minimum up/down-time.\n", + "First, the user-provided probability distributions `n_units` and `n_periods` are sampled to determine the number of generating units and the number of time steps, respectively. Then, for each unit, the probabilities `max_power` and `min_power` are sampled to determine the unit's maximum and minimum power output. To make it easier to generate valid ranges, `min_power` is not specified as the absolute power level in MW, but rather as a multiplier of `max_power`; for example, if `max_power` samples to 100 and `min_power` samples to 0.5, then the unit's power range is set to `[50,100]`. Then, the distributions `cost_startup`, `cost_prod`, `cost_prod_quad` and `cost_fixed` are sampled to determine the unit's startup, linear variable, quadratic variable, and fixed costs, while the distributions `min_uptime` and `min_downtime` are sampled to determine its minimum up/down-time.\n", "\n", "After parameters for the units have been generated, the class then generates a periodic demand curve, with a peak every 12 time steps, in the range $(0.4C, 0.8C)$, where $C$ is the sum of all units' maximum power output. Finally, all costs and demand values are perturbed by random scaling factors independently sampled from the distributions `cost_jitter` and `demand_jitter`, respectively.\n", "\n", @@ -1309,7 +1307,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 2, "id": "6217da7c", "metadata": { "ExecuteTime": { @@ -1333,78 +1331,94 @@ "min_power[0] [117.79 245.85 271.85 207.7 81.38]\n", "cost_startup[0] [3042.42 5247.56 4319.45 2912.29 6118.53]\n", "cost_prod[0] [ 6.97 14.61 18.32 22.8 39.26]\n", - "cost_fixed[0] [199.67 514.23 592.41 46.45 607.54]\n", + "cost_prod_quad[0] [0.02 0.0514 0.0592 0.0046 0.0608]\n", + "cost_fixed[0] [170.52 65.05 948.89 965.63 808.4 ]\n", "demand[0]\n", - " [ 905.06 915.41 1166.52 1212.29 1127.81 953.52 905.06 796.21 783.78\n", - " 866.23 768.62 899.59 905.06 946.23 1087.61 1004.24 1048.36 992.03\n", - " 905.06 750.82 691.48 606.15 658.5 809.95]\n", + " [ 869.31 897.58 1212.29 1124.08 930.48 1012.62 869.31 606.15]\n", "\n", "min_power[1] [117.79 245.85 271.85 207.7 81.38]\n", "max_power[1] [218.54 477.82 379.4 319.4 120.21]\n", "min_uptime[1] [7 6 3 5 7]\n", "min_downtime[1] [7 3 5 6 2]\n", "min_power[1] [117.79 245.85 271.85 207.7 81.38]\n", - "cost_startup[1] [2458.08 6200.26 4585.74 2666.05 4783.34]\n", - "cost_prod[1] [ 6.31 13.33 20.42 24.37 46.86]\n", - "cost_fixed[1] [196.9 416.42 655.57 52.51 626.15]\n", + "cost_startup[1] [3710.99 6283.5 4530.89 3526.6 4859.62]\n", + "cost_prod[1] [ 5.91 11.29 16.72 21.53 34.77]\n", + "cost_prod_quad[1] [0.0233 0.0477 0.0527 0.0047 0.0499]\n", + "cost_fixed[1] [ 196.29 51.21 1179.89 1097.07 686.62]\n", "demand[1]\n", - " [ 981.42 840.07 1095.59 1102.03 1088.41 932.29 863.67 848.56 761.33\n", - " 828.28 775.18 834.99 959.76 865.72 1193.52 1058.92 985.19 893.92\n", - " 962.16 781.88 723.15 639.04 602.4 787.02]\n", + " [ 827.37 926.76 1166.64 1128.59 939.17 948.8 950.95 639.5 ]\n", "\n", "Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (linux64 - \"Ubuntu 24.04.3 LTS\")\n", "\n", "CPU model: AMD Ryzen 9 3950X 16-Core Processor, instruction set [SSE2|AVX|AVX2]\n", "Thread count: 16 physical cores, 32 logical processors, using up to 32 threads\n", "\n", - "Optimize a model with 578 rows, 360 columns and 2128 nonzeros\n", - "Model fingerprint: 0x4dc1c661\n", - "Variable types: 120 continuous, 240 integer (240 binary)\n", + "Optimize a model with 162 rows, 120 columns and 512 nonzeros\n", + "Model fingerprint: 0x1e3651da\n", + "Model has 40 quadratic objective terms\n", + "Variable types: 40 continuous, 80 integer (80 binary)\n", "Coefficient statistics:\n", " Matrix range [1e+00, 5e+02]\n", " Objective range [7e+00, 6e+03]\n", + " QObjective range [9e-03, 1e-01]\n", " Bounds range [1e+00, 1e+00]\n", " RHS range [1e+00, 1e+03]\n", - "Presolve removed 341 rows and 133 columns\n", - "Presolve time: 0.01s\n", - "Presolved: 237 rows, 227 columns, 725 nonzeros\n", - "Variable types: 114 continuous, 113 integer (113 binary)\n", - "Found heuristic solution: objective 475243.89360\n", + "Found heuristic solution: objective 282371.35206\n", + "Presolve removed 61 rows and 40 columns\n", + "Presolve time: 0.00s\n", + "Presolved: 137 rows, 98 columns, 362 nonzeros\n", + "Presolved model has 40 quadratic objective terms\n", + "Variable types: 58 continuous, 40 integer (40 binary)\n", "\n", - "Root relaxation: objective 3.361348e+05, 96 iterations, 0.00 seconds (0.00 work units)\n", + "Root relaxation: objective 1.995341e+05, 126 iterations, 0.00 seconds (0.00 work units)\n", "\n", " Nodes | Current Node | Objective Bounds | Work\n", " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", "\n", - " 0 0 336134.820 0 18 475243.894 336134.820 29.3% - 0s\n", - "H 0 0 471441.37480 336134.820 28.7% - 0s\n", - "H 0 0 410679.27820 336134.820 18.2% - 0s\n", - "H 0 0 391706.31610 336134.820 14.2% - 0s\n", - "H 0 0 374515.31390 336134.820 10.2% - 0s\n", - "H 0 0 369369.87450 336134.820 9.00% - 0s\n", - "H 0 0 368600.14450 344055.048 6.66% - 0s\n", - "H 0 0 368180.65796 364657.488 0.96% - 0s\n", - "H 0 0 364721.76610 364657.488 0.02% - 0s\n", - " 0 0 364721.766 0 6 364721.766 364721.766 0.00% - 0s\n", + " 0 0 199534.058 0 8 282371.352 199534.058 29.3% - 0s\n", + "H 0 0 212978.80635 199534.058 6.31% - 0s\n", + "H 0 0 208079.69920 199534.058 4.11% - 0s\n", + " 0 0 203097.772 0 16 208079.699 203097.772 2.39% - 0s\n", + " 0 0 203097.772 0 4 208079.699 203097.772 2.39% - 0s\n", + " 0 0 203097.772 0 4 208079.699 203097.772 2.39% - 0s\n", + " 0 0 203097.772 0 3 208079.699 203097.772 2.39% - 0s\n", + " 0 0 203097.772 0 3 208079.699 203097.772 2.39% - 0s\n", + " 0 0 203097.772 0 3 208079.699 203097.772 2.39% - 0s\n", + " 0 0 205275.299 0 - 208079.699 205275.299 1.35% - 0s\n", + " 0 0 205777.846 0 2 208079.699 205777.846 1.11% - 0s\n", + " 0 0 205789.407 0 - 208079.699 205789.407 1.10% - 0s\n", + " 0 0 postponed 0 208079.699 205789.407 1.10% - 0s\n", + " 0 0 postponed 0 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 4 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 4 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 1 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 - 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 - 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 - 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 - 208079.699 205789.408 1.10% - 0s\n", + " 0 0 205789.408 0 2 208079.699 205789.408 1.10% - 0s\n", + "H 0 0 207525.63560 205789.408 0.84% - 0s\n", + " 0 0 205789.408 0 2 207525.636 205789.408 0.84% - 0s\n", + " 0 2 205789.408 0 2 207525.636 205789.408 0.84% - 0s\n", + "* 9 0 5 205789.40812 205789.408 0.00% 0.0 0s\n", "\n", "Cutting planes:\n", - " Gomory: 2\n", - " Cover: 7\n", - " Implied bound: 1\n", - " Clique: 19\n", - " MIR: 3\n", - " RLT: 10\n", - " Relax-and-lift: 1\n", - "\n", - "Explored 1 nodes (181 simplex iterations) in 0.03 seconds (0.01 work units)\n", + " Cover: 1\n", + " Implied bound: 6\n", + " MIR: 2\n", + " Flow cover: 2\n", + " RLT: 2\n", + "\n", + "Explored 11 nodes (621 simplex iterations) in 0.32 seconds (0.19 work units)\n", "Thread count was 32 (of 32 available processors)\n", "\n", - "Solution count 7: 364722 368600 374515 ... 475244\n", + "Solution count 5: 205789 207526 208080 ... 282371\n", + "No other solutions better than 205789\n", "\n", "Optimal solution found (tolerance 1.00e-04)\n", - "Best objective 3.647217661000e+05, best bound 3.647217661000e+05, gap 0.0000%\n", + "Best objective 2.057894080889e+05, best bound 2.057894081187e+05, gap 0.0000%\n", "\n", - "User-callback calls 816, time in user-callback 0.00 sec\n" + "User-callback calls 647, time in user-callback 0.00 sec\n" ] } ], @@ -1418,14 +1432,15 @@ "random.seed(42)\n", "np.random.seed(42)\n", "\n", - "# Generate a random instance with 5 generators and 24 time steps\n", + "# Generate a random instance with 5 generators and 8 time steps\n", "data = UnitCommitmentGenerator(\n", " n_units=randint(low=5, high=6),\n", - " n_periods=randint(low=24, high=25),\n", + " n_periods=randint(low=8, high=9),\n", " max_power=uniform(loc=50, scale=450),\n", " min_power=uniform(loc=0.5, scale=0.25),\n", " cost_startup=uniform(loc=0, scale=10_000),\n", " cost_prod=uniform(loc=0, scale=50),\n", + " cost_prod_quad=uniform(loc=0, scale=0.1),\n", " cost_fixed=uniform(loc=0, scale=1_000),\n", " min_uptime=randint(low=2, high=8),\n", " min_downtime=randint(low=2, high=8),\n", @@ -1443,6 +1458,7 @@ " print(f\"min_power[{i}]\", data[i].min_power)\n", " print(f\"cost_startup[{i}]\", data[i].cost_startup)\n", " print(f\"cost_prod[{i}]\", data[i].cost_prod)\n", + " print(f\"cost_prod_quad[{i}]\", data[i].cost_prod_quad)\n", " print(f\"cost_fixed[{i}]\", data[i].cost_fixed)\n", " print(f\"demand[{i}]\\n\", data[i].demand)\n", " print()\n", diff --git a/miplearn/problems/uc.py b/miplearn/problems/uc.py index 63f4178..3731e36 100644 --- a/miplearn/problems/uc.py +++ b/miplearn/problems/uc.py @@ -25,6 +25,7 @@ class UnitCommitmentData: min_downtime: np.ndarray cost_startup: np.ndarray cost_prod: np.ndarray + cost_prod_quad: np.ndarray cost_fixed: np.ndarray @@ -37,6 +38,7 @@ class UnitCommitmentGenerator: min_power: rv_frozen = uniform(loc=0.5, scale=0.25), cost_startup: rv_frozen = uniform(loc=0, scale=10_000), cost_prod: rv_frozen = uniform(loc=0, scale=50), + cost_prod_quad: rv_frozen = uniform(loc=0, scale=0), cost_fixed: rv_frozen = uniform(loc=0, scale=1_000), min_uptime: rv_frozen = randint(low=2, high=8), min_downtime: rv_frozen = randint(low=2, high=8), @@ -50,6 +52,7 @@ class UnitCommitmentGenerator: self.min_power = min_power self.cost_startup = cost_startup self.cost_prod = cost_prod + self.cost_prod_quad = cost_prod_quad self.cost_fixed = cost_fixed self.min_uptime = min_uptime self.min_downtime = min_downtime @@ -72,6 +75,7 @@ class UnitCommitmentGenerator: min_downtime = self.min_downtime.rvs(G) cost_startup = self.cost_startup.rvs(G) cost_prod = self.cost_prod.rvs(G) + cost_prod_quad = self.cost_prod_quad.rvs(G) cost_fixed = self.cost_fixed.rvs(G) capacity = max_power.sum() @@ -91,6 +95,7 @@ class UnitCommitmentGenerator: min_downtime = self.ref_data.min_downtime cost_startup = self.ref_data.cost_startup * self.cost_jitter.rvs(G) cost_prod = self.ref_data.cost_prod * self.cost_jitter.rvs(G) + cost_prod_quad = self.ref_data.cost_prod_quad * self.cost_jitter.rvs(G) cost_fixed = self.ref_data.cost_fixed * self.cost_jitter.rvs(G) data = UnitCommitmentData( @@ -101,6 +106,7 @@ class UnitCommitmentGenerator: min_downtime, cost_startup.round(2), cost_prod.round(2), + cost_prod_quad.round(4), cost_fixed.round(2), ) @@ -143,6 +149,7 @@ def build_uc_model_gurobipy(data: Union[str, UnitCommitmentData]) -> GurobiModel is_on[g, t] * data.cost_fixed[g] + switch_on[g, t] * data.cost_startup[g] + prod[g, t] * data.cost_prod[g] + + prod[g, t] * prod[g, t] * data.cost_prod_quad[g] for g in range(G) for t in range(T) )