# Solving model detection problems

Exact synthesis of model detection filters.`emdsyn`

Approximate synthesis of model detection filters.`amdsyn`

`FaultDetectionTools.emdsyn`

— Function```
emdsyn(sysm::MDMModel; rdim, MDSelect, nullspace = true, simple = false, minimal = true,
emtest = false, normalize = false, fast = true,
sdeg, smarg, poles, HDesign, MDtol, MDGainTol, MDfreq,
tcond, offset, atol, atol1, atol2, atol3, rtol)
-> (Q::MDFilter, R::MDFilterIF, info)
```

Solve the *exact model detection problem* (EMDP) for a given multiple synthesis model `sysm::MDMModel`

. The computed stable and proper filter objects `Q`

and `R`

contain the model detection filter, representing the solution of the EMDP, and its internal form, respectively.

The returned named tuple `info`

, with the components `info.tcond`

, `info.degs`

, `info.MDperf`

, and `info.HDesign`

, contains additional synthesis related information (see below).

For a multiple synthesis model `sysm`

containing `N`

stable component models, the `i`

-th component system `sysm.sys[i]`

has a descriptor realization of the form `sysm.sys[i] = (Ai-λEi,Bi,Ci,Di)`

and the corresponding input-output form is

` yi = Gui(λ)*u + Gdi(λ)*di + Gwi(λ)*wi + Gvi(λ)*vi`

with the Laplace- or Z-transformed plant outputs `yi`

, control inputs `u`

, disturbance inputs `di`

, noise inputs `wi`

, and auxiliary inputs `vi`

, and with `Gui(λ)`

, `Gdi(λ)`

, `Gwi(λ)`

and `Gvi(λ)`

, the corresponding transfer function matrices.

The model detection filter object `Q`

, contains in `Q.sys`

the resulting bank of `N`

filters. The `i`

-th filter `Q.sys[i]`

is in a standard state-space form and generates `r_i`

, the `i`

-th component (scalar or vector) of the overall residual vector `r := [r_1; r_2; ...; r_N]`

. The corresponding input-output (implementation) form of the `i`

-th filter is

` r_i = Qyi(λ)*y + Qui(λ)*u ,`

where `Qyi(λ)`

and `Qui(λ)`

are the transfer function matrices from the output and control inputs to the `i`

-th residual component. The dimensions of output and control inputs are contained in the integers `Q.ny`

and `Q.mu`

, respectively.

The model detection filter internal form object `R`

, contains `R.sys`

, the resulting array of `N × N`

internal form of the filters. The `(i,j)`

-th component filter `R.sys[i,j]`

is in a standard state-space form, generates the residual signal `r_ij`

, and corresponds to the input-output form

` r_ij = Ruij(λ)*u + Rdij(λ)*dj + Rwij(λ)*wj + Rvij(λ)*vj ,`

where

```
| Ruij(λ) Rdij(λ) Rwij(λ) Raij(λ) | := |Qyi(λ) Qui(λ)]*| Guj(λ) Gdj(λ) Gwj(λ) Gaj(λ) |.
| I 0 0 0 |
```

The solution of the EMDP ensures that for the `i`

-th filter, `Ruii(λ) = 0`

, `Rdii(λ) = 0`

, and `[Ruij(λ) Rdij(λ)]`

$\neq$ `0`

for `j`

$\neq$ `i`

.

Various user options can be specified via keyword arguments as follows:

`MDSelect = ifilt`

specifies in the vector `ifilt`

the indices of the desired filters to be designed (default: `ifilt = 1:N`

)

If `minimal = true`

(default), least order filter synthesis is performed to determine each of the component filters `Q.sys[i]`

for `i = 1, ...,N`

and `R.sys[i,j]`

for `i, j = 1, ...,N`

while with `minimal = false`

full order synthesis is performed.

If `HDesign = H`

, then `H`

is an `N`

-dimensional array of full row rank or empty design matrices `H = [H_1, ..., H_N]`

, where `H_i`

is the design matrix employed for the synthesis of the `i`

-th component filter (default: `HDesign = missing`

)

`rdim = q`

specifies the vector `q`

, whose `i`

-th component `q[i]`

specifies the number of residual outputs for the `i`

-th component filter `Q.sys[i]`

. If `q`

is a scalar, then a vector `rdim`

with all components equal to `q`

is assumed. The default value of `q[i]`

is chosen as follows: if `HDesign = missing`

or `H_i`

is empty then `q[i] = 1`

, if `minimal = true`

, or `q[i]`

is the number of the nullspace basis vectors used for the synthesis of `Q.sys [i]`

