flexibility using scipy Odeint in equation for ipopt

Question

I am in the process of validating a couple of models that I have written in pyoptsparse and scipy. I want to be able to test the same code in Gekko to see if the results are consistent with the more symbolic Gekko version that I originally made. To do this I have tried to use Odeint and integrating with my own code instead of using the built in Gekko capabilities. Is it possible to do it this way? is the issue with the optimizers ability to take the gradient of the Odeint?

import numpy as np
from scipy.integrate import odeint
from gekko import GEKKO

nhrs = 24
time = np.arange(nhrs)
m = GEKKO(remote=True)

gen = m.Array(m.Param, nhrs)
for geni in gen:
    geni.value = 1

tes_store = m.Array(m.Var, nhrs)
for tesi in tes_store:
    tesi.value = 0
    tesi.lower = -1
    tesi.upper = 1

T = m.Array(m.SV, nhrs)
for testi in T:
    testi.upper = 24
    testi.lower = 0

def tes_T_dt(t, T, gen, tes_store):
    return gen*tes_store

def tes_T(gen, tes_store, t):
    T_hist = [0]
    T_new = 0
    for i in range(len(t)):
        step = odeint(tes_T_dt, T_new, [0, 1], args=(gen[i].value, tes_store[i].value))
        T_new = step[1][0]
        T_hist.append(T_new)
    return T_hist

Temp = tes_T(gen, tes_store, time)[1:]

equations = []
m.Equation(T == tes_T(gen, tes_store, time)[1:])

m.Obj(-T[23])
Temp = tes_T(gen, tes_store, time)[1:]

m.options.SOLVER = 3
m.solve()

for i in range(len(T)):
    print(f'{i} gen={gen[i].value}, tes_store={tes_store[i].value}, T={T[i].value}, Temp={Temp[i]}')

I am getting the following error

 @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 false
 STOPPING...
Traceback (most recent call last):
  File "gekkopyoptsimpleTEST.py", line 54, in <module>
    m.solve()
  File "/home/prism/miniconda3/lib/python3.7/site-packages/gekko/gekko.py", line 2103, in solve
    raise Exception(response)
Exception:  @error: Equation Definition
 Equation without an equality (=) or inequality (>,<)
 false
 STOPPING...

Answer 1

All of the expressions in Gekko must use Gekko variables so that the model can be compiled to byte-code. If you open the run folder and gk_model0.apm, you will see the model that you've created with this script:

Model
Parameters
    p1 = 1
    p2 = 1
    p3 = 1
    p4 = 1
    ! removed p5-p21
    p22 = 1
    p23 = 1
    p24 = 1
End Parameters
Variables
    v1 = 0, <= 1, >= -1
    v2 = 0, <= 1, >= -1
    v3 = 0, <= 1, >= -1
    v4 = 0, <= 1, >= -1
    ! removed v5-v45
    v46 = 0, <= 24, >= 0
    v47 = 0, <= 24, >= 0
    v48 = 0, <= 24, >= 0
End Variables
Equations
    False
    minimize (-v48)
End Equations

End Model

As you can see the Equations - End Equations section, the Equation:

m.Equation(T == tes_T(gen, tes_store, time)[1:])

evaluated to False because it is comparing an array to another array that is not equal. You would need to have a Gekko expression there instead such as m.Equation(T[i] == tes_T(gen, tes_store, time)[i+1]). However, tes_T is not callable by gekko so the right side of the equation would just be a number in the model that would never update. A black-box function is not currently supported in gekko.

Gekko can integrate the differential equation. Is this close to what you are looking for to maximize the final value of T by adjusting gen and subject to the differential equation dT/dt==gen*tes_store?

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

nhrs = 24
m = GEKKO(remote=False)
m.time = np.arange(nhrs)
gen = m.Param(value=1)
tes_store = m.Var(value=0,lb=-1,ub=1)
T = m.SV(lb=0,ub=24)
m.Equation(T.dt()==gen*tes_store)
p = np.zeros(nhrs); p[-1]=1
final = m.Param(p)
m.Maximize(T*final)
m.options.SOLVER = 3
m.options.NODES = 2
m.options.IMODE = 6
m.solve()

plt.plot(m.time,T.value,'k-',label='T')
plt.plot(m.time,gen.value,'b.-',label='gen')
plt.plot(m.time,tes_store,'r--',label='tes_store')
plt.legend()
plt.show()