Derivative-Free Observability Analysis for Sensor Placement Optimization of Bioinspired Flexible Flapping Wing System

Abstract

Observability analysis of a bioinspired flexible flapping wing system provides a measure of how well the states of flexible flapping wing micro-aerial vehicles can be estimated from real-time measurements during high-speed flight. However, the traditional observability analysis approaches have trouble in terms of lack of quantitative analysis index, high computational complexity, low accuracy, and unavailability in stochastic systems with memory, including bioinspired flexible flapping wing systems. Therefore, a novel derivative-free observability analysis method is proposed here based on the generalized polynomial chaos expansion. By formulating a surrogate model to represent the relationship between the cumulative measurement and the random initial state, the observability coefficient matrix is calculated and the observability rank condition is stated. Consequently, several observability indices are proposed to quantity the observability of the system. Altogether, the proposed method avoids the disadvantages of the traditional approaches, especially in assessing the observability degree of each state and the effect of stochastic noise on observability. The validation of the proposed method is first provided by demonstrating the equivalence between the traditional and proposed methods and subsequently by comparing the observability of the Lorenz system calculated via three different approaches. Finally, the proposed method is applied on a bioinspired flexible wing system to optimize the placement of sensors, which is consistent with the natural configuration of campaniform sensilla on the wing of the hawkmoth.

1. Introduction

Due to the high maneuverability and portable structure of birds and insects, bioinspired flexible flapping wing micro-aerial vehicles (FWMAVs) have attracted considerable attention in the last decade [1,2,3]. A prerequisite to designing such agile FWMAVs is the simultaneous and accurate estimation of vehicle states, such as position and rotation rates, through the noisy information measured by the sensors on the wing. To solve this state estimation problem, the observability analysis of a bioinspired flexible flapping wing system is necessary in order to evaluate whether the system states can be inferred from the measurements. Thereby, the guidelines for measurement selection, sensor placement, and model optimization can be provided.

For a bioinspired flexible flapping wing system, deformation inevitably happens because of aerodynamic loads and inertial forces, and negatively influences the estimation of rotation rates based on the measurements from sensors on the wing. In nature, insects have sensing mechanisms that enable them to understand their flying states. For example, the strain information of a flying insect is measured and saved by campaniform sensilla on the wing, followed by conversion to action potentials (voltage spikes) in the sensory neurons [4,5]. After the delivery of these neural-encoded spikes to the central nervous system, the rotation rates are detectable. Such sensing and signal processing mechanics ensure insects’ fast and accurate sensory feedback and subsequent high agility [6]. Inspired by these biological mechanisms, the neural encoding processing of strain data can be constructed as neural encoding measurement model [7]. Unlike the raw strain measurement model used previously [8], the neural encoding measurement model has strong nonlinearity, the same as the flexible wing flapping dynamics model. In addition, the measurement at each time contains the current and previous state information. Altogether, the strong nonlinearity and memory function of the flexible flapping wing system make its observability analysis a significant challenge.

For a linear time-invariant system, the traditional observability rank condition is the simplest and most effective criterion to analyze whether the system is observable or not [9]. However, it fails in nonlinear systems, and even in linear time-variant systems. Hence, several observability analysis approaches for nonlinear systems have been developed. Rouhani et al. [10] applied the Lie-derivative-based observability analysis method to power systems with the realistic synchronous generator model, and Seo et al. [11] used the same method to analyze the observability of a missile interception system with line-of-sight angle measurement. This method is feasible only for simple and small-scale nonlinear systems, as it requires considerable derivation calculations. The empirical observability Gramian method, which is developed from the observability Gramian method through interval approximation, is more prevalent in practical applications because it does not depend explicitly on the system models and reduces the computational complexity [12]. For instance, Hinson et al. [8] utilized this method to determine the optimal sensor placement on the wing of the hawkmoth. Hodzic et al. [13] and Hinson et al. [14] respectively applied the empirical observability Gramian method to the sensor placement problem for generic wide-body aircraft and the path planning problem for under-sensed vehicles. However, the existing approaches are either derivative-based, which leads to high computational complexity, or approximation-based, which may reduce the analytical accuracy. Furthermore, these approaches are unable to measure the influence of stochastic noise on observability.

Recently, a derivative-free approach based on the generalized Polynomial Chaos (gPC) expansion has been proposed by Zheng et al. [15] to address the aforementioned problems. They further applied the derivative-free approach to a power system using the stochastic dynamic model [16]. While this method has been proven feasible and computationally efficient for nonlinear systems without control input and memory, it is impracticable for systems with memory, such as a flexible flapping wing system.

In this paper, we propose a modified gPC-based observability analysis approach for both linear and nonlinear systems incorporating memory and applying it to the analysis of a flexible flapping wing system. The observability indices obtained by the proposed approach are used to determine the optimal placement of sensors configured on the wing. The main contributions are as follows.

(1) The traditional observability analysis approaches are extended to linear and nonlinear deterministic systems with memory. By extending the definition of observability to deterministic systems with memory, a new relationship between the cumulative measurement and the initial state is formulated. Then, according to the implicit function theorem, the traditional observability rank condition for systems with memory is provided.

(2) A modified derivative-free approach based on the generalized Polynomial Chaos expansion for linear and nonlinear stochastic systems with memory is proposed, significantly reducing the computational cost and resolving the stochasticity of the system. Because the definition of observability is extended to stochastic systems with memory, it is necessary to represent this new relationship through a surrogate model. Thus, a different mapping between the stochastic input variable and standard normally distributed variable is formulated, resulting in a more appropriate surrogate model. Based on this model, the observability rank condition can subsequently be provided and its equivalence with the traditional approaches proven mathematically.

(3) Based on [15], several different observability indices, such as the first contribution rate and interference binary value, are proposed to describe the observability degree of systems, including the observability of each system state and the effect of measurement noise on the system observability.

(4) Finally, the proposed method is applied to observability analysis of a flexible flapping wing model inspired by the hawkmoth. Based on the proposed observability indices, the optimal sensor placement on the bioinspired wing is discussed and compared with the natural configuration of campaniform sensilla on the wing of the hawkmoth.

The remainder of the work is organized as follows. Section 2 formulates the flexible wing flapping dynamics model and neural encoding measurement model. Section 3 extends the traditional observability analysis approaches to linear and nonlinear systems with memory and discusses the shortcomings of the traditional approaches. Section 4 presents a derivative-free observability analysis approach and the observability indices for the linear and nonlinear system with memory, followed by the equivalence between the traditional and proposed approaches. Section 5 applies the proposed approach to the Lorenz system with memory and the bioinspired flexible flapping wing system. Section 6 discusses the optimal sensor placement based on the observability analysis. Finally, our conclusions are provided in Section 7.

2. System Model

This paper mainly focuses on the observability analysis problem of the bioinspired flexible flapping wing model based on neural encoded strain measurements. In this section, inspired by the North American hawkmoth, Manduca sexta, a low-order model of flexible wing flapping dynamics is first derived, followed by a neural encoding model of strain measurement.

2.1. Flexible Wing Flapping Dynamics

Here, the flexible flapping wing is considered as a thin and flexible cantilevered plate $\mathcal{A}$ in a rotating coordinate system, as depicted in Figure 1. This coordinate system is called the wing body frame, where $x_p$ are the positive axis points from the leading edge to the trailing edge along the wing root, the positive axis $y_p$ is perpendicular to $x_p$ and points from root to tip in the wing plane, and the positive axis $z_p$ is perpendicular to the wing plane. The cantilevered plate is assumed to deform out-of-plane without in-plane stretching or extension.

Hence, the out-of-plane plate deformation $w(x,y,t)$ is described by a finite number of orthonormal spatial modes, that is,

$$w(x,y,t)=\sum _{i=0}^{n_m}{\phi }_i(x,y){\eta }_i(t)$$

(1)

where ${\phi }_i(x,y)$ are the free-vibration mode shapes, ${\eta }_i(t)$ are the modal coordinates, and $n_m$ is the number of chosen modes used to represent the spatial deformation.

Here, two main wing mode shapes, i.e., the first bending mode and the first torsion mode, are chosen to capture the principle wing deformations [8], as illustrated in Figure 2.

We define the unit vectors in each axis direction as ${\mathit{i}}_p$, ${\mathit{j}}_p$, and ${\mathit{k}}_p$, the position of a mass element $d\mathcal{A}$ on the plate in the wing body frame is

$$\mathit{r}(x,y,t)=x{\mathit{i}}_p+y{\mathit{j}}_p+w(x,y,t){\mathit{k}}_p$$

(2)

and the velocity of $d\mathcal{A}$ is provided by

$$\begin{array}{cc}\hfill \mathit{v}& ={\mathit{v}}_0+\frac{\partial }{\partial t}\mathit{r}+{\omega }_0\times \mathit{r}\hfill \\ \hfill & =(U+Qw-Ry){\mathit{i}}_p+(V+Rx-Pw){\mathit{j}}_p+(W+\dot{w}+Py-Qx){\mathit{k}}_p\hfill \end{array}$$

