Gekko differentiation with Coolprop

Question

I am trying to solve a differential problem with Gekko for the first time. I have defined the parameters with CoolProp as suggested here: Gekko and CoolProp I would like to have a nice plot of T over time.

import CoolProp.CoolProp as CP 
from gekko import GEKKO
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

m = GEKKO(remote=False) 
m.time = np.arange(0,100,1)
#variables
T = m.Var(value=78+273)
p = m.Var()
d = m.Var()
c = m.Var()
# build cubic spline
T_data = np.arange(77.5+273,98.5+273,1)
p_data = [CP.PropsSI('P','T',Ti,'Q',1,"Water") for Ti in T_data]
d_data = [CP.PropsSI('D','T',Ti,'Q',1,"Water") for Ti in T_data]
c_data = [CP.PropsSI('O','T',Ti,'Q',1,"Water") for Ti in T_data]

m.cspline(T,p,T_data,p_data)
m.cspline(T,d,T_data,d_data)
m.cspline(T,c,T_data,c_data)

#equations
def x(d,c,T):
    return (50 + 0.1/d*T)/(1000*c)

m.Equation( T.dt() == x(d,c,T) )

m.options.IMODE = 4
m.options.SOLVER = 3 # solver (IPOPT) 
#m.options.coldstart = 1

m.solve(disp=False)

fig, ax = plt.subplots()
ax.plot(m.time, y.value)
ax.set(xlabel='time (s)', ylabel='T',
       title='gekko test')
ax.grid()
plt.show()

however, I get an error

Exception: @error: Equation Definition Equation without an equality (=) or inequality (>,<) (80.0/v21)(80.5/v21)(81.0/v21)(81.5/v21)(82.0/v21)(82.5/v21) STOPPING..

The same can be solved with scipy as follows:

def x(t,T):
    p = CP.PropsSI('P','T',T,'Q',1,"Water")
    d = CP.PropsSI('D','T',T,'Q',1,"Water")
    c  = CP.PropsSI('O','T',T,'Q',1,"Water")[1.0]).c
    dTdt = (50 + 0.1/d*T)/(1000*c)
    return dTdt
y0 = [78+273]
sol = solve_ivp(x, (0,100), y0, method='LSODA', rtol=1e-6)

fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0])
ax.set(xlabel='time (s)', ylabel='T',
       title='scipy test')
ax.grid()
plt.show()

The infeasibilities.txt shows (I believe) that there is something wrong with the spline?

***** POSSIBLE INFEASBILE EQUATIONS ************

EQ Number	Lower	Residual	Upper	Infeas.	Name
1	0.0000E+00	8.0315E+02	0.0000E+00	-8.0315E+02	p(1).n(2).cspline1.cspline: y = p_0 + p_1 (x-x_i) + p_2 (x-x_i)^2 + p_3 (x-x_i)^3
Variable	Lower	Value	Upper	$Value	Name
1	-1.2346E+20	7.8000E+01	1.2346E+20	0.0000E+00	p(1).n(2).v1
3	-1.2346E+20	0.0000E+00	1.2346E+20	0.0000E+00	p(1).n(2).v3

Any idea what causes the error and how I can get the model solved successfully?

Additionally, how do I include CoolProp data with two variables in an efficient and successful way? This doesn't seem to work properly:

T_data = np.arange(77.5+273,98.5+273,1)
D_data = np.arange(3,1004,1)
p_data=[]
for i in T_data:
    row = []
    for n in D_data:
        row.append(CP.PropsSI('P','T',i,'D',n,"Water"))
    p_data.append(row)

m.bspline(T1,D1,p1,T_data,D_data,p_data)

it results in en error

Exception: @error: Bspline Config bspline in bspline1_info.csv needs 3 lines: kx,ky,sf STOPPING...

It is a test run and familiarization with Gekko for the future more complex DAE problem, which cannot be solved with scipy. Any help is much appreciated.

EDIT 22.08 Thank you, Mr. Hedengren, for the answer below! Sadly, I still get an error when I try to do the same thing with other data, such as p or h. Simplified example:

m = GEKKO(remote=False)
u1 = m.Param(value=-81e3)
D1 = m.Param(value=808)
h1 = m.Var()
D_data = np.array([807.5, 808.0, 808.5, 809.0,809.5])
u_data = np.array([-121000, -111000, -101000, -91000,-81000])
h_data = [] 
for n in D_data:
    row = []     
    for i in u_data:         
        row.append(CP.PropsSI('H','U',i,'D',n,"Nitrogen"))
    h_data.append(row) 
# Define the B-spline interpolation
m.bspline(D1, u1, h1, D_data, u_data, h_data)

