diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 9cf2cf2..70878b4 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -32,6 +32,13 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v1
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MathOptInterface", rev="bl/qp_block_data"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:
diff --git a/ext/NLoptMathOptInterfaceExt.jl b/ext/NLoptMathOptInterfaceExt.jl
index 050b86c..7047bcf 100644
--- a/ext/NLoptMathOptInterfaceExt.jl
+++ b/ext/NLoptMathOptInterfaceExt.jl
@@ -16,11 +16,6 @@ function __init__()
     return
 end
 
-mutable struct _ConstraintInfo{F,S}
-    func::F
-    set::S
-end
-
 """
     Optimizer()
 
@@ -28,39 +23,19 @@ Create a new Optimizer object.
 """
 mutable struct Optimizer <: MOI.AbstractOptimizer
     inner::Union{NLopt.Opt,Nothing}
-    variables::MOI.Utilities.VariablesContainer{Float64}
+    # The variables, the affine and quadratic objective and constraints, and
+    # the inner nonlinear model.
+    model::MOI.Nonlinear.ModelWithQuad{Float64,MOI.Nonlinear.Model}
     starting_values::Vector{Union{Nothing,Float64}}
     nlp_data::MOI.NLPBlockData
-    nlp_model::Union{Nothing,MOI.Nonlinear.Model}
+    # Whether `nlp_data` was set through the legacy `MOI.NLPBlock` API, in
+    # which case it must not be rebuilt from the inner nonlinear model.
+    uses_nlp_block::Bool
     ad_backend::MOI.Nonlinear.AbstractAutomaticDifferentiation
     sense::Union{Nothing,MOI.OptimizationSense}
-    objective::Union{
-        MOI.VariableIndex,
-        MOI.ScalarAffineFunction{Float64},
-        MOI.ScalarQuadraticFunction{Float64},
-        Nothing,
-    }
-    linear_le_constraints::Vector{
-        _ConstraintInfo{
-            MOI.ScalarAffineFunction{Float64},
-            MOI.LessThan{Float64},
-        },
-    }
-    linear_eq_constraints::Vector{
-        _ConstraintInfo{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}},
-    }
-    quadratic_le_constraints::Vector{
-        _ConstraintInfo{
-            MOI.ScalarQuadraticFunction{Float64},
-            MOI.LessThan{Float64},
-        },
-    }
-    quadratic_eq_constraints::Vector{
-        _ConstraintInfo{
-            MOI.ScalarQuadraticFunction{Float64},
-            MOI.EqualTo{Float64},
-        },
-    }
+    # Whether an objective was explicitly set, because the QP block defaults
+    # to a zero objective.
+    has_objective::Bool
     # Parameters.
     silent::Bool
     options::Dict{String,Any}
