1. Knapsack
import numpy as np
import gurobipy as gp
from gurobipy import GRB
def generate_knapsack(num_items):
# Fix seed value
rng = np.random.default_rng(seed=0)
# Item values, weights
values = rng.uniform(low=1, high=25, size=num_items)
weights = rng.uniform(low=5, high=100, size=num_items)
# Knapsack capacity
capacity = 0.7 * weights.sum()
return values, weights, capacity
def solve_knapsack_model(values, weights, capacity):
num_items = len(values)
# Turn values and weights numpy arrays to dict
values_dict = {i: value for i, value in enumerate(values)}
weights_dict = {i: weight for i, weight in enumerate(weights)}
with gp.Env() as env:
with gp.Model(name="knapsack", env=env) as model:
# Define decision variables using the Model.addVars() method
x = model.addVars(num_items, vtype=GRB.BINARY, name="x")
# Define objective function using the Model.setObjective() method
# Build the LinExpr using the tupledict.prod() method
model.setObjective(x.prod(values_dict), sense=GRB.MAXIMIZE)
# Define capacity constraint using the Model.addConstr() method
model.addConstr(x.prod(weights_dict) <= capacity, name="knapsack_capacity")
model.optimize()
data = generate_knapsack(10000)
solve_knapsack_model(*data)Here is the log of the execution of this program:
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.2.0 24C101)
CPU model: Apple M1 Pro
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 1 rows, 10000 columns and 10000 nonzeros
Model fingerprint: 0xa2a6b3bf
Variable types: 0 continuous, 10000 integer (10000 binary)
Coefficient statistics:
Matrix range [5e+00, 1e+02]
Objective range [1e+00, 2e+01]
Bounds range [1e+00, 1e+00]
RHS range [4e+05, 4e+05]
Found heuristic solution: objective 90594.480679
Presolve time: 0.01s
Presolved: 1 rows, 10000 columns, 10000 nonzeros
Variable types: 0 continuous, 10000 integer (10000 binary)
Root relaxation: objective 1.190206e+05, 1 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 119020.617 0 1 90594.4807 119020.617 31.4% - 0s
H 0 0 119005.88462 119020.617 0.01% - 0s
H 0 0 119015.22254 119020.617 0.00% - 0s
Explored 1 nodes (1 simplex iterations) in 0.03 seconds (0.02 work units)
Thread count was 8 (of 8 available processors)
Solution count 3: 119015 119006 90594.5
Optimal solution found (tolerance 1.00e-04)
Best objective 1.190152225370e+05, best bound 1.190206174693e+05, gap 0.0045%
2. Portfolio Optimization
import json
import pandas as pd
import numpy as np
import gurobipy as gp
from gurobipy import GRB
with open("data/portfolio-example.json", "r") as f:
data = json.load(f)
n = data["num_assets"]
sigma = np.array(data["covariance"])
mu = np.array(data["expected_return"])
mu_0 = data["target_return"]
k = data["portfolio_max_size"]
with gp.Model("portfolio") as model:
x = model.addVars(n, ub=1, name="x")
y = model.addVars(n, vtype=GRB.BINARY, name="y")
risk = gp.quicksum(x[i] * sigma[i, j] * x[j] for i in range(n) for j in range(n))
model.setObjective(risk)
expected_return = gp.quicksum(mu[i] * x[i] for i in range(n))
model.addConstr(expected_return >= mu_0, name="return")
model.addConstr(x.sum() == 1, name="budget")
model.addConstr(y.sum() <= k, name="cardinality")
model.addConstrs((x[i] <= y[i] for i in range(n)), name="is_allocated")
model.optimize()
# Write the solution into a DataFrame
portfolio = [var.X for var in model.getVars() if "x" in var.VarName]
risk = model.ObjVal
expected_return = model.getRow(model.getConstrByName("return")).getValue()
df = pd.DataFrame(
data=portfolio + [risk, expected_return],
index=[f"asset_{i}" for i in range(n)] + ["risk", "return"],
columns=["Portfolio"],
)
print(df)Here is the log of the execution of this program:
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.2.0 24C101)
CPU model: Apple M1 Pro
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 23 rows, 40 columns and 100 nonzeros
Model fingerprint: 0x755216b2
Model has 210 quadratic objective terms
Variable types: 20 continuous, 20 integer (20 binary)
Coefficient statistics:
Matrix range [1e-04, 1e+00]
Objective range [0e+00, 0e+00]
QObjective range [7e-06, 4e-03]
Bounds range [1e+00, 1e+00]
RHS range [1e-04, 2e+01]
Found heuristic solution: objective 0.0000960
Presolve time: 0.00s
Presolved: 23 rows, 40 columns, 100 nonzeros
Presolved model has 210 quadratic objective terms
Variable types: 20 continuous, 20 integer (20 binary)
Root relaxation: objective 6.741224e-05, 38 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 0.00007 0 11 0.00010 0.00007 29.8% - 0s
H 0 0 0.0000674 0.00007 0.00% - 0s
0 0 0.00007 0 11 0.00007 0.00007 0.00% - 0s
Explored 1 nodes (38 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)
Solution count 2: 6.74122e-05 9.60284e-05
Optimal solution found (tolerance 1.00e-04)
Best objective 6.741224163915e-05, best bound 6.741224163915e-05, gap 0.0000%
Portfolio
asset_0 0.000000
asset_1 0.000000
asset_2 0.000000
asset_3 0.000000
asset_4 0.006812
asset_5 0.034594
asset_6 0.000000
asset_7 0.000000
asset_8 0.409840
asset_9 0.026759
asset_10 0.000000
asset_11 0.027722
asset_12 0.039911
asset_13 0.030221
asset_14 0.038511
asset_15 0.000000
asset_16 0.184718
asset_17 0.000000
asset_18 0.126673
asset_19 0.074241
risk 0.000067
return 0.000104
3. Termination Criteria
from functools import partial
import gurobipy as gp
from gurobipy import GRB
class CallbackData:
def __init__(self):
self.last_gap_change_time = -GRB.INFINITY
self.last_gap = GRB.INFINITY
def callback(model, where, *, cbdata):
if where != GRB.Callback.MIP:
return
if model.cbGet(GRB.Callback.MIP_SOLCNT) == 0:
return
best = model.cbGet(GRB.Callback.MIP_OBJBST)
bound = model.cbGet(GRB.Callback.MIP_OBJBND)
gap = abs((bound - best) / best)
time = model.cbGet(GRB.Callback.RUNTIME)
if gap < cbdata.last_gap - epsilon_to_compare_gap:
# The gap has decreased, store new gap and reset the time
cbdata.last_gap = gap
cbdata.last_gap_change_time = time
print(f"Storing new gap: {gap}, {time}")
return
if time - cbdata.last_gap_change_time > max_time_between_gap_updates:
# No change since too long
print("It's been too long...")
model.terminate()
with gp.read("data/mkp.mps.bz2") as model:
# Global variables used in the callback function
max_time_between_gap_updates = 15
epsilon_to_compare_gap = 1e-4
# Initialize data passed to the callback function
callback_data = CallbackData()
callback_func = partial(callback, cbdata=callback_data)
model.optimize(callback_func)Here is the log of the execution of this program:
Read MPS format model from file data/mkp.mps.bz2
Reading time = 0.04 seconds
mkp_bip: 30 rows, 500 columns, 15000 nonzeros
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.2.0 24C101)
CPU model: Apple M1 Pro
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 30 rows, 500 columns and 15000 nonzeros
Model fingerprint: 0xcbe37219
Variable types: 0 continuous, 500 integer (500 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+03]
Objective range [4e+02, 1e+03]
Bounds range [1e+00, 1e+00]
RHS range [1e+05, 1e+05]
Found heuristic solution: objective 177106.00000
Presolve time: 0.02s
Presolved: 30 rows, 500 columns, 15000 nonzeros
Variable types: 0 continuous, 500 integer (500 binary)
Storing new gap: 1.0921143270132012, 0.016685009002685547
Root relaxation: objective 2.164011e+05, 112 iterations, 0.00 seconds (0.01 work units)
Storing new gap: 0.22187277675516356, 0.019383907318115234
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 216401.065 0 25 177106.000 216401.065 22.2% - 0s
H 0 0 209908.00000 216401.065 3.09% - 0s
H 0 0 213709.00000 216401.065 1.26% - 0s
H 0 0 215225.00000 216401.065 0.55% - 0s
H 0 0 215436.00000 216401.065 0.45% - 0s
H 0 0 215510.00000 216401.065 0.41% - 0s
H 0 0 215624.00000 216401.065 0.36% - 0s
Storing new gap: 0.0035756687567246688, 0.10831308364868164
0 0 216395.201 0 25 215624.000 216395.201 0.36% - 0s
H 0 0 215635.00000 216395.201 0.35% - 0s
H 0 2 215645.00000 216395.201 0.35% - 0s
0 2 216395.201 0 25 215645.000 216395.201 0.35% - 0s
H 85 88 215666.00000 216395.201 0.34% 8.8 0s
Storing new gap: 0.003380226832231321, 0.29430699348449707
H 162 166 215707.00000 216395.201 0.32% 9.4 0s
Storing new gap: 0.0031895116987394937, 0.33049607276916504
H 2099 1932 215727.00000 216389.474 0.31% 7.0 0s
Storing new gap: 0.0030686933021828516, 0.9242510795593262
H20129 16350 215731.00000 216372.312 0.30% 8.0 2s
Storing new gap: 0.0029620221479527745, 2.383863925933838
41136 31672 216209.291 61 22 215731.000 216367.041 0.29% 8.1 5s
H50272 34121 215821.00000 216367.041 0.25% 8.1 5s
Storing new gap: 0.0025298742939750996, 5.620049953460693
137609 94670 216027.180 88 21 215821.000 216348.997 0.24% 8.5 10s
Storing new gap: 0.002427937967111634, 11.783638954162598
232901 178512 216148.976 76 22 215821.000 216341.790 0.24% 8.8 15s
329701 262798 216036.548 99 18 215821.000 216337.190 0.24% 9.0 20s
429310 349889 216128.018 73 23 215821.000 216333.650 0.24% 9.1 25s
It's been too long...
Cutting planes:
Gomory: 7
Explored 469770 nodes (4286693 simplex iterations) in 26.80 seconds (74.85 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: 215821 215731 215727 ... 215436
Solve interrupted
Best objective 2.158210000000e+05, best bound 2.163320000000e+05, gap 0.2368%
User-callback calls 942664, time in user-callback 1.06 sec
4. Unit Commitment Problem
import gurobipy as gp
from gurobipy import GRB
# 24 Hour Load Forecast (MW)
load_forecast = [
4, 4, 4, 4, 4, 4, 6, 6,
12, 12, 12, 12, 12, 4, 4, 4,
4, 16, 16, 16, 16, 6.5, 6.5, 6.5,
]
# solar energy forecast (MW)
solar_forecast = [
0, 0, 0, 0, 0, 0, 0.5, 1.0,
1.5, 2.0, 2.5, 3.5, 3.5, 2.5, 2.0, 1.5,
1.0, 0.5, 0, 0, 0, 0, 0, 0,
]
# global number of time intervals
nTimeIntervals = len(load_forecast)
# thermal units
thermal_units = ["gen1", "gen2", "gen3"]
# thermal units' costs (a + b*p + c*p^2), (startup and shutdown costs)
thermal_units_cost, a, b, c, sup_cost, sdn_cost = gp.multidict(
{
"gen1": [5.0, 0.5, 1.0, 2, 1],
"gen2": [5.0, 0.5, 0.5, 2, 1],
"gen3": [5.0, 3.0, 2.0, 2, 1],
}
)
# thernal units operating limits
thermal_units_limits, pmin, pmax = gp.multidict(
{"gen1": [1.5, 5.0], "gen2": [2.5, 10.0], "gen3": [1.0, 3.0]}
)
# thermal units dynamic data (initial commitment status)
thermal_units_dyn_data, init_status = gp.multidict(
{"gen1": [0], "gen2": [0], "gen3": [0]}
)
def show_results():
obj_val_s = model.ObjVal
print(f" OverAll Cost = {round(obj_val_s, 2)} ")
print("\n")
print("%5s" % "time", end=" ")
for t in range(nTimeIntervals):
print("%4s" % t, end=" ")
print("\n")
for g in thermal_units:
print("%5s" % g, end=" ")
for t in range(nTimeIntervals):
print("%4.1f" % thermal_units_out_power[g, t].X, end=" ")
print("\n")
print("%5s" % "Solar", end=" ")
for t in range(nTimeIntervals):
print("%4.1f" % solar_forecast[t], end=" ")
print("\n")
print("%5s" % "Load", end=" ")
for t in range(nTimeIntervals):
print("%4.1f" % load_forecast[t], end=" ")
print("\n")
with gp.Env() as env, gp.Model(env=env) as model:
# add variables for thermal units (power and statuses for commitment, startup and shutdown)
thermal_units_out_power = model.addVars(
thermal_units, range(nTimeIntervals), name="thermal_units_out_power"
)
thermal_units_startup_status = model.addVars(
thermal_units,
range(nTimeIntervals),
vtype=GRB.BINARY,
name="thermal_unit_startup_status",
)
thermal_units_shutdown_status = model.addVars(
thermal_units,
range(nTimeIntervals),
vtype=GRB.BINARY,
name="thermal_unit_shutdown_status",
)
thermal_units_comm_status = model.addVars(
thermal_units, range(nTimeIntervals), vtype=GRB.BINARY, name="thermal_unit_comm_status"
)
# define objective function as an empty quadratic construct and add terms
obj_fun_expr = gp.QuadExpr(0)
for t in range(nTimeIntervals):
for g in thermal_units:
obj_fun_expr.add(
a[g] * thermal_units_comm_status[g, t]
+ b[g] * thermal_units_out_power[g, t]
+ c[g] * thermal_units_out_power[g, t] * thermal_units_out_power[g, t]
+ sup_cost[g] * thermal_units_startup_status[g, t]
+ sdn_cost[g] * thermal_units_shutdown_status[g, t]
)
model.setObjective(obj_fun_expr)
# Power balance equations
for t in range(nTimeIntervals):
model.addConstr(
gp.quicksum(thermal_units_out_power[g, t] for g in thermal_units)
== load_forecast[t] - solar_forecast[t],
name="power_balance_" + str(t),
)
# Thermal units logical constraints
for t in range(nTimeIntervals):
for g in thermal_units:
if t == 0:
model.addConstr(
thermal_units_comm_status[g, t] - init_status[g]
== thermal_units_startup_status[g, t]
- thermal_units_shutdown_status[g, t],
name="logical1_" + g + "_" + str(t),
)
else:
model.addConstr(
thermal_units_comm_status[g, t]
- thermal_units_comm_status[g, t - 1]
== thermal_units_startup_status[g, t]
- thermal_units_shutdown_status[g, t],
name="logical1_" + g + "_" + str(t),
)
model.addConstr(
thermal_units_startup_status[g, t] + thermal_units_shutdown_status[g, t]
<= 1,
name="logical2_" + g + "_" + str(t),
)
# Thermal units physical constraints, using indicator constraints
for t in range(nTimeIntervals):
for g in thermal_units:
model.addConstr(
(thermal_units_comm_status[g, t] == 1) >>
(thermal_units_out_power[g, t] >= pmin[g])
)
model.addConstr(
(thermal_units_comm_status[g, t] == 1) >>
(thermal_units_out_power[g, t] <= pmax[g])
)
model.addConstr(
(thermal_units_comm_status[g, t] == 0) >>
(thermal_units_out_power[g, t] == 0.0)
)
model.optimize()
show_results()Here is the log of the execution of this program:
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.2.0 24C101)
CPU model: Apple M1 Pro
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 168 rows, 288 columns and 501 nonzeros
Model fingerprint: 0xa1375406
Model has 72 quadratic objective terms
Model has 216 simple general constraints
216 INDICATOR
Variable types: 72 continuous, 216 integer (216 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [5e-01, 5e+00]
QObjective range [1e+00, 4e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+01]
GenCon rhs range [1e+00, 1e+01]
GenCon coe range [1e+00, 1e+00]
Presolve added 1 rows and 0 columns
Presolve removed 0 rows and 125 columns
Presolve time: 0.01s
Presolved: 185 rows, 171 columns, 473 nonzeros
Presolved model has 43 quadratic objective terms
Variable types: 51 continuous, 120 integer (120 binary)
Found heuristic solution: objective 1007.6964286
Found heuristic solution: objective 907.5714286
Root relaxation: objective 8.586741e+02, 287 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 858.67410 0 44 907.57143 858.67410 5.39% - 0s
H 0 0 894.0714286 858.67410 3.96% - 0s
H 0 0 888.7380952 858.67410 3.38% - 0s
H 0 0 887.9047619 858.67410 3.29% - 0s
0 0 861.09989 0 48 887.90476 861.09989 3.02% - 0s
0 0 861.09989 0 35 887.90476 861.09989 3.02% - 0s
0 0 865.42295 0 37 887.90476 865.42295 2.53% - 0s
0 0 866.16295 0 37 887.90476 866.16295 2.45% - 0s
0 0 868.53847 0 41 887.90476 868.53847 2.18% - 0s
0 0 869.53572 0 42 887.90476 869.53572 2.07% - 0s
0 0 871.70239 0 42 887.90476 871.70239 1.82% - 0s
0 0 871.70239 0 44 887.90476 871.70239 1.82% - 0s
0 0 871.70239 0 44 887.90476 871.70239 1.82% - 0s
0 0 871.70239 0 44 887.90476 871.70239 1.82% - 0s
0 0 871.70239 0 44 887.90476 871.70239 1.82% - 0s
0 0 871.70239 0 46 887.90476 871.70239 1.82% - 0s
0 2 871.70239 0 46 887.90476 871.70239 1.82% - 0s
Cutting planes:
Cover: 2
Implied bound: 1
Clique: 1
MIR: 1
Flow cover: 7
Explored 52 nodes (1017 simplex iterations) in 0.05 seconds (0.03 work units)
Thread count was 8 (of 8 available processors)
Solution count 5: 887.905 888.738 894.071 ... 1007.7
Optimal solution found (tolerance 1.00e-04)
Best objective 8.879047619048e+02, best bound 8.879047619048e+02, gap 0.0000%
OverAll Cost = 887.9
time 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
gen1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.5 3.3 3.2 2.8 2.8 1.5 2.0 0.0 0.0 4.6 4.8 4.8 4.8 2.2 2.2 2.2
gen2 4.0 4.0 4.0 4.0 4.0 4.0 5.5 5.0 7.0 6.7 6.3 5.7 5.7 0.0 0.0 2.5 3.0 9.2 9.5 9.5 9.5 4.3 4.3 4.3
gen3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7 1.8 1.8 1.8 0.0 0.0 0.0
Solar 0.0 0.0 0.0 0.0 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.5 3.5 2.5 2.0 1.5 1.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0
Load 4.0 4.0 4.0 4.0 4.0 4.0 6.0 6.0 12.0 12.0 12.0 12.0 12.0 4.0 4.0 4.0 4.0 16.0 16.0 16.0 16.0 6.5 6.5 6.5
5. Unit Commitment Problem (matrix API)
import gurobipy as gp
from gurobipy import GRB
import numpy as np
# 24 Hour Load Forecast (MW)
load_forecast = [
4, 4, 4, 4, 4, 4, 6, 6,
12, 12, 12, 12, 12, 4, 4, 4,
4, 16, 16, 16, 16, 6.5, 6.5, 6.5,
]
# solar energy forecast (MW)
solar_forecast = [
0, 0, 0, 0, 0, 0, 0.5, 1.0,
1.5, 2.0, 2.5, 3.5, 3.5, 2.5, 2.0, 1.5,
1.0, 0.5, 0, 0, 0, 0, 0, 0,
]
# global number of time intervals
nTimeIntervals = len(load_forecast)
# thermal units
thermal_units = ["gen1", "gen2", "gen3"]
# thermal units' costs (a + b*p + c*p^2), (startup and shutdown costs)
thermal_units_cost, a, b, c, sup_cost, sdn_cost = gp.multidict(
{
"gen1": [5.0, 0.5, 1.0, 2, 1],
"gen2": [5.0, 0.5, 0.5, 2, 1],
"gen3": [5.0, 3.0, 2.0, 2, 1],
}
)
# thernal units operating limits
thermal_units_limits, pmin, pmax = gp.multidict(
{"gen1": [1.5, 5.0], "gen2": [2.5, 10.0], "gen3": [1.0, 3.0]}
)
# thermal units dynamic data (initial commitment status)
thermal_units_dyn_data, init_status = gp.multidict(
{"gen1": [0], "gen2": [0], "gen3": [0]}
)
# We need np.array instances instead of dict instances to use
# the matrix API
def dict_to_array(keys, d):
return np.array([d[k] for k in keys])
a = dict_to_array(thermal_units, a)
b = dict_to_array(thermal_units, b)
c = dict_to_array(thermal_units, c)
sup_cost = dict_to_array(thermal_units, sup_cost)
sdn_cost = dict_to_array(thermal_units, sdn_cost)
pmin = dict_to_array(thermal_units, pmin)
pmax = dict_to_array(thermal_units, pmax)
init_status = dict_to_array(thermal_units, init_status)
# We also need to turn to integer indices instead of strings
units_to_index = {u: idx for idx, u in enumerate(thermal_units)}
nb_units = len(thermal_units)
def show_results():
obj_val_s = model.ObjVal
print(f" OverAll Cost = {round(obj_val_s, 2)} ")
print("\n")
print("%5s" % "time", end=" ")
for t in range(nTimeIntervals):
print("%4s" % t, end=" ")
print("\n")
for g in thermal_units:
print("%5s" % g, end=" ")
for t in range(nTimeIntervals):
print("%4.1f" % thermal_units_out_power[units_to_index[g], t].X, end=" ")
print("\n")
print("%5s" % "Solar", end=" ")
for t in range(nTimeIntervals):
print("%4.1f" % solar_forecast[t], end=" ")
print("\n")
print("%5s" % "Load", end=" ")
for t in range(nTimeIntervals):
print("%4.1f" % load_forecast[t], end=" ")
print("\n")
with gp.Env() as env, gp.Model(env=env) as model:
# add variables for thermal units (power and statuses for commitment, startup and shutdown)
thermal_units_out_power = model.addMVar(
(nb_units, nTimeIntervals), name="thermal_units_out_power"
)
thermal_units_startup_status = model.addMVar(
(nb_units, nTimeIntervals),
vtype=GRB.BINARY,
name="thermal_unit_startup_status",
)
thermal_units_shutdown_status = model.addMVar(
(nb_units, nTimeIntervals),
vtype=GRB.BINARY,
name="thermal_unit_shutdown_status",
)
thermal_units_comm_status = model.addMVar(
(nb_units, nTimeIntervals), vtype=GRB.BINARY, name="thermal_unit_comm_status"
)
# define objective function as an empty quadratic construct and add terms
model.setObjective((a[:, None] * thermal_units_comm_status).sum() +
(b[:, None] * thermal_units_out_power).sum() +
(c[:, None] * (thermal_units_out_power * thermal_units_out_power)).sum() +
(sup_cost[:, None] * thermal_units_startup_status).sum() +
(sdn_cost[:, None] * thermal_units_shutdown_status).sum())
# Power balance equations
required_power = np.array([l - s for l, s in zip(load_forecast, solar_forecast) ])
model.addConstr(
thermal_units_out_power.sum(axis=0) == required_power,
name="power_balance"
)
# Thermal units logical constraints
# Initial condition (t = 0)
model.addConstr(
thermal_units_comm_status[:, 0] - init_status ==
thermal_units_startup_status[:, 0] - thermal_units_shutdown_status[:, 0],
name="logical1_initial"
)
# Remaining time intervals (t > 0)
model.addConstr(
thermal_units_comm_status[:, 1:] - thermal_units_comm_status[:, :-1]
== thermal_units_startup_status[:, 1:] - thermal_units_shutdown_status[:, 1:],
name="logical1_remaining"
)
model.addConstr(
thermal_units_startup_status + thermal_units_shutdown_status <= 1,
name="logical2"
)
# Thermal units physical constraints, using indicator constraints
for g in range(len(thermal_units)):
model.addGenConstrIndicator(thermal_units_comm_status[g, :], True,
thermal_units_out_power[g, :] >= pmin[g],
name="physical_min")
model.addGenConstrIndicator(thermal_units_comm_status[g, :], True,
thermal_units_out_power[g, :] <= pmax[g],
name="physical_max")
model.addGenConstrIndicator(thermal_units_comm_status[g, :], False,
thermal_units_out_power[g, :] == 0.0,
name="physical_off")
model.optimize()
show_results()Here is the log of the execution of this program:
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (mac64[arm] - Darwin 24.2.0 24C101)
CPU model: Apple M1 Pro
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 168 rows, 288 columns and 501 nonzeros
Model fingerprint: 0x3b9885a9
Model has 72 quadratic objective terms
Model has 216 simple general constraints
216 INDICATOR
Variable types: 72 continuous, 216 integer (216 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [5e-01, 5e+00]
QObjective range [1e+00, 4e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+01]
GenCon rhs range [1e+00, 1e+01]
GenCon coe range [1e+00, 1e+00]
Presolve added 1 rows and 0 columns
Presolve removed 0 rows and 125 columns
Presolve time: 0.01s
Presolved: 183 rows, 170 columns, 468 nonzeros
Presolved model has 43 quadratic objective terms
Variable types: 50 continuous, 120 integer (120 binary)
Found heuristic solution: objective 972.4464286
Found heuristic solution: objective 949.4464286
Root relaxation: objective 8.562274e+02, 297 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 856.22741 0 47 949.44643 856.22741 9.82% - 0s
H 0 0 888.6964286 856.22741 3.65% - 0s
H 0 0 887.9047619 856.22741 3.57% - 0s
0 0 859.55567 0 49 887.90476 859.55567 3.19% - 0s
0 0 859.55567 0 40 887.90476 859.55567 3.19% - 0s
0 0 860.62679 0 43 887.90476 860.62679 3.07% - 0s
0 0 861.04305 0 43 887.90476 861.04305 3.03% - 0s
0 0 863.11702 0 40 887.90476 863.11702 2.79% - 0s
0 0 863.60301 0 40 887.90476 863.60301 2.74% - 0s
0 0 864.38939 0 42 887.90476 864.38939 2.65% - 0s
0 0 865.38905 0 42 887.90476 865.38905 2.54% - 0s
0 0 865.80613 0 42 887.90476 865.80613 2.49% - 0s
0 0 865.80613 0 42 887.90476 865.80613 2.49% - 0s
0 0 865.81082 0 42 887.90476 865.81082 2.49% - 0s
0 2 865.81082 0 42 887.90476 865.81082 2.49% - 0s
Cutting planes:
Learned: 1
Cover: 2
MIR: 1
Flow cover: 2
Explored 147 nodes (1455 simplex iterations) in 0.05 seconds (0.03 work units)
Thread count was 8 (of 8 available processors)
Solution count 4: 887.905 888.696 949.446 972.446
Optimal solution found (tolerance 1.00e-04)
Best objective 8.879047619048e+02, best bound 8.879047606600e+02, gap 0.0000%
OverAll Cost = 887.9
time 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
gen1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.5 3.3 3.2 2.8 2.8 1.5 2.0 0.0 0.0 4.6 4.8 4.8 4.8 2.2 2.2 2.2
gen2 4.0 4.0 4.0 4.0 4.0 4.0 5.5 5.0 7.0 6.7 6.3 5.7 5.7 0.0 0.0 2.5 3.0 9.2 9.5 9.5 9.5 4.3 4.3 4.3
gen3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7 1.8 1.8 1.8 0.0 0.0 0.0
Solar 0.0 0.0 0.0 0.0 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.5 3.5 2.5 2.0 1.5 1.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0
Load 4.0 4.0 4.0 4.0 4.0 4.0 6.0 6.0 12.0 12.0 12.0 12.0 12.0 4.0 4.0 4.0 4.0 16.0 16.0 16.0 16.0 6.5 6.5 6.5