时间戳

首页 > 解决方案 > 求解四个非线性方程组

optimization - 求解四个非线性方程组

问题描述

我正在尝试使用模块中的 fsolve 求解四个非线性方程组scipy.optimize。雅可比矩阵 ( fjac) 用 填充nan。我不明白我做错了什么。我一直在寻找类似的问题，但没有一个与我的问题相似，这使我得出结论，我的尝试存在根本性错误。

我的代码：

import numpy as np
from scipy.optimize import fsolve
import math

#Constants
Cn = 1.0
Cf = 0.8
Nn = 3.0
gamma = 1.4
R = 287.1 #J/kg*K
dn = 0.5 * 25.4 #mm
df = 2.1516 * 25.4 #mm
dv = 1.8553 * 25.4 #mm
dbt = 0.89 * 25.4 #mm
Tb = (40 + 273.15) #K
Pr = (13 + 1.01325) * 10 ** 5 #N/m2

guess_massflow = 0.3 # kg/s
guess_Pc = 12 * 10 ** 5 #N/m2
guess_Pb = 10 * 10 ** 5 #N/m2
guess_Ps = 14 * 10 ** 5 #N/m2

def f(p):

    massflow, Pb, Pc, Ps = p

    def pi_mod(d):
        return (math.pi * d ** 2) / 4 * 10 ** (-6)
    major_mod = Nn * Cn * Pb / (Tb ** 0.5)
    radical = (gamma + 1)/(gamma - 1)
    minor_mod = (gamma / R * (2 / (gamma + 1)) ** radical) ** 0.5

    radical2 = (gamma - 1) / gamma

    def main_mod(P1, P2):
        return P2 * (1/R/Tb * 2/radical2 * (P1/P2) ** (2/gamma) * (1 - (P1/P2) ** radical2)) ** 0.5

    f1 = massflow - major_mod * pi_mod(dn) * minor_mod
    f2 = massflow - Cf * pi_mod(df) * main_mod(Pb,Pc)
    f3 = massflow - Cf * pi_mod(dbt) * main_mod(Pc,Ps)
    f4 = massflow - Cf * pi_mod(dv) * main_mod(Ps,Pr)

    return (f1, f2, f3, f4)

solution = fsolve(f, (guess_massflow, guess_Pb, guess_Pc, guess_Ps), full_output = True)

输出：

(array([  3.00000000e-01,   1.10000000e+06,   1.10000000e+06,
     1.10000000e+06]), {'nfev': 19, 'fjac': array([[ nan,  nan,  nan,  nan],
   [ nan,  nan,  nan,  nan],
   [ nan,  nan,  nan,  nan],
   [ nan,  nan,  nan,  nan]]), 'r': array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan]), 'qtf': array([ nan,  nan,  nan,  nan]), 'fvec': array([-0.65463805,  0.3       ,  0.3       , -3.45205928])}, 5, 'The iteration is not making good progress, as measured by the \n  improvement from the last ten iterations.')

我试图用 Mathcad 解决这个问题。解决方案收敛到，这是现实的：

[0.717255, 8.264713*10^5, 8.344225*10^5, 1.392793*10^6]

标签： optimizationscipy

解决方案

我认为问题出在求解器调用main_mod(P1, P2)时P1>P2，所以在这种情况下我强制函数返回零......我不知道它是否在物理上正确......但它似乎工作： 

import numpy as np
from scipy.optimize import fsolve
import math

#Constants
Cn = 1.0
Cf = 0.8
Nn = 3.0
gamma = 1.4
R = 287.1 # J/kg/K  
dn = 0.5 * 25.4 # mm
df = 2.1516 * 25.4 # mm
dv = 1.8553 * 25.4 # mm
dbt = 0.89 * 25.4 # mm
Tb = (40 + 273.15) # K
Pr = (13 + 1.01325) * 1e5 # N/m2

def pi_mod(d):
    return (math.pi * d ** 2) / 4e6  # <- mm2 to m2?

radical = (gamma + 1)/(gamma - 1)
radical2 = (gamma - 1) / gamma
minor_mod = np.sqrt( gamma/R*( 2/(gamma + 1) )**radical )

def main_mod(P1, P2):
    if P1 < P2:
        return P2 * np.sqrt( 1/R/Tb * 2/radical2 * (P1/P2)**(2/gamma) * (1 - (P1/P2)**radical2) )
    else:
        return 0

def f(p):
    massflow, Pb, Pc, Ps = p

    major_mod = Nn*Cn/np.sqrt(Tb) * Pb 

    f1 = massflow - major_mod * pi_mod(dn) * minor_mod
    f2 = massflow - Cf * pi_mod(df) * main_mod(Pb, Pc)
    f3 = massflow - Cf * pi_mod(dbt) * main_mod(Pc, Ps)
    f4 = massflow - Cf * pi_mod(dv) * main_mod(Ps, Pr)

    return (f1, f2, f3, f4)

guess_massflow = 0.3 # kg/s
guess_Pc = 12e5 # N/m2
guess_Pb = 10e5 # N/m2
guess_Ps = 14e5 # N/m2

p_zero = (guess_massflow, guess_Pb, guess_Pc, guess_Ps)

solution = fsolve(f, p_zero, full_output=True)

我还更改了单位转换pi_mod（它是磁盘的区域吗？）并将常量值和函数放在函数f之外

推荐阅读