scipy 带约束的非线性方程组

enxuqcxy  于 2023-08-05  发布在  其他
关注(0)|答案(1)|浏览(154)

大家好,我想解一个热力学方程组。不幸的是,我是Python的初学者,需要一些帮助。我试着用最小二乘法来解决它,但结果总是得到我的初始向量。
希望有人能帮忙!

from scipy.optimize import least_squares
import numpy as np

def equations(x):
    eq1 = x[7] - 1 / (np.exp((5.24677 - 1598.673 / (60.0 - -46.424)) * np.log(10.0)))
    eq2 = x[4] - 1 / (np.exp((5.08354 - 1663.125 / (60.0 - -45.622)) * np.log(10.0)))
    eq3 = 0.0 - (0.5 * 0.4) + x[2] * x[9] - x[14] * (0.02 * 0.1)
    eq4 = 0.0 - (0.5 * 0.6) + x[2] * x[5] - x[12] * (0.02 * 0.1)
    eq5 = 0.0 - (0.4 * 0.6) + x[3] * x[0] + x[14] * (0.02 * 0.1)
    eq6 = 0.0 - (0.4 * 0.4) + x[3] * x[15] + x[12] * (0.02 * 0.1)
    eq7 = 0.0 - (x[10] - x[7] * x[13])
    eq8 = 0.0 - (x[8] - x[4] * x[11])
    eq9 = x[14] - (1.0 * x[6] * (x[13] - x[9]) + x[9] * (x[14] + x[12]))
    eq10 = x[12] - (1.0 * x[6] * (x[11] - x[5]) + x[5] * (x[14] + x[12]))
    eq11 = x[14] - (0.01 * x[1] * (x[0] - x[10]) + x[0] * (x[14] + x[12]))
    eq12 = x[12] - (0.01 * x[1] * (x[15] - x[8]) + x[15] * (x[14] + x[12]))
    eq13 = 1.0 - (x[9] + x[5])
    eq14 = 1.0 - (x[0] + x[15])
    eq15 = 1.0 - (x[13] + x[11])
    eq16 = 1.0 - (x[10] + x[8])
    return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16]

# constraints
bounds = (np.zeros(16), np.array([1, np.inf, np.inf, np.inf, 1, np.inf, np.inf, 1, 1, 1, 1, 1, np.inf, 1, np.inf, 1]))

# inital 
x0 = np.zeros(16)
# x0 = np.array([1, 1e-4, 1, 1, 1, 1, 1e-4, 1, 1, 1, 1, 1, 1, 1, 1, 1])

res = least_squares(equations, x0, bounds=bounds)

print(res.x)

字符串

j0pj023g

j0pj023g1#

一个可能的问题是,第一和第二RHS的规模是巨大的相比,其他。x4和x7看起来是固定的,但您没有在初始值中反映这一点。事实上,x4和x7的上限是1,这在第一个和第二个方程中是完全没有意义的。所以这个方程组需要修正,才是合理的。
我在下面展示了忽略前两个方程(重新排序为方程4和7)会产生一个看似合理的解决方案。

from scipy.optimize import least_squares
import numpy as np

def sides(x: np.ndarray) -> np.ndarray:
    # x order as original, y reordered so that x indices are closer to their use
    return np.array((
        (x[14], 0.01*x[1]*(x[ 0] - x[10]) + x[ 0]*(x[14] + x[12])),
        (x[12], 0.01*x[1]*(x[15] - x[ 8]) + x[15]*(x[14] + x[12])),
        (0,     -0.5*0.6 + x[2]* x[ 5] - x[12]*(0.02 * 0.1)),
        (0,     -0.5*0.4 + x[2]* x[ 9] - x[14]*(0.02 * 0.1)),
        (x[4],  1 / np.exp((5.08354 - 1663.125/(60.0 - -45.622)) * np.log(10))),
        (x[12], x[6]*(x[11] - x[5]) + x[5]*(x[14] + x[12])),
        (x[14], x[6]*(x[13] - x[9]) + x[9]*(x[14] + x[12])),
        (x[7],  1 / np.exp((5.24677 - 1598.673/(60.0 - -46.424)) * np.log(10))),
        (1,     x[10] + x[ 8]),
        (1,     x[ 9] + x[ 5]),
        (0,     x[10] - x[7]*x[13]),
        (0,     x[ 8] - x[4]*x[11]),
        (0,     -0.4*0.4 + x[3]* x[15] + x[12]*(0.02 * 0.1)),
        (1,     x[13] + x[11]),
        (0,     -0.4*0.6 + x[3]* x[ 0] + x[14]*(0.02 * 0.1)),
        (1,     x[ 0] + x[15]),
    ))

def equations(x: np.ndarray) -> np.ndarray:
    left, right = sides(x)[[
        0, 1, 2, 3,
        5, 6,
        8, 9, 10, 11, 12, 13, 14, 15,
    ], :].T
    return left - right

# Ordered as original
x0 = np.array((
    1, 1e-4, 1, 1,
    10,  # 1 / np.exp((5.08354 - 1663.125/(60.0 - -45.622)) * np.log(10)),
    1, 1e-4,
    10,  # 1 / np.exp((5.24677 - 1598.673/(60.0 - -46.424)) * np.log(10)),
    1, 1, 1, 1, 1, 1, 1, 1,
))

# Ordered as original, with x[4] and x[7] bounds changed to inf (not fixed)
bounds = np.array((
    np.zeros(16),
    (
        1, np.inf, np.inf, np.inf,
        np.inf,  # 1,
        np.inf, np.inf,
        np.inf,  # 1,
        1, 1, 1, 1, np.inf, 1, np.inf, 1,
    ),
))

res = least_squares(fun=equations, x0=x0, bounds=bounds)
assert res.success, res.message

for i, (x, (left, right)) in enumerate(zip(res.x, sides(res.x))):
    print(f'y{i:<2} = {left:8.3g} ~ {right:9.3g}, '
          f'err={right-left:8.1e}, '
          f'x{i:<2} = {x:.3g}')

个字符
这个助手程序将向您展示一些选项,其中两个有问题的方程中的值可以被替换以在原始问题中产生合理的解决方案:

import numpy as np
from sympy import Symbol, Eq, solveset

for y, values in (
    (1.690, (5.08354, 1663.125, 60, -45.622)),
    (0.138, (5.24677, 1598.673, 60, -46.424)),
):
    a = Symbol('a', real=True, infinite=False)
    b = Symbol('b', real=True, infinite=False)
    c = Symbol('c', real=True, infinite=False)
    d = Symbol('d', real=True, infinite=False)
    variables = a, b, c, d

    eq = Eq(-np.log10(y), a - b/(c - d))

    print('Do one of:')
    for target, old_val in zip(variables, values):
        seq = eq.subs({
            s: v for s, v in zip(variables, values)
            if s is not target
        })
        replacement, = solveset(seq, target)

        print(f'Replace {target}={old_val} with {replacement:.1f}')
    print()
Do one of:
Replace a=5.08354 with 15.5
Replace b=1663.125 with 561.0
Replace c=60 with 267.5
Replace d=-45.622 with -253.1

Do one of:
Replace a=5.24677 with 15.9
Replace b=1598.673 with 466.8
Replace c=60 with 318.0
Replace d=-46.424 with -304.4

的字符串

相关问题