Sympy nsolve 变量作为初始猜测答案

【问题标题】：Sympy nsolve variables as initial guessSympy nsolve 变量作为初始猜测
【发布时间】：2021-11-25 10:01:06
【问题描述】：

我正在使用 Sympy nsolve 来求解非线性方程组。求解器运行良好，但我需要初始猜测的灵活性，因此我希望可以插入变量作为初始猜测。

例如，对于固定的初始猜测，一切正常：

import numpy as np
from sympy import Symbol, nsolve, cos, sin, asin, acos, S

# Definition of parameters #
tracking_angle_deg = 55
pv_area = 2.174  # pv panel area in m^2
mirror_area = pv_area / cos(0.25 * np.pi)  # mirrors area in m^2
pv_length = 1  # pv panel length in m
mirror_length = pv_length / cos(0.25 * np.pi)  # mirror length in m\
k = 0.027  # thermal conduction of air
pv_eps = 0.9  # pv emissivity
mirror_eps = 0.8  # mirrors emissivity
gravity = 9.81
sigma = 5.67e-8  # boltzman constant
nu = 15.69e-6  # kinematic viscosity of air
alpha = 22.9e-6  # thermal diffusivity of air
Pr = nu / alpha  # Prandtl of air
F_m_sky = 0.5  # shape factor mirror-sky
F_pv_sky = 1 - ((2 ** 0.5) / 2)  # shape factor pv-sky
F_pv_m = 1 - F_pv_sky  # shape factor mirror-pv
e1 = (1 - pv_eps) / (pv_eps * pv_area)
e2 = (1 - mirror_eps) / (mirror_eps * mirror_area)
f12 = 1 / (pv_area * F_pv_m)
f13 = 1 / (pv_area * F_pv_sky)
f23 = 1 / (mirror_area * F_m_sky)
Epower = 240  # electric power output from PV
Ta_kelvin = 297
Tsky = Ta_kelvin - 5
PV_flux = 535.5
m1_flux = 86.1  # absorbed flux on left/west mirror
m2_flux = 86.1  # absorbed flux on right/east mirror
teta_deg = 23
teta_rad = (teta_deg / 180) * np.pi

# angles of mirrors
phi1_rad = ((teta_deg + 45) / 180) * np.pi
phi2_rad = ((abs(teta_deg - 45)) / 180) * np.pi

J_PV_1 = Symbol('J_PV_1', real=True)
J_PV_2 = Symbol('J_PV_2', real=True)
J_m1 = Symbol('J_m1', real=True)
J_m2 = Symbol('J_m2', real=True)
Tm1 = Symbol('Tm1', real=True)
Tm2 = Symbol('Tm2', real=True)
Tpv = Symbol('Tpv', real=True)

E1 = sigma * (Tpv ** 4)
Em1 = sigma * (Tm1 ** 4)
Em2 = sigma * (Tm2 ** 4)
E3 = sigma * (Tsky ** 4)
J3 = E3

# Natural convection PV #
beta_pv = 1 / ((Tpv + Ta_kelvin) / 2)
Ra_pv = (gravity * beta_pv * (Tpv - Ta_kelvin) * (pv_length ** 3)) / (nu * alpha)
Ra_pv_crit = (2e9) / np.exp(2.94 * abs(teta_rad))
Nusselt_pv_down = 0.59 * ((Ra_pv * cos(teta_rad)) ** 0.25)
h_pv_down = (Nusselt_pv_down * k) / pv_length
Nusselt_pv_up = (0.1 + ((0.05 * abs(teta_rad)) / np.pi)) * (Ra_pv ** (1 / 3))  # turbulent flow correlation
h_pv_up = (Nusselt_pv_up * k) / pv_length

# mirror1 left west #
beta_m1 = 1 / ((Tm1 + Ta_kelvin) / 2)
Ra_m1 = abs((gravity * beta_m1 * (Tm1 - Ta_kelvin) * (mirror_length ** 3)) / (nu * alpha))
Nusselt_m1_down = 0.59 * ((Ra_m1 * cos(phi1_rad)) ** 0.25)
Nusselt_m1_up = (0.1 + ((0.05 * abs(phi1_rad)) / np.pi)) * (Ra_m1 ** (1 / 3))
h_m1_down = (Nusselt_m1_down * k) / mirror_length
h_m1_up = (Nusselt_m1_up * k) / mirror_length
h_m1 = h_m1_down + h_m1_up

