added: unit="s" keyword argument in all SimModel constructors

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
 ```
 

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,16 +63,17 @@ struct LinModel{NT<:Real} <: SimModelODE{NT}
             nu, nx, ny, nd, nk̄,
             uop, yop, dop, xop, fop,
             uname, yname, dname, xname,
+            unit,
             xs_0,
             buffer
         )
     end
 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]; i_u=1:size(sys,2), i_d=Int[], unit="s")
 
 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; i_u=1:size(sys,2), i_d=Int[], unit="s")
 
 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."))

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,23 +206,23 @@ 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)
     linfunc! = get_linearization_func(
         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."

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))")