, if `minimal = false`

; if `H_i`

specifies a full row rank design matrix, then `q[i]`

is the row dimension of `H_i`

.

If `emdtest = false`

(default), only the input channels are used for model detectability tests. If `emdtest = true`

, extended model detectability tests are performed using both control and disturbance input channels.

`MDfreq = ω`

specifies a vector of real frequency values or a scalar real frequency value for strong model detectability checks (default: `MDfreq = missing`

).

If `nullspace = false`

(default), a full-order observer based nullspace basis is used for the synthesis of the `i`

-th filter whenever the `i`

-th component system has no disturbance inputs. Otherwise, minimal proper nullspace bases are used. If `nullspace = true`

, minimal proper nullspace bases are used for the synthesis of the model detection filters.

If `simple = true`

, simple proper nullspace bases are emplyed for synthesis. The orders of the basis vectors employed for the synthesis of the `i`

-th filter are provided in `info.deg[i]`

. If `simple = false`

(default), then no simple bases are computed.

`offset = β`

specifies the boundary offset `β`

to assess the stability of poles. Accordingly, for the stability of a continuous-time system all real parts of poles must be at most `-β`

, while for the stability of a discrete-time system all moduli of poles must be at most `1-β`

. The default value used for `β`

is `sqrt(ϵ)`

, where `ϵ`

is the working machine precision.

`smarg = α`

specifies the stability margin which defines the stability domain `Cs`

of poles, as follows: for a continuous-time system, `Cs`

is the set of complex numbers with real parts at most `α`

, while for a discrete-time system, `Cs`

is the set of complex numbers with moduli at most `α < 1`

(i.e., the interior of a disc of radius `α`

centered in the origin). If `smarg`

is missing, then the employed default values are `α = -β`

for a continuous-time system and `α = 1-β`

for a discrete-time system, where `β`

is the boundary offset specified by the keyword argument `offset = β`

.

`sdeg = γ`

is the prescribed stability degree for the poles of the filters `Q`

and `R`

(default: `γ = -0.05`

for the real parts of poles for a continuous-time system and `γ = 0.95`

for the magnitudes of poles for a discrete-time system).

`poles = v`

specifies a complex vector `v`

containing a complex conjugate set of desired poles within the stability domain `Cs`

to be assigned for the filters `Q`

and `R`

(default: `poles = missing`

).

`tcond = tcmax`

specifies the maximum alowed condition number `tcmax`

of the employed non-orthogonal transformations (default: `tcmax = 1.e4`

).

`MDtol = tol1`

specifies the threshold `tol1`

for model detectability checks (default: `tol1 = 0.0001`

).

`MDGainTol = tol2`

specifies the threshold `tol2`

for strong model detectability checks (default: `tol2 = 0.01`

).

If `normalize = false`

(default), the `i`

-th component filter `Q.sys[i]`

is scaled such that the minimum gain of `R.sys[i,j]`

for `j = 1, ..., N`

, `j`

$\neq$ `i`

, is equal to one. If `normalize = true`

, the standard normalization of component filters is performed to ensure equal gains for `R.sys[1,j]`

and `R.sys[j,1]`

. This option is only possible if `ifilt[1] = 1`

(see `MDSelect`

).

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if `fast = true`

(default) or the more reliable SVD-decompositions if `fast = false`

.

The keyword arguments `atol1`

, `atol2`

and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of the state-space matrices `Ai`

, `Bi`

, `Ci`

and `Di`

, the absolute tolerance for the nonzero elements of `Ei`

, and the relative tolerance for the nonzero elements of all above matrices. The default relative tolerance is `ni*ϵ`

, where `ϵ`

is the working machine epsilon and `ni`

is the orders of the system `sysm.sys[i]`

. The keyword argument `atol3`

is an absolute tolerance for observability tests (default: internally determined value). The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

and `atol3 = atol`

.

The resulting named tuple `info`

contains `(tcond, degs, MDperf, HDesign)`

, where:

`info.tcond`

is the maximum of the condition numbers of the employed non-orthogonal transformation matrices; a warning is issued if `info.tcond >= tcmax`

;

`info.degs`

is an `N`

-dimensional vector, whose `i`

-th component is an integer vector containing the degrees of the basis vectors of the employed simple nullspace basis for the synthesis of the `i`

-th filter component, if `simple = true`

, and the degrees of the basis vectors of an equivalent polynomial nullspace basis, if `simple = false`

;

`info.MDperf`

is an `N × N`

array containing the resulting distance-mapping performance, representing the peak gains of the associated internal representations of the model detection filters. If the `(i,j)`

