arraysperformanceoptimizationgekko

How can I speed up my optimization with Gekko?

My program is optimizing the charging and decharging of a home battery to minimize the cost of electricity at the end of the year. The electricity usage of homes is measured each 15 minutes, so I have 96 measurement point in 1 day. I want to optimilize the charging and decharging of the battery for 2 days, so that day 1 takes the usage of day 2 into account. I wrote the following code and it works.

from gekko import GEKKO
import numpy as np
import pandas as pd
import time
import math

# ------------------------ Import and read input data ------------------------

file = r'D:\Bedrijfseconomie\MP Thuisbatterijen\Spyder - Gekko\Data Sim 1.xlsx'

data = pd.read_excel(file, sheet_name='Input', na_values='NaN')

dataRead = pd.DataFrame(data, columns= ['Timestep','Verbruik woning (kWh)','Prijs afname (€/kWh)',
                                            'Capaciteit batterij (kW)','Capaciteit batterij (kWh)',
                                            'Rendement (%)','Verbruikersprofiel'])

timestep = dataRead['Timestep'].to_numpy()                                 
usage_home = dataRead['Verbruik woning (kWh)'].to_numpy()
price = dataRead['Prijs afname (€/kWh)'].to_numpy()
cap_batt_kW = dataRead['Capaciteit batterij (kW)'].iloc[0]              
cap_batt_kWh = dataRead['Capaciteit batterij (kWh)'].iloc[0]
efficiency = dataRead['Rendement (%)'].iloc[0]
usersprofile = dataRead['Verbruikersprofiel'].iloc[0]

# ---------------------------- Optimization model ----------------------------
# Initialise model
m = GEKKO()

# Global options
m.options.SOLVER = 1

# Constants
snelheid_laden = cap_batt_kW/4
T = len(timestep)

loss_charging = m.Const(value = (1-efficiency)/2) 
max_cap_batt = m.Const(value = cap_batt_kWh)
min_cap_batt = m.Const(value = 0)
max_charge = m.Const(value = snelheid_laden)                                    # max battery can charge in 15min
max_decharge = m.Const(value = -snelheid_laden)                                 # max battery can decharge in 15min

# Parameters
dummy = np.array(np.ones([T]))

# Variables
e_batt = m.Array(m.Var, (T), lb = min_cap_batt, ub = max_cap_batt)              # energy in battery
usage_net = m.Array(m.Var, (T))                                                 # usage home & charge/decharge battery
price_paid = m.Array(m.Var, (T))                                                 # price paid each 15min
charging = m.Array(m.Var, (T), lb = max_decharge, ub = max_charge)              # amount charge/decharge each 15min

# Intermediates
e_batt[0] = m.Intermediate(charging[0])
for t in range(T):
    e_batt[t] = m.Intermediate(m.sum([charging[i]*(1-loss_charging) for i in range(t)]))
usage_net = [m.Intermediate(usage_home[t] + charging[t]) for t in range(T)]
price_paid = [m.Intermediate(usage_net[t] * price[t] / 100) for t in range(T)]
total_price = m.Intermediate(m.sum([price_paid[t] for t in range(T)]))

# Equations (constraints)
m.Equation([min_cap_batt*dummy[t] <= e_batt[t] for t in range(T)])
m.Equation([max_cap_batt*dummy[t] >= e_batt[t] for t in range(T)])
m.Equation([max_charge*dummy[t] >= charging[t] for t in range(T)])
m.Equation([max_decharge*dummy[t] <= charging[t] for t in range(T)])
m.Equation([min_cap_batt*dummy[t] <= usage_net[t] for t in range(T)])
m.Equation([(-1*charging[t]) <= (1-loss_charging)*e_batt[t] for t in range(T)])

# Objective 
m.Minimize(total_price)

# Solve problem
m.solve()

My code is running and it works but despite that it gives a Solution time of 10 seconds, the total time for it to run is around 8 minutes. Does anyone know a way I can speed it up?

