I'm trying to solve a long block of equations from an EES implementation using the scipy.optimze.fsolve. But in this block of equations there are CoolProp calls that have a range of validation, and sometimes it yields ValueError. I want to know if there is a strategy to avoid ValueError and let fsolve try another guesses.
This is my code:
def block1(x):
def cp_gas(Ti, Tj):
return (1000/(Tj - Ti)*(x[6]*1.25 + x[1]*(0.45 *(((Tj + 273.15)/1000)
-((Ti + 273.15)/1000)) + 1.67*(((Tj + 273.15)/1000)**2 - ((Ti + 273.15)/1000)**2)/2
- 1.27*(((Tj + 273.15)/1000)**3 - ((Ti + 273.15)/1000)**3)/3
+ 0.39*(((Tj + 273.15)/1000)**4 - ((Ti + 273.15)/1000)**4)/4)
+ x[2]*(1.79 *(((Tj + 273.15)/1000) - ((Ti + 273.15)/1000))
+ 0.107*(((Tj + 273.15)/1000)**2 - ((Ti + 273.15)/1000)**2)/2
+ 0.586*(((Tj + 273.15)/1000)**3 - ((Ti + 273.15)/1000)**3)/3
- 0.2*(((Tj + 273.15)/1000)**4 -((Ti + 273.15)/1000)**4)/4)
+ x[3]*(1.11*(((Tj + 273.15)/1000) - ((Ti + 273.15)/1000))
- 0.48*(((Tj + 273.15)/1000)**2 - ((Ti + 273.15)/1000)**2)/2
+ 0.96*(((Tj + 273.15)/1000)**3 - ((Ti + 273.15)/1000)**3)/3
- 0.42*(((Tj + 273.15)/1000)**4 - ((Ti + 273.15)/1000)**4)/4)
+ x[4]*(0.88*(((Tj + 273.15)/1000) - ((Ti + 273.15)/1000))
- 0.0001*(((Tj + 273.15)/1000)**2 - ((Ti + 273.15)/1000)**2)/2
+ 0.54*(((Tj + 273.15)/1000)**3 - ((Ti + 273.15)/1000)**3)/3
- 0.33*(((Tj + 273.15)/1000)**4 - ((Ti + 273.15)/1000)**4)/4)
+ x[5]*(0.37*(((Tj + 273.15)/1000) - ((Ti + 273.15)/1000))
+ 1.05*(((Tj + 273.15)/1000)**2 - ((Ti + 273.15)/1000)**2)/2
- 0.77*(((Tj + 273.15)/1000)**3 - ((Ti + 273.15)/1000)**3)/3
+ 0.21*(((Tj + 273.15)/1000)**4 - ((Ti + 273.15)/1000)**4)/4)))
f = np.zeros(26)
# x[24] = T_out_vent
f[0] = x[0] - cp_gas(T0, Tgas5)
f[1] = m_gas_teoria_conferindo*x[8] - x[9]*0.8 - x[7]
f[2] = x[10] + x[8] - (x[7] + x[9]*0.8)
f[3] = x[9] - x[8]*Z
f[4] = x[12] + x[13] + x[14] + x[15] + x[16] + x[17] - x[11]
f[5] = x[12] - M_CO2*x[8]/x[7]
f[6] = x[13] - M_H2O*x[8]/x[7]
f[7] = x[14] - M_N2*x[8]/x[7]
f[8] = x[15] - M_O2*x[8]/x[7]
f[9] = x[16] - M_SO2*x[8]/x[7]
f[10] = x[17] - (M_Cz*x[8] - 0.8*x[9])/x[7]
f[11] = x[18] - (e*a*((1-omega_ar) + 3.76*(1-omega_ar) + omega_ar)*(MM_ar_CBG)/(MM_CBG)*x[19])
f[12] = x[1] - ((m_gas5-x[7])*FM_g_CO2+x[7]*x[12])/(x[7]*x[11]+(m_gas5-x[7])*FM_g)
f[13] = x[2] - ((m_gas5-x[7])*FM_g_H2O+x[7]*x[13])/(x[7]*x[11]+(m_gas5-x[7])*FM_g)
f[14] = x[3] - ((m_gas5-x[7])*FM_g_N2+x[7]*x[14])/(x[7]*x[11]+(m_gas5-x[7])*FM_g)
f[15] = x[4] - ((m_gas5-x[7])*FM_g_O2+x[7]*x[15])/(x[7]*x[11]+(m_gas5-x[7])*FM_g)
f[16] = x[5] - (x[7]*x[16])/(x[7]*x[11]+(m_gas5-x[7])*FM_g)
f[17] = x[6] - (x[7]*x[17])/(x[7]*x[11]+(m_gas5-x[7])*FM_g)
f[18] = x[20] - x[21]/rho_ar_in
f[19] = (1/3600)*x[21] - (x[10]+x[18])
f[20] = ((x[10]+x[18])*h_in_vent + x[22]) - (x[10] + x[18])*x[23]
f[21] = x[23] - HAPropsSI('H', 'T', x[24] + 273.15, 'P', P_out_vent*1e3, 'W', omega_ar)/1e3
f[22] = x[22] - (0.000012523*x[20] + 0.054570445)
f[23] = x[25] - HAPropsSI('C', 'T', x[24] + 273.15, 'P', P_out_vent*1e3, 'W', omega_ar)/1e3
f[24] = m_gas5 - (x[7]+x[19]+x[18])
f[25] = eta_total - ((m_gas5*x[0]*(Tgas5-T0) - (x[10]+x[18])*x[25]*(x[24]-T0))
/(x[8]*PCI_RSU + x[19]*PCI_CBG))
return f
x = fsolve(block1, np.ones(26))
The code yields ValueError depending on constant values that are previously defined.
ValueError example:
ValueError: The output for key (8) with value (-nan) is outside the range of validity: (0) to (0.94145) :: inputs were:"H","T",1.3025950731911414e+02,"P",2.0132500000000000e+05,"W",1.0890000000000000e-02
If anyone can help me I will be grateful. Thank in advance
The function you are running doesn't handle NaN values.
You can use try/except blocks to deal with it. Or change the NaN values to a 0 (or any suitable number of your choice).
Here is a toy example to help you fix your code. You have to decide what should be the correct behavior and use one of the proposed strategies to deal with NaNs.
import bumpy as np
def f(x):
if np.isnan(x):
raise ValueError('NaN is not supported')
return x**x
test_cases = [1, 2, 3, 4, np.nan, 6, 7]
print('skip in case of error')
for x in test_cases:
try:
print(f(x))
except ValueError:
pass
print()
print('fix X in case of NaN')
for x in test_cases:
if np.isnan(x):
x = 0
print(f(x))
Output:
skip in case of error
1
4
27
256
46656
823543
fix X in case of NaN
1
4
27
256
1
46656
823543