A Trisection Algorithm for Estimating Distance Measures for Strong Observability and Strong Detectability

For reliably analyzing the properties of strong observability and strong detectability of a system, continuous distance measures can be used. When calculating these measures, it is necessary to find the global minimum of a nonconvex target function. The main contribution of this article is an optimization algorithm that guarantees to find this global minimum in a fast and efficient way by exploiting the special structure of the optimization problem. Using this optimization algorithm, the distance measures can be reliably calculated. The numerical properties and the usefulness of the algorithm in practical applications are illustrated by means of a numerical example.


I. INTRODUCTION
In real-world control applications, unknown external disturbances are often present. The correct estimation of internal system states in presence of disturbances and the reconstruction of these disturbances is of particular importance in applications, like fault detection or robust control. Many external disturbances can be modeled as unknown inputs. For state estimation in presence of unknown inputs, strong observability or strong detectability, as defined in [1], are essential system properties. Standard criteria for these properties have binary character and their reliable numerical evaluation is difficult, but important, for example, in sensor placement [2]- [4]. Therefore, measures for reliable evaluation of strong observability and strong detectability are of major interest.
For the evaluation of standard observability, or its dual controllability, of a system, a similar problem has been extensively studied. Several continuous measures for observability and controllability have been proposed, see [5] for an overview. These continuous measures can be categorized into modal, energy, and distance measures. Modal measures have been introduced in [6], and give information about the controllability or observability of every mode separately. Energy measures analyze how much energy is necessary to bring the system into a desired state, see [7]. Distance measures have been introduced in [8], and yield the smallest perturbation in system matrices rendering the system uncontrollable or unobservable.
In [9], the idea of distance measures was extended to systems with unknown inputs and distance measures, for strong observability and strong detectability, was introduced. These measures give the distance to the closest non-strongly observable or non-strongly detectable system, which is calculated by minimizing the smallest singular value of the Rosenbrock matrix. The target function of this minimization problem is nonconvex and can have multiple local minima. However, for the correct determination of the measures, knowledge of the global minimum is absolutely necessary.
For the calculation of the distance measures for standard controllability and observability, special algorithms have been proposed. In [10], an iterative algorithm is presented, where estimations for the smallest singular value and the corresponding singular vector were refined in every iteration step. This algorithm only finds the local minimum and the results depend on the initial conditions. In [11], an algorithm for finding the minimum along a straight line using a bisection algorithm was proposed. In combination with an appropriate grid, this algorithm is able to find the global minimum. The algorithm in [12] uses a similar approach, but directly determines the global minimum by a trisection algorithm with the so-called Gu's test as condition. In [13], a quasi-Newton optimization is presented that finds a minimum and uses Gu's test to check if it is the global minimum.
For the calculation of the distance measures for strong observability and strong detectability, as presented in [9], special algorithms have not been studied yet. Several standard optimization algorithms can be used, for example, grid search [14], quasi-Newton optimization [13], [15], or genetic algorithms [16]. However, these optimization algorithms do not guarantee to find the global minimum, and high computational costs are necessary for increasing the degree of confidence in a solution. A concept closely related to the distance measures for strong observability and strong detectability is the minimum phase radius, as defined in [17]. In [18], an algorithm for computing this radius is presented, which also does not guarantee to find the global minimum of the considered target function.
In this article, an optimization algorithm for the reliable and efficient calculation of measures for strong observability and strong detectability is presented, which guarantees to find the global minimum of the target function. This algorithm is based on a trisection approach in combination with Gu's test, as presented in [12] and [19], for standard controllability, and extends these ideas to the strong observability case. The properties of the algorithm are discussed and a numerical example is given, in order to show its practical usefulness.
This article is structured as follows. Section II gives preliminaries. Section III introduces the theoretical foundation for the development of the optimization algorithm. Section IV gives the optimization algorithms based on the previous theorems. In Section V, a numerical example is presented. with states x ∈ R n , unknown inputs w ∈ R q , outputs y ∈ R p , and A, E, C, and F being constant matrices of appropriate dimensions. It is assumed that no prior knowledge about the unknown inputs w is available, that rank E = q, and that more measurements than unknown inputs are available (p ≥ q). Without loss of generality, any known input is omitted. The matrices describing system Σ are combined in the matrix The operators M H , M , and σ min (M) denote the complex conjugate transpose, the spectral norm, and the smallest singular value of a matrix M, respectively. The complex conjugate of c ∈ C is denoted byc and its real part by Re(c). Furthermore, the complex right open-half-plane is denoted by C + = {c ∈ C|Re(c) > 0}.

