i'm currently working on a model of a distillation flask for a university project, the phisical problem is described by a DAE system, and i'm trying to solve it using GEKKO.
I'm facing a problem with list handling:
In this case i built a function that outputs the compressibility factor of a mixture, and it requires as inputs 3 gekko variables T1,x,y (x,y arrays)
zv1 = m.Param(value=ZCALC(n,comps,R0,p,T1.value,x,y))
m = GEKKO()
y = m.Array(m.Var,n,value=0.)
x = m.Array(m.Var,n,value=0.)
for i in range(n):
y[i].value = y0[i]
x[i].value = x0[i]
T1 = m.Var(value=3.31513478e+02, lb=300, ub=900)
If i leave the 3 values as they are i recieve some errors like:
File "F:\Codice_GEKKO\D86_GEKKO.py", line 113, in <module>
zv1 = m.Param(value=ZCALC(n,comps,R0,p,T1.value,x,y))
File "F:\Codice_GEKKO\compressibilityfactor.py", line 48, in ZCALC
zv=np.max(np.roots([1,-1,(Av-Bv-Bv**2),-Av*Bv]))
File "<__array_function__ internals>", line 6, in roots
File "C:\Users\verci\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\lib\polynomial.py", line 222, in roots
non_zero = NX.nonzero(NX.ravel(p))[0]
File "<__array_function__ internals>", line 6, in nonzero
File "C:\Users\verci\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\core\fromnumeric.py", line 1908, in nonzero
return _wrapfunc(a, 'nonzero')
File "C:\Users\verci\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\core\fromnumeric.py", line 67, in _wrapfunc
return _wrapit(obj, method, *args, **kwds)
File "C:\Users\verci\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\core\fromnumeric.py", line 44, in _wrapit
result = getattr(asarray(obj), method)(*args, **kwds)
File "C:\Users\verci\AppData\Local\Programs\Python\Python36\lib\site-packages\gekko\gk_operators.py", line 25, in __len__
return len(self.value)
File "C:\Users\verci\AppData\Local\Programs\Python\Python36\lib\site-packages\gekko\gk_operators.py", line 144, in __len__
return len(self.value)
TypeError: object of type 'int' has no len()```
Traceback (most recent call last):
File "F:\Codice_GEKKO\D86_GEKKO.py", line 113, in <module>
zv1 = m.Param(value=ZCALC(n,comps,R0,p,T1,x,y))
File "F:\Codice_GEKKO\compressibilityfactor.py", line 27, in ZCALC
(1-np.sqrt(t/tc[ii])))**2
TypeError: loop of ufunc does not support argument 0 of type GK_Operators which has no callable sqrt method
The first error is biven by the fact that x and y are not list, but they are GEKKO arrays and the second error is due to T1 not being a float (t=T1)
I found out that by using T1.value i can avoid the second error but still i have the first error
I have read the gekko documentation but i haven't been able to find a method to obtain a "standard" python list from a GEKKO array
Thank you in advance for your help
There are two different methods for obtained the value of zv.
Option 1: Initialization Calculation
The first method is to use floating point numbers to obtain a single calculation that can be used for initialization of a parameter. This first method allows any type of functions such as np.roots() or np.sqrt(). The function ZCALC() returns a floating point number. Even though Gekko variables are used as an input, the floating point number is accessed from a scalar variable with T1.value or from an array variable with x[i].value.
def ZCALC(n,comps,R0,p,T1,x,y):
# using the initialized values
t = T1.value
i = 0 # select values from x,y arrays
x1 = x[i].value
y1 = y[i].value
print('t,x[0],y[0] initialized values')
print(t,x1,y1)
# include equations for compressibility factor
w = (1-np.sqrt(t/300))**2
z = np.max(np.roots([1,-1,x1**2,x1*y1]))
# original equations from question
#(1-np.sqrt(t/tc[ii])))**2
#zv=np.max(np.roots([1,-1,(Av-Bv-Bv**2),-Av*Bv]))
return z
zv1 = m.Param(value=ZCALC(n,comps,R0,p,T1,x,y))
Option 2: Implicit Calculation
If the compressibility factor needs to change as T1 and x,y change then use Gekko variables so that the model is compiled with that dependency. The functions are only called during problem initialization. Gekko needs the equations with specific Gekko functions to enable automatic differentiation to provide gradients to the solvers.
def ZCALC2(n,comps,R0,p,T1,x,y):
# using gekko variables
t = T1
i = 0
x1 = x[i] # use index to x array
y1 = y[i] # use index to y array
# use Gekko equations, not Numpy
w = (x1/y1)*(1-m.sqrt(t/300))**2
# set lower bound to get the maximum root
zv = m.Var(value=ZCALC(n,comps,R0,p,T1,x,y),lb=10)
# solve for roots of eq with gekko, not with np.roots
eq = 1-zv+x1**2*zv+x1*y1*zv**3
m.Equation(eq==0)
return zv
zv2 = ZCALC2(n,comps,R0,p,T1,x,y)
Here is a script that shows the two methods:
import numpy as np
m=GEKKO(remote=False)
def ZCALC(n,comps,R0,p,T1,x,y):
# using the initialized values
t = T1.value
i = 0 # select values from x,y arrays
x1 = x[i].value
y1 = y[i].value
print('t,x[0],y[0] initialized values')
print(t,x1,y1)
# include equations for compressibility factor
w = (1-np.sqrt(t/300))**2
z = np.max(np.roots([1,-1,x1**2,x1*y1]))
# original equations from question
#(1-np.sqrt(t/tc[ii])))**2
#zv=np.max(np.roots([1,-1,(Av-Bv-Bv**2),-Av*Bv]))
return z
def ZCALC2(n,comps,R0,p,T1,x,y):
# using gekko variables
t = T1
i = 0
x1 = x[i] # use index to x array
y1 = y[i] # use index to y array
# use Gekko equations, not Numpy
w = (x1/y1)*(1-m.sqrt(t/300))**2
# set lower bound to get the maximum root
zv = m.Var(value=ZCALC(n,comps,R0,p,T1,x,y),lb=10)
# solve for roots of eq with gekko, not with np.roots
eq = 1-zv+x1**2*zv+x1*y1*zv**3
m.Equation(eq==0)
return zv
n = 3
y = m.Array(m.Var,n)
x = m.Array(m.Var,n)
x0 = [0.1,0.2,0.3]
y0 = [0.15,0.25,0.35]
for i in range(n):
y[i].value = y0[i]
x[i].value = x0[i]
T1 = m.Var(value=331, lb=300, ub=900)
comps = ['C2=','C3=','C2H8']; R0 = 8.314; p=10
# define Zv from initialized values (fixed parameter)
zv1 = m.Param(value=ZCALC(n,comps,R0,p,T1,x,y))
# define Zv from Gekko variables (updates with T1,x,y changes)
zv2 = ZCALC2(n,comps,R0,p,T1,x,y)
# initialized value of zv1 does not update with changes in T1,x,y
# initialized value of zv2 does update with changes in T1,x,y
print('initialized value of zv1, zv2')
print(zv1.value,zv2.value)
If the compressibility factor correlations can't be expressed as Gekko equations then try the cspline for 1D or bspline for 2D functions to create an approximation. You may be able to use the bspline function if compressibility can depend on just 2 variables T and x (replace y with an explicit calculation of x).