Solving fault detection and isolation problems

  • efdsyn Exact synthesis of fault detection filters.
  • efdisyn Exact synthesis of fault detection and isolation filters.
  • afdsyn Approximate synthesis of fault detection filters.
  • afdisyn Approximate synthesis of fault detection and isolation filters.
  • emmsyn Exact model-matching based synthesis of fault detection filters.
  • ammsyn Approximate model-matching based synthesis of fault detection filters.
FaultDetectionTools.efdsyn — Function
efdsyn(sysf::FDIModel; rdim, nullspace = true, simple = false, minimal = true, 
                       sdeg, smarg, poles, HDesign, FDtol, FDGainTol, FDfreq, 
                       tcond, offset, atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the exact fault detection problem (EFDP) for a given synthesis model sysf with additive faults. The computed stable and proper filter objects Q and R contain the fault detection filter, representing the solution of the EFDP, and its internal form, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.S, and info.HDesign, contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The fault detection filter object Q, contains in Q.sys the resulting filter in a standard state-space form, which generates the residual signal r. The corresponding input-output (implementation) form is

        r = Qy(λ)*y + Qu(λ)*u

where Qy(λ) and Qu(λ) are the transfer function matrices from the output and control inputs to the residual. The indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

The fault detection filter in internal form object R, contains R.sys, the resulting internal form of the filter in a standard state-space form, which generates the residual signal r, and corresponds to the input-output form

   r = Ru(λ)*u + Rd(λ)*d + Rf(λ)*f + Rw(λ)*w + Ra(λ)*aux ,

where 

   | Ru(λ) Rd(λ) Rf(λ) Rw(λ) Ra(λ) | = |Qy(λ) Qu(λ)|*| Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |. 
                                                     |  I     0     0     0     0    |

The solution of the EFDP ensures that Ru(λ) = 0, Rd(λ) = 0, and Rf(λ) has all its columns nonzero. The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls (void), R.disturbances (void), R.faults, R.noise and R.aux, respectively.

The resulting filters Q.sys and R.sys have observable state-space realizations (AQ,BQ,CQ,DQ) and (AQ,BR,CQ,DR), respectively, and thus share the observable pairs (AQ,CQ).

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

If minimal = true (default), a least order filter synthesis is performed, while with minimal = false a full order synthesis is performed.

If HDesign = H, a full row rank design matrix H is used to build rdim = q linear combinations of the left nullspace basis vectors (default: HDesign = missing)

rdim = q specifies the desired number q of residual outputs for Q and R. The default value of q is chosen as follows: if HDesign = missing, then q = 1, if minimal = true, or q is the number of the nullspace basis vectors used for the initial synthesis, if minimal = false; if HDesign = H specifies a full row rank design matrix H, then q is the row dimension of H.

FDfreq = freq specifies a vector of real frequency values or a scalar real frequency value for strong detectability checks (default: FDfreq = missing).

If nullspace = true (default), a minimal proper nullspace basis is used for the synthesis of the fault detection filter. If nullspace = false, a full-order observer based nullspace basis is used. This option can be only used for a proper system without disturbance inputs.

If simple = true, a simple proper nullspace basis is emplyed for synthesis. The orders of the basis vectors are provided in info.deg. If simple = false (default), then no simple basis is 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).

FDtol = tol1 specifies the threshold tol1 for fault detectability checks (default: tol1 = 0.0001).

FDGainTol = tol2 specifies the threshold tol2 for strong fault detectability checks (default: tol2 = 0.01).

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, S, 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 integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G(λ) := [ Gu(λ) Gd(λ); I 0] (also the left Kronecker indices of G(λ)), if the state-space realization of [Gu(λ) Gd(λ)] is minimal;

info.S is the binary structure matrix corresponding to the computed left nullspace basis;

info.HDesign is the design matrix H employed for the synthesis of the fault detection filter.

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

References:

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

[2] A. Varga, On computing least order fault detectors using rational nullspace bases. IFAC SAFEPROCESS'03 Symposium, Washington DC, USA, 2003.

efdsyn(sysf::FDIModel, S; rdim, nullspace = true, simple = false, minimal = true, 
                       sdeg, smarg, poles, HDesign, FDtol, FDGainTol, FDfreq, 
                       tcond, offset, atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the exact fault detection and isolation problem (EFDIP) for a given synthesis model sysf with additive faults and a given binary structure vector S. The computed stable and proper filter objects Q and R contain the fault detection filter, representing the solution of the EFDIP, and its internal form, respectively, and are determined such that R.sys[:,faults] has its j-th column nonzero if S[j] = 1 and the j-th column is zero if S[j] = 0. For the description of the keyword parameters see the function efdsyn.

FaultDetectionTools.efdisyn — Function
efdisyn(sysf::FDIModel, SFDI; rdim, nullspace = true, simple = false, minimal = true, separate = false,
                       sdeg, smarg, poles, HDesign, FDtol, FDGainTol, FDfreq, 
                       tcond, offset, atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the exact fault detection and isolation problem (EFDIP) for a given synthesis model sysf with additive faults and a given binary structure matrix SFDI with nb rows (specifications). The computed stable and proper filter objects Q and R contain the fault detection and isolation filter, representing the solution of the EFDIP, and its internal form, respectively.

The returned named tuple info, with the components info.tcond, info.degs and info.HDesign, contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The fault detection and isolation filter object Q, contains in Q.sys the resulting bank of nb 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_nb]. 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 indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

The fault detection and isolation filter in internal form object R, contains R.sys, the resulting bank of nb internal form of the filters. The i-th filter R.sys[i] is in a standard state-space form, which generates the residual signal r_i, and corresponds to the input-output form

   r_i = Rui(λ)*u + Rdi(λ)*d + Rfi(λ)*f + Rwi(λ)*w + Rai(λ)*aux ,

where 

   | Rui(λ) Rdi(λ) Rfi(λ) Rwi(λ) Rai(λ) | := |Qyi(λ) Qui(λ)]*| Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |. 
                                                             |   I     0     0     0     0   |

The solution of the EFDIP ensures that for the i-th filter, Rui(λ) = 0, Rdi(λ) = 0, and Rfi(λ) has its j-th column nonzero if the (i,j)-th element of SFDI is nonzero. The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls (void), R.disturbances (void), R.faults, R.noise and R.aux, respectively.

The resulting component filters Q.sys[i] and R.sys[i] have observable state-space realizations (AQi,BQi,CQi,DQi) and (AQi,BRi,CQi,DRi), respectively, and thus share the observable pairs (AQi,CQi).

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

If minimal = true (default), least order filter synthesis is performed to determine each of the component filters Q.sys[i] and R.sys[i] for i = 1, ...,nb, while with minimal = false full order synthesis is performed.

If HDesign = H, then H is an nb-dimensional array of full row rank or empty design matrices H = [H_1, ..., H_nb], 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 filters Q.sys[i] and R.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] and R.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.

FDfreq = freq specifies a vector of real frequency values or a scalar real frequency value for strong detectability checks (default: FDfreq = missing).

If nullspace = true (default), a minimal proper nullspace basis is used at the initial reduction step, if separate = false, or at all synthesis steps, if separate = true. If nullspace = false, a full-order observer based nullspace basis is used at the initial reduction step, if separate = false, or at all synthesis steps, if separate = true. This option can only be used for a proper system without disturbance inputs.

