Solving a constrained linear system of equations using Scipy

I am trying to solve a linear equation (AX=b) where A is an 8x8 matrix, and X and b are 8x1 vectors. X can be expressed by X=[x1, y1, x2, y2, x3, y3, x4, y4]. However, I have the following 3 constraints that I need to respect:

1) 0.5*(y1 + y2) = 0
2) 0.5*(x3 + x4) = 0
3) 0.5*(y3 + y4) = 0

First, I tried to integrate the constraints into matrix A (which became of shape 11x8) and multiplied both parts of the equation by the pseudo-inverse of A to get the solution X, but the solution I found with this method still did not satisfy the constraints. So, I created an optimzation problem in python using scipy. However, I think I am doing something wrong (probably in my objective function) since the solver message states that the optimization was terminated successful and the found solution does indeed satisfy the constraints, but when I multiply the matrix A by the found solutions for the vector X, I do not get the values of vector b. So, the system is not being solved.

I tried the following code :

A = np.array([
    [-261.60,  11.26,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [   4.07, -12.75,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,  -158.63,  -5.65,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,    -2.81, -12.14,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0, -265.99,  19.29,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,   12.59, -12.34,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0, -166.25, -12.63],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0,   -8.40, -11.14]
])

b = np.array([
     -6.95,
     16.35,
     -0.96,
     16.35,
     19.19,
    -15.85,
    -12.36,
    -15.63]).reshape(-1, 1)

def objective_function(x):
    return np.sum((np.dot(A, x) - b)**2)

def constraints(x):
    return np.array([
        0.5*(x[1] + x[3]),  # 0.5*(y1 + y2) = 0
        0.5*(x[4] + x[6]),  # 0.5*(x3 + x4) = 0
        0.5*(x[5] + x[7])   # 0.5*(y3 + y4) = 0
    ])


cons = {'type': 'eq', 'fun': constraints}
options = {'disp': True, 'verbose': 4}

res = optimize.minimize(objective_function, np.zeros(8), method='SLSQP', constraints=cons, options={'disp': True})

x_optimized = res.x

# Checking if solution was found
print(np.matmul(A, x_optimized))

At the end, the output of np.matmul(A, x_optimized) is :

[ 0.01723619  0.00041891  0.01781721 -0.00034223  0.0113301  -0.00848171
  0.01008367  0.00781168]

However, I was expecting the same values as vector b if the system was actually solved :

[ -6.95  16.35  -0.96  16.35  19.19 -15.85 -12.36 -15.63]

Since your constraints are linear you can put them all in a linear equation.

You can use least squares solution, if the system has exact solution you will get it, otherwise you will get as close as possible.

AC = np.zeros([3, A.shape[1]])
bC = np.zeros((3, 1))
AC[0][[1, 3]] = 0.5 # 0.5 * (y1 + y2) = 0
AC[1][[4, 6]] = 0.5 # 0.5 * (x3 + x4) = 0
AC[2][[5, 7]] = 0.5 # 0.5 * (y3 + y4) = 0

# Solve the augmented system
A = np.vstack([A, AC])
b = np.vstack([b, bC])
x, _, _, _ = np.linalg.lstsq(A, b, rcond=None)