how to create a simple capacity expansion model while maximizing the revenue using PYOMO?
I have create a capacity expansion model using PYOMO, where there are 2 main generators gen1 and gen2 and one backup generator lost_load. The logic is pretty simple, gen1 and gen2 will run to fulfill the demand profile and they will get the money for that and there is some negative cost associated with gen1 and gen2 capacity. Here is the code:
import datetime
import pandas as pd
import numpy as np
from pyomo.environ import *
model = ConcreteModel()
np.random.seed(24)
load_profile = np.random.randint(90, 120, 24)
model.m_index = Set(initialize=list(range(len(load_profile))))
model.grid = Var(model.m_index, domain=NonNegativeReals)
# Gen1 variable
model.gen1_cap = Var(domain=NonNegativeReals)
model.gen1_use = Var(model.m_index, domain=NonNegativeReals)
# Gen2 variables
model.gen2_cap = Var(domain=NonNegativeReals)
model.gen2_use = Var(model.m_index, domain=NonNegativeReals)
# Load profile
model.load_profile = Param(model.m_index, initialize=dict(zip(model.m_index, load_profile)))
model.lost_load = Var(model.m_index, domain=NonNegativeReals)
# Objective function
def revenue(model):
total_revenue = sum(
model.gen1_use[m] * 5.2 +
model.gen2_use[m] * 6.1 +
model.lost_load[m] * -100
for m in model.m_index)
total_fixed_cost = model.gen1_cap * -45 + model.gen2_cap * -50
total_cost = total_revenue + total_fixed_cost
return total_cost
model.obj = Objective(rule=revenue, sense=maximize)
# When i=0
def energy_balance1(model, m):
return model.grid[m] <= model.gen1_use[m] + model.gen2_use[m] + model.lost_load[m]
model.energy_balance1 = Constraint(model.m_index, rule=energy_balance1)
def grid_limit(model, m):
return model.grid[m] == model.load_profile[m]
model.grid_limit = Constraint(model.m_index, rule=grid_limit)
def max_gen1(model, m):
eq = model.gen1_use[m] <= model.gen1_cap
return eq
model.max_gen1 = Constraint(model.m_index, rule=max_gen1)
def max_gen2(model, m):
eq = model.gen2_use[m] <= model.gen2_cap
return eq
model.max_gen2 = Constraint(model.m_index, rule=max_gen2)
Solver = SolverFactory('gurobi')
Solver.options['LogFile'] = "gurobiLog"
# Solver.options['MIPGap'] = 0.50
print('\nConnecting to Gurobi Server...')
results = Solver.solve(model)
if (results.solver.status == SolverStatus.ok):
if (results.solver.termination_condition == TerminationCondition.optimal):
print("\n\n***Optimal solution found***")
print('obj returned:', round(value(model.obj), 2))
else:
print("\n\n***No optimal solution found***")
if (results.solver.termination_condition == TerminationCondition.infeasible):
print("Infeasible solution")
exit()
else:
print("\n\n***Solver terminated abnormally***")
exit()
grid_use = []
gen1 = []
gen2 = []
lost_load = []
load = []
for i in range(len(load_profile)):
grid_use.append(value(model.grid[i]))
gen1.append(value(model.gen1_use[i]))
gen2.append(value(model.gen2_use[i]))
lost_load.append(value(model.lost_load[i]))
load.append(value(model.load_profile[i]))
print('gen1 capacity: ', value(model.gen1_cap))
print('gen2 capacity: ', value(model.gen2_cap))
pd.DataFrame({
'Grid': grid_use,
'Gen1': gen1,
'Gen2': gen2,
'Shortfall': lost_load,
'Load': load
}).to_excel('capacity expansion.xlsx')
This model will work if :
def energy_balance1(model, m):
return model.grid[m] == model.gen1_use[m] + model.gen2_use[m] + model.lost_load[m]
model.energy_balance1 = Constraint(model.m_index, rule=energy_balance1)
But if I set it to (<=) greater than equal to, it gives me an infeasible error. I want this condition because in actuality I want gen1 and gen2 to run at 100% capacity at all 24 instances. After all, in the objective function, we can see the positive cost associated with gen1_use and gen2_use. So to maximize the revenue they should run at 100% capacity even if the load profile is fulfilled. Here is the output I am getting:
optimal gen1 capacity: 9MW
optimal gen2 capacity: 108MW
Although I want the dispatch something like this :
Solution
Try something like this to get started. I avoided making capacity a variable by including several distinct generator choices, which seems very reasonable.
This also incorporates a binary "use" variable, which I think you need from the described problem.
CODE:
import pyomo.environ as pyo
demand = [100, 200, 50, 22, 300]
solar_wind = [20, 100, 100, 0, 100]
net_demand = [d - s for (d, s) in zip(demand, solar_wind)]
print(net_demand)
# generator options...
gen_capacity = {'gen1': 25, 'gen2': 25, 'gen3': 50, 'gen4': 100}
running_cost = {'gen1': 12, 'gen2': 12, 'gen3': 30, 'gen4': 50}
shortage_cost = 100
m = pyo.ConcreteModel()
### SETS
m.T = pyo.Set(initialize=range(len(demand))) # time periods
m.G = pyo.Set(initialize=gen_capacity.keys()) # generators
### VARS
m.install = pyo.Var(m.G, domain=pyo.Binary)
m.use = pyo.Var(m.G, m.T, domain=pyo.Binary)
m.shortage = pyo.Var(m.T, domain=pyo.NonNegativeReals)
### PARAMS
m.use_cost = pyo.Param(m.G, initialize=running_cost)
m.capacity = pyo.Param(m.G, initialize=gen_capacity)
### OBJ
m.cost = pyo.Objective(expr=sum(m.use[g, t] * m.use_cost[g] for g in m.G for t in m.T) +
sum(m.shortage[t] * shortage_cost for t in m.T),
sense=pyo.minimize)
### CONSTRAINTS
def demand(m, t):
return sum(m.use[g, t] * m.capacity[g] for g in m.G) + m.shortage[t] >= net_demand[t]
m.C1 = pyo.Constraint(m.T, rule=demand, doc='meet demand')
def use_if_installed(m, g, t):
return m.use[g, t] <= m.install[g]
m.C2 = pyo.Constraint(m.G, m.T, rule=use_if_installed, doc='can only use if installed')
m.C3 = pyo.Constraint(expr=sum(m.install[g] for g in m.G) <= 2, doc='install at most 2 gens')
solver = pyo.SolverFactory('cbc')
result = solver.solve(m)
print(result)
m.display()
Output (a bit long)
[80, 100, -50, 22, 200]
Problem:
- Name: unknown
Lower bound: 5210.0
Upper bound: 5210.0
Number of objectives: 1
Number of constraints: 25
Number of variables: 28
Number of binary variables: 24
Number of integer variables: 24
Number of nonzeros: 24
Sense: minimize
Solver:
- Status: ok
User time: -1.0
System time: 0.0
Wallclock time: 0.0
Termination condition: optimal
Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Black box:
Number of iterations: 0
Error rc: 0
Time: 0.012147188186645508
Solution:
- number of solutions: 0
number of solutions displayed: 0
Model unknown
Variables:
install : Size=4, Index=G
Key : Lower : Value : Upper : Fixed : Stale : Domain
gen1 : 0 : 0.0 : 1 : False : False : Binary
gen2 : 0 : 0.0 : 1 : False : False : Binary
gen3 : 0 : 1.0 : 1 : False : False : Binary
gen4 : 0 : 1.0 : 1 : False : False : Binary
use : Size=20, Index=use_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
('gen1', 0) : 0 : 0.0 : 1 : False : False : Binary
('gen1', 1) : 0 : 0.0 : 1 : False : False : Binary
('gen1', 2) : 0 : 0.0 : 1 : False : False : Binary
('gen1', 3) : 0 : 0.0 : 1 : False : False : Binary
('gen1', 4) : 0 : 0.0 : 1 : False : False : Binary
('gen2', 0) : 0 : 0.0 : 1 : False : False : Binary
('gen2', 1) : 0 : 0.0 : 1 : False : False : Binary
('gen2', 2) : 0 : 0.0 : 1 : False : False : Binary
('gen2', 3) : 0 : 0.0 : 1 : False : False : Binary
('gen2', 4) : 0 : 0.0 : 1 : False : False : Binary
('gen3', 0) : 0 : 0.0 : 1 : False : False : Binary
('gen3', 1) : 0 : 0.0 : 1 : False : False : Binary
('gen3', 2) : 0 : 0.0 : 1 : False : False : Binary
('gen3', 3) : 0 : 1.0 : 1 : False : False : Binary
('gen3', 4) : 0 : 1.0 : 1 : False : False : Binary
('gen4', 0) : 0 : 1.0 : 1 : False : False : Binary
('gen4', 1) : 0 : 1.0 : 1 : False : False : Binary
('gen4', 2) : 0 : 0.0 : 1 : False : False : Binary
('gen4', 3) : 0 : 0.0 : 1 : False : False : Binary
('gen4', 4) : 0 : 1.0 : 1 : False : False : Binary
shortage : Size=5, Index=T
Key : Lower : Value : Upper : Fixed : Stale : Domain
0 : 0 : 0.0 : None : False : False : NonNegativeReals
1 : 0 : 0.0 : None : False : False : NonNegativeReals
2 : 0 : 0.0 : None : False : False : NonNegativeReals
3 : 0 : 0.0 : None : False : False : NonNegativeReals
4 : 0 : 50.0 : None : False : False : NonNegativeReals
Objectives:
cost : Size=1, Index=None, Active=True
Key : Active : Value
None : True : 5210.0
Constraints:
C1 : Size=5
Key : Lower : Body : Upper
0 : 80.0 : 100.0 : None
1 : 100.0 : 100.0 : None
2 : -50.0 : 0.0 : None
3 : 22.0 : 50.0 : None
4 : 200.0 : 200.0 : None
C2 : Size=20
Key : Lower : Body : Upper
('gen1', 0) : None : 0.0 : 0.0
('gen1', 1) : None : 0.0 : 0.0
('gen1', 2) : None : 0.0 : 0.0
('gen1', 3) : None : 0.0 : 0.0
('gen1', 4) : None : 0.0 : 0.0
('gen2', 0) : None : 0.0 : 0.0
('gen2', 1) : None : 0.0 : 0.0
('gen2', 2) : None : 0.0 : 0.0
('gen2', 3) : None : 0.0 : 0.0
('gen2', 4) : None : 0.0 : 0.0
('gen3', 0) : None : -1.0 : 0.0
('gen3', 1) : None : -1.0 : 0.0
('gen3', 2) : None : -1.0 : 0.0
('gen3', 3) : None : 0.0 : 0.0
('gen3', 4) : None : 0.0 : 0.0
('gen4', 0) : None : 0.0 : 0.0
('gen4', 1) : None : 0.0 : 0.0
('gen4', 2) : None : -1.0 : 0.0
('gen4', 3) : None : -1.0 : 0.0
('gen4', 4) : None : 0.0 : 0.0
C3 : Size=1
Key : Lower : Body : Upper
None : None : 2.0 : 2.0
- how to average in a specific dimension with numpy.mean?
- Removed Realm, but still getting this error: module importing failed: invalid token (rlm_lldb.py, line 37) File "temp.py", line 1,
- Bool object does not support item assignment
- Conda env vs venv / pyenv / virtualenv / poetry / docker, etc
- How to import values from excel with pandas into tkinter faster?
- Speed up concatenation of excel files with Pandas
- Faster way to read Excel files to pandas dataframe
- Create Azure TimerTrigger Durable Function in python
- Read same file from src and test in PyCharm
- How to use Cython to compile Python 3 into C
- Outputting variable from one thread to other
- How to POST a JSON having a single body parameter in FastAPI?
- touchscreen raspberry pi phantom touch (reset fixes)
- How to determine the encoding of text
- Specify extras_require with pip install -e
- Plotting a second scaled y axis in matplotlib from one set of data
- How can I call scikit-learn classifiers from Java?
- How does Pandas.Series.nbytes work for strings? Results don't seem to match expectations
- How to install ursina on pydroid?
- auditlog with Django and DRF
- Python: How to RECURSIVELY remove None values from a NESTED data structure (lists and dictionaries)?
- How to Mark Repeated Entries as True Starting from the Second Occurrence Using NumPy?
- Python3.5 send object over socket (Pygame cam image)
- Value based partial slicing with non-existing keys is now deprecated
- Camelot: DeprecationError: PdfFileReader is deprecated
- How can I partially fill a numpy array with a value given a range?
- Conda wont install pdfplumber
- tomli-w, adding arrays from a python script
- Override value in pydantic model with environment variable
- Polars-Python: Compare list columns