Linear Programming in 'scipy.linprog' with non-standard functions

Question

I am trying to learn about implementation of linear programming (LP) problems in scipy.linprog. I understand how it works with basic functions, for example:

max 2x+3y
st. 2x-y <= 0
    5x+y >= -10

from scipy.optimize import linprog
c = [-2, -3]
A = [[2, -1], [-5, -1]]
b = [0, 10]

Now, when it comes to more complex cases, like this there are quite a few things I don't understand.

First, I don't understand how and why Minimize 1*abs(x_0) + 1*abs(x_1) + ... + 1*abs(x_n)) is changed to Minimize 1*y_0 + 1*y_1 + ... + 1*y_n by adding the following constraint -y_i <= x_i <= y_i // for i in n. EDIT: I understand now that y_i is a slack variable, which also is the reason that the A_aux matrix is added. But I don't understand that shape of the Matrix.

Second, here they use the constraint to c = np.hstack((np.zeros(N), np.ones(N_AUX))). In my example above, there are two variables, which gives a list of the length 2. But here there are N variables, but c is a list of N*2? In my example it is easy to understand the values, but here they are [0,0,0,0,0,1,1,1,1,1] and I don't see any hint in the function to why these values are chosen.

Third, I understand the A_orig matrix.

A_orig = [[0, 1, -1, 0, 0, 0, 0, 0, 0, 0],   # orig constraint 1
          [0, 0, 1, -1, 0, 0, 0, 0, 0, 0],   # orig constraint 2
          [-1, -1, 0, 0, 0, 0, 0, 0, 0, 0],  # more interesting problem
          [0, -1, -1, 0, 0, 0, 0, 0, 0, 0]]  # ""   ""          ""

Each row simply represents the number of x_i and y_i for i in range(1,6) for every constraint.

But then a second matrix A_aux is added:

A_aux = [[-1, 0, 0, 0, 0, -1, 0, 0, 0, 0],
         [0, -1, 0, 0, 0, 0, -1, 0, 0, 0],
         [0, 0, -1, 0, 0, 0, 0, -1, 0, 0],
         [0, 0, 0, -1, 0, 0, 0, 0, -1, 0],
         [0, 0, 0, 0, -1, 0, 0, 0, 0, -1],
         [1, 0, 0, 0, 0, -1, 0, 0, 0, 0],
         [0, 1, 0, 0, 0, 0, -1, 0, 0, 0],
         [0, 0, 1, 0, 0, 0, 0, -1, 0, 0],
         [0, 0, 0, 1, 0, 0, 0, 0, -1, 0],
         [0, 0, 0, 0, 1, 0, 0, 0, 0, -1]]

Which is a skew-symmetric matrix. And I have no idea and can't find any information about why this is added. I don't understand why it is 10*10 either.

Finally, bounds are added. bnds = [(0, 50) for i in range(N)] + [(None, None) for i in range(N_AUX)] # some custom bounds. Here I find it strange that they set (0,50) for the first 5 rows, and then no bounds for the last 5 rows. I don't see any indication to why it should be bounded by 50?

There are quite a lot of questions here. But if someone could point me in the right direction I would be very grateful.

Thank you.

Answer 1

This is rather a math question than a programming question. Anyway, to answer your other questions:

In my example above, there are two variables, which gives a list of the length 2. But here there are N variables, but c is a list of N*2?

You need to introduce an additional variable y_i for each abs(x_i). You have N terms abs(x_i), therefore you'll have N+N = 2*N variables after linearizing/reformulating all the absolute values. Each variable x doesn't occurs in the objective in the reformulation, therefore each coefficient of x_i is zero in the objective. And obviously, each one of y_i is 1.

Which is a skew-symmetric matrix. And I have no idea and can't find any information about why this is added.

You need the matrix A_aux to model the constraints -y_i <= x_i <= y_i. Both y_i and x_i are optimization variables, so you can't pass this constraints as variable bounds. This only works if the bounds are given constants like -1 <= x_i = 2.

Let's say you have the constraints -y_1 <= x_1 <= y_1 and -y_2 <= x_2 <= y_2. This is the same as:

-x_1 - y_1 <= 0
 x_1 - y_1 <= 0
-x_2 - y_2 <= 0
 x_2 - y_2 <= 0

Adding the hidden zeros, this is the same as

-1*x_1 + 0*x_2 - 1*y_1 + 0*y_2 <= 0
 1*x_1 + 0*x_2 - 1*y_1 + 0*y_2 <= 0
 0*x_1 - 1*x_2 + 0*y_1 - 1*y_2 <= 0
 0*x_0 + 1*x_2 + 0*y_1 - 1*y_2 <= 0

Can you see it now? It's just a simply matrix vector product:

(-1  0  -1   0 ) (x_1) <= 0
( 1  0  -1   0 ) (x_2) <= 0
( 0 -1   0  -1 ) (y_1) <= 0
( 0  1   0  -1 ) (y_2) <= 0

Here, the left matrix is just A_aux. If you change the row order, you'll get exactly the same matrix as in your linked answer.

Finally, bounds are added. bnds = [(0, 50) for i in range(N)] + [(None, None) for i in range(N_AUX)] # some custom bounds. Here I find it strange that they set (0,50) for the first 5 rows, and then no bounds for the last 5 rows. I don't see any indication to why it should be bounded by 50?

The first N variables (x) are non-negative after the reformulation. The N addition variables y are unbounded, i.e. can be negative or positive, there there are no bounds. However, there's no real need to bound the x variables with 50 from above. It's only meant to be a sharper bound of the optimization variables based on the right-hand side of the other linear constraints in the example.

PS: I can second Sascha's comment in your linked answer, that you should prefer a library like cvxpy if you don't want to do all these transformations/reformulations on your own.