Issue
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)
Solution
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)
Answered By - John Hedengren
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.