-th component filter `R.sys[i,j]`

has the input-output form

` rij = Ruij(λ)*u + Rdij(λ)*dj + Rwij(λ)*wj + Rvij(λ)*vj ,`

then, the `(i,j)`

-th performance gain is evaluated as `info.MDperf[i,j] = ||Ruij(λ)||∞`

if `emdtest = false`

(default) or `info.MDperf[i,j] = ||[Ruij(λ) Rdij(λ)]||∞`

if `emdtest = true`

. If `MDfreq = ω`

, then each gain `info.MDperf[i,j]`

represents the maximum of 2-norm pointwise gains evaluated for all frequencies in `ω`

.

`info.HDesign`

is an `N`

-dimensional vector, whose `i`

-th component is the design matrix `H_i`

employed for the synthesis of the `i`

-th model detection filter.

*Method:* The Procedure EMD from [1] is implemented to solve the exact model detection problem. For more details on the least order synthesis of model detection filters see [2].

*References:*

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017; sec. 5.2.

[2] A. Varga, Least order fault and model detection using multi-models. IEEE CDC'09, Shanghai, China, 2009.

`FaultDetectionTools.amdsyn`

— Function```
amdsyn(sysm::MDMModel; rdim, MDSelect, nullspace = true, simple = false, minimal = true,
emtest = false, normalize = false, fast = true,
sdeg, smarg, poles, HDesign, MDtol, MDGainTol, MDfreq,
epsreg = 0.1, sdegzer, nonstd = 1, rtolinf= 0.0001,
tcond, offset, atol, atol1, atol2, atol3, rtol)
-> (Q::MDFilter, R::MDFilterIF, info)
```

Solve the *approximate model detection problem* (AMDP) for a given multiple synthesis model `sysm::MDMModel`

. The computed stable and proper filter objects `Q`

and `R`

contain the model detection filter, representing the solution of the AMDP, and its internal form, respectively.

The returned named tuple `info`

, with the components `info.tcond`

, `info.degs`

, `info.MDperf`

, `info.MDgap`

, `info.HDesign`

and `info.nonstandard`

, contains additional synthesis related information (see below).

For a multiple synthesis model `sysm`

containing `N`

stable component models, the `i`

-th component system `sysm.sys[i]`

has a descriptor realization of the form `sysm.sys[i] = (Ai-λEi,Bi,Ci,Di)`

and the corresponding input-output form is

` yi = Gui(λ)*u + Gdi(λ)*di + Gwi(λ)*wi + Gvi(λ)*vi`

with the Laplace- or Z-transformed plant outputs `yi`

, control inputs `u`

, disturbance inputs `di`

, noise inputs `wi`

, and auxiliary inputs `vi`

, and with `Gui(λ)`

, `Gdi(λ)`

, `Gwi(λ)`

and `Gvi(λ)`

, the corresponding transfer function matrices.

The model detection filter object `Q`

, contains in `Q.sys`

the resulting bank of `N`

filters. The `i`

-th filter `Q.sys[i]`

is in a standard state-space form and generates `r_i`

, the `i`

-th component (scalar or vector) of the overall residual vector `r := [r_1; r_2; ...; r_N]`

. The corresponding input-output (implementation) form of the `i`

-th filter is

` r_i = Qyi(λ)*y + Qui(λ)*u ,`

where `Qyi(λ)`

and `Qui(λ)`

are the transfer function matrices from the measured outputs and control inputs to the `i`

-th residual component. The dimensions of measured outputs and control inputs are contained in the integers `Q.ny`

and `Q.mu`

, respectively.

The model detection filter internal form object `R`

, contains `R.sys`

, the resulting array of `N × N`

internal form filters. The `(i,j)`

-th component filter `R.sys[i,j]`

is in a standard state-space form, generates the residual signal `r_ij`

, and corresponds to the input-output form

` r_ij = Ruij(λ)*u + Rdij(λ)*dj + Rwij(λ)*wj + Rvij(λ)*vj ,`

where

```
| Ruij(λ) Rdij(λ) Rwij(λ) Raij(λ) | := |Qyi(λ) Qui(λ)]*| Guj(λ) Gdj(λ) Gwj(λ) Gaj(λ) |.
| I 0 0 0 |
```

The solution of the AMDP ensures that for the `i`

-th filter, `Ruii(λ) = 0`

, `Rdii(λ) = 0`

, `[Ruij(λ) Rdij(λ)]`

$\neq$ `0`

for `j`

$\neq$ `i`

, and, additionally, the maximization of the achieved gaps (see description of `info.MDgap`

).

