diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl
index bd86ef9..4d4b623 100644
--- a/src/moi_nlp_model.jl
+++ b/src/moi_nlp_model.jl
@@ -4,7 +4,9 @@ mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}}
   meta::NLPModelMeta{Float64, Vector{Float64}}
   eval::MOI.Nonlinear.Evaluator
   lincon::LinearConstraints
+  quadcon::QuadraticConstraints
   nlcon::NonLinearStructure
+  λ::Vector{Float64}
   obj::Objective
   counters::Counters
 end
@@ -27,10 +29,11 @@ end
 
 function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = "Generic")
   index_map, nvar, lvar, uvar, x0 = parser_variables(moimodel)
-  nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel, index_map)
+  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, index_map, nvar)
 
   nlp_data = _nlp_block(moimodel)
   nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian)
+  λ = zeros(nnln)  # Lagrange multipliers for hess_coord! and hprod! without y
 
   if nlp_data.has_objective
     obj = Objective("NONLINEAR", 0.0, spzeros(Float64, nvar), COO(), 0)
@@ -38,11 +41,11 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     obj = parser_objective_MOI(moimodel, nvar, index_map)
   end
 
-  ncon = nlin + nnln
-  lcon = vcat(lin_lcon, nl_lcon)
-  ucon = vcat(lin_ucon, nl_ucon)
-  nnzj = lincon.nnzj + nlcon.nnzj
-  nnzh = obj.nnzh + nlcon.nnzh
+  ncon = nlin + quadcon.nquad + nnln
+  lcon = vcat(lin_lcon, quad_lcon, nl_lcon)
+  ucon = vcat(lin_ucon, quad_ucon, nl_ucon)
+  nnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj
+  nnzh = obj.nnzh + quadcon.nnzh + nlcon.nnzh
 
   meta = NLPModelMeta(
     nvar,
@@ -57,13 +60,13 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     nnzh = nnzh,
     lin = collect(1:nlin),
     lin_nnzj = lincon.nnzj,
-    nln_nnzj = nlcon.nnzj,
+    nln_nnzj = quadcon.nnzj + nlcon.nnzj,
     minimize = MOI.get(moimodel, MOI.ObjectiveSense()) == MOI.MIN_SENSE,
-    islp = (obj.type == "LINEAR") && (nnln == 0),
+    islp = (obj.type == "LINEAR") && (nnln == 0) && (quadcon.nquad == 0),
     name = name,
   )
 
-  return MathOptNLPModel(meta, nlp_data.evaluator, lincon, nlcon, obj, Counters()), index_map
+  return MathOptNLPModel(meta, nlp_data.evaluator, lincon, quadcon, nlcon, λ, obj, Counters()), index_map
 end
 
 function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector)
@@ -106,7 +109,13 @@ end
 
 function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVector)
   increment!(nlp, :neval_cons_nln)
-  MOI.eval_constraint(nlp.eval, c, x)
+  for i = 1:(nlp.quadcon.nquad)
+    qcon = nlp.quadcon.constraints[i]
+    c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
+  end
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    MOI.eval_constraint(nlp.eval, view(c, (nlp.quadcon.nquad + 1):(nlp.meta.nnln)), x)
+  end
   return c
 end
 
@@ -125,8 +134,20 @@ function NLPModels.jac_nln_structure!(
   rows::AbstractVector{<:Integer},
   cols::AbstractVector{<:Integer},
 )
-  view(rows, 1:(nlp.nlcon.nnzj)) .= nlp.nlcon.jac_rows
-  view(cols, 1:(nlp.nlcon.nnzj)) .= nlp.nlcon.jac_cols
+  if nlp.quadcon.nquad > 0
+    index = 0
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      view(rows, index+1:index+qcon.nnzg) .= i
+      view(cols, index+1:index+qcon.nnzg) .= qcon.g
+      index += qcon.nnzg
+    end
+  end
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj)
+    view(rows, ind_nnln) .= nlp.quadcon.nquad .+ nlp.nlcon.jac_rows
+    view(cols, ind_nnln) .= nlp.nlcon.jac_cols
+  end
   return rows, cols
 end
 
@@ -138,7 +159,33 @@ end
 
 function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::AbstractVector)
   increment!(nlp, :neval_jac_nln)
-  MOI.eval_constraint_jacobian(nlp.eval, view(vals, 1:(nlp.nlcon.nnzj)), x)
+  if nlp.quadcon.nquad > 0
+    index = 0
+    view(vals, 1:nlp.quadcon.nnzj) .= 0.0
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      for (j, ind) in enumerate(qcon.b.nzind)
+        k = qcon.dg[ind]
+        vals[index+k] += qcon.b.nzval[j]
+      end
+      for j = 1:qcon.nnzh
+        row = qcon.A.rows[j]
+        col = qcon.A.cols[j]
+        val = qcon.A.vals[j]
+        k1 = qcon.dg[row]
+        vals[index+k1] += val * x[col]
+        if row != col
+          k2 = qcon.dg[col]
+          vals[index+k2] += val * x[row]
+        end
+      end
+      index += qcon.nnzg
+    end
+  end
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    ind_nnln = (nlp.quadcon.nnzj + 1):(nlp.quadcon.nnzj + nlp.nlcon.nnzj)
+    MOI.eval_constraint_jacobian(nlp.eval, view(vals, ind_nnln), x)
+  end
   return vals
 end
 
@@ -166,7 +213,18 @@ function NLPModels.jprod_nln!(
   Jv::AbstractVector,
 )
   increment!(nlp, :neval_jprod_nln)
