diff --git a/Project.toml b/Project.toml
index 2d467c9fc..bf0edac79 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,6 +28,7 @@ PrecompileTools = "1"
 ProgressLogging = "0.1"
 Random = "1.6"
 RecipesBase = "1"
+TestItemRunner = "1.1"
 julia = "1.10"
 
 [extras]
@@ -35,6 +36,7 @@ DAQP = "c47d62df-3981-49c8-9651-128b1cd08617"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test", "Documenter", "Plots", "DAQP"]
+test = ["Test", "TestItemRunner", "Documenter", "Plots", "DAQP"]
diff --git a/test/runtests.jl b/test/runtests.jl
index ac55e1346..c2efcd985 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,9 +7,26 @@ using LinearAlgebra
 using Random: randn
 using JuMP, OSQP, Ipopt, DAQP, ForwardDiff
 using Plots
-using Test
+using Test, TestItemRunner
+
+@run_package_tests 
+
+#@testitem "ModelPredictiveControl.jl" begin
+
+@testsnippet SetupMPCtests begin
+    using ControlSystemsBase, LinearAlgebra
+    using Random: randn
+    using JuMP, OSQP, Ipopt, DAQP, ForwardDiff
+    using Plots
+    Ts = 400.0
+    sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
+            tf(-0.74,[800.0,1])   tf(0.74,[800.0,1])    tf(-0.74,[800.0,1])   ] 
+    sys_ss = minreal(ss(sys))
+    Gss = c2d(sys_ss[:,1:2], Ts, :zoh)
+    Gss2 = c2d(sys_ss[:,1:2], 0.5Ts, :zoh)
+end
+
 
-@testset "ModelPredictiveControl.jl" begin
 include("test_sim_model.jl")
 include("test_state_estim.jl")
 include("test_predictive_control.jl")
@@ -29,6 +46,6 @@ DocMeta.setdocmeta!(
 doctest(ModelPredictiveControl, testset="DocTest")
 ENV["JULIA_DEBUG"] = old_debug_level
 
-end;
+#end;
 
 nothing
\ No newline at end of file
diff --git a/test/test_plot_sim.jl b/test/test_plot_sim.jl
index 3450b7bba..dac5e9710 100644
--- a/test/test_plot_sim.jl
+++ b/test/test_plot_sim.jl
@@ -1,8 +1,4 @@
-Ts = 400.0
-sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
-        tf(-0.74,[800.0,1])   tf(0.74,[800.0,1])    tf(-0.74,[800.0,1])   ]
-
-@testset "SimModel quick simulation" begin
+@testitem "SimModel quick simulation" setup=[SetupMPCtests] begin
     model = LinModel(sys, Ts, i_d=[3])
     res = sim!(model, 15)
     display(res)
@@ -21,7 +17,7 @@ sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
     @test_throws ArgumentError SimResult(model, [res.U_data model.uop], res.Y_data, res.D_data)
 end
 
-@testset "SimModel Plots" begin
+@testitem "SimModel Plots" setup=[SetupMPCtests] begin
     model = LinModel(sys, Ts, i_d=[3])
     res = sim!(model, 15, [1, 3], [-10])
     p = plot(res, plotx=true)
@@ -45,7 +41,7 @@ end
     @test p[3][1][:y] ≈ res.X_data[4, :]
 end
 
-@testset "StateEstimator quick simulation" begin
+@testitem "StateEstimator quick simulation" setup=[SetupMPCtests] begin
     estim = SteadyKalmanFilter(LinModel(sys, Ts, i_d=[3]))
     res = sim!(estim, 15)
     @test isa(res.obj, SteadyKalmanFilter)
@@ -67,7 +63,7 @@ end
     @test res_man.X̂_data ≈ res.X̂_data
 end
 
-@testset "StateEstimator Plots" begin
+@testitem "StateEstimator Plots" setup=[SetupMPCtests] begin
     estim = MovingHorizonEstimator(LinModel(sys, Ts, i_d=[3]), He=5)
     estim = setconstraint!(estim, x̂min=[-100,-101,-102,-103,-Inf,-Inf])
     estim = setconstraint!(estim, x̂max=[+100,+101,+102,+103,+Inf,+Inf])
@@ -124,7 +120,7 @@ end
     @test all(p7[end-2][3][:y] .≈ +103)
 end
 
