diff --git a/src/MultiObjectiveAlgorithms.jl b/src/MultiObjectiveAlgorithms.jl
index 3a9a68d..4447c85 100644
--- a/src/MultiObjectiveAlgorithms.jl
+++ b/src/MultiObjectiveAlgorithms.jl
@@ -34,8 +34,10 @@ function dominates(sense, a::SolutionPoint, b::SolutionPoint)
     end
 end
 
+_sort!(solutions::Vector{SolutionPoint}) = sort!(solutions; by = x -> x.y)
+
 function filter_nondominated(sense, solutions::Vector{SolutionPoint})
-    solutions = sort(solutions; by = x -> x.y)
+    _sort!(solutions)
     nondominated_solutions = SolutionPoint[]
     for candidate in solutions
         if any(test -> dominates(sense, test, candidate), solutions)
@@ -594,6 +596,7 @@ function MOI.optimize!(model::Optimizer)
     model.termination_status = status
     if solutions !== nothing
         model.solutions = solutions
+        _sort!(model.solutions)
     end
     if MOI.get(model, ComputeIdealPoint())
         _compute_ideal_point(model, start_time)
diff --git a/test/algorithms/Chalmet.jl b/test/algorithms/Chalmet.jl
index 2e54f52..8ffefb7 100644
--- a/test/algorithms/Chalmet.jl
+++ b/test/algorithms/Chalmet.jl
@@ -53,22 +53,17 @@ function test_knapsack_min()
     MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
     MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
     MOI.optimize!(model)
