Gekko Optimization doesn't give me a unique answer although there is a unique answer for my model

Viewed 79

my optimization model works but with different initial values for the main variable (Pre), it gives a different answer! not the optimal one! but this should have one answer. I do not understand why!

try:
    from pip import main as pipmain
except:
   from pip._internal import main as pipmain
pipmain(['install','gekko'])

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


#Initialize Model
m = GEKKO(remote=False)

#define parameter
df=pd.read_excel (r'C:\Users\....')
Pw=pd.DataFrame(df).values
eta = m.Const(value=0.6)
Pre=m.Var(lb=20, ub=30)
Pre.value=23


def f(Pw,Pre):
    Dplus=m.Var(value=0)
    Dminus=m.Var(value=0)
    for i in range(744):
       D=float(Pw[i])-Pre.value
       if D>=0:
          Dplus.value=Dplus.value+D*eta.value
       elif D<0:
          Dminus.value=Dminus.value+D
    return Dplus+Dminus

#constraint:
m.Equation(f(Pw,Pre)>=0)

#Objective:
m.Minimize(f(Pw,Pre))

#Set global options: m.options.IMODE = 2 #steady state optimization

#Solve simulation: m.solve()

1 Answers

You shouldn't use .value to build model equations because it only references the initial guess value, not the variable value. Use the m.if3() (preferred) or m.if2() functions to use conditional statements in the model. Here is an example of m.if3():

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=False)
p = m.Param()
y = m.if3(p-4,p**2,p+1)

# solve with condition<0
p.value = 3
m.solve(disp=False)
print(y.value)

# solve with condition>=0
p.value = 5
m.solve(disp=False)
print(y.value)

An even better way to incorporate conditional statements is to use slack variables. It appears that your application is for energy storage where there is inefficiency in the storage process given by eta. Here is a simple energy storage problem (see problem 4) where the loss from storage and retrieval is embedded in the optimization.

energy storage

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