Nomenclature

$A$

state matrix

$B$

input matrix

$C$

output matrix

$D$

feedforward matrix

$I$

identity matrix

$K$

constant feedback gain matrix

$\ell $

$n - r \gt 0$

${L_i}$

auxiliary matrix

$m$

dimension of $\boldsymbol{u}$

$n$

dimension of $\boldsymbol{x}$

$r$

dimension of $\boldsymbol{y}$

$t$

time

$u$

input (or control) vector

$U$

auxiliary matrix

$V$

auxiliary matrix

$\boldsymbol{x}$

state vector

$\boldsymbol{y}$

output vector

$\boldsymbol{z}$

transformed state vector

Greek Symbol

${\lambda _i}$

generic eigenvalue

${\boldsymbol{\nu}_i}$

generic eigenvector

Subscripts

a

augmented

Cl

closed loop

ℓ

of dimension equal to ℓ

r

of dimension equal to r

Superscripts

a

achievable

d

desired

†

Moore-Penrose inverse

∼

transformed

.

time derivative

2.0 Problem statement and procedure architecture

A flat truss structure, entirely made of high-stiffness composite material for typical aerospace use, is considered as a test sample and is shown in Fig. 1. Three of the 46 rods are active and equipped with sensors (to monitor the outputs) and actuators (to apply the input forces). Although an in-depth study of the properties of the sensor/actuator system is beyond the scope of this paper, piezoelectric elements of ‘1–3 rod composite’ [Reference Chee, Tong and Steven23] are considered optimal for the present purpose. Indeed, like all piezoelectric devices [Reference Preumont24], they can act simultaneously as sensors and actuators due to their ability to translate mechanical force into electrical energy and vice versa. In addition, their constituent material is a high-stiffness composite, with properties very similar to those of passive rods, so as to authorise their full assimilation with the remaining members when in the passive phase. Actuators, strategically located, are supposed to excite only the structural modes of interest, while sensors are supposed to measure the corresponding modes.

Figure 1. Schematic view of a space truss with three piezoelectric sensors/actuators.

An external random excitation is considered to simulate the action of a thruster, in the hypothesis of measuring its force by sensors.

The various phases of the procedure are accomplished by means of a combined and integrated use of FEM, Python scripting and Fortran programming. Indeed, the Python language is used to develop a unique code, divided into several scripts, both to perform pre- and post-processing tasks within the Abaqus 2022 FE software and to implement system identification and feedback control algorithms; see the block diagram in Fig. 2. In addition, a dedicated subroutine was created in Fortran language to process the control law in FEA. This specific approach is made possible by the well-established FE modeling technique through Python scripting [25], the communication between FE processing and Fortran subroutines [26], the use of specific Python modules for scientific [Reference Virtanen and Gommers27] and array [Reference Harris and Millman28] calculation, and for control system development [Reference Boni, Bassetto, Quarta and Mengali29–Reference Fuller, Greiner and Moore31].

Figure 2. Block diagram of the proposed procedure.

Two different versions of Python are used, namely the older Python 2.7, which is the only one supported by the Abaqus software, and the newer Python 3.11, due to the integration of its libraries with control system techniques. The use of Pickle files enables total communication between scripts written in different Python versions. Indeed, a Pickle file can serialises Python object structures, by converting an object in memory into a stream of bytes to be stored as a binary file on disk. Once re-loaded into a Python program, this binary file can be de-serialised into a Python object.

3.0 System identification

Over the years, many system identification techniques have been developed for space structure applications to cope with their extremely low damping [Reference Lew, Juang and Longman32]. In this context, Markov parameters [Reference Juang, Phan, Horta and Longman33] arise in a realisation theory that finds a state-space representation of a linear time-invariant system from its pulse response. The disadvantage of solving the Markov parameters of the system directly in the time domain from the input and output data is the need to invert an input matrix that becomes very large for slightly damped systems. In this respect, one possible solution to this problem is to obtain, again from input-output data of the system to be identified, the Markov parameters of an asymptotically stable observer. This method, chosen in the present work, is called the OKID [Reference Phan, Horta, Juang and Longman34] and consists of a procedure in which the state-space model and the corresponding observer are determined simultaneously.