A. Strong Observability and Strong Detectability
Strong observability or strong detectability were introduced in [1], and are necessary system properties for state estimation in the presence of unknown inputs. These properties are connected to the transmission zeros of a system. Transmission zeros are the values λ 0 ∈ C, for which the matrix satisfies rank(R(S, λ 0 )) < n + q, as defined in [20]. According to [1], a system Σ is strongly detectable if and only if (iff) it has no transmission zeros λ 0 in the positive complex half plane C + , i.e., iff Re(λ 0 ) < 0 holds for all λ 0 , where R(S, λ 0 ) is rank deficient. A system Σ is strongly observable iff it has no transmission zeros at all.

B. Distance Measures for Strong Observability and Strong Detectability
The rank condition for strong observability and strong detectability has a binary character, and therefore, its reliable numerical evaluation is difficult. In particular, numerical errors while checking the conditions or small parametric perturbations of the system matrices could change the result. In order to circumvent these issues, the continuous measures and their calculation are introduced in [9, Theorem 2]. The measures are defined as the minimum norm parametric perturbation Δ that renders the system non-strongly observable or non-strongly detectable. This perturbation is calculated as the smallest singular value of R(S, λ) in the corresponding part of the complex plane and can be considered the distance from system S to a perturbed system with a transmission zero at λ.
In the distance measures calculation, the minimization problem (5) or (7) needs to be solved. The corresponding target function can have multiple local minima, because multiple systems with different transmission zeros λ but equal distance to Σ may exist. Therefore, f (λ) is a nonconvex function, which can include local minima that give only upper bounds for the global minimum. For the correct calculation of the distance measures, it is therefore necessary to have an optimization algorithm, which reliably finds the global minimum of f (λ).

III. THEORETICAL FOUNDATION FOR A GLOBAL OPTIMIZATION ALGORITHM
In order to reliably find the global minimum of f (λ) despite its nonconvexity, an algorithm is proposed to verify the existence of successively tighter, lower, or upper bounds in every iteration step. In this section, first a theorem for the existence of such bounds is given, then a proposition for their verification is derived, and finally a corollary of this proposition is presented for a lower dimension case.

A. Fundamental Theorem
The following theorem allows to check if a pair of values (q, q − δ) contains a lower or upper bound of the global minimum.
denote the minimum of f (λ) in D γ and the minimum of f (λ) along the border ∂D γ , respectively. Assume that holds.
The following implications are obtained from Theorem 1 and bounds for μ * can be derived, provided that q ≤ μ γ holds: If λ int exist, for which (10) holds, then yields an upper bound. If no such point exists, then the assumptions of Theorem 1 are contradicted, and yields a lower bound. Remark: Note that Theorem 1 is formulated such that bounds for either μ SO or μ SD can be deduced by means of a parameter γ. For the determination of μ SO , set γ = −∞ to ensure that ∂D γ includes C. For the determination of μ SD , set γ = 0 to restrict D γ to C + . Proof: A full proof of Theorem 1 is given in Appendix A. In the following, a short summary of the proof is presented. Fig. 1 shows the geometrical relations of the proof as a visual aid.
Due to the continuity of f (λ), it is possible to find a closed curve ∂D γ , where f (λ) = q holds with the point of minimum P * , at which f (λ * ) = μ * (S) holds, in its interior. Consider two points P 1 and P 2 located on ∂D γ with the same imaginary part as P * . There exists a lower bound for the distances η 1 and η 2 from the two points P 1 and P 2 to P * . The lower bound can be derived from the fact that the change of function value of f (λ) is bounded by f (λ * + η) − f (λ * ) ≤ |η| as can be shown by using Weyl's inequality [21]. By inserting the function values for f (λ * + η) and f (λ * ), the lower bounds are given by η 1 ≥ q − μ SO (S) and η 2 ≥ q − μ SO (S), respectively. When ∂G is shifted left and right by δ, new versions ∂G L and ∂G R are obtained, respectively. There exist points of intersection between ∂G L and ∂G R at P int = (α int , β int ), where (10) holds since δ ≤ q − μ SO (S) holds. The rightshifted version ∂G R is strictly in D γ , with Re(∂G R ) ≥ (γ + δ), as shown in Fig. 1, and therefore any point of intersection P int is in D γ with α int > (γ + δ).

