result using coolprop seems to be wrong

[fatsen]

When calculating with coolprop, the result seems to be wrong. The scripts below is simplified to show the problem.

from CoolProp.CoolProp import PropsSI   #coolprop's source code is C++, and it has a python wrapper.
from gekko import GEKKO

m = GEKKO(remote=False)

a = m.Var(value=1, name='a')
t = m.Var(value=20, name='t')

fluid_name = 'Water'
# m.Equation(a + PropsSI('P', 'Q', 1, 'T', t + 273.15, fluid_name) / 1000 == 0) 
# with the above equation, [TypeError: must be real number, not GK_Operators] will be raised.
m.Equation(a + PropsSI('P', 'Q', 1, 'T', t.value + 273.15, fluid_name) / 1000 == 0)
m.Equation(a + t == 0)
m.solve()

for i in m._variables:
    print(i.name, i.value)

print(a.value[0] + PropsSI('P', 'Q', 1, 'T', t.value[0] + 273.15, fluid_name) / 1000)  # should be zero
print(a.value[0]+t.value[0])  # should be zero

the outputs are:

a [-2.3393181883]
t [2.3393181883]
-1.6160227470100634
0.0

The output of [-1.616] should be 0. Other equations not using coolprop are OK. Is coolprop not compatible with gekko? Thank you very much.

[APMonitor]

This is a great question, but it would be better to post on StackOverflow. See a related question on CoolProp and gekko: https://stackoverflow.com/questions/62832981/gekko-and-coolprop

The issue you're encountering is that GEKKO cannot directly handle calls to external libraries like CoolProp within the equation-building process. GEKKO expects explicit mathematical expressions that it can compile with automatic differentiation for optimization, but CoolProp's PropsSI function doesn't provide this.

To fix this, you can decouple CoolProp calculations from the GEKKO equations by calculating the thermodynamic properties outside of the GEKKO model (before defining the equations). One workaround is to update the variables dynamically after solving the model or use a cubic spline, as demonstrated below. Here’s how you can modify your script to avoid direct CoolProp calls in GEKKO's equations:

from CoolProp.CoolProp import PropsSI
from gekko import GEKKO
import numpy as np

# Create the GEKKO model
m = GEKKO(remote=False)

# Define the variables
Tsat = m.Var(lb=273.15,ub=373.15,name='Tsat')  # Tempearture (°C)
Psat = m.Var(value=101.325, name='Psat')  # Psat (kPa)

# Define the fluid name for CoolProp
fluid_name = 'Water'

# Generate data for the cubic spline
n = 100  # Number of data points
T_vals = np.linspace(273.15, 373.15, n)  # Temperature range in °C
P_vals = [PropsSI('P', 'Q', 1, 'T', T, fluid_name) / 1000 for T in T_vals]  # Psat in kPa

print(T_vals,P_vals)

# Build the cubic spline in GEKKO
m.cspline(Tsat, Psat, T_vals, P_vals)

# Define the equations using the spline
m.Equation(Psat==101.325)

# Solve the GEKKO model
m.solve()

# Print the results
for i in m._variables:
    print(i.name, i.value)

This produces the output to calculate Tsat at 1 atm:

tsat [373.12408224]
psat [101.325]