diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
deleted file mode 100644
index 144c91b1f..000000000
--- a/.github/workflows/CI.yml
+++ /dev/null
@@ -1,49 +0,0 @@
-name: CI
-on:
-  push:
-    branches:
-      - main
-    tags: ['*']
-  pull_request:
-concurrency:
-  # Skip intermediate builds: always.
-  # Cancel intermediate builds: only if it is a pull request build.
-  group: ${{ github.workflow }}-${{ github.ref }}
-  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
-jobs:
-  test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
-    runs-on: ${{ matrix.os }}
-    permissions: # needed for julia-actions/cache delete old caches that it has created
-      actions: write
-      contents: read
-    continue-on-error: ${{ matrix.version == 'pre' }}
-    strategy:
-      fail-fast: false
-      matrix:
-        version:
-          - 'lts' # long-term support release
-          - '1'   # latest stable 1.x release
-          - 'pre' # latest stable prerelease
-          # - 'nightly' # commented since noisy + 'pre' allows testing upcoming versions
-        os:
-          - ubuntu-latest
-        arch:
-          - x64
-    steps:
-      - name: Set JULIA_DEBUG environment variable if applicable
-        if: ${{ runner.debug == '1' }}
-        run: echo "JULIA_DEBUG=ModelPredictiveControl" >> $GITHUB_ENV
-      - uses: actions/checkout@v7
-      - uses: julia-actions/setup-julia@v3
-        with:
-          version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v3
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-runtest@v1
-      - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v7
-        with:
-          token: ${{ secrets.CODECOV_TOKEN }}
-          fail_ci_if_error: false 
\ No newline at end of file
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
deleted file mode 100644
index 3042569cb..000000000
--- a/.github/workflows/TagBot.yml
+++ /dev/null
@@ -1,33 +0,0 @@
-name: TagBot
-on:
-  issue_comment:
-    types:
-      - created
-  workflow_dispatch:
-    inputs:
-      lookback:
-        default: 3
-permissions:
-  actions: read
-  checks: read
-  contents: write
-  deployments: read
-  issues: read
-  discussions: read
-  packages: read
-  pages: read
-  pull-requests: read
-  repository-projects: read
-  security-events: read
-  statuses: read
-jobs:
-  TagBot:
-    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
-    runs-on: ubuntu-latest
-    steps:
-      - uses: JuliaRegistries/TagBot@v1
-        with:
-          token: ${{ secrets.GITHUB_TOKEN }}
-          # Edit the following line to reflect the actual name of the GitHub Secret containing your private key
-          ssh: ${{ secrets.DOCUMENTER_KEY }}
-          # ssh: ${{ secrets.NAME_OF_MY_SSH_PRIVATE_KEY_SECRET }}
\ No newline at end of file
diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index 49c57a06e..821c0aad6 100644
--- a/.github/workflows/benchmark.yml
+++ b/.github/workflows/benchmark.yml
@@ -12,4 +12,5 @@ jobs:
     steps:
       - uses: MilesCranmer/AirspeedVelocity.jl@action-v1
         with:
-          julia-version: '1'
\ No newline at end of file
+          julia-version: '1'
+          job-summary: 'true'
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
deleted file mode 100644
index 70a059074..000000000
--- a/.github/workflows/documentation.yml
+++ /dev/null
@@ -1,26 +0,0 @@
-name: Documentation
-
-on:
-  push:
-    branches:
-      - main # update to match your development branch (master, main, dev, trunk, ...)
-    tags: '*'
-  pull_request:
-
-jobs:
-  build:
-    permissions:
-      contents: write
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v7
-      - uses: julia-actions/setup-julia@v3
-        with:
-          version: '1'
-      - name: Install dependencies
-        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
-      - name: Build and deploy
-        env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
\ No newline at end of file
diff --git a/.github/workflows/juliaci.yml b/.github/workflows/juliaci.yml
new file mode 100644
index 000000000..ea09a233e
--- /dev/null
+++ b/.github/workflows/juliaci.yml
@@ -0,0 +1,25 @@
+name: Julia CI
+
+on:
+  push: {branches: [main,master]}
+  pull_request: {types: [opened,synchronize,reopened,ready_for_review,converted_to_draft]}
+  issue_comment: {types: [created]}
+  workflow_dispatch:
+    inputs:
+      feature:
+        type: choice
+        description: What to run
+        options: [DocDeploy, LintAndTest, TagBot]
+
+jobs:
+  julia-ci:
+    uses: julia-testitems/testitem-workflow/.github/workflows/juliaci.yml@v2
+    with:
+      testitem-timeout: 3600
+      include-rc-versions: true
+      include-windows-x86: false
+      include-linux-x86: false
+      include-macos-x64: false
+    permissions: write-all
+    secrets:
+      codecov_token: ${{ secrets.CODECOV_TOKEN }}
\ No newline at end of file
diff --git a/src/estimator/construct.jl b/src/estimator/construct.jl
index 328593021..c7f187787 100644
--- a/src/estimator/construct.jl
+++ b/src/estimator/construct.jl
@@ -359,11 +359,8 @@ and unstable `model` (see Examples).
 ```jldoctest
 julia> model = LinModel(append(tf(3, [10, 1]), tf(2, [1, 0]), tf(4,[-5, 1])), 1.0);
 
-julia> nint_ym = default_nint(model)
-3-element Vector{Int64}:
- 1
- 0
- 1
+julia> default_nint(model) |> println
+[1, 0, 1]
 ```
 """
 function default_nint(model::LinModel, i_ym=1:model.ny, nint_u=0)
diff --git a/src/sim_model.jl b/src/sim_model.jl
index 15ed1595e..ffdb9f6e7 100644
--- a/src/sim_model.jl
+++ b/src/sim_model.jl
@@ -289,7 +289,7 @@ function savetime!(model::SimModel)
     return model.t[]
 end
 
-"""
+@doc raw"""
     periodsleep(model::SimModel, busywait=false) -> nothing
 
 Sleep for `model.Ts` s minus the time elapsed since the last call to [`savetime!`](@ref).
@@ -305,14 +305,14 @@ simple soft real-time simulations, see the example below.
     restart it at a convenient time e.g.: just before calling `periodsleep`.
 
 # Examples
-```jldoctest
+```jldoctest; filter = r"(\d+\.\d)\d*" => s"\1"
 julia> model = LinModel(tf(2, [0.3, 1]), 0.25);
 
 julia> function sim_realtime!(model)
            t_0 = time()
            for i=1:3
                t = savetime!(model)      # first function called
-               println(round(t - t_0, digits=3))
+               println(round(t - t_0, digits=2))
                updatestate!(model, [1])
                periodsleep(model, true)  # last function called
            end
diff --git a/test/1_test_sim_model.jl b/test/1_test_sim_model.jl
index e2eabb750..8a618b864 100644
--- a/test/1_test_sim_model.jl
+++ b/test/1_test_sim_model.jl
@@ -144,7 +144,7 @@ end
         updatestate!(linmodel1, [1])
         periodsleep(linmodel1)
     end
-    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.025))
+    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.05))
     linmodel2 = LinModel(tf(2, [0.1, 1]), 0.25)
     times2 = zeros(5)
     for i=1:5
@@ -152,7 +152,7 @@ end
         updatestate!(linmodel2, [1])
         periodsleep(linmodel2, true)
     end
-    @test all(isapprox.(diff(times2[2:end]), 0.25, atol=0.0001))
+    @test all(isapprox.(diff(times2[2:end]), 0.25, atol=0.05))
 end
 
 @testitem "NonLinModel construction" setup=[SetupMPCtests] begin
@@ -406,7 +406,7 @@ end
         updatestate!(nonlinmodel1, [1])
         periodsleep(nonlinmodel1)
     end
-    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.025))
+    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.05))
     linmodel2 = LinModel(tf(2, [0.1, 1]), 0.25)
     nonlinmodel2 = NonLinModel(
         (x,u,_,_)->linmodel2.A*x + linmodel2.Bu*u,
@@ -419,5 +419,5 @@ end
         updatestate!(nonlinmodel2, [1])
         periodsleep(nonlinmodel2, true)
     end
-    @test all(isapprox.(diff(times2[2:end]), 0.25, atol=0.0001))
+    @test all(isapprox.(diff(times2[2:end]), 0.25, atol=0.05))
 end
\ No newline at end of file
diff --git a/test/2_test_state_estim.jl b/test/2_test_state_estim.jl
index 3c98ebc2c..b920e0d06 100644
--- a/test/2_test_state_estim.jl
+++ b/test/2_test_state_estim.jl
@@ -149,7 +149,7 @@ end
         updatestate!(kalmanfilter1, [1], [1])
         periodsleep(kalmanfilter1, true)
     end
-    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.01))
+    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.05))
 end
     
 @testitem "KF construction" setup=[SetupMPCtests] begin
@@ -1031,7 +1031,7 @@ end
     )
 end
 
-@testitem "MHE estimation and getinfo (LinModel)" setup=[SetupMPCtests] begin
+@testitem "MHE estim. & getinfo (LinModel, SS)" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, ForwardDiff
     using JuMP, DAQP
     linmodel = LinModel(sys,Ts,i_u=[1,2], i_d=[3])
@@ -1123,6 +1123,14 @@ end
         updatestate!(mhe5, [10, 50], [50, NaN], [5])
     )
     
+end
+
+@testitem "MHE estim. & getinfo (LinModel, MS)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, ForwardDiff
+    using JuMP, DAQP
+    linmodel = LinModel(sys,Ts,i_u=[1,2], i_d=[3])
+    linmodel = setop!(linmodel, uop=[10,50], yop=[50,30], dop=[5])
+
     mhe5 = MovingHorizonEstimator(linmodel, He=2, direct=true, transcription=MultipleShooting())
     for i in 1:40
         preparestate!(mhe5, [51, 32], [5])
@@ -1141,7 +1149,7 @@ end
 
 end
 
-@testitem "MHE estimation and getinfo (NonLinModel)" setup=[SetupMPCtests] begin
+@testitem "MHE estim. & getinfo (NonLinModel, SS)" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, ForwardDiff
     using JuMP, Ipopt, DifferentiationInterface, SparseMatrixColorings, SparseConnectivityTracer
     import ForwardDiff
@@ -1229,6 +1237,60 @@ end
     preparestate!(mhe2, [50, 30], [5])
     @test mhe2([5]) ≈ [50, 30] atol=5e-3
 
+    Q̂ = diagm([1/4, 1/4, 1/4, 1/4].^2) 
+    R̂ = diagm([1, 1].^2)
+    optim = Model(Ipopt.Optimizer)
+    covestim = ExtendedKalmanFilter(nonlinmodel, 1:2, 0, 0, Q̂, Q̂, R̂)
+    mhe5 = MovingHorizonEstimator(nonlinmodel, 1, 1:2, 0, 0, Q̂, Q̂, R̂; optim, covestim)
+    preparestate!(mhe5, [50, 30], [5])
+    x̂ = updatestate!(mhe5, [10, 50], [50, 30], [5])
+    @test x̂ ≈ zeros(4) atol=1e-9
+    @test mhe5.x̂0 ≈ zeros(4) atol=1e-9
+    preparestate!(mhe5, [50, 30], [5])
+    @test evaloutput(mhe5, [5]) ≈ mhe5([5]) ≈ [50, 30]
+    info = getinfo(mhe5)
+    @test info[:x̂] ≈ x̂ atol=1e-9
+    @test info[:Ŷ][end-1:end] ≈ [50, 30] atol=1e-9  
+
+    # coverage of the branch with error termination status (with an infeasible problem):
+    mhe_infeas = MovingHorizonEstimator(nonlinmodel, He=1, Cwt=Inf)
+    mhe_infeas = setconstraint!(mhe_infeas, v̂min=[1, 1], v̂max=[-1, -1]) 
+    @test_logs(
+        (:error, "MHE terminated without solution: estimation in open-loop "*
+                 "(more info in debug log)"), 
+        preparestate!(mhe_infeas, [0, 0], [0])
+    )
+
+    @test_nowarn ModelPredictiveControl.info2debugstr(info)
+    @test_throws ErrorException setstate!(mhe1, [1,2,3,4,5,6], diagm(.1:.1:.6))
+
+    hessian = AutoSparse(
+        AutoForwardDiff();
+        sparsity_detector=TracerSparsityDetector(),
+        coloring_algorithm=GreedyColoringAlgorithm(),
+    )
+
+    mhe8 = MovingHorizonEstimator(nonlinmodel; He=2, hessian)
+    @test_logs(
+        (:warn, "NaN values in the MHE measurements ym: ignoring them in the objective"),
+        preparestate!(mhe8, [50, NaN], [5])
+    )
+    @test mhe8.x̂0 ≈ zeros(6) atol=1e-9
+
+end
+
+@testitem "MHE estim. & getinfo (NonLinModel, MS)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, ForwardDiff
+    using JuMP, Ipopt, DifferentiationInterface, SparseMatrixColorings, SparseConnectivityTracer
+    import ForwardDiff
+
+    linmodel = LinModel(sys,Ts,i_u=[1,2], i_d=[3])
+    linmodel = setop!(linmodel, uop=[10,50], yop=[50,30], dop=[5])
+    f(x,u,d,model) = model.A*x + model.Bu*u + model.Bd*d
+    h(x,d,model)   = model.C*x + model.Dd*d
+    nonlinmodel = NonLinModel(f, h, Ts, 2, 4, 2, 1, solver=nothing, p=linmodel)
+    nonlinmodel = setop!(nonlinmodel, uop=[10,50], yop=[50,30], dop=[5])
+
     mhe3 = MovingHorizonEstimator(
         nonlinmodel, He=3, direct=false, transcription=MultipleShooting(f_threads=true)
     )
@@ -1244,22 +1306,13 @@ end
     updatestate!(mhe3, [10, 50], [50, 30], [5])
     info = getinfo(mhe3) # test getinfo when Nk<He 
     @test info[:Ŵ] ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-    @test info[:V̂] ≈ [0.0, 0.0]    
+    @test info[:V̂] ≈ [0.0, 0.0]
+end
 
-    Q̂ = diagm([1/4, 1/4, 1/4, 1/4].^2) 
-    R̂ = diagm([1, 1].^2)
-    optim = Model(Ipopt.Optimizer)
-    covestim = ExtendedKalmanFilter(nonlinmodel, 1:2, 0, 0, Q̂, Q̂, R̂)
-    mhe5 = MovingHorizonEstimator(nonlinmodel, 1, 1:2, 0, 0, Q̂, Q̂, R̂; optim, covestim)
-    preparestate!(mhe5, [50, 30], [5])
-    x̂ = updatestate!(mhe5, [10, 50], [50, 30], [5])
-    @test x̂ ≈ zeros(4) atol=1e-9
-    @test mhe5.x̂0 ≈ zeros(4) atol=1e-9
-    preparestate!(mhe5, [50, 30], [5])
-    @test evaloutput(mhe5, [5]) ≈ mhe5([5]) ≈ [50, 30]
-    info = getinfo(mhe5)
-    @test info[:x̂] ≈ x̂ atol=1e-9
-    @test info[:Ŷ][end-1:end] ≈ [50, 30] atol=1e-9
+@testitem "MHE estim. & getinfo (NonLinModel, TC)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, ForwardDiff
+    using JuMP, Ipopt, DifferentiationInterface, SparseMatrixColorings, SparseConnectivityTracer
+    import ForwardDiff
 
     f! = (ẋ,x,u,_,_) -> ẋ .= -0.001x .+ u 
     h! = (y,x,_,_)   -> y .= x 
@@ -1292,6 +1345,16 @@ end
     info = getinfo(mhe7) # test getinfo when Nk<He 
     @test info[:Ŵ] ≈ [0.0, 0.0]
     @test info[:V̂] ≈ [0.0]
+end
+
+@testitem "MHE estim. & getinfo (NonLinModel, OC)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, ForwardDiff
+    using JuMP, Ipopt, DifferentiationInterface, SparseMatrixColorings, SparseConnectivityTracer
+    import ForwardDiff
+
+    f! = (ẋ,x,u,_,_) -> ẋ .= -0.001x .+ u 
+    h! = (y,x,_,_)   -> y .= x 
+    nonlinmodel_c = NonLinModel(f!, h!, 500, 1, 1, 1)
 
     transcription = OrthogonalCollocation(f_threads=true, h_threads=true)
     mhe10 = MovingHorizonEstimator(
@@ -1313,41 +1376,10 @@ end
         updatestate!(mhe11, [-6], [13])
     end
     preparestate!(mhe11, [13])
-    @test mhe11() ≈ [13] atol=5e-3  
-
-    # coverage of the branch with error termination status (with an infeasible problem):
-    mhe_infeas = MovingHorizonEstimator(nonlinmodel, He=1, Cwt=Inf)
-    mhe_infeas = setconstraint!(mhe_infeas, v̂min=[1, 1], v̂max=[-1, -1]) 
-    @test_logs(
-        (:error, "MHE terminated without solution: estimation in open-loop "*
-                 "(more info in debug log)"), 
-        preparestate!(mhe_infeas, [0, 0], [0])
-    )
-
-    # for coverage of NLP functions, the univariate syntax of JuMP.@operator
-    mhe7 = MovingHorizonEstimator(nonlinmodel, He=1, Cwt=Inf)
-    setconstraint!(mhe7, v̂min=[-51,-52], v̂max=[53,54])
-    x̂ = preparestate!(mhe7, [50, 30], [5])
-    @test x̂ ≈ zeros(6) atol=1e-9
-    @test_nowarn ModelPredictiveControl.info2debugstr(info)
-    @test_throws ErrorException setstate!(mhe1, [1,2,3,4,5,6], diagm(.1:.1:.6))
-
-    hessian = AutoSparse(
-        AutoForwardDiff();
-        sparsity_detector=TracerSparsityDetector(),
-        coloring_algorithm=GreedyColoringAlgorithm(),
-    )
-
-    mhe8 = MovingHorizonEstimator(nonlinmodel; He=2, hessian)
-    @test_logs(
-        (:warn, "NaN values in the MHE measurements ym: ignoring them in the objective"),
-        preparestate!(mhe8, [50, NaN], [5])
-    )
-    @test mhe8.x̂0 ≈ zeros(6) atol=1e-9
-
+    @test mhe11() ≈ [13] atol=5e-3
 end
 
-@testitem "MHE estimation with unfilled window" setup=[SetupMPCtests] begin
+@testitem "MHE estim. with unfilled window" setup=[SetupMPCtests] begin
     f(x,u,_,_) = 0.5x + u
     h(x,_,_) = x
     model = NonLinModel(f, h, 10.0, 1, 1, 1, solver=nothing)
@@ -1373,7 +1405,7 @@ end
     
 end
 
-@testitem "MHE fallbacks for arrival covariance estimation" setup=[SetupMPCtests] begin
+@testitem "MHE fallbacks for arrival covariance" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
     linmodel = setop!(LinModel(sys,Ts,i_u=[1,2], i_d=[3]), uop=[10,50], yop=[50,30], dop=[5])
     f(x,u,d,model) = model.A*x + model.Bu*u + model.Bd*d
@@ -1592,8 +1624,10 @@ end
     linmodel2 = LinModel(sys, Ts, i_u=[1,2], i_d=[3])
     linmodel2 = setop!(linmodel2, uop=[10,50], yop=[50,30], dop=[5])
     function gclv!(LHS, X̂e, _, _, _, _, _, _, _, nx̂, _ )
-        for i in 1:div(length(X̂e), nx̂)
-            LHS[(i-1)+1] = 0.5 - X̂e[(i-1)*nx̂ + 1]  # First state >= 0.5
+        nc_Nk = length(X̂e) ÷ nx̂
+        LHS .= 0 # for when Nk < He, fill unused values with 0s
+        for i in 1:nc_Nk
+            LHS[i] = 0.5 - X̂e[(i-1)*nx̂ + 1]  # First state >= 0.5
         end
         return nothing
     end
@@ -1678,8 +1712,10 @@ end
     nonlinmodel2 = NonLinModel(f, h, Ts, 2, 4, 2, 1, solver=nothing, p=linmodel2)
     nonlinmodel2 = setop!(nonlinmodel2, uop=[10,50], yop=[50,30], dop=[5])
     function gcnlv!(LHS, X̂e, _, _, _, _, _, _, _, nx̂, _)
-        for i in 1:div(length(X̂e), nx̂)
-            LHS[(i-1)+1] = 0.5 - X̂e[(i-1)*nx̂ + 1]  # First state >= 0.5
+        nc_Nk = length(X̂e) ÷ nx̂
+        LHS .= 0 # for when Nk < He, fill unused values with 0s
+        for i in 1:nc_Nk
+            LHS[i] = 0.5 - X̂e[(i-1)*nx̂ + 1]  # First state >= 0.5
         end
         return nothing
     end
@@ -1800,7 +1836,7 @@ end
         updatestate!(mhe, [11, 50], y, [25])
         updatestate!(kf,  [11, 50], y, [25])
     end
-    @test X̂_mhe ≈ X̂_kf atol=1e-6 rtol=1e-6
+    @test X̂_mhe ≈ X̂_kf atol=1e-3 rtol=1e-3
     kf  = KalmanFilter(linmodel, nint_ym=0, direct=true)
     # recuperate P̂(-1|-1) exact value using the Kalman filter:
     preparestate!(kf, [50, 30], [20])
@@ -1818,7 +1854,7 @@ end
         updatestate!(mhe, [11, 50], y, [25])
         updatestate!(kf,  [11, 50], y, [25])
     end
-    @test X̂_mhe ≈ X̂_kf atol=1e-6 rtol=1e-6
+    @test X̂_mhe ≈ X̂_kf atol=1e-3 rtol=1e-3
 
     f = (x,u,d,model) -> model.A*x + model.Bu*u + model.Bd*d
     h = (x,d,model)   -> model.C*x + model.Dd*d
@@ -1842,8 +1878,8 @@ end
         updatestate!(ukf, [11, 50], y, [25])
         updatestate!(ekf, [11, 50], y, [25])
     end
-    @test X̂_mhe ≈ X̂_ukf atol=1e-6 rtol=1e-6
-    @test X̂_mhe ≈ X̂_ekf atol=1e-6 rtol=1e-6
+    @test X̂_mhe ≈ X̂_ukf atol=1e-3 rtol=1e-3
+    @test X̂_mhe ≈ X̂_ekf atol=1e-3 rtol=1e-3
     
     ukf = UnscentedKalmanFilter(nonlinmodel, nint_ym=0, direct=true)
     ekf = ExtendedKalmanFilter(nonlinmodel, nint_ym=0, direct=true)
@@ -1869,8 +1905,8 @@ end
         updatestate!(ukf, [11, 50], y, [25])
         updatestate!(ekf, [11, 50], y, [25])
     end
-    @test X̂_mhe ≈ X̂_ukf atol=1e-6 rtol=1e-6
-    @test X̂_mhe ≈ X̂_ekf atol=1e-6 rtol=1e-6 
+    @test X̂_mhe ≈ X̂_ukf atol=1e-3 rtol=1e-3
+    @test X̂_mhe ≈ X̂_ekf atol=1e-3 rtol=1e-3 
 end
 
 @testitem "MHE LinModel v.s. NonLinModel" setup=[SetupMPCtests] begin
diff --git a/test/3_test_predictive_control.jl b/test/3_test_predictive_control.jl
index 6f4f636e9..83ba626e7 100644
--- a/test/3_test_predictive_control.jl
+++ b/test/3_test_predictive_control.jl
@@ -89,7 +89,7 @@
     @test_throws DimensionMismatch LinMPC(model, Wy=ones(2, model.ny), Wu=ones(3, model.nu))
 end
 
-@testitem "LinMPC moves and getinfo" setup=[SetupMPCtests] begin
+@testitem "LinMPC moves & getinfo" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
     linmodel = setop!(LinModel(tf(5, [2, 1]), 3), yop=[10])
@@ -581,7 +581,7 @@ end
 
 @testitem "LinMPC real-time simulations" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
-    linmodel1 = LinModel(tf(2, [10, 1]), 0.1)
+    linmodel1 = LinModel(tf(2, [10, 1]), 0.25)
     mpc1 = LinMPC(linmodel1)
     times1 = zeros(5)
     for i=1:5
@@ -590,7 +590,7 @@ end
         updatestate!(mpc1, [1], [1])
         periodsleep(mpc1, true)
     end
-    @test all(isapprox.(diff(times1[2:end]), 0.1, atol=0.01))
+    @test all(isapprox.(diff(times1[2:end]), 0.25, atol=0.05))
 end
 
 @testitem "ExplicitMPC construction" setup=[SetupMPCtests] begin
@@ -637,7 +637,7 @@ end
 
 end
 
-@testitem "ExplicitMPC moves and getinfo" setup=[SetupMPCtests] begin
+@testitem "ExplicitMPC moves & getinfo" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
     model = LinModel(tf(5, [2, 1]), 3)
     mpc1 = ExplicitMPC(model, Nwt=[0], Hp=1000, Hc=1)
@@ -905,7 +905,7 @@ end
     @test_throws ArgumentError OrthogonalCollocation(2)
 end
 
-@testitem "NonLinMPC moves and getinfo (LinModel)" setup=[SetupMPCtests] begin
+@testitem "NonLinMPC moves & getinfo (LinModel, SS)" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
     using DifferentiationInterface
     import FiniteDiff
@@ -948,14 +948,6 @@ end
     preparestate!(nmpc6, [0])
     @test moveinput!(nmpc6, [0]) ≈ [0.0] atol=5e-2
     
-    nmpc9 = NonLinMPC(linmodel, Nwt=[0], Hp=100, Hc=1, transcription=MultipleShooting())
-    preparestate!(nmpc9, [10])
-    u = moveinput!(nmpc9, [20])
-    @test u ≈ [2] atol=5e-2
-    info = getinfo(nmpc9)
-    @test info[:u] ≈ u
-    @test info[:Ŷ][end] ≈ 20 atol=5e-2
-    
     # coverage of the branch with error termination status (with an infeasible problem):
     nmpc_infeas = NonLinMPC(linmodel, Hp=1, Hc=1, Cwt=Inf)
     nmpc_infeas = setconstraint!(nmpc_infeas, umin=[+1], umax=[-1])
@@ -969,7 +961,22 @@ end
     @test_nowarn ModelPredictiveControl.info2debugstr(info)
 end
 
-@testitem "NonLinMPC moves and getinfo (NonLinModel)" setup=[SetupMPCtests] begin
+@testitem "NonLinMPC moves & getinfo (LinModel, MS)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
+    using DifferentiationInterface
+    import FiniteDiff
+    linmodel = setop!(LinModel(tf(5, [2000, 1]), 3000.0), yop=[10])
+    
+    nmpc9 = NonLinMPC(linmodel, Nwt=[0], Hp=100, Hc=1, transcription=MultipleShooting())
+    preparestate!(nmpc9, [10])
+    u = moveinput!(nmpc9, [20])
+    @test u ≈ [2] atol=5e-2
+    info = getinfo(nmpc9)
+    @test info[:u] ≈ u
+    @test info[:Ŷ][end] ≈ 20 atol=5e-2
+end
+
+@testitem "NonLinMPC moves & getinfo (NonLinModel, SS)" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
     using DifferentiationInterface
     import FiniteDiff
@@ -977,10 +984,6 @@ end
     f = (x,u,d,model) -> model.A*x + model.Bu*u + model.Bd*d
     h = (x,d,model)   -> model.C*x + model.Dd*d
     nonlinmodel = NonLinModel(f, h, 3000.0, 1, 2, 1, 1, solver=nothing, p=linmodel2)
-
-    f! = (ẋ,x,u,_,_) -> ẋ .= -0.001x .+ u 
-    h! = (y,x,_,_) -> y .= x 
-    nonlinmodel_c = NonLinModel(f!, h!, 500, 1, 1, 1)
     
     nmpc1 = NonLinMPC(nonlinmodel, Nwt=[0], Hp=100, Hc=1)
     preparestate!(nmpc1, [0], [0])
@@ -1004,6 +1007,43 @@ end
     u = moveinput!(nmpc3, [0], d, R̂u=fill(12, nmpc3.Hp))
     @test u ≈ [12] atol=5e-2
 
+    nonlinmodel2 = NonLinModel{Float32}(f, h, 3000.0, 1, 2, 1, 1, solver=nothing, p=linmodel2)
+    nmpc7  = NonLinMPC(nonlinmodel2, Hp=10)
+    y = similar(nonlinmodel2.yop)
+    ModelPredictiveControl.h!(y, nonlinmodel2, Float32[0,0], Float32[0], nonlinmodel2.p)
+    preparestate!(nmpc7, [0], [0])
+    @test moveinput!(nmpc7, [0], [0]) ≈ [0.0] atol=5e-2
+
+    nmpc10 = setconstraint!(NonLinMPC(
+        nonlinmodel, Nwt=[0], Hp=100, Hc=1, 
+        gradient=AutoFiniteDiff(),
+        jacobian=AutoFiniteDiff(),
+        hessian=AutoFiniteDiff()
+        ), ymax=[100], ymin=[-100]
+    )
+    preparestate!(nmpc10, [0], [0])
+    u = moveinput!(nmpc10, [10], [0])
+    @test u ≈ [2] atol=5e-2
+    info = getinfo(nmpc10)
+    @test info[:u] ≈ u
+    @test info[:Ŷ][end] ≈ 10 atol=5e-2
+    
+    nmpc11 = NonLinMPC(nonlinmodel, Hp=10, Hc=[1, 2, 3, 4], Nwt=[10])
+    preparestate!(nmpc11, y, [0])
+    moveinput!(nmpc11, [10], [0])
+    ΔU_diff = diff(getinfo(nmpc11)[:U])
+    @test ΔU_diff[[2, 4, 5, 7, 8, 9]] ≈ zeros(6) atol=1e-9
+end
+
+@testitem "NonLinMPC moves & getinfo (NonLinModel, MS)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
+    using DifferentiationInterface
+    import FiniteDiff
+    linmodel2 = LinModel([tf(5, [2000, 1]) tf(7, [8000,1])], 3000.0, i_d=[2])
+    f = (x,u,d,model) -> model.A*x + model.Bu*u + model.Bd*d
+    h = (x,d,model)   -> model.C*x + model.Dd*d
+    nonlinmodel = NonLinModel(f, h, 3000.0, 1, 2, 1, 1, solver=nothing, p=linmodel2)
+
     transcription = MultipleShooting()
     nmpc4 = NonLinMPC(nonlinmodel; Nwt=[0], Hp=100, Hc=1, transcription)
     preparestate!(nmpc4, [0], [0])
@@ -1022,6 +1062,20 @@ end
     info = getinfo(nmpc4t)
     @test info[:u] ≈ u
     @test info[:Ŷ][end] ≈ 10 atol=5e-2
+end
+
+@testitem "NonLinMPC moves & getinfo (NonLinModel, TC)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
+    using DifferentiationInterface
+    import FiniteDiff
+    linmodel2 = LinModel([tf(5, [2000, 1]) tf(7, [8000,1])], 3000.0, i_d=[2])
+    f = (x,u,d,model) -> model.A*x + model.Bu*u + model.Bd*d
+    h = (x,d,model)   -> model.C*x + model.Dd*d
+    nonlinmodel = NonLinModel(f, h, 3000.0, 1, 2, 1, 1, solver=nothing, p=linmodel2)
+
+    f! = (ẋ,x,u,_,_) -> ẋ .= -0.001x .+ u 
+    h! = (y,x,_,_) -> y .= x 
+    nonlinmodel_c = NonLinModel(f!, h!, 500, 1, 1, 1)
 
     transcription = TrapezoidalCollocation(0, f_threads=true, h_threads=true)
     nmpc5 = NonLinMPC(nonlinmodel_c; Nwt=[0], Hp=100, Hc=1, transcription)
@@ -1034,6 +1088,20 @@ end
     preparestate!(nmpc5_1, [0.0])
     u = moveinput!(nmpc5_1, [1/0.001])
     @test u ≈ [1.0] atol=5e-2
+end
+
+@testitem "NonLinMPC moves & getinfo (NonLinModel, OC)" setup=[SetupMPCtests] begin
+    using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
+    using DifferentiationInterface
+    import FiniteDiff
+    linmodel2 = LinModel([tf(5, [2000, 1]) tf(7, [8000,1])], 3000.0, i_d=[2])
+    f = (x,u,d,model) -> model.A*x + model.Bu*u + model.Bd*d
+    h = (x,d,model)   -> model.C*x + model.Dd*d
+    nonlinmodel = NonLinModel(f, h, 3000.0, 1, 2, 1, 1, solver=nothing, p=linmodel2)
+
+    f! = (ẋ,x,u,_,_) -> ẋ .= -0.001x .+ u 
+    h! = (y,x,_,_) -> y .= x 
+    nonlinmodel_c = NonLinModel(f!, h!, 500, 1, 1, 1)
 
     transcription = OrthogonalCollocation(0, 4)
     nmpc6 = NonLinMPC(nonlinmodel_c; Nwt=[0], Hp=100, Hc=1, transcription)
@@ -1049,33 +1117,6 @@ end
     setstate!(nmpc6_1, [5.0])
     moveinput!(nmpc6_1)
     @test nmpc6_1.con.FS ≈ nmpc6_1.con.KS*nmpc6_1.estim.x̂0 # OC + IMC: special lin. eq. method
-    
-    nonlinmodel2 = NonLinModel{Float32}(f, h, 3000.0, 1, 2, 1, 1, solver=nothing, p=linmodel2)
-    nmpc7  = NonLinMPC(nonlinmodel2, Hp=10)
-    y = similar(nonlinmodel2.yop)
-    ModelPredictiveControl.h!(y, nonlinmodel2, Float32[0,0], Float32[0], nonlinmodel2.p)
-    preparestate!(nmpc7, [0], [0])
-    @test moveinput!(nmpc7, [0], [0]) ≈ [0.0] atol=5e-2
-
-    nmpc10 = setconstraint!(NonLinMPC(
-        nonlinmodel, Nwt=[0], Hp=100, Hc=1, 
-        gradient=AutoFiniteDiff(),
-        jacobian=AutoFiniteDiff(),
-        hessian=AutoFiniteDiff()
-        ), ymax=[100], ymin=[-100]
-    )
-    preparestate!(nmpc10, [0], [0])
-    u = moveinput!(nmpc10, [10], [0])
-    @test u ≈ [2] atol=5e-2
-    info = getinfo(nmpc10)
-    @test info[:u] ≈ u
-    @test info[:Ŷ][end] ≈ 10 atol=5e-2
-    
-    nmpc11 = NonLinMPC(nonlinmodel, Hp=10, Hc=[1, 2, 3, 4], Nwt=[10])
-    preparestate!(nmpc11, y, [0])
-    moveinput!(nmpc11, [10], [0])
-    ΔU_diff = diff(getinfo(nmpc11)[:U])
-    @test ΔU_diff[[2, 4, 5, 7, 8, 9]] ≈ zeros(6) atol=1e-9
 end
 
 @testitem "NonLinMPC step disturbance rejection" setup=[SetupMPCtests] begin
diff --git a/test/5_test_extensions.jl b/test/5_test_extensions.jl
index 1384c0bd5..ee4ec83be 100644
--- a/test/5_test_extensions.jl
+++ b/test/5_test_extensions.jl
@@ -1,3 +1,4 @@
+# TODO: add `skip=(Sys.ARCH==:x86)` because bug in DAQP.jl, and activate Linux x86 CI job in juliaci.yml
 @testitem "LinearMPCext general" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra, JuMP, DAQP
     import LinearMPC
diff --git a/test/7_test_aqua.jl b/test/7_test_aqua.jl
index 2eca61e33..fb907c80b 100644
--- a/test/7_test_aqua.jl
+++ b/test/7_test_aqua.jl
@@ -1,9 +1,50 @@
-@testitem "Aqua" begin
-    using Aqua
-    # All the functions defined inside `ModelPredictiveControl.get_nonlinobj_op` with the 
-    # `Vararg`s have unbound type parameters. This is necessary for the splatting syntax of
-    # `JuMP.@operator`, and JuMP will never call these functions with 0 argument, so
-    # defining zero-argument methods would be useless, so I disable this check here.
-    unbound_args = false
-    Aqua.test_all(ModelPredictiveControl; unbound_args)
+@testitem "Aqua ambiguities" begin
+    using Aqua
+    Aqua.test_ambiguities(ModelPredictiveControl)
+end
+
+### All the functions defined inside `ModelPredictiveControl.get_nonlinobj_op` with the 
+### `Vararg`s have unbound type parameters. This is necessary for the splatting syntax of
+### `JuMP.@operator`, and JuMP will never call these functions with 0 argument, so
+### defining zero-argument methods would be useless, so I disable this check here.
+# @testitem "Aqua unbound args" begin
+#     using Aqua
+#     Aqua.test_unbound_args(ModelPredictiveControl)
+# end
+
+@testitem "Aqua undefined exports" begin
+    using Aqua
+    Aqua.test_undefined_exports(ModelPredictiveControl)
+end
+
+@testitem "Aqua project extras" begin
+    using Aqua
+    Aqua.test_project_extras(ModelPredictiveControl)
+end
+
+@testitem "Aqua stale deps" begin
+    using Aqua
+    Aqua.test_stale_deps(ModelPredictiveControl)
+end
+
+@testitem "Aqua deps compat" begin
+    using Aqua
+    Aqua.test_deps_compat(ModelPredictiveControl)
+end
+
+@testitem "Aqua piracies" begin
+    using Aqua
+    Aqua.test_piracies(ModelPredictiveControl)
+end
+
+### This package does not use a __init__ function, and will very likely never use one,
+### so the persistent task test is useless (and slow, because of the long precompilation)
+# @testitem "Aqua persistent tasks" begin
+#     using Aqua
+#     Aqua.test_persistent_tasks(ModelPredictiveControl)
+# end
+
+@testitem "Aqua undocumented names" begin
+    using Aqua
+    Aqua.test_undocumented_names(ModelPredictiveControl)
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 5a707242e..ca1e49c84 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,8 +1,10 @@
 using ModelPredictiveControl
 using Test, TestItems, TestItemRunner
 
-@run_package_tests(verbose=true)
+@run_package_tests
 
+# Not needed for TestItems discovery, but including the files means `Pkg.test` parses them 
+# and catches syntax errors early:
 include("0_test_module.jl")
 include("1_test_sim_model.jl")
 include("2_test_state_estim.jl")