Optimize linear mathematical model SCIPY

Question

I'm having trouble to understand how should I translate my mathematical model to code such that I can use Scipy to optimize it.

The model is the following:

[from Budak, Gercek & Kara, Imdat & Ic, Yusuf & Kasimbeyli, Refail. (2019). New mathematical models for team formation of sports clubs before the match. Central European Journal of Operations Research. 27. 10.1007/s10100-017-0491-x.] [1]: https://i.stack.imgur.com/6U9YS.png

As you can see the model is a triple sumation so its kinda tricky for me to translate that to code as the exmples I've seen are fairly simple.

I dont, have much trouble with the constraints tho. Any help would be appreciated.

#Edit

Ok so i manage to use pulp to kind of solve the problem, however. I can't get one of my constraints to work.

Here is my code

opt_mod = plp.LpProblem("MIP Model", plp.LpMaximize)
x_vars = plp.LpVariable.dicts("x", pairs, cat=plp.LpBinary)

for i in set_i:
    opt_mod+=plp.lpSum(x_vars[i,j] for j in set_j)<=1

for j in set_j:
    opt_mod+=plp.lpSum(x_vars[i,j] for i in set_i)==1

for j in set_j:
    for y in set_y:
        opt_mod += plp.lpSum(float(SSiy[i,y]) * x_vars[i,j] 

for i in set_i) >= float(THjy[j,y])

for i in set_i:
    for j in set_j:
        opt_mod += x_vars[i,j] <= int(PMij[i,j])

opt_mod += plp.lpSum(SSiy[i,y]*SMyj[y,j]*MSj[j]* x_vars[i,j]
    for i in set_i for y in set_y for j in set_j)

However the first constraint is not working properly. It should force the solution to not have repeating values

Nonetheless:

opt_mod.solve() for x in x_vars: if x_vars[x].value()==1: print(x)

I get (0,1) (0,4) (0,5) (0,8) (0,9) (1,0)

etc.

Im not shure, why this is not working.