B. Verification of Intersection
The challenge when using Theorem 1 for finding bounds on μ SO and μ SD is now to verify if a point of intersection P int does exist in order to show that (10) holds. For this verification the following proposition is introduced.
Remark: For the verification of the existence of P int and consequently λ int satisfying (10) in Theorem 1, the following steps are necessary. First, check if the matrix pair Π(δ), Γ has a purely real-valued eigenvalue α i , for which α i ≥ (γ + δ) holds. Second, for every α i , calculate the eigenvalues of M(α i + δ) and M(α i − δ). If they share a common purely imaginary eigenvalue jβ int for a α int = α i , then there exists a point P int = (α int , β int ) and λ int satisfying (10).

C. Minimization Along the Border
Applying Theorem 1 requires the knowledge of μ γ , the minimum value of f (λ) on the border ∂D γ located at P * γ = (α * γ + jβ * γ ). It is necessary to find an algorithm for minimizing f (λ) along ∂D γ , i.e., a line with λ = γ + jβ with β ∈ R. For the construction of such an algorithm, a corollary for an upper bound of μ γ can be deduced from Proposition 2 by inserting α int = γ and δ = 0. In case γ = −∞ is chosen, μ γ = σ min (F) holds, as given in Appendix A1. Otherwise, the following corollary can be used.
with eigenvectors w and v of M(γ) for the eigenvalues jβ and −jβ, respectively.
Since δ = 0 holds, (10) reduces to and can be used as a bound for a bisection algorithm.

IV. ALGORITHMS
In this section, a trisection algorithm for calculation of μ SO and μ SD is presented. The steps of the algorithm are described, the conditions for the bounds and their values are given, and the initial condition of the algorithm is discussed. In addition, a bisection algorithm for calculation of μ γ is given.

A. Trisection Algorithm for the Calculation of Strong Observability and Strong Detectability
Using the results presented in Section III, an algorithm for finding the global minimum μ * of f (λ), given by Algorithm 1, can be constructed by the following steps. First, a possible range μ * ∈ [L 0 , U 0 ] is initialized, such that it includes the global minimum. The definition of L 0 and U 0 is given in the next paragraph. Afterwards, in every iteration step i, the parameters are chosen, such that the resulting bounds split the range of the previous step in thirds, hence, the name trisection algorithm. By verifying if a point of intersection exists using Proposition 2, either (11) or (12) can be used for adapting L i or U i accordingly. The stopping criterion for the algorithm is given by U i − L i < . The range [L i , U i ] always includes μ * , and therefore it can be guaranteed to find μ * with a precision of .
The values L 0 and U 0 are chosen as L 0 = 0 and U 0 = μ γ , respectively, because 0 ≤ μ * ≤ μ γ always holds. Furthermore, selecting U 0 = μ γ ensures that the condition q ≤ μ γ in Theorem 1 holds in every iteration step of the algorithm.
The possible range for μ * is reduced by a factor of 2/3 in every iteration step, and can be calculated as U i − L i = (U 0 − L 0 )(2/3) i . Starting from this equation and inserting L 0 = 0, the necessary number of iteration steps for finding the minimum μ * with precision is given by Remark: A reduction of computation cost is possible by reformulating the generalized eigenvalue problem of size (2n) 2 × (2n) 2 in (51). Similar to the considerations in [19], (51) can be reduced to a 2n × 2n standard eigenvalue problem by splitting X up into submatrices, and further exploiting symmetry properties. This modification is given in Appendix B1.

