基于投资组合问题的凸二次规划模型及求解——Gurobi求解器+高阶牛顿法（python）_gurobi求解大规模问题

本学期上了最优化课程，学习到了一些凸优化的知识，想要尝试将这些优化理论应用到实际中。

这里基于均值-方差分析方法，针对以确定收益水平下的最小化风险为目标的投资组合问题建立了具有等式约束的凸二次规划模型，并分别利用GUROBI求解器和一种高阶牛顿算法对其进行求解。

模型描述如下（公式太多的地方就直接贴图片了）：

分别利用Gurobi求解器和高阶牛顿法求解：

Gurobi：Gurobi是美国Gurobi Optimization公司开发的大规模优化求解器，目前已经成为综合能力全球领先的数学规划求解器，作为一个全局优化器，Gurobi支持的模型类型包括：连续和混合整数线性问题、凸目标或约束连续和混合整数二次问题、含有指对数、三角函数、高阶多项式目标和约束以及任何形式的分段约束的非线性问题等。另外，Gurobi提供了方便轻巧的API接口，支持C++、Java、Python等多种语言，同时支持包括Windows、Linux以及OS在内的多种平台。

基于绝对值光华逼近的高阶牛顿法（公式真的多我是真的懒，依旧贴图片）：

代码如下：

from gurobipy import GRB, GurobiError
import gurobipy as gp
from numpy import *
import math
import pandas as pd
 
c=array([[0],[0],[0]])
Q=array([[0.2,0.3,-0.01],[0.3,2.4,0.5],[-0.01,0.5,1.1]])
A=array([[0.06,0.12,0.09],[1,1,1]])
b=array([[0.06],[1]])
c=array([[0],[0],[0]])
x=array([[0],[0],[0]])
y=array([[0],[0]])
k=1000
e=0.000001
 
def guroby_solver(Q,A,b,c):
    try:
 
        m = gp.Model("QP_guroby")
        var = []
        constr = 0
        obj = 0
 
        for i in range(len(A[0])):
            var.append(m.addVar(lb=0, vtype=GRB.CONTINUOUS, name='var[%d]' % i))  #定义变量
 
        for i in range(len(A)):
            constr = 0
            for j in range(len(A[0])):
                if A[i][j] != 0:
                    constr = constr + A[i][j] * var[j]
            m.addConstr(constr == b[i][0], name='constraint %d' % i)     #添加约束
 
        for i in range(len(A[0])):
            for j in range(len(A[0])):
                obj = obj + Q[i][j] * var[i] * var[j]
        for i in range(len(A[0])):
            obj = obj + c[i][0] * var[i]
        m.setObjective(obj, GRB.MINIMIZE)    #添加目标函数
 
        m.optimize()
        variable=[]
        for v in m.getVars():
            print(v.varName, v.x)
            variable.append(v.x)
        print('Obj:', m.objVal)
 
    except GurobiError:
        print('Error reported')
    return variable,m.objVal          #GUROBI求解器求解函数
 
def newton(x,y,e,k,Q,A,b,c):
 
    def cal_pk(m):
        x_pk=zeros((len(m),1))
        for i in range(len(m)):
            x_pk[i][0]=math.sqrt(m[i][0]**2+1/(k**2))
        return x_pk
 
    def cal_Fk(X,m):
        x_split=X[:m,:]
        y_split=X[m:,:]
        Fk_up=dot(A,cal_pk(x_split)-x_split)-b
        Fk_down=cal_pk(x_split)+x_split+ \
                dot(A.T,y_split)+dot(Q,x_split)-dot(Q,cal_pk(x_split))-c
        Fk=vstack((Fk_up,Fk_down))
        return Fk
 
    def cal_Fk_grad(X,m,n):
        x_split=X[:m,:]
        Dk=zeros((m,m))
        for i in range(m):
            Dk[i][i]=x_split[i][0]/math.sqrt(x_split[i][0]**2+1/(k**2))
        O=zeros((n,n))
        Fk_grad_up=hstack((dot(A,Dk-identity(3)),O))
        Fk_grad_down=hstack((Dk+identity(3)+dot(Q,identity(3)-Dk),A.T))
        Fk_grad=vstack((Fk_grad_up,Fk_grad_down))
        return Fk_grad
 
    X=vstack((x,y))
    while True:
        if dot(cal_Fk(X,len(x)).T,cal_Fk(X,len(x)))<e:
            break
        U=X-0.5*dot(linalg.inv(cal_Fk_grad(X,len(x),len(y))),cal_Fk(X,len(x)))
        V=X-dot(linalg.inv(cal_Fk_grad(U,len(x),len(y))),cal_Fk(X,len(x)))
        X=V+dot(linalg.inv(cal_Fk_grad(X,len(x),len(y))) \
                -2*linalg.inv(cal_Fk_grad(U,len(x),len(y))),cal_Fk(V,len(x)))
    re_x=X[:len(x),:]
    z=zeros((len(x),1))
    for i in range(len(x)):
        z[i][0]=abs(re_x[i][0])-re_x[i][0]
    mid=dot(Q,z)
    obj_fuc=dot(z.T,mid)+dot(c.T,z)
    print('高阶牛顿法的最优解为：',z)
    print('目标函数值为：',obj_fuc)
 
    return z,obj_fuc              #牛顿法求解函数
 
df=pd.DataFrame(index=range(13),columns=['rp','z1_g','z2_g', \
                                         'z3_g','obj_g','z1_n', \
                                         'z2_n','z3_n','obj_n'])
for i in range(13):
    b[0][0]=0.06+i*0.005
    m,n=guroby_solver(Q,A,b,c)
    p,q=newton(x,y,e,k,Q,A,b,c)
    df['rp'][i]=b[0][0]
    df['z1_g'][i]=m[0]
    df['z2_g'][i]=m[1]
    df['z3_g'][i]=m[2]
    df['obj_g'][i]=n
    df['z1_n'][i]=p[0][0]
    df['z2_n'][i]=p[1][0]
    df['z3_n'][i]=p[2][0]
    df['obj_n'][i]=q[0][0]
df.to_csv('Result.csv')     #将计算结果导入excel

参考文献：

雍龙泉,贾伟,黎延海.基于光滑逼近函数的高阶牛顿法求解凸二次规划[J].科学技术与工程,2021,21(06):2151-2156.