Python GEKKO: how to insert a function into the solver?

Question

I have quite a big project I want to solve with GEKKO. It consists of quite a large number of partial differential equations, and I have a function that uses an iterative process to calculate steady state "leak". However, GEKKO runs the function only during initialization. I want GEKKO to solve this task by taking into account that function. It would be really hard to write this function in GEKKOs equations. But without GEKKO it would be really hard to solve the Partial Differential Equations. So I am stuck, I would appreciate any help.

Here is a simple example I want to implement

a = 0.01

def CalculateLeak(P,a):
    a = a - P*0.1
    print("Leak is calculated")
    return a

m = GEKKO(remote = False)
tf = 10
nt = int(tf/1) + 1
m.time = np.linspace(0,tf,nt)
P = m.Var(0.1)
m.Equation(P.dt() == P*0.1 - CalculateLeak(P,a))
m.options.IMODE = 7
m.solve(disp = False)

print("Finished")
print(a)

Below is the function I actually want to add to my project. It calculates the adsorption amount based on P (pressure), T (temperature), y1,y2,y3,y4,y5,y6 are molar fractions. All these variables should come from the partial differential equations solved in time by GEKKO (probably except temperature which can be assumed constant some time). Every time iteration this function would calculate the amount of gas adsorbed and return to Partial Diff. Eq. as a source term.

import numpy as np 
import matplotlib.pyplot as plt


