Adding constraint to limit battery cycling to 1 charge and 1 discharge every 24 hours

Question

This enquiry is an extension to the question found in : '@Error: Solution not found' being returned when using gekko for optimization. "ind_1" and "ind_2" are lists of length 8760 containing 0s/1s. Certain hours of the year may earn additional revenue, so these indicator lists are used to distinguish those hours (further used in the maximization function I am trying to build onto this model by limiting the battery cycle to at MOST 1 charge and discharge every 24 hours. As an initial simplistic approach, I am attempting to sum up the battery command signals for each 24 hour segment and limiting it to at most 8000 kWh. You can find my approach below:

m = Gekko(remote=False)
#variables
e_battery = m.Var(lb=0, ub=4000, value=2000) #energy in battery at time t, battery size 4 MWh, initial value is 2MWh
command = m.Var(lb=-1000, ub=1000) #command power -1 to 1 (in MW)
e_price = m.Param(value = price) #price is a list of 8760 values
ind_1 = m.Param(value = ind_1) 
ind_2 = m.Param(value = ind_2)
peak_list = m.Param(value = peak_load_list) #list of the monthly peaks (an array of length 8760)
load_list = m.Param(value = load) #hourly electric load
m.time = np.linspace(0,8759, 8760)
m.Equation(e_battery.dt() == command) 

#The next 2 constraints are to ensure that the new load (original load + battery operation) is greater than 0, but less than the peak load for that month
m.Equation(load_list + command >= 0)
m.Equation(load_list + command <= peak_list)

#Here is the code to limit the cycling. "abs(command)" is used since "command" can be negative (discharge) or positive (charge), and a full charge and full discharge will equate to 8000 kWh. 
daily_sum=0
for i in range(8760):
    daily_sum += abs(command)
    if i%24==0 and i!=0: #when i=0, it's the beginning of the first day so we can skip it
        m.Equation(daily_sum <= 8000)
        daily_sum = 0 #reset to 0 in preparation for the first hour of the next day
m.Maximize((-command)*(e_price + ind_1*ind1_price + ind_2*ind2_price))
m.options.IMODE = 6
m.solve()

When adding the cycling constraint, the following output is returned:

 --------- APM Model Size ------------
 Each time step contains
   Objects      :  0
   Constants    :  0
   Variables    :  373
   Intermediates:  0
   Connections  :  0
   Equations    :  368
   Residuals    :  368
Error: At line 1545 of file apm.f90
Traceback: not available, compile with -ftrace=frame or -ftrace=full
Fortran runtime error: Out of memory

Does this particular implementation work using gekko's framework? Would I have to initialize a different type of variable for "command"? Also, I haven't been able to find many relevant examples of using for loops for the equations, so I'm very aware that my implementation might be well off. Would love to hear anyone's thoughts and/or suggestions, thanks.

Answer 1

Binary variables indicate when a destination has been reached (e_battery>3999 or e_battery<1). Integrating those binary variables gives an indication of how many times in a day the limit has been reached. One possible solution is to limit the integral of the binary variable to be less than the day count.

Below are two examples with soft constraints and hard constraints. The number of time points is reduced from 8760 to 120 (5 days) for testing.

from gekko import Gekko
import numpy as np
m = Gekko(remote=False)

n = 120 # hours
price=np.ones(n)
e_battery = m.Var(lb=0, ub=4000, value=2000) #energy in battery at time t
# battery size 4 MWh, initial value is 2MWh
command = m.Var(lb=-1000, ub=1000) #command power -1 to 1 (in MW)
e_price = m.Param(value = price) #price is a list of 8760 values
ind_1=1; ind_2=1
ind1_price=1; ind2_price=1
ind_1 = m.Param(value = ind_1) 
ind_2 = m.Param(value = ind_2)
m.time = np.linspace(0,n-1,n)
m.Equation(e_battery.dt() == command)

day = 24
discharge = m.Intermediate(m.integral(m.if3(e_battery+1,1,0)))
charge    = m.Intermediate(m.integral(m.if3(e_battery-3999,0,1)))
x         = np.ones_like(m.time)
for i in range(1,n):
    if i%day==0:
        x[i] = x[i-1] + 1
    else:
        x[i] = x[i-1]
limit = m.Param(x)

soft_constraints = True
if soft_constraints:
    derr = m.CV(value=0)
    m.Equation(derr==limit-discharge)
    derr.STATUS = 1
    derr.SPHI = 1; derr.WSPHI = 1000
    derr.SPLO = 0; derr.WSPLO = 1000

    cerr = m.CV(value=0)
    m.Equation(cerr==limit-charge)
    cerr.STATUS = 1
    cerr.SPHI = 1; cerr.WSPHI = 1000
    cerr.SPLO = 0; cerr.WSPLO = 1000
else:
    # Hard Constraints
    m.Equation(charge<=limit)
    m.Equation(charge>=limit-1)

    m.Equation(discharge<=limit)
    m.Equation(discharge>=limit-1)

m.Minimize(command*(e_price + ind_1*ind1_price + ind_2*ind2_price))
m.options.IMODE = 6
m.solve()

import matplotlib.pyplot as plt
plt.figure(figsize=(10,6))
plt.subplot(3,1,1)
plt.plot(m.time,limit.value,'g-.',label='Limit')
plt.plot(m.time,discharge.value,'b:',label='Discharge')
plt.plot(m.time,charge.value,'r--',label='Charge')
plt.legend(); plt.xlabel('Time'); plt.ylabel('Cycles'); plt.grid()
plt.subplot(3,1,2)
plt.plot(m.time,command.value,'k-',label='command')
plt.legend(); plt.xlabel('Time'); plt.ylabel('Command'); plt.grid()
plt.subplot(3,1,3)
plt.plot(m.time,e_battery.value,'g-.',label='Battery Charge')
plt.legend(); plt.xlabel('Time'); plt.ylabel('Battery'); plt.grid()
plt.show()

The application in the original question runs out of memory because 8760 equations are each simultaneously integrated over 8760 time steps. Try posing equations that are written once but valid over the entire frame. The current objective function is to minimize electricity usage. You may need to include constraints or an objective function to meet demand. Otherwise, the solution is to never use electricity because it is minimized (e.g. maximize(-command)). Here are similar Grid Energy Benchmark problems that may help.