*Article* **Dynamic Deformation Reconstruction of Variable Section WING with Fiber Bragg Grating Sensors**

#### **Zhen Fu 1, Yong Zhao 2, Hong Bao 1,\* and Feifei Zhao <sup>1</sup>**


Received: 27 June 2019; Accepted: 26 July 2019; Published: 30 July 2019

**Abstract:** In order to monitor the variable-section wing deformation in real-time, this paper proposes a dynamic reconstruction algorithm based on the inverse finite element method and fuzzy network to sense the deformation of the variable-section beam structure. Firstly, based on Timoshenko beam theory and inverse finite element framework, a deformation reconstruction model of variable-section beam element was established. Then, considering the installation error of the fiber Bragg grating (FBG) sensor and the dynamic un-modeled error caused by the difference between the static model and dynamic model, the real-time measured strain was corrected using a solidified fuzzy network. The parameters of the fuzzy network were learned using support vector machines to enhance the generalization ability of the fuzzy network. The loading deformation experiment shows that the deformation of the variable section wing can be reconstructed with the proposed algorithm in high precision.

**Keywords:** deformation reconstruction; fuzzy system; inverse finite element method; strain modification; variable section wing

#### **1. Introduction**

Through the real-time perception of the deformation of the wing structure, the flight attitude is adjusted in real time by the actuation and control systems to ensure the safety of the aircraft. This is one of the key technologies for the next generation of intelligent space vehicles. Using the strain sensors arranged on the wing structure to obtain the strain data of the structure in real time, and reconstructing the deformation information of the structure using the specific strain displacement relation [1], that is, shape sensing, it has become a research hotspot of the wing structure deformation. Due to the advantages of small size, anti-electromagnetic interference and the ability to form a large-scale quasi-distributed sensing network [2], the FBG strain sensors are widely used in the field of shape sensing.

The key to shape sensing is constructing a relationship between the structural deformation and the strain measurement. There are many research methods of strain-based deformation reconstruction proposed by domestic and foreign scholars, such as the modal transformation method, piecewise linearization method, and finite element method. The modal transformation method can accurately reconstruct the deformation of plate and beam structure [3], but it needs accurate finite element model. Another problem is that the algorithm is difficult to apply to the large deformation of geometric nonlinear structures because it is based on the principle of linear superposition [4]. The Ko method is based on the classical Euler Bernoulli beam theory [5]. By integrating the discrete surface strain measurements with piecewise continuous polynomials, high-precision reconstruction of the beam element in a one-dimensional direction can be achieved. The Ko method can be used for deformation

reconstruction of cantilever beam structure, and ground experiment for deformation of wing structure has been realized [5,6]. The scheme is only suitable for one-dimensional deformation of the beam structure, but it is difficult to estimate the element deformation under multidimensional complex loads because the scheme requires a large number of strain sensors.

In order to meet the requirements of deformation monitoring, the deformation reconstruction algorithm should be able to adapt to complex topology and boundary conditions, and maintain accuracy and stability under a wide range of loads or elastic inertia changes of materials. Tessler and Spangler have proposed an inverse finite element method (iFEM) [7], which is based on the least square variational principle and uses the measured strain data to estimate the deformation of plate structure [8–10]. Furthermore, they have applied it to shape sensing of plate and shell structure undergoing large displacements [11]. Based on inverse finite element theory, Gherlone realized static deformation reconstruction of wing plate structure [12]. Furthermore, the inverse finite beam element, which is suitable for beam/frame structure deformation, is constructed by combining iFEM with Timoshenko beam theory, and the three-dimensional deformation reconstruction of frame structure under static condition is realized [4]. Based on the theory of ZigZag and the idea of inverse finite element, Cerracchio constructed the deformation reconstruction model of plate element for composite sandwich structure [13]. Yong Zhao et al. have studied the effect of sensor distribution on the reconstruction accuracy in the process of inverse finite beam element reconstruction, and realized the 3D deformational reconstruction of the wing like frame structure [14–16]. Xinglin Pan et al. have deeply studied the effect of measurement strain error on reconstruction accuracy in the process of static deformation reconstruction of beam structure, and put forward the error correction using fuzzy network [17]. When the wing is regarded as a constant-section beam structure, the above research studies the deformation reconstruction of the wing structure with high precision [18,19]. For some variable-section wings, when the element division is dense enough, the partial element can be regarded as a constant-section beam, but the number of sensors used is increased, and furthermore, the error caused by element superposition is so big that it cannot be neglected [20]. Guanghong Chuan et.al analyzed the stiffness matrix error of the variable-section Timoshenko Beam for the different rates of the section areas [21]. With the above conclusion, the reconstruction model of the cross-section wing is constructed in this paper. Meanwhile, the above studies only focus on the structural deformation reconstruction under static load, but not on the structural deformation reconstruction under dynamic excitation. Furthermore, it is found that the reconstruction accuracy and stability of inverse finite element method will be affected by the accuracy of the deformation reconstruction model and the sensor placement.

For the sake of sensing the deformation of the cross-section wing accurately, a variable cross-section inverse finite beam element model is proposed in this paper. Meanwhile, this paper derives the relationship between strain and deformation of the cross-section wing based on inverse finite element method. The deformation reconstruction of a cross-section wing is accomplished with the strain measurement of the structure surface. Furthermore, in order to remove the influence of the strain sensors installation and dynamic un-modeled error for the dynamic deformation reconstruction, this paper proposes a self-structuring linear support vector regression algorithm fuzzy network (SSILSVRFN) to modify the strain measurements. The inverse finite element method is employed to reconstruct the deformation displacement with the modified strain, and the accuracy of the deformation reconstruction is improved. The content of this paper is divided into the following parts: firstly, based on the inverse finite element framework, the inverse finite element reconstruction equation for the variable cross-section beam is constructed for the deformation reconstruction of glass fiber wing model with variable cross section. Then, a fuzzy network based on iterative linear support vector regression is constructed to correct the real-time error of dynamic strain measurement. Finally, the vibration experiments on the fiberglass wing model show that the proposed algorithm can effectively improve the accuracy and stability of wing deformation reconstruction.

#### **2. A Deformation–Reconstruction Model for Wing with Variable Cross-Section**

#### *2.1. Inverse Finite Element Model for Variable Cross-Section Beam*