-@testset "PredictiveController quick simulation" begin
+@testitem "PredictiveController quick simulation" setup=[SetupMPCtests] begin
     mpc1 = LinMPC(LinModel(sys, Ts, i_d=[3]))
     res = sim!(mpc1, 15)
     @test isa(res.obj, LinMPC)
@@ -158,7 +154,7 @@ end
     @test res_man.Ry_data ≈ res.Ry_data
 end
 
-@testset "PredictiveController Plots" begin
+@testitem "PredictiveController Plots" setup=[SetupMPCtests] begin
     estim = MovingHorizonEstimator(LinModel(sys, Ts, i_d=[3]), He=5)
     estim = setconstraint!(estim, x̂min=[-100,-101,-102,-103,-104,-105])
     estim = setconstraint!(estim, x̂max=[+100,+101,+102,+103,+104,+105])
diff --git a/test/test_predictive_control.jl b/test/test_predictive_control.jl
index 278082370..ed6145aab 100644
--- a/test/test_predictive_control.jl
+++ b/test/test_predictive_control.jl
@@ -1,8 +1,4 @@
-Ts = 400.0
-sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
-        tf(-0.74,[800.0,1])   tf(0.74,[800.0,1])    tf(-0.74,[800.0,1])   ]
-
-@testset "LinMPC construction" begin
+@testitem "LinMPC construction" setup=[SetupMPCtests] begin
     model = LinModel(sys, Ts, i_d=[3])
     mpc1 = LinMPC(model, Hp=15)
     @test isa(mpc1.estim, SteadyKalmanFilter)
@@ -58,7 +54,7 @@ sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
     @test_throws ArgumentError LinMPC(model, Cwt=-1)
 end
 
-@testset "LinMPC moves and getinfo" begin
+@testitem "LinMPC moves and getinfo" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(tf(5, [2, 1]), 3), yop=[10])
     mpc1 = LinMPC(linmodel, Nwt=[0], Hp=1000, Hc=1)
     r = [15]
@@ -91,7 +87,7 @@ end
     @test_throws DimensionMismatch moveinput!(mpc3; R̂u = fill(0, mpc1.Hp+1))
 end
 
-@testset "LinMPC step disturbance rejection" begin
+@testitem "LinMPC step disturbance rejection" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(tf(5, [2, 1]), 3.0), yop=[10])
     r = [15]
     outdist = [5]
@@ -99,9 +95,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(mpc_im, ym)
-        u = moveinput!(mpc_im, r)
+        global u = moveinput!(mpc_im, r)
         updatestate!(mpc_im, u, ym)
         updatestate!(linmodel, u)
     end
@@ -111,9 +107,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(mpc_nint_u, ym)
-        u = moveinput!(mpc_nint_u, r)
+        global u = moveinput!(mpc_nint_u, r)
         updatestate!(mpc_nint_u, u, ym)
         updatestate!(linmodel, u)
     end
@@ -123,9 +119,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(mpc_nint_ym, ym)
-        u = moveinput!(mpc_nint_ym, r)
+        global u = moveinput!(mpc_nint_ym, r)
         updatestate!(mpc_nint_ym, u, ym)
         updatestate!(linmodel, u)
     end
@@ -133,7 +129,7 @@ end
     @test ym ≈ r   atol=1e-2 
 end
 
-@testset "LinMPC other methods" begin
+@testitem "LinMPC other methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     mpc1 = LinMPC(linmodel1)
     @test initstate!(mpc1, [10, 50], [50, 30+1]) ≈ [zeros(3); [1]]
@@ -151,7 +147,7 @@ end
     @test_throws ArgumentError updatestate!(mpc1, [0,0])
 end
 
-@testset "LinMPC set constraints" begin
+@testitem "LinMPC set constraints" setup=[SetupMPCtests] begin
     model = LinModel(sys, Ts, i_d=[3])
     mpc = LinMPC(model, Hp=1, Hc=1)
 
@@ -218,7 +214,7 @@ end
     @test_throws ArgumentError setconstraint!(mpc3, c_ymax=[1, 1])
 end
 
-@testset "LinMPC constraint violation" begin
+@testitem "LinMPC constraint violation" setup=[SetupMPCtests] begin
     model = LinModel(tf([2], [10, 1]), 3.0)
     mpc = LinMPC(model, Hp=50, Hc=5)
 
