5.4.4. C++的SDP建模与优化¶
在本节中，我们将使用 MindOpt C++ 语言的 API 来建模以及求解 半定规划问题示例 中的问题。
5.4.4.1. SDP示例1¶
首先，需要引入头文件：
23#include "MindoptCpp.h"
第一步：创建优化模型：
先建立一空的MindOpt模型。
47 MdoModel model;
第二步：SDP模型输入
接着，我们将目标函数改为maximization，并新增 \(\mathbf{X}\) 对称矩阵变量。
51 /*------------------------------------------------------------------*/
52 /* Step 2. Input model. */
53 /*------------------------------------------------------------------*/
54 /* Change to maximization problem. */
55 model.setIntAttr(MDO_INT_ATTR::MIN_SENSE, MDO_NO);
56
57 /* Add matrix variable. */
58 model.addSymMat(dim_mat, mat_name);
我们利用 mindopt::MdoModel::replaceSymMatObjs() 来输入目标函数中的 \(\mathbf{C}\) 的非零元素。
其中，第一个参数为 矩阵的索引值 (此处为0)，第二个参数为 矩阵中非零元素的个数，第三个和第四个参数分别代表 矩阵中非零元素其对应的行和
列之索引 ，而最后一个参数则为 非零元数值的索引。
60 /* Input objective coefficients. */
61 model.replaceSymMatObjs(0, C_size, C_row_indices, C_col_indices, C_values);
接着，我们输入第一个约束。我们首先输一空约束；接着再利用 mindopt::MdoModel::replaceSymMatElements() 输入约束中的
\(\mathbf{A}\) 的非零元素。其中，第一和第二个参数分别为： 约束的索引值 (此处为0)以及 矩阵的索引值 (此处为0)，第三个参数为 矩阵的非零元素的个数，第四个和第五个参数代表 矩阵中非零元素其对应的行和列之索引，而最后一个参数则为 非零元的数值之索引。
63 /* Input first constraint. */
64 model.addCons(1.0, 1.0, "c0");
65 model.replaceSymMatElements(
66 0, 0, A_size, A_row_indices_in_cone, A_col_indices_in_cone, A_values_in_cone);
第三步：求解SDP模型
模型输入后，我们接着用 mindopt::MdoModel::solveProb() 来求解问题，并用 mindopt::MdoModel::displayResults() 呈现求解结果。
71 /* Solve the problem. */
72 model.solveProb();
73 model.displayResults();
第四步: 取得SDP模型的解
我们利用泛用函数 mindopt::MdoModel::getRealAttr() 来取得最优的的目标函数值
82 std::cout << "Optimizer terminated with an OPTIMAL status" << std::endl;
83 std::cout << std::setprecision(6) << std::cout.width(8) << std::fixed << std::scientific <<
84 " - Primal objective : " << model.getRealAttr(MDO_REAL_ATTR::PRIMAL_OBJ_VAL) << std::endl;
最后，我们利用 mindopt::MdoModel::getRealAttrSymMat() 来取得 \(\mathbf{X}\) 的变量值。其中，
第一个参数定义 属性值 (附注: mindopt::MdoModel::getRealAttrSymMat() 为泛用函数，因此必须配合
属性值来定义函数的使用目的；MDO_REAL_ATTR::SYM_MAT_PRIMAL_SOLN定义目的为
获取(对称矩阵)变量的数值)，第二个参数为： 矩阵的索引值 (此处为0)，第三个参数为 矩阵的尺寸。
85 auto soln = model.getRealAttrSymMat(
86 MDO_REAL_ATTR::SYM_MAT_PRIMAL_SOLN, 0, dim_mat * dim_mat, NULL, NULL);
文件链接 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 * Matrix
16 * C = [ -3 0 1 ] A = [ 3 0 1 ]
17 * [ 0 -2 0 ] [ 0 4 0 ]
18 * [ 1 0 -3 ] [ 1 0 5 ]
19 * End
20 */
21#include <iostream>
22#include <iomanip>
23#include "MindoptCpp.h"
24
25using namespace mindopt;
26
27int main(void)
28{
29 const std::string mat_name = "X";
30 const int num_mats = 1;
31 const int dim_mat = 3; /* Dimension of the matrix variables. */
32
33 const int C_size = 4;
34 const int C_row_indices[] = { 0, 0, 1, 2 }; /* Row index of a matrix variable. */
35 const int C_col_indices[] = { 0, 2, 1, 2 }; /* Column index of a matrix variable. */
36 const double C_values[] = { -3.0, 1.0, -2.0, -3.0 }; /* Values of a matrix variable. */
37
38 const int A_size = 4;
39 const int A_row_indices_in_cone[] = { 0, 2, 1, 2 }; /* Row index in a matrix variable. */
40 const int A_col_indices_in_cone[] = { 0, 0, 1, 2 }; /* Column index in a matrix variable. */
41 const double A_values_in_cone[] = { 3.0, 1.0, 4.0, 5.0 }; /* Values of a matrix variable. */
42
43 /*------------------------------------------------------------------*/
44 /* Step 1. Create a model and change the parameters. */
45 /*------------------------------------------------------------------*/
46 /* Create an empty model. */
47 MdoModel model;
48
49 try
50 {
51 /*------------------------------------------------------------------*/
52 /* Step 2. Input model. */
53 /*------------------------------------------------------------------*/
54 /* Change to maximization problem. */
55 model.setIntAttr(MDO_INT_ATTR::MIN_SENSE, MDO_NO);
56
57 /* Add matrix variable. */
58 model.addSymMat(dim_mat, mat_name);
59
60 /* Input objective coefficients. */
61 model.replaceSymMatObjs(0, C_size, C_row_indices, C_col_indices, C_values);
62
63 /* Input first constraint. */
64 model.addCons(1.0, 1.0, "c0");
65 model.replaceSymMatElements(
66 0, 0, A_size, A_row_indices_in_cone, A_col_indices_in_cone, A_values_in_cone);
67
68 /*------------------------------------------------------------------*/
69 /* Step 3. Solve the problem and populate the result. */
70 /*------------------------------------------------------------------*/
71 /* Solve the problem. */
72 model.solveProb();
73 model.displayResults();
74
75 switch (model.getStatus())
76 {
77 case MDO_UNKNOWN:
78 std::cout <<"Optimizer terminated with an UNKNOWN status." << std::endl;
79 break;
80 case MDO_OPTIMAL:
81 {
82 std::cout << "Optimizer terminated with an OPTIMAL status" << std::endl;
83 std::cout << std::setprecision(6) << std::cout.width(8) << std::fixed << std::scientific <<
84 " - Primal objective : " << model.getRealAttr(MDO_REAL_ATTR::PRIMAL_OBJ_VAL) << std::endl;
85 auto soln = model.getRealAttrSymMat(
86 MDO_REAL_ATTR::SYM_MAT_PRIMAL_SOLN, 0, dim_mat * dim_mat, NULL, NULL);
87 std::cout << "X = " << std::endl;
88 for (auto i = 0; i < dim_mat; ++i)
89 {
90 std::cout << " (";
91 for (auto j = 0; j < dim_mat; ++j)
92 {
93 std::cout << " " << soln[i * dim_mat + j];
94 }
95 std::cout << " )" << std::endl;
96 }
97 break;
98 }
99 case MDO_INFEASIBLE:
100 std::cout << "Optimizer terminated with an INFEASIBLE status." << std::endl;
101 break;
102 case MDO_UNBOUNDED:
103 std::cout << "Optimizer terminated with an UNBOUNDED status." << std::endl;
104 break;
105 case MDO_INF_OR_UBD:
106 std::cout << "Optimizer terminated with an INFEASIBLE or UNBOUNDED status." << std::endl;
107 break;
108 }
109 }
110 catch (MdoException & e)
111 {
112 std::cerr << "===================================" << std::endl;
113 std::cerr << "Error : code <" << e.getResult() << ">" << std::endl;
114 std::cerr << "Reason : " << model.explainResult(e.getResult()) << std::endl;
115 std::cerr << "===================================" << std::endl;
116
117 return static_cast<int>(e.getResult());
118 }
119
120 return static_cast<int>(MDO_OKAY);
121}
5.4.4.2. SDP示例2¶
首先，需要引入头文件：
30#include "MindoptCpp.h"
第一步：创建MindOpt模型
首先，我们必须先建立一空的MindOpt模型。
67 /* Create an empty model. */
68 MdoModel model;
第二步：SDP模型输入
接着，我们将目标函数改为maximization，并新增 \(x_0\) 和 \(x_1\) 两个变量，以及 \(\mathbf{X}_0\) 和 \(\mathbf{X}_1\) 两个对称矩阵变量。
72 /*------------------------------------------------------------------*/
73 /* Step 2. Input model. */
74 /*------------------------------------------------------------------*/
75 /* Change to maximization problem. */
76 model.setIntAttr(MDO_INT_ATTR::MIN_SENSE, MDO_NO);
77
78 /* Add variables. */
79 std::vector<MdoVar> x;
80 x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, "x0", MDO_NO));
81 x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, "x1", MDO_NO));
82
83 /* Add matrix variables. */
84 model.addSymMats(num_mats, dim_mats, mat_names);
我们利用 mindopt::MdoModel::replaceSymMatObjs() 来输入目标函数中的 \(\mathbf{C}_0\) 的非零元素。
其中，第一个参数为 矩阵的索引值 (此处为C0_mat_idx)，第二个参数为 矩阵的非零元素个数， 第三和第四个参数分别代表 矩阵中非零元素其对应的行和列之索引 ，
而最后一个参数则为 非零元的数值索引 。
86 /* Input objective coefficients. */
87 model.replaceSymMatObjs(C0_mat_idx, C0_size, C0_row_indices, C0_col_indices, C0_values);
同理，我们再接着输入目标函数中的 \(\mathbf{C}_1\) 的非零元素。
88 model.replaceSymMatObjs(C1_mat_idx, C1_size, C1_row_indices, C1_col_indices, C1_values);
接着，我们输入第一个约束。我们首先输入空约束，接着再利用 mindopt::MdoModel::replaceSymMatElements() 输入约束中的 \(\mathbf{A}_{00}\) 的非零元素。其中，第一和第二个参数分别为： 约束的索引值 (此处为0)以及 矩阵的索引值 (此处为 A00_mat_idx )，第三个参数为 矩阵的非零元素的个数，第四个和第五个参数代表 矩阵中非零元素其对应的行和列之索引，而最后一个参数则为 非零元的数值之索引。
90 /* Input first constraint. */
91 model.addCons(1.0 * x[0] == 1.0, "c0");
92 model.replaceSymMatElements(
93 0, A00_mat_idx, A00_size, A00_row_indices_in_cone, A00_col_indices_in_cone, A00_values_in_cone);
同理，我们再接着输入第二个约束。
95 /* Input second constraint. */
96 model.addCons(1.0 * x[1] == 2.0, "c1");
97 model.replaceSymMatElements(
98 1, A11_mat_idx, A11_size, A11_row_indices_in_cone, A11_col_indices_in_cone, A11_values_in_cone);
第三步：求解SDP模型
模型输入后，我们接着用 mindopt::MdoModel::solveProb() 来求解问题，并用 mindopt::MdoModel::displayResults() 呈现求解结果。
100 /*------------------------------------------------------------------*/
101 /* Step 3. Solve the problem and populate the result. */
102 /*------------------------------------------------------------------*/
103 /* Solve the problem. */
104 model.solveProb();
105 model.displayResults();
第四步: 取得SDP模型的解
我们利用泛用函数 mindopt::MdoModel::getRealAttr() 来
取得最优的的目标函数值，并用泛用函数 mindopt::MdoModel::getRealAttrArray() 来
取得 \(x_0\) 和 \(x_1\) 两个变量值；此处第一个参数定义了 属性值 (附注: mindopt::MdoModel::getRealAttrArray() 为泛用函数，
因此必须配合属性值来定义函数的使用目的；MDO_REAL_ATTR::PRIMAL_SOLN定义目的为获取(线性)变量的数值)，
第二和第三个参数分别代表第一个索引和最后一个索引 加上1，
也就是 \(x_j,\) 其中 \(0 \leq j < 2\)。
115 std::cout << std::setprecision(6) << std::cout.width(8) << std::fixed << std::scientific <<
116 " - Primal objective : " << model.getRealAttr(MDO_REAL_ATTR::PRIMAL_OBJ_VAL) << std::endl;
最后，我们利用 mindopt::MdoModel::getRealAttrSymMat() 来取得 \(\mathbf{X}_0\) 和 \(\mathbf{X}_1\) 的变量值。其中，
第一个参数定义 属性值 (附注: mindopt::MdoModel::getRealAttrSymMat() 为泛用函数，因此必须配合
属性值来定义函数的使用目的；MDO_REAL_ATTR::SYM_MAT_PRIMAL_SOLN定义目的为
获取(对称矩阵)变量的数值)，第二个参数为： 矩阵的索引值 (此处为b)，第三个参数为 矩阵的尺寸。
124 auto soln = model.getRealAttrSymMat(
125 MDO_REAL_ATTR::SYM_MAT_PRIMAL_SOLN, b, dim_mats[b] * dim_mats[b], NULL, NULL);
文件链接 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 * Matrix
20 * C0 = [ 2 1 ] A00 = [ 3 1 ]
21 * [ 1 2 ] [ 1 3 ]
22 *
23 * C0 = [ 3 0 1 ] A00 = [ 3 0 1 ]
24 * [ 0 2 0 ] [ 0 4 0 ]
25 * [ 1 0 3 ] [ 1 0 5 ]
26 * End
27 */
28#include <iostream>
29#include <iomanip>
30#include "MindoptCpp.h"
31
32using namespace mindopt;
33
34int main(void)
35{
36 const int num_mats = 2;
37 const int dim_mats[] = { 2, 3 }; /* Dimension of the matrix variables. */
38 const std::string mat_names[] = { "X0", "X1" };
39
40 const int C0_size = 3;
41 const int C0_mat_idx = 0;
42 const int C0_row_indices[] = { 0, 1, 1 }; /* Row index of a matrix variable. */
43 const int C0_col_indices[] = { 0, 0, 1 }; /* Column index of a matrix variable. */
44 const double C0_values[] = { 2.0, 1.0, 2.0 }; /* Values of a matrix variable. */
45
46 const int C1_size = 4;
47 const int C1_mat_idx = 1;
48 const int C1_row_indices[] = { 0, 0, 1, 2 }; /* Row index of a matrix variable. */
49 const int C1_col_indices[] = { 0, 2, 1, 2 }; /* Column index of a matrix variable. */
50 const double C1_values[] = { 3.0, 1.0, 2.0, 3.0 }; /* Values of a matrix variable. */
51
52 const int A00_size = 3;
53 const int A00_mat_idx = 0;
54 const int A00_row_indices_in_cone[] = { 0, 1, 1 }; /* Row index in a matrix variable. */
55 const int A00_col_indices_in_cone[] = { 0, 0, 1 }; /* Column index in a matrix variable. */
56 const double A00_values_in_cone[] = { 3.0, 1.0, 3.0 }; /* Values of a matrix variable. */
57
58 const int A11_size = 4;
59 const int A11_mat_idx = 1;
60 const int A11_row_indices_in_cone[] = { 0, 2, 1, 2 }; /* Row index in a matrix variable. */
61 const int A11_col_indices_in_cone[] = { 0, 0, 1, 2 }; /* Column index in a matrix variable. */
62 const double A11_values_in_cone[] = { 3.0, 1.0, 4.0, 5.0 }; /* Values of a matrix variable. */
63
64 /*------------------------------------------------------------------*/
65 /* Step 1. Create a model and change the parameters. */
66 /*------------------------------------------------------------------*/
67 /* Create an empty model. */
68 MdoModel model;
69
70 try
71 {
72 /*------------------------------------------------------------------*/
73 /* Step 2. Input model. */
74 /*------------------------------------------------------------------*/
75 /* Change to maximization problem. */
76 model.setIntAttr(MDO_INT_ATTR::MIN_SENSE, MDO_NO);
77
78 /* Add variables. */
79 std::vector<MdoVar> x;
80 x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, "x0", MDO_NO));
81 x.push_back(model.addVar(0.0, MDO_INFINITY, 0.0, "x1", MDO_NO));
82
83 /* Add matrix variables. */
84 model.addSymMats(num_mats, dim_mats, mat_names);
85
86 /* Input objective coefficients. */
87 model.replaceSymMatObjs(C0_mat_idx, C0_size, C0_row_indices, C0_col_indices, C0_values);
88 model.replaceSymMatObjs(C1_mat_idx, C1_size, C1_row_indices, C1_col_indices, C1_values);
89
90 /* Input first constraint. */
91 model.addCons(1.0 * x[0] == 1.0, "c0");
92 model.replaceSymMatElements(
93 0, A00_mat_idx, A00_size, A00_row_indices_in_cone, A00_col_indices_in_cone, A00_values_in_cone);
94
95 /* Input second constraint. */
96 model.addCons(1.0 * x[1] == 2.0, "c1");
97 model.replaceSymMatElements(
98 1, A11_mat_idx, A11_size, A11_row_indices_in_cone, A11_col_indices_in_cone, A11_values_in_cone);
99
100 /*------------------------------------------------------------------*/
101 /* Step 3. Solve the problem and populate the result. */
102 /*------------------------------------------------------------------*/
103 /* Solve the problem. */
104 model.solveProb();
105 model.displayResults();
106
107 switch (model.getStatus())
108 {
109 case MDO_UNKNOWN:
110 std::cout <<"Optimizer terminated with an UNKNOWN status." << std::endl;
111 break;
112 case MDO_OPTIMAL:
113 {
114 std::cout << "Optimizer terminated with an OPTIMAL status" << std::endl;
115 std::cout << std::setprecision(6) << std::cout.width(8) << std::fixed << std::scientific <<
116 " - Primal objective : " << model.getRealAttr(MDO_REAL_ATTR::PRIMAL_OBJ_VAL) << std::endl;
117 auto soln = model.getRealAttrArray(MDO_REAL_ATTR::PRIMAL_SOLN, 0, 2);
118 for (auto j = 0; j < 2; ++j)
119 {
120 std::cout << "x[" << j << "] = " << soln[j] << std::endl;
121 }
122 for (auto b = 0; b < num_mats; ++b)
123 {
124 auto soln = model.getRealAttrSymMat(
125 MDO_REAL_ATTR::SYM_MAT_PRIMAL_SOLN, b, dim_mats[b] * dim_mats[b], NULL, NULL);
126 std::cout << "X[" << b<< "] = " << std::endl;
127 for (auto i = 0; i < dim_mats[b]; ++i)
128 {
129 std::cout << " (";
130 for (auto j = 0; j < dim_mats[b]; ++j)
131 {
132 std::cout << " " << soln[i * dim_mats[b] + j];
133 }
134 std::cout << " )" << std::endl;
135 }
136 }
137 break;
138 }
139 case MDO_INFEASIBLE:
140 std::cout << "Optimizer terminated with an INFEASIBLE status." << std::endl;
141 break;
142 case MDO_UNBOUNDED:
143 std::cout << "Optimizer terminated with an UNBOUNDED status." << std::endl;
144 break;
145 case MDO_INF_OR_UBD:
146 std::cout << "Optimizer terminated with an INFEASIBLE or UNBOUNDED status." << std::endl;
147 break;
148 }
149 }
150 catch (MdoException & e)
151 {
152 std::cerr << "===================================" << std::endl;
153 std::cerr << "Error : code <" << e.getResult() << ">" << std::endl;
154 std::cerr << "Reason : " << model.explainResult(e.getResult()) << std::endl;
155 std::cerr << "===================================" << std::endl;
156
157 return static_cast<int>(e.getResult());
158 }
159
160 return static_cast<int>(MDO_OKAY);
161}