-  MOI.eval_constraint_jacobian_product(nlp.eval, Jv, x, v)
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
+    MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v)
+  end
+  (nlp.meta.nnln == nlp.quadcon.nquad) && (Jv .= 0.0)
+  if nlp.quadcon.nquad > 0
+    for i = 1:(nlp.quadcon.nquad)
+      # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
+      qcon = nlp.quadcon.constraints[i]
+      Jv[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
+    end
+  end
   return Jv
 end
 
@@ -194,7 +252,19 @@ function NLPModels.jtprod_nln!(
   Jtv::AbstractVector,
 )
   increment!(nlp, :neval_jtprod_nln)
-  MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, v)
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
+    MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln))
+  end
+  (nlp.meta.nnln == nlp.quadcon.nquad) && (Jtv .= 0.0)
+  if nlp.quadcon.nquad > 0
+    for i = 1:(nlp.quadcon.nquad)
+      # Jtv += v[i] * (Aᵢ * x + bᵢ)
+      qcon = nlp.quadcon.constraints[i]
+      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i])
+      Jtv .+= v[i] .* qcon.b
+    end
+  end
   return Jtv
 end
 
@@ -207,9 +277,18 @@ function NLPModels.hess_structure!(
     view(rows, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.rows
     view(cols, 1:(nlp.obj.nnzh)) .= nlp.obj.hessian.cols
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
-    view(rows, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows
-    view(cols, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols
+  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
+    view(rows, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_rows
+    view(cols, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)) .= nlp.nlcon.hess_cols
+  end
+  if nlp.quadcon.nquad > 0
+    index = nlp.obj.nnzh
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows
+      view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols
+      index += qcon.nnzh
+    end
   end
   return rows, cols
 end
@@ -225,15 +304,23 @@ function NLPModels.hess_coord!(
   if nlp.obj.type == "QUADRATIC"
     view(vals, 1:(nlp.obj.nnzh)) .= obj_weight .* nlp.obj.hessian.vals
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
-    MOI.eval_hessian_lagrangian(
-      nlp.eval,
-      view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)),
+  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
+    λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
+    MOI.eval_hessian_lagrangian(nlp.eval,
+      view(vals, (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)),
       x,
       obj_weight,
-      view(y, nlp.meta.nln),
+      λ
     )
   end
+  if nlp.quadcon.nquad > 0
+    index = nlp.obj.nnzh
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      view(vals, (index + 1):(index + qcon.nnzh)) .= y[nlp.meta.nlin + i] .* qcon.A.vals
+      index += qcon.nnzh
+    end
+  end
   return vals
 end
 
@@ -252,7 +339,9 @@ function NLPModels.hess_coord!(
     view(vals, (nlp.obj.nnzh + 1):(nlp.meta.nnzh)) .= 0.0
   end
   if nlp.obj.type == "NONLINEAR"
-    MOI.eval_hessian_lagrangian(nlp.eval, vals, x, obj_weight, zeros(nlp.meta.nnln))
+    view(vals, 1:(nlp.obj.nnzh + nlp.quadcon.nnzh)) .= 0.0
+    ind_nnln = (nlp.obj.nnzh + nlp.quadcon.nnzh + 1):(nlp.meta.nnzh)
+    MOI.eval_hessian_lagrangian(nlp.eval, view(vals, ind_nnln), x, obj_weight, nlp.λ)
   end
   return vals
 end
@@ -269,19 +358,20 @@ function NLPModels.hprod!(
   if (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == 0)
     hv .= 0.0
   end
-  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > 0)
-    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, view(y, nlp.meta.nln))
+  if (nlp.obj.type == "NONLINEAR") || (nlp.meta.nnln > nlp.quadcon.nquad)
+    λ = view(y, (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon))
+    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, λ)
   end
   if nlp.obj.type == "QUADRATIC"
-    nlp.meta.nnln == 0 && (hv .= 0.0)
-    coo_sym_add_mul!(
-      nlp.obj.hessian.rows,
-      nlp.obj.hessian.cols,
-      nlp.obj.hessian.vals,
-      v,
-      hv,
-      obj_weight,
-    )
+    (nlp.meta.nnln == nlp.quadcon.nquad) && (hv .= 0.0)
+    coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight)
+  end
+  if nlp.quadcon.nquad > 0
+    (nlp.obj.type == "LINEAR") && (hv .= 0.0)
+    for i = 1:(nlp.quadcon.nquad)
+      qcon = nlp.quadcon.constraints[i]
+      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nlp.meta.nlin + i])
+    end
   end
   return hv
 end
@@ -302,7 +392,7 @@ function NLPModels.hprod!(
     coo_sym_add_mul!(nlp.obj.hessian.rows, nlp.obj.hessian.cols, nlp.obj.hessian.vals, v, hv, obj_weight)
   end
   if nlp.obj.type == "NONLINEAR"
-    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, zeros(nlp.meta.nnln))
+    MOI.eval_hessian_lagrangian_product(nlp.eval, hv, x, v, obj_weight, nlp.λ)
   end
   return hv
 end
diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl
index 0c1cd8d..c2f7650 100644
--- a/src/moi_nls_model.jl
+++ b/src/moi_nls_model.jl
@@ -9,7 +9,9 @@ mutable struct MathOptNLSModel <: AbstractNLSModel{Float64, Vector{Float64}}
   linequ::LinearEquations
   nlequ::NonLinearStructure
   lincon::LinearConstraints
+  quadcon::QuadraticConstraints
   nlcon::NonLinearStructure
+  λ::Vector{Float64}
   counters::NLSCounters
 end
 
@@ -30,20 +32,21 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
 
   _nlp_sync!(cmodel)
   moimodel = backend(cmodel)