(3)

where ${\mathit{v}}_0=U{\mathit{i}}_p+V{\mathit{j}}_p+W{\mathit{k}}_p$ is the velocity of the origin of wing body frame and ${\omega }_0=P{\mathit{i}}_p+Q{\mathit{j}}_p+R{\mathit{k}}_p$ is the angular velocity of wing body frame.

We define the mass density of the wing as $\rho (x,y)$, the thickness as $h(x,y)$, the kinetic energy as $T_e$, and the potential energy $U_e$ of the wing $\mathcal{A}$ as follows:

$$\begin{array}{cc}\hfill T_e(t)& =\frac{1}{2}{\int }_{\mathcal{A}}\rho (x,y)h(x,y)\mathit{v}(x,y,t)·\mathit{v}(x,y,t)d\mathcal{A}\hfill \\ \hfill U_e(t)& =\frac{1}{2}{\int }_{\mathcal{A}}\frac{h{(x,y)}^3}{12}\mathsf{\Lambda }{(x,y,t)}^T\mathit{M}\mathsf{\Lambda }(x,y,t)d\mathcal{A}\hfill \end{array}$$

(4)

where $\mathit{M}$ is the matrix of the material constants and $\mathsf{\Lambda }$ is the strain vector:

$$\mathsf{\Lambda }={\left[\begin{array}{ccc}\frac{{\partial }^{2}w}{\partial x^2}& \frac{{\partial }^{2}w}{\partial y^2}& \frac{2{\partial }^{2}w}{\partial x\partial y}\end{array}\right]}^T$$

(5)

Because the generalized coordinates describing the structural configuration are modal coordinates $\eta$, Lagrange’s equation is formulated as

$$\frac{\mathrm{d}}{\mathrm{d}t}(\frac{\partial T_e}{\partial {\dot{\eta }}_i})-\frac{\partial T_e}{\partial {\eta }_i}+\frac{\partial U_e}{\partial {\eta }_i}=Q_i,i=1,\cdots ,n_m$$

(6)

where $Q_i$ are the exogenous non-conservative generalized forces.

Substituting Equations (1), (3), and (4) into Lagrange’s equation, the equations of motion for the flexible wing in the wing body frame is

$$\ddot{\mathit{\eta }}+\mathbf{\Omega }\mathit{\eta }+{\mathit{M}}_a(\left[\begin{array}{c}\dot{W}\\ \dot{P}\\ \dot{Q}\end{array}\right]-{\mathit{C}}_a)=\mathit{Q}$$

(7)

where the modal coordinate vector is $\eta ={\left[{\eta }_1{\eta }_2\right]}^T$, the total aerodynamic force is $\mathit{Q}={\left[Q_1{Q}_2\right]}^T$, and the stiff matrix $\mathbf{\Omega }$, applied acceleration mass matrix ${\mathit{M}}_a$, and Coriolis force ${\mathit{C}}_a$ are provided by Equations (8), (9), and (10), respectively:

$$\mathbf{\Omega }=\left[\begin{array}{cc}{\omega }_1^2-P^2-Q^2& 0\\ 0& {\omega }_2^2-P^2-Q^2\end{array}\right]$$

(8)

where ${\omega }_i=2\pi f_i$ and $f_i$ denote the corresponding mode frequency of the ith mode shape.

$${\mathit{M}}_a={\int }_{\mathcal{A}}\left[\begin{array}{ccc}\rho h{\phi }_1& \rho h{\phi }_{1}y& -\rho h{\phi }_{1}x\\ \rho h{\phi }_2& \rho h{\phi }_{2}y& -\rho h{\phi }_{2}x\end{array}\right]d\mathcal{A}$$

(9)

$${\mathit{C}}_a=\left[\begin{array}{c}QU-PV\\ -QR\\ PR\end{array}\right]$$

(10)

Let the position of the feathering rotation axis relative to the origin of the wing body frame be ${\mathit{r}}_r=x_r{\mathit{i}}_p+0{\mathit{j}}_p+0{\mathit{k}}_p$; then, the velocity and acceleration of wing body frame are computed by

$$\begin{array}{cc}\hfill {\mathit{v}}_o& =-{\mathit{\omega }}_0\times {\mathit{r}}_r\hfill \\ \hfill {\dot{\mathit{v}}}_o& =-{\dot{\mathit{\omega }}}_0\times {\mathit{r}}_r-{\mathit{\omega }}_0\times ({\mathit{\omega }}_0\times {\mathit{r}}_r)\hfill \end{array}$$

(11)

Substituting Equation (11) into (7), the motion equation can be rewritten as

$$\ddot{\mathit{\eta }}+\mathbf{\Omega }\mathit{\eta }+{\mathit{M}}_a\left[\begin{array}{c}{x}_r(\dot{Q}-2PR)\\ \dot{P}+QR\\ \dot{Q}-PR\end{array}\right]=\mathit{Q}$$

(12)

The aerodynamics force $\mathit{Q}$ is modeled as the aerodynamics model provided in [8], that is,

$$\mathit{Q}={\mathit{Q}}_{t0}+{\mathit{Q}}_{a0}+Q_{t\eta }\eta +Q_{a\dot{\eta }}\dot{\eta }+Q_{a\ddot{\eta }}\ddot{\eta }$$

(13)

where ${\mathit{Q}}_{t0}$ and $Q_{t\eta }$ denote the translational generalized force vector and translational force stiffness matrix, respectively, and ${\mathit{Q}}_{a0}$, $Q_{a\dot{\eta }}$, and $Q_{a\ddot{\eta }}$ represent the added mass force vector, added mass damping matrix, and added mass acceleration matrix, respectively. The detailed equations and corresponding derivations are omitted because the aerodynamics model is not the focus of this work.

Let the state vector be $\mathit{x}={\left[\mathit{\eta }\dot{\mathit{\eta }}PQR\right]}^T$ and the control input vector be $\mathit{u}={\left[\dot{P}\dot{Q}\dot{R}\right]}^T$; then, the flexible wing flapping dynamics model is formulated by substituting Equation (13) into (12):

$$\dot{\mathit{x}}=\left[\begin{array}{c}\dot{\mathit{\eta }}\\ {(I-Q_{a\ddot{\eta }})}^{-1}[Q_{a\dot{\eta }}\dot{\mathit{\eta }}-(\mathbf{\Omega }-Q_{t\eta })\mathit{\eta }+{\mathit{Q}}_{t0}+{\mathit{Q}}_{a0}+{\mathit{M}}_{a_1}{x}_r(2PR-\dot{Q})+{\mathit{M}}_{a_2}(-QR-\dot{P})+{\mathit{M}}_{a_3}(PR-\dot{Q})]\\ \dot{P}\\ \dot{Q}\\ \dot{R}\end{array}\right]$$

(14)

Define the feathering angle as $\alpha$, the elevation angle as $\theta$, and the position angle as $\zeta$, as shown in Figure 3; these three Euler angles represent the motion of the wing body frame with respect to the inertial frame.

Assuming that the 3-1-2 rotation sequence is used, ${\mathit{\omega }}_0$ and the Euler angles have the following relationship:

$${\mathit{\omega }}_0=\left[\begin{array}{ccc}cos\alpha & 0& -sin\alpha cos\theta \\ 0& 1& sin\theta \\ sin\alpha & 0& cos\alpha cos\theta \end{array}\right]\left[\begin{array}{c}\dot{\theta }\\ \dot{\alpha }\\ \dot{\zeta }\end{array}\right]$$

(15)

where $\dot{\theta }$, $\dot{\alpha }$, and $\dot{\zeta }$ denote the feathering, elevation, and position angular velocity, respectively. The control input $\mathit{u}$ can be computed by taking the derivative of Equation (15) with respect to time t.

2.2. Neural Encoding Measurement Model

Assuming that a strain sensor on the wing plane is located at ${\mathit{r}}_s=x_s\overrightarrow{{\mathit{i}}_p}+y_s\overrightarrow{{\mathit{j}}_p}+z_s\overrightarrow{{\mathit{k}}_p}$, the strain measurement $ϵ(x_s,y_s,t)={\left[{ϵ}_{xx}(x_s,y_s,t){ϵ}_{yy}(x_s,y_s,t){ϵ}_{zz}(x_s,y_s,t)\right]}^T$ is provided by

$$ϵ(x_s,y_s,t)=-\sum _{i=1}^{n_m}{z}_s\left[\begin{array}{c}\frac{{\partial }^2{\phi }_i(x_s,y_s)}{\partial X^2}\\ \frac{{\partial }^2{\phi }_i(x_s,y_s)}{\partial y^2}\\ \frac{2{\partial }^2{\phi }_i(x_s,y_s)}{\partial x\partial y}\end{array}\right]{\eta }_i(t)$$

(16)

where ${ϵ}_{xx}$, ${ϵ}_{yy}$, and ${ϵ}_{xy}$ mean the chord-wise bending strain, span-wise bending strain, and shear strain, respectively. Because the assumed structural modes do not have chord-wise bending, ${ϵ}_{xx}$ is effectively zero, and only ${ϵ}_{yy}$ and ${ϵ}_{xy}$ are assumed to be measured by the sensors on the wing.

