diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
index 48913f5..9127294 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/tutorial.md
@@ -23,7 +23,7 @@ Let's define the famous Rosenbrock function
 ```math
 f(x) = (x_1 - 1)^2 + 100(x_2 - x_1^2)^2,
 ```
-with starting point ``x^0 = (-1.2,1.0)``.
+with starting point ``x_0 = (-1.2, 1.0)``.
 
 ```@example jumpnlp
 using NLPModels, NLPModelsJuMP, JuMP
@@ -31,7 +31,7 @@ using NLPModels, NLPModelsJuMP, JuMP
 x0 = [-1.2; 1.0]
 model = Model() # No solver is required
 @variable(model, x[i=1:2], start=x0[i])
-@NLobjective(model, Min, (x[1] - 1)^2 + 100 * (x[2] - x[1]^2)^2)
+@objective(model, Min, (x[1] - 1)^2 + 100 * (x[2] - x[1]^2)^2)
 
 nlp = MathOptNLPModel(model)
 ```
@@ -148,9 +148,9 @@ using NLPModels, NLPModelsJuMP, JuMP
 model = Model()
 x0 = [-1.2; 1.0]
 @variable(model, x[i=1:2] >= 0.0, start=x0[i])
-@NLobjective(model, Min, (x[1] - 1)^2 + 100 * (x[2] - x[1]^2)^2)
+@objective(model, Min, (x[1] - 1)^2 + 100 * (x[2] - x[1]^2)^2)
 @constraint(model, x[1] + x[2] == 3.0)
-@NLconstraint(model, x[1] * x[2] >= 1.0)
+@constraint(model, x[1] * x[2] >= 1.0)
 
 nlp = MathOptNLPModel(model)
 
@@ -184,7 +184,7 @@ using NLPModels, NLPModelsJuMP, JuMP
 model = Model()
 x0 = [-1.2; 1.0]
 @variable(model, x[i=1:2], start=x0[i])
-@NLexpression(model, F1, x[1] - 1)
+@expression(model, F1, x[1] - 1)
 @NLexpression(model, F2, 10 * (x[2] - x[1]^2))
 
 nls = MathOptNLSModel(model, [F1, F2], name="rosen-nls")
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index 9ec6d2c..199e79e 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -86,6 +86,12 @@ MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{<:ALS}) = true
 MOI.supports_constraint(::Optimizer, ::Type{VQF}, ::Type{<:VLS}) = true
 MOI.supports_constraint(::Optimizer, ::Type{SNF}, ::Type{<:ALS}) = true
 
+function _add_identity_constraints!(index_map, src, F, S)
+  for ci in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
+    index_map[ci] = ci
+  end
+end
+
 function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
   if !haskey(dest.options, "solver")
     error(
@@ -94,6 +100,10 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
   end
   dest.nlp, index_map = nlp_model(src)
   dest.solver = dest.options["solver"](dest.nlp)
+  # Add identity constraint mappings to the index map
+  for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
+    _add_identity_constraints!(index_map, src, F, S)
+  end
   return index_map
 end
 
diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl
index a4975aa..11168d0 100644
--- a/src/moi_nlp_model.jl
+++ b/src/moi_nlp_model.jl
@@ -3,6 +3,8 @@ export MathOptNLPModel
 mutable struct MathOptNLPModel <: AbstractNLPModel{Float64, Vector{Float64}}
   meta::NLPModelMeta{Float64, Vector{Float64}}
   eval::MOI.Nonlinear.Evaluator
+  jump_variables::Dict{String,Int}
+  jump_constraints::Dict{String,Int}
   lincon::LinearConstraints
   quadcon::QuadraticConstraints
   nlcon::NonLinearStructure
@@ -26,15 +28,15 @@ function MathOptNLPModel(jmodel::JuMP.Model; kws...)
 end
 
 function MathOptNLPModel(moimodel::MOI.ModelLike; kws...)
-  return nlp_model(moimodel; kws...)[1]
+  nlp, _ = nlp_model(moimodel; kws...)
+  return nlp
 end
 
 function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String = "Generic")
-  index_map, nvar, lvar, uvar, x0 = parser_variables(moimodel)
-  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon =
-    parser_MOI(moimodel, index_map, nvar)
+  jump_variables, variables, nvar, lvar, uvar, x0, index_map = parser_variables(moimodel)
+  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon, jump_constraints_linear, jump_constraints_quadratic, valid_label = parser_MOI(moimodel, variables, index_map)
 
