二次规划&二次约束规划建模求解_二次约束二次规划

二次规划&二次约束

• 二次规划（Quadratic Programming, QP）:是运筹学中的一种优化问题，它涉及求解一个目标函数为二次函数的最优化问题，同时受到线性等式或不等式的约束。目标二次，约束一次。
• 二次约束规划：目标函数是凸二次函数，约束条件是凸二次约束和线性约束。（Quadratically Constrained Programming, QCP）

COPT４.０新增凸 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 bounds​x1​+2x2​+3x3​−0.5(33x12​+22x22​+11x32​−12x1​x2​−23x2​x3​)−x1​+x2​+x3​≤20x1​−3x2​+x3​≤300≤x1​400≤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);
    }
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

使用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 bounds​2x12​−4x1​x2​+4x22​−6x1​−3x2​x1​+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 的值
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

使用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 to​x+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.')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

参考：
https://guide.coap.online/copt/zh-doc/modeling.html#qp