使用Gurobi optimizer求解二次规划(QP)问题…
m
i
n
i
m
i
z
e
:
x
+
y
+
x
2
+
x
y
+
y
2
+
y
z
+
z
2
minimize:x + y + x^2 + xy + y^2 + yz + z^2
minimize:x+y+x2+xy+y2+yz+z2
s.t.
x
+
2
y
+
3
z
≥
4
x + 2y + 3z \geq 4
x+2y+3z≥4
x
+
y
≥
1
x + y \geq 1
x+y≥1
x
,
y
,
z
≥
0
x, y, z \geq 0
x,y,z≥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)>]