diff --git a/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl b/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl
index b6c1601..f21214f 100644
--- a/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl
+++ b/ext/MultiObjectiveAlgorithmsPolyhedraExt.jl
@@ -37,12 +37,14 @@ end
 function MOA.minimize_multiobjective!(
     algorithm::MOA.Sandwiching,
     model::MOA.Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
 )
-    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+    @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
     solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}()
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    n = MOI.output_dimension(model.f)
-    scalars = MOI.Utilities.scalarize(model.f)
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
+    n = MOI.output_dimension(f)
+    scalars = MOI.Utilities.scalarize(f)
     status = MOI.OPTIMAL
     δ_OPS_optimizer = MOI.instantiate(model.optimizer_factory)
     if MOI.supports(δ_OPS_optimizer, MOI.Silent())
@@ -52,26 +54,26 @@ function MOA.minimize_multiobjective!(
     anchors = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}()
     yI, yUB = zeros(n), zeros(n)
     for (i, f_i) in enumerate(scalars)
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
         MOA.optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !MOA._is_scalar_status_optimal(model)
             return status, nothing
         end
-        X, Y = MOA._compute_point(model, variables, model.f)
+        X, Y = MOA._compute_point(model, variables, f)
         model.ideal_point[i] = Y[i]
         yI[i] = Y[i]
         anchors[Y] = X
