zoukankan      html  css  js  c++  java
  • 在Python中利用CVXOPT求解二次规划问题

    工作中需要用到cvxopt,cvxopt安装有坑,大家注意下.
    1.首先一定要卸载numpy,无论是直接安装的,还是anaconda安装的,主要是必须用whl安装numpy才不会有包的冲突
    2.二次规划包的使用
    二次规划的标准形式如下

    Minquad frac{1}{2}x^{t}Px+q^{t}x

    s.t. quad Gxleq h

    Ax=b

    Python 代码如下

    from cvxopt import matrix
    import cvxopt.solvers as sol
    result = sol.qp(P, Q, G, h, A, b)

    问题描述:
        在实际生活中,我们经常会遇到一些优化问题,简单的线性规划可以作图求解,但是对于目标函数包含二次项时,则需要另觅它法在金融实践中,马科维茨均方差模型就有实际的二次优化需求

    作为金融实践中常用的方法,本篇将对CVXOPT中求解二次规划的问题进行举例详细说明,关于该方法在均方差优化中的实践应用,参见后续发帖

    1、二次规划问题的标准形式

    min12xTPx+qTx
    s.t.Gx≤h
    Ax=b

    上式中,x为所要求解的列向量,xT表示x的转置
    接下来,按步骤对上式进行相关说明:
        上式表明,任何二次规划问题都可以转化为上式的结构,事实上用cvxopt的第一步就是将实际的二次规划问题转换为上式的结构,写出对应的P、q、G、h、A、b目标函数若为求max,可以通过乘以−1,将最大化问题转换为最小化问题Gx≤b表示的是所有的不等式约束,同样,若存在诸如x≥0的限制条件,也可以通过乘以−1转换为"≤"的形式Ax=b表示所有的等式约束

    2、以一个标准的例子进行过程说明

    min(x,y)12x2+3x+4y
    s.t.x,y≥0
    x+3y≥15
    2x+5y≤100
    3x+4y≤80

    例子中,需要求解的是x,y,我们可以把它写成向量的形式,同时,也需要将限制条件按照上述标准形式进行调整,用矩阵形式表示,如下所示:

    min(x,y)12[xy]T[100][xy]+[34]T[xy]
    [−10 0−1-1−3 25 34][xy]≤[0-1510080]

        如上所示,目标函数和限制条件均转化成了二次规划的标准形式,这是第一步,也是最难的一步,接下来的事情就简单了对比上式和标准形式,不难得出:P=[100],q=[34],G=[−10 0−1-1−3 25 34],h=[0-1510080]
    接下来就是几行简单的代码,目的是告诉计算机上面的参数具体是什么

    from cvxopt import solvers, matrix
    P = matrix([[1.0,0.0],[0.0,0.0]]) # matrix里区分int和double,所以数字后面都需要加小数点
    q = matrix([3.0,4.0])
    G = matrix([[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.0,4.0]])
    h = matrix([0.0,0.0,-15.0,100.0,80.0])
    sol = solvers.qp(P,q,G,h) # 调用优化函数solvers.qp求解
    print sol['x'] # 打印结果,sol里面还有很多其他属性,读者可以自行了解


    pcost dcost gap pres dres

    0: 1.0780e+02 -7.6366e+02 9e+02 1e-16 4e+01
    1: 9.3245e+01 9.7637e+00 8e+01 1e-16 3e+00
    2: 6.7311e+01 3.2553e+01 3e+01 6e-17 1e+00
    3: 2.6071e+01 1.5068e+01 1e+01 2e-16 7e-01
    4: 3.7092e+01 2.3152e+01 1e+01 2e-16 4e-01
    5: 2.5352e+01 1.8652e+01 7e+00 8e-17 3e-16
    6: 2.0062e+01 1.9974e+01 9e-02 6e-17 3e-16
    7: 2.0001e+01 2.0000e+01 9e-04 6e-17 3e-16
    8: 2.0000e+01 2.0000e+01 9e-06 9e-17 2e-16
    Optimal solution found.
    [ 7.13e-07]
    [ 5.00e+00]


        看了上面的代码,是不是觉得很简单。因为难点不在代码,而是在于将实际优化问题转化为标准形式的过程在上面的例子中,并没有出现等号,当出现等式约束时,过程一样,找到A,b,然后运行代码 sol = solvers.qp(P,q,G,h,A,b) 即可求解
    扩展:上述定义各个矩阵参数用的是最直接的方式,其实也可以结合Numpy来定义上述矩阵

    from cvxopt import solvers, matrix
    import numpy as np
    P = matrix(np.diag([1.0,0])) # 对于一些特殊矩阵,用numpy创建会方便很多(在本例中可能感受不大)
    q = matrix(np.array([3.0,4]))
    G = matrix(np.array([[-1.0,0],[0,-1],[-1,-3],[2,5],[3,4]]))
    h = matrix(np.array([0.0,0,-15,100,80]))
    sol = solvers.qp(P,q,G,h)



    pcost dcost gap pres dres

    0: 1.0780e+02 -7.6366e+02 9e+02 1e-16 4e+01
    1: 9.3245e+01 9.7637e+00 8e+01 1e-16 3e+00
    2: 6.7311e+01 3.2553e+01 3e+01 6e-17 1e+00
    3: 2.6071e+01 1.5068e+01 1e+01 2e-16 7e-01
    4: 3.7092e+01 2.3152e+01 1e+01 2e-16 4e-01
    5: 2.5352e+01 1.8652e+01 7e+00 8e-17 3e-16
    6: 2.0062e+01 1.9974e+01 9e-02 6e-17 3e-16
    7: 2.0001e+01 2.0000e+01 9e-04 6e-17 3e-16
    8: 2.0000e+01 2.0000e+01 9e-06 9e-17 2e-16
    Optimal solution found.

    先写到这吧,关于二次规划在均方差优化中的实践应用,参见后续发帖,欢迎交流~~出处
    发布于 2018-05-11

  • 相关阅读:
    System.currentTimeMillis();
    Spark Core源代码分析: Spark任务模型
    CocoaPods on Xcode 6 and Yosemite
    leetcode第一刷_Binary Tree Inorder Traversal
    struts2讲义----建立一个struts2工程
    Java实现 蓝桥杯VIP 算法提高 班级排名
    Java实现 蓝桥杯VIP 算法提高 种树
    Java实现 蓝桥杯VIP 算法提高 种树
    Java实现 蓝桥杯VIP 算法提高 种树
    Java实现 蓝桥杯VIP 算法提高 种树
  • 原文地址:https://www.cnblogs.com/Rvin/p/9639308.html
Copyright © 2011-2022 走看看