@@ -276,7 +272,7 @@ end
     setconstraint!(mpc, x̂min=[-1e6,-Inf], x̂max=[+1e6,+Inf])
 end
 
-@testset "LinMPC terminal cost" begin
+@testitem "LinMPC terminal cost" setup=[SetupMPCtests] begin
     model = LinModel(ss([0.5 -0.4;0.6 0.5], [1 0;0 1], [1 0; 0 1], 0, 1))
     K = lqr(Discrete, model.A, model.Bu, I, 0.5I)
     M_end = ControlSystemsBase.are(Discrete, model.A, model.Bu, I, 0.5I)
@@ -296,14 +292,14 @@ end
     X_lqr = zeros(2,20)
     x=[1,1]
     for i=1:20
-        u = -K*x
+        global u = -K*x
         X_lqr[:,i] = x
-        x = model.A*x + model.Bu*u
+        global x = model.A*x + model.Bu*u
     end
     @test all(isapprox.(X_mpc, X_lqr, atol=1e-5))
 end
 
-@testset "LinMPC set model" begin
+@testitem "LinMPC set model" setup=[SetupMPCtests] begin
     estim = KalmanFilter(setop!(LinModel(tf(5, [2, 1]), 3), yop=[10], uop=[1]))
     mpc = LinMPC(estim, Nwt=[0], Cwt=1e4, Hp=1000, Hc=1)
     mpc = setconstraint!(mpc, umin=[-24], umax=[26])
@@ -342,7 +338,7 @@ end
     @test mpc.weights.L_Hp ≈ diagm(1.1:1000.1)
 end
 
-@testset "LinMPC real-time simulations" begin
+@testitem "LinMPC real-time simulations" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(tf(2, [10, 1]), 0.1)
     mpc1 = LinMPC(linmodel1)
     times1 = zeros(5)
@@ -355,7 +351,7 @@ end
     @test all(isapprox.(diff(times1[2:end]), 0.1, atol=0.01))
 end
 
-@testset "ExplicitMPC construction" begin
+@testitem "ExplicitMPC construction" setup=[SetupMPCtests] begin
     model = LinModel(sys, Ts, i_d=[3])
     mpc1 = ExplicitMPC(model, Hp=15)
     @test isa(mpc1.estim, SteadyKalmanFilter)
@@ -385,7 +381,7 @@ end
     @test_throws ArgumentError LinMPC(model, Hp=0)
 end
 
-@testset "ExplicitMPC moves and getinfo" begin
+@testitem "ExplicitMPC moves and getinfo" setup=[SetupMPCtests] begin
     mpc1 = ExplicitMPC(LinModel(tf(5, [2, 1]), 3), Nwt=[0], Hp=1000, Hc=1)
     r = [5]
     preparestate!(mpc1, [0])
@@ -412,7 +408,7 @@ end
 end
 
 
-@testset "ExplicitMPC step disturbance rejection" begin
+@testitem "ExplicitMPC step disturbance rejection" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(tf(5, [2, 1]), 3.0), yop=[10])
     r = [15]
     outdist = [5]
@@ -420,9 +416,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(mpc_im, ym)
-        u = moveinput!(mpc_im, r)
+        global u = moveinput!(mpc_im, r)
         updatestate!(mpc_im, u, ym)
         updatestate!(linmodel, u)
     end
@@ -432,9 +428,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(mpc_nint_u, ym)
-        u = moveinput!(mpc_nint_u, r)
+        global u = moveinput!(mpc_nint_u, r)
         updatestate!(mpc_nint_u, u, ym)
         updatestate!(linmodel, u)
     end
@@ -444,9 +440,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(mpc_nint_ym, ym)
-        u = moveinput!(mpc_nint_ym, r)
+        global u = moveinput!(mpc_nint_ym, r)
         updatestate!(mpc_nint_ym, u, ym)
         updatestate!(linmodel, u)
     end
@@ -454,7 +450,7 @@ end
     @test ym ≈ r   atol=1e-2 
 end
 
-@testset "ExplicitMPC other methods" begin
+@testitem "ExplicitMPC other methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     mpc1 = ExplicitMPC(linmodel1)
     @test initstate!(mpc1, [10, 50], [50, 30+1]) ≈ [zeros(3); [1]]
@@ -468,13 +464,13 @@ end
     @test_throws ArgumentError updatestate!(mpc1, [0,0])
 end
 