-  nlin, lincon, lin_lcon, lin_ucon = parser_MOI(moimodel, index_map)
+  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon = parser_MOI(moimodel, index_map, nvar)
 
   nlp_data = _nlp_block(moimodel)
   nnln, nlcon, nl_lcon, nl_ucon = parser_NL(nlp_data, hessian = hessian)
+  λ = zeros(nnln)  # Lagrange multipliers for hess_coord! and hprod! without y
 
   nequ = nlinequ + nnlnequ
   Fnnzj = linequ.nnzj + nlequ.nnzj
   Fnnzh = nlequ.nnzh
 
-  ncon = nlin + nnln
-  lcon = vcat(lin_lcon, nl_lcon)
-  ucon = vcat(lin_ucon, nl_ucon)
-  cnnzj = lincon.nnzj + nlcon.nnzj
-  cnnzh = lls.nnzh + nlcon.nnzh
+  ncon = nlin + quadcon.nquad + nnln
+  lcon = vcat(lin_lcon, quad_lcon, nl_lcon)
+  ucon = vcat(lin_ucon, quad_ucon, nl_ucon)
+  cnnzj = lincon.nnzj + quadcon.nnzj + nlcon.nnzj
+  cnnzh = lls.nnzh + quadcon.nnzh + nlcon.nnzh
 
   meta = NLPModelMeta(
     nvar,
@@ -58,7 +61,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
     nnzh = cnnzh,
     lin = collect(1:nlin),
     lin_nnzj = lincon.nnzj,
-    nln_nnzj = nlcon.nnzj,
+    nln_nnzj = quadcon.nnzj + nlcon.nnzj,
     minimize = objective_sense(cmodel) == MOI.MIN_SENSE,
     islp = false,
     name = name,
@@ -73,7 +76,9 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
     linequ,
     nlequ,
     lincon,
+    quadcon,
     nlcon,
+    λ,
     NLSCounters(),
   )
 end
@@ -253,7 +258,13 @@ end
 
 function NLPModels.cons_nln!(nls::MathOptNLSModel, x::AbstractVector, c::AbstractVector)
   increment!(nls, :neval_cons_nln)
-  MOI.eval_constraint(nls.ceval, c, x)
+  for i = 1:(nls.quadcon.nquad)
+    qcon = nls.quadcon.constraints[i]
+    c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
+  end
+  if nls.meta.nnln > nls.quadcon.nquad
+    MOI.eval_constraint(nls.ceval, view(c, (nls.quadcon.nquad + 1):(nls.meta.nnln)), x)
+  end
   return c
 end
 
@@ -272,8 +283,20 @@ function NLPModels.jac_nln_structure!(
   rows::AbstractVector{<:Integer},
   cols::AbstractVector{<:Integer},
 )
-  view(rows, 1:(nls.nlcon.nnzj)) .= nls.nlcon.jac_rows
-  view(cols, 1:(nls.nlcon.nnzj)) .= nls.nlcon.jac_cols
+  if nls.quadcon.nquad > 0
+    index = 0
+    for i = 1:(nls.quadcon.nquad)
+      qcon = nls.quadcon.constraints[i]
+      view(rows, index+1:index+qcon.nnzg) .= i
+      view(cols, index+1:index+qcon.nnzg) .= qcon.g
+      index += qcon.nnzg
+    end
+  end
+  if nls.meta.nnln > nls.quadcon.nquad
+    ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj)
+    view(rows, ind_nnln) .= nls.quadcon.nquad .+ nls.nlcon.jac_rows
+    view(cols, ind_nnln) .= nls.nlcon.jac_cols
+  end
   return rows, cols
 end
 
@@ -285,7 +308,33 @@ end
 
 function NLPModels.jac_nln_coord!(nls::MathOptNLSModel, x::AbstractVector, vals::AbstractVector)
   increment!(nls, :neval_jac_nln)
-  MOI.eval_constraint_jacobian(nls.ceval, view(vals, 1:(nls.nlcon.nnzj)), x)
+  if nls.quadcon.nquad > 0
+    index = 0
+    view(vals, 1:nls.quadcon.nnzj) .= 0.0
+    for i = 1:(nls.quadcon.nquad)
+      qcon = nls.quadcon.constraints[i]
+      for (j, ind) in enumerate(qcon.b.nzind)
+        k = qcon.dg[ind]
+        vals[index+k] += qcon.b.nzval[j]
+      end
+      for j = 1:qcon.nnzh
+        row = qcon.A.rows[j]
+        col = qcon.A.cols[j]
+        val = qcon.A.vals[j]
+        k1 = qcon.dg[row]
+        vals[index+k1] += val * x[col]
+        if row != col
+          k2 = qcon.dg[col]
+          vals[index+k2] += val * x[row]
+        end
+      end
+      index += qcon.nnzg
+    end
+  end
+  if nls.meta.nnln > nls.quadcon.nquad
+    ind_nnln = (nls.quadcon.nnzj + 1):(nls.quadcon.nnzj + nls.nlcon.nnzj)
+    MOI.eval_constraint_jacobian(nls.ceval, view(vals, ind_nnln), x)
+  end
   return vals
 end
 
@@ -313,7 +362,18 @@ function NLPModels.jprod_nln!(
   Jv::AbstractVector,
 )
   increment!(nls, :neval_jprod_nln)
-  MOI.eval_constraint_jacobian_product(nls.ceval, Jv, x, v)
+  if nls.meta.nnln > nls.quadcon.nquad
+    ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln)
+    MOI.eval_constraint_jacobian_product(nls.ceval, view(Jv, ind_nnln), x, v)
+  end
+  (nls.meta.nnln == nls.quadcon.nquad) && (Jv .= 0.0)
+  if nls.quadcon.nquad > 0
+    for i = 1:(nls.quadcon.nquad)
+      # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
+      qcon = nls.quadcon.constraints[i]
+      Jv[i] += coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
+    end
+  end
   return Jv
 end
 