def FastIAST(P_gas,T,y1,y2,y3,y4,y5,y6):
    IP1_CH4                               = 0.0016     #kmol/kg
    IP2_CH4                               = 0          #1/bar
    IP3_CH4                               = 4.2E-05    #1/bar
    IP4_CH4                               = 2922.78    #K

    IP1_C2H6                              = 0.0027     #kmol/kg
    IP2_C2H6                              = 0.0        #1/bar 
    IP3_C2H6                              = 2.66E-04   #1/bar
    IP4_C2H6                              = 2833.77    # K

    IP1_C3H8                              = 0.0062     #kmol/kg
    IP2_C3H8                              = 0.0        #1/bar 
    IP3_C3H8                              = 3.75E-04   #1/bar
    IP4_C3H8                              = 2795.28    #K

    IP1_C4H10                             = 0.007      #kmol/kg
    IP2_C4H10                             = 0.0        #1/bar 
    IP3_C4H10                             = 0.0015     #1/bar
    IP4_C4H10                             = 2600       #K

    IP1_CO2                               = 0.0028     #kmol/kg
    IP2_CO2                               = 0.0        #(kmol/kg)/bar 
    IP3_CO2                               = 0.000748   #1/bar
    IP4_CO2                               = 2084.44    #K

    IP1_N2                                = 0.0075     #kmol/kg
    IP2_N2                                = 0.0        #(kmol/kg)/bar 
    IP3_N2                                = 0.00099    #1/bar
    IP4_N2                                = 935.77     #K
    
    Q1 = IP1_CH4 - IP2_CH4*T # Isotherm max capacity CH4
    Q2 = IP1_C2H6 - IP2_C2H6*T # Isotherm max capacity C2H6
    Q3 = IP1_C3H8 - IP2_C3H8*T # Isotherm max capacity C3H8
    Q4 = IP1_C4H10 - IP2_C4H10*T # Isotherm max capacity C4H10
    Q5 = IP1_CO2 - IP2_CO2*T # Isotherm max capacity CO2
    Q6 = IP1_N2 - IP2_N2*T  # Isotherm max capacity N2


    b1 = IP3_CH4*np.exp(IP4_CH4/T) # Isotherm affinity coeff. CH4
    b2 = IP3_C2H6*np.exp(IP4_C2H6/T) # Isotherm affinity coeff. C2H6
    b3 = IP3_C3H8*np.exp(IP4_C3H8/T) # Isotherm affinity coeff. C3H8
    b4 = IP3_C4H10*np.exp(IP4_C4H10/T) # Isotherm affinity coeff. C4H10
    b5 = IP3_CO2*np.exp(IP4_CO2/T) # Isotherm affinity coeff. CO2
    b6 = IP3_N2*np.exp(IP4_N2/T) # Isotherm affinity coeff. N2
    error = 0 # 1 - there was an error in the programm, 0 - OK
    N = 6 # Number of components
    #Langmuir Isotherm
    SingleComponentCapacity = np.array([Q1,Q2,Q3,Q4,Q5,Q6]) #Langmuir Isotherm capacity of every component
    AffinityCoefficient = np.array([b1,b2,b3,b4,b5,b6]) #Langmuir Affinity Coefficient of every component
    fractionGas = np.array([y1,y2,y3,y4,y5,y6]) #Gas fraction of every component

    #Initialization
    fastiastGraphConcentration = np.zeros(N)
    fastiastGraphFraction = np.zeros(N)
    fastiastPressure = 0
    adsorbedFraction = np.zeros(N)
    adsorbedConcentration = np.zeros(N)

    #Checking..
    if (len(fractionGas) < N or len(SingleComponentCapacity) < N or len(AffinityCoefficient) < N):
        print("You have the incorrect number of components")
        error = 1
    if np.sum(fractionGas) < 0.95 or np.sum(fractionGas) > 1.05:
        error = 1
        print("The molar fractions sum is not equal to 1")

    ###Calculation###

    kappa_old = np.zeros(N)
    delta_kappa = np.ones(N)
    kappa = np.zeros(N)
    CmuT = 0


    partialPressureComponents = fractionGas*P_gas
    for k in range(N):
        CmuT += SingleComponentCapacity[k]*AffinityCoefficient[k]*partialPressureComponents[k]

    for k in range(N):
        kappa[k] = CmuT/(SingleComponentCapacity[k])

    i = 0
    while np.any((delta_kappa) > 1e-4):
        f = np.zeros(N)
        fDerivative = np.zeros(N)
        g = np.zeros(N)
        sigma = np.zeros(N)
        phi = np.zeros((N,N))
        phi = np.matrix(phi)

        for k in range(N):
            f[k] = SingleComponentCapacity[k]*(np.log(1+kappa[k]))
            fDerivative[k] = SingleComponentCapacity[k]*(1/(1+kappa[k]))

        for k in range(N-1):
            g[k] = f[k] - f[k+1]
        for k in range(N):
            g[N-1] += AffinityCoefficient[k]*partialPressureComponents[k]/kappa[k]
        g[N-1] = g[N-1] - 1 


        for k in range(N-1):
            phi[k,k] = fDerivative[k]
            phi[k,k+1] = -fDerivative[k+1]
        for k in range(0,N):
            phi[N-1,k] = - (AffinityCoefficient[k]*partialPressureComponents[k]/(kappa[k]**2))

        sigma = np.linalg.solve(phi, g)
        kappa_old = kappa
        kappa = kappa_old - sigma

        delta_kappa = np.abs(kappa-kappa_old)


        i += 1
        if i > 20 or np.any(kappa<0):
            print("No convergence")
            error = 1
            break
        if np.any(kappa < 0):
            print("No convergence")
            error = 1
            break
        adsorbedFraction = partialPressureComponents*AffinityCoefficient/kappa
        adsorbedConcentrationPure = SingleComponentCapacity*(kappa
                                      /(1+kappa))  
        
        C_total = 0
        for k in range(0,N):
            C_total += ( (adsorbedFraction[k]) / adsorbedConcentrationPure[k])
        C_total = 1/C_total
        adsorbedConcentration =  C_total*adsorbedFraction   

    fastiastGraphConcentration=np.vstack((fastiastGraphConcentration, adsorbedConcentration))
    fastiastGraphFraction=np.vstack((fastiastGraphFraction, adsorbedFraction))
    fastiastPressure=np.vstack((fastiastPressure, P_gas))
    if error == 0:
    ###Result###
        return(fastiastGraphConcentration[1,:])
    else:
        return(0)

 
FastIAST(2,298,1,0,0,0,0,0)

Answer 1

There is no current method to call external "black-box" functions with Gekko. One of the reasons that Gekko performs well is that it compiles the functions to byte-code as if they were written in FORTRAN or C++ and it uses automatic differentiation to provide sparse 1st and 2nd derivatives to the solvers. One work-around is to use a c-spline (1D) or b-spline (2D) to approximate the function if there are only one or two independent variables. The simple problem would qualify but the FastIAST has 8 independent variables so that approach wouldn't work. There is also the deep learning library in Gekko to approximate functions of any dimension, but it may be more difficult to control the approximation error. There are new developments coming that may allow external function calls and interfaces to other machine learning libraries that would allow function approximations. As of Gekko v1.0.4, external black-box function calls aren't possible. Python function calls are allowed such as:

from gekko import GEKKO
m = GEKKO()
def f(x,c):
    y = m.sum([(xi-c)**2 for xi in x])
    return y
x1 = m.Array(m.Var,5)
p  = 2.1
m.Minimize(f(x1,p))
m.Equation(f(x1,0)<=10)
m.solve()
print(x1)