diff --git a/src/Nonlinear/Nonlinear.jl b/src/Nonlinear/Nonlinear.jl
index 9e8f0fc334..20e6a0807d 100644
--- a/src/Nonlinear/Nonlinear.jl
+++ b/src/Nonlinear/Nonlinear.jl
@@ -48,5 +48,6 @@ include("model.jl")
 include("evaluator.jl")
 
 include("ReverseAD/ReverseAD.jl")
+include("SymbolicAD/SymbolicAD.jl")
 
 end  # module
diff --git a/src/Nonlinear/SymbolicAD/SymbolicAD.jl b/src/Nonlinear/SymbolicAD/SymbolicAD.jl
new file mode 100644
index 0000000000..c18aa5b40d
--- /dev/null
+++ b/src/Nonlinear/SymbolicAD/SymbolicAD.jl
@@ -0,0 +1,225 @@
+# 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 SymbolicAD
+
+import MathOptInterface as MOI
+
+"""
+    simplify(f)
+
+Return a simplified version of the function `f`.
+
+!!! warning
+    This function is not type stable by design.
+"""
+simplify(f) = f
+
+function simplify(f::MOI.ScalarAffineFunction{T}) where {T}
+    f = MOI.Utilities.canonical(f)
+    if isempty(f.terms)
+        return f.constant
+    end
+    return f
+end
+
+function simplify(f::MOI.ScalarQuadraticFunction{T}) where {T}
+    f = MOI.Utilities.canonical(f)
+    if isempty(f.quadratic_terms)
+        return simplify(MOI.ScalarAffineFunction(f.affine_terms, f.constant))
+    end
+    return f
+end
+
+function simplify(f::MOI.ScalarNonlinearFunction)
+    for i in 1:length(f.args)
+        f.args[i] = simplify(f.args[i])
+    end
+    return _eval_if_constant(simplify(Val(f.head), f))
+end
+
+function simplify(f::MOI.VectorAffineFunction{T}) where {T}
+    f = MOI.Utilities.canonical(f)
+    if isempty(f.terms)
+        return f.constants
+    end
+    return f
+end
+
+function simplify(f::MOI.VectorQuadraticFunction{T}) where {T}
+    f = MOI.Utilities.canonical(f)
+    if isempty(f.quadratic_terms)
+        return simplify(MOI.VectorAffineFunction(f.affine_terms, f.constants))
+    end
+    return f
+end
+
+function simplify(f::MOI.VectorNonlinearFunction)
+    return MOI.VectorNonlinearFunction(simplify.(f.rows))
+end
+
+# If a ScalarNonlinearFunction has only constant arguments, we should return
+# the vaålue.
+
+_isnum(::Any) = false
+
+_isnum(::Union{Bool,Integer,Float64}) = true
+
+function _eval_if_constant(f::MOI.ScalarNonlinearFunction)
+    if all(_isnum, f.args) && hasproperty(Base, f.head)
+        return getproperty(Base, f.head)(f.args...)
+    end
+    return f
+end
+
+_eval_if_constant(f) = f
+
+_iszero(x::Any) = _isnum(x) && iszero(x)
+
+_isone(x::Any) = _isnum(x) && isone(x)
+
+"""
+    _isexpr(f::Any, head::Symbol[, n::Int])
+
+Return `true` if `f` is a `ScalarNonlinearFunction` with head `head` and, if
+specified, `n` arguments.
+"""
+_isexpr(::Any, ::Symbol, n::Int = 0) = false
+
+_isexpr(f::MOI.ScalarNonlinearFunction, head::Symbol) = f.head == head
+
+function _isexpr(f::MOI.ScalarNonlinearFunction, head::Symbol, n::Int)
+    return _isexpr(f, head) && length(f.args) == n
+end
+
+"""
+    simplify(::Val{head}, f::MOI.ScalarNonlinearFunction)
+
+Return a simplified version of `f` where the head of `f` is `head`.
+
+Implementing this method enables custom simplification rules for different
+operators without needing a giant switch statement.
+"""
+simplify(::Val, f::MOI.ScalarNonlinearFunction) = f
+
+function simplify(::Val{:*}, f::MOI.ScalarNonlinearFunction)
+    new_args = Any[]
+    first_constant = 0
+    for arg in f.args
+        if _isexpr(arg, :*)
+            # If the child is a :*, lift its arguments to the parent
+            append!(new_args, arg.args)
+        elseif _iszero(arg)
+            # If any argument is zero, the entire expression must be false
+            return false
+        elseif _isone(arg)
+            # Skip any arguments that are one
+        elseif arg isa Real
+            # Collect all constant arguments into a single value
+            if first_constant == 0
+                push!(new_args, arg)
+                first_constant = length(new_args)
+            else
+                new_args[first_constant] *= arg
+            end
+        else
+            push!(new_args, arg)
+        end
+    end
+    if isempty(new_args)
+        return true
+    elseif length(new_args) == 1
+        return only(new_args)
+    end
+    return MOI.ScalarNonlinearFunction(:*, new_args)
+end
+
+function simplify(::Val{:+}, f::MOI.ScalarNonlinearFunction)
+    if length(f.args) == 1
+        # +(x) -> x
+        return only(f.args)
+    elseif length(f.args) == 2 && _isexpr(f.args[2], :-, 1)
+        # +(x, -y) -> -(x, y)
+        return MOI.ScalarNonlinearFunction(
+            :-,
+            Any[f.args[1], f.args[2].args[1]],
+        )
+    end
+    new_args = Any[]
+    first_constant = 0
+    for arg in f.args
+        if _isexpr(arg, :+)
+            # If a child is a :+, lift its arguments to the parent
+            append!(new_args, arg.args)
+        elseif _iszero(arg)
+            # Skip any zero arguments
+        elseif arg isa Real
+            # Collect all constant arguments into a single value
+            if first_constant == 0
+                push!(new_args, arg)
+                first_constant = length(new_args)
+            else
+                new_args[first_constant] += arg
+            end
+        else
+            push!(new_args, arg)
+        end
+    end
+    if isempty(new_args)
+        # +() -> false
+        return false
+    elseif length(new_args) == 1
+        # +(x) -> x
+        return only(new_args)
+    end
+    return MOI.ScalarNonlinearFunction(:+, new_args)
+end
+
+function simplify(::Val{:-}, f::MOI.ScalarNonlinearFunction)
+    if length(f.args) == 1
+        if _isexpr(f.args[1], :-, 1)
+            # -(-(x)) => x
+            return f.args[1].args[1]
+        end
+    elseif length(f.args) == 2
+        if _iszero(f.args[1])
+            # 0 - x => -x
+            return MOI.ScalarNonlinearFunction(:-, Any[f.args[2]])
+        elseif _iszero(f.args[2])
+            # x - 0 => x
+            return f.args[1]
+        elseif f.args[1] == f.args[2]
+            # x - x => 0
+            return false
+        elseif _isexpr(f.args[2], :-, 1)
+            # x - -(y) => x + y
+            return MOI.ScalarNonlinearFunction(
+                :+,
+                Any[f.args[1], f.args[2].args[1]],
+            )
+        end
+    end
+    return f
+end
+
+function simplify(::Val{:^}, f::MOI.ScalarNonlinearFunction)
+    if _iszero(f.args[2])
+        # x^0 => 1
+        return true
+    elseif _isone(f.args[2])
+        # x^1 => x
+        return f.args[1]
+    elseif _iszero(f.args[1])
+        # 0^x => 0
+        return false
+    elseif _isone(f.args[1])
+        # 1^x => 1
+        return true
+    end
+    return f
+end
+
+end  # module
diff --git a/src/Utilities/functions.jl b/src/Utilities/functions.jl
index 2d558f92e4..9f238913da 100644
--- a/src/Utilities/functions.jl
+++ b/src/Utilities/functions.jl
@@ -1071,6 +1071,10 @@ function canonicalize!(
     return f
 end
 
+function canonical(f::MOI.ScalarNonlinearFunction)::MOI.ScalarNonlinearFunction
+    return MOI.Nonlinear.SymbolicAD.simplify(f)
+end
+
 function canonicalize!(f::MOI.ScalarNonlinearFunction)
     for (i, arg) in enumerate(f.args)
         if !is_canonical(arg)
@@ -1080,6 +1084,11 @@ function canonicalize!(f::MOI.ScalarNonlinearFunction)
     return f
 end
 
+function canonical(f::MOI.VectorNonlinearFunction)
+    rows = MOI.Nonlinear.SymbolicAD.simplify.(f.rows)
+    return MOI.VectorNonlinearFunction(rows)
+end
+
 function canonicalize!(f::MOI.VectorNonlinearFunction)
     for (i, fi) in enumerate(f.rows)
         f.rows[i] = canonicalize!(fi)
diff --git a/src/functions.jl b/src/functions.jl
index bc658737c4..3bf7f08723 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -926,11 +926,23 @@ function _is_approx(x::AbstractArray, y::AbstractArray; kwargs...)
            all(z -> _is_approx(z[1], z[2]; kwargs...), zip(x, y))
 end
 
+function _is_univariate_plus(f)
+    if f.head == :+ && length(f.args) == 1
+        return only(f.args) isa ScalarNonlinearFunction
+    end
+    return false
+end
+
 function Base.isapprox(
     f::ScalarNonlinearFunction,
     g::ScalarNonlinearFunction;
     kwargs...,
 )
+    if _is_univariate_plus(f)
+        return isapprox(only(f.args), g.args; kwargs...)
+    elseif _is_univariate_plus(g)
+        return isapprox(f, only(g.args); kwargs...)
+    end
     if f.head != g.head || length(f.args) != length(g.args)
         return false
     end
@@ -1127,22 +1139,40 @@ _order(x::VariableIndex, y::Real, z::VariableIndex) = (y, x, z)
 _order(x::VariableIndex, y::VariableIndex, z::Real) = (z, x, y)
 _order(x, y, z) = nothing
 
+_order_quad(x, y) = nothing
+_order_quad(x::VariableIndex, y::VariableIndex) = (x, y)
+
 function Base.convert(
     ::Type{ScalarQuadraticTerm{T}},
     f::ScalarNonlinearFunction,
 ) where {T}
-    if f.head != :* || length(f.args) != 3
+    if f.head != :*
         throw(InexactError(:convert, ScalarQuadraticTerm, f))
+    elseif length(f.args) == 2
+        # Deal with *(x, y)
+        ret_2 = _order_quad(f.args[1], f.args[2])
+        if ret_2 === nothing
+            throw(InexactError(:convert, ScalarQuadraticTerm, f))
+        end
+        coef = one(T)
+        if ret_2[1] == ret_2[2]
+            coef *= 2
+        end
+        return ScalarQuadraticTerm(coef, ret_2[1], ret_2[2])
+    elseif length(f.args) == 3
+        # *(constant, x, y)
+        ret = _order(f.args[1], f.args[2], f.args[3])
+        if ret === nothing
+            throw(InexactError(:convert, ScalarQuadraticTerm, f))
+        end
+        coef = convert(T, ret[1])
+        if ret[2] == ret[3]
+            coef *= 2
+        end
+        return ScalarQuadraticTerm(coef, ret[2], ret[3])
+    else
+        return throw(InexactError(:convert, ScalarQuadraticTerm, f))
     end
-    ret = _order(f.args[1], f.args[2], f.args[3])
-    if ret === nothing
-        throw(InexactError(:convert, ScalarQuadraticTerm, f))
-    end
-    coef = convert(T, ret[1])
-    if ret[2] == ret[3]
-        coef *= 2
-    end
-    return ScalarQuadraticTerm(coef, ret[2], ret[3])
 end
 
 function _add_to_function(
@@ -1157,7 +1187,11 @@ function _add_to_function(
     arg::ScalarNonlinearFunction,
 ) where {T}
     if arg.head == :* && length(arg.args) == 2
-        push!(f.affine_terms, convert(ScalarAffineTerm{T}, arg))
+        if _order_quad(arg.args[1], arg.args[2]) === nothing
+            push!(f.affine_terms, convert(ScalarAffineTerm{T}, arg))
+        else
+            push!(f.quadratic_terms, convert(ScalarQuadraticTerm{T}, arg))
+        end
     elseif arg.head == :* && length(arg.args) == 3
         push!(f.quadratic_terms, convert(ScalarQuadraticTerm{T}, arg))
     else
@@ -1174,15 +1208,12 @@ function Base.convert(
     f::ScalarNonlinearFunction,
 ) where {T}
     if f.head == :*
-        if length(f.args) == 2
-            quad_terms = ScalarQuadraticTerm{T}[]
-            affine_terms = [convert(ScalarAffineTerm{T}, f)]
-            return ScalarQuadraticFunction{T}(quad_terms, affine_terms, zero(T))
-        elseif length(f.args) == 3
-            quad_terms = [convert(ScalarQuadraticTerm{T}, f)]
-            affine_terms = ScalarAffineTerm{T}[]
-            return ScalarQuadraticFunction{T}(quad_terms, affine_terms, zero(T))
-        end
+        g = ScalarQuadraticFunction{T}(
+            ScalarQuadraticTerm{T}[],
+            ScalarAffineTerm{T}[],
+            zero(T),
+        )
+        return _add_to_function(g, f)
     elseif f.head == :^ && length(f.args) == 2 && f.args[2] == 2
         return convert(
             ScalarQuadraticFunction{T},
diff --git a/test/Bridges/Constraint/NormInfinityBridge.jl b/test/Bridges/Constraint/NormInfinityBridge.jl
index f61f482572..4f37077e81 100644
--- a/test/Bridges/Constraint/NormInfinityBridge.jl
+++ b/test/Bridges/Constraint/NormInfinityBridge.jl
@@ -647,9 +647,9 @@ function test_NormInfinity_VectorNonlinearFunction()
     g = MOI.VectorNonlinearFunction([
         MOI.ScalarNonlinearFunction(
             :+,
-            Any[MOI.ScalarNonlinearFunction(:-, Any[v_sin]), u_p],
+            Any[MOI.ScalarNonlinearFunction(:-, Any[v_sin]), u],
         ),
-        MOI.ScalarNonlinearFunction(:+, Any[v_sin, u_p]),
+        MOI.ScalarNonlinearFunction(:+, Any[v_sin, u]),
     ])
     @test ≈(MOI.get(inner, MOI.ConstraintFunction(), indices[1]), g)
     h = MOI.VectorNonlinearFunction([
@@ -689,10 +689,7 @@ function test_NormOne_VectorNonlinearFunction()
     u, v, w = inner_variables
     v_sin = MOI.ScalarNonlinearFunction(:sin, Any[v])
     g = MOI.VectorNonlinearFunction([
-        MOI.ScalarNonlinearFunction(
-            :-,
-            Any[MOI.ScalarNonlinearFunction(:+, Any[u]), 0.0+1.0*w],
-        ),
+        MOI.ScalarNonlinearFunction(:-, Any[u, 0.0+1.0*w]),
         MOI.ScalarNonlinearFunction(
             :+,
             Any[MOI.ScalarNonlinearFunction(:-, Any[v_sin]), w],
diff --git a/test/Bridges/Constraint/SquareBridge.jl b/test/Bridges/Constraint/SquareBridge.jl
index e9e1c0036a..12526956f4 100644
--- a/test/Bridges/Constraint/SquareBridge.jl
+++ b/test/Bridges/Constraint/SquareBridge.jl
@@ -295,21 +295,19 @@ function test_VectorNonlinearFunction_mixed_type()
     @test length(indices) == 1
     g = MOI.get(inner, MOI.ConstraintFunction(), indices[1])
     y = MOI.get(inner, MOI.ListOfVariableIndices())
-    gis = vcat(
-        Any[MOI.ScalarNonlinearFunction(:log, Any[y[i]]) for i in 1:2],
-        1.0 * y[3] + 2.0,
-        y[4],
-    )
+    gis = [
+        MOI.ScalarNonlinearFunction(:log, Any[y[1]]),
+        MOI.ScalarNonlinearFunction(:log, Any[y[2]]),
+        MOI.ScalarNonlinearFunction(:+, Any[y[3], 2.0]),
+        MOI.ScalarNonlinearFunction(:+, Any[y[4]]),
+    ]
     @test g ≈ MOI.VectorNonlinearFunction(gis[[1, 3, 4]])
     F, S = MOI.ScalarNonlinearFunction, MOI.EqualTo{Float64}
     indices = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}())
     @test length(indices) == 1
     @test ≈(
         MOI.get(inner, MOI.ConstraintFunction(), indices[1]),
-        MOI.ScalarNonlinearFunction(
-            :-,
-            Any[convert(MOI.ScalarNonlinearFunction, gis[3]), gis[2]],
-        ),
+        MOI.ScalarNonlinearFunction(:-, Any[gis[3], gis[2]]),
     )
     return
 end
diff --git a/test/Nonlinear/Nonlinear.jl b/test/Nonlinear/Nonlinear.jl
index 743ed350ed..9045d93589 100644
--- a/test/Nonlinear/Nonlinear.jl
+++ b/test/Nonlinear/Nonlinear.jl
@@ -1219,3 +1219,4 @@ end  # TestNonlinear
 TestNonlinear.runtests()
 
 include("ReverseAD.jl")
+include("SymbolicAD.jl")
diff --git a/test/Nonlinear/SymbolicAD.jl b/test/Nonlinear/SymbolicAD.jl
new file mode 100644
index 0000000000..d8f768d2cc
--- /dev/null
+++ b/test/Nonlinear/SymbolicAD.jl
@@ -0,0 +1,262 @@
+# 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 TestSymbolicAD
+
+using Test
+import MathOptInterface as MOI
+import MathOptInterface.Nonlinear: SymbolicAD
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+function test_simplify()
+    x = MOI.VariableIndex(1)
+    @test SymbolicAD.simplify(x) === x
+    @test SymbolicAD.simplify(1.0) === 1.0
+    return
+end
+
+function test_simplify_ScalarAffineFunction()
+    f = zero(MOI.ScalarAffineFunction{Float64})
+    @test SymbolicAD.simplify(f) == 0.0
+    f = MOI.ScalarAffineFunction{Float64}(MOI.ScalarAffineTerm{Float64}[], 2.0)
+    @test SymbolicAD.simplify(f) == 2.0
+    x = MOI.VariableIndex(1)
+    @test SymbolicAD.simplify(1.0 * x + 1.0) ≈ 1.0 * x + 1.0
+    @test SymbolicAD.simplify(1.0 * x + 2.0 * x + 1.0) ≈ 3.0 * x + 1.0
+    return
+end
+
+function test_simplify_ScalarQuadraticFunction()
+    x = MOI.VariableIndex(1)
+    f = MOI.ScalarQuadraticFunction(
+        MOI.ScalarQuadraticTerm{Float64}[],
+        [MOI.ScalarAffineTerm{Float64}(1.0, x)],
+        1.0,
+    )
+    @test SymbolicAD.simplify(f) ≈ 1.0 * x + 1.0
+    @test SymbolicAD.simplify(1.0 * x * x + 1.0) ≈ 1.0 * x * x + 1.0
+    g = 1.0 * x * x + 2.0 * x * x + 1.0
+    @test SymbolicAD.simplify(g) ≈ 3.0 * x * x + 1.0
+    return
+end
+
+function test_simplify_ScalarNonlinearFunction()
+    x = MOI.VariableIndex(1)
+    # sin(3 * (x^0)) -> sin(3)
+    f = MOI.ScalarNonlinearFunction(:^, Any[x, 0])
+    g = MOI.ScalarNonlinearFunction(:*, Any[3, f])
+    h = MOI.ScalarNonlinearFunction(:sin, Any[g])
+    @test SymbolicAD.simplify(h) ≈ sin(3)
+    # sin(log(x)) -> sin(log(x))
+    f = MOI.ScalarNonlinearFunction(:log, Any[x])
+    g = MOI.ScalarNonlinearFunction(:sin, Any[f])
+    @test SymbolicAD.simplify(g) ≈ g
+    @test MOI.Utilities.canonical(g) ≈ g
+    return
+end
+
+# simplify(::Val{:*}, f::MOI.ScalarNonlinearFunction)
+function test_simplify_ScalarNonlinearFunction_multiplication()
+    x, y, z = MOI.VariableIndex.(1:3)
+    # *(x, *(y, z)) -> *(x, y, z)
+    @test ≈(
+        SymbolicAD.simplify(
+            MOI.ScalarNonlinearFunction(
+                :*,
+                Any[x, MOI.ScalarNonlinearFunction(:*, Any[y, z])],
+            ),
+        ),
+        MOI.ScalarNonlinearFunction(:*, Any[x, y, z]),
+    )
+    # *(x, *(y, z, *(x, 2))) -> *(x, y, z, x, 2)
+    f = MOI.ScalarNonlinearFunction(:*, Any[x, 2])
+    @test ≈(
+        SymbolicAD.simplify(
+            MOI.ScalarNonlinearFunction(
+                :*,
+                Any[x, MOI.ScalarNonlinearFunction(:*, Any[y, z, f])],
+            ),
+        ),
+        MOI.ScalarNonlinearFunction(:*, Any[x, y, z, x, 2]),
+    )
+    # *(x, 3, 2) -> *(x, 6)
+    @test ≈(
+        SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:*, Any[x, 3, 2])),
+        MOI.ScalarNonlinearFunction(:*, Any[x, 6]),
+    )
+    # *(3, x, 2) -> *(6, x)
+    @test ≈(
+        SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:*, Any[3, x, 2])),
+        MOI.ScalarNonlinearFunction(:*, Any[6, x]),
+    )
+    # *(x, 1) -> x
+    @test ≈(SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:*, Any[x, 1])), x)
+    # *(x, 0) -> 0
+    @test ≈(SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:*, Any[x, 0])), 0)
+    # *(-(x, x), 1) -> 0
+    f = MOI.ScalarNonlinearFunction(:-, Any[x, x])
+    @test ≈(SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:*, Any[f, 1])), 0)
+    # *() -> true
+    @test ≈(SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:*, Any[])), 1)
+    return
+end
+
+# simplify(::Val{:+}, f::MOI.ScalarNonlinearFunction)
+function test_simplify_ScalarNonlinearFunction_addition()
+    x, y, z = MOI.VariableIndex.(1:3)
+    # (+(x, +(y, z)))=>(+(x, y, z)),
+    @test ≈(
+        SymbolicAD.simplify(
+            MOI.ScalarNonlinearFunction(
+                :+,
+                Any[x, MOI.ScalarNonlinearFunction(:+, Any[y, z])],
+            ),
+        ),
+        MOI.ScalarNonlinearFunction(:+, Any[x, y, z]),
+    )
+    # +(sin(x), -cos(x))=>sin(x)-cos(x),
+    sinx = MOI.ScalarNonlinearFunction(:sin, Any[x])
+    cosx = MOI.ScalarNonlinearFunction(:cos, Any[x])
+    @test ≈(
+        SymbolicAD.simplify(
+            MOI.ScalarNonlinearFunction(
+                :+,
+                Any[sinx, MOI.ScalarNonlinearFunction(:-, Any[cosx])],
+            ),
+        ),
+        MOI.ScalarNonlinearFunction(:-, Any[sinx, cosx]),
+    )
+    # (+(x, 1, 2))=>(+(x, 3)),
+    @test ≈(
+        SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:+, Any[x, 1, 2])),
+        MOI.ScalarNonlinearFunction(:+, Any[x, 3]),
+    )
+    # (+(1, x, 2))=>(+(3, x)),
+    @test ≈(
+        SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:+, Any[1, x, 2])),
+        MOI.ScalarNonlinearFunction(:+, Any[3, x]),
+    )
+    # +(x, 0) -> x
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:+, Any[x, 0])) ≈ x
+    # +(0, x) -> x
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:+, Any[0, x])) ≈ x
+    # +(-(x, x), 0) -> 0
+    f = MOI.ScalarNonlinearFunction(
+        :+,
+        Any[MOI.ScalarNonlinearFunction(:-, Any[x, x]), 0],
+    )
+    @test SymbolicAD.simplify(f) === false
+    return
+end
+
+# simplify(::Val{:-}, f::MOI.ScalarNonlinearFunction)
+function test_simplify_ScalarNonlinearFunction_subtraction()
+    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
+    f = MOI.ScalarNonlinearFunction(:-, Any[x])
+    # -x -> -x
+    @test SymbolicAD.simplify(f) ≈ f
+    # -(-(x)) -> x
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:-, Any[f])) ≈ x
+    # -(x, 0) -> x
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:-, Any[x, 0])) ≈ x
+    # -(0, x) -> -x
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:-, Any[0, x])) ≈ f
+    # -(x, x) -> 0
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:-, Any[x, x])) ≈ 0
+    # -(x, -y) -> +(x, y)
+    f = MOI.ScalarNonlinearFunction(
+        :-,
+        Any[x, MOI.ScalarNonlinearFunction(:-, Any[y])],
+    )
+    @test SymbolicAD.simplify(f) ≈ MOI.ScalarNonlinearFunction(:+, Any[x, y])
+    # -(x, y) -> -(x, y)
+    f = MOI.ScalarNonlinearFunction(:-, Any[x, y])
+    @test SymbolicAD.simplify(f) ≈ f
+    return
+end
+
+# simplify(::Val{:^}, f::MOI.ScalarNonlinearFunction)
+function test_simplify_ScalarNonlinearFunction_power()
+    x, y = MOI.VariableIndex(1), MOI.VariableIndex(2)
+    # x^0 -> 1
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:^, Any[x, 0])) == 1
+    # x^1 -> x
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:^, Any[x, 1])) == x
+    # 0^x -> 0
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:^, Any[0, x])) == 0
+    # 1^x -> 1
+    @test SymbolicAD.simplify(MOI.ScalarNonlinearFunction(:^, Any[1, x])) == 1
+    # x^y -> x^y
+    f = MOI.ScalarNonlinearFunction(:^, Any[x, y])
+    @test SymbolicAD.simplify(f) ≈ f
+    return
+end
+
+function test_simplify_VectorAffineFunction()
+    f = MOI.VectorAffineFunction{Float64}(
+        MOI.VectorAffineTerm{Float64}[],
+        [0.0, 1.0, 2.0],
+    )
+    @test SymbolicAD.simplify(f) == [0.0, 1.0, 2.0]
+    x = MOI.VariableIndex(1)
+    f = MOI.Utilities.operate(vcat, Float64, 1.0, x, 2.0 * x + 1.0 * x)
+    @test SymbolicAD.simplify(f) ≈ f
+    return
+end
+
+function test_simplify_VectorQuadraticFunction()
+    f = MOI.VectorQuadraticFunction{Float64}(
+        MOI.VectorQuadraticTerm{Float64}[],
+        MOI.VectorAffineTerm{Float64}[],
+        [0.0, 1.0, 2.0],
+    )
+    @test SymbolicAD.simplify(f) == [0.0, 1.0, 2.0]
+    x = MOI.VariableIndex(1)
+    f = MOI.VectorQuadraticFunction{Float64}(
+        MOI.VectorQuadraticTerm{Float64}[],
+        [MOI.VectorAffineTerm{Float64}(2, MOI.ScalarAffineTerm(3.0, x))],
+        [1.0, 0.0],
+    )
+    g = MOI.Utilities.operate(vcat, Float64, 1.0, 3.0 * x)
+    @test SymbolicAD.simplify(f) ≈ g
+    f = MOI.Utilities.operate(vcat, Float64, 1.0, 2.0 * x * x)
+    @test SymbolicAD.simplify(f) ≈ f
+    return
+end
+
+function test_simplify_VectorNonlinearFunction()
+    x = MOI.VariableIndex.(1:3)
+    y = MOI.ScalarNonlinearFunction(
+        :+,
+        Any[MOI.ScalarNonlinearFunction(:^, Any[xi, 2]) for xi in x],
+    )
+    x_plus = [MOI.ScalarNonlinearFunction(:+, Any[xi]) for xi in x]
+    function wrap(f)
+        return MOI.ScalarNonlinearFunction(
+            :+,
+            Any[MOI.ScalarNonlinearFunction(:-, Any[f, 0.0]), 0.0],
+        )
+    end
+    f = MOI.VectorNonlinearFunction(wrap.([y; x_plus]))
+    g = MOI.VectorNonlinearFunction([y; x_plus])
+    @test SymbolicAD.simplify(f) ≈ g
+    @test isapprox(MOI.Utilities.canonical(f), g)
+    return
+end
+
+end  # module
+
+TestSymbolicAD.runtests()
diff --git a/test/functions.jl b/test/functions.jl
index 503dc50da5..44d5cd586a 100644
--- a/test/functions.jl
+++ b/test/functions.jl
@@ -396,10 +396,13 @@ function test_convert_ScalarNonlinearFunction_ScalarQuadraticTerm()
     @test convert(MOI.ScalarQuadraticTerm{Float64}, g) == f
     @test convert(MOI.ScalarQuadraticTerm{Float64}, h) == f
     @test convert(MOI.ScalarQuadraticTerm{Float64}, i) == f
+    j = MOI.ScalarNonlinearFunction(:*, Any[x, x])
+    @test convert(MOI.ScalarQuadraticTerm{Float64}, j) == f
     for f_error in (
         MOI.ScalarNonlinearFunction(:*, Any[1.0, x, 2.0]),
         MOI.ScalarNonlinearFunction(:+, Any[1.0, x, x]),
-        MOI.ScalarNonlinearFunction(:*, Any[x, x]),
+        MOI.ScalarNonlinearFunction(:*, Any[1.0, x]),
+        MOI.ScalarNonlinearFunction(:*, Any[x]),
     )
         @test_throws(
             InexactError,