Sensitivity matrices as keys to local structural system properties of large-scale nonlinear systems

Sensitivities are shown to play a key role in a highly efficient algorithm, presented in this paper, to establish three fundamental structural system properties: local structural identifiability, local observability, and local strong accessibility. Sensitivities have the advantageous property to be governed by linear dynamics, also if the system itself is nonlinear. By integrating their linear dynamics over a short time period, and by sampling the result, a sensitivity matrix is obtained. If this sensitivity matrix satisfies a rank condition, then the local structural system property under investigation holds. This rank condition will be referred to in this paper as the sensitivity rank condition (SERC). Applying a singular value decomposition (SVD) to the sensitivity matrix not only determines its rank but also pinpoints exactly the system components causing a possible failure to satisfy the local structural system property. The algorithm is highly efficient because integration of linear systems over short time-periods and computation of an SVD are computationally cheap. Therefore, it allows for the handling of large-scale systems in the order of seconds, as opposed to conventional algorithms that mostly rely on Lie series expansions and a corresponding Lie algebraic rank condition (LARC). The SERC and LARC algorithms are mathematically and computationally compared through a series of examples.


Introduction
In establishing local structural identifiability of nonlinear dynamical systems [1][2][3][4], we found that parametric output sensitivity matrices turned out to be the key to a simple and efficient algorithm that determines local structural identifiability, even for large-scale nonlinear systems. Moreover, the signatures obtained from this algorithm pinpoint exactly to the parameters involved in total correlations that, in fact, cause the lack of identifiability. As compared to alternative algorithms, based on, e.g., Lie brackets [5][6][7][8][9][10] or directed graphs [11][12][13], this algorithm puts hardly any restrictions on system size or system structure. Moreover, in terms of efficiency, the algorithm is highly competitive. Together, this opens up a large field of applications, e.g., in mechanical engineering, systems biology, and chemical engineering where use is made of nonlinear, possibly large-scale dynamic systems described by sets of ordinary differential equations.
The beneficial properties of the algorithm are largely due to the fact that parametric output sensitivities are governed by linear dynamics, irrespective of the possibly nonlinear dynamics of the underlying system. They also relate to the fact that the linear dynamics need only be integrated over a small timeinterval and the fact that, after sampling of the result, an SVD is the only additional computation required. Very recently, in an extensive overview and comparison of computational methods for identifiability [14], this algorithm was mentioned to be incredibly fast. In the earlier comparisons [15,16], the computational complexity of establishing identifiability had already been noted. This motivates us to extend this algorithm, developed in [1][2][3] for local structural identifiability, to local observability and local strong accessibility, as presented in this paper. All three properties are highly fundamental structural properties of nonlinear systems.
Computational methods for local identifiability of parameters, often also including initial conditions, have received significant attention in the last decade, as can be seen from the comparison and large number of references in [14]. For nonlinear systems, a Lie algebraic rank condition (LARC) still plays a central role in almost all of them and is also used to check local observability [9]. LARC is also central to algorithms detecting local strong accessibility, which is closely related to local controllability [7,17]. The computation of LARC requires repeated symbolic Lie differentiation or Lie Bracketing which are computationally expensive, especially with growing system dimensions. But even for small-scale systems, the computation may be very expensive [18]. In addition, no known stopping condition may exist as long as LARC is not satisfied [19]. Because local observability is almost dual to local strong accessibility [5,7,20], it is remarkable that a similar development of algorithms detecting local strong accessibility seems not to have taken place [17,21,22]. A possible reason for this has recently been discussed by us [18]. To extend the algorithm, we will further investigate and exploit duality in this paper.
The algorithm based on sensitivities culminates in computing what we refer to as the sensitivity rank condition (SERC). The shorthand SERC will also be used to refer to the algorithm itself. We also compare SERC with a conventional algorithm establishing the same three properties by means of a Lie algebraic rank condition (LARC) [5][6][7][8][9][10]. This conventional algorithm, also denoted by the shorthand LARC, relies on Lie series expansion, differential geometry and symbolic computation. From a practical perspective, SERC is shown to heavily outperform LARC in terms of computational efficiency and feasibility for largescale systems. These results are illustrated with a series of examples, including both small-and largescale systems.
The outline of this paper is as follows. Starting from the local structural identifiability algorithm, Sect. 2 presents the extensions to assess local observability and local strong accessibility revealing that all three properties can be treated within one single framework. Section 3 presents a series of examples concerning both small-and large-scale systems, as well as details concerning implementation of the algorithm. One large-scale example demonstrates how the amount by which SERC outperforms LARC in computational efficiency, heavily grows with increasing system dimensions. Section 3 concerns an extensive comparison and discussion of SERC versus LARC. For local strong accessibility and local observability, generic equivalence of SERC and LARC is proved in appendix 2. A similar proof seems possible for local structural identifiability. But the additional technical issues involved are considered to be beyond the scope of this paper, that focusses on application of SERC. Section 5 presents conclusions, an important one stating that SERC reveals how linear dynamics are fundamental in determining local structural properties of dynamic systems, also if the systems themselves are nonlinear.
2 Algorithms establishing local structural identifiability, local observability and local strong accessibility

