Run out of memory issue in Gekko python

Question

Now, I'm trying to solve the optimization problem as above and I made it in a code as belows.

However for large N (such as 300), the code failed to work because of run out of the memory.

So, I tried m = GEKKO(remote=True) instead of m = GEKKO(remote=False), but it was not done even after 10-12 hours.

I'm now finding a way to solve this problem while using m = GEKKO(remote=False).

Code

# Import package
from gekko import GEKKO
import numpy as np

# Define parameters
P_CO = 600                      # $/tonCO
beta_CO2 = 1                    # no unit
P_CO2 = 60                      # $/tonCO2eq
E_ref = 3.1022616               # tonCO2eq/tonCO
E_dir = -1.600570692            # tonCO2eq/tonCO
E_indir_others = 0.3339226804   # tonCO2eq/tonCO
E_indir_elec_cons = 18.46607256 # GJ/tonCO
C1_CAPEX = 285695               # no unit
C2_CAPEX = 188.42               # no unit
C1_FOX = 82282                  # no unit
C2_FOX = 24.094                 # no unit
C1_ROX = 4471.5                 # no unit
C2_ROX = 96.034                 # no unit
C1_UOX = 1983.7                 # no unit
C2_UOX = 249.79                 # no unit
r = 0.08                        # discount rate
N = 300                         # number of scenarios
T = 30                          # total time period
GWP_init = 0.338723235          # 2020 Electricity GWP in EU 27 countries
theta_max = 1600000             # Max capacity

# Function to make GWP_EU matrix (TxN matrix)
def Electricity_GWP(GWP_init, n_years, num_episodes):

    GWP_mean = 0.36258224*np.exp(-0.16395611*np.arange(1, n_years+2)) + 0.03091272
    GWP_mean = GWP_mean.reshape(-1,1)
    GWP_Yearly = np.tile(GWP_mean, num_episodes) 

    noise = np.zeros((n_years+1, num_episodes))
    stdev2050 = GWP_mean[-1] * 0.25 
    stdev = np.arange(0, stdev2050 * (1 + 1/n_years), stdev2050/n_years)

    for i in range(n_years+1):
        noise[i,:] = np.random.normal(0, stdev[i], num_episodes) 

    GWP_forecast = GWP_Yearly + noise 

    return GWP_forecast

GWP_EU = Electricity_GWP(GWP_init, T, N) # (T+1)*N matrix
GWP_EU = GWP_EU[1:,:] # T*N matrix

print(np.shape(GWP_EU))

# Build Gekko model
m = GEKKO(remote=False)
theta = m.Array(m.Var, N, lb=0, ub=theta_max)
demand = np.ones((T,1))
demand[0] = 8031887.589
for k in range(1,11):
    demand[k] = demand[k-1] * 1.026 
for k in range(11,21):
    demand[k] = demand[k-1] * 1.016
for k in range(21,T):
    demand[k] = demand[k-1] * 1.011 
demand = 0.12 * demand
demand = np.tile(demand, N) # T*N matrix

print(np.shape(demand))

obj = m.sum([m.sum([((1/(1+r))**(t+1))*((P_CO*m.min3(demand[t,s], theta[s])) \
            + (beta_CO2*P_CO2*m.min3(demand[t,s], theta[s])*(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s])) \
            - (C1_CAPEX+C2_CAPEX*theta[s]+C1_FOX+C2_FOX*theta[s])-(C1_ROX+C2_ROX*m.min3(demand[t,s], theta[s])+C1_UOX+C2_UOX*m.min3(demand[t,s], theta[s]))) for t in range(T)]) for s in range(N)])
m.Maximize(obj/N)
m.solve() 

Answer 1

The multiple m.min3() expressions lead to many additional variables. Try defining it once and substituting the value m3 into the objective expression.

m3 = [[m.min3(demand[t,s],theta[s]) for t in range(T)] for s in range(N)]

obj = m.sum([sum([((1/(1+r))**(t+1))*((P_CO*m3[s][t]) \
            + (beta_CO2*P_CO2*m3[s][t]\
            *(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s])) \
            - (C1_CAPEX+C2_CAPEX*theta[s]+C1_FOX+C2_FOX*theta[s])\
            - (C1_ROX+C2_ROX*m3[s][t]\
            +C1_UOX+C2_UOX*m3[s][t])) \
                for t in range(T)]) for s in range(N)])