@@ -341,7 +401,19 @@ function NLPModels.jtprod_nln!(
   Jtv::AbstractVector,
 )
   increment!(nls, :neval_jtprod_nln)
-  MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, v)
+  if nls.meta.nnln > nls.quadcon.nquad
+    ind_nnln = (nls.quadcon.nquad + 1):(nls.meta.nnln)
+    MOI.eval_constraint_jacobian_transpose_product(nls.ceval, Jtv, x, view(v, ind_nnln))
+  end
+  (nls.meta.nnln == nls.quadcon.nquad) && (Jtv .= 0.0)
+  if nls.quadcon.nquad > 0
+    for i = 1:(nls.quadcon.nquad)
+      # Jtv += v[i] * (Aᵢ * x + bᵢ)
+      qcon = nls.quadcon.constraints[i]
+      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i])
+      Jtv .+= v[i] .* qcon.b
+    end
+  end
   return Jtv
 end
 
@@ -354,9 +426,18 @@ function NLPModels.hess_structure!(
     view(rows, 1:(nls.lls.nnzh)) .= nls.lls.hessian.rows
     view(cols, 1:(nls.lls.nnzh)) .= nls.lls.hessian.cols
   end
-  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > 0)
-    view(rows, (nls.lls.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows
-    view(cols, (nls.lls.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols
+  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad)
+    view(rows, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_rows
+    view(cols, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= nls.nlcon.hess_cols
+  end
+  if nls.quadcon.nquad > 0
+    index = nls.lls.nnzh
+    for i = 1:(nls.quadcon.nquad)
+      qcon = nls.quadcon.constraints[i]
+      view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows
+      view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols
+      index += qcon.nnzh
+    end
   end
   return rows, cols
 end
@@ -372,15 +453,24 @@ function NLPModels.hess_coord!(
   if nls.nls_meta.nlin > 0
     view(vals, 1:(nls.lls.nnzh)) .= obj_weight .* nls.lls.hessian.vals
   end
-  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > 0)
+  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad)
+    λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon))
     MOI.eval_hessian_lagrangian(
       nls.ceval,
-      view(vals, (nls.lls.nnzh + 1):(nls.meta.nnzh)),
+      view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)),
       x,
       obj_weight,
-      view(y, nls.meta.nln),
+      λ
     )
   end
+  if nls.quadcon.nquad > 0
+    index = nls.lls.nnzh
+    for i = 1:(nls.quadcon.nquad)
+      qcon = nls.quadcon.constraints[i]
+      view(vals, (index + 1):(index + qcon.nnzh)) .= y[nls.meta.nlin + i] .* qcon.A.vals
+      index += qcon.nnzh
+    end
+  end
   return vals
 end
 
@@ -394,16 +484,14 @@ function NLPModels.hess_coord!(
   if nls.nls_meta.nlin > 0
     view(vals, 1:(nls.lls.nnzh)) .= obj_weight .* nls.lls.hessian.vals
   end
+  view(vals, (nls.lls.nnzh + 1):(nls.lls.nnzh + nls.quadcon.nnzh)) .= 0.0
   if nls.nls_meta.nnln > 0
-    MOI.eval_hessian_lagrangian(
-      nls.ceval,
-      view(vals, (nls.lls.nnzh + 1):(nls.meta.nnzh)),
-      x,
-      obj_weight,
-      zeros(nls.meta.nnln),
-    )
+    ind_nnln = (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)
+    MOI.eval_hessian_lagrangian(nls.ceval, view(vals, ind_nnln), x, obj_weight, nls.λ)
   else
-    view(vals, (nls.lls.nnzh + 1):(nls.meta.nnzh)) .= 0.0
+    if nls.meta.nnln > nls.quadcon.nquad
+      view(vals, (nls.lls.nnzh + nls.quadcon.nnzh + 1):(nls.meta.nnzh)) .= 0.0
+    end
   end
   return vals
 end
@@ -417,11 +505,12 @@ function NLPModels.hprod!(
   obj_weight::Float64 = 1.0,
 )
   increment!(nls, :neval_hprod)
-  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > 0)
-    MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, view(y, nls.meta.nln))
+  if (nls.nls_meta.nnln > 0) || (nls.meta.nnln > nls.quadcon.nquad)
+    λ = view(y, (nls.meta.nlin + nls.quadcon.nquad + 1):(nls.meta.ncon))
+    MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, λ)
   end
+  (nls.nls_meta.nnln == 0) && (nls.meta.nnln == nls.quadcon.nquad) && (hv .= 0.0)
   if nls.nls_meta.nlin > 0
-    (nls.nls_meta.nnln == 0) && (nls.meta.nnln == 0) && (hv .= 0.0)
     coo_sym_add_mul!(
       nls.lls.hessian.rows,
       nls.lls.hessian.cols,
@@ -431,6 +520,12 @@ function NLPModels.hprod!(
       obj_weight,
     )
   end
+  if nls.quadcon.nquad > 0
+    for i = 1:(nls.quadcon.nquad)
+      qcon = nls.quadcon.constraints[i]
+      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nls.meta.nlin + i])
+    end
+  end
   return hv
 end
 
@@ -443,7 +538,7 @@ function NLPModels.hprod!(
 )
   increment!(nls, :neval_hprod)
   if nls.nls_meta.nnln > 0
-    MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, zeros(nls.meta.nnln))
+    MOI.eval_hessian_lagrangian_product(nls.ceval, hv, x, v, obj_weight, nls.λ)
   end
   if nls.nls_meta.nlin > 0
     (nls.nls_meta.nnln == 0) && (hv .= 0.0)