-        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+        MOI.set(inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
         MOA.optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !MOA._is_scalar_status_optimal(model)
             MOA._warn_on_nonfinite_anti_ideal(algorithm, MOI.MIN_SENSE, i)
             return status, nothing
         end
         _, Y = MOA._compute_point(model, variables, f_i)
         yUB[i] = Y
-        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+        MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
         e_i = Float64.(1:n .== i)
         MOI.add_constraint(
             δ_OPS_optimizer,
@@ -86,14 +88,13 @@ function MOA.minimize_multiobjective!(
     end
     IPS = [yUB, keys(anchors)...]
     merge!(solutions, anchors)
-    u = MOI.add_variables(model.inner, n)
+    u = MOI.add_variables(inner, n)
     u_constraints = [ # u_i >= 0 for all i = 1:n
-        MOI.add_constraint(model.inner, u_i, MOI.GreaterThan{Float64}(0))
-        for u_i in u
+        MOI.add_constraint(inner, u_i, MOI.GreaterThan{Float64}(0)) for u_i in u
     ]
     f_constraints = [ # f_i + u_i <= yUB_i for all i = 1:n
         MOI.Utilities.normalize_and_add_constraint(
-            model.inner,
+            inner,
             scalars[i] + u[i],
             MOI.LessThan(yUB[i]),
         ) for i in 1:n
@@ -113,14 +114,14 @@ function MOA.minimize_multiobjective!(
         end
         # would not terminate when precision is set to 0
         new_f = sum(w[i] * (scalars[i] + u[i]) for i in 1:n) # w' * (f(x) + u)
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
         MOA.optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !MOA._is_scalar_status_optimal(model)
             return status, nothing
         end
-        β̄ = MOI.get(model.inner, MOI.ObjectiveValue())
-        X, Y = MOA._compute_point(model, variables, model.f)
+        β̄ = MOI.get(inner, MOI.ObjectiveValue())
+        X, Y = MOA._compute_point(model, variables, f)
         solutions[Y] = X
         MOI.add_constraint(
             δ_OPS_optimizer,
@@ -130,9 +131,9 @@ function MOA.minimize_multiobjective!(
         IPS = push!(IPS, Y)
         H = _halfspaces(IPS)
     end
-    MOI.delete.(model.inner, f_constraints)
-    MOI.delete.(model.inner, u_constraints)
-    MOI.delete.(model.inner, u)
+    MOI.delete.(inner, f_constraints)
+    MOI.delete.(inner, u_constraints)
+    MOI.delete.(inner, u)
     return status, [MOA.SolutionPoint(X, Y) for (Y, X) in solutions]
 end
 
diff --git a/src/MultiObjectiveAlgorithms.jl b/src/MultiObjectiveAlgorithms.jl
index a4dccd3..f444e90 100644
--- a/src/MultiObjectiveAlgorithms.jl
+++ b/src/MultiObjectiveAlgorithms.jl
@@ -701,7 +701,9 @@ This function is part of the public developer API. You should not call it from
 user-facing code. You may use it when implementing new algorithms in third-party
 packages.
 """
-function minimize_multiobjective! end
+function minimize_multiobjective!(alg::AbstractAlgorithm, model::Optimizer)
+    return minimize_multiobjective!(alg, model, model.inner, model.f)
+end
 
 """
     optimize_multiobjective!(
diff --git a/src/algorithms/Chalmet.jl b/src/algorithms/Chalmet.jl
index 33e2c4b..401190c 100644
--- a/src/algorithms/Chalmet.jl
+++ b/src/algorithms/Chalmet.jl
@@ -26,66 +26,73 @@ This algorithm is restricted to problems with:
 struct Chalmet <: AbstractAlgorithm end
 
 function _solve_constrained_model(
-    model::Optimizer,
     ::Chalmet,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
     rhs::Vector{Float64},
 )
-    f = MOI.Utilities.scalarize(model.f)
-    g = sum(1.0 * fi for fi in f)
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(g)}(), g)
+    f_scalars = MOI.Utilities.scalarize(model.f)
+    g = MOI.Utilities.operate(+, Float64, f_scalars...)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(g)}(), g)
     sets = MOI.LessThan.(rhs .- 1)
-    c = MOI.Utilities.normalize_and_add_constraint.(model.inner, f, sets)
+    c = MOI.Utilities.normalize_and_add_constraint.(inner, f_scalars, sets)
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if !_is_scalar_status_optimal(status)
         _log_subproblem_solve(model, "subproblem not optimal")
         MOI.delete.(model, c)
         return status, nothing
     end
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    X, Y = _compute_point(model, variables, model.f)
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
+    X, Y = _compute_point(model, variables, f)
     _log_subproblem_solve(model, Y)
     MOI.delete.(model, c)
     return status, SolutionPoint(X, Y)
 end
 
-function minimize_multiobjective!(algorithm::Chalmet, model::Optimizer)
-    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
-    if MOI.output_dimension(model.f) != 2
+function minimize_multiobjective!(
+    algorithm::Chalmet,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
+    @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+    if MOI.output_dimension(f) != 2
         error("Chalmet requires exactly two objectives")
     end
     solutions = SolutionPoint[]
     E = Tuple{Int,Int}[]
     Q = Tuple{Int,Int}[]
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    f1, f2 = MOI.Utilities.scalarize(model.f)
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
+    f1, f2 = MOI.Utilities.scalarize(f)
     y1, y2 = zeros(2), zeros(2)
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if !_is_scalar_status_optimal(status)
         return status, solutions
     end
     _, y1[2] = _compute_point(model, variables, f2)
     _log_subproblem_solve(model, variables)
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f1)}(), f1)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(f1)}(), f1)
     y1_constraint = MOI.Utilities.normalize_and_add_constraint(
-        model.inner,
+        inner,
         f2,
         MOI.LessThan(y1[2]),
     )
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if !_is_scalar_status_optimal(status)
         return status, solutions
     end
     x1, y1[1] = _compute_point(model, variables, f1)
     _log_subproblem_solve(model, y1)
-    MOI.delete(model.inner, y1_constraint)
+    MOI.delete(inner, y1_constraint)
     push!(solutions, SolutionPoint(x1, y1))
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f1)}(), f1)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(f1)}(), f1)
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if !_is_scalar_status_optimal(status)
         return status, solutions
     end
