Gekko - Pumped Hydro electric problem - python

Question

I am trying to model a pumped hydro electric problem using GEKKO and python, and I am running into some errors trying to avoid charging and discharging at the same time. The idea is using price as input to Pump/Release water and obtain the highest possible benefit, having a maximum storage of 4 hours of water release. My problem is that I can't seem to work out binary variables for the problem to work:

import numpy as np
from gekko import GEKKO

horas = 4
año = 2018
Eff = 0.75
PMD_calculo = np.array([1,1.5,3,0.5,0.2])
#PMD_calculo = np.array(PMD_año[año]['Precio'])
Num_horas = len(PMD_calculo)

m = GEKKO(remote=False)
m.options.SOLVER = 1
m.solver_options = ['minlp_integer_tol 0.0']
#variables

Carga = m.Var(Num_horas, lb=0, ub=horas) # Storage (from 0 to 4 hours of energy)

Tur = m.Array(m.Var, Num_horas, lb=0, ub=1)  # Generation
Bom = m.Array(m.Var, Num_horas, lb=0, ub=1/Eff)  # Pumping

T = m.Array(m.Var, Num_horas, lb=0, ub=1, integer = True)
B = m.Array(m.Var, Num_horas, lb=0, ub=1, integer = True)

m.Equation(Carga[i] <= Bom[i-1]*B[i-1] - Tur[i-1]*T[i-1] for i in range(1,len(PMD_calculo) ))# charge / discharge
m.Equation(B[i]+T[i] <= 1 for i in range(len(PMD_calculo))) # Avoid charge and discharge

m.Maximize(Tur*T*PMD_calculo - Bom*PMD_calculo*B)

m.solve()

With this I get the following error:

----------------------------------------------------------------
 APMonitor, Version 1.0.0
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 @error: Inequality Definition
 invalid inequalities: z > x < y
 <generatorobject<genexpr>at0x000001b29a296eb0>
 STOPPING . . .

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-45-26c2eb591bd7> in <module>
     21 m.Maximize(Tur*T*PMD_calculo - Bom*PMD_calculo*B)
     22 
---> 23 m.solve()

C:\ProgramData\Anaconda3\lib\site-packages\gekko\gekko.py in solve(self, disp, debug, GUI, **kwargs)
   2138                 print("Error:", errs)
   2139             if (debug >= 1) and record_error:
-> 2140                 raise Exception(apm_error)
   2141 
   2142         else: #solve on APM server

Exception: @error: Inequality Definition
 invalid inequalities: z > x < y
 <generatorobject<genexpr>at0x000001b29a296eb0>
 STOPPING . . .

I THINK that i am getting the correct output when not using binary variables (charge-discharge either cancel each other out in the same hour or in multiple hours) but it still makes the process hard to process

Answer 1

The syntax problem is with the equation and objective definitions. Use Equations instead of Equation when there is a list of equations. Use m.sum() when there are multiple objective terms.

# charge / discharge
m.Equations([Carga <= Bom[i-1]*B[i-1] - Tur[i-1]*T[i-1] 
             for i in range(1,Num_horas)])

# Avoid charge and discharge
m.Equations([B[i]+T[i] <= 1 
             for i in range(Num_horas)])

m.Maximize(m.sum(Tur*T*PMD_calculo - Bom*PMD_calculo*B))

The complete script is:

import numpy as np
from gekko import GEKKO

horas = 4
año = 2018
Eff = 0.75
PMD_calculo = np.array([1,1.5,3,0.5,0.2])
Num_horas = len(PMD_calculo)

m = GEKKO(remote=False)
m.options.SOLVER = 1
m.solver_options = ['minlp_integer_tol 0.0']

# Storage (from 0 to 4 hours of energy)
Carga = m.Var(Num_horas, lb=0, ub=horas)

Tur = m.Array(m.Var, Num_horas, lb=0, ub=1)  # Generation
Bom = m.Array(m.Var, Num_horas, lb=0, ub=1/Eff)  # Pumping

T = m.Array(m.Var, Num_horas, lb=0, ub=1, integer = True)
B = m.Array(m.Var, Num_horas, lb=0, ub=1, integer = True)

