diff --git a/docs/src/background/algebraic_modeling_languages.md b/docs/src/background/algebraic_modeling_languages.md
index 2256b849238..6ff34d1d948 100644
--- a/docs/src/background/algebraic_modeling_languages.md
+++ b/docs/src/background/algebraic_modeling_languages.md
@@ -138,6 +138,9 @@ julia> function algebraic_knapsack(c, w, b)
            @objective(model, Max, sum(c[i] * x[i] for i = 1:n))
            @constraint(model, sum(w[i] * x[i] for i = 1:n) <= b)
            optimize!(model)
+           if termination_status(model) != OPTIMAL
+               error("Not solved correctly")
+           end
            return value.(x)
        end
 algebraic_knapsack (generic function with 1 method)
@@ -179,6 +182,9 @@ julia> function nonalgebraic_knapsack(c, w, b)
            con = build_constraint(error, lhs, MOI.LessThan(b))
            add_constraint(model, con)
            optimize!(model)
+           if termination_status(model) != OPTIMAL
+               error("Not solved correctly")
+           end
            return value.(x)
        end
 nonalgebraic_knapsack (generic function with 1 method)
@@ -219,6 +225,9 @@ julia> function mathoptinterface_knapsack(optimizer, c, w, b)
                MOI.LessThan(b),
            )
            MOI.optimize!(model)
+           if MOI.get(model, MOI.TerminationStatus()) != MOI.OPTIMAL
+               error("Not solved correctly")
+           end
            return MOI.get.(model, MOI.VariablePrimal(), x)
        end
 mathoptinterface_knapsack (generic function with 1 method)
@@ -257,6 +266,9 @@ julia> function highs_knapsack(c, w, b)
                w,
            )
            Highs_run(model)
+           if Highs_getModelStatus(model) != kHighsModelStatusOptimal
+               error("Not solved correctly")
+           end
            x = fill(NaN, 2)
            Highs_getSolution(model, x, C_NULL, C_NULL, C_NULL)
            Highs_destroy(model)
diff --git a/docs/src/manual/solutions.md b/docs/src/manual/solutions.md
index 823f9c22734..99d83b4d164 100644
--- a/docs/src/manual/solutions.md
+++ b/docs/src/manual/solutions.md
@@ -32,6 +32,34 @@ Subject to
  y[b] ≤ 1
 ```
 
+## Check if an optimal solution exists
+
+Use [`is_solved_and_feasible`](@ref) to check if the solver found an optimal
+solution:
+```jldoctest solutions
+julia> is_solved_and_feasible(model)
+true
+```
+
+By default, [`is_solved_and_feasible`](@ref) returns `true` for both global and
+local optima. Pass `allow_local = false` to check if the solver found a globally
+optimal solution:
+```jldoctest solutions
+julia> is_solved_and_feasible(model; allow_local = false)
+true
+```
+
+Pass `dual = true` to check if the solver found an optimal dual solution in
+addition to an optimal primal solution:
+```jldoctest solutions
+julia> is_solved_and_feasible(model; dual = true)
+true
+```
+
+If this function returns `false`, use the functions mentioned below like
+[`solution_summary`](@ref), [`termination_status`](@ref), [`primal_status`](@ref),
+and [`dual_status`](@ref) to understand what solution (if any) the solver found.
+
 ## Solutions summary
 
 [`solution_summary`](@ref) can be used for checking the summary of the
@@ -267,25 +295,101 @@ And data, a 2-element Vector{Float64}:
 
 ## Recommended workflow
 
-The recommended workflow for solving a model and querying the solution is
-something like the following:
+You should always check whether the solver found a solution before calling
+solution functions like [`value`](@ref) or [`objective_value`](@ref).
+
+A simple approach for small scripts and notebooks is to use
+[`is_solved_and_feasible`](@ref):
+
 ```jldoctest solutions
-julia> begin
-           if termination_status(model) == OPTIMAL
+julia> function solve_and_print_solution(model)
+           optimize!(model)
+           if !is_solved_and_feasible(model; dual = true)
+               error(
+                   """
+                   The model was not solved correctly:
+                   termination_status : $(termination_status(model))
+                   primal_status      : $(primal_status(model))
+                   dual_status        : $(dual_status(model))
+                   raw_status         : $(raw_status(model))
+                   """,
+               )
+           end
+           println("Solution is optimal")
+           println("  objective value = ", objective_value(model))
+           println("  primal solution: x = ", value(x))
+           println("  dual solution: c1 = ", dual(c1))
+           return
+       end
+solve_and_print_solution (generic function with 1 method)
+
+julia> solve_and_print_solution(model)
+Solution is optimal
+  objective value = -205.14285714285714
+  primal solution: x = 15.428571428571429
+  dual solution: c1 = 1.7142857142857142
+```
+
+For code like libraries that should be more robust to the range of possible
+termination and result statuses, do some variation of the following:
+```jldoctest solutions
+julia> function solve_and_print_solution(model)
+           status = termination_status(model)
+           if status in (OPTIMAL, LOCALLY_SOLVED)
                println("Solution is optimal")
-           elseif termination_status(model) == TIME_LIMIT && has_values(model)
-               println("Solution is suboptimal due to a time limit, but a primal solution is available")
+           elseif status in (ALMOST_OPTIMAL, ALMOST_LOCALLY_SOLVED)
+               println("Solution is optimal to a relaxed tolerance")
+           elseif status == TIME_LIMIT
+               println(
+                   "Solver stopped due to a time limit. If a solution is available, " *
+                   "it may be suboptimal."
+               )
+           elseif status in (
+               ITERATION_LIMIT, NODE_LIMIT, SOLUTION_LIMIT, MEMORY_LIMIT,
+               OBJECTIVE_LIMIT, NORM_LIMIT, OTHER_LIMIT,
+           )
+               println(
+                   "Solver stopped due to a limit. If a solution is available, it " *
+                   "may be suboptimal."
+               )
+           elseif status in (INFEASIBLE, LOCALLY_INFEASIBLE)
+               println("The problem is primal infeasible")
+           elseif status == DUAL_INFEASIBLE
+               println(
+                   "The problem is dual infeasible. If a primal feasible solution " *
+                   "exists, the problem is unbounded. To check, set the objective " *
+                   "to `@objective(model, Min, 0)` and re-solve. If the problem is " *
+                   "feasible, the primal is unbounded. If the problem is " *
+                   "infeasible, both the primal and dual are infeasible.",
+               )
+           elseif status == INFEASIBLE_OR_UNBOUNDED
+               println(
+                   "The model is either infeasible or unbounded. Set the objective " *
+                   "to `@objective(model, Min, 0)` and re-solve to disambiguate. If " *
+                   "the problem was infeasible, it will still be infeasible. If the " *
+                   "problem was unbounded, it will now have a finite optimal solution.",
+               )
            else
-               error("The model was not solved correctly.")
+               println(
+                   "The model was not solved correctly. The termination status is $status",
+               )
            end
-           println("  objective value = ", objective_value(model))
-           if primal_status(model) == FEASIBLE_POINT
+           if primal_status(model) in (FEASIBLE_POINT, NEARLY_FEASIBLE_POINT)
+               println("  objective value = ", objective_value(model))
                println("  primal solution: x = ", value(x))
+           elseif primal_status(model) == INFEASIBILITY_CERTIFICATE
+               println("  primal certificate: x = ", value(x))
            end
-           if dual_status(model) == FEASIBLE_POINT
+           if dual_status(model) in (FEASIBLE_POINT, NEARLY_FEASIBLE_POINT)
                println("  dual solution: c1 = ", dual(c1))
+           elseif dual_status(model) == INFEASIBILITY_CERTIFICATE
+               println("  dual certificate: c1 = ", dual(c1))
            end
+           return
        end
+solve_and_print_solution (generic function with 1 method)
+
+julia> solve_and_print_solution(model)
 Solution is optimal
   objective value = -205.14285714285714
   primal solution: x = 15.428571428571429
diff --git a/docs/src/tutorials/algorithms/benders_decomposition.jl b/docs/src/tutorials/algorithms/benders_decomposition.jl
index 7ca882bedee..fcab4850e91 100644
--- a/docs/src/tutorials/algorithms/benders_decomposition.jl
+++ b/docs/src/tutorials/algorithms/benders_decomposition.jl
@@ -164,7 +164,7 @@ function solve_subproblem(x)
     con = @constraint(model, A_2 * y .<= b - A_1 * x)
     @objective(model, Min, c_2' * y)
     optimize!(model)
-    @assert termination_status(model) == OPTIMAL
+    @assert is_solved_and_feasible(model; dual = true)
     return (obj = objective_value(model), y = value.(y), π = dual.(con))
 end
 
@@ -194,6 +194,7 @@ ABSOLUTE_OPTIMALITY_GAP = 1e-6
 println("Iteration  Lower Bound  Upper Bound          Gap")
 for k in 1:MAXIMUM_ITERATIONS
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     lower_bound = objective_value(model)
     x_k = value.(x)
     ret = solve_subproblem(x_k)
@@ -211,6 +212,7 @@ end
 # Finally, we can obtain the optimal solution
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test value.(x) == [0.0, 1.0]  #src
 x_optimal = value.(x)
 
@@ -268,6 +270,7 @@ set_attribute(lazy_model, MOI.LazyConstraintCallback(), my_callback)
 # Now when we optimize!, our callback is run:
 
 optimize!(lazy_model)
+@assert is_solved_and_feasible(lazy_model)
 
 # For this model, the callback algorithm required more solves of the subproblem:
 
@@ -323,7 +326,7 @@ print(subproblem)
 function solve_subproblem(model, x)
     fix.(model[:x_copy], x)
     optimize!(model)
-    @assert termination_status(model) == OPTIMAL
+    @assert is_solved_and_feasible(model; dual = true)
     return (
         obj = objective_value(model),
         y = value.(model[:y]),
@@ -341,6 +344,7 @@ end
 println("Iteration  Lower Bound  Upper Bound          Gap")
 for k in 1:MAXIMUM_ITERATIONS
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     lower_bound = objective_value(model)
     x_k = value.(x)
     ret = solve_subproblem(subproblem, x_k)
@@ -358,6 +362,7 @@ end
 # Finally, we can obtain the optimal solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test value.(x) == [0.0, 1.0]  #src
 x_optimal = value.(x)
 
diff --git a/docs/src/tutorials/algorithms/cutting_stock_column_generation.jl b/docs/src/tutorials/algorithms/cutting_stock_column_generation.jl
index 93e2623e2d1..9e7809772b4 100644
--- a/docs/src/tutorials/algorithms/cutting_stock_column_generation.jl
+++ b/docs/src/tutorials/algorithms/cutting_stock_column_generation.jl
@@ -235,6 +235,7 @@ set_silent(model)
 @objective(model, Min, sum(x))
 @constraint(model, demand[i in 1:I], patterns[i]' * x >= data.pieces[i].d)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # This solution requires 421 rolls. This solution is sub-optimal because the
@@ -252,6 +253,7 @@ solution_summary(model)
 
 unset_integer.(x)
 optimize!(model)
+@assert is_solved_and_feasible(model; dual = true)
 π_13 = dual(demand[13])
 
 # Using the economic interpretation of the dual variable, we can say that a one
@@ -282,6 +284,7 @@ function solve_pricing(data::Data, π::Vector{Float64})
     @constraint(model, sum(data.pieces[i].w * y[i] for i in 1:I) <= data.W)
     @objective(model, Max, sum(π[i] * y[i] for i in 1:I))
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     number_of_rolls_saved = objective_value(model)
     if number_of_rolls_saved > 1 + 1e-8
         ## Benefit of pattern is more than the cost of a new roll plus some
@@ -312,6 +315,7 @@ solve_pricing(data, zeros(I))
 while true
     ## Solve the linear relaxation
     optimize!(model)
+    @assert is_solved_and_feasible(model; dual = true)
     ## Obtain a new dual vector
     π = dual.(demand)
     ## Solve the pricing problem
@@ -362,6 +366,7 @@ sum(ceil.(Int, solution.rolls))
 
 set_integer.(x)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution = DataFrames.DataFrame([
     (pattern = p, rolls = value(x_p)) for (p, x_p) in enumerate(x)
 ])
diff --git a/docs/src/tutorials/algorithms/parallelism.md b/docs/src/tutorials/algorithms/parallelism.md
index 2af638522cf..978760b5bd0 100644
--- a/docs/src/tutorials/algorithms/parallelism.md
+++ b/docs/src/tutorials/algorithms/parallelism.md
@@ -215,6 +215,7 @@ my_lock = Threads.ReentrantLock()
 Threads.@threads for i in 1:10
     set_lower_bound(x, i)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     Threads.lock(my_lock) do
         push!(solutions, i => objective_value(model))
     end
@@ -251,6 +252,7 @@ julia> Threads.@threads for i in 1:10
            @objective(model, Min, x)
            set_lower_bound(x, i)
            optimize!(model)
+           @assert is_solved_and_feasible(sudoku)
            Threads.lock(my_lock) do
                push!(solutions, i => objective_value(model))
            end
@@ -295,6 +297,7 @@ julia> Distributed.@everywhere begin
                @objective(model, Min, x)
                set_lower_bound(x, i)
                optimize!(model)
+               @assert is_solved_and_feasible(sudoku)
                return objective_value(model)
            end
        end
diff --git a/docs/src/tutorials/algorithms/tsp_lazy_constraints.jl b/docs/src/tutorials/algorithms/tsp_lazy_constraints.jl
index 62eb3d94978..51add8abeb0 100644
--- a/docs/src/tutorials/algorithms/tsp_lazy_constraints.jl
+++ b/docs/src/tutorials/algorithms/tsp_lazy_constraints.jl
@@ -199,6 +199,7 @@ subtour(x::AbstractMatrix{VariableRef}) = subtour(value.(x))
 
 iterative_model = build_tsp_model(d, n)
 optimize!(iterative_model)
+@assert is_solved_and_feasible(iterative_model)
 time_iterated = solve_time(iterative_model)
 cycle = subtour(iterative_model[:x])
 while 1 < length(cycle) < n
@@ -209,6 +210,7 @@ while 1 < length(cycle) < n
         sum(iterative_model[:x][i, j] for (i, j) in S) <= length(cycle) - 1,
     )
     optimize!(iterative_model)
+    @assert is_solved_and_feasible(iterative_model)
     global time_iterated += solve_time(iterative_model)
     global cycle = subtour(iterative_model[:x])
 end
@@ -262,6 +264,7 @@ set_attribute(
     subtour_elimination_callback,
 )
 optimize!(lazy_model)
+@assert is_solved_and_feasible(lazy_model)
 objective_value(lazy_model)
 
 # This finds the same optimal tour:
diff --git a/docs/src/tutorials/applications/optimal_power_flow.jl b/docs/src/tutorials/applications/optimal_power_flow.jl
index d23b09034ee..a54bbafa338 100644
--- a/docs/src/tutorials/applications/optimal_power_flow.jl
+++ b/docs/src/tutorials/applications/optimal_power_flow.jl
@@ -41,7 +41,7 @@ import DataFrames
 import Ipopt
 import LinearAlgebra
 import SparseArrays
-import Test #src
+import Test
 
 # ## Initial formulation
 
@@ -137,6 +137,7 @@ println("Objective value (basic lower bound) : $basic_lower_bound")
 
 @constraint(model, sum(P_G) >= sum(P_Demand))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 better_lower_bound = round(objective_value(model); digits = 2)
 println("Objective value (better lower bound): $better_lower_bound")
 
@@ -280,6 +281,7 @@ P_G = real(S_G)
 # We're finally ready to solve our nonlinear AC-OPF problem:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test isapprox(objective_value(model), 3087.84; atol = 1e-2)  #src
 solution_summary(model)
 
@@ -418,9 +420,8 @@ optimize!(model)
 
 #-
 
+Test.@test is_solved_and_feasible(model; allow_almost = true)
 sdp_relaxation_lower_bound = round(objective_value(model); digits = 2)
-Test.@test termination_status(model) in (OPTIMAL, ALMOST_OPTIMAL)           #src
-Test.@test primal_status(model) in (FEASIBLE_POINT, NEARLY_FEASIBLE_POINT)  #src
 Test.@test isapprox(sdp_relaxation_lower_bound, 2753.04; rtol = 1e-3)     #src
 println(
     "Objective value (W & V relax. lower bound): $sdp_relaxation_lower_bound",
diff --git a/docs/src/tutorials/applications/power_systems.jl b/docs/src/tutorials/applications/power_systems.jl
index b5e9a5d289a..9cd58397271 100644
--- a/docs/src/tutorials/applications/power_systems.jl
+++ b/docs/src/tutorials/applications/power_systems.jl
@@ -115,6 +115,7 @@ function solve_economic_dispatch(generators::Vector, wind, scenario)
     @constraint(model, sum(g[i] for i in 1:N) + w == scenario.demand)
     ## Solve statement
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     ## return the optimal value of the objective function and its minimizers
     return (
         g = value.(g),
@@ -216,6 +217,7 @@ function solve_economic_dispatch_inplace(
             wind.variable_cost * w,
         )
         optimize!(model)
+        @assert is_solved_and_feasible(model)
         push!(obj_out, objective_value(model))
         push!(w_out, value(w))
         push!(g1_out, value(g[1]))
@@ -382,6 +384,7 @@ function solve_unit_commitment(generators::Vector, wind, scenario)
     if status != OPTIMAL
         return (status = status,)
     end
+    @assert primal_status(model) == FEASIBLE_POINT
     return (
         status = status,
         g = value.(g),
@@ -525,6 +528,7 @@ function solve_nonlinear_economic_dispatch(
     )
     @constraint(model, sum(g[i] for i in 1:N) + sqrt(w) == scenario.demand)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return (
         g = value.(g),
         w = value(w),
diff --git a/docs/src/tutorials/applications/two_stage_stochastic.jl b/docs/src/tutorials/applications/two_stage_stochastic.jl
index 8b87bd745fb..ac7a8bdce70 100644
--- a/docs/src/tutorials/applications/two_stage_stochastic.jl
+++ b/docs/src/tutorials/applications/two_stage_stochastic.jl
@@ -86,6 +86,7 @@ set_silent(model)
 @expression(model, z[ω in Ω], 5y[ω] - 0.1 * (x - y[ω]))
 @objective(model, Max, -2x + sum(P[ω] * z[ω] for ω in Ω))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # The optimal number of pies to make is:
@@ -158,6 +159,7 @@ function CVaR(Z::Vector{Float64}, P::Vector{Float64}; γ::Float64)
     @constraint(model, [i in 1:N], z[i] >= ξ - Z[i])
     @objective(model, Max, ξ - 1 / γ * sum(P[i] * z[i] for i in 1:N))
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return objective_value(model)
 end
 
@@ -216,6 +218,7 @@ set_silent(model)
 @constraint(model, [ω in Ω], z[ω] >= ξ - Z[ω])
 @objective(model, Max, -2x + ξ - 1 / γ * sum(P[ω] * z[ω] for ω in Ω))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # When ``\gamma = 0.4``, the optimal number of pies to bake is:
 
diff --git a/docs/src/tutorials/applications/web_app.jl b/docs/src/tutorials/applications/web_app.jl
index 7481899871e..4d90fe619e1 100644
--- a/docs/src/tutorials/applications/web_app.jl
+++ b/docs/src/tutorials/applications/web_app.jl
@@ -57,7 +57,7 @@ function endpoint_solve(params::Dict{String,Any})
         "primal_status" => primal_status(model),
     )
     ## Only include the `x` key if it has a value.
-    if has_values(model)
+    if primal_status(model) == FEASIBLE_POINT
         ret["x"] = value(x)
     end
     return ret
diff --git a/docs/src/tutorials/conic/arbitrary_precision.jl b/docs/src/tutorials/conic/arbitrary_precision.jl
index c06a3def722..2323d25bbb5 100644
--- a/docs/src/tutorials/conic/arbitrary_precision.jl
+++ b/docs/src/tutorials/conic/arbitrary_precision.jl
@@ -76,6 +76,7 @@ print(model)
 # Let's solve and inspect the solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model; dual = true)
 solution_summary(model)
 
 # The value of each decision variable is a `BigFloat`:
@@ -100,6 +101,7 @@ value.(x) .- [3 // 7, 3 // 14]
 set_attribute(model, "tol_gap_abs", 1e-32)
 set_attribute(model, "tol_gap_rel", 1e-32)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value.(x) .- [3 // 7, 3 // 14]
 
 # ## Rational arithmetic
@@ -141,6 +143,7 @@ print(model)
 # Let's solve and inspect the solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # The optimal values are given in exact rational arithmetic:
diff --git a/docs/src/tutorials/conic/dualization.jl b/docs/src/tutorials/conic/dualization.jl
index 149c2d3f5f2..8a2fc18fde2 100644
--- a/docs/src/tutorials/conic/dualization.jl
+++ b/docs/src/tutorials/conic/dualization.jl
@@ -130,6 +130,7 @@ print(model_dual)
 
 set_optimizer(model_primal, SCS.Optimizer)
 optimize!(model_primal)
+@assert is_solved_and_feasible(model_primal; dual = true)
 
 # (There are five rows in the constraint matrix because SCS expects problems in
 # geometric conic form, and so JuMP has reformulated the `X, PSD` variable
@@ -153,6 +154,7 @@ objective_value(model_primal)
 
 set_optimizer(model_dual, SCS.Optimizer)
 optimize!(model_dual)
+@assert is_solved_and_feasible(model_dual; dual = true)
 
 # and the solution we obtain is:
 
@@ -182,6 +184,7 @@ objective_value(model_dual)
 
 set_optimizer(model_primal, Dualization.dual_optimizer(SCS.Optimizer))
 optimize!(model_primal)
+@assert is_solved_and_feasible(model_primal; dual = true)
 
 # The performance is the same as if we solved `model_dual`, and the correct
 # solution is returned to `X`:
@@ -197,6 +200,7 @@ dual.(primal_c)
 
 set_optimizer(model_dual, Dualization.dual_optimizer(SCS.Optimizer))
 optimize!(model_dual)
+@assert is_solved_and_feasible(model_dual; dual = true)
 
 #-
 
diff --git a/docs/src/tutorials/conic/ellipse_approx.jl b/docs/src/tutorials/conic/ellipse_approx.jl
index 6ebff779540..2b08720ca47 100644
--- a/docs/src/tutorials/conic/ellipse_approx.jl
+++ b/docs/src/tutorials/conic/ellipse_approx.jl
@@ -49,7 +49,7 @@ import LinearAlgebra
 import Plots
 import Random
 import SCS
-import Test  #src
+import Test
 
 # ## Data
 
@@ -110,8 +110,7 @@ m, n = size(S)
 @constraint(model, [t; vec(Z)] in MOI.RootDetConeSquare(n))
 @objective(model, Max, t)
 optimize!(model)
-Test.@test termination_status(model) == OPTIMAL    #src
-Test.@test primal_status(model) == FEASIBLE_POINT  #src
+Test.@test is_solved_and_feasible(model)
 solution_summary(model)
 
 # ## Results
@@ -211,6 +210,7 @@ f = [1 - S[i, :]' * Z * S[i, :] + 2 * S[i, :]' * z - s for i in 1:m]
 ## The former @objective(model, Max, t)
 @objective(model, Max, 1 * t + 0)
 optimize!(model)
+Test.@test is_solved_and_feasible(model)
 Test.@test isapprox(D, value.(Z); atol = 1e-6)  #src
 solve_time_1 = solve_time(model)
 
@@ -233,6 +233,7 @@ print_active_bridges(model)
 
 remove_bridge(model, MOI.Bridges.Constraint.GeoMeanToPowerBridge)
 optimize!(model)
+Test.@test is_solved_and_feasible(model)
 
 # This time, the solve took:
 
diff --git a/docs/src/tutorials/conic/experiment_design.jl b/docs/src/tutorials/conic/experiment_design.jl
index 02c07b15fab..75680c92dbd 100644
--- a/docs/src/tutorials/conic/experiment_design.jl
+++ b/docs/src/tutorials/conic/experiment_design.jl
@@ -141,6 +141,7 @@ for i in 1:q
 end
 @objective(aOpt, Min, sum(u))
 optimize!(aOpt)
+@assert is_solved_and_feasible(aOpt)
 objective_value(aOpt)
 
 #-
@@ -182,6 +183,7 @@ set_silent(eOpt)
 @constraint(eOpt, sum(np) <= n)
 @objective(eOpt, Max, t)
 optimize!(eOpt)
+@assert is_solved_and_feasible(eOpt)
 objective_value(eOpt)
 #-
 value.(np)
@@ -212,6 +214,7 @@ set_silent(dOpt)
 E = V * LinearAlgebra.diagm(0 => np ./ n) * V'
 @constraint(dOpt, [t; 1; triangle_vec(E)] in MOI.LogDetConeTriangle(q))
 optimize!(dOpt)
+@assert is_solved_and_feasible(dOpt)
 objective_value(dOpt)
 #-
 value.(np)
diff --git a/docs/src/tutorials/conic/logistic_regression.jl b/docs/src/tutorials/conic/logistic_regression.jl
index cdfa33fe74a..890bde86267 100644
--- a/docs/src/tutorials/conic/logistic_regression.jl
+++ b/docs/src/tutorials/conic/logistic_regression.jl
@@ -195,11 +195,12 @@ X, y = generate_dataset(n, p; shift = 10.0);
 model = build_logit_model(X, y, λ)
 set_optimizer(model, SCS.Optimizer)
 set_silent(model)
-JuMP.optimize!(model)
+optimize!(model)
+@assert is_solved_and_feasible(model)
 
 #-
 
-θ♯ = JuMP.value.(model[:θ])
+θ♯ = value.(model[:θ])
 
 # It appears that the speed of convergence is not that impacted by the correlation
 # of the dataset, nor by the penalty $\lambda$.
@@ -237,11 +238,12 @@ count_nonzero(v::Vector; tol = 1e-6) = sum(abs.(v) .>= tol)
 sparse_model = build_sparse_logit_model(X, y, λ)
 set_optimizer(sparse_model, SCS.Optimizer)
 set_silent(sparse_model)
-JuMP.optimize!(sparse_model)
+optimize!(sparse_model)
+@assert is_solved_and_feasible(sparse_model)
 
 #-
 
-θ♯ = JuMP.value.(sparse_model[:θ])
+θ♯ = value.(sparse_model[:θ])
 println(
     "Number of non-zero components: ",
     count_nonzero(θ♯),
diff --git a/docs/src/tutorials/conic/min_ellipse.jl b/docs/src/tutorials/conic/min_ellipse.jl
index 1fc33027a85..d6a6f4e74fb 100644
--- a/docs/src/tutorials/conic/min_ellipse.jl
+++ b/docs/src/tutorials/conic/min_ellipse.jl
@@ -45,7 +45,7 @@ using JuMP
 import LinearAlgebra
 import Plots
 import SCS
-import Test  #src
+import Test
 
 # ## Data
 
@@ -125,8 +125,7 @@ end
 # Now, solve the program:
 
 optimize!(model)
-Test.@test termination_status(model) == OPTIMAL    #src
-Test.@test primal_status(model) == FEASIBLE_POINT  #src
+Test.@test is_solved_and_feasible(model)
 solution_summary(model)
 
 # ## Results
diff --git a/docs/src/tutorials/conic/quantum_discrimination.jl b/docs/src/tutorials/conic/quantum_discrimination.jl
index d743cc0b005..252fd95d7f8 100644
--- a/docs/src/tutorials/conic/quantum_discrimination.jl
+++ b/docs/src/tutorials/conic/quantum_discrimination.jl
@@ -98,6 +98,7 @@ E = [@variable(model, [1:d, 1:d] in HermitianPSDCone()) for i in 1:N]
 # Now we optimize:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # The probability of guessing correctly is:
@@ -140,6 +141,7 @@ push!(E, E_N)
 # Then we can check that we get the same solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 #-
diff --git a/docs/src/tutorials/conic/simple_examples.jl b/docs/src/tutorials/conic/simple_examples.jl
index a89e8de8c79..4d5951f15ce 100644
--- a/docs/src/tutorials/conic/simple_examples.jl
+++ b/docs/src/tutorials/conic/simple_examples.jl
@@ -64,6 +64,7 @@ function solve_max_cut_sdp(weights)
     @objective(model, Max, 0.25 * LinearAlgebra.dot(L, X))
     @constraint(model, LinearAlgebra.diag(X) .== 1)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     V = svd_cholesky(value(X))
     Random.seed!(N)
     r = rand(N)
@@ -133,6 +134,7 @@ function example_k_means_clustering()
     @constraint(model, [i = 1:m], sum(Z[i, :]) .== 1)
     @constraint(model, LinearAlgebra.tr(Z) == num_clusters)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     Z_val = value.(Z)
     current_cluster, visited = 0, Set{Int}()
     solution = [1, 1, 2, 1, 2, 2]  #src
@@ -185,10 +187,12 @@ function example_correlation_problem()
     @constraint(model, 0.4 <= ρ["B", "C"] <= 0.5)
     @objective(model, Max, ρ["A", "C"])
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     println("An upper bound for ρ_AC is $(value(ρ["A", "C"]))")
     Test.@test value(ρ["A", "C"]) ≈ 0.87195 atol = 1e-4  #src
     @objective(model, Min, ρ["A", "C"])
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     println("A lower bound for ρ_AC is $(value(ρ["A", "C"]))")
     Test.@test value(ρ["A", "C"]) ≈ -0.978 atol = 1e-3  #src
     return
@@ -265,8 +269,7 @@ function example_minimum_distortion()
     fix(Q[1, 1], 0)
     @objective(model, Min, c²)
     optimize!(model)
-    Test.@test termination_status(model) == OPTIMAL
-    Test.@test primal_status(model) == FEASIBLE_POINT
+    Test.@test is_solved_and_feasible(model)
     Test.@test objective_value(model) ≈ 4 / 3 atol = 1e-4
     ## Recover the minimal distorted embedding:
     X = [zeros(3) sqrt(value.(Q)[2:end, 2:end])]
@@ -350,8 +353,7 @@ function example_theta_problem()
     J = ones(Int, 5, 5)
     @objective(model, Max, LinearAlgebra.dot(J, X))
     optimize!(model)
-    Test.@test termination_status(model) == OPTIMAL
-    Test.@test primal_status(model) == FEASIBLE_POINT
+    Test.@test is_solved_and_feasible(model)
     Test.@test objective_value(model) ≈ sqrt(5) rtol = 1e-4
     println("The Lovász number is: $(objective_value(model))")
     return
@@ -382,6 +384,7 @@ function example_robust_uncertainty_sets()
     @constraint(model, [((1-ɛ)/ɛ) (u - μ)'; (u-μ) Σ] >= 0, PSDCone())
     @objective(model, Max, c' * u)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     exact =
         μhat' * c +
         Γ1(𝛿 / 2, N) * LinearAlgebra.norm(c) +
diff --git a/docs/src/tutorials/conic/start_values.jl b/docs/src/tutorials/conic/start_values.jl
index d632fc62045..0ebc7d62f4c 100644
--- a/docs/src/tutorials/conic/start_values.jl
+++ b/docs/src/tutorials/conic/start_values.jl
@@ -67,6 +67,7 @@ model = Model(SCS.Optimizer)
 @constraint(model, sum(x) <= 1)
 @objective(model, Max, sum(i * x[i] for i in 1:3))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # By looking at the log, we can see that SCS took 75 iterations to find the optimal
 # solution. Now we set the optimal solution as our starting point:
diff --git a/docs/src/tutorials/conic/tips_and_tricks.jl b/docs/src/tutorials/conic/tips_and_tricks.jl
index 6d558a2116c..9de21fb0463 100644
--- a/docs/src/tutorials/conic/tips_and_tricks.jl
+++ b/docs/src/tutorials/conic/tips_and_tricks.jl
@@ -98,6 +98,7 @@ set_silent(model)
 @constraint(model, [t; x] in SecondOrderCone())
 @objective(model, Min, t)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(t), value.(x)
 
 # ## Rotated Second-Order Cone
@@ -120,6 +121,7 @@ set_silent(model)
 @constraint(model, [t; 0.5; residuals] in RotatedSecondOrderCone())
 @objective(model, Min, t)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(θ), value(t)
 
 # ## Exponential Cone
@@ -142,6 +144,7 @@ set_silent(model)
 @objective(model, Min, z)
 @constraint(model, [x, 1, z] in MOI.ExponentialCone())
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(z), exp(1.5)
 
 # ### Logarithm
@@ -155,6 +158,7 @@ set_silent(model)
 @objective(model, Max, x)
 @constraint(model, [x, 1, z] in MOI.ExponentialCone())
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(x), log(1.5)
 
 # ### Log-sum-exp
@@ -216,6 +220,7 @@ set_silent(model)
 @constraint(model, A * x .<= b)
 @constraint(model, [i = 1:n], [t[i], x[i], 1] in MOI.ExponentialCone())
 optimize!(model)
+@assert is_solved_and_feasible(model)
 objective_value(model)
 
 # The [`MOI.ExponentialCone`](@ref) has a dual, the [`MOI.DualExponentialCone`](@ref),
@@ -234,6 +239,7 @@ set_silent(model)
 @constraint(model, A * x .<= b)
 @constraint(model, [t; ones(n); x] in MOI.RelativeEntropyCone(2n + 1))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 objective_value(model)
 
 # ## PowerCone
@@ -255,6 +261,7 @@ set_silent(model)
 @constraint(model, [t, 1, x] in MOI.PowerCone(1 / 3))
 @objective(model, Min, t)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(t), value(x)
 
 # The [`MOI.PowerCone`](@ref) has a dual, the [`MOI.DualPowerCone`](@ref),
@@ -277,6 +284,7 @@ function p_norm(x::Vector, p)
     @constraint(model, sum(r) == t)
     @objective(model, Min, t)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value(t)
 end
 
@@ -320,6 +328,7 @@ set_silent(model)
 @objective(model, Min, t)
 @constraint(model, t .* I - A in PSDCone())
 optimize!(model)
+@assert is_solved_and_feasible(model)
 objective_value(model)
 
 # ## GeometricMeanCone
@@ -353,6 +362,7 @@ set_silent(model)
 @constraint(model, [t; vec(X)] in MOI.RootDetConeSquare(2))
 @constraint(model, X .== [2 1; 1 3])
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(t), sqrt(LinearAlgebra.det(value.(X)))
 
 # If `X` is symmetric, then you can use [`MOI.RootDetConeTriangle`](@ref)
@@ -390,6 +400,7 @@ set_silent(model)
 @constraint(model, X .== [2 1; 1 3])
 @constraint(model, u == 0.5)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))
 
 # If `X` is symmetric, then you can use [`MOI.LogDetConeTriangle`](@ref)
@@ -410,6 +421,7 @@ set_silent(model)
 @constraint(model, X .== [2 1; 1 3])
 @constraint(model, u == 0.5)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5))
 
 # ## Other Cones and Functions
diff --git a/docs/src/tutorials/getting_started/debugging.jl b/docs/src/tutorials/getting_started/debugging.jl
index a86a89d1d3f..8a6be8d6841 100644
--- a/docs/src/tutorials/getting_started/debugging.jl
+++ b/docs/src/tutorials/getting_started/debugging.jl
@@ -346,6 +346,7 @@ set_silent(model)
 # for variables with large positive or negative values in the optimal solution.
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 for var in all_variables(model)
     if var == objective
         continue
diff --git a/docs/src/tutorials/getting_started/design_patterns_for_larger_models.jl b/docs/src/tutorials/getting_started/design_patterns_for_larger_models.jl
index ea9530beccb..634365c2469 100644
--- a/docs/src/tutorials/getting_started/design_patterns_for_larger_models.jl
+++ b/docs/src/tutorials/getting_started/design_patterns_for_larger_models.jl
@@ -55,6 +55,7 @@ model = Model(HiGHS.Optimizer)
 @objective(model, Max, sum(profit[i] * x[i] for i in 1:N))
 @constraint(model, sum(weight[i] * x[i] for i in 1:N) <= capacity)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value.(x)
 
 # The benefits of this approach are:
@@ -87,6 +88,7 @@ function solve_knapsack_1(profit::Vector, weight::Vector, capacity::Real)
     @objective(model, Max, sum(profit[i] * x[i] for i in 1:N))
     @constraint(model, sum(weight[i] * x[i] for i in 1:N) <= capacity)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value.(x)
 end
 
@@ -159,6 +161,7 @@ function solve_knapsack_2(data::KnapsackData)
         sum(v.weight * x[k] for (k, v) in data.objects) <= data.capacity,
     )
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value.(x)
 end
 
@@ -230,6 +233,7 @@ function solve_knapsack_3(data::KnapsackData; binary_knapsack::Bool)
         sum(v.weight * x[k] for (k, v) in data.objects) <= data.capacity,
     )
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value.(x)
 end
 
@@ -272,6 +276,7 @@ function solve_knapsack_4(data::KnapsackData, config::AbstractConfiguration)
         sum(v.weight * x[k] for (k, v) in data.objects) <= data.capacity,
     )
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value.(x)
 end
 
@@ -359,6 +364,7 @@ function solve_knapsack_5(data::KnapsackData, config::AbstractConfiguration)
     add_knapsack_constraints(model, data, config)
     add_knapsack_objective(model, data, config)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value.(model[:x])
 end
 
@@ -384,7 +390,7 @@ function solve_knapsack_6(
     add_knapsack_constraints(model, data, config)
     add_knapsack_objective(model, data, config)
     optimize!(model)
-    if termination_status(model) != OPTIMAL
+    if !is_solved_and_feasible(model)
         @warn("Model not solved to optimality")
         return nothing
     end
@@ -519,7 +525,7 @@ function _solve_knapsack(
     _add_knapsack_constraints(model, data, config)
     _add_knapsack_objective(model, data, config)
     JuMP.optimize!(model)
-    if JuMP.termination_status(model) != JuMP.OPTIMAL
+    if !JuMP.is_solved_and_feasible(model)
         @warn("Model not solved to optimality")
         return nothing
     end
diff --git a/docs/src/tutorials/getting_started/getting_started_with_JuMP.jl b/docs/src/tutorials/getting_started/getting_started_with_JuMP.jl
index 8d04bfebab1..9814bdaa54d 100644
--- a/docs/src/tutorials/getting_started/getting_started_with_JuMP.jl
+++ b/docs/src/tutorials/getting_started/getting_started_with_JuMP.jl
@@ -178,7 +178,13 @@ optimize!(model)
 #     Julia has a convention that functions which mutate their arguments should
 #     end in `!`. A common example is `push!`.
 
-# Now let's see what information we can query about the solution.
+# Now let's see what information we can query about the solution,
+# starting with [`is_solved_and_feasible`](@ref):
+
+is_solved_and_feasible(model)
+
+# We can get more information about the solution by querying the three types of
+# statuses.
 
 # [`termination_status`](@ref) tells us why the solver stopped:
 
@@ -212,6 +218,17 @@ value(y)
 shadow_price(c1)
 shadow_price(c2)
 
+# !!! warning
+#     You should always check whether the solver found a solution before calling
+#     solution functions like [`value`](@ref) or [`objective_value`](@ref). A
+#     common workflow is:
+#     ```julia
+#     optimize!(model)
+#     if !is_solved_and_feasible(model)
+#         error("Solver did not find an optimal solution")
+#     end
+#     ```
+
 # That's it for our simple model. In the rest of this tutorial, we expand on
 # some of the basic JuMP operations.
 
@@ -269,7 +286,7 @@ function solve_infeasible()
     @constraint(model, x + y >= 3)
     @objective(model, Max, x + 2y)
     optimize!(model)
-    if termination_status(model) != OPTIMAL
+    if !is_solved_and_feasible(model)
         @warn("The model was not solved correctly.")
         return
     end
@@ -506,4 +523,5 @@ c = [1, 3, 5, 2]
 @constraint(vector_model, A * x .== b)
 @objective(vector_model, Min, c' * x)
 optimize!(vector_model)
+@assert is_solved_and_feasible(vector_model)
 objective_value(vector_model)
diff --git a/docs/src/tutorials/getting_started/getting_started_with_data_and_plotting.jl b/docs/src/tutorials/getting_started/getting_started_with_data_and_plotting.jl
index c977d57c5b9..e617457cc8a 100644
--- a/docs/src/tutorials/getting_started/getting_started_with_data_and_plotting.jl
+++ b/docs/src/tutorials/getting_started/getting_started_with_data_and_plotting.jl
@@ -364,6 +364,10 @@ optimize!(model)
 
 solution_summary(model)
 
+# Just to be sure, check that the solver found an optimal solution:
+
+@assert is_solved_and_feasible(model)
+
 # ### Solution
 
 # Let's  have a look at the solution in more detail:
diff --git a/docs/src/tutorials/linear/callbacks.jl b/docs/src/tutorials/linear/callbacks.jl
index f4517c8f374..b647337a507 100644
--- a/docs/src/tutorials/linear/callbacks.jl
+++ b/docs/src/tutorials/linear/callbacks.jl
@@ -13,7 +13,7 @@
 using JuMP
 import GLPK
 import Random
-import Test  #src
+import Test
 
 # !!! info
 #     This tutorial uses the [MathOptInterface](@ref moi_documentation) API.
@@ -60,11 +60,10 @@ function example_lazy_constraint()
     end
     set_attribute(model, MOI.LazyConstraintCallback(), my_callback_function)
     optimize!(model)
-    Test.@test termination_status(model) == OPTIMAL    #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
-    Test.@test lazy_called    #src
-    Test.@test value(x) == 1  #src
-    Test.@test value(y) == 2  #src
+    Test.@test is_solved_and_feasible(model)
+    Test.@test lazy_called
+    Test.@test value(x) == 1
+    Test.@test value(y) == 2
     println("Optimal solution (x, y) = ($(value(x)), $(value(y)))")
     return
 end
@@ -100,9 +99,8 @@ function example_user_cut_constraint()
     end
     set_attribute(model, MOI.UserCutCallback(), my_callback_function)
     optimize!(model)
-    Test.@test termination_status(model) == OPTIMAL  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
-    Test.@test callback_called  #src
+    Test.@test is_solved_and_feasible(model)
+    Test.@test callback_called
     @show callback_called
     return
 end
@@ -128,16 +126,15 @@ function example_heuristic_solution()
         ret =
             MOI.submit(model, MOI.HeuristicSolution(cb_data), x, floor.(x_vals))
         println("Heuristic solution status = $(ret)")
-        Test.@test ret in (                     #src
-            MOI.HEURISTIC_SOLUTION_ACCEPTED,    #src
-            MOI.HEURISTIC_SOLUTION_REJECTED,    #src
-        )                                       #src
+        Test.@test ret in (
+            MOI.HEURISTIC_SOLUTION_ACCEPTED,
+            MOI.HEURISTIC_SOLUTION_REJECTED,
+        )
     end
     set_attribute(model, MOI.HeuristicCallback(), my_callback_function)
     optimize!(model)
-    Test.@test termination_status(model) == OPTIMAL  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
-    Test.@test callback_called  #src
+    Test.@test is_solved_and_feasible(model)
+    Test.@test callback_called
     return
 end
 
@@ -174,11 +171,10 @@ function example_solver_dependent_callback()
     end
     set_attribute(model, GLPK.CallbackFunction(), my_callback_function)
     optimize!(model)
-    Test.@test termination_status(model) == OPTIMAL  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
-    Test.@test lazy_called  #src
-    Test.@test value(x) == 1  #src
-    Test.@test value(y) == 2  #src
+    Test.@test is_solved_and_feasible(model)
+    Test.@test lazy_called
+    Test.@test value(x) == 1
+    Test.@test value(y) == 2
     return
 end
 
diff --git a/docs/src/tutorials/linear/cannery.jl b/docs/src/tutorials/linear/cannery.jl
index aaf3c3b344e..13fae00c5a1 100644
--- a/docs/src/tutorials/linear/cannery.jl
+++ b/docs/src/tutorials/linear/cannery.jl
@@ -21,19 +21,19 @@
 using JuMP
 import HiGHS
 import JSON
-import Test  #src
+import Test
 
 # ## Formulation
 
-# The cannery problem assumes we are optimizing the shipment of cases of 
+# The cannery problem assumes we are optimizing the shipment of cases of
 # cans from production plants ``p \in P`` to markets ``m \in M``.
 
 # Each production plant ``p`` has a capacity ``c_p``, and each market ``m``
-# has a demand ``d_m``. The shipping cost per case of cans from plant ``p`` 
+# has a demand ``d_m``. The shipping cost per case of cans from plant ``p``
 # to market ``m`` is ``d_{p,m}``.
 
 # We wish to find the distribution plan ``x_{p,m}``, the number of cases of cans
-# to ship from plant ``p`` to market ``m``, for ``p \in P`` and ``m \in M``  
+# to ship from plant ``p`` to market ``m``, for ``p \in P`` and ``m \in M``
 # that minimizes the shipping costs. We can formulate our problem as the
 # following linear program:
 # ```math
@@ -121,10 +121,8 @@ solution_summary(model)
 
 # What's the optimal shipment?
 
+Test.@test is_solved_and_feasible(model)
+Test.@test isapprox(objective_value(model), 1_680.0, atol = 1e-6)  #src
 for p in P, m in M
     println(p, " => ", m, ": ", value(x[p, m]))
 end
-
-Test.@test termination_status(model) == OPTIMAL    #src
-Test.@test primal_status(model) == FEASIBLE_POINT  #src
-Test.@test objective_value(model) == 1_680.0           #src
diff --git a/docs/src/tutorials/linear/constraint_programming.jl b/docs/src/tutorials/linear/constraint_programming.jl
index ed47265c03f..a420186effa 100644
--- a/docs/src/tutorials/linear/constraint_programming.jl
+++ b/docs/src/tutorials/linear/constraint_programming.jl
@@ -29,6 +29,7 @@ set_silent(model)
 @variable(model, 1 <= x[1:4] <= 4, Int)
 @constraint(model, x in MOI.AllDifferent(4))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value.(x)
 
 # ## BinPacking
@@ -44,6 +45,7 @@ set_silent(model)
 @variable(model, 1 <= x[1:length(weights)] <= number_of_bins, Int)
 @constraint(model, x in MOI.BinPacking(capacity, weights))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value.(x)
 
 # Here, the value of `x[i]` is the bin that item `i` was placed into.
@@ -59,6 +61,7 @@ set_silent(model)
 @variable(model, x[1:4], Int)
 @constraint(model, x in MOI.Circuit(4))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # Let's see what tour was found, starting at node number `1`:
 y = round.(Int, value.(x))
@@ -112,6 +115,7 @@ n = 1
 # Let's check that we found a valid solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value.(x)
 
 # ## CountBelongs
@@ -130,6 +134,7 @@ set_silent(model)
 set = Set([2, 3])
 @constraint(model, [n; x] in MOI.CountBelongs(1 + length(x), set))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(n), value.(x)
 
 # ## CountDistinct
@@ -144,6 +149,7 @@ set_silent(model)
 @objective(model, Max, sum(x))
 @constraint(model, [n; x] in MOI.CountDistinct(1 + length(x)))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(n), value.(x)
 
 # ## CountGreaterThan
@@ -163,6 +169,7 @@ set_silent(model)
 @objective(model, Max, sum(x))
 @constraint(model, [n; y; x] in MOI.CountGreaterThan(1 + 1 + length(x)))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value(n), value(y), value.(x)
 
 # Here `n` is strictly greater than the count, and there is no limit on how
@@ -187,4 +194,5 @@ set_silent(model)
 @variable(model, x[i = 1:3], Int)
 @constraint(model, x in MOI.Table(table))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 value.(x)
diff --git a/docs/src/tutorials/linear/diet.jl b/docs/src/tutorials/linear/diet.jl
index 6e9ef81915a..413608f942e 100644
--- a/docs/src/tutorials/linear/diet.jl
+++ b/docs/src/tutorials/linear/diet.jl
@@ -16,7 +16,7 @@ using JuMP
 import CSV
 import DataFrames
 import HiGHS
-import Test  #hide
+import Test
 
 # ## Formulation
 
@@ -145,7 +145,7 @@ print(model)
 # Let's optimize and take a look at the solution:
 
 optimize!(model)
-Test.@test primal_status(model) == FEASIBLE_POINT        #hide
+@assert is_solved_and_feasible(model)
 Test.@test objective_value(model) ≈ 11.8288 atol = 1e-4  #hide
 solution_summary(model)
 
@@ -178,8 +178,9 @@ dairy_foods = ["milk", "ice cream"]
 is_dairy = map(name -> name in dairy_foods, foods.name)
 dairy_constraint = @constraint(model, sum(foods[is_dairy, :x]) <= 6)
 optimize!(model)
-Test.@test termination_status(model) == INFEASIBLE  #hide
-Test.@test primal_status(model) == NO_SOLUTION      #hide
+Test.@test !is_solved_and_feasible(model)
+Test.@test termination_status(model) == INFEASIBLE
+Test.@test primal_status(model) == NO_SOLUTION
 solution_summary(model)
 
 # There exists no feasible solution to our problem. Looks like we're stuck
diff --git a/docs/src/tutorials/linear/facility_location.jl b/docs/src/tutorials/linear/facility_location.jl
index b435f45f0fc..8b87955f072 100644
--- a/docs/src/tutorials/linear/facility_location.jl
+++ b/docs/src/tutorials/linear/facility_location.jl
@@ -130,6 +130,7 @@ set_silent(model)
 # Solve the uncapacitated facility location problem with HiGHS
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 println("Optimal value: ", objective_value(model))
 
 # ### Visualizing the solution
@@ -256,6 +257,7 @@ set_silent(model)
 # Solve the problem
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 println("Optimal value: ", objective_value(model))
 
 # ### Visualizing the solution
diff --git a/docs/src/tutorials/linear/factory_schedule.jl b/docs/src/tutorials/linear/factory_schedule.jl
index f034a93532b..7e2e96f3def 100644
--- a/docs/src/tutorials/linear/factory_schedule.jl
+++ b/docs/src/tutorials/linear/factory_schedule.jl
@@ -186,6 +186,7 @@ function solve_factory_scheduling(
         )
     )
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     schedules = Dict{Symbol,Vector{Float64}}(
         Symbol(f) => value.(production[:, f]) for f in factories
     )
diff --git a/docs/src/tutorials/linear/finance.jl b/docs/src/tutorials/linear/finance.jl
index 2e376dfd308..ae95dc92336 100644
--- a/docs/src/tutorials/linear/finance.jl
+++ b/docs/src/tutorials/linear/finance.jl
@@ -92,7 +92,7 @@ end)
 )
 
 optimize!(financing)
-
+@assert is_solved_and_feasible(financing)
 objective_value(financing)
 
 # ## Combinatorial auctions
@@ -136,9 +136,8 @@ auction = Model(HiGHS.Optimizer)
 for i in 1:6
     @constraint(auction, sum(y[j] for j in 1:6 if i in bid_items[j]) <= 1)
 end
-
 optimize!(auction)
-
+@assert is_solved_and_feasible(auction)
 objective_value(auction)
 
 #-
diff --git a/docs/src/tutorials/linear/geographic_clustering.jl b/docs/src/tutorials/linear/geographic_clustering.jl
index 448ae2e53df..681f7cf9155 100644
--- a/docs/src/tutorials/linear/geographic_clustering.jl
+++ b/docs/src/tutorials/linear/geographic_clustering.jl
@@ -151,6 +151,7 @@ end
 # We can then call `optimize!` and review the results.
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # ### Reviewing the Results
 
diff --git a/docs/src/tutorials/linear/knapsack.jl b/docs/src/tutorials/linear/knapsack.jl
index 6e85075a29a..24077e416e6 100644
--- a/docs/src/tutorials/linear/knapsack.jl
+++ b/docs/src/tutorials/linear/knapsack.jl
@@ -96,6 +96,7 @@ print(model)
 # We can now solve the optimization problem and inspect the results.
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # The items chosen are
@@ -124,8 +125,7 @@ function solve_knapsack_problem(;
     @objective(model, Max, profit' * x)
     @constraint(model, weight' * x <= capacity)
     optimize!(model)
-    @assert termination_status(model) == OPTIMAL
-    @assert primal_status(model) == FEASIBLE_POINT
+    @assert is_solved_and_feasible(model)
     println("Objective is: ", objective_value(model))
     println("Solution is:")
     for i in 1:n
diff --git a/docs/src/tutorials/linear/lp_sensitivity.jl b/docs/src/tutorials/linear/lp_sensitivity.jl
index 5d9e3deecf1..0bc608cf19a 100644
--- a/docs/src/tutorials/linear/lp_sensitivity.jl
+++ b/docs/src/tutorials/linear/lp_sensitivity.jl
@@ -39,6 +39,7 @@ model = Model(HiGHS.Optimizer)
 @constraint(model, c2, 7x + 12y >= 120)
 @constraint(model, c3, x + y <= 20)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model; verbose = true)
 
 # Can you identify:
diff --git a/docs/src/tutorials/linear/mip_duality.jl b/docs/src/tutorials/linear/mip_duality.jl
index cb03cacf41b..b82dd1ab379 100644
--- a/docs/src/tutorials/linear/mip_duality.jl
+++ b/docs/src/tutorials/linear/mip_duality.jl
@@ -58,6 +58,7 @@ print(model)
 # If we optimize this model, we obtain a [`dual_status`](@ref) of [`NO_SOLUTION`](@ref):
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 dual_status(model)
 
 # This is because HiGHS cannot compute the duals of a mixed-integer program. We
@@ -73,6 +74,7 @@ print(model)
 # dual:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 dual_status(model)
 
 # and a marginal price of electricity of \$100/MWh:
@@ -94,6 +96,7 @@ print(model)
 # the [`fix_discrete_variables`](@ref) function:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 dual_status(model)
 
 #-
@@ -113,6 +116,7 @@ print(model)
 #-
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 dual_status(model)
 
 #-
diff --git a/docs/src/tutorials/linear/multi.jl b/docs/src/tutorials/linear/multi.jl
index 91cef976613..3916ac2cb34 100644
--- a/docs/src/tutorials/linear/multi.jl
+++ b/docs/src/tutorials/linear/multi.jl
@@ -24,7 +24,7 @@ import DataFrames
 import HiGHS
 import SQLite
 import Tables
-import Test  #src
+import Test
 
 const DBInterface = SQLite.DBInterface
 
@@ -177,8 +177,7 @@ end
 # Finally, we can optimize the model:
 
 optimize!(model)
-Test.@test termination_status(model) == OPTIMAL     #src
-Test.@test primal_status(model) == FEASIBLE_POINT   #src
+Test.@test is_solved_and_feasible(model)
 Test.@test objective_value(model) == 225_700.0      #src
 solution_summary(model)
 
diff --git a/docs/src/tutorials/linear/multi_commodity_network.jl b/docs/src/tutorials/linear/multi_commodity_network.jl
index 3f6f234a76b..124d1e9af66 100644
--- a/docs/src/tutorials/linear/multi_commodity_network.jl
+++ b/docs/src/tutorials/linear/multi_commodity_network.jl
@@ -20,7 +20,7 @@ import DataFrames
 import HiGHS
 import SQLite
 import SQLite.DBInterface
-import Test  #src
+import Test
 
 # ## Formulation
 
@@ -201,8 +201,7 @@ df = DataFrames.leftjoin(
 # Finally, we can optimize the model:
 
 optimize!(model)
-Test.@test termination_status(model) == OPTIMAL     #src
-Test.@test primal_status(model) == FEASIBLE_POINT   #src
+Test.@test is_solved_and_feasible(model)
 solution_summary(model)
 
 # update the solution in the DataFrames:
diff --git a/docs/src/tutorials/linear/multi_objective_examples.jl b/docs/src/tutorials/linear/multi_objective_examples.jl
index b4d9c6fa170..cbca784b18e 100644
--- a/docs/src/tutorials/linear/multi_objective_examples.jl
+++ b/docs/src/tutorials/linear/multi_objective_examples.jl
@@ -37,6 +37,7 @@ solution_summary(model)
 #-
 
 for i in 1:result_count(model)
+    @assert is_solved_and_feasible(model; result = i)
     print(i, ": z = ", round.(Int, objective_value(model; result = i)), " | ")
     println("x = ", value.([x1, x2]; result = i))
 end
@@ -65,6 +66,7 @@ solution_summary(model)
 #-
 
 for i in 1:result_count(model)
+    @assert is_solved_and_feasible(model; result = i)
     print(i, ": z = ", round.(Int, objective_value(model; result = i)), " | ")
     println("x = ", round.(Int, value.(x; result = i)))
 end
@@ -109,6 +111,7 @@ solution_summary(model)
 #-
 
 for i in 1:result_count(model)
+    @assert is_solved_and_feasible(model; result = i)
     print(i, ": z = ", round.(Int, objective_value(model; result = i)), " | ")
     X = round.(Int, value.(x; result = i))
     print("Path:")
diff --git a/docs/src/tutorials/linear/multi_objective_knapsack.jl b/docs/src/tutorials/linear/multi_objective_knapsack.jl
index e8f81f92528..a95b98497cf 100644
--- a/docs/src/tutorials/linear/multi_objective_knapsack.jl
+++ b/docs/src/tutorials/linear/multi_objective_knapsack.jl
@@ -121,6 +121,7 @@ set_attribute(model, MOA.Algorithm(), MOA.EpsilonConstraint())
 # Let's solve the problem and see the solution
 
 optimize!(model)
+@assert termination_status(model) == OPTIMAL
 solution_summary(model)
 
 # There are 9 solutions available. We can also use [`result_count`](@ref) to see
@@ -137,6 +138,14 @@ solution_summary(model; result = 5)
 
 #-
 
+@assert primal_status(model; result = 5) == FEASIBLE_POINT
+
+#-
+
+@assert is_solved_and_feasible(model; result = 5)
+
+#-
+
 objective_value(model; result = 5)
 
 # Note that because we set a vector of two objective functions, the objective
diff --git a/docs/src/tutorials/linear/n-queens.jl b/docs/src/tutorials/linear/n-queens.jl
index 36fba85372f..c75ff6d824c 100644
--- a/docs/src/tutorials/linear/n-queens.jl
+++ b/docs/src/tutorials/linear/n-queens.jl
@@ -66,6 +66,7 @@ end
 # a feasible solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # We can now review the solution that our model found:
 
diff --git a/docs/src/tutorials/linear/network_flows.jl b/docs/src/tutorials/linear/network_flows.jl
index b87030a028e..89c3e98391f 100644
--- a/docs/src/tutorials/linear/network_flows.jl
+++ b/docs/src/tutorials/linear/network_flows.jl
@@ -79,6 +79,7 @@ set_silent(shortest_path)
 @constraint(shortest_path, [i = 1:n], sum(x[i, :]) - sum(x[:, i]) == b[i],)
 @objective(shortest_path, Min, sum(G .* x))
 optimize!(shortest_path)
+@assert is_solved_and_feasible(shortest_path)
 objective_value(shortest_path)
 #-
 value.(x)
@@ -123,6 +124,7 @@ set_silent(assignment)
 @constraint(assignment, [j = 1:n], sum(y[j, :]) == 1)
 @objective(assignment, Max, sum(G .* y))
 optimize!(assignment)
+@assert is_solved_and_feasible(assignment)
 objective_value(assignment)
 #-
 value.(y)
@@ -163,6 +165,7 @@ max_flow = Model(HiGHS.Optimizer)
 @constraint(max_flow, [i = 1:n; i != 1 && i != 8], sum(f[i, :]) == sum(f[:, i]))
 @objective(max_flow, Max, sum(f[1, :]))
 optimize!(max_flow)
+@assert is_solved_and_feasible(max_flow)
 objective_value(max_flow)
 #-
 value.(f)
diff --git a/docs/src/tutorials/linear/piecewise_linear.jl b/docs/src/tutorials/linear/piecewise_linear.jl
index 493604437c0..6ad6fff89ee 100644
--- a/docs/src/tutorials/linear/piecewise_linear.jl
+++ b/docs/src/tutorials/linear/piecewise_linear.jl
@@ -52,6 +52,7 @@ function outer_approximate_x_squared(x̄)
     @objective(model, Min, y)
     @constraint(model, x == x̄)  # <-- a trivial constraint just for testing.
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value(y)
 end
 
@@ -102,6 +103,7 @@ function outer_approximate_log(x̄)
     @objective(model, Max, y)
     @constraint(model, x == x̄)  # <-- a trivial constraint just for testing.
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value(y)
 end
 
@@ -167,6 +169,7 @@ function inner_approximate_x_squared(x̄)
     @objective(model, Min, y)
     @constraint(model, x == x̄)  # <-- a trivial constraint just for testing.
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value(y)
 end
 
@@ -209,6 +212,7 @@ function inner_approximate_log(x̄)
     @objective(model, Max, y)
     @constraint(model, x == x̄)  # <-- a trivial constraint just for testing.
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value(y)
 end
 
@@ -262,6 +266,7 @@ function piecewise_linear_sin(x̄)
     end)
     @constraint(model, x == x̄)  # <-- a trivial constraint just for testing.
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     return value(y)
 end
 
diff --git a/docs/src/tutorials/linear/sudoku.jl b/docs/src/tutorials/linear/sudoku.jl
index 8e2c207188c..b13b6d2389f 100644
--- a/docs/src/tutorials/linear/sudoku.jl
+++ b/docs/src/tutorials/linear/sudoku.jl
@@ -134,6 +134,7 @@ end
 
 # solve problem
 optimize!(sudoku)
+@assert is_solved_and_feasible(sudoku)
 
 # Extract the values of x
 x_val = value.(x);
@@ -202,6 +203,7 @@ for i in 1:9, j in 1:9
 end
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # Display the solution
 
diff --git a/docs/src/tutorials/linear/transp.jl b/docs/src/tutorials/linear/transp.jl
index dd5ab0e253e..a62810ca3e8 100644
--- a/docs/src/tutorials/linear/transp.jl
+++ b/docs/src/tutorials/linear/transp.jl
@@ -120,6 +120,7 @@ function solve_transportation_problem(data::Containers.DenseAxisArray)
     @constraint(model, [o in O], sum(x[o, :]) <= data[o, "SUPPLY"])
     @constraint(model, [d in D], sum(x[:, d]) == data["DEMAND", d])
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     ## Pretty print the solution in the format of the input
     print("    ", join(lpad.(D, 7, ' ')))
     for o in O
diff --git a/docs/src/tutorials/nonlinear/classifiers.jl b/docs/src/tutorials/nonlinear/classifiers.jl
index 73a47342e1a..351742b74c6 100644
--- a/docs/src/tutorials/nonlinear/classifiers.jl
+++ b/docs/src/tutorials/nonlinear/classifiers.jl
@@ -25,7 +25,7 @@ import Ipopt
 import LinearAlgebra
 import Plots
 import Random
-import Test #src
+import Test
 
 # ## Data and visualisation
 
@@ -127,8 +127,7 @@ function solve_SVM_classifier(P::Matrix, labels::Vector; C::Float64 = C_0)
     D = LinearAlgebra.Diagonal(labels)
     @constraint(model, D * (P * w .- g) .+ y .>= 1)
     optimize!(model)
-    Test.@test termination_status(model) == LOCALLY_SOLVED  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
+    Test.@test is_solved_and_feasible(model)
     slack = extrema(value.(y))
     println("Minimum slack: ", slack[1], "\nMaximum slack: ", slack[2])
     classifier(x) = line(x; w = value.(w), g = value(g))
@@ -234,8 +233,7 @@ function solve_dual_SVM_classifier(P::Matrix, labels::Vector; C::Float64 = C_0)
     @objective(model, Min, 1 / 2 * u' * D * P * P' * D * u - sum(u))
     @constraint(model, con, sum(D * u) == 0)
     optimize!(model)
-    Test.@test termination_status(model) == LOCALLY_SOLVED  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
+    Test.@test is_solved_and_feasible(model)
     w = P' * D * value.(u)
     g = dual(con)
     classifier(x) = line(x; w = w, g = g)
@@ -322,8 +320,7 @@ function solve_kernel_SVM_classifier(
     con = @constraint(model, sum(D * u) == 0)
     @objective(model, Min, 1 / 2 * u' * D * K * D * u - sum(u))
     optimize!(model)
-    Test.@test termination_status(model) == LOCALLY_SOLVED  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
+    Test.@test is_solved_and_feasible(model)
     u_sol, g_sol = value.(u), dual(con)
     function classifier(v::Vector)
         return sum(
diff --git a/docs/src/tutorials/nonlinear/complementarity.jl b/docs/src/tutorials/nonlinear/complementarity.jl
index 3b5024a018a..3dde3f22350 100644
--- a/docs/src/tutorials/nonlinear/complementarity.jl
+++ b/docs/src/tutorials/nonlinear/complementarity.jl
@@ -47,6 +47,7 @@ set_silent(model)
 @variable(model, 0 <= x[1:4] <= 10, start = 0)
 @constraint(model, M * x + q ⟂ x)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test value.(x) ≈ [2.8, 0.0, 0.8, 1.2]  #src
 value.(x)
 
@@ -67,6 +68,7 @@ set_silent(model)
 @constraint(model, w + 2x - 2y + 4z - 6 ⟂ z)
 @constraint(model, w - x + 2y - 2z - 2 ⟂ y)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test value.([w, x, y, z]) ≈ [2.8, 0.0, 0.8, 1.2]  #src
 value.([w, x, y, z])
 
@@ -102,6 +104,7 @@ set_silent(model)
     end
 )
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test isapprox(value(p["new-york"]), 0.225; atol = 1e-3)  #src
 value.(p)
 
@@ -137,6 +140,7 @@ set_silent(model)
     end
 )
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test isapprox(value(C_G), 0.996; atol = 1e-3)  #src
 value(K)
 
@@ -189,6 +193,7 @@ set_silent(model)
 ## Production does not exceed capacity
 @constraint(model, [ω = 1:5], x - Y[ω] ⟂ μ[ω])
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # An equilibrium solution is to build 389 MW:
diff --git a/docs/src/tutorials/nonlinear/nested_problems.jl b/docs/src/tutorials/nonlinear/nested_problems.jl
index 5ee010472e3..27b069469e9 100644
--- a/docs/src/tutorials/nonlinear/nested_problems.jl
+++ b/docs/src/tutorials/nonlinear/nested_problems.jl
@@ -86,7 +86,7 @@ function solve_lower_level(x...)
     )
     @constraint(model, (y[1] - 10)^2 + (y[2] - 10)^2 <= 25)
     optimize!(model)
-    @assert termination_status(model) == LOCALLY_SOLVED
+    @assert is_solved_and_feasible(model)
     return objective_value(model), value.(y)
 end
 
@@ -149,6 +149,7 @@ model = Model(Ipopt.Optimizer)
 @operator(model, op_V, 2, V, ∇V, ∇²V)
 @objective(model, Min, x[1]^2 + x[2]^2 + op_V(x[1], x[2]))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # The optimal objective value is:
@@ -228,6 +229,7 @@ cache = Cache(Float64[], NaN, Float64[])
 )
 @objective(model, Min, x[1]^2 + x[2]^2 + op_cached_f(x[1], x[2]))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # an we can check we get the same objective value:
diff --git a/docs/src/tutorials/nonlinear/portfolio.jl b/docs/src/tutorials/nonlinear/portfolio.jl
index 666d92b8ab9..420c13bacf0 100644
--- a/docs/src/tutorials/nonlinear/portfolio.jl
+++ b/docs/src/tutorials/nonlinear/portfolio.jl
@@ -158,6 +158,7 @@ set_silent(model)
 @constraint(model, sum(x) <= 1000)
 @constraint(model, r' * x >= 50)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # The optimal allocation of our assets is:
@@ -209,6 +210,7 @@ set_optimizer_attribute(model, MOA.SolutionLimit(), 50)
 ## a single objective sense `Min`, and negate any `Max` objectives:
 @objective(model, Min, [variance, -expected_return])
 optimize!(model)
+@assert termination_status(model) == OPTIMAL
 solution_summary(model)
 
 # The algorithm found 50 different solutions. Let's plot them to see how they
diff --git a/docs/src/tutorials/nonlinear/querying_hessians.jl b/docs/src/tutorials/nonlinear/querying_hessians.jl
index d4a44ee72cc..e740a33e829 100644
--- a/docs/src/tutorials/nonlinear/querying_hessians.jl
+++ b/docs/src/tutorials/nonlinear/querying_hessians.jl
@@ -71,6 +71,7 @@ set_silent(model)
 @constraint(model, g_2, (x[1] + x[2])^2 <= 2)
 @objective(model, Min, (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2)
 optimize!(model)
+@assert is_solved_and_feasible(model)
 
 # ## The analytic solution
 
diff --git a/docs/src/tutorials/nonlinear/rocket_control.jl b/docs/src/tutorials/nonlinear/rocket_control.jl
index 33b3895b834..478e8319c87 100644
--- a/docs/src/tutorials/nonlinear/rocket_control.jl
+++ b/docs/src/tutorials/nonlinear/rocket_control.jl
@@ -127,6 +127,7 @@ ddt(x::Vector, t::Int) = (x[t] - x[t-1]) / Δt
 # Now we optimize the model and check that we found a solution:
 
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model)
 
 # Finally, we plot the solution:
diff --git a/docs/src/tutorials/nonlinear/simple_examples.jl b/docs/src/tutorials/nonlinear/simple_examples.jl
index 3a8ae666cba..1eade4f5462 100644
--- a/docs/src/tutorials/nonlinear/simple_examples.jl
+++ b/docs/src/tutorials/nonlinear/simple_examples.jl
@@ -25,8 +25,7 @@ function example_rosenbrock()
     @variable(model, y)
     @objective(model, Min, (1 - x)^2 + 100 * (y - x^2)^2)
     optimize!(model)
-    Test.@test termination_status(model) == LOCALLY_SOLVED
-    Test.@test primal_status(model) == FEASIBLE_POINT
+    Test.@test is_solved_and_feasible(model)
     Test.@test objective_value(model) ≈ 0.0 atol = 1e-10
     Test.@test value(x) ≈ 1.0
     Test.@test value(y) ≈ 1.0
@@ -87,8 +86,7 @@ function example_clnlbeam()
     primal_status      = $(primal_status(model))
     objective_value    = $(objective_value(model))
     """)
-    Test.@test termination_status(model) == LOCALLY_SOLVED  #src
-    Test.@test primal_status(model) == FEASIBLE_POINT  #src
+    Test.@test is_solved_and_feasible(model)
     Test.@test objective_value(model) ≈ 350.0  #src
     return
 end
@@ -116,6 +114,7 @@ function example_mle()
         sum((data[i] - μ)^2 for i in 1:n) / (2 * σ^2)
     )
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     println("μ             = ", value(μ))
     println("mean(data)    = ", Statistics.mean(data))
     println("σ^2           = ", value(σ)^2)
@@ -126,6 +125,7 @@ function example_mle()
     ## You can even do constrained MLE!
     @constraint(model, μ == σ^2)
     optimize!(model)
+    @assert is_solved_and_feasible(model)
     Test.@test value(μ) ≈ value(σ)^2
     println()
     println("With constraint μ == σ^2:")
@@ -153,12 +153,11 @@ function example_qcp()
     @constraint(model, x * x + y * y - z * z <= 0)
     @constraint(model, x * x - y * z <= 0)
     optimize!(model)
+    Test.@test is_solved_and_feasible(model)
     print(model)
     println("Objective value: ", objective_value(model))
     println("x = ", value(x))
     println("y = ", value(y))
-    Test.@test termination_status(model) == LOCALLY_SOLVED
-    Test.@test primal_status(model) == FEASIBLE_POINT
     Test.@test objective_value(model) ≈ 0.32699 atol = 1e-5
     Test.@test value(x) ≈ 0.32699 atol = 1e-5
     Test.@test value(y) ≈ 0.25707 atol = 1e-5
diff --git a/docs/src/tutorials/nonlinear/space_shuttle_reentry_trajectory.jl b/docs/src/tutorials/nonlinear/space_shuttle_reentry_trajectory.jl
index a775ee2c509..c6725f9ad32 100644
--- a/docs/src/tutorials/nonlinear/space_shuttle_reentry_trajectory.jl
+++ b/docs/src/tutorials/nonlinear/space_shuttle_reentry_trajectory.jl
@@ -303,7 +303,7 @@ end
 
 set_silent(model)  # Hide solver's verbose output
 optimize!(model)  # Solve for the control and state
-@assert termination_status(model) == LOCALLY_SOLVED
+@assert is_solved_and_feasible(model)
 
 ## Show final cross-range of the solution
 println(
diff --git a/docs/src/tutorials/nonlinear/tips_and_tricks.jl b/docs/src/tutorials/nonlinear/tips_and_tricks.jl
index 26f96d0d97f..a51f5e6108a 100644
--- a/docs/src/tutorials/nonlinear/tips_and_tricks.jl
+++ b/docs/src/tutorials/nonlinear/tips_and_tricks.jl
@@ -52,6 +52,7 @@ set_silent(model)
 @constraint(model, op_foo_2(x[1], x[2]) <= 2)
 function_calls = 0
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test objective_value(model) ≈ √3 atol = 1e-4
 Test.@test value.(x) ≈ [1.0, 1.0] atol = 1e-4
 println("Naive approach: function calls = $(function_calls)")
@@ -120,6 +121,7 @@ set_silent(model)
 @constraint(model, op_foo_2(x[1], x[2]) <= 2)
 function_calls = 0
 optimize!(model)
+@assert is_solved_and_feasible(model)
 Test.@test objective_value(model) ≈ √3 atol = 1e-4
 Test.@test value.(x) ≈ [1.0, 1.0] atol = 1e-4
 println("Memoized approach: function_calls = $(function_calls)")
diff --git a/docs/src/tutorials/nonlinear/user_defined_hessians.jl b/docs/src/tutorials/nonlinear/user_defined_hessians.jl
index b7f75e42d20..f2e272229d8 100644
--- a/docs/src/tutorials/nonlinear/user_defined_hessians.jl
+++ b/docs/src/tutorials/nonlinear/user_defined_hessians.jl
@@ -72,4 +72,5 @@ model = Model(Ipopt.Optimizer)
 @operator(model, op_rosenbrock, 2, rosenbrock, ∇rosenbrock, ∇²rosenbrock)
 @objective(model, Min, op_rosenbrock(x[1], x[2]))
 optimize!(model)
+@assert is_solved_and_feasible(model)
 solution_summary(model; verbose = true)
diff --git a/src/optimizer_interface.jl b/src/optimizer_interface.jl
index 533b4ed0c5b..aa42f18186b 100644
--- a/src/optimizer_interface.jl
+++ b/src/optimizer_interface.jl
@@ -582,6 +582,63 @@ function dual_status(model::GenericModel; result::Int = 1)
     return MOI.get(model, MOI.DualStatus(result))::MOI.ResultStatusCode
 end
 
+"""
+    is_solved_and_feasible(
+        model::GenericModel;
+        allow_local::Bool = true,
+        allow_almost::Bool = false,
+        dual::Bool = false,
+        result::Int = 1,
+    )
+
+Return `true` if the model has a feasible primal solution associated with result
+index `result` and the [`termination_status`](@ref) is [`OPTIMAL`](@ref) (the
+solver found a global optimum) or [`LOCALLY_SOLVED`](@ref) (the solver found a
+local optimum, which may also be the global optimum, but the solver could not
+prove so).
+
+If `allow_local = false`, then this function returns `true` only if the
+[`termination_status`](@ref) is [`OPTIMAL`](@ref).
+
+If `allow_almost = true`, then the [`termination_status`](@ref) may additionally
+be [`ALMOST_OPTIMAL`](@ref) or [`ALMOST_LOCALLY_SOLVED`](@ref) (if `allow_local`),
+and the [`primal_status`](@ref) and [`dual_status`](@ref) may additionally be
+[`NEARLY_FEASIBLE_POINT`](@ref).
+
+If `dual`, additionally check that an optimal dual solution is available.
+
+If this function returns `false`, use [`termination_status`](@ref),
+[`result_count`](@ref), [`primal_status`](@ref) and [`dual_status`](@ref) to
+understand what solutions are available (if any).
+"""
+function is_solved_and_feasible(
+    model::GenericModel;
+    dual::Bool = false,
+    allow_local::Bool = true,
+    allow_almost::Bool = false,
+    result::Int = 1,
+)
+    status = termination_status(model)
+    ret =
+        (status == OPTIMAL) ||
+        (allow_local && (status == LOCALLY_SOLVED)) ||
+        (allow_almost && (status == ALMOST_OPTIMAL)) ||
+        (allow_almost && allow_local && (status == ALMOST_LOCALLY_SOLVED))
+    if ret
+        primal = primal_status(model; result)
+        ret &=
+            (primal == FEASIBLE_POINT) ||
+            (allow_almost && (primal == NEARLY_FEASIBLE_POINT))
+    end
+    if ret && dual
+        dual_stat = dual_status(model; result)
+        ret &=
+            (dual_stat == FEASIBLE_POINT) ||
+            (allow_almost && (dual_stat == NEARLY_FEASIBLE_POINT))
+    end
+    return ret
+end
+
 """
     solve_time(model::GenericModel)
 
diff --git a/test/test_model.jl b/test/test_model.jl
index 20bcb345d6e..c9134eab2ad 100644
--- a/test/test_model.jl
+++ b/test/test_model.jl
@@ -1244,4 +1244,73 @@ function test_caching_mps_model()
     return
 end
 
+function test_is_solved_and_feasible()
+    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
+    model = direct_model(mock)
+    for term in [
+        MOI.OPTIMAL,
+        MOI.LOCALLY_SOLVED,
+        MOI.ALMOST_OPTIMAL,
+        MOI.ALMOST_LOCALLY_SOLVED,
+        MOI.TIME_LIMIT,
+    ]
+        _global = term == MOI.OPTIMAL
+        has_local = _global || (term == MOI.LOCALLY_SOLVED)
+        _almost_global = _global || (term == MOI.ALMOST_OPTIMAL)
+        _almost_local =
+            has_local || _almost_global || (term == MOI.ALMOST_LOCALLY_SOLVED)
+        for primal in
+            [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT, MOI.NO_SOLUTION]
+            _primal = primal == MOI.FEASIBLE_POINT
+            _almost_primal = _primal || primal == MOI.NEARLY_FEASIBLE_POINT
+            for dual in
+                [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT, MOI.NO_SOLUTION]
+                _dual = dual == MOI.FEASIBLE_POINT
+                _almost_dual = _dual || dual == MOI.NEARLY_FEASIBLE_POINT
+                MOI.set(mock, MOI.TerminationStatus(), term)
+                MOI.set(mock, MOI.PrimalStatus(), primal)
+                MOI.set(mock, MOI.DualStatus(), dual)
+                @test is_solved_and_feasible(model) == (has_local && _primal)
+                @test is_solved_and_feasible(model; dual = true) ==
+                      (has_local && _primal && _dual)
+                @test is_solved_and_feasible(model; allow_local = false) ==
+                      (_global && _primal)
+                @test is_solved_and_feasible(
+                    model;
+                    dual = true,
+                    allow_local = false,
+                ) == (_global && _primal && _dual)
+                @test is_solved_and_feasible(model; allow_almost = true) ==
+                      (_almost_local && _almost_primal)
+                @test is_solved_and_feasible(
+                    model;
+                    dual = true,
+                    allow_almost = true,
+                ) == (_almost_local && _almost_primal && _almost_dual)
+                @test is_solved_and_feasible(
+                    model;
+                    allow_local = false,
+                    allow_almost = true,
+                ) == (_almost_global && _almost_primal)
+                @test is_solved_and_feasible(
+                    model;
+                    dual = true,
+                    allow_local = false,
+                    allow_almost = true,
+                ) == (_almost_global && _almost_primal && _almost_dual)
+                MOI.set(mock, MOI.ResultCount(), 3)
+                MOI.set(mock, MOI.PrimalStatus(3), primal)
+                MOI.set(mock, MOI.DualStatus(3), dual)
+                @test !is_solved_and_feasible(model; result = 2)
+                @test !is_solved_and_feasible(model; dual = true, result = 2)
+                @test is_solved_and_feasible(model; result = 3) ==
+                      (has_local && _primal)
+                @test is_solved_and_feasible(model; dual = true, result = 3) ==
+                      (has_local && _primal && _dual)
+            end
+        end
+    end
+    return
+end
+
 end  # module TestModels