@@ -73,29 +48,13 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     function Optimizer()
         return new(
             nothing,
-            MOI.Utilities.VariablesContainer{Float64}(),
+            MOI.Nonlinear.ModelWithQuad(MOI.Nonlinear.Model()),
             Union{Nothing,Float64}[],
             MOI.NLPBlockData([], _EmptyNLPEvaluator(), false),
-            nothing,
+            false,
             MOI.Nonlinear.SparseReverseMode(),
             nothing,
-            nothing,
-            _ConstraintInfo{
-                MOI.ScalarAffineFunction{Float64},
-                MOI.LessThan{Float64},
-            }[],
-            _ConstraintInfo{
-                MOI.ScalarAffineFunction{Float64},
-                MOI.EqualTo{Float64},
-            }[],
-            _ConstraintInfo{
-                MOI.ScalarQuadraticFunction{Float64},
-                MOI.LessThan{Float64},
-            }[],
-            _ConstraintInfo{
-                MOI.ScalarQuadraticFunction{Float64},
-                MOI.EqualTo{Float64},
-            }[],
+            false,
             false,
             copy(_DEFAULT_OPTIONS),
             NaN,
@@ -110,36 +69,34 @@ struct _EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
 
 MOI.initialize(::_EmptyNLPEvaluator, ::Vector{Symbol}) = nothing
 
+MOI.features_available(::_EmptyNLPEvaluator) = [:Grad, :Jac]
+
+MOI.jacobian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
+
 MOI.eval_constraint(::_EmptyNLPEvaluator, g, x) = nothing
 
 MOI.eval_constraint_jacobian(::_EmptyNLPEvaluator, J, x) = nothing
 
 function MOI.empty!(model::Optimizer)
     model.inner = nothing
-    MOI.empty!(model.variables)
+    model.model = MOI.Nonlinear.ModelWithQuad(MOI.Nonlinear.Model())
     empty!(model.starting_values)
     model.nlp_data = MOI.NLPBlockData([], _EmptyNLPEvaluator(), false)
-    model.nlp_model = nothing
+    model.uses_nlp_block = false
     model.sense = nothing
-    model.objective = nothing
-    empty!(model.linear_le_constraints)
-    empty!(model.linear_eq_constraints)
-    empty!(model.quadratic_le_constraints)
-    empty!(model.quadratic_eq_constraints)
+    model.has_objective = false
     model.status = :NOT_CALLED
     return
 end
 
 function MOI.is_empty(model::Optimizer)
-    return MOI.is_empty(model.variables) &&
+    return MOI.is_empty(model.model.variables) &&
            isempty(model.starting_values) &&
-           model.nlp_data.evaluator isa _EmptyNLPEvaluator &&
-           model.nlp_model === nothing &&
+           !model.uses_nlp_block &&
+           MOI.is_empty(model.model.inner) &&
            model.sense == nothing &&
-           isempty(model.linear_le_constraints) &&
-           isempty(model.linear_eq_constraints) &&
-           isempty(model.quadratic_le_constraints) &&
-           isempty(model.quadratic_eq_constraints)
+           !model.has_objective &&
+           length(model.model) == 0
 end
 
 function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
@@ -147,7 +104,7 @@ function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
     if model.sense !== nothing
         push!(ret, MOI.ObjectiveSense())
     end
-    if model.objective !== nothing
+    if model.has_objective
         F = MOI.get(model, MOI.ObjectiveFunctionType())
         push!(ret, MOI.ObjectiveFunction{F}())
     end
@@ -164,38 +121,6 @@ MOI.get(::Optimizer, ::MOI.SolverName) = "NLopt"
 
 MOI.get(::Optimizer, ::MOI.SolverVersion) = "$(NLopt.version())"
 
-function _constraints(
-    model,
-    ::Type{<:MOI.ScalarAffineFunction},
-    ::Type{<:MOI.LessThan},
-)
-    return model.linear_le_constraints
-end
-
-function _constraints(
-    model,
-    ::Type{<:MOI.ScalarAffineFunction},
-    ::Type{<:MOI.EqualTo},
-)
-    return model.linear_eq_constraints
-end
-
-function _constraints(
-    model,
-    ::Type{<:MOI.ScalarQuadraticFunction},
-    ::Type{<:MOI.LessThan},
-)
-    return model.quadratic_le_constraints
-end
-
-function _constraints(
-    model,
-    ::Type{<:MOI.ScalarQuadraticFunction},
-    ::Type{<:MOI.EqualTo},
-)
-    return model.quadratic_eq_constraints
-end
-
 function MOI.supports_constraint(
     ::Optimizer,
     ::Type{
@@ -219,20 +144,12 @@ function MOI.get(
     },
     S<:Union{MOI.LessThan{Float64},MOI.EqualTo{Float64}},
 }
-    return length(_constraints(model, F, S))
+    return MOI.get(model.model, MOI.NumberOfConstraints{F,S}())
 end
 
 function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintTypesPresent)
-    constraints = MOI.get(model.variables, attr)
-    function _check(model, F, S)
-        if !isempty(_constraints(model, F, S))
-            push!(constraints, (F, S))
-        end
-    end
-    _check(model, MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})
-    _check(model, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})
-    _check(model, MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})
-    _check(model, MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64})
+    constraints = MOI.get(model.model.variables, attr)
+    append!(constraints, MOI.get(model.model, attr))
     return constraints
 end
 
@@ -246,7 +163,7 @@ function MOI.get(
     },
     S<:Union{MOI.LessThan{Float64},MOI.EqualTo{Float64}},
 }
-    return MOI.ConstraintIndex{F,S}.(eachindex(_constraints(model, F, S)))
+    return MOI.get(model.model, MOI.ListOfConstraintIndices{F,S}())
 end
 
 function MOI.get(
@@ -260,7 +177,7 @@ function MOI.get(
     },
     S<:Union{MOI.LessThan{Float64},MOI.EqualTo{Float64}},
 }
