CODEWITHSUNDEEP
pythonnumpygekko
Thomas Conrad
Thomas Conrad

Reputation: 57

Max equation length limit in Gekko constraint

I am working on optimizing a long list of equations under a constraint. The objective function can be split into several objective functions, but i am unable to split my constraint function in any way that i can find. I have tried to solve this problem in several different ways over the last month, using both sympy and scipy.optimize, but both take ages to run. In this example i have limited the length to one of the shortest where i get the error, but i would like to be able to expand this as much as i want by increasing constraintmonth, nummonthstodo and adding to db1permonth. How can i get around the length limit error:

Exception:  @error: Max Equation Length
 Error with line number:           67
 (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((
 [...lots of numbers here..]
 0/(1.0+(2.718281828459045^(-(((v14-476.0))/(120)))))))))+((-2.02635827)*((200.0
 /(1.0+(2.718281

 APM model error: string >       15000  characters
 Consider breaking up the line into multiple equations

 The may also be due to only using newline character CR
   instead of CR LF (for Windows) or LF (for MacOS/Linux)
 To fix this problem, save APM file with appropriate newline characters

 STOPPING...

In this code:

import numpy as np
from gekko import GEKKO    
from math import e


m = GEKKO()

nummonthstodo = 20
constraintmonth = 15
markets=3
fudgefactor = 1.18
constraint = -100000
ig = np.tile([1169.00,476.00,191.00],(nummonthstodo))
dynamic_scale = np.tile(np.array([100,120,80]),nummonthstodo)
turning_point =np.tile( np.array([1169.00,476.00,191.00]),nummonthstodo)
eo14_scale = np.tile(np.array([72.00,200.00,500.00]),nummonthstodo)
db1life = np.tile(np.array([2131.64415,  927.62553,  376.79813]),nummonthstodo)

db1permonth = np.array([-7.69469125e+01, -8.27515903e+01, -8.23699131e+01,  4.45215469e+02,
        3.69705278e+02,  1.58977477e+02, -3.88706011e+00, -3.02636074e+00,
       -2.79420475e+00, -7.67567831e-01, -3.03178864e+00, -2.79921627e+00,
        6.82114433e+01,  5.11866268e+01,  7.41661986e+01, -3.22847133e+00,
       -2.28929507e+00, -1.28846317e+00, -7.10922109e-01, -2.29483815e+00,
       -1.29158293e+00,  5.33895986e+01,  3.53279803e+01,  4.70884129e+01,
       -2.98307479e+00, -2.02448766e+00, -8.09906719e-01, -8.69398271e-01,
       -2.02635827e+00, -8.10655068e-01,  4.40599327e+01,  2.74298296e+01,
        3.40114832e+01, -2.83729295e+00, -1.90180770e+00, -5.82931770e-01,
       -4.06643364e+01, -2.85105835e+01, -8.73890925e+00,  2.05687539e+02,
        1.15368035e+02,  2.54844135e+01, -1.93419651e+00, -1.17779492e+00,
       -4.36383453e-01, -3.54686842e-01, -1.17891974e+00, -4.36800210e-01,
        3.28820412e+01,  1.73557733e+01,  2.03414999e+01, -1.85294103e+00,
       -1.09751072e+00, -3.47833721e-01, -4.44622167e-01, -1.09856091e+00,
       -3.48166558e-01,  2.90051018e+01,  1.46659871e+01,  1.69382174e+01])



x = m.Array(m.Var,(len(ig)))

# lower bounds
for i,xi in enumerate(x):
    xi.value = ig[i]
    xi.lower = 0
    xi.upper = 1500


def eo14(a,b):
    return np.divide(np.tile(eo14_scale[0:3],int((b-a)/3)),(np.ones(b-a)+e**(-np.divide((x[a:b]-np.tile(turning_point[0:3],int((b-a)/3))),np.tile(dynamic_scale[0:3],int((b-a)/3))))))
def db1(a,b):
    return db1permonth[a:b]
def cac(a,b):
    return x[a:b]

eo14s = [eo14(i*markets,(i+1)*markets) for i in range(nummonthstodo)]
db1s = [db1(i*markets,(i+1)*markets) for i in range(nummonthstodo)]
cacs = [cac(i*markets,(i+1)*markets) for i in range(nummonthstodo)]

db1month = 0

for i in range(constraintmonth):
    for j in range(constraintmonth-i):
        db1month = db1month + np.multiply(db1s[j],eo14s[i])

spend = 0

for i in range(constraintmonth):
        spend = spend + np.multiply(cacs[i],eo14s[i])

alleo14 = eo14(0,len(ig))

m.Equation(m.sum(db1month - spend) >= constraint)

for i in range(len(ig)):
   m.Obj(-(np.multiply(db1life-x,alleo14))[i])

m.solve()

print(x)

Thanks in advance for any help you can give, as a student and a beginner in data science this has really been giving me headaches for some time now.

Upvotes: 2

Views: 609

Answers (1)

John Hedengren
John Hedengren

Reputation: 14346

You can overcome the maximum equation length by declaring your products as a list and then use m.sum(). Otherwise, it builds a very large expression that exceeds the limit of 15,000 characters per equation. You could also use m.Intermediate() variables to reduce the expression size. Here is a new version of your program that solves successfully.

import numpy as np
from gekko import GEKKO    
from math import e

m = GEKKO(remote=False)

nummonthstodo = 20
constraintmonth = 15
markets=3
fudgefactor = 1.18
constraint = -100000
ig = np.tile([1169.00,476.00,191.00],(nummonthstodo))
dynamic_scale = np.tile(np.array([100,120,80]),nummonthstodo)
turning_point =np.tile( np.array([1169.00,476.00,191.00]),nummonthstodo)
eo14_scale = np.tile(np.array([72.00,200.00,500.00]),nummonthstodo)
db1life = np.tile(np.array([2131.64415,  927.62553,  376.79813]),nummonthstodo)

db1permonth = np.array([-7.69469125e+01, -8.27515903e+01, -8.23699131e+01,  4.45215469e+02,
        3.69705278e+02,  1.58977477e+02, -3.88706011e+00, -3.02636074e+00,
       -2.79420475e+00, -7.67567831e-01, -3.03178864e+00, -2.79921627e+00,
        6.82114433e+01,  5.11866268e+01,  7.41661986e+01, -3.22847133e+00,
       -2.28929507e+00, -1.28846317e+00, -7.10922109e-01, -2.29483815e+00,
       -1.29158293e+00,  5.33895986e+01,  3.53279803e+01,  4.70884129e+01,
       -2.98307479e+00, -2.02448766e+00, -8.09906719e-01, -8.69398271e-01,
       -2.02635827e+00, -8.10655068e-01,  4.40599327e+01,  2.74298296e+01,
        3.40114832e+01, -2.83729295e+00, -1.90180770e+00, -5.82931770e-01,
       -4.06643364e+01, -2.85105835e+01, -8.73890925e+00,  2.05687539e+02,
        1.15368035e+02,  2.54844135e+01, -1.93419651e+00, -1.17779492e+00,
       -4.36383453e-01, -3.54686842e-01, -1.17891974e+00, -4.36800210e-01,
        3.28820412e+01,  1.73557733e+01,  2.03414999e+01, -1.85294103e+00,
       -1.09751072e+00, -3.47833721e-01, -4.44622167e-01, -1.09856091e+00,
       -3.48166558e-01,  2.90051018e+01,  1.46659871e+01,  1.69382174e+01])

x = m.Array(m.Var,len(ig),lb=0,ub=1500)
for i,xi in enumerate(x):
    xi.value = ig[i]

def eo14(a,b):
    return np.divide(np.tile(eo14_scale[0:3],int((b-a)/3)),\
                     (np.ones(b-a)+e**(-np.divide((x[a:b]\
                        -np.tile(turning_point[0:3],int((b-a)/3))),\
                            np.tile(dynamic_scale[0:3],int((b-a)/3))))))
def db1(a,b):
    return db1permonth[a:b]
def cac(a,b):
    return x[a:b]

eo14s = [eo14(i*markets,(i+1)*markets) for i in range(nummonthstodo)]
db1s = [db1(i*markets,(i+1)*markets) for i in range(nummonthstodo)]
cacs = [cac(i*markets,(i+1)*markets) for i in range(nummonthstodo)]

prod1 = [db1s[j]*eo14s[i] for j in range(constraintmonth-i) \
                         for i in range(constraintmonth)]
db1month = m.sum(prod1)

prod2 = [cacs[i]*eo14s[i] for i in range(constraintmonth)]
spend = m.sum(prod2)

alleo14 = eo14(0,len(ig))

m.Equation(db1month - spend >= constraint)

for i in range(len(ig)):
   m.Maximize((np.multiply(db1life-x,alleo14))[i])

m.solve()

print(x)

This gives the solver output:

Number of Iterations....: 9

                                   (scaled)                 (unscaled)
Objective...............: -1.7974777841412724e+003  -2.8845129361943263e+006
Dual infeasibility......:  1.3421365394256974e-009   2.1538014234551661e-006
Constraint violation....:  0.0000000000000000e+000   0.0000000000000000e+000
Complementarity.........:  2.1861799117816599e-010   3.5082849379387629e-007
Overall NLP error.......:  1.3421365394256974e-009   2.1538014234551661e-006


Number of objective function evaluations             = 10
Number of objective gradient evaluations             = 10
Number of equality constraint evaluations            = 10
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 10
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 9
Total CPU secs in IPOPT (w/o function evaluations)   =      0.066
Total CPU secs in NLP function evaluations           =      0.001

EXIT: Optimal Solution Found.

 The solution was found.

 The final value of the objective function is  -2884512.9361943263

 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :  0.127 sec
 Objective      :  -2884512.9361943263
 Successful solution
 ---------------------------------------------------

And the solution:

[[1359.5263445] [561.97176237] [203.39330742] [1359.5263445]
 [561.97176237] [203.39330742] [1359.5263445] [561.97176237]
 [203.39330742] [1359.5263445] [561.97176237] [203.39330742]
 [1359.5263445] [561.97176237] [203.39330742] [1359.5263445]
 [561.97176237] [203.39330742] [1359.5263445] [561.97176237]
 [203.39330742] [1359.5263445] [561.97176237] [203.39330742]
 [1359.5263445] [561.97176237] [203.39330742] [1359.5263445]
 [561.97176237] [203.39330742] [1359.5263445] [561.97176237]
 [203.39330742] [1359.5263445] [561.97176237] [203.39330742]
 [1359.5263445] [561.97176237] [203.39330742] [1359.5263445]
 [561.97176237] [203.39330742] [1359.5263445] [561.97176237]
 [203.39330742] [1359.5263445] [561.97176237] [203.39330742]
 [1359.5263445] [561.97176237] [203.39330742] [1359.5263445]
 [561.97176237] [203.39330742] [1359.5263445] [561.97176237]
 [203.39330742] [1359.5263445] [561.97176237] [203.39330742]]

Upvotes: 2

Related Questions