B. Bisection Algorithm for Computation of Initial Bound
From Corollary 3, it is possible to construct an algorithm for finding the minimum μ γ given by Algorithm 2. First, a possible range μ γ ∈ [L 0 , U 0 ] is defined, with L 0 = 0 and U 0 = σ min (F) to ensure that μ γ is included, as given in Appendix A1. Afterwards, in every iteration step i a parameter is chosen, splitting the range of the previous step in half. For deriving the bounds for the bisection algorithm, the following set of equivalences can be drawn from Corollary 3: If M(γ) has a purely imaginary eigenvalue, then U i = q ≥ μ γ is a new upper bound for μ γ . Conversely, L i = q < μ γ is a new lower bound for μ γ , if M(γ) does not have a purely imaginary eigenvalue. The stopping criterion for the algorithm is given by U i − L i < . Therefore, it can be guaranteed to find μ γ with a precision of .
The possible range for μ γ is reduced by a factor of 1/2 in every iteration step. Following the same considerations as in the previous section, the necessary number of iteration steps for finding the minimum along the imaginary axis μ γ with precision is given by

V. NUMERICAL EXAMPLE
In order to show the practical usefulness of the presented algorithms, they are implemented in MATLAB and tested for a numerical example. For comparison with a state-of-the-art global optimization algorithm, the genetic algorithm ga, included in the MATLAB Global Optimization Toolbox, is used. In order to find the global minimum with acceptable confidence, an initial population size of 10 3 and a minimum improvement of the fitness function of 10 −10 are chosen for the ga algorithm.
As a practical application, the strong observability and strong detectability of the steel beam, shown in Fig. 2, are analyzed for different sensor configurations. One end of the beam is fixed, and on the free end an unknown force F is acting. At two locations z 1 and z 2 , the displacement, the acceleration, or both can be measured.
A linear modal model of the cantilevered beam, as presented for example in [22], is given by where Ω = diag(ω 0,1 , ω 0,2 , . . . , ω 0,k ) is a diagonal matrix containing the resonance frequencies ω 0,j , (j = 1, . . . , k), Γ =    , γ 2 , . . . , γ k ) is a diagonal matrix containing the corresponding modal damping ratios γ j , (j = 1, . . . , k), and k is the number of modes considered in the model. The matrices Φ w ∈ R k×1 and T ∈ R 2×k collect the mode shape vectors at the points z 1 and z 2 , respectively. The mode shapes describe how sensitive the system is for excitation and measurement at a specific point. The beam is modeled by the first k = 8 modes and the resulting system is of order n = 16. The modal parameters for the system (25) are determined using the finite element simulation program Ansys Academic Student 2019 R3 and are given in Table I.
As an illustration of the iterative calculation process for the beam with acceleration and position measurement at z 1 , Fig. 3 shows the mesh plot of f (λ) and the closed curves with f (λ) = q i . It can be observed how the closed curves first include multiple local minima but get closer and closer to the global minimum in every iteration step.
Using algorithms presented in Section IV with a precision of = 10 −10 , the distance measures for strong observability and strong detectability are calculated. The results are given in Table II. If only an acceleration sensor is used, the system is neither strongly observable nor strongly detectable. If an acceleration sensor and a displacement sensor are used at the same position z 1 , then the system is not strongly observable but strongly detectable. If acceleration sensors and displacement sensors are used at different locations, the system is strongly observable and strongly detectable. All these properties are correctly indicated by the corresponding values of μ SO and μ SD .
For the calculation of μ SO , 70 iteration steps are necessary in every sensor configuration, as can be determined by (22). The number of iteration steps for the calculation of μ SD ranges from 28 to 43 iteration steps because it depends on the start value given by the minimum along the imaginary axis. This minimum is first calculated by using the bisection algorithm in 42 iteration steps, as can be calculated using (24). In Table II, the ratio between the calculation times on off-the-shelf hardware (Notebook, Intel Core i5-8265 U, 1.80 GHz) of the trisection algorithm, with reduced eigenvalue problem t tri and of the genetic algorithm t gen , is given for different sensor configurations. It can be seen that the trisection algorithm finds the minimum roughly two to five times faster than the genetic algorithm.

VI. CONCLUSION
This article presented an algorithm for reliable and efficient calculation of measures for strong observability and strong detectability of linear and time-invariant systems. The challenge in calculating these measures was to solve a minimization problem with nonconvex target function. The main contribution of this article was an algorithm, which guarantees to find the global minimum of the target function for a given precision within a fixed number of iteration steps. The algorithm was constructed by exploiting the special structure of the target function. Using even more sophisticated numerical techniques, further reduction of the computational load may be possible. The computationally most expensive operation was finding the real eigenvalues of (51) or (59). In [19], a divide and conquer approach for efficient calculation of real eigenvalues was presented, which may be used in the presented algorithm as well.