Using a Kalman filter as an observer overcomes the problem of signals disturbed by noise. In the matrix formulation developed by Phan et al. [Reference Phan, Horta, Juang and Longman34], the Markov parameters of the system can be uniquely obtained from those of the observer/Kalman filter. Once observer and system Markov parameters are determined, the ERA [Reference Juang and Pappa35] is used to realize the state-space model. Thus, the method referred to as OKID/ERA provides an identification of the undamaged structure (before) and of the damaged structure (after) with a state-space model in the following form

(1) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\dot{\boldsymbol{x}} = A{\rm{\;}}\boldsymbol{x} + B{\rm{\;}}\boldsymbol{u}}\\[5pt]{\boldsymbol{y} = C{\rm{\;}}\boldsymbol{x} + D{\rm{\;}}\boldsymbol{u}}\end{array}} \right.\end{align}

where $\boldsymbol{x} \in {\mathbb{R}^n}$ denotes the system state, $\boldsymbol{u} \in {\mathbb{R}^m}$ is the input (or control) variable, and $\boldsymbol{y} \in {\mathbb{R}^r}$ is the output variable. Note that $\boldsymbol{x}$ , $\boldsymbol{u}$ , and $\boldsymbol{y}$ are all functions of time $t$ . Moreover, without loss of generality, it is assumed that $m \lt r \lt n$ , that is, the number of measured outputs exceeds the number of inputs. The physical nature of the identified state is not known a priori because OKID algorithm selects only the set of matrices $\left\{ {A,{\rm{\;}}B,{\rm{\;}}C,{\rm{\;}}D} \right\}$ that best correlate input with output data. A state transformation can anyway help in giving a clear meaning to $\boldsymbol{x}$ . In particular, if the system order is equal to the number of outputs (i.e. if $r = n$ ), the matrix $C$ is square and can be inverted. Thus, by introducing the transformation

(2) \begin{align}\boldsymbol{z} = C{\rm{\;}}\boldsymbol{x}\end{align}

the following state-space system is obtained

(3) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\dot{\boldsymbol{z}} = \tilde A{\rm{\;}}\boldsymbol{z} + \tilde B\;\boldsymbol{u}}\\[5pt]{\boldsymbol{y} = \tilde C{\rm{\;}}\boldsymbol{z}}\end{array}} \right.\end{align}

where

(4) \begin{align}\tilde A = C{\rm{\;}}A{\rm{\;}}{C^{ - 1}},\quad {\tilde{\rm{B}}} = C{\rm{\;}}B,\quad {\tilde{\rm{C}}} = I\end{align}

According to this transformation, the state variables acquire a clear physical meaning: system outputs and state variables are all coincident with nodal displacements.

However, in the case where the order of the system is greater than that of the measured outputs ( $n \gt r$ ), it turns out that $C$ is rectangular, thus not invertible. To proceed as above, firstly the state vector is partitioned as

\begin{align*}\boldsymbol{x} = \left[ {\begin{array}{*{20}{l}}{{\boldsymbol{x}_r}}\\[2pt]{{\boldsymbol{x}_\ell }}\end{array}} \right]\end{align*}

with ${\boldsymbol{x}_r} \in {\mathbb{R}^r},$ ${\boldsymbol{x}_\ell } \in {\mathbb{R}^\ell }$ and $\ell = n - r$ in such a way that

\begin{align*}\boldsymbol{y} = \left[ {\begin{array}{l@{\quad}l}{{C_1}} {} {}& {{C_2}}\end{array}} \right]\left[ {\begin{array}{*{20}{l}}{{\boldsymbol{x}_r}}\\[2pt]{{\boldsymbol{x}_\ell }}\end{array}} \right]\end{align*}

Then, an augmented output vector is introduced as

(5) \begin{align}{\boldsymbol{y}_a} = \left[ {\begin{array}{*{20}{l}}\boldsymbol{y}\\[2pt]{{\boldsymbol{x}_\ell }}\end{array}} \right] = {C_a}\left[ {\begin{array}{*{20}{l}}{{\boldsymbol{x}_r}}\\[2pt]{{\boldsymbol{x}_\ell }}\end{array}} \right]\end{align}

where

(6) \begin{align}{C_a} = \left[ {\begin{array}{l@{\quad}l}{{C_1}} {}& {{C_2}}\\[3pt]0 {} {}& {{I_\ell }}\end{array}} \right]\end{align}

and ${I_\ell } \in {\mathbb{R}^{\ell \times \ell }}$ is the identity matrix. After this artificial extension of the outputs, the state-space model can be transformed as follows:

(7) \begin{align}\boldsymbol{z} = {C_a}\boldsymbol{x}\end{align}

which now plays the role of Equation (2). Following the above procedure, the system representation becomes

(8) \begin{align}\dot{\boldsymbol{z}} = \tilde A\boldsymbol{x} + \tilde B\boldsymbol{u}\end{align}

(9) \begin{align}\boldsymbol{y} = \tilde C\boldsymbol{z}\end{align}

where

(10) \begin{align}\tilde A = {C_a}AC_a^{ - 1}\end{align}

(11) \begin{align}\tilde B = {C_a}B\end{align}

(12) \begin{align}\tilde C = \left[ {\begin{array}{*{20}{l}}{{I_r}} {}\quad 0\end{array}} \right]\end{align}

The transformation of the state variables of Equation (7) provides clear physical meaning to the identified state space systems because the same type of transformation gives identical meaning to the state variables of the undamaged and damaged structure models, allowing them to be compared. Because of this comparability, system identification techniques enable health monitoring of the structure throughout its life. This can be done by applying independent, time-varying forces through the piezoelectric actuators as well as a disturbance external to the structure. The first identification, in undamaged conditions, represents a reference for those subsequently done: changes in the truss low frequency modes may indicate that damages have crept in the structure and would lead to its replace or refurbish.

4.0 Control law

Once the presence of damage is detected, the problem is how to restore, at least partially, the dynamic behaviour of the structure. In this work, the damaged truss is forced to behave like the intact one by means of a static output feedback controller. Moreover, due to the state transformation in Equation (12), the output feedback behaves as a state feedback designed to match the identifiable eigenvalues and corresponding right eigenvectors of the intact with those of the damaged structure. Accordingly, the control law is in the form

(13) \begin{align}\boldsymbol{u} = K\;\boldsymbol{y}\end{align}

where $K$ is the constant feedback gain matrix and the state matrix in closed loop becomes

(14) \begin{align}{\tilde A_{cl}} = \tilde A + \tilde B{\rm{\;}}K{\tilde{\rm{C}}}\end{align}

As shown by Srinathkumar [Reference Srinathkumar36], ${\rm{max}}\left( {m,r} \right)$ closed loop eigenvalues can be arbitrarily assigned with ${\rm{min}}\left( {m,r} \right)$ components of the corresponding right eigenvectors, provided the system is controllable and observable. The gain matrix can be calculated according to the method developed by Andry et al. [Reference Andry, Shapiro and Chung37], which accounts for the concept of achievable eigenvectors. Indeed, as explained in Ref. (Reference Srinathkumar36), it is not possible to completely assign the right eigenvectors of the controlled system, but it is possible to guarantee that the desired eigenvectors are “best” approximated in a least square sense.

According to this procedure, if ${\lambda _i}$ is a generic eigenvalue and $\boldsymbol{\nu}_i^d$ is the desired right eigenvector, the corresponding achievable eigenvector is

(15) \begin{align}\boldsymbol\nu _i^a = {L_i}{\rm{\;}}L_i^\dagger {\rm{\;}}\boldsymbol\nu _i^d\end{align}

where

(16) \begin{align}{L_i} = {\big( {{\lambda _i}{\rm{\;}}I - \tilde A} \big)^{ - 1}}\tilde{\rm{B}}\end{align}

and $L_i^\dagger $ is the Moore-Penrose inverse of ${L_i}$ . The final expression of the gain matrix discussed in Ref. (Reference Andry, Shapiro and Chung37) is obtained only for a special state transformation addressed to simplify the $B$ matrix structure. In this work, the original procedure has been modified accordingly to be compliant with the state transformation in Equation (3). The obtained expression of the feedback matrix is

(17) \begin{align}K = U{\rm{\;}}{( {\tilde C{\rm{\;}}V} )^{ - 1}} = {\tilde B^{ - 1}}{\rm{\;}}U{\left( V \right)^{ - 1}}\end{align}

where the columns of matrix $V$ are the achievable eigenvectors and the columns of matrix $U$ are related to the achievable eigenvalues and eigenvectors. Note that, as pointed out in Ref. (Reference Andry, Shapiro and Chung37), when the number of poles to be placed is smaller than the dimension of the state matrix the feedback control can adversely affect some of the unconstrained eigenvalues such that system instability can result. In the problem at hand, this phenomenon has been avoided by making $r = n$ with a balanced reduction of the identified systems.

