0%

【漫漫科研路\C&C++】CPLEX解SOCP问题

IBM CPLEX可以解SOCP问题,但是需要先将这个SOCP问题化为指定的格式。本文首先介绍SOCP问题,然后举例介绍如何将SOCP问题转化为CPLEX认可的输入格式并求解。

SOCP的介绍

关于SOCP问题的介绍,可以参考Boyd等人写的Convex Optimization 或者是维基百科的SOCP词条 。这里摘录Convex Optimization一书中关于SOCP的定义:

图1

CPLEX求解SOCP问题

使用CPLEX求解SOCP问题,一般需要将问题转化为CPLEX可以识别的格式。CPLEX的例子ilosocpex1(位于安装目录的examples文件夹内,例如:C:\Program Files\IBM\ILOG\CPLEX_Enterprise_Server129\CPLEX_Studio\cplex\examples)给予了格式说明:
图2
下面我们首先给一个直接可以使用CPLEX求解的例子,然后在此基础上考虑一个更一般的例子(需要变量替换来符合格式)。

一个简单的例子

图3
如上图所示,q1, q2可以直接转化为前面提到的CPLEX认可格式。下面给出源代码如下(注意项目需要预先配置好,配置请见上一篇博文):

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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include <cmath>
#include <string>
#include <iostream>
#include <ilcplex/ilocplex.h>

ILOSTLBEGIN // import namespace std

// Tolerance for testing KKT conditions.
#define TESTTOL 1e-9
// Tolerance for barrier convergence.
#define CONVTOL 1e-9

// A Simple Example
// Minimize
// obj: x1 + x2 + x3 + x4 + x5 + x6
// Subject To
// c1: x1 + x2 + x5 = 8
// c2: x3 + x5 + x6 = 10
// q1: x1 >= |(x2, x3)| ---->>>> q1: [ -x1^2 + x2^2 + x3^2 ] <= 0 and x1 >=0
// q2: x4 >= |x5| ---->>>> q2: [ -x4^2 + x5^2 ] <= 0 and x4 >=0
// Bounds
// x2 Free
// x3 Free
// x5 Free
// x6 Free
// End

static void
createmodel(IloModel& model, IloObjective &obj, IloNumVarArray &x,
IloRangeArray &rngs, IloIntArray& cone)
{
IloEnv env = model.getEnv();

// Create variables.
x.add(IloNumVar(env, 0, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, 0, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));

// Create objective function and immediately store it in return value.
// obj = IloMinimize(env, x1 + x2 + x3 + x4 + x5 + x6);
obj = IloMinimize(env, x[0]+x[1] + x[2] + x[3] + x[4] + x[5]);
// Setup model.
model.add(x[0] + x[1] + x[4]==8);
model.add(x[2] + x[4] + x[5]==10);
model.add(-x[0] * x[0] + x[1] * x[1] + x[2] * x[2]<=0);

//equal to model.add(-x[3] * x[3] + x[4] * x[4]<=0), useful for lots of variables
double a[] = {0,0,0, -1,1 };
IloExpr temp(env);
for (IloInt i = 3; i < 5; i++)
{
temp += a[i] * x[i] * x[i];
}
model.add(temp <= 0);
temp.end();

model.add(obj);
}

int
main(void)
{
IloEnv env;
int retval = -1;

try {
// Create the model.
IloModel model(env);
IloCplex cplex(env);
IloObjective obj(env);
IloNumVarArray vars(env);
IloRangeArray rngs(env);
IloIntArray cone(env);
createmodel(model, obj, vars, rngs, cone);

// Extract model.
cplex.extract(model);

// Solve the problem. If we cannot find an _optimal_ solution then
// there is no point in checking the KKT conditions and we throw an
// exception.
cplex.setParam(IloCplex::Param::Barrier::QCPConvergeTol, CONVTOL);
if (!cplex.solve() || cplex.getStatus() != IloAlgorithm::Optimal)
throw string("Failed to solve problem to optimality");

IloNumArray vals_x(env);
env.out() << "Solution status = " << cplex.getStatus() << endl;
env.out() << "Solution value = " << cplex.getObjValue() << endl;
cplex.getValues(vals_x, vars);
env.out() << "Values = " << vals_x << endl;
env.end();
}
catch (IloException &e) {
cerr << "IloException: " << e << endl;
if (env.getImpl())
env.end();
::abort();
}
catch (string& e) {
cerr << e << endl;
if (env.getImpl())
env.end();
::abort();
}
return retval;
}