To determine the rotation rates, flying insects such as the hawkmoth Manduca sexta deliver the strain data from the campaniform sensilla on the wing to the central nervous system via sensory neurons, where the strain data are converted into action potentials (voltage spikes). This transformation process is called neural encoding. Inspired by this neural encoding mechanism, we constructed a neural encoding measurement model instead of using unprocessed strain measurements (16).

Several experiments have verified that the neural encoding model of hawkmoth wings is based on the two main functions [6,17], namely, the spike-triggered average (STA) and a nonlinear activation (NLA) functions.

The STA function is a measure that relates continuous strain stimulus to the spike, which represents the average stimuli taken during spike occurrences. This is approximated as an exponentially decaying sinusoidal function with a delay a

$$STA(t)=cos(2\pi f_{STA}(-t+a))exp(\frac{-{(-t+a)}^2}{b^2})$$

(17)

where b is the width and $f_{STA}$ is the STA frequency.

The firing rate of a neuron $\varkappa (x_s,y_s,t)$ can be estimated by convolution of the strain stimulus and the STA:

$$\varkappa (x_s,y_s,t)=\frac{1}{C_{\varkappa }}{\int }_0^{t_M}ϵ(x_s,y_s,t-\tau )STA(\tau )d\tau$$

(18)

where $C_{\varkappa }$ is a normalization constant and $t_M$ is the memory length.

Furthermore, a saturation function is added to properly reflect the neuron’s nonlinear behaviour. Hence, by substituting $\varkappa (x_s,y_s,t)$ into the STA function, the probability of firing $P_{fire}(x_s,y_s,t)$, which implies the probability of firing an action potential at the coordinate $(x_s,y_s)$ on the wing, is provided by

$$\begin{array}{cc}\hfill P_{fire}(x_s,y_s,t)& =NLA(\varkappa (x_s,y_s,t))\hfill \\ \hfill & =\frac{1}{1+exp(-c(\varkappa (x_s,y_s,t)-d))}\hfill \end{array}$$

(19)

where c is the slope, d represents the half-maximum position of the NLA function, $P_{fire}(x_s,y_s,t)$ is regarded as the neural encoding strain measurement, and the neural encoding measurement model is provided by Equation (19).

3. Traditional Observability Analysis for Deterministic System with Memory

In this section, the traditional deterministic observability analysis methods are extended to both linear and nonlinear systems with memory; the limitations of the existing methods are then discussed.

3.1. Linear System with Memory

Consider the following general linear discrete-time time-invariant system with memory:

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+1}& =\mathit{A}{\mathit{x}}_k+\mathit{B}{\mathit{u}}_k\hfill \\ \hfill {\mathit{y}}_k& =\sum _{\tau =0}^N{\mathit{C}}_{\tau }{\mathit{x}}_{k-\tau }\hfill \end{array}$$

(20)

where ${\mathit{x}}_k\in {\mathbb{R}}^{n\times 1}$, ${\mathit{u}}_k\in {\mathbb{R}}^{n_u\times 1}$, and ${\mathit{y}}_k\in {\mathbb{R}}^{m\times 1}$ are the state, control input, and measurement vector at time k, respectively. Here, N is the memory length, which implies that ${\mathit{y}}_k$ explicitly hinges on the information of states from time $(k-N)$ to time k and that the measurement information can be obtained if and only if $k\ge N$.

Hence, the definition of observability can be extended to the discrete-time dynamics system with memory (20) as follows:

Definition 1.

The system is (locally) observable over the interval $[k_1,k_2]$ if the initial state ${\mathit{x}}_{k_1-N}$ can be uniquely determined from ${\mathit{y}}_k$, $k\in [k_1,k_2]$.

The observability rank condition is determined by the following theorem:

Theorem 1.

System (20) is (locally) observable if and only if the observability matrix

$$\mathit{O}=\left[\begin{array}{c}\overline{\mathit{C}}\\ \overline{\mathit{C}}\mathit{A}\\ ⋮\\ \overline{\mathit{C}}{\mathit{A}}^{n-1}\end{array}\right]$$

(21)

is a full column matrix where $\overline{\mathit{C}}$ is defined as

$$\overline{\mathit{C}}=\sum _{\tau =0}^N{\mathit{C}}_{\tau }{\mathit{A}}^{N-\tau }$$

(22)

Proof of Theorem 1.

Define the cumulative measurement vector ${\mathcal{Y}}_k\in {\mathbb{R}}^{mn\times 1}$ as

$${\mathcal{Y}}_k={\left[\begin{array}{cccc}{\mathit{y}}_k& {\mathit{y}}_{k+1}& \cdots & {\mathit{y}}_{k+n-1}\end{array}\right]}^T$$

(23)

Substituting Equation (20) into Equation (23), the relationship between the arbitrary initial state ${\mathit{x}}_{k-N}$ and its corresponding measurements ${\mathcal{Y}}_k$ is provided as

$${\mathcal{Y}}_k=\left[\begin{array}{c}\overline{\mathit{C}}\\ \overline{\mathit{C}}\mathit{A}\\ ⋮\\ \overline{\mathit{C}}{\mathit{A}}^{n-1}\end{array}\right]{\mathit{x}}_{k-N}+\mathit{D}\left[\begin{array}{c}{\mathit{u}}_{k-1}\\ {\mathit{u}}_{k-2}\\ ⋮\\ {\mathit{u}}_{k-N}\end{array}\right]$$

(24)

where the $(i,j)$ element of the feedthrough matrix $\mathit{D}\in {\mathbb{R}}^{N\times (n-1)}$ is

$$\mathit{D}(i,j)=\sum _{\tau =1}^i{\mathit{C}}_{\tau -1}{\mathit{A}}^{N+j-1-\tau }\mathit{B}$$

(25)

According to the implicit function theorem [18], the initial state ${\mathit{x}}_{k-N}$ can be uniquely determined from measurements ${\mathcal{Y}}_k$ if and only if the Jacobian matrix

$$\mathit{O}=\frac{\partial {\mathcal{Y}}_k}{\partial {\mathit{x}}_{k-N}}$$

(26)

is a full column matrix.    □

The Gramian observability method provides an equivalent observability condition: the system (20) is (locally) observable if and only if the observability Gramian matrix

$${\mathit{W}}_{o_k}=\sum _{\tau =0}^k{({\mathit{A}}^T)}^{\tau }{\overline{\mathit{C}}}^T\overline{\mathit{C}}{\mathit{A}}^{\tau }$$

(27)

is nonsingular.

These observability rank conditions make the observability analysis of a linear system simple and efficient. However, observability analysis become a significant challenge with regard to nonlinear systems and even to linear time-varying systems.

3.2. Nonlinear System with Memory

Now, consider a general continuous-time nonlinear system with memory

$$\begin{array}{cc}\hfill \dot{\mathit{x}}(t)& =\mathit{f}(\mathit{x}(t),\mathit{u}(t))\hfill \end{array}$$

(28a)

$$\begin{array}{cc}\hfill \mathit{y}(t)& =\mathit{h}(\mathit{x}(t),\mathit{x}(t-d_1),\cdots ,\mathit{x}(t-d_N))\hfill \end{array}$$

(28b)

where $\mathit{x}(t)\in {\mathbb{R}}^{n\times 1}$, $\mathit{u}(t)\in {\mathbb{R}}^{n_u\times 1}$ and $\mathit{y}(t)\in {\mathbb{R}}^{m\times 1}$ are the state, control input, and measurement vector at time t($t\ge N$), respectively, $d_i=i{D}_t$ and $D_t$ is regarded as the time step, and $\mathit{f}$ and $\mathit{h}$ are vector-valued functions.

The definition of observability of the continuous-time dynamic system with memory (28) is then as follows.

Definition 2.

The system is (locally) observable over the interval $[t_1,t_2]$ if the initial state $\mathit{x}(t_1-d_N)$ can be uniquely determined from $\mathit{y}(t)$, $t\in [t_1,t_2]$.

According to Equation (28a), there obviously exists a function between $\mathit{x}(t-\tau ),\tau =d_0,d_1,\cdots ,d_N$ and $\mathit{x}(t-d_N)$; hence, the measurement (28b) can be rewritten as

$$\mathit{y}(t)=\overline{\mathit{h}}(\mathit{x}(t-d_N),\overline{\mathit{u}}(t))$$

(29)

where $\overline{\mathit{u}}(t)$ is a combination of known control inputs $\mathit{u}$ from time $(t-d_N)$ to time t, denoted as $\overline{\mathit{u}}(t)=[\mathit{u}(t-d_N)\mathit{u}(t-d_{N-1})\cdots \mathit{u}(t)]$. For simplicity of notation, we write $\overline{\mathit{h}}(\mathit{x}(t-d_N))$ instead of Equation (29).

In this case, the observability rank condition is declared by the following theorem.

Theorem 2.

