From d6a1d48a461bd06427c8240a3ae1f1de42430059 Mon Sep 17 00:00:00 2001
From: franckgaga <franckgaga2@gmail.com>
Date: Tue, 29 Sep 2026 15:53:35 -0400
Subject: [PATCH 1/5] added: `unit="s"` kwarg in all `SimModel` constructors

It modifies the suffix of the x-axis label in the plot (no underlying time data conversion), and also the unit of time shown when pretty-printing objects.

This is not strictly necessary for plots since we can modify the x-axis label after a is generated, but we still needs to this on all the generated plots, thus simplifying the API in the end.
---
 src/model/linmodel.jl       | 37 +++++++++++++++++++++----------------
 src/model/nonlinmodel.jl    | 18 +++++++++++-------
 src/model/nonlinmodeldae.jl | 30 ++++++++++++++++++------------
 src/predictive_control.jl   |  3 ++-
 src/sim_model.jl            |  3 ++-
 src/state_estim.jl          |  3 ++-
 6 files changed, 56 insertions(+), 38 deletions(-)

diff --git a/src/model/linmodel.jl b/src/model/linmodel.jl
index d0d5c1f66..a38a3e518 100644
--- a/src/model/linmodel.jl
+++ b/src/model/linmodel.jl
@@ -22,9 +22,10 @@ struct LinModel{NT<:Real} <: SimModelODE{NT}
     yname::Vector{String}
     dname::Vector{String}
     xname::Vector{String}
+    unit::String
     xs_0::Vector{NT}
     buffer::SimModelBuffer{NT}
-    function LinModel{NT}(A, Bu, C, Bd, Dd, Ts) where {NT<:Real}
+    function LinModel{NT}(A, Bu, C, Bd, Dd, Ts, unit="s") where {NT<:Real}
         A, Bu = to_mat(A, 1, 1), to_mat(Bu, 1, 1)
         nu, nx = size(Bu, 2), size(A, 2)
         (C == I) && (C = Matrix{NT}(I, nx, nx))
@@ -62,6 +63,7 @@ struct LinModel{NT<:Real} <: SimModelODE{NT}
             nu, nx, ny, nd, nk̄,
             uop, yop, dop, xop, fop,
             uname, yname, dname, xname,
+            unit,
             xs_0,
             buffer
         )
@@ -69,9 +71,9 @@ struct LinModel{NT<:Real} <: SimModelODE{NT}
 end
 
 @doc raw"""
-    LinModel(sys::StateSpace[, Ts]; i_u=1:size(sys,2), i_d=Int[])
+    LinModel(sys::StateSpace[, Ts]; i_u=1:size(sys,2), i_d=Int[], unit="s")
 
-Construct a linear model from state-space model `sys` with sampling time `Ts` in seconds.
+Construct a linear model from state-space model `sys` with sampling time `Ts`.
 
 The system `sys` can be continuous or discrete-time (`Ts` can be omitted for the latter).
 For continuous dynamics, its state-space equations are (discrete case in Extended Help):
@@ -83,11 +85,13 @@ For continuous dynamics, its state-space equations are (discrete case in Extende
 ```
 with the state ``\mathbf{x}`` and output ``\mathbf{y}`` vectors. The ``\mathbf{s}`` vector 
 comprises the manipulated inputs ``\mathbf{u}`` and measured disturbances ``\mathbf{d}``, 
-in any order. `i_u` provides the indices of ``\mathbf{s}`` that are manipulated, and `i_d`, 
-the measured disturbances. The constructor automatically discretizes continuous systems,
-resamples discrete ones if `Ts ≠ sys.Ts`, computes a new balancing and minimal state-space
-realization, and separates the ``\mathbf{s}`` terms in two parts (details in Extended Help). 
-The rest of the documentation assumes discrete models since all systems end up in this form.
+in any order. The keyword argument `i_u` provides the indices of ``\mathbf{s}`` that are
+manipulated, and `i_d`, the measured disturbances. The string `unit` sets the x-label suffix
+in the plots (no underlying time data conversion). The constructor automatically discretizes
+continuous systems, resamples discrete ones if `Ts ≠ sys.Ts`, computes a new balancing and
+minimal state-space realization and separates the ``\mathbf{s}`` terms in two parts (details
+in Extended Help). The rest of the documentation assumes discrete models since all systems
+end up in this form.
 
 See also [`ss`](@extref ControlSystemsBase.ss)
 
@@ -134,7 +138,7 @@ LinModel with a sample time Ts = 0.1 s:
         \mathbf{y}(k)   &=  \mathbf{C x}(k) + \mathbf{D_d d}(k)
     \end{aligned}
     ```
-    Use the syntax [`LinModel{NT}(A, Bu, C, Bd, Dd, Ts)`](@ref) to force a specific
+    Use the syntax [`LinModel{NT}(A, Bu, C, Bd, Dd, Ts, unit)`](@ref) to force a specific
     state-space representation. 
     
     It is assumed that ``\mathbf{D_u=0}`` (or `sys` is strictly proper) since otherwise the
@@ -152,7 +156,8 @@ function LinModel(
     sys::StateSpace{E, NT},
     Ts::Union{Real,Nothing} = nothing;
     i_u::AbstractVector{Int} = 1:size(sys,2),
-    i_d::AbstractVector{Int} = Int[]
+    i_d::AbstractVector{Int} = Int[],
+    unit="s"
 ) where {E, NT<:Real}
     if !isempty(i_d)
         # common indexes in i_u and i_d are interpreted as measured disturbances d :
@@ -198,12 +203,12 @@ function LinModel(
     Bd  = sys_dis.B[:,nu+1:end]
     C   = sys_dis.C
     Dd  = sys_dis.D[:,nu+1:end]
-    return LinModel{NT}(A, Bu, C, Bd, Dd, Ts)
+    return LinModel{NT}(A, Bu, C, Bd, Dd, Ts, unit)
 end
 
 
 @doc raw"""