-@testset "ExplicitMPC constraints" begin
+@testitem "ExplicitMPC constraints" setup=[SetupMPCtests] begin
     model = LinModel(sys, Ts, i_d=[3])
     mpc = ExplicitMPC(model, Hp=1, Hc=1)
     @test_throws ErrorException setconstraint!(mpc, umin=[0.0, 0.0])
 end
 
-@testset "ExplicitMPC set model" begin
+@testitem "ExplicitMPC set model" setup=[SetupMPCtests] begin
     estim = KalmanFilter(setop!(LinModel(tf(5, [2, 1]), 3), yop=[10], uop=[1]))
     mpc = ExplicitMPC(estim, Nwt=[0], Hp=1000, Hc=1)
     @test mpc.Yop ≈ fill(10.0, 1000)
@@ -503,7 +499,7 @@ end
     @test mpc.weights.L_Hp ≈ diagm(1.1:1000.1)
 end
 
-@testset "NonLinMPC construction" begin
+@testitem "NonLinMPC construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_d=[3])
     nmpc0 = NonLinMPC(linmodel1, Hp=15)
     @test isa(nmpc0.estim, SteadyKalmanFilter)
@@ -566,7 +562,7 @@ end
     @test_logs (:warn, Regex(".*")) NonLinMPC(nonlinmodel, Hp=15, gc=(Ue,_,_,_,_)->Ue, nc=0)    
 end
 
-@testset "NonLinMPC moves and getinfo" begin
+@testitem "NonLinMPC moves and getinfo" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(tf(5, [2000, 1]), 3000.0), yop=[10])
     Hp = 1000
     nmpc_lin = NonLinMPC(linmodel, Nwt=[0], Hp=Hp, Hc=1)
@@ -638,7 +634,7 @@ end
     @test_nowarn ModelPredictiveControl.info2debugstr(info)
 end
 
-@testset "NonLinMPC step disturbance rejection" begin
+@testitem "NonLinMPC step disturbance rejection" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(tf(5, [2000, 1]), 3000.0), yop=[10])
     r = [15]
     outdist = [5]
@@ -646,9 +642,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(nmpc_im, ym)
-        u = moveinput!(nmpc_im, r)
+        global u = moveinput!(nmpc_im, r)
         updatestate!(nmpc_im, u, ym)
         updatestate!(linmodel, u)
     end
@@ -658,9 +654,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(nmpc_nint_u, ym)
-        u = moveinput!(nmpc_nint_u, r)
+        global u = moveinput!(nmpc_nint_u, r)
         updatestate!(nmpc_nint_u, u, ym)
         updatestate!(linmodel, u)
     end
@@ -670,9 +666,9 @@ end
     linmodel.x0 .= 0
     ym, u = linmodel() - outdist, [0.0]
     for i=1:25
-        ym = linmodel() - outdist
+        global ym = linmodel() - outdist
         preparestate!(nmpc_nint_ym, ym)
-        u = moveinput!(nmpc_nint_ym, r)
+        global u = moveinput!(nmpc_nint_ym, r)
         updatestate!(nmpc_nint_ym, u, ym)
         updatestate!(linmodel, u)
     end
@@ -680,7 +676,7 @@ end
     @test ym ≈ r   atol=1e-2
 end
 
-@testset "NonLinMPC other methods" begin
+@testitem "NonLinMPC other methods" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     f = (x,u,_,_) -> linmodel.A*x + linmodel.Bu*u
     h = (x,_,_)   -> linmodel.C*x
@@ -697,7 +693,7 @@ end
     @test nmpc1.estim.x̂0 ≈ [0,0,0,0] atol=1e-6
 end
 
-@testset "NonLinMPC set constraints" begin
+@testitem "NonLinMPC set constraints" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_d=[3])
     nmpc_lin = NonLinMPC(linmodel1, Hp=1, Hc=1)
 
@@ -732,7 +728,7 @@ end
 
 end
 
-@testset "NonLinMPC constraint violation" begin
+@testitem "NonLinMPC constraint violation" setup=[SetupMPCtests] begin
     gc(Ue, Ŷe, _ ,p , ϵ) = [p[1]*(Ue[1:end-1] .- 4.2 .- ϵ); p[2]*(Ŷe[2:end] .- 3.14 .- ϵ)]
     Hp=50
 
@@ -876,7 +872,7 @@ end
 
 end
 