运行结果如下图:

图4


一个更一般的例子

在前面例子的基础上,我们只是改变了约束q1,使其更一般化,如下图所示:
图5

为了将q1转化为合适的格式,我们使用变量替换 $x_7=x_1+x_2 $ 。因此只需在前面源代码中更改 createmodel 函数中的部分代码。为保持代码完整性,我们依旧给出完整的代码:

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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include <cmath>
#include <string>
#include <iostream>
#include <ilcplex/ilocplex.h>

ILOSTLBEGIN // import namespace std

// Tolerance for testing KKT conditions.
#define TESTTOL 1e-9
// Tolerance for barrier convergence.
#define CONVTOL 1e-9

// A Simple Example
// Minimize
// obj: x1 + x2 + x3 + x4 + x5 + x6
// Subject To
// c1: x1 + x2 + x5 = 8
// c2: x3 + x5 + x6 = 10
// q1: x1 + x2 >= |(x1, x2, x3)| ---->>>> q1:[-x7^2+x1^2+x2^2+x3^2]<=0 and x7=x1+x2>=0
// q2: x4 >= |x5| ---->>>> q2: [ -x4^2 + x5^2 ] <= 0 and x4 >=0
// Bounds
// x1 Free
// x2 Free
// x3 Free
// x5 Free
// x6 Free
// End

static void
createmodel(IloModel& model, IloObjective &obj, IloNumVarArray &x,
IloRangeArray &rngs, IloIntArray& cone)
{
IloEnv env = model.getEnv();

// Create variables.
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, 0, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, -IloInfinity, IloInfinity));
x.add(IloNumVar(env, 0, IloInfinity));// x7

// Create objective function and immediately store it in return value.
// obj = IloMinimize(env, x1 + x2 + x3 + x4 + x5 + x6);
obj = IloMinimize(env, x[0]+x[1] + x[2] + x[3] + x[4] + x[5]);
// Setup model.
model.add(x[0] + x[1] + x[4]==8);
model.add(x[2] + x[4] + x[5]==10);
model.add(x[6] - x[0] - x[1] == 0); // x7 = x1 + x2
model.add(-x[6] * x[6]+ x[0] * x[0] + x[1] * x[1] + x[2] * x[2]<=0); //[-x7^2+x1^2+x2^2+x3^2]<=0

//equal to model.add(-x[3] * x[3] + x[4] * x[4]<=0), useful for lots of variables
double a[] = {0,0,0, -1,1 };
IloExpr temp(env);
for (IloInt i = 3; i < 5; i++)
{
temp += a[i] * x[i] * x[i];
}
model.add(temp <= 0);
temp.end();

model.add(obj);
}

int
main(void)
{
IloEnv env;
int retval = -1;

try {
// Create the model.
IloModel model(env);
IloCplex cplex(env);
IloObjective obj(env);
IloNumVarArray vars(env);
IloRangeArray rngs(env);
IloIntArray cone(env);
createmodel(model, obj, vars, rngs, cone);

// Extract model.
cplex.extract(model);

// Solve the problem. If we cannot find an _optimal_ solution then
// there is no point in checking the KKT conditions and we throw an
// exception.
cplex.setParam(IloCplex::Param::Barrier::QCPConvergeTol, CONVTOL);
if (!cplex.solve() || cplex.getStatus() != IloAlgorithm::Optimal)
throw string("Failed to solve problem to optimality");

IloNumArray vals_x(env);
env.out() << "Solution status = " << cplex.getStatus() << endl;
env.out() << "Solution value = " << cplex.getObjValue() << endl;
cplex.getValues(vals_x, vars);
env.out() << "Values = " << vals_x << endl;
env.end();
}
catch (IloException &e) {
cerr << "IloException: " << e << endl;
if (env.getImpl())
env.end();
::abort();
}
catch (string& e) {
cerr << e << endl;
if (env.getImpl())
env.end();
::abort();
}
return retval;
}

运行结果如下图:
图6

坚持原创技术分享,您的支持将鼓励我继续创作!