-    LinModel(sys::TransferFunction[, Ts]; i_u=1:size(sys,2), i_d=Int[])
+    LinModel(sys::TransferFunction[, Ts]; unit="s", i_u=1:size(sys,2), i_d=Int[])
 
 Convert to minimal realization state-space when `sys` is a transfer function.
 
@@ -230,7 +235,7 @@ end
 
 
 """
-    LinModel(sys::DelayLtiSystem, Ts; i_u=1:size(sys,2), i_d=Int[])
+    LinModel(sys::DelayLtiSystem, Ts; unit="s", i_u=1:size(sys,2), i_d=Int[])
 
 Discretize with zero-order hold when `sys` is a continuous system with delays.
 
@@ -242,7 +247,7 @@ function LinModel(sys::DelayLtiSystem, Ts::Real; kwargs...)
 end
 
 @doc raw"""
-    LinModel{NT}(A, Bu, C, Bd, Dd, Ts)
+    LinModel{NT}(A, Bu, C, Bd, Dd, Ts, unit="s")
 
 Construct the model from the discrete state-space matrices `A, Bu, C, Bd, Dd` directly.
 
@@ -252,8 +257,8 @@ syntax do not modify the state-space representation provided in argument (`minre
 called). Care must be taken to ensure that the model is controllable and observable. The 
 optional parameter `NT` explicitly set the number type of vectors (default to `Float64`).
 """
-LinModel{NT}(A, Bu, C, Bd, Dd, Ts) where NT<:Real
-LinModel(A, Bu, C, Bd, Dd, Ts) = LinModel{Float64}(A, Bu, C, Bd, Dd, Ts)
+LinModel{NT}(A, Bu, C, Bd, Dd, Ts, unit) where NT<:Real
+LinModel(A, Bu, C, Bd, Dd, Ts, unit="s") = LinModel{Float64}(A, Bu, C, Bd, Dd, Ts, unit)
 
 function validate_transcription(::LinModel, ::CollocationMethod)
     throw(ArgumentError("Collocation methods are not supported for LinModel."))
diff --git a/src/model/nonlinmodel.jl b/src/model/nonlinmodel.jl
index e9f6b2ab7..035b2f2af 100644
--- a/src/model/nonlinmodel.jl
+++ b/src/model/nonlinmodel.jl
@@ -44,13 +44,15 @@ struct NonLinModel{
     yname::Vector{String}
     dname::Vector{String}
     xname::Vector{String}
+    unit::String
     xs_0::Vector{NT}
     jacobian::JB
     linfunc!::LF
     buffer::SimModelBuffer{NT}
     function NonLinModel{NT}(
-        solver::DS, f!::F, h!::H, Ts, nu, nx, ny, nd, 
-        p::PT, jacobian::JB, linfunc!::LF
+        solver::DS, f!::F, h!::H, Ts, nu, nx, ny, nd,  
+        p::PT, jacobian::JB, linfunc!::LF,
+        unit="s"
     ) where {
             NT<:Real, 
             DS<:DiffSolver,
@@ -85,6 +87,7 @@ struct NonLinModel{
             nu, nx, ny, nd, nk̄, 
             uop, yop, dop, xop, fop,
             uname, yname, dname, xname,
+            unit,
             xs_0,
             jacobian, linfunc!,
             buffer
@@ -138,7 +141,7 @@ See also [`LinModel`](@ref), and [`NonLinModelDAE`](@ref) to include algebraic e
 # Arguments
 - `f::Function` or `f!`: state function of the model.
 - `h::Function` or `h!`: output function of the model.
-- `Ts`: sampling time of the model in seconds.
+- `Ts`: sampling time of the model (see `unit` below).
 - `nu`: number of manipulated inputs.
 - `nx`: number of states.
 - `ny`: number of outputs.
@@ -148,6 +151,7 @@ See also [`LinModel`](@ref), and [`NonLinModelDAE`](@ref) to include algebraic e
   dynamics, use `nothing` for discrete-time models (default to 4th order [`RungeKutta`](@ref)).
 - `jacobian=AutoForwardDiff()`: an `AbstractADType` backend when [`linearize`](@ref) is
    called, see [`DifferentiationInterface` doc](@extref DifferentiationInterface List).
+- `unit="s"`: time unit shown in plots x-label suffix (no underlying time data conversion).
 
 # Examples
 ```jldoctest
@@ -202,7 +206,7 @@ NonLinModel with a sample time Ts = 2.0 s:
 """
 function NonLinModel{NT}(
     f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0;
-    p=NT[], solver=RungeKutta(4), jacobian=AutoForwardDiff()
+    p=NT[], solver=RungeKutta(4), jacobian=AutoForwardDiff(), unit="s"
 ) where {NT<:Real}
     isnothing(solver) && (solver=EmptySolver())
     f!, h! = get_mutating_functions(NT, f, h)
@@ -210,15 +214,15 @@ function NonLinModel{NT}(
         NT, f!, h!, Ts, nu, nx, ny, nd, p, solver, jacobian
     )
     return NonLinModel{NT}(
-        solver, f!, h!, Ts, nu, nx, ny, nd, p, jacobian, linfunc!
+        solver, f!, h!, Ts, nu, nx, ny, nd, p, jacobian, linfunc!, unit
     )
 end
 
 function NonLinModel(
     f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0;
-    p=Float64[], solver=RungeKutta(4), jacobian=AutoForwardDiff()
+    p=Float64[], solver=RungeKutta(4), jacobian=AutoForwardDiff(), unit="s"
 )
-    return NonLinModel{Float64}(f, h, Ts, nu, nx, ny, nd; p, solver, jacobian)
+    return NonLinModel{Float64}(f, h, Ts, nu, nx, ny, nd; p, solver, jacobian, unit)
 end
 
 "Get the mutating functions `f!` and `h!` from the provided functions in argument."
diff --git a/src/model/nonlinmodeldae.jl b/src/model/nonlinmodeldae.jl
index 317387e7b..5fab8733f 100644
--- a/src/model/nonlinmodeldae.jl
+++ b/src/model/nonlinmodeldae.jl
@@ -51,6 +51,7 @@ struct NonLinModelDAE{
     yname::Vector{String}
     dname::Vector{String}
     xname::Vector{String}
+    unit::String
     xs_0::Vector{NT}
     as_0::Vector{NT}
     x0_optim::Vector{NT}
@@ -64,9 +65,9 @@ struct NonLinModelDAE{
         xs_0,
         as_0,
         transcription::TM, 
-        optim_state::JMS,
-        optim_output::JMO,
-        jacobian::JB, hessian::HB
+        optim_state::JMS, optim_output::JMO,
+        jacobian::JB, hessian::HB,
+        unit="s"
     ) where {
             NT<:Real, 
             TM<:CollocationMethod,
@@ -123,6 +124,7 @@ struct NonLinModelDAE{
             nu, nx, na, ny, nd, 
             uop, yop, dop, xop, fop,
             uname, yname, dname, xname,
+            unit,
             xs_0, as_0,
             x0_optim, u0_optim, d0_optim,
             iszero_Ha,
@@ -182,7 +184,7 @@ See also [`NonLinModel`](@ref) for ODEs.
 # Arguments
 - `fq::Function` or `fq!`: combined state and algebraic function of the model.
 - `h::Function` or `h!`: output function of the model.
-- `Ts`: sampling time of the model in seconds.
+- `Ts`: sampling time of the model (see `unit` below).
 - `nu`: number of manipulated inputs.
 - `nx`: number of states.
 - `na`: number of algebraic variables.
@@ -191,17 +193,18 @@ See also [`NonLinModel`](@ref) for ODEs.
 - `p=[]`: parameters of the model (any type).
 - `xs_0=zeros(nx)`: initial guess (optimization warm-start) for the states.
 - `as_0=zeros(na)`: initial guess (optimization warm-start) for the algebraic variables.
-- `transcription=OrthogonalCollocation()` : a [`TrapezoidalCollocation`](@ref) or 
+- `transcription=OrthogonalCollocation()`: a [`TrapezoidalCollocation`](@ref) or 
    [`OrthogonalCollocation`](@ref) instance for open-loop simulations.
-- `optim_state=JuMP.Model(Ipopt.Optimizer)` : nonlinear optimizer for [`updatestate!`](@ref),
+- `optim_state=JuMP.Model(Ipopt.Optimizer)`: nonlinear optimizer for [`updatestate!`](@ref),
    provided as a [`JuMP.Model`](@extref) object (default to [`Ipopt`](https://github.com/jump-dev/Ipopt.jl) optimizer).
-- `optim_output=JuMP.Model(Ipopt.Optimizer)` : nonlinear optimizer for [`evaloutput`](@ref),
+- `optim_output=JuMP.Model(Ipopt.Optimizer)`: nonlinear optimizer for [`evaloutput`](@ref),
    provided as a [`JuMP.Model`](@extref) object (default to [`Ipopt`](https://github.com/jump-dev/Ipopt.jl) optimizer).
-- `jacobian=AutoForwardDiff()` : an `AbstractADType` backend for the Jacobian of the
+- `jacobian=AutoForwardDiff()`: an `AbstractADType` backend for the Jacobian of the
    nonlinear constraints, see [`DifferentiationInterface` doc](@extref DifferentiationInterface List)
-- `hessian=false` : an `AbstractADType` backend or `Bool` for the Hessian of the Lagrangian, 
+- `hessian=false`: an `AbstractADType` backend or `Bool` for the Hessian of the Lagrangian, 
    see `jacobian` above for the options. The default `false` skip it and use the
     quasi-Newton method of `optim` (see Extended Help).
+- `unit="s"`: time unit shown in plots x-label suffix (no underlying time data conversion).
 
 # Examples
 ```jldoctest
@@ -251,12 +254,13 @@ function NonLinModelDAE{NT}(
     optim_output  = JuMP.Model(DEFAULT_NLP_OPTIMIZER, add_bridges=false),
     jacobian = DEFAULT_JACDENSE,
     hessian = false,
+    unit = "s"
 ) where {NT<:Real}
     fq!, h! = get_mutating_functions_dae(NT, fq, h)
     hessian = validate_hessian(hessian, DEFAULT_NONLINDAE_HESSIAN)
     return NonLinModelDAE{NT}(
         fq!, h!, Ts, nu, nx, na, ny, nd, p, xs_0, as_0,
-        transcription, optim_state, optim_output, jacobian, hessian
+        transcription, optim_state, optim_output, jacobian, hessian, unit
     )
 end
 
@@ -271,10 +275,11 @@ function NonLinModelDAE(
     optim_output  = JuMP.Model(DEFAULT_NLP_OPTIMIZER, add_bridges=false),
     jacobian = DEFAULT_JACDENSE,
     hessian = false,
+    unit = "s"
 )
     return NonLinModelDAE{Float64}(
         fq, h, Ts, nu, nx, na, ny, nd; 
-        p, xs_0, as_0, transcription, optim_state, optim_output, jacobian, hessian
+        p, xs_0, as_0, transcription, optim_state, optim_output, jacobian, hessian, unit
     )
 end
 
@@ -933,11 +938,12 @@ function getinfo(model::NonLinModelDAE{NT}) where NT<:Real
 end
 
 function Base.show(io::IO, model::NonLinModelDAE)
+    Ts, unit = model.Ts, model.unit
     nu, nd = model.nu, model.nd
     nx, ny = model.nx, model.ny
     na = model.na
     n = maximum(ndigits.((nu, nx, ny, nd))) + 1
-    println(io, "$(nameof(typeof(model))) with a sample time Ts = $(model.Ts) s:")
+    println(io, "$(nameof(typeof(model))) with a sample time Ts = $Ts $unit:")
     println(io, "├ state optimizer: $(JuMP.solver_name(model.optim_state))")
     println(io, "├ output optimizer: $(JuMP.solver_name(model.optim_output))")
     println(io, "├ transcription: $(transcription_str(model.transcription))")
diff --git a/src/predictive_control.jl b/src/predictive_control.jl
index 40bd550db..4fd426a3f 100644
--- a/src/predictive_control.jl
+++ b/src/predictive_control.jl
@@ -30,12 +30,13 @@ include("controller/transcription.jl")
 
 function Base.show(io::IO, mpc::PredictiveController)
     estim, model = mpc.estim, mpc.estim.model
+    Ts, unit = model.Ts, model.unit
     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)
     n = maximum(ndigits.((Hp, Hc, nu, nx̂, nym, nyu, nd, other_dims...))) + 1
-    println(io, "$(nameof(typeof(mpc))) controller with a sample time Ts = $(model.Ts) s:")
+    println(io, "$(nameof(typeof(mpc))) controller with a sample time Ts = $Ts $unit:")
     println(io, "├ estimator: $(nameof(typeof(mpc.estim)))")
     println(io, "├ model: $(nameof(typeof(model)))")
     println(io, "├ optimizer: $(JuMP.solver_name(mpc.optim))")
diff --git a/src/sim_model.jl b/src/sim_model.jl
index 33dbc27d6..3d4f4d9e8 100644
--- a/src/sim_model.jl
+++ b/src/sim_model.jl
@@ -387,10 +387,11 @@ include("model/nonlinmodel.jl")
 include("model/nonlinmodeldae.jl")
 
 function Base.show(io::IO, model::SimModel)
+    Ts, unit = model.Ts, model.unit
     nu, nd = model.nu, model.nd
     nx, ny = model.nx, model.ny
     n = maximum(ndigits.((nu, nx, ny, nd))) + 1
-    println(io, "$(nameof(typeof(model))) with a sample time Ts = $(model.Ts) s:")
+    println(io, "$(nameof(typeof(model))) with a sample time Ts = $Ts $unit:")
     print_details(io, model)
     println(io, "└ dimensions:")
     println(io, "  ├$(lpad(nu, n)) manipulated inputs u")
diff --git a/src/state_estim.jl b/src/state_estim.jl
index 49744a7e9..3b38b122e 100644
--- a/src/state_estim.jl
+++ b/src/state_estim.jl
@@ -32,11 +32,12 @@ include("estimator/manual.jl")
 
 function Base.show(io::IO, estim::StateEstimator)
     model = estim.model
+    Ts, unit = model.Ts, model.unit
     nu, nd = model.nu, model.nd
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
     other_dims = get_other_dims(estim)
     n = maximum(ndigits.((nu, nx̂, nym, nyu, nd, other_dims...))) + 1
-    println(io, "$(nameof(typeof(estim))) estimator with a sample time Ts = $(model.Ts) s:")
+    println(io, "$(nameof(typeof(estim))) estimator with a sample time Ts = $Ts $unit:")
     println(io, "├ model: $(nameof(typeof(estim.model)))")
     print_details(io, estim)
     println(io, "└ dimensions:")

From 67e51baea87d53e9831b4b97283b5aaf4fe0e6b1 Mon Sep 17 00:00:00 2001
From: franckgaga <franckgaga2@gmail.com>
Date: Tue, 29 Sep 2026 16:07:49 -0400
Subject: [PATCH 2/5] test: simple construction tests with `unit` kwarg

---
 test/1_test_sim_model.jl | 26 +++++++++++++++-----------
 1 file changed, 15 insertions(+), 11 deletions(-)

diff --git a/test/1_test_sim_model.jl b/test/1_test_sim_model.jl
index 559792804..c8d6a1fb8 100644
--- a/test/1_test_sim_model.jl
+++ b/test/1_test_sim_model.jl
@@ -1,6 +1,6 @@
 @testitem "LinModel construction" setup=[SetupMPCtests] begin
     using .SetupMPCtests, ControlSystemsBase, LinearAlgebra
-    linmodel1 = LinModel(sys, Ts, i_u=1:2)
+    linmodel1 = LinModel(sys, Ts, i_u=1:2, unit="h")
     @test linmodel1.nx == 2
     @test linmodel1.nu == 2
     @test linmodel1.nd == 0
@@ -10,6 +10,7 @@
     @test linmodel1.Bd  ≈ zeros(2,0)
     @test linmodel1.C   ≈ Gss.C
     @test linmodel1.Dd  ≈ zeros(2,0)
+    @test linmodel1.unit == "h"
 
     linmodel2 = LinModel(Gss)
     setop!(linmodel2, uop=[10,50], yop=[50,30])
@@ -77,8 +78,9 @@
     linmodel11 = LinModel(Gss.A, Gss.B, I, 0, 0, Ts)
     @test linmodel11.ny == linmodel11.nx
 
-    linmodel12 = LinModel{Float32}(Gss.A, Gss.B, Gss.C, zeros(2, 0), zeros(2, 0), Ts)
+    linmodel12 = LinModel{Float32}(Gss.A, Gss.B, Gss.C, zeros(2, 0), zeros(2, 0), Ts, "h")
     @test isa(linmodel12, LinModel{Float32})
+    @test linmodel12.unit == "h"
 
     linmodel13 = LinModel(sys,Ts,i_d=[3])
     linmodel13 = setname!(linmodel13, 
@@ -164,15 +166,16 @@ end
     linmodel1 = LinModel(sys,Ts,i_u=[1,2])
     f1!(x,u,_,model) = model.A*x + model.Bu*u
     h1!(x,_,model)   = model.C*x
-    daemodel = NonLinModel(f1!,h1!,Ts,2,2,2,solver=nothing,p=linmodel1)
-    @test daemodel.nx == 2
-    @test daemodel.nu == 2
-    @test daemodel.nd == 0
-    @test daemodel.ny == 2
-    ẋ, y = daemodel.buffer.x, daemodel.buffer.y
-    daemodel.f!(ẋ, [0,0],[0,0],[1],daemodel.p)
+    nonlinmodel1 = NonLinModel(f1!,h1!,Ts,2,2,2, solver=nothing, p=linmodel1, unit="h")
+    @test nonlinmodel1.nx == 2
+    @test nonlinmodel1.nu == 2
+    @test nonlinmodel1.nd == 0
+    @test nonlinmodel1.ny == 2
+    @test nonlinmodel1.unit == "h"
+    ẋ, y = nonlinmodel1.buffer.x, nonlinmodel1.buffer.y
+    nonlinmodel1.f!(ẋ, [0,0],[0,0],[1],nonlinmodel1.p)
     @test ẋ ≈ zeros(2,)
-    daemodel.h!(y,[0,0],[1],daemodel.p)
+    nonlinmodel1.h!(y,[0,0],[1],nonlinmodel1.p)
     @test y ≈ zeros(2,)
 
     linmodel2 = LinModel(sys,Ts,i_d=[3])
@@ -451,12 +454,13 @@ end
     Ts = 1
     p = [0.5]
 
-    dae = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny; p)
+    dae = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny; p, unit="h")
     @test dae.nx == nx
     @test dae.na == na
     @test dae.nu == nu
     @test dae.nd == 0
     @test dae.ny == ny
+    @test dae.unit == "h"
     @test dae.iszero_Ha == false 
     ẋ, q, y = dae.buffer.x, dae.buffer.a, dae.buffer.y
     dae.fq!(ẋ, q, [0], [0], [0], [0], dae.p)

From 6f045d4e12c75f2ae0d9ab0fb2c4a7ecaf0b7fb7 Mon Sep 17 00:00:00 2001
From: franckgaga <franckgaga2@gmail.com>
Date: Tue, 29 Sep 2026 16:12:38 -0400
Subject: [PATCH 3/5] added: show correct time unit in the plots

---
 src/plot_sim.jl | 74 ++++++++++++++++++++++++++-----------------------
 1 file changed, 40 insertions(+), 34 deletions(-)

diff --git a/src/plot_sim.jl b/src/plot_sim.jl
index 467e19856..313dd13d3 100644
--- a/src/plot_sim.jl
+++ b/src/plot_sim.jl
@@ -389,6 +389,8 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:SimModel}) = nothing
     nd = length(indices_d)
     nx = length(indices_x)
 
+    unit = model.unit
+
     layout_mat = Matrix{Tuple{Int64, Int64}}(undef, 1, 0)
     ny > 0 && (layout_mat = [layout_mat (ny, 1)])
     nu > 0 && (layout_mat = [layout_mat (nu, 1)])
@@ -401,7 +403,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:SimModel}) = nothing
     for i in 1:ny
         i_y = indices_y[i]
         @series begin
-            i == ny && (xguide --> "Time (s)")
+            i == ny && (xguide --> "Time ($unit)")
             yguide  --> yname[i_y]
             color   --> 1
             subplot --> subplot_base + i
@@ -415,7 +417,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:SimModel}) = nothing
     for i in 1:nu
         i_u = indices_u[i]
         @series begin
-            i == nu && (xguide --> "Time (s)")
+            i == nu && (xguide --> "Time ($unit)")
             yguide     --> uname[i_u]
             color      --> 1
             subplot    --> subplot_base + i
@@ -430,7 +432,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:SimModel}) = nothing
     for i in 1:nd
         i_d = indices_d[i]
         @series begin
-            i == nd && (xguide --> "Time (s)")
+            i == nd && (xguide --> "Time ($unit)")
             yguide  --> dname[i_d]
             color   --> 1
             subplot --> subplot_base + i
@@ -444,7 +446,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:SimModel}) = nothing
     for i in 1:nx
         i_x = indices_x[i]
         @series begin
-            i == nx && (xguide --> "Time (s)")
+            i == nx && (xguide --> "Time ($unit)")
             yguide     --> xname[i_x]
             color      --> 1
             subplot    --> subplot_base + i
@@ -532,6 +534,8 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
     nx  = length(indices_x)
     nx̂  = length(indices_x̂)
     nxx̂ = length(indices_xx̂)
+
+    unit = model.unit
     
     layout_mat = Matrix{Tuple{Int64, Int64}}(undef, 1, 0)
     ny  ≠ 0 && (layout_mat = [layout_mat (ny, 1)])
@@ -547,7 +551,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
     for i in 1:ny
         i_y = indices_y[i]
         @series begin
-            i == ny && (xguide --> "Time (s)")
+            i == ny && (xguide --> "Time ($unit)")
             yguide  --> yname[i_y]
             color   --> 1
             subplot --> subplot_base + i
@@ -557,7 +561,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
         end
         if i_y in indices_ŷ
             @series begin
-                i == ny && (xguide --> "Time (s)")
+                i == ny && (xguide --> "Time ($unit)")
                 yguide  --> yname[i_y]
                 color     --> 2
                 subplot   --> subplot_base + i
@@ -574,7 +578,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
     for i in 1:nu
         i_u = indices_u[i]
         @series begin
-            i == nu && (xguide --> "Time (s)")
+            i == nu && (xguide --> "Time ($unit)")
             yguide     --> uname[i_u]
             color      --> 1
             subplot    --> subplot_base + i
@@ -589,7 +593,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
     for i in 1:nd
         i_d = indices_d[i]
         @series begin
-            i == nd && (xguide --> "Time (s)")
+            i == nd && (xguide --> "Time ($unit)")
             yguide  --> dname[i_d]
             color   --> 1
             subplot --> subplot_base + i
@@ -603,7 +607,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
     for i in 1:nx
         i_x = indices_x[i]
         @series begin
-            i == nx && (xguide --> "Time (s)")
+            i == nx && (xguide --> "Time ($unit)")
             yguide     --> xname[i_x]
             color      --> 1
             subplot    --> subplot_base + i
@@ -617,7 +621,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
     for i in 1:nx̂
         i_x̂ = indices_x̂[i]
         @series begin
-            i == nx̂ && (xguide --> "Time (s)")
+            i == nx̂ && (xguide --> "Time ($unit)")
             yguide     --> x̂name[i_x̂]
             color      --> 2
             subplot    --> subplot_base + i
@@ -630,7 +634,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
         x̂min_i, x̂max_i = X̂min[end-2*estim.nx̂+i_x̂], X̂max[end-2*estim.nx̂+i_x̂]
         if i_x̂ in indices_x̂min && !isinf(x̂min_i)
             @series begin
-                i == nx̂ && (xguide --> "Time (s)")
+                i == nx̂ && (xguide --> "Time ($unit)")
                 yguide     --> x̂name[i_x̂]
                 color      --> 4
                 subplot    --> subplot_base + i
@@ -643,7 +647,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
         end
         if i_x̂ in indices_x̂max && !isinf(x̂max_i)
             @series begin
-                i == nx̂ && (xguide --> "Time (s)")
+                i == nx̂ && (xguide --> "Time ($unit)")
                 yguide     --> x̂name[i_x̂]
                 color      --> 5
                 subplot    --> subplot_base + i
@@ -663,7 +667,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
         isplotted_x̂ = i_xx̂ ≤ size(res.X̂_data, 1)
         if isplotted_x
             @series begin
-                i == nxx̂ && (xguide --> "Time (s)")
+                i == nxx̂ && (xguide --> "Time ($unit)")
                 yguide     --> xx̂name[i_xx̂]
                 color      --> 1
                 subplot    --> subplot_base + i
@@ -674,7 +678,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
         end
         if isplotted_x̂
             @series begin
-                i == nxx̂ && (xguide --> "Time (s)")
+                i == nxx̂ && (xguide --> "Time ($unit)")
                 yguide     --> xx̂name[i_xx̂]
                 color      --> 2
                 subplot    --> subplot_base + i
@@ -687,7 +691,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
             x̂min_i, x̂max_i = X̂min[end-2*estim.nx̂+i_xx̂], X̂max[end-2*estim.nx̂+i_xx̂]
             if i_xx̂ in indices_x̂min && !isinf(x̂min_i)
                 @series begin
-                    i == nxx̂ && (xguide --> "Time (s)")
+                    i == nxx̂ && (xguide --> "Time ($unit)")
                     yguide     --> xx̂name[i_xx̂]
                     color      --> 4
                     subplot    --> subplot_base + i
@@ -700,7 +704,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:StateEstimator}) = nothing
             end
             if i_xx̂ in indices_x̂max && !isinf(x̂max_i)
                 @series begin
-                    i == nxx̂ && (xguide --> "Time (s)")
+                    i == nxx̂ && (xguide --> "Time ($unit)")
                     yguide     --> xx̂name[i_xx̂]
                     color      --> 5
                     subplot    --> subplot_base + i
@@ -815,6 +819,8 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
     nx̂  = length(indices_x̂)
     nxx̂ = length(indices_xx̂)
 
+    unit = model.unit
+
     layout_mat = Matrix{Tuple{Int64, Int64}}(undef, 1, 0)
     ny  ≠ 0 && (layout_mat = [layout_mat (ny, 1)])
     nu  ≠ 0 && (layout_mat = [layout_mat (nu, 1)])
@@ -829,7 +835,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
     for i in 1:ny
         i_y = indices_y[i]
         @series begin
-            i == ny && (xguide --> "Time (s)")
+            i == ny && (xguide --> "Time ($unit)")
             yguide  --> yname[i_y]
             color   --> 1
             subplot --> subplot_base + i
@@ -839,7 +845,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         end
         if i_y in indices_ŷ
             @series begin
-                i == ny && (xguide --> "Time (s)")
+                i == ny && (xguide --> "Time ($unit)")
                 yguide  --> yname[i_y]
                 color     --> 2
                 subplot   --> subplot_base + i
@@ -853,7 +859,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         M_Hp_i = mpc.weights.M_Hp[i_y, i_y]
         if i_y in indices_ry && !iszero(M_Hp_i)
             @series begin
-                i == ny && (xguide --> "Time (s)")
+                i == ny && (xguide --> "Time ($unit)")
                 yguide    --> yname[i_y]
                 color     --> 3
                 subplot   --> subplot_base + i
@@ -867,7 +873,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         ymin_i, ymax_i = Ymin[i_y], Ymax[i_y]
         if i_y in indices_ymin && !isinf(ymin_i)
             @series begin
-                i == ny && (xguide --> "Time (s)")
+                i == ny && (xguide --> "Time ($unit)")
                 yguide    --> yname[i_y]
                 color     --> 4
                 subplot   --> subplot_base + i
@@ -880,7 +886,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         end
         if i_y in indices_ymax && !isinf(ymax_i)
             @series begin
-                i == ny && (xguide --> "Time (s)")
+                i == ny && (xguide --> "Time ($unit)")
                 yguide    --> yname[i_y]
                 color     --> 5
                 subplot   --> subplot_base + i
@@ -897,7 +903,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
     for i in 1:nu
         i_u = indices_u[i]
         @series begin
-            i == nu && (xguide --> "Time (s)")
+            i == nu && (xguide --> "Time ($unit)")
             yguide     --> uname[i_u]
             color      --> 1
             subplot    --> subplot_base + i
@@ -909,7 +915,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         L_Hp_i = mpc.weights.L_Hp[i_u, i_u]
         if i_u in indices_ru && !iszero(L_Hp_i)
             @series begin
-                i == nu && (xguide --> "Time (s)")
+                i == nu && (xguide --> "Time ($unit)")
                 yguide    --> uname[i_u]
                 color     --> 3
                 subplot   --> subplot_base + i
@@ -923,7 +929,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         umin_i, umax_i = Umin[i_u], Umax[i_u]
         if i_u in indices_umin && !isinf(umin_i)
             @series begin
-                i == nu && (xguide --> "Time (s)")
+                i == nu && (xguide --> "Time ($unit)")
                 yguide    --> uname[i_u]
                 color     --> 4
                 subplot   --> subplot_base + i
@@ -936,7 +942,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         end
         if i_u in indices_umax && !isinf(umax_i)
             @series begin
-                i == nu && (xguide --> "Time (s)")
+                i == nu && (xguide --> "Time ($unit)")
                 yguide    --> uname[i_u]
                 color     --> 5
                 subplot   --> subplot_base + i
@@ -953,7 +959,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
     for i in 1:nd
         i_d = indices_d[i]
         @series begin
-            i == nd && (xguide --> "Time (s)")
+            i == nd && (xguide --> "Time ($unit)")
             yguide  --> dname[i_d]
             color   --> 1
             subplot --> subplot_base + i
@@ -967,7 +973,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
     for i in 1:nx
         i_x = indices_x[i]
         @series begin
-            i == nx && (xguide --> "Time (s)")
+            i == nx && (xguide --> "Time ($unit)")
             yguide     --> xname[i_x]
             color      --> 1
             subplot    --> subplot_base + i
@@ -981,7 +987,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
     for i in 1:nx̂
         i_x̂ = indices_x̂[i]
         @series begin
-            i == nx̂ && (xguide --> "Time (s)")
+            i == nx̂ && (xguide --> "Time ($unit)")
             yguide     --> x̂name[i_x̂]
             color      --> 2
             subplot    --> subplot_base + i
@@ -994,7 +1000,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         x̂min_i, x̂max_i = X̂min[end-2*estim.nx̂+i_x̂], X̂max[end-2*estim.nx̂+i_x̂]
         if i_x̂ in indices_x̂min && !isinf(x̂min_i)
             @series begin
-                i == nx̂ && (xguide --> "Time (s)")
+                i == nx̂ && (xguide --> "Time ($unit)")
                 yguide     --> x̂name[i_x̂]
                 color      --> 4
                 subplot    --> subplot_base + i
@@ -1007,7 +1013,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         end
         if i_x̂ in indices_x̂max && !isinf(x̂max_i)
             @series begin
-                i == nx̂ && (xguide --> "Time (s)")
+                i == nx̂ && (xguide --> "Time ($unit)")
                 yguide     --> x̂name[i_x̂]
                 color      --> 5
                 subplot    --> subplot_base + i
@@ -1027,7 +1033,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         isplotted_x̂ = i_xx̂ ≤ size(res.X̂_data, 1)
         if isplotted_x
             @series begin
-                i == nxx̂ && (xguide --> "Time (s)")
+                i == nxx̂ && (xguide --> "Time ($unit)")
                 yguide     --> xx̂name[i_xx̂]
                 color      --> 1
                 subplot    --> subplot_base + i
@@ -1038,7 +1044,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
         end
         if isplotted_x̂
             @series begin
-                i == nxx̂ && (xguide --> "Time (s)")
+                i == nxx̂ && (xguide --> "Time ($unit)")
                 yguide     --> xx̂name[i_xx̂]
                 color      --> 2
                 subplot    --> subplot_base + i
@@ -1051,7 +1057,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
             x̂min_i, x̂max_i = X̂min[end-2*estim.nx̂+i_xx̂], X̂max[end-2*estim.nx̂+i_xx̂]
             if i_xx̂ in indices_x̂min && !isinf(x̂min_i)
                 @series begin
-                    i == nxx̂ && (xguide --> "Time (s)")
+                    i == nxx̂ && (xguide --> "Time ($unit)")
                     yguide     --> xx̂name[i_xx̂]
                     color      --> 4
                     subplot    --> subplot_base + i
@@ -1064,7 +1070,7 @@ plot_recipe(::Nothing, ::SimResult{<:Real, <:PredictiveController}) = nothing
             end
             if i_xx̂ in indices_x̂max && !isinf(x̂max_i)
                 @series begin
-                    i == nxx̂ && (xguide --> "Time (s)")
+                    i == nxx̂ && (xguide --> "Time ($unit)")
                     yguide     --> xx̂name[i_xx̂]
                     color      --> 5
                     subplot    --> subplot_base + i

From eb61bec4862acc7dae8976f011b94b9a251460f1 Mon Sep 17 00:00:00 2001
From: franckgaga <franckgaga2@gmail.com>
Date: Tue, 29 Sep 2026 16:23:28 -0400
Subject: [PATCH 4/5] doc: use the new `unit` argument in the pH tutorial

---
 docs/src/manual/nonlinmpc2.md | 21 +++++++++++++--------
 1 file changed, 13 insertions(+), 8 deletions(-)

diff --git a/docs/src/manual/nonlinmpc2.md b/docs/src/manual/nonlinmpc2.md
index f62776e99..e82ed1acd 100644
--- a/docs/src/manual/nonlinmpc2.md
+++ b/docs/src/manual/nonlinmpc2.md
@@ -184,7 +184,8 @@ Kw = 1.0e-14    # water dissociation constant [mol^2/L^2]
 Ka = 1.75e-5    # acid dissociation constant [mol/L]
 V = 1000.0      # reactor volume [L]
 
-Ts = 0.5        # Sample time [h]
+Ts = 0.5        # sample time
+unit = "h"      # time unit shown on plot x-label suffixes
 nu, nx, na, ny, nd = 1, 2, 1, 1, 1
 p = [c_Ain, c_Bin, Kw, Ka, V]
 
@@ -194,7 +195,10 @@ vx, vy = [raw"$c_A$ (mol/L)", raw"$c_B$ (mol/L)"], [raw"$\mathrm{pH}$"]
 transcription = TrapezoidalCollocation()
 xs_0, as_0 = [0.025, 0.025], [7]
 
-plant = NonLinModelDAE(fq!, h!, Ts, nu, nx, na, ny, nd; p=p, xs_0, as_0, transcription)
+plant = NonLinModelDAE(
+    fq!, h!, Ts, nu, nx, na, ny, nd; 
+    p=p, xs_0, as_0, transcription, unit
+)
 plant = setname!(plant, u=vu, x=vx, y=vy, d=vd)
 ```
 
@@ -236,12 +240,11 @@ N = 81
 res = simDAE(plant, N; x_0)
 ```
 
-We plot the results by modifying the x-axis label to substitute the default time units
-to hours:
+We can now plot the result:
 
 ```@example 1
 using Plots
-plot(res, plotu=true, plotd=true, xlabel="Time (h)")
+plot(res, plotu=true, plotd=true)
 savefig("plot1_DAEpH.svg"); nothing # hide
 ```
 
@@ -280,7 +283,10 @@ p̂ = [c_Bin, Kw, Ka, V]
 nx̂ = nx + 1
 vx̂ = [vx; raw"$c_{Ain}$ (mol/L)"]
 x̂s_0 = [xs_0; 0.1]
-model = NonLinModelDAE(f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; p=p̂, xs_0=x̂s_0, as_0, transcription)
+model = NonLinModelDAE(
+    f̂q!, ĥ!, Ts, nu, nx̂, na, ny, nd; 
+    p=p̂, xs_0=x̂s_0, as_0, transcription, unit
+)
 model = setname!(model, u=vu, x=vx̂, y=vy, d=vd)
 ```
 
@@ -346,8 +352,7 @@ function simMHE(mhe, plant, N; x_0, x̂_0)
 end
 x̂_0 = [0.025, 0.025, c_Ain]
 res = simMHE(mhe, plant, N; x_0, x̂_0)
-p = plot(res, plotd=false, plotu=false, plotxwithx̂=true, plotx̂min=false, xlabel="Time (h)")
-xlabel!(p[2], ""); xlabel!(p[3], "") # remove xlabel on c_A and c_B plots
+p = plot(res, plotd=false, plotu=false, plotxwithx̂=true, plotx̂min=false)
 savefig(p, "plot2_DAEpH.svg"); nothing # hide
 ```
 

From 1b11030e057366e56a57f5a73e7fdd492ee05159 Mon Sep 17 00:00:00 2001
From: franckgaga <franckgaga2@gmail.com>
Date: Tue, 29 Sep 2026 18:42:58 -0400
Subject: [PATCH 5/5] doc: minor correction

---
 src/model/linmodel.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/model/linmodel.jl b/src/model/linmodel.jl
index a38a3e518..8fff8feb6 100644
--- a/src/model/linmodel.jl
+++ b/src/model/linmodel.jl
@@ -208,7 +208,7 @@ end
 
 
 @doc raw"""
-    LinModel(sys::TransferFunction[, Ts]; unit="s", i_u=1:size(sys,2), i_d=Int[])
+    LinModel(sys::TransferFunction[, Ts]; i_u=1:size(sys,2), i_d=Int[], unit="s")
 
 Convert to minimal realization state-space when `sys` is a transfer function.
 
@@ -235,7 +235,7 @@ end
 
 
 """
-    LinModel(sys::DelayLtiSystem, Ts; unit="s", i_u=1:size(sys,2), i_d=Int[])
+    LinModel(sys::DelayLtiSystem, Ts; i_u=1:size(sys,2), i_d=Int[], unit="s")
 
 Discretize with zero-order hold when `sys` is a continuous system with delays.