-  nlp_data = _nlp_block(moimodel)
+  nlp_data, valid_label2, jump_constraints_nonlinear = _nlp_block(moimodel)
   nlcon = parser_NL(nlp_data, hessian = hessian)
   oracles = parser_oracles(moimodel)
   counters = Counters()
@@ -44,7 +46,7 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
   if nlp_data.has_objective
     obj = Objective("NONLINEAR", 0.0, spzeros(Float64, nvar), COO(), 0)
   else
-    obj = parser_objective_MOI(moimodel, nvar, index_map)
+    obj = parser_objective_MOI(moimodel, variables, index_map)
   end
 
   # Total counts
@@ -74,9 +76,26 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     hess_available = hessian && oracles.hessian_oracles_supported,
   )
 
-  return MathOptNLPModel(
+  # Label of the constraints
+  jump_constraints = Dict{String, Int}()
+  if valid_label && valid_label2 && (oracles.ncon == 0)
+    sizehint!(jump_constraints, ncon)
+    for (key, val) in jump_constraints_linear
+      jump_constraints[key] = val
+    end
+    for (key, val) in jump_constraints_quadratic
+      jump_constraints[key] = val + nlin
+    end
+    for (key, val) in jump_constraints_nonlinear
+      jump_constraints[key] = val + nlin + quadcon.nquad
+    end
+  end
+
+  nlp = MathOptNLPModel(
     meta,
     nlp_data.evaluator,
+    jump_variables,
+    jump_constraints,
     lincon,
     quadcon,
     nlcon,
@@ -85,8 +104,8 @@ function nlp_model(moimodel::MOI.ModelLike; hessian::Bool = true, name::String =
     hv,
     obj,
     counters,
-  ),
-  index_map
+  )
+  return nlp, index_map
 end
 
 function NLPModels.obj(nlp::MathOptNLPModel, x::AbstractVector)
diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl
index eb11cb6..e96823e 100644
--- a/src/moi_nls_model.jl
+++ b/src/moi_nls_model.jl
@@ -5,6 +5,8 @@ mutable struct MathOptNLSModel <: AbstractNLSModel{Float64, Vector{Float64}}
   nls_meta::NLSMeta{Float64, Vector{Float64}}
   Feval::MOI.Nonlinear.Evaluator
   ceval::MOI.Nonlinear.Evaluator
+  jump_variables::Dict{String,Int}
+  jump_constraints::Dict{String,Int}
   lls::Objective
   linequ::LinearEquations
   nlequ::NonLinearStructure
@@ -26,17 +28,16 @@ Construct a `MathOptNLSModel` from a `JuMP` model and a container of JuMP
 """
 function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::String = "Generic")
   moimodel = backend(cmodel)
-  index_map, nvar, lvar, uvar, x0 = parser_variables(moimodel)
+  jump_variables, variables, nvar, lvar, uvar, x0, index_map = parser_variables(moimodel)
 
-  lls, linequ, nlinequ = parser_linear_expression(cmodel, nvar, index_map, F)
-  Feval, nlequ, nnlnequ = parser_nonlinear_expression(cmodel, nvar, F, hessian = hessian)
+  lls, linequ, nlinequ = parser_linear_expression(cmodel, variables, F)
+  Feval, nlequ, nnlnequ = parser_nonlinear_expression(cmodel, variables, F, hessian = hessian)
 
   _nlp_sync!(cmodel)
   moimodel = backend(cmodel)
-  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon =
-    parser_MOI(moimodel, index_map, nvar)
+  nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon, jump_constraints_linear, jump_constraints_quadratic, valid_label = parser_MOI(moimodel, variables, index_map)
 
-  nlp_data = _nlp_block(moimodel)
+  nlp_data, valid_label2, jump_constraints_nonlinear = _nlp_block(moimodel)
   nlcon = parser_NL(nlp_data, hessian = hessian)
   oracles = parser_oracles(moimodel)
   nls_counters = NLSCounters()
@@ -75,11 +76,28 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
 
   nls_meta = NLSMeta(nequ, nvar, nnzj = Fnnzj, nnzh = Fnnzh, lin = collect(1:nlinequ))
 
+  # Label of the constraints
+  jump_constraints = Dict{String, Int}()
+  if valid_label && valid_label2 && (oracles.ncon == 0)
+    sizehint!(jump_constraints, ncon)
+    for (key, val) in jump_constraints_linear
+      jump_constraints[key] = val
+    end
+    for (key, val) in jump_constraints_quadratic
+      jump_constraints[key] = val + nlin
+    end
+    for (key, val) in jump_constraints_nonlinear
+      jump_constraints[key] = val + nlin + quadcon.nquad
+    end
+  end
+
   return MathOptNLSModel(
     meta,
     nls_meta,
     Feval,
     nlp_data.evaluator,
+    jump_variables,
+    jump_constraints,
     lls,
     linequ,
     nlequ,
diff --git a/src/utils.jl b/src/utils.jl
index 9013716..599c68f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -252,7 +252,7 @@ function coo_sym_dot(
 end
 
 """
-    parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon, index_map)
+    parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
 
 Parse a `ScalarAffineFunction` fun with its associated set.
 `linrows`, `lincols`, `linvals`, `lin_lcon` and `lin_ucon` are updated.
