Как избежать ошибки значения в fsolve из scipy.optimize - PullRequest
Новая
       114

Как избежать ошибки значения в fsolve из scipy.optimize

0 голосов
/ 20 июня 2020

Я пытаюсь решить длинный блок уравнений из реализации EES, используя scipy.optimze.fsolve. Но в этом блоке уравнений есть вызовы CoolProp, которые имеют диапазон проверки, и иногда он дает ValueError. Я хочу знать, есть ли стратегия избежать ValueError и позволить fsolve попробовать другие догадки.

Это мой код: 

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))

Код возвращает ValueError в зависимости от значений констант, которые были определены ранее.

Пример ValueError: 

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

Если кто поможет мне буду признателен. Заранее благодарю

1 Ответ

0 голосов
/ 20 июня 2020

Функция, которую вы выполняете, не обрабатывает значения NaN.

Вы можете использовать блоки try / except, чтобы справиться с этим. Или измените значения NaN на 0 (или любое подходящее число по вашему выбору).

Вот игрушечный пример, который поможет вам исправить ваш код. Вы должны решить, какое поведение должно быть правильным, и использовать одну из предложенных стратегий для работы с NaN. 

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))

Вывод: 

skip in case of error
1
4
27
256
46656
823543

fix X in case of NaN
1
4
27
256
1
46656
823543
...