5.0 Finite element models

The FE model of the nine-bay truss shown in Fig. 1 is executed entirely by Python scripting. The truss structure is made of 46 graphite/epoxy rods with hollow circular cross section of mean radius $340{\rm{\;mm}}$ and thickness $2{\rm{\;mm}}$ . The length of each non-diagonal element is $1.5{\rm{\;m}}$ . The material is assumed to be homogeneous with density $1700{\rm{\;kg}}/{{\rm{m}}^2}$ and isotropic, according to the one-dimensional nature of rod elements, with Young’s modulus equal to $90{\rm{\;GPa}}$ , coincident with the E11 value of the actual composite material. No type of damping is assigned to this material. The 2-node linear 2-D truss elements of type T2D2 are used for the rods, while axial connectors elements are adopted to model the piezoelectric sensor/actuators. Each connector, which has the same elastic stiffness as the rods, behaves as a generic truss element if no axial displacements or forces are externally imposed to promote its function as an actuator. The structure is clamped at the left end, while an external excitation normal to the longitudinal axis of the truss is applied to the right upper end node. A Gaussian white noise is added to the excitation to simulate the presence of an external disturbance. Three actuators are supposed to be located in the horizontal, vertical and diagonal direction as shown in Fig. 1. The location of the actuators is based upon simplified controllability considerations (singular values analysis).

The FE model of the undamaged truss originates three different types of analysis: (i) to calculate the lowest three natural frequencies and modes of vibrations of the structure; (ii) to enable system identification; (iii) and to record the response to external excitation. The first two analyses are grouped and cascaded as a frequency step followed by one of dynamic implicit type. In the last analysis, a single dynamic implicit step is processed. In each case, the simulation lasts 5 seconds and is based on the implicit operator defined by Hilber et al. [Reference Hilber and Hughes38] for the temporal integration of the dynamic problem. The set of input data used to identify the system is shown in Fig. 3: axial forces are applied to the actuators and a transversal force at the upper left node of the truss, all following a time history based on white noise. The white noise amplitude is generated with zero mean value, unitary standard deviation, and $2 \times {10^3}$ frequency. The vertical component of displacement and velocity are recorded at the ends of the actuators and at the upper right node as output data, as shown in Fig. 3. The system to be identified has 4 inputs and 14 outputs. The last analysis is addressed to record the vertical displacement of the left upper node when the external transversal force with double-step amplitude shown in Fig. 4 acts on it.

Figure 3. Scheme of the inputs for system identification.

Figure 4. Time history of the external force.

The FE model of the damaged structure is obtained by simply removing some passive members: in the example shown in Fig. 5, two horizontal rods are removed from the top of the truss. All FEAs performed on the undamaged truss are repeated for the damaged structure. A final analysis is then added to simulate the control law developed to restore the dynamic behaviour of the truss. This analysis exploits the possibility of defining a load history through an external subroutine written in Fortran.

Figure 5. Introduction of damage as member removal.

In addition, the Abaqus software allows sensors to be placed at strategic nodal points that can transmit the desired output variables to the subroutine at each time increment. Once the control gain matrix is defined according to the procedure described in Section 4, its coefficients are copied into the subroutine and used to define the loads on the actuators as a linear combination of their end displacements, continuously updated by the sensors. In the latter analysis, a dynamic implicit step is implemented, in which the external load and the control forces act simultaneously.

6.0 Numerican results

6.1 System identification

The Python script that implements the OKID/ERA algorithm, according to the procedure described in Section 3, provides the four matrices $\left\{ {A,{\rm{\;}}B,{\rm{\;}}C,{\rm{\;}}D} \right\}$ necessary to define the discrete space-state model of the truss. The identified 14th-order model with a time interval of 0.005 seconds is assembled and then transformed into continuous form in a dedicated Python script. The response of the system is simulated to both input loads used for its identification and to the external double step force. All these results are obtained for both undamaged and damaged structure. In view of the control strategy development, for the damaged truss, a balanced reduction is applied to the system, in order to decrease its order from 14 to 6, which is double the number of complex poles to be allocated.

The proposed methodology depends heavily on the effectiveness of the OKID/ERA algorithm in identifying a dynamic system with realistic behaviour. Because of the importance of an experimental reference, the results provided by FEA are used, in this paper, as reference values. As an example, for the undamaged structure, Fig. 6 shows the displacement in vertical direction of the upper right node of the truss in case of white noise input used for identification. A perfect superposition between the two data sets is clearly visible.

