PBH test for observability of linear systems

Author: franckgaga

src/model/linmodel.jl

@@ -182,18 +182,20 @@ end
 
 Evaluate the steady-state vector when `model` is a [`LinModel`](@ref).
 
-Omitting the operating points, the method evaluates the equilibrium ``\mathbf{x}(∞)`` from:
+Following [`setop!`](@ref) notation, the method evaluates the equilibrium ``\mathbf{x}(∞)``
+from:
 ```math
-    \mathbf{x}(∞) = \mathbf{(I - A)^{-1}(B_u u + B_d d)}
+    \mathbf{x}(∞) = \mathbf{(I - A)^{-1}(B_u u_0 + B_d d_0)}
 ```
-with the manipulated inputs held constant at ``\mathbf{u}`` and, the measured disturbances, 
-at ``\mathbf{d}``. The Moore-Penrose pseudo-inverse computes ``\mathbf{(I - A)^{-1}}``
-to support integrating `model` (integrator states will be 0).
+with the manipulated inputs held constant at ``\mathbf{u_0}`` and, the measured
+disturbances, at ``\mathbf{d_0}``. The Moore-Penrose pseudo-inverse computes 
+``\mathbf{(I - A)^{-1}}`` to support integrating `model` (integrator states will be 0).
 """
 function steadystate(model::LinModel, u, d=Float64[])
     return pinv(I - model.A)*(model.Bu*(u - model.uop) + model.Bd*(d - model.dop))
 end
 
+
 """
     f(model::LinModel, x, u, d)
 

src/state_estim.jl

@@ -154,6 +154,15 @@ function augment_model(model::LinModel, As, Cs)
     return Â, B̂u, Ĉ, B̂d, D̂d
 end
 
+"""
+    isobservable(A, C)
+
+Verify if the state-space `A` and `C` matrices of are observable.
+
+Based on the [Popov-Belevitch-Hautus (PBH) test](https://en.wikipedia.org/wiki/Hautus_lemma).
+"""
+isobservable(A, C) = all(map(λ->rank([λ * I - A; C]) == size(A, 1), eigvals(A)))
+
 @doc raw"""
     f̂(estim::StateEstimator, x̂, u, d)
 
@@ -293,11 +302,12 @@ function updatestate!(estim::StateEstimator, u, ym, d=Float64[])
     return estim.x̂
 end
 
-include("estimator/kalman.jl")
-include("estimator/luenberger.jl")
-include("estimator/internal_model.jl")
-
 "Get [`InternalModel`](@ref) output `ŷ` from current measured outputs `ym` and dist. `d`."
 evalŷ(estim::InternalModel, ym, d) = evaloutput(estim,ym, d)
 "Other [`StateEstimator`](@ref) ignores `ym` to evaluate `ŷ`."
-evalŷ(estim::StateEstimator, _, d) = evaloutput(estim, d)
+evalŷ(estim::StateEstimator, _, d) = evaloutput(estim, d)
+
+
+include("estimator/kalman.jl")
+include("estimator/luenberger.jl")
+include("estimator/internal_model.jl")