Find Fraction using LP

Question

I have two polygons BP and GP described by the set of inequalities constraints -x+y<=1 and x+y<= 5 and x-y<=3 and -y <= 0 for the Black Polygon and -1<=x<=4 and 0 <= y <= 3 for the Green Polygon.

My goal here is to use an LP to find the optimal solution of a fraction problem: given that B is in GP what is the maximum value lambda such that

B = lambda*B_BP + (1-lambda)*B_GP

In other words I would like to find the largest fraction of Bthat is inside the polygon in the above sense. For that I am struggling to write a LP program, I think that if we write BP as a matrix inequality condition we get that every B_BP is such that M_BP*B_BP <= C were C is a the vector (1,5,3,0) and M_BP is the matrix ((-1,1),(1,1),(1,-1),(0,-1)). So I think it should go with something like, given B = x_1+x_2

maximize lambda

subject to M_BP*L_BP <= C_B

and B_BP >= 0

Where I suppose (this is all my attempt, might be all very wrong) that L_BP = (x,y) vector and lambda = (x+y)/normalization and also that C_B relates somehow to the vector B.

Sorry if my first question is too messy, I'm starting here.

Answer 1

If I understand your problem correctly, I think it is possible to solve this exact, in a mathematical sense. Let me explain. Since the objective function is linear in lambda (as Superkogito points out), a maximum (or minimum) is always achieved in one of the corner points. By using this, you can calculate lambda.

Let me start with some simple examples. For any point within the black polygon it is clear that lambda is 1: you can just put B = B_BP. Now let's take B = (-1, 3). If you would take any black point other than B_BP = (-1, 0) and have lambda > 0, then with any green point within the square your x-coordinate will be higher than -1. So the best you can do is put B_BP = (-1, 0). Then the green point should be B_GP = (-1, 3), therefore lambda = 0.

Following the same logic, you can see that on the edge defined by the endpoints (-1, 0) and (-1, 3) you would always use B_BP = (-1, 0) and B_GP = (-1, 3). Going up this edge, lambda will decrease from 1 to 0. In (-1, 1) lambda = 2/3, in (-1, 2) lambda = 1/3. Same thing for the upper edge between (-1, 3) and (2, 3): in (0, 3) lambda = 1/3 and so forth. For the green triangle with corner (4, 3) you have to use (4, 3) as endpoint. Then in (3, 3) for example lambda = 1/2.

The interesting issue is obviously in the inner parts of the triangles. Here also B_GP is one of the corners, and B_BP is on the black line that is the side of the triangle. You can show this by assuming that there is another solution where this is not the case, and then proving that it is possible to increase lambda, by shifting B_GP or B_BP a little bit to the left or to the right. A complete mathematical proof would be too long for here, I guess. Now let's take the triangle on the left, with green corner (-1, 3) and the black side between (-1, 0) and (2, 3). Then for maximal lambda, B_GP = (-1, 3) and B_BP is on the black side. Because B = lambda * B_BP + (1 - lambda) * B_GP, this gives you two equations. And because B_BP is on the line y = x + 1, that gives you a third equation. From these you can solve the x- and y-coordinates of B_BP and lambda (given a point B).

I have done this and arrive at lambda = (Bx - By + 4) / 3. The coordinates of B_BP are then B_BPx = (Bx + 1 + lambda) / lambda and B_BPy = (By - 3 + 3 * lambda) / lambda (just fill in the lambda). For example, for the point (0, 2) you would have lambda = 2/3. You can do exactly the same for the other two triangles, I have done that as well.

I'll summarize:

For the left triangle: lambda = (Bx - By + 4) / 3

For the upper right triangle: lambda = (-By - Bx + 7) / 2

For the lower right triangle: lambda = By - Bx + 4

Programming this is now trivial. If you want, I can give you the corresponding coordinates of B_BP for the other two triangles, please let me know. By the way, you can arrive at them by drawing a straight line through the corner of the triangle and B.

Of course this only works when the objective function is linear. In other cases you would have to use an approach as suggested by Superkogito. I hope I have understood your question correctly, please forgive me if not! At least I found it a nice mathematical challenge.

Answer 2

The problem needs a better formulation in my opinion. I am not sure if this solves your issue but hopefully it helps. So I suggest using scipy.optimize.minimize to solve this optimization problem and just by inverting the sign or using the inverse you can transform your maximization into a minimization.

Also, since you are basing your code on random points from BP, GP and the random point B, you should feed those also into your input vector. From the input vector you can compute the lambda coefficients (I named this variable k in my code). This approach will return the values of the input vector that fulfills the constraints with the smallest output of the objective function fun aka the largest kx and largest ky.

The previously explained approach can be implemented as follows:

import numpy as np
import scipy.optimize


# compute k
def k(x):
    x_bp, y_bp = x[0], x[1]
    x_gp, y_gp = x[2], x[3]
    bx  , by   = x[4], x[5]

    # compute k
    kx = (bx - x_gp) / (x_bp - x_gp)
    ky = (by - y_gp) / (y_bp - y_gp)
    return kx, ky


# define objctive/cost function
def fun(x):
    kx, ky = k(x)
    return 1 / kx + 1 / ky

# define constraints (bp raandom point constraints)
cons = ({'type': 'ineq', 'fun': lambda x:  x[0] + x[1] - 5},
        {'type': 'ineq', 'fun': lambda x: -x[0] - x[1] - 3})

# init bounds 
bnds = ((None, None), (0, None), 
        (-1., 4.), (0., 3.), 
        (-1., 4.), (0., 3.))

# init vars vec
#x = [x_bp, y_bp, x_gp, y_gp, bx,   by]
x0 = [ 0.1,  0.2,  0.3,  0.5, 0.6, 0.5]

# optimize 
res = scipy.optimize.minimize(fun, 
                              x0, 
                              bounds=bnds, 
                              constraints=cons)

# print result
print("optimal x: ", np.round(res.x, 3))

# print optimal k 
print("optimal k: ", np.round(k(res.x), 3))

Note that you might need to tweak the code a bit and play with the initial value x0 and bounds but this should do the trick. The posted snippet results in the following output:

optimal x:  [0.1 0.2 0.3 0.5 0.6 0.5]
optimal k:  [-1.5 -0. ]