Solution

  • There are a few ways to speed up the Gekko code:

    For additional diagnosis, try setting m.options.DIAGLEVEL=1 to get a detailed timing report of how long it takes to compile the model and perform each function, 1st derivative, and 2nd derivative calculation. It also gives a detailed view of the solver versus model time during the solution phase.

    Update with Data File Testing

    Thanks for sending the data file. The run directory shows that the model file is 58,682 lines long. It takes a while to compile a model that size. Here is the solution from the files you sent:

     --------- APM Model Size ------------
 Each time step contains
   Objects      :          193
   Constants    :            5
   Variables    :        20641
   Intermediates:          578
   Connections  :        18721
   Equations    :        20259
   Residuals    :        19681
 
 Number of state variables:          20641
 Number of total equations: -        19873
 Number of slack variables: -         1152
 ---------------------------------------
 Degrees of freedom       :           -384
 
 * Warning: DOF <= 0
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  3.37044E+01  5.00000E+00
    1  2.81987E+01  1.00000E-10
    2  2.81811E+01  5.22529E-12
    3  2.81811E+01  2.10942E-15
    4  2.81811E+01  2.10942E-15
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    10.5119999999879      sec
 Objective      :    28.1811214884047     
 Successful solution
 ---------------------------------------------------

    Here is a version that uses IMODE=6 instead. You define the variables and equations once and let Gekko handle the time discretization. It makes a much more efficient model because there is no unnecessary duplication of equations.

    Results

    from gekko import GEKKO
import numpy as np
import pandas as pd
import time
import math

# ------------------------ Import and read input data ------------------------

file = r'Data Sim 1.xlsx'

data = pd.read_excel(file, sheet_name='Input', na_values='NaN')

dataRead = pd.DataFrame(data, columns= ['Timestep','Verbruik woning (kWh)','Prijs afname (€/kWh)',
                                            'Capaciteit batterij (kW)','Capaciteit batterij (kWh)',
                                            'Rendement (%)','Verbruikersprofiel'])

timestep = dataRead['Timestep'].to_numpy()                                 
usage_home = dataRead['Verbruik woning (kWh)'].to_numpy()
price = dataRead['Prijs afname (€/kWh)'].to_numpy()
cap_batt_kW = dataRead['Capaciteit batterij (kW)'].iloc[0]              
cap_batt_kWh = dataRead['Capaciteit batterij (kWh)'].iloc[0]
efficiency = dataRead['Rendement (%)'].iloc[0]
usersprofile = dataRead['Verbruikersprofiel'].iloc[0]

# ---------------------------- Optimization model ----------------------------
# Initialise model
m = GEKKO()
m.open_folder()

# Global options
m.options.SOLVER = 1
m.options.IMODE  = 6

# Constants
snelheid_laden = cap_batt_kW/4
m.time = timestep

loss_charging = m.Const(value = (1-efficiency)/2) 
max_cap_batt = m.Const(value = cap_batt_kWh)
min_cap_batt = m.Const(value = 0)
max_charge = m.Const(value = snelheid_laden)                                    # max battery can charge in 15min
max_decharge = m.Const(value = -snelheid_laden)                                 # max battery can decharge in 15min

# Parameters
usage_home = m.Param(usage_home)
price = m.Param(price)

# Variables
e_batt = m.Var(value=0, lb = min_cap_batt, ub = max_cap_batt)  # energy in battery
price_paid = m.Var()                                              # price paid each 15min
charging = m.Var(lb = max_decharge, ub = max_charge)              # amount charge/decharge each 15min
usage_net = m.Var(lb=min_cap_batt)

# Equations
m.Equation(e_batt==m.integral(charging*(1-loss_charging)))
m.Equation(usage_net==usage_home + charging)
price_paid = m.Intermediate(usage_net * price / 100)
m.Equation(-charging <= (1-loss_charging)*e_batt)

# Objective 
m.Minimize(price_paid)

# Solve problem
m.solve()

