Optimizing Summation function with bounds from binary selector - GEKKO

Question

I am trying to optimize a summation function. From a similar question, I found half an answer for my problem(I have made my objective function similar to the one in that question).

Problem :

Attempt :

from gekko import GEKKO
import numpy as np
from scipy.special import factorial

m = GEKKO(remote=False)
x = m.sos1([1,2,3])
yb = m.Array(m.Var,3,lb=0,ub=1,integer=True)
m.Equation(m.sum(yb)==1)
y = m.sum([yb[i]*(i+1) for i in range(3)])
m.Equation(m.sum([yb[i]*(i+1) for i in range(3)])<2)
yf = factorial(np.linspace(0,len(yb),len(yb)+1))
obj = m.exp(x)
for j in range(0,len(yb)):
    obj += x**j/yf[j]
    m.Maximize(obj*yb[j])
m.solve()
print('x='+str(x.value[0]))
print('y='+str(y.value[0]))
print('Objective='+str(-m.options.objfcnval))

This works but it only works when I add m.Maximize inside the for loop. Running the same without the m.Maximize inside the loop produces wrong result(like below). The problem is the binary selector yb. When yb = [0,1,0] => y=2, the summation bound would be i=0 to 2. But I am not able to get this behaviour. I cannot run the for loop with range(y) because it is a GKVariable and also it changes the equation size as variable values change.

obj = m.exp(x) + x**0/yf[0]
for j in range(1,len(yb)+1):
    obj += x**j/yf[j]*yb[j-1] # It non-zero only when yb[j-1]=1 and misses the values before that.
#Without it I'm computing the entire bound everytime and changing yb has no impact.
m.Maximize(obj)

In the working solution, how does Gekko give me the correct result? How to reproduce this behaviour, when the function is a constraint as below (i.e. without using m.maximize inside the for loop)? Thank you.

Edit: (Modified the above python code). The below figure illustrates the problem. Consider yb=[0,1,0] => y=2. So the summation should to stop at i=0 to 1

I can do this easily if x and yb were normal lists then I could get the expected results from below.

import numpy as np
from scipy.special import factorial

x = 3
yb = [0,1,0]
y = np.sum([yb[i]*(i+1) for i in range(3)])
yf = factorial(np.linspace(0,len(yb),len(yb)+1))
obj = np.exp(x)
for j in range(0,y):
    obj += x**j/yf[j] 
print(obj)

If I don't want the final equation size to change with the variable, I can do the below.

x = 3
yb = [0,1,0]
y = np.sum([yb[i]*(i+1) for i in range(3)])
yf = factorial(np.linspace(0,len(yb),len(yb)+1))
obj = np.exp(x)
for i in range(0,y):
    yb[i] = 1
for j in range(0,len(yb)):
    obj += (x**j/yf[j]) * yb[j]
print(obj)

I cannot do both the option with Gekko because x and yb are GKVariable (design variable) and cannot be interpreted as an integer (i.e) range(y) doesn't work. Is there any workaround for this ? Thank you.

Answer 1

As there aren't any responses, I'm sharing the workaround that I came up with. It probably is not an efficient way of doing it but it works.

m = GEKKO(remote=False)
x = m.sos1([1,2,3])
yb = m.Array(m.Var,3,lb=0,ub=1,integer=True)
m.Equation(m.sum(yb)==1)
y = m.sum([yb[i]*(i+1) for i in range(3)])
#m.Equation(m.sum([yb[i]*(i+1) for i in range(3)])<2)
yf = factorial(np.linspace(0,len(yb),len(yb)+1))
obj = m.exp(x) + x**0/yf[0]
obj_temp=[]
for j in range(1,len(yb)+1):
    obj_temp.append(x**j/yf[j])
    for k in range(j):
        obj+= obj_temp[k] * yb[j-1]
m.Maximize(obj)        
m.solve()
print('x='+str(x.value[0]))
print('y='+str(y.value[0]))
print('Objective='+str(-m.options.objfcnval))

However, I don't know how Gekko replicates this additional for loop internally and why it only works when the summation function is an objective (m.Maximze) and fails when it is a constraint(m.Equation).