From 2d894db0a7ab946ea45d534143316f903667522c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 13:42:46 +0000
Subject: [PATCH 1/7] Add ModelWithQuad and privatize the product functions of
 QPBlockData

Making a first release with public add_constraint_jacobian_product,
add_constraint_jacobian_transpose_product, and
add_hessian_lagrangian_product commits MOI to three new function names.
Instead, rename them _add_... and add ModelWithQuad/EvaluatorWithQuad:
the evaluator implements the standard MOI.AbstractNLPEvaluator
interface (including the product callbacks, composed inner-first so
that store-semantics inner implementations such as ReverseAD do not
clobber the QP contribution), so consumers only rely on existing
generic functions and the QPBlockData internals stay private.

ModelWithQuad routes affine and quadratic objectives and constraints
to a QPBlockData and everything else to an inner model, forwards the
MOI attribute queries of the QP block, and tracks where the objective
lives. EvaluatorWithQuad remaps the variables of the QP block to their
consecutive index in ordered_variables during MOI.initialize (variables
absent from ordered_variables are parameters and keep their index,
with the parameters dictionary aliased so per-solve value syncs are
visible), and MOI.NLPBlockData(evaluator) assembles the combined
constraint bounds.
---
 docs/src/submodules/Nonlinear/reference.md |   5 +-
 src/Nonlinear/Nonlinear.jl                 |   1 +
 src/Nonlinear/model_with_quad.jl           | 478 +++++++++++++++++++++
 src/Nonlinear/qp_block_data.jl             |  12 +-
 test/Nonlinear/test_model_with_quad.jl     | 309 +++++++++++++
 5 files changed, 796 insertions(+), 9 deletions(-)
 create mode 100644 src/Nonlinear/model_with_quad.jl
 create mode 100644 test/Nonlinear/test_model_with_quad.jl

diff --git a/docs/src/submodules/Nonlinear/reference.md b/docs/src/submodules/Nonlinear/reference.md
index 468f44cdfe..f5c3746069 100644
--- a/docs/src/submodules/Nonlinear/reference.md
+++ b/docs/src/submodules/Nonlinear/reference.md
@@ -73,9 +73,8 @@ Nonlinear.ExprGraphOnly
 Nonlinear.SparseReverseMode
 Nonlinear.SymbolicMode
 Nonlinear.QPBlockData