System (28) is (locally) observable if and only if the observability matrix

$$\mathit{O}(t)={\left[\begin{array}{cccc}\frac{\partial {\mathit{L}}_{\mathit{f}}^0\overline{\mathit{h}}(\mathit{x}(t-d_N))}{\partial \mathit{x}(t-d_N)}& \frac{\partial {\mathit{L}}_{\mathit{f}}^1\overline{\mathit{h}}(\mathit{x}(t-d_N))}{\partial \mathit{x}(t-d_N)}& \cdots & \frac{\partial {\mathit{L}}_{\mathit{f}}^{n-1}\overline{\mathit{h}}(\mathit{x}(t-d_N))}{\partial \mathit{x}(t-d_N)}\end{array}\right]}^T$$

(30)

is a full column matrix in which the Lie derivative is defined by

$$\begin{array}{cc}\hfill {\mathit{L}}_{\mathit{f}}^0\overline{\mathit{h}}(\mathit{x}(t-d_N))& =\overline{\mathit{h}}(\mathit{x}(t-d_N))\hfill \\ \hfill {\mathit{L}}_{\mathit{f}}^k\overline{\mathit{h}}(\mathit{x}(t-d_N))& =\frac{\partial {\mathit{L}}_{\mathit{f}}^{k-1}\overline{\mathit{h}}(\mathit{x}(t-d_N))}{\partial \mathit{x}(t-d_N)}\mathit{f}(\mathit{x}(t-d_N))\hfill \end{array}$$

(31)

Proof of Theorem 2.

Define the cumulative measurement vector $\mathcal{Y}(t)\in {\mathbb{R}}^{mn\times 1}$ as

$$\mathcal{Y}(t)={\left[\begin{array}{cccc}\mathit{y}(t)& \dot{\mathit{y}}(t)& \cdots & {\mathit{y}}^{(n-1)}(t)\end{array}\right]}^T$$

(32)

According to the implicit function theorem [18], the initial state $\mathit{x}(t-d_N)$ can be uniquely determined from measurements $\mathcal{Y}(t)$ if and only if

$$rank(\frac{\mathrm{d}\mathcal{Y}(t)}{\mathrm{d}\mathit{x}(t-d_N)})=n$$

(33)

With the chain rule, it is simple to obtain

$${\overline{\mathit{y}}}^{(k)}(t)={\mathit{L}}_{\mathit{f}}^k\overline{\mathit{h}}(\mathit{x}(t-d_N))$$

(34)

Hence, the system is locally observable if condition (30) holds.    □

This differential geometric method provides a derivative-based analytical observability analysis approach for nonlinear systems called Lie derivative-based observability analysis. However, it usually encounters computational difficulties when solving the higher-order Lie derivatives of complex nonlinear systems, which makes it impracticable in real applications.

In analytically intractable cases [13], the most widely-used approach is the empirical observability Gramian approach [12], which provides a numerical way to approximate the traditional observability Gramian (27) by adding small perturbations in each initial state and comparing the output for each perturbation.

Define ${\mathit{y}}^{+i}$ as the output with the ith nominal initial state ${\mathit{x}}_{0,i}$ affected by a positive disturbance $ϵ$ and let ${\mathit{y}}^{-i}$ be the output with ${\mathit{x}}_{0,i}$ affected by a negative disturbance $-ϵ$. The change in the measurement caused by perturbations in each initial state is denoted as $\Delta {\mathit{y}}_i={\mathit{y}}^{+i}-{\mathit{y}}^{-i}$, and the empirical observability Gramian can be computed by

$${\tilde{\mathit{W}}}_o(t)=\frac{1}{4{ϵ}^2}{\int }_0^t\left[\begin{array}{c}\Delta {\mathit{y}}_1^T\\ \Delta {\mathit{y}}_2^T\\ ⋮\\ \Delta {\mathit{y}}_r^T\end{array}\right]\left[\begin{array}{cccc}\Delta {\mathit{y}}_1& \Delta {\mathit{y}}_2& \cdots & \Delta {\mathit{y}}_r\end{array}\right]d\tau$$

(35)

where r is the number of states of interest, which implies that the empirical observability Gramian is computed by simulating the system $2r$ times in total.

3.3. Limitations of Traditional Observability Analysis

Although these traditional methods are effective at providing a binary answer to the observability question, they have disadvantages which limit their practicality.

First, due to lack of a unified assessment metric, different systems correspond to different observability rank conditions, leading to confusion in the understanding and use of observability analysis. Second, because multiple derivative calculations are generally required in the analytical method for nonlinear systems, applying them to complex models, such as our bioinspired flexible wing model, is impractical. Furthermore, the traditional methods fail to precisely quantify each state’s observability. However, in observability analysis it is crucial to determine which states are observable and which are not. Finally, the existing methods are formulated in a deterministic framework, and hence, are not able to assess the observability of stochastic systems.

Therefore, a derivative-free observability analysis method is proposed here to overcome the above weaknesses. It is able to remarkably reduce the computational burden, analyze the observability degree of each state, and evaluate the effect of noise on observability.

4. gPC-Based Observability Analysis for Stochastic System with Memory

In this section, based on a brief review of the gPC expansion, the surrogate model of cumulative measurement is formulated, then the observability-coefficient matrix is provided for qualitative analysis of observability. Consequently, observability indices are proposed to quantify the observability of the system. Finally, the equivalence between the traditional and gPC-based methods for linear systems with memory is proven.

4.1. The Generalized Polynomial Chaos Expansion

The generalized Polynomial Chaos expansion is an efficient uncertainty propagation method that represents the stochastic output with a weighted sum of orthogonal polynomial chaos basis functions of random input variables.

Let y and $\mathit{\xi }=\left[\begin{array}{cccc}{\xi }_1& {\xi }_2& \cdots & {\xi }_n\end{array}\right]$ denote the output and random input variables, respectively, following a known probability distribution. The stochastic output can be expressed as

$$y=\sum _{i=0}^{n_p}{\gamma }_i{\varphi }_i(\mathit{\xi })$$

(36)

where ${\varphi }_i$ is a polynomial chaos basis function, ${\gamma }_i$ is the ith polynomial chaos coefficient, $n_p=\frac{(n+p)!}{n!p!}-1$, and p is the maximum order of the polynomial chaos basis functions.

The mean $\mu$ and variance ${\sigma }^2$ of output y can be directly obtained as

$$\begin{array}{cc}\hfill \mu & ={\gamma }_0\hfill \\ \hfill {\sigma }^2& =\sum _{i=1}^{n_p}{\gamma }_i^2\hfill \end{array}$$

(37)

In general, higher orders yield higher accuracy. However, the number of unknown coefficients and the computational burden increase with the order. Many tests have proven that the increase in order has a negligible effect on the accuracy improvement when the order is larger than 2 [19]. Hence, the second order truncated gPC expansion is adopted in this paper, which is

$$y={\gamma }_0{\varphi }_0+\sum _{i=1}^n{\gamma }_i{\varphi }_1({\xi }_i)+\sum _{i=1}^n{\gamma }_{i,i}{\varphi }_2({\xi }_i^2)$$

(38)

where ${\varphi }_0$, ${\varphi }_1({\xi }_i)$, and ${\varphi }_2({\xi }_i^2)$ denote the zeroth-order, first-order, and second-order polynomial chaos bases, respectively, and ${\gamma }_0$, ${\gamma }_i$, and ${\gamma }_{i,i}$ represent the corresponding polynomial chaos coefficients. Note that only $2n+1$ polynomial chaos coefficients need to be determined.

Assuming that the random variable ${\xi }_i$ follows a standard normal distribution, the corresponding polynomial chaos bases are as follows [19]:

$$\begin{array}{cc}\hfill & {\varphi }_0=1\hfill \\ \hfill & {\varphi }_1(\xi )=\xi \hfill \\ \hfill & {\varphi }_2(\xi )=\frac{{\xi }^2-1}{\sqrt{2}}\hfill \end{array}$$

(39)

4.2. Observability Rank Condition

Consider a general nonlinear stochastic system with memory (28), where $\mathit{x}=[x_1{x}_2\cdots x_n]$ and the ith state random variable $x_i$ follows the known distribution.

Distinct from [15], the definition of observability of stochastic systems with memory is extended as follows.

Definition 3.

The system is (locally) observable over the interval $[t_1,t_2]$ if the initial state $\mathit{x}(t_1-N)$ can be inferred from $\mathit{y}(t)$, $t\in [t_1,t_2]$ and its solution satisfies a $95\%$ confidence interval.

In order to resolve probabilistic uncertainty propagation in such a system, a surrogate model is formulated by the following steps to represent the above relationship between the random initial state $\mathit{x}(t-d_N)$ and its corresponding cumulative measurement vector $\mathcal{Y}(t)$ (32).

First, we determine the mapping relationship between the ith element of the stochastic input variable $\mathit{x}(t-d_N)$ and a random variable ${\xi }_i$ following a standard normal distribution, which is

$$x_i(t-d_N)=F_i^{-1}(T_i({\xi }_i))$$

