1

我正在尝试做一个优化问题，该问题需要计算受实现中变量影响的新协方差矩阵。

我可以通过 scipy 优化在我的目标函数中使用 numpy.cov 最小化来做到这一点。但是，由于我需要整数约束，我无法想出一个解决方案来解决我的 cvxpy、gekko 问题，因为大多数在线优化问题都有一个固定的协方差矩阵。

以下是我的 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})

如果有人有解决方案将不胜感激！谢谢你。

4

1 回答 1

0

xi和的协方差yi是用 显式计算的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))

该函数np.cov(xi,yi)返回一个 2x2 对称矩阵

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

[[0.80333333 2.26666667]
 [2.26666667 6.66666667]]

Gekko 需要基于梯度的优化器的协方差公式的符号形式。下面是一个cov()使用 Gekko 变量创建符号协方差计算的函数。

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 和 Numpy 对固定值xi和yi值产生相同的结果：

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

既然cov()函数已经验证过了，就可以切换x和y被计算整数值了，比如：

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

要获得整数解，请在命令m.options.SOLVER=1前切换到 (APOT) 求解器。m.solve()

于 2022-02-10T14:37:13.350 回答