Search code examples
gekko

GEKKO: Array size as a model variable


I'm quite new to Gekko. Is it possible to vary the size of a model array as part of an optimization? I am running a simple problem where various numbers of torsional springs engage at different angles, and I would like to allow the model to change the number of engagement angles. Each spring has several component variables, which I am also attempting to define as arrays of variables. However, the size definition of the array theta_engage, below, has not accepted int(n_engage.value). I get the following error:

TypeError: int() argument must be a string, a bytes-like object or a number, not 'GK_Value'

Relevant code:

n_engage = m.Var(2, lb=1, ub=10, integer=True)

theta_engage = m.Array(m.Var, (int(n_engage.value)))
theta_engage[0].value = 0.0
theta_engage[0].lower = 0.0
theta_engage[0].upper = 85.0
theta_engage[1].value = 15.0
theta_engage[1].lower = 0.0
theta_engage[1].upper = 85.0

If I try to define the size of theta_engage only by n_engage.value, I get this error:

TypeError: expected sequence object with len >= 0 or a single integer

I suppose I could define the array at the maximum size I am willing to accept and allow the number of springs to have a lower bound of 0, but I would have to enforce a minimum number of total springs somehow in the constraints. If Gekko is capable of varying the size of the arrays this way it seems to me the more elegant solution.

Any help is much appreciated.


Solution

  • The problem structure can't be changed iteration-to-iteration. However, it is easy to define a binary variable b that either activates or deactivates those parts of the model that should be included or excluded.

    from gekko import GEKKO
    import numpy as np
    m = GEKKO()
    # number of springs
    n = 10
    # number of engaged springs (1-10)
    nb = m.Var(2, lb=1, ub=n, integer=True)
    # engaged springs (binary, 0-1)
    b = m.Array(m.Var,n,lb=0,ub=1,integer=True)
    # angle of engaged springs
    θ = m.Array(m.Param,n,lb=0,ub=85)
    # initialize values
    t0 = [0,15,20,25,30,15,30,25,10,50]
    for i,ti in enumerate(t0):
        θ[i].value = ti
    # contributing spring forces
    F = [m.Intermediate(b[i]*m.cos((np.pi/180.0)*θ[i])) \
         for i in range(10)]
    # force constraint
    m.Equation(m.sum(F)>=3)
    # engaged springs
    m.Equation(nb==m.sum(b))
    # minimize engaged springs
    m.Minimize(nb)
    # optimize with APOPT solver
    m.options.SOLVER=1
    m.solve()
    # print solution
    print(b)
    

    This gives a solution in 0.079 sec that springs 1, 3, 9, and 10 should be engaged. It selects the minimum number of springs (4) to achieve the required force that is equivalent to 3 springs at 0 angle.

     Successful solution
     ---------------------------------------------------
     Solver         :  APOPT (v1.0)
     Solution time  :   7.959999999729916E-002 sec
     Objective      :    4.00000000000000     
     Successful solution
     ---------------------------------------------------
    [[1.0] [0.0] [1.0] [0.0] [0.0] [0.0] [0.0] [0.0] [1.0] [1.0]]