Gekko - Pumped Hydro electric problem - python

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