# mirror2 right east #
beta_m2 = 1 / ((Tm2 + Ta_kelvin) / 2)
Ra_m2 = abs((gravity * beta_m2 * (Tm2 - Ta_kelvin) * (mirror_length ** 3)) / (nu * alpha))
Nusselt_m2_down = 0.59 * ((Ra_m2 * cos(phi2_rad)) ** 0.25)
Nusselt_m2_up = (0.1 + ((0.05 * abs(phi2_rad)) / np.pi)) * (Ra_m2 ** (1 / 3))
h_m2_down = (Nusselt_m2_down * k) / mirror_length
h_m2_up = (Nusselt_m2_up * k) / mirror_length
h_m2 = h_m2_down + h_m2_up

# equations #
eq1 = ((E1 - J_PV_1) / e1) + ((J_m1 - J_PV_1) / f12) + ((J3 - J_PV_1) / f13)
eq2 = ((E1 - J_PV_2) / e1) + ((J_m2 - J_PV_2) / f12) + ((J3 - J_PV_2) / f13)
eq3 = ((Em1 - J_m1) / e2) + ((J_PV_1 - J_m1) / f12) + ((J3 - J_m1) / f23)
eq4 = ((Em2 - J_m2) / e2) + ((J_PV_2 - J_m2) / f12) + ((J3 - J_m2) / f23)
eq5 = (PV_flux * pv_area) - Epower - (h_pv_up * pv_area * (Tpv - Ta_kelvin)) - (
        h_pv_down * pv_area * (Tpv - Ta_kelvin)) - ((E1 - J_PV_1) / e1) - ((E1 - J_PV_2) / e1)
eq6 = (m1_flux * mirror_area) - (h_m1 * mirror_area * (Tm1 - Ta_kelvin)) - (
        mirror_eps * sigma * mirror_area * ((Tm1 ** 4) - (Tsky ** 4))) - ((Em1 - J_m1) / e2)
eq7 = (m2_flux * mirror_area) - (h_m2 * mirror_area * (Tm2 - Ta_kelvin)) - (
        mirror_eps * sigma * mirror_area * ((Tm2 ** 4) - (Tsky ** 4))) - ((Em2 - J_m2) / e2)

initial_guess = [PV_flux, PV_flux, PV_flux, PV_flux, Ta_kelvin, Ta_kelvin, Ta_kelvin + 15]

print(nsolve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], [J_PV_1, J_PV_2, J_m1, J_m2, Tm1, Tm2, Tpv],
             [500, 500, 500, 500, 293, 293, 310], domain=S.Reals))

输出：

Matrix([[593.171007328550], [593.142501153849], [491.462389303786], [491.059251115046], [304.687760815851], [304.609734541627], [321.670012318818]])

但是我不能以任何我尝试的方式在初始猜测中使用变量（有限数字），例如：对于变量：

print([PV_flux, Ta_kelvin])
[535.5, 297]

当我在 nsolve 中使用它们时：

print(nsolve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], [J_PV_1, J_PV_2, J_m1, J_m2, Tm1, Tm2, Tpv],
             [PV_flux, PV_flux, PV_flux, PV_flux, Ta_kelvin, Ta_kelvin, Ta_kelvin], domain=S.Reals))

加注：

Traceback (most recent call last):
  File "C:/Users/", line 154, in <module>
    a=temp(535.5,240.1,86.1, 86.1,23,24)
  File "C:/Users/", line 105, in temp
    return(nsolve([eq1,eq2,eq3,eq4,eq5,eq6,eq7], [J_PV_1,J_PV_2,J_m1,J_m2,Tm1,Tm2,Tpv], [PV_flux,PV_flux,PV_flux,PV_flux,Ta_kelvin,Ta_kelvin,Ta_kelvin], domain = S.Reals))
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\sympy\utilities\decorator.py", line 92, in func_wrapper
    return func(*args, **kwargs)
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\sympy\solvers\solvers.py", line 3026, in nsolve
    x = findroot(f, x0, J=J, **kwargs)
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\mpmath\calculus\optimization.py", line 960, in findroot
    for x, error in iterations:
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\mpmath\calculus\optimization.py", line 658, in __iter__
    s = self.ctx.lu_solve(Jx, fxn)
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\mpmath\matrices\linalg.py", line 227, in lu_solve
    A, p = ctx.LU_decomp(A)
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\mpmath\matrices\linalg.py", line 143, in LU_decomp
    ctx.swap_row(A, j, p[j])
  File "C:\Users\bents\anaconda3\envs\ben\lib\site-packages\mpmath\matrices\matrices.py", line 868, in swap_row
    A[i,k], A[j,k] = A[j,k], A[i,k]
  File "C:\Users\\anaconda3\envs\ben\lib\site-packages\mpmath\matrices\matrices.py", line 491, in __getitem__
    if key[0] >= self.__rows or key[1] >= self.__cols:
TypeError: '>=' not supported between instances of 'NoneType' and 'int'

谢谢

【问题讨论】：

这些变量是什么类型的？它们都必须是有限数。你试过print([PV_flux,PV_flux,PV_flux,PV_flux,Ta_kelvin,Ta_kelvin,Ta_kelvin]) 来检查它们吗？
它们是有限数。 nsolve 命令在函数内部，这些变量是函数的输入。
@JohanC 已编辑。谢谢
你能展示实际的可运行代码吗？
@OscarBenjamin 我编辑了这篇文章。参数多，有点长

标签： sympy

【解决方案1】：

您选择的特定值对应于方程的奇异性，例如：

In [29]: eq5
Out[29]: 
                                                                                0.25               
                                                4                    ⎛Tpv - 297⎞                   
19.566⋅J_PV_1 + 19.566⋅J_PV_2 - 2.2187844e-6⋅Tpv  - 13.7890542929877⋅⎜─────────⎟    ⋅(Tpv - 297) - 
                                                                     ⎜Tpv   297⎟                   
                                                                     ⎜─── + ───⎟                   
                                                                     ⎝ 2     2 ⎠                   

                            0.333333333333333                      
                 ⎛Tpv - 297⎞                                       
18.8042682352567⋅⎜─────────⎟                 ⋅(Tpv - 297) + 924.177
                 ⎜Tpv   297⎟                                       
                 ⎜─── + ───⎟                                       
                 ⎝ 2     2 ⎠

如果您使用Tpv = 297，那么您将遇到数字问题。确实有雅可比行列式中的 nan：

In [5]: eqs = [eq1, eq2, eq3, eq4, eq5, eq6, eq7]

In [6]: syms = [J_PV_1, J_PV_2, J_m1, J_m2, Tm1, Tm2, Tpv]

In [7]: vals = [535.5, 535.5, 535.5, 535.5, 297, 297, 297]

In [8]: rep = dict(zip(syms, vals))

In [9]: Matrix(eqs).jacobian(syms).subs(rep)
Out[9]: 
⎡     -21.74              0          1.53725014229955           0                 0                 0          116.255751364922 ⎤
⎢                                                                                                                               ⎥
⎢       0               -21.74               0          1.53725014229955          0                 0          116.255751364922 ⎥
⎢                                                                                                                               ⎥
⎢1.53725014229955         0          -15.3725014229955          0          73.0713156818435         0                  0        ⎥
⎢                                                                                                                               ⎥
⎢       0          1.53725014229955          0          -15.3725014229955         0          73.0713156818435          0        ⎥
⎢                                                                                                                               ⎥
⎢     19.566            19.566               0                  0                 0                 0          -232.511502729845⎥
⎢                                                                                                                               ⎥
⎢       0                 0          12.2980011383964           0                nan                0                  0        ⎥
⎢                                                                                                                               ⎥
⎣       0                 0                  0          12.2980011383964          0                nan                 0        ⎦

如果您将 297 更改为 298，那么 nsolve 将成功：

In [20]: vals = [535.5, 535.5, 535.5, 535.5, 298, 298, 298]

In [21]: nsolve(eqs, syms, vals)
Out[21]: 
⎡593.17100732855 ⎤
⎢                ⎥
⎢593.142501153849⎥
⎢                ⎥
⎢491.462389303786⎥
⎢                ⎥
⎢491.059251115046⎥
⎢                ⎥
⎢304.687760815851⎥
⎢                ⎥
⎢304.609734541627⎥
⎢                ⎥
⎣321.670012318818⎦

【讨论】：