List handling in GEKKO python

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).