Below is a test script that increases the size of N from 10 to 100.

It shows the compile time and solution time for each of the cases. The number of variables for N=10 to N=100 is shown as the x-axis.

# Import package
from gekko import GEKKO
import numpy as np
import time
import matplotlib.pyplot as plt

# Define parameters
P_CO = 600                      # $/tonCO
beta_CO2 = 1                    # no unit
P_CO2 = 60                      # $/tonCO2eq
E_ref = 3.1022616               # tonCO2eq/tonCO
E_dir = -1.600570692            # tonCO2eq/tonCO
E_indir_others = 0.3339226804   # tonCO2eq/tonCO
E_indir_elec_cons = 18.46607256 # GJ/tonCO
C1_CAPEX = 285695               # no unit
C2_CAPEX = 188.42               # no unit
C1_FOX = 82282                  # no unit
C2_FOX = 24.094                 # no unit
C1_ROX = 4471.5                 # no unit
C2_ROX = 96.034                 # no unit
C1_UOX = 1983.7                 # no unit
C2_UOX = 249.79                 # no unit
r = 0.08                        # discount rate
T = 30                          # total time period
GWP_init = 0.338723235          # 2020 Electricity GWP in EU 27 countries
theta_max = 1600000             # Max capacity

# Function to make GWP_EU matrix (TxN matrix)
def Electricity_GWP(GWP_init, n_years, num_episodes):

    GWP_mean = 0.36258224*np.exp(-0.16395611*np.arange(1, n_years+2)) + 0.03091272
    GWP_mean = GWP_mean.reshape(-1,1)
    GWP_Yearly = np.tile(GWP_mean, num_episodes) 

    noise = np.zeros((n_years+1, num_episodes))
    stdev2050 = GWP_mean[-1] * 0.25 
    stdev = np.arange(0, stdev2050 * (1 + 1/n_years), stdev2050/n_years)

    for i in range(n_years+1):
        noise[i,:] = np.random.normal(0, stdev[i], num_episodes) 

    GWP_forecast = GWP_Yearly + noise 

    return GWP_forecast

Nx = np.array([10,20,30,40,60,80,100])
tv = 122*Nx+1
tt = np.zeros_like(Nx)
tsolve = np.zeros_like(Nx)
for i,N in enumerate(Nx):
    # N=number of scenarios
    print(N)
    ts = time.time()
    GWP_EU = Electricity_GWP(GWP_init, T, N) # (T+1)*N matrix
    GWP_EU = GWP_EU[1:,:] # T*N matrix

    # Build Gekko model
    m = GEKKO(remote=False)
    theta = m.Array(m.Var, N, lb=0, ub=theta_max)
    demand = np.ones((T,1))
    demand[0] = 8031887.589
    for k in range(1,11):
        demand[k] = demand[k-1] * 1.026 
    for k in range(11,21):
        demand[k] = demand[k-1] * 1.016
    for k in range(21,T):
        demand[k] = demand[k-1] * 1.011 
    demand = 0.12 * demand
    demand = np.tile(demand, N) # T*N matrix

    m3 = [[m.min3(demand[t,s],theta[s]) for t in range(T)] for s in range(N)]

    obj = m.sum([sum([((1/(1+r))**(t+1))*((P_CO*m3[s][t]) \
                + (beta_CO2*P_CO2*m3[s][t]\
                  *(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s])) \
                - (C1_CAPEX+C2_CAPEX*theta[s]+C1_FOX+C2_FOX*theta[s])\
                - (C1_ROX+C2_ROX*m3[s][t]\
                   +C1_UOX+C2_UOX*m3[s][t])) \
                        for t in range(T)]) for s in range(N)])
    m.Maximize(obj/N)
    m.solve(disp=True)
    tsolve[i] = m.options.SOLVETIME
    tt[i] = time.time()-ts-tsolve[i]
    
plt.figure(figsize=(8,5))
plt.plot(tv,tt,'ro-',label='Compile time')
plt.plot(tv,tsolve,'b--',label='Solve time')
plt.legend(); plt.grid()
plt.ylabel('Time (sec)'); plt.xlabel('Problem Size')
plt.savefig('Solution_time.png',dpi=300)
plt.show()

One other thing that improved the solution time is to use sum() instead of m.sum() for the inner summation. It is faster to solve this way, but other problems are faster with m.sum().