-    return copy(_constraints(model, F, S)[c.value].func)
+    return MOI.get(model.model, MOI.ConstraintFunction(), c)
 end
 
 function MOI.get(
@@ -274,7 +191,7 @@ function MOI.get(
     },
     S<:Union{MOI.LessThan{Float64},MOI.EqualTo{Float64}},
 }
-    return _constraints(model, F, S)[c.value].set
+    return MOI.get(model.model, MOI.ConstraintSet(), c)
 end
 
 # ObjectiveSense
@@ -372,7 +289,7 @@ end
 
 function MOI.add_variable(model::Optimizer)
     push!(model.starting_values, nothing)
-    return MOI.add_variable(model.variables)
+    return MOI.add_variable(model.model.variables)
 end
 
 function MOI.supports_constraint(
@@ -399,7 +316,7 @@ function MOI.get(
         MOI.ListOfConstraintIndices{MOI.VariableIndex},
     },
 )
-    return MOI.get(model.variables, attr)
+    return MOI.get(model.model.variables, attr)
 end
 
 function MOI.get(
@@ -407,14 +324,14 @@ function MOI.get(
     attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
     ci::MOI.ConstraintIndex{MOI.VariableIndex},
 )
-    return MOI.get(model.variables, attr, ci)
+    return MOI.get(model.model.variables, attr, ci)
 end
 
 function MOI.is_valid(
     model::Optimizer,
     index::Union{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex}},
 )
-    return MOI.is_valid(model.variables, index)
+    return MOI.is_valid(model.model.variables, index)
 end
 
 function MOI.add_constraint(
@@ -427,7 +344,7 @@ function MOI.add_constraint(
         MOI.Interval{Float64},
     },
 )
-    return MOI.add_constraint(model.variables, vi, set)
+    return MOI.add_constraint(model.model.variables, vi, set)
 end
 
 function MOI.set(
@@ -436,14 +353,14 @@ function MOI.set(
     ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
     set::S,
 ) where {S}
-    return MOI.set(model.variables, attr, ci, set)
+    return MOI.set(model.model.variables, attr, ci, set)
 end
 
 function MOI.delete(
     model::Optimizer,
     ci::MOI.ConstraintIndex{MOI.VariableIndex},
 )
-    return MOI.delete(model.variables, ci)
+    return MOI.delete(model.model.variables, ci)
 end
 
 # constraints
@@ -462,7 +379,7 @@ function MOI.is_valid(
     },
     S<:Union{MOI.LessThan{Float64},MOI.EqualTo{Float64}},
 }
-    return 1 <= ci.value <= length(_constraints(model, F, S))
+    return MOI.is_valid(model.model, ci)
 end
 
 function _check_inbounds(model, f::MOI.VariableIndex)
@@ -499,9 +416,7 @@ function MOI.add_constraint(
     S<:Union{MOI.LessThan{Float64},MOI.EqualTo{Float64}},
 }
     _check_inbounds(model, func)
-    constraints = _constraints(model, F, S)
-    push!(constraints, _ConstraintInfo(func, set))
-    return MOI.ConstraintIndex{F,S}(length(constraints))
+    return MOI.add_constraint(model.model, func, set)
 end
 
 # MOI.VariablePrimalStart
@@ -539,10 +454,11 @@ end
 MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
 
 function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
-    if model.nlp_model !== nothing
+    if !MOI.is_empty(model.model.inner)
         error("Cannot mix the new and legacy nonlinear APIs")
     end
     model.nlp_data = nlp_data
+    model.uses_nlp_block = !(nlp_data.evaluator isa _EmptyNLPEvaluator)
     return
 end
 
