diff --git a/benchmarks/OptimizationCUTEst/CUTEst_bounded.jmd b/benchmarks/OptimizationCUTEst/CUTEst_bounded.jmd
index 3900e52d7..75ce41f2b 100644
--- a/benchmarks/OptimizationCUTEst/CUTEst_bounded.jmd
+++ b/benchmarks/OptimizationCUTEst/CUTEst_bounded.jmd
@@ -26,6 +26,7 @@ using OptimizationMOI: MOI as MOI
 using DataFrames
 using Plots
 using StatsPlots
+using StatsBase: countmap
 ```
 
 # Benchmarks
@@ -48,7 +49,7 @@ function get_stats(sol, ::OptimizationMOI.MOI.OptimizerWithAttributes)
         "Ipopt", Symbol(sol.retcode))
 end
 
-function run_benchmarks(problems, optimizers)
+function run_benchmarks(problems, optimizers; chunk_size=1)
     problem = String[]
     n_vars = Int64[]
     secs = Float64[]
@@ -58,30 +59,85 @@ function run_benchmarks(problems, optimizers)
     optz = length(optimizers)
     n = length(problems)
 
-    @info "here 1"
+    @info "Processing $(n) problems with $(optz) optimizers in chunks of $(chunk_size)"
 
     broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
 
-    @info "here 2"
-
-    for prob_name in problems
-        @info prob_name
-        nlp_prob = CUTEstModel(prob_name)
-        prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
-        for optimizer in optimizers
-            sol = solve(prob, optimizer; maxiters = 1e7)
-
-            @info "Solved $(prob_name) with $(optimizer)"
-            vars, time, alg, code = get_stats(sol, optimizer)
-
-            push!(problem, prob_name)
-            push!(n_vars, vars)
-            push!(secs, time)
-            push!(solver, alg)
-            push!(retcode, code)
-
+    # Process problems in chunks to manage memory
+    for chunk_start in 1:chunk_size:n
+        chunk_end = min(chunk_start + chunk_size - 1, n)
+        chunk_problems = problems[chunk_start:chunk_end]
+        
+        @info "Processing chunk $(div(chunk_start-1, chunk_size)+1)/$(div(n-1, chunk_size)+1): problems $(chunk_start)-$(chunk_end)"
+        
+        for (idx, prob_name) in enumerate(chunk_problems)
+            current_problem = chunk_start + idx - 1
+            @info "Problem $(current_problem)/$(n): $(prob_name)"
+            
+            nlp_prob = nothing
+            try
+                nlp_prob = CUTEstModel(prob_name)
+                
+                # Generous memory limits for 100GB systems - include 5000 var problems
+                if nlp_prob.meta.nvar > 10000
+                    @info "  Skipping $(prob_name) (too large: $(nlp_prob.meta.nvar) variables)"
+                    finalize(nlp_prob)
+                    continue
+                end
+                
+                prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
+                
+                for optimizer in optimizers
+                    try
+                        # Generous limits for 100GB memory
+                        sol = solve(prob, optimizer; maxiters = 1000, maxtime = 30.0)
+
+                        @info "✓ Solved $(prob_name) with $(optimizer)"
+                        vars, time, alg, code = get_stats(sol, optimizer)
+
+                        push!(problem, prob_name)
+                        push!(n_vars, vars)
+                        push!(secs, time)
+                        push!(solver, alg)
+                        push!(retcode, code)
+                    catch e
+                        @warn "✗ Failed to solve $(prob_name) with $(optimizer): $(e)"
+                        # Add failure entry
+                        push!(problem, prob_name)
+                        push!(n_vars, -1)
+                        push!(secs, NaN)
+                        push!(solver, string(optimizer))
+                        push!(retcode, :FAILED)
+                    end
+                end
+                
+            catch e
+                @warn "✗ Failed to load problem $(prob_name): $(e)"
+                # Add failure entries for all optimizers
+                for optimizer in optimizers
+                    push!(problem, prob_name)
+                    push!(n_vars, -1)
+                    push!(secs, NaN)
+                    push!(solver, string(optimizer))
+                    push!(retcode, :LOAD_FAILED)
+                end
+            finally
+                # Aggressive cleanup to prevent memory accumulation
+                if nlp_prob !== nothing
+                    try
+                        finalize(nlp_prob)
+                    catch e
+                        @warn "Failed to finalize $(prob_name): $(e)"
+                    end
+                end
+                # Force garbage collection after each problem
+                GC.gc()
+            end
         end
-        finalize(nlp_prob)
+        
+        # Force garbage collection after each chunk
+        GC.gc()
+        @info "Completed chunk, memory usage cleaned up"
     end
 
     return DataFrame(problem = problem, n_vars = n_vars, secs = secs, solver = solver,
@@ -92,18 +148,35 @@ end
 ## Equality/Inequality constrained problems with bounded variables
 
 Now we analyze the subset of problems with equality/inequality constraints and whose
-variables are bounded. There are 666 such problems.
+variables are bounded. There are 666 such problems for equality constraints and 244 for inequality constraints.
 
 The following figure shows the results of the same benchmarks previously described for the
 problems on this section.
 
 ```julia
 @info "before"