-@testset "NonLinMPC set model" begin
+@testitem "NonLinMPC set model" setup=[SetupMPCtests] begin
     estim = KalmanFilter(setop!(LinModel(tf(5, [200, 1]), 300), yop=[10], uop=[1]))
     mpc = NonLinMPC(estim, Nwt=[0], Cwt=1e4, Hp=1000, Hc=1)
     mpc = setconstraint!(mpc, umin=[-24], umax=[26])
@@ -928,7 +924,7 @@ end
     @test_throws ErrorException setmodel!(nmpc, deepcopy(nonlinmodel))
 end
 
-@testset "LinMPC v.s. NonLinMPC" begin
+@testitem "LinMPC v.s. NonLinMPC" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(sys,Ts,i_d=[3]), uop=[10,50], yop=[50,30], dop=[20])
     f = (x,u,d,_) -> linmodel.A*x + linmodel.Bu*u + linmodel.Bd*d
     h = (x,d,_)   -> linmodel.C*x + linmodel.Dd*d
diff --git a/test/test_sim_model.jl b/test/test_sim_model.jl
index 87a47888c..24fc35861 100644
--- a/test/test_sim_model.jl
+++ b/test/test_sim_model.jl
@@ -1,11 +1,4 @@
-Ts = 4.0
-sys = [ tf(1.90,[18.0,1])   tf(1.90,[18.0,1])   tf(1.90,[18.0,1]);
-        tf(-0.74,[8.0,1])   tf(0.74,[8.0,1])    tf(-0.74,[8.0,1])   ] 
-sys_ss = minreal(ss(sys))
-Gss = c2d(sys_ss[:,1:2], Ts, :zoh)
-Gss2 = c2d(sys_ss[:,1:2], 0.5Ts, :zoh)
-
-@testset "LinModel construction" begin
+@testitem "LinModel construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys, Ts, i_u=1:2)
     @test linmodel1.nx == 2
     @test linmodel1.nu == 2
@@ -29,7 +22,7 @@ Gss2 = c2d(sys_ss[:,1:2], 0.5Ts, :zoh)
     @test linmodel2.dop ≈ zeros(0,1)
 
     linmodel3 = LinModel(Gss, 0.5Ts)
-    @test linmodel3.Ts == 2.0
+    @test linmodel3.Ts == 200.0
     @test linmodel3.A ≈ Gss2.A
     @test linmodel3.C ≈ Gss2.C
 
@@ -55,7 +48,7 @@ Gss2 = c2d(sys_ss[:,1:2], 0.5Ts, :zoh)
     @test linmodel5.yop ≈ [50,30]
     @test linmodel5.dop ≈ [20]
 
-    linmodel6 = LinModel([delay(4) delay(4)]*sys,Ts,i_d=[3])
+    linmodel6 = LinModel([delay(Ts) delay(Ts)]*sys,Ts,i_d=[3])
     @test linmodel6.nx == 3
     @test sum(eigvals(linmodel6.A) .≈ 0) == 1
 
@@ -103,7 +96,7 @@ Gss2 = c2d(sys_ss[:,1:2], 0.5Ts, :zoh)
     @test_throws ErrorException LinModel(sys_ss,Ts)
 end
 
-@testset "LinModel sim methods" begin
+@testitem "LinModel sim methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(Gss), uop=[10,50], yop=[50,30])
     @test updatestate!(linmodel1, [10, 50]) ≈ zeros(2)
     @test updatestate!(linmodel1, [10, 50], Float64[]) ≈ zeros(2)
@@ -123,7 +116,7 @@ end
     @test_throws DimensionMismatch evaloutput(linmodel1, zeros(1))
 end
 
-@testset "LinModel real time simulations" begin
+@testitem "LinModel real time simulations" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(tf(2, [10, 1]), 0.1)
     times1 = zeros(5)
     for i=1:5
@@ -142,7 +135,7 @@ end
     @test all(isapprox.(diff(times2[2:end]), 0.001, atol=0.0001))
 end
 
-@testset "NonLinModel construction" begin
+@testitem "NonLinModel construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     f1(x,u,_,model) = model.A*x + model.Bu*u
     h1(x,_,model)   = model.C*x
@@ -249,7 +242,7 @@ end
         (x,_)->linmodel1.C*x, Ts, 2, 4, 2, 1, solver=nothing)
 end
 
