L

#### Luis Enriquez-Contreras

##### Guest

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:

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:

Code:

```
===============================================================================
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:

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()
```

<p>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.</p>

<p>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.</p>

<p>Figure 1: Near Ideal Battery Operation from the Load Following Algorithm

<a href="https://i.sstatic.net/53OBhFMH.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/53OBhFMH.png" alt="" /></a>

Figure 2: Battery Operation from the Convex Algorithm

<a href="https://i.sstatic.net/UJU3cnED.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/UJU3cnED.png" alt="" /></a>

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.</p>

<p>Thank you for your assistance.</p>

<p>Convex Optimization CVXPy Code:</p>

<pre><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

*== 'on':*

mult = 0.1079

elif pk

mult = 0.1079

elif pk

*== 'mid':*

mult = 0.0874

elif pkmult = 0.0874

elif pk

*== 'off':*

mult = 0.0755

alpha = np.append(alpha, (mult * np.ones(lgmult = 0.0755

alpha = np.append(alpha, (mult * np.ones(lg

*)))*

# 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]# 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]

elif pk[j] == 'mid':

mult = val[1]elif pk[j] == 'mid':

mult = val[1]

elif pk[j] == 'off':

mult = val[2]elif pk[j] == 'off':

mult = val[2]

beta_i = np.append(beta_i, (mult * np.ones(lg[j])))

betabeta_i = np.append(beta_i, (mult * np.ones(lg[j])))

beta

*= 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 # Power drawn from the grid at t

E_B = cp.Variable # Energy in the Battery

P_B_plus = cp.Variable # Battery charging power at t

P_B_minus = cp.Variable # Battery discharging power at t

P_SL = cp.Variable # 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:

# 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()

</code></pre>

<p>Convex Optimization CVXPy Output:</p>

<pre><code>===============================================================================

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

</code></pre>

<p>Load Following Algorithm Code:</p>

<pre><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 ((solarbeta_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 # Power drawn from the grid at t

E_B = cp.Variable # Energy in the Battery

P_B_plus = cp.Variable # Battery charging power at t

P_B_minus = cp.Variable # Battery discharging power at t

P_SL = cp.Variable # 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:

# 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()

</code></pre>

<p>Convex Optimization CVXPy Output:</p>

<pre><code>===============================================================================

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

</code></pre>

<p>Load Following Algorithm Code:</p>

<pre><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

*- load**) >= 0):*

# Battery can only charge up to 100 kW

if (solar# Battery can only charge up to 100 kW

if (solar

*- load**) > 200:*

output.append(-200)

else:

output.append(loadoutput.append(-200)

else:

output.append(load

*- solar**)*

# With non depleted charged battery and excessive load: Discharge the batteries and send power to the gtid

elif (battery > (200 + (load# With non depleted charged battery and excessive load: Discharge the batteries and send power to the gtid

elif (battery > (200 + (load

*/4))) and ((solar**- load**) < 0):*

# Battery can only discharge up to 100 kW

if (solar# Battery can only discharge up to 100 kW

if (solar

*- load**) < -200:*

output.append(200)

else:

output.append(loadoutput.append(200)

else:

output.append(load

*- solar**)*

else:

output.append(0)

battery += (-0.25 * outputelse:

output.append(0)

battery += (-0.25 * output

*)*

soc.append(battery / 1000)

net_load.append(solarsoc.append(battery / 1000)

net_load.append(solar

*- load**+ output**)*

# 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()

</code></pre># 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()

</code></pre>