Systems and trajectories
In this paper, we consider dynamical systems represented by ordinary differential equations. These can be written in the state-space format Here, variable t denotes continuous time, x t ð Þ is a vector containing the state variables, u t ð Þ is a vector containing the input variables of the system, h represents the vector of parameters and f is an analytic vector function also referred to as the system dynamics. Furthermore, y t ð Þ is a vector containing the output variables available for observation, identification and control and h is an analytic vector function that determines how the output variables depend on the state variables and parameters. Analytic vector functions f , h may be nonlinear and f need not be affine in u as is required in most conventional controllability, observability and structural identifiability analyses [5][6][7][8][9][10]. The algorithms to be presented will perform computations along an a-priori specified trajectory governed by Eqs.
This trajectory is fixed if we fully specify the initial conditions, the parameter values and input: All results obtained from our algorithms are therefore conditioned on a choice of (2.3). This choice will be further discussed in Sect. 3.

Local structural identifiability algorithm
Our starting point is the algorithm presented in [1] that addresses local structural identifiability. The so-called parametric state and output sensitivity matrix functions ox oh t ð Þ, oy oh t ð Þ associated with system (2.1), (2.2) play a key role. For notational convenience, these matrix functions will be written as x h t ð Þ and y h t ð Þ having dimensions n Â p and m Â p, respectively. At each time t, x h and y h indicate the sensitivity of state vector x t ð Þ and output vector y t ð Þ to the parameter vector h. Straightforward differentiation of Eqs. (2.1), (2.2) with respect to parameter vector h, and interchanging differentiation with respect to different, independent parameters, reveals that the dynamics of these sensitivity matrix functions is governed by Simultaneous integration of differential Eqs. (2.1), (2.4) from the initial conditions x t 0 ð Þ in (2.3) and, Sensitivity matrix function y h t ð Þ having p independent columns is a necessary and sufficient condition to uniquely determine h from y t ð Þ of system (2.1), (2.2), given a choice of (2.3) [23]. One way to check this rank condition is to evaluate y h t ð Þ; at discrete times t 0 \t 1 \::\t N , and to compute the rank of the resulting parametric output sensitivity matrix: To avoid obvious rank deficiency of Y h , the number of discrete time points N should be chosen such that, In principle, the sampling instants t 0 \t 1 \::\t N in (2.7) may be chosen arbitrarily close, thus requiring only a very small integration interval t 0 ; t N ½ : In practice, this interval should be selected large enough for all timescales of the system, especially the larger ones, to become manifest. To verify this, one may increase the length of the integration interval and check whether this increases the rank of sensitivity matrix Y h .
The rank condition, referred to as SERC, now reads as rank Y h ð Þ ¼ p. It is a necessary and sufficient condition to uniquely determine h from y t ð Þ of system (2.1), (2.2), given a choice of (2.3). To determine SERC, it is convenient to compute a singular value decomposition (SVD) of Y h . If this matrix is rank deficient, its SVD yields one or more zero singular values. Corresponding to each zero singular value, a singular vector is obtained from the SVD. The nonzero elements of such a singular vector indicate the parameters that are totally correlated and therefore cannot be identified uniquely. A graphical representation of these singular vectors is called a signature [1].

Local observability algorithm
By considering the initial state x t 0 ð Þ as the to be identified system parameters h, local structural identifiability turns into local observability of state x t 0 ð Þ, i.e., local structural identifiability of the initial state x t 0 ð Þ from measurements y t ð Þ, t ! t 0 . In this case, state vector x t 0 ð Þ acts as parameter vector h in Eqs. (2.1)-(2.5). For notational convenience, let x 0 denote x t 0 ð Þ. Because in Eqs. (2.4), (2.5), f and h do not depend on x 0 , these equations simplify to 9Þ where x x 0 t ð Þ, y x 0 t ð Þ represent the parametric sensitivity matrix functions of x t ð Þ and y t ð Þ with respect to the initial state x 0 . Simultaneous integration of differential Eqs. (2.1), (2.9) from the initial condition (2.3) and initial condition 11Þ In Eq. (2.11), I n denotes the n Â n identity matrix, being the sensitivity of x t 0 ð Þ ¼ x 0 with respect to itself. Substitution of x t ð Þ, x x 0 t ð Þ in (2.9), (2.10) then yields y t ð Þ, y x 0 t ð Þ. Sensitivity matrix function y x 0 t ð Þ having n independent columns is necessary and sufficient to uniquely obtain x 0 from y t ð Þ of system (2.1), (2.2), given a choice of (2.3) [23]. As in (2.7), to check this rank condition, we concatenate evaluations of sensitivity matrix function y x 0 t ð Þ at times t 0 \t 1 \::\t N to obtain the following sensitivity matrix: To avoid obvious rank deficiency of Y x 0 , the number of discrete time points N should now be chosen such that N þ 1 ð Þm ! n: ð2:13Þ Similar as for local structural identifiability, SERC for local observability now reads as rank Y x 0 ð Þ ¼ n. It is a necessary and sufficient condition to uniquely obtain x 0 from y t ð Þ of system (2.1), (2.2), given a choice of (2.3). If Y x 0 does not have full rank n, the nonzero elements of singular vectors corresponding to zero singular values indicate state variables that are totally correlated and therefore cannot be identified uniquely [1].

