GEKKO: count ocurrences in an 2D array as condition

I have a grid of 9 points in the xy coordinate system with following coordinates (1,1), (1,2), (1,3), (2,1), (2,2), (2,3), (3,1), (3,2), (3,3).

The grid looks like the following:

I would like to solve follwoing non-linear integer optimization problem using GEKKO: where i=1,2,3 and j=1,2,3. In other words, the pairs (xi, yi) can only appear once in each point of the grid.

So far I wrote this piece of code, but I don't know how to define the constraints in GEKKO.

from gekko import GEKKO

m = GEKKO() # create GEKKO model
w13 = ((x1-x3)**2+(y1-y3)**2)**0.5
w57 = ((x5-x7)**2+(y5-y7)**2)**0.5
obj_fun = 2*w13 + 3*w57
m.Minimize(obj_fun)
m.options.SOLVER = 1
m.solve()

Edit: Inspired by John's comment, I've tried out this code

from gekko import GEKKO
import numpy as np
m = GEKKO()
m.options.SOLVER = 1

b = m.Array(m.Var, (3,3), lb=0, ub=1, integer=True)
k = m.Array(m.Var, (3))
k[0].value = 1
k[1].value = 2
k[2].value = 3
x1 = m.sum(b.dot(k))
x3 = m.sum(b.dot(k))
x5 = m.sum(b.dot(k))
x7 = m.sum(b.dot(k))
y1 = m.sum(b.dot(k))
y3 = m.sum(b.dot(k))
y5 = m.sum(b.dot(k))
y7 = m.sum(b.dot(k))

w13 = ((x1-x3)**2+(y1-y3)**2)**0.5
w57 = ((x5-x7)**2+(y5-y7)**2)**0.5
obj_fun = 2*w13 + 3*w57
m.Equation(m.sum(b) <= 1)
m.Equation(m.sum(b) >= 1)
m.Minimize(obj_fun)
m.solve()

,but unfortunately I got an error message:

Traceback (most recent call last):
  File "C:\Users\Arman.Kohnic\AppData\Local\Programs\Python\Python310\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 23, in <module>
  File "C:\Users\...\Local\Programs\Python\Python310\lib\site-packages\gekko\gk_operators.py", line 25, in __len__
    return len(self.value)
  File "C:\Users\...\Local\Programs\Python\Python310\lib\site-packages\gekko\gk_operators.py", line 144, in __len__
    return len(self.value)
TypeError: object of type 'int' has no len()