@@ -94,20 +101,20 @@ function minimize_multiobjective!(algorithm::Chalmet, model::Optimizer)
     if y2[1] ≈ solutions[1].y[1]
         return MOI.OPTIMAL, solutions
     end
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
     y2_constraint = MOI.Utilities.normalize_and_add_constraint(
-        model.inner,
+        inner,
         f1,
         MOI.LessThan(y2[1]),
     )
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if !_is_scalar_status_optimal(status)
         return status, solutions
     end
     x2, y2[2] = _compute_point(model, variables, f2)
     _log_subproblem_solve(model, y2)
-    MOI.delete(model.inner, y2_constraint)
+    MOI.delete(inner, y2_constraint)
     push!(solutions, SolutionPoint(x2, y2))
     push!(Q, (1, 2))
     t = 3
@@ -118,7 +125,8 @@ function minimize_multiobjective!(algorithm::Chalmet, model::Optimizer)
         r, s = pop!(Q)
         yr, ys = solutions[r].y, solutions[s].y
         rhs = [max(yr[1], ys[1]), max(yr[2], ys[2])]
-        status, solution = _solve_constrained_model(model, algorithm, rhs)
+        status, solution =
+            _solve_constrained_model(algorithm, model, inner, f, rhs)
         if !_is_scalar_status_optimal(status)
             push!(E, (r, s))
             continue
diff --git a/src/algorithms/Dichotomy.jl b/src/algorithms/Dichotomy.jl
index 7332c86..c6bd523 100644
--- a/src/algorithms/Dichotomy.jl
+++ b/src/algorithms/Dichotomy.jl
@@ -61,37 +61,49 @@ function MOI.get(alg::Dichotomy, attr::SolutionLimit)
 end
 
 function _solve_weighted_sum(
-    model::Optimizer,
     ::Dichotomy,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
     weights::Vector{Float64},
 )
-    f = _scalarise(model.f, weights)
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f)}(), f)
+    f_scalar = _scalarise(f, weights)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(f_scalar)}(), f_scalar)
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if !_is_scalar_status_optimal(status)
         return status, nothing
     end
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    X, Y = _compute_point(model, variables, model.f)
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
+    X, Y = _compute_point(model, variables, f)
     _log_subproblem_solve(model, Y)
     return status, SolutionPoint(X, Y)
 end
 
 function optimize_multiobjective!(algorithm::Dichotomy, model::Optimizer)