Linking observability of linear and nonlinear systems through sensitivities
Here, using sensitivities, we will show how observability of linear systems can be used to asses local observability of system (2.1), (2.2), that may be nonlinear. Recall that sensitivity matrix functions x x 0 t ð Þ; y x 0 t ð Þ have the advantageous property of being governed by the dynamics (2.9), (2.10) that are linear. This linearity is easily understood from the fact that the sensitivity of x t ð Þ and y t ð Þ with respect to x t 0 ð Þ ¼ x 0 concerns the propagation of infinitesimal deviations dx t 0 ð Þ;dy t 0 ð Þ from x t 0 ð Þ;y t 0 ð Þ along a trajectory of the nonlinear system (2.1), (2.2) specified by (2.3). The corresponding dynamic equations of dx t ð Þ, dy t ð Þ are therefore linear dynamic equations evaluated along the trajectory under consideration causing these linear dynamics to be time-varying. Associated with time-varying linear systems is the state transition matrix U t; t 0 ð Þ, which for system (2.14) is defined in the following way [24]: By substituting (2.15) in (2.14), one immediately sees that the state transition matrix U t; t 0 ð Þ also satisfies matrix differential equation (2.9) if we identify x x 0 t ð Þ with U t; t 0 ð Þ. Moreover, the initial condition (2.11) of sensitivity matrix function x x 0 t ð Þ is precisely the initial condition of state transition matrix U t; t 0 ð Þ. Therefore, sensitivity matrix function x x 0 t ð Þ is nothing but the state transition matrix U t; t 0 ð Þ of the linear dynamics (2.14). The linear system (2.14) is observable if matrix function oh ox t ð ÞU t; t 0 ð Þhas n independent columns [24]. Note that, in [24], this type of observability is called Aobservability [25]. From Eqs. (2.10) and (2.16), we find that Therefore, the condition for observability of the linear dynamics (2.14) is identical to the condition that sensitivity matrix function y x 0 t ð Þ; t 0 t t N has n independent columns, as required for local observability of x 0 from the trajectory. Thus, we arrive at the following important conclusion: For nonlinear system (2.1), (2.2), local observability of x 0 from a trajectory specified by (2.3) is equivalent with observability of the linear dynamics (2.14) along that trajectory.
Observability of linear systems is very well developed as part of linear control system design [24][25][26][27]. Specifically, observability is an important property for the successful application of linear output perturbation feedback controller designs, to control nonlinear systems along optimal state and control trajectories [27]. In this context, however, the linear dynamics represent an approximation of the nonlinear dynamics (2.1), (2.2) close to the trajectory, whereas in our context, they are an exact description of the propagation of sensitivities along the trajectory.