For a typical variable-section beam (Figure 1), the displacement of point B on the beam surface can be expressed with six kinematic variables  *u*, *v*, *w*, θ*x*, θ*y*, θ*<sup>z</sup>* on point A [4]

$$\begin{aligned} u\_x(\mathbf{x}, y, z) &= u(\mathbf{x}) + z \partial\_y(\mathbf{x}) - y \partial\_z(\mathbf{x}) \\ u\_y(\mathbf{x}, y, z) &= v(\mathbf{x}) - z \partial\_x(\mathbf{x}) \\ u\_z(\mathbf{x}, y, z) &= w(\mathbf{x}) + y \partial\_x(\mathbf{x}) \end{aligned} \tag{1}$$

where *ux*, *uy*, and *uz* are the displacements along the *x*, *y*, and *z* axes, respectively, with *u*(*x*), *v*(*x*), and *w*(*x*) denoting the displacements at *y* = *z* = 0; θ*x*(*x*), θ*y*(*x*), and θ*z*(*x*) are the rotations about the three coordinate axes; positive orientations for the displacements and rotations are depicted in Figure 1. These kinematic assumptions neglect the effect of axial warping due to torsion, i.e., each cross-section remains flat and rigid with respect to thickness-stretch deformations along the *y* and *z*-axes. The six kinematic variables can be grouped in vector form as

$$\mathbf{u}(\mathbf{x}) = \begin{bmatrix} \mu\_r \upsilon\_r w\_r \theta\_{\mathbf{x}\_r} \theta\_{\mathbf{y}\_r} \theta\_z \end{bmatrix}^T \tag{2}$$

In the finite element framework, the deformation of any point along the centroid-axis of the beam element can be obtained by interpolation of the determined shape function *N*(*x*) and the node degrees of freedom *u*ε.

$$
u(\mathbf{x}) = \mathcal{N}(\mathbf{x})\mathfrak{u}^{\mathfrak{x}} \tag{3}$$

Based on the small strain hypothesis, the relationship between the section strain at any point along the centroid-axis and the node degree of freedom is

$$e(\mu) = B(\mathbf{x})\mu^{\varepsilon} \tag{4}$$

where the matrix *B*(*x*) (see Appendix A) contains the derivatives of the shape functions *N*(*x*), and the section strain vector *e*(*u*) is expressed as

$$e(\mathbf{u}) \equiv \begin{bmatrix} e\_{1'}e\_{2'}e\_3, e\_4, e\_5, e\_6 \end{bmatrix}^T \tag{5}$$
 
$$e\_1(\mathbf{x}) \equiv \boldsymbol{u}\_{,x}(\mathbf{x}) \quad e\_4(\mathbf{x}) \equiv \boldsymbol{w}\_{,x}(\mathbf{x}) + \boldsymbol{\theta}\_y(\mathbf{x})$$
 
$$e\_2(\mathbf{x}) \equiv \boldsymbol{\theta}\_{y,x}(\mathbf{x}) \quad e\_5(\mathbf{x}) \equiv \boldsymbol{v}\_{,x}(\mathbf{x}) - \boldsymbol{\theta}\_z(\mathbf{x})$$
 
$$e\_3(\mathbf{x}) \equiv -\boldsymbol{\theta}\_{z,x}(\mathbf{x}) \quad e\_6(\mathbf{x}) \equiv \boldsymbol{\theta}\_{x,x}(\mathbf{x})$$

And the surface strain vector [ ε*x*(*x*, *y*, *z*), γ*xz*(*x*, *y*), γ*xy*(*x*, *y*) ] *<sup>T</sup>* of the beam element can be expressed with the above section strain vector *e*(*u*) in theory:

$$\begin{aligned} \varepsilon\_{\mathbf{x}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) &= \varepsilon\_{1}(\mathbf{x}) + \varepsilon \varepsilon\_{2}(\mathbf{x}) + \mathcal{y} \varepsilon\_{3}(\mathbf{x}) \\ \gamma\_{\mathbf{x}\varepsilon}(\mathbf{x}, \mathbf{y}) &= \varepsilon\_{4}(\mathbf{x}) + \mathcal{y} \varepsilon\_{6}(\mathbf{x}) \\ \gamma\_{\mathbf{x}\mathbf{y}}(\mathbf{x}, \mathbf{y}) &= \varepsilon\_{5}(\mathbf{x}) - \varepsilon \varepsilon\_{6}(\mathbf{x}) \end{aligned} \tag{6}$$

While, the section strain vector *e*(*u*) is calculated from the kinematic variables and cannot be directly measured. In the inverse finite element method, the least square error between the section strain *e*(*u*) and the section strain *e*<sup>ε</sup> obtained from the surface measured strain is constructed as following

$$\mathcal{Q}(\mathbf{u}) = \left\| \mathbf{e}(\mathbf{u}) - \mathbf{e}^{\varepsilon} \right\|^2 \tag{7}$$

when ∅(u) obtains the minimum value, the section strain *e*(*u*) is replaced with the section strain *e*<sup>ε</sup> to form a relationship model between the beam element node degree of freedom and the section strain

$$k^{\mu}u^{\mu} = f^{\mu} \tag{8}$$

where *ke* = *<sup>L</sup> n* ∗ *<sup>n</sup> i*=1 *BT*[*xi*]*B*[*xi*] , *f* <sup>ε</sup> = *<sup>L</sup> n* ∗ *<sup>n</sup> i*=1 *BT*[*xi*]*e*ε*<sup>i</sup>* . Note that *<sup>k</sup><sup>e</sup>* resembles an element stiffness matrix of the direct finite element method and *f* <sup>ε</sup> resembles the load vector; *L* is the element length (Figure 1); *n* and *xi* (0 ≤ *xi* ≤ *L*) are, respectively, the number and the axial coordinate of the locations where the section strains are evaluated, and the superscript ε*i* is used to denote the section strains computed from the surface measured strain.

**Figure 1.** Schematic diagram of a variable section beam.

#### *2.2. Calculation of Section Strain of Variable Section Beam Element*

Correctly solving the section strain is the key to the element deformation reconstruction. For the constitutive cross-section beam (Figure 2), the relationship between the section strain *e*<sup>ε</sup> and the external load force and moment can be expressed as [4]