diff --git a/src/utils.jl b/src/utils.jl
index e1cf1c6..35d3b33 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -6,10 +6,18 @@ import NLPModels.increment!, NLPModels.decrement!
 using JuMP, MathOptInterface
 const MOI = MathOptInterface
 
+# VariableIndex
+const VI = MOI.VariableIndex  # VariableIndex(value)
+
 # ScalarAffineFunctions and VectorAffineFunctions
-const SAF = MOI.ScalarAffineFunction{Float64}
-const VAF = MOI.VectorAffineFunction{Float64}
-const AF = Union{SAF, VAF}
+const SAF = MOI.ScalarAffineFunction{Float64}  # ScalarAffineFunction{T}(terms, constant)
+const VAF = MOI.VectorAffineFunction{Float64}  # VectorAffineFunction{T}(terms, constants)
+const AF  = Union{SAF, VAF}
+
+# ScalarQuadraticFunctions and VectorQuadraticFunctions
+const SQF = MOI.ScalarQuadraticFunction{Float64}  # ScalarQuadraticFunction{T}(affine_terms, quadratic_terms, constant)
+const VQF = MOI.VectorQuadraticFunction{Float64}  # VectorQuadraticFunction{T}(affine_terms, quadratic_terms, constants)
+const QF  = Union{SQF, VQF}
 
 # AffLinSets and VecLinSets
 const ALS = Union{
@@ -21,10 +29,6 @@ const ALS = Union{
 const VLS = Union{MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros}
 const LS = Union{ALS, VLS}
 
-const VI = MOI.VariableIndex
-const SQF = MOI.ScalarQuadraticFunction{Float64}
-const LinQuad = Union{VI, SAF, SQF}
-
 # Expressions
 const VF = VariableRef
 const AE = GenericAffExpr{Float64, VariableRef}
@@ -32,6 +36,9 @@ const LE = Union{VF, AE}
 const QE = GenericQuadExpr{Float64, VariableRef}
 const NLE = NonlinearExpression
 
+const LinQuad = Union{VI, SAF, SQF}
+
+# Sparse matrix in coordinate format
 mutable struct COO
   rows::Vector{Int}
   cols::Vector{Int}
@@ -45,6 +52,23 @@ mutable struct LinearConstraints
   nnzj::Int
 end
 
+# xᵀAx + bᵀx
+mutable struct QuadraticConstraint
+  A::COO
+  b::SparseVector{Float64}
+  g::Vector{Int}
+  dg::Dict{Int,Int}
+  nnzg::Int
+  nnzh::Int
+end
+
+mutable struct QuadraticConstraints
+  nquad::Int
+  constraints::Vector{QuadraticConstraint}
+  nnzj::Int
+  nnzh::Int
+end
+
 mutable struct NonLinearStructure
   jac_rows::Vector{Int}
   jac_cols::Vector{Int}
@@ -200,11 +224,140 @@ function parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_uco
 end
 
 """
-    parser_MOI(moimodel, index_map)
+    parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+
+Parse a `ScalarQuadraticFunction` fun with its associated set.
+`qcons`, `quad_lcon`, `quad_ucon` are updated.
+"""
+function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+  _index(v::MOI.VariableIndex) = index_map[v].value
+
+  b = spzeros(Float64, nvar)
+  rows = Int[]
+  cols = Int[]
+  vals = Float64[]
+
+  # Parse a ScalarAffineTerm{Float64}(coefficient, variable_index)
+  for term in fun.affine_terms
+    b[_index(term.variable)] = term.coefficient
+  end
+
+  # Parse a ScalarQuadraticTerm{Float64}(coefficient, variable_index_1, variable_index_2)
+  for term in fun.quadratic_terms
+    i = _index(term.variable_1)
+    j = _index(term.variable_2)
+    if i ≥ j
+      push!(rows, i)
+      push!(cols, j)
+    else
+      push!(rows, j)
+      push!(cols, i)
+    end
+    push!(vals, term.coefficient)
+  end
+
+  if typeof(set) in (MOI.Interval{Float64}, MOI.GreaterThan{Float64})
+    push!(quad_lcon, -fun.constant + set.lower)
+  elseif typeof(set) == MOI.EqualTo{Float64}
+    push!(quad_lcon, -fun.constant + set.value)
+  else
+    push!(quad_lcon, -Inf)
+  end
+
+  if typeof(set) in (MOI.Interval{Float64}, MOI.LessThan{Float64})
+    push!(quad_ucon, -fun.constant + set.upper)
+  elseif typeof(set) == MOI.EqualTo{Float64}
+    push!(quad_ucon, -fun.constant + set.value)
+  else
+    push!(quad_ucon, Inf)
+  end
+
+  A = COO(rows, cols, vals)
+  g = unique(vcat(rows, cols, b.nzind))  # sparsity pattern of Ax + b
+  nnzg = length(g)
+  # dg is a dictionary where:
+  # - The key `r` specifies a row index in the vector Ax + b.
+  # - The value `dg[r]` is a position in the vector (of length nnzg)
+  # where the non-zero entries of the Jacobian for row `r` are stored.
+  dg = Dict{Int,Int}(g[p] => p for p = 1:nnzg)
+  nnzh = length(vals)
+  qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh)
+  push!(qcons, qcon)
+end
+
+"""
+    parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+
+Parse a `VectorQuadraticFunction` fun with its associated set.
+`qcons`, `quad_lcon`, `quad_ucon` are updated.
+"""
+function parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+  _index(v::MOI.VariableIndex) = index_map[v].value
+
+  ncon = length(fun.constants)
+  for k = 1:ncon
+    b = spzeros(Float64, nvar)
+    rows = Int[]
+    cols = Int[]
+    vals = Float64[]
+
+    # Parse a VectorAffineTerm{Float64}(output_index, scalar_term)
+    for affine_term in fun.affine_terms
+      if affine_term.output_index == k
+        b[_index(affine_term.scalar_term.variable)] = affine_term.scalar_term.coefficient
+      end
+    end
+
+    # Parse a VectorQuadraticTerm{Float64}(output_index, scalar_term)
+    for quadratic_term in fun.quadratic_terms
+      if quadratic_term.output_index == k
+        i = _index(quadratic_term.scalar_term.variable_1)
+        j = _index(quadratic_term.scalar_term.variable_2)
+        if i ≥ j
+          push!(rows, i)
+          push!(cols, j)
+        else
+          push!(rows, j)
+          push!(cols, i)
+        end
+        push!(vals, quadratic_term.scalar_term.coefficient)
+      end
+    end
+
+    constant = fun.constants[k]
+
+    if typeof(set) in (MOI.Nonnegatives, MOI.Zeros)
+      append!(quad_lcon, constant)
+    else
+      append!(quad_lcon, -Inf)
+    end
+
+    if typeof(set) in (MOI.Nonpositives, MOI.Zeros)
+      append!(quad_ucon, -constant)
+    else
+      append!(quad_ucon, Inf)
+    end
+
+    A = COO(rows, cols, vals)
+    g = unique(vcat(rows, cols, b.nzind))  # sparsity pattern of Ax + b
+    nnzg = length(g)
+    # dg is a dictionary where:
+    # - The key `r` specifies a row index in the vector Ax + b.
+    # - The value `dg[r]` is a position in the vector (of length nnzg)
+    # where the non-zero entries of the Jacobian for row `r` are stored.
+    dg = Dict{Int,Int}(g[p] => p for p = 1:nnzg)
+    nnzh = length(vals)
+    qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh)
+    push!(qcons, qcon)
+  end
+end
+
+"""
+    parser_MOI(moimodel, index_map, nvar)
 
 Parse linear constraints of a `MOI.ModelLike`.
 """
