GEKKO optimizer @TypeError: 'GK_Operators' object is not callable

Question

I am doing weather data simulation and have four equations in four variables. There are lot of other parameters but their value is provided. I am using GEKKO optimizer and getting this error. I am new to programming. Please help. I have tried to provide as much information as possible.

# all libraries
from scipy.integrate import quad
import numpy as np
# all input parameters are as below
h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 #  thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)
# calculation of incoming water vapor flux
Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
    return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)
def f2(y):
    return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)
from gekko import GEKKO
m = GEKKO()
s = m.Var(value = 0.5)
qm = m.Var(value = 1)
tg = m.Var(value = 250)
thetam = m.Var(value = 350)
m.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))-(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))==0)
m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)-LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
m.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)*int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
m.solve()

Answer 1

The error points to a specific equation in the model.

m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*\
          (ph/ps)**(rd/cp)))**(1/7)*sigma*\
          (thetam*(ph/ps)**(rd/cp))**(4))(1-(0.75*((2/3*qm*ps/g*\
          (1-(ph/ps)**(3/2)))**(1/7))))+(sigma*(thetam**(4))*m*A*\
          ((2/3*qm*ps/g)**(m))*int2)-(sigma*(tg**(4)))-((c1*uz\
          + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp)\
          -LHV*((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
          *((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)

Traceback (most recent call last):
  File "C:\Users\johnh\Desktop\test.py", line 59, in <module>
    m.Equation(RS*(1-(0.17-0.085*s))+(1.24*((qm*ph)/(0.622*thetam*\
TypeError: 'GK_Operators' object is not callable

The equation is missing a missing at least 1 operator such as on line 3 of that equation that should include a missing operator between the middle )( parenthesis.

          (thetam*(ph/ps)**(rd/cp))**(4))*(1-(0.75*((2/3*qm*ps/g*\

One other issue is that the variable m=1/7 is redefined as a gekko model with m=GEKKO(). You can use a different name for the model such as gk=GEKKO() to avoid conflict.

A good way to troubleshoot the model and make it more readable is to break down the individual pieces of the model into sub-equations such as by using an Intermediate variable with gk.Intermediate(). Another option is to define the pieces outside the equation definition. I've shown an example below:

m.Equation(a+b*c+d**e)

Simplified:

p1 = m.Intermediate(b*c) # define sub-equation with Intermediate
p2 = d**e                # define sub-equation without Intermediate
m.Equation(a+p1+p2)

The advantage of using m.Intermediate() is that the model generally solves faster for larger models. Using this approach revealed another problem. The integration with quad needs to return just the first element. By default the error is also returned.

int2 = quad(f2, yh, 1)[0]

With all of these fixes, the model is now complete:

# all libraries
from scipy.integrate import quad
import numpy as np
# all input parameters are as below
h = 1000 # height of air slab in meters
ps = 100000 # surface presuure in pascals
ph = 88000 # pressure at slab top in pascals
atop = 0.2 # entrainment parameter
b = 0.3 # moistening parameter
zh = 0.2 # hydrologically active soil depth in meters
zt = 0.4 #  thermally active soil depth in m
n = 0.25 # soil porosity
vmin = 0.5 # Mineral volume fraction of soil
vorg = 0.25 # Organic volume fraction of soil
c = 1 # Exponent in evaporation efficiency
eta = 1 # Coefficient in runoff ratio
r = 2 # Exponent in runoff ratio
uz = 4 # Mixed layer wind speed in m/s
l = 500000 # Length scale of region
qin = 8/1000 #Effective humidity of incoming air
g = 9.81 # gravitational accelaration in m/s**2
rd = 287.058 # gas constant for dry air in SI uits
cp = 1005 # dry air specific heat at constant pressure in SI units
sigma = 5.670374419*10**(-8)
A = 0.75 # integration term
m = 1/7 # integration term
c1 = 0.001 # coefficient for forced velocity
c2 = 0.00025 # coefficient for buoyancy velocity
LHV = 22.5*10**5 # latent heat of vaporization of water in SI units
rv = 461 # gas consatant for water vapor in SI units
rho = 1.225 # air density in SI units
RS = 1000 # incoming solar radiation in SI units
yh = (ph/ps)**(3/2)
# calculation of incoming water vapor flux
Qin = ((ps-ph)/g)/l*uz*qin
def f1(y):
    return (y**(8/3*rd/cp))*((y-yh)**(m-1))
int1 = quad(f1, yh, 1)[0]
def f2(y):
    return (y**(8/3*rd/cp))*((1-y)**(m-1))
int2 = quad(f2, yh, 1)[0]
from gekko import GEKKO
gk = GEKKO(remote=False)
s = gk.Var(value = 0.5)
qm = gk.Var(value = 1)
tg = gk.Var(value = 250)
thetam = gk.Var(value = 350)
gk.Equation(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
           *((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))\
           -(((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
           *((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))))*\
           eta*(s**r))-((s**c)*((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))\
           *((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm)*rho))==0)
gk.Equation(Qin-(((ps-ph)/g)/l*uz*qm)+((s**c)*((c1*uz + c2*((g/thetam*h*\
          (thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))\
          -qm)*rho))-((1-b)*(Qin+((s**c)*((c1*uz + c2*((g/thetam*h*\
          (thetam-tg))**0.5))*((0.622*611/ps*(2.718**(-(1/273.15-1/tg)\
          *LHV/rv)))-qm)*rho))))==0)
p1 = gk.Intermediate(1-(0.17-0.085*s))
p2 = gk.Intermediate(((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp))))
p3 = gk.Intermediate((1.24*p2**(1/7)*sigma*(thetam*(ph/ps)**(rd/cp))**(4)))
p4 = gk.Intermediate((1-(0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7)))))
p5 = gk.Intermediate(sigma*(thetam**(4))*m*A*((2/3*qm*ps/g)**(m))*int2)
p6 = gk.Intermediate((sigma*(tg**(4))))
p7 = gk.Intermediate((c1*uz+ c2*((g/thetam*h*(thetam-tg))**0.5)))
p8 = gk.Intermediate((p7*(tg-thetam)*rho*cp))
p9 = gk.Intermediate((c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5)))
p10 = gk.Intermediate(((0.622*611/ps*(2.718**(-(1/273.15-1/tg)*LHV/rv)))-qm))
p11 = gk.Intermediate(LHV*((s**c)*(p9*p10*rho)))
gk.Equation(RS*p1+p3*p4+p5-p6-p8-p11==0)
gk.Equation(((1.24*((qm*ph)/(0.622*thetam*(ph/ps)**(rd/cp)))**(1/7)*\
          sigma*(thetam*(ph/ps)**(rd/cp))**(4))+(sigma*(tg**(4))))*\
          0.75*((2/3*qm*ps/g*(1-(ph/ps)**(3/2)))**(1/7))-sigma*(thetam**(4))\
          *m*A*((2/3*qm*ps/g)**(m))*int2-sigma*thetam**(4)*m*A*(2/3*qm*ps/g)**(m)\
          *int1+1.2*(c1*uz + c2*((g/thetam*h*(thetam-tg))**0.5))*(tg-thetam)*rho*cp==0)
gk.options.SOLVER = 1
gk.open_folder()
gk.solve()

However, the solver is unable to find a solution to the problem:

 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
 
 Iter    Objective  Convergence
    0  2.99819E-17  3.88868E-01
    1  4.00000E-20  1.35463E-01
    2  4.00000E-20  7.81250E-03
    3  4.00000E-20  0.00000E+00
 No feasible solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.0291 sec
 Objective      :  0.
 Unsuccessful with error code  0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found

The run folder opens with gk.open_folder(). The infeasibilities.txt file shows that all of the equations are evaluating to NaN. That means that at least one of the variables is diverging to infinity. I recommend that you check the model, possibly include upper and lower bounds, or other troubleshooting methods.