A. Proof of Theorem 1
The singular value decomposition is continuous in the change of matrix entries, and therefore f (λ) is continuous in λ. Due to the fact that holds, it follows that there exists a λ 1 , such that f (λ 1 ) = q.
Since f (λ) is continuous and (26) holds, there exist closed curves that satisfy f (λ) = q. The point of minimum P * = (α * , β * ) in D γ , at which f (λ * ) = f (α * + jβ * ) = μ * holds, must lie within one of these closed curves. Let ∂G denotes the closed curve with the smallest area G with P * in its interior. From (26), it follows that all points on ∂G are in D γ , as shown in Fig. 1. On a line parallel to the real axis with P * on it, there exist two points P 1 = (α * − η 1 , β * ) and P 2 = (α * + η 2 , β * ) on ∂G with η 1 , η 2 ∈ R, η 1 > 0, and η 2 > 0. If the line crosses ∂G multiple times, then P 1 and P 2 are chosen as the points closest to P * . For the points P 1 and P 2 holds. The point P * lies within G and also on the line segment between P 1 and P 2 . Therefore, any point on this line segment lies within G. The change of singular values of a matrix T ∈ R m×n is bounded by the norm of the perturbation U ∈ R m×n as described by Weyl's inequality, see for example [21]. If with η ∈ R is inserted in (28) and the smallest singular value is considered, then |σ min (R(S, λ + η)) − σ min (R(S, λ))| ≤ ηI 0 0 0 (30) holds, and the change of function value is bounded by the norm of change of argument |η|. Combining (31) with the relation (27) for points P 1 and P 2 , it is possible to give lower bounds for their distances η 1 and η 2 to the point of minimum P * as If the curve ∂G is shifted left and right, i.e., in parallel to the imaginary axis, by a distance δ, two continuous closed curves are obtained. On the curve shifted to the left ∂G L , there exist two points P 1L = (α * − η 1 − δ, β * ) and P 2L = (α * + η 2 − δ, β * ). On the curve shifted to the right ∂G R , there exist two points P 1R = (α * − η 1 + δ, β * ) and P 2R = (α * + η 2 + δ, β * ). The point P 1L is left of P 1R since their real parts satisfy and therefore P 1L is in the exterior of ∂G R . Furthermore, P 2L is to the right of P 1R , because holds due to the fact that δ ≤ q − μ SO ≤ η 1 and, analogously, δ ≤ η 2 . From (34) and (35), it follows that two points P 1L and P 2L exist on ∂G L , one in the interior and one in the exterior of ∂G R , respectively. Since ∂G L and ∂G R are continuous closed curves, there exists at least one point of intersection P int = (α int , β int ) of the two curves ∂G L and ∂G R .
Since ∂G is strictly in D γ , also the right-shifted version ∂G R is strictly in D γ , with Re(∂G R ) ≥ (γ + δ), as shown in Fig. 1. Therefore, any point of intersection P int between ∂G L and ∂G R is in D γ with α int > (γ + δ).
From the definition of the shifted curves (33), it follows that λ int + δ and λ int − δ are on ∂G and solutions for (10). Therefore, (10) has at least one solution λ int in D γ with α int > (γ + δ).

B. Limit of f
The operations of finding the limit of the function and calculating the singular values can be interchanged, and it follows that lim |λ|→∞ f (λ) = σ min (Q ∞ (S, λ)) = min{c, σ min (F)} = σ min (F) holds for the limit of f (λ).

C. Proof of Proposition 2
From the definition of the singular value decomposition it follows that if f (λ) = σ min (R(S, λ)) = q holds, then there exists a vector z such that R(S, λ) H R(S, λ)z = q 2 z holds. It is possible to split (39) up into the following two equations: From the second line of (40a), w = Cz 1 + Fz 2 can be calculated and inserted into the second line of (40b), giving E T v+F T Cz 1 +F T Fz 2 = q 2 z 2 where G = (F T F − q 2 I) −1 . Since q satisfies q < σ min (F) this inverse exists. Inserting w and z 2 into the first lines of (40a) and (40b), one inserted into first and fourth lines, to get vec(X 11 ) vec(X 22 ) = 2α int vec(X 11 ) vec(X 22 ) .