Trying to solve this non linear optimization using GEKKO, getting this error

Question

@Error: setting an array element with a sequence

I am trying to mninimize the downside risk.

I have a two dimensional array of returns shape(1000, 10), and the portfolio starts with $100. Compound that 10 times by each return in a row. Do that for all the rows. Compare that last cell's value for each row with mean of last column's values. Keep the value if it's less than mean or else zero. So we will have an array of (1000, 1). At the end I am finding the standard deviation of that.

Objective is to minimize the standard deviation. Constraints: weights need to be less than 1

the expected return i.e. wt*ret should be equal to a value like 7%. I have to do that for couple of values like 7%, 8% or 10%.

wt = np.array([0.4, 0.3, 0.3])

cov = array([[0.00026566, 0.00016167, 0.00011949],
   [0.00016167, 0.00065866, 0.00021662],
   [0.00011949, 0.00021662, 0.00043748]])

ret =[.098, 0.0620,.0720]
iterations = 10000

return_sim = np.random.multivariate_normal(ret, cov, iterations)

def simulations(wt):

    downside =[]
    fund_ret =np.zeros((1000,10))
    prt_ret = np.dot(return_sim , wt)

    re_ret = np.array(prt_ret).reshape(1000, 10) #10 years

    for m in range(len(re_ret)):
        fund_ret[m][0] = 100 * (1 + re_ret[m][0])  #start with $100
        for n in range(9):

            fund_ret[m][n+1] = fund_ret[m][n]* (1 + re_ret[m][n+1])

    mean = np.mean(fund_ret[:,-1])  #just need the last column and all rows

    for i in range(1000): 
        downside.append(np.maximum((mean - fund_ret[i,-1]), 0))

    return np.std(downside)

b = GEKKO()
w = b.Array(b.Var,3,value=0.33,lb=1e-5, ub=1)
b.Equation(b.sum(w)<=1)
b.Equation(np.dot(w,ret) == .07) 
b.Minimize(simulations(w))
b.solve(disp=False)

#simulations(wt)

If you comment out the gekko section and call the simulation function at the bottom, it works fine

Answer 1

In this case, you would want to consider a different optimizer such as scipy.minimize.optimize. The function np.std() is not currently supported in Gekko. Gekko compiles the model into byte-code for automatic differentiation so you need to fit the problem into a form that is supported. Gekko's approach has several advantages, especially for large-scale or non-linear problems. For small problems with fewer than 100 variables and nearly linear constraints, an optimizer such as scipy.minimize.optimize is often a viable option. Here is your problem with a solution:

import numpy as np
from scipy.optimize import minimize

wt = np.array([0.4, 0.3, 0.3])

cov = np.array([[0.00026566, 0.00016167, 0.00011949],
   [0.00016167, 0.00065866, 0.00021662],
   [0.00011949, 0.00021662, 0.00043748]])

ret =[.098, 0.0620,.0720]
iterations = 10000
return_sim = np.random.multivariate_normal(ret, cov, iterations)

def simulations(wt):
    downside =[]
    fund_ret =np.zeros((1000,10))
    prt_ret = np.dot(return_sim , wt)
    re_ret = np.array(prt_ret).reshape(1000, 10) #10 years

    for m in range(len(re_ret)):
        fund_ret[m][0] = 100 * (1 + re_ret[m][0]) #start with $100
        for n in range(9):
            fund_ret[m][n+1] = fund_ret[m][n]* (1+re_ret[m][n+1])

    #just need the last column and all rows
    mean = np.mean(fund_ret[:,-1])  

    for i in range(1000): 
        downside.append(np.maximum((mean - fund_ret[i,-1]), 0))

    return np.std(downside)

b = (1e-5,1); bnds=(b,b,b)
cons = ({'type': 'ineq', 'fun': lambda x:  sum(x)-1},\
        {'type': 'eq',   'fun': lambda x:  np.dot(x,ret)-.07})
sol = minimize(simulations,wt,bounds=bnds,constraints=cons)
w = sol.x
print(w)

This produces the solution sol with optimal values w=sol.x:

     fun: 6.139162309118155
     jac: array([ 8.02691203, 10.04863131,  9.49171901])
 message: 'Optimization terminated successfully.'
    nfev: 33
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([0.09741111, 0.45326888, 0.44932001])