Various user options can be specified via keyword arguments as follows:

`MDSelect = ifilt`

specifies in the vector `ifilt`

the indices of the desired filters to be designed (default: `ifilt = 1:N`

)

If `minimal = true`

(default), least order filter synthesis is performed to determine each of the component filters `Q.sys[i]`

for `i = 1, ...,N`

and `R.sys[i,j]`

for `i, j = 1, ...,N`

while with `minimal = false`

full order synthesis is performed.

If `HDesign = H`

, then `H`

is an `N`

-dimensional array of full row rank or empty design matrices `H = [H_1, ..., H_N]`

, where `H_i`

is the design matrix employed for the synthesis of the `i`

-th component filter (default: `HDesign = missing`

)

`rdim = q`

specifies the vector `q`

, whose `i`

-th component `q[i]`

specifies the number of residual outputs for the `i`

-th component filter `Q.sys[i]`

. If `q`

is a scalar, then a vector `rdim`

with all components equal to `q`

is assumed. The default value of `q[i]`

is chosen as follows: if `HDesign = missing`

or `H_i`

is empty then `q[i] = 1`

, if `minimal = true`

, or `q[i]`

is the number of the nullspace basis vectors used for the synthesis of `Q.sys [i]`

, if `minimal = false`

; if `H_i`

specifies a full row rank design matrix, then `q[i]`

is the row dimension of `H_i`

.

If `emdtest = false`

(default), only the input channels are used for model detectability tests. If `emdtest = true`

, extended model detectability tests are performed using both control and disturbance input channels.

`MDfreq = ω`

specifies a vector of real frequency values or a scalar real frequency value for strong model detectability checks (default: `MDfreq = missing`

).

If `nullspace = false`

(default), a full-order observer based nullspace basis is used for the synthesis of the `i`

-th filter whenever the `i`

-th component system has no disturbance inputs. Otherwise, minimal proper nullspace bases are used. If `nullspace = true`

, minimal proper nullspace bases are used for the synthesis of the model detection filters.

If `simple = true`

, simple proper nullspace bases are emplyed for synthesis. The orders of the basis vectors employed for the synthesis of the `i`

-th filter are provided in `info.deg[i]`

. If `simple = false`

(default), then no simple bases are computed.

`offset = β`

specifies the boundary offset `β`

to assess the stability of poles. Accordingly, for the stability of a continuous-time system all real parts of poles must be at most `-β`

, while for the stability of a discrete-time system all moduli of poles must be at most `1-β`

. The default value used for `β`

is `sqrt(ϵ)`

, where `ϵ`

is the working machine precision.

`smarg = α`

specifies the stability margin which defines the stability domain `Cs`

of poles, as follows: for a continuous-time system, `Cs`

is the set of complex numbers with real parts at most `α`

, while for a discrete-time system, `Cs`

is the set of complex numbers with moduli at most `α < 1`

(i.e., the interior of a disc of radius `α`

centered in the origin). If `smarg`

is missing, then the employed default values are `α = -β`

for a continuous-time system and `α = 1-β`

for a discrete-time system, where `β`

is the boundary offset specified by the keyword argument `offset = β`

.

`sdeg = γ`

is the prescribed stability degree for the poles of the filters `Q`

and `R`

(default: `γ = -0.05`

for the real parts of poles for a continuous-time system and `γ = 0.95`

for the magnitudes of poles for a discrete-time system).

`poles = v`

specifies a complex vector `v`

containing a complex conjugate set of desired poles within the stability domain `Cs`

to be assigned for the filters `Q`

and `R`

(default: `poles = missing`

).

`tcond = tcmax`

specifies the maximum alowed condition number `tcmax`

of the employed non-orthogonal transformations (default: `tcmax = 1.e4`

).

`MDtol = tol1`

specifies the threshold `tol1`

for model detectability checks (default: `tol1 = 0.0001`

).

`MDGainTol = tol2`

specifies the threshold `tol2`

for strong model detectability checks (default: `tol2 = 0.01`

).

If `normalize = false`

(default), the `i`

-th component filter `Q.sys[i]`

is scaled such that the minimum gain of `R.sys[i,j]`

for `j = 1, ..., N`

, `j`

$\neq$ `i`

, is equal to one. If `normalize = true`

, the standard normalization of component filters is performed to ensure equal gains for `R.sys[1,j]`

and `R.sys[j,1]`

. This option is only possible if `ifilt[1] = 1`

(see `MDSelect`

).

`epsreg = ϵ`

specifies the value of the regularization parameter `ϵ`

(default: `ϵ = 0.1`

)

