using GEKKO to optimize with a large number of variables, bounds and constraints

Question

I'm tyring to solve an optimizatio problem using GEKKO. The original problem is one where you have a linear objective function with thousands of variables, constraints (some of which are non-linear) and boundaries. I was kindly suggested to use GEKKO, and since I don't fully understand the mechanics of it I'm having some problem in implementing it. I keep getting object of type 'int' has no len() error. I've simplified the problem so now it's only twenty varialbes and there's no constraints but 10 bounds that have to be respected. This way I think I get to isolate and pindown the source of error. As follows is the code:

import pandas as pd
import numpy as np
import copy as cp
from gekko import GEKKO  
from scipy.optimize import minimize

df_uni = pd.read_csv(r'D:\US\IIT\lectures_textbooks\QIS\week_4\universe_sets.csv')
df_uni["begin_wt"]=np.zeros(10)

#df_holding=df_uni
df_holding=cp.deepcopy(df_uni)

x_holding=np.array([0.1,0.25,0.2,0.2,0.05,0.05,0.1,0.15,-0.05,-0.05])

df_holding["begin_wt"]=x_holding

df_holding.loc[df_holding['begin_wt'] >0, "up_limit"] = df_holding['begin_wt']
df_holding.loc[df_holding['begin_wt'] <0, "up_limit"] = 0

Alpha_pickup=3
df_holding.loc[df_holding['begin_wt'] >0, "Alpha"] = df_holding['Alpha']+Alpha_pickup
df_holding.loc[df_holding['begin_wt'] <0, "Alpha"] = df_holding['Alpha']-Alpha_pickup


df_holding.loc[df_holding['begin_wt'] >0, "low_limit"] = 0
df_holding.loc[df_holding['begin_wt'] <0,"low_limit"]=df_holding['begin_wt']


df_holding=df_holding.drop("begin_w",axis=1)
df_uni=df_uni.drop("begin_w",axis=1)

sect_offset=0.1
lncap_offset=0.1

sect1=sum(df_uni.loc[df_holding['Sector_1'] ==1]['ben_wt'])
sect2=sum(df_uni.loc[df_holding['Sector_2'] ==1]['ben_wt'])

lncap1=sum(df_uni.loc[df_holding['Sector_1'] ==1]['lncap'])
lncap2=sum(df_uni.loc[df_holding['Sector_2'] ==1]['lncap'])


list_uni_alpha=list(df_uni['Alpha'])
list_holding_alpha=list(df_holding['Alpha'])
bind_list_alpha=list_uni_alpha+list_holding_alpha

#x=[1 for i in range(20)]

def objective(x):
    l=0
    sum_of_Alpha=0
    for i in bind_list_alpha:
        sum_of_Alpha=sum_of_Alpha+x[l]*i
        print(sum_of_Alpha)
        l=l+1
    return sum_of_Alpha
# constraints always writing them in terms of f(x)>0
# consolidated weights are bound 
# security offsets
uni_begin=list(df_uni['begin_wt'])
holding_begin=list(df_holding['begin_wt'])
#initial guess
ig=cp.deepcopy(uni_begin+holding_begin)

m=GEKKO()
x = m.Array(m.Var,(20)) 
x_holding=x[10:20]
i=0
#bounds
for xi in x_holding:

    xi.value = x_holding[i]
    xi.lower = df_holding['low_limit'][i]
    xi.upper = df_holding['up_limit'][i]
    i = i + 1

m.Obj(objective(x))
m.solve()
print(x)

Forgive me for including the block of code that doesn't seem to be relevant. But to give context to thoso who are familar with portfolio construction, I'm actually trying to construct a portfolio of stocks. The obejctve function is the linear combination of the stocks' alphas. "The holding" means the stock you're currently holding while "the universe" is the large pool of stocks you get to invest in. I'm doing active management so I tend to overweigh stocks whose prospect I think are good and underweigh those prosepct I don't think are good. But of course I don't want my portfolio to look vastly different than the benchmark because that would make the portfolio bear a lot of systematic risk. Hence those constraints you'll see towards the end of the code. I've been looking for an optimizer that can accommodate constraints written in the form of aX=b, where both a and b are array-like and X is a matrix. But for now, I think this particular optimizer will do me just as good!

Thank you!

score 3 · Answer 1 · answered Dec 20 '19 at 16:19

3

You may have a problem with the line xi.value = x_holding[i] where gekko needs an initial guess value that is a number, not a gekko variable. Here is a simplified version of your problem:

from gekko import GEKKO 
import numpy as np

def objective(x):
    return m.sum(x)

m=GEKKO()
x = m.Array(m.Var,(20)) 
for i,xi in enumerate(x[0:10]):
    xi.value = 0.5
    xi.lower = 0
    xi.upper = i
for i,xi in enumerate(x[10:]):
    xi.value = 0.5
    xi.lower = 0
    xi.upper = i

m.Maximize(objective(x))
m.solve()
print(x)

This gives a solution [0,1,...,8,9,0,1,...,8,9] because all variables go to upper limits with m.Maximize. They all go to the minimum with m.Minimize. I agree with the other answer that it will be more efficient to solve with the m.axb and m.qobj functions with sparse matrices for your portfolio optimization, especially for a large-scale problem.

answered Dec 20 '19 at 16:19

John Hedengren

12,068
1
21
25