-    X_E = Float64[
-        0 0 1 1 1 0 1 1 1 1
-        1 0 1 1 1 0 1 1 0 1
-        0 1 1 1 1 0 1 0 1 1
+    results = [
+        [1, 0, 1, 1, 1, 0, 1, 1, 0, 1] => [-3394, -3817],
+        [0, 1, 1, 1, 1, 0, 1, 0, 1, 1] => [-3042, -4627],
+        [0, 0, 1, 1, 1, 0, 1, 1, 1, 1] => [-2854, -4636],
     ]
-    Y_N = Float64[
-        -2854 -4636
-        -3394 -3817
-        -3042 -4627
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(x_sol, X_E'; atol = 1e-6)
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)
-    @test isapprox(y_sol, Y_N'; atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
+    @test MOI.get(model, MOI.ResultCount()) == length(results)
+    for (i, (x_sol, y_sol)) in enumerate(results)
+        @test ≈(x_sol, MOI.get(model, MOI.VariablePrimal(i), x); atol = 1e-6)
+        @test ≈(y_sol, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ [-3394, -4636]
     return
 end
 
@@ -103,22 +98,17 @@ function test_knapsack_max()
     MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
     MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
     MOI.optimize!(model)
-    X_E = Float64[
-        0 0 1 1 1 0 1 1 1 1
-        1 0 1 1 1 0 1 1 0 1
-        0 1 1 1 1 0 1 0 1 1
-    ]
-    Y_N = Float64[
-        2855 4636
-        3395 3817
-        3043 4627
+    results = [
+        [0, 0, 1, 1, 1, 0, 1, 1, 1, 1] => [2855, 4636],
+        [0, 1, 1, 1, 1, 0, 1, 0, 1, 1] => [3043, 4627],
+        [1, 0, 1, 1, 1, 0, 1, 1, 0, 1] => [3395, 3817],
     ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(x_sol, X_E'; atol = 1e-6)
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)
-    @test isapprox(y_sol, Y_N'; atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
+    @test MOI.get(model, MOI.ResultCount()) == length(results)
+    for (i, (x_sol, y_sol)) in enumerate(results)
+        @test ≈(x_sol, MOI.get(model, MOI.VariablePrimal(i), x); atol = 1e-6)
+        @test ≈(y_sol, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ [3395, 4636]
     return
 end
 
@@ -228,6 +218,6 @@ function test_too_many_objectives()
     return
 end
 
-end
+end  # module TestChalmet
 
 TestChalmet.run_tests()
diff --git a/test/algorithms/Dichotomy.jl b/test/algorithms/Dichotomy.jl
index 9d4d41f..eeb486b 100644
--- a/test/algorithms/Dichotomy.jl
+++ b/test/algorithms/Dichotomy.jl
@@ -45,13 +45,13 @@ function test_moi_bolp_1()
     MOI.Utilities.loadfromstring!(
         model,
         """
-variables: x, y
-minobjective: [2 * x + y + 1, x + 3 * y]
-c1: x + y >= 1.0
-c2: 0.5 * x + y >= 0.75
-c3: x >= 0.0
-c4: y >= 0.25
-""",
+        variables: x, y
+        minobjective: [2 * x + y + 1, x + 3 * y]
+        c1: x + y >= 1.0
+        c2: 0.5 * x + y >= 0.75
+        c3: x >= 0.0
+        c4: y >= 0.25
+        """,
     )
     x = MOI.get(model, MOI.VariableIndex, "x")
     y = MOI.get(model, MOI.VariableIndex, "y")
@@ -81,21 +81,21 @@ function test_moi_bolp_1_maximize()
     MOI.Utilities.loadfromstring!(
         model,
         """
-variables: x, y
-maxobjective: [-2.0 * x + -1.0 * y, -1.0 * x + -3.0 * y + 0.5]
-c1: x + y >= 1.0
-c2: 0.5 * x + y >= 0.75
-c3: x >= 0.0
-c4: y >= 0.25
-""",
+        variables: x, y
+        maxobjective: [-2.0 * x + -1.0 * y, -1.0 * x + -3.0 * y + 0.5]
+        c1: x + y >= 1.0
+        c2: 0.5 * x + y >= 0.75
+        c3: x >= 0.0
+        c4: y >= 0.25
+        """,
     )
     x = MOI.get(model, MOI.VariableIndex, "x")
     y = MOI.get(model, MOI.VariableIndex, "y")
     MOI.optimize!(model)
     @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
     @test MOI.get(model, MOI.ResultCount()) == 3
-    X = [[0.0, 1.0], [0.5, 0.5], [1.0, 0.25]]
-    Y = [-[1.0, 2.5], -[1.5, 1.5], -[2.25, 1.25]]
+    X = [[1.0, 0.25], [0.5, 0.5], [0.0, 1.0]]
+    Y = [[-2.25, -1.25], [-1.5, -1.5], [-1.0, -2.5]]
     for i in 1:3
         @test MOI.get(model, MOI.PrimalStatus(i)) == MOI.FEASIBLE_POINT
         @test MOI.get(model, MOI.DualStatus(i)) == MOI.NO_SOLUTION
@@ -117,13 +117,13 @@ function test_moi_bolp_1_reversed()
     MOI.Utilities.loadfromstring!(
         model,
         """
-variables: x, y
-minobjective: [x + 3 * y, 2 * x + y]
-c1: x + y >= 1.0
-c2: 0.5 * x + y >= 0.75
-c3: x >= 0.0
-c4: y >= 0.25
-""",
+        variables: x, y
+        minobjective: [x + 3 * y, 2 * x + y]
+        c1: x + y >= 1.0
+        c2: 0.5 * x + y >= 0.75
+        c3: x >= 0.0
+        c4: y >= 0.25
+        """,
     )
     x = MOI.get(model, MOI.VariableIndex, "x")
     y = MOI.get(model, MOI.VariableIndex, "y")
@@ -153,13 +153,13 @@ function test_moi_bolp_1_scalar()
     MOI.Utilities.loadfromstring!(
         model,
         """
-variables: x, y
-minobjective: [2 * x + y, x + 3 * y]
-c1: x + y >= 1.0
-c2: 0.5 * x + y >= 0.75
-c3: x >= 0.0
-c4: y >= 0.25
-""",
+        variables: x, y
+        minobjective: [2 * x + y, x + 3 * y]
+        c1: x + y >= 1.0
+        c2: 0.5 * x + y >= 0.75
+        c3: x >= 0.0
+        c4: y >= 0.25
+        """,
     )
     x = MOI.get(model, MOI.VariableIndex, "x")
     y = MOI.get(model, MOI.VariableIndex, "y")
@@ -219,17 +219,16 @@ function test_biobjective_knapsack()
         MOI.LessThan(900.0),
     )
     MOI.optimize!(model)
-    results = Dict(
-        [955.0, 906.0] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
-        [948.0, 939.0] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
-        [934.0, 971.0] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
+    results = [
         [918.0, 983.0] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
-    )
+        [934.0, 971.0] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
+        [948.0, 939.0] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
+        [955.0, 906.0] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
+    ]
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
-        X = findall(elt -> elt > 0.9, x_sol)
-        Y = MOI.get(model, MOI.ObjectiveValue(i))
-        @test results[Y] == X
+        @test results[i][2] == findall(elt -> elt > 0.9, x_sol)
+        @test results[i][1] ≈ MOI.get(model, MOI.ObjectiveValue(i))
     end
     return
 end
@@ -353,9 +352,9 @@ function test_three_objective()
     MOI.Utilities.loadfromstring!(
         model,
         """
-variables: x
-maxobjective: [1.0 * x, -1.0 * x, 2.0 * x + 2.0]
-""",
+        variables: x
+        maxobjective: [1.0 * x, -1.0 * x, 2.0 * x + 2.0]
+        """,
     )
     @test_throws(
         ErrorException("Only scalar or bi-objective problems supported."),
@@ -408,6 +407,6 @@ function test_vector_of_variables_objective()
     return
 end
 
-end
+end  # module TestDichotomy
 
 TestDichotomy.run_tests()
diff --git a/test/algorithms/DominguezRios.jl b/test/algorithms/DominguezRios.jl
index f79efed..fab67d1 100644
--- a/test/algorithms/DominguezRios.jl
+++ b/test/algorithms/DominguezRios.jl
@@ -11,6 +11,8 @@ import HiGHS
 import MultiObjectiveAlgorithms as MOA
 import MultiObjectiveAlgorithms: MOI
 
+include(joinpath(dirname(@__DIR__), "problems.jl"))
+
 function run_tests()
     if Sys.WORD_SIZE == 32
         return  # Skip on 32-bit because HiGHS fails
@@ -25,478 +27,11 @@ function run_tests()
     return
 end
 
-function test_knapsack_min_p3()
-    p = 3
-    n = 10
-    W = 2137.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-    ]
-    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
+function test_problems()
     model = MOA.Optimizer(HiGHS.Optimizer)
     MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
     MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
-            for i in 1:p for j in 1:n
-        ],
-        ones(p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 1 1 1 0 1 0 1 1
-        0 0 1 1 1 0 1 1 1 1
-        0 1 1 1 1 0 0 1 0 1
-        1 0 1 1 1 0 1 1 0 1
-        1 1 1 1 1 0 0 0 1 0
-        1 1 1 1 1 0 0 0 0 1
-        1 0 1 1 1 0 0 1 1 0
-    ]
-    Y_N = Float64[
-        -3042 -4627 -3189
-        -2854 -4636 -3076
-        -2854 -3570 -3714
-        -3394 -3817 -3408
-        -2706 -3857 -3304
-        -2997 -3539 -3509
-        -2518 -3866 -3191
-    ]
-    Y_N .+= 1
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_max_p3()
-    p = 3
-    n = 10
-    W = 2137.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-    ]
-    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
-            i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 1 1 1 0 1 0 1 1
-        0 0 1 1 1 0 1 1 1 1
-        0 1 1 1 1 0 0 1 0 1
-        1 0 1 1 1 0 1 1 0 1
-        1 1 1 1 1 0 0 0 1 0
-        1 1 1 1 1 0 0 0 0 1
-        1 0 1 1 1 0 0 1 1 0
-    ]
-    Y_N = Float64[
-        3042 4627 3189
-        2854 4636 3076
-        2854 3570 3714
-        3394 3817 3408
-        2706 3857 3304
-        2997 3539 3509
-        2518 3866 3191
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_min_p4()
-    p = 4
-    n = 10
-    W = 2653.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-        557 898 148 63 78 964 246 662 386 272
-    ]
-    w = Float64[979 448 355 955 426 229 9 695 322 889]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
-            for i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 1 0 1 1 1 1 1 0
-        0 1 1 1 1 0 1 0 1 0
-        0 1 1 0 1 0 1 0 1 1
-        0 0 1 0 1 1 1 0 1 1
-        1 1 1 0 1 1 1 0 0 0
-        0 0 1 1 0 1 1 1 1 0
-        0 1 1 1 1 1 1 0 0 0
-        0 1 1 1 0 1 1 0 1 0
-        0 1 0 0 0 1 1 1 1 1
-        0 0 1 1 1 1 1 0 1 0
-    ]
-    Y_N = Float64[
-        -3152 -3232 -3596 -3382
-        -2725 -4064 -2652 -1819
-        -2862 -3648 -3049 -2028
-        -2435 -3618 -2282 -2094
-        -3269 -2320 -3059 -2891
-        -1904 -3253 -2530 -2469
-        -2883 -3237 -2535 -2397
-        -2092 -3244 -2643 -2705
-        -2146 -1944 -2947 -3428
-        -2298 -4034 -1885 -1885
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_max_p4()
-    p = 4
-    n = 10
-    W = 2653.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-        557 898 148 63 78 964 246 662 386 272
-    ]
-    w = Float64[979 448 355 955 426 229 9 695 322 889]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
-            i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 1 0 1 1 1 1 1 0
-        0 1 1 1 1 0 1 0 1 0
-        0 1 1 0 1 0 1 0 1 1
-        0 0 1 0 1 1 1 0 1 1
-        1 1 1 0 1 1 1 0 0 0
-        0 0 1 1 0 1 1 1 1 0
-        0 1 1 1 1 1 1 0 0 0
-        0 1 1 1 0 1 1 0 1 0
-        0 1 0 0 0 1 1 1 1 1
-        0 0 1 1 1 1 1 0 1 0
-    ]
-    Y_N = Float64[
-        3152 3232 3596 3382
-        2725 4064 2652 1819
-        2862 3648 3049 2028
-        2435 3618 2282 2094
-        3269 2320 3059 2891
-        1904 3253 2530 2469
-        2883 3237 2535 2397
-        2092 3244 2643 2705
-        2146 1944 2947 3428
-        2298 4034 1885 1885
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
-    return
-end
-
-function test_assignment_min_p3()
-    p = 3
-    n = 5
-    C = Float64[
-        6 1 20 2 3
-        2 6 9 10 18
-        1 6 20 5 9
-        6 8 6 9 6
-        7 10 10 6 2
-        17 20 8 8 20
-        10 13 1 10 15
-        4 11 1 13 1
-        19 13 7 18 17
-        15 3 5 1 11
-        10 7 1 19 12
-        2 15 12 10 3
-        11 20 16 12 9
-        10 15 20 11 7
-        1 9 20 7 6
-    ]
-    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
-    MOI.set(model, MOI.Silent(), true)
-    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    for i in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    for j in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0
-        1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
-        0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0
-    ]
-    Y_N = Float64[
-        23 43 44
-        38 33 53
-        40 47 37
-        20 52 54
-        45 33 34
-        43 51 31
-        28 33 58
-        29 29 59
-        35 49 39
-        50 40 32
-        16 61 47
-        37 55 36
-        39 43 41
-        18 47 67
-        24 39 45
-        28 66 39
-        34 60 42
-        22 37 63
-        22 54 47
-        17 43 71
-        35 38 56
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol =
-        hcat([MOI.get(model, MOI.VariablePrimal(i), vec(x)) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    return
-end
-
-function test_assignment_max_p3()
-    p = 3
-    n = 5
-    C = Float64[
-        6 1 20 2 3
-        2 6 9 10 18
-        1 6 20 5 9
-        6 8 6 9 6
-        7 10 10 6 2
-        17 20 8 8 20
-        10 13 1 10 15
-        4 11 1 13 1
-        19 13 7 18 17
-        15 3 5 1 11
-        10 7 1 19 12
-        2 15 12 10 3
-        11 20 16 12 9
-        10 15 20 11 7
-        1 9 20 7 6
-    ]
-    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
-    MOI.set(model, MOI.Silent(), true)
-    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    for i in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    for j in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(-C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0
-        1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
-        0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0
-    ]
-    Y_N = Float64[
-        -23 -43 -44
-        -38 -33 -53
-        -40 -47 -37
-        -20 -52 -54
-        -45 -33 -34
-        -43 -51 -31
-        -28 -33 -58
-        -29 -29 -59
-        -35 -49 -39
-        -50 -40 -32
-        -16 -61 -47
-        -37 -55 -36
-        -39 -43 -41
-        -18 -47 -67
-        -24 -39 -45
-        -28 -66 -39
-        -34 -60 -42
-        -22 -37 -63
-        -22 -54 -47
-        -17 -43 -71
-        -35 -38 -56
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol =
-        hcat([MOI.get(model, MOI.VariablePrimal(i), vec(x)) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    Problems.run_tests(model)
     return
 end
 
@@ -604,52 +139,6 @@ function test_vector_of_variables_objective()
     return
 end
 
-function test_issue_105()
-    cost = [100.0, 120.0, 150.0, 110.0, 200.0, 170.0]
-    time = [8.0, 3.0, 4.0, 2.0, 5.0, 4.0]
-    capacity = [10.0, 8.0]
-    demand = [5.0, 8.0, 5.0]
-    m, n = 2, 3
-    model = MOI.instantiate(; with_bridge_type = Float64) do
-        return MOA.Optimizer(HiGHS.Optimizer)
-    end
-    MOI.set(model, MOA.Algorithm(), MOA.DominguezRios())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, m * n)
-    MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
-    MOI.add_constraint.(model, x, MOI.Integer())
-    X = reshape(x, m, n)
-    for i in 1:m
-        f_i = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[i, :]), 0.0)
-        MOI.add_constraint(model, f_i, MOI.LessThan(capacity[i]))
-    end
-    for j in 1:n
-        f_j = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[:, j]), 0.0)
-        MOI.add_constraint(model, f_j, MOI.EqualTo(demand[j]))
-    end
-    f = MOI.Utilities.operate(
-        vcat,
-        Float64,
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(cost, x), 0.0),
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(time, x), 0.0),
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    @test MOI.get(model, MOI.ResultCount()) == 6
-    y = [MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:6]
-    @test sort!(y) == [
-        [2380.0, 81.0, 18.0],
-        [2440.0, 78.0, 18.0],
-        [2500.0, 75.0, 18.0],
-        [2560.0, 72.0, 18.0],
-        [2620.0, 69.0, 18.0],
-        [2680.0, 66.0, 18.0],
-    ]
-    return
-end
-
-end
+end  # module TestDominguezRios
 
 TestDominguezRios.run_tests()
diff --git a/test/algorithms/EpsilonConstraint.jl b/test/algorithms/EpsilonConstraint.jl
index b733c37..d8f978e 100644
--- a/test/algorithms/EpsilonConstraint.jl
+++ b/test/algorithms/EpsilonConstraint.jl
@@ -47,23 +47,23 @@ function test_biobjective_knapsack()
         MOI.LessThan(900.0),
     )
     MOI.optimize!(model)
-    results = Dict(
-        [956, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
-        [950, 915] => [1, 2, 5, 6, 8, 9, 10, 11, 15, 16, 17],
-        [949, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
-        [944, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
-        [937, 942] => [1, 2, 3, 5, 6, 10, 11, 12, 15, 16, 17],
-        [936, 947] => [2, 5, 6, 8, 9, 10, 11, 12, 15, 16, 17],
-        [935, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
-        [928, 972] => [2, 3, 5, 6, 8, 9, 10, 11, 12, 16, 17],
+    results = [
         [919, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
-    )
+        [928, 972] => [2, 3, 5, 6, 8, 9, 10, 11, 12, 16, 17],
+        [935, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
+        [936, 947] => [2, 5, 6, 8, 9, 10, 11, 12, 15, 16, 17],
+        [937, 942] => [1, 2, 3, 5, 6, 10, 11, 12, 15, 16, 17],
+        [944, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
+        [949, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
+        [950, 915] => [1, 2, 5, 6, 8, 9, 10, 11, 15, 16, 17],
+        [956, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
+    ]
     @test MOI.get(model, MOI.ResultCount()) == 9
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
         Y = MOI.get(model, MOI.ObjectiveValue(i))
-        @test results[round.(Int, Y)] == X
+        @test results[i] == (round.(Int, Y) => X)
     end
     @test MOI.get(model, MOI.ObjectiveBound()) == [956.0, 983.0]
     return
@@ -91,23 +91,23 @@ function test_biobjective_knapsack_atol()
         MOI.LessThan(900.0),
     )
     MOI.optimize!(model)
-    results = Dict(
-        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
-        [949, 915] => [1, 2, 5, 6, 8, 9, 10, 11, 15, 16, 17],
-        [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
-        [943, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
-        [936, 942] => [1, 2, 3, 5, 6, 10, 11, 12, 15, 16, 17],
-        [935, 947] => [2, 5, 6, 8, 9, 10, 11, 12, 15, 16, 17],
-        [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
-        [927, 972] => [2, 3, 5, 6, 8, 9, 10, 11, 12, 16, 17],
+    results = [
         [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
-    )
+        [927, 972] => [2, 3, 5, 6, 8, 9, 10, 11, 12, 16, 17],
+        [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
+        [935, 947] => [2, 5, 6, 8, 9, 10, 11, 12, 15, 16, 17],
+        [936, 942] => [1, 2, 3, 5, 6, 10, 11, 12, 15, 16, 17],
+        [943, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
+        [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
+        [949, 915] => [1, 2, 5, 6, 8, 9, 10, 11, 15, 16, 17],
+        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
+    ]
     @test MOI.get(model, MOI.ResultCount()) == 9
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
         Y = MOI.get(model, MOI.ObjectiveValue(i))
-        @test results[round.(Int, Y)] == X
+        @test results[i] == (round.(Int, Y) => X)
     end
     @test MOI.get(model, MOI.ObjectiveBound()) == [955.0, 983.0]
     return
@@ -137,18 +137,18 @@ function test_biobjective_knapsack_atol_large()
         MOI.LessThan(900.0),
     )
     MOI.optimize!(model)
-    results = Dict(
-        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
-        [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
-        [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
+    results = [
         [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
-    )
+        [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
+        [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
+        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
+    ]
     @test MOI.get(model, MOI.ResultCount()) == 4
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
         Y = MOI.get(model, MOI.ObjectiveValue(i))
-        @test results[round.(Int, Y)] == X
+        @test results[i] == (round.(Int, Y) => X)
     end
     return
 end
@@ -176,7 +176,7 @@ function test_biobjective_knapsack_min()
         MOI.LessThan(900.0),
     )
     MOI.optimize!(model)
-    results = Dict(
+    results = [
         [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
         [949, 915] => [1, 2, 5, 6, 8, 9, 10, 11, 15, 16, 17],
         [948, 939] => [1, 2, 3, 5, 6, 8, 10, 11, 15, 16, 17],
@@ -186,13 +186,13 @@ function test_biobjective_knapsack_min()
         [934, 971] => [2, 3, 5, 6, 8, 10, 11, 12, 15, 16, 17],
         [927, 972] => [2, 3, 5, 6, 8, 9, 10, 11, 12, 16, 17],
         [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
-    )
+    ]
     @test MOI.get(model, MOI.ResultCount()) == 9
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
         Y = MOI.get(model, MOI.ObjectiveValue(i))
-        @test results[-round.(Int, Y)] == X
+        @test results[i] == (-round.(Int, Y) => X)
     end
     @test MOI.get(model, MOI.ObjectiveBound()) == [-955.0, -983.0]
     return
@@ -222,17 +222,17 @@ function test_biobjective_knapsack_min_solution_limit()
         MOI.LessThan(900.0),
     )
     MOI.optimize!(model)
-    results = Dict(
-        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
-        [943, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
+    results = [
         [918, 983] => [2, 3, 4, 5, 6, 8, 10, 11, 12, 16, 17],
-    )
+        [943, 940] => [2, 3, 5, 6, 8, 9, 10, 11, 15, 16, 17],
+        [955, 906] => [2, 3, 5, 6, 9, 10, 11, 14, 15, 16, 17],
+    ]
     @test MOI.get(model, MOI.ResultCount()) == 3
     for i in 1:MOI.get(model, MOI.ResultCount())
         x_sol = MOI.get(model, MOI.VariablePrimal(i), x)
         X = findall(elt -> elt > 0.9, x_sol)
         Y = MOI.get(model, MOI.ObjectiveValue(i))
-        @test results[round.(Int, Y)] == X
+        @test results[i] == (round.(Int, Y) => X)
     end
     @test MOI.get(model, MOI.ObjectiveBound()) == [955.0, 983.0]
     return
@@ -491,6 +491,6 @@ function test_too_many_objectives()
     return
 end
 
-end
+end  # module TestEpsilonConstraint
 
 TestEpsilonConstraint.run_tests()
diff --git a/test/algorithms/KirlikSayin.jl b/test/algorithms/KirlikSayin.jl
index 5242b4c..077d57c 100644
--- a/test/algorithms/KirlikSayin.jl
+++ b/test/algorithms/KirlikSayin.jl
@@ -11,6 +11,8 @@ import HiGHS
 import MultiObjectiveAlgorithms as MOA
 import MultiObjectiveAlgorithms: MOI
 
+include(joinpath(dirname(@__DIR__), "problems.jl"))
+
 function run_tests()
     for name in names(@__MODULE__; all = true)
         if startswith("$name", "test_")
@@ -22,484 +24,11 @@ function run_tests()
     return
 end
 
-function test_knapsack_min_p3()
-    p = 3
-    n = 10
-    W = 2137.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-    ]
-    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
+function test_problems()
     model = MOA.Optimizer(HiGHS.Optimizer)
     MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
     MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
-            for i in 1:p for j in 1:n
-        ],
-        ones(p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 0 1 1 1 0 1 1 0 1
-        0 1 1 1 1 0 1 0 1 1
-        1 1 1 1 1 0 0 0 0 1
-        0 1 1 1 1 0 0 1 0 1
-        1 1 1 1 1 0 0 0 1 0
-        1 0 1 1 1 0 0 1 1 0
-        0 0 1 1 1 0 1 1 1 1
-    ]
-    Y_N = Float64[
-        -3394 -3817 -3408
-        -3042 -4627 -3189
-        -2997 -3539 -3509
-        -2854 -3570 -3714
-        -2706 -3857 -3304
-        -2518 -3866 -3191
-        -2854 -4636 -3076
-    ]
-    Y_N .+= 1
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_max_p3()
-    p = 3
-    n = 10
-    W = 2137.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-    ]
-    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
-            i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 0 1 1 1 0 1 1 0 1
-        0 1 1 1 1 0 1 0 1 1
-        1 1 1 1 1 0 0 0 0 1
-        0 1 1 1 1 0 0 1 0 1
-        1 1 1 1 1 0 0 0 1 0
-        1 0 1 1 1 0 0 1 1 0
-        0 0 1 1 1 0 1 1 1 1
-    ]
-    Y_N = Float64[
-        3394 3817 3408
-        3042 4627 3189
-        2997 3539 3509
-        2854 3570 3714
-        2706 3857 3304
-        2518 3866 3191
-        2854 4636 3076
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_min_p4()
-    p = 4
-    n = 10
-    W = 2653.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-        557 898 148 63 78 964 246 662 386 272
-    ]
-    w = Float64[979 448 355 955 426 229 9 695 322 889]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
-            for i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 1 1 0 1 1 1 0 0 0
-        0 1 1 0 1 1 1 1 1 0
-        0 1 1 1 1 1 1 0 0 0
-        0 1 1 0 1 0 1 0 1 1
-        0 1 1 1 1 0 1 0 1 0
-        0 0 1 0 1 1 1 0 1 1
-        0 1 1 1 0 1 1 0 1 0
-        0 0 1 1 0 1 1 1 1 0
-        0 0 1 1 1 1 1 0 1 0
-        0 1 0 0 0 1 1 1 1 1
-    ]
-    Y_N = Float64[
-        -3269 -2320 -3059 -2891
-        -3152 -3232 -3596 -3382
-        -2883 -3237 -2535 -2397
-        -2862 -3648 -3049 -2028
-        -2725 -4064 -2652 -1819
-        -2435 -3618 -2282 -2094
-        -2092 -3244 -2643 -2705
-        -1904 -3253 -2530 -2469
-        -2298 -4034 -1885 -1885
-        -2146 -1944 -2947 -3428
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_max_p4()
-    p = 4
-    n = 10
-    W = 2653.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-        557 898 148 63 78 964 246 662 386 272
-    ]
-    w = Float64[979 448 355 955 426 229 9 695 322 889]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
-            i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 1 1 0 1 1 1 0 0 0
-        0 1 1 0 1 1 1 1 1 0
-        0 1 1 1 1 1 1 0 0 0
-        0 1 1 0 1 0 1 0 1 1
-        0 1 1 1 1 0 1 0 1 0
-        0 0 1 0 1 1 1 0 1 1
-        0 1 1 1 0 1 1 0 1 0
-        0 0 1 1 0 1 1 1 1 0
-        0 0 1 1 1 1 1 0 1 0
-        0 1 0 0 0 1 1 1 1 1
-    ]
-    Y_N = Float64[
-        3269 2320 3059 2891
-        3152 3232 3596 3382
-        2883 3237 2535 2397
-        2862 3648 3049 2028
-        2725 4064 2652 1819
-        2435 3618 2282 2094
-        2092 3244 2643 2705
-        1904 3253 2530 2469
-        2298 4034 1885 1885
-        2146 1944 2947 3428
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
-    return
-end
-
-function test_assignment_min_p3()
-    if Sys.WORD_SIZE == 32
-        return  # Skip on 32-bit because HiGHS fails
-    end
-    p = 3
-    n = 5
-    C = Float64[
-        6 1 20 2 3
-        2 6 9 10 18
-        1 6 20 5 9
-        6 8 6 9 6
-        7 10 10 6 2
-        17 20 8 8 20
-        10 13 1 10 15
-        4 11 1 13 1
-        19 13 7 18 17
-        15 3 5 1 11
-        10 7 1 19 12
-        2 15 12 10 3
-        11 20 16 12 9
-        10 15 20 11 7
-        1 9 20 7 6
-    ]
-    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
-    MOI.set(model, MOI.Silent(), true)
-    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    for i in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    for j in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
-        0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0
-        0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0
-    ]
-    Y_N = Float64[
-        16 61 47
-        17 43 71
-        18 47 67
-        20 52 54
-        22 37 63
-        23 43 44
-        22 54 47
-        28 66 39
-        34 60 42
-        24 39 45
-        35 49 39
-        37 55 36
-        28 33 58
-        35 38 56
-        39 43 41
-        38 33 53
-        45 33 34
-        43 51 31
-        40 47 37
-        29 29 59
-        50 40 32
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol =
-        hcat([MOI.get(model, MOI.VariablePrimal(i), vec(x)) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
-    return
-end
-
-function test_assignment_max_p3()
-    if Sys.WORD_SIZE == 32
-        return  # Skip on 32-bit because HiGHS fails
-    end
-    p = 3
-    n = 5
-    C = Float64[
-        6 1 20 2 3
-        2 6 9 10 18
-        1 6 20 5 9
-        6 8 6 9 6
-        7 10 10 6 2
-        17 20 8 8 20
-        10 13 1 10 15
-        4 11 1 13 1
-        19 13 7 18 17
-        15 3 5 1 11
-        10 7 1 19 12
-        2 15 12 10 3
-        11 20 16 12 9
-        10 15 20 11 7
-        1 9 20 7 6
-    ]
-    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
-    MOI.set(model, MOI.Silent(), true)
-    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    for i in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    for j in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(-C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
-        0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0
-        0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0
-    ]
-    Y_N = Float64[
-        -16 -61 -47
-        -17 -43 -71
-        -18 -47 -67
-        -20 -52 -54
-        -22 -37 -63
-        -23 -43 -44
-        -22 -54 -47
-        -28 -66 -39
-        -34 -60 -42
-        -24 -39 -45
-        -35 -49 -39
-        -37 -55 -36
-        -28 -33 -58
-        -35 -38 -56
-        -39 -43 -41
-        -38 -33 -53
-        -45 -33 -34
-        -43 -51 -31
-        -40 -47 -37
-        -29 -29 -59
-        -50 -40 -32
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol =
-        hcat([MOI.get(model, MOI.VariablePrimal(i), vec(x)) for i in 1:N]...)
-    @test isapprox(sort(x_sol; dims = 1), sort(X_E'; dims = 1); atol = 1e-6)
-    y_sol = vcat([MOI.get(model, MOI.ObjectiveValue(i))' for i in 1:N]...)
-    @test isapprox(sort(y_sol; dims = 1), sort(Y_N; dims = 1); atol = 1e-6)
+    Problems.run_tests(model)
     return
 end
 
@@ -608,52 +137,6 @@ function test_vector_of_variables_objective()
     return
 end
 
-function test_issue_105()
-    cost = [100.0, 120.0, 150.0, 110.0, 200.0, 170.0]
-    time = [8.0, 3.0, 4.0, 2.0, 5.0, 4.0]
-    capacity = [10.0, 8.0]
-    demand = [5.0, 8.0, 5.0]
-    m, n = 2, 3
-    model = MOI.instantiate(; with_bridge_type = Float64) do
-        return MOA.Optimizer(HiGHS.Optimizer)
-    end
-    MOI.set(model, MOA.Algorithm(), MOA.KirlikSayin())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, m * n)
-    MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
-    MOI.add_constraint.(model, x, MOI.Integer())
-    X = reshape(x, m, n)
-    for i in 1:m
-        f_i = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[i, :]), 0.0)
-        MOI.add_constraint(model, f_i, MOI.LessThan(capacity[i]))
-    end
-    for j in 1:n
-        f_j = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[:, j]), 0.0)
-        MOI.add_constraint(model, f_j, MOI.EqualTo(demand[j]))
-    end
-    f = MOI.Utilities.operate(
-        vcat,
-        Float64,
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(cost, x), 0.0),
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(time, x), 0.0),
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    @test MOI.get(model, MOI.ResultCount()) == 6
-    y = [MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:6]
-    @test sort!(y) == [
-        [2380.0, 81.0, 18.0],
-        [2440.0, 78.0, 18.0],
-        [2500.0, 75.0, 18.0],
-        [2560.0, 72.0, 18.0],
-        [2620.0, 69.0, 18.0],
-        [2680.0, 66.0, 18.0],
-    ]
-    return
-end
-
-end
+end  # TestKirlikSayin
 
 TestKirlikSayin.run_tests()
diff --git a/test/algorithms/Lexicographic.jl b/test/algorithms/Lexicographic.jl
index cf99f0a..14993e9 100644
--- a/test/algorithms/Lexicographic.jl
+++ b/test/algorithms/Lexicographic.jl
@@ -61,16 +61,16 @@ function test_knapsack_default()
     MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
     MOI.add_constraint(model, sum(1.0 * x[i] for i in 1:4), MOI.LessThan(2.0))
     MOI.optimize!(model)
-    results = Dict(
+    results = [
         [0, 1, 1] => [0, 1, 0, 1],
         [1, 0, 1] => [1, 0, 0, 1],
         [1, 1, 0] => [1, 1, 0, 0],
-    )
+    ]
     @test MOI.get(model, MOI.ResultCount()) == 3
     for i in 1:MOI.get(model, MOI.ResultCount())
         X = round.(Int, MOI.get(model, MOI.VariablePrimal(i), x))
         Y = round.(Int, MOI.get(model, MOI.ObjectiveValue(i)))
-        @test results[Y] == X
+        @test results[i] == (Y => X)
     end
     return
 end
@@ -197,6 +197,6 @@ function test_knapsack_time_limit()
     return
 end
 
-end
+end  # module TestLexicographic
 
 TestLexicographic.run_tests()
diff --git a/test/algorithms/TambyVanderpooten.jl b/test/algorithms/TambyVanderpooten.jl
index 67b0a5e..231fc4c 100644
--- a/test/algorithms/TambyVanderpooten.jl
+++ b/test/algorithms/TambyVanderpooten.jl
@@ -11,6 +11,8 @@ import HiGHS
 import MultiObjectiveAlgorithms as MOA
 import MultiObjectiveAlgorithms: MOI
 
+include(joinpath(dirname(@__DIR__), "problems.jl"))
+
 function run_tests()
     for name in names(@__MODULE__; all = true)
         if startswith("$name", "test_")
@@ -22,502 +24,11 @@ function run_tests()
     return
 end
 
-function test_knapsack_min_p3()
-    p = 3
-    n = 10
-    W = 2137.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-    ]
-    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
+function test_problems()
     model = MOA.Optimizer(HiGHS.Optimizer)
     MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
     MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
-            for i in 1:p for j in 1:n
-        ],
-        ones(p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 0 1 1 1 0 1 1 0 1
-        0 1 1 1 1 0 1 0 1 1
-        1 1 1 1 1 0 0 0 0 1
-        0 1 1 1 1 0 0 1 0 1
-        1 1 1 1 1 0 0 0 1 0
-        1 0 1 1 1 0 0 1 1 0
-        0 0 1 1 1 0 1 1 1 1
-    ]
-    Y_N = Float64[
-        -3394 -3817 -3408
-        -3042 -4627 -3189
-        -2997 -3539 -3509
-        -2854 -3570 -3714
-        -2706 -3857 -3304
-        -2518 -3866 -3191
-        -2854 -4636 -3076
-    ]
-    Y_N .+= 1
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)'
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)'
-    y_sol, x_sol = y_sol[sortperm(collect(eachrow(y_sol))), :],
-    x_sol[sortperm(collect(eachrow(y_sol))), :]
-    Y_N, X_E = Y_N[sortperm(collect(eachrow(Y_N))), :],
-    X_E[sortperm(collect(eachrow(Y_N))), :]
-    @test isapprox(x_sol, X_E; atol = 1e-6)
-    @test isapprox(y_sol, Y_N; atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_max_p3()
-    p = 3
-    n = 10
-    W = 2137.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-    ]
-    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
-            i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 0 1 1 1 0 1 1 0 1
-        0 1 1 1 1 0 1 0 1 1
-        1 1 1 1 1 0 0 0 0 1
-        0 1 1 1 1 0 0 1 0 1
-        1 1 1 1 1 0 0 0 1 0
-        1 0 1 1 1 0 0 1 1 0
-        0 0 1 1 1 0 1 1 1 1
-    ]
-    Y_N = Float64[
-        3394 3817 3408
-        3042 4627 3189
-        2997 3539 3509
-        2854 3570 3714
-        2706 3857 3304
-        2518 3866 3191
-        2854 4636 3076
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)'
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)'
-    y_sol, x_sol = y_sol[sortperm(collect(eachrow(y_sol))), :],
-    x_sol[sortperm(collect(eachrow(y_sol))), :]
-    Y_N, X_E = Y_N[sortperm(collect(eachrow(Y_N))), :],
-    X_E[sortperm(collect(eachrow(Y_N))), :]
-    @test isapprox(x_sol, X_E; atol = 1e-6)
-    @test isapprox(y_sol, Y_N; atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_min_p4()
-    p = 4
-    n = 10
-    W = 2653.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-        557 898 148 63 78 964 246 662 386 272
-    ]
-    w = Float64[979 448 355 955 426 229 9 695 322 889]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
-            for i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 1 1 0 1 1 1 0 0 0
-        0 1 1 0 1 1 1 1 1 0
-        0 1 1 1 1 1 1 0 0 0
-        0 1 1 0 1 0 1 0 1 1
-        0 1 1 1 1 0 1 0 1 0
-        0 0 1 0 1 1 1 0 1 1
-        0 1 1 1 0 1 1 0 1 0
-        0 0 1 1 0 1 1 1 1 0
-        0 0 1 1 1 1 1 0 1 0
-        0 1 0 0 0 1 1 1 1 1
-    ]
-    Y_N = Float64[
-        -3269 -2320 -3059 -2891
-        -3152 -3232 -3596 -3382
-        -2883 -3237 -2535 -2397
-        -2862 -3648 -3049 -2028
-        -2725 -4064 -2652 -1819
-        -2435 -3618 -2282 -2094
-        -2092 -3244 -2643 -2705
-        -1904 -3253 -2530 -2469
-        -2298 -4034 -1885 -1885
-        -2146 -1944 -2947 -3428
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)'
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)'
-    y_sol, x_sol = y_sol[sortperm(collect(eachrow(y_sol))), :],
-    x_sol[sortperm(collect(eachrow(y_sol))), :]
-    Y_N, X_E = Y_N[sortperm(collect(eachrow(Y_N))), :],
-    X_E[sortperm(collect(eachrow(Y_N))), :]
-    @test isapprox(x_sol, X_E; atol = 1e-6)
-    @test isapprox(y_sol, Y_N; atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(minimum(Y_N; dims = 1))
-    return
-end
-
-function test_knapsack_max_p4()
-    p = 4
-    n = 10
-    W = 2653.0
-    C = Float64[
-        566 611 506 180 817 184 585 423 26 317
-        62 84 977 979 874 54 269 93 881 563
-        664 982 962 140 224 215 12 869 332 537
-        557 898 148 63 78 964 246 662 386 272
-    ]
-    w = Float64[979 448 355 955 426 229 9 695 322 889]
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, n)
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    MOI.add_constraint(
-        model,
-        MOI.ScalarAffineFunction(
-            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
-            0.0,
-        ),
-        MOI.LessThan(W),
-    )
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
-            i in 1:p for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        1 1 1 0 1 1 1 0 0 0
-        0 1 1 0 1 1 1 1 1 0
-        0 1 1 1 1 1 1 0 0 0
-        0 1 1 0 1 0 1 0 1 1
-        0 1 1 1 1 0 1 0 1 0
-        0 0 1 0 1 1 1 0 1 1
-        0 1 1 1 0 1 1 0 1 0
-        0 0 1 1 0 1 1 1 1 0
-        0 0 1 1 1 1 1 0 1 0
-        0 1 0 0 0 1 1 1 1 1
-    ]
-    Y_N = Float64[
-        3269 2320 3059 2891
-        3152 3232 3596 3382
-        2883 3237 2535 2397
-        2862 3648 3049 2028
-        2725 4064 2652 1819
-        2435 3618 2282 2094
-        2092 3244 2643 2705
-        1904 3253 2530 2469
-        2298 4034 1885 1885
-        2146 1944 2947 3428
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol = hcat([MOI.get(model, MOI.VariablePrimal(i), x) for i in 1:N]...)'
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)'
-    y_sol, x_sol = y_sol[sortperm(collect(eachrow(y_sol))), :],
-    x_sol[sortperm(collect(eachrow(y_sol))), :]
-    Y_N, X_E = Y_N[sortperm(collect(eachrow(Y_N))), :],
-    X_E[sortperm(collect(eachrow(Y_N))), :]
-    @test isapprox(x_sol, X_E; atol = 1e-6)
-    @test isapprox(y_sol, Y_N; atol = 1e-6)
-    @test MOI.get(model, MOI.ObjectiveBound()) ≈ vec(maximum(Y_N; dims = 1))
-    return
-end
-
-function test_assignment_min_p3()
-    p = 3
-    n = 5
-    C = Float64[
-        6 1 20 2 3
-        2 6 9 10 18
-        1 6 20 5 9
-        6 8 6 9 6
-        7 10 10 6 2
-        17 20 8 8 20
-        10 13 1 10 15
-        4 11 1 13 1
-        19 13 7 18 17
-        15 3 5 1 11
-        10 7 1 19 12
-        2 15 12 10 3
-        11 20 16 12 9
-        10 15 20 11 7
-        1 9 20 7 6
-    ]
-    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
-    MOI.set(model, MOI.Silent(), true)
-    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    for i in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    for j in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
-        0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0
-        0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0
-    ]
-    Y_N = Float64[
-        16 61 47
-        17 43 71
-        18 47 67
-        20 52 54
-        22 37 63
-        23 43 44
-        22 54 47
-        28 66 39
-        34 60 42
-        24 39 45
-        35 49 39
-        37 55 36
-        28 33 58
-        35 38 56
-        39 43 41
-        38 33 53
-        45 33 34
-        43 51 31
-        40 47 37
-        29 29 59
-        50 40 32
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol =
-        hcat([MOI.get(model, MOI.VariablePrimal(i), vec(x)) for i in 1:N]...)'
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)'
-    y_sol, x_sol = y_sol[sortperm(collect(eachrow(y_sol))), :],
-    x_sol[sortperm(collect(eachrow(y_sol))), :]
-    Y_N, X_E = Y_N[sortperm(collect(eachrow(Y_N))), :],
-    X_E[sortperm(collect(eachrow(Y_N))), :]
-    @test isapprox(x_sol, X_E; atol = 1e-6)
-    @test isapprox(y_sol, Y_N; atol = 1e-6)
-    return
-end
-
-function test_assignment_max_p3()
-    p = 3
-    n = 5
-    C = Float64[
-        6 1 20 2 3
-        2 6 9 10 18
-        1 6 20 5 9
-        6 8 6 9 6
-        7 10 10 6 2
-        17 20 8 8 20
-        10 13 1 10 15
-        4 11 1 13 1
-        19 13 7 18 17
-        15 3 5 1 11
-        10 7 1 19 12
-        2 15 12 10 3
-        11 20 16 12 9
-        10 15 20 11 7
-        1 9 20 7 6
-    ]
-    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
-    model = MOA.Optimizer(HiGHS.Optimizer)
-    MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
-    MOI.set(model, MOI.Silent(), true)
-    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
-    MOI.add_constraint.(model, x, MOI.ZeroOne())
-    for i in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    for j in 1:n
-        MOI.add_constraint(
-            model,
-            MOI.ScalarAffineFunction(
-                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
-                0.0,
-            ),
-            MOI.EqualTo(1.0),
-        )
-    end
-    f = MOI.VectorAffineFunction(
-        [
-            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(-C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
-        ],
-        fill(0.0, p),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    X_E = Float64[
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1
-        0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1
-        0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
-        0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
-        0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0
-        0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0
-        0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0
-        0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0
-        0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0
-    ]
-    Y_N = Float64[
-        -16 -61 -47
-        -17 -43 -71
-        -18 -47 -67
-        -20 -52 -54
-        -22 -37 -63
-        -23 -43 -44
-        -22 -54 -47
-        -28 -66 -39
-        -34 -60 -42
-        -24 -39 -45
-        -35 -49 -39
-        -37 -55 -36
-        -28 -33 -58
-        -35 -38 -56
-        -39 -43 -41
-        -38 -33 -53
-        -45 -33 -34
-        -43 -51 -31
-        -40 -47 -37
-        -29 -29 -59
-        -50 -40 -32
-    ]
-    N = MOI.get(model, MOI.ResultCount())
-    x_sol =
-        hcat([MOI.get(model, MOI.VariablePrimal(i), vec(x)) for i in 1:N]...)'
-    y_sol = hcat([MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]...)'
-    y_sol, x_sol = y_sol[sortperm(collect(eachrow(y_sol))), :],
-    x_sol[sortperm(collect(eachrow(y_sol))), :]
-    Y_N, X_E = Y_N[sortperm(collect(eachrow(Y_N))), :],
-    X_E[sortperm(collect(eachrow(Y_N))), :]
-    @test isapprox(x_sol, X_E; atol = 1e-6)
-    @test isapprox(y_sol, Y_N; atol = 1e-6)
+    Problems.run_tests(model)
     return
 end
 
@@ -625,52 +136,6 @@ function test_vector_of_variables_objective()
     return
 end
 
-function test_issue_105()
-    cost = [100.0, 120.0, 150.0, 110.0, 200.0, 170.0]
-    time = [8.0, 3.0, 4.0, 2.0, 5.0, 4.0]
-    capacity = [10.0, 8.0]
-    demand = [5.0, 8.0, 5.0]
-    m, n = 2, 3
-    model = MOI.instantiate(; with_bridge_type = Float64) do
-        return MOA.Optimizer(HiGHS.Optimizer)
-    end
-    MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
-    MOI.set(model, MOI.Silent(), true)
-    x = MOI.add_variables(model, m * n)
-    MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
-    MOI.add_constraint.(model, x, MOI.Integer())
-    X = reshape(x, m, n)
-    for i in 1:m
-        f_i = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[i, :]), 0.0)
-        MOI.add_constraint(model, f_i, MOI.LessThan(capacity[i]))
-    end
-    for j in 1:n
-        f_j = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[:, j]), 0.0)
-        MOI.add_constraint(model, f_j, MOI.EqualTo(demand[j]))
-    end
-    f = MOI.Utilities.operate(
-        vcat,
-        Float64,
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(cost, x), 0.0),
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(time, x), 0.0),
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
-    )
-    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
-    MOI.optimize!(model)
-    @test MOI.get(model, MOI.ResultCount()) == 6
-    y = [MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:6]
-    @test sort!(y) == [
-        [2380.0, 81.0, 18.0],
-        [2440.0, 78.0, 18.0],
-        [2500.0, 75.0, 18.0],
-        [2560.0, 72.0, 18.0],
-        [2620.0, 69.0, 18.0],
-        [2680.0, 66.0, 18.0],
-    ]
-    return
-end
-
-end
+end  # module TestTambyVanderpooten
 
 TestTambyVanderpooten.run_tests()
diff --git a/test/problems.jl b/test/problems.jl
new file mode 100644
index 0000000..eb83ce6
--- /dev/null
+++ b/test/problems.jl
@@ -0,0 +1,444 @@
+#  Copyright 2019, Oscar Dowson and contributors
+#  This Source Code Form is subject to the terms of the Mozilla Public License,
+#  v.2.0. If a copy of the MPL was not distributed with this file, You can
+#  obtain one at http://mozilla.org/MPL/2.0/.
+
+module Problems
+
+using Test
+import MathOptInterface as MOI
+
+function run_tests(model::MOI.ModelLike)
+    for name in names(@__MODULE__; all = true)
+        if startswith("$name", "test_")
+            @testset "$name" begin
+                MOI.empty!(model)
+                getfield(@__MODULE__, name)(model)
+            end
+        end
+    end
+    return
+end
+
+function test_problem_knapsack_min_p3(model)
+    p = 3
+    n = 10
+    W = 2137.0
+    C = Float64[
+        566 611 506 180 817 184 585 423 26 317
+        62 84 977 979 874 54 269 93 881 563
+        664 982 962 140 224 215 12 869 332 537
+    ]
+    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
+    x = MOI.add_variables(model, n)
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
+            0.0,
+        ),
+        MOI.LessThan(W),
+    )
+    f = MOI.VectorAffineFunction(
+        [
+            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
+            for i in 1:p for j in 1:n
+        ],
+        ones(p),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    results = [
+        [1, 0, 1, 1, 1, 0, 1, 1, 0, 1] => [-3394, -3817, -3408],
+        [0, 1, 1, 1, 1, 0, 1, 0, 1, 1] => [-3042, -4627, -3189],
+        [1, 1, 1, 1, 1, 0, 0, 0, 0, 1] => [-2997, -3539, -3509],
+        [0, 0, 1, 1, 1, 0, 1, 1, 1, 1] => [-2854, -4636, -3076],
+        [0, 1, 1, 1, 1, 0, 0, 1, 0, 1] => [-2854, -3570, -3714],
+        [1, 1, 1, 1, 1, 0, 0, 0, 1, 0] => [-2706, -3857, -3304],
+        [1, 0, 1, 1, 1, 0, 0, 1, 1, 0] => [-2518, -3866, -3191],
+    ]
+    @assert MOI.get(model, MOI.ResultCount()) == length(results)
+    for i in 1:length(results)
+        X = MOI.get(model, MOI.VariablePrimal(i), x)
+        Y = MOI.get(model, MOI.ObjectiveValue(i))
+        @test isapprox(results[i][1], X; atol = 1e-6)
+        @test isapprox(results[i][2] .+ 1, Y; atol = 1e-6)
+    end
+    @test ≈(
+        MOI.get(model, MOI.ObjectiveBound()),
+        vec(minimum(mapreduce(last, hcat, results); dims = 2)) .+ 1,
+    )
+    return
+end
+
+function test_problem_knapsack_max_p3(model)
+    p = 3
+    n = 10
+    W = 2137.0
+    C = Float64[
+        566 611 506 180 817 184 585 423 26 317
+        62 84 977 979 874 54 269 93 881 563
+        664 982 962 140 224 215 12 869 332 537
+    ]
+    w = Float64[557, 898, 148, 63, 78, 964, 246, 662, 386, 272]
+    x = MOI.add_variables(model, n)
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
+            0.0,
+        ),
+        MOI.LessThan(W),
+    )
+    f = MOI.VectorAffineFunction(
+        [
+            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
+            i in 1:p for j in 1:n
+        ],
+        fill(0.0, p),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    results = [
+        [1, 0, 1, 1, 1, 0, 0, 1, 1, 0] => [2518, 3866, 3191],
+        [1, 1, 1, 1, 1, 0, 0, 0, 1, 0] => [2706, 3857, 3304],
+        [0, 1, 1, 1, 1, 0, 0, 1, 0, 1] => [2854, 3570, 3714],
+        [0, 0, 1, 1, 1, 0, 1, 1, 1, 1] => [2854, 4636, 3076],
+        [1, 1, 1, 1, 1, 0, 0, 0, 0, 1] => [2997, 3539, 3509],
+        [0, 1, 1, 1, 1, 0, 1, 0, 1, 1] => [3042, 4627, 3189],
+        [1, 0, 1, 1, 1, 0, 1, 1, 0, 1] => [3394, 3817, 3408],
+    ]
+    N = MOI.get(model, MOI.ResultCount())
+    @assert N == length(results)
+    for i in 1:length(results)
+        X = MOI.get(model, MOI.VariablePrimal(i), x)
+        Y = MOI.get(model, MOI.ObjectiveValue(i))
+        @test isapprox(results[i][1], X; atol = 1e-6)
+        @test isapprox(results[i][2], Y; atol = 1e-6)
+    end
+    @test ≈(
+        MOI.get(model, MOI.ObjectiveBound()),
+        vec(maximum(mapreduce(last, hcat, results); dims = 2)),
+    )
+    return
+end
+
+function test_problem_knapsack_min_p4(model)
+    p = 4
+    n = 10
+    W = 2653.0
+    C = Float64[
+        566 611 506 180 817 184 585 423 26 317
+        62 84 977 979 874 54 269 93 881 563
+        664 982 962 140 224 215 12 869 332 537
+        557 898 148 63 78 964 246 662 386 272
+    ]
+    w = Float64[979 448 355 955 426 229 9 695 322 889]
+    x = MOI.add_variables(model, n)
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
+            0.0,
+        ),
+        MOI.LessThan(W),
+    )
+    f = MOI.VectorAffineFunction(
+        [
+            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(-C[i, j], x[j]))
+            for i in 1:p for j in 1:n
+        ],
+        fill(0.0, p),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    results = reverse([
+        [0, 0, 1, 1, 0, 1, 1, 1, 1, 0] => [1904, 3253, 2530, 2469],
+        [0, 1, 1, 1, 0, 1, 1, 0, 1, 0] => [2092, 3244, 2643, 2705],
+        [0, 1, 0, 0, 0, 1, 1, 1, 1, 1] => [2146, 1944, 2947, 3428],
+        [0, 0, 1, 1, 1, 1, 1, 0, 1, 0] => [2298, 4034, 1885, 1885],
+        [0, 0, 1, 0, 1, 1, 1, 0, 1, 1] => [2435, 3618, 2282, 2094],
+        [0, 1, 1, 1, 1, 0, 1, 0, 1, 0] => [2725, 4064, 2652, 1819],
+        [0, 1, 1, 0, 1, 0, 1, 0, 1, 1] => [2862, 3648, 3049, 2028],
+        [0, 1, 1, 1, 1, 1, 1, 0, 0, 0] => [2883, 3237, 2535, 2397],
+        [0, 1, 1, 0, 1, 1, 1, 1, 1, 0] => [3152, 3232, 3596, 3382],
+        [1, 1, 1, 0, 1, 1, 1, 0, 0, 0] => [3269, 2320, 3059, 2891],
+    ])
+    @test MOI.get(model, MOI.ResultCount()) == length(results)
+    for (i, (x_sol, y_sol)) in enumerate(results)
+        @test ≈(x_sol, MOI.get(model, MOI.VariablePrimal(i), x); atol = 1e-6)
+        @test ≈(-y_sol, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ -[3269, 4064, 3596, 3428]
+    return
+end
+
+function test_problem_knapsack_max_p4(model)
+    p = 4
+    n = 10
+    W = 2653.0
+    C = Float64[
+        566 611 506 180 817 184 585 423 26 317
+        62 84 977 979 874 54 269 93 881 563
+        664 982 962 140 224 215 12 869 332 537
+        557 898 148 63 78 964 246 662 386 272
+    ]
+    w = Float64[979 448 355 955 426 229 9 695 322 889]
+    x = MOI.add_variables(model, n)
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    MOI.add_constraint(
+        model,
+        MOI.ScalarAffineFunction(
+            [MOI.ScalarAffineTerm(w[j], x[j]) for j in 1:n],
+            0.0,
+        ),
+        MOI.LessThan(W),
+    )
+    f = MOI.VectorAffineFunction(
+        [
+            MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(C[i, j], x[j])) for
+            i in 1:p for j in 1:n
+        ],
+        fill(0.0, p),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    results = [
+        [0, 0, 1, 1, 0, 1, 1, 1, 1, 0] => [1904, 3253, 2530, 2469],
+        [0, 1, 1, 1, 0, 1, 1, 0, 1, 0] => [2092, 3244, 2643, 2705],
+        [0, 1, 0, 0, 0, 1, 1, 1, 1, 1] => [2146, 1944, 2947, 3428],
+        [0, 0, 1, 1, 1, 1, 1, 0, 1, 0] => [2298, 4034, 1885, 1885],
+        [0, 0, 1, 0, 1, 1, 1, 0, 1, 1] => [2435, 3618, 2282, 2094],
+        [0, 1, 1, 1, 1, 0, 1, 0, 1, 0] => [2725, 4064, 2652, 1819],
+        [0, 1, 1, 0, 1, 0, 1, 0, 1, 1] => [2862, 3648, 3049, 2028],
+        [0, 1, 1, 1, 1, 1, 1, 0, 0, 0] => [2883, 3237, 2535, 2397],
+        [0, 1, 1, 0, 1, 1, 1, 1, 1, 0] => [3152, 3232, 3596, 3382],
+        [1, 1, 1, 0, 1, 1, 1, 0, 0, 0] => [3269, 2320, 3059, 2891],
+    ]
+    @test MOI.get(model, MOI.ResultCount()) == length(results)
+    for (i, (x_sol, y_sol)) in enumerate(results)
+        @test ≈(x_sol, MOI.get(model, MOI.VariablePrimal(i), x); atol = 1e-6)
+        @test ≈(y_sol, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    @test MOI.get(model, MOI.ObjectiveBound()) ≈ [3269, 4064, 3596, 3428]
+    return
+end
+
+function test_problem_assignment_min_p3(model)
+    p = 3
+    n = 5
+    C = Float64[
+        6 1 20 2 3
+        2 6 9 10 18
+        1 6 20 5 9
+        6 8 6 9 6
+        7 10 10 6 2
+        17 20 8 8 20
+        10 13 1 10 15
+        4 11 1 13 1
+        19 13 7 18 17
+        15 3 5 1 11
+        10 7 1 19 12
+        2 15 12 10 3
+        11 20 16 12 9
+        10 15 20 11 7
+        1 9 20 7 6
+    ]
+    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
+    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    for i in 1:n
+        MOI.add_constraint(
+            model,
+            MOI.ScalarAffineFunction(
+                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
+                0.0,
+            ),
+            MOI.EqualTo(1.0),
+        )
+    end
+    for j in 1:n
+        MOI.add_constraint(
+            model,
+            MOI.ScalarAffineFunction(
+                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
+                0.0,
+            ),
+            MOI.EqualTo(1.0),
+        )
+    end
+    f = MOI.VectorAffineFunction(
+        [
+            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
+        ],
+        fill(0.0, p),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    results = [
+        [0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1] => [16, 61, 47],
+        [0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1] => [17, 43, 71],
+        [0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1] => [18, 47, 67],
+        [0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1] => [20, 52, 54],
+        [0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1] => [22, 37, 63],
+        [0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1] => [22, 54, 47],
+        [0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0] => [23, 43, 44],
+        [0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0] => [24, 39, 45],
+        [0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0] => [28, 33, 58],
+        [0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0] => [28, 66, 39],
+        [0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0] => [29, 29, 59],
+        [0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1] => [34, 60, 42],
+        [0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0] => [35, 38, 56],
+        [0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0] => [35, 49, 39],
+        [0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1] => [37, 55, 36],
+        [1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0] => [38, 33, 53],
+        [0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0] => [39, 43, 41],
+        [0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0] => [40, 47, 37],
+        [0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0] => [43, 51, 31],
+        [0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0] => [45, 33, 34],
+        [0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0] => [50, 40, 32],
+    ]
+    @test MOI.get(model, MOI.ResultCount()) == length(results)
+    for (i, (x_sol, y_sol)) in enumerate(results)
+        x_primal = MOI.get(model, MOI.VariablePrimal(i), vec(x))
+        @test ≈(vec(x_sol), x_primal; atol = 1e-6)
+        @test ≈(y_sol, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    return
+end
+
+function test_problem_assignment_max_p3(model)
+    p = 3
+    n = 5
+    C = Float64[
+        6 1 20 2 3
+        2 6 9 10 18
+        1 6 20 5 9
+        6 8 6 9 6
+        7 10 10 6 2
+        17 20 8 8 20
+        10 13 1 10 15
+        4 11 1 13 1
+        19 13 7 18 17
+        15 3 5 1 11
+        10 7 1 19 12
+        2 15 12 10 3
+        11 20 16 12 9
+        10 15 20 11 7
+        1 9 20 7 6
+    ]
+    C = permutedims(reshape(C, (n, p, n)), [2, 1, 3])
+    x = [MOI.add_variable(model) for i in 1:n, j in 1:n]
+    MOI.add_constraint.(model, x, MOI.ZeroOne())
+    for i in 1:n
+        MOI.add_constraint(
+            model,
+            MOI.ScalarAffineFunction(
+                [MOI.ScalarAffineTerm(1.0, x[i, j]) for j in 1:n],
+                0.0,
+            ),
+            MOI.EqualTo(1.0),
+        )
+    end
+    for j in 1:n
+        MOI.add_constraint(
+            model,
+            MOI.ScalarAffineFunction(
+                [MOI.ScalarAffineTerm(1.0, x[i, j]) for i in 1:n],
+                0.0,
+            ),
+            MOI.EqualTo(1.0),
+        )
+    end
+    f = MOI.VectorAffineFunction(
+        [
+            MOI.VectorAffineTerm(k, MOI.ScalarAffineTerm(-C[k, i, j], x[i, j])) for k in 1:p for i in 1:n for j in 1:n
+        ],
+        fill(0.0, p),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    results = reverse([
+        [0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1] => [16, 61, 47],
+        [0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1] => [17, 43, 71],
+        [0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1] => [18, 47, 67],
+        [0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1] => [20, 52, 54],
+        [0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1] => [22, 37, 63],
+        [0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1] => [22, 54, 47],
+        [0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0] => [23, 43, 44],
+        [0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0] => [24, 39, 45],
+        [0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0] => [28, 33, 58],
+        [0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0] => [28, 66, 39],
+        [0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0] => [29, 29, 59],
+        [0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1] => [34, 60, 42],
+        [0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0] => [35, 38, 56],
+        [0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0] => [35, 49, 39],
+        [0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1] => [37, 55, 36],
+        [1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0] => [38, 33, 53],
+        [0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0] => [39, 43, 41],
+        [0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0] => [40, 47, 37],
+        [0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0] => [43, 51, 31],
+        [0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0] => [45, 33, 34],
+        [0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0] => [50, 40, 32],
+    ])
+    @test MOI.get(model, MOI.ResultCount()) == length(results)
+    for (i, (x_sol, y_sol)) in enumerate(results)
+        x_primal = MOI.get(model, MOI.VariablePrimal(i), vec(x))
+        @test ≈(vec(x_sol), x_primal; atol = 1e-6)
+        @test ≈(-y_sol, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    return
+end
+
+function test_problem_issue_105(model)
+    cost = [100.0, 120.0, 150.0, 110.0, 200.0, 170.0]
+    time = [8.0, 3.0, 4.0, 2.0, 5.0, 4.0]
+    capacity = [10.0, 8.0]
+    demand = [5.0, 8.0, 5.0]
+    m, n = 2, 3
+    x = MOI.add_variables(model, m * n)
+    MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
+    MOI.add_constraint.(model, x, MOI.Integer())
+    X = reshape(x, m, n)
+    for i in 1:m
+        f_i = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[i, :]), 0.0)
+        MOI.add_constraint(model, f_i, MOI.LessThan(capacity[i]))
+    end
+    for j in 1:n
+        f_j = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, X[:, j]), 0.0)
+        MOI.add_constraint(model, f_j, MOI.EqualTo(demand[j]))
+    end
+    f = MOI.Utilities.operate(
+        vcat,
+        Float64,
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(cost, x), 0.0),
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(time, x), 0.0),
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
+    )
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+    MOI.optimize!(model)
+    @test MOI.get(model, MOI.ResultCount()) == 6
+    for (i, y) in enumerate([
+        [2380.0, 81.0, 18.0],
+        [2440.0, 78.0, 18.0],
+        [2500.0, 75.0, 18.0],
+        [2560.0, 72.0, 18.0],
+        [2620.0, 69.0, 18.0],
+        [2680.0, 66.0, 18.0],
+    ])
+        @test ≈(y, MOI.get(model, MOI.ObjectiveValue(i)); atol = 1e-6)
+    end
+    return
+end
+
+end  # module Problems
diff --git a/test/test_utilities.jl b/test/test_utilities.jl
index fc86fbc..b65be5f 100644
--- a/test/test_utilities.jl
+++ b/test/test_utilities.jl
@@ -95,11 +95,11 @@ function test_filter_nondominated_triple()
             MOA.SolutionPoint(x, [1, 1, 1][p]),
         ]
         # The permutation can change the ordering of the solutions that are
-        # returnned, so we can't use `@test min_sol == solutions[1:2]`
-        min_sol = MOA.filter_nondominated(MOI.MIN_SENSE, solutions)
+        # returned, so we can't use `@test min_sol == solutions[1:2]`
+        min_sol = MOA.filter_nondominated(MOI.MIN_SENSE, copy(solutions))
         @test solutions[1] in min_sol && solutions[2] in min_sol
         @test length(min_sol) == 2
-        max_sol = MOA.filter_nondominated(MOI.MAX_SENSE, solutions)
+        max_sol = MOA.filter_nondominated(MOI.MAX_SENSE, copy(solutions))
         @test solutions[2] in max_sol && solutions[3] in max_sol
         @test length(max_sol) == 2
     end