breaking change: continuous dynamics by default for NonLinModel

Author: franckgaga

docs/src/manual/nonlinmpc.md

@@ -53,13 +53,13 @@ const par = (9.8, 0.4, 1.2, 0.3)
 f(x, u, _ ) = pendulum(par, x, u)
 h(x, _ )    = [180/π*x[1]]  # [°]
 Ts, nu, nx, ny = 0.1, 1, 2, 1
-solver = RungeKutta()
-model = NonLinModel(f, h, Ts, nu, nx, ny; solver)
+model = NonLinModel(f, h, Ts, nu, nx, ny)
 ```
 
 The output function ``\mathbf{h}`` converts the ``θ`` angle to degrees. Note that special
 characters like ``θ`` can be typed in the Julia REPL or VS Code by typing `\theta` and
 pressing the `<TAB>` key. The tuple `par` is constant here to improve the [performance](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-untyped-global-variables).
+Note that a 4th order [`RungeKutta`](@ref) differential equation solver is used by default.
 It is good practice to first simulate `model` using [`sim!`](@ref) as a quick sanity check:
 
 ```@example 1
@@ -92,7 +92,7 @@ estimator tuning is tested on a plant with a 25 % larger friction coefficient ``
 ```@example 1
 const par_plant = (par[1], par[2], 1.25*par[3], par[4])
 f_plant(x, u, _ ) = pendulum(par_plant, x, u)
-plant = NonLinModel(f_plant, h, Ts, nu, nx, ny; solver)
+plant = NonLinModel(f_plant, h, Ts, nu, nx, ny)
 res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
 plot(res, plotu=false, plotxwithx̂=true)
 savefig(ans, "plot2_NonLinMPC.svg"); nothing # hide
@@ -170,8 +170,8 @@ Kalman Filter similar to the previous one (``\mathbf{y^m} = θ`` and ``\mathbf{y
 ```@example 1
 h2(x, _ ) = [180/π*x[1], x[2]]
 nu, nx, ny = 1, 2, 2
-model2 = NonLinModel(f      , h2, Ts, nu, nx, ny; solver)
-plant2 = NonLinModel(f_plant, h2, Ts, nu, nx, ny; solver)
+model2 = NonLinModel(f      , h2, Ts, nu, nx, ny)
+plant2 = NonLinModel(f_plant, h2, Ts, nu, nx, ny)
 estim2 = UnscentedKalmanFilter(model2; σQ, σR, nint_u, σQint_u, i_ym=[1])
 ```
 
@@ -184,7 +184,7 @@ function JE(UE, ŶE, _ )
     τ, ω = UE[1:end-1], ŶE[2:2:end-1]
     return Ts*sum(τ.*ω)
 end
-empc = NonLinMPC(estim2, Hp=20, Hc=2, Mwt=[0.5, 0], Nwt=[2.5], Ewt=4e3, JE=JE)
+empc = NonLinMPC(estim2, Hp=20, Hc=2, Mwt=[0.5, 0], Nwt=[2.5], Ewt=3e9, JE=JE)
 empc = setconstraint!(empc, umin=[-1.5], umax=[+1.5])
 ```
 

src/controller/nonlinmpc.jl

@@ -145,7 +145,7 @@ This method uses the default state estimator :
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1);
+julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> mpc = NonLinMPC(model, Hp=20, Hc=1, Cwt=1e6)
 NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedKalmanFilter estimator and:
