Optimization integer programming with covariance matrix

Question

I am trying to do a optimization problem which requires the calculation of a new covariance matrix affected by the variable within the implementation.

I am able to do so with scipy optimization Minimize using numpy.cov within my objective function. However, as I need to have integer constraints, I am not able to think of a solution which tackles my issue with cvxpy, gekko since most of the optimization problem online have a fixed covariance matrix.

Below is my code for scipy:

room_revpar = np.array(df.iloc[:,1:10])
nla = np.array([753.2,1077.6, 1278.6,1463.9,1657.0,1990.6,2404.9,2754.6,3464.72])
min_nla = 270517.16
max_nla = 271270.359995

def objective(x, room_revpar,nla,sign = -1.0):
    room_revenue = room_revpar * x
    avg_revenue = np.mean(room_revenue, axis = 0)
    
    total_revenue = sum(avg_revenue)
    
    cov_matrix = np.cov(room_revenue.T)
    
    total_nla = np.matmul(x.T, nla)
    weights = x * nla / total_nla
    
    portfolio_sd = np.sqrt(np.matmul(np.matmul(weights.T, cov_matrix), weights))
    
    adj_risk = total_revenue / portfolio_sd
    return sign * adj_risk

def constraint1(x, nla, min_nla):
    total_nla = np.matmul(x.T, nla)
    return total_nla - min_nla

def constraint2(x, nla, max_nla):
    total_nla = np.matmul(x.T, nla)
    return max_nla - total_nla 

con1 = {'type': 'ineq', 'fun': constraint1, 'args': (nla, min_nla)}
con2 = {'type': 'ineq', 'fun': constraint2, 'args': (nla, max_nla)}

from scipy.optimize import minimize
x = np.ones(9)
sol = minimize(objective,x0 = x, args = (room_revpar, nla), constraints = (con1,con2), options = {'maxiter': 100000})

Would appreciate if anybody has a solution! Thank you.

Answer 1

The covariance of xi and yi is calculated explicitly with np.cov().

import numpy as np
xi = [2.1,2.5,3.6,4.0]
yi = [8,10,12,14]
print(np.cov(xi,yi))

The function np.cov(xi,yi) returns a 2x2 symmetric matrix

[[cov[xi,xi],cov[xi,yi]],
 [cov[xi,yi],cov[yi,yi]]]

[[0.80333333 2.26666667]
 [2.26666667 6.66666667]]

Gekko needs a symbolic form of the covariance formula for the gradient-based optimizer. Below is a function cov() that creates the symbolic covariance calculation with Gekko variables.

import numpy as np
from gekko import GEKKO

def cov(m,x,y,ddof=1):
    ''' Calculate the covariance matrix of x, y
    Inputs:
      m: Gekko model
      x: x vector of equal length to y
      y: y vector of equal length to x
      [ddof=1]: delta degrees of freedom
    Returns:
      c: covariance as a Gekko variable
    '''
    nx = len(x); ny = len(y)  # length of x, y
    if nx!=ny:
        print('Error: mismatch of x and y')
    xm = m.sum(x)/nx  # mean of x
    ym = m.sum(y)/ny  # mean of y
    c = m.Var()       # covariance
    m.Equation(c==(m.sum([(x[i]-xm)*(y[i]-ym) \
                     for i in range(nx)]))/(nx-ddof))
    return c

m = GEKKO()

n = 4
x = m.Array(m.Param,n)
y = m.Array(m.Param,n)
xi = [2.1,2.5,3.6,4.0]
yi = [8,10,12,14]
for i in range(n):
    x[i].value = xi[i]
    y[i].value = yi[i]

c0 = cov(m,x,y,ddof=0)
c1 = cov(m,x,y)

m.solve(disp=False)

print('Covariance (Numpy) population cov: ', np.cov(xi,yi,ddof=0)[0,1])
print('Covariance (Numpy) sample cov: ', np.cov(xi,yi)[0,1])
print('Covariance (Gekko) population cov: ', c0.value[0])
print('Covariance (Gekko) sample cov: ', c1.value[0])

Gekko and Numpy produce the same results for the fixed xi and yi values:

Covariance (Numpy) population cov:  1.7
Covariance (Numpy) sample cov:  2.2666666666666666
Covariance (Gekko) population cov:  1.7
Covariance (Gekko) sample cov:  2.2666666667

Now that the cov() function is verified, you can switch x and y to be calculated integer values such as:

x = m.Array(m.Var,n,lb=0,ub=10,integer=True)
y = m.Array(m.Var,n,lb=0,ub=5,integer=True)

To obtain an integer solution, switch to m.options.SOLVER=1 (APOPT) solver before the m.solve() command.