(40)

where $F_i$ is the cumulative probability function of $x_i(t-d_N)$, $F_i^{-1}$ is the inverse function of $F_i$, and $T_i$ denotes the cumulative probability function of ${\xi }_i$.

Next, collocations points (CPs) are selected to construct the surrogate model. CPs are a finite sample set of $\mathit{\xi }$, and the ith combination of CPs is expressed as ${\mathit{\xi }}_i=[{\xi }_{i,1}{\xi }_{i,2}\cdots {\xi }_{i,n}]$. The element values of the CPs are generated by the roots of one higher-order one-dimensional Hermite polynomial [19]. For example, three roots of the third order Hermite polynomial ${\varphi }_3(\xi )=\frac{{\xi }^3-3\xi }{\sqrt{6}}$, which are $\{-\sqrt{3},0,\sqrt{3}\}$, are used to to comprise the CPs, as a second order polynomial chaos expansion is adopted in this paper. Because there are $2n+1$ unknown polynomial chaos coefficients, $2n+1$ independent combinations of CPs are chosen randomly from the $3^n$ possible combinations, which results in the matrix $\mathsf{\Xi }\in {\mathbb{R}}^{(2n+1)\times n}$, as follows:

$$\mathsf{\Xi }=\left[\begin{array}{cccc}{\xi }_{1,1}& {\xi }_{1,2}& \cdots & {\xi }_{1,n}\\ {\xi }_{2,1}& {\xi }_{2,2}& \cdots & {\xi }_{2,n}\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_{2n+1,1}& {\xi }_{2n+1,2}& \cdots & {\xi }_{2n+1,n}\end{array}\right]$$

(41)

which is the full rank where the value of the ith element of the sth CP ${\xi }_{s,i}$ is chosen randomly from $\{-\sqrt{3},0,\sqrt{3}\}$.

Then, $2n+1$ samples of stochastic input variables are transformed from $2n+1$ CPs based on Equation (40). Substituting them into Equation (32) in place of the original input variables $\mathit{x}(t-d_N)$, the output matrix $\mathit{Y}(t)\in {\mathbb{R}}^{(2n+1)\times mn}$ is obtained, where $\mathit{Y}(t)={[{\mathcal{Y}}_1(t){\mathcal{Y}}_2(t)\cdots {\mathcal{Y}}_{2n+1}(t)]}^T$ and ${\mathcal{Y}}_i$ represents the measurement output from the ith sample.

Consequently, according to Equation (38), the surrogate model is provided by

$$\mathit{Y}(t)=\mathit{H}\mathbf{\Gamma }(t)$$

(42)

where the basis matrix of polynomial chaos bases $\mathit{H}\in {\mathbb{R}}^{(2n+1)\times (2n+1)}$ is provided by Equation (43) and the coefficient matrix of polynomial chaos coefficients $\mathbf{\Gamma }(t)\in {\mathbb{R}}^{(2n+1)\times mn}$ is Equation (44); here, ${\gamma }_0^l$, ${\gamma }_i^l$, and ${\gamma }_{i,i}^l$ represent the polynomial chaos coefficients with respect to the ith state and lth measurement:

$$\mathit{H}=\left[\begin{array}{ccccccccc}{\varphi }_0& {\varphi }_1({\xi }_{1,1})& {\varphi }_1({\xi }_{1,2})& \cdots & {\varphi }_1({\xi }_{1,n})& {\varphi }_2({\xi }_{1,1}^2)& {\varphi }_2({\xi }_{1,2}^2)& \cdots & {\varphi }_2({\xi }_{1,n}^2)\\ {\varphi }_0& {\varphi }_1({\xi }_{2,1})& {\varphi }_1({\xi }_{2,2})& \cdots & {\varphi }_1({\xi }_{2,n})& {\varphi }_2({\xi }_{2,1}^2)& {\varphi }_2({\xi }_{2,2}^2)& \cdots & {\varphi }_2({\xi }_{2,n}^2)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ {\varphi }_0& {\varphi }_1({\xi }_{2n+1,1})& {\varphi }_1({\xi }_{2n+1,2})& \cdots & {\varphi }_1({\xi }_{2n+1,n})& {\varphi }_2({\xi }_{2n+1,1}^2)& {\varphi }_2({\xi }_{2n+1,2}^2)& \cdots & {\varphi }_2({\xi }_{2n+1,n}^2)\end{array}\right]$$

(43)

$$\mathbf{\Gamma }(t)=\left[\begin{array}{cccc}{\gamma }_0^1& {\gamma }_0^2& \cdots & {\gamma }_0^{mn}\\ {\gamma }_1^1& {\gamma }_1^2& \cdots & {\gamma }_1^{mn}\\ ⋮& ⋮& \ddots & ⋮\\ {\gamma }_n^1& {\gamma }_n^2& \cdots & {\gamma }_n^{mn}\\ {\gamma }_{1,1}^1& {\gamma }_{1,1}^2& \cdots & {\gamma }_{1,1}^{mn}\\ ⋮& ⋮& \ddots & ⋮\\ {\gamma }_{n,n}^1& {\gamma }_{n,n}^2& \cdots & {\gamma }_{n,n}^{mn}\end{array}\right]$$

(44)

Because $\mathsf{\Xi }$ is full rank, it is obvious that $\mathit{H}$ is nonsingular and its inverse matrix ${\mathit{H}}^{-1}$ exists. Hence, the unknown polynomial chaos coefficients can be obtained by

$$\mathbf{\Gamma }(t)={\mathit{H}}^{-1}\mathit{Y}(t)$$

(45)

According to Equation (37) and the orthogonal property [20], the zero-order polynomial chaos coefficient ${\gamma }_0^l$ denotes the mean of the estimated value for the lth measurement and the first and second order polynomial chaos coefficients ${\gamma }_i^l$ and ${\gamma }_{i,i}^l$ represent the contribution of the ith state to the uncertainty of the lth measurement. Supposing that the contribution of the ith state to the uncertainty of the lth measurement is zero, the lth measurement is unaffected by the change in this state, which implies that the ith state cannot be uniquely determined from the lth measurement. In order to uniquely infer n states from the given measurements, n valid measurements, for which the contributions of the states to the uncertainties of the estimated values for these measurements are linearly independent, are necessary.

Therefore, we define the observability coefficient matrix $\mathbf{\Phi }(t)\in {\mathbb{R}}^{2n\times mn}$ as follows:

$$\mathbf{\Phi }(t)=\left[\begin{array}{cccc}{\gamma }_1^1& {\gamma }_1^2& \cdots & {\gamma }_1^{mn}\\ ⋮& ⋮& \ddots & ⋮\\ {\gamma }_n^1& {\gamma }_n^2& \cdots & {\gamma }_n^{mn}\\ {\gamma }_{1,1}^1& {\gamma }_{1,1}^2& \cdots & {\gamma }_{1,1}^{mn}\\ ⋮& ⋮& \ddots & ⋮\\ {\gamma }_{n,n}^1& {\gamma }_{n,n}^2& \cdots & {\gamma }_{n,n}^{mn}\end{array}\right]$$

(46)

and the observability rank condition is stated by the following theorem.

Theorem 3.

System (28) is (locally) observable if and only if the observability coefficient matrix $\mathbf{\Phi }(t)$ has n linearly independent column vectors and the ith and $(n+i)$th rows are not zero vectors.

Note that by using a surrogate model to represent ${\mathcal{Y}}_k$ in Equation (23), the proposed gPC-based observability analysis method can solve the observability problem of linear and nonlinear discrete-time systems with memory.

The detailed gPC-based observability analysis procedure is summarized in Algorithm 1.

Algorithm 1 gPC-based Observability Analysis Procedure

1: Determine the mapping between ith random variable $x_i(t-d_N)$ and a normal variable ${\xi }_i$ via Equation (40); 2: Select CPs (41) based on the linear independence method and substitute them into basis matrix (43); 3: Transform all CPs into the samples of stochastic input variables based on the mapping (40); 4: Compute output matrix $\mathit{Y}(t)$; 5: Calculate the coefficient matrix (45) and the observability-coefficient matrix $\mathbf{\Phi }(t)$; 6: Analyze Observability according to the observability rank condition Theorem 3.

4.3. Degree of Observability

The observability rank condition can only draw a binary conclusion as to whether the system is observable. However, it is limited to quantifying the observability of the system. Hence, three more quantitative observability indices are proposed to further analyze the observability degree.

1.

Condition Number: The condition number $\kappa (\mathbf{\Phi })$ is the ratio of the maximum singular value to the minimum singular value, which is utilized to show the numerical stability of the system states, and is provided by

$$\kappa (\mathbf{\Phi })=\frac{{\sigma }_{max}(\mathbf{\Phi })}{{\sigma }_{min}(\mathbf{\Phi })}$$

(47)

where ${\sigma }_{max}(\mathbf{\Phi })$ and ${\sigma }_{min}(\mathbf{\Phi })$ are the maximal and minimal singular values of $\mathbf{\Phi }$, respectively.