Figure 6. Comparison between FEM and predicted outputs for undamaged truss with white noise input.

For the damaged truss, the FE results are compared with the responses predicted by both the full and reduced order system. As an example, Fig. 7 shows the displacement in vertical direction of the upper left node of the truss for the case of external double step force. The graph in figure highlights the excellent ability of the identified system, even in its reduced form, to capture the structural response.

Figure 7. Comparison between FEM and predicted outputs for damaged truss with external double step force.

In Fig. 8, the eigenvalues extracted from FEAs are compared with the 14 of the identified systems for the undamaged truss case. The frequency step implemented in Abaqus allows to calculate the eigenvalues of the symmetrised system, that is, neglecting damping or any source of asymmetry in the stiffness matrix. The symmetrised system has real squared eigenvalues ${\lambda ^2}$ and real eigenvectors only: if the stiffness matrix is positive semidefinite there are only pairs of conjugate imaginary eigenvalues $ \pm \lambda i$ , where $\lambda $ is the circular frequency [39]. The eigenvalues identified by FEAs are accordingly aligned on the imaginary axis, symmetrically distributed with respect to the real axis. In contrast, the eigenvalues of the identified system depend on the input/output data of the dynamic implicit step and can capture any source of damping.

Figure 8. Comparison between FEM and predicted eigenvalues for undamaged truss.

Despite the absence of intrinsic material damping, the operator adopted to integrate the dynamic equations is characterised by user-defined numerical damping, which varies from a minimum in transient fidelity analyses to a maximum in steady-state ones. The nature of artificial numerical damping is to damp out unwanted high-frequency modes to reduce the simulation time up to steady state conditions. When performing transient fidelity analysis, as in the present study, lower frequencies are not altered by artificial damping, while the higher the frequencies, the more they are damped. This phenomenon is evident in Fig. 8, where the predicted poles lie on a pseudo-ellipse arc tangent to the imaginary axis: the three lowest frequencies of FEAs, related to the poles highlighted in red, are totally coincident with the predicted ones. In Fig. 9, a similar comparison between the predicted and FE eigenvalues is made for the undamaged truss case. The predicted poles, compared with the previous case, are slightly more damped, consistently with the dynamic behaviour of artificial damping, which varies according to the convergence performance of the specific FE model.

Figure 9. Comparison between FEM and predicted eigenvalues for damaged truss.

6.2 Control strategy

In the control strategy outlined in Section 4, the application of an output-feedback control law was simulated in such a way that the first three complex eigenvalues and the components of the right eigenvectors of the healthy structure can be assigned to the damaged structure to restore, at least in part, the dynamic behaviour of the truss.

In Fig. 10, it is possible to directly compare the first eigenvalues of the two systems and observe the effect of the controller. As is evident from the figure, only the imaginary part of the poles was considered in the control strategy. The vertical displacement of the upper right node due to the double step force is shown in Fig. 11 for the undamaged, damaged and controlled structure. From the figure, a substantial overlap between the structural response of the undamaged truss and the controlled truss can be seen: the controller, designed to simply recover the healthy eigenstructure, also shows good trajectory tracking capability in the presence of an external force.

Figure 10. Comparison between eigenvalues of undamaged and damaged truss.

Figure 11. Comparison between structural response of undamaged, damaged and controlled truss.

The effects of control on the structure is finally studied by simulating it in a dedicated FEA, as explained in Section 5. The overlay plot tool, available in Abaqus, allows to superimpose the deformation of the whole truss in undamaged and damaged-controlled conditions. In particular, Fig. 12 shows the displacement field at an instant of maximum vibration amplitude. From this figure, it is possible to detect the substantial effectiveness of the controller in recovering the healthy structural behaviour along the whole truss axis and not only at the point of external force application.

Figure 12. Overlay plots of truss deformations in undamaged and damaged-controlled conditions.

Another important aspect is the possibility of designing the controller not only to recover but also to improve the dynamic response of the structure. In fact, it is possible to add damping to the vibrations simply by placing poles with negative real part. For example, looking at Fig. 13, the control restores frequencies to their original values and, at the same time, adds damping by placing poles with negative real part of 0.25.

Figure 13. Comparison between structural response of undamaged, damaged and controlled truss with added damping.