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()
- RuntimeError: Trying to backward through the graph a second time on loss tensor
- ERROR: post() got an unexpected keyword argument 'proxies'
- How do I automatically fix an invalid JSON string?
- Python generating a list of dates between two dates
- ModuleNotFoundError: No module named 'json'
- Check if string is None, empty or has spaces only
- invalid command name ".!canvas"
- SqlAlchemy: How to make a LONGBLOB column in mysql?
- Validate if input number is prime
- OpenAI remove key access while using AAD authentication
- Converting a Text File to a String in Python
- Problems with repeated draws of Sobol sequence in scipy
- Render NumPy array in FastAPI
- Any faster and memory-efficient alternative of torch.autograd.functional.jacobian(model.decoder, latent_l)?
- Convert integer to string in Python
- starting container process caused: exec: "fastapi": executable file not found in $PATH: unknown
- How to change version of Python in PyCharm
- convert image from CV_64F to CV_8U
- Pyspark new column when otherwise results in "should be a column" error
- FastAPI as a Windows service
- Regex to match a whole number not ending in some digits
- How can I measure the width of a string rendering via tkFont without creating a window first?
- The tilde operator in Python
- 2D stacking and optimizing
- Compute percentage of positive rows in a group_by polars DataFrame
- Extract domain name from URL in Python
- How to use AssertRaisesMessage() in Django tests
- Langchain ChatOllama always produces invalid format: expected "json" error
- Python client library for Language Server Protocol (LSP)
- How can I convert a python urandom to a string?