@@ -285,7 +285,7 @@ function parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_uco
 end
 
 """
-    parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon, index_map)
+    parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon)
 
 Parse a `VectorAffineFunction` fun with its associated set.
 `linrows`, `lincols`, `linvals`, `lin_lcon` and `lin_ucon` are updated.
@@ -314,7 +314,7 @@ function parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_uco
 end
 
 """
-    parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+    parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon)
 
 Parse a `ScalarQuadraticFunction` fun with its associated set.
 `qcons`, `quad_lcon`, `quad_ucon` are updated.
@@ -376,7 +376,7 @@ function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
 end
 
 """
-    parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+    parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon)
 
 Parse a `VectorQuadraticFunction` fun with its associated set.
 `qcons`, `quad_lcon`, `quad_ucon` are updated.
@@ -443,11 +443,19 @@ function parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
 end
 
 """
-    parser_MOI(moimodel, index_map, nvar)
+    parser_MOI(moimodel, variables)
 
 Parse linear constraints of a `MOI.ModelLike`.
 """
-function parser_MOI(moimodel, index_map, nvar)
+function parser_MOI(moimodel, variables, index_map)
+
+  # Number of variables
+  nvar = length(variables)
+
+  # Ensure that each constraint has a valid label
+  valid_label = true
+  jump_constraints_linear = Dict{String, Int}()
+  jump_constraints_quadratic = Dict{String, Int}()
 
   # Variables associated to linear constraints
   nlin = 0
@@ -475,31 +483,39 @@ function parser_MOI(moimodel, index_map, nvar)
     F <: AF || F <: QF || F == SNF || F == VI || error("Function $F is not supported.")
     S <: LS || error("Set $S is not supported.")
 
+    (F == VI) && continue
     conindices = MOI.get(moimodel, MOI.ListOfConstraintIndices{F, S}())
     for cidx in conindices
+      cname = MOI.get(moimodel, MOI.ConstraintName(), cidx)
       fun = MOI.get(moimodel, MOI.ConstraintFunction(), cidx)
-      if F == VI
-        index_map[cidx] = MOI.ConstraintIndex{F, S}(fun.value)
-        continue
-      else
-        index_map[cidx] = MOI.ConstraintIndex{F, S}(nlin)
-      end
       set = MOI.get(moimodel, MOI.ConstraintSet(), cidx)
       if typeof(fun) <: SAF
         parser_SAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon, index_map)
         nlin += 1
+        if valid_label && (cname != "")
+          jump_constraints_linear[cname] = nlin
+        else
+          valid_label = false
+        end
       end
       if typeof(fun) <: VAF
         parser_VAF(fun, set, linrows, lincols, linvals, nlin, lin_lcon, lin_ucon, index_map)
         nlin += set.dimension
+        valid_label = false
       end
       if typeof(fun) <: SQF
         parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
         nquad += 1
+        if valid_label && (cname != "")
+          jump_constraints_quadratic[cname] = nquad
+        else
+          valid_label = false
+        end
       end
       if typeof(fun) <: VQF
         parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
         nquad += set.dimension
+        valid_label = false
       end
     end
   end
@@ -514,15 +530,28 @@ function parser_MOI(moimodel, index_map, nvar)
   end
   quadcon = QuadraticConstraints(nquad, qcons, quad_nnzj, quad_nnzh)
 
-  return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon
+  return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon, jump_constraints_linear, jump_constraints_quadratic, valid_label
 end
 
 # Affine or quadratic, nothing to do
-_nlp_model(::MOI.Nonlinear.Model, ::MOI.ModelLike, ::Type, ::Type) = false
+function _nlp_model(::MOI.Nonlinear.Model, ::MOI.ModelLike, ::Dict{String, Int}, ::Type, ::Type)
+  has_nonlinear = false
+  valid_label = true
+  return has_nonlinear, valid_label
+end
 
-function _nlp_model(dest::MOI.Nonlinear.Model, src::MOI.ModelLike, F::Type{SNF}, S::Type)
+function _nlp_model(dest::MOI.Nonlinear.Model, src::MOI.ModelLike, jump_constraints_nonlinear::Dict{String, Int}, F::Type{SNF}, S::Type)
   has_nonlinear = false
+  valid_label = true
+  ncon = 0
   for ci in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
+    cname = MOI.get(src, MOI.ConstraintName(), ci)
+    ncon += 1
+    if valid_label && (cname != "")
+      jump_constraints_nonlinear[cname] = ncon
+    else
+      valid_label = false
+    end
     MOI.Nonlinear.add_constraint(
       dest,
       MOI.get(src, MOI.ConstraintFunction(), ci),
@@ -530,12 +559,14 @@ function _nlp_model(dest::MOI.Nonlinear.Model, src::MOI.ModelLike, F::Type{SNF},
     )
     has_nonlinear = true
   end
-  return has_nonlinear
+  return has_nonlinear, valid_label
 end
 
-function _nlp_model(model::MOI.ModelLike)::Union{Nothing, MOI.Nonlinear.Model}
+function _nlp_model(model::MOI.ModelLike)
   nlp_model = MOI.Nonlinear.Model()
   has_nonlinear = false
+  valid_label = true
+  jump_constraints_nonlinear = Dict{String, Int}()
   for attr in MOI.get(model, MOI.ListOfModelAttributesSet())
     if attr isa MOI.UserDefinedFunction
       has_nonlinear = true
@@ -544,7 +575,9 @@ function _nlp_model(model::MOI.ModelLike)::Union{Nothing, MOI.Nonlinear.Model}
     end
   end
   for (F, S) in MOI.get(model, MOI.ListOfConstraintTypesPresent())
-    has_nonlinear |= _nlp_model(nlp_model, model, F, S)
+    has_nonlinear_constraint, valid_label_constraint = _nlp_model(nlp_model, model, jump_constraints_nonlinear, F, S)
+    has_nonlinear |= has_nonlinear_constraint
+    valid_label = valid_label && valid_label_constraint
   end
   F = MOI.get(model, MOI.ObjectiveFunctionType())
   if F <: SNF
@@ -552,16 +585,16 @@ function _nlp_model(model::MOI.ModelLike)::Union{Nothing, MOI.Nonlinear.Model}
     has_nonlinear = true
   end
   if !has_nonlinear
-    return nothing
+    return nothing, valid_label, jump_constraints_nonlinear
   end
-  return nlp_model
+  return nlp_model, valid_label, jump_constraints_nonlinear
 end
 
 function _nlp_block(model::MOI.ModelLike)
   # Old interface with `@NL...`
   nlp_data = MOI.get(model, MOI.NLPBlock())
   # New interface with `@constraint` and `@objective`
-  nlp_model = _nlp_model(model)
+  nlp_model, valid_label, jump_constraints_nonlinear = _nlp_model(model)
   vars = MOI.get(model, MOI.ListOfVariableIndices())
   if isnothing(nlp_data)
     if isnothing(nlp_model)
@@ -574,6 +607,7 @@ function _nlp_block(model::MOI.ModelLike)
       nlp_data = MOI.NLPBlockData(evaluator)
     end
   else
+    valid_label = false
     if !isnothing(nlp_model)
       error(
         "Cannot optimize a model which contains the features from " *
@@ -583,7 +617,7 @@ function _nlp_block(model::MOI.ModelLike)
       )
     end
   end
-  return nlp_data
+  return nlp_data, valid_label, jump_constraints_nonlinear
 end
 
 """
