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!

Answer 1

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.

Answer 2

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.