How to initialize an 3D array variable in Gekko?

Question

I'm trying to solve a step from a three-dimensional master timetabling model, which involvs periods(5), courses(19) and locations(8).

So I have a problem to initialize these variables with an 3D array in Gekko. Without this initialization the algorithm doesn't converge, after more than 15 minutes run and 1000 iterations.

When I try initialize, this error appears:

" raise Exception(response) Exception: @error: Equation Definition Equation without an equality (=) or inequality (>,<) true STOPPING... "

How can I fix this problem? Follows a version of my code:

import numpy as np
from gekko import GEKKO

# Input data

# Schedule of periods and courses
sched = np.array([ [0,  1,  0,  0,  1], [0, 0,  1,  1,  0], [0, 0,  1,  1,  0], \
[0, 0,  0,  0,  1], [1, 0,  0,  0,  1], [0, 0,  0,  1,  1], [0, 1,  1,  0,  0], \
[1, 0,  0,  1,  0], [0, 1,  0,  0,  1], [1, 1,  0,  0,  0], [0, 1,  1,  0,  0], \
[0, 1,  1,  0,  0], [1, 0,  0,  1,  0], [1, 0,  0,  1,  0], [0, 0,  1,  0,  1], \
[1, 0,  1,  0,  0], [0, 1,  0,  1,  0], [0, 0,  1,  1,  0], [0, 1,  0,  0,  1] ], dtype=np.int64)