-    if MOI.output_dimension(model.f) > 2
-        error("Only scalar or bi-objective problems supported.")
-    end
-    if MOI.output_dimension(model.f) == 1
+    return optimize_multiobjective!(algorithm, model, model.inner, model.f)
+end
+
+function optimize_multiobjective!(
+    algorithm::Dichotomy,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
+    if MOI.output_dimension(f) == 1
         if (ret = _check_premature_termination(model)) !== nothing
             return ret, nothing
         end
-        status, solution = _solve_weighted_sum(model, algorithm, [1.0])
+        status, solution =
+            _solve_weighted_sum(algorithm, model, inner, f, [1.0])
         return status, [solution]
+    elseif MOI.output_dimension(f) > 2
+        error("Only scalar or bi-objective problems supported.")
     end
     solutions = Dict{Float64,SolutionPoint}()
     for (i, w) in (1 => 1.0, 2 => 0.0)
-        status, solution = _solve_weighted_sum(model, algorithm, [w, 1.0 - w])
+        status, solution =
+            _solve_weighted_sum(algorithm, model, inner, f, [w, 1.0 - w])
         if !_is_scalar_status_optimal(status)
             return status, nothing
         end
@@ -113,7 +125,8 @@ function optimize_multiobjective!(algorithm::Dichotomy, model::Optimizer)
         (a, b) = popfirst!(queue)
         y_d = solutions[a].y .- solutions[b].y
         w = y_d[2] / (y_d[2] - y_d[1])
-        status, solution = _solve_weighted_sum(model, algorithm, [w, 1.0 - w])
+        status, solution =
+            _solve_weighted_sum(algorithm, model, inner, f, [w, 1.0 - w])
         if !_is_scalar_status_optimal(status)
             break # Exit the solve with some error.
         elseif solution ≈ solutions[a] || solution ≈ solutions[b]
diff --git a/src/algorithms/DominguezRios.jl b/src/algorithms/DominguezRios.jl
index 64ec5a8..78aca8b 100644
--- a/src/algorithms/DominguezRios.jl
+++ b/src/algorithms/DominguezRios.jl
@@ -152,19 +152,24 @@ _isapprox(::Nothing, ::_DominguezRiosBox) = false
 
 _isapprox(A::_DominguezRiosBox, B::_DominguezRiosBox) = A.l ≈ B.l && A.u ≈ B.u
 
-function minimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
-    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+function minimize_multiobjective!(
+    algorithm::DominguezRios,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
+    @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
     n = MOI.output_dimension(model.f)
     L = [_DominguezRiosBox[] for i in 1:n]
     scalars = MOI.Utilities.scalarize(model.f)
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
     yI, yN = zeros(n), zeros(n)
     # Ideal and Nadir point estimation
     for (i, f_i) in enumerate(scalars)
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
-        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
+        MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
         optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !_is_scalar_status_optimal(status)
             return status, nothing
         end
@@ -172,9 +177,9 @@ function minimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
         _log_subproblem_solve(model, variables)
         yI[i] = Y
         model.ideal_point[i] = Y
-        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+        MOI.set(inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
         optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !_is_scalar_status_optimal(status)
             _warn_on_nonfinite_anti_ideal(algorithm, MOI.MIN_SENSE, i)
             return status, nothing
@@ -183,14 +188,14 @@ function minimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
         _log_subproblem_solve(model, variables)
         yN[i] = Y + 1
     end
-    MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
     ϵ = 1 / (2 * n * (maximum(yN - yI) - 1))
     # If ϵ is small, then the scalar objectives can contain terms that fall
     # below the tolerance level of the solver. To fix this, we rescale the
     # objective so that the coefficients have magnitude `1e+00` or larger.
     scale = max(1.0, 1 / ϵ)
     push!(L[1], _DominguezRiosBox(yI, yN, 0.0))
-    t_max = MOI.add_variable(model.inner)
+    t_max = MOI.add_variable(inner)
     solutions = SolutionPoint[]
     k = 0
     status = MOI.OPTIMAL
@@ -215,7 +220,7 @@ function minimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
         w = scale ./ max.(1, B.u - yI)
         constraints = [
             MOI.Utilities.normalize_and_add_constraint(
-                model.inner,
+                inner,
                 # `w` is the scaled version here. This epigraph constraint will
                 # make `t_max` similarly scaled.
                 t_max - (w[i] * (scalars[i] - yI[i])),
@@ -226,12 +231,12 @@ function minimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
         #  * t_max is already covered in `constraints`
         #  * the `ϵ` term is already covered in `w`
         new_f = t_max + ϵ * sum(w[i] * (scalars[i] - yI[i]) for i in 1:n)
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
         optimize_inner!(model)
         if _is_scalar_status_optimal(model)
             X, Y = _compute_point(model, variables, model.f)
             _log_subproblem_solve(model, Y)
-            obj = MOI.get(model.inner, MOI.ObjectiveValue())
+            obj = MOI.get(inner, MOI.ObjectiveValue())
             # We need to undo the scaling of the scalar objective. There's no
             # need to unscale `Y` because we have evaluated this explicitly from
             # the modified `model.f`.
@@ -251,8 +256,8 @@ function minimize_multiobjective!(algorithm::DominguezRios, model::Optimizer)
             _log_subproblem_solve(model, "subproblem not optimal")
             deleteat!(L[k], i)
         end
-        MOI.delete.(model.inner, constraints)
+        MOI.delete.(inner, constraints)
     end
-    MOI.delete(model.inner, t_max)
+    MOI.delete(inner, t_max)
     return status, solutions
 end
diff --git a/src/algorithms/EpsilonConstraint.jl b/src/algorithms/EpsilonConstraint.jl
index a8bad83..6442cea 100644
--- a/src/algorithms/EpsilonConstraint.jl
+++ b/src/algorithms/EpsilonConstraint.jl
@@ -74,20 +74,22 @@ end
 function minimize_multiobjective!(
     algorithm::EpsilonConstraint,
     model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
 )
-    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
-    if MOI.output_dimension(model.f) != 2
+    @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+    if MOI.output_dimension(f) != 2
         error("EpsilonConstraint requires exactly two objectives")
     end
     # Compute the bounding box of the objectives using Hierarchical().
     alg = Hierarchical()
     MOI.set.(Ref(alg), ObjectivePriority.(1:2), [1, 0])
-    status, solution_1 = optimize_multiobjective!(alg, model)
+    status, solution_1 = minimize_multiobjective!(alg, model, inner, f)
     if !_is_scalar_status_optimal(status)
         return status, nothing
     end
     MOI.set(alg, ObjectivePriority(2), 2)
-    status, solution_2 = optimize_multiobjective!(alg, model)
+    status, solution_2 = minimize_multiobjective!(alg, model, inner, f)
     if !_is_scalar_status_optimal(status)
         return status, nothing
     end
@@ -101,10 +103,10 @@ function minimize_multiobjective!(
         ε = abs(right - left) / (n_points - 1)
     end
     solutions = SolutionPoint[only(solution_1), only(solution_2)]
-    f1, f2 = MOI.Utilities.eachscalar(model.f)
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
+    f1, f2 = MOI.Utilities.eachscalar(f)
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(f2)}(), f2)
     # Add epsilon constraint
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
     bound = right - ε
     constant = MOI.constant(f1, Float64)
     ci = MOI.Utilities.normalize_and_add_constraint(
@@ -125,7 +127,7 @@ function minimize_multiobjective!(
         if !_is_scalar_status_optimal(model)
             break
         end
-        X, Y = _compute_point(model, variables, model.f)
+        X, Y = _compute_point(model, variables, f)
         _log_subproblem_solve(model, Y)
         if isempty(solutions) || !(Y ≈ solutions[end].y)
             push!(solutions, SolutionPoint(X, Y))
diff --git a/src/algorithms/Hierarchical.jl b/src/algorithms/Hierarchical.jl
index 33aabc4..ca18b5b 100644
--- a/src/algorithms/Hierarchical.jl
+++ b/src/algorithms/Hierarchical.jl
@@ -89,11 +89,16 @@ function _sorted_priorities(priorities::Vector{Int})
     return [findall(isequal(u), priorities) for u in unique_priorities]
 end
 
-function minimize_multiobjective!(algorithm::Hierarchical, model::Optimizer)
-    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
-    objectives = MOI.Utilities.eachscalar(model.f)
+function minimize_multiobjective!(
+    algorithm::Hierarchical,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
+    @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+    objectives = MOI.Utilities.eachscalar(f)
     N = length(objectives)
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
     # Find list of objectives with same priority
     constraints = Any[]
     priorities = [MOI.get(algorithm, ObjectivePriority(i)) for i in 1:N]
@@ -103,9 +108,9 @@ function minimize_multiobjective!(algorithm::Hierarchical, model::Optimizer)
         # Solve weighted sum
         new_vector_f = objectives[indices]
         new_f = _scalarise(new_vector_f, weights[indices])
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(new_f)}(), new_f)
         optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !_is_scalar_status_optimal(status)
             return status, nothing
         end
@@ -113,7 +118,7 @@ function minimize_multiobjective!(algorithm::Hierarchical, model::Optimizer)
             break
         end
         if !model.silent
-            X, Y = _compute_point(model, variables, model.f)
+            X, Y = _compute_point(model, variables, f)
             _log_subproblem_solve(model, Y)
         end
         # Add tolerance constraints
@@ -125,7 +130,7 @@ function minimize_multiobjective!(algorithm::Hierarchical, model::Optimizer)
             push!(constraints, ci)
         end
     end
-    X, Y = _compute_point(model, variables, model.f)
+    X, Y = _compute_point(model, variables, f)
     _log_subproblem_solve(model, Y)
     # Remove tolerance constraints
     for c in constraints
diff --git a/src/algorithms/KirlikSayin.jl b/src/algorithms/KirlikSayin.jl
index 5e7ee0d..0637f03 100644
--- a/src/algorithms/KirlikSayin.jl
+++ b/src/algorithms/KirlikSayin.jl
@@ -81,26 +81,31 @@ function _update_list(L::Vector{_Rectangle}, f::Vector{Float64})
     return L_new
 end
 
-function minimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
-    @assert MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+function minimize_multiobjective!(
+    algorithm::KirlikSayin,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
+    @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
     solutions = SolutionPoint[]
     # Problem with p objectives.
     # Set k = 1, meaning the nondominated points will get projected
     # down to the objective {2, 3, ..., p}
     k = 1
     YN = Vector{Float64}[]
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    n = MOI.output_dimension(model.f)
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
+    n = MOI.output_dimension(f)
     yI, yN = zeros(n), zeros(n)
     # This tolerance is really important!
     δ = 1.0
-    scalars = MOI.Utilities.scalarize(model.f)
+    scalars = MOI.Utilities.scalarize(f)
     # Ideal and Nadir point estimation
     for (i, f_i) in enumerate(scalars)
         # Ideal point
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
         optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !_is_scalar_status_optimal(status)
             return status, nothing
         end
@@ -108,19 +113,19 @@ function minimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
         _log_subproblem_solve(model, variables)
         model.ideal_point[i] = yI[i] = Y
         # Nadir point
-        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+        MOI.set(inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
         optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
         if !_is_scalar_status_optimal(status)
             # Repair ObjectiveSense before exiting
-            MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+            MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
             _warn_on_nonfinite_anti_ideal(algorithm, MOI.MIN_SENSE, i)
             return status, nothing
         end
         _, Y = _compute_point(model, variables, f_i)
         _log_subproblem_solve(model, variables)
         yN[i] = Y + δ
-        MOI.set(model.inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+        MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
     end
     L = [_Rectangle(_project(yI, k), _project(yN, k))]
     status = MOI.OPTIMAL
@@ -135,18 +140,14 @@ function minimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
         #   minimize: f_1(x)
         #       s.t.: f_i(x) <= u_i - δ
         @assert k == 1
-        MOI.set(
-            model.inner,
-            MOI.ObjectiveFunction{typeof(scalars[k])}(),
-            scalars[k],
-        )
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(scalars[k])}(), scalars[k])
         ε_constraints = Any[]
         for (i, f_i) in enumerate(scalars)
             if i == k
                 continue
             end
             ci = MOI.Utilities.normalize_and_add_constraint(
-                model.inner,
+                inner,
                 f_i,
                 MOI.LessThan{Float64}(uᵢ[i-1] - δ),
             )
@@ -161,14 +162,14 @@ function minimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
             MOI.delete.(model, ε_constraints)
             continue
         end
-        zₖ = MOI.get(model.inner, MOI.ObjectiveValue())
+        zₖ = MOI.get(inner, MOI.ObjectiveValue())
         # Solving the second stage model: Q_k(ε, zₖ)
-        # Set objective sum(model.f)
+        # Set objective sum(f)
         sum_f = MOI.Utilities.operate(+, Float64, scalars...)
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(sum_f)}(), sum_f)
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(sum_f)}(), sum_f)
         # Constraint to eliminate weak dominance
         zₖ_constraint = MOI.Utilities.normalize_and_add_constraint(
-            model.inner,
+            inner,
             scalars[k],
             MOI.EqualTo(zₖ),
         )
@@ -182,7 +183,7 @@ function minimize_multiobjective!(algorithm::KirlikSayin, model::Optimizer)
             _remove_rectangle(L, _Rectangle(_project(yI, k), uᵢ))
             continue
         end
-        X, Y = _compute_point(model, variables, model.f)
+        X, Y = _compute_point(model, variables, f)
         _log_subproblem_solve(model, Y)
         Y_proj = _project(Y, k)
         if !(Y in YN)
diff --git a/src/algorithms/Lexicographic.jl b/src/algorithms/Lexicographic.jl
index 821827a..0edad0a 100644
--- a/src/algorithms/Lexicographic.jl
+++ b/src/algorithms/Lexicographic.jl
@@ -69,10 +69,19 @@ function MOI.set(alg::Lexicographic, ::LexicographicAllPermutations, val::Bool)
 end
 
 function optimize_multiobjective!(algorithm::Lexicographic, model::Optimizer)
+    return optimize_multiobjective!(algorithm, model, model.inner, model.f)
+end
+
+function optimize_multiobjective!(
+    algorithm::Lexicographic,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
     sequence = 1:MOI.output_dimension(model.f)
     perm = MOI.get(algorithm, LexicographicAllPermutations())
     if !something(perm, _default(LexicographicAllPermutations()))
-        return _solve_in_sequence(algorithm, model, sequence)
+        return _solve_in_sequence(algorithm, model, inner, f, sequence)
     end
     if perm === nothing && length(sequence) >= 5
         o, n = length(sequence), factorial(length(sequence))
@@ -102,7 +111,8 @@ function optimize_multiobjective!(algorithm::Lexicographic, model::Optimizer)
     solutions = SolutionPoint[]
     status = MOI.OPTIMAL
     for sequence in Combinatorics.permutations(sequence)
-        status, solution = _solve_in_sequence(algorithm, model, sequence)
+        status, solution =
+            _solve_in_sequence(algorithm, model, inner, f, sequence)
         if !isempty(solution)
             push!(solutions, solution[1])
         end
@@ -116,11 +126,13 @@ end
 function _solve_in_sequence(
     algorithm::Lexicographic,
     model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
     sequence::AbstractVector{Int},
 )
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
     constraints = Any[]
-    scalars = MOI.Utilities.eachscalar(model.f)
+    scalars = MOI.Utilities.eachscalar(f)
     solution = SolutionPoint[]
     status = MOI.OPTIMAL
     for i in sequence
@@ -128,27 +140,27 @@ function _solve_in_sequence(
             status = ret
             break
         end
-        f = scalars[i]
-        MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f)}(), f)
+        scalar_f = scalars[i]
+        MOI.set(inner, MOI.ObjectiveFunction{typeof(scalar_f)}(), scalar_f)
         optimize_inner!(model)
-        status = MOI.get(model.inner, MOI.TerminationStatus())
-        primal_status = MOI.get(model.inner, MOI.PrimalStatus())
+        status = MOI.get(inner, MOI.TerminationStatus())
+        primal_status = MOI.get(inner, MOI.PrimalStatus())
         if _is_scalar_status_feasible_point(primal_status)
-            X, Y = _compute_point(model, variables, model.f)
+            X, Y = _compute_point(model, variables, f)
             _log_subproblem_solve(model, Y)
             solution = [SolutionPoint(X, Y)]
         end
         if !_is_scalar_status_optimal(status)
             break
         end
-        X, Y = _compute_point(model, variables, f)
+        X, Y = _compute_point(model, variables, scalar_f)
         rtol = MOI.get(algorithm, ObjectiveRelativeTolerance(i))
-        set = if MOI.get(model.inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
+        set = if MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
             MOI.LessThan(Y + rtol * abs(Y))
         else
             MOI.GreaterThan(Y - rtol * abs(Y))
         end
-        ci = MOI.Utilities.normalize_and_add_constraint(model, f, set)
+        ci = MOI.Utilities.normalize_and_add_constraint(model, scalar_f, set)
         push!(constraints, ci)
     end
     for c in constraints
diff --git a/src/algorithms/RandomWeighting.jl b/src/algorithms/RandomWeighting.jl
index d9d287d..9f5ab47 100644
--- a/src/algorithms/RandomWeighting.jl
+++ b/src/algorithms/RandomWeighting.jl
@@ -39,20 +39,29 @@ function MOI.get(alg::RandomWeighting, attr::SolutionLimit)
 end
 
 function optimize_multiobjective!(algorithm::RandomWeighting, model::Optimizer)
+    return optimize_multiobjective!(algorithm, model, model.inner, model.f)
+end
+
+function optimize_multiobjective!(
+    algorithm::RandomWeighting,
+    model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
+)
     if MOI.get(model, MOI.TimeLimitSec()) === nothing &&
        algorithm.solution_limit === nothing
         error("At least `MOI.TimeLimitSec` or `MOA.SolutionLimit` must be set")
     end
     solutions = SolutionPoint[]
     sense = MOI.get(model, MOI.ObjectiveSense())
-    P = MOI.output_dimension(model.f)
-    variables = MOI.get(model.inner, MOI.ListOfVariableIndices())
-    f = _scalarise(model.f, ones(P))
-    MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f)}(), f)
+    P = MOI.output_dimension(f)
+    variables = MOI.get(inner, MOI.ListOfVariableIndices())
+    scalar_f = _scalarise(f, ones(P))
+    MOI.set(inner, MOI.ObjectiveFunction{typeof(scalar_f)}(), scalar_f)
     optimize_inner!(model)
