solving Ax =b with constraints for a non-square matrix A, using python

Question

I solved this problem in Mathcad, but I don't know how to transfer it to Python.

Mathcad:

I tried this code:

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)

Result: {(0, 1.66644368166375x4 - 2.70780339743064x5, 3.99892914904867x4 - 6.49785771642014x5, x4, x5)}

it's fine if x4 and x5 are not found, but the values of each x should be > 0 and < 1. I do not know how to do this with e.g. numpy/scipy.

Answer 1

Here is a solution with scipy.optimize.minimize() that follows this answer to a related question and that tries to replicate your Mathcad solution:

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]

The result is pretty close to what you found in Mathcad, although the way to get there was admittedly a bit more tedious here. Note: just like your Mathcad solution but unlike your sympy solution, the approach with scipy.optimize.minimize() picks a single solution from the solution space.

Answer 2

There is not one solution. There is an infinite number of solutions. Sympy has (poor, but present) support for inequality solution, so use it:

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}]')

F1 = 1.11022302462516e-16*F4
F2 = 1.66644368166375*F4 - 2.70780339743065*F5
F3 = 3.99892914904867*F4 - 6.49785771642014*F5
0.001 < F4
0.001 < F5
F4 < 1.0
F5 < 1.0
F4 > 1.62489943538159*F5 + 0.000250066946106784
F4 > 1.62489943538159*F5 + 0.000600080285342505
F4 < 1.62489943538159*F5 + 0.250066946106784
F4 < 1.62489943538159*F5 + 0.600080285342505

Example solutions:
F=[0.0000, 0.0010, 0.0024, 0.0022, 0.0010] eqs=[+1.8e-15, +1.4e-14, -7.1e-15]
F=[0.0000, 0.0010, 0.0024, 0.0022, 0.0010] eqs=[+1.8e-15, +1.4e-14, -7.1e-15]
F=[0.0000, 0.1049, 0.2518, 0.0646, 0.0010] eqs=[ 0.     , +1.6e-13, +2.9e-13]
F=[0.0000, 0.0010, 0.0024, 0.0646, 0.0394] eqs=[+2.8e-13, +2.3e-13, +2.3e-13]
F=[0.0000, 0.2089, 0.5012, 0.1270, 0.0010] eqs=[+4.5e-13, +3.9e-13, +5.2e-13]
F=[0.0000, 0.0010, 0.0024, 0.1270, 0.0778] eqs=[+2.5e-13, +9.1e-13,  0.     ]
F=[0.0000, 0.3128, 0.7506, 0.1893, 0.0010] eqs=[ 0.     , +3.9e-13, +5.2e-13]
F=[0.0000, 0.0010, 0.0024, 0.1893, 0.1161] eqs=[+9.4e-14, +9.1e-13,  0.     ]
F=[0.0000, 0.4167, 1.0000, 0.2517, 0.0010] eqs=[+1.8e-12, +1.3e-12, +5.2e-13]
F=[0.0000, 0.0010, 0.0024, 0.2517, 0.1545] eqs=[+1.6e-12, +9.1e-13,  0.     ]

This means that, so far as all of the equalities and inequalities above are satisfied, any value for F4 and F5 is acceptable.

Also note that the system is not satisfiable if F0 is constrained to > 1e-3.