# Set solver options
m.options.SOLVER = 3
m.options.DIAGLEVEL = 10

# Solve the optimization problem
m.solve()

Could someone perhaps help me understand why the solution doesn't work in this case?

John Hedengren · Accepted Answer · 2023-08-11T20:48:43.770

The infeasible solution is because d and c are initialized to a value of 0 and they appear in the denominator of the equation. Setting d.value=1 & c.value=1 or initializing with m.Var(1) resolves the problem.

import CoolProp.CoolProp as CP 
from gekko import GEKKO
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

def x(t,T):
    p = CP.PropsSI('P','T',T,'Q',1,"Water")
    d = CP.PropsSI('D','T',T,'Q',1,"Water")
    c = CP.PropsSI('O','T',T,'Q',1,"Water")
    dTdt = (50 + (0.1/d)*T)/(1000*c)
    return dTdt
y0 = [78+273]
sol = solve_ivp(x, (0,100), y0, method='LSODA', rtol=1e-6)

plt.figure(figsize=(8,3))
plt.subplot(1,2,1)
ax = plt.gca()
ax.plot(sol.t, sol.y[0])
ax.set(xlabel='time (s)', ylabel='T',
       title='scipy test')
ax.grid()


m = GEKKO(remote=False) 
m.time = np.arange(0,101,1)

#variables
T = m.Var(value=78+273)
p = m.Var(1)
d = m.Var(1)
c = m.Var(1)

# build cubic spline
T_data = np.arange(77.5+273,98.5+273,1)
#T.lower=min(T_data); T.upper=max(T_data)
p_data = [CP.PropsSI('P','T',Ti,'Q',1,"Water") for Ti in T_data]
d_data = [CP.PropsSI('D','T',Ti,'Q',1,"Water") for Ti in T_data]
c_data = [CP.PropsSI('O','T',Ti,'Q',1,"Water") for Ti in T_data]

m.cspline(T,p,T_data,p_data)
m.cspline(T,d,T_data,d_data)
m.cspline(T,c,T_data,c_data)

#equations
def x(d,c,T):
    return (50 + (0.1/d)*T)/(1000*c)

m.Equation( T.dt() == x(d,c,T) )

m.options.IMODE = 4
m.options.SOLVER = 3 # solver (IPOPT) 

m.solve(disp=False)

plt.subplot(1,2,2)
ax = plt.gca()
ax.plot(m.time, T.value)
ax.set(xlabel='time (s)',title='gekko test')
ax.grid()
plt.tight_layout()
plt.savefig('results.png',dpi=300)
plt.show()

The documentation has an example of using bspline().

from gekko import GEKKO
import numpy as np

m = GEKKO(remote=False)

xgrid = np.linspace(-1, 1, 20)
ygrid = xgrid
xg,yg = np.meshgrid(xgrid,ygrid)
z_data = xg*yg

x = m.Var(0.5,-1,1)
y = m.Var(0.5,-1,1)
z = m.Var(2)

m.bspline(x,y,z,xgrid,ygrid,z_data)

m.Minimize(z)

m.solve()

The xgrid and ygrid values are vectors while the z_data value is a matrix. In your example, the number of column characters with D_data exceeds the maximum. Use D_data as the first entry to create something with fewer column characters and more rows. The alternative is to use fewer mesh points such as D_data = np.arange(3,1004,10) to also improve the speed of building the b-spline.

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

m = GEKKO()
T_data = np.arange(77.5+273,98.5+273,1)
D_data = np.arange(3,1004,1)
p_data=[]
T1,D1 = m.Array(m.Param,2,value=100)
p1 = m.Var()
for n in D_data:
    row = []
    for i in T_data:    
        row.append(CP.PropsSI('P','T',i,'D',n,"Water"))
    p_data.append(row)

m.bspline(D1,T1,p1,D_data,T_data,p_data)

m.solve()

This solves successfully:

Number of nonzeros in equality constraint Jacobian...:        1
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        1
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 4.26e+04 0.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 0.00e+00 0.00e+00 -11.0 4.26e+04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

 The solution was found.

 The final value of the objective function is   0.000000000000000E+000

 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   5.700000000160799E-003 sec
 Objective      :   0.000000000000000E+000
 Successful solution
 ---------------------------------------------------

Hi Mr. Hedengren, thank you so much for your answer! Could you perhaps look at the edit above? The solution doesn't seem to work in that case — Beginner2020, Aug 22 '23 at 14:43
Thanks for the follow-up question. Could you create a new question and reference this one? — John Hedengren, Aug 22 '23 at 18:51

Gekko differentiation with Coolprop

1 Answers1