How to fix Python Gekko Max Equation Length error

from gekko import GEKKO import numpy as np m = GEKKO() m.options.SOLVER=1 # optional solver settings with APOPT m.solver_options = ['minlp_maximum_iterations max', \ # minlp iterations with integer solution 'minlp_max_iter_with_int_sol max', \ # treat minlp as nlp 'minlp_as_nlp min', \ # nlp sub-problem max iterations 'nlp_maximum_iterations 5000', \ # 1 = depth first, 2 = breadth first 'minlp_branch_method 1', \ # maximum deviation from whole number 'minlp_integer_tol 0.0005', \ # covergence tolerance 'minlp_gap_tol 0.0001' ,\ 'max_cpu_time 10min'] n = 128 p = 8 q = 16 x = m.Array(m.Var,(n,n),integer=True) for i in range(n): for j in range(n): x[i,j].value = 1 x[i,j].lower = 0 x[i,j].upper = 1 s = 60 a = m.Array(m.Param,(x.shape[0],)) a_ = [172,282,10, 264, 287, 442, 393, 428, 484, 444, 344, 250, 293, 34, 473, 110, 338, 62, 250, 205, 81, 336, 249, 199, 328, 447, 408, 82, 357, 334, 181, 133, 378, 79, 292, 103, 494, 382, 10, 477, 237, 267, 337, 313, 395, 110, 114, 381, 52, 232, 457, 69, 167, 344, 363, 284, 136 ,240, 258, 449, 119, 317, 370, 404, 197, 102, 428, 238, 321, 103, 282, 37, 41, 86, 496 ,31, 148, 245, 78, 219, 37, 115, 378, 129, 37, 448, 415, 368, 238, 406, 408, 100, 112 ,289 ,205, 329, 338, 396, 494, 145, 355, 45, 5, 220, 411, 167, 85, 369, 426, 442, 406 ,217, 57, 176, 12, 368, 444, 227, 337, 63, 267, 216, 227, 182, 408, 116, 385, 140] for i in range(len(a)): a[i] = a_[i] O1 = np.array(range(n)) O_= np.transpose(np.nonzero(O1+1)) O = np.zeros((x.shape[0],x.shape[1]),dtype=int) for i in range(n): for j in range(n): O[i,O_[i]] = 1 d = m.Array(m.Param,(p,q)) for i in range(p): for j in range(q): d[i,j] = (round(p/2)-int(i/2))*4-(j%4) #condition m.Equation(m.sum(abs(x[:,:]-O[:,:]))<=2*s) for i in range(n): m.Equation(m.sum(x[:,i])==1) for i in range(n): m.Equation(m.sum(x[i,:])==1) #Objective m.Obj(np.sum(a.T*np.sum(x*d.reshape(1,n),axis=1))) #Set global options m.options.IMODE = 3 #Solve simulation m.solve(disp=False,debug=True)

Below is a version that should work, if you let it solve long enough although a solution is unlikely because you have 128*128 = 16384 binary variables. There are 2^16384 solutions and a Mixed Integer Nonlinear Programming solver such as APOPT is unlikely to find a solution quickly if the number of integer variables is greater than a few hundred. It can solve continuous variables with many thousands of variables but mixed integer is more difficult because it must solve them with a branch and bound method. You may want to try an MILP solver like CPLEX or Gurobi if you can make your problem linear (substitute slack variables for abs function).

Here are a few suggestions on your model:

Do not use numpy for model expressions. Functions such as np.abs and np.sum do not provide the needed information for the solvers. You should use m.abs2 (or m.abs3) and m.sum instead.
The arrays a and d can be a numpy arrays of constants instead of Gekko parameters.
You can break up at least one of your expressions into Intermediates to help with solution times.
The objective function is a scalar result of a.T x[:,:] d. You can use multiple m.Obj functions from a[i] x[i,j] d[j] instead of one large expression that results from the matrix multiplication.
I recommend that you start with a smaller problem and then scale-up to a larger problem as you gain confidence in the model and solution.

from gekko import GEKKO
import numpy as np

m = GEKKO()
m.options.SOLVER=1
# optional solver settings with APOPT
m.solver_options = ['minlp_maximum_iterations 1000', \
                    # minlp iterations with integer solution
                    'minlp_max_iter_with_int_sol 100', \
                    # treat minlp as nlp
                    'minlp_as_nlp 1', \
                    # nlp sub-problem max iterations
                    'nlp_maximum_iterations 5000', \
                    # 1 = depth first, 2 = breadth first
                    'minlp_branch_method 1', \
                    # maximum deviation from whole number
                    'minlp_integer_tol 0.0005', \
                    # covergence tolerance
                    'minlp_gap_tol 0.0001']

n = 128 
p = 8 
q = 16 

x = m.Array(m.Var,(n,n),integer=True)

for i in range(n):
    for j in range(n):
        x[i,j].value = 1
        x[i,j].lower = 0
        x[i,j].upper = 1


s = 60

# don't need a gekko array here, constants are okay
#a = m.Array(m.Param,(x.shape[0],))
a = np.array([172,282,10, 264, 287, 442, 393, 428, 484, 444, 344, 250, 293,  34, 473, 110, 338,  62,
 250, 205,  81, 336, 249, 199, 328, 447, 408,  82, 357, 334, 181, 133, 378,  79, 292, 103,
 494, 382,  10, 477, 237, 267, 337, 313, 395, 110, 114, 381,  52, 232, 457,  69, 167, 344,
 363, 284, 136 ,240, 258, 449, 119, 317, 370, 404, 197, 102, 428, 238, 321, 103, 282,  37,
  41,  86, 496  ,31, 148, 245,  78, 219,  37, 115, 378, 129,  37, 448, 415, 368, 238, 406,
 408, 100, 112 ,289 ,205, 329, 338, 396, 494, 145, 355,  45,   5, 220, 411, 167,  85, 369,
 426, 442, 406 ,217,  57, 176,  12, 368, 444, 227, 337,  63, 267, 216, 227, 182, 408, 116,
 385, 140])

# use .value to assign a value instead of over-writing the gekko variable
#   with a floating point number
#for i in range(len(a)):
#    a[i].value = a_[i]

O1 = np.array(range(n))
O_= np.transpose(np.nonzero(O1+1))

O = np.zeros((x.shape[0],x.shape[1]),dtype=int)
for i in range(n):
    for j in range(n):
        O[i,O_[i]] = 1

d = np.empty((p,q)) # doesn't need to be a gekko array
for i in range(p):
    for j in range(q):
        d[i,j] = (round(p/2)-int(i/2))*4-(j%4)

#condition - use m.abs2 or m.abs3 for versions of ABS that
# have continuous first and second derivatives
xabs = [[m.Intermediate(m.abs3(x[i,j]-O[i,j])) for j in range(n)] for i in range(n)]
# matrix summation
msum = m.sum([m.sum([xabs[i][j] for i in range(n)]) for j in range(n)])
m.Equation(msum<=2*s)

for i in range(n):
    m.Equation(m.sum(x[:,i])==1)

for i in range(n):
    m.Equation(m.sum(x[i,:])==1)

#Objective - can add multiple objectives
e = d.reshape(n)
for i in range(n):
    for j in range(n):
        m.Obj(a[i]*x[i,j]*e[j])
#m.Obj(m.sum(a.T*m.sum(x*d.reshape(1,n))))

#Set global options
m.options.IMODE = 3

#Solve simulation
m.solve(disp=True,debug=True)

Recommended topics

Hot tags