diff --git a/Project.toml b/Project.toml
index 376738f0a..7c2b7a18b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ModelPredictiveControl"
 uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c"
-version = "2.4.3"
+version = "2.4.4"
 authors = ["Francis Gagnon"]
 
 [deps]
diff --git a/src/controller/construct.jl b/src/controller/construct.jl
index d9002cd3c..2338cfbe9 100644
--- a/src/controller/construct.jl
+++ b/src/controller/construct.jl
@@ -248,17 +248,20 @@ julia> mpc = setconstraint!(mpc, umin=[0], umax=[100], Δumin=[-10], Δumax=[+10
 LinMPC controller with a sample time Ts = 4.0 s:
 ├ estimator: SteadyKalmanFilter
 ├ model: LinModel
-├ optimizer: OSQP
+├ optimizer: OSQP 
 ├ transcription: SingleShooting
 └ dimensions:
-  ├ 10 prediction steps Hp
-  ├  2 control steps Hc
-  ├  1 slack variable ϵ (control constraints)
-  ├  1 manipulated inputs u (0 integrating states)
-  ├  2 estimated states x̂
-  ├  1 measured outputs ym (1 integrating states)
-  ├  0 unmeasured outputs yu
-  └  0 measured disturbances d
+  │ ├ 10 prediction steps Hp
+  │ ├  2 control steps Hc
+  │ ├  1 manipulated inputs u (0 integrating states)
+  │ ├  2 estimated states x̂
+  │ ├  1 measured outputs ym (1 integrating states)
+  │ ├  0 unmeasured outputs yu
+  │ └  0 measured disturbances d
+  └ optimization:
+    ├  3 decision variables Z̃ (1 slack variable)
+    ├ 25 linear inequality constraints A (0 custom)
+    └  0 linear equality constraints Aeq
 ```
 
 # Extended Help
diff --git a/src/controller/explicitmpc.jl b/src/controller/explicitmpc.jl
index eff16ef88..cd264cd8c 100644
--- a/src/controller/explicitmpc.jl
+++ b/src/controller/explicitmpc.jl
@@ -122,6 +122,8 @@ ExplicitMPC controller with a sample time Ts = 4.0 s:
 ├ estimator: SteadyKalmanFilter
 ├ model: LinModel
 └ dimensions:
+  ├ 30 prediction steps Hp
+  ├  1 control steps Hc
   ├  1 manipulated inputs u (0 integrating states)
   ├  4 estimated states x̂
   ├  2 measured outputs ym (2 integrating states)
@@ -188,6 +190,8 @@ function Base.show(io::IO, mpc::ExplicitMPC)
     println(io, "├ estimator: $(nameof(typeof(mpc.estim)))")
     println(io, "├ model: $(nameof(typeof(model)))")
     println(io, "└ dimensions:")
+    println(io, "  ├$(lpad(Hp, n)) prediction steps Hp")
+    println(io, "  ├$(lpad(Hc, n)) control steps Hc")
     print_estim_dim(io, mpc.estim, n)
 end
 
diff --git a/src/controller/linmpc.jl b/src/controller/linmpc.jl
index 12a90d7c9..7b7329963 100644
--- a/src/controller/linmpc.jl
+++ b/src/controller/linmpc.jl
@@ -180,17 +180,20 @@ julia> mpc = LinMPC(model, Mwt=[0, 1], Nwt=[0.5], Hp=30, Hc=1)
 LinMPC controller with a sample time Ts = 4.0 s:
 ├ estimator: SteadyKalmanFilter
 ├ model: LinModel
-├ optimizer: OSQP
+├ optimizer: OSQP 
 ├ transcription: SingleShooting
 └ dimensions:
-  ├ 30 prediction steps Hp
-  ├  1 control steps Hc
-  ├  1 slack variable ϵ (control constraints)
-  ├  1 manipulated inputs u (0 integrating states)
-  ├  4 estimated states x̂
-  ├  2 measured outputs ym (2 integrating states)
-  ├  0 unmeasured outputs yu
-  └  0 measured disturbances d
+  │ ├ 30 prediction steps Hp
+  │ ├  1 control steps Hc
+  │ ├  1 manipulated inputs u (0 integrating states)
+  │ ├  4 estimated states x̂
+  │ ├  2 measured outputs ym (2 integrating states)
+  │ ├  0 unmeasured outputs yu
+  │ └  0 measured disturbances d
+  └ optimization:
+    ├ 2 decision variables Z̃ (1 slack variable)
+    ├ 1 linear inequality constraints A (0 custom)
+    └ 0 linear equality constraints Aeq
 ```
 
 # Extended Help
@@ -267,17 +270,20 @@ julia> mpc = LinMPC(estim, Mwt=[0, 1], Nwt=[0.5], Hp=30, Hc=1)
 LinMPC controller with a sample time Ts = 4.0 s:
 ├ estimator: KalmanFilter
 ├ model: LinModel
-├ optimizer: OSQP
+├ optimizer: OSQP 
 ├ transcription: SingleShooting
 └ dimensions:
-  ├ 30 prediction steps Hp
-  ├  1 control steps Hc
-  ├  1 slack variable ϵ (control constraints)
-  ├  1 manipulated inputs u (0 integrating states)
-  ├  3 estimated states x̂
-  ├  1 measured outputs ym (1 integrating states)
-  ├  1 unmeasured outputs yu
-  └  0 measured disturbances d
+  │ ├ 30 prediction steps Hp
+  │ ├  1 control steps Hc
+  │ ├  1 manipulated inputs u (0 integrating states)
+  │ ├  3 estimated states x̂
+  │ ├  1 measured outputs ym (1 integrating states)
+  │ ├  1 unmeasured outputs yu
+  │ └  0 measured disturbances d
+  └ optimization:
+    ├ 2 decision variables Z̃ (1 slack variable)
+    ├ 1 linear inequality constraints A (0 custom)
+    └ 0 linear equality constraints Aeq
 ```
 """
 function LinMPC(
@@ -331,4 +337,15 @@ function init_optimization!(mpc::LinMPC, model::LinModel, optim::JuMP.GenericMod
     @constraint(optim, linconstrainteq, Aeq*Z̃var .== beq)
     set_objective_hessian!(mpc, model, Z̃var)
     return nothing
+end
+
+"Print the decision variable and linear constraint dimensions for `LinMPC`."
+function print_optim_dim(io::IO, mpc::LinMPC)
+    nZ̃, nϵ = length(mpc.Z̃), mpc.nϵ
+    nA, nW, nAeq = sum(mpc.con.i_b) , mpc.con.nw*(mpc.Hp + 1), size(mpc.con.Aeq, 1)
+    m = maximum(ndigits.((nZ̃, nA, nAeq))) + 1
+    println(io, "  └ optimization:")
+    println(io, "    ├$(lpad(nZ̃, m)) decision variables Z̃ ($nϵ slack variable)")
+    println(io, "    ├$(lpad(nA, m)) linear inequality constraints A ($nW custom)")
+    print(io,   "    └$(lpad(nAeq, m)) linear equality constraints Aeq")
 end
\ No newline at end of file
diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl
index 8c675bc27..b60e86de0 100644
--- a/src/controller/nonlinmpc.jl
+++ b/src/controller/nonlinmpc.jl
@@ -253,14 +253,19 @@ NonLinMPC controller with a sample time Ts = 10.0 s:
 ├ jacobian: AutoSparse (AutoForwardDiff, TracerSparsityDetector, GreedyColoringAlgorithm)
 ├ hessian: nothing
 └ dimensions:
-  ├ 20 prediction steps Hp
-  ├ 10 control steps Hc
-  ├  1 slack variable ϵ (control constraints)
-  ├  1 manipulated inputs u (0 integrating states)
-  ├  2 estimated states x̂
-  ├  1 measured outputs ym (1 integrating states)
-  ├  0 unmeasured outputs yu
-  └  0 measured disturbances d
+  │ ├ 20 prediction steps Hp
+  │ ├ 10 control steps Hc
+  │ ├  1 manipulated inputs u (0 integrating states)
+  │ ├  2 estimated states x̂
+  │ ├  1 measured outputs ym (1 integrating states)
+  │ ├  0 unmeasured outputs yu
+  │ └  0 measured disturbances d
+  └ optimization:
+    ├ 51 decision variables Z̃ (1 slack variable)
+    ├  1 linear inequality constraints A (0 custom)
+    ├ 20 linear equality constraints Aeq
+    ├  0 nonlinear inequality constraints g (0 custom)
+    └ 20 nonlinear equality constraints geq
 ```
 
 # Extended Help
@@ -396,14 +401,19 @@ NonLinMPC controller with a sample time Ts = 10.0 s:
 ├ jacobian: AutoForwardDiff
 ├ hessian: nothing
 └ dimensions:
-  ├ 20 prediction steps Hp
-  ├  2 control steps Hc
-  ├  1 slack variable ϵ (control constraints)
-  ├  1 manipulated inputs u (0 integrating states)
-  ├  2 estimated states x̂
-  ├  1 measured outputs ym (1 integrating states)
-  ├  0 unmeasured outputs yu
-  └  0 measured disturbances d
+  │ ├ 20 prediction steps Hp
+  │ ├  2 control steps Hc
+  │ ├  1 manipulated inputs u (0 integrating states)
+  │ ├  2 estimated states x̂
+  │ ├  1 measured outputs ym (1 integrating states)
+  │ ├  0 unmeasured outputs yu
+  │ └  0 measured disturbances d
+  └ optimization:
+    ├ 3 decision variables Z̃ (1 slack variable)
+    ├ 1 linear inequality constraints A (0 custom)
+    ├ 0 linear equality constraints Aeq
+    ├ 0 nonlinear inequality constraints g (0 custom)
+    └ 0 nonlinear equality constraints geq
 ```
 """
 function NonLinMPC(
@@ -1161,3 +1171,17 @@ function print_backends(io::IO, mpc::NonLinMPC)
     println(io, "├ jacobian: $(backend_str(mpc.jacobian))")
     println(io, "├ hessian: $(backend_str(mpc.hessian))")
 end
+
+"Print the decision variable, linear and nonlinear constraint dimensions for `NonLinMPC`."
+function print_optim_dim(io::IO, mpc::NonLinMPC)
+    nZ̃, nϵ = length(mpc.Z̃), mpc.nϵ
+    nA, nW, nAeq = sum(mpc.con.i_b) , mpc.con.nw*(mpc.Hp + 1), size(mpc.con.Aeq, 1)
+    ng, nc, neq = sum(mpc.con.i_g), mpc.con.nc, mpc.con.neq
+    m = maximum(ndigits.((nZ̃, nA, nAeq, ng, neq))) + 1
+    println(io, "  └ optimization:")
+    println(io, "    ├$(lpad(nZ̃, m)) decision variables Z̃ ($nϵ slack variable)")
+    println(io, "    ├$(lpad(nA, m)) linear inequality constraints A ($nW custom)")
+    println(io, "    ├$(lpad(nAeq, m)) linear equality constraints Aeq")
+    println(io, "    ├$(lpad(ng, m)) nonlinear inequality constraints g ($nc custom)")
+    print(io,   "    └$(lpad(neq, m)) nonlinear equality constraints geq")
+end
\ No newline at end of file
diff --git a/src/estimator/internal_model.jl b/src/estimator/internal_model.jl
index 4079e1713..1ff3c4245 100644
--- a/src/estimator/internal_model.jl
+++ b/src/estimator/internal_model.jl
@@ -358,13 +358,14 @@ function evaloutput(estim::InternalModel, d)
 end
 
 "Print InternalModel information without i/o integrators."
-function print_estim_dim(io::IO, estim::InternalModel, n)
+function print_estim_dim(io::IO, estim::InternalModel, n; firstchars=" ")
     nu, nd = estim.model.nu, estim.model.nd
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
-    println(io, "  ├$(lpad(nu, n)) manipulated inputs u")
-    println(io, "  ├$(lpad(nx̂, n)) estimated states x̂")
-    println(io, "  ├$(lpad(nym, n)) measured outputs ym")
-    println(io, "  ├$(lpad(nyu, n)) unmeasured outputs yu")
-    print(io,   "  └$(lpad(nd, n)) measured disturbances d")
+    f = firstchars
+    println(io, "$f ├$(lpad(nu, n)) manipulated inputs u")
+    println(io, "$f ├$(lpad(nx̂, n)) estimated states x̂")
+    println(io, "$f ├$(lpad(nym, n)) measured outputs ym")
+    println(io, "$f ├$(lpad(nyu, n)) unmeasured outputs yu")
+    print(io,   "$f └$(lpad(nd, n)) measured disturbances d")
 end
 
diff --git a/src/estimator/mhe.jl b/src/estimator/mhe.jl
index c5ecac97d..9122148ec 100644
--- a/src/estimator/mhe.jl
+++ b/src/estimator/mhe.jl
@@ -22,16 +22,26 @@ end
 print_backends(::IO, ::MovingHorizonEstimator, ::LinModel) = nothing
 
 "Print the overall dimensions of the MHE `estim` with left padding `n`."
-function print_estim_dim(io::IO, estim::MovingHorizonEstimator, n)
+function print_estim_dim(io::IO, estim::MovingHorizonEstimator, n; firstchars=nothing)
     nu, nd = estim.model.nu, estim.model.nd
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
     He, nε = estim.He, estim.nε
     niu, niym = sum(estim.nint_u), sum(estim.nint_ym)
-    println(io, "  ├$(lpad(He, n)) estimation steps He")
-    println(io, "  ├$(lpad(nε, n)) slack variable ε (estimation constraints)")
-    println(io, "  ├$(lpad(nu, n)) manipulated inputs u ($niu integrating states)")
-    println(io, "  ├$(lpad(nx̂, n)) estimated states x̂")
-    println(io, "  ├$(lpad(nym, n)) measured outputs ym ($niym integrating states)")
-    println(io, "  ├$(lpad(nyu, n)) unmeasured outputs yu")
-    print(io,   "  └$(lpad(nd, n)) measured disturbances d")
+    println(io, "  │ ├$(lpad(He, n)) estimation steps He")
+    println(io, "  │ ├$(lpad(nu, n)) manipulated inputs u ($niu integrating states)")
+    println(io, "  │ ├$(lpad(nx̂, n)) estimated states x̂")
+    println(io, "  │ ├$(lpad(nym, n)) measured outputs ym ($niym integrating states)")
+    println(io, "  │ ├$(lpad(nyu, n)) unmeasured outputs yu")
+    print(io,   "  │ └$(lpad(nd, n)) measured disturbances d")
+    if isnothing(firstchars) # the user prints the MHE object itself, not a controller:
+        nZ̃, nε = length(estim.Z̃), estim.nε
+        nA = sum(estim.con.i_b)
+        ng, nc = sum(estim.con.i_g), estim.con.nc 
+        m = maximum(ndigits.((nZ̃, nA, ng))) + 1
+        println(io)
+        println(io, "  └ optimization:")
+        println(io, "    ├$(lpad(nZ̃, m)) decision variables Z̃ ($nε slack variable)")
+        println(io, "    ├$(lpad(nA, m)) linear inequality constraints A")
+        print(io,   "    └$(lpad(ng, m)) nonlinear inequality constraints g ($nc custom)")
+    end
 end
\ No newline at end of file
diff --git a/src/estimator/mhe/construct.jl b/src/estimator/mhe/construct.jl
index 7128b9f07..73abc47f2 100644
--- a/src/estimator/mhe/construct.jl
+++ b/src/estimator/mhe/construct.jl
@@ -332,13 +332,16 @@ MovingHorizonEstimator estimator with a sample time Ts = 10.0 s:
 ├ arrival covariance: UnscentedKalmanFilter 
 ├ direct: true
 └ dimensions:
-  ├ 5 estimation steps He
-  ├ 0 slack variable ε (estimation constraints)
-  ├ 1 manipulated inputs u (0 integrating states)
-  ├ 2 estimated states x̂
-  ├ 1 measured outputs ym (1 integrating states)
-  ├ 0 unmeasured outputs yu
-  └ 0 measured disturbances d
+  │ ├ 5 estimation steps He
+  │ ├ 1 manipulated inputs u (0 integrating states)
+  │ ├ 2 estimated states x̂
+  │ ├ 1 measured outputs ym (1 integrating states)
+  │ ├ 0 unmeasured outputs yu
+  │ └ 0 measured disturbances d
+  └ optimization:
+    ├ 12 decision variables Z̃ (0 slack variable)
+    ├  0 linear inequality constraints A
+    └  0 nonlinear inequality constraints g (0 custom)
 ```
 
 # Extended Help
@@ -772,17 +775,20 @@ julia> estim = MovingHorizonEstimator(LinModel(ss(0.5,1,1,0,1)), He=3);
 julia> estim = setconstraint!(estim, x̂min=[-50, -50], x̂max=[50, 50])
 MovingHorizonEstimator estimator with a sample time Ts = 1.0 s:
 ├ model: LinModel
-├ optimizer: OSQP
-├ arrival covariance: KalmanFilter
+├ optimizer: OSQP 
+├ arrival covariance: KalmanFilter 
 ├ direct: true
 └ dimensions:
-  ├ 3 estimation steps He
-  ├ 0 slack variable ε (estimation constraints)
-  ├ 1 manipulated inputs u (0 integrating states)
-  ├ 2 estimated states x̂
-  ├ 1 measured outputs ym (1 integrating states)
-  ├ 0 unmeasured outputs yu
-  └ 0 measured disturbances d
+  │ ├ 3 estimation steps He
+  │ ├ 1 manipulated inputs u (0 integrating states)
+  │ ├ 2 estimated states x̂
+  │ ├ 1 measured outputs ym (1 integrating states)
+  │ ├ 0 unmeasured outputs yu
+  │ └ 0 measured disturbances d
+  └ optimization:
+    ├  8 decision variables Z̃ (0 slack variable)
+    ├ 16 linear inequality constraints A
+    └  0 nonlinear inequality constraints g (0 custom)
 ```
 
 # Extended Help
diff --git a/src/predictive_control.jl b/src/predictive_control.jl
index 11691b1e5..964ce1ff7 100644
--- a/src/predictive_control.jl
+++ b/src/predictive_control.jl
@@ -28,7 +28,7 @@ include("controller/nonlinmpc.jl")
 
 function Base.show(io::IO, mpc::PredictiveController)
     estim, model = mpc.estim, mpc.estim.model
-    Hp, Hc, nϵ = mpc.Hp, mpc.Hc, mpc.nϵ
+    Hp, Hc = mpc.Hp, mpc.Hc
     nu, nd = model.nu, model.nd
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
     other_dims = get_other_dims(estim)
@@ -40,15 +40,19 @@ function Base.show(io::IO, mpc::PredictiveController)
     println(io, "├ transcription: $(nameof(typeof(mpc.transcription)))")
     print_backends(io, mpc)
     println(io, "└ dimensions:")
-    println(io, "  ├$(lpad(Hp, n)) prediction steps Hp")
-    println(io, "  ├$(lpad(Hc, n)) control steps Hc")
-    println(io, "  ├$(lpad(nϵ, n)) slack variable ϵ (control constraints)")
-    print_estim_dim(io, mpc.estim, n)
+    println(io, "  │ ├$(lpad(Hp, n)) prediction steps Hp")
+    println(io, "  │ ├$(lpad(Hc, n)) control steps Hc")
+    print_estim_dim(io, mpc.estim, n, firstchars="  │")
+    println(io) # add a linebreak since `print_estim_dim` ends with a `print` (w/o ln) 
+    print_optim_dim(io, mpc)
 end
 
 "No differentiation backends to print for a `PredictiveController` by default."
 print_backends(::IO, ::PredictiveController) = nothing
 
+"No dimensions related to the optimization problem by default."
+print_optim_dim(io::IO, ::PredictiveController) =  println(io, "  └ optimization: nothing")
+
 "Functor allowing callable `PredictiveController` object as an alias for `moveinput!`."
 function (mpc::PredictiveController)(
     ry::AbstractVector = mpc.estim.model.yop, 
diff --git a/src/state_estim.jl b/src/state_estim.jl
index ce7ee5bcf..c4beb6f89 100644
--- a/src/state_estim.jl
+++ b/src/state_estim.jl
@@ -50,13 +50,14 @@ function print_details(io::IO, estim::StateEstimator)
 end
 
 "Print the overall dimensions of the state estimator `estim` with left padding `n`."
-function print_estim_dim(io::IO, estim::StateEstimator, n)
+function print_estim_dim(io::IO, estim::StateEstimator, n; firstchars=" ")
     nu, nd = estim.model.nu, estim.model.nd
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
     niu, niym = sum(estim.nint_u), sum(estim.nint_ym)
-    println(io, "  ├$(lpad(nu, n)) manipulated inputs u ($niu integrating states)")
-    println(io, "  ├$(lpad(nx̂, n)) estimated states x̂")
-    println(io, "  ├$(lpad(nym, n)) measured outputs ym ($niym integrating states)")
-    println(io, "  ├$(lpad(nyu, n)) unmeasured outputs yu")
-    print(io,   "  └$(lpad(nd, n)) measured disturbances d")
+    f = firstchars
+    println(io, "$f ├$(lpad(nu, n)) manipulated inputs u ($niu integrating states)")
+    println(io, "$f ├$(lpad(nx̂, n)) estimated states x̂")
+    println(io, "$f ├$(lpad(nym, n)) measured outputs ym ($niym integrating states)")
+    println(io, "$f ├$(lpad(nyu, n)) unmeasured outputs yu")
+    print(io,   "$f └$(lpad(nd, n)) measured disturbances d")
 end
\ No newline at end of file