I cannot get the optimal solution properly

Question

#pip install gekko
from gekko import GEKKO
import numpy as np
import pandas as pd 

lista1 = [0,1,2,3,4,3,2,1,2]
lista2 = [1,0,1,2,3,3,2,2,1]
lista3 = [2,1,0,1,2,3,3,3,2]
lista4 = [3,2,1,0,1,2,1,2,2]
lista5 = [4,3,2,1,0,1,2,3,3]
lista6 = [3,3,3,2,1,0,1,2,2]
lista7 = [2,2,3,1,2,1,0,1,1]
lista8 = [1,2,3,2,3,2,1,0,2]
lista9 = [2,1,2,2,3,2,1,2,0]
df_distancias = pd.DataFrame(list(zip(lista1,lista2,lista3,lista4,lista5,lista6,lista7,lista8,lista9)),columns=['col1','col2','col3','col4','col5','col6','col7','col8','col9'])
distance_matrix = df_distancias.to_numpy()
m = GEKKO() 

## Set variables as a matrix shaped (9,6)
q = m.Array(m.Var,(9,6),lb=0,ub=1,integer=True) ## Get matrix product  

producto_matricial = np.dot(distance_matrix,q) ## Set objective function

def objective_function(q):
    rdo = sum(q.flatten() * producto_matricial.flatten())
    return rdo 

## Set the params to constraint our variables matrix

vector_columna = m.Param(value=[3,5,4,4,4,3])
vector_fila =  m.Param(value=[4,3,2,4,2,2,2,2,2])
vector_fila ## Set columns and rows constraints

for j in range(6):
    m.Equation(np.sum(q[:,j])==vector_columna[j]) 
for i in range(9):
    m.Equation(np.sum(q[i,:])==vector_fila[i]) 

## Minimize objective function
m.Minimize(objective_function(q)) ## Select APOPT solver because we´re looking for binary solution 

m.options.SOLVER = 1 m.solve(disp=False) 

And the error I get is :

Exception Traceback (most recent call last)<ipython-input-18-c3170ab08189> in <module>
----> 1 m.solve(disp=False) /usr/local/lib/python3.9/dist-packages/gekko/gk_write_files.py in _write_csv(self)
    159                 length = np.size(np.array(vp.value).flatten())
    160                 if self.options.IMODE in (1,3) and length > 1:
--> 161                     raise Exception('This steady-state IMODE only allows scalar values.')
    162                 elif self.options.IMODE == 2 and length == 1:
    163                     #in MPU, the first vp checked could be a scalar value (FV, param, var initial guess
Exception: This steady-state IMODE only allows scalar valueshas context menu

My objective function consists of two steps:

1. First step, matrix product of a symmetric matrix by the variable matrix
2. The next step is the sum of the scalar product of the result of the first step through the array of variables

Answer 1

The error is because vector_columna and vector_fila need to be scalars instead of arrays as Gekko parameters. It is easier to just leave them as a Python list for this problem.

vector_columna = [3,5,4,4,4,3]
vector_fila =  [4,3,2,4,2,2,2,2,2]

Using m.sum() instead of np.sum() also improves the efficiency of the solution. There were a few syntax errors that are now fixed.

import numpy as np
import pandas as pd 
from gekko import GEKKO

lista1 = [0,1,2,3,4,3,2,1,2]
lista2 = [1,0,1,2,3,3,2,2,1]
lista3 = [2,1,0,1,2,3,3,3,2]
lista4 = [3,2,1,0,1,2,1,2,2]
lista5 = [4,3,2,1,0,1,2,3,3]
lista6 = [3,3,3,2,1,0,1,2,2]
lista7 = [2,2,3,1,2,1,0,1,1]
lista8 = [1,2,3,2,3,2,1,0,2]
lista9 = [2,1,2,2,3,2,1,2,0]
df_distancias = pd.DataFrame(list(zip(lista1,lista2,lista3,lista4,
                                      lista5,lista6,lista7,lista8,
                                      lista9)),
                             columns=['col1','col2','col3','col4',
                                      'col5','col6','col7','col8',
                                      'col9'])
distance_matrix = df_distancias.to_numpy()

m = GEKKO() 

## Set variables as a matrix shaped (9,6)
q = m.Array(m.Var,(9,6),lb=0,ub=1,integer=True)

producto_matricial = np.dot(distance_matrix,q)

def objective_function(q):
    rdo = sum(q.flatten() * producto_matricial.flatten())
    return rdo 

## Set the params to constraint our variables matrix
vector_columna = [3,5,4,4,4,3]
vector_fila =  [4,3,2,4,2,2,2,2,2]
vector_fila ## Set columns and rows constraints

for j in range(6):
    m.Equation(m.sum(q[:,j])==vector_columna[j]) 
for i in range(9):
    m.Equation(m.sum(q[i,:])==vector_fila[i]) 

## Minimize objective function
m.Minimize(objective_function(q))

m.options.SOLVER = 1
m.solve(disp=False) 

print(q)

Solution is non-zero now:

[[[1.0] [0.0] [1.0] [1.0] [0.0] [1.0]]
 [[0.0] [0.0] [1.0] [1.0] [1.0] [0.0]]
 [[0.0] [0.0] [0.0] [1.0] [1.0] [0.0]]
 [[0.0] [1.0] [1.0] [1.0] [1.0] [0.0]]
 [[0.0] [1.0] [0.0] [0.0] [1.0] [0.0]]
 [[0.0] [1.0] [0.0] [0.0] [0.0] [1.0]]
 [[1.0] [1.0] [0.0] [0.0] [0.0] [0.0]]
 [[1.0] [0.0] [0.0] [0.0] [0.0] [1.0]]
 [[0.0] [1.0] [1.0] [0.0] [0.0] [0.0]]]