-@testset "NonLinModel sim methods" begin
+@testitem "NonLinModel sim methods" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     f1(x,u,_,model) = model.A*x + model.Bu*u
     h1(x,_,model)   = model.C*x
@@ -268,7 +261,7 @@ end
     @test_throws DimensionMismatch evaloutput(nonlinmodel, zeros(1))
 end
 
-@testset "NonLinModel linearization" begin
+@testitem "NonLinModel linearization" setup=[SetupMPCtests] begin
     Ts = 1.0
     f1(x,u,d,_) = x.^5 + u.^4 + d.^3
     h1(x,d,_)   = x.^2 + d
@@ -313,17 +306,17 @@ end
     linmodel3 = linearize(nonlinmodel3; x, u, d)
     for i=1:N
         ynl = nonlinmodel3(d)
-        yl  = linmodel3(d)
+        global yl  = linmodel3(d)
         Ynl[i] = ynl[1]
         Yl[i]  = yl[1]
-        linmodel3 = linearize(nonlinmodel3; u, d)
+        global linmodel3 = linearize(nonlinmodel3; u, d)
         updatestate!(nonlinmodel3, u, d)
         updatestate!(linmodel3, u, d)
     end
     @test all(isapprox.(Ynl, Yl, atol=1e-6))
 end
 
-@testset "NonLinModel real time simulations" begin
+@testitem "NonLinModel real time simulations" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(tf(2, [10, 1]), 0.1)
     nonlinmodel1 = NonLinModel(
         (x,u,_,_)->linmodel1.A*x + linmodel1.Bu*u,
diff --git a/test/test_state_estim.jl b/test/test_state_estim.jl
index d130e450c..c56a2f9e7 100644
--- a/test/test_state_estim.jl
+++ b/test/test_state_estim.jl
@@ -1,8 +1,4 @@
-Ts = 400.0
-sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
-        tf(-0.74,[800.0,1])   tf(0.74,[800.0,1])    tf(-0.74,[800.0,1])   ] 
-
-@testset "SteadyKalmanFilter construction" begin
+@testitem "SteadyKalmanFilter construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     skalmanfilter1 = SteadyKalmanFilter(linmodel1)
     @test skalmanfilter1.nym == 2
@@ -64,7 +60,7 @@ sys = [ tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1])   tf(1.90,[1800.0,1]);
     @test_throws ErrorException SteadyKalmanFilter(linmodel1, nint_u=[1,1], nint_ym=[1,1])
 end
 
-@testset "SteadyKalmanFilter estimator methods" begin
+@testitem "SteadyKalmanFilter estimator methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     skalmanfilter1 = SteadyKalmanFilter(linmodel1, nint_ym=[1, 1])
     preparestate!(skalmanfilter1, [50, 30])
@@ -115,7 +111,7 @@ end
     @test_throws ArgumentError updatestate!(skalmanfilter1, [10, 50])
 end 
 
-@testset "SteadyKalmanFilter set model" begin
+@testitem "SteadyKalmanFilter set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     skalmanfilter = SteadyKalmanFilter(linmodel, nint_ym=0)
@@ -125,7 +121,7 @@ end
     @test_throws ErrorException setmodel!(skalmanfilter, linmodel, R̂=[0.01])
 end
 
-@testset "SteadyKalmanFilter real-time simulations" begin
+@testitem "SteadyKalmanFilter real-time simulations" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(tf(2, [10, 1]), 0.1)
     skalmanfilter1 = SteadyKalmanFilter(linmodel1)
     times1 = zeros(5)
@@ -138,7 +134,7 @@ end
     @test all(isapprox.(diff(times1[2:end]), 0.1, atol=0.01))
 end
     
-@testset "KalmanFilter construction" begin
+@testitem "KalmanFilter construction" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     kalmanfilter1 = KalmanFilter(linmodel1)
     @test kalmanfilter1.nym == 2
@@ -189,7 +185,7 @@ end
     @test_throws ErrorException KalmanFilter(linmodel1, nint_ym=0, σP_0=[1])
 end
 
-@testset "KalmanFilter estimator methods" begin
+@testitem "KalmanFilter estimator methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     kalmanfilter1 = KalmanFilter(linmodel1)
     preparestate!(kalmanfilter1, [50, 30])
@@ -235,7 +231,7 @@ end
     @test_throws ArgumentError updatestate!(kalmanfilter1, [10, 50])
 end
 