If simple = true, simple proper nullspace bases are emplyed for synthesis. The orders of the basis vectors employed for the synthesis of i-th filter are provided in info.deg[i]. If simple = false (default), then no simple bases are computed.

If separate = false (default), a two-step synthesis procedure is employed, where a minimal proper nullspace basis is used at the initial reduction step. If separate = true, the filter components are separately determined by solving appropriately formulated fault detection problems.

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).

FDtol = tol1 specifies the threshold tol1 for fault detectability checks (default: tol1 = 0.0001).

FDGainTol = tol2 specifies the threshold tol2 for strong fault detectability checks (default: tol2 = 0.01).

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, HDesign), where:

info.tcond is an nb-dimensional vector, whose i-th component is the maximum of the condition numbers of the employed non-orthogonal transformation matrices employed for the synthesis of the i-th filter component; a warning is issued if any info.tcond[i] >= tcmax;

info.degs is an nb-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, ifsimple = false`;

info.HDesign is an nb-dimensional vector, whose i-th component is the design matrix H_i employed for the synthesis of the i-th fault detection filter.

Method: The Procedure EFDI from [1] is implemented to solve the exact fault detection and isolation problem. This procedure relies on the nullspace-based synthesis method proposed in [2]. For more details on the least order synthesis of fault detection filters see [3].

References:

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

[2] A. Varga, On designing least order residual generators for fault detection and isolation. 16th International Conference on Control Systems and Computer Science, Bucharest, Romania, 2007.

[3] A. Varga, On computing least order fault detectors using rational nullspace bases. IFAC SAFEPROCESS'03 Symposium, Washington DC, USA, 2003.

FaultDetectionTools.afdsyn — Function
afdsyn(sysf::FDIModel; rdim, nullspace = true, simple = false, minimal = true, exact = false, 
                       gamma = 1, epsreg = 0.1, sdegzer, nonstd = 1, freq, sdeg, smarg, poles, 
                       HDesign, HDesign2, scale2, FDtol, FDGainTol, FDfreq, 
                       tcond, offset, atol, atol1, atol2, atol3, rtol, fast = true) 
                          -> (Q::FDFilter, R::FDFilterIF, info)

Solve the approximate fault detection problem (AFDP) for a given synthesis model sysf with additive faults. The computed stable and proper filter objects Q and R contain the fault detection filter, representing the solution of the AFDP, and its internal form, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.degs2, info.S, info.S2, info.HDesign, info.HDesign2, info.freq and info.gap contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The fault detection filter object Q, contains in Q.sys the resulting filter in a standard state-space form, which generates the residual signal r. The corresponding input-output (implementation) form is

        r = Qy(λ)*y + Qu(λ)*u  ,

where Qy(λ) and Qu(λ) are the transfer function matrices from the output and control inputs to the residual. The indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

The fault detection filter in internal form object R, contains R.sys, the resulting internal form of the filter in a standard state-space form, which generates the residual signal r, and corresponds to the input-output form

   r = Ru(λ)*u + Rd(λ)*d + Rf(λ)*f + Rw(λ)*w + Ra(λ)*aux ,

where 

   | Ru(λ) Rd(λ) Rf(λ) Rw(λ) Ra(λ) | = |Qy(λ) Qu(λ)|*| Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |. 
                                                     |  I     0     0     0     0    |

The solution of the AFDP ensures that Ru(λ) = 0, Rd(λ) = 0, Rf(λ) has all its columns nonzero and the H∞-norm of Rw(λ) satisfies ||Rw(λ)||∞ < γ, where the bound γ is specified via the keyword argument gamma. The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls (void), R.disturbances (void), R.faults, R.noise and R.aux, respectively.

The transfer function matrices Q(λ) = [ Qy(λ) Qu(λ) ] and R(λ) = [ Ru(λ) Rd(λ) Rf(λ) Rw(λ) Ra(λ) ] of the resulting filters Q.sys and R.sys, respectively, have, in general, the partitioned forms

 Q(λ) = [ Q1(λ) ] ,   R(λ) = [ R1(λ) ] ,                      (1)
        [ Q2(λ) ]            [ R2(λ) ]

where the filters Q1(λ) and R1(λ) with q1 outputs are the solution of an AFDP, while the filters Q2(λ) and R2(λ) with q2 outputs are the solution of an exact fault detection problem formulated for a reduced system obtained by decoupling the control and disturbance inputs from the residuals (see [4]). The overall resulting filters Q.sys and R.sys have observable state-space realizations (AQ,BQ,CQ,DQ) and (AQ,BR,CQ,DR), respectively, and thus share the observable pairs (AQ,CQ).

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

If minimal = true (default), a least order filter synthesis is performed, while with minimal = false a full order synthesis is performed.

If exact = true, an exact synthesis (without optimization) is performed, while with exact = false (default), an approximate synthesis is performed.

If HDesign = H1, a design matrix H1 of full row rank q1 is used to build q1 linear combinations of the left nullspace basis vectors of G1(λ) := [ Gu(λ) Gd(λ); I 0]. H1 is used in the synthesis of the components Q1(λ) and R1(λ) in (1) (default: HDesign = missing).

If HDesign2 = H2, a design matrix H2 of full row rank q2 is used to build q2 linear combinations of the left nullspace basis vectors of G2(λ) := [ Gu(λ) Gd(λ) Gw(λ); I 0 0]. H2 is used in the synthesis of the components Q2(λ) and R2(λ) in (1) (default: HDesign2 = missing)

rdim = q specifies the desired number q of residual outputs for Q and R. If rdim = missing, the default value of q is chosen as q = q1 + q2, where the default values of q1 and q2 are chosen taking into account the rank rw of Rw(λ) in the reduced system (see [4]), as follows: if HDesign = missing, then q1 = min(1,rw), if minimal = true, or q1 = rw, if minimal = false; if HDesign = H, then q1 is the row dimension of the design matrix H2. if HDesign2 = missing, then q2 = 1-min(1,rw), if minimal = true, or q2 = nvec-rw, if minimal = false, where nvec is the number of the nullspace basis vectors used for the initial synthesis (see [1]); if HDesign2 = H2, then q2 is the row dimension of the design matrix H2.

FDfreq = freq specifies a vector of real frequency values or a scalar real frequency value for strong detectability checks (default: FDfreq = missing).

If nullspace = true (default), a minimal proper nullspace basis is used for the initial synthesis of the fault detection filter. If nullspace = false, a full-order observer based nullspace basis is used. This option can be only used for a proper system without disturbance inputs.

If simple = true, a simple proper nullspace basis is emplyed for synthesis. The orders of the basis vectors employed for the synthesis are provided in info.deg and info.deg2 (see below). If simple = false (default), then no simple basis is 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).

scale2 = σ2 specifies the scaling factor σ2 to be employed for the components Q2(λ) and R2(λ) in (1), i.e., use σ2*Q2(λ) and σ2*R2(λ) instead of Q2(λ) and R2(λ). (default: σ2 is chosen to ensure the minimum gap provided by Q1(λ))

freq = val specifies the values of a test frequency to be employed to check the full row rank admissibility condition (default: randomly generated in the interval (0,1)).

tcond = tcmax specifies the maximum alowed condition number tcmax of the employed non-orthogonal transformations (default: tcmax = 1.e4).

FDtol = tol1 specifies the threshold tol1 for fault detectability checks (default: tol1 = 0.0001).

FDGainTol = tol2 specifies the threshold tol2 for strong fault detectability checks (default: tol2 = 0.01).

gamma = γ specifies the allowed upper bound for ∥Rw(λ)∥∞ (default: γ = 1).

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 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, degs2, S, S2, HDesign, HDesign2, freq, gap), 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 integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G1(λ) := [ Gu(λ) Gd(λ); I 0] (also the left Kronecker indices of G1(λ)), if the state-space realization of [Gu(λ) Gd(λ)] is minimal. This information has been used in the case minimal = true to determine the least order of components Q1(λ) and R1(λ) in (1).

info.degs2 is an integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G2(λ) := [ Gu(λ) Gd(λ) Gw(λ); I 0 0] (also the left Kronecker indices of G2(λ)), if the state-space realization of [Gu(λ) Gd(λ) Gw(λ)] is minimal. This information has been used in the case minimal = true to determine the least order of components Q2(λ) and R2(λ) in (1).

info.S is the binary structure matrix of the reduced system corresponding to the computed left nullspace basis of G1(λ) := [ Gu(λ) Gd(λ); I 0];

info.S2 is the binary structure matrix of the reduced system corresponding to the computed left nullspace basis of G2(λ) := [ Gu(λ) Gd(λ) Gw(λ); I 0 0];

info.HDesign is the design matrix H1 employed for the synthesis of the components Q1(λ) and R1(λ) in (1) of the fault detection filter;

info.HDesign2 is the design matrix H2 employed for the synthesis of the components Q2(λ) and R2(λ) in (1) of the fault detection filter;

info.freq is the frequency value employed to check the full row rank admissibility condition;

info.gap is the achieved gap ∥Rf(λ)∥∞−/∥Rw(λ)∥∞, where the H−minus index is computed over the whole frequency range, if FDfreq = missing, or over the frequency values contained in freq if FDfreq = freq.

Method: An extension of the Procedure AFD from [1] is implemented to solve the approximate fault detection problem (see also [2] and Remark 5.10 of [1]). The employed regularization approach, based on the modified co-outer-co-inner factorization, is discussed in [3], see also Remark 5.8 of [1]. For the details of the implemented method, see the documentation of the afdsyn function in [4].

References:

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

[2] A. Varga, General computational approach for optimal fault detection. Proc. IFAC Symposium SAFEPROCESS, Barcelona, Spain, pp. 107–112, 2009.

[3] K. Glover and A. Varga, On solving non-standard H-/H_2/inf fault detection problems. Proc. IEEE CDC, Orlando, FL, USA, pp. 891–896, 2011.

[4] A. Varga, Fault Detection and Isolation Tools (FDITOOLS) User's Guide, arXiv:1703.08480. 

afdsyn(sysf::FDIModel, SFDI; rdim, nullspace = true, simple = false, minimal = true, exact = false, 
                       gamma = 1, epsreg = 0.1, sdegzer, nonstd = 1, freq, sdeg, smarg, poles, 
                       HDesign, HDesign2, scale2, FDtol, FDGainTol, FDfreq, 
                       tcond, offset, atol, atol1, atol2, atol3, rtol, fast = true)  
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the approximate fault detection and isolation problem (AFDIP) for a given synthesis model sysf with additive faults and a given binary structure vector SFDI. The computed stable and proper filter objects Q and R contain the fault detection filter, representing the solution of the AFDIP, and its internal form, respectively, and are determined such that the transfer function matrix of R.sys[:,faults] has its j-th column nonzero if SFDI[j] = 1. If the solution of a strong AFDIP is feasible, then the j-th column is zero if SFDI[j] = 0. If only a the solution of a weak AFDIP is feasible, then the j-th column may be nonzero if SFDI[j] = 0. For the description of the keyword parameters see the function afdsyn.

FaultDetectionTools.afdisyn — Function
afdisyn(sysf::FDIModel, SFDI; rdim, nullspace = true, simple = false, minimal = true, separate = false,
                       gamma = 1, epsreg = 0.1, sdegzer, nonstd = 1, freq, sdeg, smarg, poles, 
                       HDesign, HDesign2, scale2, FDtol, FDGainTol, FDfreq, 
                       tcond, offset, atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the approximate fault detection and isolation problem (AFDIP) for a given synthesis model sysf with additive faults and a given binary structure matrix SFDI with nb rows (specifications). The computed stable and proper filter objects Q and R contain the fault detection and isolation filter, representing the solution of the AFDIP, and its internal form, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.degs2, info.HDesign, info.HDesign2 andinfo.gap contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The fault detection and isolation filter object Q, contains in Q.sys the resulting bank of nb 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_nb]. 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 indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

The fault detection and isolation filter in internal form object R, contains R.sys, the resulting bank of nb internal form of the filters. The i-th filter R.sys[i] is in a standard state-space form, which generates the residual signal r_i, and corresponds to the input-output form

   r_i = Rui(λ)*u + Rdi(λ)*d + Rfi(λ)*f + Rwi(λ)*w + Rai(λ)*aux ,

where 

   | Rui(λ) Rdi(λ) Rfi(λ) Rwi(λ) Rai(λ) | := |Qyi(λ) Qui(λ)]*| Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |. 
                                                             |   I     0     0     0     0   |

The solution of the AFDIP ensures that for the i-th filter, Rui(λ) = 0, Rdi(λ) = 0, Rfi(λ) has its j-th column nonzero if the (i,j)-th element of SFDI is nonzero, and the H∞-norm of Rwi(λ) satisfies ||Rwi(λ)||∞ < γ, where the bound γ is specified via the keyword argument gamma. The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls (void), R.disturbances (void), R.faults, R.noise and R.aux, respectively.

The transfer function matrices Qi(λ) := [ Qyi(λ) Qui(λ) ] and Ri(λ) := [ Rui(λ) Rdi(λ) Rfi(λ) Rwi(λ) Rai(λ) ] of the i-th components of the resulting filters Q.sys and R.sys, respectively, have, in general, the partitioned forms

 Qi(λ) = [ Q1i(λ) ] ,   Ri(λ) = [ R1i(λ) ] ,                      (1)
         [ Q2i(λ) ]             [ R2i(λ) ]

where the filters Q1i(λ) and R1i(λ) with q1i outputs are the solution of an AFDP, while the filters Q2i(λ) and R2i(λ) with q2i outputs are the solution of an exact fault detection problem formulated for a reduced system obtained by decoupling the control and disturbance inputs from the residuals (see [4]). The overall resulting component filters Q.sys[i] and R.sys[i] have observable state-space realizations (AQi,BQi,CQi,DQi) and (AQi,BRi,CQi,DRi), respectively, and thus share the observable pairs (AQi,CQi).

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

If minimal = true (default), least order filter synthesis is performed to determine each of the component filters Q.sys[i] and R.sys[i] for i = 1, ...,nb, while with minimal = false full order synthesis is performed.

If exact = true, an exact synthesis (without optimization) is performed, while with exact = false (default), an approximate synthesis is performed.

If HDesign = H1, then H1 is an nb-dimensional array of full row rank or empty design matrices, where H1[i] is the design matrix employed for the synthesis of the components Q1i(λ) and R1i(λ) in (1) of the i-th filter (default: HDesign = missing)

If HDesign2 = H2, then H2 is an nb-dimensional array of full row rank or empty design matrices, where H2[i] is the design matrix employed for the synthesis of the components Q2i(λ) and R2i(λ) in (1) of the i-th filter (default: HDesign2 = missing)

rdim = q specifies the vector q, whose i-th component q[i] specifies the number of residual outputs for the i-th component filters Q.sys[i] and R.sys[i]. If q is a scalar, then a vector rdim with all components equal to q is assumed. If rdim = missing, the default value of q[i] is chosen as q[i] = q1i + q2i, where the default values of q1i and q2i are chosen taking into account the rank rwi of Rwi(λ) in the reduced system (see [2]), as follows: if HDesign = missing, then q1i = min(1,rwi), if minimal = true, or q1i = rwi, if minimal = false; if HDesign = H1, then q1i is the row dimension of the nonemty design matrix H1[i], or if H1[i] is empty, the above choice for HDesign = missing is employed; if HDesign2 = missing, then q2i = 1-min(1,rwi), if minimal = true, or q2i is set to its maximum achievable value, if minimal = false (see [1]); if HDesign2 = H2, then q2i is the row dimension of the nonemty design matrix H2[i], or if H2[i] is empty, the above choice for HDesign2 = missing is employed.

FDfreq = freq specifies a vector of real frequency values or a scalar real frequency value for strong detectability checks (default: FDfreq = missing).

If nullspace = true (default), a minimal proper nullspace basis is used at the initial reduction step, if separate = false, or at all synthesis steps, if separate = true. If nullspace = false, a full-order observer based nullspace basis is used at the initial reduction step, if separate = false, or at all synthesis steps, if separate = true. This option can only be used for a proper system without disturbance inputs.

If simple = true, simple proper nullspace bases are emplyed for synthesis. If simple = false (default), then no simple bases are computed.

If separate = false (default), a two-step synthesis procedure is employed, where a minimal proper nullspace basis is used at the initial reduction step. If separate = true, the filter components are separately determined by solving appropriately formulated fault detection problems.

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).

FDtol = tol1 specifies the threshold tol1 for fault detectability checks (default: tol1 = 0.0001).

FDGainTol = tol2 specifies the threshold tol2 for strong fault detectability checks (default: tol2 = 0.01).

scale2 = σ2 specifies the vector of scaling factors σ2 to be employed for the components Q2i(λ) and R2i(λ) in (1), i.e., use σ2[i]*Q2i(λ) and σ2[i]*R2i(λ) instead of Q2i(λ) and R2i(λ). (default: For scale2 = missing, each σ2[i] is chosen to ensure the minimum gap provided by Q1i(λ))

gamma = γ specifies the allowed upper bound for the resulting ∥Rwi(λ)∥∞ (default: γ = 1).

epsreg = ϵ specifies the value of the regularization parameter ϵ used in afdsyn (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 in afdsyn, 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 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, HDesign, HDesign2, freq, gap), where:

info.tcond is an nb-dimensional vector, whose i-th component is the maximum of the condition numbers of the employed non-orthogonal transformation matrices employed for the synthesis of the i-th filter component; a warning is issued if any info.tcond[i] >= tcmax;

info.HDesign = H1 is an nb-dimensional vector of design matrices, whose i-th component H1[i] is the design matrix to be employed for the synthesis of the components Q1i(λ) and R1i(λ) in (1) of the i-th fault detection filter.

info.HDesign2 = H2 is an nb-dimensional vector of design matrices, whose i-th component H2[i] is the design matrix to be employed for the synthesis of the components Q2i(λ) and R2i(λ) in (1) of the i-th fault detection filter.

info.freq is the frequency value employed to check the full row rank admissibility condition.

info.gap is an nb-dimensional vector, whose i-th component is the achieved gap for the synthesis of the i-th filter component.

Method: The Procedure AFDI from [1] is implemented to solve the approximate fault detection and isolation problem. For implementation details, see [2].

References:

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

[2] A. Varga, Fault Detection and Isolation Tools (FDITOOLS) User's Guide, arXiv:1703.08480.

FaultDetectionTools.emmsyn — Function
emmsyn(sysf::FDIModel, sysr::FDFilterIF; nullspace = true, simple = false, minimal = true, regmin = true, normalize = "gain", 
                       sdeg, smarg, poles, freq, HDesign, tcond, offset, 
                       atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the exact model-matching problem (EMMP) for a given synthesis model sysf::FDIModel with additive faults and a given stable reference filter sysr::FDFilterIF. The computed stable and proper filter objects Q and R contain the fault detection filter, representing the solution of the EMMP, and its internal form, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.M, info.freq and info.HDesign, contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The continuous- or discrete-time reference filter sysr.sys is in a standard or descriptor state-space form sysr.sys = (Ar-λEr,Br,Cr,Dr), which corresponds to the input-output form 

   yr = Mru(λ)*u + Mrd(λ)*d + Mrf(λ)*f + Mrw(λ)*w + Mra(λ)*aux,

with the Laplace- or Z-transformed reference filter outputs yr, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Mru(λ), Mrd(λ), Mrf(λ), Mrw(λ), and Mra(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysr.controls, sysr.disturbances, sysr.faults, sysr.noise and sysr.aux, respectively. If any of the above vectors is void, then the corresponding transfer function matrix is considered null.

The fault detection filter object Q, contains in Q.sys the resulting filter in a standard state-space form, which generates the residual signal r. The corresponding input-output (implementation) form is

        r = Qy(λ)*y + Qu(λ)*u

where Qy(λ) and Qu(λ) are the transfer function matrices from the output and control inputs to the residual. The indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

The fault detection filter internal form object R, contains R.sys, the resulting internal form of the filter in a standard state-space form, which generates the residual signal r, and corresponds to the input-output form

   r = Ru(λ)*u + Rd(λ)*d + Rf(λ)*f + Rw(λ)*w + Ra(λ)*aux ,

where 

   | Ru(λ) Rd(λ) Rf(λ) Rw(λ) Ra(λ) | = |Qy(λ) Qu(λ)|*| Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |. 
                                                     |  I     0     0     0     0    |

The solution of the standard EMMP is computed if sysr.noise and sysr.aux are void and ensures that Ru(λ) = M(λ)*Mru(λ), Rd(λ) = M(λ)*Mrd(λ) and Rf(λ) = M(λ)*Mrf(λ), where M(λ) is the transfer function matrix of a stable, diagonal and invertible updating filter returned in info.M. This filter is determined to guarantee the stability of resulting filters Q and R. If sysr.noise and sysr.aux are not both void, then the extended EMMP is solved which additionally ensures Rw(λ) = M(λ)*Mrw(λ) and Ra(λ) = M(λ)*Mra(λ). The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls, R.disturbances, R.faults, R.noise and R.aux, respectively.

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

If minimal = true (default), a least order filter synthesis is performed, while with minimal = false no least order synthesis is performed.

If regmin = true (default), the regularization (see [1]) is performed for the case when sysr.controls and sysr.disturbances are void with the selection of a least order left annihilator Nl(λ) such that Nl(λ)*[Gu(λ) Gd(λ); I 0 ]. If regmin = false, the regularization is performed by choosing Nl(λ) a minimal left nullspace basis of G(λ) = [Gu(λ) Gd(λ); I 0 ].

If HDesign = H is a full row rank design matrix, then H*Nl(λ) is used instead Nl(λ) (default: HDesign = missing).

An initial reduction step is performed using the nullspace-based approach (see [1]) if sysr.controls, sysr.disturbances, sysr.noise and sysr.aux are void and minimal = false. In this case, if nullspace = true (default), a minimal proper nullspace basis is used at the initial reduction step, while, if nullspace = false, a full-order observer based nullspace basis is used at the initial reduction step. This later option can only be used for a proper system without disturbance inputs. The nullspace option is ignored if any of sysr.controls, sysr.disturbances, sysr.noise or sysr.aux is non-void or if minimal = true

If simple = true, a simple proper nullspace basis Nl(λ) is emplyed as left annihilator for synthesis. The orders of the basis vectors are provided in info.deg. If simple = false (default), then a minimal proper nullspace basis is 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).

freq = val specifies the value of a test frequency to be employed to check the full column rank (i.e., left-invertibility) solvability condition (default: randomly generated in the interval (0,1)). The employed value of freq is returned in info.freq.

normalize = job specifies the option for the normalization of the diagonal elements of the updating matrix M(λ) as follows:

  job = "gain"    – scale with the gains of the zero-pole-gain representation (default);
  job = "dcgain"  – scale with the DC-gains;
  job = "infnorm" – scale with the values of infinity-norms.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, M, freq, 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 integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G(λ) := [ Gu(λ) Gd(λ); I 0] (also the left Kronecker indices of G(λ)), if the state-space realization of [Gu(λ) Gd(λ)] is minimal;

info.M is the employed stable and invertible updating filter used to solve the EMMP, with a diagonal transfer function matrix M(λ);

info.freq is the employed frequency used to check left invertibility (set to missing if no frequency-based left invertibility check was performed)

info.HDesign is the design matrix H employed for the synthesis of the fault detection filter Q; H = missing if no design matrix was involved.

Method: The synthesis Procedures EMM and EMMS from [1] are implemented. Procedure EMM relies on the model-matching synthesis method proposed in [2], while Procedure EMMS uses the inversion-based method proposed in [3]. Procedure EMM is generally employed, unless a strong exact fault detection and isolation problem (strong EFDIP) is solved, in which case Procedure EMMS is used.

The strong EFDIP corresponds to the choice of the reference filter sysr such that Mru(λ) = 0, Mrd(λ) = 0, Mrf(λ) is invertible, Mrw(λ) = 0 and Mra(λ) = 0. In this case, only the indices of fault inputs sysr.faults must be specified and the indices of the rest of inputs must be void. The solution of a fault estimation problem can be targeted by choosing Mrf(λ) = I and checking that the resulting info.M = I.

References:

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

[2] A. Varga, New computational approach for the design of fault detection and isolation filters. In M. Voicu (Ed.), "Advances in Automatic Control", vol. 754 of The Kluwer International Series in Engineering and Computer Science, Kluwer Academic Publishers, pp. 367-381, 2003.

[3] A. Varga. New computational paradigms in solving fault detection and isolation problems. Annual Reviews in Control, 37:25–42, 2013. 

emmsyn(sysf::FDIModel, sysr::FDIFilterIF; nullspace = true, simple = false, minimal = true, regmin = true, normalize = "gain", 
                       sdeg, smarg, poles, freq, HDesign, tcond, offset, 
                       atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDIFilter, R::FDIFilterIF, info)

Solve the exact model-matching problem (EMMP) for a given synthesis model sysf::FDIModel with additive faults and a given bank of stable reference filters sysr::FDIFilterIF. The computed stable and proper filter objects Q and R contain the bank of fault detection filters, representing the component-wise solution of the EMMP, and their internal forms, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.M, info.freq and info.HDesign, contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The continuous- or discrete-time reference filters packed in sysr are in a standard or descriptor state-space form, where the i-th filter sysr.sys[i] = (Ari-λEri,Bri,Cri,Dri) corresponds to the input-output form 

   yri = Mrui(λ)*u + Mrdi(λ)*d + Mrfi(λ)*f + Mrwi(λ)*w + Mrai(λ)*aux,

with the Laplace- or Z-transformed reference filter outputs yri, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Mrui(λ), Mrdi(λ), Mrfi(λ), Mrwi(λ), and Mrai(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysr.controls, sysr.disturbances, sysr.faults, sysr.noise and sysr.aux, respectively. If any of the above vectors is void, then the corresponding transfer function matrices are considered null.

The fault detection and isolation filter object Q, contains in its i-th component Q.sys[i] the resulting filter in a standard state-space form, which generates the i-th component ri of the residual signal. The corresponding input-output (implementation) form is

        ri = Qyi(λ)*y + Qui(λ)*u

where Qyi(λ) and Qui(λ) are the transfer function matrices from the output and control inputs to the residual. The indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

The fault detection and isolation filter internal form object R, contains in its i-th component R.sys[i], the resulting internal form of the filter in a standard state-space form, which generates the i-th component ri of residual signal , and corresponds to the input-output form

   ri = Rui(λ)*u + Rdi(λ)*d + Rfi(λ)*f + Rwi(λ)*w + Rai(λ)*aux ,

where 

   | Rui(λ) Rdi(λ) Rfi(λ) Rwi(λ) Rai(λ) | = |Qyi(λ) Qui(λ)|*| Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |. 
                                                            |  I     0     0     0     0    |

The component-wise solution of the standard EMMP is computed if sysr.noise and sysr.aux are void and ensures that Rui(λ) = Mi(λ)*Mrui(λ), Rdi(λ) = Mi(λ)*Mrdi(λ) and Rfi(λ) = Mi(λ)*Mrfi(λ), where Mi(λ) is the transfer function matrix of a stable, diagonal and invertible updating filter returned in the i-th component of the vector info.M. This filter is determined to guarantee the stability of the i-th components of resulting filters Q and R. If sysr.noise and sysr.aux are not both void, then the extended EMMP is component-wise solved which additionally ensures Rwi(λ) = Mi(λ)*Mrwi(λ) and Rai(λ) = Mi(λ)*Mrai(λ). The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls, R.disturbances, R.faults, R.noise and R.aux, respectively.

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

If minimal = true (default), least order filter syntheses are performed, while with minimal = false no least order synthesis are performed.

If regmin = true (default), the regularization (see [1]) is performed for the case when sysr.controls and sysr.disturbances are void with the selection of a least order left annihilator Nl(λ) such that Nl(λ)*[Gu(λ) Gd(λ); I 0 ]. If regmin = false, the regularization is performed by choosing Nl(λ) a minimal left nullspace basis of G(λ) = [Gu(λ) Gd(λ); I 0 ].

If HDesign = H is a vector of full row rank design matrices, then H[i]*Nl(λ) is used instead Nl(λ) for the synthesis of the i-th filter (default: HDesign = missing).

An initial reduction step is performed using the nullspace-based approach (see [1]) if sysr.controls, sysr.disturbances, sysr.noise and sysr.aux are void and minimal = false. In this case, if nullspace = true (default), a minimal proper nullspace basis is used at the initial reduction step, while, if nullspace = false, a full-order observer based nullspace basis is used at the initial reduction step. This later option can only be used for a proper system without disturbance inputs. The nullspace option is ignored if any of sysr.controls, sysr.disturbances, sysr.noise or sysr.aux is non-void or if minimal = true

If simple = true, a simple proper nullspace basis Nl(λ) is emplyed as left annihilator for synthesis. The orders of the basis vectors are provided in info.deg. If simple = false (default), then a minimal proper nullspace basis is 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).

freq = val specifies the value of a test frequency to be employed to check the full column rank (i.e., left-invertibility) solvability condition (default: randomly generated in the interval (0,1)). The employed value of freq is returned in info.freq.

normalize = job specifies the option for the normalization of the diagonal elements of the updating matrices Mi(λ) as follows:

  job = "gain"    – scale with the gains of the zero-pole-gain representation (default);
  job = "dcgain"  – scale with the DC-gains;
  job = "infnorm" – scale with the values of infinity-norms.

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, M, freq, 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 integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G(λ) := [ Gu(λ) Gd(λ); I 0] (also the left Kronecker indices of G(λ)), if the state-space realization of [Gu(λ) Gd(λ)] is minimal;

info.M is a vector of descriptor systems, whose i-th system info.M[i] contains the employed stable and invertible updating filter used to solve the i-th EMMP, with a diagonal transfer function matrix Mi(λ);

info.freq is the employed frequency used to check left invertibility (set to missing if no frequency-based left invertibility check was performed)

info.HDesign is a vector of design matrices H, where H[i] is the design matrix employed for the synthesis of the i-th component of the fault detection filter Q; H[i] is an empty matrix if no design matrix was involved.

Method: The synthesis Procedures EMM and EMMS from [1] are used to determine the component filters. Procedure EMM relies on the model-matching synthesis method proposed in [2], while Procedure EMMS uses the inversion-based method proposed in [3]. Procedure EMM is generally employed, unless a strong exact fault detection and isolation problem (strong EFDIP) is solved, in which case Procedure EMMS is used.

The strong EFDIP corresponds to the choice of each component of the bank of reference filters sysr such that Mrui(λ) = 0, Mrdi(λ) = 0, Mrfi(λ) is invertible, Mrwi(λ) = 0 and Mrai(λ) = 0. In this case, only the indices of fault inputs sysr.faults must be specified and the indices of the rest of inputs must be void. The solution of a fault estimation problem can be targeted by choosing Mrfi(λ) = ei, where ei is the i-th row of the appropriate identity matrix, and checking that the resulting info.M = 1.

References:

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

[2] A. Varga, New computational approach for the design of fault detection and isolation filters. In M. Voicu (Ed.), "Advances in Automatic Control", vol. 754 of The Kluwer International Series in Engineering and Computer Science, Kluwer Academic Publishers, pp. 367-381, 2003.

[3] A. Varga. New computational paradigms in solving fault detection and isolation problems. Annual Reviews in Control, 37:25–42, 2013.

FaultDetectionTools.ammsyn — Function
ammsyn(sysf::FDIModel, sysr::FDFilterIF; nullspace = true, simple = false, mindeg = false, 
                       regmin = true, normalize = "infnorm", H2syn = false, reltol = 1.e-4, 
                       sdeg, smarg, poles, freq, HDesign, tcond, offset, 
                       atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the approximate model-matching problem (AMMP) for a given synthesis model sysf::FDIModel with additive faults and a given stable reference filter sysr::FDFilterIF. The computed stable and proper filter objects Q and R contain the fault detection filter, representing a solution of the AMMP, and its internal form, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.M, info.freq, info.HDesign, info.nonstandard, info.gammaopt0, info.gammaopt and info.gammasub, contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The continuous- or discrete-time reference filter sysr.sys is in a standard or descriptor state-space form sysr.sys = (Ar-λEr,Br,Cr,Dr), which corresponds to the input-output form 

   yr = Mru(λ)*u + Mrd(λ)*d + Mrf(λ)*f + Mrw(λ)*w + Mra(λ)*aux,

with the Laplace- or Z-transformed reference filter outputs yr, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Mru(λ), Mrd(λ), Mrf(λ), Mrw(λ), and Mra(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysr.controls, sysr.disturbances, sysr.faults, sysr.noise and sysr.aux, respectively. If any of the above vectors is void, then the corresponding transfer function matrix is considered null.

The fault detection filter object Q, contains in Q.sys the resulting filter in a standard state-space form, which generates the residual signal r. The corresponding input-output (implementation) form is

        r = Qy(λ)*y + Qu(λ)*u

where Qy(λ) and Qu(λ) are the transfer function matrices from the output and control inputs to the residual. The indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

Let define 

  Ge(λ) = | Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |,  Mr(λ) = | Mru(λ) Mrd(λ) Mrf(λ) Mrw(λ) Mra(λ) | . 
          |  I     0     0     0     0    |

In the standard case, Ge(λ) has no zeros on the boundary of the stability domain, and the resulting stable filter Q(λ) := |Qy(λ) Qu(λ)| is Q(λ) = Q0(λ), where Q0(λ) is the optimal solution of the H∞- or H2-norm error minimization problem

gammaopt0 = ||Q0(λ)*Ge(λ)-M(λ)*Mr(λ)|| = min,         (1)

where M(λ) = M0(λ) is an updating factor chosen as follows: M0(λ) = I in the case of emplyoing the H∞ norm, while in the case of employing the H2 norm, M0(λ) = I for a discrete-time system or, for a continuous-time system, M0(λ) is determined a stable, diagonal, and invertible transfer function matrix, which ensures the existence of a finite H2-norm.

In the non-standard case, Ge(λ) has zeros on the boundary of the stability domain, and the resulting optimal filter Q0(λ), which solves the H∞- or H2-norm error minimization problem (1) is a possibly unstable or improper. A second updating factor M1(λ) is determined, with the same properties as M0(λ), which ensures that the computed stable and proper filter Q(λ) := M1(λ)*Q0(λ) represents a suboptimal solution of an updated H∞- or H2-norm error minimization problem, for which the achieved suboptimal model-matching performance is

gammasub = ||Q(λ)*Ge(λ)-M(λ)*Mr(λ)|| ,            (2)

where M(λ) := M1(λ)*M0(λ). The optimal solution Qt(λ) of the updated H∞- or H2-norm error minimization problem 

gammaopt = ||Qt(λ)*Ge(λ)-M(λ)*Mr(λ)|| = min ,     (3)

is still possibly unstable or improper. The values of gammaopt0, gammaopt and gammasub are returned in info.gammaopt0, info.gammaopt and info.gammasub, respectively.

The fault detection filter internal form object R, contains R.sys, the resulting internal form of the filter in a standard state-space form, which generates the residual signal r, and corresponds to the input-output form

   r = Ru(λ)*u + Rd(λ)*d + Rf(λ)*f + Rw(λ)*w + Ra(λ)*aux ,

where 

   | Ru(λ) Rd(λ) Rf(λ) Rw(λ) Ra(λ) | = Q(λ)*Ge(λ).

The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls, R.disturbances, R.faults, R.noise and R.aux, respectively. The state-space realization of the resulting M(λ) is returned in info.M.

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

If H2syn = false (default), a H∞-norm based synthesis is performed, while if H2syn = true, a H2-norm based synthesis is performed.

reltol = tol specifies the relative tolerance tol for the desired accuracy of γ-iteration (default: tol = 1.e-4).

If mindeg = true, a least order filter synthesis is performed, if possible, while with minimal = false (default) no least order synthesis is performed.

If regmin = true (default), the regularization (see [1]) is performed for the case when sysr.controls and/or sysr.disturbances are void with the selection of a least order left annihilator Nl(λ) of G(λ) = [Gu(λ) Gd(λ); I 0 ]. If regmin = false, the regularization is performed by choosing a left annihilator Nl(λ) as a minimal left nullspace basis of G(λ).

If HDesign = H is a full row rank design matrix, then H*Nl(λ) is used as left annihilator instead Nl(λ) (default: HDesign = missing).

If nullspace = true (default) and sysr.controls and/or sysr.disturbances are void, a minimal proper nullspace basis is used at the initial reduction step. If nullspace = false and sysr.controls and/or sysr.disturbances are void, a full-order observer based nullspace basis is used at the initial reduction step. This option can only be used for a proper system without disturbance inputs. The nullspace option is ignored if both sysr.controls and sysr.disturbances are non-void.

If simple = true, a simple proper nullspace basis is emplyed as left annihilator for synthesis. The orders of the basis vectors are provided in info.deg. If simple = false (default), then a minimal proper nullspace basis is 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).

freq = val specifies the value of a test frequency to be employed to check the full column rank (i.e., left-invertibility) solvability condition (default: randomly generated in the interval (0,1)). The employed value of freq is returned in info.freq.

normalize = job specifies the option for the normalization of the diagonal elements of the updating matrix M(λ) as follows:

  job = "gain"    – scale with the gains of the zero-pole-gain representation;
  job = "dcgain"  – scale with the DC-gains;
  job = "infnorm" – scale with the values of infinity-norms (default).

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, M, freq, HDesign, gammaopt0, gammaopt, gammasub, 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 integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G(λ) := [ Gu(λ) Gd(λ); I 0] (also the left Kronecker indices of G(λ)), if the state-space realization of [Gu(λ) Gd(λ)] is minimal;

info.M is the employed stable and invertible updating filter used to solve the AMMP, with a stable, diagonal and invertible transfer function matrix M(λ);

info.freq is the employed frequency used to check left invertibility (set to missing if no frequency-based left invertibility check was performed)

info.HDesign is the design matrix H employed for the synthesis of the fault detection filter Q; H = missing if no design matrix was involved;

info.gammaopt0 is the optimal performance gammaopt0 for the original problem (1);

info.gammaopt is the optimal performance gammaopt for the updated problem (3);

info.gammasub is the suboptimal performance gammasub in (2);

info.nonstandard is set to true for a non-standard problem (i.e., Ge(λ) has zeros on the boundary of the stability domain), and set to false for a standard problem (i.e., Ge(λ) has no zeros on the boundary of the stability domain).

Method: The synthesis Procedure AMMS from [1] is implemented. The Procedure AMMS relies on the approximate model-matching synthesis method proposed in [2]. For more details on computational aspects see [3].

References:

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

[2] A. Varga, Integrated computational algorithm for solving H_inf-optimal FDI problems. In Proc. of the IFAC World Congress, Milano, Italy, pp. 10187–10192, 2011.

[3] A. Varga. Descriptor system techniques in solving H_2/H-Inf-optimal fault detection and isolation problems". In L. T. Biegler, S. L. Campbell, and V. Mehrmann (Eds.), Control and Optimization with Differential-Algebraic Constraints, vol. 23 of Advances in Design and Control, pp. 105–125. SIAM, 2012. 

ammsyn(sysf::FDIModel, sysr::FDIFilterIF; nullspace = true, simple = false, mindeg = false, 
                       regmin = true, normalize = "infnorm", H2syn = false, reltol = 1.e-4, 
                       sdeg, smarg, poles, freq, HDesign, tcond, offset, 
                       atol, atol1, atol2, atol3, rtol, fast = true) 
                       -> (Q::FDFilter, R::FDFilterIF, info)

Solve the approximate model-matching problem (AMMP) for a given synthesis model sysf::FDIModel with additive faults and a given bank of stable reference filters sysr::FDIFilterIF. The computed stable and proper filter objects Q and R contain the bank of fault detection filters, representing the component-wise solution of the AMMP, and their internal forms, respectively.

The returned named tuple info, with the components info.tcond, info.degs, info.M, info.freq, info.HDesign, info.nonstandard, info.gammaopt0, info.gammaopt and info.gammasub, contains additional synthesis related information (see below).

The continuous- or discrete-time system sysf.sys is in a standard or descriptor state-space form sysf.sys = (A-λE,B,C,D), which corresponds to the input-output form 

   y = Gu(λ)*u + Gd(λ)*d + Gf(λ)*f + Gw(λ)*w + Ga(λ)*aux,

with the Laplace- or Z-transformed plant outputs y, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Gu(λ), Gd(λ), Gf(λ), Gw(λ), and Ga(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysf.controls, sysf.disturbances, sysf.faults, sysf.noise and sysf.aux, respectively.

The continuous- or discrete-time reference filters packed in sysr are in a standard or descriptor state-space form, where the i-th filter sysr.sys[i] = (Ari-λEri,Bri,Cri,Dri) corresponds to the input-output form 

   yri = Mrui(λ)*u + Mrdi(λ)*d + Mrfi(λ)*f + Mrwi(λ)*w + Mrai(λ)*aux,

with the Laplace- or Z-transformed reference filter outputs yri, control inputs u, disturbance inputs d, fault inputs f, noise inputs w and auxiliary inputs aux, and with Mrui(λ), Mrdi(λ), Mrfi(λ), Mrwi(λ), and Mrai(λ) the corresponding transfer-function matrices. The indices of control, disturbance, fault, noise and auxiliary inputs are contained in the associated integer vectors sysr.controls, sysr.disturbances, sysr.faults, sysr.noise and sysr.aux, respectively. If any of the above vectors is void, then the corresponding transfer function matrices are considered null.

The fault detection and isolation filter object Q, contains in its i-th component Q.sys[i] the resulting filter in a standard state-space form, which generates the i-th component ri of the residual signal. The corresponding input-output (implementation) form is

        ri = Qyi(λ)*y + Qui(λ)*u

where Qyi(λ) and Qui(λ) are the transfer function matrices from the output and control inputs to the residual. The indices of output and control inputs are contained in the integer vectors Q.outputs and Q.controls, respectively.

Let define 

  Ge(λ) = | Gu(λ) Gd(λ) Gf(λ) Gw(λ) Ga(λ) |,  Mri(λ) = | Mrui(λ) Mrdi(λ) Mrfi(λ) Mrwi(λ) Mrai(λ) | . 
          |  I     0     0     0     0    |

In the standard case, Ge(λ) has no zeros on the boundary of the stability domain, and each resulting component stable filter Qi(λ) := |Qyi(λ) Qui(λ)| is Qi(λ) = Q0i(λ), where Q0i(λ) is the optimal solution of the H∞- or H2-norm error minimization problem

gammaopt0i = ||Q0i(λ)*Ge(λ)-Mi(λ)*Mri(λ)|| = min,         (1)

where Mi(λ) = M0i(λ) is an updating factor chosen as follows: M0i(λ) = I in the case of emplyoing the H∞ norm, while in the case of employing the H2 norm, M0i(λ) = I for a discrete-time system or, for a continuous-time system, M0i(λ) is determined a stable, diagonal, and invertible transfer function matrix, which ensures the existence of a finite H2-norm.

In the non-standard case, Ge(λ) has zeros on the boundary of the stability domain, and each resulting optimal component filter Q0i(λ), which solves the H∞- or H2-norm error minimization problem (1) is a possibly unstable or improper. A second updating factor M1i(λ) is determined, with the same properties as M0i(λ), which ensures that the computed stable and proper filter Qi(λ) := M1i(λ)*Q0i(λ) represents a suboptimal solution of an updated H∞- or H2-norm error minimization problem, for which the achieved suboptimal model-matching performance is

gammasubi = ||Qi(λ)*Ge(λ)-Mi(λ)*Mri(λ)|| ,            (2)

where Mi(λ) := M1i(λ)*M0i(λ). The optimal solutions Qti(λ) of the updated H∞- or H2-norm error minimization problem 

gammaopti = ||Qti(λ)*Ge(λ)-Mi(λ)*Mri(λ)|| = min ,     (3)

is still possibly unstable or improper. The values of gammaopt0i, gammaopti and gammasubi are returned in the i-th components of the vectors info.gammaopt0, info.gammaopt and info.gammasub, respectively.

The fault detection and isolation filter internal form object R, contains in its i-th component R.sys[i], the resulting internal form of the filter in a standard state-space form, which generates the i-th component ri of residual signal , and corresponds to the input-output form

   ri = Rui(λ)*u + Rdi(λ)*d + Rfi(λ)*f + Rwi(λ)*w + Rai(λ)*aux ,

where 

   | Rui(λ) Rdi(λ) Rfi(λ) Rwi(λ) Rai(λ) | = Qi(λ)*Ge(λ).

The indices of the inputs u, d, f, w and aux of the resulting filter R.sys are contained in the integer vectors R.controls, R.disturbances, R.faults, R.noise and R.aux, respectively. The state-space realization of the resulting Mi(λ) is returned in the i-th component of the vector info.M.

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

If H2syn = false (default), a H∞-norm based synthesis is performed, while if H2syn = true, a H2-norm based synthesis is performed.

reltol = tol specifies the relative tolerance tol for the desired accuracy of γ-iteration (default: tol = 1.e-4).

If mindeg = true, least order filter syntheses are performed, if possible, while with minimal = false (default) no least order synthesis are performed.

If regmin = true (default), the regularization (see [1]) is performed for the case when sysr.controls and/or sysr.disturbances are void with the selection of a least order left annihilator Nl(λ) of G(λ) = [Gu(λ) Gd(λ); I 0 ]. If regmin = false, the regularization is performed by choosing a left annihilator Nl(λ) as a minimal left nullspace basis of G(λ).

If HDesign = H is a vector of full row rank design matrices, then H[i]*Nl(λ) is used as left annihilator instead Nl(λ) for the synthesis of the i-th filter (default: HDesign = missing).

If nullspace = true (default) and sysr.controls and/or sysr.disturbances are void, a minimal proper nullspace basis is used at the initial reduction step. If nullspace = false and sysr.controls and/or sysr.disturbances are void, a full-order observer based nullspace basis is used at the initial reduction step. This option can only be used for a proper system without disturbance inputs. The nullspace option is ignored if both sysr.controls and sysr.disturbances are non-void.

If simple = true, a simple proper nullspace basis is emplyed as left annihilator for synthesis. The orders of the basis vectors are provided in info.deg. If simple = false (default), then a minimal proper nullspace basis is 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).

freq = val specifies the value of a test frequency to be employed to check the full column rank (i.e., left-invertibility) solvability condition (default: randomly generated in the interval (0,1)). The employed value of freq is returned in info.freq.

normalize = job specifies the option for the normalization of the diagonal elements of the updating matrices Mi(λ) as follows:

  job = "gain"    – scale with the gains of the zero-pole-gain representation;
  job = "dcgain"  – scale with the DC-gains;
  job = "infnorm" – scale with the values of infinity-norms (default).

The rank determinations in the performed reductions are based on rank revealing QR-decompositions with column pivoting if fast = true 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 A, B, C, D, the absolute tolerance for the nonzero elements of E, and the relative tolerance for the nonzero elements of A, B, C, D and E. The default relative tolerance is n*ϵ, where ϵ is the working machine epsilon and n is the order of the system sysf.sys. 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, M, freq, HDesign, gammaopt0, gammaopt, gammasub, 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 integer vector containing the increasingly ordered degrees of a left minimal polynomial nullspace basis of G(λ) := [ Gu(λ) Gd(λ); I 0] (also the left Kronecker indices of G(λ)), if the state-space realization of [Gu(λ) Gd(λ)] is minimal;

info.M is a vector of descriptor systems, whose i-th system info.M[i] contains the employed stable and invertible updating filter used to solve the i-th AMMP, with a diagonal transfer function matrix Mi(λ);

info.freq is the employed frequency used to check left invertibility (set to missing if no frequency-based left invertibility check was performed)

info.HDesign is a vector of design matrices H, where H[i] is the design matrix employed for the synthesis of the i-th component of the fault detection filter Q; H[i] is an empty matrix if no design matrix was involved.

info.gammaopt0 is a vector whose i-th component is the optimal performance gammaopt0i for the i-th original problem (1);

info.gammaopt is a vector whose i-th component is the optimal performance gammaopti for the i-th updated problem (3);

info.gammasub is a vector whose i-th component is the suboptimal performance gammasubi in (2);

info.nonstandard is set to true for a non-standard problem (i.e., Ge(λ) has zeros on the boundary of the stability domain), and set to false for a standard problem (i.e., Ge(λ) has no zeros on the boundary of the stability domain).

Method: The synthesis Procedure AMMS from [1] is used to determine the component filters. The Procedure AMMS relies on the approximate model-matching synthesis method proposed in [2]. For more details on computational aspects see [3].

References:

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

[2] A. Varga, Integrated computational algorithm for solving H_inf-optimal FDI problems. In Proc. of the IFAC World Congress, Milano, Italy, pp. 10187–10192, 2011.

[3] A. Varga. Descriptor system techniques in solving H_2/H-Inf-optimal fault detection and isolation problems". In L. T. Biegler, S. L. Campbell, and V. Mehrmann (Eds.), Control and Optimization with Differential-Algebraic Constraints, vol. 23 of Advances in Design and Control, pp. 105–125. SIAM, 2012.