-    status = MOI.get(model.inner, MOI.TerminationStatus())
+    status = MOI.get(inner, MOI.TerminationStatus())
     if _is_scalar_status_optimal(status)
-        X, Y = _compute_point(model, variables, model.f)
+        X, Y = _compute_point(model, variables, f)
         _log_subproblem_solve(model, Y)
         push!(solutions, SolutionPoint(X, Y))
     else
@@ -69,12 +78,12 @@ function optimize_multiobjective!(algorithm::RandomWeighting, model::Optimizer)
                 return ret, filter_nondominated(sense, solutions)
             end
             weights = rand(P)
-            f = _scalarise(model.f, weights)
-            MOI.set(model.inner, MOI.ObjectiveFunction{typeof(f)}(), f)
+            scalar_f = _scalarise(f, weights)
+            MOI.set(inner, MOI.ObjectiveFunction{typeof(scalar_f)}(), scalar_f)
             optimize_inner!(model)
-            status = MOI.get(model.inner, MOI.TerminationStatus())
+            status = MOI.get(inner, MOI.TerminationStatus())
             if _is_scalar_status_optimal(status)
-                X, Y = _compute_point(model, variables, model.f)
+                X, Y = _compute_point(model, variables, f)
                 _log_subproblem_solve(model, Y)
                 push!(solutions, SolutionPoint(X, Y))
             end
diff --git a/src/algorithms/TambyVanderpooten.jl b/src/algorithms/TambyVanderpooten.jl
index d0436d5..c6d33a5 100644
--- a/src/algorithms/TambyVanderpooten.jl
+++ b/src/algorithms/TambyVanderpooten.jl
@@ -123,13 +123,15 @@ end
 function minimize_multiobjective!(
     algorithm::TambyVanderpooten,
     model::Optimizer,
+    inner::MOI.ModelLike,
+    f::MOI.AbstractVectorFunction,
 )
     solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}()
     status = _minimize_multiobjective!(
         algorithm,
         model,
-        model.inner,
-        model.f,
+        inner,
+        f,
         solutions,
     )::MOI.TerminationStatusCode
     return status, SolutionPoint[SolutionPoint(X, Y) for (Y, X) in solutions]