The numpy array stores gekko objects at each row and column location. To access the gekko value for variable x
, you need to get x.value[0]
. You can convert the q
matrix to a numpy matrix with
# convert to matrix form
qr = np.array([[q[i,j].value[0] for j in range(4)] for i in range(4)])
Here is a complete version of your code that checks that the row and column constraints are satisfied.
import numpy as np
import scipy.optimize as opt
from gekko import GEKKO
p= np.array([4, 5, 6.65, 12]) #p = prices
pmx = np.triu(p - p[:, np.newaxis]) #pmx = price matrix, upper triangular
m = GEKKO(remote=False)
q = m.Array(m.Var,(4,4),lb=0,ub=10)
# only upper triangular can change
for i in range(4):
for j in range(4):
if j<=i:
q[i,j].upper=0 # set upper bound = 0
def profit(q):
profit = np.sum(q.flatten() * pmx.flatten())
return profit
for i in range(4):
m.Equation(np.sum(q[i,:])<=10)
m.Equation(np.sum(q[:,i])<=8)
m.Maximize(profit(q))
m.solve()
print(q)
# convert to matrix form
qr = np.array([[q[i,j].value[0] for j in range(4)] for i in range(4)])
for i in range(4):
rs = qr[i,:].sum()
print('Row sum ' + str(i) + ' = ' + str(rs))
cs = qr[:,i].sum()
print('Col sum ' + str(i) + ' = ' + str(cs))
The results of this script are:
[[[0.0] [2.5432017412] [3.7228765674] [3.7339217013]]
[[0.0] [0.0] [4.2771234426] [4.2660783187]]
[[0.0] [0.0] [0.0] [0.0]]
[[0.0] [0.0] [0.0] [0.0]]]
Row sum 0 = 10.000000009899999
Col sum 0 = 0.0
Row sum 1 = 8.5432017613
Col sum 1 = 2.5432017412
Row sum 2 = 0.0
Col sum 2 = 8.00000001
Row sum 3 = 0.0
Col sum 3 = 8.00000002
There is a slight violation of the constraints because of the solver tolerance. You can adjust the solver tolerance with m.options.ATOL=1e-6
to something that is more like 1e-8
if you'd like to have a more precise answer.