CODEWITHSUNDEEP
drake
Yuki
Yuki

Reputation: 139

drake gurobi cannot handle quadratic constraint

I am currently working to implement a simple quadratic constraint on an optimization problem. Gurobi's website says that quadratic constraints can be implemented. Does drake not have an interface with this constraint using Gurobi?

The code is as follows.

#include <Eigen/Dense>
#include <math.h>
#include <iostream>

#include "drake/solvers/mathematical_program.h"
#include "drake/solvers/solve.h"
#include "drake/solvers/gurobi_solver.h"
#include "drake/solvers/scs_solver.h"

using namespace std;
using namespace Eigen;
using namespace drake;

int main(){
    solvers::MathematicalProgram prog; 

    // Instantiate the decision variables
    solvers::VectorXDecisionVariable x = prog.NewContinuousVariables(2, "x");

    // Define constraints 
    for(int i = 0; i < 2; i++){
      prog.AddConstraint(x[i]*x[i] <= 2); //Replacing this with a linear constraint 
                                          //such as prog.AddLinearConstraint(-5 <= x[i] && x[i] <= 5); 
                                          //will work
    }

    // Define the cost
    MatrixXd Q(2,2);
    Q << 2, 1,
         1, 4;
    VectorXd b(2);
    b << 0, 3;
    double c = 1.0;
    prog.AddQuadraticCost(Q, b, c, x);

    solvers::GurobiSolver solver;
    cout << "Gurobi available? " << solver.is_enabled() << endl;

    auto result = solver.Solve(prog);
    cout << "Is optimization successful?" << result.is_success() << endl;
    cout << "Optimal x: " << result.GetSolution().transpose() << endl;
    cout << "solver is: " << result.get_solver_id().name() << endl;
    cout << "computation time is: " << result.get_solver_details<solvers::GurobiSolver>().optimizer_time;

    return 0;
}

Upvotes: 5

Views: 218

Answers (2)

Hongkai Dai
Hongkai Dai

Reputation: 2766

Drake currently (intentionally) doesn't support the quadratic constraint. For better solver performance, I would recommend to reformulate the optimization problem without quadratic constraint, but with second-order cone constraints.

First for a quadratic constraint

xᵀx <= a²

This can be regarded as a Lorentz cone constraint

(y, x) is in the Lorentz cone
y = a

In Drake, you could add this constraint as

prog.AddLorentzConeConstraint((drake::VectorX<drake::symbolic::Expression>(x.rows() + 1) << a, x).finished());

For the more general quadratic constraint 0.5xᵀQx + bᵀx + c<=0, we also provide AddQuadraticAsRotatedLorentzConeConstraint which imposes a rotated Lorentz cone constraint.

prog.AddQuadraticAsRotatedLorentzConeConstraint(Q, b, c, x);

For the quadratic cost

min 0.5xᵀQx + bᵀx + c

you could also reformulate it as a linear cost with a rotated Lorentz cone constraint

min z
z >= 0.5xᵀQx + bᵀx + c

In Drake, the code is

z = prog.NewContinuousVariables<1>()(0);
prog.AddLinearCost(z);
prog.AddRotatedLorentzConeConstraint(symbolic::Expression(z), symbolic::Expression(1), 0.5 * x.dot(Q * x) + b.dot(x) + c);

The reason why we would prefer Lorentz cone constraint, instead of a quadratic-constrained-quadratic-program (QCQP), is that with Lorentz cone constraint, you end up with a second-order conic optimization problem (SOCP), this optimization problem has a lot of nice features, in terms of its dual form. You could read the paper Alizadeh's paper for more details. And Mosek also recommends SOCP over QCQP https://docs.mosek.com/9.2/rmosek/prob-def-quadratic.html#a-recommendation

Upvotes: 3

jwnimmer-tri
jwnimmer-tri

Reputation: 2449

Possibly the (as-of-yet unimplemented) feature request at https://github.com/RobotLocomotion/drake/issues/6341 is what you're looking for?

Upvotes: 0

Related Questions