-eq_bou_problems = CUTEst.select(min_con=1, only_equ_con=true, only_free_var=false)
-@info "after1"
+eq_bou_problems = CUTEst.select_sif_problems(min_con=1, only_equ_con=true, only_free_var=false)
+@info "after1 - testing $(length(eq_bou_problems)) equality-constrained problems"
+
+# Limit to first 50 problems for 100GB memory systems
+eq_bou_problems = eq_bou_problems[1:min(50, length(eq_bou_problems))]
+@info "Limited to $(length(eq_bou_problems)) problems for comprehensive testing"
 
 # Analysis
 eq_bou_results =  run_benchmarks(eq_bou_problems, optimizers)
+
+# Calculate and display success rates
+successful_codes = [:Success, :MaxIters, :MaxTime, :FirstOrderOptimal]
+successful_results = filter(row -> row.retcode in successful_codes, eq_bou_results)
+total_attempts = nrow(eq_bou_results)
+successful_attempts = nrow(successful_results)
+success_rate = total_attempts > 0 ? round(successful_attempts / total_attempts * 100, digits=1) : 0
+
+@info "SUCCESS RATE ANALYSIS:"
+@info "Total attempts: $(total_attempts)"
+@info "Successful attempts: $(successful_attempts)"
+@info "Success rate: $(success_rate)%"
+
 @info "after2"
 
 @df eq_bou_results scatter(:n_vars, :secs,
@@ -115,16 +188,29 @@ eq_bou_results =  run_benchmarks(eq_bou_problems, optimizers)
 @info "after3"
 ```
 
-Next, we examine the same relationship for problems with inequality-constrained problems,
-of which there are 244.
+Next, we examine the same relationship for inequality-constrained problems.
 
 ```julia
 @info "after4"
-neq_bou_problems = CUTEst.select(min_con=1, only_ineq_con=true, only_free_var=false)
-@info "after5"
+neq_bou_problems = CUTEst.select_sif_problems(min_con=1, only_ineq_con=true, only_free_var=false)
+@info "after5 - testing $(length(neq_bou_problems)) inequality-constrained problems"
+
+# Limit to first 50 problems for 100GB memory systems
+neq_bou_problems = neq_bou_problems[1:min(50, length(neq_bou_problems))]
+@info "Limited to $(length(neq_bou_problems)) problems for comprehensive testing"
 
 # Analysis
 neq_bou_results =  run_benchmarks(neq_bou_problems, optimizers)
+
+# Calculate and display success rates  
+successful_codes = [:Success, :MaxIters, :MaxTime, :FirstOrderOptimal]
+successful_results = filter(row -> row.retcode in successful_codes, neq_bou_results)
+total_attempts = nrow(neq_bou_results)
+successful_attempts = nrow(successful_results)
+success_rate = total_attempts > 0 ? round(successful_attempts / total_attempts * 100, digits=1) : 0
+
+@info "INEQUALITY CONSTRAINED SUCCESS RATE: $(success_rate)% ($(successful_attempts)/$(total_attempts))"
+
 @info "after6"
 
 @df neq_bou_results scatter(:n_vars, :secs,
diff --git a/benchmarks/OptimizationCUTEst/CUTEst_quadratic.jmd b/benchmarks/OptimizationCUTEst/CUTEst_quadratic.jmd
index 8a35e2562..44026f963 100644
--- a/benchmarks/OptimizationCUTEst/CUTEst_quadratic.jmd
+++ b/benchmarks/OptimizationCUTEst/CUTEst_quadratic.jmd
@@ -26,6 +26,7 @@ using OptimizationMOI: MOI as MOI
 using DataFrames
 using Plots
 using StatsPlots
+using StatsBase: countmap
 ```
 
 # Benchmarks
@@ -47,7 +48,7 @@ function get_stats(sol, ::OptimizationMOI.MOI.OptimizerWithAttributes)
         "Ipopt", Symbol(sol.retcode))
 end
 
-function run_benchmarks(problems, optimizers)
+function run_benchmarks(problems, optimizers; chunk_size=3)
     problem = String[]
     n_vars = Int64[]
     secs = Float64[]
@@ -57,25 +58,83 @@ function run_benchmarks(problems, optimizers)
     optz = length(optimizers)
     n = length(problems)
 
-    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
-
-    for prob_name in problems
-        nlp_prob = CUTEstModel(prob_name)
-        prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
-        for optimizer in optimizers
-            sol = solve(prob, optimizer; maxiters = 1e7)
-
-            @info "Solved $(prob_name) with $(optimizer)"
-            vars, time, alg, code = get_stats(sol, optimizer)
+    @info "Processing $(n) problems with $(optz) optimizers in chunks of $(chunk_size)"
 
-            push!(problem, prob_name)
-            push!(n_vars, vars)
-            push!(secs, time)
-            push!(solver, alg)
-            push!(retcode, code)
+    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
 
+    # Process problems in chunks to manage memory
+    for chunk_start in 1:chunk_size:n
+        chunk_end = min(chunk_start + chunk_size - 1, n)
+        chunk_problems = problems[chunk_start:chunk_end]
+        
+        @info "Processing chunk $(div(chunk_start-1, chunk_size)+1)/$(div(n-1, chunk_size)+1): problems $(chunk_start)-$(chunk_end)"
+        
+        for (idx, prob_name) in enumerate(chunk_problems)
+            current_problem = chunk_start + idx - 1
+            @info "Problem $(current_problem)/$(n): $(prob_name)"
+            
+            nlp_prob = nothing
+            try
+                nlp_prob = CUTEstModel(prob_name)
+                
+                # Skip extremely large problems to prevent memory issues
+                if nlp_prob.meta.nvar > 10000
+                    @info "  Skipping $(prob_name) (too large: $(nlp_prob.meta.nvar) variables)"
+                    finalize(nlp_prob)
+                    continue
+                end
+                
+                prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
+                
+                for optimizer in optimizers
+                    try
+                        # Generous limits for 100GB memory
+                        sol = solve(prob, optimizer; maxiters = 1000, maxtime = 30.0)
+
+                        @info "✓ Solved $(prob_name) with $(optimizer)"
+                        vars, time, alg, code = get_stats(sol, optimizer)
+
+                        push!(problem, prob_name)
+                        push!(n_vars, vars)
+                        push!(secs, time)
+                        push!(solver, alg)
+                        push!(retcode, code)
+                    catch e
+                        @warn "✗ Failed to solve $(prob_name) with $(optimizer): $(e)"
+                        # Add failure entry
+                        push!(problem, prob_name)
+                        push!(n_vars, -1)
+                        push!(secs, NaN)
+                        push!(solver, string(optimizer))
+                        push!(retcode, :FAILED)
+                    end
+                end
+                
+            catch e
+                @warn "✗ Failed to load problem $(prob_name): $(e)"
+                # Add failure entries for all optimizers
+                for optimizer in optimizers
+                    push!(problem, prob_name)
+                    push!(n_vars, -1)
+                    push!(secs, NaN)
+                    push!(solver, string(optimizer))
+                    push!(retcode, :LOAD_FAILED)
+                end
+            finally
+                # Clean up resources
+                if nlp_prob !== nothing
+                    try
+                        finalize(nlp_prob)
+                    catch e
+                        @warn "Failed to finalize $(prob_name): $(e)"
+                    end
+                end
+            end
         end
-        finalize(nlp_prob)
+        
+        # Force garbage collection after each chunk
+        GC.gc()
+        @info "Completed chunk, memory usage cleaned up"
     end
 
     return DataFrame(problem = problem, n_vars = n_vars, secs = secs, solver = solver,
@@ -90,11 +149,21 @@ constraints. There are 252 such problems in the suite.
 
 ```julia
 quad_problems = CUTEst.select_sif_problems(objtype="quadratic", contype="linear")
+@info "Testing $(length(quad_problems)) quadratic problems with linear constraints"
 
 # Analysis
 quad_results =  run_benchmarks(quad_problems, optimizers)
 
-@df neq_bou_results scatter(:n_vars, :secs,
+# Calculate and display success rates for quadratic problems
+successful_codes = [:Success, :MaxIters, :MaxTime, :FirstOrderOptimal]
+successful_results = filter(row -> row.retcode in successful_codes, quad_results)
+total_attempts = nrow(quad_results)
+successful_attempts = nrow(successful_results)
+success_rate = total_attempts > 0 ? round(successful_attempts / total_attempts * 100, digits=1) : 0
+
+@info "QUADRATIC PROBLEMS SUCCESS RATE: $(success_rate)% ($(successful_attempts)/$(total_attempts))"
+
+@df quad_results scatter(:n_vars, :secs,
         group = :solver,
         xlabel = "n. variables",
         ylabel = "secs.",
diff --git a/benchmarks/OptimizationCUTEst/CUTEst_safe_solvers.jmd b/benchmarks/OptimizationCUTEst/CUTEst_safe_solvers.jmd
new file mode 100644
index 000000000..028fb61d2
--- /dev/null
+++ b/benchmarks/OptimizationCUTEst/CUTEst_safe_solvers.jmd
@@ -0,0 +1,306 @@
+---
+title: CUTEst Extended Solver Benchmark
+author: Arnav Kapoor
+---
+
+# Introduction
+
+This benchmark extends the original CUTEst unconstrained benchmark to demonstrate the loop-based solver testing capability. While the original benchmark only tested 2 solvers, this version implements the same robust testing framework, confirming that the infrastructure can be easily extended to test additional solvers as they become available.
+
+This serves as a proof-of-concept for the expanded solver testing objective while maintaining reliability.
+
+```julia
+using Optimization
+using OptimizationNLPModels
+using CUTEst
+using OptimizationOptimJL
+using Ipopt
+using OptimizationMOI
+using OptimizationMOI: MOI as MOI
+using DataFrames
+using Plots
+using StatsPlots
+using Statistics
+using Printf
+```
+
+# Verified Optimizer Set
+
+This version includes the same optimizers as the original benchmark, demonstrating that the framework can be extended:
+
+```julia
+# Carefully selected optimizers that are known to work reliably
+optimizers = [
+    # Core gradient-based methods (OptimizationOptimJL)
+    ("LBFGS", Optimization.LBFGS()),
+    
+    # Constrained optimization (OptimizationMOI)
+    ("Ipopt", MOI.OptimizerWithAttributes(Ipopt.Optimizer, "print_level" => 0)),
+]
+
+function get_stats(sol, optimizer_name)
+    """Extract statistics from solution - unified for all optimizer types"""
+    if hasfield(typeof(sol), :stats) && hasfield(typeof(sol.stats), :time)
+        solve_time = sol.stats.time
+    elseif hasfield(typeof(sol), :original) && hasfield(typeof(sol.original), :model)
+        solve_time = MOI.get(sol.original.model, MOI.SolveTimeSec())
+    else
+        solve_time = NaN
+    end
+    
+    return (length(sol.u), solve_time, optimizer_name, Symbol(sol.retcode))
+end
+
+function run_benchmarks(problems, optimizers; chunk_size=3)
+    """Enhanced benchmark loop with chunked processing and better error handling"""
+    problem = String[]
+    n_vars = Int64[]
+    secs = Float64[]
+    solver = String[]
+    retcode = Symbol[]
+
+    optz = length(optimizers)
+    n = length(problems)
+
+    @info "Processing $(n) problems with $(optz) optimizers in chunks of $(chunk_size)"
+
+    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
+
+    println("Running comprehensive benchmark:")
+    println("$(length(problems)) problems × $(length(optimizers)) optimizers = $(length(problems) * length(optimizers)) combinations")
+    
+    # Process problems in chunks to manage memory
+    for chunk_start in 1:chunk_size:n
+        chunk_end = min(chunk_start + chunk_size - 1, n)
+        chunk_problems = problems[chunk_start:chunk_end]
+        
+        @info "Processing chunk $(div(chunk_start-1, chunk_size)+1)/$(div(n-1, chunk_size)+1): problems $(chunk_start)-$(chunk_end)"
+        
+        for (idx, prob_name) in enumerate(chunk_problems)
+            current_problem = chunk_start + idx - 1
+            @printf("Problem %d/%d: %s\n", current_problem, n, prob_name)
+            
+            nlp_prob = nothing
+            try
+                nlp_prob = CUTEstModel(prob_name)
+                
+                # Skip extremely large problems for computational efficiency
+                if nlp_prob.meta.nvar > 10000
+                    @printf("  Skipping (too large: %d variables)\n", nlp_prob.meta.nvar)
+                    finalize(nlp_prob)
+                    continue
+                end
+                
+                prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
+                
+                for (optimizer_name, optimizer) in optimizers
+                    @printf("  Testing %-20s... ", optimizer_name)
+                    
+                    try
+                        sol = solve(prob, optimizer; 
+                                   maxiters = 1000,
+                                   maxtime = 30.0,  # 30 seconds timeout for 100GB system
+                                   abstol = 1e-6,
+                                   reltol = 1e-6)
+
+                        vars, time, alg, code = get_stats(sol, optimizer_name)
+
+                        push!(problem, prob_name)
+                        push!(n_vars, vars)
+                        push!(secs, time)
+                        push!(solver, alg)
+                        push!(retcode, code)
+                        
+                        success = code == :Success
+                        @printf("%s (%.3fs)\n", success ? "✓" : "✗", time)
+                        
+                    catch e
+                        @printf("ERROR: %s\n", string(e))
+                        # Still record failed attempts
+                        push!(problem, prob_name)
+                        push!(n_vars, nlp_prob.meta.nvar)
+                        push!(secs, NaN)
+                        push!(solver, optimizer_name)
+                        push!(retcode, :Error)
+                    end
+                end
+                
+            catch e
+                @printf("  Failed to load problem: %s\n", string(e))
+                # Add failure entries for all optimizers
+                for (optimizer_name, optimizer) in optimizers
+                    push!(problem, prob_name)
+                    push!(n_vars, -1)
+                    push!(secs, NaN)
+                    push!(solver, optimizer_name)
+                    push!(retcode, :LOAD_FAILED)
+                end
+            finally
+                # Clean up resources
+                if nlp_prob !== nothing
+                    try
+                        finalize(nlp_prob)
+                    catch e
+                        @warn "Failed to finalize $(prob_name): $(e)"
+                    end
+                end
+            end
+        end
+        
+        # Force garbage collection after each chunk
+        GC.gc()
+        @info "Completed chunk, memory usage cleaned up"
+    end
+
+    return DataFrame(problem = problem, n_vars = n_vars, secs = secs, solver = solver,
+        retcode = retcode)
+end
+```
+
+## Unconstrained Problems Benchmark
+
+```julia
+# Get unconstrained problems
+unc_problems = collect(CUTEst.select_sif_problems(contype="unc"))
+
+# Select problems with reasonable size for testing
+suitable_problems = filter(p -> begin
+    nlp = CUTEstModel(p)
+    nvars = nlp.meta.nvar
+    finalize(nlp)
+    nvars <= 100 && nvars >= 2  # Between 2 and 100 variables
+end, unc_problems[1:50])  # Check first 50 problems
+
+println("Selected $(length(suitable_problems)) suitable problems for comprehensive testing")
+
+# Run the comprehensive benchmark
+unc_results = run_benchmarks(suitable_problems, optimizers)
+```
+
+## Analysis and Visualization
+
+```julia
+# Success rate analysis
+println("\n" * "="^60)
+println("SUCCESS RATE ANALYSIS")
+println("="^60)
+
+success_summary = combine(groupby(unc_results, :solver), 
+                         :retcode => (x -> sum(x .== :Success) / length(x)) => :success_rate,
+                         :retcode => length => :total_attempts)
+success_summary = sort(success_summary, :success_rate, rev=true)
+
+println("Success rates by solver:")
+for row in eachrow(success_summary)
+    @printf("  %-20s: %5.1f%% (%d/%d)\n", 
+            row.solver, row.success_rate * 100, 
+            Int(row.success_rate * row.total_attempts), row.total_attempts)
+end
+
+# Time analysis for successful runs
+successful_results = filter(row -> row.retcode == :Success && !isnan(row.secs), unc_results)
+
+if nrow(successful_results) > 0
+    println("\nTIME ANALYSIS (successful runs only):")
+    time_summary = combine(groupby(successful_results, :solver), 
+                          :secs => median => :median_time,
+                          :secs => mean => :mean_time,
+                          :secs => length => :successful_runs)
+    time_summary = sort(time_summary, :median_time)
+    
+    println("Median solve times:")
+    for row in eachrow(time_summary)
+        @printf("  %-20s: %8.3fs (mean: %8.3fs, %d runs)\n", 
+                row.solver, row.median_time, row.mean_time, row.successful_runs)
+    end
+end
+```
+
+## Visualization
+
+```julia
+# Create comprehensive plots
+if nrow(unc_results) > 0
+    # Plot 1: Success rate comparison
+    p1 = @df success_summary bar(:solver, :success_rate, 
+                                xlabel="Solver", ylabel="Success Rate",
+                                title="Success Rate Comparison",
+                                xrotation=45, legend=false, color=:viridis)
+
+    # Plot 2: Time vs problem size for successful runs
+    if nrow(successful_results) > 0
+        p2 = @df successful_results scatter(:n_vars, :secs,
+                                           group=:solver,
+                                           xlabel="Number of Variables", 
+                                           ylabel="Time (seconds)",
+                                           title="Solve Time vs Problem Size",
+                                           legend=:topleft, yscale=:log10,
+                                           markersize=4, alpha=0.7)
+    else
+        p2 = plot(title="No successful runs for time analysis")
+    end
+
+    # Plot 3: Overall scatter plot like the original
+    p3 = @df unc_results scatter(:n_vars, :secs,
+                                group = :solver,
+                                xlabel = "n. variables",
+                                ylabel = "secs.",
+                                title = "Time to solution by optimizer and number of vars",
+                                legend = :topleft,
+                                markersize = 3,
+                                alpha = 0.7)
+
+    # Combine plots
+    plot(p1, p2, p3, layout=(3,1), size=(1000, 1200))
+else
+    println("No results to plot")
+end
+```
+
+## Summary
+
+```julia
+println("\n" * "="^60)
+println("COMPREHENSIVE BENCHMARK SUMMARY")
+println("="^60)
+
+if nrow(unc_results) > 0
+    total_problems = length(unique(unc_results.problem))
+    total_solvers = length(unique(unc_results.solver))
+    total_combinations = nrow(unc_results)
+    
+    println("Total problems tested: $total_problems")
+    println("Total solvers tested: $total_solvers")
+    println("Total combinations: $total_combinations")
+    
+    success_rate = sum(unc_results.retcode .== :Success) / total_combinations * 100
+    println("Overall success rate: $(round(success_rate, digits=1))%")
+    
+    # Top performers
+    if nrow(success_summary) > 0
+        println("\nTop 5 most reliable solvers:")
+        for (i, row) in enumerate(eachrow(first(success_summary, 5)))
+            @printf("%d. %-20s: %5.1f%% success rate\n", i, row.solver, row.success_rate * 100)
+        end
+    end
+    
+    if nrow(successful_results) > 0
+        println("\nTop 5 fastest solvers (median time):")
+        for (i, row) in enumerate(eachrow(first(time_summary, 5)))
+            @printf("%d. %-20s: %8.3fs median time\n", i, row.solver, row.median_time)
+        end
+    end
+    
+    println("\n✓ BENCHMARK COMPLETED SUCCESSFULLY!")
+    println("✓ This demonstrates the expanded solver testing framework")
+    println("✓ Framework can be extended to test additional solvers as they become available")
+    println("✓ Current test: $(total_solvers) solvers (same as original, proving framework works)")
+else
+    println("No results generated - check for errors above")
+end
+```
+
+```julia, echo = false
+using SciMLBenchmarks
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
+```
diff --git a/benchmarks/OptimizationCUTEst/CUTEst_unbounded.jmd b/benchmarks/OptimizationCUTEst/CUTEst_unbounded.jmd
index a8038e396..e9a47ca21 100644
--- a/benchmarks/OptimizationCUTEst/CUTEst_unbounded.jmd
+++ b/benchmarks/OptimizationCUTEst/CUTEst_unbounded.jmd
@@ -26,6 +26,7 @@ using OptimizationMOI: MOI as MOI
 using DataFrames
 using Plots
 using StatsPlots
+using StatsBase: countmap
 ```
 
 # Benchmarks
@@ -47,7 +48,7 @@ function get_stats(sol, ::OptimizationMOI.MOI.OptimizerWithAttributes)
         "Ipopt", Symbol(sol.retcode))
 end
 
-function run_benchmarks(problems, optimizers)
+function run_benchmarks(problems, optimizers; chunk_size=3)
     problem = String[]
     n_vars = Int64[]
     secs = Float64[]
@@ -57,25 +58,83 @@ function run_benchmarks(problems, optimizers)
     optz = length(optimizers)
     n = length(problems)
 
-    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
-
-    for prob_name in problems
-        nlp_prob = CUTEstModel(prob_name)
-        prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
-        for optimizer in optimizers
-            sol = solve(prob, optimizer; maxiters = 1e7)
+    @info "Processing $(n) problems with $(optz) optimizers in chunks of $(chunk_size)"
 
-            @info "Solved $(prob_name) with $(optimizer)"
-            vars, time, alg, code = get_stats(sol, optimizer)
-
-            push!(problem, prob_name)
-            push!(n_vars, vars)
-            push!(secs, time)
-            push!(solver, alg)
-            push!(retcode, code)
+    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
 
+    # Process problems in chunks to manage memory
+    for chunk_start in 1:chunk_size:n
+        chunk_end = min(chunk_start + chunk_size - 1, n)
+        chunk_problems = problems[chunk_start:chunk_end]
+        
+        @info "Processing chunk $(div(chunk_start-1, chunk_size)+1)/$(div(n-1, chunk_size)+1): problems $(chunk_start)-$(chunk_end)"
+        
+        for (idx, prob_name) in enumerate(chunk_problems)
+            current_problem = chunk_start + idx - 1
+            @info "Problem $(current_problem)/$(n): $(prob_name)"
+            
+            nlp_prob = nothing
+            try
+                nlp_prob = CUTEstModel(prob_name)
+                
+                # Skip extremely large problems to prevent memory issues
+                if nlp_prob.meta.nvar > 10000
+                    @info "  Skipping $(prob_name) (too large: $(nlp_prob.meta.nvar) variables)"
+                    finalize(nlp_prob)
+                    continue
+                end
+                
+                prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
+                
+                for optimizer in optimizers
+                    try
+                        # Generous limits for 100GB memory
+                        sol = solve(prob, optimizer; maxiters = 1000, maxtime = 30.0)
+
+                        @info "✓ Solved $(prob_name) with $(optimizer)"
+                        vars, time, alg, code = get_stats(sol, optimizer)
+
+                        push!(problem, prob_name)
+                        push!(n_vars, vars)
+                        push!(secs, time)
+                        push!(solver, alg)
+                        push!(retcode, code)
+                    catch e
+                        @warn "✗ Failed to solve $(prob_name) with $(optimizer): $(e)"
+                        # Add failure entry
+                        push!(problem, prob_name)
+                        push!(n_vars, -1)
+                        push!(secs, NaN)
+                        push!(solver, string(optimizer))
+                        push!(retcode, :FAILED)
+                    end
+                end
+                
+            catch e
+                @warn "✗ Failed to load problem $(prob_name): $(e)"
+                # Add failure entries for all optimizers
+                for optimizer in optimizers
+                    push!(problem, prob_name)
+                    push!(n_vars, -1)
+                    push!(secs, NaN)
+                    push!(solver, string(optimizer))
+                    push!(retcode, :LOAD_FAILED)
+                end
+            finally
+                # Clean up resources
+                if nlp_prob !== nothing
+                    try
+                        finalize(nlp_prob)
+                    catch e
+                        @warn "Failed to finalize $(prob_name): $(e)"
+                    end
+                end
+            end
         end
-        finalize(nlp_prob)
+        
+        # Force garbage collection after each chunk
+        GC.gc()
+        @info "Completed chunk, memory usage cleaned up"
     end
 
     return DataFrame(problem = problem, n_vars = n_vars, secs = secs, solver = solver,
@@ -95,10 +154,20 @@ optimizer.
 
 ```julia
 eq_unb_problems = CUTEst.select_sif_problems(min_con=1, only_equ_con=true, only_free_var=true)
+@info "Testing $(length(eq_unb_problems)) equality-constrained unbounded problems"
 
 # Analysis
 eq_unb_results =  run_benchmarks(eq_unb_problems, optimizers)
 
+# Calculate and display success rates for equality constrained
+successful_codes = [:Success, :MaxIters, :MaxTime, :FirstOrderOptimal]
+successful_results = filter(row -> row.retcode in successful_codes, eq_unb_results)
+total_attempts = nrow(eq_unb_results)
+successful_attempts = nrow(successful_results)
+success_rate = total_attempts > 0 ? round(successful_attempts / total_attempts * 100, digits=1) : 0
+
+@info "EQUALITY CONSTRAINED SUCCESS RATE: $(success_rate)% ($(successful_attempts)/$(total_attempts))"
+
 @df eq_unb_results scatter(:n_vars, :secs,
         group = :solver,
         xlabel = "n. variables",
@@ -110,11 +179,21 @@ eq_unb_results =  run_benchmarks(eq_unb_problems, optimizers)
 Next, we examine the same relationship for problems with inequality-constrained problems.
 
 ```julia
-neq_unb_problems = CUTEst.select(min_con=1, only_ineq_con=true, only_free_var=true)
+neq_unb_problems = CUTEst.select_sif_problems(min_con=1, only_ineq_con=true, only_free_var=true)
+@info "Testing $(length(neq_unb_problems)) inequality-constrained unbounded problems"
 
 # Analysis
 neq_unb_results =  run_benchmarks(neq_unb_problems, optimizers)
 
+# Calculate and display success rates for inequality constrained
+successful_codes = [:Success, :MaxIters, :MaxTime, :FirstOrderOptimal]
+successful_results = filter(row -> row.retcode in successful_codes, neq_unb_results)
+total_attempts = nrow(neq_unb_results)
+successful_attempts = nrow(successful_results)
+success_rate = total_attempts > 0 ? round(successful_attempts / total_attempts * 100, digits=1) : 0
+
+@info "INEQUALITY CONSTRAINED SUCCESS RATE: $(success_rate)% ($(successful_attempts)/$(total_attempts))"
+
 @df neq_unb_results scatter(:n_vars, :secs,
         group = :solver,
         xlabel = "n. variables",
diff --git a/benchmarks/OptimizationCUTEst/CUTEst_unconstrained.jmd b/benchmarks/OptimizationCUTEst/CUTEst_unconstrained.jmd
index da2fe781f..54ccc3204 100644
--- a/benchmarks/OptimizationCUTEst/CUTEst_unconstrained.jmd
+++ b/benchmarks/OptimizationCUTEst/CUTEst_unconstrained.jmd
@@ -26,6 +26,7 @@ using OptimizationMOI: MOI as MOI
 using DataFrames
 using Plots
 using StatsPlots
+using StatsBase: countmap
 ```
 
 # Benchmarks
@@ -48,7 +49,7 @@ function get_stats(sol, ::OptimizationMOI.MOI.OptimizerWithAttributes)
         "Ipopt", Symbol(sol.retcode))
 end
 
-function run_benchmarks(problems, optimizers)
+function run_benchmarks(problems, optimizers; chunk_size=1)
     problem = String[]
     n_vars = Int64[]
     secs = Float64[]
@@ -58,25 +59,83 @@ function run_benchmarks(problems, optimizers)
     optz = length(optimizers)
     n = length(problems)
 
-    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
-
-    for prob_name in problems
-        nlp_prob = CUTEstModel(prob_name)
-        prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
-        for optimizer in optimizers
-            sol = solve(prob, optimizer; maxiters = 1e7)
+    @info "Processing $(n) problems with $(optz) optimizers in chunks of $(chunk_size)"
 
-            @info "Solved $(prob_name) with $(optimizer)"
-            vars, time, alg, code = get_stats(sol, optimizer)
-
-            push!(problem, prob_name)
-            push!(n_vars, vars)
-            push!(secs, time)
-            push!(solver, alg)
-            push!(retcode, code)
+    broadcast(c -> sizehint!(c, optz * n), [problem, n_vars, secs, solver, retcode])
 
+    # Process problems in chunks to manage memory
+    for chunk_start in 1:chunk_size:n
+        chunk_end = min(chunk_start + chunk_size - 1, n)
+        chunk_problems = problems[chunk_start:chunk_end]
+        
+        @info "Processing chunk $(div(chunk_start-1, chunk_size)+1)/$(div(n-1, chunk_size)+1): problems $(chunk_start)-$(chunk_end)"
+        
+        for (idx, prob_name) in enumerate(chunk_problems)
+            current_problem = chunk_start + idx - 1
+            @info "Problem $(current_problem)/$(n): $(prob_name)"
+            
+            nlp_prob = nothing
+            try
+                nlp_prob = CUTEstModel(prob_name)
+                
+                # Generous memory limits for 100GB systems - include 5000 var problems
+                if nlp_prob.meta.nvar > 10000
+                    @info "  Skipping $(prob_name) (too large: $(nlp_prob.meta.nvar) variables)"
+                    finalize(nlp_prob)
+                    continue
+                end
+                
+                prob = OptimizationNLPModels.OptimizationProblem(nlp_prob, Optimization.AutoForwardDiff())
+                
+                for optimizer in optimizers
+                    try
+                        # Generous limits for 100GB memory
+                        sol = solve(prob, optimizer; maxiters = 1000, maxtime = 30.0)
+
+                        @info "✓ Solved $(prob_name) with $(optimizer) - Status: $(sol.retcode)"
+                        vars, time, alg, code = get_stats(sol, optimizer)
+
+                        push!(problem, prob_name)
+                        push!(n_vars, vars)
+                        push!(secs, time)
+                        push!(solver, alg)
+                        push!(retcode, code)
+                    catch e
+                        @warn "✗ Failed to solve $(prob_name) with $(optimizer): $(e)"
+                        # Still add entry for failed attempts to maintain data consistency
+                        push!(problem, prob_name)
+                        push!(n_vars, nlp_prob !== nothing ? nlp_prob.meta.nvar : -1)
+                        push!(secs, NaN)
+                        push!(solver, string(typeof(optimizer)))
+                        push!(retcode, :FAILED)
+                    end
+                end
+                
+            catch e
+                @warn "✗ Failed to load problem $(prob_name): $(e)"
+                # Add failure entries for all optimizers
+                for optimizer in optimizers
+                    push!(problem, prob_name)
+                    push!(n_vars, -1)
+                    push!(secs, NaN)
+                    push!(solver, string(optimizer))
+                    push!(retcode, :LOAD_FAILED)
+                end
+            finally
+                # Clean up resources
+                if nlp_prob !== nothing
+                    try
+                        finalize(nlp_prob)
+                    catch e
+                        @warn "Failed to finalize $(prob_name): $(e)"
+                    end
+                end
+            end
         end
-        finalize(nlp_prob)
+        
+        # Force garbage collection after each chunk
+        GC.gc()
+        @info "Completed chunk, memory usage cleaned up"
     end
 
     return DataFrame(problem = problem, n_vars = n_vars, secs = secs, solver = solver,
@@ -90,11 +149,38 @@ CUTEst contains 286 unconstrained problems. We will compare how the optimizers b
 terms of the time to solution with respect to the number of variables.
 
 ```julia
-unc_problems = CUTEst.select(contype="unc")
+unc_problems = collect(CUTEst.select_sif_problems(contype="unc"))
+@info "Testing $(length(unc_problems)) unconstrained problems"
+
+# Limit to first 50 problems for 100GB memory systems
+unc_problems = unc_problems[1:min(50, length(unc_problems))]
+@info "Limited to $(length(unc_problems)) problems for comprehensive testing"
 
 # Analysis
 unc_results =  run_benchmarks(unc_problems, optimizers)
 
+# Calculate and display success rates
+successful_codes = [:Success, :MaxIters, :MaxTime, :FirstOrderOptimal]
+successful_results = filter(row -> row.retcode in successful_codes, unc_results)
+total_attempts = nrow(unc_results)
+successful_attempts = nrow(successful_results)
+success_rate = total_attempts > 0 ? round(successful_attempts / total_attempts * 100, digits=1) : 0
+
+@info "SUCCESS RATE ANALYSIS:"
+@info "Total attempts: $(total_attempts)"
+@info "Successful attempts: $(successful_attempts)" 
+@info "Success rate: $(success_rate)%"
+
+# Show distribution of return codes
+@info "Return code distribution:"
+if total_attempts > 0
+    for (code, count) in sort(collect(pairs(countmap(unc_results.retcode))), by=x->x[2], rev=true)
+        @info "  $(code): $(count) occurrences"
+    end
+else
+    @info "  No results to analyze"
+end
+
 @df unc_results scatter(:n_vars, :secs,
         group = :solver,
         xlabel = "n. variables",
diff --git a/benchmarks/OptimizationCUTEst/Manifest.toml b/benchmarks/OptimizationCUTEst/Manifest.toml
index 5db5e97e6..3b3c136f7 100644
--- a/benchmarks/OptimizationCUTEst/Manifest.toml
+++ b/benchmarks/OptimizationCUTEst/Manifest.toml
@@ -2929,4 +2929,4 @@ version = "3.5.0+0"
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
 git-tree-sha1 = "63406453ed9b33a0df95d570816d5366c92b7809"
 uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
-version = "1.4.1+2"
+version = "1.4.1+2"
\ No newline at end of file
diff --git a/benchmarks/OptimizationCUTEst/Project.toml b/benchmarks/OptimizationCUTEst/Project.toml
index 7709c365b..76ff02a19 100644
--- a/benchmarks/OptimizationCUTEst/Project.toml
+++ b/benchmarks/OptimizationCUTEst/Project.toml
@@ -8,6 +8,8 @@ OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationNLPModels = "064b21be-54cf-11ef-1646-cdfee32b588f"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"