@@ -633,7 +667,7 @@ function parser_oracles(moimodel)
     MOI.ListOfConstraintIndices{MOI.VectorOfVariables, MOI.VectorNonlinearOracle{Float64}}(),
   )
     f = MOI.get(moimodel, MOI.ConstraintFunction(), ci)  # ::MOI.VectorOfVariables
-    set = MOI.get(moimodel, MOI.ConstraintSet(), ci)       # ::MOI.VectorNonlinearOracle{Float64}
+    set = MOI.get(moimodel, MOI.ConstraintSet(), ci)     # ::MOI.VectorNonlinearOracle{Float64}
 
     cache = _VectorNonlinearOracleCache(set)
     push!(oracles, (f, cache))
@@ -668,40 +702,49 @@ Parse variables informations of a `MOI.ModelLike`.
 """
 function parser_variables(model::MOI.ModelLike)
   # Number of variables and bounds constraints
-  vars = MOI.get(model, MOI.ListOfVariableIndices())
-  nvar = length(vars)
+  variables = MOI.get(model, MOI.ListOfVariableIndices())
+  nvar = length(variables)
   lvar = zeros(nvar)
   uvar = zeros(nvar)
+
   # Initial solution
   x0 = zeros(nvar)
   has_start = MOI.VariablePrimalStart() in MOI.get(model, MOI.ListOfVariableAttributesSet())
 
+  jump_variables = Dict{String,Int}()
+  sizehint!(jump_variables, nvar)
+
+  # Build an IndexMap from source MOI.VariableIndex to contiguous 1:n.
+  # After variable deletion, vi.value may not be contiguous.
   index_map = MOI.Utilities.IndexMap()
-  for (i, vi) in enumerate(vars)
-    index_map[vi] = MOI.VariableIndex(i)
-  end
+  for (idx, vi) in enumerate(variables)
+    index_map[vi] = MOI.VariableIndex(idx)
+    name = MOI.get(model, MOI.VariableName(), vi)
+    jump_variables[name] = idx
 
-  for (i, vi) in enumerate(vars)
-    lvar[i], uvar[i] = MOI.Utilities.get_bounds(model, Float64, vi)
+    lvar[idx], uvar[idx] = MOI.Utilities.get_bounds(model, Float64, vi)
     if has_start
       val = MOI.get(model, MOI.VariablePrimalStart(), vi)
       if val !== nothing
-        x0[i] = val
+        x0[idx] = val
       end
     end
   end
 
-  return index_map, nvar, lvar, uvar, x0
+  return jump_variables, variables, nvar, lvar, uvar, x0, index_map
 end
 
 """
