CGAL QP solver versus Matlab quadprog()

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

CGAL QP solver versus Matlab quadprog()

DaniCamps
Dear all,

I am finding that CGAL returns INFEASIBLE for a problem that Matlab
quadprog() can solve. I do not know if this is normal, or I am making a
mistake somewhere. I leave here the input of the program and would
appreciate if anyone can give me a hint.

BR

Daniel

Input of the program
=============

D
==
D[0][0]=0.000050
D[0][1]=0.000050
D[0][2]=0.000000
D[0][3]=0.000000
D[1][0]=0.000050
D[1][1]=0.000050
D[1][2]=0.000000
D[1][3]=0.000000
D[2][0]=0.000000
D[2][1]=0.000000
D[2][2]=0.000050
D[2][3]=0.000050
D[3][0]=0.000000
D[3][1]=0.000000
D[3][2]=0.000050
D[3][3]=0.000050

c:
==
c[0][0]=-0.002359
c[1][0]=-0.002359
c[2][0]=-0.002399
c[3][0]=-0.002399

A:
==
A[0][0]=1.000000
A[0][1]=0.000000
A[0][2]=1.000000
A[0][3]=0.000000
A[1][0]=0.000000
A[1][1]=1.000000
A[1][2]=0.000000
A[1][3]=1.000000

b:
==
b[0][0]=1.000000
b[1][0]=1.000000

ub:
==
ub[0][0]=0.914808
ub[1][0]=0.933087
ub[2][0]=1.061726
ub[3][0]=1.028364

lb:
==
lb[0][0]=0.000000
lb[1][0]=0.000000
lb[2][0]=0.000000
lb[3][0]=0.000000

The problem has equality constraints. Matlab quadprog() returns the
following solution:

Cost function
========
fval:[0][0]=-0.004717

x:
=
x[0][0]=0.293237
x[1][0]=0.306496
x[2][0]=0.706763
x[3][0]=0.693504

I configure the program in CGAL as follows:

// choose exact integral type
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpzf.h>
typedef CGAL::Gmpzf ET;
#else
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
#endif

// program and solution types
typedef CGAL::Quadratic_program<double> Program;
typedef CGAL::Quadratic_program_solution<ET> Solution;

Solution s = CGAL::solve_quadratic_program(qp, ET());







--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: CGAL QP solver versus Matlab quadprog()

DaniCamps
Hi,

After some tweaking I get now the following result:

status:          OPTIMAL
objective value: -1.27827e-26/2.71051e-24
variable values:
  0: 0/1
  1: 6e-05/0.0001
  2: 1.35525e-20/1.35525e-20
  3: 4e-05/0.0001

It turns out that the objective value is correct (the same I get through
Matlab), i.e. -0.00471, but the variable values are wrong. In fact if you
evaluate the objective function (x'*D*x + c'*x) at the x = (0; 6e-9; 1;
4e-9) you get a value of -0.0024, not corresponding to the optimal one.
Instead the solution delivered by Matlab quadprog() does correspond to the
optimal value of the cost function.

Any help is appreciated.

BR

Daniel

PS: I post here my program

// example: construct a quadratic program from data
// the QP below is the first quadratic program example in the user manual
#include <iostream>
#include <cassert>
#include <CGAL/basic.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>

// choose exact integral type
#ifdef CGAL_USE_GMP
#include <CGAL/Gmpzf.h>
typedef CGAL::Gmpzf ET;
#else
#include <CGAL/MP_Float.h>
typedef CGAL::MP_Float ET;
#endif

// program and solution types
typedef CGAL::Quadratic_program<double> Program;
typedef CGAL::Quadratic_program_solution<ET> Solution;

int main() {
  // by default, we have a nonnegative QP with Ax <= b
  Program qp (CGAL::EQUAL, true, 0, true, 0);

  //Setting matrix A
  qp.set_a(0,0,1);
  qp.set_a(1,0,0);
  qp.set_a(2,0,1);
  qp.set_a(3,0,0);
  qp.set_a(0,1,0);
  qp.set_a(1,1,1);
  qp.set_a(2,1,0);
  qp.set_a(3,1,1);

  //Setting constraints vector b
  qp.set_b(0,1);
  qp.set_b(1,1);

  // Setting matrix D
  qp.set_d(0,0,0.00005);
  qp.set_d(1,0,0.00005);
  qp.set_d(1,1,0.00005);
  qp.set_d(2,0,0);
  qp.set_d(2,1,0);
  qp.set_d(2,2,0.00005);
  qp.set_d(3,0,0);
  qp.set_d(3,1,0);
  qp.set_d(3,2,0.00005);
  qp.set_d(3,3,0.00005);

  //Setting vector c
  qp.set_c(0,-0.002359);
  qp.set_c(1,-0.002359);
  qp.set_c(2,-0.002399);
  qp.set_c(3,-0.002399);

  //Setting upper bounds
  qp.set_u(0, true, 0.914808);
  qp.set_u(1, true, 0.933087);
  qp.set_u(2, true, 1.061726);
  qp.set_u(3, true, 1.028364);

  //Setting lower bounds
  qp.set_l(0, true, 0);
  qp.set_l(1, true, 0);
  qp.set_l(2, true, 0);
  qp.set_l(3, true, 0);

  //Setting euality constraints
  qp.set_r(0, CGAL::EQUAL);
  qp.set_r(1, CGAL::EQUAL);


  Solution s = CGAL::solve_quadratic_program(qp, ET());
  assert (s.solves_quadratic_program(qp));

  // output solution
  std::cout << s;
  return 0;
}




--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss