I have formulated the following problem that tries to minimize the cost of imported energy by controlling a battery:
import pyomo.environ as pyo
def self_consumption(pv, demand, prices):
pv.index = np.arange(1, len(pv) + 1)
demand.index = np.arange(1, len(demand) + 1)
prices.index = np.arange(1, len(prices) + 1)
# Define the model
model = pyo.ConcreteModel()
# Define the set of timesteps
model.timesteps = pyo.Set(initialize = pyo.RangeSet(len(pv)), ordered = True)
# Define the inputs of the model
model.b_efficiency = pyo.Param(initialize = ETTA)
model.b_cap = pyo.Param(initialize = BATTERY_CAPACITY)
model.b_min_soc = pyo.Param(initialize = BATTERY_SOC_MIN)
model.b_max_soc = pyo.Param(initialize = BATTERY_SOC_MAX)
model.b_charging_rate = pyo.Param(initialize = BATTERY_CHARGE_RATE)
model.Ppv = pyo.Param(model.timesteps, initialize = pv.to_dict()['value'], within = pyo.Any)
model.Pdemand = pyo.Param(model.timesteps, initialize = demand.to_dict()['value'], within = pyo.Any)
model.day_ahead_prices = pyo.Param(model.timesteps, initialize = prices.to_dict(), within = pyo.Any)
# Define the decision variables of the model
model.Pbat_ch = pyo.Var(model.timesteps, within = pyo.NonNegativeReals, bounds = (0, model.b_charging_rate))
model.Pbat_dis = pyo.Var(model.timesteps, within = pyo.NonNegativeReals, bounds = (0, model.b_charging_rate))
model.Ebat = pyo.Var(model.timesteps, within = pyo.NonNegativeReals, bounds = (model.b_min_soc * model.b_cap, model.b_max_soc * model.b_cap))
model.Pgrid = pyo.Var(model.timesteps)
model.is_charging = pyo.Var(model.timesteps, within = pyo.Binary)
model.AbsPgrid = pyo.Var(model.timesteps, within=pyo.NonNegativeReals)
# Define the constraints of the model
def BatEnergyRule(model, t):
if t == 1:
return model.Ebat[t] == model.b_cap/2 # battery initialization at half of the capacity (assumption)
else:
return model.Ebat[t] == model.Ebat[t-1] + (model.b_efficiency * model.Pbat_ch[t] - model.Pbat_dis[t]/model.b_efficiency)
model.cons1 = pyo.Constraint(model.timesteps, rule = BatEnergyRule)
def PowerBalanceRule(model, t):
return model.Pgrid[t] == model.Pdemand[t] - model.Ppv[t] + model.Pbat_ch[t] - model.Pbat_dis[t]
model.cons2 = pyo.Constraint(model.timesteps, rule = PowerBalanceRule)
# Charge/Discharge constraint
def charge_discharge_rule(model, t):
return model.Pbat_ch[t] <= model.is_charging[t] * model.b_charging_rate
model.charge_constraint = pyo.Constraint(model.timesteps, rule=charge_discharge_rule)
def discharge_charge_rule(model, t):
return model.Pbat_dis[t] <= (1 - model.is_charging[t]) * model.b_charging_rate
model.discharge_constraint = pyo.Constraint(model.timesteps, rule=discharge_charge_rule)
def AbsPgridRule1(model, t):
return model.AbsPgrid[t] >= model.Pgrid[t]
model.cons6 = pyo.Constraint(model.timesteps, rule=AbsPgridRule1)
def AbsPgridRule2(model, t):
return model.AbsPgrid[t] >= -model.Pgrid[t]
model.cons7 = pyo.Constraint(model.timesteps, rule=AbsPgridRule2)
# Define the objective function
def ObjRule(model):
return sum(model.day_ahead_prices[t] * model.AbsPgrid[t] for t in model.timesteps)
model.obj = pyo.Objective(rule=ObjRule, sense=pyo.minimize)
# Choose a solver
opt = pyo.SolverFactory('glpk')
# Solve the optimization problem
result = opt.solve(model, tee = True)
solution = {
"Pbat_ch": {t: model.Pbat_ch[t].value for t in model.timesteps},
"Pbat_dis": {t: model.Pbat_dis[t].value for t in model.timesteps},
"Ebat": {t: model.Ebat[t].value for t in model.timesteps},
"Pgrid": {t: model.Pgrid[t].value for t in model.timesteps},
}
return solution
However, by inspecting the solution I notice that the Ebat at the first timestep is not as I have defined it (namely b_cap/2 = 8.3/2 = 4.15, but 0.415). Do you know why is this happening?
The code to run the model is the following:
import numpy as np
import pandas as pd
# battery constants
BATTERY_CAPACITY = 8.3 # Energy capacity of the battery in kWh
BATTERY_CHARGE_RATE = 2.6 # Charging/discharging rate of battery in kW
BATTERY_SOC_MAX = 1 # Maximum amount of state of charge in %
BATTERY_SOC_MIN = 0.05 # Minimum amount of state of charge in %
ETTA = 0.96 # Battery efficiency in %
# Create 'pv' DataFrame
pv = pd.DataFrame({'value': np.random.uniform(0, 5, 24)})
# Create 'demand' DataFrame
demand = pd.DataFrame({'value': np.random.uniform(0, 5, 24)})
# Create 'prices' Series
prices = pd.Series(np.random.uniform(0.1, 0.5, 24), name='prices')
optimization_results = self_consumption(pv, demand, prices)
Here is your code as posted with 2 small tweaks that should be inconsequential (used rich to format the dictionary printing and cbc for solver). Below that are the results from a run.
import pyomo.environ as pyo
import numpy as np
import pandas as pd
from rich import print
# battery constants
BATTERY_CAPACITY = 8.3 # Energy capacity of the battery in kWh
BATTERY_CHARGE_RATE = 2.6 # Charging/discharging rate of battery in kW
BATTERY_SOC_MAX = 1 # Maximum amount of state of charge in %
BATTERY_SOC_MIN = 0.05 # Minimum amount of state of charge in %
ETTA = 0.96 # Battery efficiency in %
# Create 'pv' DataFrame
pv = pd.DataFrame({'value': np.random.uniform(0, 5, 24)})
# Create 'demand' DataFrame
demand = pd.DataFrame({'value': np.random.uniform(0, 5, 24)})
# Create 'prices' Series
prices = pd.Series(np.random.uniform(0.1, 0.5, 24), name='prices')
def self_consumption(pv, demand, prices):
pv.index = np.arange(1, len(pv) + 1)
demand.index = np.arange(1, len(demand) + 1)
prices.index = np.arange(1, len(prices) + 1)
# Define the model
model = pyo.ConcreteModel()
# Define the set of timesteps
model.timesteps = pyo.Set(initialize=pyo.RangeSet(len(pv)), ordered=True)
# Define the inputs of the model
model.b_efficiency = pyo.Param(initialize=ETTA)
model.b_cap = pyo.Param(initialize=BATTERY_CAPACITY)
model.b_min_soc = pyo.Param(initialize=BATTERY_SOC_MIN)
model.b_max_soc = pyo.Param(initialize=BATTERY_SOC_MAX)
model.b_charging_rate = pyo.Param(initialize=BATTERY_CHARGE_RATE)
model.Ppv = pyo.Param(model.timesteps, initialize=pv.to_dict()['value'], within=pyo.Any)
model.Pdemand = pyo.Param(model.timesteps, initialize=demand.to_dict()['value'], within=pyo.Any)
model.day_ahead_prices = pyo.Param(model.timesteps, initialize=prices.to_dict(), within=pyo.Any)
# Define the decision variables of the model
model.Pbat_ch = pyo.Var(model.timesteps, within=pyo.NonNegativeReals, bounds=(0, model.b_charging_rate))
model.Pbat_dis = pyo.Var(model.timesteps, within=pyo.NonNegativeReals, bounds=(0, model.b_charging_rate))
model.Ebat = pyo.Var(model.timesteps, within=pyo.NonNegativeReals,
bounds=(model.b_min_soc * model.b_cap, model.b_max_soc * model.b_cap))
model.Pgrid = pyo.Var(model.timesteps)
model.is_charging = pyo.Var(model.timesteps, within=pyo.Binary)
model.AbsPgrid = pyo.Var(model.timesteps, within=pyo.NonNegativeReals)
# Define the constraints of the model
def BatEnergyRule(model, t):
if t == 1:
return model.Ebat[t] == model.b_cap / 2 # battery initialization at half of the capacity (assumption)
else:
return model.Ebat[t] == model.Ebat[t - 1] + (
model.b_efficiency * model.Pbat_ch[t] - model.Pbat_dis[t] / model.b_efficiency)
model.cons1 = pyo.Constraint(model.timesteps, rule=BatEnergyRule)
def PowerBalanceRule(model, t):
return model.Pgrid[t] == model.Pdemand[t] - model.Ppv[t] + model.Pbat_ch[t] - model.Pbat_dis[t]
model.cons2 = pyo.Constraint(model.timesteps, rule=PowerBalanceRule)
# Charge/Discharge constraint
def charge_discharge_rule(model, t):
return model.Pbat_ch[t] <= model.is_charging[t] * model.b_charging_rate
model.charge_constraint = pyo.Constraint(model.timesteps, rule=charge_discharge_rule)
def discharge_charge_rule(model, t):
return model.Pbat_dis[t] <= (1 - model.is_charging[t]) * model.b_charging_rate
model.discharge_constraint = pyo.Constraint(model.timesteps, rule=discharge_charge_rule)
def AbsPgridRule1(model, t):
return model.AbsPgrid[t] >= model.Pgrid[t]
model.cons6 = pyo.Constraint(model.timesteps, rule=AbsPgridRule1)
def AbsPgridRule2(model, t):
return model.AbsPgrid[t] >= -model.Pgrid[t]
model.cons7 = pyo.Constraint(model.timesteps, rule=AbsPgridRule2)
# Define the objective function
def ObjRule(model):
return sum(model.day_ahead_prices[t] * model.AbsPgrid[t] for t in model.timesteps)
model.obj = pyo.Objective(rule=ObjRule, sense=pyo.minimize)
# Choose a solver
opt = pyo.SolverFactory('cbc')
# Solve the optimization problem
result = opt.solve(model, tee=True)
solution = {
"Pbat_ch": {t: model.Pbat_ch[t].value for t in model.timesteps},
"Pbat_dis": {t: model.Pbat_dis[t].value for t in model.timesteps},
"Ebat": {t: model.Ebat[t].value for t in model.timesteps},
"Pgrid": {t: model.Pgrid[t].value for t in model.timesteps},
}
return solution
optimization_results = self_consumption(pv, demand, prices)
print(optimization_results)
Welcome to the CBC MILP Solver
Version: 2.10.7
Build Date: Aug 2 2023
command line - /opt/homebrew/opt/cbc/bin/cbc -printingOptions all -import /var/folders/7l/f196n6c974x3yjx5s37t69dc0000gn/T/tmpq1dkxtns.pyomo.lp -stat=1 -solve -solu /var/folders/7l/f196n6c974x3yjx5s37t69dc0000gn/T/tmpq1dkxtns.pyomo.soln (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 117 (-27) rows, 117 (-27) columns and 325 (-32) elements
Statistics for presolved model
Original problem has 24 integers (24 of which binary)
Presolved problem has 24 integers (24 of which binary)
==== 93 zero objective 25 different
==== absolute objective values 25 different
==== for integers 24 zero objective 1 different
24 variables have objective of 0
==== for integers absolute objective values 1 different
24 variables have objective of 0
===== end objective counts
Problem has 117 rows, 117 columns (24 with objective) and 325 elements
There are 1 singletons with objective
Column breakdown:
24 of type 0.0->inf, 48 of type 0.0->up, 0 of type lo->inf,
21 of type lo->up, 0 of type free, 0 of type fixed,
0 of type -inf->0.0, 0 of type -inf->up, 24 of type 0.0->1.0
Row breakdown:
20 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0,
1 of type E other, 0 of type G 0.0, 0 of type G 1.0,
0 of type G other, 24 of type L 0.0, 0 of type L 1.0,
71 of type L other, 0 of type Range 0.0->1.0, 1 of type Range other,
0 of type Free
Continuous objective value is 1.68999 - 0.00 seconds
Cgl0003I 0 fixed, 0 tightened bounds, 23 strengthened rows, 0 substitutions
Cgl0004I processed model has 117 rows, 117 columns (24 integer (24 of which binary)) and 348 elements
Cbc0038I Initial state - 5 integers unsatisfied sum - 1.53704
Cbc0038I Pass 1: suminf. 0.00000 (0) obj. 1.84167 iterations 6
Cbc0038I Solution found of 1.84167
Cbc0038I Relaxing continuous gives 1.84167
Cbc0038I Before mini branch and bound, 19 integers at bound fixed and 39 continuous
Cbc0038I Full problem 117 rows 117 columns, reduced to 25 rows 23 columns
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I Round again with cutoff of 1.8399
Cbc0038I Pass 2: suminf. 0.04772 (1) obj. 1.8399 iterations 1
Cbc0038I Pass 3: suminf. 0.08272 (1) obj. 1.8399 iterations 2
Cbc0038I Pass 4: suminf. 0.39138 (2) obj. 1.8399 iterations 3
Cbc0038I Pass 5: suminf. 1.04429 (4) obj. 1.8399 iterations 6
Cbc0038I Pass 6: suminf. 0.39138 (2) obj. 1.8399 iterations 3
Cbc0038I Pass 7: suminf. 0.04772 (1) obj. 1.8399 iterations 1
Cbc0038I Pass 8: suminf. 0.08272 (1) obj. 1.8399 iterations 5
Cbc0038I Pass 9: suminf. 0.18461 (2) obj. 1.8399 iterations 6
Cbc0038I Pass 10: suminf. 0.18461 (2) obj. 1.8399 iterations 1
Cbc0038I Pass 11: suminf. 0.38576 (3) obj. 1.8399 iterations 2
Cbc0038I Pass 12: suminf. 0.18461 (2) obj. 1.8399 iterations 2
Cbc0038I Pass 13: suminf. 0.94249 (4) obj. 1.8399 iterations 3
Cbc0038I Pass 14: suminf. 0.94249 (4) obj. 1.8399 iterations 1
Cbc0038I Pass 15: suminf. 0.94249 (4) obj. 1.8399 iterations 0
Cbc0038I Pass 16: suminf. 0.52221 (3) obj. 1.8399 iterations 1
Cbc0038I Pass 17: suminf. 0.52221 (3) obj. 1.8399 iterations 0
Cbc0038I Pass 18: suminf. 0.04772 (1) obj. 1.8399 iterations 4
Cbc0038I Pass 19: suminf. 0.08272 (1) obj. 1.8399 iterations 3
Cbc0038I Pass 20: suminf. 0.57319 (3) obj. 1.8399 iterations 6
Cbc0038I Pass 21: suminf. 0.48983 (1) obj. 1.8399 iterations 6
Cbc0038I Pass 22: suminf. 0.61653 (2) obj. 1.8399 iterations 1
Cbc0038I Pass 23: suminf. 0.98912 (3) obj. 1.8399 iterations 2
Cbc0038I Pass 24: suminf. 0.98912 (3) obj. 1.8399 iterations 0
Cbc0038I Pass 25: suminf. 0.61653 (2) obj. 1.8399 iterations 4
Cbc0038I Pass 26: suminf. 0.83138 (3) obj. 1.8399 iterations 5
Cbc0038I Pass 27: suminf. 0.48983 (1) obj. 1.8399 iterations 2
Cbc0038I Pass 28: suminf. 0.48983 (1) obj. 1.8399 iterations 0
Cbc0038I Pass 29: suminf. 0.94249 (4) obj. 1.8399 iterations 4
Cbc0038I Pass 30: suminf. 0.94249 (4) obj. 1.8399 iterations 0
Cbc0038I Pass 31: suminf. 0.60490 (3) obj. 1.8399 iterations 1
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 19 integers at bound fixed and 39 continuous
Cbc0038I Full problem 117 rows 117 columns, reduced to 25 rows 23 columns
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 1.84167 - took 0.01 seconds
Cbc0012I Integer solution of 1.8416688 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0038I Full problem 117 rows 117 columns, reduced to 68 rows 69 columns
Cbc0031I 4 added rows had average density of 18.75
Cbc0013I At root node, 19 cuts changed objective from 1.8240671 to 1.8416688 in 1 passes
Cbc0014I Cut generator 0 (Probing) - 6 row cuts average 3.0 elements, 2 column cuts (2 active) in 0.000 seconds - new frequency is 1
Cbc0014I Cut generator 1 (Gomory) - 5 row cuts average 16.6 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 3 row cuts average 2.3 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1
Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100
Cbc0014I Cut generator 6 (TwoMirCuts) - 5 row cuts average 15.8 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100
Cbc0001I Search completed - best objective 1.841668825018124, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 1.82407 to 1.84167
Probing was tried 1 times and created 8 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 1 times and created 5 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 1 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 1 times and created 5 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Result - Optimal solution found
Objective value: 1.84166883
Enumerated nodes: 0
Total iterations: 0
Time (CPU seconds): 0.01
Time (Wallclock seconds): 0.02
Total time (CPU seconds): 0.01 (Wallclock seconds): 0.02
{
'Pbat_ch': {
1: 0.0,
2: 0.72073244,
3: 1.1548358,
4: 0.0,
5: 0.0,
6: 1.4225377,
7: 1.4145129,
8: 1.1800824,
9: 1.401624,
10: 2.3134876,
11: 0.0,
12: 0.0,
13: 0.0,
14: 0.0,
15: 0.0,
16: 2.6,
17: 0.0,
18: 0.0,
19: 0.0,
20: 1.9031677,
21: 0.072007854,
22: 1.68232,
23: 1.0115796,
24: 2.1131942
},
'Pbat_dis': {
1: 0.2385651,
2: 0.0,
3: 0.0,
4: 2.6,
5: 2.2705603,
6: 0.0,
7: 0.0,
8: 0.0,
9: 0.0,
10: 0.0,
11: 2.1272945,
12: 0.51614516,
13: 0.40099184,
14: 1.5481068,
15: 0.076779037,
16: 0.0,
17: 1.7763798,
18: 0.24885847,
19: 1.9521439,
20: 0.0,
21: 0.0,
22: 0.0,
23: 0.0,
24: 0.0
},
'Ebat': {
1: 4.15,
2: 4.8419031,
3: 5.9505455,
4: 3.2422122,
5: 0.8770452,
6: 2.2426814,
7: 3.6006138,
8: 4.7334929,
9: 6.0790519,
10: 8.3,
11: 6.0840682,
12: 5.546417,
13: 5.1287172,
14: 3.516106,
15: 3.4361278,
16: 5.9321278,
17: 4.0817322,
18: 3.8225046,
19: 1.7890214,
20: 3.6160624,
21: 3.6851899,
22: 5.3002171,
23: 6.2713335,
24: 8.3
},
'Pgrid': {
1: 0.0,
2: 0.0,
3: 0.0,
4: -4.9994836,
5: -4.5241992,
6: 0.0,
7: 0.0,
8: 0.0,
9: 0.0,
10: 0.0,
11: -3.0453041,
12: 0.0,
13: 0.0,
14: 0.0,
15: 0.0,
16: -0.1398704,
17: 0.0,
18: 0.0,
19: 0.0,
20: 0.0,
21: 0.0,
22: 0.0,
23: 0.0,
24: 0.0
}
}
Process finished with exit code 0