Thank you! I don't understand how the way I intialized the guess is wrong though. I literally copied your demo code from another thread here, https://stackoverflow.com/questions/52944970/how-to-use-arrays-in-gekko-optimizer-for-python I also notice that in this reply you just provided, you did not even initialize the guess. Lastly, I have once again spent some time browsing through all the examples suppliled in the website of gekko and just couldn't find an example having a large number of constraints of both non-linear and linear type. Would you nudge me in the right direction? – Chen Lizi Dec 21 '19 at 01:27
The guess for x[i] is initialized after the variable is created with `xi.value = 0.5`. You can also initialize when you declare the variable by using an optional argument to set `value = 0.5`. There are additional examples here: https://github.com/BYU-PRISM/GEKKO/tree/master/examples although most of the problems are examples, not the full-scale problems such as https://github.com/BYU-PRISM/USTAR-Artificial-Lift and https://github.com/BYU-PRISM/hale-trajectory that have 10,000+ variables and nonlinear equations. – John Hedengren Dec 21 '19 at 05:33
Thank you! I'm looking into the three links you posted here. – Chen Lizi Dec 22 '19 at 02:39
Would you care to show what the optimization problems they try to solve are in here? https://github.com/BYU-PRISM/GEKKO/tree/master/examples As I was trying to understand how GEKKO works by looking at any of these examples, I found no explanation or statement of the optimization problem it tried to solve. Is there a book in which you can easily find these optimization? Thank you again! – Chen Lizi Dec 23 '19 at 01:55
Many of these are on the dynamic optimization course at https://apmonitor.com/do or in the APMonitor documentation at https://apmonitor.com/wiki – John Hedengren Dec 23 '19 at 04:14

score 1 · Answer 2 · answered Dec 20 '19 at 15:29

Gekko has functions to load dense or sparse matrices for linear programming problems of the form:

min c x
s.t. A1 x = b1
     A2 x < b2

If you have a very large-scale problem with many zeros in the matrix then the sparse form may be most efficient. Here is the way you are writing models equations:

from gekko import GEKKO
m = GEKKO()
x1 = m.Var(lb=0, ub=5) # Product 1
x2 = m.Var(lb=0, ub=4) # Product 2
m.Maximize(100*x1+125*x2) # Profit function
m.Equation(3*x1+6*x2<=30) # Units of A
m.Equation(8*x1+4*x2<=44) # Units of B
m.solve(disp=False)
p1 = x1.value[0]; p2 = x2.value[0]
print ('Product 1 (x1): ' + str(p1))
print ('Product 2 (x2): ' + str(p2))
print ('Profit        : ' + str(100*p1+125*p2))

If you would like to use the built-in linear equations and quadratic objective models of Gekko, in dense matrix form it is:

from gekko import GEKKO
m = GEKKO(remote=False)
c = [100, 125]
A = [[3, 6], [8, 4]]
b = [30, 44]
x = m.qobj(c,otype='max')
m.axb(A,b,x=x,etype='<')
x[0].lower=0; x[0].upper=5
x[1].lower=0; x[1].upper=4
m.options.solver = 1
m.solve(disp=True)
print ('Product 1 (x1): ' + str(x[0].value[0]))
print ('Product 2 (x2): ' + str(x[1].value[0]))
print ('Profit        : ' + str(m.options.objfcnval))

In sparse matrix form it is:

# solve with GEKKO and sparse matrices
import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
# [[row indices],[column indices],[values]]
A_sparse = [[1,1,2,2],[1,2,1,2],[3,6,8,4]]
# [[row indices],[values]]
b_sparse = [[1,2],[30,44]]
x = m.axb(A_sparse,b_sparse,etype='<',sparse=True)
# [[row indices],[values]]
c_sparse = [[1,2],[100,125]]
m.qobj(c_sparse,x=x,otype='max',sparse=True)
x[0].lower=0; x[0].upper=5
x[1].lower=0; x[1].upper=4
m.solve(disp=True)
print(m.options.OBJFCNVAL)
print('x: ' + str(x))

The sparse matrices are stored in coordinate list (COO) form with [rows,columns,values]. I prefer matrices in compressed sparse row (CSR) form but COO is okay if the problem is not massive and approaching the memory limitations of your computer. You may also want to look at solvers CPLEX, Gurobi, or Xpress Mosel if your problem is linear because these are dedicated linear solvers versus Gekko that uses Mixed Integer Nonlinear Programming solvers. They should give the same answers but the Mixed Integer Linear Programming solvers will be faster.

Thank you! But my real problem is one with a linear objective function and thousands of constraints of both linear types and nonlinear types, which is why I wonder if there is a way of using m.Equation and the format of m.axb together, one for nonlinear and the other for linear? — Chen Lizi, Dec 21 '19 at 01:31
Yes, you can include the m.axb m.qobj models with other equations and objective functions. — John Hedengren, Dec 23 '19 at 04:16
@John Hedengren Thank you very much for your reply! You're saying I could actually mix them up. I'm going to give it a try. I'm genuinely surprized though that there don't seem to be as many non-linear optimization packages catering to real world portfolio construction. I've just talked to my supervisor and he told me he used to do everything either in SAS or in C++ or in Matlab. Btw, when I say real world construction, I mean the way they tend to construct a portfolio of equity trying to beat a benchmark. — Chen Lizi, Dec 24 '19 at 07:28
Nonlinear optimization is somewhat unique. Many use linear models because they are much faster for large-scale problems. One of the unique things about Gekko is the ability to mix linear / nonlinear with time-series empirical or fundamentals-based equations. — John Hedengren, Dec 25 '19 at 20:54

using GEKKO to optimize with a large number of variables, bounds and constraints

2 Answers2