-@testset "KalmanFilter set model" begin
+@testitem "KalmanFilter set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     kalmanfilter = KalmanFilter(linmodel, nint_ym=0)
@@ -263,7 +259,7 @@ end
     @test kalmanfilter.R̂ ≈ [1e-6]
 end
 
-@testset "Luenberger construction" begin
+@testitem "Luenberger construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     lo1 = Luenberger(linmodel1)
     @test lo1.nym == 2
@@ -303,7 +299,7 @@ end
     @test_throws ErrorException Luenberger(LinModel(tf(1,[1, 0]),0.1), poles=[0.5,0.6])
 end
     
-@testset "Luenberger estimator methods" begin
+@testitem "Luenberger estimator methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     lo1 = Luenberger(linmodel1, nint_ym=[1, 1])
     preparestate!(lo1, [50, 30])
@@ -348,7 +344,7 @@ end
     @test isa(x̂, Vector{Float32})
 end
 
-@testset "Luenberger set model" begin
+@testitem "Luenberger set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     lo = Luenberger(linmodel, nint_ym=0)
@@ -356,7 +352,7 @@ end
     @test_throws ErrorException setmodel!(lo, deepcopy(linmodel))
 end
 
-@testset "InternalModel construction" begin
+@testitem "InternalModel construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     internalmodel1 = InternalModel(linmodel1)
     @test internalmodel1.nym == 2
@@ -425,7 +421,7 @@ end
     @test_throws ErrorException InternalModel(linmodel1, stoch_ym=ss(1,1,1,0,Ts).*I(2))
 end    
     
-@testset "InternalModel estimator methods" begin
+@testitem "InternalModel estimator methods" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]) , uop=[10,50], yop=[50,30])
     internalmodel1 = InternalModel(linmodel1)
     preparestate!(internalmodel1, [50, 30] .+ 1)
@@ -462,7 +458,7 @@ end
     @test isa(x̂, Vector{Float32})
 end
 
-@testset "InternalModel set model" begin
+@testitem "InternalModel set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     internalmodel = InternalModel(linmodel)
@@ -487,7 +483,7 @@ end
     @test internalmodel.x̂0 ≈ [3.0 - 8.0]
 end
  
-@testset "UnscentedKalmanFilter construction" begin
+@testitem "UnscentedKalmanFilter construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_d=[3])
     f(x,u,d,_) = linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
     h(x,d,_)   = linmodel1.C*x + linmodel1.Du*d
@@ -544,7 +540,7 @@ end
     @test isa(ukf10, UnscentedKalmanFilter{Float32})
 end
 
-@testset "UnscentedKalmanFilter estimator methods" begin
+@testitem "UnscentedKalmanFilter estimator methods" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     f(x,u,_,_) = linmodel1.A*x + linmodel1.Bu*u
     h(x,_,_)   = linmodel1.C*x
@@ -592,7 +588,7 @@ end
     @test isa(x̂, Vector{Float32})
 end
 
-@testset "UnscentedKalmanFilter set model" begin
+@testitem "UnscentedKalmanFilter set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     ukf1 = UnscentedKalmanFilter(linmodel, nint_ym=0)
@@ -628,7 +624,7 @@ end
     @test_throws ErrorException setmodel!(ukf2, deepcopy(nonlinmodel))
 end
 
-@testset "ExtendedKalmanFilter construction" begin
+@testitem "ExtendedKalmanFilter construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_d=[3])
     f(x,u,d,_) = linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
     h(x,d,_)   = linmodel1.C*x + linmodel1.Du*d
@@ -681,7 +677,7 @@ end
     @test isa(ekf8, ExtendedKalmanFilter{Float32})
 end
 
-@testset "ExtendedKalmanFilter estimator methods" begin
+@testitem "ExtendedKalmanFilter estimator methods" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     f(x,u,_,_) = linmodel1.A*x + linmodel1.Bu*u
     h(x,_,_)   = linmodel1.C*x
@@ -729,7 +725,7 @@ end
     @test isa(x̂, Vector{Float32})
 end
 
-@testset "ExtendedKalmanFilter set model" begin
+@testitem "ExtendedKalmanFilter set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     ekf1 = ExtendedKalmanFilter(linmodel, nint_ym=0)
@@ -765,7 +761,7 @@ end
     @test_throws ErrorException setmodel!(ekf2, deepcopy(nonlinmodel))
 end
 