-    parser_objective_MOI(moimodel, nvar, index_map)
+    parser_objective_MOI(moimodel, variables)
 
 Parse linear and quadratic objective of a `MOI.ModelLike`.
 """
-function parser_objective_MOI(moimodel, nvar, index_map)
+function parser_objective_MOI(moimodel, variables, index_map)
   _index(v::MOI.VariableIndex) = index_map[v].value
 
+  # Number of variables
+  nvar = length(variables)
+
   # Variables associated to linear and quadratic objective
   type = "UNKNOWN"
   constant = 0.0
@@ -751,11 +794,14 @@ function parser_objective_MOI(moimodel, nvar, index_map)
 end
 
 """
-    parser_linear_expression(cmodel, nvar, index_map, F)
+    parser_linear_expression(cmodel, variables, F)
 
 Parse linear expressions of type `VariableRef` and `GenericAffExpr{Float64,VariableRef}`.
 """
-function parser_linear_expression(cmodel, nvar, index_map, F)
+function parser_linear_expression(cmodel, variables, F)
+
+  # Number of variables
+  nvar = length(variables)
 
   # Variables associated to linear expressions
   rows = Int[]
@@ -810,7 +856,11 @@ function parser_linear_expression(cmodel, nvar, index_map, F)
     end
   end
   moimodel = backend(cmodel)
-  lls = parser_objective_MOI(moimodel, nvar, index_map)
+  index_map = MOI.Utilities.IndexMap()
+  for (idx, vi) in enumerate(variables)
+    index_map[vi] = MOI.VariableIndex(idx)
+  end
+  lls = parser_objective_MOI(moimodel, variables, index_map)
   return lls, LinearEquations(COO(rows, cols, vals), constants, length(vals)), nlinequ
 end
 
@@ -844,11 +894,14 @@ function add_constraint_model(Fmodel, Fi::AbstractArray)
 end
 
 """
-    parser_nonlinear_expression(cmodel, nvar, F; hessian)
+    parser_nonlinear_expression(cmodel, variables, F; hessian)
 
 Parse nonlinear expressions of type `NonlinearExpression`.
 """
-function parser_nonlinear_expression(cmodel, nvar, F; hessian::Bool = true)
+function parser_nonlinear_expression(cmodel, variables, F; hessian::Bool = true)
+
+  # Number of variables
+  nvar = length(variables)
 
   # Nonlinear least squares model
   F_is_array_of_containers = F isa Array{<:AbstractArray}