import numpy as np
from scipy.optimize import minimize

def objective(x):
    sum = 0
    for i in range(4):
        for j in range(3):
            sum = sum + x[i][j]*q[i][j]
    return -1*sum

def constraint1(i):
    def g(x):
        con = 20
        for k in range(3):
            con = con - x[i][k]
        return con
    return g

def constraint2(k):
    def h(x):
        sum_eq = 20
        for i in range(4):
            sum_eq = sum_eq - x[i][k]
        return sum_eq
    return h

# initial guesses
x0 = np.array([[5,5,5],[5,5,5],[5,5,5],[5,5,5]])
q = np.array([[5,4,3],[5,4,3],[5,4,3],[5,4,3]])

# show initial objective
print('Initial Objective: ' + str(objective(x0)))

# optimize
cons = []
for i in range(4):
    cons.append({'type': 'ineq', 'fun': constraint1(i)})

for k in range(3):
    cons.append({'type': 'ineq', 'fun': constraint2(k)})


solution = minimize(objective,x0,method='SLSQP',\
                    bounds=None,constraints=cons)
x = solution.x

# show final objective
print('Final Objective: ' + str(objective(x)))

# print solution
print('Solution')
print('x00 = ' + str(x[0][0]))
print('x01 = ' + str(x[0][1]))
print('x02 = ' + str(x[0][2]))
print('x10 = ' + str(x[1][0]))
print('x11 = ' + str(x[1][1]))
print('x12 = ' + str(x[1][2]))
print('x20 = ' + str(x[1][0]))
print('x21 = ' + str(x[2][1]))
print('x22 = ' + str(x[2][2]))
print('x23 = ' + str(x[2][3]))
print('x30 = ' + str(x[3][0]))
print('x31 = ' + str(x[3][1]))
print('x32 = ' + str(x[3][2]))