Find Fraction using LP

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. ]

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.