Merge pull request #117 from JuliaControl/mtk_obsv_error

Added: details in observability error message, modify CI for new release and MTK exemple debug

Author: franckgaga

.github/workflows/CI.yml

@@ -19,8 +19,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
           - '1'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest

docs/Project.toml

@@ -16,4 +16,4 @@ JuMP = "1"
 LinearAlgebra = "1.6"
 Logging = "1.6"
 Plots = "1"
-ModelingToolkit = "9"
+ModelingToolkit = "~9.41"

docs/src/manual/mtk.md

@@ -22,6 +22,9 @@ the last section.
     as a basic starting template to combine both packages. There is no guarantee that it
     will work for all corner cases.
 
+!!! compat
+    The example only work with `ModelingToolkit.jl` v9.41 releases.
+
 We first construct and instantiate the pendulum model:
 
 ```@example 1
@@ -76,14 +79,18 @@ function generate_f_h(model, inputs, outputs)
         end
         return nothing
     end
-    h_ = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs)
+    return_inplace = true
+    (_, h_ip) = ModelingToolkit.build_explicit_observed_function(
+        io_sys, outputs; inputs, return_inplace
+    )
+    println(h_ip)
     u_nothing = fill(nothing, nu)
     function h!(y, x, _ , p)
-        y .= try
+        try
             # MTK.jl supports a `u` argument in `h_` function but not this package. We set
             # `u` as a vector of nothing and `h_` function will presumably throw an
             # MethodError it this argument is used inside the function
-            h_(x, u_nothing, p, nothing)
+            h_ip(y, x, u_nothing, p, nothing)
         catch err
             if err isa MethodError
                 error("NonLinModel only support strictly proper systems (no manipulated "*

src/estimator/construct.jl

@@ -156,7 +156,9 @@ function augment_model(model::LinModel{NT}, As, Cs_u, Cs_y; verify_obsv=true) wh
     if verify_obsv && !ControlSystemsBase.observability(Â, Ĉ)[:isobservable]
         error("The augmented model is unobservable. You may try to use 0 integrator on "*
               "model integrating outputs with nint_ym parameter. Adding integrators at both "*
-              "inputs (nint_u) and outputs (nint_ym) can also violate observability.")
+              "inputs (nint_u) and outputs (nint_ym) can also violate observability. If the "*
+              "model is still unobservable without any integrators, you may need to call "*
+              "sminreal or minreal on your system.")
     end
     x̂op, f̂op = [model.xop; zeros(nxs)], [model.fop; zeros(nxs)]
     return Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op

src/estimator/kalman.jl

@@ -38,11 +38,15 @@ struct SteadyKalmanFilter{NT<:Real, SM<:LinModel} <: StateEstimator{NT}
         Â, B̂u, Ĉ, B̂d, D̂d, x̂op, f̂op = augment_model(model, As, Cs_u, Cs_y)
         Ĉm, D̂dm = Ĉ[i_ym, :], D̂d[i_ym, :]
         validate_kfcov(nym, nx̂, Q̂, R̂)
+        if ny == nym
+            R̂_y = R̂
+        else
+            R̂_y = zeros(NT, ny, ny)
+            R̂_y[i_ym, i_ym] = R̂
+            R̂_y = Hermitian(R̂_y, :L)
+        end
         K̂ = try
-            Q̂_kalman = Matrix(Q̂) # Matrix() required for Julia 1.6
-            R̂_kalman = zeros(NT, ny, ny)
-            R̂_kalman[i_ym, i_ym] = R̂
-            ControlSystemsBase.kalman(Discrete, Â, Ĉ, Q̂_kalman, R̂_kalman; direct)[:, i_ym] 
+            ControlSystemsBase.kalman(Discrete, Â, Ĉ, Q̂, R̂_y; direct)[:, i_ym]
         catch my_error
             if isa(my_error, ErrorException)
                 error("Cannot compute the optimal Kalman gain K̂ for the "* 

src/model/solver.jl

@@ -47,13 +47,11 @@ function get_solver_functions(NT::DataType, solver::RungeKutta, fc!, hc!, Ts, _
     k4_cache::DiffCache{Vector{NT}, Vector{NT}}   = DiffCache(zeros(NT, nx), Nc)
     f! = function inner_solver_f!(xnext, x, u, d, p)
         CT = promote_type(eltype(x), eltype(u), eltype(d))
-        # dummy variable for get_tmp, necessary for PreallocationTools + Julia 1.6 :
-        var::CT = 0
-        xcur = get_tmp(xcur_cache, var)
-        k1   = get_tmp(k1_cache, var)
-        k2   = get_tmp(k2_cache, var)
-        k3   = get_tmp(k3_cache, var)
-        k4   = get_tmp(k4_cache, var)
+        xcur = get_tmp(xcur_cache, CT)
+        k1   = get_tmp(k1_cache, CT)
+        k2   = get_tmp(k2_cache, CT)
+        k3   = get_tmp(k3_cache, CT)
+        k4   = get_tmp(k4_cache, CT)
         xterm = xnext
         @. xcur = x
         for i=1:solver.supersample