@@ -562,14 +478,16 @@ function MOI.supports(
 end
 
 function MOI.get(model::Optimizer, ::MOI.ObjectiveFunctionType)
-    if model.nlp_model !== nothing && model.nlp_model.objective !== nothing
+    if model.model.inner.objective !== nothing
         return MOI.ScalarNonlinearFunction
+    elseif !model.has_objective
+        return Nothing
     end
-    return typeof(model.objective)
+    return MOI.get(model.model, MOI.ObjectiveFunctionType())
 end
 
-function MOI.get(model::Optimizer, ::MOI.ObjectiveFunction{F}) where {F}
-    return convert(F, model.objective)::F
+function MOI.get(model::Optimizer, attr::MOI.ObjectiveFunction{F}) where {F}
+    return MOI.get(model.model, attr)::F
 end
 
 function MOI.set(
@@ -584,21 +502,18 @@ function MOI.set(
     },
 }
     _check_inbounds(model, func)
-    model.objective = func
-    if model.nlp_model !== nothing
-        MOI.Nonlinear.set_objective(model.nlp_model, nothing)
-    end
+    model.has_objective = true
+    MOI.Nonlinear.set_objective(model.model, func)
     return
 end
 
 # ScalarNonlinearFunction
 
-function _init_nlp_model(model)
-    if model.nlp_model === nothing
-        if !(model.nlp_data.evaluator isa _EmptyNLPEvaluator)
-            error("Cannot mix the new and legacy nonlinear APIs")
-        end
-        model.nlp_model = MOI.Nonlinear.Model()
+# The nonlinear model is `model.model.inner` and always exists; this guard
+# only rejects mixing it with the legacy `MOI.NLPBlock` API.
+function _check_no_nlp_block(model)
+    if model.uses_nlp_block
+        error("Cannot mix the new and legacy nonlinear APIs")
     end
     return
 end
@@ -614,11 +529,8 @@ function MOI.is_valid(
         MOI.Interval{Float64},
     },
 }
-    if model.nlp_model === nothing
-        return false
-    end
     index = MOI.Nonlinear.ConstraintIndex(ci.value)
-    return MOI.is_valid(model.nlp_model, index)
+    return MOI.is_valid(model.model.inner, index)
 end
 
 function MOI.supports_constraint(
@@ -646,8 +558,8 @@ function MOI.add_constraint(
         MOI.Interval{Float64},
     },
 )
-    _init_nlp_model(model)
-    index = MOI.Nonlinear.add_constraint(model.nlp_model, f, set)
+    _check_no_nlp_block(model)
+    index = MOI.Nonlinear.add_constraint(model.model.inner, f, set)
     return MOI.ConstraintIndex{typeof(f),typeof(set)}(index.value)
 end
 
@@ -663,8 +575,8 @@ function MOI.set(
     attr::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction},
     func::MOI.ScalarNonlinearFunction,
 )
-    _init_nlp_model(model)
-    MOI.Nonlinear.set_objective(model.nlp_model, func)
+    _check_no_nlp_block(model)
+    MOI.Nonlinear.set_objective(model.model, func)
     return
 end
 
@@ -690,9 +602,9 @@ end
 MOI.supports(model::Optimizer, ::MOI.UserDefinedFunction) = true
 
 function MOI.set(model::Optimizer, attr::MOI.UserDefinedFunction, args)