Local strong accessibility algorithm
By considering sensitivities, in section 1.3, we found the important equivalence between observability of the linear dynamics along trajectories and observability of the initial state x 0 from such trajectories. Since for linear systems, observability and controllability are dual properties, this suggests a similar equivalence between local controllability and controllability of linear dynamics along trajectories. Based on this idea, in [20,18], we already presented some preliminary results. In this section, the dual results are presented. We will use these in section 3 to argue that the equivalence holds if local controllability is replaced with local strong accessibility, which is a slightly weaker property [7]. Consider the linear system describing the propagation of infinitesimal perturbations dx,du along trajectories of system (2.1). To pursue duality, the following substitutions must be made in section 1.2 [20,18,[23][24][25]: As dual equivalents of Eqs. (2.9)-(2.11), we now obtain The dual equivalent of sensitivity matrix (2.12) is containing samples of u x N t ð Þ. The dual equivalent of (2.17) becomes From (2.25), u x N t ð Þ having n independent columns is the dual equivalent of U t N ; t N À t ð Þ of ou having n independent rows which is precisely the condition for the linear dynamics (2.18) to be what is called Acontrollable in [24]. To avoid obvious rank deficiency in (2.23), N should now be chosen such that, Þr ! n: ð2:26Þ

Implementation
The algorithms described in the previous sections rely on numerical integration of the dynamics (2.1), (2.2), (2.4), (2.5) for local structural identifiability, and (2.1), (2.2), (2.9), (2.10) for local observability along a trajectory determined by a choice of (2.3). In section 3, we will argue that, in general, selecting the input u t ð Þ, 3) arbitrarily, results in sufficient excitation of the linear sensitivity dynamics required for SERC to produce the same rank as LARC. The sampling instants at which the sensitivities are evaluated are taken to be equidistant here, i.e., To satisfy (2.8), we select N ¼ 2p, and to satisfy (2.13), (2.26) N ¼ 2n.
The algorithm for local structural identifiability has been used earlier in [1,3], on both small-and largescale systems to assess local structural identifiability of h, sometimes with the initial condition of the state variables x 0 included as additional parameters. As explained in section 1.2, identifiability of the initial state variables x 0 alone is equivalent with local observability of x 0 . Since in [1][2][3], identifiability of parameters and initial states has already been covered extensively, most examples in this paper concern local strong accessibility. These examples require numerical integration of Eqs. (2.1), (2.20), (2.21) backwards in time. We will find that in all examples, numerical integration stays tractable and is efficient, even for large-scale systems, since we can restrict integration to relatively short time intervals. The only requirement for the length of this time interval is, of course, that all timescales of the system become manifest. The partial derivatives in the sensitivity equations are computed using automatic differentiation [28], since pure symbolic differentiation may become intractable for largescale systems [29]. Alternatively, one may use complex-valued algebra for a very accurate determination of the Jacobi matrices that appear in the sensitivity equations [30]. We used MATLAB R2019A for our computations on a PC having a 3.10 gigahertz Intel Core i5-8600 processor and running on Windows 10. For automatic differentiation, we used a method obtained from MATLAB Central, supplied by Martin de la Gorce, that supports sparse matrices. The latter improves the algorithm efficiency, especially for large systems. For numerical integration, we used MATLAB function ode45 with default settings. We deliberately do not select an integrator for stiff systems, since these may approximate fast transient behavior corresponding to small time scales, that are important for detecting structural system properties.
For local structural identifiability and local observability, a singular value decomposition (SVD) of sensitivity matrix Y x 0 in (2.7), (2.12) plays a key role in the analysis. Similarly, for local strong accessibility, an SVD of sensitivity matrix U x N in (2.23) plays this key role. The SVD of this matrix can be represented by This decomposition is a sum of equally sized matrices u i v T i , weighted by their corresponding singular values r i ! 0. If one or more r i are (numerically) zero, U x N is rank deficient. Column vectors v i corresponding to zero singular values r i are basis vectors of the null space of U x N . Given that N þ 1 ð Þr ! n, these vectors determine directions in state space along which state x N cannot be controlled locally. In the examples presented in the next two sections, the column vectors v i will be graphically presented as to clearly show these directions. As in [1], we call these graphs the strong accessibility signature. Of course, if the system satisfies SERC, then nonzero singular values are present and the strong accessibility signature is empty.

