二次规划(QP)求解器只依赖于NumPy/SciPy?

b1uwtaje  于 9个月前  发布在  其他
关注(0)|答案(6)|浏览(121)

我希望学生在作业中解决二次规划,而无需安装额外的软件,如cvxopt等。是否有只依赖于NumPy/SciPy的Python实现?

hec6srdp

hec6srdp1#

我对二次规划不是很熟悉,但我认为你可以用scipy.optimize的约束最小化算法来解决这类问题。下面是一个例子:

import numpy as np
from scipy import optimize
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

# minimize
#     F = x[1]^2 + 4x[2]^2 -32x[2] + 64

# subject to:
#      x[1] + x[2] <= 7
#     -x[1] + 2x[2] <= 4
#      x[1] >= 0
#      x[2] >= 0
#      x[2] <= 4

# in matrix notation:
#     F = (1/2)*x.T*H*x + c*x + c0

# subject to:
#     Ax <= b

# where:
#     H = [[2, 0],
#          [0, 8]]

#     c = [0, -32]

#     c0 = 64

#     A = [[ 1, 1],
#          [-1, 2],
#          [-1, 0],
#          [0, -1],
#          [0,  1]]

#     b = [7,4,0,0,4]

H = np.array([[2., 0.],
              [0., 8.]])

c = np.array([0, -32])

c0 = 64

A = np.array([[ 1., 1.],
              [-1., 2.],
              [-1., 0.],
              [0., -1.],
              [0.,  1.]])

b = np.array([7., 4., 0., 0., 4.])

x0 = np.random.randn(2)

def loss(x, sign=1.):
    return sign * (0.5 * np.dot(x.T, np.dot(H, x))+ np.dot(c, x) + c0)

def jac(x, sign=1.):
    return sign * (np.dot(x.T, H) + c)

cons = {'type':'ineq',
        'fun':lambda x: b - np.dot(A,x),
        'jac':lambda x: -A}

opt = {'disp':False}

def solve():

    res_cons = optimize.minimize(loss, x0, jac=jac,constraints=cons,
                                 method='SLSQP', options=opt)

    res_uncons = optimize.minimize(loss, x0, jac=jac, method='SLSQP',
                                   options=opt)

    print '\nConstrained:'
    print res_cons

    print '\nUnconstrained:'
    print res_uncons

    x1, x2 = res_cons['x']
    f = res_cons['fun']

    x1_unc, x2_unc = res_uncons['x']
    f_unc = res_uncons['fun']

    # plotting
    xgrid = np.mgrid[-2:4:0.1, 1.5:5.5:0.1]
    xvec = xgrid.reshape(2, -1).T
    F = np.vstack([loss(xi) for xi in xvec]).reshape(xgrid.shape[1:])

    ax = plt.axes(projection='3d')
    ax.hold(True)
    ax.plot_surface(xgrid[0], xgrid[1], F, rstride=1, cstride=1,
                    cmap=plt.cm.jet, shade=True, alpha=0.9, linewidth=0)
    ax.plot3D([x1], [x2], [f], 'og', mec='w', label='Constrained minimum')
    ax.plot3D([x1_unc], [x2_unc], [f_unc], 'oy', mec='w',
              label='Unconstrained minimum')
    ax.legend(fancybox=True, numpoints=1)
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    ax.set_zlabel('F')

字符串
输出量:

Constrained:
  status: 0
 success: True
    njev: 4
    nfev: 4
     fun: 7.9999999999997584
       x: array([ 2.,  3.])
 message: 'Optimization terminated successfully.'
     jac: array([ 4., -8.,  0.])
     nit: 4

Unconstrained:
  status: 0
 success: True
    njev: 3
    nfev: 5
     fun: 0.0
       x: array([ -2.66453526e-15,   4.00000000e+00])
 message: 'Optimization terminated successfully.'
     jac: array([ -5.32907052e-15,  -3.55271368e-15,   0.00000000e+00])
     nit: 3


的数据

jmo0nnb3

jmo0nnb32#

这可能是一个迟来的答案,但我发现CVXOPT-http://cvxopt.org/-作为Quadratic Programming常用的免费python库。然而,它不容易安装,因为它需要安装其他依赖项。

igetnqfo

igetnqfo3#