-    _init_nlp_model(model)
+    _check_no_nlp_block(model)
     MOI.Nonlinear.register_operator(
-        model.nlp_model,
+        model.model.inner,
         attr.name,
         attr.arity,
         args...,
@@ -703,106 +615,34 @@ end
 ### MOI.ListOfSupportedNonlinearOperators
 
 function MOI.get(model::Optimizer, attr::MOI.ListOfSupportedNonlinearOperators)
-    _init_nlp_model(model)
-    return MOI.get(model.nlp_model, attr)
+    _check_no_nlp_block(model)
+    return MOI.get(model.model.inner, attr)
 end
 
 # optimize!
 
-function _fill_gradient(grad, x, f::MOI.VariableIndex)
-    grad[f.value] = 1.0
-    return
-end
-
-function _fill_gradient(grad, x, f::MOI.ScalarAffineFunction{Float64})
-    for term in f.terms
-        grad[term.variable.value] += term.coefficient
-    end
-    return
-end
-
-function _fill_gradient(grad, x, f::MOI.ScalarQuadraticFunction{Float64})
-    for term in f.affine_terms
-        grad[term.variable.value] += term.coefficient
-    end
-    for term in f.quadratic_terms
-        i, j = term.variable_1.value, term.variable_2.value
-        grad[i] += term.coefficient * x[j]
-        if i != j
-            grad[j] += term.coefficient * x[i]
-        end
-    end
-    return
-end
-
-function _fill_result(result::Vector, x, offset, constraints::Vector)
-    for (i, constraint) in enumerate(constraints)
-        lhs = MOI.Utilities.eval_variables(vi -> x[vi.value], constraint.func)
-        result[offset+i] = lhs - MOI.constant(constraint.set)
-    end
-    return
-end
-
-function _fill_jacobian(jac, x, offset, term::MOI.ScalarAffineTerm)
-    jac[term.variable.value, offset] += term.coefficient
-    return
-end
-
-function _fill_jacobian(jac, x, offset, term::MOI.ScalarQuadraticTerm)
-    i, j = term.variable_1.value, term.variable_2.value
-    jac[i, offset] += term.coefficient * x[j]
-    if i != j
-        jac[j, offset] += term.coefficient * x[i]
-    end
-    return
-end
-
-function _fill_jacobian(jac, x, offset, f::MOI.ScalarAffineFunction)
-    for term in f.terms
-        _fill_jacobian(jac, x, offset, term)
-    end
-    return
-end
-
-function _fill_jacobian(jac, x, offset, f::MOI.ScalarQuadraticFunction)
-    for term in f.affine_terms
-        _fill_jacobian(jac, x, offset, term)
-    end
-    for q_term in f.quadratic_terms
-        _fill_jacobian(jac, x, offset, q_term)
-    end
-    return
-end
-
-function _fill_jacobian(jac, x, offset, constraints::Vector)
-    for (i, constraint) in enumerate(constraints)
-        _fill_jacobian(jac, x, offset + i, constraint.func)
-    end
-    return
-end
-
-function objective_fn(model::Optimizer, x::Vector, grad::Vector)
-    # The order of the conditions is important. NLP objectives override regular
-    # objectives.
+function objective_fn(model::Optimizer, evaluator, x::Vector, grad::Vector)
+    # The objective sink of `evaluator` routes between the QP block and the
+    # inner nonlinear model, and evaluates to zero if no objective is set. A
+    # legacy NLPBlock objective is evaluated directly, because the sink does
+    # not know about it.
+    legacy_nlp = model.uses_nlp_block && model.nlp_data.has_objective
     if length(grad) > 0
         fill!(grad, 0.0)
         if model.sense == MOI.FEASIBILITY_SENSE
             # nothing
-        elseif model.nlp_data.has_objective
+        elseif legacy_nlp
             MOI.eval_objective_gradient(model.nlp_data.evaluator, grad, x)
-        elseif model.objective !== nothing
-            _fill_gradient(grad, x, model.objective)
+        else
+            MOI.eval_objective_gradient(evaluator, grad, x)
         end
     end
     if model.sense == MOI.FEASIBILITY_SENSE
         return 0.0
-    elseif model.nlp_data.has_objective
+    elseif legacy_nlp
         return MOI.eval_objective(model.nlp_data.evaluator, x)
-    elseif model.objective !== nothing
-        return MOI.Utilities.eval_variables(vi -> x[vi.value], model.objective)
     end
-    # No ObjectiveFunction is set, but ObjectiveSense is?
-    return 0.0
+    return MOI.eval_objective(evaluator, x)
 end
 
 function _initialize_options!(model::Optimizer)
@@ -847,91 +687,91 @@ function MOI.optimize!(model::Optimizer)
     num_variables = length(model.starting_values)
     model.inner = NLopt.Opt(model.options["algorithm"], num_variables)
     _initialize_options!(model)
-    if model.nlp_model !== nothing
-        vars = MOI.VariableIndex.(1:num_variables)
+    vars = MOI.VariableIndex.(1:num_variables)
+    # Rebuild even when the inner model is empty: a previous solve may have
+    # left a stale `nlp_data` (for example, a nonlinear objective replaced by
+    # a quadratic one).
+    if !model.uses_nlp_block
         model.nlp_data = MOI.NLPBlockData(
-            MOI.Nonlinear.Evaluator(model.nlp_model, model.ad_backend, vars),
+            MOI.Nonlinear.Evaluator(model.model.inner, model.ad_backend, vars),
         )
     end
-    NLopt.lower_bounds!(model.inner, model.variables.lower)
-    NLopt.upper_bounds!(model.inner, model.variables.upper)
-    nonlinear_equality_indices = findall(
-        bound -> bound.lower == bound.upper,
-        model.nlp_data.constraint_bounds,
-    )
-    nonlinear_inequality_indices = findall(
-        bound -> bound.lower != bound.upper,
+    # One evaluator for the whole model: the rows of the QP block come first,
+    # followed by the rows of `model.nlp_data`.
+    evaluator =
+        MOI.Nonlinear.EvaluatorWithQuad(model.model, model.nlp_data.evaluator)
+    constraint_bounds = vcat(
+        MOI.NLPBoundsPair.(model.model.qp.g_L, model.model.qp.g_U),
         model.nlp_data.constraint_bounds,
     )
-    num_nlpblock_constraints = length(model.nlp_data.constraint_bounds)
+    num_constraints = length(constraint_bounds)
+    NLopt.lower_bounds!(model.inner, model.model.variables.lower)
+    NLopt.upper_bounds!(model.inner, model.model.variables.upper)
+    equality_indices =
+        findall(bound -> bound.lower == bound.upper, constraint_bounds)
+    inequality_indices =
+        findall(bound -> bound.lower != bound.upper, constraint_bounds)
     # map from eqidx/ineqidx to index in equalities/inequalities
-    constrmap = zeros(Int, num_nlpblock_constraints)
-    for (i, k) in enumerate(nonlinear_equality_indices)
+    constrmap = zeros(Int, num_constraints)
+    for (i, k) in enumerate(equality_indices)
         constrmap[k] = i
     end
-    num_nlpblock_inequalities = 0
-    for (i, k) in enumerate(nonlinear_inequality_indices)
-        num_nlpblock_inequalities += 1
-        constrmap[k] = num_nlpblock_inequalities
-        bounds = model.nlp_data.constraint_bounds[k]
+    num_inequalities = 0
+    for k in inequality_indices
+        num_inequalities += 1
+        constrmap[k] = num_inequalities
+        bounds = constraint_bounds[k]
         if !isinf(bounds.lower) && !isinf(bounds.upper)
             # constraint has bounds on both sides, keep room for it
-            num_nlpblock_inequalities += 1
+            num_inequalities += 1
         end
     end
     if string(model.options["algorithm"])[2] == 'N'
         # Derivative free optimizer chosen
-        MOI.initialize(model.nlp_data.evaluator, Symbol[])
-    elseif num_nlpblock_constraints > 0
-        MOI.initialize(model.nlp_data.evaluator, [:Grad, :Jac])
+        MOI.initialize(evaluator, Symbol[])
+    elseif num_constraints > 0
+        MOI.initialize(evaluator, [:Grad, :Jac])
     else
-        MOI.initialize(model.nlp_data.evaluator, [:Grad])
+        MOI.initialize(evaluator, [:Grad])
     end
     if model.sense == MOI.MAX_SENSE
-        NLopt.max_objective!(model.inner, (x, g) -> objective_fn(model, x, g))
+        NLopt.max_objective!(
+            model.inner,
+            (x, g) -> objective_fn(model, evaluator, x, g),
+        )
     else
-        NLopt.min_objective!(model.inner, (x, g) -> objective_fn(model, x, g))
+        NLopt.min_objective!(
+            model.inner,
+            (x, g) -> objective_fn(model, evaluator, x, g),
+        )
     end
     Jac_IJ = Tuple{Int,Int}[]
-    if num_nlpblock_constraints > 0
-        append!(Jac_IJ, MOI.jacobian_structure(model.nlp_data.evaluator))
+    if num_constraints > 0
+        append!(Jac_IJ, MOI.jacobian_structure(evaluator))
     end
     Jac_val = zeros(length(Jac_IJ))
-    g_vec = zeros(num_nlpblock_constraints)
+    g_vec = zeros(num_constraints)
     function equality_constraint_fn(result::Vector, x::Vector, jac::Matrix)
         if length(jac) > 0
             fill!(jac, 0.0)
-            MOI.eval_constraint_jacobian(model.nlp_data.evaluator, Jac_val, x)
+            MOI.eval_constraint_jacobian(evaluator, Jac_val, x)
             for ((row, col), val) in zip(Jac_IJ, Jac_val)
-                bounds = model.nlp_data.constraint_bounds[row]
+                bounds = constraint_bounds[row]
                 if bounds.lower == bounds.upper
                     jac[col, constrmap[row]] += val
                 end
             end
-            offset = length(nonlinear_equality_indices)
-            _fill_jacobian(jac, x, offset, model.linear_eq_constraints)
-            offset += length(model.linear_eq_constraints)
-            _fill_jacobian(jac, x, offset, model.quadratic_eq_constraints)
         end
-        MOI.eval_constraint(model.nlp_data.evaluator, g_vec, x)
-        for (i, index) in enumerate(nonlinear_equality_indices)
-            bounds = model.nlp_data.constraint_bounds[index]
-            result[i] = g_vec[index] - bounds.upper
+        MOI.eval_constraint(evaluator, g_vec, x)
+        for (i, index) in enumerate(equality_indices)
+            result[i] = g_vec[index] - constraint_bounds[index].upper
         end
-        offset = length(nonlinear_equality_indices)
-        _fill_result(result, x, offset, model.linear_eq_constraints)
-        offset += length(model.linear_eq_constraints)
-        _fill_result(result, x, offset, model.quadratic_eq_constraints)
         return
     end
-    num_equality_constraints =
-        length(nonlinear_equality_indices) +
-        length(model.linear_eq_constraints) +
-        length(model.quadratic_eq_constraints)
-    if num_equality_constraints > 0
+    if length(equality_indices) > 0
         NLopt.equality_constraint!(
             model.inner,
-            num_equality_constraints,
+            length(equality_indices),
             equality_constraint_fn,
             model.options["constrtol_abs"],
         )
@@ -942,9 +782,9 @@ function MOI.optimize!(model::Optimizer)
     function inequality_constraint_fn(result::Vector, x::Vector, jac::Matrix)
         if length(jac) > 0
             fill!(jac, 0.0)
-            MOI.eval_constraint_jacobian(model.nlp_data.evaluator, Jac_val, x)
+            MOI.eval_constraint_jacobian(evaluator, Jac_val, x)
             for ((row, col), val) in zip(Jac_IJ, Jac_val)
-                bounds = model.nlp_data.constraint_bounds[row]
+                bounds = constraint_bounds[row]
                 if bounds.lower == bounds.upper
                     continue  # This is an equality constraint
                 elseif isinf(bounds.lower)  # upper bound
@@ -956,17 +796,11 @@ function MOI.optimize!(model::Optimizer)
                     jac[col, constrmap[row]+1] -= val
                 end
             end
-            offset = num_nlpblock_inequalities
-            _fill_jacobian(jac, x, offset, model.linear_le_constraints)
-            offset += length(model.linear_le_constraints)
-            _fill_jacobian(jac, x, offset, model.quadratic_le_constraints)
         end
-        # Fill in the result. The first entries are from NLPBlock, and the value
-        # of g(x) is placed in g_vec.
-        MOI.eval_constraint(model.nlp_data.evaluator, g_vec, x)
-        for row in 1:num_nlpblock_constraints
+        MOI.eval_constraint(evaluator, g_vec, x)
+        for row in 1:num_constraints
             index = constrmap[row]
-            bounds = model.nlp_data.constraint_bounds[row]
+            bounds = constraint_bounds[row]
             if bounds.lower == bounds.upper
                 continue  # This is an equality constraint
             elseif isinf(bounds.lower)  # g(x) <= u --> g(x) - u <= 0
@@ -978,20 +812,12 @@ function MOI.optimize!(model::Optimizer)
                 result[index+1] = bounds.lower - g_vec[row]
             end
         end
-        offset = num_nlpblock_inequalities
-        _fill_result(result, x, offset, model.linear_le_constraints)
-        offset += length(model.linear_le_constraints)
-        _fill_result(result, x, offset, model.quadratic_le_constraints)
         return
     end
-    num_inequality_constraints =
-        num_nlpblock_inequalities +
-        length(model.linear_le_constraints) +
-        length(model.quadratic_le_constraints)
-    if num_inequality_constraints > 0
+    if num_inequalities > 0
         NLopt.inequality_constraint!(
             model.inner,
-            num_inequality_constraints,
+            num_inequalities,
             inequality_constraint_fn,
             model.options["constrtol_abs"],
         )
@@ -999,7 +825,7 @@ function MOI.optimize!(model::Optimizer)
     # Set MOI.VariablePrimalStart, clamping to bound nearest 0 if not given.
     model.solution = something.(
         model.starting_values,
-        clamp.(0.0, model.variables.lower, model.variables.upper),
+        clamp.(0.0, model.model.variables.lower, model.model.variables.upper),
     )
     start_time = time()
     model.objective_value, _, model.status =