I have a question about PV+ESS cost minimize by using python GEKKO

Question

I would like to minimize the electricity bill for my home by using self-consumption PV+ESS. I would like to use the power generated through PV at home, store the excess power in ESS, and if the power bill is relatively low, I would like to charge the ESS directly instead of using it at home to minimize the daily power bill. At this point, should we use conditional statements? Currently, there is an error in my code.

from gekko import GEKKO
import numpy as np       

m = GEKKO()          

hour = 24
Num_ESS = 1
x = m.Array(m.Var,(hour,Num_ESS))
ESS_c = m.Array(m.Var,(hour,Num_ESS))
ESS_d = m.Array(m.Var,(hour,Num_ESS))
SOC_t = m.Array(m.Var,(hour,Num_ESS))

## cost
TOU = [57.7,57.7,57.7,57.7,57.7,57.7,57.7,57.7,57.7,
       108.9,131.4,131.4,108.9,131.4,131.4,131.4,131.4,108.9,108.9,108.9,
       108.9,108.9,108.9,57.7]

## ESS charging / discharging and SOC 
for tt in range(0,hour):
    ESS_c[tt,0].lower = 0; ESS_d[tt,0].lower = 0; SOC_t[tt,0].lower = 1; x[tt,0].lower = 0
    ESS_c[tt,0].upper = 2; ESS_d[tt,0].upper = 2; SOC_t[tt,0].upper = 4; x[tt,0].upper = 2

## PV data
PV_a = [0,0,0,0,0,0.5,0.7,1,1.3,1.5,2,2.5,3,1.3,1,0.7,0.5,0,0,0,0,0,0,0]

## load data
people_a = [0.3,0.3,0.4,0.5,0.7,1.18,1.4,1.6,1.6,1.3,0.8,0.8,1.1,1.7,1.7,1.5,1.1,0.9,0.4,0.2,0.2,0.1,0.1,0.1]

#ESS parameter
ess_ini = 1
ess_max = 4
ess_min = 0
power = 3
C_ess = 0.95
D_ess = 0.95

# charing to ESS

eq_charging = np.zeros((hour,1))
eq_charging = list(eq_charging)

for tt in range(0,hour):
    eq_charging[tt] = PV_a[tt] == ESS_c[tt,0]

m.Equations(eq_charging)

## Supply/Demand Match
eq_total = np.zeros((hour,1))
eq_total = list(eq_total)

for tt in range(0,hour):
    eq_total[tt] = (PV_a[tt] + ESS_d[tt,0] + x[tt,0]) == people_a[tt] +ESS_c[tt,0]

m.Equations(eq_total)

## SCO 
SOC_t[0,0] = 1

eq_ESS_SOC = np.zeros((hour))
eq_ESS_SOC = list(eq_ESS_SOC)

for tt in range(0,hour):
    eq_ESS_SOC[tt] = SOC_t[tt-1,0] + (ESS_c[tt,0]*C_ess -ESS_d[tt,0]*D_ess) == SOC_t[tt,0]
    
m.Equations(eq_ESS_SOC)


# object function
F = np.zeros((hour*Num_ESS))
F= F.tolist()

for tt in range(0,hour):
    for i in range(0,Num_ESS):
        F[i+tt*Num_ESS] = TOU[tt]*x[tt,0]

    F_Obj = np.sum(F)

m.options.IMODE = 3
m.options.SOLVER = 1
m.solve(disp=False)

apm 203.249.127.36_gk_model0 <br><pre> ----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 Warning: there is insufficient data in CSV file 203.249.127.36_gk_model0.csv
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :           24
   Constants    :            0
   Variables    :          696
   Intermediates:            0
   Connections  :          600
   Equations    :          649
   Residuals    :          649
 
 Number of state variables:            696
 Number of total equations: -          672
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :             24
 
 ----------------------------------------------

 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[11], line 2
      1 m.Obj(F_Obj)
----> 2 m.solve(disp=True)

File c:\Users\GCU\Desktop\test\test_1\Lib\site-packages\gekko\gekko.py:2185, in GEKKO.solve(self, disp, debug, GUI, **kwargs)
   2183 #print APM error message and die
   2184 if (debug >= 1) and ('@error' in response):
