Quadratic Program (QP)

使用Gurobi optimizer求解二次规划(QP)问题…

$$ minimize:x + y + x^2 + xy + y^2 + yz + z^2 $$
s.t.
$$ x + 2y + 3z \geq 4 $$
$$ x + y \geq 1 $$
$$ x, y, z \geq 0 $$

1 2	import numpy as np from gurobipy import *

VAR_NUM = 3  # 变量数目
A = np.array([[1, 2, 3],[1, 1, 0]])              # 约束系数
b = np.array([4, 1])                             # 资源系数
c = np.array([1, 1, 0])                          # 目标函数系数
Q = np.array([[1, 1, 0], [0, 1, 1], [0, 0, 1]])  #

1
2
3

lb = [0,0,0]   # 变量取值下限
ub = [GRB.INFINITY, GRB.INFINITY, GRB.INFINITY]           # 变量取值上限
vtype = [GRB.CONTINUOUS, GRB.CONTINUOUS, GRB.CONTINUOUS]  # 变量类型

1	sense = [GRB.GREATER_EQUAL, GRB.GREATER_EQUAL] # 约束条件类型

实例化模型

1	model = Model()

Academic license - for non-commercial use only

创建决策变量

# 方式一：for + addVar()
vars = []
for i in range(VAR_NUM):
    # 添加变量
    x = model.addVar(lb=lb[i], ub=ub[i], vtype=vtype[i])  
    vars.append(x)

设置目标函数

# 木匾函数表达式
obj_func = QuadExpr()
for i in range(Q.shape[0]):
    for j in range(Q.shape[1]):
        if Q[i,j] != 0:
            obj_func += Q[i,j]*vars[i]*vars[j]
for k in range(c.shape[0]):
    if c[k] != 0:
        obj_func += c[k]*vars[k]
# 设置目标函数
model.setObjective(obj_func)

加入约束条件

# 方式一：for + addConstr()
rows,cols = A.shape
for i in range(rows):
    expr = LinExpr()  
    for j in range(cols):
        expr += A[i,j]*vars[j]
    # 添加约束条件
    model.addConstr(expr, sense[i], b[i])

模型优化

1	model.optimize()

Optimize a model with 2 rows, 3 columns and 5 nonzeros
Model has 5 quadratic objective terms
Coefficient statistics:
  Matrix range     [1e+00, 3e+00]
  Objective range  [1e+00, 1e+00]
  QObjective range [2e+00, 2e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 4e+00]
Presolve time: 0.02s
Presolved: 2 rows, 3 columns, 5 nonzeros
Presolved model has 5 quadratic objective terms
Ordering time: 0.00s

Barrier statistics:
 Free vars  : 2
 AA' NZ     : 6.000e+00
 Factor NZ  : 1.000e+01
 Factor Ops : 3.000e+01 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   1.68862705e+05 -1.66862803e+05  1.50e+03 4.63e-07  9.99e+05     0s
   1   3.32226526e+05 -3.31124005e+05  1.50e-03 4.55e-13  1.33e+05     0s
   2   4.88052576e+04 -4.83822339e+04  1.50e-09 2.84e-14  1.94e+04     0s
   3   7.20061857e+03 -7.03783959e+03  2.84e-14 0.00e+00  2.85e+03     0s
   4   1.07437613e+03 -1.01125049e+03  1.78e-14 1.78e-15  4.17e+02     0s
   5   1.65062979e+02 -1.40052218e+02  1.07e-14 3.55e-15  6.10e+01     0s
   6   2.73860883e+01 -1.68291419e+01  8.88e-16 2.66e-15  8.84e+00     0s
   7   5.67667771e+00 -1.96124352e-01  8.88e-16 2.22e-16  1.17e+00     0s
   8   2.86519627e+00  2.67284732e+00  1.11e-16 3.33e-16  3.85e-02     0s
   9   2.85719105e+00  2.85649015e+00  0.00e+00 0.00e+00  1.40e-04     0s
  10   2.85714291e+00  2.85714220e+00  6.66e-16 0.00e+00  1.40e-07     0s
  11   2.85714286e+00  2.85714286e+00  1.11e-16 1.46e-16  1.40e-10     0s

Barrier solved model in 11 iterations and 0.07 seconds
Optimal objective 2.85714286e+00

1 2	x = model.getAttr('x', vars) x

[0.5714285707563889, 0.42857142925598307, 0.8571428569222507]

1 2	R = model.getConstrs() R

[<gurobi.Constr R0>, <gurobi.Constr R1>]

1 2	func = model.getObjective() func

<gurobi.QuadExpr: C0 + C1 + [ C0 ^ 2 + C0 * C1 + C1 ^ 2 + C1 * C2 + C2 ^ 2 ]>

1	model.getVars()

[<gurobi.Var C0 (value 0.5714285707563889)>,
 <gurobi.Var C1 (value 0.42857142925598307)>,
 <gurobi.Var C2 (value 0.8571428569222507)>]