-Nonlinear.add_constraint_jacobian_product
-Nonlinear.add_constraint_jacobian_transpose_product
-Nonlinear.add_hessian_lagrangian_product
+Nonlinear.ModelWithQuad
+Nonlinear.EvaluatorWithQuad
 
 ```
 
diff --git a/src/Nonlinear/Nonlinear.jl b/src/Nonlinear/Nonlinear.jl
index e057052029..2475fba33c 100644
--- a/src/Nonlinear/Nonlinear.jl
+++ b/src/Nonlinear/Nonlinear.jl
@@ -43,5 +43,6 @@ include("ReverseAD/ReverseAD.jl")
 include("SymbolicAD/SymbolicAD.jl")
 
 include("qp_block_data.jl")
+include("model_with_quad.jl")
 
 end  # module
diff --git a/src/Nonlinear/model_with_quad.jl b/src/Nonlinear/model_with_quad.jl
new file mode 100644
index 0000000000..517e2ac4b8
--- /dev/null
+++ b/src/Nonlinear/model_with_quad.jl
@@ -0,0 +1,478 @@
+# Copyright (c) 2017: Miles Lubin and contributors
+# Copyright (c) 2017: Google Inc.
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+"""
+    ModelWithQuad{T,M}(
+        qp::QPBlockData{T},
+        inner::M;
+        objective_sink::Symbol = :none,
+    ) where {T,M}
+
+A model layer that stores affine and quadratic objectives and constraints in
+a [`QPBlockData`](@ref) and forwards everything else to the `inner` model,
+typically a [`Model`](@ref). Solvers with their own storage can pass an
+existing `qp` and set `inner` to `nothing` when they manage the nonlinear
+model themselves.
+
+`ModelWithQuad(inner)` and `ModelWithQuad{T}(inner)` create an empty
+[`QPBlockData`](@ref), with `T` defaulting to `Float64`.
+
+Add constraints with [`add_constraint`](@ref) or `MOI.add_constraint`, and
+set the objective with [`set_objective`](@ref) or `MOI.set`: affine and
+quadratic functions are routed to the QP block, everything else to the inner
+model. `objective_sink` records where the objective currently lives (`:none`,
+`:quad` or `:inner`).
+
+Create the corresponding evaluator, [`EvaluatorWithQuad`](@ref), with
+`Evaluator(model, backend, ordered_variables)`, or construct it directly from
+an inner `MOI.AbstractNLPEvaluator`. The rows of the QP block come first,
+followed by the rows of the inner evaluator.
+
+Functions added to the QP block may contain parameters, following the
+convention documented in [`QPBlockData`](@ref): a variable is a parameter if
+and only if its index is a key of `qp.parameters`.
+"""
+mutable struct ModelWithQuad{T,M}
+    qp::QPBlockData{T}
+    inner::M
+    objective_sink::Symbol # :none, :quad or :inner
+
+    function ModelWithQuad{T}(
+        qp::QPBlockData{T},
+        inner::M;
+        objective_sink::Symbol = :none,
+    ) where {T,M}
+        return new{T,M}(qp, inner, objective_sink)
+    end
+end
+
+function ModelWithQuad{T}(inner) where {T}
+    return ModelWithQuad{T}(QPBlockData{T}(), inner)
+end
+
+ModelWithQuad(inner) = ModelWithQuad{Float64}(inner)
+
+"""
+    Base.length(model::ModelWithQuad)
+
+The number of affine and quadratic constraint rows of `model`, which come
+before the rows of the inner model in the corresponding evaluator.
+"""
+Base.length(model::ModelWithQuad) = length(model.qp)
+
+# Methods forwarded to the inner model.
+
+function add_parameter(model::ModelWithQuad, value::Real)
+    return add_parameter(model.inner, value)
+end
+
+add_expression(model::ModelWithQuad, expr) = add_expression(model.inner, expr)
+
+Base.getindex(model::ModelWithQuad, index::ExpressionIndex) = model.inner[index]
+
+function register_operator(
+    model::ModelWithQuad,
+    op::Symbol,
+    nargs::Int,
+    f::Function...,
+)
+    return register_operator(model.inner, op, nargs, f...)
+end
+
+function MOI.is_valid(model::ModelWithQuad, index::ConstraintIndex)
+    return MOI.is_valid(model.inner, index)
+end
+
+function MOI.get(
+    model::ModelWithQuad,
+    attr::MOI.ListOfSupportedNonlinearOperators,
+)
+    return MOI.get(model.inner, attr)
+end
+
+# The objective.
+
+# Support for models whose inner model is managed by the solver.
+set_objective(::Nothing, ::Nothing) = nothing
+
+function set_objective(
+    model::ModelWithQuad{T},
+    obj::Union{
+        MOI.VariableIndex,
+        MOI.ScalarAffineFunction{T},
+        MOI.ScalarQuadraticFunction{T},
+    },
+) where {T}
+    MOI.set(model.qp, MOI.ObjectiveFunction{typeof(obj)}(), obj)
+    set_objective(model.inner, nothing)
+    model.objective_sink = :quad
+    return
+end
+
+function set_objective(model::ModelWithQuad{T}, obj) where {T}
+    F = MOI.ScalarAffineFunction{T}
+    MOI.set(model.qp, MOI.ObjectiveFunction{F}(), zero(F))
+    set_objective(model.inner, obj)
+    model.objective_sink = obj === nothing ? :none : :inner
+    return
+end
+
+function MOI.set(
+    model::ModelWithQuad{T},
+    ::MOI.ObjectiveFunction{F},
+    func::F,
+) where {
+    T,
+    F<:Union{
+        MOI.VariableIndex,
+        MOI.ScalarAffineFunction{T},
+        MOI.ScalarQuadraticFunction{T},
+    },
+}
+    set_objective(model, func)
+    return
+end
+
+function MOI.get(model::ModelWithQuad, attr::MOI.ObjectiveFunctionType)
+    return MOI.get(model.qp, attr)
+end
+
+function MOI.get(model::ModelWithQuad, attr::MOI.ObjectiveFunction{F}) where {F}
+    return MOI.get(model.qp, attr)
+end
+
+# The affine and quadratic constraints. The MOI attribute methods are
+# forwarded to the QP block, which implements them.
+
+const _QPFunction{T} =
+    Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}
+
+const _QPSet{T} =
+    Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}}
+
+function add_constraint(
+    model::ModelWithQuad{T},
+    func::_QPFunction{T},
+    set::_QPSet{T},
+) where {T}
+    return MOI.add_constraint(model.qp, func, set)
+end
+
+function add_constraint(model::ModelWithQuad, func, set)
+    return add_constraint(model.inner, func, set)
+end
+
+function MOI.add_constraint(
+    model::ModelWithQuad{T},
+    func::_QPFunction{T},
+    set::_QPSet{T},
+) where {T}
+    return MOI.add_constraint(model.qp, func, set)
+end
+
+function MOI.get(model::ModelWithQuad, attr::MOI.ListOfConstraintTypesPresent)
+    return MOI.get(model.qp, attr)
+end
+
+function MOI.is_valid(
+    model::ModelWithQuad{T},
+    ci::MOI.ConstraintIndex{F,S},
+) where {T,F<:_QPFunction{T},S<:_QPSet{T}}
+    return MOI.is_valid(model.qp, ci)
+end
+
+function MOI.get(
+    model::ModelWithQuad{T},
+    attr::Union{MOI.ListOfConstraintIndices{F,S},MOI.NumberOfConstraints{F,S}},
+) where {T,F<:_QPFunction{T},S<:_QPSet{T}}
+    return MOI.get(model.qp, attr)
+end
+
+function MOI.get(
+    model::ModelWithQuad{T},
+    attr::Union{
+        MOI.ConstraintFunction,
+        MOI.ConstraintSet,
+        MOI.ConstraintDualStart,
+    },
+    ci::MOI.ConstraintIndex{F,S},
+) where {T,F<:_QPFunction{T},S<:_QPSet{T}}
+    return MOI.get(model.qp, attr, ci)
+end
+
+function MOI.set(
+    model::ModelWithQuad{T},
+    attr::MOI.ConstraintSet,
+    ci::MOI.ConstraintIndex{F,S},
+    set::S,
+) where {T,F<:_QPFunction{T},S<:_QPSet{T}}
+    return MOI.set(model.qp, attr, ci, set)
+end
+
+function MOI.set(
+    model::ModelWithQuad{T},
+    attr::MOI.ConstraintDualStart,
+    ci::MOI.ConstraintIndex{F,S},
+    value,
+) where {T,F<:_QPFunction{T},S<:_QPSet{T}}
+    return MOI.set(model.qp, attr, ci, value)
+end
+
+"""
+    EvaluatorWithQuad(
+        model::ModelWithQuad,
+        inner::MOI.AbstractNLPEvaluator,
+        ordered_variables::Vector{MOI.VariableIndex},
+    ) <: MOI.AbstractNLPEvaluator
+
+The evaluator of a [`ModelWithQuad`](@ref) layer. It implements the
+[`MOI.AbstractNLPEvaluator`](@ref) interface: the rows of the QP block come
+first, followed by the rows of `inner`, and the Jacobian and Hessian product
+callbacks compose the contributions of the two blocks.
+
+Create it with `Evaluator(model::ModelWithQuad, backend, ordered_variables)`,
+which recursively creates the evaluator of the inner model, or construct it
+directly from an existing inner evaluator.
+
+During [`MOI.initialize`](@ref), the variables of the QP block are mapped to
+their consecutive 1-based index in `ordered_variables`; variables absent from
+`ordered_variables` are parameters and keep their index.
+"""
+mutable struct EvaluatorWithQuad{T,M,E<:MOI.AbstractNLPEvaluator} <:
+               MOI.AbstractNLPEvaluator
+    model::ModelWithQuad{T,M}
+    inner::E
+    ordered_variables::Vector{MOI.VariableIndex}
+    # A copy of `model.qp` with the variables mapped to their consecutive
+    # 1-based index in `ordered_variables`, and the number of entries of its
+    # Jacobian and of its Hessian of the Lagrangian. Rebuilt during
+    # `MOI.initialize`.
+    qp::QPBlockData{T}
+    qp_nnzj::Int
+    qp_nnzh::Int
+
+    function EvaluatorWithQuad(
+        model::ModelWithQuad{T,M},
+        inner::E,
+        ordered_variables::Vector{MOI.VariableIndex},
+    ) where {T,M,E<:MOI.AbstractNLPEvaluator}
+        return new{T,M,E}(
+            model,
+            inner,
+            ordered_variables,
+            QPBlockData{T}(),
+            0,
+            0,
+        )
+    end
+end
+
+function Evaluator(
+    model::ModelWithQuad,
+    backend::AbstractAutomaticDifferentiation,
+    ordered_variables::Vector{MOI.VariableIndex},
+)
+    inner = Evaluator(model.inner, backend, ordered_variables)
+    return EvaluatorWithQuad(model, inner, ordered_variables)
+end
+
+function MOI.features_available(d::EvaluatorWithQuad)
+    features = MOI.features_available(d.inner)
+    return filter(f -> f in (:Grad, :Jac, :JacVec, :Hess, :HessVec), features)
+end
+
+function MOI.initialize(
+    d::EvaluatorWithQuad{T},
+    features::Vector{Symbol},
+) where {T}
+    index_map = Dict{MOI.VariableIndex,MOI.VariableIndex}(
+        x => MOI.VariableIndex(i) for (i, x) in enumerate(d.ordered_variables)
+    )
+    # Variables absent from `ordered_variables` are parameters: they keep
+    # their index, which the `parameters` dictionary uses as key.
+    fmap = v::MOI.VariableIndex -> get(index_map, v, v)
+    src = d.model.qp
+    qp = QPBlockData{T}()
+    qp.objective = MOI.Utilities.map_indices(fmap, src.objective)
+    qp.objective_function_type = src.objective_function_type
+    for f in src.constraints
+        push!(qp.constraints, MOI.Utilities.map_indices(fmap, f))
+    end
+    append!(qp.g_L, src.g_L)
+    append!(qp.g_U, src.g_U)
+    append!(qp.mult_g, src.mult_g)
+    append!(qp.function_type, src.function_type)
+    append!(qp.bound_type, src.bound_type)
+    # Alias, do not copy: parameter values updated in the model between
+    # solves must be visible to the evaluator without re-initializing.
+    qp.parameters = src.parameters
+    d.qp = qp
+    d.qp_nnzj = length(MOI.jacobian_structure(qp))
+    d.qp_nnzh = length(MOI.hessian_lagrangian_structure(qp))
+    MOI.initialize(d.inner, features)
+    return
+end
+
+function MOI.eval_objective(d::EvaluatorWithQuad{T}, x) where {T}
+    sink = d.model.objective_sink
+    if sink == :quad
+        return MOI.eval_objective(d.qp, x)
+    elseif sink == :inner
+        return MOI.eval_objective(d.inner, x)
+    else
+        return zero(T)
+    end
+end
+
+function MOI.eval_objective_gradient(d::EvaluatorWithQuad{T}, grad, x) where {T}
+    sink = d.model.objective_sink
+    if sink == :quad
+        MOI.eval_objective_gradient(d.qp, grad, x)
+    elseif sink == :inner
+        MOI.eval_objective_gradient(d.inner, grad, x)
+    else
+        grad .= zero(T)
+    end
+    return
+end
+
+function MOI.eval_constraint(d::EvaluatorWithQuad, g, x)
+    m = length(d.qp)
+    MOI.eval_constraint(d.qp, view(g, 1:m), x)
+    MOI.eval_constraint(d.inner, view(g, (m+1):length(g)), x)
+    return
+end
+
+function MOI.jacobian_structure(d::EvaluatorWithQuad)
+    J = MOI.jacobian_structure(d.qp)
+    offset = length(d.qp)
+    # An evaluator is only required to implement `jacobian_structure` if it
+    # supports `:Jac`. If the inner evaluator does not (it then must not have
+    # any rows for the stack to be usable), append nothing.
+    if :Jac in MOI.features_available(d.inner)
+        for (row, col) in MOI.jacobian_structure(d.inner)
+            push!(J, (row + offset, col))
+        end
+    end
+    return J
+end
+
+function MOI.eval_constraint_jacobian(d::EvaluatorWithQuad, J, x)
+    MOI.eval_constraint_jacobian(d.qp, J, x)
+    MOI.eval_constraint_jacobian(d.inner, view(J, (d.qp_nnzj+1):length(J)), x)
+    return
+end
+
+function MOI.hessian_lagrangian_structure(d::EvaluatorWithQuad)
+    H = MOI.hessian_lagrangian_structure(d.qp)
+    if :Hess in MOI.features_available(d.inner)
+        append!(H, MOI.hessian_lagrangian_structure(d.inner))
+    end
+    return H
+end
+
+function MOI.eval_hessian_lagrangian(d::EvaluatorWithQuad, H, x, σ, μ)
+    m = length(d.qp)
+    # If the objective is not in the QP block, `d.qp.objective` is zero, so
+    # passing `σ` is harmless; and vice versa for the inner evaluator.
+    MOI.eval_hessian_lagrangian(d.qp, H, x, σ, view(μ, 1:m))
+    MOI.eval_hessian_lagrangian(
+        d.inner,
+        view(H, (d.qp_nnzh+1):length(H)),
+        x,
+        σ,
+        view(μ, (m+1):length(μ)),
+    )
+    return
+end
+
+# The rows of the two blocks are disjoint, so zero everything and let each
+# block write its own rows.
+function MOI.eval_constraint_jacobian_product(d::EvaluatorWithQuad, y, x, w)
+    fill!(y, zero(eltype(y)))
+    m = length(d.qp)
+    MOI.eval_constraint_jacobian_product(
+        d.inner,
+        view(y, (m+1):length(y)),
+        x,
+        w,
+    )
+    _add_constraint_jacobian_product(d.qp, y, x, w)
+    return
+end
+
+# Both blocks accumulate into the same variable-dimensional output. Call the
+# inner evaluator FIRST because implementations are allowed to overwrite the
+# output, and accumulate the QP block afterwards.
+function MOI.eval_constraint_jacobian_transpose_product(
+    d::EvaluatorWithQuad,
+    y,
+    x,
+    w,
+)
+    fill!(y, zero(eltype(y)))
+    m = length(d.qp)
+    MOI.eval_constraint_jacobian_transpose_product(
+        d.inner,
+        y,
+        x,
+        view(w, (m+1):length(w)),
+    )
+    _add_constraint_jacobian_transpose_product(d.qp, y, x, view(w, 1:m))
+    return
+end
+
+function MOI.eval_hessian_lagrangian_product(
+    d::EvaluatorWithQuad,
+    H,
+    x,
+    v,
+    σ,
+    μ,
+)
+    fill!(H, zero(eltype(H)))
+    m = length(d.qp)
+    MOI.eval_hessian_lagrangian_product(
+        d.inner,
+        H,
+        x,
+        v,
+        σ,
+        view(μ, (m+1):length(μ)),
+    )
+    _add_hessian_lagrangian_product(d.qp, H, x, v, σ, view(μ, 1:m))
+    return
+end
+
+# The lower and upper bounds of each constraint row, in the row order of the
+# evaluator. Solvers that use their own inner evaluator type can add a method
+# for it so that `MOI.NLPBlockData(::EvaluatorWithQuad)` works.
+function _constraint_bounds(evaluator::Evaluator)
+    return MOI.NLPBoundsPair[
+        _bound(c.set) for (_, c) in evaluator.model.constraints
+    ]
+end
+
+function _constraint_bounds(d::EvaluatorWithQuad)
+    bounds = MOI.NLPBoundsPair[
+        MOI.NLPBoundsPair(l, u) for
+        (l, u) in zip(d.model.qp.g_L, d.model.qp.g_U)
+    ]
+    return append!(bounds, _constraint_bounds(d.inner))
+end
+
+_has_objective(d::Evaluator) = d.model.objective !== nothing
+
+function _has_objective(d::EvaluatorWithQuad)
+    if d.model.objective_sink == :quad
+        return true
+    end
+    return _has_objective(d.inner)
+end
+
+function MOI.NLPBlockData(d::EvaluatorWithQuad)
+    return MOI.NLPBlockData(_constraint_bounds(d), d, _has_objective(d))
+end
diff --git a/src/Nonlinear/qp_block_data.jl b/src/Nonlinear/qp_block_data.jl
index 85e6ad8ce4..22881cef71 100644
--- a/src/Nonlinear/qp_block_data.jl
+++ b/src/Nonlinear/qp_block_data.jl
@@ -705,7 +705,7 @@ end
 # These are used to add the QP contribution on top of the NL contribution.
 
 """
-    add_constraint_jacobian_product(
+    _add_constraint_jacobian_product(
         block::QPBlockData{T},
         y::AbstractVector{T},
         x::AbstractVector{T},
@@ -720,7 +720,7 @@ accumulates into `y` instead of storing the result, so that the contributions
 of several blocks can be composed: the caller is responsible for zeroing `y`
 before the first contribution.
 """
-function add_constraint_jacobian_product(
+function _add_constraint_jacobian_product(
     block::QPBlockData{T},
     y::AbstractVector{T},
     x::AbstractVector{T},
@@ -733,7 +733,7 @@ function add_constraint_jacobian_product(
 end
 
 """
-    add_constraint_jacobian_transpose_product(
+    _add_constraint_jacobian_transpose_product(
         block::QPBlockData{T},
         y::AbstractVector{T},
         x::AbstractVector{T},
@@ -748,7 +748,7 @@ function accumulates into `y` instead of storing the result, so that the
 contributions of several blocks can be composed: the caller is responsible
 for zeroing `y` before the first contribution.
 """
-function add_constraint_jacobian_transpose_product(
+function _add_constraint_jacobian_transpose_product(
     block::QPBlockData{T},
     y::AbstractVector{T},
     x::AbstractVector{T},
@@ -761,7 +761,7 @@ function add_constraint_jacobian_transpose_product(
 end
 
 """
-    add_hessian_lagrangian_product(
+    _add_hessian_lagrangian_product(
         block::QPBlockData{T},
         H::AbstractVector{T},
         x::AbstractVector{T},
@@ -778,7 +778,7 @@ accumulates into `H` instead of storing the result, so that the contributions
 of several blocks can be composed: the caller is responsible for zeroing `H`
 before the first contribution.
 """
-function add_hessian_lagrangian_product(
+function _add_hessian_lagrangian_product(
     block::QPBlockData{T},
     H::AbstractVector{T},
     x::AbstractVector{T},
diff --git a/test/Nonlinear/test_model_with_quad.jl b/test/Nonlinear/test_model_with_quad.jl
new file mode 100644
index 0000000000..a3a8967dc9
--- /dev/null
+++ b/test/Nonlinear/test_model_with_quad.jl
@@ -0,0 +1,309 @@
+# Copyright (c) 2017: Miles Lubin and contributors
+# Copyright (c) 2017: Google Inc.
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+module TestNonlinearModelWithQuad
+
+using Test
+import MathOptInterface as MOI
+
+import MathOptInterface.Nonlinear
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+# A model with, in row order:
+#   row 1 (quad layer, linear):    2x + 3y <= 4
+#   row 2 (quad layer, quadratic): x^2 + xy + y in [0, 1]
+#   row 3 (inner nlp):             sin(x) <= 0.5
+# and the objective x^2 in the quad layer.
+function _test_model(x, y)
+    model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    Nonlinear.set_objective(
+        model,
+        MOI.ScalarQuadraticFunction(
+            [MOI.ScalarQuadraticTerm(2.0, x, x)],
+            MOI.ScalarAffineTerm{Float64}[],
+            0.0,
+        ),
+    )
+    c1 = MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(2.0, x), MOI.ScalarAffineTerm(3.0, y)],
+            0.0,
+        ),
+        MOI.LessThan(4.0),
+    )
+    @test c1 isa MOI.ConstraintIndex{
+        MOI.ScalarAffineFunction{Float64},
+        MOI.LessThan{Float64},
+    }
+    c2 = Nonlinear.add_constraint(
+        model,
+        MOI.ScalarQuadraticFunction(
+            [
+                MOI.ScalarQuadraticTerm(2.0, x, x),
+                MOI.ScalarQuadraticTerm(1.0, x, y),
+            ],
+            [MOI.ScalarAffineTerm(1.0, y)],
+            0.0,
+        ),
+        MOI.Interval(0.0, 1.0),
+    )
+    c3 = Nonlinear.add_constraint(model, :(sin($x)), MOI.LessThan(0.5))
+    @test c3 isa Nonlinear.ConstraintIndex
+    @test length(model) == 2
+    return model
+end
+
+function test_evaluator_with_quad()
+    # Use non-consecutive variable indices to test that the evaluator remaps
+    # the variables to their consecutive index in `ordered_variables`.
+    x, y = MOI.VariableIndex(5), MOI.VariableIndex(9)
+    model = _test_model(x, y)
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
+    @test d isa Nonlinear.EvaluatorWithQuad
+    @test d.inner isa Nonlinear.Evaluator
+    @test MOI.features_available(d) == [:Grad, :Jac, :JacVec, :Hess, :HessVec]
+    MOI.initialize(d, [:Grad, :Jac, :Hess])
+    xv = [1.0, 2.0]  # x = 1, y = 2
+    @test MOI.eval_objective(d, xv) == 1.0
+    grad = fill(NaN, 2)
+    MOI.eval_objective_gradient(d, grad, xv)
+    @test grad == [2.0, 0.0]
+    g = fill(NaN, 3)
+    MOI.eval_constraint(d, g, xv)
+    @test g ≈ [8.0, 5.0, sin(1.0)]
+    # Jacobian: accumulate the sparse entries into a dense matrix.
+    J_structure = MOI.jacobian_structure(d)
+    J_values = fill(NaN, length(J_structure))
+    MOI.eval_constraint_jacobian(d, J_values, xv)
+    J = zeros(3, 2)
+    for ((row, col), value) in zip(J_structure, J_values)
+        J[row, col] += value
+    end
+    @test J ≈ [
+        2.0 3.0
+        4.0 2.0
+        cos(1.0) 0.0
+    ]
+    # Hessian of the Lagrangian: accumulate into a dense matrix.
+    H_structure = MOI.hessian_lagrangian_structure(d)
+    σ, μ = 2.0, [10.0, 100.0, 1_000.0]
+    H_values = fill(NaN, length(H_structure))
+    MOI.eval_hessian_lagrangian(d, H_values, xv, σ, μ)
+    H = zeros(2, 2)
+    for ((row, col), value) in zip(H_structure, H_values)
+        H[row, col] += value
+    end
+    # σ * ∇²(x^2) + μ₂ * ∇²(x^2 + xy) + μ₃ * ∇²(sin(x))
+    @test H[1, 1] ≈ 2σ + 2 * μ[2] - sin(1.0) * μ[3]
+    @test H[1, 2] + H[2, 1] ≈ μ[2]
+    @test H[2, 2] ≈ 0.0
+    block = MOI.NLPBlockData(d)
+    @test block.has_objective
+    @test block.constraint_bounds == [
+        MOI.NLPBoundsPair(-Inf, 4.0),
+        MOI.NLPBoundsPair(0.0, 1.0),
+        MOI.NLPBoundsPair(-Inf, 0.5),
+    ]
+    return
+end
+
+function test_evaluator_products()
+    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
+    model = _test_model(x, y)
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
+    MOI.initialize(d, [:Grad, :Jac, :JacVec, :Hess, :HessVec])
+    xv = [1.0, 2.0]
+    # Dense Jacobian from the sparse callback, as the reference.
+    J_structure = MOI.jacobian_structure(d)
+    J_values = fill(NaN, length(J_structure))
+    MOI.eval_constraint_jacobian(d, J_values, xv)
+    J = zeros(3, 2)
+    for ((row, col), value) in zip(J_structure, J_values)
+        J[row, col] += value
+    end
+    w = [1.0, -2.0]
+    Jv = fill(NaN, 3)
+    MOI.eval_constraint_jacobian_product(d, Jv, xv, w)
+    @test Jv ≈ J * w
+    u = [1.0, -1.0, 2.0]
+    Jtv = fill(NaN, 2)
+    MOI.eval_constraint_jacobian_transpose_product(d, Jtv, xv, u)
+    @test Jtv ≈ J' * u
+    # Dense Hessian of the Lagrangian, as the reference.
+    H_structure = MOI.hessian_lagrangian_structure(d)
+    σ, μ = 2.0, [10.0, 100.0, 1_000.0]
+    H_values = fill(NaN, length(H_structure))
+    MOI.eval_hessian_lagrangian(d, H_values, xv, σ, μ)
+    H = zeros(2, 2)
+    for ((row, col), value) in zip(H_structure, H_values)
+        H[row, col] += value
+        if row != col
+            H[col, row] += value
+        end
+    end
+    v = [1.0, -3.0]
+    Hv = fill(NaN, 2)
+    MOI.eval_hessian_lagrangian_product(d, Hv, xv, v, σ, μ)
+    @test Hv ≈ H * v
+    return
+end
+
+function test_objective_sink_switching()
+    x = MOI.VariableIndex(1)
+    model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    @test model.objective_sink == :none
+    f = MOI.ScalarQuadraticFunction(
+        [MOI.ScalarQuadraticTerm(2.0, x, x)],
+        MOI.ScalarAffineTerm{Float64}[],
+        0.0,
+    )
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    @test model.objective_sink == :quad
+    @test MOI.get(model, MOI.ObjectiveFunctionType()) ==
+          MOI.ScalarQuadraticFunction{Float64}
+    @test MOI.get(model, MOI.ObjectiveFunction{typeof(f)}()) ≈ f
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac, :Hess])
+    @test MOI.eval_objective(d, [3.0]) == 9.0
+    @test MOI.NLPBlockData(d).has_objective
+    # Switch to a nonlinear objective: the quadratic objective must be
+    # cleared, including its Hessian entries.
+    Nonlinear.set_objective(model, :(sin($x)))
+    @test model.objective_sink == :inner
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac, :Hess])
+    @test MOI.eval_objective(d, [3.0]) == sin(3.0)
+    @test MOI.NLPBlockData(d).has_objective
+    H_structure = MOI.hessian_lagrangian_structure(d)
+    H = fill(NaN, length(H_structure))
+    MOI.eval_hessian_lagrangian(d, H, [3.0], 1.0, Float64[])
+    @test sum(H) ≈ -sin(3.0)
+    # Switch to a linear objective, and then remove it.
+    g = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0)
+    Nonlinear.set_objective(model, g)
+    @test model.objective_sink == :quad
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac])
+    @test MOI.eval_objective(d, [3.0]) == 7.0
+    Nonlinear.set_objective(model, nothing)
+    @test model.objective_sink == :none
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac])
+    @test MOI.eval_objective(d, [3.0]) == 0.0
+    grad = fill(NaN, 1)
+    MOI.eval_objective_gradient(d, grad, [3.0])
+    @test grad == [0.0]
+    @test !MOI.NLPBlockData(d).has_objective
+    return
+end
+
+function test_quad_parameters()
+    x, p = MOI.VariableIndex(1), MOI.VariableIndex(42)
+    model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    # Register the parameter before any structure query. Its value may be
+    # updated between evaluations.
+    model.qp.parameters[p.value] = 5.0
+    Nonlinear.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(2.0, x), MOI.ScalarAffineTerm(3.0, p)],
+            0.0,
+        ),
+        MOI.LessThan(10.0),
+    )
+    Nonlinear.add_constraint(
+        model,
+        MOI.ScalarQuadraticFunction(
+            [MOI.ScalarQuadraticTerm(1.0, p, x)],
+            MOI.ScalarAffineTerm{Float64}[],
+            0.0,
+        ),
+        MOI.LessThan(10.0),
+    )
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac, :Hess])
+    g = fill(NaN, 2)
+    MOI.eval_constraint(d, g, [1.0])
+    @test g == [2.0 * 1.0 + 3.0 * 5.0, 5.0 * 1.0]
+    # Parameters never appear in the Jacobian or Hessian structure.
+    @test MOI.jacobian_structure(d) == [(1, 1), (2, 1)]
+    J = fill(NaN, 2)
+    MOI.eval_constraint_jacobian(d, J, [1.0])
+    @test J == [2.0, 5.0]
+    @test isempty(MOI.hessian_lagrangian_structure(d))
+    # Updating the parameter value must be visible without re-initializing.
+    model.qp.parameters[p.value] = 7.0
+    MOI.eval_constraint(d, g, [1.0])
+    @test g == [2.0 * 1.0 + 3.0 * 7.0, 7.0 * 1.0]
+    return
+end
+
+function test_attribute_forwarding()
+    x = MOI.VariableIndex(1)
+    model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    F, S = MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}
+    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
+    ci = MOI.add_constraint(model, f, MOI.GreaterThan(1.0))
+    @test MOI.is_valid(model, ci)
+    @test !MOI.is_valid(model, typeof(ci)(ci.value + 1))
+    @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 1
+    @test MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) == [ci]
+    @test (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
+    @test MOI.get(model, MOI.ConstraintFunction(), ci) ≈ f
+    @test MOI.get(model, MOI.ConstraintSet(), ci) == MOI.GreaterThan(1.0)
+    MOI.set(model, MOI.ConstraintSet(), ci, MOI.GreaterThan(2.0))
+    @test MOI.get(model, MOI.ConstraintSet(), ci) == MOI.GreaterThan(2.0)
+    @test MOI.get(model, MOI.ConstraintDualStart(), ci) === nothing
+    MOI.set(model, MOI.ConstraintDualStart(), ci, 1.5)
+    @test MOI.get(model, MOI.ConstraintDualStart(), ci) == 1.5
+    # Nonlinear-model forwarding
+    p = Nonlinear.add_parameter(model, 2.0)
+    @test p isa Nonlinear.ParameterIndex
+    ex = Nonlinear.add_expression(model, :($p * $x))
+    @test model[ex] isa Nonlinear.Expression
+    Nonlinear.register_operator(model, :my_square, 1, z -> z^2)
+    c = Nonlinear.add_constraint(model, :(my_square($ex)), MOI.LessThan(1.0))
+    @test MOI.is_valid(model, c)
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac])
+    g = fill(NaN, 2)
+    MOI.eval_constraint(d, g, [3.0])
+    @test g == [3.0, 36.0]
+    return
+end
+
+function test_quad_only_with_empty_inner()
+    x = MOI.VariableIndex(1)
+    model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    Nonlinear.add_constraint(
+        model,
+        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
+        MOI.GreaterThan(1.0),
+    )
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    MOI.initialize(d, [:Grad, :Jac, :Hess])
+    g = fill(NaN, 1)
+    MOI.eval_constraint(d, g, [1.5])
+    @test g == [1.5]
+    @test isempty(MOI.hessian_lagrangian_structure(d))
+    @test MOI.NLPBlockData(d).constraint_bounds == [MOI.NLPBoundsPair(1.0, Inf)]
+    return
+end
+
+end  # module
+
+TestNonlinearModelWithQuad.runtests()

From 729379991c0bfde7694c716c38230a0c222d9f68 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 14:45:05 +0000
Subject: [PATCH 2/7] Remove the unused nothing-inner support of ModelWithQuad

No solver passes `inner = nothing`: the consumers all wrap an eager inner
`Nonlinear.Model`, so drop the `set_objective(::Nothing, ::Nothing)` hook
and the docstring sentence advertising it. A solver that manages its own
nonlinear storage can still define `set_objective` for its inner type.
---
 src/Nonlinear/model_with_quad.jl | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/Nonlinear/model_with_quad.jl b/src/Nonlinear/model_with_quad.jl
index 517e2ac4b8..f4c40b3d34 100644
--- a/src/Nonlinear/model_with_quad.jl
+++ b/src/Nonlinear/model_with_quad.jl
@@ -13,9 +13,7 @@
 
 A model layer that stores affine and quadratic objectives and constraints in
 a [`QPBlockData`](@ref) and forwards everything else to the `inner` model,
-typically a [`Model`](@ref). Solvers with their own storage can pass an
-existing `qp` and set `inner` to `nothing` when they manage the nonlinear
-model themselves.
+typically a [`Model`](@ref).
 
 `ModelWithQuad(inner)` and `ModelWithQuad{T}(inner)` create an empty
 [`QPBlockData`](@ref), with `T` defaulting to `Float64`.
@@ -95,9 +93,6 @@ end
 
 # The objective.
 
-# Support for models whose inner model is managed by the solver.
-set_objective(::Nothing, ::Nothing) = nothing
-
 function set_objective(
     model::ModelWithQuad{T},
     obj::Union{

From e655a282df349c1cfdb1f5eae3e7ae6f98555248 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 15:20:31 +0000
Subject: [PATCH 3/7] Remove MOI.set for the objective of ModelWithQuad

Nonlinear.set_objective is the single way to set the objective: it does
everything the MOI.set method did, and it also accepts nothing to clear the
objective, which the MOI attribute interface cannot express. The MOI getters
stay since the function API has no counterpart for them.
---
 src/Nonlinear/model_with_quad.jl       | 22 +++-------------------
 test/Nonlinear/test_model_with_quad.jl |  2 +-
 2 files changed, 4 insertions(+), 20 deletions(-)

diff --git a/src/Nonlinear/model_with_quad.jl b/src/Nonlinear/model_with_quad.jl
index f4c40b3d34..6b704e4196 100644
--- a/src/Nonlinear/model_with_quad.jl
+++ b/src/Nonlinear/model_with_quad.jl
@@ -19,9 +19,9 @@ typically a [`Model`](@ref).
 [`QPBlockData`](@ref), with `T` defaulting to `Float64`.
 
 Add constraints with [`add_constraint`](@ref) or `MOI.add_constraint`, and
-set the objective with [`set_objective`](@ref) or `MOI.set`: affine and
-quadratic functions are routed to the QP block, everything else to the inner
-model. `objective_sink` records where the objective currently lives (`:none`,
+set the objective with [`set_objective`](@ref): affine and quadratic
+functions are routed to the QP block, everything else to the inner model.
+`objective_sink` records where the objective currently lives (`:none`,
 `:quad` or `:inner`).
 
 Create the corresponding evaluator, [`EvaluatorWithQuad`](@ref), with
@@ -115,22 +115,6 @@ function set_objective(model::ModelWithQuad{T}, obj) where {T}
     return
 end
 
-function MOI.set(
-    model::ModelWithQuad{T},
-    ::MOI.ObjectiveFunction{F},
-    func::F,
-) where {
-    T,
-    F<:Union{
-        MOI.VariableIndex,
-        MOI.ScalarAffineFunction{T},
-        MOI.ScalarQuadraticFunction{T},
-    },
-}
-    set_objective(model, func)
-    return
-end
-
 function MOI.get(model::ModelWithQuad, attr::MOI.ObjectiveFunctionType)
     return MOI.get(model.qp, attr)
 end
diff --git a/test/Nonlinear/test_model_with_quad.jl b/test/Nonlinear/test_model_with_quad.jl
index a3a8967dc9..58eb872920 100644
--- a/test/Nonlinear/test_model_with_quad.jl
+++ b/test/Nonlinear/test_model_with_quad.jl
@@ -171,7 +171,7 @@ function test_objective_sink_switching()
         MOI.ScalarAffineTerm{Float64}[],
         0.0,
     )
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    Nonlinear.set_objective(model, f)
     @test model.objective_sink == :quad
     @test MOI.get(model, MOI.ObjectiveFunctionType()) ==
           MOI.ScalarQuadraticFunction{Float64}

From 00f3978b42c993b0dab6face7a4ce51d55653722 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 16:14:31 +0000
Subject: [PATCH 4/7] Move the variables and the parameter convention into
 ModelWithQuad

ModelWithQuad now owns a MOI.Utilities.VariablesContainer and implements
MOI.add_variable, guaranteeing variable indices 1:n like
MOI.Utilities.MatrixOfConstraints, so the EvaluatorWithQuad no longer remaps
the QP block. Parameters are added with MOI.add_constrained_variable: their
indices are offset by _PARAMETER_OFFSET (moved here from Ipopt), QPBlockData
is back to the simple offset-based _is_parameter instead of the parameters
dictionary, and its parameters vector aliases the parameter storage of the
inner Nonlinear.Model, so solvers no longer sync parameter values before a
solve.
---
 src/Nonlinear/model_with_quad.jl       | 212 +++++++++++++++----------
 src/Nonlinear/qp_block_data.jl         | 136 ++++++++--------
 test/Nonlinear/test_model_with_quad.jl |  57 ++++---
 3 files changed, 234 insertions(+), 171 deletions(-)

diff --git a/src/Nonlinear/model_with_quad.jl b/src/Nonlinear/model_with_quad.jl
index 6b704e4196..4295c6c0ad 100644
--- a/src/Nonlinear/model_with_quad.jl
+++ b/src/Nonlinear/model_with_quad.jl
@@ -11,13 +11,20 @@
         objective_sink::Symbol = :none,
     ) where {T,M}
 
-A model layer that stores affine and quadratic objectives and constraints in
-a [`QPBlockData`](@ref) and forwards everything else to the `inner` model,
-typically a [`Model`](@ref).
+A model layer that owns the variables of the model, stores affine and
+quadratic objectives and constraints in a [`QPBlockData`](@ref), and forwards
+everything else to the `inner` model, typically a [`Model`](@ref).
 
 `ModelWithQuad(inner)` and `ModelWithQuad{T}(inner)` create an empty
 [`QPBlockData`](@ref), with `T` defaulting to `Float64`.
 
+Add variables with `MOI.add_variable`: the layer guarantees that the variable
+indices are `1:n`, like `MOI.Utilities.MatrixOfConstraints`. Add parameters
+with `MOI.add_constrained_variable(model, ::MOI.Parameter)`: parameters get
+indices offset by [`_PARAMETER_OFFSET`](@ref), their values are stored in the
+inner model through [`add_parameter`](@ref), and `qp.parameters` aliases that
+storage, so a parameter update is visible to both blocks.
+
 Add constraints with [`add_constraint`](@ref) or `MOI.add_constraint`, and
 set the objective with [`set_objective`](@ref): affine and quadratic
 functions are routed to the QP block, everything else to the inner model.
@@ -25,15 +32,12 @@ functions are routed to the QP block, everything else to the inner model.
 `:quad` or `:inner`).
 
 Create the corresponding evaluator, [`EvaluatorWithQuad`](@ref), with
-`Evaluator(model, backend, ordered_variables)`, or construct it directly from
-an inner `MOI.AbstractNLPEvaluator`. The rows of the QP block come first,
-followed by the rows of the inner evaluator.
-
-Functions added to the QP block may contain parameters, following the
-convention documented in [`QPBlockData`](@ref): a variable is a parameter if
-and only if its index is a key of `qp.parameters`.
+`Evaluator(model, backend)`, or construct it directly from an inner
+`MOI.AbstractNLPEvaluator`. The rows of the QP block come first, followed by
+the rows of the inner evaluator.
 """
 mutable struct ModelWithQuad{T,M}
+    variables::MOI.Utilities.VariablesContainer{T}
     qp::QPBlockData{T}
     inner::M
     objective_sink::Symbol # :none, :quad or :inner
@@ -43,7 +47,14 @@ mutable struct ModelWithQuad{T,M}
         inner::M;
         objective_sink::Symbol = :none,
     ) where {T,M}
-        return new{T,M}(qp, inner, objective_sink)
+        model = new{T,M}(
+            MOI.Utilities.VariablesContainer{T}(),
+            qp,
+            inner,
+            objective_sink,
+        )
+        _share_parameters(model.qp, model.inner)
+        return model
     end
 end
 
@@ -53,6 +64,83 @@ end
 
 ModelWithQuad(inner) = ModelWithQuad{Float64}(inner)
 
+# The QP block reads the parameter values from the inner model's storage.
+_share_parameters(::QPBlockData, ::Any) = nothing
+
+function _share_parameters(qp::QPBlockData{Float64}, inner::Model)
+    qp.parameters = inner.parameters
+    return
+end
+
+# The variables and the parameters.
+
+MOI.add_variable(model::ModelWithQuad) = MOI.add_variable(model.variables)
+
+function MOI.add_constrained_variable(
+    model::ModelWithQuad{T},
+    set::MOI.Parameter{T},
+) where {T}
+    p = add_parameter(model.inner, set.value)
+    x = MOI.VariableIndex(_PARAMETER_OFFSET + p.value)
+    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(x.value)
+    return x, ci
+end
+
+function MOI.is_valid(model::ModelWithQuad, x::MOI.VariableIndex)
+    if _is_parameter(x)
+        return 1 <= x.value - _PARAMETER_OFFSET <= length(model.qp.parameters)
+    end
+    return MOI.is_valid(model.variables, x)
+end
+
+function MOI.is_valid(
+    model::ModelWithQuad{T},
+    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}},
+) where {T}
+    return MOI.is_valid(model, MOI.VariableIndex(ci.value))
+end
+
+function MOI.get(
+    model::ModelWithQuad{T},
+    ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.Parameter{T}},
+) where {T}
+    return length(model.qp.parameters)
+end
+
+function MOI.get(
+    model::ModelWithQuad{T},
+    ::MOI.ListOfConstraintIndices{F,S},
+) where {T,F<:MOI.VariableIndex,S<:MOI.Parameter{T}}
+    n = length(model.qp.parameters)
+    return MOI.ConstraintIndex{F,S}.(_PARAMETER_OFFSET .+ (1:n))
+end
+
+function MOI.get(
+    model::ModelWithQuad{T},
+    ::MOI.ConstraintFunction,
+    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}},
+) where {T}
+    return MOI.VariableIndex(ci.value)
+end
+
+function MOI.get(
+    model::ModelWithQuad{T},
+    ::MOI.ConstraintSet,
+    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}},
+) where {T}
+    return MOI.Parameter(model.qp.parameters[ci.value-_PARAMETER_OFFSET])
+end
+
+function MOI.set(
+    model::ModelWithQuad{T},
+    ::MOI.ConstraintSet,
+    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}},
+    set::MOI.Parameter{T},
+) where {T}
+    model.qp.parameters[ci.value-_PARAMETER_OFFSET] = set.value
+    return
+end
+
 """
     Base.length(model::ModelWithQuad)
 
@@ -204,7 +292,6 @@ end
     EvaluatorWithQuad(
         model::ModelWithQuad,
         inner::MOI.AbstractNLPEvaluator,
-        ordered_variables::Vector{MOI.VariableIndex},
     ) <: MOI.AbstractNLPEvaluator
 
 The evaluator of a [`ModelWithQuad`](@ref) layer. It implements the
@@ -212,50 +299,38 @@ The evaluator of a [`ModelWithQuad`](@ref) layer. It implements the
 first, followed by the rows of `inner`, and the Jacobian and Hessian product
 callbacks compose the contributions of the two blocks.
 
-Create it with `Evaluator(model::ModelWithQuad, backend, ordered_variables)`,
-which recursively creates the evaluator of the inner model, or construct it
-directly from an existing inner evaluator.
+Create it with `Evaluator(model::ModelWithQuad, backend)`, which recursively
+creates the evaluator of the inner model, or construct it directly from an
+existing inner evaluator.
 
-During [`MOI.initialize`](@ref), the variables of the QP block are mapped to
-their consecutive 1-based index in `ordered_variables`; variables absent from
-`ordered_variables` are parameters and keep their index.
+The QP block is evaluated as stored: [`ModelWithQuad`](@ref) owns the
+variables of the model, so their indices are the columns `1:n` and no
+remapping is needed.
 """
 mutable struct EvaluatorWithQuad{T,M,E<:MOI.AbstractNLPEvaluator} <:
                MOI.AbstractNLPEvaluator
     model::ModelWithQuad{T,M}
     inner::E
-    ordered_variables::Vector{MOI.VariableIndex}
-    # A copy of `model.qp` with the variables mapped to their consecutive
-    # 1-based index in `ordered_variables`, and the number of entries of its
-    # Jacobian and of its Hessian of the Lagrangian. Rebuilt during
-    # `MOI.initialize`.
-    qp::QPBlockData{T}
+    # The number of entries of the Jacobian and of the Hessian of the
+    # Lagrangian of the QP block, computed during `MOI.initialize`.
     qp_nnzj::Int
     qp_nnzh::Int
 
     function EvaluatorWithQuad(
         model::ModelWithQuad{T,M},
         inner::E,
-        ordered_variables::Vector{MOI.VariableIndex},
     ) where {T,M,E<:MOI.AbstractNLPEvaluator}
-        return new{T,M,E}(
-            model,
-            inner,
-            ordered_variables,
-            QPBlockData{T}(),
-            0,
-            0,
-        )
+        return new{T,M,E}(model, inner, 0, 0)
     end
 end
 
 function Evaluator(
     model::ModelWithQuad,
     backend::AbstractAutomaticDifferentiation,
-    ordered_variables::Vector{MOI.VariableIndex},
 )
-    inner = Evaluator(model.inner, backend, ordered_variables)
-    return EvaluatorWithQuad(model, inner, ordered_variables)
+    vars = MOI.get(model.variables, MOI.ListOfVariableIndices())
+    inner = Evaluator(model.inner, backend, vars)
+    return EvaluatorWithQuad(model, inner)
 end
 
 function MOI.features_available(d::EvaluatorWithQuad)
@@ -263,34 +338,9 @@ function MOI.features_available(d::EvaluatorWithQuad)
     return filter(f -> f in (:Grad, :Jac, :JacVec, :Hess, :HessVec), features)
 end
 
-function MOI.initialize(
-    d::EvaluatorWithQuad{T},
-    features::Vector{Symbol},
-) where {T}
-    index_map = Dict{MOI.VariableIndex,MOI.VariableIndex}(
-        x => MOI.VariableIndex(i) for (i, x) in enumerate(d.ordered_variables)
-    )
-    # Variables absent from `ordered_variables` are parameters: they keep
-    # their index, which the `parameters` dictionary uses as key.
-    fmap = v::MOI.VariableIndex -> get(index_map, v, v)
-    src = d.model.qp
-    qp = QPBlockData{T}()
-    qp.objective = MOI.Utilities.map_indices(fmap, src.objective)
-    qp.objective_function_type = src.objective_function_type
-    for f in src.constraints
-        push!(qp.constraints, MOI.Utilities.map_indices(fmap, f))
-    end
-    append!(qp.g_L, src.g_L)
-    append!(qp.g_U, src.g_U)
-    append!(qp.mult_g, src.mult_g)
-    append!(qp.function_type, src.function_type)
-    append!(qp.bound_type, src.bound_type)
-    # Alias, do not copy: parameter values updated in the model between
-    # solves must be visible to the evaluator without re-initializing.
-    qp.parameters = src.parameters
-    d.qp = qp
-    d.qp_nnzj = length(MOI.jacobian_structure(qp))
-    d.qp_nnzh = length(MOI.hessian_lagrangian_structure(qp))
+function MOI.initialize(d::EvaluatorWithQuad, features::Vector{Symbol})
+    d.qp_nnzj = length(MOI.jacobian_structure(d.model.qp))
+    d.qp_nnzh = length(MOI.hessian_lagrangian_structure(d.model.qp))
     MOI.initialize(d.inner, features)
     return
 end
@@ -298,7 +348,7 @@ end
 function MOI.eval_objective(d::EvaluatorWithQuad{T}, x) where {T}
     sink = d.model.objective_sink
     if sink == :quad
-        return MOI.eval_objective(d.qp, x)
+        return MOI.eval_objective(d.model.qp, x)
     elseif sink == :inner
         return MOI.eval_objective(d.inner, x)
     else
@@ -309,7 +359,7 @@ end
 function MOI.eval_objective_gradient(d::EvaluatorWithQuad{T}, grad, x) where {T}
     sink = d.model.objective_sink
     if sink == :quad
-        MOI.eval_objective_gradient(d.qp, grad, x)
+        MOI.eval_objective_gradient(d.model.qp, grad, x)
     elseif sink == :inner
         MOI.eval_objective_gradient(d.inner, grad, x)
     else
@@ -319,15 +369,15 @@ function MOI.eval_objective_gradient(d::EvaluatorWithQuad{T}, grad, x) where {T}
 end
 
 function MOI.eval_constraint(d::EvaluatorWithQuad, g, x)
-    m = length(d.qp)
-    MOI.eval_constraint(d.qp, view(g, 1:m), x)
+    m = length(d.model.qp)
+    MOI.eval_constraint(d.model.qp, view(g, 1:m), x)
     MOI.eval_constraint(d.inner, view(g, (m+1):length(g)), x)
     return
 end
 
 function MOI.jacobian_structure(d::EvaluatorWithQuad)
-    J = MOI.jacobian_structure(d.qp)
-    offset = length(d.qp)
+    J = MOI.jacobian_structure(d.model.qp)
+    offset = length(d.model.qp)
     # An evaluator is only required to implement `jacobian_structure` if it
     # supports `:Jac`. If the inner evaluator does not (it then must not have
     # any rows for the stack to be usable), append nothing.
@@ -340,13 +390,13 @@ function MOI.jacobian_structure(d::EvaluatorWithQuad)
 end
 
 function MOI.eval_constraint_jacobian(d::EvaluatorWithQuad, J, x)
-    MOI.eval_constraint_jacobian(d.qp, J, x)
+    MOI.eval_constraint_jacobian(d.model.qp, J, x)
     MOI.eval_constraint_jacobian(d.inner, view(J, (d.qp_nnzj+1):length(J)), x)
     return
 end
 
 function MOI.hessian_lagrangian_structure(d::EvaluatorWithQuad)
-    H = MOI.hessian_lagrangian_structure(d.qp)
+    H = MOI.hessian_lagrangian_structure(d.model.qp)
     if :Hess in MOI.features_available(d.inner)
         append!(H, MOI.hessian_lagrangian_structure(d.inner))
     end
@@ -354,10 +404,10 @@ function MOI.hessian_lagrangian_structure(d::EvaluatorWithQuad)
 end
 
 function MOI.eval_hessian_lagrangian(d::EvaluatorWithQuad, H, x, σ, μ)
-    m = length(d.qp)
-    # If the objective is not in the QP block, `d.qp.objective` is zero, so
+    m = length(d.model.qp)
+    # If the objective is not in the QP block, `d.model.qp.objective` is zero, so
     # passing `σ` is harmless; and vice versa for the inner evaluator.
-    MOI.eval_hessian_lagrangian(d.qp, H, x, σ, view(μ, 1:m))
+    MOI.eval_hessian_lagrangian(d.model.qp, H, x, σ, view(μ, 1:m))
     MOI.eval_hessian_lagrangian(
         d.inner,
         view(H, (d.qp_nnzh+1):length(H)),
@@ -372,14 +422,14 @@ end
 # block write its own rows.
 function MOI.eval_constraint_jacobian_product(d::EvaluatorWithQuad, y, x, w)
     fill!(y, zero(eltype(y)))
-    m = length(d.qp)
+    m = length(d.model.qp)
     MOI.eval_constraint_jacobian_product(
         d.inner,
         view(y, (m+1):length(y)),
         x,
         w,
     )
-    _add_constraint_jacobian_product(d.qp, y, x, w)
+    _add_constraint_jacobian_product(d.model.qp, y, x, w)
     return
 end
 
@@ -393,14 +443,14 @@ function MOI.eval_constraint_jacobian_transpose_product(
     w,
 )
     fill!(y, zero(eltype(y)))
-    m = length(d.qp)
+    m = length(d.model.qp)
     MOI.eval_constraint_jacobian_transpose_product(
         d.inner,
         y,
         x,
         view(w, (m+1):length(w)),
     )
-    _add_constraint_jacobian_transpose_product(d.qp, y, x, view(w, 1:m))
+    _add_constraint_jacobian_transpose_product(d.model.qp, y, x, view(w, 1:m))
     return
 end
 
@@ -413,7 +463,7 @@ function MOI.eval_hessian_lagrangian_product(
     μ,
 )
     fill!(H, zero(eltype(H)))
-    m = length(d.qp)
+    m = length(d.model.qp)
     MOI.eval_hessian_lagrangian_product(
         d.inner,
         H,
@@ -422,7 +472,7 @@ function MOI.eval_hessian_lagrangian_product(
         σ,
         view(μ, (m+1):length(μ)),
     )
-    _add_hessian_lagrangian_product(d.qp, H, x, v, σ, view(μ, 1:m))
+    _add_hessian_lagrangian_product(d.model.qp, H, x, v, σ, view(μ, 1:m))
     return
 end
 
diff --git a/src/Nonlinear/qp_block_data.jl b/src/Nonlinear/qp_block_data.jl
index 22881cef71..1c7b3e65b8 100644
--- a/src/Nonlinear/qp_block_data.jl
+++ b/src/Nonlinear/qp_block_data.jl
@@ -4,12 +4,6 @@
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
 # This file is adapted from `Ipopt.jl/ext/IpoptMathOptInterfaceExt/utils.jl`.
-#
-# Unlike the Ipopt version, a variable is treated as a parameter if and only
-# if its index is a key of the `parameters` dictionary, instead of an
-# index-offset convention. Parameters must therefore be registered in
-# `parameters` before any structure query, but their values may be updated
-# freely between function evaluations.
 
 @enum(
     _FunctionType,
@@ -72,11 +66,12 @@ the solver through the same callbacks as an [`MOI.AbstractNLPEvaluator`](@ref)
 
 ## Parameters
 
-A variable is treated as a parameter if and only if its index is a key of the
-`parameters` dictionary, which maps the raw `MOI.VariableIndex` value of the
-parameter to its current value. Register every parameter in `parameters`
-before querying any structure; the values may be updated freely between
-function evaluations.
+A variable is treated as a parameter if and only if its index is offset by
+[`_PARAMETER_OFFSET`](@ref); see [`_is_parameter`](@ref). The value of the
+parameter `x` is `parameters[x.value - _PARAMETER_OFFSET]`, following the
+indexing of [`ParameterIndex`](@ref), so that `parameters` can alias the
+parameter storage of a [`Model`](@ref). The values may be updated freely
+between function evaluations.
 """
 mutable struct QPBlockData{T}
     objective::Union{MOI.ScalarAffineFunction{T},MOI.ScalarQuadraticFunction{T}}
@@ -89,7 +84,7 @@ mutable struct QPBlockData{T}
     mult_g::Vector{Union{Nothing,T}}
     function_type::Vector{_FunctionType}
     bound_type::Vector{_BoundType}
-    parameters::Dict{Int64,T}
+    parameters::Vector{T}
 
     function QPBlockData{T}() where {T}
         return new(
@@ -101,21 +96,40 @@ mutable struct QPBlockData{T}
             Union{Nothing,T}[],
             _FunctionType[],
             _BoundType[],
-            Dict{Int64,T}(),
+            T[],
         )
     end
 end
 
-_is_parameter(v::MOI.VariableIndex, p::Dict) = haskey(p, v.value)
+"""
+    _PARAMETER_OFFSET
+
+The offset of the `MOI.VariableIndex` value of a parameter: the variable
+`x` is a parameter if and only if `x.value >= _PARAMETER_OFFSET`, and
+`x.value - _PARAMETER_OFFSET` is the value of the corresponding
+[`ParameterIndex`](@ref).
+"""
+const _PARAMETER_OFFSET = 0x00f0000000000000
+
+"""
+    _is_parameter(x::MOI.VariableIndex)
 
-function _value(v::MOI.VariableIndex, x, p::Dict)
-    return _is_parameter(v, p) ? p[v.value] : x[v.value]
+Return whether `x` is a parameter, following the [`_PARAMETER_OFFSET`](@ref)
+convention.
+"""
+_is_parameter(x::MOI.VariableIndex) = x.value >= _PARAMETER_OFFSET
+
+function _value(v::MOI.VariableIndex, x, p::Vector)
+    if _is_parameter(v)
+        return p[v.value-_PARAMETER_OFFSET]
+    end
+    return x[v.value]
 end
 
 function _eval_function(
     f::MOI.ScalarQuadraticFunction{T},
     x::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::T where {T}
     y = f.constant
     for term in f.affine_terms
@@ -136,7 +150,7 @@ end
 function _eval_function(
     f::MOI.ScalarAffineFunction{T},
     x::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::T where {T}
     y = f.constant
     for term in f.terms
@@ -149,20 +163,19 @@ function _eval_dense_gradient(
     ∇f::AbstractVector{T},
     f::MOI.ScalarQuadraticFunction{T},
     x::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Nothing where {T}
     for term in f.affine_terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             ∇f[term.variable.value] += term.coefficient
         end
     end
     for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1, p)
+        if !_is_parameter(term.variable_1)
             v = _value(term.variable_2, x, p)
             ∇f[term.variable_1.value] += term.coefficient * v
         end
-        if term.variable_1 != term.variable_2 &&
-           !_is_parameter(term.variable_2, p)
+        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
             v = _value(term.variable_1, x, p)
             ∇f[term.variable_2.value] += term.coefficient * v
         end
@@ -174,10 +187,10 @@ function _eval_dense_gradient(
     ∇f::AbstractVector{T},
     f::MOI.ScalarAffineFunction{T},
     x::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Nothing where {T}
     for term in f.terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             ∇f[term.variable.value] += term.coefficient
         end
     end
@@ -188,19 +201,18 @@ function _append_sparse_gradient_structure!(
     f::MOI.ScalarQuadraticFunction,
     J,
     row,
-    p::Dict,
+    p::Vector,
 )
     for term in f.affine_terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             push!(J, (row, term.variable.value))
         end
     end
     for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1, p)
+        if !_is_parameter(term.variable_1)
             push!(J, (row, term.variable_1.value))
         end
-        if term.variable_1 != term.variable_2 &&
-           !_is_parameter(term.variable_2, p)
+        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
             push!(J, (row, term.variable_2.value))
         end
     end
@@ -211,10 +223,10 @@ function _append_sparse_gradient_structure!(
     f::MOI.ScalarAffineFunction,
     J,
     row,
-    p::Dict,
+    p::Vector,
 )
     for term in f.terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             push!(J, (row, term.variable.value))
         end
     end
@@ -225,23 +237,22 @@ function _eval_sparse_gradient(
     ∇f::AbstractVector{T},
     f::MOI.ScalarQuadraticFunction{T},
     x::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Int where {T}
     i = 0
     for term in f.affine_terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             i += 1
             ∇f[i] = term.coefficient
         end
     end
     for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1, p)
+        if !_is_parameter(term.variable_1)
             v = _value(term.variable_2, x, p)
             i += 1
             ∇f[i] = term.coefficient * v
         end
-        if term.variable_1 != term.variable_2 &&
-           !_is_parameter(term.variable_2, p)
+        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
             v = _value(term.variable_1, x, p)
             i += 1
             ∇f[i] = term.coefficient * v
@@ -254,11 +265,11 @@ function _eval_sparse_gradient(
     ∇f::AbstractVector{T},
     f::MOI.ScalarAffineFunction{T},
     x::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Int where {T}
     i = 0
     for term in f.terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             i += 1
             ∇f[i] = term.coefficient
         end
@@ -269,11 +280,10 @@ end
 function _append_sparse_hessian_structure!(
     f::MOI.ScalarQuadraticFunction,
     H,
-    p::Dict,
+    p::Vector,
 )
     for term in f.quadratic_terms
-        if _is_parameter(term.variable_1, p) ||
-           _is_parameter(term.variable_2, p)
+        if _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
             continue
         end
         push!(H, (term.variable_1.value, term.variable_2.value))
@@ -284,7 +294,7 @@ end
 function _append_sparse_hessian_structure!(
     ::MOI.ScalarAffineFunction,
     H,
-    ::Dict,
+    ::Vector,
 )
     return nothing
 end
@@ -293,12 +303,11 @@ function _eval_sparse_hessian(
     ∇²f::AbstractVector{T},
     f::MOI.ScalarQuadraticFunction{T},
     σ::T,
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Int where {T}
     i = 0
     for term in f.quadratic_terms
-        if _is_parameter(term.variable_1, p) ||
-           _is_parameter(term.variable_2, p)
+        if _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
             continue
         end
         i += 1
@@ -311,7 +320,7 @@ function _eval_sparse_hessian(
     ∇²f::AbstractVector{T},
     f::MOI.ScalarAffineFunction{T},
     σ::T,
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Int where {T}
     return 0
 end
@@ -588,11 +597,11 @@ function _add_Jv_product(
     y::AbstractVector{T},
     x::AbstractVector{T},
     w::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
     i::Int,
 )::Nothing where {T}
     for term in f.terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             y[i] += term.coefficient * w[term.variable.value]
         end
     end
@@ -604,21 +613,20 @@ function _add_Jv_product(
     y::AbstractVector{T},
     x::AbstractVector{T},
     w::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
     i::Int,
 )::Nothing where {T}
     for term in f.affine_terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             y[i] += term.coefficient * w[term.variable.value]
         end
     end
     for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1, p)
+        if !_is_parameter(term.variable_1)
             v = _value(term.variable_2, x, p)
             y[i] += term.coefficient * v * w[term.variable_1.value]
         end
-        if term.variable_1 != term.variable_2 &&
-           !_is_parameter(term.variable_2, p)
+        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
             v = _value(term.variable_1, x, p)
             y[i] += term.coefficient * v * w[term.variable_2.value]
         end
@@ -631,11 +639,11 @@ function _add_Jtv_product(
     y::AbstractVector{T},
     x::AbstractVector{T},
     w::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
     i::Int,
 )::Nothing where {T}
     for term in f.terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             y[term.variable.value] += term.coefficient * w[i]
         end
     end
@@ -647,21 +655,20 @@ function _add_Jtv_product(
     y::AbstractVector{T},
     x::AbstractVector{T},
     w::AbstractVector{T},
-    p::Dict{Int64,T},
+    p::Vector{T},
     i::Int,
 )::Nothing where {T}
     for term in f.affine_terms
-        if !_is_parameter(term.variable, p)
+        if !_is_parameter(term.variable)
             y[term.variable.value] += term.coefficient * w[i]
         end
     end
     for term in f.quadratic_terms
-        if !_is_parameter(term.variable_1, p)
+        if !_is_parameter(term.variable_1)
             v = _value(term.variable_2, x, p)
             y[term.variable_1.value] += term.coefficient * v * w[i]
         end
-        if term.variable_1 != term.variable_2 &&
-           !_is_parameter(term.variable_2, p)
+        if term.variable_1 != term.variable_2 && !_is_parameter(term.variable_2)
             v = _value(term.variable_1, x, p)
             y[term.variable_2.value] += term.coefficient * v * w[i]
         end
@@ -675,11 +682,10 @@ function _add_Hv_product(
     x::AbstractVector{T},
     v::AbstractVector{T},
     λ::T,
-    p::Dict{Int64,T},
+    p::Vector{T},
 )::Nothing where {T}
     for term in f.quadratic_terms
-        if _is_parameter(term.variable_1, p) ||
-           _is_parameter(term.variable_2, p)
+        if _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
             continue
         end
         i, j = term.variable_1.value, term.variable_2.value
@@ -697,7 +703,7 @@ function _add_Hv_product(
     x::AbstractVector{T},
     v::AbstractVector{T},
     λ::T,
-    p::Dict{Int64,T},
+    p::Vector{T},
 ) where {T}
     return nothing
 end
diff --git a/test/Nonlinear/test_model_with_quad.jl b/test/Nonlinear/test_model_with_quad.jl
index 58eb872920..02bcbaaebf 100644
--- a/test/Nonlinear/test_model_with_quad.jl
+++ b/test/Nonlinear/test_model_with_quad.jl
@@ -27,8 +27,12 @@ end
 #   row 2 (quad layer, quadratic): x^2 + xy + y in [0, 1]
 #   row 3 (inner nlp):             sin(x) <= 0.5
 # and the objective x^2 in the quad layer.
-function _test_model(x, y)
+function _test_model()
     model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    x = MOI.add_variable(model)
+    y = MOI.add_variable(model)
+    @test (x, y) == (MOI.VariableIndex(1), MOI.VariableIndex(2))
+    @test MOI.is_valid(model, x) && !MOI.is_valid(model, MOI.VariableIndex(3))
     Nonlinear.set_objective(
         model,
         MOI.ScalarQuadraticFunction(
@@ -64,15 +68,12 @@ function _test_model(x, y)
     c3 = Nonlinear.add_constraint(model, :(sin($x)), MOI.LessThan(0.5))
     @test c3 isa Nonlinear.ConstraintIndex
     @test length(model) == 2
-    return model
+    return model, x, y
 end
 
 function test_evaluator_with_quad()
-    # Use non-consecutive variable indices to test that the evaluator remaps
-    # the variables to their consecutive index in `ordered_variables`.
-    x, y = MOI.VariableIndex(5), MOI.VariableIndex(9)
-    model = _test_model(x, y)
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
+    model, x, y = _test_model()
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     @test d isa Nonlinear.EvaluatorWithQuad
     @test d.inner isa Nonlinear.Evaluator
     @test MOI.features_available(d) == [:Grad, :Jac, :JacVec, :Hess, :HessVec]
@@ -122,9 +123,8 @@ function test_evaluator_with_quad()
 end
 
 function test_evaluator_products()
-    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
-    model = _test_model(x, y)
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x, y])
+    model, x, y = _test_model()
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac, :JacVec, :Hess, :HessVec])
     xv = [1.0, 2.0]
     # Dense Jacobian from the sparse callback, as the reference.
@@ -163,8 +163,8 @@ function test_evaluator_products()
 end
 
 function test_objective_sink_switching()
-    x = MOI.VariableIndex(1)
     model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    x = MOI.add_variable(model)
     @test model.objective_sink == :none
     f = MOI.ScalarQuadraticFunction(
         [MOI.ScalarQuadraticTerm(2.0, x, x)],
@@ -176,7 +176,7 @@ function test_objective_sink_switching()
     @test MOI.get(model, MOI.ObjectiveFunctionType()) ==
           MOI.ScalarQuadraticFunction{Float64}
     @test MOI.get(model, MOI.ObjectiveFunction{typeof(f)}()) ≈ f
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(d, [3.0]) == 9.0
     @test MOI.NLPBlockData(d).has_objective
@@ -184,7 +184,7 @@ function test_objective_sink_switching()
     # cleared, including its Hessian entries.
     Nonlinear.set_objective(model, :(sin($x)))
     @test model.objective_sink == :inner
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac, :Hess])
     @test MOI.eval_objective(d, [3.0]) == sin(3.0)
     @test MOI.NLPBlockData(d).has_objective
@@ -196,12 +196,12 @@ function test_objective_sink_switching()
     g = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0)
     Nonlinear.set_objective(model, g)
     @test model.objective_sink == :quad
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac])
     @test MOI.eval_objective(d, [3.0]) == 7.0
     Nonlinear.set_objective(model, nothing)
     @test model.objective_sink == :none
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac])
     @test MOI.eval_objective(d, [3.0]) == 0.0
     grad = fill(NaN, 1)
@@ -212,11 +212,18 @@ function test_objective_sink_switching()
 end
 
 function test_quad_parameters()
-    x, p = MOI.VariableIndex(1), MOI.VariableIndex(42)
     model = Nonlinear.ModelWithQuad(Nonlinear.Model())
-    # Register the parameter before any structure query. Its value may be
-    # updated between evaluations.
-    model.qp.parameters[p.value] = 5.0
+    x = MOI.add_variable(model)
+    p, cp = MOI.add_constrained_variable(model, MOI.Parameter(5.0))
+    @test p.value == Nonlinear._PARAMETER_OFFSET + 1
+    @test MOI.is_valid(model, p) && MOI.is_valid(model, cp)
+    @test MOI.get(model, MOI.ConstraintFunction(), cp) == p
+    @test MOI.get(model, MOI.ConstraintSet(), cp) == MOI.Parameter(5.0)
+    F, S = MOI.VariableIndex, MOI.Parameter{Float64}
+    @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 1
+    @test MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) == [cp]
+    # The value is stored in the inner model, aliased by the QP block.
+    @test model.qp.parameters === model.inner.parameters
     Nonlinear.add_constraint(
         model,
         MOI.ScalarAffineFunction(
@@ -234,7 +241,7 @@ function test_quad_parameters()
         ),
         MOI.LessThan(10.0),
     )
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac, :Hess])
     g = fill(NaN, 2)
     MOI.eval_constraint(d, g, [1.0])
@@ -246,15 +253,15 @@ function test_quad_parameters()
     @test J == [2.0, 5.0]
     @test isempty(MOI.hessian_lagrangian_structure(d))
     # Updating the parameter value must be visible without re-initializing.
-    model.qp.parameters[p.value] = 7.0
+    MOI.set(model, MOI.ConstraintSet(), cp, MOI.Parameter(7.0))
     MOI.eval_constraint(d, g, [1.0])
     @test g == [2.0 * 1.0 + 3.0 * 7.0, 7.0 * 1.0]
     return
 end
 
 function test_attribute_forwarding()
-    x = MOI.VariableIndex(1)
     model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    x = MOI.add_variable(model)
     F, S = MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}
     f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
     ci = MOI.add_constraint(model, f, MOI.GreaterThan(1.0))
@@ -278,7 +285,7 @@ function test_attribute_forwarding()
     Nonlinear.register_operator(model, :my_square, 1, z -> z^2)
     c = Nonlinear.add_constraint(model, :(my_square($ex)), MOI.LessThan(1.0))
     @test MOI.is_valid(model, c)
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac])
     g = fill(NaN, 2)
     MOI.eval_constraint(d, g, [3.0])
@@ -287,14 +294,14 @@ function test_attribute_forwarding()
 end
 
 function test_quad_only_with_empty_inner()
-    x = MOI.VariableIndex(1)
     model = Nonlinear.ModelWithQuad(Nonlinear.Model())
+    x = MOI.add_variable(model)
     Nonlinear.add_constraint(
         model,
         MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
         MOI.GreaterThan(1.0),
     )
-    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode(), [x])
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
     MOI.initialize(d, [:Grad, :Jac, :Hess])
     g = fill(NaN, 1)
     MOI.eval_constraint(d, g, [1.5])

From 12bba957ab4e26d21d1ebbdd11729f6e67dbe086 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 16:41:48 +0000
Subject: [PATCH 5/7] Always alias the parameter storage of the inner model

Sharing only when the inner model is a Nonlinear.Model silently left
qp.parameters empty for any other inner model type. Assume instead that the
inner model exposes its parameter values as parameters::Vector{T}, like
Nonlinear.Model does, and always alias it; an inner model without that field
fails loudly at construction.
---
 src/Nonlinear/model_with_quad.jl | 20 ++++++++------------
 1 file changed, 8 insertions(+), 12 deletions(-)

diff --git a/src/Nonlinear/model_with_quad.jl b/src/Nonlinear/model_with_quad.jl
index 4295c6c0ad..4b2e034907 100644
--- a/src/Nonlinear/model_with_quad.jl
+++ b/src/Nonlinear/model_with_quad.jl
@@ -21,9 +21,10 @@ everything else to the `inner` model, typically a [`Model`](@ref).
 Add variables with `MOI.add_variable`: the layer guarantees that the variable
 indices are `1:n`, like `MOI.Utilities.MatrixOfConstraints`. Add parameters
 with `MOI.add_constrained_variable(model, ::MOI.Parameter)`: parameters get
-indices offset by [`_PARAMETER_OFFSET`](@ref), their values are stored in the
-inner model through [`add_parameter`](@ref), and `qp.parameters` aliases that
-storage, so a parameter update is visible to both blocks.
+indices offset by [`_PARAMETER_OFFSET`](@ref), and their values are stored in
+the inner model through [`add_parameter`](@ref). The inner model must expose
+that storage as `parameters::Vector{T}`, like [`Model`](@ref) does:
+`qp.parameters` aliases it, so a parameter update is visible to both blocks.
 
 Add constraints with [`add_constraint`](@ref) or `MOI.add_constraint`, and
 set the objective with [`set_objective`](@ref): affine and quadratic
@@ -53,7 +54,10 @@ mutable struct ModelWithQuad{T,M}
             inner,
             objective_sink,
         )
-        _share_parameters(model.qp, model.inner)
+        # The QP block reads the parameter values from the storage of the
+        # inner model, which must expose them as `parameters::Vector{T}`,
+        # like [`Model`](@ref) does.
+        model.qp.parameters = inner.parameters
         return model
     end
 end
@@ -64,14 +68,6 @@ end
 
 ModelWithQuad(inner) = ModelWithQuad{Float64}(inner)
 
-# The QP block reads the parameter values from the inner model's storage.
-_share_parameters(::QPBlockData, ::Any) = nothing
-
-function _share_parameters(qp::QPBlockData{Float64}, inner::Model)
-    qp.parameters = inner.parameters
-    return
-end
-
 # The variables and the parameters.
 
 MOI.add_variable(model::ModelWithQuad) = MOI.add_variable(model.variables)

From 1df877d9bf1c35527e3e27536c5507b29635a017 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 16:47:40 +0000
Subject: [PATCH 6/7] Define _is_parameter on the affine and quadratic terms

They belong next to the convention; Ipopt and MadNLP defined them on their
local alias of the function, which pirates it.
---
 src/Nonlinear/qp_block_data.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/Nonlinear/qp_block_data.jl b/src/Nonlinear/qp_block_data.jl
index 1c7b3e65b8..a7c5c7ebe7 100644
--- a/src/Nonlinear/qp_block_data.jl
+++ b/src/Nonlinear/qp_block_data.jl
@@ -119,6 +119,12 @@ convention.
 """
 _is_parameter(x::MOI.VariableIndex) = x.value >= _PARAMETER_OFFSET
 
+_is_parameter(term::MOI.ScalarAffineTerm) = _is_parameter(term.variable)
+
+function _is_parameter(term::MOI.ScalarQuadraticTerm)
+    return _is_parameter(term.variable_1) || _is_parameter(term.variable_2)
+end
+
 function _value(v::MOI.VariableIndex, x, p::Vector)
     if _is_parameter(v)
         return p[v.value-_PARAMETER_OFFSET]

From f04f2c4c3c3ac5e0fb1308c0b9cc47f4c52bfe81 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Sat, 15 Aug 2026 17:26:51 +0000
Subject: [PATCH 7/7] Move the parameter substitution into ModelWithQuad

The substitution of the offset parameter indices by ParameterIndex before
the inner model parses a function is index arithmetic tied to the parameter
convention of the layer, so it belongs here: the generic add_constraint and
set_objective now perform it, and the solvers just forward.
---
 src/Nonlinear/model_with_quad.jl       | 42 ++++++++++++++++++++++++++
 test/Nonlinear/test_model_with_quad.jl | 14 +++++++++
 2 files changed, 56 insertions(+)

diff --git a/src/Nonlinear/model_with_quad.jl b/src/Nonlinear/model_with_quad.jl
index 4b2e034907..8176616abf 100644
--- a/src/Nonlinear/model_with_quad.jl
+++ b/src/Nonlinear/model_with_quad.jl
@@ -145,6 +145,42 @@ before the rows of the inner model in the corresponding evaluator.
 """
 Base.length(model::ModelWithQuad) = length(model.qp)
 
+# Replace the parameters of `f`, encoded as `MOI.VariableIndex`es offset by
+# [`_PARAMETER_OFFSET`](@ref), by the corresponding [`ParameterIndex`](@ref),
+# which the inner model understands. An affine or quadratic function that
+# contains a parameter is converted to `MOI.ScalarNonlinearFunction`, because
+# the inner model parses such functions with their variable indices verbatim.
+_replace_parameters(f) = f
+
+function _replace_parameters(f::MOI.VariableIndex)
+    if _is_parameter(f)
+        return ParameterIndex(f.value - _PARAMETER_OFFSET)
+    end
+    return f
+end
+
+function _replace_parameters(f::MOI.ScalarAffineFunction)
+    if any(_is_parameter, f.terms)
+        return _replace_parameters(convert(MOI.ScalarNonlinearFunction, f))
+    end
+    return f
+end
+
+function _replace_parameters(f::MOI.ScalarQuadraticFunction)
+    if any(_is_parameter, f.affine_terms) ||
+       any(_is_parameter, f.quadratic_terms)
+        return _replace_parameters(convert(MOI.ScalarNonlinearFunction, f))
+    end
+    return f
+end
+
+function _replace_parameters(f::MOI.ScalarNonlinearFunction)
+    for (i, arg) in enumerate(f.args)
+        f.args[i] = _replace_parameters(arg)
+    end
+    return f
+end
+
 # Methods forwarded to the inner model.
 
 function add_parameter(model::ModelWithQuad, value::Real)
@@ -194,6 +230,9 @@ end
 function set_objective(model::ModelWithQuad{T}, obj) where {T}
     F = MOI.ScalarAffineFunction{T}
     MOI.set(model.qp, MOI.ObjectiveFunction{F}(), zero(F))
+    if !isempty(model.qp.parameters)
+        obj = _replace_parameters(obj)
+    end
     set_objective(model.inner, obj)
     model.objective_sink = obj === nothing ? :none : :inner
     return
@@ -225,6 +264,9 @@ function add_constraint(
 end
 
 function add_constraint(model::ModelWithQuad, func, set)
+    if !isempty(model.qp.parameters)
+        func = _replace_parameters(func)
+    end
     return add_constraint(model.inner, func, set)
 end
 
diff --git a/test/Nonlinear/test_model_with_quad.jl b/test/Nonlinear/test_model_with_quad.jl
index 02bcbaaebf..d061a1b1c8 100644
--- a/test/Nonlinear/test_model_with_quad.jl
+++ b/test/Nonlinear/test_model_with_quad.jl
@@ -256,6 +256,20 @@ function test_quad_parameters()
     MOI.set(model, MOI.ConstraintSet(), cp, MOI.Parameter(7.0))
     MOI.eval_constraint(d, g, [1.0])
     @test g == [2.0 * 1.0 + 3.0 * 7.0, 7.0 * 1.0]
+    # A nonlinear constraint with the parameter in an embedded affine
+    # subfunction: the layer substitutes the parameter before the inner model
+    # parses the function.
+    aff = MOI.ScalarAffineFunction(
+        [MOI.ScalarAffineTerm(3.0, p), MOI.ScalarAffineTerm(1.0, x)],
+        0.0,
+    )
+    snf = MOI.ScalarNonlinearFunction(:sqrt, Any[aff])
+    Nonlinear.add_constraint(model, snf, MOI.LessThan(10.0))
+    d = Nonlinear.Evaluator(model, Nonlinear.SparseReverseMode())
+    MOI.initialize(d, [:Grad, :Jac])
+    g = fill(NaN, 3)
+    MOI.eval_constraint(d, g, [1.0])
+    @test g ≈ [2.0 + 3.0 * 7.0, 7.0, sqrt(3.0 * 7.0 + 1.0)]
     return
 end