doc : add docstring to construct LinModel with a specific state-space representation

Author: franckgaga

docs/src/internals/predictive_control.md

@@ -9,7 +9,7 @@ The prediction methodology of this module is mainly based on Maciejowski textboo
 [^1]: Maciejowski, J. 2000, "Predictive control : with constraints", 1st ed., Prentice Hall,
      ISBN 978-0201398236.
 
-## Controller Initialization
+## Controller Construction
 
 ```@docs
 ModelPredictiveControl.init_predmat

docs/src/internals/state_estim.md

@@ -4,7 +4,7 @@
 Pages = ["state_estim.md"]
 ```
 
-## Estimator Initialization
+## Estimator Construction
 
 ```@docs
 ModelPredictiveControl.init_estimstoch
@@ -27,13 +27,22 @@ ModelPredictiveControl.ĥ
 ModelPredictiveControl.remove_op!
 ```
 
-## Update Estimate
+## Init Estimate
 
 !!! info
     All these methods assume that the operating points are already removed in `u`, `ym`
     and `d` arguments. Strickly speaking, the arguments should be called `u0`, `ym0` and
     `d0`, following [`setop!`](@ref) notation. The `0` is dropped to simplify the notation.
 
+```@docs
+ModelPredictiveControl.init_estimate!
+```
+
+## Update Estimate
+
+!!! info
+    Same as above: the arguments should be called `u0`, `ym0` and `d0`, strickly speaking.
+
 ```@docs
 ModelPredictiveControl.update_estimate!
 ```

src/estimator/internal_model.jl

@@ -214,20 +214,11 @@ function update_estimate!(estim::InternalModel, u, ym, d=Float64[])
     return x̂d
 end
 
-@doc raw"""
-    initstate!(estim::InternalModel, u, ym, d=Float64[]) -> x̂d
-
-Init `estim.x̂d` / `x̂s` states from current inputs `u`, meas. outputs `ym` and disturb. `d`.
-
-The deterministic state `estim.x̂d` initialization method is identical to 
-[`initstate!(::StateEstimator)`](@ref). The stochastic states `estim.x̂s` are init at 0. 
-"""
-function initstate!(estim::InternalModel, u, ym, d=Float64[])
-    model = estim.model
-    init_estimate!(estim, model, u, ym, d)
+"Init the stochastic states `estim.x̂s` of [`InternalModel`](@ref) at 0. " 
+function initstate_post!(estim::InternalModel)
     # TODO: best method to initialize internal model stochastic states ? not sure...
     estim.x̂s[:] = zeros(estim.nxs)
