cvxopt求解二次型规划

最新推荐文章于 2023-12-31 23:57:54 发布

原创最新推荐文章于 2023-12-31 23:57:54 发布 · 4.3k 阅读

22 ·

CC 4.0 BY-SA版权

本文深入解析二次型规划的数学模型及应用实例，介绍如何利用Python的cvxopt库解决二次规划问题，涵盖矩阵表达和求解过程。

QP definition

参考文献见：https://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf

二次型规划可划归为以下模型

$\begin{aligned} min_{x} \quad \frac{1}{2}x^TPx+q^Tx \\ subject \quad to \quad Gx \le h \\ \quad \quad \quad Ax=b \\ \end{aligned}$

一个例子

原问题
$\begin{aligned} min_{x,y} \quad \frac{1}{2}x^2+3x+4y \\ subject \quad to\quad x,y\ge0 \\ \quad \quad \quad x+3y\ge15 \\ \quad \quad \quad 2x+5y\le100 \\ \quad \quad \quad 3x+4y\le80 \\ \end{aligned}$
矩阵表达,公式

$\begin{aligned}min_{x,y} \quad \frac{1}{2} \begin{bmatrix}x \\y \end{bmatrix}^T \begin{bmatrix}1 &0 \\0&0 \end{bmatrix}\begin{bmatrix}x \\y \end{bmatrix}+ \begin{bmatrix}3 \\4 \end{bmatrix}\begin{bmatrix}x \\y \end{bmatrix} \\\begin{bmatrix}-1 & 0 \\ 0 & -1 \\ -1 & -3 \\2 &5 \\3 & 4\end{bmatrix}\begin{bmatrix}x \\y \end{bmatrix} \le\begin{bmatrix} 0 \\ 0 \\ -15 \\ 100 \\80\end{bmatrix} \end{aligned}$

求解

import numpy as np
import cvxopt
from cvxopt import matrix,solvers

使用numpy 和 matrix结合

P = matrix(np.diag([1,0]), tc='d')
q = matrix(np.array([3,4]), tc='d')
G = matrix(np.array([[-1,0],[0,-1],[-1,-3],[2,5],[3,4]]),tc='d')
h = matrix(np.array([0,0,-15,100,80]), tc='d')
sol = solvers.qp(P,q,G,h)
#sol = solvers.qp(P,q,G,h,A,b)
print(sol['x'])
print(sol['primal objective'])

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

20.00000617311241

使用 matrix直接求解

这里matrix的输入是转置的，比如G

Glist = [[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.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]])
print(Glist)
print(G)

[[-1.0, 0.0, -1.0, 2.0, 3.0], [0.0, -1.0, -3.0, 5.0, 4.0]]
[-1.00e+00  0.00e+00]
[ 0.00e+00 -1.00e+00]
[-1.00e+00 -3.00e+00]
[ 2.00e+00  5.00e+00]
[ 3.00e+00  4.00e+00]

P = matrix(np.diag([1,0]), tc='d')
q = matrix(np.array([3,4]), tc='d')
G = matrix(np.array([[-1,0],[0,-1],[-1,-3],[2,5],[3,4]]),tc='d')
h = matrix(np.array([0,0,-15,100,80]), tc='d')
sol = solvers.qp(P,q,G,h)
#sol = solvers.qp(P,q,G,h,A,b)
print(sol['x'])
print(sol['primal objective'])

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

20.00000617311241