5.4.3. C++的SDP建模与优化¶
在本节中，我们将使用 MindOpt 的 C++ API 来建模以及求解 半定规划问题示例 中的问题。
5.4.3.1. SDP示例1¶
首先，需要引入头文件：
25#include "MindoptCpp.h"
第一步：创建优化模型：
先建立一空的MindOpt模型。
46 MDOEnv env = MDOEnv();
47 MDOModel m = MDOModel(env);
第二步：SDP模型输入
我们利用 MDOModel::addPsdVar() 创建一个新的半正定矩阵变量，并同时设定其对应的目标函数中的矩阵系数。其中，第一个参数为 矩阵系数，这里它是通过 MDOMatrix::coo() 方法创建的一个 MDOMatrix 类的实例。注意到该系数矩阵的维度应与矩阵变量相同。第二个参数为该变量的名称。
54 /* Add variables. */
55 MDOPsdVar psd_var = m.addPsdVar(MDOMatrix::coo(dim_mat, dim_mat, C_nnz, C_nz_indices, C_nz_values), "X0");
接着，我们输入第一个约束。我们利用 MDOModel::addPsdConstr() 建立一个带半正定矩阵变量的约束。其中，第一个参数分别为该约束中的 半正定表达式，即 \(\langle \mathbf{A},\mathbf{X} \rangle\)。
我们通过 MDOPsdExpr 类的初始化方法 MDOPsdExpr::MDOPsdExpr() 创建。第二个系数为该约束的属性，目前我们只支持 MDO_EQUAL (‘=’)，代表等式约束。第三个参数为该约束的右侧值 (此处为1.0)。最后一个参数为该约束的名称。
57 /* Add constraints. */
58 m.addPsdConstr(MDOPsdExpr(psd_var, MDOMatrix::coo(dim_mat, dim_mat, A_nnz, A_nz_indices, A_nz_values)), MDO_EQUAL, 1.0, "C0");
最后，我们利用 MDOModel::setObjective() 设定优化函数中的线性部分 (此处为0)，并将优化方向改为maximization (-1)。
60 /* Set objective function. */
61 MDOLinExpr objective = 0;
62 m.setObjective(objective, -1);
第三步：求解SDP模型
模型输入后，我们接着用 MDOModel::optimize() 来求解问题。
67 // Optimize the model
68 m.optimize();
第四步: 取得SDP模型的解
我们利用泛用函数 MDOModel::get() 来取得最优的的目标函数值，即 PrimalObjVal 属性。
69 if (m.get(MDO_IntAttr_Status) == MDO_OPTIMAL)
70 {
71 cout << "Primal objective val: " << m.get(MDO_DoubleAttr_PrimalObjVal) << endl;
72 }
73 else
74 {
75 cout << "No feasible solution." << endl;
76 }
文件链接 MdoSdoEx1.cpp 提供了完整源代码：
1/**
2 * Description
3 * -----------
4 *
5 * Semidefinite optimization (row-wise input).
6 *
7 * Formulation
8 * -----------
9 *
10 * Maximize
11 * obj: tr(C X)
12 *
13 * Subject To
14 * c0 : tr(A X) = 1
15 * X is p.s.d.
16 *
17 * Matrix
18 * C = [ -3 0 1 ] A = [ 3 0 1 ]
19 * [ 0 -2 0 ] [ 0 4 0 ]
20 * [ 1 0 -3 ] [ 1 0 5 ]
21 * End
22 */
23
24#include <iostream>
25#include "MindoptCpp.h"
26
27using namespace std;
28
29int main(void)
30{
31 /* Model data. */
32 int num_mats = 1;
33 int dim_mat = 3; /* Dimension of the matrix variables. */
34
35 int C_nnz = 4;
36 int C_nz_indices[] = { 0, 2, 4, 8 }; /* Nonzero vectorized index of obj coeff. */
37 double C_nz_values[] = { -3.0, 1.0, -2.0, -3.0 }; /* Nonzero values of obj coeff. */
38
39 int A_nnz = 4;
40 int A_nz_indices[] = { 0, 2, 4, 8 }; /* Nonzero vectorized index of constr coeff. */
41 double A_nz_values[] = { 3.0, 1.0, 4.0, 5.0 }; /* Nonzero values of constr coeff. */
42
43 /*------------------------------------------------------------------*/
44 /* Step 1. Create environment and model. */
45 /*------------------------------------------------------------------*/
46 MDOEnv env = MDOEnv();
47 MDOModel m = MDOModel(env);
48
49 try
50 {
51 /*------------------------------------------------------------------*/
52 /* Step 2. Input model. */
53 /*------------------------------------------------------------------*/
54 /* Add variables. */
55 MDOPsdVar psd_var = m.addPsdVar(MDOMatrix::coo(dim_mat, dim_mat, C_nnz, C_nz_indices, C_nz_values), "X0");
56
57 /* Add constraints. */
58 m.addPsdConstr(MDOPsdExpr(psd_var, MDOMatrix::coo(dim_mat, dim_mat, A_nnz, A_nz_indices, A_nz_values)), MDO_EQUAL, 1.0, "C0");
59
60 /* Set objective function. */
61 MDOLinExpr objective = 0;
62 m.setObjective(objective, -1);
63
64 /*------------------------------------------------------------------*/
65 /* Step 3. Solve the problem and populate optimization result. */
66 /*------------------------------------------------------------------*/
67 // Optimize the model
68 m.optimize();
69 if (m.get(MDO_IntAttr_Status) == MDO_OPTIMAL)
70 {
71 cout << "Primal objective val: " << m.get(MDO_DoubleAttr_PrimalObjVal) << endl;
72 }
73 else
74 {
75 cout << "No feasible solution." << endl;
76 }
77 }
78 catch (MDOException &e)
79 {
80 cout << "Error code = " << e.getErrorCode() << endl;
81 cout << e.getMessage() << endl;
82 }
83 catch (...)
84 {
85 cout << "Error during optimization." << endl;
86 }
87
88 return static_cast<int>(MDO_OKAY);
89}
5.4.3.2. SDP示例2¶
首先，需要引入头文件：
32#include "MindoptCpp.h"
第一步：创建MindOpt模型
首先，我们必须先建立一空的MindOpt模型。
64 MDOEnv env = MDOEnv();
65 MDOModel m = MDOModel(env);
第二步：SDP模型输入
接着，我们利用 MDOModel::addPsdVar() 创建两个新的半正定矩阵变量，并同时设定他们对应的目标函数中的矩阵系数。其中，第一个参数为 矩阵系数，这里它是通过 MDOMatrix::coo() 方法创建的一个 MDOMatrix 类的实例。注意到该系数矩阵的维度应与矩阵变量相同。第二个参数为该变量的名称。
73 MDOPsdVar psd_var0 = m.addPsdVar(MDOMatrix::coo(dim_mat[0], dim_mat[0], C0_nnz, C0_nz_indices, C0_nz_values), "X0");
74 MDOPsdVar psd_var1 = m.addPsdVar(MDOMatrix::coo(dim_mat[1], dim_mat[1], C1_nnz, C1_nz_indices, C1_nz_values), "X1");
我们再利用 MDOModel::addVar() 创建两个线性变量。其中，第一个和第二个参数为 变量下界与上界。第三个参数为目标函数中该变量对应的系数。第四个参数为该变量的类型，在这里是连续变量类型 (‘C’)。最后一个参数是该变量名称。
75 MDOVar var0 = m.addVar(0.0, MDO_INFINITY, 0.0, 'C', "x0");
76 MDOVar var1 = m.addVar(0.0, MDO_INFINITY, 0.0, 'C', "x1");
下一步，我们创建约束。我们利用 MDOModel::addPsdConstr() 建立带半正定矩阵变量的约束。其中，第一个参数分别为该约束中的 半正定表达式。由于在该例子中，约束同时包含半正定变量与线性变量，我们分两步构造其半正定表达式。我们首先通过 MDOPsdExpr 类的初始化方法 MDOPsdExpr::MDOPsdExpr() 创建半正定部分，再通过 MDOPsdExpr::addTerms() 添加剩余的线性部分， 得到最终的该约束对应的半正定表达式。第二个系数为该约束的属性，目前我们只支持 MDO_EQUAL (‘=’)，代表等式约束。第三个参数为该约束的右侧值 (此处为1.0)。最后一个参数为该约束的名称。
78 /* Add constraints. */
79 MDOPsdExpr psd_expr0 = MDOPsdExpr(psd_var0, MDOMatrix::coo(dim_mat[0], dim_mat[0], A0_nnz, A0_nz_indices, A0_nz_values));
80 const MDOVar row0_vars[] = { var0 };
81 psd_expr0.addTerms(row0_values, row0_vars, 1);
82 m.addPsdConstr(psd_expr0, MDO_EQUAL, 1.0, "C0");
83
84 MDOPsdExpr psd_expr1 = MDOPsdExpr(psd_var1, MDOMatrix::coo(dim_mat[0], dim_mat[0], A1_nnz, A1_nz_indices, A1_nz_values));
85 const MDOVar row1_vars[] = { var1 };
86 psd_expr1.addTerms(row1_values, row1_vars, 1);
87 m.addPsdConstr(psd_expr1, MDO_EQUAL, 2.0, "C1");
最后，我们利用 MDOModel::setObjective() 设定优化函数中的线性部分 (此处为0)，并将优化方向改为maximization (-1)。
89 /* Set objective function. */
90 MDOLinExpr objective = 0;
91 m.setObjective(objective, -1);
第三步：求解SDP模型
模型输入后，我们接着用 MDOModel::optimize() 来求解问题。
96 // Optimize the model
97 m.optimize();
第四步: 取得SDP模型的解
我们利用泛用函数 MDOModel::get() 来取得最优的的目标函数值，即 PrimalObjVal 属性。
98 if (m.get(MDO_IntAttr_Status) == MDO_OPTIMAL)
99 {
100 cout << "Primal objective val: " << m.get(MDO_DoubleAttr_PrimalObjVal) << endl;
101 }
102 else
103 {
104 cout << "No feasible solution." << endl;
105 }
文件链接 MdoSdoEx2.cpp 提供了完整源代码：
1/**
2 * Description
3 * -----------
4 *
5 * Semidefinite optimization (row-wise input).
6 *
7 * Formulation
8 * -----------
9 *
10 * Maximize
11 * obj: tr(C0 X0) + tr(C0 X1) + 0 x0 + 0 x1
12 *
13 * Subject To
14 * c0 : tr(A00 X0) + 1 x0 = 1
15 * c1 : tr(A00 X1) + 1 x1 = 2
16 * Bounds
17 * 0 <= x0
18 * 0 <= x1
19 * X1, X2 are p.s.d.
20 *
21 * Matrix
22 * C0 = [ 2 1 ] A00 = [ 3 1 ]
23 * [ 1 2 ] [ 1 3 ]
24 *
25 * C0 = [ 3 0 1 ] A00 = [ 3 0 1 ]
26 * [ 0 2 0 ] [ 0 4 0 ]
27 * [ 1 0 3 ] [ 1 0 5 ]
28 * End
29 */
30
31#include <iostream>
32#include "MindoptCpp.h"
33
34using namespace std;
35
36int main(void)
37{
38 /* Model data. */
39 int num_mats = 1;
40 int dim_mat[] = { 2, 3 }; /* Dimension of the matrix variables. */
41
42 int C0_nnz = 3;
43 int C0_nz_indices[] = { 0, 1, 3 }; /* Nonzero vectorized index of obj coeff. */
44 double C0_nz_values[] = { 2.0, 1.0, 2.0 }; /* Nonzero values of obj coeff. */
45
46 int C1_nnz = 4;
47 int C1_nz_indices[] = { 0, 2, 4, 8 }; /* Nonzero vectorized index of obj coeff. */
48 double C1_nz_values[] = { 3.0, 1.0, 2.0, 3.0 }; /* Nonzero values of obj coeff. */
49
50 int A0_nnz = 3;
51 int A0_nz_indices[] = { 0, 1, 3 }; /* Nonzero vectorized index of constr coeff. */
52 double A0_nz_values[] = { 3.0, 1.0, 3.0 }; /* Nonzero values of constr coeff. */
53
54 int A1_nnz = 4;
55 int A1_nz_indices[] = { 0, 2, 4, 8 }; /* Nonzero vectorized index of constr coeff. */
56 double A1_nz_values[] = { 3.0, 1.0, 4.0, 5.0 }; /* Nonzero values of constr coeff. */
57
58 const double row0_values[] = { 1.0 };
59 const double row1_values[] = { 1.0 };
60
61 /*------------------------------------------------------------------*/
62 /* Step 1. Create environment and model. */
63 /*------------------------------------------------------------------*/
64 MDOEnv env = MDOEnv();
65 MDOModel m = MDOModel(env);
66
67 try
68 {
69 /*------------------------------------------------------------------*/
70 /* Step 2. Input model. */
71 /*------------------------------------------------------------------*/
72 /* Add variables. */
73 MDOPsdVar psd_var0 = m.addPsdVar(MDOMatrix::coo(dim_mat[0], dim_mat[0], C0_nnz, C0_nz_indices, C0_nz_values), "X0");
74 MDOPsdVar psd_var1 = m.addPsdVar(MDOMatrix::coo(dim_mat[1], dim_mat[1], C1_nnz, C1_nz_indices, C1_nz_values), "X1");
75 MDOVar var0 = m.addVar(0.0, MDO_INFINITY, 0.0, 'C', "x0");
76 MDOVar var1 = m.addVar(0.0, MDO_INFINITY, 0.0, 'C', "x1");
77
78 /* Add constraints. */
79 MDOPsdExpr psd_expr0 = MDOPsdExpr(psd_var0, MDOMatrix::coo(dim_mat[0], dim_mat[0], A0_nnz, A0_nz_indices, A0_nz_values));
80 const MDOVar row0_vars[] = { var0 };
81 psd_expr0.addTerms(row0_values, row0_vars, 1);
82 m.addPsdConstr(psd_expr0, MDO_EQUAL, 1.0, "C0");
83
84 MDOPsdExpr psd_expr1 = MDOPsdExpr(psd_var1, MDOMatrix::coo(dim_mat[0], dim_mat[0], A1_nnz, A1_nz_indices, A1_nz_values));
85 const MDOVar row1_vars[] = { var1 };
86 psd_expr1.addTerms(row1_values, row1_vars, 1);
87 m.addPsdConstr(psd_expr1, MDO_EQUAL, 2.0, "C1");
88
89 /* Set objective function. */
90 MDOLinExpr objective = 0;
91 m.setObjective(objective, -1);
92
93 /*------------------------------------------------------------------*/
94 /* Step 3. Solve the problem and populate optimization result. */
95 /*------------------------------------------------------------------*/
96 // Optimize the model
97 m.optimize();
98 if (m.get(MDO_IntAttr_Status) == MDO_OPTIMAL)
99 {
100 cout << "Primal objective val: " << m.get(MDO_DoubleAttr_PrimalObjVal) << endl;
101 }
102 else
103 {
104 cout << "No feasible solution." << endl;
105 }
106 }
107 catch (MDOException &e)
108 {
109 cout << "Error code = " << e.getErrorCode() << endl;
110 cout << e.getMessage() << endl;
111 }
112 catch (...)
113 {
114 cout << "Error during optimization." << endl;
115 }
116
117 return static_cast<int>(MDO_OKAY);
118}