我遇到了一个很好的解决方案,并希望将其发布。在NICTA的ELEFANT机器学习工具包中有一个LOQO的Python实现(截至本文发布时为http://elefant.forge.nicta.com.au)。看看optimization.intpointsolver。这是由Alex Smola编写的,我已经使用了相同代码的C版本并取得了巨大的成功。

sg24os4d

sg24os4d4#

qpsolvers软件包似乎也符合要求。它只依赖于NumPy,可以通过pip install qpsolvers安装。然后,您可以执行以下操作:

from numpy import array, dot
from qpsolvers import solve_qp

M = array([[1., 2., 0.], [-8., 3., 2.], [0., 1., 1.]])
P = dot(M.T, M)  # quick way to build a symmetric matrix
q = dot(array([3., 2., 3.]), M).reshape((3,))
G = array([[1., 2., 1.], [2., 0., 1.], [-1., 2., -1.]])
h = array([3., 2., -2.]).reshape((3,))

# min. 1/2 x^T P x + q^T x with G x <= h
print "QP solution:", solve_qp(P, q, G, h)

字符串
您还可以通过更改solver关键字参数(例如solver='cvxopt'solver='osqp')来尝试不同的QP解算器(例如Curious提到的CVXOPT)。

mepcadol

mepcadol5#

mystic提供了一个纯Python实现的非线性/非凸优化算法,具有通常只在QP求解器中才能找到的高级约束功能。mystic实际上提供了比大多数QP求解器更强大的约束。但是,如果你正在寻找优化算法速度,那么下面的代码不适合你。mystic并不慢,但它是纯python,而不是python绑定到C。如果您正在寻找非线性求解器中的灵活性和QP约束功能,那么您可能会感兴趣。

"""
Maximize: f = 2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2

Subject to: -2*x[0] + 2*x[1] <= -2
             2*x[0] - 4*x[1] <= 0
               x[0]**3 -x[1] == 0

where: 0 <= x[0] <= inf
       1 <= x[1] <= inf
"""
import numpy as np
import mystic.symbolic as ms
import mystic.solvers as my
import mystic.math as mm

# generate constraints and penalty for a nonlinear system of equations 
ieqn = '''
   -2*x0 + 2*x1 <= -2
    2*x0 - 4*x1 <= 0'''
eqn = '''
     x0**3 - x1 == 0'''
cons = ms.generate_constraint(ms.generate_solvers(ms.simplify(eqn,target='x1')))
pens = ms.generate_penalty(ms.generate_conditions(ieqn), k=1e3)
bounds = [(0., None), (1., None)]

# get the objective
def objective(x, sign=1):
  x = np.asarray(x)
  return sign * (2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2)

# solve    
x0 = np.random.rand(2)
sol = my.fmin_powell(objective, x0, constraint=cons, penalty=pens, disp=True,
                     bounds=bounds, gtol=3, ftol=1e-6, full_output=True,
                     args=(-1,))

print 'x* = %s; f(x*) = %s' % (sol[0], -sol[1])

字符串
需要注意的是,mystic可以通用地将LP、QP和高阶等式和不等式约束应用于任何给定的优化器,而不仅仅是一个特殊的QP求解器。其次,mystic可以消化符号数学,因此定义/输入约束的容易程度比处理函数的矩阵和导数要好一点。mystic依赖于numpy,如果安装了scipy,则将使用scipy(但是,scipy不是必需的)。mystic利用sympy来处理符号约束,但通常优化也不需要它。
输出量:

Optimization terminated successfully.
         Current function value: -2.000000
         Iterations: 3
         Function evaluations: 103

x* = [ 2.  1.]; f(x*) = 2.0


在这里获取mystichttps://github.com/uqfoundation

f3temu5u

f3temu5u6#

可以尝试python模块quadprogQP-tasksOLS和尝试example

import numpy as np
import quadprog
from scipy import optimize
# e.g. least-squares problem
#print("    min. || G * x - a ||^2")
#print("    s.t. C * x <= b")
#print('\n')  

def solve_qp_scipy(G, a, C, b, meq):
    def f(x):
        return 0.5 * np.dot(x, G).dot(x) - np.dot(a, x)

    constraints = []
    if C is not None:
        constraints = [{
            'type': 'eq' if i < meq else 'ineq',
            'fun': lambda x, C=C, b=b, i=i: (np.dot(C.T, x) - b)[i]
        } for i in range(C.shape[1])]

    result = optimize.minimize(
        f, x0=np.zeros(len(G)), method='SLSQP', constraints=constraints,
        tol=1e-10, options={'maxiter': 2000})
    return result

# data: 
M = np.array([[1.0, 2.0, 0.0], [-8.0, 3.0, 2.0], [0.0, 1.0, 1.0]])
G = np.dot(M.T, M)  # this is a positive definite matrix
a = np.dot(np.array([3.0, 2.0, 3.0]), M)   # @ M
C = np.array([[1.0, 2.0, 1.0], [2.0, 0.0, 1.0], [-1.0, 2.0, -1.0]])
b = np.array([3.0, 2.0, -2.0])     #.reshape((3,))
meq = b.size

# for dense problem from https://github.com/qpsolvers/qpsolvers/blob/main/examples/quadratic_program.py
print("    min. 1/2 x^T P x + a^T x")
print("    s.t. G * x <= b")

# quadprog
xf, f, xu, iters, lagr, iact = quadprog.solve_qp(G, a, C, b, meq)
print(xf, '\n', f, '\n')

# SciPy
result = solve_qp_scipy(G, a, C, b, meq)
print(result.x, '\n', result.fun, '\n')

np.testing.assert_array_almost_equal(result.x, xf)
np.testing.assert_array_almost_equal(result.fun, f)

字符串
P.S.其他qpsolvers,我不能pip installmsys32\mingw32 gcc 9.3.0vs_BT2019 需要,我想
P.P.S.有价值的参考资料here-“如何不与统计数据撒谎:总结基准测试结果的正确方法”

相关问题