I am trying to port an application from Matlab to C++. This application contains a linear optimization problem of the form:
min. F * x
s.t. A * x < Rhs
s.t. Lb <= x <= Ub
which needs to be solved.
The matlab API of CPLEX supports writing all coefficients as vectors (F, Rhs, Lb, and Ub) or matrices (A) respectively. The C++ interface however does not. Instead the problem needs to be constructed by row or by column.
Here I am stuck since I can't seem to geht the same result in Matlab and C++ even for a small toy problem.
t_v_l = [1 0 0 1 0 0 1 0 0]';
t_v_m = [0 1 0 0 1 0 0 1 0]';
t_v_r = [0 0 1 0 0 1 0 0 1]';
t_h_o = [1 1 1 0 0 0 0 0 0]';
t_h_m = [0 0 0 1 1 1 0 0 0]';
t_h_u = [0 0 0 0 0 0 1 1 1]';
target = [0 1 0 1 1 1 0 1 0];
dict = [t_v_l, t_v_m, t_v_r, t_h_o, t_h_m, t_h_u];
M = size(dict, 1);
N = size(dict, 2);
F = zeros(N + M, 1);
F(N + 1 : N + M) = ones(M,1);
Lb = sparse(zeros(N + M, 1));
Lb(N + 1 : N + M) = 0;
Ub = ones(N + M, 1);
Ub(N + 1 : N + M) = 999;
Rhs = zeros(2 * M, 1);
Rhs( 1 : M) = -target;
Rhs(M + 1 : 2 * M) = target;
Lhs = -999 * ones(2 * M, 1);
A = zeros(2 * M, N + M);
A( 1 : M, 1 : N ) = -dict;
A(M + 1 : 2 * M, 1 : N ) = dict;
A( 1 : M, N + 1 : N + M) = -eye(M);
A(M + 1 : 2 * M, N + 1 : N + M) = -eye(M);
plex = Cplex('myModel');
plex.addCols(F, [], Lb, Ub);
plex.addRows(Lhs, A, Rhs);
plex.writeModel('ex.lp'); % Optional
plex.solve();
X = plex.Solution.x;
fprintf("Solution:\n")
disp(X');
Conceptually I am trying to reconstruct a binary target vector from different binary templates. The first 6 entries in the solution show which templates have to be chosen (in this case template 2 and 5). The remaining entries correspond to slack variables where the reconstruction is not entirely accurate due to overlap or noise in the target. The code above works and gives me the correct result.
However, when I try to follow the examples (like Ilodiet.cpp) to produce the same results in the Concert API the optimization does not work anymore.
#ifndef IL_STD
#define IL_STD
#endif
#include <cstring>
#include <ilcplex/ilocplex.h>
#include <vector>
ILOSTLBEGIN
void main()
{
// Input Data
std::vector<vector<float>> dict =
{
{ 1, 0, 0, 1, 0, 0, 1, 0, 0 },
{ 0, 1, 0, 0, 1, 0, 0, 1, 0 },
{ 0, 0, 1, 0, 0, 1, 0, 0, 1 },
{ 1, 1, 1, 0, 0, 0, 0, 0, 0 },
{ 0, 0, 0, 1, 1, 1, 0, 0, 0 },
{ 0, 0, 0, 0, 0, 0, 1, 1, 1 }
};
std::vector<float> target =
{ 0, 1, 0, 1, 1, 1, 0, 1, 0 };
const int M = target.size();
const int N = dict.size();
// Build Matrix A
std::vector<std::vector<float>> A;
for(int r = 0; r < M; ++r)
{
std::vector<float> row(N + M, 0);
for(int c = 0; c < N; ++c)
{
row[c] = -dict[c][r];
}
row[N + r] = -1;
A.push_back(row);
}
for (int r = M; r < 2 * M; ++r)
{
std::vector<float> row(N + M, 0);
for (int c = 0; c < N; ++c)
{
row[c] = dict[c][r - M];
}
row[N + r - M] = -1;
A.push_back(row);
}
// Build Model
IloEnv env;
IloModel mod(env);
IloNumArray F(env, N + M);
IloNumArray Lb(env, N + M);
IloNumArray Ub(env, N + M);
for (int i = 0; i < N; ++i)
{
F[i] = 0.f;
Lb[i] = 0.f;
Ub[i] = 1.f;
}
for (int i = N; i < N + M; ++i)
{
Lb[i] = 0.f;
Ub[i] = 999.f;// IloInfinity;
F[i] = 1.f;
}
IloNumArray Rhs(env, 2 * M);
IloNumArray Lhs(env, 2 * M);
for (int i = 0; i < M; ++i)
{
Rhs[i] = -target[i];
Lhs[i] = -999.f;// -IloInfinity;
Rhs[M + i] = target[i];
Lhs[M + i] = -999.f;// -IloInfinity;
}
IloNumVarArray x(env, Lb, Ub, ILOFLOAT);
IloObjective minimize = IloMinimize(env, IloScalProd(x, F));
mod.add(IloExtractable(minimize));
for (int r = 0; r < A.size(); ++r)
{
IloExpr expr(env);
for (int c = 0; c < A[r].size(); ++c) {
expr += F[c] * A[r][c];
}
mod.add(expr <= Rhs[r]);
expr.end();
}
IloCplex cplex(mod);
cplex.exportModel("ex.lp");
cplex.solve();
std::cout << "Solution status = " << cplex.getStatus() << endl;
std::cout << "Cost = " << cplex.getObjValue() << endl;
std::vector<int> ret;
for (int i = 0; i < N; i++)
{
try
{
std::cout << "Template i: " << cplex.getValue(x[i]) << std::endl;
}
catch(...)
{
std::cout << "Template i: 0 (optimized away)" << std::endl;
}
}
for (int i = N; i < N + M; i++)
{
try
{
std::cout << "Slack i: " << cplex.getValue(x[i]) << std::endl;
}
catch (...)
{
std::cout << "Slack i: 0 (optimized away)" << std::endl;
}
}
env.end();
}
The output tells me that the presolver eliminated all rows and columns and my solution is empty. The exported problem files also look completely different. If I take a look at the individual arrays however, they look exactly the same. So I guess I am doing something wrong while passing my information to the Concert API.
Where is my mistake? Any help is appreciated!
The code for printing the variables (since IloArrays are hard to look at in the debugger) in Matlab:
fprintf("F:\n");
disp(F');
fprintf("Lb:\n");
disp(Lb');
fprintf("Ub:\n");
disp(Ub');
fprintf("Lhs:\n");
disp(Lhs');
fprintf("Rhs:\n");
disp(Rhs');
fprintf("A:\n");
disp(A);
And in C++:
std::cout << "F:" << std::endl;
for (int i = 0; i < N + M; ++i)
{
std::cout << F[i] << ", ";
}
std::cout << "\n\nLb:" << std::endl;
for (int i = 0; i < N + M; ++i)
{
std::cout << Lb[i] << ", ";
}
std::cout << "\n\nUb:" << std::endl;
for (int i = 0; i < N + M; ++i)
{
std::cout << Ub[i] << ", ";
}
std::cout << "\n\nLhs:" << std::endl;
for (int i = 0; i < 2 * M; ++i)
{
std::cout << Lhs[i] << ", ";
}
std::cout << "\n\nRhs:" << std::endl;
for (int i = 0; i < 2 * M; ++i)
{
std::cout << Rhs[i] << ", ";
}
std::cout << "\n\nA:" << std::endl;
for (int r = 0; r < 2 * M; ++r)
{
for(int c = 0; c < N + M; ++c)
{
std::cout.width(2);
std::cout.fill(' ');
std::cout << A[r][c] << ", ";
}
std::cout << "\n";
}
