diff --git a/src/Bridges/Constraint/bridges/quad_to_soc.jl b/src/Bridges/Constraint/bridges/quad_to_soc.jl
index 465ae4f5c9..b7a081288d 100644
--- a/src/Bridges/Constraint/bridges/quad_to_soc.jl
+++ b/src/Bridges/Constraint/bridges/quad_to_soc.jl
@@ -292,7 +292,8 @@ function MOI.set(
     # |        U * x        |
     # where we compute `x'Qx/2` and `U * x` using the starting values of the variable.
     soc = MOI.get(model, MOI.ConstraintFunction(), bridge.soc)
-    Ux = MOI.Utilities.eval_variables(MOI.Utilities.eachscalar(soc)[3:end]) do v
+    f = MOI.Utilities.eachscalar(soc)[3:end]
+    Ux = MOI.Utilities.eval_variables(model, f) do v
         return _primal_start_or_error(model, attr, v)
     end
     if bridge.less_than
diff --git a/src/Bridges/Objective/bridges/slack.jl b/src/Bridges/Objective/bridges/slack.jl
index 739624b4ff..5b45599967 100644
--- a/src/Bridges/Objective/bridges/slack.jl
+++ b/src/Bridges/Objective/bridges/slack.jl
@@ -224,7 +224,7 @@ function MOI.set(
     end
     MOI.set(model, MOI.VariablePrimalStart(), b.slack, zero(T))
     f = MOI.get(model, MOI.ConstraintFunction(), b.constraint)
-    f_val = MOI.Utilities.eval_variables(f) do v
+    f_val = MOI.Utilities.eval_variables(model, f) do v
         return MOI.get(model, MOI.VariablePrimalStart(), v)
     end
     f_val -= MOI.constant(MOI.get(model, MOI.ConstraintSet(), b.constraint))
diff --git a/src/Nonlinear/operators.jl b/src/Nonlinear/operators.jl
index 7e51f108fc..ee49768111 100644
--- a/src/Nonlinear/operators.jl
+++ b/src/Nonlinear/operators.jl
@@ -895,3 +895,72 @@ function eval_comparison_function(
         return lhs > rhs
     end
 end
+
+# This method is implmented here because it needs the OperatorRegistry type from
+# the Nonlinear, which doesn't exist when the Utilities submodule is defined.
+
+function MOI.Utilities.eval_variables(
+    value_fn::F,
+    model::MOI.ModelLike,
+    f::MOI.ScalarNonlinearFunction,
+) where {F}
+    registry = OperatorRegistry()
+    return _evaluate_expr(registry, value_fn, model, f)
+end
+
+function _evaluate_expr(
+    ::OperatorRegistry,
+    value_fn::Function,
+    model::MOI.ModelLike,
+    f::MOI.AbstractFunction,
+)
+    return MOI.Utilities.eval_variables(value_fn, model, f)
+end
+
+function _evaluate_expr(
+    ::OperatorRegistry,
+    ::Function,
+    ::MOI.ModelLike,
+    f::Number,
+)
+    return f
+end
+
+function _evaluate_expr(
+    registry::OperatorRegistry,
+    value_fn::Function,
+    model::MOI.ModelLike,
+    expr::MOI.ScalarNonlinearFunction,
+)
+    op = expr.head
+    if !_is_registered(registry, op, length(expr.args))
+        udf = MOI.get(model, MOI.UserDefinedFunction(op, length(expr.args)))
+        if udf === nothing
+            throw(MOI.UnsupportedNonlinearOperator(op))
+        end
+        args = map(expr.args) do arg
+            return _evaluate_expr(registry, value_fn, model, arg)
+        end
+        return first(udf)(args...)
+    end
+    if length(expr.args) == 1 && haskey(registry.univariate_operator_to_id, op)
+        arg = _evaluate_expr(registry, value_fn, model, expr.args[1])
+        return eval_univariate_function(registry, op, arg)
+    elseif haskey(registry.multivariate_operator_to_id, op)
+        args = map(expr.args) do arg
+            return _evaluate_expr(registry, value_fn, model, arg)
+        end
+        return eval_multivariate_function(registry, op, args)
+    elseif haskey(registry.logic_operator_to_id, op)
+        @assert length(expr.args) == 2
+        x = _evaluate_expr(registry, value_fn, model, expr.args[1])
+        y = _evaluate_expr(registry, value_fn, model, expr.args[2])
+        return eval_logic_function(registry, op, x, y)
+    else
+        @assert haskey(registry.comparison_operator_to_id, op)
+        @assert length(expr.args) == 2
+        x = _evaluate_expr(registry, value_fn, model, expr.args[1])
+        y = _evaluate_expr(registry, value_fn, model, expr.args[2])
+        return eval_comparison_function(registry, op, x, y)
+    end
+end
diff --git a/src/Utilities/functions.jl b/src/Utilities/functions.jl
index fe87c8bce1..80d618b741 100644
--- a/src/Utilities/functions.jl
+++ b/src/Utilities/functions.jl
@@ -78,7 +78,7 @@ Returns the output type that results if a function of type `F` is evaluated
 using variables with numeric type `T`.
 
 In other words, this is the return type for
-`MOI.Utilities.eval_variables(value_fn::Function, f::F)` for a function
+`MOI.Utilities.eval_variables(value_fn::Function, model, f::F)` for a function
 `value_fn(::MOI.VariableIndex)::T`.
 """
 function value_type end
@@ -107,7 +107,13 @@ Returns the value of function `f` if each variable index `vi` is evaluated as
 
 Note that `value_fn` must return a Number. See [`substitute_variables`](@ref)
 for a similar function where `value_fn` returns an
-    [`MOI.AbstractScalarFunction`](@ref).
+[`MOI.AbstractScalarFunction`](@ref).
+
+!!! warning
+    The two-argument version of `eval_variables` is deprecated and may be
+    removed in MOI v2.0.0. Use the three-argument method
+    `eval_variables(::Function, ::MOI.ModelLike, ::MOI.AbstractFunction)`
+    instead.
 """
 function eval_variables end
 
@@ -164,6 +170,31 @@ function eval_variables(value_fn::Function, f::MOI.VectorQuadraticFunction)
     return out
 end
 
+"""
+    eval_variables(
+        value_fn::Function,
+        model::MOI.ModelLike,
+        f::MOI.AbstractFunction,
+    )
+
+Returns the value of function `f` if each variable index `vi` is evaluated as
+`value_fn(vi)`.
+
+Note that `value_fn` must return a Number. See [`substitute_variables`](@ref)
+for a similar function where `value_fn` returns an
+[`MOI.AbstractScalarFunction`](@ref).
+"""
+function eval_variables(
+    value_fn::F,
+    model::MOI.ModelLike,
+    f::MOI.AbstractFunction,
+) where {F}
+    return eval_variables(value_fn, f)
+end
+
+# The `eval_variables(::F, ::MOI.ModelLike, ::MOI.ScalarNonlinearFunction)`
+# method is defined in the MOI.Nonlinear submodule.
+
 """
     map_indices(index_map::Function, attr::MOI.AnyAttribute, x::X)::X where {X}
 
diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl
index 0ea130dc72..b56820d9b1 100644
--- a/src/Utilities/results.jl
+++ b/src/Utilities/results.jl
@@ -26,10 +26,9 @@ function get_fallback(model::MOI.ModelLike, attr::MOI.ObjectiveValue)
     MOI.check_result_index_bounds(model, attr)
     F = MOI.get(model, MOI.ObjectiveFunctionType())
     f = MOI.get(model, MOI.ObjectiveFunction{F}())
-    obj = eval_variables(
-        vi -> MOI.get(model, MOI.VariablePrimal(attr.result_index), vi),
-        f,
-    )
+    obj = eval_variables(model, f) do vi
+        return MOI.get(model, MOI.VariablePrimal(attr.result_index), vi)
+    end
     if is_ray(MOI.get(model, MOI.PrimalStatus()))
         # Dual infeasibiltiy certificates do not include the primal objective
         # constant.
@@ -187,7 +186,7 @@ function get_fallback(
 )
     MOI.check_result_index_bounds(model, attr)
     f = MOI.get(model, MOI.ConstraintFunction(), idx)
-    c = eval_variables(f) do vi
+    c = eval_variables(model, f) do vi
         return MOI.get(model, MOI.VariablePrimal(attr.result_index), vi)
     end
     if is_ray(MOI.get(model, MOI.PrimalStatus()))
diff --git a/test/Utilities/functions.jl b/test/Utilities/functions.jl
index 77cd8fae4a..79b214f5fc 100644
--- a/test/Utilities/functions.jl
+++ b/test/Utilities/functions.jl
@@ -247,6 +247,42 @@ function test_eval_variables()
     return
 end
 
+function test_eval_variables_scalar_nonlinear_function()
+    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
+    x = MOI.add_variable(model)
+    f = MOI.ScalarNonlinearFunction(:log, Any[x])
+    @test MOI.Utilities.eval_variables(xi -> 0.5, model, f) ≈ log(0.5)
+    f = MOI.ScalarNonlinearFunction(:*, Any[x, 2.0*x, 1.0*x+2.0])
+    @test MOI.Utilities.eval_variables(xi -> 0.5, model, f) ≈
+          0.5 * (2.0 * 0.5) * (1.0 * 0.5 + 2.0)
+    f = MOI.ScalarNonlinearFunction(
+        :ifelse,
+        Any[MOI.ScalarNonlinearFunction(:<, Any[x, 1.0]), 0.0, x],
+    )
+    @test MOI.Utilities.eval_variables(xi -> 0.5, model, f) ≈ 0.0
+    @test MOI.Utilities.eval_variables(xi -> 1.5, model, f) ≈ 1.5
+    f = MOI.ScalarNonlinearFunction(
+        :||,
+        Any[
+            MOI.ScalarNonlinearFunction(:<, Any[x, 0.0]),
+            MOI.ScalarNonlinearFunction(:>, Any[x, 1.0]),
+        ],
+    )
+    @test MOI.Utilities.eval_variables(xi -> 0.5, model, f) ≈ 0.0
+    @test MOI.Utilities.eval_variables(xi -> 1.5, model, f) ≈ 1.0
+    @test MOI.Utilities.eval_variables(xi -> -0.5, model, f) ≈ 1.0
+    my_square(x, y) = (x - y)^2
+    MOI.set(model, MOI.UserDefinedFunction(:my_square, 2), (my_square,))
+    f = MOI.ScalarNonlinearFunction(:my_square, Any[x, 1.0])
+    @test MOI.Utilities.eval_variables(xi -> 0.5, model, f) ≈ (0.5 - 1.0)^2
+    f = MOI.ScalarNonlinearFunction(:bad_f, Any[x, 1.0])
+    @test_throws(
+        MOI.UnsupportedNonlinearOperator(:bad_f),
+        MOI.Utilities.eval_variables(xi -> 0.5, model, f)
+    )
+    return
+end
+
 function test_substitute_variables()
     # We do tests twice to make sure the function is not modified
     subs = Dict(w => 1.0y + 1.0z, x => 2.0y + 1.0, y => 1.0y, z => -1.0w)