Small-scale system examples
This artificial example is used in [31] to show how initial conditions x t 0 ð Þ can influence strong accessibility of systems which in turn may influence identifiability. The system is represented by ( The system (2.1), (3.3) is local strong accessible from all states except for states satisfying irrespective of the input u t ð Þ, which causes the system to be not local strong accessible from these states. This in turn causes a lack of local identifiability of system (2.1), (2.2), (3.3), (3.4) that does not occur otherwise [31]. We will use our algorithms to confirm these results. Taking arbitrarily, our algorithm for local strong accessibility produces the singular values r 1 ¼ 9:3806e À 01, r 2 ¼ 2:1674e À 02 of U x N in (2.23), confirming local strong accessibility. Next, we take (3.5) except for From (3.3), the choice (3.6) implies x 2 t ð Þ ¼ h 4 =h 3 , irrespective of the input u t ð Þ: Then, our algorithm produces the singular values r 1 ¼ 9:3522e À 01, r 2 ¼ 0. Clearly, the zero singular value indicates that local strong accessibility does not hold. The corresponding right singular vector in (3.2) is found to be m 2 ¼ 0; 1 ½ T . This confirms that only the second state is not local strong accessible.
Suppose we want to estimate all p ¼ 4 parameters of system (2.1), (2.2), (3.3), (3.4). Then, given (3.5), with x t N ð Þ replaced by x t 0 ð Þ, our identifiability algorithm produces singular values 1.7372e + 00, 7.2723e À 02, 1.6213e À 02, 5.6044e À 04, indicating local identifiability. Next, take (3.5) except for, Again, from (3.3), the choice (3.7) implies x 2 t ð Þ ¼ h 4 =h 3 , irrespective of the input u t ð Þ. Furthermore, observe that (3.7) represents a very special case in which the initial state x t 0 ð Þ depends on the parameters h. This requires a modification of the initial conditions (2.6) (initiating the sensitivity dynamics) of the identifiability algorithm Then, our algorithm produces the singular values 2.2712e + 00, 1.6340e À 01, 3.5852e À 02, 0 indicating that local identifiability is lost again. To the singular value r 4 ¼ 0 in (3.2), the right singular vector m 4 ¼ 0; 0; À1; 0 ½ T corresponds. It indicates that only the third parameter is unidentifiable. All these results completely match the ones found in [31]. It is important to note that, of these examples, the one with the largest computation time took only 0:04 seconds.
Example 2 (Example 6.4 from [8] where local strong accessibility is called local reachability).
This example concerns the classic car parking problem where the dynamics of the car are represented by The dynamics show that in this case, there is no drift vector field, but instead two control vector fields (see section 3 for this terminology) can be manipulated.
We arbitrarily selected u t ð Þ ¼ sin 2pt x t N ð Þ ¼ 1.0129; 0:9605 ; 0:8504; 0:5950 The singular values of sensitivity matrix U x N in (2.23) are then computed to be 3:4230, 3:0000, 0:1001, 0:01782. Since no numerically zero singular values are found, local strong accessibility is established, as in [8]. The computation time required by the algorithm was 0.049 s. Selecting  From [8], we know that system (2.1), (3.10) has two modes that are not locally strong accessible. We  Fig. 1 in descending order. Note the logarithmic scale and the large gap between the singular values clearly indicating that two singular values are numerically zero. They reveal the two modes that are not locally strong accessible as already found in [8]. The computation time required by the algorithm was 0.1 s.
The right panel of Fig. 1 shows what we call the strong accessibility signature showing all elements of the two column vectors v i in (3.2) corresponding to the two numerically zero singular values r 3 , r 4 . Element j of each singular vector corresponds to state variable j, j ¼ 1; 2; ::; n. Elements of the one-but-last singular vector m 3 are denoted by 'o' and those of the last singular vector m 4 by 'x'. From this strong accessibility signature on the right, we conclude that state variable 1 is locally strong accessible while the two modes that are not locally strong accessible involve state variables 2, 3 and 4. The computation time required by the algorithm is 0.343 s.

Large-scale system examples
To demonstrate the high efficiency of our algorithm, and its applicability to large-scale systems, the first two examples presented in this section are those with the highest state dimension appearing in [32][33][34]. For one of these examples, it is very easy to extend the state dimension, as we will do. Furthermore, we will extend the analysis performed in [32][33][34] by checking local strong accessibility for different combinations of control inputs.
In this example, we consider a platoon consisting of one robot that is the leader and several other robots following the leader. The leader is the only robot with full information of the terrain and mission. The robots that follow maintain the same heading as the leading robot but have full control over their individual inputs. In the planar state-space model of the platoon, the leader and each follower are represented by five state variables. State variables 1 and 3 represent planar x; y coordinates of the robot. State variables 2 and 4 and their derivatives model damping present in each robot. State variable 5 represents angular position with respect to a planar reference frame. The dynamics of the leader is given by The followers are numbered by j ¼ 1; 2; ::; M. The dynamics of the j th follower is governed by Since no singular value is numerically zero, this confirms local strong accessibility in each case. Numerical integration was performed starting at a terminal state x N , having uniformly distributed components between À1 and 1. All controls were taken to be constant during backwards integration from terminal time t N ¼ 1 to time t 0 ¼ 0. The components of the control vector were also uniformly distributed random numbers between À1 and 1: The corresponding computation times required by our algorithm are presented in Table 1. To show the high efficiency of the algorithm for large-scale systems, a comparison is made with the Mathematica package called ProPac [8]. Computation times to execute the function ControlDistribution within Mathematica 11.2 on the very same PC are also recorded in Table 1. This function uses Lie brackets to establish local strong accessibility (called local reachability in [8]). The factor by which our algorithm outperforms the one using Lie brackets is seen to dramatically increase with the state dimension.
The numerical efficiency of the algorithm allows extended analysis of large systems, e.g., we may easily remove some control inputs and check the effect on local strong accessibility. For the case M ¼ 4, n ¼ 25 ð Þafter removing control input u 4 , Fig. 3 shows the singular values of U x N and the corresponding strong accessibility signature. These reveal that in this case, one mode is not locally strong accessible. Its direction involves state variables 3, 8 and 10. Thus, the other state variables are still locally strong accessible. If we remove u 3 , instead of u 4 , four modes appear to be not locally strong accessible, as can be seen from Fig. 4. These modes involve state variables 6, 7, 8 and 9 that all belong to the first robot that follows. The other state variables and robots are therefore still locally strong accessible.
Finally, we consider the case M ¼ 4 with u 6 removed. Figure 5 shows the resulting singular values and strong accessibility signature. As expected from symmetry arguments in [32], the effect of removing u 6 is similar to removing u 4 . This can also be concluded from the structure of equation (3.12). By comparing Fig. 5 with Fig. 3, we observe that one mode is obtained that is not locally strong accessible in both cases. Moreover, when u 6 is removed, the strong accessibility signature is shifted five indices to the  Example 5 A combined heat and power system of a hotel building (BCHP) [33,34].
In this example, we consider the controllability of temperature in different rooms and parts of a Hotel building by adjusting flow rates in the system. The state-space representation of this system is given in appendix 1. The system has 14 state variables all representing different temperatures in the system. It has 12 control inputs that represent different flow rates in the system. In [33], it is reported that the system is locally strong accessible. Deactivating sets of control inputs and applying our algorithm revealed that local strong accessibility is still preserved when C 1 ¼ u 1 ; u 5 ; u 6 ; u 8 ; u 9 f g are set to zero. The same holds for C 2 ¼ u 1 ; u 6 ; u 8 ; u 9 ; u 12 f g . These results follow from Fig. 6 where C 0 indicates the case that no controls are excluded. Numerical integration was performed starting from terminal state x N having randomly generated components between 0 and 1. All controls were taken to be constant during integration from terminal time t N ¼ 1 to time t 0 ¼ 0. The components of the control vector were generated randomly between À1 and 1. To the components of parameter vector h, we assigned random values between 1 and 2. The algorithm execution time varied between 0.361 and 0.409 s. Figures 7 and 8 show the singular values and strong accessibility signatures for four other sets of excluded control inputs, namely C 3 ¼ u 1 ; u 5 ; u 6 ; u 8 ; u 9 ; u 12 f g , C 4 ¼ u 1 ; u 5 ; u 6 ; u 7 ; u 9 f g , C 5 ¼ u 1 ; u 7 ; u 9 ; u 11 ; u 12 f g , Since all of them are nonzero numerically this confirms local strong accessibility in each case and C 6 ¼ u 1 ; u 2 ; u 6 ; u 8 ; u 9 ; u 12 f g . Observe that for the excluded sets C 3 À C 5 , one mode occurs that is not locally strong accessible. For C 3 , this mode coincides with state variable 11, for C 4 with state variable 6, and for C 5 with state variable 12. All other state variables remain locally strong accessible. For C 6 , two modes occur that are not locally strong accessible, coinciding with state variables 1 and 2.
The next and final example is a challenging one since it contains parameters that are ill-identifiable, i.e., practically very hard to identify. In addition, the system is also ill-conditioned in some cases. The example concerns large-scale mass-spring-damper systems that contain nonlinearities and in some cases flutter. Flutter implies stiffness of the system, which reduces the efficiency of numerical integration. Still, with some modification, our algorithm still comes up  with correct solutions in the order of seconds up to a minute.
Example 6: Identifiability of parameters of a massspring-damper system.
The system consists of n M unit masses m i ¼ 1; i ¼ 1; 2; :; n M connected in series with nonlinear springs and dampers. At each end of the series, a nonlinear spring and linear damper connect to a fixed wall, see Fig. 9. The state variables of this system are x i ; i ¼ 1; 2; :; n M representing displacements of the masses from their equilibrium position and x n M þi , i ¼ 1; 2; :; n M representing their velocities, with positive values pointing to the right. In Fig. 9, s i ; r i ; d i ; i ¼ 1; 2; ::; n M þ 1 represent ordinary spring constants acting linearly, nonlinear spring constants, and damping coefficients acting linearly. The system dynamics are represented by (2.1) with Assuming the displacement of only one mass is measured, we study the identifiability of parameters s i ; r i ; d i ;i ¼ 1; 2; :; n M þ 1 for different values of n M .
We selected n M in between 1 and 15 and random parameter values s i ; d i in between 0.2 and 1 and r i ¼ s i =2, i ¼ 1; 2; :; n M þ 1. The initial conditions x i , i ¼ 1; 2; :; 2n M were also selected randomly between À1 and 1. Then, in each case, all 3 n M þ 1 ð Þparameters turned out to be identifiable regardless of which mass displacement is measured. On the other hand, parameters of springs and dampers close to the right wall turned out to be ill-identifiable from measurement of the position x 1 of mass m 1 on the left, especially for larger values of n M .
Here, we only present results of some interesting, difficult cases. For the large-scale example with n M ¼ 15, while measuring the displacement x 1 of the first mass m 1 ; we provide the singular values obtained from one trajectory in the left panel in Fig. 10 and of five trajectories in the right panel. In the latter case, the five sensitivity matrices (2.7) obtained from each trajectory, having identical parameters but different initial conditions, are stacked on top of one another resulting in a concatenated sensitivity matrix to which we apply the singular value decomposition (3.2). As explained in [3], [35], this concatenation can greatly improve the visibility of a possible gap in the singular values indicating a lack of identifiability. Since in this case all parameters are identifiable, the stacking only reduces the range of the singular values significantly. Note that, the ratio of singular values on the left is limited by the machine precision of 2:2 Â 10 À16 which may prevent locating a possible gap. This problem is thus resolved by concatenating the sensitivity matrices of five trajectories on the right. The execution time is approximately 3.3 s for each single trajectory.
As a final case, we present the results for n M ¼ 5 and n M ¼ 15 with flutter introduced by enlarging the spring constants s 1 ; r 1 ; s 16 ; r 16 with a factor 1000 and at the same time reducing d 1 ; d 16 with a factor 1000 turning the system into a stiff one. Moreover, perfect symmetric motions are realized on each side of the middle mass m 3 for n M ¼ 5, and m 8 for n M ¼ 15. In the latter case, this is achieved by taking s 8þi ¼ s 8Ài ; x 0;iþ7 ¼ Àx 0;iÀ7 , i ¼ 1; 2; ::; 7, x 0;i ¼ 0;i ¼ 8; 16; 17; 18; ::; 30;  Fig. 11 are singular trajectories, as further explained and discussed in the next section.

Is SERC a free lunch?
The benefits of SERC as compared to LARC have been clearly demonstrated in the previous section. They mostly relate to an enormous gain in efficiency, enabling the handling of large-scale systems, for which LARC tends to be very inefficient or to break down. Furthermore, LARC is restricted to nonlinear systems (2.1), (2.2) that are affine in the input. These are represented by (2.1), (2.2) with where f i x; h ð Þ, i ¼ 0; 1; ::; r are analytic vector functions. Note that, the property local strong accessibility, which is verified by both SERC and LARC, is a slightly weaker property than local controllability. If the so-called drift vector field f 0 x; h ð Þ ¼ 0, local strong accessibility becomes equal to local controllability.
The vector fields f i x; h ð Þ, i ¼ 1; 2; ::; r that multiply u i in (4.1), are called control vector fields [7]. SERC performs computations along the a-priori selected trajectory specified by parameters h, terminal state x t N ð Þ and input u t ð Þ; t 2 t 0 ; t N ½ in (2.3). Since system (2.1), (2.2) is assumed analytic, results concerning system structure obtained for terminal state x t N ð Þ on the trajectory, hold for all states x t ð Þ, t 2 t 0 ; t N ½ on the trajectory. Establishing the linear dynamics by which sensitivities propagate requires differentiation of the nonlinear system dynamics. This can be very accurately handled by automatic differentiation [28]. After equidistant sampling over time, a numerical sensitivity matrix is obtained of which the numerical rank must be determined. Apart from automatic differentiation, these numerical computations may be prone to numerical errors. As the examples in the previous section show, with default settings for numerical integration, and taking a still short, but long enough trajectory, leads generally to a clear gap in the singular values, if any, to properly determine the numerical rank. This gap can be increased, and the range of singular values decreased, by concatenating sensitivity matrices of different trajectories, as illustrated by Example 6. From this example, it is clear that when the range of singular values approaches the machine precision, our algorithm cannot locate a possible gap in the lower range. Since the singular values are actually measures of identifiability, when being close to the machine precision, they indicate extreme ill-identifiability. This is valuable information to engineers concerning practical identifiability but may prevent the algorithm from correctly establishing structural identifiability [36].
LARC on the other hand restricts analysis to a single state [7], denoted by x L in Fig. 12, not requiring specification of a trajectory. It verifies local strong accessibility for state x L , complying with local definitions of this property [7], [9]. LARC does not require specification of a trajectory since it locally considers all possible trajectories containing x L . This is realized by considering a Lie series expansion of the nonlinear  .2) around state x L , requiring Lie bracketing [5][6][7][8][9][10]. A disadvantage of LARC is that the number of Lie Brackets that are needed to reach the maximum rank is generally not known a-priori and may grow very large [19]. This causes the outcome to be uncertain, whenever the rank condition is not met and may cause the analysis to become infeasible [14][15][16], [18]. On the other hand, LARC determines structural identifiability and does not suffer from stiffness and other types of numerical ill-conditioning, such as ill-identifiability.
Only in exceptional cases, the difference between SERC and LARC may result in different outcomes. One example is the classic car parking problem (Example 2 in section 2.2, Example 1 in [18]) if u t ð Þ ¼ 0, t 2 t 0 ; t N ½ . Then, the trajectories are steadystate trajectories with corresponding sensitivity matrices that are rank deficient, instead of full rank. Although SERC correctly indicates a failure of controllability of the linear dynamics (2.14), other trajectories containing x L will not indicate this failure, see Fig. 12. This shows the importance of choosing the trajectory appropriately. In [37], singular trajectories are defined and recognized on which the equivalence between controllability of the linear dynamics (2.14) and local strong accessibility fails. This equivalence concerns precisely the one between SERC and LARC in this paper, when their rank is full. Moreover, in [37], trajectories are proved to be generically nonsingular. In Fig. 12, we have therefore adopted the same terminology. In appendix 2, using arguments from [37] and [5], the generic equivalence between SERC and LARC for local strong accessibility is also proved when LARC does not give full rank.
So far in this section, the discussion focused on local strong accessibility. Using similar arguments as for local strong accessibility, we also prove generic equivalence between LARC and SERC for local observability in appendix 2. Most likely, the proof for local structural identifiability can be obtained in a similar manner. But since this involves additional technicalities, it is considered out of the scope of this paper that focuses on application of SERC. Experimenting with SERC on different nonlinear systems and different trajectories, we found that if one selects the input u t ð Þ, t 2 t 0 ; t N ½ in (2.3) arbitrarily, SERC equaled LARC in all cases. This is in line with all examples in section 2 and [1][2][3]. In the case of local structural identifiability, the selection of a trajectory through an arbitrary choice of the input u t ð Þ, t 2 t 0 ; t N ½ in (2.3) may be identified with selecting an input that is ''sufficiently exciting'' for the output sensitivities, meaning that the input excites all relevant dynamics of the linear sensitivity system (2.4), (2.5) [38]. If one is in doubt about the nonsingularity of the trajectory containing x L in Fig. 12, one may of course test with several trajectories containing x L in Fig. 12 and select the one for which SERC gives the highest rank.