@@ -166,8 +166,8 @@ NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedK
 
     The optimization relies on [`JuMP`](https://github.com/jump-dev/JuMP.jl) automatic 
     differentiation (AD) to compute the objective and constraint derivatives. Optimizers 
-    generally benefit from exact derivatives like AD. However, the [`NonLinModel`](@ref) `f` 
-    and `h` functions must be compatible with this feature. See [Automatic differentiation](https://jump.dev/JuMP.jl/stable/manual/nlp/#Automatic-differentiation)
+    generally benefit from exact derivatives like AD. However, the [`NonLinModel`](@ref) 
+    state-space functions must be compatible with this feature. See [Automatic differentiation](https://jump.dev/JuMP.jl/stable/manual/nlp/#Automatic-differentiation)
     for common mistakes when writing these functions.
 
     Note that if `Cwt≠Inf`, the attribute `nlp_scaling_max_gradient` of `Ipopt` is set to 
@@ -221,7 +221,7 @@ Use custom state estimator `estim` to construct `NonLinMPC`.
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1);
+julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> estim = UnscentedKalmanFilter(model, σQint_ym=[0.05]);
 

src/estimator/kalman.jl

@@ -434,7 +434,7 @@ represents the measured outputs of ``\mathbf{ĥ}`` function (and unmeasured one
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1);
+julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> estim = UnscentedKalmanFilter(model, σR=[1], nint_ym=[2], σP0int_ym=[1, 1])
 UnscentedKalmanFilter estimator with a sample time Ts = 10.0 s, NonLinModel and:
@@ -685,7 +685,7 @@ automatic differentiation.
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.2x+u, (x,_)->-3x, 5.0, 1, 1, 1);
+julia> model = NonLinModel((x,u,_)->0.2x+u, (x,_)->-3x, 5.0, 1, 1, 1, solver=nothing);
 
 julia> estim = ExtendedKalmanFilter(model, σQ=[2], σQint_ym=[2], σP0=[0.1], σP0int_ym=[0.1])
 ExtendedKalmanFilter estimator with a sample time Ts = 5.0 s, NonLinModel and:

src/estimator/mhe/construct.jl

@@ -206,7 +206,7 @@ matrix ``\mathbf{P̂}_{k-N_k}(k-N_k+1)`` is estimated with an [`ExtendedKalmanFi
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1);
+julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> estim = MovingHorizonEstimator(model, He=5, σR=[1], σP0=[0.01])
 MovingHorizonEstimator estimator with a sample time Ts = 10.0 s, Ipopt optimizer, NonLinModel and:

src/model/linearization.jl

@@ -22,7 +22,7 @@ Jacobians of ``\mathbf{f}`` and ``\mathbf{h}`` functions are automatically compu
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->x.^3 + u, (x,_)->x, 0.1, 1, 1, 1);
+julia> model = NonLinModel((x,u,_)->x.^3 + u, (x,_)->x, 0.1, 1, 1, 1, solver=nothing);
 
 julia> linmodel = linearize(model, x=[10.0], u=[0.0]); 
 

src/model/linmodel.jl

@@ -35,26 +35,35 @@ end
 
 Construct a linear model from state-space model `sys` with sampling time `Ts` in second.
 
-`Ts` can be omitted when `sys` is discrete-time. Its state-space matrices are:
+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):
 ```math
 \begin{aligned}
-    \mathbf{x}(k+1) &= \mathbf{A x}(k) + \mathbf{B z}(k) \\
-    \mathbf{y}(k)   &= \mathbf{C x}(k) + \mathbf{D z}(k)
+    \mathbf{ẋ}(t) &= \mathbf{A x}(t) + \mathbf{B z}(t) \\
+    \mathbf{y}(t) &= \mathbf{C x}(t) + \mathbf{D z}(t)
 \end{aligned}
 ```
 with the state ``\mathbf{x}`` and output ``\mathbf{y}`` vectors. The ``\mathbf{z}`` vector 
 comprises the manipulated inputs ``\mathbf{u}`` and measured disturbances ``\mathbf{d}``, 
 in any order. `i_u` provides the indices of ``\mathbf{z}`` that are manipulated, and `i_d`, 
-the measured disturbances. See Extended Help if `sys` is continuous-time, or discrete-time
-with `Ts ≠ sys.Ts`.
+the measured disturbances. The constructor automatically discretizes continuous systems,
+resamples discrete ones if `Ts ≠ sys.Ts`, compute a new realization if not minimal, and 
+separate the ``\mathbf{z}`` terms in two parts (details in Extended Help). The rest of the
+documentation assumes discrete dynamics since all systems end up in this form.
 
-See also [`ss`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.ss),
-[`tf`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.tf).
+See also [`ss`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.ss)
 
 # Examples
 ```jldoctest
-julia> model = LinModel(ss(0.4, 0.2, 0.3, 0, 0.1))
-Discrete-time linear model with a sample time Ts = 0.1 s and:
+julia> model1 = LinModel(ss(-0.1, 1.0, 1.0, 0), 2.0) # continuous-time StateSpace
+LinModel with a sample time Ts = 2.0 s and:
+ 1 manipulated inputs u
+ 1 states x
+ 1 outputs y
+ 0 measured disturbances d
+
+julia> model2 = LinModel(ss(0.4, 0.2, 0.3, 0, 0.1)) # discrete-time StateSpace
+LinModel with a sample time Ts = 0.1 s and:
  1 manipulated inputs u
  1 states x
  1 outputs y
@@ -63,9 +72,9 @@ Discrete-time linear model with a sample time Ts = 0.1 s and:
 
 # Extended Help
 !!! details "Extended Help"
-    State-space matrices are similar if `sys` is continuous (replace ``\mathbf{x}(k+1)``
-    with ``\mathbf{ẋ}(t)`` and ``k`` with ``t`` on the LHS). In such a case, it's 
-    discretized with [`c2d`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.c2d)
+    State-space equations are similar if `sys` is discrete-time (replace ``\mathbf{ẋ}(t)``
+    with ``\mathbf{x}(k+1)`` and ``k`` with ``t`` on the LHS). Continuous dynamics are 
+    internally discretized using [`c2d`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.c2d)
     and `:zoh` for manipulated inputs, and `:tustin`, for measured disturbances. Lastly, if 
     `sys` is discrete and the provided argument `Ts ≠ sys.Ts`, the system is resampled by
     using the aforementioned discretization methods.
@@ -124,8 +133,8 @@ function LinModel(
             sysd_dis = sysd
         end     
     end
-    # minreal to merge common poles if possible and ensure observability
-    sys_dis = minreal([sysu_dis sysd_dis])
+    # minreal to merge common poles if possible and ensure controllability/observability:
+    sys_dis = minreal([sysu_dis sysd_dis]) # same realization if already minimal
     nu  = length(i_u)
     A   = sys_dis.A
     Bu  = sys_dis.B[:,1:nu]
@@ -144,10 +153,12 @@ Convert to minimal realization state-space when `sys` is a transfer function.
 `sys` is equal to ``\frac{\mathbf{y}(s)}{\mathbf{z}(s)}`` for continuous-time, and 
 ``\frac{\mathbf{y}(z)}{\mathbf{z}(z)}``, for discrete-time.
 
+See also [`tf`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.tf)
+
 # Examples
 ```jldoctest
 julia> model = LinModel([tf(3, [30, 1]) tf(-2, [5, 1])], 0.5, i_d=[2])
-Discrete-time linear model with a sample time Ts = 0.5 s and:
+LinModel with a sample time Ts = 0.5 s and:
  1 manipulated inputs u
  2 states x
  1 outputs y
@@ -245,6 +256,4 @@ function h!(y, model::LinModel, x, d)
     mul!(y, model.C,  x, 1, 1)
     mul!(y, model.Dd, d, 1, 1)
     return nothing
-end
-
-typestr(model::LinModel) = "linear"
+end

src/model/nonlinmodel.jl

@@ -24,24 +24,25 @@ struct NonLinModel{NT<:Real, F<:Function, H<:Function, DS<:DiffSolver} <: SimMod
 end
 
 @doc raw"""
-    NonLinModel{NT}(f::Function,  h::Function,  Ts, nu, nx, ny, nd=0; solver=nothing)
-    NonLinModel{NT}(f!::Function, h!::Function, Ts, nu, nx, ny, nd=0; solver=nothing)
+    NonLinModel{NT}(f::Function,  h::Function,  Ts, nu, nx, ny, nd=0; solver=RungeKutta(4))
+    NonLinModel{NT}(f!::Function, h!::Function, Ts, nu, nx, ny, nd=0; solver=RungeKutta(4))
 
 Construct a nonlinear model from state-space functions `f`/`f!` and `h`/`h!`.
 
-Default arguments assume discrete-time dynamics, in which the state update ``\mathbf{f}``
-and output ``\mathbf{h}`` functions are defined as :
+Both continuous and discrete-time models are supported. The default arguments assume 
+continuous dynamics. Use `solver=nothing` for the discrete case (see Extended Help). The
+functions are defined as:
 ```math
-    \begin{aligned}
-    \mathbf{x}(k+1) &= \mathbf{f}\Big( \mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k) \Big) \\
-    \mathbf{y}(k)   &= \mathbf{h}\Big( \mathbf{x}(k), \mathbf{d}(k) \Big)
-    \end{aligned}
+\begin{aligned}
+    \mathbf{ẋ}(t) &= \mathbf{f}\Big( \mathbf{x}(t), \mathbf{u}(t), \mathbf{d}(t) \Big) \\
+    \mathbf{y}(t) &= \mathbf{h}\Big( \mathbf{x}(t), \mathbf{d}(t) \Big)
+\end{aligned}
 ```
-Denoting ``\mathbf{x}(k+1)`` as `xnext`, they can be implemented in two possible ways:
+They can be implemented in two possible ways:
 
-- Non-mutating functions (out-of-place): they must be defined as `f(x, u, d) -> xnext` and
+- Non-mutating functions (out-of-place): defined them as `f(x, u, d) -> ẋ` and
   `h(x, d) -> y`. This syntax is simple and intuitive but it allocates more memory.
-- Mutating functions (in-place): they must be defined as `f!(xnext, x, u, d) -> nothing` and
+- Mutating functions (in-place): defined them as `f!(ẋ, x, u, d) -> nothing` and
   `h!(y, x, d) -> nothing`. This syntax reduces the allocations and potentially the 
   computational burden as well.
 
@@ -50,10 +51,11 @@ manipulated inputs, states, outputs and measured disturbances.
 
 !!! tip
     Replace the `d` argument with `_` if `nd = 0` (see Examples below).
-
-Specifying a differential equation solver with the the `solver` keyword argument implies a
-continuous-time model (see Extended Help for details). The optional parameter `NT`
-explicitly set the number type of vectors (default to `Float64`).
+    
+A 4th order [`RungeKutta`](@ref) solver discretizes the differential equations by default. 
+The rest of the documentation assumes discrete dynamics since all models end up in this 
+form. The optional parameter `NT` explicitly set the number type of vectors (default to 
+`Float64`).
 
 !!! warning
     The two functions must be in pure Julia to use the model in [`NonLinMPC`](@ref),
@@ -63,19 +65,23 @@ See also [`LinModel`](@ref).
 
 # Examples
 ```jldoctest
-julia> model1 = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1)
-Discrete-time nonlinear model with a sample time Ts = 10.0 s, empty solver and:
+julia> f!(ẋ, x, u, _ ) = (ẋ .= -0.2x .+ u; nothing);
+
+julia> h!(y, x, _ ) = (y .= 0.1x; nothing);
+
+julia> model1 = NonLinModel(f!, h!, 5.0, 1, 1, 1)               # continuous dynamics
+NonLinModel with a sample time Ts = 5.0 s, RungeKutta solver and:
  1 manipulated inputs u
  1 states x
  1 outputs y
  0 measured disturbances d
 
-julia> f!(ẋ, x, u, _ ) = (ẋ .= -0.1x .+ u; nothing);
+julia> f(x, u, _ ) = 0.1x + u;
 
-julia> h!(y, x, _ ) = (y .= 2x; nothing);
+julia> h(x, _ ) = 2x;
 
-julia> model2 = NonLinModel(f!, h!, 5.0, 1, 1, 1, solver=RungeKutta())
-Discrete-time nonlinear model with a sample time Ts = 5.0 s, RungeKutta solver and:
+julia> model2 = NonLinModel(f, h, 2.0, 1, 1, 1, solver=nothing) # discrete dynamics
+NonLinModel with a sample time Ts = 2.0 s, empty solver and:
  1 manipulated inputs u
  1 states x
  1 outputs y
@@ -84,17 +90,22 @@ Discrete-time nonlinear model with a sample time Ts = 5.0 s, RungeKutta solver a
 
 # Extended Help
 !!! details "Extended Help"
-    State-space equations are similar for continuous-time models (replace ``\mathbf{x}(k+1)``
-    with ``\mathbf{ẋ}(t)`` and ``k`` with ``t`` on the LHS), with also two possible 
-    implementations (second one to reduce the allocations):
-
-    - Non-mutating functions (out-of-place): they must be defined as `f(x, u, d) -> ẋ` and
-      `h(x, d) -> y`.
-    - Mutating functions (in-place): they must be defined as `f!(ẋ, x, u, d) -> nothing` and
-      `h!(y, x, d) -> nothing`. 
+    The state-space functions are similar for discrete dynamics:
+    ```math
+    \begin{aligned}
+        \mathbf{x}(k+1) &= \mathbf{f}\Big( \mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k) \Big) \\
+        \mathbf{y}(k)   &= \mathbf{h}\Big( \mathbf{x}(k), \mathbf{d}(k) \Big)
+    \end{aligned}
+    ```
+    with two possible implementations as well:
+
+    - Non-mutating functions: defined them as `f(x, u, d) -> xnext` and `h(x, d) -> y`.
+    - Mutating functions: defined them as `f!(xnext, x, u, d) -> nothing` and
+      `h!(y, x, d) -> nothing`.
 """
 function NonLinModel{NT}(
-    f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0; solver=nothing
+    f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0; 
+    solver=RungeKutta(4)
 ) where {NT<:Real}
     isnothing(solver) && (solver=EmptySolver())
     ismutating_f = validate_f(NT, f)
@@ -107,7 +118,8 @@ function NonLinModel{NT}(
 end
 
 function NonLinModel(
-    f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0; solver=nothing
+    f::Function, h::Function, Ts::Real, nu::Int, nx::Int, ny::Int, nd::Int=0; 
+    solver=RungeKutta(4)
 )
     return NonLinModel{Float64}(f, h, Ts, nu, nx, ny, nd; solver)
 end
@@ -158,6 +170,5 @@ f!(xnext, model::NonLinModel, x, u, d) = model.f!(xnext, x, u, d)
 "Call `h!(y, x, d)` with `model.h` method for [`NonLinModel`](@ref)."
 h!(y, model::NonLinModel, x, d) = model.h!(y, x, d)
 
-typestr(model::NonLinModel) = "nonlinear"
 detailstr(model::NonLinModel) = ", $(typeof(model.solver).name.name) solver"
 detailstr(::NonLinModel{<:Real, <:Function, <:Function, <:EmptySolver}) = ", empty solver"

src/model/solver.jl

@@ -24,11 +24,11 @@ struct RungeKutta <: DiffSolver
 end
 
 """
-    RungeKutta(order::Int=4; supersample::Int=1)
+    RungeKutta(order=4; supersample=1)
 
 Create a Runge-Kutta solver with optional super-sampling.
 
-Only the fourth order Runge-Kutta is supported for now. The keyword argument `supersample`
+Only the 4th order Runge-Kutta is supported for now. The keyword argument `supersample`
 provides the number of internal steps (default to 1 step).
 """
 function RungeKutta(order::Int=4; supersample::Int=1)

src/plot_sim.jl

@@ -101,7 +101,7 @@ on them (see Examples below). Note that the method mutates `plant` internal stat
 
 # Examples
 ```julia-repl
-julia> plant = NonLinModel((x,u,d)->0.1x+u+d, (x,_)->2x, 10.0, 1, 1, 1, 1);
+julia> plant = NonLinModel((x,u,d)->0.1x+u+d, (x,_)->2x, 10.0, 1, 1, 1, 1, solver=nothing);
 
 julia> res = sim!(plant, 15, [0], [0], x0=[1])
 Simulation results of NonLinModel with 15 time steps.

src/precompile.jl

@@ -67,7 +67,7 @@ sim!(exmpc, 3, [55, 30])
 f(x,u,_) = model.A*x + model.Bu*u
 h(x,_) = model.C*x
 
-nlmodel = setop!(NonLinModel(f, h, Ts, 2, 2, 2), uop=[10, 10], yop=[50, 30])
+nlmodel = setop!(NonLinModel(f, h, Ts, 2, 2, 2, solver=nothing), uop=[10, 10], yop=[50, 30])
 y = nlmodel()
 nmpc_im = setconstraint!(NonLinMPC(InternalModel(nlmodel), Hp=10, Cwt=Inf), ymin=[45, -Inf])
 initstate!(nmpc_im, nlmodel.uop, y)