import matplotlib.pyplot as plt
plt.plot(m.time,e_batt.value,label='Battery Charge')
plt.plot(m.time,charging.value,label='Charging')
plt.plot(m.time,price_paid.value,label='Price')
plt.plot(m.time,usage_net.value,label='Net Usage')
plt.xlabel('Time'); plt.grid(); plt.legend(); plt.show()

    The model is only 31 lines long (see gk0_model.apm) and it solves much faster (a couple seconds total).

     --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            5
   Variables    :            8
   Intermediates:            1
   Connections  :            0
   Equations    :            6
   Residuals    :            5
 
 Number of state variables:           1337
 Number of total equations: -          955
 Number of slack variables: -          191
 ---------------------------------------
 Degrees of freedom       :            191
 
 ----------------------------------------------
 Dynamic Control with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  3.46205E+01  3.00000E-01
    1  3.30649E+01  4.41141E-10
    2  3.12774E+01  1.98558E-11
    3  3.03148E+01  1.77636E-15
    4  2.96824E+01  3.99680E-15
    5  2.82700E+01  8.88178E-16
    6  2.82039E+01  1.77636E-15
    7  2.81334E+01  8.88178E-16
    8  2.81085E+01  1.33227E-15
    9  2.81039E+01  8.88178E-16
 
 Iter    Objective  Convergence
   10  2.81005E+01  8.88178E-16
   11  2.80999E+01  1.77636E-15
   12  2.80996E+01  8.88178E-16
   13  2.80996E+01  8.88178E-16
   14  2.80996E+01  8.88178E-16
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :   0.527499999996508      sec
 Objective      :    28.0995878585948     
 Successful solution
 ---------------------------------------------------

    There is no long compile time. Also, the solver time is reduced from 10 sec to 0.5 sec. The objective function is nearly the same (28.18 versus 28.10).

    Here is a complete version without the data file dependency (in case the data file isn't available in the future).

    from gekko import GEKKO
import numpy as np
import pandas as pd
import time
import math

# ------------------------ Import and read input data ------------------------
timestep = np.arange(1,193)                                 
usage_home = np.array([0.05,0.07,0.09,0.07,0.05,0.07,0.07,0.07,0.06,
                          0.05,0.07,0.07,0.09,0.07,0.06,0.07,0.07,
                          0.07,0.16,0.12,0.17,0.08,0.10,0.11,0.06,
                          0.06,0.06,0.06,0.06,0.07,0.07,0.07,0.08,
                          0.08,0.06,0.07,0.07,0.07,0.07,0.05,0.07,
                          0.07,0.07,0.07,0.21,0.08,0.07,0.08,0.27,
                          0.12,0.09,0.10,0.11,0.09,0.09,0.08,0.08,
                          0.12,0.15,0.08,0.10,0.08,0.10,0.09,0.10,
                          0.09,0.08,0.10,0.12,0.10,0.10,0.10,0.11,
                          0.10,0.10,0.11,0.13,0.21,0.12,0.10,0.10,
                          0.11,0.10,0.11,0.12,0.12,0.10,0.11,0.10,
                          0.10,0.10,0.11,0.10,0.10,0.09,0.08,0.12,
                          0.10,0.11,0.11,0.10,0.06,0.05,0.06,0.06,
                          0.06,0.07,0.06,0.06,0.05,0.06,0.05,0.06,
                          0.05,0.06,0.05,0.06,0.07,0.06,0.09,0.10,
                          0.10,0.22,0.08,0.06,0.05,0.06,0.08,0.08,
                          0.07,0.08,0.07,0.07,0.16,0.21,0.08,0.08,
                          0.09,0.09,0.10,0.09,0.09,0.08,0.12,0.24,
                          0.09,0.08,0.09,0.08,0.10,0.24,0.08,0.09,
                          0.09,0.08,0.08,0.07,0.06,0.05,0.06,0.07,
                          0.07,0.05,0.05,0.06,0.05,0.28,0.11,0.20,
                          0.10,0.09,0.28,0.10,0.15,0.09,0.10,0.18,
                          0.12,0.13,0.30,0.10,0.11,0.10,0.10,0.11,
                          0.10,0.21,0.10,0.10,0.12,0.10,0.08])
price = np.array([209.40,209.40,209.40,209.40,193.00,193.00,193.00,
                    193.00,182.75,182.75,182.75,182.75,161.60,161.60,
                    161.60,161.60,154.25,154.25,154.25,154.25,150.70,
                    150.70,150.70,150.70,150.85,150.85,150.85,150.85,
                    150.00,150.00,150.00,150.00,153.25,153.25,153.25,
                    153.25,153.25,153.25,153.25,153.25,151.35,151.35,
                    151.35,151.35,151.70,151.70,151.70,151.70,154.95,
                    154.95,154.95,154.95,150.20,150.20,150.20,150.20,
                    153.75,153.75,153.75,153.75,160.55,160.55,160.55,
                    160.55,179.90,179.90,179.90,179.90,202.00,202.00,
                    202.00,202.00,220.25,220.25,220.25,220.25,245.75,
                    245.75,245.75,245.75,222.90,222.90,222.90,222.90,
                    203.40,203.40,203.40,203.40,205.30,205.30,205.30,
                    205.30,192.80,192.80,192.80,192.80,177.00,177.00,
                    177.00,177.00,159.90,159.90,159.90,159.90,152.50,
                    152.50,152.50,152.50,143.95,143.95,143.95,143.95,
                    142.10,142.10,142.10,142.10,143.75,143.75,143.75,
                    143.75,170.80,170.80,170.80,170.80,210.35,210.35,
                    210.35,210.35,224.45,224.45,224.45,224.45,226.30,
                    226.30,226.30,226.30,227.85,227.85,227.85,227.85,
                    225.45,225.45,225.45,225.45,225.80,225.80,225.80,
                    225.80,224.50,224.50,224.50,224.50,220.30,220.30,
                    220.30,220.30,220.00,220.00,220.00,220.00,221.90,
                    221.90,221.90,221.90,230.25,230.25,230.25,230.25,
                    233.60,233.60,233.60,233.60,225.20,225.20,225.20,
                    225.20,179.85,179.85,179.85,179.85,171.85,171.85,
                    171.85,171.85,162.90,162.90,162.90,162.90,158.85,
                    158.85,158.85,158.85])
cap_batt_kW = 3.00             
cap_batt_kWh = 5.00
efficiency = 0.95
usersprofile = 1

# ---------------------------- Optimization model ----------------------------
# Initialise model
m = GEKKO()
#m.open_folder()

# Global options
m.options.SOLVER = 1
m.options.IMODE  = 6

# Constants
snelheid_laden = cap_batt_kW/4
m.time = timestep

loss_charging = m.Const(value = (1-efficiency)/2) 
max_cap_batt = m.Const(value = cap_batt_kWh)
min_cap_batt = m.Const(value = 0)
max_charge = m.Const(value = snelheid_laden)                                    # max battery can charge in 15min
max_decharge = m.Const(value = -snelheid_laden)                                 # max battery can decharge in 15min

# Parameters
usage_home = m.Param(usage_home)
price = m.Param(price)

# Variables
e_batt = m.Var(value=0, lb = min_cap_batt, ub = max_cap_batt)  # energy in battery
price_paid = m.Var()                                              # price paid each 15min
charging = m.Var(lb = max_decharge, ub = max_charge)              # amount charge/decharge each 15min
usage_net = m.Var(lb=min_cap_batt)

# Equations
m.Equation(e_batt==m.integral(charging*(1-loss_charging)))
m.Equation(usage_net==usage_home + charging)
price_paid = m.Intermediate(usage_net * price / 100)
m.Equation(-charging <= (1-loss_charging)*e_batt)

# Objective 
m.Minimize(price_paid)

# Solve problem
m.solve()

import matplotlib.pyplot as plt
plt.plot(m.time,e_batt.value,label='Battery Charge')
plt.plot(m.time,charging.value,label='Charging')
plt.plot(m.time,price_paid.value,label='Price')
plt.plot(m.time,usage_net.value,label='Net Usage')
plt.xlabel('Time'); plt.grid(); plt.legend(); plt.show()