$$\begin{aligned} N &= A\_{\ge} \mathfrak{e}\_1{}^{\varepsilon} & M\_{\ge} &= J\_{\ge} \mathfrak{e}\_6{}^{\varepsilon} \\ Q\_y &= G\_y \mathfrak{e}\_5{}^{\varepsilon} & M\_y &= D\_y \mathfrak{e}\_2{}^{\varepsilon} \\ Q\_{\overline{z}} &= G\_{\overline{z}} \mathfrak{e}\_4{}^{\varepsilon} & M\_{\overline{z}} &= D\_{\overline{z}} \mathfrak{e}\_3{}^{\varepsilon} \end{aligned} \tag{9}$$

where *Ax* <sup>≡</sup> *EA* is the axial rigidity; *Gy* <sup>≡</sup> *<sup>k</sup>*<sup>2</sup> *yGA* and *Gz* <sup>≡</sup> *<sup>k</sup>*<sup>2</sup> *zGA* are the shear rigidities, with *k*<sup>2</sup> *<sup>y</sup>* and *k*<sup>2</sup> *z* denoting the shear correction factors; *Jx* ≡ *GIp* is the torsional rigidity; *Dy* ≡ *EIy* and *Dz* ≡ *EIz* are the bending rigidities; *A* is the section area; *E* (Young's modulus) and *G* (shear modulus) are the elastic constants. *M* is a moment acting in a certain direction on the unit, and its unit is *N*\**m*; the unit of *Q* is *N*.

**Figure 2.** The force form of the beam.

Once the type of the external load is known, the format of the section strain *e*<sup>ε</sup> can be determined with using Equation (9). For example, the concentrated load, bending moment and torsion performed on the beam element, the in-plane force *N*, *Qy*, *Qz* and torsion *Mx* of the element are constant along the axis direction, and the bending moment *My* and *Mz* change along the *x*-axis. The conclusion is valid for both constant-section beam and variable-section beam [21]. The section strain *e*<sup>ε</sup> <sup>1</sup> can be expressed as

$$e\_1^x(\mathbf{x}) = \frac{\mathbb{C}\_1}{A\_\mathbf{x}(\mathbf{x})} \tag{10}$$

where *C*<sup>1</sup> is an unknown constant and *Ax*(*x*) is the tensile stiffness of the element section. Similarly, the section strains *e*<sup>ε</sup> <sup>4</sup>, *<sup>e</sup>*<sup>ε</sup> <sup>5</sup>, and *<sup>e</sup>*<sup>ε</sup> <sup>6</sup> can be written as follows

$$\begin{array}{l} \mathfrak{e}\_{4}^{\,\,\ell}(\mathbf{x}) = \frac{\underset{\mathbf{C}\_{4}}{\mathbf{C}\_{\mathbf{f}\_{2}}(\mathbf{x})}}{\mathbf{C}\_{\mathbf{f}\_{2}}(\mathbf{x})}\\ \mathfrak{e}\_{5}^{\,\,\ell}(\mathbf{x}) = \frac{\underset{\mathbf{C}\_{\mathbf{f}}(\mathbf{x})}{\mathbf{C}\_{\mathbf{f}}(\mathbf{x})}}{\mathbf{C}\_{\mathbf{f}}(\mathbf{x})} \end{array} \tag{11}$$

where *Gy*(*x*), *Gz*(*x*), and *Ip*(*x*) are the bending stiffness and torsional stiffness of the variable section beam element [21]; *C*4, *C*5, and *C*<sup>6</sup> are unknown constants; *C*<sup>4</sup> and *C*<sup>5</sup> can be computed with the following equation in [4]:

$$\begin{array}{l} \mathcal{C}\_{4} = Q\_{\mathcal{Y}} = \frac{\partial D\_{z} \mathcal{C}\_{3}}{\frac{\partial \mathcal{X}}{\partial x}}\\ \mathcal{C}\_{5} = Q\_{z} = \frac{\partial D\_{y} \mathcal{C}\_{2}}{\partial x} \end{array} \tag{12}$$

Since the changes of bending moment *My* and *Mz* along the centroid-axis of the beam element are linear, *My* and *Mz* can be assumed as

$$\begin{aligned} M\_{\mathcal{Y}}(\mathbf{x}) &= \mathbb{C}\_{5}\mathbf{x} + \mathbb{C}\_{2} \\ M\_{z}(\mathbf{x}) &= \mathbb{C}\_{4}\mathbf{x} + \mathbb{C}\_{3} \end{aligned} \tag{13}$$

where *C*<sup>2</sup> and *C*<sup>3</sup> are unknown constants. Submitting the Equation (13) into Equation (9), *e*<sup>ε</sup> <sup>2</sup> and *<sup>e</sup>*<sup>ε</sup> <sup>3</sup> can be expressed as

$$c\_2^{\varepsilon}(\mathbf{x}) = \frac{\mathbb{C}\_5 \mathbf{x} + \mathbb{C}\_2}{D\_{\mathcal{Y}}(\mathbf{x})} \; , \; c\_3^{\varepsilon}(\mathbf{x}) = \frac{\mathbb{C}\_4 \mathbf{x} + \mathbb{C}\_3}{D\_{\mathcal{Z}}(\mathbf{x})} \tag{14}$$

For the variable-section wing (Figure 3), the relationship between the surface strain measurement ε∗ and the section strain vector is ascertained with using the strain-tensor transformation from the (θ, *X*,*r*) to (*X*, *Y*, *Z*) [4,22,23]:

$$
\varepsilon^\* = \varepsilon\_x \left( \cos^2 \beta - \nu \sin^2 \beta \right) + \gamma\_{x0} \cos \beta \sin \beta \tag{15}
$$

**Figure 3.** Orthogonal and cylindrical coordinate systems of the wing section.

Substituting *r* = *Ri* in Equation (6), yields

$$
\varepsilon\_x = \varepsilon\_1 + \varepsilon\_2 R\_i \sin \theta + \varepsilon\_3 R\_i \cos \theta
$$

$$
\gamma\_{x0} = \varepsilon\_4 \cos \theta - \varepsilon\_5 \sin \theta + \varepsilon\_6 R\_i
\tag{16}
$$

**Figure 4.** Location and coordinate system of a FBG sensor placed on the wing external surface.

Substituting (10), (11), (14), and (16) into (15) gives the relationship between the strain measurements and the section strain for the variable section beam element as (Figure 4)

$$\begin{array}{ll} \varepsilon^\*(\mathbf{x}\_i, \boldsymbol{\theta}\_i, \boldsymbol{\beta}\_i) &= \frac{\mathbb{C}\_1}{A\_{\mathbf{z}}(\mathbf{x}\_i)} (\cos^2 \boldsymbol{\beta}\_i - \boldsymbol{\upsilon} \sin^2 \boldsymbol{\beta}\_i) + \frac{\mathbb{C}\_3 \mathbf{x}\_i + \mathbb{C}\_2}{D\_{\mathbf{y}}(\mathbf{x}\_i)} (\cos^2 \boldsymbol{\beta}\_i - \boldsymbol{\upsilon} \sin^2 \boldsymbol{\beta}\_i) \mathbb{R}\_i \sin \boldsymbol{\theta}\_i \\ &+ \frac{\mathbb{C}\_4 \mathbf{x}\_i + \mathbb{C}\_3}{D\_{\mathbf{z}}(\mathbf{x}\_i)} (\cos^2 \boldsymbol{\beta}\_i - \boldsymbol{\upsilon} \sin^2 \boldsymbol{\beta}\_i) \mathbb{R}\_i \cos \boldsymbol{\theta}\_i + \frac{\mathbb{C}\_4}{G\_{\mathbf{z}}(\mathbf{x}\_i)} \cos \boldsymbol{\beta}\_i \sin \boldsymbol{\beta}\_i \cos \boldsymbol{\theta}\_i \\ &- \frac{\mathbb{C}\_5}{G\_{\mathbf{y}}(\mathbf{x}\_i)} \cos \boldsymbol{\beta}\_i \sin \boldsymbol{\beta}\_i \sin \boldsymbol{\theta}\_i + \frac{\mathbb{C}\_6}{I\_{\mathbf{y}}(\mathbf{x}\_i)} R\_i \cos \boldsymbol{\beta}\_i \sin \boldsymbol{\beta}\_i \end{array} \tag{17}$$

where R*i*= <sup>2</sup> *y*2 *<sup>i</sup>* <sup>+</sup> *<sup>z</sup>*<sup>2</sup> *<sup>i</sup>* ; ν is Poisson's ratio; *Ri* is the polar radius of the beam section; *xi*, θ and β are the positions of the strain measurements on the surface of the beam element.

Therefore, the unknown parameters *C*1, *C*2, *C*3, *C*4, *C*5, *C*<sup>6</sup> can be solved from six different surface measurement strain values *xi*, θ*i*, β*i*,(*i* = 1, ... , 6) with using the Equation (17), and the section strains *e*<sup>ε</sup> will be determined with Equations (10), (11), and (14).

#### **3. Strain Error Correction**

Because of the strain measurement system error and the model error, there is a big deviation between the actual displacement and the displacement computed from the strain measurement data with iFEM. The strain measurement system error consist of the location error resulted from the FBG sensors attachment process, measurement error of strain measurement instrument et al. The model error, including the dynamic un-modeled error caused by the difference between the static model and dynamic model, the connection and segmentation of the structure elements. For removing the above errors, the basic modification strategy is: (1) the actual strain value εˆ is computed from actual deformation captured from the third-party measurement instrument with reconstruction model, Equation (8). (2) the fuzzy network of the self-structuring iterative linear support vector regression fuzzy network (SSILSVRFN) algorithm [24] is trained and fixed with the above actual strain value εˆ and the corresponding measured strain values ε. (3) The actual deformation of the structure for any loading cases can be computed from the actual strain value εˆ *<sup>A</sup>* computed from the corresponding strain measurement ε*<sup>A</sup>* with the above fixed fuzzy network.

For the deformation reconstruction based on strain measurement, there are many error factors in the measurement of structural surface strain. Thus, the difference exists between the measured strain and the actual strain; and the corresponding relationship between the measured strain and the actual strain is difficult to be accurately described with mathematical expressions. The self-structuring fuzzy network (SSFN) is attempted to approach this unknown relationship because the fuzzy network has the characteristic of infinite approximation to any function mapping relation. SSFN algorithm constructs the network from zero rule, and the adjustment of the rules number and structure is based on the training data as current network input. The training network is difficult to achieve the desired application effect in the testing phase because network structure is optimal for the current data. This paper proposes a SSILSVRFN using the support vector regression theory and cluster idea. The SSILSVRFN system is divided into the structure training phase and the parameter learning phase.

The structure training phase adopts the SSRG (self-structuring rule generation) algorithm for automatic fuzzy rule generation and initialization. According to the spatial distribution of the input data and the analysis of the overall data, a reasonable distribution of fuzzy sets center and width is achieved. Not only the size of the whole network is effectively reduced, but also the impact on the network structure of the data order is avoided.

The SSILSVRFN proposed in this paper is divided into five layers, as shown in Figure 5.

**Figure 5.** Structure of the self-structuring linear support vector regression algorithm fuzzy network (SSILSVRFN).

In parameter learning phase, the support vector regression (SVR) is based on the structural risk minimization principle; and SVR learning algorithm has a better performance in generalization and prevention of over-fitting. Therefore, SSILSVRFN adopts iterative linear SVR (ILSVR) to iteratively adjust parameter of fuzzy rules.

The SSILSVRFN training and learning system diagram is shown in Figure 6.

**Figure 6.** The SSILSVRFN training.

#### *3.1. SSILSVRFN Structure Learning*

The structure learning determines the fuzzy rule number and initial fuzzy set center *mij* and width σ*ij*. Once the number of fuzzy rules is determined, the numbers of nodes in the layer 2, 3 and 4 will be accordingly determined, which are *nr*, *r* and *r*, respectively. In the structure learning phase, a fuzzy rule is regarded as a cluster that corresponds to a rule node of layer 3 in the input space. A SSRG algorithm is proposed to determine the suitable number of rules.

The flow chart (Figure 7) is as follows

**Figure 7.** The flow chart of self-structuring rule generation SSRG.

The algorithm steps are presented as follows

Step 1: classify the training samples, and set the number of iteration *Nite* = 0. Assign each input training data *x* to cluster P, which is calculated as follows

$$\mathbf{P} = \arg\min\_{1 \le p \le r} d(\mathbf{x}\_{\prime}, m\_p) \tag{18}$$

$$d\left(\mathbf{x}, m\_p\right) = \sum\_{j=1}^{n} \left|\mathbf{x}^j - m\_{pj}\right|\tag{19}$$

Step 2: calculate the maximum variance. For the first cluster, the variance is defined as follows

$$\mathbf{I} = \arg\max\_{1 \le i \le r} \overline{\sigma}\_i^2 \tag{20a}$$

$$
\overline{\sigma}\_i^2 = \sum\_{j=1}^n \overline{\sigma}\_{ij}^2 \tag{20b}
$$

$$\nabla\_{ij}^{2} = \frac{1}{\nabla\_{i}} \sum\_{\mathbf{x} \in \mathbb{C}\_{i}} \left(\mathbf{x}^{j} - m\_{ij}\right)^{2} \tag{20c}$$

Step 3: if the maximum value of cluster variance σ<sup>2</sup> *<sup>I</sup>* is less than the threshold value *vth* and the current iteration number is less than 20, stop; otherwise, the cluster with the largest difference will be split to generate a new cluster. The center point vectors of the new clusters are, respectively, set as *mI* + ε and *mI* − ε, where ε is a constant vector with a small value. This paper sets the *i*th component of ε to be 1% of the domain of input variable *xi* .

Step 4: update the centers of clusters. Recalculate the center of cluster **P** as

$$m\mathbf{p} \leftarrow m\mathbf{p} + \frac{1}{N\_P} \sum\_{k=1}^{N\_P} \left(\mathbf{x}\_{\mathbf{k}} - m\mathbf{p}\right) \tag{21}$$

where *Np* is the number of samples in cluster **P**. The iteration number is updated as *Nite* = *Nite* + 1.

#### *3.2. SSILSVRFN Parameter Learning*

The iterative linear SVR (ILSVR) learning algorithm is used to adjust the antecedent and consequent parameters of fuzzy rules in the SSILSVRFN algorithm. In ILSVR algorithm, the purpose of learning parameters is to optimize antecedent and consequent parameters of fuzzy rules based on the cost function.

#### 3.2.1. Consequent Parameter Learning

In this section, the fourth layer represented by the structure diagram in Figure 5 is used to calculate the output value of each rule. The calculation expression is as follows

$$f^i = \mu^i(\mathbf{x}) \cdot \left(a\_{i0} + \sum\_{j=1}^n a\_{ij} \mathbf{x}^j \right) = \mu^i(\mathbf{x}) \cdot \sum\_{j=0}^n a\_{ij} \mathbf{x}^j, \mathbf{x}^0 \Delta \ \mathbf{1} \tag{22}$$

where *ai*<sup>0</sup> <sup>+</sup> *<sup>n</sup> j*=0 *aijx<sup>j</sup>* is the corresponding consequent value of the node; and μ*<sup>i</sup>* (*x*) is the corresponding firing strength.

The fifth layer represents the output variable of the output layer, which calculates as follows

$$y^\* = \sum\_{i=1}^r f^i + b \tag{23}$$

where *b* is the compensation constant and *f<sup>i</sup>* is the output value of the node.

By substituting (22) into (23), the output of SSILSVRFN can be transformed as follows

$$y^\*(\mathbf{x}) = \sum\_{i=1}^r \sum\_{j=0}^n a\_{ij} \mu^i \mathbf{x}^j + b \tag{24}$$

Equation (24) shows that the output *y*<sup>∗</sup> is a linear function of μ*<sup>i</sup> x<sup>j</sup>* with weights *aij*. Therefore, linear SVR can be employed to learn parameters *aij*. After structure learning, the number of rules *r* and the initial rule antecedent part parameters are determined. The input data *xk* is transformed to the following vector

$$\begin{array}{lcl}\mathfrak{sp}(\mathbf{x\_{k}}) &= \left[\mathfrak{q}^{1}(\mathbf{x\_{k}}) \; , \; \dots \; \; \; \; \mathfrak{q}^{r(n+1)}(\mathbf{x\_{k}})\right] \\ &= \left[\mathfrak{\mu}^{1}\mathbf{x\_{k}}^{0} \; , \dots \; \; \; \; \mu^{1}\mathbf{x\_{k}}^{n} \; , \dots \; \; \; \mu^{r}\mathbf{x\_{k}}^{0} \; , \dots \; \; \; \mu^{r}\mathbf{x\_{k}}^{n}\right] \in \mathfrak{R}^{r(n+1)} \end{array} \tag{25}$$

The vector ϕ as input into a linear SVR, and the training data pairs are represented as follows

$$\begin{array}{lcl} \mathbf{S} &= \left\{ (\mathbf{x}\_1, y\_1), \ (\mathbf{x}\_2, y\_2), \ \dots, \ (\mathbf{x}\_{\mathbf{N}}, y\_{\mathbf{N}}) \ \right\} \\ &= \left\{ \{\boldsymbol{\varphi}(\mathbf{x}\_1), y\_1\}, \ \{\boldsymbol{\varphi}(\mathbf{x}\_2), y\_2\}, \ \dots, \ \{\boldsymbol{\varphi}(\mathbf{x}\_{\mathbf{N}}), y\_{\mathbf{N}}\} \end{array} \tag{26}$$

The linear regression function *y*∗ (*x*) is given as

$$y^\*(\mathbf{x}) = \sum\_{k=1}^N (\alpha\_k - \beta\_k) \mathbf{x}\_k \mathbf{x} + b \tag{27}$$

According to Equations (26) and (27), the optimal linear regression function *y*∗ (*x*) is given as

$$y^\*(\mathbf{x}) = \sum\_{k=1}^N \left( a\_k - \hat{a}\_k \right) \left\langle \ \varphi(\mathbf{x}), \ \ \left. \ \varphi(\mathbf{x}\_k) \right\rangle \right\rangle + b \tag{28}$$

where α*<sup>k</sup>* c and αˆ *<sup>k</sup>* are calculated by SVR software, i.e., Library for Support Vector Machines (LIBSVM) in this section. Based on Equation (25), Equation (28) can be represented as follows

$$\begin{array}{lcl} y^\*(\mathbf{x}) &= \sum\_{k=1}^N \left( \alpha\_k - \hat{\alpha}\_k \right) \sum\_{m=1}^{r(n+1)} q^m(\mathbf{x}) q^m(\mathbf{x}\_k) + b \\ &= \sum\_{m=1}^{r(n+1)} \left[ \sum\_{k=1}^N \left( \alpha\_k - \hat{\alpha}\_k \right) q^m(\mathbf{x}\_k) \right] q^m(\mathbf{x}) + b \\ &= \sum\_{i=1}^r \sum\_{j=0}^n \left[ \sum\_{k=1}^N \left( \alpha\_k - \hat{\alpha}\_k \right) \mu^i \mathbf{x}\_k^j \right] \mu^j \mathbf{x}^j + b \end{array} \tag{29}$$

Since Equation (29) is equivalent to Equation (24), by comparing these two expressions, the mathematical expression of the parameters *aij* can be obtained.

$$\mathbf{a}\_{ij} = \sum\_{k=1}^{N} (\mathbf{a}\_k - \mathbf{\hat{a}}\_k) \mu^i \mathbf{x}\_{\mathbf{k}'}^j \mathbf{i} = 1, \dots, r, j = 0, \dots, n \tag{30}$$

#### 3.2.2. Antecedent Parameter Learning

The initial parameters *mij* and σ*ij* in SSILSVRFN are determined by the SSRG algorithm. Further, these parameters are tuned based on the minimization of the cost function. The output function of SSILSVRFN in Equation (24) can be written as follows

$$\mathbf{y}^\*(\mathbf{x}) = \sum\_{i=1}^r \left\{ \exp \left[ -\sum\_{j=1}^n \left[ \frac{\mathbf{x}^j - m\_{ij}}{a\_{ij}} \right]^2 \right] \right\} \cdot \left( \sum\_{j=0}^n a\_{ij} \mathbf{x}^j \right) + b \tag{31}$$

The above equation shows that the antecedent parameters *mij* and σ*ij* are not the linear combination coefficients of the fuzzy network output *y*∗ . Therefore, the linear SVR cannot be directly applied. In response to this problem, this section uses Taylor series expansion to linearize it. Since the parameters *mij* and σ*ij* are independent of each other, the fuzzy network output *y*<sup>∗</sup> is expanded as shown below.

$$\begin{aligned} y^\*(\mathbf{p}) &= y^\*(\mathbf{p}(0)) + \sum\_{i=1}^r \sum\_{j=1}^n \Delta m\_{ij} (m\_{ij} - m\_{ij}(0)) \\ &+ \sum\_{i=1}^r \sum\_{j=1}^n \Delta \sigma\_{ij} (\sigma\_{ij} - \sigma\_{ij}(0)) + y^\*\_h \end{aligned} \tag{32}$$

where *y*∗ *<sup>h</sup>* is the remainder of the Taylor expansion, <sup>P</sup> = [*m*11, ... , *mrn*, <sup>σ</sup>11, ...σ*rn*] denotes the antecedent parameters vector. The first partial derivative Δ*mij* and Δσ*ij* are, respectively, written as follows

$$
\Delta m\_{i\bar{j}} = \frac{\partial y^\*}{\partial m\_{i\bar{j}}} = f^i \cdot 2 \cdot \frac{x\_{\bar{j}} - m\_{i\bar{j}}}{\left(\sigma\_{i\bar{j}}\right)^2} \tag{33}
$$

$$
\Delta \sigma\_{i\bar{j}} = \frac{\partial y^\*}{\partial \sigma\_{i\bar{j}}} = f^i \cdot 2 \cdot \frac{\left(\mathbf{x}\_{\bar{j}} - m\_{i\bar{j}}\right)^2}{\left(\sigma\_{i\bar{j}}\right)^3} \tag{34}
$$

Let *y*ˆ∗ = *y*∗ (**p**) − *y*<sup>∗</sup> (**p**(0)). Equation (32) can be further transformed as follows

$$\mathcal{G}^\* \approx \sum\_{i=1}^r \sum\_{j=1}^n \left[ \Delta m\_{ij} \left[ m\_{ij} - m\_{ij}(0) \right] + \Delta \sigma\_{ij} \left[ \sigma\_{ij} - \sigma\_{ij}(0) \right] \right] \tag{35}$$

Equation (35) shows that *y*ˆ<sup>∗</sup> is a linear function with linear combination coefficients *mij* − *mij*(0) and σ*ij* − σ*ij*(0). Therefore, linear SVR can be used to optimize tuning parameters *mij* − *mij*(0) and

σ*ij* − σ*ij*(0). Being similar to the derivation of consequent parameter learning, the parameters *mij* and σ*ij* can be solved.

$$\begin{aligned} m\_{i\bar{j}} &= m\_{i\bar{j}}(0) + \sum\_{k=1}^{N} \left[ [a\_k - \hat{\alpha}\_k] \Delta m\_{i\bar{j}}(\mathbf{x}\_k) \right] \\ \sigma\_{i\bar{j}} &= \sigma\_{i\bar{j}}(0) + \sum\_{k=1}^{N} \left[ [a\_k - \hat{\alpha}\_k] \Delta \sigma\_{i\bar{j}}(\mathbf{x}\_k) \right] i = 1, \dots, r, j = 0, \dots, n \end{aligned} \tag{36}$$

The linearized function *y*ˆ∗ in Equation (35) approximates the output of SSILSVRFN *y*∗ with the assumption that the updated values of *mij* and σ*ij* are very close to the original *mij*(0) and σ*ij*(0). After the optimization of the linear SVR, the updated parameters may be too large to meet this assumption. Thus, ILSVR is used to solve this problem. After several linear SVR learning iterations, the parameter learning tends to converge. The updated values of *mij* and σ*ij* can meet the assumption in a Taylor expansion because the change in the two parameters tends to zero.

Throughout the learning process of parameters, the steps of the entire ILSVR algorithm can be summarized as follows

Step 1: set the maximal iteration number *Tite* and set the iteration index *l* to zero.

Step 2: calculate consequent parameters *aij*(*l*), *i* = 1, ... ,*r*, *j* = 0, ... , *n* by using the linear SVR and Equation (30).

Step 3: based on the learned parameters *aij*(*l*) in Step 2, calculate the antecedent parameters *mij*(*l*) and σ*ij*(*l*) using the linear SVR and Equation (36).

Step 4: if *l* ≥ *Tite*, this algorithm stops; otherwise, update the iteration index *l* = *l* + 1 and go back to step 2.

The overall structure block diagram of the variable section wing dynamic deformation reconstruction method is shown in the Figure 8.

**Figure 8.** Structural block diagram of dynamic deformation reconstruction method for variable section wing.

#### **4. Verifications through Simulations and Experimentation**

In order to assess the accuracy and effectiveness of the measurement scheme proposed in this paper, simulations and model experimentation are performed. In Section 4.1, a finite element model of a variable-section wing is constructed through direct finite element (FE) analyses software (ABAQUS), used to assess the accuracy of the reconstructing model of the variable-section wing. In Section 4.2, a physical wing model is tested under different dynamic loads, in order to assess the feasibility of measurement scheme proposed in this paper.

#### *4.1. Simulation Test*

The FE model of the whole variable section wing (Figure 9A) is directly obtained from the CAD model of the wing (Figure 9B) with using ABAQUS. The whole wing is combined the skin with the framework, and the framework mainly consists of glass fiber rib plate and glass fiber sheet. The single wingspan is 1500 mm long. The whole wing structure is divided into two segments for deformation reconstruction: the first section is the wing root part which is 600 mm; the second section is the aileron which is 900 mm. For the largest cross-sectional area of the wing root, the length and the width are 500 mm and 70 mm, respectively. For the smallest cross-sectional area of the wing tip, the length is 200 mm and the width is 28 mm (Figure 10). The cross section area decreases along the span direction, and the shape of each section is the similar (Figure 10B). For the wing skin, the Young's modulus is *E* = 8.3 GPa, the Poisson ratio is *v* = 0.22, and the density is ρ = 4500 kg/m3. For the glass fiber rib plate and the glass fiber sheet, the Young's modulus is *E* = 7.3 GPa, the Poisson ratio is *v* = 0.22, and the density is ρ = 5000 kg/m3. The wing is connected to the airframe with two thin-walled carbon fiber beams. For the carbon fiber beam, the Young's modulus is *E* = 210 GPa, the Poisson ratio is *v* = 0.307, and the density is ρ = 4000 kg/m3. The longer beam length is 757 mm, the external radius is 18 mm and the thickness is 2 mm, and the shorter beam length is 345 mm, the external radius is 7.5 mm and the thickness is 1 mm.

**Figure 9.** The finite element (FE) model and the CAD model of the variable section wing. (**A**) The FE model; (**B**) The CAD model.

**Figure 10.** Wing structure. (**A**) Framework of the wing; (**B**) wing cross section.

Two different loads are performed on the FE model. Herein, the point C (Figure 11) is selected as the target point to assess the accuracy of the reconstructing model. In Figure 11, the hexahedron element is used to divide the wing structure, and the total number of the element is 4546. The comparison between the deformation computed from the strain values with using the reconstruction model (Equation (7)) and the deformation analyzed with using the FE model for the point C is shown in Table 1. The strain values are obtained from the FE analysis.

**Figure 11.** The contour plot of the wing deformation.

The comparisons of Table 1 demonstrate that the reconstructing model presents higher accuracy, for the main deformation *wz* under the two loading cases, the percent errors remain is below 6.0%. Although the percentage errors are relatively high for the displacements in the other two directions and the rotations for three directions, the absolute errors are very small. During the flight, the main deformation of the wing is the displacement along the vertical direction, Z (Figure 11).

**Table 1.** The comparison between the FE analysis and the reconstructing using inverse finite element method (iFEM). The displacements are expressed in millimeter and rotations are expressed in radian.


#### *4.2. Physical Model Test*

For the aim of assessing the effectiveness of the measurement scheme proposed in this paper, a dynamic loading test is performed on the variable-section wing physical model. The size of the physical wing (Figure 12) is same to the CAD model of the wing (Figure 9B).

**Figure 12.** Variable-section wing physical model.

In the experiment, the strain data are obtained from the strain measurement system composed of FBG strain sensors (Fiber Bragg Grating| os1100, Micron Optics, Atlanta, GA, USA) and the FBG interrogator (Optical Sensing Instrument| Si 155, Micron Optics, Atlanta, GA, USA). Twelve FBG strain sensors (the range of initial wavelength is (1527 nm, 1564 nm)) are placed at different locations along the wing and used to capture the surface strain. The placement of the strain sensors is shown in Figure 13.

**Figure 13.** Finite element model of wing with variable cross section.

In order to evaluate the accuracy of the deformation reconstructed from the algorithm proposed in this paper, the actual deformation of the wing structure is captured from the 3D optical measurement instrument (see Figure 14A, NDI Optrotrak Certus) which is abbreviated as NDI. Several position sensors which send the infrared lights to the CCD cameras of the NDI to reflect the deformation of the wing structure. The accuracy of NDI is 0.1 mm in its measurement range. The whole experiment system and the coordinate can be seen in Figure 14B.

**Figure 14.** Wing tests. (**A**) NDI Optrotrak Certus; (**B**) experiment system.

In the experiment, twenty different static loadings are performed on the end of the wing, which used for the fuzzy network training. The total weight of the loading is 20 kg. The dynamic load is caused by the sudden removal of the load applied on the wing tip (Figure 15).

**Figure 15.** Dynamic load performed on the variable section wing model.

Under the working condition, the aircraft is mainly subjected to the vertical lifting force. Therefore, in this experiment, the load along the *Y*-axis is mainly performed on the wing model. When the kinematic variable *u<sup>e</sup>* of the wing is solved by the measured strain through the Equation (8), any deformation of the wing surface can be calculated by using the Equation (1). To verify the accuracy of the reconstruction, the root mean square error RMS (Root Mean Square) is applied.

$$RMS = \sqrt[2]{\sum\_{i=1}^{j} \left(\text{disp}^{\text{NDI}}(\mathbf{x}\_{i}) - \text{disp}^{\text{iFEM}}(\mathbf{x}\_{i})\right)^{2}/j} \tag{37}$$

where, *disp*(*xi*) is the displacement of one node along the wing in one direction; the superscript 'NDI' refers to the deformation values captured from NDI; 'iFEM' refers to the predicted values computed from strain data with iFEM; and *j* is the number of the nodes used to describe the wing deformation. In the test, the number of nodes used to describe the beam deformation is 6, and these nodes are placed on the surface of the wing (Figure 14B). Meanwhile, the node of maximum deformation, point C (Figure 13) is taken account. The tracking of point C for dynamic loading cases can be seen in Figure 16, and the comparison between the deformations of the point C computed from the strain data are shown in Tables 2 and 3. The maximum percentage error among the individual nodal displacements is shown in Table 4. The comparison between the deformations of the whole wing computed from the strain data are shown in Table 5.

**Figure 16.** Comparison between deformation displacements of point C for the different time. (**A**) *X*-axis; (**B**) *Z*-axis; and (**C**) *Y*-axis.

**Table 2.** Comparison of deformation along the y-axis for different time by removing the load of 20 kg loaded at the wing tip, the measured strain is obtained from the fiber Bragg grating (FBG) sensor measurement, and the actual strain is calculated from the NDI measurement with iFEM. *CU <sup>Y</sup>* and *CM Y* are, respectively, the deformation computed from the unmodified and modified strain measurements with iFEM at point C.



**Table 3.** Comparison of deformation along the *x*-axis and *z*-axis for different time by removing the load of 20 kg loaded at the wing tip.

**Table 4.** The maximum percentage error among the individual nodal displacements, *dispU <sup>Y</sup>* and *dispM Y* are, respectively, the displacement computed form the unmodified and modified strain measurements with iFEM.


It is found that, the deformation compute from the modified strain data with iFEM on point C is close to the result captured from the NDI. From the Table 2, it is shown that the minimum value of the percentage error of the C point displacement is 3.9%, and the maximum value is 6.7% along the *y*-axis. Although the percentage errors are relatively high for the displacements in the other two directions, deformations along the *x*-axis and *z*-axis are very small from the Table 3. From the Table 4, it is shown that the maximum percentage error among the individual nodal displacements is 6.7%, and the maximum value is 4.1%.

From the Table 5, it is shown that the accuracy of the deformation reconstruction with iFEM is increased when the strain measurements are modified with using SSILSVRFN algorithm proposed in this paper. The RMS of the reconstruction for the maximum deformation is 12.01 mm when the strain measurements are unmodified, while RMS is 3.97 mm when the strain measurements are modified with using SSILSVRFN algorithm. The reduction percentage of RMS is 66.8% at least.

**Table 5.** Comparison of deformation along the *y*-axis for different time, *RMSU* and *RMSM* are, respectively, the accuracy of the deformation computed form the unmodified and modified strain measurements with iFEM.


During the flight, the main deformation of the wing is the displacement along the vertical direction (Figures 11 and 12). From the Tables 1 and 2, numerical and experimental studies show that:


The dynamic deformation reconstruction algorithm proposed in this paper shows high precision for the wing structure with variable section.

#### **5. Conclusions**

In view of the effects of the dynamic un-modeled error and the sensor placement error on the cross-section wing deformation sensing, a dynamic deformation reconstruction algorithm for wing structure with variable section is proposed in this paper. For the sake of reducing the influence of the strain measurement error, the in situ strain data are modified with using SSILSVRFN in real time. The test performed on the wing shows that, the dynamic deformation of the variable cross-section wing can be accurately reconstructed from the strain measurement with the scheme proposed in this paper. The measurement scheme proposed in this paper can be applied to the dynamic deformation reconstruction of the structure from the variable section beam. If the measurement objects are different, the reconstruction model and the correction network need to be re-established. Nevertheless, the dynamic loading in this paper is caused by sudden removal of the load applied on the wing tip. The further work is performing a dynamic experiment on the shaking platform to simulate the dynamic deformation caused by different airflow excitations, to assess the effectiveness and accuracy of the reconstruction algorithm.

**Author Contributions:** Z.F. carried out the majority of work on the paper, proposing a deformation-reconstruction model for wing with variable cross-section, fuzzy correction method and writing the article; H.B. analyzed the method; Y.Z. and F.Z. gave advice regarding the writing of the article.

**Funding:** This research was funded by [Key Technology Research on Steerable 110-m Aperture Radio Telescope] grant number [2015CB857100], [Xinjiang Astronomical Observatory] grant number [2014KL012], and [Shanghai Aerospace Science and Technology Innovation Fund] grant number [SAST201413].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The interpolation shape function in [4] is only applicable to constant-section beam elements. This paper also uses interdependent interpolations (Tessler and Dong, 1981) to derive the interpolation shape function suitable for the variable section beam element.

$$\begin{aligned} u(\xi) &= \sum\_{i=1,\mathcal{I},2} L\_i^{(2)}(\xi) u\_i \\ v(\xi) &= \sum\_{j=1,\mathcal{J},s,2} L\_j^{(3)}(\xi) v\_i - \sum\_{j=1,\mathcal{J},s,2} N\_j^{(4)} \theta\_{zj} \\ w(\xi) &= \sum\_{j=1,\mathcal{J},s,2} L\_j^{(3)}(\xi) w\_i + \sum\_{j=1,\mathcal{J},s,2} N\_j^{(4)} \theta\_{yj} \\ \theta\_x(\xi) &= \sum\_{i=1,\mathcal{I},2} L\_i^{(2)}(\xi) \theta\_{xi} \\ \theta\_y(\xi) &= \sum\_{j=1,\mathcal{J},s,2} L\_j^{(3)} \theta\_{yj} \\ \theta\_z(\xi) &= \sum\_{j=1,\mathcal{J},s,2} L\_j^{(3)} \theta\_{zj} \end{aligned} \tag{A1}$$

The shape function is expressed as follows

$$\begin{array}{l} [L\_1^{(2)}, L\_1^{(2)}, L\_2^{(2)}] = \frac{1}{2} [\xi(\xi - 1), 2(1 - \xi^2), \xi(\xi + 1)] \\ [L\_1^{(1)}, L\_2^{(2)}] = \frac{1}{16} ((9\xi^2 - 1)) [(1 - \xi), (1 + \xi)] \\ [L\_q^{(3)}, L\_s^{(3)}] = \frac{9}{16} (\xi^2 - 1) [(3\xi - 1), -(3\xi + 1)] \\ [N\_1^{(4)}, N\_2^{(4)}] = \frac{1}{128} [(9\xi^2 - 1)(\xi^2 - 1), -(9\xi^2 - 1)(\xi^2 - 1)] \\ [N\_q^{(4)}, N\_s^{(4)}] = \frac{31\epsilon}{128} [-(9\xi^2 - 1)(\xi^2 - 1), (9\xi^2 - 1)(\xi^2 - 1)] \end{array} \tag{A2}$$

An element is referred to a local axial coordinate *x* ∈ [0, *L*]*,* where *L* denotes the element length. Furthermore, a non-dimensional coordinate ξ ≡ 2*x <sup>L</sup>* − 1 ∈ [−1, 1] is used to define the element shape function (Figure A1).

**Figure A1.** Inverse finite element geometry and nodal topology.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