`sdegzer = δ`

specifies the prescribed stability degree `δ`

for zeros shifting (default: `δ = −0.05`

for a continuous-time system `sysf.sys`

and `δ = 0.95`

for a discrete-time system `sysf.sys`

).

`nonstd = job`

specifies the option to handle nonstandard optimization problems as follows:

```
job = 1 – use the quasi-co-outer–co-inner factorization (default);
job = 2 – use the modified co-outer–co-inner factorization with the
regularization parameter `ϵ`;
job = 3 – use the Wiener-Hopf type co-outer–co-inner factorization;
job = 4 – use the Wiener-Hopf type co-outer-co-inner factorization with
zero shifting of the non-minimum phase factor using the
stabilization parameter `δ`;
job = 5 – use the Wiener-Hopf type co-outer-co-inner factorization with
the regularization of the non-minimum phase factor using the
regularization parameter `ϵ`.
```

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if `fast = true`

(default) or the more reliable SVD-decompositions if `fast = false`

.

The keyword arguments `atol1`

, `atol2`

and `rtol`

, specify, respectively, the absolute tolerance for the nonzero elements of the state-space matrices `Ai`

, `Bi`

, `Ci`

and `Di`

, the absolute tolerance for the nonzero elements of `Ei`

, and the relative tolerance for the nonzero elements of all above matrices. The default relative tolerance is `ni*ϵ`

, where `ϵ`

is the working machine epsilon and `ni`

is the orders of the system `sysm.sys[i]`

. The keyword argument `atol3`

is an absolute tolerance for observability tests (default: internally determined value). The keyword argument `atol`

can be used to simultaneously set `atol1 = atol`

, `atol2 = atol`

and `atol3 = atol`

.

The resulting named tuple `info`

contains `(tcond, degs, MDperf, MDgap, HDesign, nonstandard)`

, where:

`info.tcond`

is the maximum of the condition numbers of the employed non-orthogonal transformation matrices; a warning is issued if `info.tcond >= tcmax`

;

`info.degs`

is an `N`

-dimensional vector, whose `i`

-th component is an integer vector containing the degrees of the basis vectors of the employed simple nullspace basis for the synthesis of the `i`

-th filter component, if `simple = true`

, and the degrees of the basis vectors of an equivalent polynomial nullspace basis, if `simple = false`

;

`info.MDperf`

is an `N × N`

array containing the resulting distance-mapping performance, representing the peak gains of the associated internal representations of the model detection filters. If the `(i,j)`

-th component filter `R.sys[i,j]`

has the input-output form

` rij = Ruij(λ)*u + Rdij(λ)*dj + Rwij(λ)*wj + Rvij(λ)*vj ,`

then, the `(i,j)`

-th performance gain is evaluated as `info.MDperf[i,j] = ||Ruij(λ)||∞`

if `emdtest = false`

(default) or `info.MDperf[i,j] = ||[Ruij(λ) Rdij(λ)]||∞`

if `emdtest = true`

. If `MDfreq = ω`

, then each gain `info.MDperf[i,j]`

represents the maximum of 2-norm pointwise gains evaluated for all frequencies in `ω`

.

`info.MDgap`

is an `N`

-dimensional vector of achieved noise gaps, where the `i`

-th gap is computed as `info.MDgap[i] = β[i]/γ[i]`

, where `β[i] = min(||Rij(λ)||∞`

for `i`

$\neq$ `j`

) and `γ[i] = ||Rwii(λ)||∞`

. `Rij(λ)`

is defined as `Rij(λ) = Ruij(λ)`

if `emtest = false`

(default) or `Rij(λ) = [Ruij(λ) Rdij(λ)]`

if `emtest = true`

. If `MDfreq = ω`

, where `ω`

is a given vector of real frequency values, then each gain `β[i]`

represents the minimum of the maximum of 2-norm pointwise gains evaluated for all frequencies in `ω`

.

`info.HDesign`

is an `N`

-dimensional vector, whose `i`

-th component is the design matrix `H_i`

employed for the synthesis of the `i`

-th model detection filter.

`info.nonstandard`

is an `N`

-dimensional vector, whose `i`

-th component is `true`

if the synthesis of the `i`

-th filter involved a nonstandard optimization problem and is `false`

otherwise.

*Method:* An extension of the Procedure AMD from [1] is implemented to solve the approximate model detection problem. This procedure relies on the nullspace-based synthesis method of model detection filters proposed in [2].

*References:*

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017; sec. 5.2.

[2] A. Varga, Least order fault and model detection using multi-models. IEEE CDC'09, Shanghai, China, 2009.