What is the proper way to use factorial and range functions when doing optimization with Gekko?

Question

I have been trying to solve a mixed integer nonlinear programming problem using Gekko but I am having trouble with the usage of range and factorial functions.

My model has one variable and it aims to minimize c while keeping Wq on a certain limit. When I execute the code below;

m = GEKKO() 
m.options.SOLVER=1 
def ineq1(c):
    arr_rate = 60
    ser_rate = 25
    p = arr_rate/(ser_rate)
    eps = 0
    for m in range(c):
        eps += p**m/factorial(m)
    Po = 1/(eps+p**c/(factorial(c)*(1-p/(c))))
    Lq = Po*p**(c+1)/(c*factorial(c)*(1-p/c)**2)
    Wq = Lq/arr_rate
    return Wq
x1 = m.Var(value=2,lb=1,ub=10,integer=True)
m.Equation(ineq1(x1)<=0.005)
m.Obj(x1) 
m.solve(disp=False)

I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-427-6b79ae34d974> in <module>
     16 x1 = m.Var(value=1,lb=1,ub=5,integer=True)
     17 # Equations
---> 18 m.Equation(ineq1(x1)<=0.5)
     19 m.Obj(x1) # Objective
     20 m.solve(disp=False) # Solve

<ipython-input-427-6b79ae34d974> in ineq1(c1)
      8     p = arr_rate/(ser_rate)
      9     eps = 0
---> 10     for m in range(c):
     11         eps += p**m/factorial(m)
     12     Po = 1/(eps+p**c/(factorial(c)*(1-p/(c))))

TypeError: 'GKVariable' object cannot be interpreted as an integer

Apparently, Gekko does not want me to force the variable to be an integer but I have to use range and factorial functions in order to get my model properly worked. I would appreciate any kind of advice, thanks.

Answer 1

This isn't a good application for Mixed Integer Programming because some of the trial solutions are non-integer. How about just evaluating x1 with higher numbers until the constraint is satisfied:

from math import factorial

def ineq1(c):
    arr_rate = 60
    ser_rate = 25
    p = arr_rate/(ser_rate)
    eps = 0
    for m in range(c):
        eps += p**m/factorial(m)
    Po = 1/(eps+p**c/(factorial(c)*(1-p/(c))))
    Lq = Po*p**(c+1)/(c*factorial(c)*(1-p/c)**2)
    Wq = Lq/arr_rate
    return Wq

for x1 in range(1,15):
    if ineq1(x1)<=0.005:
        print('Constraint ineq1(x1)<0.005: ' + '{0:10.5f}'.format(ineq1(x1)) \
              + ' satisfied with minimum x1=' + str(x1))
        break

Stop at the first value that satisfies the constraint and report it as the optimal solution.

Constraint ineq1(x1)<0.005:   -0.06857 satisfied with minimum x1=1