Convex Function in CVXPy does not appear to be producing the optimal solution for a battery in an electircal system
I am working on optimizing a battery system(1 MWh) that, in conjunction with solar power(200 kW nameplate capacity), aims to reduce electricity costs for a commercial building. For those unfamiliar with commercial electricity pricing, here's a brief overview: charges typically include energy charges (based on total energy consumption over a month) and demand charges (based on the peak energy usage within a 15-minute interval, referred to as peak demand). These rates vary throughout the day based on Time of Use (TOU).
The objective is to minimize monthly energy consumption to lower energy charges, while also ensuring enough energy is stored in the battery to mitigate sudden spikes in consumption to reduce demand charges. The battery should achieve this by charging when solar generation exceeds building consumption and discharging when consumption exceeds solar production. This should be straightforward with a basic load-following algorithm when sufficient solar and battery storage are available.
I tested this approach with data that successfully optimized battery operations (shown in Figure 1). However, using convex optimization resulted in significantly poorer performance (Figure 2). The optimized solution from the convex solver increased energy consumption and worsened demand charges compared to not using a battery at all. Despite optimizing for TOU rates, the solver's output falls short of an ideal solution. I have thoroughly reviewed my code, objectives, and constraints, and they appear correct to me. My hypothesis is that the solver's algorithm might prioritize sending excess power to the grid (resulting in positive peaks), potentially in an an attempt negative peaks. Maybe that is why there is a random peak on the last data point.
Figure 1: Near Ideal Battery Operation from the Load Following Algorithm
Figure 2: Battery Operation from the Convex Algorithm
Ideally, I aim to minimize both energy and demand charges, except when it's economical to store excess power for anticipated high-demand periods. Any insights or suggestions on refining this approach would be greatly appreciated.
Thank you for your assistance.
Convex Optimization CVXPy Code:
# Import libraries needed
import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
# One day of 15-minute load data
load = [ 36, 42, 40, 42, 40, 44, 42, 42, 40, 32, 32, 32, 32,
30, 34, 30, 32, 30, 32, 32, 32, 32, 32, 32, 30, 32,
32, 34, 54, 62, 66, 66, 76, 76, 80, 78, 80, 80, 82,
78, 46, 104, 78, 76, 74, 78, 82, 88, 96, 84, 94, 92,
92, 92, 92, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 86, 86, 82, 66, 72, 56, 56, 54, 48, 48, 42,
50, 42, 46, 46, 46, 42, 42, 42, 44, 44, 36, 34, 32,
34, 32, 34, 32, 32]
# One day of 15-minute solar data
solar = [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,
2, 6, 14, 26, 46, 66, 86, 104, 120, 138, 154, 168, 180,
190, 166, 152, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 190, 178, 164, 148, 132, 114, 96,
76, 58, 40, 22, 4, 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]
# Define alpha matrix which are the TOU energy charges for one day
lg = [31, 16, 25, 20, 4] # Length of each TOU period in 15 minute intervals
pk = ['off', 'mid', 'on', 'mid', 'off'] # Classifcation of each TOU period
alpha = np.array([])
for i in range(len(lg)):
if pk[i] == 'on':
mult = 0.1079
elif pk[i] == 'mid':
mult = 0.0874
elif pk[i] == 'off':
mult = 0.0755
alpha = np.append(alpha, (mult * np.ones(lg[i])))
# Define beta matricies which are the TOU demand charges for one day
val = [[0.1709, 0, 0], [0, 0.0874, 0], [0, 0, 0.0755]]
beta = {}
for i in range(len(val)):
beta_i = np.array([])
for j in range(len(lg)):
if pk[j] == 'on':
mult = val[0][i]
elif pk[j] == 'mid':
mult = val[1][i]
elif pk[j] == 'off':
mult = val[2][i]
beta_i = np.append(beta_i, (mult * np.ones(lg[j])))
beta[i] = beta_i
beta_ON = np.zeros((96, 96))
np.fill_diagonal(beta_ON, beta[0])
beta_MID = np.zeros((96, 96))
np.fill_diagonal(beta_MID, beta[1])
beta_OFF = np.zeros((96, 96))
np.fill_diagonal(beta_OFF, beta[2])
# Declare Parameters
eta_plus=0.96 # charging efficiency
eta_minus=0.96 # discharging efficiency
Emax=900 # SOC upper limit
Emin=200 # SOC lower limit
E_init=500 # initial state of charge
P_B_plus_max=200 # charging power limit
P_B_minus_max=200 # discharging power limit
opt_load=load #declaring optimal load
n=96 #declaring number of timestpes for each optimization
del_t=1/4 #time delta
d = 1 # int(len(load) / n ) # number of days
# Declare the arrays for the data outputs
pg = np.array([])
psl = np.array([])
eb = np.array([])
pbp = np.array([])
pbn = np.array([])
for i in range(d):
# Declare constraints List
cons = []
# Pull solar and load data for nth day
P_S = solar[int(n*i) : int(n*i + n)]
P_L = load[int(n*i) : int(n*i + n)]
# Declare variables
P_G = cp.Variable(n) # Power drawn from the grid at t
E_B = cp.Variable(n) # Energy in the Battery
P_B_plus = cp.Variable(n) # Battery charging power at t
P_B_minus = cp.Variable(n) # Battery discharging power at t
P_SL = cp.Variable(n) # Solar power fed to load at t
obj = cp.Minimize(cp.sum(cp.matmul(alpha, P_G) * del_t) + cp.max(cp.matmul(beta_OFF, P_G)) + cp.max(cp.matmul(beta_MID, P_G)) + cp.max(cp.matmul(beta_ON, P_G)))
for t in range(n):
# First iteration of constraints has an inital amount of energy for the battery.
if t == 0:
cons_temp = [
E_B[t] == E_init,
E_B[t] >= Emin,
E_B[t] <= Emax,
P_B_plus[t] >= 0,
P_B_plus[t] <= P_B_plus_max,
P_B_minus[t] >= 0,
P_B_minus[t] <= P_B_minus_max,
P_SL[t] + P_B_plus[t]/eta_plus == P_S[t],
P_SL[t] + P_G[t] + P_B_minus[t]*eta_minus == P_L[t],
P_SL[t] >= 0
]
# Subsequent iterations use have the amount of energy from the battery calcuated from the previous constraint
else:
cons_temp = [
E_B[t] == E_B[t - 1] + del_t*(P_B_plus[t - 1] - P_B_minus[t - 1]),
E_B[t] >= Emin,
E_B[t] <= Emax,
P_B_plus[t] >= 0,
P_B_plus[t] <= P_B_plus_max,
P_B_minus[t] >= 0,
P_B_minus[t] <= P_B_minus_max,
P_SL[t] + P_B_plus[t]/eta_plus == P_S[t],
P_SL[t] + P_G[t] + P_B_minus[t]*eta_minus == P_L[t],
P_SL[t] >= 0
]
cons += cons_temp
# Solve CVX Problem
prob = cp.Problem(obj, cons)
prob.solve(solver=cp.CBC, verbose = True, qcp = True)
# Store solution
pg = np.append(pg, P_G.value)
psl = np.append(psl, P_SL.value)
eb = np.append(eb, E_B.value)
pbp = np.append(pbp, P_B_plus.value)
pbn = np.append(pbn, P_B_minus.value)
# Update energy stored in battery for next iteration
E_init = E_B[n - 1]
# Plot Output
time = np.arange(0, 24, 0.25) # 24 hours, 15-minute intervals
plt.figure(figsize=(10, 6))
plt.plot(time, solar, label='Solar')
plt.plot(time, [i * -1 for i in load], label='Load before Optimization')
plt.plot(time, [i * -1 for i in pg], label='Load after Optimization')
plt.plot(time, pbn - pbp, label='Battery Operation')
# Adding labels and title
plt.xlabel('Time')
plt.ylabel('Demand (kW)')
plt.title('Battery Optimization Output')
# Adding legend
plt.legend()
# Display the plot
plt.grid(True)
plt.show()
Convex Optimization CVXPy Output:
===============================================================================
CVXPY
v1.3.2
===============================================================================
(CVXPY) Jun 22 03:24:36 PM: Your problem has 480 variables, 960 constraints, and 0 parameters.
(CVXPY) Jun 22 03:24:36 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Jun 22 03:24:36 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Jun 22 03:24:36 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
(CVXPY) Jun 22 03:24:36 PM: Compiling problem (target solver=CBC).
(CVXPY) Jun 22 03:24:36 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CBC
(CVXPY) Jun 22 03:24:36 PM: Applying reduction Dcp2Cone
(CVXPY) Jun 22 03:24:36 PM: Applying reduction CvxAttr2Constr
(CVXPY) Jun 22 03:24:36 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Jun 22 03:24:36 PM: Applying reduction CBC
(CVXPY) Jun 22 03:24:37 PM: Finished problem compilation (took 8.116e-01 seconds).
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Jun 22 03:24:37 PM: Invoking solver CBC to obtain a solution.
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
(CVXPY) Jun 22 03:24:37 PM: Problem status: optimal
(CVXPY) Jun 22 03:24:37 PM: Optimal value: -4.894e+01
(CVXPY) Jun 22 03:24:37 PM: Compilation took 8.116e-01 seconds
(CVXPY) Jun 22 03:24:37 PM: Solver (including time spent in interface) took 5.628e-03 seconds
Load Following Algorithm Code:
# Import libraries needed
import matplotlib.pyplot as plt
# One day of 15-minute load data
load = [ 36, 42, 40, 42, 40, 44, 42, 42, 40, 32, 32, 32, 32,
30, 34, 30, 32, 30, 32, 32, 32, 32, 32, 32, 30, 32,
32, 34, 54, 62, 66, 66, 76, 76, 80, 78, 80, 80, 82,
78, 46, 104, 78, 76, 74, 78, 82, 88, 96, 84, 94, 92,
92, 92, 92, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 86, 86, 82, 66, 72, 56, 56, 54, 48, 48, 42,
50, 42, 46, 46, 46, 42, 42, 42, 44, 44, 36, 34, 32,
34, 32, 34, 32, 32]
# One day of 15-minute solar data
solar = [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,
2, 6, 14, 26, 46, 66, 86, 104, 120, 138, 154, 168, 180,
190, 166, 152, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 190, 178, 164, 148, 132, 114, 96,
76, 58, 40, 22, 4, 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]
battery = 500
output = [] #
soc = [] # State of Charge of Battery
net_load = [] # "Optimized Load"
for i in range(96):
# With non fully charged battery and excess solar: Pull power from the solar panels to charge the batteries
if (battery < 900) and ((solar[i] - load[i]) >= 0):
# Battery can only charge up to 100 kW
if (solar[i] - load[i]) > 200:
output.append(-200)
else:
output.append(load[i] - solar[i])
# With non depleted charged battery and excessive load: Discharge the batteries and send power to the gtid
elif (battery > (200 + (load[i]/4))) and ((solar[i] - load[i]) < 0):
# Battery can only discharge up to 100 kW
if (solar[i] - load[i]) < -200:
output.append(200)
else:
output.append(load[i] - solar[i])
else:
output.append(0)
battery += (-0.25 * output[i])
soc.append(battery / 1000)
net_load.append(solar[i] - load[i] + output[i])
# Plot Output
time = np.arange(0, 24, 0.25) # 24 hours, 15-minute intervals
plt.figure(figsize=(10, 6))
plt.plot(time, solar, label='Solar')
plt.plot(time, [i * -1 for i in load], label='Load before Optimization')
plt.plot(time, net_load, label='Load after Optimization')
plt.plot(time, output, label='Battery Operation')
# Adding labels and title
plt.xlabel('Time')
plt.ylabel('Demand (kW)')
plt.title('Battery Optimization Output')
# Adding legend
plt.legend()
# Display the plot
plt.grid(True)
plt.show()
Solution
A couple caveats:
- I didn't look at the load following stuff, but you should compute the costs for that and compare with the optimization
- I think the way you are constructing the cost vector with all of the numpy/diagonals/ones/zeros/etc. is pretty confusing, but that is just an aside. I plotted the costs (
alpha) and it looks fine.
First thing: You are asking CVXPY solve to use CBC solver in your solve statement. CBC is NOT a nonlinear solver and you have max() functions in your objective, which is nonlinear, unless CVXPY is doing some enormously complex linearization under the hood (unlikely) you are getting a junk result by using that solver, unless you reformulate and get rid of the max() stuff. (Another caveat: I don't know what you're doing with the max stuff in the obj, but it isn't super relevant). So try just letting CVXPY choose by omitting that flag (as I show below)
After that, the results look ... well ... more credible. Over to you. I plotted the cost (x1000) so it would show, added the battery state, and flipped your battery usage line such that positive == charging (makes more sense to me.) The battery charges when rates are low and dumps what it has at peak rate, so that looks credible.
Your code (w/ minor tweaks mentioned):
import sys
# Import libraries needed
import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
# One day of 15-minute load data
load = [36, 42, 40, 42, 40, 44, 42, 42, 40, 32, 32, 32, 32,
30, 34, 30, 32, 30, 32, 32, 32, 32, 32, 32, 30, 32,
32, 34, 54, 62, 66, 66, 76, 76, 80, 78, 80, 80, 82,
78, 46, 104, 78, 76, 74, 78, 82, 88, 96, 84, 94, 92,
92, 92, 92, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100,
100, 100, 86, 86, 82, 66, 72, 56, 56, 54, 48, 48, 42,
50, 42, 46, 46, 46, 42, 42, 42, 44, 44, 36, 34, 32,
34, 32, 34, 32, 32]
# One day of 15-minute solar data
solar = [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,
2, 6, 14, 26, 46, 66, 86, 104, 120, 138, 154, 168, 180,
190, 166, 152, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
200, 200, 200, 200, 200, 200, 190, 178, 164, 148, 132, 114, 96,
76, 58, 40, 22, 4, 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]
# Define alpha matrix which are the TOU energy charges for one day
lg = [31, 16, 25, 20, 4] # Length of each TOU period in 15 minute intervals
pk = ['off', 'mid', 'on', 'mid', 'off'] # Classifcation of each TOU period
alpha = np.array([])
for i in range(len(lg)):
if pk[i] == 'on':
mult = 0.1079
elif pk[i] == 'mid':
mult = 0.0874
elif pk[i] == 'off':
mult = 0.0755
alpha = np.append(alpha, (mult * np.ones(lg[i])))
# Define beta matricies which are the TOU demand charges for one day
val = [[0.1709, 0, 0], [0, 0.0874, 0], [0, 0, 0.0755]]
beta = {}
for i in range(len(val)):
beta_i = np.array([])
for j in range(len(lg)):
if pk[j] == 'on':
mult = val[0][i]
elif pk[j] == 'mid':
mult = val[1][i]
elif pk[j] == 'off':
mult = val[2][i]
beta_i = np.append(beta_i, (mult * np.ones(lg[j])))
beta[i] = beta_i
beta_ON = np.zeros((96, 96))
np.fill_diagonal(beta_ON, beta[0])
beta_MID = np.zeros((96, 96))
np.fill_diagonal(beta_MID, beta[1])
beta_OFF = np.zeros((96, 96))
np.fill_diagonal(beta_OFF, beta[2])
# Declare Parameters
eta_plus = 0.96 # charging efficiency
eta_minus = 0.96 # discharging efficiency
Emax = 900 # SOC upper limit
Emin = 200 # SOC lower limit
E_init = 500 # initial state of charge
P_B_plus_max = 200 # charging power limit
P_B_minus_max = 200 # discharging power limit
opt_load = load # declaring optimal load
n = 96 # declaring number of timestpes for each optimization
del_t = 1 / 4 # time delta
d = 1 # int(len(load) / n ) # number of days
# Declare the arrays for the data outputs
pg = np.array([])
psl = np.array([])
eb = np.array([])
pbp = np.array([])
pbn = np.array([])
print(alpha)
for i in range(d):
# Declare constraints List
cons = []
# Pull solar and load data for nth day
P_S = solar[int(n * i): int(n * i + n)]
P_L = load[int(n * i): int(n * i + n)]
# Declare variables
P_G = cp.Variable(n) # Power drawn from the grid at t
E_B = cp.Variable(n) # Energy in the Battery
P_B_plus = cp.Variable(n) # Battery charging power at t
P_B_minus = cp.Variable(n) # Battery discharging power at t
P_SL = cp.Variable(n) # Solar power fed to load at t
obj = cp.Minimize(
cp.sum(cp.matmul(alpha, P_G) * del_t)
+ cp.max(cp.matmul(beta_OFF, P_G)) + cp.max(
cp.matmul(beta_MID, P_G)) + cp.max(cp.matmul(beta_ON, P_G)))
print(obj)
for t in range(n):
# First iteration of constraints has an inital amount of energy for the battery.
if t == 0:
cons_temp = [
E_B[t] == E_init,
E_B[t] >= Emin,
E_B[t] <= Emax,
P_B_plus[t] >= 0,
P_B_plus[t] <= P_B_plus_max,
P_B_minus[t] >= 0,
P_B_minus[t] <= P_B_minus_max,
P_SL[t] + P_B_plus[t] / eta_plus == P_S[t],
P_SL[t] + P_G[t] + P_B_minus[t] * eta_minus == P_L[t],
P_SL[t] >= 0
]
# Subsequent iterations use have the amount of energy from the battery calcuated from the previous constraint
else:
cons_temp = [
E_B[t] == E_B[t - 1] + del_t * (P_B_plus[t - 1] - P_B_minus[t - 1]),
E_B[t] >= Emin,
E_B[t] <= Emax,
P_B_plus[t] >= 0,
P_B_plus[t] <= P_B_plus_max,
P_B_minus[t] >= 0,
P_B_minus[t] <= P_B_minus_max,
P_SL[t] + P_B_plus[t] / eta_plus == P_S[t],
P_SL[t] + P_G[t] + P_B_minus[t] * eta_minus == P_L[t],
P_SL[t] >= 0
]
cons += cons_temp
# Solve CVX Problem
prob = cp.Problem(obj, cons)
prob.solve(verbose=True, qcp=True)
# Store solution
pg = np.append(pg, P_G.value)
psl = np.append(psl, P_SL.value)
eb = np.append(eb, E_B.value)
pbp = np.append(pbp, P_B_plus.value)
pbn = np.append(pbn, P_B_minus.value)
# Update energy stored in battery for next iteration
E_init = E_B[n - 1]
# Plot Output
time = np.arange(0, 24, 0.25) # 24 hours, 15-minute intervals
plt.figure(figsize=(10, 6))
plt.plot(time, solar, label='Solar')
plt.plot(time, [i * -1 for i in load], label='Load before Optimization')
plt.plot(time, [i * -1 for i in pg], label='Load after Optimization')
plt.plot(time, pbp - pbn, label='Battery Operation')
plt.plot(time, [t*10000 for t in alpha], label='cost[x10000]')
plt.plot(time, eb, label='Battery state')
# Adding labels and title
plt.xlabel('Time')
plt.ylabel('Demand (kW)')
plt.title('Battery Optimization Output')
# Adding legend
plt.legend()
# Display the plot
plt.grid(True)
plt.show()
- AttributeError: install_layout when attempting to install a package in a virtual environment
- Python list comprehension - want to avoid repeated evaluation
- Hash algorithm for dynamic growing/streaming data?
- matplotlib - making labels for violin plots
- Python How to I check if last element has been reached in iterator tool chain?
- Polars and the Lazy API: How to drop columns that contain only null values?
- Why are my Mean, Var, and Std outputs from NumPy different from what the online grader expects?
- Correlation dataframe convertion from results from pl.corr
- Polars DataFrame transformation
- Discord rate limiting while only sending 1 request per minute
- Check if column contains (/,-,_, *or~) and split in another column - Pandas
- How to draw a rectangle at (x,y) in a PyQt GraphicsView?
- how to calculate correlation between ten columns with polars
- How to set class attribute with await in __init__
- Detect hindi encoding, response received from Facebook API in Python
- Is it possible to write a horizontal if statement with a multi-line body?
- Max length of items in list
- Cannot subclass multiprocessing Queue in Python 3.5
- How can I get notified of updates to Python packages in a unified way?
- Using python AST to traverse code and extract return statements
- merge groups of columns in a polars dataframe to single columns
- Group Pandas DataFrame by Continuous Date Ranges
- Flask login @login_required not working
- Odoo: one2many and many2one? KeyError:'___'
- merge some columns in a Polars dataframe and duplicate the others
- Python: Create table from string mixed with separators using FOR loops
- How do I type hint a method with the type of the enclosing class?
- How can I verify an emails DKIM signature in Python?
- Writing a class that accepts a callback in Python?
- Python Paramiko channel.exec_command not returning output intermittently