How does the Time Delay (m.delay) method work in Gekko

Question

I read the Time Delay method in APMonitor, and used it along with the m.Intermediate method in Gekko as given in the code below. I have two questions and appreciate your help: First) According to m.delay(B_SOC_initial, B_SOC_Prev, 1), I expect that B_SOC_Prev[1] is equal to B_SOC_initial[0], but it is equal to zero (Please see the attached image). I do not know where I am wrong. I appreciate your advice. Second) I defined a simple objective function as c[i] = Generation which maximizes the sum of c[i] values. I need to have the optimum value of the objective function, but the result shows the sum of all values, not the optimum value. The expected objective function value is 103000, but it shows 309000 (Please the attached image).In other words, the objective function gets back the the sum of c[i]'s in all iterations, not the optimum run's case. Thank you so much for your help.

from gekko import GEKKO
import numpy as np

#Battery Specifications
Battery_capacity=589200 #kW-h
Battery_Initial_capacity=117840

n=3
m = GEKKO(remote=False)
m.time = np.linspace(0, n - 1, n)

Demand =m.Param([0,46258,46952])
Generation =m.Param([0,53000,50000])

B_ch = m.Var(lb=0)
B_disch = m.Var(lb=0)

B_SOC_Prev= m.Var(value=Battery_Initial_capacity)
B_SOC_initial=m.Intermediate(B_SOC_Prev+B_ch-B_disch)
m.delay(B_SOC_initial, B_SOC_Prev, 1)

m.Equations([Generation + B_disch - B_ch == Demand,
             B_disch * B_ch <= 0,
             ]
            )

c = [None]*(n)
for i in range(n):
    c[i] = Generation

m.Maximize(sum(c))

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 1',\
                    'minlp_integer_tol 1e-8', \
                    'minlp_integer_leaves 2']
m.options.SOLVER   = 1
m.options.IMODE    = 6

# m.open_folder()
m.solve(disp=False)

print("Obj. Value=",-m.options.OBJFCNVAL)

for i in range(n):
    print("B_ch[{}]={} B_disch[{}]={} B_SOC_initial[{}]={}  B_SOC_Prev[{}]={}".format(i,B_ch[i],i,B_disch[i], i,B_SOC_initial[i],i,B_SOC_Prev[i]))

Answer 1

It looks like there is an initialization issue with using an intermediate variable and the time-delay model. Try using m.Var() types for both:

B_SOC_Prev= m.Var(value=Battery_Initial_capacity)
B_SOC_initial=m.Var(value=Battery_Initial_capacity)
m.Equation(B_SOC_initial==B_SOC_Prev+B_ch-B_disch)
m.delay(B_SOC_initial, B_SOC_Prev, 1)

The B_SOC_Prev[1] is no longer zero.

B_ch[0]=0.0 B_disch[0]=0.0 B_SOC_initial[0]=117840.0  B_SOC_Prev[0]=117840.0
B_ch[1]=6742.0 B_disch[1]=0.0 B_SOC_initial[1]=124582.0  B_SOC_Prev[1]=117840.0
B_ch[2]=3048.0 B_disch[2]=0.0 B_SOC_initial[2]=127630.0  B_SOC_Prev[2]=124582.0

Additional information on the time delay model is in the documentation.

For issue #2, use m.Maximize(Generation) instead of m.Maximize(sum(c)). This gives the correct objective function of 103000 because gekko automatically sums the values across the time horizon for that parameter. Note that it is currently a m.Param() type so the optimizer does not adjust the Generation values.

Here is a complete script:

from gekko import GEKKO
import numpy as np

#Battery Specifications
Battery_capacity=589200 #kW-h
Battery_Initial_capacity=117840

n=3
m = GEKKO(remote=False)
m.time = np.linspace(0, n - 1, n)

Demand =m.Param([0,46258,46952])
Generation =m.Param([0,53000,50000])

B_ch = m.Var(lb=0)
B_disch = m.Var(lb=0)

B_SOC_Prev= m.Var(value=Battery_Initial_capacity)
B_SOC_initial=m.Var(value=Battery_Initial_capacity)
m.Equation(B_SOC_initial==B_SOC_Prev+B_ch-B_disch)
m.delay(B_SOC_initial, B_SOC_Prev, 1)

m.Equations([Generation + B_disch - B_ch == Demand,
             B_disch * B_ch <= 0,
             ]
            )

m.Maximize(Generation)

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 1',\
                    'minlp_integer_tol 1e-8', \
                    'minlp_integer_leaves 2']
m.options.SOLVER   = 1
m.options.IMODE    = 6

# m.open_folder()
m.solve(disp=True)

print("Obj. Value=",-m.options.OBJFCNVAL)

for i in range(n):
    print(f"B_ch[{i}]={B_ch[i]} " + 
          f"B_disch[{i}]={B_disch[i]} " +
          f"B_SOC_initial[{i}]={B_SOC_initial[i]} " +
          f"B_SOC_Prev[{i}]={B_SOC_Prev[i]}")