diff --git a/README.md b/README.md
index b35894ded..d160d35ea 100644
--- a/README.md
+++ b/README.md
@@ -51,8 +51,7 @@ Our goal is controlling the first output $y_1$, but the second one $y_2$ should
 35:
 
 ```julia
-mhe = MovingHorizonEstimator(model)
-mpc = LinMPC(mhe, Mwt=[1, 0], Nwt=[0.1])
+mpc = LinMPC(model, Mwt=[1, 0], Nwt=[0.1])
 mpc = setconstraint!(mpc, ymax=[Inf, 35])
 ```
 
diff --git a/src/ModelPredictiveControl.jl b/src/ModelPredictiveControl.jl
index 423475c5f..8bb73005f 100644
--- a/src/ModelPredictiveControl.jl
+++ b/src/ModelPredictiveControl.jl
@@ -14,6 +14,8 @@ using DifferentiationInterface: hessian!, value_gradient_and_hessian!, prepare_h
 using DifferentiationInterface: Constant, Cache
 using SparseConnectivityTracer: TracerSparsityDetector
 using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern
+using SparseMatrixColorings: NaturalOrder, LargestFirst, SmallestLast
+using SparseMatrixColorings: IncidenceDegree, DynamicLargestFirst
 
 import ProgressLogging
 
diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl
index e622e9965..7ac5ddc75 100644
--- a/src/controller/nonlinmpc.jl
+++ b/src/controller/nonlinmpc.jl
@@ -5,7 +5,7 @@ const DEFAULT_NONLINMPC_JACDENSE  = AutoForwardDiff()
 const DEFAULT_NONLINMPC_JACSPARSE = AutoSparse(
     AutoForwardDiff();
     sparsity_detector=TracerSparsityDetector(),
-    coloring_algorithm=GreedyColoringAlgorithm(),
+    coloring_algorithm=GreedyColoringAlgorithm(ALL_COLORING_ORDERS),
 )
 const DEFAULT_NONLINMPC_HESSIAN = DEFAULT_NONLINMPC_JACSPARSE
 
@@ -291,10 +291,17 @@ NonLinMPC controller with a sample time Ts = 10.0 s:
     AutoSparse(
         AutoForwardDiff(); 
         sparsity_detector  = TracerSparsityDetector(), 
-        coloring_algorithm = GreedyColoringAlgorithm()
+        coloring_algorithm = GreedyColoringAlgorithm((
+            NaturalOrder(),
+            LargestFirst(),
+            SmallestLast(),
+            IncidenceDegree(),
+            DynamicLargestFirst()
+        ))
     )
     ```
-    This is also the sparse backend selected for the Hessian of the Lagrangian function if 
+    that is, it will test many coloring orders at preparation and keep the best. This is
+    also the sparse backend selected for the Hessian of the Lagrangian function if 
     `oracle=true` and `hessian=true`, which is the second exception. Second order 
     derivatives are only supported with `oracle=true` option.
     
diff --git a/src/estimator/mhe/construct.jl b/src/estimator/mhe/construct.jl
index 6069b6ea4..5d7c9120a 100644
--- a/src/estimator/mhe/construct.jl
+++ b/src/estimator/mhe/construct.jl
@@ -5,7 +5,7 @@ const DEFAULT_NONLINMHE_JACOBIAN  = AutoForwardDiff()
 const DEFAULT_NONLINMHE_HESSIAN   = AutoSparse(
     AutoForwardDiff();
     sparsity_detector=TracerSparsityDetector(),
-    coloring_algorithm=GreedyColoringAlgorithm(),
+    coloring_algorithm=GreedyColoringAlgorithm(ALL_COLORING_ORDERS),
 ) 
 
 @doc raw"""
@@ -382,9 +382,17 @@ MovingHorizonEstimator estimator with a sample time Ts = 10.0 s:
     AutoSparse(
         AutoForwardDiff(); 
         sparsity_detector  = TracerSparsityDetector(), 
-        coloring_algorithm = GreedyColoringAlgorithm()
+        coloring_algorithm = GreedyColoringAlgorithm((
+            NaturalOrder(),
+            LargestFirst(),
+            SmallestLast(),
+            IncidenceDegree(),
+            DynamicLargestFirst()
+        ))
     )
     ```
+    that is, it will test many coloring orders at preparation and keep the best. 
+    
     The slack variable ``ε`` relaxes the constraints if enabled, see [`setconstraint!`](@ref). 
     It is disabled by default for the MHE (from `Cwt=Inf`) but it should be activated for
     problems with two or more types of bounds, to ensure feasibility (e.g. on the estimated
diff --git a/src/general.jl b/src/general.jl
index 3972d6c32..4ce717392 100644
--- a/src/general.jl
+++ b/src/general.jl
@@ -6,6 +6,15 @@ const DEFAULT_LWT = 0.0
 const DEFAULT_CWT = 1e5
 const DEFAULT_EWT = 0.0
 
+"All deterministic algorithms for matrix coloring order in `SparseMatrixColoring.jl`."
+const ALL_COLORING_ORDERS = (
+    NaturalOrder(),
+    LargestFirst(),
+    SmallestLast(),
+    IncidenceDegree(),
+    DynamicLargestFirst(),
+)
+
 "Termination status that means 'no solution available'."
 const ERROR_STATUSES = (
     JuMP.INFEASIBLE, JuMP.DUAL_INFEASIBLE, JuMP.LOCALLY_INFEASIBLE,