-    return estim.x̂d
+    return nothing
 end
 
 @doc raw"""

src/estimator/kalman.jl

@@ -725,7 +725,7 @@ function update_estimate!(estim::ExtendedKalmanFilter, u, ym, d=Float64[])
 end
 
 "Initialize the covariance estimate `P̂` for the time-varying Kalman Filters" 
-function initstate_cov!(
+function initstate_post!(
     estim::Union{KalmanFilter, UnscentedKalmanFilter, ExtendedKalmanFilter}
 ) 
     estim.P̂.data[:] = estim.P̂0

src/model/linmodel.jl

@@ -68,15 +68,17 @@ and `:zoh` for manipulated inputs, and `:tustin`, for measured disturbances. Las
 the aforementioned discretization methods.
 
 Note that the constructor transforms the system to its minimal realization using [`minreal`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.minreal)
-to favor observability. As a consequence, the final state-space representation may be 
-different from the one provided in `sys`. It is also converted into a more practical form
-(``\mathbf{D_u=0}`` because of the zero-order hold):
+to favor controllability and observability. As a consequence, the final state-space
+representation may be different from the one provided in `sys`. It is also converted into a
+more practical form (``\mathbf{D_u=0}`` because of the zero-order hold):
 ```math
 \begin{aligned}
     \mathbf{x}(k+1) &=  \mathbf{A x}(k) + \mathbf{B_u u}(k) + \mathbf{B_d d}(k) \\
     \mathbf{y}(k)   &=  \mathbf{C x}(k) + \mathbf{D_d d}(k)
 \end{aligned}
 ```
+Use the syntax [`LinModel(A, Bu, C, Bd, Dd, Ts, nu, nx, ny, nd)`](@ref) to force a specific
+state-space representation.
 """
 function LinModel(
     sys::StateSpace,
@@ -180,6 +182,15 @@ function LinModel(sys::DelayLtiSystem, Ts::Real; kwargs...)
     return LinModel(sys_dis, Ts; kwargs...)
 end
 
+@doc raw"""
+    LinModel(A, Bu, C, Bd, Dd, Ts, nu, nx, ny, nd)
+
+Construct the model from the discrete state-space matrices `A, Bu, C, Bd, Dd` directly.
+
+This syntax do not modify the state-space representation provided in argument ([`minreal`](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/constructors/#ControlSystemsBase.minreal)
+is not called). Care must be taken to ensure that the model is controllable and observable.
+"""
+LinModel(A, Bu, C, Bd, Dd, Ts, nu, nx, ny, nd)
 
 """
     f(model::LinModel, x, u, d)

src/state_estim.jl

@@ -283,74 +283,76 @@ end
 
 Init `estim.x̂` states from current inputs `u`, measured outputs `ym` and disturbances `d`.
 
-The method set the error covariance to `estim.P = estim.P̂0` (if applicable) and tries to 
-find a good steady-state to initialize `estim.x̂` estimate. 
+The method removes the operating points with [`remove_op!`](@ref) and call 
+[`init_estimate!`](@ref):
 
 - If `estim.model` is a [`LinModel`](@ref), it finds the steady-state of the augmented model
   using `u` and `d` arguments, and uses the `ym` argument to enforce that ``\mathbf{ŷ^m} = 
   \mathbf{y^m}``. For control applications, this solution produces a bumpless manual to 
-  automatic transfer. See Extended Help for details.
+  automatic transfer. See [`init_estimate!`](@ref) for details.
 - Else, `estim.x̂` is left unchanged. Use [`setstate!`](@ref) to manually modify it.
 
+If applicable, it also sets the error covariance to `estim.P = estim.P̂0`.
+
 # Examples
 ```jldoctest
 julia> estim = SteadyKalmanFilter(LinModel(tf(3, [10, 1]), 0.5), nint_ym=[2]);
 
-julia> u = [1]; ym = [3 - 0.1]; x̂ = round.(initstate!(estim, u, ym), digits=3)
+julia> u = [1]; y = [3 - 0.1]; x̂ = round.(initstate!(estim, u, y), digits=3)
 3-element Vector{Float64}:
   5.0
   0.0
  -0.1
 
-julia> x̂ ≈ updatestate!(estim, u, ym)
+julia> x̂ ≈ updatestate!(estim, u, y)
 true
 
-julia> evaloutput(estim) ≈ ym
+julia> evaloutput(estim) ≈ y
 true
 ```
 
-# Extended Help
-Based on [`setop!`](@ref) notation, the resulting system to solve for [`LinModel`](@ref) is:
+"""
+function initstate!(estim::StateEstimator, u, ym, d=Float64[])
+    # --- init state estimate ----
+    u0, d0, ym0 = remove_op!(estim, u, d, ym)
+    init_estimate!(estim, estim.model, u0, ym0, d0)
+    # --- init covariance error estimate, if applicable ---
+    initstate_post!(estim)
+    return estim.x̂
+end
+
+@doc raw"""
+    init_estimate!(estim::StateEstimator, model::LinModel, u, ym, d)
+
+Init `estim.x̂` estimate with the steady-state solution if `model` is a [`LinModel`](@ref).
+
+Using `u`, `ym` and `d` arguments, the steady-state problem combined to the equality 
+constraint ``\mathbf{ŷ^m} = \mathbf{y^m}`` engenders the following system to solve :
 ```math
 \begin{bmatrix}
     \mathbf{I} - \mathbf{Â}             \\
     \mathbf{Ĉ}
 \end{bmatrix} \mathbf{x̂} =
 \begin{bmatrix}
-    \mathbf{B̂_u u_0} + \mathbf{B̂_d d_0} \\
-    \mathbf{y_0} - \mathbf{D̂_d d_0}
+    \mathbf{B̂_u u} + \mathbf{B̂_d d} \\
+    \mathbf{y} - \mathbf{D̂_d d}
 \end{bmatrix}
 ```
-with ``\mathbf{u_0 = u - u_{op}}`` and ``\mathbf{d_0 = d - d_{op}}``. The vector
-``\mathbf{y_0}`` comprises the measured ``\mathbf{y_0^m = y^m - y_{op}^m}`` and unmeasured
-``\mathbf{y_0^u = 0}`` outputs.
+in which ``\mathbf{y}`` comprises the measured ``\mathbf{y^m}`` and unmeasured
+``\mathbf{y^u = 0}`` outputs.
 """
-function initstate!(estim::StateEstimator, u, ym, d=Float64[])
-    model = estim.model
-    # --- init covariance error estimate (if applicable) ---
-    initstate_cov!(estim)
-    # --- init lastu0 for PredictiveControllers ---
-    estim.lastu0[:] = u - model.uop
-    # --- init state estimate ----
-    init_estimate!(estim, model, u, ym, d)
-    return estim.x̂
-end
-
-"By default, state estimators do not need initialization of covariance estimate."
-initstate_cov!(::StateEstimator) = nothing
-
-"Init estimate x̂ with the steady-state solution for if `model` is a [`LinModel`](@ref)."
 function init_estimate!(estim::StateEstimator, model::LinModel, u, ym, d)
     Â, B̂u, Ĉ, B̂d, D̂d = estim.Â, estim.B̂u, estim.Ĉ, estim.B̂d, estim.D̂d
-    y0 = zeros(model.ny)
-    y0[estim.i_ym] = ym - model.yop[estim.i_ym]
-    u0 = u - model.uop
-    d0 = d - model.dop
-    estim.x̂[:] = [(I - Â); Ĉ]\[B̂u*u0 + B̂d*d0; y0 - D̂d*d0]
+    y = zeros(model.ny)
+    y[estim.i_ym] = ym
+    estim.x̂[:] = [(I - Â); Ĉ]\[B̂u*u + B̂d*d; y - D̂d*d]
 end
-"Do nothing if `model` is not a [`LinModel`](@ref)."
+"Left `estim.x̂` estimate unchanged if `model` is not a [`LinModel`](@ref)."
 init_estimate!(::StateEstimator, ::SimModel, _ , _ , _ ) = nothing
 
+"By default, do nothing at the end of `initstate!` (used to init the covariance estimate)."
+initstate_post!(::StateEstimator) = nothing
+
 @doc raw"""
     evaloutput(estim::StateEstimator, d=Float64[]) -> ŷ
 

test/test_sim_model.jl

@@ -57,7 +57,6 @@ Gss2 = c2d(sys_ss[:,1:2], 0.5Ts, :zoh)
 
     linmodel6 = LinModel([delay(4) delay(4)]*sys,Ts,i_d=[3])
     @test linmodel6.nx == 3
-    println(eigvals(linmodel6.A))
     @test sum(eigvals(linmodel6.A) .≈ 0) == 1
 
     linmodel7 = LinModel(