-function parser_MOI(moimodel, index_map)
+function parser_MOI(moimodel, index_map, nvar)
 
   # Variables associated to linear constraints
   nlin = 0
@@ -214,10 +367,16 @@ function parser_MOI(moimodel, index_map)
   lin_lcon = Float64[]
   lin_ucon = Float64[]
 
+  # Variables associated to quadratic constraints
+  nquad = 0
+  qcons = QuadraticConstraint[]
+  quad_lcon = Float64[]
+  quad_ucon = Float64[]
+
   contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent())
   for (F, S) in contypes
-    F <: AF || F == MOI.ScalarNonlinearFunction || F == VI || @warn("Function $F is not supported.")
-    S <: LS || @warn("Set $S is not supported.")
+    F <: AF || F <: QF || F == MOI.ScalarNonlinearFunction || F == VI || error("Function $F is not supported.")
+    S <: LS || error("Set $S is not supported.")
 
     conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}())
     for cidx in conindices
@@ -237,13 +396,29 @@ function parser_MOI(moimodel, index_map)
         parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon, index_map)
         nlin += set.dimension
       end
+      if typeof(fun) <: SQF
+        parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+        nquad += 1
+      end
+      if typeof(fun) <: VQF
+        parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+        nquad += set.dimension
+      end
     end
   end
   coo = COO(linrows, lincols, linvals)
-  nnzj = length(linvals)
-  lincon = LinearConstraints(coo, nnzj)
+  lin_nnzj = length(linvals)
+  lincon = LinearConstraints(coo, lin_nnzj)
+
+  quad_nnzj = 0
+  quad_nnzh = 0
+  for i = 1:nquad
+    quad_nnzj += qcons[i].nnzg
+    quad_nnzh += qcons[i].nnzh
+  end
+  quadcon = QuadraticConstraints(nquad, qcons, quad_nnzj, quad_nnzh)
 
-  return nlin, lincon, lin_lcon, lin_ucon
+  return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon
 end
 
 # Affine or quadratic, nothing to do
@@ -314,9 +489,9 @@ function _nlp_block(model::MOI.ModelLike)
 end
 
 """
-    parser_NL(jmodel, moimodel)
+    parser_NL(nlp_data; hessian)
 
-Parse nonlinear constraints of a `MOI.Nonlinear.Evaluator`.
+Parse nonlinear constraints of an `nlp_data`.
 """
 function parser_NL(nlp_data; hessian::Bool = true)
   nnln = length(nlp_data.constraint_bounds)
@@ -373,7 +548,7 @@ function parser_variables(model::MOI.ModelLike)
 end
 
 """
-    parser_objective_MOI(moimodel, nvar)
+    parser_objective_MOI(moimodel, nvar, index_map)
 
 Parse linear and quadratic objective of a `MOI.ModelLike`.
 """
@@ -419,8 +594,8 @@ function parser_objective_MOI(moimodel, nvar, index_map)
         push!(rows, i)
         push!(cols, j)
       else
-        push!(cols, j)
-        push!(rows, i)
+        push!(rows, j)
+        push!(cols, i)
       end
       push!(vals, term.coefficient)
     end
@@ -429,7 +604,7 @@ function parser_objective_MOI(moimodel, nvar, index_map)
 end
 
 """
-    parser_linear_expression(cmodel, nvar, F)
+    parser_linear_expression(cmodel, nvar, index_map, F)
 
 Parse linear expressions of type `VariableRef` and `GenericAffExpr{Float64,VariableRef}`.
 """
@@ -522,7 +697,7 @@ function add_constraint_model(Fmodel, Fi::AbstractArray)
 end
 
 """
-    parser_nonlinear_expression(cmodel, nvar, F)
+    parser_nonlinear_expression(cmodel, nvar, F; hessian)
 
 Parse nonlinear expressions of type `NonlinearExpression`.
 """