-> 2185     raise Exception(response)
   2187 #load results
   2188 def byte2str(byte):

Exception:  @error: Solution Not Found

Answer 1

The infeasible solution is because of the constraints on the variables:

## ESS charging / discharging and SOC 
for tt in range(0,hour):
    ESS_c[tt,0].lower = 0; ESS_d[tt,0].lower = 0
    SOC_t[tt,0].lower = 1; x[tt,0].lower = 0
    ESS_c[tt,0].upper = 2; ESS_d[tt,0].upper = 2
    SOC_t[tt,0].upper = 4; x[tt,0].upper = 2

It appears that there is a problem with the lower bound and upper bound(s). Widening the limits allow the problem to solve successfully.

## ESS charging / discharging and SOC 
for tt in range(0,hour):
    ESS_c[tt,0].lower = 0; ESS_d[tt,0].lower = 0
    SOC_t[tt,0].lower = 0.8; x[tt,0].lower = 0
    relax = 1.7
    ESS_c[tt,0].upper = 2+relax; ESS_d[tt,0].upper = 2+relax
    SOC_t[tt,0].upper = 4+relax; x[tt,0].upper = 2+relax

Here is a complete version that solves successfully:

from gekko import GEKKO
import numpy as np       

m = GEKKO()          

hour = 24
Num_ESS = 1
x = m.Array(m.Var,(hour,Num_ESS))
ESS_c = m.Array(m.Var,(hour,Num_ESS))
ESS_d = m.Array(m.Var,(hour,Num_ESS))
SOC_t = m.Array(m.Var,(hour,Num_ESS))

## cost
TOU = [57.7,57.7,57.7,57.7,57.7,57.7,57.7,57.7,57.7,
       108.9,131.4,131.4,108.9,131.4,
       131.4,131.4,131.4,108.9,108.9,108.9,
       108.9,108.9,108.9,57.7]

## ESS charging / discharging and SOC 
for tt in range(0,hour):
    ESS_c[tt,0].lower = 0; ESS_d[tt,0].lower = 0
    SOC_t[tt,0].lower = 0.8; x[tt,0].lower = 0
    relax = 1.7
    ESS_c[tt,0].upper = 2+relax; ESS_d[tt,0].upper = 2+relax
    SOC_t[tt,0].upper = 4+relax; x[tt,0].upper = 2+relax

## PV data
PV_a = [0,0,0,0,0,0.5,0.7,1,1.3,1.5,
        2,2.5,3,1.3,1,0.7,0.5,0,0,0,0,0,0,0]

## load data
people_a = [0.3,0.3,0.4,0.5,0.7,
            1.18,1.4,1.6,1.6,1.3,0.8,
            0.8,1.1,1.7,1.7,1.5,1.1,
            0.9,0.4,0.2,0.2,0.1,0.1,0.1]

#ESS parameter
ess_ini = 1
ess_max = 4
ess_min = 0
power = 3
C_ess = 0.95
D_ess = 0.95

# charing to ESS

eq_charging = np.zeros((hour,1))
eq_charging = list(eq_charging)

for tt in range(0,hour):
    eq_charging[tt] = PV_a[tt] == ESS_c[tt,0]

m.Equations(eq_charging)

## Supply/Demand Match
eq_total = np.zeros((hour,1))
eq_total = list(eq_total)

for tt in range(0,hour):
    eq_total[tt] = (PV_a[tt] + ESS_d[tt,0] + x[tt,0]) \
                     == people_a[tt] +ESS_c[tt,0]

m.Equations(eq_total)

## SCO 
SOC_t[0,0] = 1

eq_ESS_SOC = np.zeros((hour))
eq_ESS_SOC = list(eq_ESS_SOC)

for tt in range(0,hour):
    eq_ESS_SOC[tt] = SOC_t[tt-1,0] + (ESS_c[tt,0]*C_ess 
                     -ESS_d[tt,0]*D_ess) == SOC_t[tt,0]
    
m.Equations(eq_ESS_SOC)