Suppose that the condition number is enormous or even goes to infinity, the observability-coefficient matrix $\mathbf{\Phi }$ is ill-conditioned, and the system is weakly observable; this implies that certain states are hard to observe. However, if the condition number is close to one, $\mathbf{\Phi }$ is well-conditioned and the system is brawny observable. Note that all system states are considered equally observable when the condition number equals one.

2.

Contribution Rate: Two different equations for the contribution rate are provided in this paper, both of which are regarded as quantitative indices of each state’s observability. The first contribution rate, ${\chi }_i^1$, is defined as the contribution of the ith state to the uncertainty of all measurements. According to the surrogate model, it is provided by

$${\chi }_i^1=\frac{{\sum }_{l=1}^{mn}{({\gamma }_i^l)}^2+{({\gamma }_{i,i}^l)}^2}{{\sum }_{i=1}^n{\sum }_{l=1}^{mn}{({\gamma }_i^l)}^2+{({\gamma }_{i,i}^l)}^2}$$

(48)

which represents the influence of a specific state on the measurements. The ith state is brawnier observable when its first contribution rate ${\chi }_i^1$ is closer to one. Instead, the ith state is weakly observable as ${\chi }_i^1$ approaches zero.

The second contribution rate, ${\chi }_i^2$, is the maximum contribution of the ith state to each measurement. We define the proportion of the contribution of the ith state to the variance of the lth measurement as ${\chi }_{i,l}^2$, that is,

$${\chi }_{i,l}^2=\frac{{({\gamma }_i^l)}^2+{({\gamma }_{i,i}^l)}^2}{{\sum }_{i=1}^n{({\gamma }_i^l)}^2+{({\gamma }_{i,i}^l)}^2}$$

(49)

where the numerator means the contribution of the ith state to the variance of the lth measurement and the denominator denotes the variance of the lth measurement.

Hence, ${\chi }_i^2$ equals the maximum of $({\chi }_{i,1}^2,{\chi }_{i,2}^2,\cdots ,{\chi }_{i,mn}^2)$. As in the first contribution rate, the larger the second contribution rate is, the more the state is brawnier observable.

3.

Interference Rate: If the contribution of the ith state to the variance of the lth measurement, ${({\gamma }_i^l)}^2+{({\gamma }_{i,i}^l)}^2$, is smaller than the measurement noise variance ${\sigma }_v^2$, the changes in measurements contain too much environmental interference, and the initial state is difficult to distinguish from noisy measurements with a high confidence level.

Therefore, we define the proportion of the measurement noise variance ${\sigma }_v^2$ to the variance of the lth measurement as the interference rate:

$$V_l=\frac{{\sigma }_v^2}{{\sum }_{i=1}^n{({\gamma }_i^l)}^2+{({\gamma }_{i,i}^l)}^2}$$

(50)

If the interference rates $V_1,V_2,\cdots ,V_{mn}$ are larger than all the corresponding contribution of any states ${\chi }_{i,1}^2,{\chi }_{i,2}^2,\cdots ,{\chi }_{i,mn}^2$, the system is considered to be weakly observable due to noisy measurements.

To simplification simulation, an equivalent index called the interference binary value $\mathsf{{\rm Y}}$ is used:

$$\mathsf{{\rm Y}}=\left\{\begin{array}{c}0\underset{i}{max}\underset{l}{min}{V}_l/{\chi }_{i,l}^2\le 1\hfill \\ 1otherwise\hfill \end{array}\right.$$

(51)

When $\mathsf{{\rm Y}}$ equals zero, the measurement noise has negligible effect on observability analysis. Inversely, the states are hard to infer from noisy measurements when $\mathsf{{\rm Y}}$ equals one.

4.4. Equivalence between the Traditional and Proposed Approaches

To further clarify the rationality of the proposed method from another perspective, the connection between the traditional and proposed methods for a general linear system with memory (20) is demonstrated as a particular case.

For a general linear system with memory (20), the relationship between the random initial state ${\mathit{x}}_{k-N}$ and the cumulative measurement ${\mathcal{Y}}_k$ at time k is provided by Equation (24). We define $\Delta {\mathit{x}}_{k-N}\in {\mathbb{R}}^{n\times 1}$ and $\Delta {\mathcal{Y}}_k\in {\mathbb{R}}^{m\times 1}$ as the state and measurement uncertainty vectors, respectively; as the control input vectors $\mathit{u}$ from time $k-N$ to $k-1$ are determined at time k, the uncertainty propagation function is

$$\Delta {\mathcal{Y}}_k=\mathit{O}\Delta {\mathit{x}}_{k-N}$$

(52)

where $\mathit{O}$ is provided by Equation (21).

Assuming that there exists an open neighborhood U of ${\mathit{x}}_{k-N}$ and that ${\mathit{x}}_{k-N}+\Delta {\mathit{x}}_{k-N}\in U$, we can find that ${\mathit{x}}_{k-N}$ and ${\mathit{x}}_{k-N}+\Delta {\mathit{x}}_{k-N}$ are distinguishable from their corresponding measurements if and only if the observability matrix $\mathit{O}$ has full column ranking based on Theorem 4, which implies that system (20) is locally observable at time k.

Theorem 4.

Consider a non-homogeneous linear system of equations as

$$\mathit{A}\mathit{x}=\mathit{b}$$

(53)

where $\mathit{A}\in {\mathbb{R}}^{m\times n}$, $\mathit{x}\in {\mathbb{R}}^{n\times 1}$, and $\mathit{b}\in {\mathbb{R}}^{m\times 1}$.

Then, there exists a unique solution if and only if the rank of the coefficient matrix $\mathit{A}$ is n.

Taking n samples of stochastic input variables into consideration, the above equation can be rewritten as

$$\Delta \mathit{Y}=\Delta \mathit{X}\mathit{O}$$

(54)

where $\Delta \mathit{Y}={[\Delta {\mathcal{Y}}_k^1,\Delta {\mathcal{Y}}_k^2,\cdots ,\Delta {\mathcal{Y}}_k^n]}^T$ and $\Delta \mathit{X}={[\Delta {\mathit{x}}_{k-N}^1,\Delta {\mathit{x}}_{k-N}^2,\cdots ,\Delta {\mathit{x}}_{k-N}^n]}^T$. $\Delta {\mathcal{Y}}_k^i$, where $\Delta {\mathit{x}}_k^i$ represents the measurement and state uncertainty vector from ith samples at time k.

In the proposed approach, the surrogate model of the relationship between ${\mathit{x}}_{k-N}$ and ${\mathcal{Y}}_k$ is expressed as

$${\mathit{Y}}_k=\mathit{H}{\mathbf{\Gamma }}_k$$

(55)

Because the second-order polynomial chaos coefficients are zeros for a linear system, the surrogate model is simplified as

$${\mathit{Y}}_k=\overline{\mathit{H}}{\overline{\mathbf{\Gamma }}}_k$$

(56)

where

$$\overline{\mathit{H}}=\left[\begin{array}{ccccc}{\varphi }_0& {\varphi }_1({\xi }_{1,1})& {\varphi }_1({\xi }_{1,2})& \cdots & {\varphi }_1({\xi }_{1,n})\\ {\varphi }_0& {\varphi }_1({\xi }_{2,1})& {\varphi }_1({\xi }_{2,2})& \cdots & {\varphi }_1({\xi }_{2,n})\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\varphi }_0& {\varphi }_1({\xi }_{n+1,1})& {\varphi }_1({\xi }_{n+1,2})& \cdots & {\varphi }_1({\xi }_{n+1,n})\end{array}\right]$$

(57)

$${\overline{\mathbf{\Gamma }}}_k=\left[\begin{array}{cccc}{\gamma }_0^1& {\gamma }_0^2& \cdots & {\gamma }_0^{mn}\\ {\gamma }_1^1& {\gamma }_1^2& \cdots & {\gamma }_1^{mn}\\ ⋮& ⋮& \ddots & ⋮\\ {\gamma }_n^1& {\gamma }_n^2& \cdots & {\gamma }_n^{mn}\end{array}\right]$$

(58)

As mentioned above, ${\gamma }_0^l$ is the mean of the lth measurement and ${\gamma }_i^l$ stands for the contribution of the ith state to the uncertainty of the lth measurement. Hence,

$$\Delta \mathit{Y}=\Delta \mathit{X}\overline{\mathbf{\Phi }}$$

(59)

where

$$\overline{\mathbf{\Phi }}=\left[\begin{array}{cccc}{\gamma }_1^1& {\gamma }_1^2& \cdots & {\gamma }_1^{mn}\\ ⋮& ⋮& \ddots & ⋮\\ {\gamma }_n^1& {\gamma }_n^2& \cdots & {\gamma }_n^{mn}\end{array}\right]$$

(60)

System (20) is locally observable if and only if $rank(\overline{\mathbf{\Phi }})=n$.

The same observability rank condition for linear systems with memory are derived from the traditional and proposed method. This convincingly proves that gPC-based observability analysis is equivalent to the traditional method.

