numpy 使用Python求解非方阵A的约束条件Ax =B

kdfy810k  于 2023-11-18  发布在  Python
关注(0)|答案(2)|浏览(126)

我在Mathcad中解决了这个问题,但我不知道如何将其转移到Python中。
Mathcad:
x1c 0d1x的数据
我试了这个代码:

from sympy import symbols, nonlinsolve
Q = np.array([230.8084,119.1916,76.943,153.8654,196.1346])
x1, x2, x3, x4, x5 = symbols('x1, x2, x3, x4, x5', real=True)
eq1 =  (Q[0]**2)*x1 - (Q[1]**2)*x2 + (Q[2]**2)*x3
eq2 =  (Q[0]**2)*x1 + (Q[3]**2)*x4 - (Q[4]**2)*x5 - (Q[2]**2)*x3
eq3 =  (Q[2]**2)*x3 - (Q[3]**2)*x4 + (Q[4]**2)*x5
q = nonlinsolve([eq1,eq2,eq3 ], [x1, x2, x3, x4, x5])
print(q)

字符串
结果:{(0,1.66644368166375x4 - 2.70780339743064x5,3.99892914904867x4 - 6.49785771642014x5,x4,x5)}
如果没有找到x4和x5也没关系,但是每个x的值应该> 0和< 1。我不知道如何使用例如numpy/scipy来做到这一点。

tkqqtvp1

tkqqtvp11#

下面是一个scipy.optimize.minimize()的解决方案,它遵循this answer的相关问题,并试图复制您的Mathcad解决方案:

import numpy as np
from scipy.optimize import minimize, Bounds

# We only ever need the squared values of Q, so we square them directly
q_sq = np.array([230.8084, 119.1916, 76.943, 153.8654, 196.1346]) ** 2

# Define the function to be optimized
def f(x):

    # Set up the system of equations ...
    rows = np.array([[q_sq[0], -q_sq[1], q_sq[2],        0,        0],
                     [q_sq[0], -q_sq[1],       0,  q_sq[3], -q_sq[4]],
                     [      0,        0, q_sq[2], -q_sq[3],  q_sq[4]]])
    y = rows @ x
    
    # ... the results of which should all be zero
    return y @ y

lower_bds = 0.001 * np.ones(5)  # Constrain all x >= 0.001
upper_bds = 0.999 * np.ones(5)  # Constrain all x <= 0.999
initial_guess = lower_bds + 1e-6  # Initial guess slightly above lower bounds

res = minimize(f, x0=initial_guess, bounds=Bounds(lb=lower_bds, ub=upper_bds))

print(f"Converged? {res.success}")
print(f"Result: {res.x}")
# >>> Converged? True
# >>> Result: [0.001      0.00416655 0.001      0.00193571 0.00103738]

字符串
结果非常接近你在Mathcad中找到的结果,尽管在这里实现的方法确实有点乏味。注意:就像你的Mathcad解决方案一样,但与你的symny解决方案不同,scipy.optimize.minimize()的方法从解决方案空间中选择一个解决方案。

pn9klfpd

pn9klfpd2#

没有一个解。有无限多的解。Sympy对不等式解有(差的,但现在)支持,所以使用它:

from sympy import linsolve, reduce_inequalities, symbols, Matrix
from sympy.core.numbers import Zero
from sympy.matrices.dense import matrix_multiply_elementwise as mul

Q = Matrix((230.8084, 119.1916, 76.943, 153.8654, 196.1346))
F = Matrix(symbols('F1, F2, F3, F4, F5', real=True, finite=True))
Q2F = mul(mul(Q, Q), F)
eqs = Matrix((
    (1, -1,  1,  0,  0),
    (1,  0, -1,  1, -1),
    (0,  0,  1, -1,  1),
)) * Q2F

sol, = linsolve(eqs, tuple(F))

ineqs = [
    sol[0] > 0
] + [
    s > 1e-3 for s in sol[1:]
] + [s < 1 for s in sol]
free_sym = next(
    sym
    for sol_entry in sol
    for sym in sol_entry.free_symbols
)
ineq_soln = reduce_inequalities(ineqs, [free_sym])

for sym, equal_to in zip(F, sol):
    if sym is not equal_to:
        print(f'{sym} = {equal_to}')
for predicate in ineq_soln.args:
    print(predicate)
print()

print('Example solutions:')
f4_min = 1.62489943538159*1e-3 + 0.000600080285342505
f4_max = 1.62489943538159*1e-3 + 0.250066946106784
for i in range(0, 101, 25):
    f4 = i/100*(f4_max - f4_min) + f4_min
    f5_min = 1e-3
    f5_max = (f4 - 0.000600080285342505)/1.62489943538159
    for f5 in (f5_min, f5_max):
        substitutions = {F[3]: f4, F[4]: f5}
        eq_errors = eqs.subs({
            sym: sol_entry.subs(substitutions)
            for sym, sol_entry in zip(F, sol)
        })
        f_desc = ', '.join(
            f'{s.subs(substitutions):.4f}'
            for s in sol
        )
        eq_desc = ', '.join(
            ' 0.     ' if isinstance(f, Zero)
            else f'{f:+.1e}'
            for f in eq_errors
        )
        print(f'F=[{f_desc}] eqs=[{eq_desc}]')

个字符
这意味着,只要满足上述所有等式和不等式,F4和F5的任何值都是可以接受的。
还应注意,如果F0被约束为> 1 e-3,则系统不可满足。

相关问题