-@testset "MovingHorizonEstimator construction" begin
+@testitem "MovingHorizonEstimator construction" setup=[SetupMPCtests] begin
     linmodel1 = LinModel(sys,Ts,i_d=[3])
     f(x,u,d,_) = linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
     h(x,d,_)   = linmodel1.C*x + linmodel1.Du*d
@@ -846,7 +842,7 @@ end
     @test_throws ArgumentError MovingHorizonEstimator(linmodel1, Cwt=-1)
 end
 
-@testset "MovingHorizonEstimator estimation and getinfo" begin
+@testitem "MovingHorizonEstimator estimation and getinfo" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2], i_d=[3]), uop=[10,50], yop=[50,30], dop=[5])
     f(x,u,d,_) = linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
     h(x,d,_)   = linmodel1.C*x + linmodel1.Dd*d
@@ -976,7 +972,7 @@ end
     @test_nowarn ModelPredictiveControl.info2debugstr(info)
 end
 
-@testset "MovingHorizonEstimator fallbacks for arrival covariance estimation" begin
+@testitem "MovingHorizonEstimator fallbacks for arrival covariance estimation" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(sys,Ts,i_u=[1,2], i_d=[3]), uop=[10,50], yop=[50,30], dop=[5])
     f(x,u,d,_) = linmodel.A*x + linmodel.Bu*u + linmodel.Bd*d
     h(x,d,_)   = linmodel.C*x + linmodel.Dd*d
@@ -1020,7 +1016,7 @@ end
     @test mhe.invP̄ ≈ invP̄_copy
 end
 
-@testset "MovingHorizonEstimator set constraints" begin
+@testitem "MovingHorizonEstimator set constraints" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     mhe1 = MovingHorizonEstimator(linmodel1, He=1, nint_ym=0, Cwt=1e3)
     setconstraint!(mhe1, x̂min=[-51,-52], x̂max=[53,54])
@@ -1101,7 +1097,7 @@ end
     @test_throws ArgumentError setconstraint!(mhe4, c_v̂max=[1,1])
 end
 
-@testset "MovingHorizonEstimator constraint violation" begin
+@testitem "MovingHorizonEstimator constraint violation" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
     mhe = MovingHorizonEstimator(linmodel1, He=1, nint_ym=0)
 
@@ -1199,7 +1195,7 @@ end
     @test info[:V̂] ≈ [-1,-1] atol=5e-2
 end
 
-@testset "MovingHorizonEstimator set model" begin
+@testitem "MovingHorizonEstimator set model" setup=[SetupMPCtests] begin
     linmodel = LinModel(ss(0.5, 0.3, 1.0, 0, 10.0))
     linmodel = setop!(linmodel, uop=[2.0], yop=[50.0], xop=[3.0], fop=[3.0])
     mhe = MovingHorizonEstimator(linmodel, He=1, nint_ym=0, direct=false)
@@ -1243,7 +1239,7 @@ end
     @test_throws ErrorException setmodel!(mhe2, deepcopy(nonlinmodel))
 end
 
-@testset "MovingHorizonEstimator v.s. Kalman filters" begin
+@testitem "MovingHorizonEstimator v.s. Kalman filters" setup=[SetupMPCtests] begin
     linmodel1 = setop!(LinModel(sys,Ts,i_d=[3]), uop=[10,50], yop=[50,30], dop=[20])
     mhe = MovingHorizonEstimator(linmodel1, He=3, nint_ym=0, direct=false)
     kf  = KalmanFilter(linmodel1, nint_ym=0, direct=false)
@@ -1319,7 +1315,7 @@ end
     @test X̂_mhe ≈ X̂_ekf atol=1e-3 rtol=1e-3 
 end
 
-@testset "MovingHorizonEstimator LinModel v.s. NonLinModel" begin
+@testitem "MovingHorizonEstimator LinModel v.s. NonLinModel" setup=[SetupMPCtests] begin
     linmodel = setop!(LinModel(sys,Ts,i_d=[3]), uop=[10,50], yop=[50,30], dop=[20])
     f = (x,u,d,_) -> linmodel.A*x + linmodel.Bu*u + linmodel.Bd*d
     h = (x,d,_)   -> linmodel.C*x + linmodel.Dd*d