5. Simulations

In this section, the observability of a Lorenz system with memory is first analyzed by the traditional and proposed methods to validate the feasibility of gPC-based observability analysis. Then, we apply the proposed method to observability analysis of a bioinspired flapping wing system based on the hawkmoth, Manduca sexta.

5.1. Lorenz System with Memory

Consider the following Lorenz system with memory:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\dot{x}}_1(t)=-10x_1(t)+10x_2(t)\hfill \\ {\dot{x}}_2(t)=28x_1(t)+x_1(t)x_3(t)-x_2(t)\hfill \\ {\dot{x}}_3(t)=x_1(t)x_2(t)-\frac{8}{3}{x}_3(t)\hfill \end{array}\right.\hfill \\ \hfill & \left\{\begin{array}{c}{y}_1(t)=x_1(t)+\frac{1}{2}{x}_1(t-d_1)+\frac{1}{4}{x}_1(t-d_2)\hfill \\ y_2(t)=x_2(t)+\frac{1}{2}{x}_2(t-d_1)+\frac{1}{4}{x}_2(t-d_2)\hfill \end{array}\right.\hfill \end{array}$$

(61)

with starting time $t_0=0s$, final time $t_f=10s$, time step $dt=0.01s$, and initial state $x_1(0)=x_2(0)=x_3(0)=1$.

5.1.1. Lie Derivative-Based Observability Analysis

The observability analysis is first performed by using the Lie derivative-based approach. The observability matrix $\mathit{O}(t)$ is calculated by Equation (32), and its rank is shown in

Table 1. The rank of observability matrices calculated via three approaches.

Observability Matrix	Rank
Lie Derivative-Based Observability Matrix	3
Empirical Observability Gramian Matrix	3
gPC-based Observability Coefficient Matrix	6
gPC-based First Order Observability Coefficient Matrix	3

. In addition, the condition number of $\mathit{O}(t)$ is shown in Figure 4 to measure the observability degree.

It can be seen that the system is observable, as the observability matrix is full rank and the value of the condition number changes drastically and reaches large values at intervals, revealing that the system is weakly observable. However, this approach can neither quantify the degree of observability for each state nor account for the effect of the observation noise, thereby reducing its reliability in practice.

5.1.2. Empirical Observability Gramian Analysis

Empirical observability Gramian is the most widely used method in practical applications due to its avoidance of the need to calculate the complicated Lie derivatives of nonlinear systems. The rank of the observability matrix ${\tilde{\mathit{W}}}_o(t)$ computed by Equation (35) is shown in Table 1, and the condition number is present in Figure 5.

The same conclusion that the system is weakly observable is drawn via empirical observability Gramian analysis, as the observability matrix maintains full rank while the condition numbers are far from one. However, the condition number is more stable than that computed by the Lie derivative-based approach, most likely due to the inaccuracy of the approximation.

5.1.3. gPC-Based Observability Analysis

Here, the proposed observability analysis approach is utilized. The observability coefficient matrix is computed through Equation (44). In the observability coefficient matrix, the first order polynomial chaos coefficients are larger and more significant than the second order polynomial chaos coefficients [16]; thus, we name the first n rows of the observability coefficient matrix as the first order observability coefficient matrix. The ranks of these two observability matrices are shown in Table 1. It can be seen that both of the observability coefficient matrices are full rank at all times, meaning that the system is observable.

Then, the first and second contribution rates of each state are shown in Figure 6.

It is evident from Figure 6 that the contribution rate of the third state is much lower than that of the first two states and varies periodically with time, which implies that the third state is weakly observable while the others are brawny observable.

This conclusion can be proven by the condition number. The condition number of the first order coefficient matrix is displayed in Figure 7a. It can be seen that it has enormous values and a similar change to the contribution rate of the third state. This is because the weak observability of the third state has a negative impact on the observability of the whole system, thereby resulting in an ill-conditioned matrix.

The coefficients of the states with brawny observability extracted from the first-order coefficient matrix constitute the active coefficient matrix. Its condition number is presented in Figure 7b, showing that the condition number of the active coefficient matrix is very close to one. This implies that the active coefficient matrix is well-conditioned and the first two states are brawny observable. Therefore, it can be concluded that the condition number is strongly affected by the weak observability of any state. When the condition number has large values, at least one state of the system is weakly observable.

In addition, the effect of the measurement noise on the observability is shown in Figure 8. Two different noise variances are assumed, $1\times {10}^{-9}$ and $1\times {10}^{-11}$. When the noise variance becomes larger, more interference binary values within the simulation time equal one, which implies that the observability of the system is weaker.

Finally, the computation complexity of the three approaches is compared by computation time in

Table 2. Computation time of the three approaches.

Approach	Computation Time (s)
Lie Derivative-Based Observability Analysis	0.7940
Empirical Observability Gramian Analysis	87.1700
gPC-based Observability Analysis	1.7980

. The computation time was computed using MATLAB R2021b on a 2.70 GHz Intel(R) Core(TM) i7 processor with 8 GB RAM.

Though Lie derivative-based observability analysis has the shortest time cost in the case where the derivation is completed, it is limited for systems with high nonlinearity in actual application because of the complicated derivative calculation. The proposed method effectively reduces the computational cost compared to empirical observability Gramian analysis.

5.2. Bioinspired Flexible Flapping Wing System

Considering the system model as Equation (14) and the measurement model as Equation (19), the proposed observability analysis approach is used to determine whether the neural encoding strain measurements are sufficient to reconstruct the wing rotation rates (P, Q, R). The observability analysis plays a vital role in the optimal placement of sensors on the wing in practice, which is further discussed in the next section.

The simulation of the full bioinspired flexible flapping wing model requires several steps. First, the wing-beat period and the total simulation time are assumed to be $T_{beat}$ and $2T_{beat}$, respectively. The time step is set as $T_{beat}/50$. The first period is used to fill the memory, and the observability analysis is based on the measurements at the second period.

Second, the control input sequence $\mathit{u}(t)={\left[\dot{P}(t)\dot{Q}(t)\dot{R}(t)\right]}^T$ is generated by derivation of Equation (15) with respect to time t. According to [21,22], the Euler angles are provided by

$$\begin{array}{cc}\hfill \theta (t)& =0\hfill \\ \hfill \alpha (t)& =\frac{\pi }{2}-A_{\alpha }tanh(\frac{\pi }{2}sin(\frac{2\pi t}{T_{beat}}))\hfill \\ \hfill \zeta (t)& =-A_{\zeta }cos(\frac{2\pi t}{T_{beat}})\hfill \end{array}$$

(62)

where $A_{\zeta }$ is the position angle amplitude and $A_{\alpha }$ is the feathering angle amplitude.

Then, the stiff matrix $\mathbf{\Omega }$, applied acceleration mass matrix ${\mathit{M}}_a$, and aerodynamics force $\mathit{Q}$ are computed using the wing structural mode shapes ${\phi }_i(x,y)$, their corresponding $f_i$, the mass density $\rho (x,y)$, and the thickness $h(x,y)$. We assume that the mass density of the wing is constant and model the thickness as an exponential decrease from root to tip and from the leading edge (LE) to the trailing edge (TE) [23], that is,

$$h(x,y)=\frac{1}{2}(t_{LE}{e}^{-a_{LE}\frac{y}{c_{wing}}}+t_{RO}{e}^{-a_{RO}\frac{x-x_{LE}(y)}{d_{wing}}})$$

(63)

where $t_{LE}$ and $t_{RO}$ are the leading edge and the root thickness, respectively, $c_{wing}$ and $d_{wing}$ are the chord and spanwise length, respectively, $a_{LE}$ and $a_{RO}$ are the respective decay rates along the two directions, $x_{LE}(y_0)$ means the x-ordinate of the intersection point at the leading edge, and $y=y_0$.

Finally, bysubstituting these pre-processing data and the detailed parameters listed in

Table 3. Detailed parameters used in the simulation.

Parameter	Symbol	Value	Unit
Wing-Beat Period	$T_{beat}$	40	ms
Memory Length	$t_M$	40	ms
Feathering Angle Amplitude	$A_{\alpha }$	45	deg
Position Angle Amplitude	$A_{\zeta }$	60	deg
Mass Density	$\rho$	220	kg/m${}^3$
Chord Length	$c_{wing}$	22	mm
Spanwise Length	$d_{wing}$	50	mm
Leading Edge Thickness	$t_{RO}$	$0.26$	mm
Root Thickness	$t_{RO}$	$0.26$	mm
Decay Rate from LE to TE	$a_{LE}$	$0.24$	−
Decay Rate from Root to Tip	$a_{RO}$	$0.24$	−
Mode Frequency of First Bending Mode	$f_1$	50	Hz
Mode Frequency of First Torsion Mode	$f_2$	55	Hz
Rotation Axis	$x_r$	$-4.68$	mm
Air Fluid Density	${\rho }_f$	$1.225$	kg/m${}^3$
Delay of STA Function	a	5	ms
Width of STA Function	b	4	ms
STA Frequency	$f_{STA}$	$159.155$	Hz
Slope of NLA Function	c	$0.15$	−
Half-Maximum position of NLA Function	d	$0.005$	−
Normalization Constant	$C_{\varkappa }$	$1.084\times {10}^{-4}$	−

into Equations (14) and (19), the full bioinspired flexible flapping wing model is constructed and the observability indices are calculated using the gPC-based observability analysis approach.

To account for the observability changes on different locations of different sensors, the wing’s surface is meshed into a $50\times 22$ grid. One shear sensor and one bending sensor are placed at each point of intersection on the wing. Here, the first contribution rate is adopted to describe the observability degree of each state, while the minimum contribution rate provides a measure of the weakest observable state.

The averaged contribution rate of the three wing rotation rates for all sensors is shown in Figure 9.

The above illustrates that R is the weakest observable state when only shear or bending strain measurements are obtained, while P and Q are brawnier observable. In addition, shear strain measurements result in brawnier observability of the whole wing model.

On the other hand, bending strain measurements provide the smaller condition number, with a mean value of 52.5629, while the averaged condition number of the shear strain measurements is 87.049, which implies that the observability of all rotation rates are more balanced when using bending strain measurements.

Then, the averaged minimum contribution rate for each sensor location and type within the simulation time is obtained and normalized by the max-min normalization method, depicted in Figure 10.

This indicates that the rotation rates can be accurately estimated when bending strain sensors are concentrated in the wing root and shear strain sensors are placed in the upper part of the trailing edge, as because the bending strain at the wing root and shear strain at the upper part of the trailing edge are the largest and contain more information about the rotation rates. In addition, bending strain sensors located in the wing root provide more measurements to encode the rotation rate compared to the shear strain sensors at their best location.

Consequently, we chose five typical locations, namely, the root, the tip, the middle of the leading edge, the middle of the trailing edge, and the center of the wing, to show the changes in observability degree during the simulation time. The minimum contribution rates of the rotation rates for the two types of sensors are shown in Figure 11 and Figure 12.

It can be seen that the bending strain sensor at the wing root results in the brawniest observability, while the shear strain sensor at the same place leads to the opposite result. For shear strain sensors, similar observability degrees can be achieved from the other four locations; in particular, the sensors at the tip and the middle of the trailing edge are capable of obtaining the most valuable measurements. However, the effective bending measurements are mainly obtained from the sensor at the root. This result is the same as that demonstrated in Figure 10.

The condition numbers of the rotation rates for each sensor location and type are shown in Figure 13 and Figure 14. It can be seen that the condition number for each sensor location and type varies dramatically over time. The observability of the states is primarily more balanced when the shear strain sensors are located at the root and the bending strain sensors are located at the middle of the trailing edge.

When the sensors are assumed to be inaccurate, the noise variances of the shear and bending strain sensors are $1\times {10}^{-18}$ and $1\times {10}^{-19}$, respectively; the interference binary values from noisy measurements are shown in Figure 15 and Figure 16.

The figures reveal that the observability based on the bending strain measurements is more sensitive to measurement noise than that based on noisy shear strain measurements. Furthermore, for sensors of the same type in different locations, the same measurement noise shows a different effect on the observability. This is because the same sensors at different locations encode different information about the wing rotation rates.

These observability indices provide the main metrics for optimizing the placement of sensors, which is discussed in the next section.

6. Observability-Based Optimal Sensor Placement

Here, all sensors are assumed to be located at the wing veins, as this is more consistent with the characteristics of the hawkmoth [24]. The optimal sensor problem is to find the set of sensors that can encode the most information about wing rotation rates based on the proposed observability indices when the possible location and type are given.

As mentioned in Section 4.3, the observability degrees of the rotation rates are more balanced when the condition number is smaller, and the system becomes more observable as the minimum contribution rate becomes larger. In addition, only when the interference binary value equals zero can the negative effect of measurement noise on observability be neglected. Thus, we define the observability-based criterion for the ith sensor as $J_i$, which is

$$J_i={\kappa }_i+{\omega }_{\varpi }\varpi ({\chi }_{min,i}^1)+{\omega }_{\mathsf{{\rm Y}}}{\mathsf{{\rm Y}}}_i$$

(64)

where the normalized condition number of the ith sensor is denoted as ${\kappa }_i$, ${\chi }_{min,i}^1$ stands for the minimum contribution rate of the ith sensor, the max-min normalization function is $\varpi (x)=(max-x)/(max-min)$, ${\mathsf{{\rm Y}}}_i$ is the interference binary value of the ith sensor, ${\omega }_{\varpi }$ and ${\omega }_{\mathsf{{\rm Y}}}$ are the weights, and $J_i$ shows the ith sensor’s ability to encode information about the wing rotation rates, which increases as $J_i$ decreases.

To obtain the optimal sensor placement, a set of possible locations of sensors is first determined by sampling the points at every 2 mm along each vein. Then, the observability-based optimal sensor placement problem can be modeled as

$$\begin{array}{cc}\hfill & \underset{{\beta }_i}{min}\sum _{i=1}^p{\beta }_i{J}_i\hfill \\ \hfill & s.t.\left\{\begin{array}{c}{\beta }_i\in 0,1\hfill \\ \sum _{i=1}^p{\beta }_i=r\hfill \end{array}\right.\hfill \end{array}$$

(65)

where ${\beta }_i$ is a binary activation function, analogous to an off/on switch, and the number of potential sensors in the set $p=336$, r is the desired number of sensors for placement on the wing.

Based on the convex optimization problem and the traversal method, r optimal locations out of all potential places on the wing veins are obtained at each time step. The optimal locations change with flexible wing motion and external noise during flight, and thus each possible point’s cumulative number for being an optimal location is counted. A larger cumulative number indicates that a strain sensor placed here is more conducive to improving the observability of the flexible wing flapping system by following the estimation accuracy improvement of rotation rates. The areas where more optimal locations are clustered suggest the most suitable places to arrange the strain sensors.

We assume that the sensors are inaccurate and that measurements inevitably contain noise, as sensors are usually imprecise in practical applications due to technical or external environmental reasons. The noise variances of the shear and bending strain measurements are set as $1\times {10}^{-18}$ and $1\times {10}^{-19}$, respectively. Let ${\omega }_{\varpi }$ be 5, ${\omega }_{\mathsf{{\rm Y}}}$ be 100, and r be 20; then, each possible location’s cumulative number for being an optimal location during the simulation is plotted in Figure 17 (when over 5).

The figure shows that the optimal locations for strain sensors on the wing are divided into three groups: one near the root, one in the center, and one on the upper part of the trailing edge. The wing root is the best location, followed by the upper part of the trailing edge, then the center.

When the weight of the minimum contribution rate ${\omega }_{\varpi }$ is reduced to 1, the numerical stability of the states becomes more critical. The cumulative number for each possible location to be optimal in this case is shown in Figure 18. Obviously, more optimal locations are concentrated at the wing root, while the number at the other locations decreases.

With $r=30$ the numbers for all three groups are increased, as illustrated in Figure 19.

When the noise variance of the shear and bending strain measurements diminishes to $1\times {10}^{-19}$ and $1\times {10}^{-20}$, respectively, the best locations on the wing are those shown in Figure 20. In this situation, the numbers for the best locations at the center and upper part of the wing increase.

Finally, process noise with the variance of $\left[0.010.010.010.01111\right]$ is considered, which is mainly caused by inaccuracy of the dynamic model and unknown control inputs. The best locations on the wing are shown in Figure 21.

The best sensor locations are again clustered into three similar areas, although now most of them are located on the wing root. This is because the observability of the wing system is the brawniest, and the least affected by process noise when the sensors are located at the wing root.

Although the optimal sensor placement varies slightly with the different parameter selection, it remains concentrated in three main areas of the wing: one at the wing root, one in the center, and one on the upper part of the trailing edge, which is consistent with the measured locations of the campaniform sensilla on the hawkmoth wing depicted in Figure 22 [24]. The similar optimal placement we obtained based on the gPC-based observability analysis approach proves its validity and practicality.

7. Discussion and Conclusions

This paper proposes a derivative-free observability analysis approach based on the generalized Polynomial Chaos expansion for stochastic systems with memory. The observability rank condition provides a binary answer to the observability question, followed by several observability indices to describe the observability degree of the system from different perspectives. For instance, the contribution rate quantifies the observability degree of each state, the condition number presents the numerical stability of all states, and the interference rate measures the influence of measurement noise on the observability. In addition, the equivalence between the traditional and proposed approach for linear systems with memory is proven. The effectiveness of the proposed approach is mathematically demonstrated by applying three different approaches to analysis of the Lorenz system and comparing the results. Sequentially, the proposed approach is utilized to analyze the observability of a flexible wing system and determine the optimal sensor placement. The results show that bending strain sensors should be located at the wing root, while shear strain sensors should be placed at the center and upper part of the wing. The optimal placement we obtained based on the proposed method is similar to the natural distribution of campaniform sensilla on the wing of the hawkmoth.