diff --git a/test/nlp_consistency.jl b/test/nlp_consistency.jl
index f8108d9..b116823 100644
--- a/test/nlp_consistency.jl
+++ b/test/nlp_consistency.jl
@@ -4,7 +4,7 @@ for problem in nlp_problems
     problem_f = eval(Symbol(lowercase(problem)))
     nlp_moi = MathOptNLPModel(problem_f())
     nlps = [nlp_manual; nlp_moi]
-    consistent_nlps(nlps, linear_api = true)
+    consistent_nlps(nlps, linear_api=true, test_slack=false)
     view_subarray_nlp(nlp_moi)
   end
 end
diff --git a/test/nlp_problems/hs10.jl b/test/nlp_problems/hs10.jl
index e6d47bd..3e657e1 100644
--- a/test/nlp_problems/hs10.jl
+++ b/test/nlp_problems/hs10.jl
@@ -8,7 +8,7 @@ function hs10()
 
   @objective(nlp, Min, x[1] - x[2])
 
-  @NLconstraint(nlp, -3 * x[1]^2 + 2 * x[1] * x[2] - x[2]^2 ≥ -1)
+  @constraint(nlp, -3 * x[1]^2 + 2 * x[1] * x[2] - x[2]^2 ≥ -1)
 
   return nlp
 end
diff --git a/test/nlp_problems/hs100.jl b/test/nlp_problems/hs100.jl
new file mode 100644
index 0000000..92156cb
--- /dev/null
+++ b/test/nlp_problems/hs100.jl
@@ -0,0 +1,24 @@
+function hs100(args...; kwargs...)
+  nlp = Model()
+  x0 = [1, 2, 0, 4, 0, 1, 1]
+  @variable(nlp, x[i = 1:7], start = x0[i])
+
+  @NLconstraint(nlp, 127 - 2 * x[1]^2 - 3 * x[2]^4 - x[3] - 4 * x[4]^2 - 5 * x[5] ≥ 0)
+  @constraint(nlp, 282 - 7 * x[1] - 3 * x[2] - 10 * x[3]^2 - x[4] + x[5] ≥ 0)
+  @constraint(nlp, -196 + 23 * x[1] + x[2]^2 + 6 * x[6]^2 - 8 * x[7] ≤ 0)
+  @constraint(nlp, -4 * x[1]^2 - x[2]^2 + 3 * x[1] * x[2] - 2 * x[3]^2 - 5 * x[6] + 11 * x[7] ≥ 0)
+
+  @NLobjective(
+    nlp,
+    Min,
+    (x[1] - 10)^2 +
+    5 * (x[2] - 12)^2 +
+    x[3]^4 +
+    3 * (x[4] - 11)^2 +
+    10 * x[5]^6 +
+    7 * x[6]^2 +
+    x[7]^4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]
+  )
+
+  return nlp
+end
diff --git a/test/nlp_problems/hs11.jl b/test/nlp_problems/hs11.jl
index 2217383..fd23caf 100644
--- a/test/nlp_problems/hs11.jl
+++ b/test/nlp_problems/hs11.jl
@@ -8,7 +8,7 @@ function hs11()
 
   @objective(nlp, Min, (x[1] - 5)^2 + x[2]^2 - 25)
 
-  @NLconstraint(nlp, -x[1]^2 + x[2] ≥ 0)
+  @constraint(nlp, -x[1]^2 + x[2] ≥ 0)
 
   return nlp
 end
diff --git a/test/nlp_problems/hs14.jl b/test/nlp_problems/hs14.jl
index 01bae89..e349d98 100644
--- a/test/nlp_problems/hs14.jl
+++ b/test/nlp_problems/hs14.jl
@@ -8,7 +8,7 @@ function hs14()
 
   @objective(nlp, Min, (x[1] - 2)^2 + (x[2] - 1)^2)
 
-  @NLconstraint(nlp, -x[1]^2 / 4 - x[2]^2 + 1 ≥ 0)
+  @constraint(nlp, -x[1]^2 / 4 - x[2]^2 + 1 ≥ 0)
 
   @constraint(nlp, x[1] - 2 * x[2] + 1 == 0)
 
diff --git a/test/nlp_problems/hs219.jl b/test/nlp_problems/hs219.jl
new file mode 100644
index 0000000..ada2147
--- /dev/null
+++ b/test/nlp_problems/hs219.jl
@@ -0,0 +1,16 @@
+function hs219(args...; kwargs...)
+  nlp = Model()
+  x0 = [10, 10, 10, 10]
+  @variable(nlp, x[i = 1:4], start = x0[i])
+
+  @constraint(nlp, x[1]^2 - x[2] - x[4]^2 == 0)
+  @NLconstraint(nlp, x[2] - x[1]^3 - x[3]^2 == 0)
+
+  @NLobjective(
+    nlp,
+    Min,
+    -x[1]
+  )
+
+  return nlp
+end
diff --git a/test/nlp_problems/hs6.jl b/test/nlp_problems/hs6.jl
index dde92d9..d229738 100644
--- a/test/nlp_problems/hs6.jl
+++ b/test/nlp_problems/hs6.jl
@@ -8,7 +8,7 @@ function hs6()
 
   @NLobjective(nlp, Min, (1 - x[1])^2)
 
-  @NLconstraint(nlp, 10 * (x[2] - x[1]^2) == 0)
+  @constraint(nlp, 10 * (x[2] - x[1]^2) == 0)
 
   return nlp
 end