Conclusions
Sensitivity matrices and an associated rank condition (SERC) were shown to be keys to a unifying approach establishing three fundamental properties that characterize nonlinear system structure: local structural identifiability, local observability and local strong accessibility. This approach resulted in a highly Fig. 12 LARC evaluated at x L , versus SERC evaluated along nonsingular state trajectories x ns;i t ð Þ, i ¼ 1; 2; 3 and the singular trajectory x s t ð Þ, all containing x L but having different inputs u t ð Þ, t 2 t 0 ; t N ½ . LARC equals SERC along nonsingular trajectories. Along singular trajectories the rank obtained from SERC is less than the one from LARC. In practice, the number of nonsingular trajectories greatly outnumber the number of singular trajectories efficient algorithm enabling the structural analysis of systems with a state dimension up to 250 within several seconds. Our algorithm opens the possibility to verify these properties for large-scale systems where other methods tend to fail or break down. The importance of sensitivity matrices has first been recognized in the context of system identification. In this paper, it was shown how it carries over to local observability and local strong accessibility.
The behavior of sensitivity matrix functions appears to be governed by time-varying linear dynamics. Integrating the linear sensitivity dynamics in parallel with the nonlinear system dynamics over short periods of time provides the sensitivity matrix functions. After appropriate sampling, a sensitivity matrix is obtained. A singular value decomposition applied to this matrix together with a sensitivity rank condition (SERC) yields enough information to decide whether the structural property holds. Moreover, the null space corresponding to zero singular values, also obtained from the singular value decomposition, reveals locally the cause of lack of identifiability, observability, or strong accessibility. It pinpoints the exact state variables and/or parameters that are involved in a total correlation. This null space can be graphically represented and is termed the identifiability/observability/strong accessibility signature.
For the examples presented in this paper and several references, the outcomes obtained along different, arbitrarily selected trajectories of the same system, are generally identical. This can be expected since local strong accessibility, local observability and local structural identifiability are properties of nonlinear systems that, in general, have a global character. Example 1 in this paper illustrated how, for a tiny set of initial conditions, structural properties may change. Last, but not least, we emphasize that the high efficiency of the algorithm based on sensitivities stems from their linear dynamics along trajectories. This holds even if the system itself is nonlinear.

Declarations
Conflict of interest The authors have no relevant financial or non-financial interests to disclose. Data will be made available on reasonable request.
Appendix 1: Dynamic model of the combined heat and power system of a hotel building [33,34] used in Example 5 f x; u; h ð Þ¼ k 9 x 10 À x 9 ð Þþ x 5 V 9 À x 9 V 9 u 9 k 10 x 9 À x 10 ð Þþ b 3 V 10 À x 10 V 10 u 10 output y e t