Quadratic Program (QP)

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

---

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

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]])  #

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

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

实例化模型

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])

模型优化

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

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

[0.5714285707563889, 0.42857142925598307, 0.8571428569222507]

R = model.getConstrs()
R

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

func = model.getObjective()
func

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

model.getVars()

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