from scipy.optimize import linprog

points = data['Points']
cost = data['Cost']
gk   = data['Position'].apply(lambda x: 1 if x == 'GK' else 0)
defe = data['Position'].apply(lambda x: 1 if x == 'DEF' else 0)
mid  = data['Position'].apply(lambda x: 1 if x == 'MID' else 0)
stri = data['Position'].apply(lambda x: 1 if x == 'STR' else 0)
xi   = np.ones(data.shape[0])

A_upperbounds = np.array([cost, defe, mid, stri])
b_upperbounds = np.array([100, 5, 5, 3])
A_equality = np.array([gk, xi])
b_equality = np.array([1, 11])
bounds = [(0, 1) for x in range(data.shape[0])]

solution = linprog(
    c=-points,
    A_ub=A_upperbounds,
    b_ub=b_upperbounds,
    A_eq=A_equality,
    b_eq=b_equality,
    bounds=bounds
)