From 146fb6b615c5ecae97f119644edf46ef40c4853d Mon Sep 17 00:00:00 2001 From: "Alinson S. Xavier" Date: Mon, 8 Dec 2025 15:31:22 -0600 Subject: [PATCH] Implement UnitCommitmentPerturber --- docs/guide/problems.ipynb | 121 +++++++++++++++++------------- miplearn/problems/uc.py | 152 +++++++++++++++++++++++++++----------- tests/problems/test_uc.py | 3 - 3 files changed, 179 insertions(+), 97 deletions(-) diff --git a/docs/guide/problems.ipynb b/docs/guide/problems.ipynb index 6b6f188..a52a7ab 100644 --- a/docs/guide/problems.ipynb +++ b/docs/guide/problems.ipynb @@ -1357,9 +1357,9 @@ "\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", + "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.\n", "\n", - "If `fix_units=True`, then the list of generators (with their respective parameters) is kept the same for all generated instances. If `cost_jitter` and `demand_jitter` are provided, the instances will still have slightly different costs and demands." + "To create multiple instances with the same unit structure but different values, you can use `UnitCommitmentPerturber`. This class takes an existing UnitCommitmentData instance and generates new instances by applying randomization factors to the existing costs and demand while keeping the unit characteristics (power limits, minimum up/down-times) fixed. More specifically, the costs are multiplied by random scaling factors sampled from the distribution `cost_jitter`, and the demand is multiplied by random scaling factors sampled from `demand_jitter`." ] }, { @@ -1380,7 +1380,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 12, "id": "6217da7c", "metadata": { "ExecuteTime": { @@ -1401,25 +1401,23 @@ "max_power[0] [218.54 477.82 379.4 319.4 120.21]\n", "min_uptime[0] [7 6 3 5 7]\n", "min_downtime[0] [7 3 5 6 2]\n", - "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_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", + "cost_startup[0] [3710.99 6283.5 4530.89 3526.6 4859.62]\n", + "cost_prod[0] [ 5.91 11.29 16.72 21.53 34.77]\n", + "cost_prod_quad[0] [0.0233 0.0477 0.0527 0.0047 0.0499]\n", + "cost_fixed[0] [ 196.29 51.21 1179.89 1097.07 686.62]\n", "demand[0]\n", - " [ 869.31 897.58 1212.29 1124.08 930.48 1012.62 869.31 606.15]\n", + " [ 827.37 926.76 1166.64 1128.59 939.17 948.8 950.95 639.5 ]\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] [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", + "cost_startup[1] [3594.78 5571.07 3954.24 2276.77 5540.27]\n", + "cost_prod[1] [ 6.36 16.29 19.58 27.21 38.71]\n", + "cost_prod_quad[1] [0.0162 0.0569 0.0669 0.0047 0.069 ]\n", + "cost_fixed[1] [169.99 65.79 914.51 736.5 649.91]\n", "demand[1]\n", - " [ 827.37 926.76 1166.64 1128.59 939.17 948.8 950.95 639.5 ]\n", + " [ 783.34 954.21 1262.44 1175.56 980.96 926.35 844.7 559.58]\n", "\n", "Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (linux64 - \"Ubuntu 22.04.5 LTS\")\n", "\n", @@ -1427,60 +1425,77 @@ "Thread count: 16 physical cores, 16 logical processors, using up to 16 threads\n", "\n", "Optimize a model with 162 rows, 120 columns and 512 nonzeros (Min)\n", - "Model fingerprint: 0x1e3651da\n", + "Model fingerprint: 0xf8c1c4e3\n", "Model has 120 linear objective coefficients\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", + " Objective range [6e+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", - "Found heuristic solution: objective 282371.35206\n", - "Presolve removed 61 rows and 40 columns\n", + "Found heuristic solution: objective 277606.67683\n", + "Presolve removed 45 rows and 32 columns\n", "Presolve time: 0.00s\n", - "Presolved: 137 rows, 98 columns, 362 nonzeros\n", + "Presolved: 147 rows, 103 columns, 416 nonzeros\n", "Presolved model has 40 quadratic objective terms\n", - "Variable types: 58 continuous, 40 integer (40 binary)\n", + "Variable types: 55 continuous, 48 integer (48 binary)\n", "\n", - "Root relaxation: objective 1.995341e+05, 126 iterations, 0.00 seconds (0.00 work units)\n", + "Root relaxation: objective 1.823746e+05, 145 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 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 2 208079.699 203097.772 2.39% - 0s\n", - " 0 0 204065.874 0 2 208079.699 204065.874 1.93% - 0s\n", - " 0 0 205410.033 0 2 208079.699 205410.033 1.28% - 0s\n", - " 0 0 205413.331 0 2 208079.699 205413.331 1.28% - 0s\n", - "H 0 0 205789.40802 205458.260 0.16% - 0s\n", - "* 0 0 0 205789.40802 205458.260 0.16% - 0s\n", - " 0 0 - 0 205789.408 205787.978 0.00% - 0s\n", + " 0 0 182374.559 0 9 277606.677 182374.559 34.3% - 0s\n", + "H 0 0 190641.39953 182374.559 4.34% - 0s\n", + "H 0 0 190346.53552 182374.559 4.19% - 0s\n", + " 0 0 187223.908 0 24 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 5 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 5 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 7 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 - 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 2 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 3 190346.536 187223.908 1.64% - 0s\n", + " 0 0 187223.908 0 - 190346.536 187223.908 1.64% - 0s\n", + " 0 0 188108.163 0 3 190346.536 188108.163 1.18% - 0s\n", + " 0 0 189470.596 0 3 190346.536 189470.596 0.46% - 0s\n", + " 0 0 189884.665 0 3 190346.536 189884.665 0.24% - 0s\n", + " 0 0 190074.608 0 3 190346.536 190074.608 0.14% - 0s\n", + " 0 0 190074.608 0 3 190346.536 190074.608 0.14% - 0s\n", + " 0 0 190088.760 0 3 190346.536 190088.760 0.14% - 0s\n", + " 0 0 190134.867 0 3 190346.536 190134.867 0.11% - 0s\n", + " 0 0 190141.852 0 3 190346.536 190141.852 0.11% - 0s\n", + " 0 0 190141.870 0 3 190346.536 190141.870 0.11% - 0s\n", + " 0 0 190144.673 0 3 190346.536 190144.673 0.11% - 0s\n", + " 0 0 190148.372 0 3 190346.536 190148.372 0.10% - 0s\n", + " 0 0 190155.514 0 3 190346.536 190155.514 0.10% - 0s\n", + " 0 0 190157.727 0 3 190346.536 190157.727 0.10% - 0s\n", + " 0 0 190158.746 0 3 190346.536 190158.746 0.10% - 0s\n", + " 0 0 190159.827 0 3 190346.536 190159.827 0.10% - 0s\n", + " 0 0 190160.144 0 3 190346.536 190160.144 0.10% - 0s\n", + " 0 0 190183.316 0 3 190346.536 190183.316 0.09% - 0s\n", + " 0 0 190190.168 0 3 190346.536 190190.168 0.08% - 0s\n", + " 0 0 190270.542 0 3 190346.536 190270.542 0.04% - 0s\n", + " 0 0 190273.871 0 3 190346.536 190273.871 0.04% - 0s\n", + " 0 0 190317.529 0 3 190346.536 190317.529 0.02% - 0s\n", + " 0 0 190317.932 0 3 190346.536 190317.932 0.02% - 0s\n", + " 0 0 - 0 190346.536 190337.053 0.00% - 0s\n", "\n", "Cutting planes:\n", - " Gomory: 1\n", - " Cover: 1\n", - " Implied bound: 1\n", - " MIR: 7\n", - " Flow cover: 3\n", + " Implied bound: 3\n", + " MIR: 5\n", " Relax-and-lift: 1\n", "\n", - "Explored 1 nodes (431 simplex iterations) in 0.02 seconds (0.01 work units)\n", + "Explored 1 nodes (1077 simplex iterations) in 0.05 seconds (0.02 work units)\n", "Thread count was 16 (of 16 available processors)\n", "\n", - "Solution count 4: 205789 208080 212979 282371 \n", + "Solution count 3: 190347 190641 277607 \n", "\n", "Optimal solution found (tolerance 1.00e-04)\n", - "Best objective 2.057894080182e+05, best bound 2.057879779509e+05, gap 0.0007%\n", + "Best objective 1.903465355189e+05, best bound 1.903370533334e+05, gap 0.0050%\n", "\n", - "User-callback calls 540, time in user-callback 0.00 sec\n" + "User-callback calls 654, time in user-callback 0.00 sec\n" ] } ], @@ -1488,14 +1503,14 @@ "import random\n", "import numpy as np\n", "from scipy.stats import uniform, randint\n", - "from miplearn.problems.uc import UnitCommitmentGenerator, build_uc_model_gurobipy\n", + "from miplearn.problems.uc import UnitCommitmentGenerator, UnitCommitmentPerturber, build_uc_model_gurobipy\n", "\n", "# Set random seed to make example reproducible\n", "random.seed(42)\n", "np.random.seed(42)\n", "\n", - "# Generate a random instance with 5 generators and 8 time steps\n", - "data = UnitCommitmentGenerator(\n", + "# Generate a reference instance with 5 generators and 8 time steps\n", + "generator = UnitCommitmentGenerator(\n", " n_units=randint(low=5, high=6),\n", " n_periods=randint(low=8, high=9),\n", " max_power=uniform(loc=50, scale=450),\n", @@ -1506,10 +1521,15 @@ " 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", + ")\n", + "reference_instance = generator.generate(1)[0]\n", + "\n", + "# Generate perturbed instances using the reference\n", + "perturber = UnitCommitmentPerturber(\n", " cost_jitter=uniform(loc=0.75, scale=0.5),\n", " demand_jitter=uniform(loc=0.9, scale=0.2),\n", - " fix_units=True,\n", - ").generate(10)\n", + ")\n", + "data = perturber.perturb(reference_instance, 10)\n", "\n", "# Print problem data for the two first instances\n", "for i in range(2):\n", @@ -1517,7 +1537,6 @@ " print(f\"max_power[{i}]\", data[i].max_power)\n", " print(f\"min_uptime[{i}]\", data[i].min_uptime)\n", " print(f\"min_downtime[{i}]\", data[i].min_downtime)\n", - " 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", diff --git a/miplearn/problems/uc.py b/miplearn/problems/uc.py index 3731e36..262b9d2 100644 --- a/miplearn/problems/uc.py +++ b/miplearn/problems/uc.py @@ -30,6 +30,12 @@ class UnitCommitmentData: class UnitCommitmentGenerator: + """Random instance generator for the Unit Commitment Problem. + + Generates instances by creating new random unit commitment problems with + parameters sampled from user-provided probability distributions. + """ + def __init__( self, n_units: rv_frozen = randint(low=1_000, high=1_001), @@ -42,10 +48,32 @@ class UnitCommitmentGenerator: 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), - cost_jitter: rv_frozen = uniform(loc=0.75, scale=0.5), - demand_jitter: rv_frozen = uniform(loc=0.9, scale=0.2), - fix_units: bool = False, ) -> None: + """Initialize the problem generator. + + Parameters + ---------- + n_units: rv_frozen + Probability distribution for number of units. + n_periods: rv_frozen + Probability distribution for number of periods. + max_power: rv_frozen + Probability distribution for maximum power output. + min_power: rv_frozen + Probability distribution for minimum power output (as fraction of max_power). + cost_startup: rv_frozen + Probability distribution for startup costs. + cost_prod: rv_frozen + Probability distribution for production costs. + cost_prod_quad: rv_frozen + Probability distribution for quadratic production costs. + cost_fixed: rv_frozen + Probability distribution for fixed costs. + min_uptime: rv_frozen + Probability distribution for minimum uptime. + min_downtime: rv_frozen + Probability distribution for minimum downtime. + """ self.n_units = n_units self.n_periods = n_periods self.max_power = max_power @@ -56,49 +84,33 @@ class UnitCommitmentGenerator: self.cost_fixed = cost_fixed self.min_uptime = min_uptime self.min_downtime = min_downtime - self.cost_jitter = cost_jitter - self.demand_jitter = demand_jitter - self.fix_units = fix_units - self.ref_data: Optional[UnitCommitmentData] = None def generate(self, n_samples: int) -> List[UnitCommitmentData]: def _sample() -> UnitCommitmentData: - if self.ref_data is None: - T = self.n_periods.rvs() - G = self.n_units.rvs() + T = self.n_periods.rvs() + G = self.n_units.rvs() - # Generate unit parameteres - max_power = self.max_power.rvs(G) - min_power = max_power * self.min_power.rvs(G) - max_power = max_power - min_uptime = self.min_uptime.rvs(G) - 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() + # Generate unit parameteres + max_power = self.max_power.rvs(G) + min_power = max_power * self.min_power.rvs(G) + max_power = max_power + min_uptime = self.min_uptime.rvs(G) + 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() - # Generate periodic demand in the range [0.4, 0.8] * capacity, with a peak every 12 hours. - demand = np.sin([i / 6 * pi for i in range(T)]) - demand *= uniform(loc=0, scale=1).rvs(T) - demand -= demand.min() - demand /= demand.max() / 0.4 - demand += 0.4 - demand *= capacity - else: - T, G = len(self.ref_data.demand), len(self.ref_data.max_power) - demand = self.ref_data.demand * self.demand_jitter.rvs(T) - min_power = self.ref_data.min_power - max_power = self.ref_data.max_power - min_uptime = self.ref_data.min_uptime - 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) + # Generate periodic demand in the range [0.4, 0.8] * capacity, with a peak every 12 hours. + demand = np.sin([i / 6 * pi for i in range(T)]) + demand *= uniform(loc=0, scale=1).rvs(T) + demand -= demand.min() + demand /= demand.max() / 0.4 + demand += 0.4 + demand *= capacity - data = UnitCommitmentData( + return UnitCommitmentData( demand.round(2), min_power.round(2), max_power.round(2), @@ -110,10 +122,64 @@ class UnitCommitmentGenerator: cost_fixed.round(2), ) - if self.ref_data is None and self.fix_units: - self.ref_data = data + return [_sample() for _ in range(n_samples)] - return data + +class UnitCommitmentPerturber: + """Perturbation generator for existing Unit Commitment instances. + + Takes an existing UnitCommitmentData instance and generates new instances + by applying randomization factors to the existing costs and demand while + keeping the unit structure fixed. + """ + + def __init__( + self, + cost_jitter: rv_frozen = uniform(loc=0.75, scale=0.5), + demand_jitter: rv_frozen = uniform(loc=0.9, scale=0.2), + ) -> None: + """Initialize the perturbation generator. + + Parameters + ---------- + cost_jitter: rv_frozen + Probability distribution for randomization factors applied to costs. + demand_jitter: rv_frozen + Probability distribution for randomization factors applied to demand. + """ + assert isinstance( + cost_jitter, rv_frozen + ), "cost_jitter should be a SciPy probability distribution" + assert isinstance( + demand_jitter, rv_frozen + ), "demand_jitter should be a SciPy probability distribution" + self.cost_jitter = cost_jitter + self.demand_jitter = demand_jitter + + def perturb( + self, + instance: UnitCommitmentData, + n_samples: int, + ) -> List[UnitCommitmentData]: + def _sample() -> UnitCommitmentData: + T, G = len(instance.demand), len(instance.max_power) + demand = instance.demand * self.demand_jitter.rvs(T) + cost_startup = instance.cost_startup * self.cost_jitter.rvs(G) + cost_prod = instance.cost_prod * self.cost_jitter.rvs(G) + cost_prod_quad = instance.cost_prod_quad * self.cost_jitter.rvs(G) + cost_fixed = instance.cost_fixed * self.cost_jitter.rvs(G) + + return UnitCommitmentData( + demand.round(2), + instance.min_power, + instance.max_power, + instance.min_uptime, + instance.min_downtime, + cost_startup.round(2), + cost_prod.round(2), + cost_prod_quad.round(4), + cost_fixed.round(2), + ) return [_sample() for _ in range(n_samples)] diff --git a/tests/problems/test_uc.py b/tests/problems/test_uc.py index e471f8f..d90633f 100644 --- a/tests/problems/test_uc.py +++ b/tests/problems/test_uc.py @@ -25,9 +25,6 @@ def test_generator() -> None: cost_fixed=uniform(loc=1, scale=1), min_uptime=randint(low=1, high=8), min_downtime=randint(low=1, high=8), - cost_jitter=uniform(loc=0.75, scale=0.5), - demand_jitter=uniform(loc=0.9, scale=0.2), - fix_units=True, ) data = gen.generate(1) assert data[0].demand.tolist() == [430.3, 511.29, 484.91, 860.61]