diff --git a/CoMISo/NSolver/BoundConstraint.hh b/CoMISo/NSolver/BoundConstraint.hh
index c238a236418f85b95e04bbf67c2481d88d92142c..becb621cf62625a7cdeb7fa667ff98d2a64db9d9 100644
--- a/CoMISo/NSolver/BoundConstraint.hh
+++ b/CoMISo/NSolver/BoundConstraint.hh
@@ -57,6 +57,7 @@ public:
   virtual void   eval_gradient   ( const double* _x, SVectorNC& _g      ) { _g.resize(n_); _g.coeffRef(idx_) = 1.0; }
   virtual void   eval_hessian    ( const double* _x, SMatrixNC& _h      ) { _h.clear(); _h.resize(n_,n_); }
 
+  virtual bool   is_linear() { return true;}
 
   // set/get values
   unsigned int& idx()   {return idx_;}
diff --git a/CoMISo/NSolver/GUROBISolver.cc b/CoMISo/NSolver/GUROBISolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..522b2df8b249f0ebcc73962409f0848beb7c29dc
--- /dev/null
+++ b/CoMISo/NSolver/GUROBISolver.cc
@@ -0,0 +1,181 @@
+//=============================================================================
+//
+//  CLASS GUROBISolver - IMPLEMENTATION
+//
+//=============================================================================
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include <CoMISo/Config/config.hh>
+#if COMISO_GUROBI_AVAILABLE
+//=============================================================================
+
+
+#include "GUROBISolver.hh"
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ========================================================== 
+
+
+//// Constructor
+//GUROBISolver::
+//GUROBISolver()
+//{
+//
+//}
+
+//-----------------------------------------------------------------------------
+
+// ********** SOLVE **************** //
+bool
+GUROBISolver::
+solve(NProblemInterface*                  _problem,
+      std::vector<NConstraintInterface*>& _constraints,
+      std::vector<PairUiV>&               _discrete_constraints,
+      const double                        _time_limit)
+{
+  try
+  {
+    //----------------------------------------------
+    // 0. set up environment
+    //----------------------------------------------
+
+    GRBEnv   env   = GRBEnv();
+    GRBModel model = GRBModel(env);
+
+    model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
+
+
+    //----------------------------------------------
+    // 1. allocate variables
+    //----------------------------------------------
+
+    // determine variable types: 0->real, 1->integer, 2->bool
+    std::vector<char> vtypes(_problem->n_unknowns(),0);
+    for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
+      switch(_discrete_constraints[i].second)
+      {
+        case Integer: vtypes[_discrete_constraints[i].first] = 1; break;
+        case Binary : vtypes[_discrete_constraints[i].first] = 2; break;
+        default     : break;
+      }
+
+    // GUROBI variables
+    std::vector<GRBVar> vars;
+    // first all
+    for( int i=0; i<_problem->n_unknowns(); ++i)
+      switch(vtypes[i])
+      {
+        case 0 : vars.push_back( model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_CONTINUOUS) ); break;
+        case 1 : vars.push_back( model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_INTEGER   ) ); break;
+        case 2 : vars.push_back( model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_BINARY    ) ); break;
+      }
+
+
+    // Integrate new variables
+    model.update();
+
+    //----------------------------------------------
+    // 2. setup constraints
+    //----------------------------------------------
+
+    // get zero vector
+    std::vector<double> x(_problem->n_unknowns(), 0.0);
+
+    for(unsigned int i=0; i<_constraints.size();  ++i)
+    {
+      if(!_constraints[i]->is_linear())
+        std::cerr << "Warning: GUROBISolver received a problem with non-linear constraints!!!" << std::endl;
+
+      GRBLinExpr lin_expr;
+      NConstraintInterface::SVectorNC gc;
+      _constraints[i]->eval_gradient(P(x), gc);
+
+      NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
+      for(; v_it; ++v_it)
+//        lin_expr += v_it.value()*vars[v_it.index()];
+        lin_expr = lin_expr + vars[v_it.index()]*v_it.value();
+
+      double b = _constraints[i]->eval_constraint(P(x));
+
+      switch(_constraints[i]->constraint_type())
+      {
+        case NConstraintInterface::NC_EQUAL         : model.addConstr(lin_expr + b == 0); break;
+        case NConstraintInterface::NC_LESS_EQUAL    : model.addConstr(lin_expr + b <= 0); break;
+        case NConstraintInterface::NC_GREATER_EQUAL : model.addConstr(lin_expr + b >= 0); break;
+      }
+    }
+    model.update();
+
+    //----------------------------------------------
+    // 3. setup energy
+    //----------------------------------------------
+
+    if(!_problem->constant_hessian())
+      std::cerr << "Warning: GUROBISolver received a problem with non-constant hessian!!!" << std::endl;
+
+    GRBQuadExpr objective;
+
+    // add quadratic part
+    NProblemInterface::SMatrixNP H;
+    _problem->eval_hessian(P(x), H);
+    for( int i=0; i<H.outerSize(); ++i)
+      for (NProblemInterface::SMatrixNP::InnerIterator it(H,i); it; ++it)
+        objective += 0.5*it.value()*vars[it.row()]*vars[it.col()];
+
+
+    // add linear part
+    std::vector<double> g(_problem->n_unknowns());
+    _problem->eval_gradient(P(x), P(g));
+    for(unsigned int i=0; i<g.size(); ++i)
+      objective += g[i]*vars[i];
+
+    // add constant part
+    objective += _problem->eval_f(P(x));
+
+    model.set(GRB_IntAttr_ModelSense, 1);
+    model.setObjective(objective);
+
+    //----------------------------------------------
+    // 4. solve problem
+    //----------------------------------------------
+
+    model.optimize();
+
+    //----------------------------------------------
+    // 5. store result
+    //----------------------------------------------
+
+    for(unsigned int i=0; i<vars.size(); ++i)
+      x[i] = vars[i].get(GRB_DoubleAttr_X);
+
+    _problem->store_result(P(x));
+
+    std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
+    return true;
+  }
+  catch(GRBException e)
+  {
+    std::cout << "Error code = " << e.getErrorCode() << std::endl;
+    std::cout << e.getMessage() << std::endl;
+    return false;
+  }
+  catch(...)
+  {
+    std::cout << "Exception during optimization" << std::endl;
+    return false;
+  }
+
+  return false;
+}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // COMISO_GUROBI_AVAILABLE
+//=============================================================================
diff --git a/CoMISo/NSolver/GUROBISolver.hh b/CoMISo/NSolver/GUROBISolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d07dce6747ac44bf63a0f094e8d3a828fbfdd340
--- /dev/null
+++ b/CoMISo/NSolver/GUROBISolver.hh
@@ -0,0 +1,83 @@
+//=============================================================================
+//
+//  CLASS GUROBISolver
+//
+//=============================================================================
+
+
+#ifndef COMISO_GUROBISOLVER_HH
+#define COMISO_GUROBISOLVER_HH
+
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include <CoMISo/Config/config.hh>
+#if COMISO_GUROBI_AVAILABLE
+
+//== INCLUDES =================================================================
+
+#include <CoMISo/Config/CoMISoDefines.hh>
+#include <vector>
+#include "NProblemInterface.hh"
+#include "NConstraintInterface.hh"
+#include "VariableType.hh"
+
+#include <gurobi_c++.h>
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+	      
+
+/** \class NewtonSolver GUROBISolver.hh
+
+    Brief Description.
+  
+    A more elaborate description follows.
+*/
+class COMISODLLEXPORT GUROBISolver
+{
+public:
+   
+  typedef std::pair<unsigned int, VariableType> PairUiV;
+
+  /// Default constructor -> set up IpOptApplication
+  GUROBISolver() {}
+ 
+  /// Destructor
+  ~GUROBISolver() {}
+
+  // ********** SOLVE **************** //
+  bool solve(NProblemInterface*                  _problem,                // problem instance
+             std::vector<NConstraintInterface*>& _constraints,            // linear constraints
+             std::vector<PairUiV>&               _discrete_constraints,   // discrete constraints
+             const double                        _time_limit = 60     ); // time limit in seconds
+
+protected:
+  double* P(std::vector<double>& _v)
+  {
+    if( !_v.empty())
+      return ((double*)&_v[0]);
+    else
+      return 0;
+  }
+
+private:
+
+};
+
+
+
+//=============================================================================
+} // namespace COMISO
+
+//=============================================================================
+#endif // COMISO_GUROBI_AVAILABLE
+//=============================================================================
+#endif // ACG_GUROBISOLVER_HH defined
+//=============================================================================
+
diff --git a/CoMISo/NSolver/LinearConstraint.hh b/CoMISo/NSolver/LinearConstraint.hh
index 7d03f6b4aaca1398682dcde8c6dd09acba5a750a..d3bb143b057b2ded14dc70fdf3c144c7dec0db43 100644
--- a/CoMISo/NSolver/LinearConstraint.hh
+++ b/CoMISo/NSolver/LinearConstraint.hh
@@ -65,6 +65,8 @@ public:
 
   virtual void eval_hessian    ( const double* _x, SMatrixNC& _h      );
 