# charge / discharge
m.Equations([Carga <= Bom[i-1]*B[i-1] - Tur[i-1]*T[i-1] 
             for i in range(1,Num_horas)])

# Avoid charge and discharge
m.Equations([B[i]+T[i] <= 1 
             for i in range(Num_horas)])

m.Maximize(m.sum(Tur*T*PMD_calculo - Bom*PMD_calculo*B))

m.solve()

 Number of state variables:             30
 Number of total equations: -            9
 Number of slack variables: -            9
 ---------------------------------------
 Degrees of freedom       :             12
 
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:      0.00 NLPi:    5 Dpth:    0 Lvs:    3 Obj: -6.67E-02 Gap:       NaN
--Integer Solution:   0.00E+00 Lowest Leaf:  -6.67E-02 Gap:   6.67E-02
Iter:     2 I:  0 Tm:      0.00 NLPi:    2 Dpth:    1 Lvs:    2 Obj:  0.00E+00 Gap:  6.67E-02
Iter:     3 I:  0 Tm:      0.00 NLPi:    2 Dpth:    1 Lvs:    4 Obj: -2.00E-01 Gap:  6.67E-02
--Integer Solution:  -2.00E-01 Lowest Leaf:  -2.00E-01 Gap:   0.00E+00
Iter:     4 I:  0 Tm:      0.00 NLPi:    2 Dpth:    2 Lvs:    4 Obj: -2.00E-01 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   1.459999999497086E-002 sec
 Objective      :  -0.200000000000000     
 Successful solution
 ---------------------------------------------------

Alternative to Binary Variables

Another way to solve this problem is with slack variables instead of binary variables. See Energy Benchmarks with storage and retrieval (Benchmark IV) from the paper:

Gates, N.S., Hill, D.C., Billings, B.W., Powell, K.M., Hedengren, J.D., Benchmarks for Grid Energy Management with Python Gekko, 60th IEEE Conference on Decision and Control (CDC), Tutorial Session: Open Source Software for Control Applications, Austin, TX, Dec 13-15, 2021.

The fourth benchmark problem models a hybrid system with a single generator with constant production constraints coupled with energy storage that together must meet an oscillating electricity demand. The goal of the problem is to minimize the required power production and use energy storage to capture excess generation serve the oscillating energy demand while keeping the base-load generator production constant. As before, the model has perfect foresight. In order to prevent the energy storage from charging and discharging simultaneously without requiring mixed-integer variables, slack variables are used to control when the storage charges and discharges, allowing it to switch modes in a way that is both continuous and differentiable. This allows the modeling language to use automatic differentiation to provide exact gradients to the solver.

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

m = GEKKO(remote=False)
m.time = np.linspace(0,1,101)

g = m.FV(); g.STATUS = 1 # production
s = m.Var(1e-2, lb=0)    # storage inventory
store = m.Var()          # store energy rate
s_in = m.Var(lb=0)       # store slack variable
recover = m.Var()        # recover energy rate
s_out = m.Var(lb=0)         # recover slack variable
eta = 0.7
d = m.Param(-2*np.sin(2*np.pi*m.time)+10)
m.periodic(s)
m.Equations([g + recover/eta - store >= d,
             g - d == s_out - s_in,
             store == g - d + s_in,
             recover == d - g + s_out,
             s.dt() == store - recover/eta,
             store * recover <= 0])
m.Minimize(g)

m.options.SOLVER   = 1
m.options.IMODE    = 6
m.options.NODES    = 3
m.solve()

plt.figure(figsize=(6,3))
plt.subplot(2,1,1)
plt.plot(m.time,d,'r-',label='Demand')
plt.plot(m.time,g,'b:',label='Prod')
plt.legend(); plt.grid(); plt.xlim([0,1])

plt.subplot(2,1,2)
plt.plot(m.time,s,'k-',label='Storage')
plt.plot(m.time,store,'g--', label='Store Rate')
plt.plot(m.time,recover,'b:', label='Recover Rate')
plt.legend(); plt.grid(); plt.xlim([0,1])
plt.show()