diff --git a/test/nlp_problems/hs61.jl b/test/nlp_problems/hs61.jl
new file mode 100644
index 0000000..5e422cc
--- /dev/null
+++ b/test/nlp_problems/hs61.jl
@@ -0,0 +1,12 @@
+"HS61 model"
+function hs61(args...; kwargs...)
+  nlp = Model()
+  @variable(nlp, x[i = 1:3], start = 0)
+
+  @constraint(nlp, 3 * x[1] - 2 * x[2]^2 - 7 == 0)
+  @constraint(nlp, 4 * x[1] - x[3]^2 - 11 == 0)
+
+  @NLobjective(nlp, Min, 4 * x[1]^2 + 2 * x[2]^2 + 2 * x[3]^2 - 33 * x[1] + 16 * x[2] - 24 * x[3])
+
+  return nlp
+end
diff --git a/test/nlp_problems/quadcon.jl b/test/nlp_problems/quadcon.jl
new file mode 100644
index 0000000..020ca1d
--- /dev/null
+++ b/test/nlp_problems/quadcon.jl
@@ -0,0 +1,13 @@
+function quadcon()
+  nlp = Model()
+
+  @variable(nlp, x[1:3])
+
+  @constraint(nlp, [x[1], x[1]^2 + x[1] * x[2] + 3.0] in MOI.Nonnegatives(2))
+  @constraint(nlp, [x[2], x[2]^2 + x[2] * x[3] + 4.0] in MOI.Nonpositives(2))
+  @constraint(nlp, [x[3], x[3]^2 + x[3] * x[1] + 5.0] in MOI.Zeros(2))
+
+  @objective(nlp, Min, x[1] * x[2] * exp(x[3]))
+
+  return nlp
+end
diff --git a/test/nls_consistency.jl b/test/nls_consistency.jl
index a7a868f..2025c52 100644
--- a/test/nls_consistency.jl
+++ b/test/nls_consistency.jl
@@ -8,7 +8,7 @@ for problem in nls_problems
     if isdefined(Main, Symbol(spc))
       push!(nlss, eval(Meta.parse(spc))())
     end
-    consistent_nlss(nlss, linear_api = true)
+    consistent_nlss(nlss, linear_api=true, test_slack=false)
     view_subarray_nls(nls_moi)
   end
 end
diff --git a/test/nls_problems/nlshs20.jl b/test/nls_problems/nlshs20.jl
index 98d500b..7b1ef17 100644
--- a/test/nls_problems/nlshs20.jl
+++ b/test/nls_problems/nlshs20.jl
@@ -6,9 +6,9 @@ function nlshs20()
   x0 = [-2.0; 1.0]
   @variable(model, lvar[i] ≤ x[i = 1:2] ≤ uvar[i], start = x0[i])
 
-  @NLconstraint(model, x[1] + x[2]^2 ≥ 0)
-  @NLconstraint(model, x[1]^2 + x[2] ≥ 0)
-  @NLconstraint(model, x[1]^2 + x[2]^2 ≥ 1)
+  @constraint(model, x[1] + x[2]^2 ≥ 0)
+  @constraint(model, x[1]^2 + x[2] ≥ 0)
+  @constraint(model, x[1]^2 + x[2]^2 ≥ 1)
 
   @expression(model, F1, 1 - x[1])
   @NLexpression(model, F2, 10 * (x[2] - x[1]^2))
diff --git a/test/nls_problems/nlsqc.jl b/test/nls_problems/nlsqc.jl
new file mode 100644
index 0000000..3c80c81
--- /dev/null
+++ b/test/nls_problems/nlsqc.jl
@@ -0,0 +1,13 @@
+function nlsqc()
+  nls = Model()
+
+  @variable(nls, x[1:3])
+
+  @constraint(nls, [x[1], x[1]^2 + x[1] * x[2] + 3.0] in MOI.Nonnegatives(2))
+  @constraint(nls, [x[2], x[2]^2 + x[2] * x[3] + 4.0] in MOI.Nonpositives(2))
+  @constraint(nls, [x[3], x[3]^2 + x[3] * x[1] + 5.0] in MOI.Zeros(2))
+
+  @NLexpression(nls, F[i = 1:3], x[i]^2 - i^2)
+
+  return MathOptNLSModel(nls, F, name = "nlsqc")
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 49d13d9..382de7e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,9 +4,11 @@ using Test, Printf
 
 nlp_problems = setdiff(NLPModelsTest.nlp_problems, ["MGH01Feas"])
 nls_problems = NLPModelsTest.nls_problems
-extra_nls_problems = ["HS30", "HS43", "MGH07", "nlsnohesspb"]
 
-for problem in lowercase.(nlp_problems ∪ ["nohesspb"])
+extra_nlp_problems = ["nohesspb", "hs61", "hs100", "hs219", "quadcon"]
+extra_nls_problems = ["nlsnohesspb", "HS30", "HS43", "MGH07", "nlsqc"]
+
+for problem in lowercase.(nlp_problems ∪ extra_nlp_problems)
   include(joinpath("nlp_problems", "$problem.jl"))
 end
 
diff --git a/test/test_moi_nlp_model.jl b/test/test_moi_nlp_model.jl
index 162694d..bfafe05 100644
--- a/test/test_moi_nlp_model.jl
+++ b/test/test_moi_nlp_model.jl
@@ -11,7 +11,7 @@ println("Testing MathOptNLPModel")
   "‖c(x₀)‖"
 )
 # Test that every problem can be instantiated.
-for prob in Symbol.(lowercase.(nlp_problems ∪ ["nohesspb"]))
+for prob in Symbol.(lowercase.(nlp_problems ∪ extra_nlp_problems))
   prob_fn = eval(prob)
   nlp = MathOptNLPModel(prob_fn(), hessian = (prob != :nohesspb))
   n = nlp.meta.nvar