+  virtual bool is_linear() { return true;}
+
   // inherited from base
 //  virtual ConstraintType  constraint_type (                                      ) { return type_; }
 
diff --git a/CoMISo/NSolver/NConstraintInterface.hh b/CoMISo/NSolver/NConstraintInterface.hh
index 3e0f115c4f7d6ad90faa89bc763cdb1469d5c2a3..5a557af3d3cfecec491d7b9bc8b7681775c6f250 100644
--- a/CoMISo/NSolver/NConstraintInterface.hh
+++ b/CoMISo/NSolver/NConstraintInterface.hh
@@ -67,6 +67,8 @@ public:
     return false;
   }
 
+  virtual bool   is_linear() { return false;}
+
   virtual double gradient_update_factor( const double* _x, double _eps = 1e-6 )
   {
     double val = eval_constraint(_x);
diff --git a/CoMISo/NSolver/VariableType.hh b/CoMISo/NSolver/VariableType.hh
new file mode 100644
index 0000000000000000000000000000000000000000..6fcff566d0c5ec9f216ea5a423d7425d0ecaa217
--- /dev/null
+++ b/CoMISo/NSolver/VariableType.hh
@@ -0,0 +1,32 @@
+//=============================================================================
+//
+//  ENUM VariableType
+//
+//=============================================================================
+
+
+#ifndef COMISO_VARIABLETYPE_HH
+#define COMISO_VARIABLETYPE_HH
+
+
+//== INCLUDES =================================================================
+
+#include <CoMISo/Config/CoMISoDefines.hh>
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO { 
+
+//== CLASS DEFINITION =========================================================
+
+enum VariableType { Real, Integer, Binary};
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+#endif // COMISO_VARIABLETYPE_HH defined
+//=============================================================================
+