# Initial allocation of all periods, courses and locations
alloc=np.array([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,\
                0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int64)

# Number of students enrolled in each course
enrol = np.array([ 60, 60,  60,   40,   40,  110,  120,   50,   60,    55,    50, \
                   55,    40,    64,    72, 50,    50,    55,    55], dtype=np.float64)

# Capacity of each location (classroom)
capac = np.array([ 60,   60,  120,   60,   80,   60,   60,   65], dtype=np.float64)

# Total costs of using each location
costs = np.array([ 9017.12, 9017.12, 12050.24, 9017.12, 9413.68,  9017.12, \
                   9017.12, 9188.96 ])

# Estimated cost of each location by period and student
ecost = np.repeat(np.array([[costs[i]*pow(enrol[j]*5,-1) for j in range(19)] for i in range(8)]), 5)

# The model construction

m = GEKKO()

# Constant arrays
x = m.Array(m.Const,(19,5))
y = m.Array(m.Const,(8,19,5))
N = m.Array(m.Const,(19))
C = m.Array(m.Const,(8))
Ec = m.Array(m.Const,(8,19,5))
Ecy = m.Array(m.Const,(8,19,5))
Alt = m.Array(m.Const,(8,19,5))

for k in range(5):
   for j in range(19):
      N[j] = enrol[j]
      x[j,k] = sched[j,k]
      for i in range(8):
         C[i] = capac[i]
         Ec[i,j,k] = ecost[k+j*5+i*19*5]
         y[i,j,k] = alloc[k+j*5+i*19*5]
         Ecy[i,j,k] = Ec[i,j,k]*y[i,j,k]
         if sched[j,k]==1:
            Alt[i,j,np.where(sched[j,:]==1)[0][0]]=-sched[j,k]*(1-sum(sched[j,:]))
            if sum(sched[j,:])==2:
               Alt[i,j,np.where(sched[j,:]==1)[0][1]]=sched[j,k]*(1-sum(sched[j,:]))
         else:
            Alt[i,j,k]=0

# Initialize the variable z with the initial value y: 
# These commented approaches produce the error.
z = m.Array(m.Var,(8,19,5),lb=0,ub=1,integer=True) 
#for i in range(8):
#   for j in range(19):
#      for k in range(5):
#         z[i,j,k] = y[i,j,k]
# nor
#z = m.Array(m.Var,(8,19,5),value=y,lb=0,ub=1,integer=True)  

# Intermediate equations
Ecz = m.Array(m.Var,(8,19,5),lb=0)
Altz = m.Array(m.Var,(8,19))
for i in range(8):
   for j in range(19):
      Altz[i,j]=m.Intermediate(m.sum(Alt[i,j,:]*z[i,j,:]))
      for k in range(5):
         Ecz[i,j,k]=m.Intermediate(Ec[i,j,k]*z[i,j,k])

# Constraints
m.Equation(m.sum(m.sum(m.sum(Ecz)))<=m.sum(m.sum(m.sum(Ecy))))
for j in range(19):
   for k in range(5):
      m.Equation(m.sum(z[:,j,k])==x[j,k])
for i in range(8):
   for k in range(5):
      m.Equation(m.sum(z[i,:,k])==m.sum(y[i,:,k]))
for i in range(8): 
   for j in range(19):
      m.Equation(m.sum((C[i]/N[j]-x[j,:])*z[i,j,:])>=0)

# Objective: to minimize the quantity of courses allocated in different locations      
# Example: with the solution y, I have 12 courses in different locations in the periods 
# print(sum([sum(Alt[i,j,:]*y[i,j,:])**2 for j in range(19) for i in range(8)])/2) 
for i in range(8): 
   for j in range(19):
      m.Obj(Altz[i,j]**2/2)

# Options and final results
m.options.SOLVER=1
m.options.IMODE=2
m.solve()
print(z)
print(m.options.OBJFCNVAL)

Note: My original problem has 20 periods, 171 courses, and 18 locations.

Answer 1

Use z[i,j,k].value = y[i,j,k] to give an initial guess for z. Using z[i,j,k] = y[i,j,k] redefines z entries as floating point numbers instead of gekko variable types.

One other issue is that the variables Ecz and Altz are defined as Variables as m.Var and then overridden as Intermediates. Instead, try allocating them and assigning them as intermediates:

Ecz = np.empty((8,19,5),dtype=object)
Altz = np.empty((8,19),dtype=object)

Use flatten() to simplify the summation of all elements of the 3 dimensional array.

m.Equation(m.sum(Ecz.flatten())<=sum(Ecy.flatten()))

The constant arrays can be defined as numpy arrays to avoid additional symbolic processing by Gekko. This speeds up the model compile time but has no effect on the final solution.

x   = np.empty((19,5))
y   = np.empty((8,19,5))
N   = np.empty((19))
C   = np.empty((8))
Ec  = np.empty((8,19,5))
Ecy = np.empty((8,19,5))
Alt = np.empty((8,19,5))

The IMODE should be 3 for optimization. IMODE=2 is for parameter regression. IMODE=2 should also work for this problem but 3 is the correct option because you aren't trying to fit to data.

m.options.IMODE=3

Try using IPOPT to obtain an initial non-integer solution and then use APOPT to find an integer solution.

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 1']

Mixed Integer Nonlinear Programming (MINLP) problems can be challenging to solve so you may need to use some of the solver options to speed up the solution. Try minlp_branch_method 1 to help the solver find an initial integer solution to do better pruning. The gap tolerance can also help to speed up the solution if a sub-optimal solution is okay. Below is the complete script. Consider using remote=False to run locally instead of using the public servers, especially for large optimization problems.

import numpy as np
from gekko import GEKKO

# Input data

# Schedule of periods and courses
sched = np.array([ [0,  1,  0,  0,  1], [0, 0,  1,  1,  0], [0, 0,  1,  1,  0], \
[0, 0,  0,  0,  1], [1, 0,  0,  0,  1], [0, 0,  0,  1,  1], [0, 1,  1,  0,  0], \
[1, 0,  0,  1,  0], [0, 1,  0,  0,  1], [1, 1,  0,  0,  0], [0, 1,  1,  0,  0], \
[0, 1,  1,  0,  0], [1, 0,  0,  1,  0], [1, 0,  0,  1,  0], [0, 0,  1,  0,  1], \
[1, 0,  1,  0,  0], [0, 1,  0,  1,  0], [0, 0,  1,  1,  0], [0, 1,  0,  0,  1] ], dtype=np.int64)

# Initial allocation of all periods, courses and locations
alloc=np.array([0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,\
                0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,\
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], dtype=np.int64)

# Number of students enrolled in each course
enrol = np.array([ 60, 60,  60,   40,   40,  110,  120,   50,   60,    55,    50, \
                   55,    40,    64,    72, 50,    50,    55,    55], dtype=np.float64)

# Capacity of each location (classroom)
capac = np.array([ 60,   60,  120,   60,   80,   60,   60,   65], dtype=np.float64)

# Total costs of using each location
costs = np.array([ 9017.12, 9017.12, 12050.24, 9017.12, 9413.68,  9017.12, \
                   9017.12, 9188.96 ])

# Estimated cost of each location by period and student
ecost = np.repeat(np.array([[costs[i]*pow(enrol[j]*5,-1) for j in range(19)] for i in range(8)]), 5)

# The model construction

m = GEKKO(remote=True)

# Constant arrays
x   = np.empty((19,5))
y   = np.empty((8,19,5))
N   = np.empty((19))
C   = np.empty((8))
Ec  = np.empty((8,19,5))
Ecy = np.empty((8,19,5))
Alt = np.empty((8,19,5))

for k in range(5):
   for j in range(19):
      N[j] = enrol[j]
      x[j,k] = sched[j,k]
      for i in range(8):
         C[i] = capac[i]
         Ec[i,j,k] = ecost[k+j*5+i*19*5]
         y[i,j,k] = alloc[k+j*5+i*19*5]
         Ecy[i,j,k] = Ec[i,j,k]*y[i,j,k]
         if sched[j,k]==1:
            Alt[i,j,np.where(sched[j,:]==1)[0][0]]=-sched[j,k]*(1-sum(sched[j,:]))
            if sum(sched[j,:])==2:
               Alt[i,j,np.where(sched[j,:]==1)[0][1]]=sched[j,k]*(1-sum(sched[j,:]))
         else:
            Alt[i,j,k]=0

# Initialize the variable z with the initial value y: 
# These commented approaches produce the error.
z = m.Array(m.Var,(8,19,5),lb=0,ub=1,integer=True) 
for i in range(8):
   for j in range(19):
      for k in range(5):
         z[i,j,k].value = y[i,j,k]
# nor
#z = m.Array(m.Var,(8,19,5),value=y,lb=0,ub=1,integer=True)  

# Intermediate equations
Ecz = np.empty((8,19,5),dtype=object)
Altz = np.empty((8,19),dtype=object)
for i in range(8):
   for j in range(19):
      Altz[i,j]=m.Intermediate(m.sum(Alt[i,j,:]*z[i,j,:]))
      for k in range(5):
         Ecz[i,j,k]=m.Intermediate(Ec[i,j,k]*z[i,j,k])

# Constraints
m.Equation(m.sum(Ecz.flatten())<=sum(Ecy.flatten()))
for j in range(19):
   for k in range(5):
      m.Equation(m.sum(z[:,j,k])==x[j,k])
for i in range(8):
   for k in range(5):
      m.Equation(m.sum(z[i,:,k])==m.sum(y[i,:,k]))
for i in range(8): 
   for j in range(19):
      m.Equation(m.sum((C[i]/N[j]-x[j,:])*z[i,j,:])>=0)

# Objective: to minimize the quantity of courses allocated in different locations      
# Example: with the solution y, I have 12 courses in different locations in the periods 
# print(sum([sum(Alt[i,j,:]*y[i,j,:])**2 for j in range(19) for i in range(8)])/2) 
for i in range(8): 
   for j in range(19):
      m.Obj(Altz[i,j]**2/2)

# Options and final results
m.options.IMODE=3

# Initialize with IPOPT
m.options.SOLVER=3
m.solve()

# Integer solution with APOPT
m.options.SOLVER=1

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 1']

m.solve()


print(z)
print(m.options.OBJFCNVAL)