How to formulate a linear minimization problem with scipy

I have a problem of the form:

min (x1 - k1) + (x2 -k2) + ...

k are constants

with linear constraints:

A*x = B

B's size is 3x1, A is 3xn and x is nx1

I also have linear bounds inequalities in the form

l1 <= x1 <= u1

l2 <= x2 <= u2 ...etc

I'm trying to resolve it in python with scipy :

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

A = \
np.array([[ 0.11,  0.1333,  0.1333,  0.01],
          [ 0.02,  6.667,  0.1333,  0.12],
          [0.0933,  0.6667,  0.6,  0.01]])

B = \
np.array([[25],
          [57],
          [28]])

l = \
np.array([[50],
          [20],
          [5],
          [50]])
u = \
np.array([[200],
          [80],
          [20],
          [200]])

def objfun(x):
    return np.linalg.norm((x - (u-l)/2))


x0 = \
np.array([[1],
          [2],
          [3],
          [4]])

bounds = Bounds(l, u)


eq_cons = {'type': 'eq', 'fun': lambda x: np.matmul(A,x)-B}
res = minimize(objfun, x0,  method='SLSQP',
               constraints=[eq_cons], options={'ftol': 1e-9, 'disp': True},
               bounds=bounds)

but it gives me the following error:

all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

Thanks for the help.

As already mentioned in the comments and Rorys answer, you should use scipy.optimize.linprog for linear problems (LPs), i.e. your objective function and the constraint functions are linear. Only use scipy.optimize.minimize for nonlinear problems. The former interfaces sophisticated high-performance solvers (like HiGHS) for LPs, and thus it will have much better performance in most cases.

That being said, it's also important to stress that minimizing the linear objective function

c^T * (x-k) = 1*(x1-k1) + 1*(x2-k2) + ...

is not equivalent to minimizing the euclidean vector norm of your objective

||c^T * (x-k)|| = sqrt( (1*(x1-k1) + 1*(x2-k2) + ...)**2 )

as it makes your objective function nonlinear and changes the fundamental properties of your optimization problem. Spoiler alert: your LP is infeasible, i.e. it doesn't have an optimal solution whilst your nonlinear problem (NLP) has one.

Here's a code snippet that solves your LP via linprog:

import numpy as np
from scipy.optimize import linprog

# Linear constraint A*x = B
A = np.array([[ 0.11,  0.1333,  0.1333,  0.01],
              [ 0.02,  6.667,  0.1333,  0.12],
              [0.0933,  0.6667,  0.6,  0.01]])

B = np.array([25, 57, 28])

# Bounds l <= x <= u
l = np.array([50, 20, 5, 50])
u =  np.array([200, 80, 20, 200])
bounds = [(ll, uu) for ll,uu in zip(l, u)]

# objective f(x) = c^T * x - c^T*k (vector notation)
c = np.ones(4)
k = (u-l)/2

# we minimize c^T*x  s.t.  A*x = B, l <= x <= u
# In order to obtain the objective value of the original LP, we
# simply subtract c^T*k from the objective value res.fun

res = linprog(c, A_ub=None, b_ub=None, A_eq=A, b_eq=B, bounds=bounds)