# object function
F = np.zeros((hour*Num_ESS))
F= F.tolist()

for tt in range(0,hour):
    for i in range(0,Num_ESS):
        F[i+tt*Num_ESS] = TOU[tt]*x[tt,0]

    F_Obj = np.sum(F)

m.options.IMODE = 3
m.options.SOLVER = 1
m.solve(disp=True)

 Number of state variables:             96
 Number of total equations: -           72
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :             24
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  3.94326E-14  3.00000E+00
    1  4.44089E-26  4.44089E-16
    2  4.44089E-26  4.44089E-16
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.490000000194414E-002 sec
 Objective      :   0.000000000000000E+000
 Successful solution
 ---------------------------------------------------

I recommend a visualization of the solution to determine the points where the problem becomes infeasible. If it isn't possible to widen the bounds, set soft constraints on the outputs (penalty in the objective function) so that only the least important constraints are satisfied. Here is an example for explicit ranking of the constraints.

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO()
m.time = np.linspace(0,10,101)

# Dynamic control options
m.options.IMODE = 6
m.options.CV_TYPE = 1
m.options.MV_TYPE = 0
m.options.SOLVER = 3
m.options.MV_STEP_HOR = 1
m.options.NODES = 3

#Define Manipulated Variables
u = m.MV(name='u')

#Define Controled Variables
y = m.CV(1,name='y')
z = m.CV(1,name='z')
s = m.CV(1,name='s')

# Environmental Constraint
#setup CV
# tau is the speed of the CV response, 0=step, 1 = 63.2# of the way
#   to the new setpoint in 1 sec, only if tr_init is 1 or 2.
# with tr_init=0, it is just a pure dead-band
# specifying the speed to get to the set point
# get to 63.2# of sp withing tau seconds
y.TAU = 5
y.STATUS = 1
y.TR_INIT = 2
y.SPHI = 5
y.SPLO = 4
y.FSTATUS = 0
y.WSPHI = 100
y.WSPLO = 100

# Operational Constraint
z.TAU = 4
z.STATUS = 1
z.TR_INIT = 2
z.SPHI = 7
z.SPLO = 6
z.FSTATUS = 0
z.WSPHI = 50
z.WSPLO = 50

# Safety Constraint
s.TAU = 10
s.STATUS = 1
s.TR_INIT = 2
s.TR_OPEN = 3
s.SPHI = 11
s.SPLO = 10
s.FSTATUS = 0
s.WSPHI = 200
s.WSPLO = 200

#setup MV (u)
u.STATUS = 1
u.DCOST = 0
u.LOWER = 0
u.UPPER = 1000
u.COST = 0

# process model
tau = 1
K = 3
m.Equation(tau*y.dt()+y==u)
m.Equation(z==y)
m.Equation(s==y)

# solve problem
m.solve(disp=True)

# get additional solution information
import json
with open(m.path+'//results.json') as f:
    results = json.load(f)

# create plot
p, ax = plt.subplots(nrows=2, ncols=1, \
                     gridspec_kw={'height_ratios':[3,1]})

ax[0].plot(m.time,results['s.tr_hi'],'r-.',lw=2)
ax[0].plot(m.time,results['y.tr_hi'],'b:',lw=2)
ax[0].plot(m.time,results['z.tr_hi'],'--',color='orange',lw=2)
ax[0].plot(m.time,results['z'],'k-',lw=3)
ax[0].legend(['Priority 1: Safety Constraint',\
            'Priority 2: Environmental Constraint',\
            'Priority 3: Economic Constraint','Response'],loc=4)
ax[0].plot(m.time,results['z.tr_lo'],'--',color='orange',lw=2)
ax[0].plot(m.time,results['y.tr_lo'],'b:',lw=2)
ax[0].plot(m.time,results['s.tr_lo'],'r-.',lw=2)
ax[0].set_ylabel('Pressure (bar)')

ax[1].plot(m.time,u.value,'b-',lw=2)
ax[1].legend(['Manipulated Variable'])
ax[1].set_ylabel('MV')
ax[1].set_xlabel('Time (min)')

plt.show()