二次规划&二次约束规划建模求解
二次规划&二次约束
- 二次规划(Quadratic Programming, QP):是运筹学中的一种优化问题,它涉及求解一个目标函数为二次函数的最优化问题,同时受到线性等式或不等式的约束。目标二次,约束一次。
- 二次约束规划:目标函数是凸二次函数,约束条件是凸二次约束和线性约束。(Quadratically Constrained Programming, QCP)
COPT4.0新增凸 QP 和凸 QCP 问题求解能力。其中二次规划问题 QP 是解决在目标函数内部有如 x 平方以及 xy 等二次项的这样的问题。二次规划问题最早在金融领域提出,用来做投资组合优化。二次约束规划问题 QCP 则是在问题约束之内有二次项。若一个规划问题,同时有二次约束以及二次目标函数,则称之为二次约束二次规划问题(英文简称 QCQP)。
支持求解二次规划问题的求解器
支持:
- Cplex
- Gurobi
- 杉树COPT
- OR-Tools:cp_model
不支持:
- XPress:支持线性规划、整数规划
- Lingo:线性
- LP_sovle:是一个混合整数线性规划(MILP)求解器
求解二次规划/二次约束规
使用Java调用Cplex求解二次规划问题
min x 1 + 2 x 2 + 3 x 3 − 0.5 ( 33 x 1 2 + 22 x 2 2 + 11 x 3 2 − 12 x 1 x 2 − 23 x 2 x 3 ) subject to − x 1 + x 2 + x 3 ≤ 20 x 1 − 3 x 2 + x 3 ≤ 30 with bounds 0 ≤ x 1 40 0 ≤ x 2 ≤ ∞ 0 ≤ x 3 ≤ ∞ \begin{align} \min \quad & x_1+2x_2+3x_3-0.5(33x_1^2+22x_2^2+11x_3^2-12x_1x_2-23x_2x_3)\\ \text {subject to} \quad & -x_{1}+x_{2}+x_{3} \leq 20 \\ & x_{1}-3 x_{2}+x_{3} \leq 30 \\ \text {with bounds} \quad & 0 \leq x_1 40 \\ & 0 \leq x_2 \leq \infty\\ & 0 \leq x_3 \leq \infty\\ \end{align} minsubject towith boundsx1+2x2+3x3−0.5(33x12+22x22+11x32−12x1x2−23x2x3)−x1+x2+x3≤20x1−3x2+x3≤300≤x1400≤x2≤∞0≤x3≤∞
package org.example;
import ilog.concert.*;
import ilog.cplex.*;
public class CplexQuadraticExample {
public static void main(String[] args) {
try {
IloCplex cplex = new IloCplex();
cplex.setOut(null);
IloLPMatrix lp = populateByRow(cplex);
if (cplex.solve()) {
double[] x = cplex.getValues(lp);
double[] dj = cplex.getReducedCosts(lp);
double[] pi = cplex.getDuals(lp);
double[] slack = cplex.getSlacks(lp);
System.out.println("Solution status = " + cplex.getStatus());
System.out.println("Solution value = " + cplex.getObjValue());
int nvars = x.length;
for (int j = 0; j < nvars; ++j) {
System.out.println("Variable " + j +
": Value = " + x[j] +
" Reduced cost = " + dj[j]);
}
int ncons = slack.length;
for (int i = 0; i < ncons; ++i) {
System.out.println("Constraint " + i +
": Slack = " + slack[i] +
" Pi = " + pi[i]);
}
cplex.exportModel("qpex1.lp");
}
cplex.end();
} catch (IloException e) {
System.err.println("Concert exception '" + e + "' caught");
}
}
static IloLPMatrix populateByRow(IloMPModeler model) throws IloException {
IloLPMatrix lp = model.addLPMatrix();
// 决策变量取值范围
double[] lb = {0.0, 0.0, 0.0};
double[] ub = {40.0, Double.MAX_VALUE, Double.MAX_VALUE};
// 声明决策变量
IloNumVar[] x = model.numVarArray(model.columnArray(lp, 3), lb, ub);
// - x0 + x1 + x2 <= 20
// x0 - 3*x1 + x2 <= 30
double[] lhs = {-Double.MAX_VALUE, -Double.MAX_VALUE};
double[] rhs = {20.0, 30.0};
double[][] val = {
{-1.0, 1.0, 1.0},
{1.0, -3.0, 1.0}
};
int[][] ind = {
{0, 1, 2},
{0, 1, 2}
};
lp.addRows(lhs, rhs, ind, val);
// Q = 0.5 ( 33*x0*x0 + 22*x1*x1 + 11*x2*x2 - 12*x0*x1 - 23*x1*x2 )
IloNumExpr x00 = model.prod(33.0, x[0], x[0]);
IloNumExpr x11 = model.prod(22.0, x[1], x[1]);
IloNumExpr x22 = model.prod(11.0, x[2], x[2]);
IloNumExpr x01 = model.prod(-12.0, x[0], x[1]);
IloNumExpr x12 = model.prod(-23.0, x[1], x[2]);
IloNumExpr Q = model.prod(0.5, model.sum(x00, x11, x22, x01, x12));
// maximize x0 + 2*x1 + 3*x2 + Q
double[] objvals = {1.0, 2.0, 3.0};
model.add(model.maximize(model.diff(model.scalProd(x, objvals), Q)));
return (lp);
}
}
使用Python调用Gurobi求解二次规划问题
min 2 x 1 2 − 4 x 1 x 2 + 4 x 2 2 − 6 x 1 − 3 x 2 subject to x 1 + x 2 ≤ 3 4 x 1 + x 2 ≤ 9 with bounds x 1 , x 2 ≥ 0 \begin{align} \min \quad & 2x_1^2-4x_1x_2+4x_2^2-6x_1-3x_2\\ \text {subject to} \quad & x_{1} + x_{2} \leq 3 \\ & 4x_{1} + x_{2} \leq 9 \\ \text {with bounds} \quad & x_1,x_2 \geq 0 \\ \end{align} minsubject towith bounds2x12−4x1x2+4x22−6x1−3x2x1+x2≤34x1+x2≤9x1,x2≥0
from gurobipy import *
# 创建模型
M_QCP = Model("QCP")
# 变量声明
x1 = M_QCP.addVar(lb=0, ub=1e4, name="x1")
x2 = M_QCP.addVar(lb=0, ub=1e4, name="x2")
# 设置目标函数
M_QCP.setObjective(2 * x1 ** 2 - 4 * x1 * x2 + 4 * x2 ** 2 - 6 * x1 - 3 * x2, GRB.MINIMIZE)
# 添加约束
M_QCP.addConstr(x1 + x2 <= 3, "Con1")
M_QCP.addConstr(4 * x1 + x2 <= 9, "Con2")
M_QCP.Params.NonConvex = 2
# Optimize model
M_QCP.optimize()
M_QCP.write("QCP.lp")
print(' =========最优解======== ')
print('Obj is :', M_QCP.ObjVal) # 输出目标值
print('x1 is :', x1.x) # 输出 x1 的值
print('x2 is :', x2.x) # 输出 x2 的值
使用Python调用OR-Tools求解二次约束规划问题
max x + y subject to x 2 − y 2 ≤ 100 \begin{align} \max \quad & x+y\\ \text {subject to} \quad & x^2 - y^2 \leq 100 \\ \end{align} maxsubject tox+yx2−y2≤100
from ortools.sat.python import cp_model
model = cp_model.CpModel()
solver = cp_model.CpSolver()
# x,y,z为0到100之间的整数
x = model.NewIntVar(-100, 100, 'x')
y = model.NewIntVar(-100, 100, 'y')
xx = model.NewIntVar(-10000, 10000, 'mm')
yy = model.NewIntVar(-10000, 10000, 'mm')
model.AddMultiplicationEquality(xx, x, x)
model.Add(xx - yy <= 100)
# 目标函数
model.Maximize(x + y)
status = solver.Solve(model)
if status == cp_model.OPTIMAL:
print('x =', solver.Value(x))
print('y =', solver.Value(y))
print('xx =', solver.Value(xx))
print('yy =', solver.Value(yy))
print('obj =', solver.ObjectiveValue())
else:
print('No solution found.')
参考:
https://guide.coap.online/copt/zh-doc/modeling.html#qp