# Adaptive Model Predictive Vibration Control of a Cantilever Beam with Real-Time Parameter Estimation

Authors: Gergely Takács; Tomáš Polóni; Boris Rohal’-Ilkiv
Journal: Shock and Vibration
 (2014)

---

## Abstract

This paper presents an adaptive-predictive vibration control system using extended Kalman filtering for the joint estimation of system states and model parameters. A fixed-free cantilever beam equipped with piezoceramic actuators serves as a test platform to validate the proposed control strategy. Deflection readings taken at the end of the beam have been used to reconstruct the position and velocity information for a second-order state-space model. In addition to the states, the dynamic system has been augmented by the unknown model parameters: stiffness, damping constant, and a voltage/force conversion constant, characterizing the actuating effect of the piezoceramic transducers. The states and parameters of this augmented system have been estimated in real time, using the hybrid extended Kalman filter. The estimated model parameters have been applied to define the continuous state-space model of the vibrating system, which in turn is discretized for the predictive controller. The model predictive control algorithm generates state predictions and dual-mode quadratic cost prediction matrices based on the updated discrete state-space models. The resulting cost function is then minimized using quadratic programming to find the sequence of optimal but constrained control inputs. The proposed active vibration control system is implemented and evaluated experimentally to investigate the viability of the control method.


---

## Body

## 1. Introduction

Undesirable mechanical and structural vibrations may often cause discomfort in humans and in certain engineering applications can even lead to catastrophic failure or other extreme consequences. Passive vibration attenuation approaches are popular in engineering practice, but the required structural changes tend to get troublesome with low frequency vibrations [1, 2]. With the advent of new actuator and sensor types and the availability of cheap computing technology, active vibration control (AVC) has become an important tool in managing excessive vibration levels [3, 4].When designing the algorithm support for AVC systems, a frequent assumption is that the controlled structure maintains its dynamic properties throughout the control procedure. This assumption enables the straightforward tuning of some controllers used in vibration attenuation like positive position feedback (PPF) [2, 4], while, in model-based algorithms such as linear quadratic (LQ) [3, 5] or model predictive control [6–8] (MPC), it allows the use of a relatively precise nominal model. However, not every AVC application fits this convenient premise. A time-varying system behavior may simply detune the controller causing it to operate with suboptimal performance, but it may also affect the stability of the closed-loop control system. Also, one of the most important properties of self-reliant structural control systems is adaptivity, implying a degree of in situ intelligence [9].One of the possible ways to handle model and parameter changes is robust controller design. The basis of the design approach is the a priori analysis of system behavior given a compact set of uncertain parameters. Within this predefined bounded set, the controller remains stable and fulfills certain tradeoff performance criteria; nevertheless, the best performance is usually achieved only with the nominal model. Robust versions of numerous well-known control methods have been considered for the control of vibrating mechanical structures, including robust pole-placement [10], robust H
∞ [11], robust LQ [12], and robust min–max LQ [1]. Nonetheless, a real adaptive vibration control system offering a degree of self-reliance requires advanced control looking beyond the boundaries of mere robustness [9].In addition to robust controller design practices, the problem of variable structural properties is often addressed by introducing a level of adaptivity into the control loop. Adaptivity enables the vibration control system to function under varying load conditions, partial actuator failure, or performance degradation, a limited change in geometric configuration or varying physical parameters. Modeling uncertainties—and to a degree nonlinearities—can be also compensated by adaptivity [13]. Thus, the advantage of adaptive algorithms in general is their superior performance, although they tend to be computationally more expensive than robust implementations. An excellent overview of adaptive structural control with a focus on self-reliance is given by Hyland and Davis in [9].Possibly one of the best known and most widely used adaptive methods applied for AVC and especially active noise control (ANC) is the least mean square (LMS) feedforward control [14–16], where the changes in system dynamics are translated into the coefficient variations in the adaptive finite impulse response filter (FIR). A frequent choice of adaptive AVC encountered in literature is model reference adaptive control (MRAC) [13], which does not require the use of explicitly defined system parameters [17]. Other adaptive vibration control schemes proposed in the past include sliding mode control [18, 19], neural networks [20, 21], or an adaptive control scheme featuring a high-gain observer for compensating model errors, which was introduced for vibrating microcantilevers by Zhang et al. [22]. Suitable online identification methods for adaptive AVC include the subspace identification or autoregressive moving average (ARMAX) models, adaptive feedforward control, while variants of the extended Kalman filter (EKF) have been used for real-time modal parameter diagnostics applications [23].In order to formulate an adaptive vibration controller, instead of identifying the coefficients in FIR filters or using the models without explicit parameterization, we apply the augmented continuous-discrete extended Kalman filter [24–26] to estimate the dynamic states as physically interpretable parameters of a vibrating mechanical system, such as the stiffness and damping coefficient. The extended Kalman filter is possibly the most widely used algorithm for the estimation of states and parameters of nonlinear dynamic systems [24]; here, it performs the online parameter identification for adaptive-predictive vibration control. The use of EKF or its alternatives to identify the parameters of vibrating systems is also common; however, most of these methods involve offline parameter estimation from measured data or online diagnostics and do not use EKF in real-time to supplement model data to model-based control systems. Examples include disturbance force estimation on a free-free beam [27], health monitoring of base-isolated structures under seismic loads [28], and stiffness and damping coefficient estimation in hydromount systems for reducing vibrations in automotive applications [29]. Recently, Szabat and Orlowska-Kowalska used the extended Kalman filter to adaptively control torsional vibration in a two-mass rotating drive system [30, 31]. The EKF was previously used as a benchmark algorithm for the open-loop state and parameter estimation of the free vibrations of the cantilever beam considered in this work, contrasting the results to moving horizon observer (MHO) in an offline simulation [32].In this work, the parameters estimated by the EKF are used to update the continuous model, which is in turn discretized, to create the state predictions for an input constrained infinite horizon dual-mode MPC algorithm [33, 34]. Regarded by many as one of the most important developments in control engineering [35], the key advantage of using a model-based predictive algorithm in any system is that it can take the effect of future state, input, or output constraints into account when computing the control move [36]. Constraints are essential in every application field, but piezoceramics often used in vibration control are particularly prone to failure and performance degradation by depolarization, which may occur if voltage limits are continually exceeded [8, 37]. The explicit constraint handling feature offered by the MPC method is unique and valuable property amongst the available control algorithms [34]. In addition to this hallmark feature, the performance of a well-designed and tuned MPC algorithm may surpass other traditional control strategies used in vibration control. This is due to the fact that model predictive control is an optimal control method. A significant drawback of using advanced optimization-based control methods such as the MPC is the computational cost of solving the problem online [6, 37].The combination of EKF and MPC resulting in the proposed adaptive EKF-MPC vibration control scheme is applied to a clamped aluminum cantilever beam which is equipped with piezoceramic transducers and its deformation is measured at the free end using a laser triangulation system (Figure1). Let us assume that the dynamics of this beam can be modeled using a second-order differential equation describing the motion of the point-mass-damper system. This assumption holds true, if the first resonance of the structure dominates the dynamic response, which is the case for a class of real-life structures [8, 38]. It is necessary to use such a simple dynamic model in order to ensure the real-time feasibility of the computation intensive quadratic programming solver required for constrained MPC. Furthermore, let us assume that the force contribution from the actuators can be modeled by the input voltage multiplied by a scalar conversion constant and that the beam is excited by an external shock-like disturbance force.Clamped cantilever with active vibration control.
(a)
Active beam assembly
(b)
Piezoceramic actuators
(c)
Laser feedback and stingerThe disturbance force is a series of impacts delivered to the beam, emulating the release test similar to transient vibration typical for aerospace structures [39]. The aim of the adaptive EKF-MPC vibration control algorithm then is to attempt to keep the beam position near its equilibrium using constrained inputs and despite the external disturbances, while continually adapting to possible changes in the dynamic behavior by observing the unknown system parameters in addition to the system states. The control system will shorten the overall settling time of the structure and concurrently adapt the algorithm properties to parameter changes.The vibration control algorithm proposed here unifies the advantages of the adaptivity enabled by the online extended Kalman filter and constrained model predictive control. The EKF-MPC method offers greater damping performance in the face of time-varying physical parameters and model uncertainties, while at the same time ensuring the fulfillment of process constraints in order to increase the safety and reliability of modern AVC systems. The advantages introduced above come at a price, which is mainly the increase in complexity and computational load.After the formerly presented introductory section, the theoretical methodology of the work is discussed in Section2. First, the modeling of the dynamics of the vibrating cantilever including the parameter augmentation is presented; then, we review the formulation of the extended Kalman filter for the given problem. The discretization procedure of the continuous system model and the constrained model predictive algorithm is followed by an overall summary of the proposed EKF-MPC vibration control algorithm. Section 3 presents the experimental hardware in more detail along with the common assumptions and utilized settings, while the different experimental scenarios are also introduced in this section. We report the results of the experiments and discuss their meaning in Section 4 and, finally, the paper is concluded in Section 5.


## 2. Methodology

### 2.1. System Dynamics Augmented by Unknown Parameters

Let us assume that the vibration dynamics of the fixed-free cantilever beam can be approximated by a second-order differential equation describing the motion of a spring-mass-damper with an external force input [3, 4] and an unknown external disturbance
(1)
m
q
¨
(
t
)
+
b
q
˙
(
t
)
+
k
q
(
t
)
=
c
u
(
t
)
+
F
e
(
t
)
,
where m (kg) is the equivalent mass of the beam and the actuators, b (Ns/m) is the equivalent viscous damping coefficient, and k (N/m) is the equivalent stiffness; let c (N/V) represent the scalar coefficient transforming input voltage to force F
(
t
)
=
c
u
(
t
) (N) and let F
e
(
t
) (N) be the disturbance [40]. The input to this system is the voltage supplied to the actuators u
(
t
) (V) and the single measured output is the position of the beam q
(
t
) (m).Such a second-order differential equation representing a 1 DOF vibrating system only takes the first resonant frequency into account; however, from the viewpoint of the control system, this is not an issue, since the first resonance dominates the response. Single degree of freedom models are routinely used in literature to identify the physical parameters of vibrating structures using the EKF and other statistical estimation methods [32, 41, 42].If we choose the position and velocity as state variables (x
1
(
t
)
=
q
(
t
)
,
x
2
(
t
)
=
q
˙
(
t
)), it is possible to express this equation in a continuous, linear, state-space form
(2)
x
˙
(
t
)
=
Φ
x
(
t
)
+
Γ
u
(
t
)
+
Θ
F
e
(
t
)
,
(3)
y
(
t
)
=
C
x
(
t
)
.MatrixΦ
∈
ℝ
2
×
2 is the dynamics matrix, Γ
∈
ℝ
2
×
1 is the input matrix, and Θ
∈
ℝ
2
×
1 is the unknown force transition matrix, which for the assumed model take the following form [40, 43]:
(4)
Φ
=
[
0
1
-
k
m
-
b
m
]
Γ
=
[
0
c
m
]
Θ
=
[
0
1
m
]
.Furthermore,C
=
[
1
0
] outputs the position y
(
t
)
=
q
(
t
) (m), measured at a single point at the beam tip, while direct input feedthrough is omitted in this case. Assuming that the equivalent mass m is known and constant but other parameters may vary with time, the system dynamics noted as Φ
(
t
) and Γ
(
t
) are time variable and are continually updated in the EKF-MPC algorithm.The continuous statex
(
t
) describing only the system dynamics as defined in (2) can be augmented by the vector of unknown parameters p
(
t
). In this case,
(5)
p
(
t
)
=
[
k
(
t
)
b
(
t
)
c
(
t
)
F
e
(
t
)
]
T
,
where p
(
t
) is expressing the unknown and potentially changing damping, stiffness, force conversion parameter, and the external disturbance. The new augmented state becomes
(6)
x
a
(
t
)
=
[
x
(
t
)
p
(
t
)
]
T
.The system dynamics, which is augmented by the unknown parameters, can be expressed as(7)
x
˙
(
t
)
=
f
~
c
(
x
(
t
)
,
p
(
t
)
,
u
(
t
)
)
+
w
s
(
t
)
,
p
˙
(
t
)
=
w
p
(
t
)
,
where w
s
(
t
) is the system process noise and w
p
(
t
) is the parameter process noise. Equation (7) can be combined to get the augmented dynamics
(8)
x
˙
a
(
t
)
=
f
c
(
x
a
(
t
)
,
u
(
t
)
)
+
w
(
t
)
,
(9)
y
(
t
)
=
h
(
x
a
(
t
)
,
u
(
t
)
)
+
v
(
t
)
=
[
1
0
0
0
0
0
]
x
a
(
t
)
+
v
(
t
)
,
where the function f
c represents the continuous augmented dynamics, h
c is the continuous measurement function, and w
(
t
)
=
[
w
s
(
t
)
w
p
(
t
)
]
T expresses system and parameter process noise. The process noise and measurement noise have the properties of white noise; in other words, they have a sequentially uncorrelated Gaussian distribution with a zero mean:
(10)
w
(
t
)
~
N
(
0
,
Q
f
t
)
,
v
(
t
)
~
N
(
0
,
R
f
)
,
where Q
f
t is a process noise covariance matrix and R
f is a measurement noise covariance matrix.Note that, by augmenting the dynamics of the mass-spring-damper system by the parameters in (8) and (9), we essentially turn the otherwise linear state-space model into a nonlinear formulation, since the parameters will be used in multiplicative operations inside the model [24, 25].Even though the real system dynamics is continuous in nature, in practice the output can be only observed at discrete sampling timest
=
k
T, k
=
1,2
,
3
,
…, where T is the sampling interval. Therefore, we will consider a numerically simulated nonlinear system [24] described and propagated by
(11)
x
a
(
k
+
1
)
=
f
(
x
a
(
k
)
,
u
(
k
)
)
+
w
(
k
)
,
y
(
k
)
=
h
(
x
a
(
k
)
,
u
(
k
)
)
+
v
(
k
)
=
[
1
0
0
0
0
0
]
x
a
(
k
)
+
v
(
k
)
.


### 2.2. Hybrid Extended Kalman Filter

The EKF algorithm is commonly used to estimate the parameters of nonlinear dynamic systems and the variant summarized here is known in literature as the continuous-discrete or hybrid EKF [24]. The EKF [24, 25] is initialized with the initial estimate of the augmented state
(12)
x
^
a
0
+
=
E
[
x
a
0
]
and the covariance matrix of the initial state error estimate
(13)
P
0
+
=
E
[
(
x
a
0
-
x
^
a
0
+
)
(
x
a
0
-
x
^
a
0
+
)
T
]
.To obtain the a priori state estimate (denoted by the- subscript) in the continuous part of the filter, the state estimate of the augmented dynamic system (8) is simulatively propagated from time t
=
(
k
-
1
) one step ahead to the next time instant t
=
(
k
)
(14)
x
^
a
(
k
)
-
=
f
(
x
^
a
(
k
-
1
)
+
,
u
(
k
)
)
,using the system dynamics of the mass-spring-damper augmented by the change in parametersp
(
t
)
=
[
k
(
t
)
b
(
t
)
c
(
t
)
F
e
(
t
)
]
T
(15)
q
˙
(
t
)
=
x
˙
1
(
t
)
=
x
2
(
t
)
q
¨
(
t
)
=
x
˙
2
(
t
)
=
-
1
m
x
3
(
t
)
x
1
(
t
)
-
1
m
x
4
(
t
)
x
2
(
t
)
+
1
m
x
5
(
t
)
u
(
t
)
+
1
m
x
6
(
t
)
k
˙
(
t
)
=
x
˙
3
(
t
)
=
0
b
˙
(
t
)
=
x
˙
4
(
t
)
=
0
c
˙
(
t
)
=
x
˙
5
(
t
)
=
0
F
˙
e
(
t
)
=
x
˙
6
(
t
)
=
0
.The time update of the a priori covariance matrix estimate [24, 25] is given by the equation
(16)
P
˙
(
t
)
=
Z
(
x
^
a
(
t
)
)
P
(
t
)
+
P
(
t
)
Z
T
(
x
^
a
(
t
)
)
+
Q
f
t
,
where the Jacobian of the state equation (15) is expressed by(17)
Z
(
x
^
a
)
=
∂
f
c
(
x
a
)
∂
x
a
|
x
a
=
x
^
a
=
[
0
1
0
0
0
0
-
1
m
x
3
(
t
)
-
1
m
x
4
(
t
)
-
1
m
x
1
(
t
)
-
1
m
x
2
(
t
)
1
m
u
(
t
)
1
m
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
]
x
a
=
x
^
a
.Concluding the continuous part of the hybrid EKF, the covariance matrix estimate is obtained by simulative propagation of (16) from time t
=
(
k
-
1
) to the sample at t
=
(
k
)
(18)
P
(
k
)
-
=
g
(
P
(
k
-
1
)
+
,
Z
(
x
^
a
(
k
-
1
)
+
)
)
.The Kalman gain matrix at sample (k) is
(19)
K
f
(
k
)
=
P
(
k
)
-
L
T
[
L
P
(
k
)
-
L
T
+
M
R
f
M
T
]
-
1
.The augmented state estimate is updated to obtain the a posteriori estimate (denoted by the+ subscript) along with the a posteriori update of the error covariance matrix
(20)
x
^
a
(
k
)
+
=
x
^
a
(
k
)
-
+
K
f
(
k
)
[
y
(
k
)
-
h
(
x
^
a
(
k
)
-
,
u
(
k
)
)
]
,
(21)
P
(
k
)
+
=
[
I
-
K
f
(
k
)
L
]
P
(
k
)
-
[
I
-
K
f
(
k
)
L
]
T
+
K
f
(
k
)
M
R
f
M
T
K
f
(
k
)
T
,
using the Jacobians of the measurement equation h
(
x
a
(
t
)
) with respect to x
a
(
t
) and v
a
(
t
)
(22)
L
(
k
)
=
∂
h
(
x
a
(
t
)
)
∂
x
a
(
t
)
|
x
a
(
t
)
=
x
^
a
(
t
)
-
=
[
1
0
0
0
0
0
]
M
(
k
)
=
∂
h
(
x
a
(
t
)
)
∂
v
a
(
t
)
|
x
a
(
t
)
=
x
^
a
(
t
)
-
=
1
.Bothx
^
a
(
k
)
+ and P
(
k
)
+ are stored and used in the next step of the a priori computation.


### 2.3. Model Updates and Discretization

The continuous model in (2) defined by Φ
(
k
) and Γ
(
k
) is updated at discrete times (
k
) based on the latest available estimates of the parameters k
(
k
), b
(
k
), and c
(
k
) by the EKF introduced in the previous section. The linear time-invariant MPC formulation assumed in this work expects a discrete-time model given by
(23)
x
(
k
+
1
)
=
A
(
k
)
x
(
k
)
+
B
(
k
)
u
(
k
)
,
(24)
y
(
k
)
=
C
x
(
k
)
.MatrixA
(
k
)
∈
ℝ
2
×
2 is the dynamics matrix, B
(
k
)
∈
ℝ
2
×
1 is the input matrix, and C
=
[
1
0
] outputs the position y
(
k
)
=
x
1
(
k
)
=
q
(
k
) (m). Using a discrete sampling time T (s), the continuous-time state and input matrices Φ
(
k
) and Γ
(
k
) are discretized [44] using
(25)
A
(
k
)
=
I
+
Φ
(
k
)
T
Ψ
(
k
)
,
(26)
B
(
k
)
=
Ψ
(
k
)
T
Γ
(
k
)
,
where the term Ψ
(
k
) is expressed by the infinite series
(27)
Ψ
(
k
)
=
I
+
Φ
(
k
)
T
2
!
+
Φ
(
k
)
T
3
!
+
Φ
(
k
)
T
4
!
+
⋯
,
which in practice can be approximated using successive iterations j of [44]
(28)
Ψ
(
k
)
≈
I
+
Φ
(
k
)
T
2
(
I
+
Φ
(
k
)
T
3
(
⋯
Φ
(
k
)
T
j
-
1
(
I
+
Φ
(
k
)
T
j
)
)
⋯
)
.


### 2.4. Model Predictive Control

The online MPC optimization problem involves minimizing the cost function subject to constraints on the voltage input. The particular MPC method considered here is known as the infinite horizon constrained dual-mode MPC algorithm [33, 34]. After choosing a horizon of n steps, we may predict the evolution of the future states—compactly denoted by x
→
(
k
)—based on the estimate of the current state x
0
(
k
)
=
x
^
0
(
k
)
+ (excluding parameters p
(
k
)) and the sequence of future inputs u
→
(
k
), at time (
k
) by the recursive substitution of states using (23) [34]. We may express this in a compact notation by [33, 34]
(29)
x
→
(
k
)
=
M
(
k
)
x
0
(
k
)
+
N
(
k
)
u
→
(
k
)
,
where matrices M
(
k
) and N
(
k
) for step (
k
) are [33](30)
M
(
k
)
=
[
A
(
k
)
0
A
(
k
)
1
…
A
(
k
)
n
-
2
A
(
k
)
n
-
1
A
(
k
)
n
]
T
(31)
N
(
k
)
=
[
0
0
…
0
0
0
B
(
k
)
0
…
0
0
0
A
(
k
)
B
(
k
)
B
(
k
)
…
0
0
0
⋮
⋮
⋱
⋮
⋮
⋮
A
(
k
)
n
-
1
B
(
k
)
A
(
k
)
n
-
2
B
(
k
)
…
A
(
k
)
2
B
(
k
)
A
(
k
)
B
(
k
)
B
(
k
)
]
.Let us define the linear, quadratic cost function, expressing the contribution of future statesx
→
(
k
) and future inputs u
→
(
k
) into the numerical control quality indicator J(-). This cost function uses the dual-mode paradigm, utilizing free inputs for the first n steps and a fixed feedback controller afterwards [45]:
(32)
J
(
k
)
=
∑
i
=
0
n
-
1
‍
(
x
(
k
+
i
)
T
Q
x
(
k
+
i
)
+
u
(
k
+
i
)
T
R
u
(
k
+
i
)
)
+
x
(
k
+
n
)
T
P
f
(
k
)
x
(
k
+
n
)
,
where Q
=
Q
T
≥
0 is the user determined penalization matrix for states and R
=
R
T
≥
0 is the penalization matrix for inputs. Furthermore, P
f
(
k
) is the solution of the unconstrained, infinite horizon quadratic regulation problem [33, 34, 45] at sample (
k
). The recursive feasibility of constraints to provide stability guarantees [46, 47] is not implemented in this work, so that the feasibility of the real-time implementation is ensured given the relatively short sampling periods of the vibration control problem [8, 37].Unlike in the case of MPC with a nominal model, here the terminal weightP
f
(
k
) must be recomputed online as the solution of the discrete-time algebraic Riccati equation (DARE) [48]:
(33)
A
(
k
)
T
P
f
(
k
)
A
(
k
)
-
P
f
(
k
)
-
A
(
k
)
T
P
f
(
k
)
B
(
k
)
(
R
+
B
(
k
)
T
P
f
(
k
)
B
(
k
)
)
-
1
B
(
k
)
P
f
(
k
)
A
(
k
)
+
Q
=
0
.For the numerical solution of this classical control problem, amongst others, iterative cost reduction algorithms [49], Schur methods [50], or eigenvalue decomposition methods [48, 51] may be applied, of which we consider the latter procedure based on the work of Pappas et al. [48]. Let us define G
(
k
) as
(34)
G
(
k
)
=
B
(
k
)
R
-
1
B
(
k
)
T
.Consider the generalized eigenvalue problem as(35)
Y
(
k
)
V
(
k
)
=
U
(
k
)
X
(
k
)
V
(
k
)
,
where U
(
k
) is the diagonal matrix containing generalized eigenvalues and the full matrix V
(
k
) corresponding to the generalized principal vectors of (35). We define X
(
k
) and Y
(
k
) as
(36)
X
(
k
)
=
[
I
G
(
k
)
0
A
(
k
)
T
]
Y
(
k
)
=
[
A
(
k
)
0
-
Q
I
]
.After solving (35), the diagonal matrix U
(
k
) can be used to identify the indices of stable eigenvalues. The eigenvectors of V
(
k
) that are corresponding to the identified stable eigenvalues can be inserted in a matrix W
(
k
) used as the basis of the stable eigenspace, which is partitioned into two submatrices
(37)
W
(
k
)
=
[
W
1
(
k
)
W
2
(
k
)
]
.Finally, the solution of the discrete-time algebraic Riccati equation in (33) can be computed as [48]
(38)
P
f
(
k
)
=
W
2
(
k
)
W
1
(
k
)
-
1
.Since the adaptive MPC algorithm considered here does not implement a recursive constraint checking horizon needed to guarantee stability on an infinite horizon [46, 47], it does not explicitly need the LQ gain itself, only the solution of the DARE. Nevertheless, this work uses the unconstrained equivalent LQ gain as a certain diagnostic indicator of the algorithm behavior. The unconstrained equivalent LQ gain can be computed by evaluating
(39)
K
(
k
)
=
(
R
+
B
(
k
)
T
P
f
(
k
)
B
(
k
)
)
-
1
B
(
k
)
T
P
f
(
k
)
A
(
k
)
.The model predictive control algorithm is formulated based on the minimization of the dual-mode cost in (32). To find the solution of the MPC problem at a given sample, one must execute the minimization [33, 34]
(40)
min
⁡
u
k
J
(
u
k
,
x
k
)
=
∑
i
=
0
n
-
1
‍
(
x
(
k
+
i
)
T
Q
x
(
k
+
i
)
+
u
(
k
+
i
)
T
R
u
(
k
+
i
)
)
∑
i
=
0
n
-
1
‍
+
x
(
k
+
n
)
T
P
f
(
k
)
x
(
k
+
n
)
,
subject to the following constraints:
(41)
u
_
≤
u
(
k
+
i
)
≤
u
¯
,
i
=
0,1
,
2
,
…
,
n
-
1
(42)
x
0
(
k
)
=
x
(
k
)
(43)
x
(
k
+
1
+
i
)
=
A
(
k
)
x
(
k
+
i
)
+
B
(
k
)
u
(
k
+
i
)
,
i
≥
0
(44)
y
(
k
+
i
)
=
C
x
(
k
+
i
)
,
i
≥
0
(45)
u
(
k
+
i
)
=
K
(
k
)
x
(
k
+
i
)
,
i
≥
n
,
where (41) expresses the predefined process constraints on the voltage levels, (42) is the observed dynamic state containing position and velocity estimates, (43) and (44) are the discrete linear state-space model updated by the parameters, and (45) is the constraint which is implicitly included in the terminal cost and, due to the receding horizon implementation, it is never realized.In order to make this optimization procedure suitable for numerical solvers, the cost and the constraints are transformed. The cost function in (32) can be compactly denoted [33, 34] by
(46)
J
(
k
)
=
u
→
(
k
)
T
H
(
k
)
u
→
(
k
)
+
2
x
0
(
k
)
T
G
(
k
)
T
u
→
(
k
)
+
x
0
(
k
)
T
F
(
k
)
x
0
(
k
)
.Since the partx
0
(
k
)
T
F
(
k
)
x
0
(
k
) does not depend on the optimization variable u
→
(
k
) and contributes only a fixed value to the cost at each iteration, we may omit it from the optimization procedure. For the variable model structure, H
(
k
) and G
(
k
) are evaluated online at each step (
k
) by [33, 34]
(47)
H
(
k
)
=
∑
i
=
0
n
-
1
‍
N
i
(
k
)
T
Q
N
i
(
k
)
+
N
n
(
k
)
T
P
f
(
k
)
N
n
(
k
)
+
ℛ
,
G
(
k
)
=
∑
i
=
0
n
-
1
‍
N
i
(
k
)
T
Q
M
i
(
k
)
+
N
n
(
k
)
T
P
f
(
k
)
M
n
(
k
)
,
where the index i denotes the ith block row, respectively, and n denotes the last block row of N
(
k
) and M
(
k
) identified, discretized, and computed for time (
k
). Matrix ℛ is the block matrix with the input penalty R on its main diagonal. Using this formulation, the constraints defined by (43)–(45) are included in the cost function itself.The upper and lower constraints in (41) which are enforced on all future input variables u
_
≤
u
→
(
k
)
<
u
¯ until the end of the horizon can be compactly expressed [34] by
(48)
[
I
-
I
]
u
→
(
k
)
≤
[
1
u
¯
-
1
u
_
]
.In this case, the constraints define limits on the allowable voltage potential supplied to the actuators.


### 2.5.  Adaptive-Predictive EKF-MPC Vibration Control Strategy

The resulting adaptive-predictive EKF-MPC control strategy can be summarized in the following algorithm [40].Algorithm 1.
Perform the following set of operations at each sampling instant(
k
): (1)
propagate the statex
^
a
(
k
)
- in (14) and covariance matrix P
˙
(
t
) in (16) in simulation to obtain the a priori estimate;
(2)
sample the actual deflectiony
(
k
) filtered by the low-pass and running mean filter;
(3)
compute the Kalman gainK
f
(
k
) using (19);
(4)
based on the measurement sample and gain, update the statex
^
a
(
k
)
+ through (20) to get a posteriori state estimates; then update the covariance matrix P
(
k
)
+ in (21);
(5)
based on the new parametersp
(
k
), reassemble the continuous system model Φ
(
k
) and Γ
(
k
) by (4);
(6)
approximateΨ
(
k
) through (28) to discretize the system matrices A
(
k
) and B
(
k
) using (25);
(7)
solve the discrete-time algebraic Riccati equation in (33) for the terminal weighting P
f
(
k
) matrix using (34)–(38);
(8)
use the discretized model to compute the state prediction matricesM
(
k
) and N
(
k
) by (30);
(9)
use the state prediction matrices to compute the cost prediction matricesH
(
k
) and G
(
k
) through (46);
(10)
minimize the quadratic MPC cost functionJ
(
k
) in (46) subject to the following input constraints: u
¯
≤
u
(
k
+
i
)
≤
u
_, i
=
0,1
,
2
,
…
,
n
-
1;
(11)
apply the first element of the vector of optimal control movesu
→
(
k
) to the controlled system;
(12)
repeat the procedure from(
1
).


## 2.1. System Dynamics Augmented by Unknown Parameters

Let us assume that the vibration dynamics of the fixed-free cantilever beam can be approximated by a second-order differential equation describing the motion of a spring-mass-damper with an external force input [3, 4] and an unknown external disturbance
(1)
m
q
¨
(
t
)
+
b
q
˙
(
t
)
+
k
q
(
t
)
=
c
u
(
t
)
+
F
e
(
t
)
,
where m (kg) is the equivalent mass of the beam and the actuators, b (Ns/m) is the equivalent viscous damping coefficient, and k (N/m) is the equivalent stiffness; let c (N/V) represent the scalar coefficient transforming input voltage to force F
(
t
)
=
c
u
(
t
) (N) and let F
e
(
t
) (N) be the disturbance [40]. The input to this system is the voltage supplied to the actuators u
(
t
) (V) and the single measured output is the position of the beam q
(
t
) (m).Such a second-order differential equation representing a 1 DOF vibrating system only takes the first resonant frequency into account; however, from the viewpoint of the control system, this is not an issue, since the first resonance dominates the response. Single degree of freedom models are routinely used in literature to identify the physical parameters of vibrating structures using the EKF and other statistical estimation methods [32, 41, 42].If we choose the position and velocity as state variables (x
1
(
t
)
=
q
(
t
)
,
x
2
(
t
)
=
q
˙
(
t
)), it is possible to express this equation in a continuous, linear, state-space form
(2)
x
˙
(
t
)
=
Φ
x
(
t
)
+
Γ
u
(
t
)
+
Θ
F
e
(
t
)
,
(3)
y
(
t
)
=
C
x
(
t
)
.MatrixΦ
∈
ℝ
2
×
2 is the dynamics matrix, Γ
∈
ℝ
2
×
1 is the input matrix, and Θ
∈
ℝ
2
×
1 is the unknown force transition matrix, which for the assumed model take the following form [40, 43]:
(4)
Φ
=
[
0
1
-
k
m
-
b
m
]
Γ
=
[
0
c
m
]
Θ
=
[
0
1
m
]
.Furthermore,C
=
[
1
0
] outputs the position y
(
t
)
=
q
(
t
) (m), measured at a single point at the beam tip, while direct input feedthrough is omitted in this case. Assuming that the equivalent mass m is known and constant but other parameters may vary with time, the system dynamics noted as Φ
(
t
) and Γ
(
t
) are time variable and are continually updated in the EKF-MPC algorithm.The continuous statex
(
t
) describing only the system dynamics as defined in (2) can be augmented by the vector of unknown parameters p
(
t
). In this case,
(5)
p
(
t
)
=
[
k
(
t
)
b
(
t
)
c
(
t
)
F
e
(
t
)
]
T
,
where p
(
t
) is expressing the unknown and potentially changing damping, stiffness, force conversion parameter, and the external disturbance. The new augmented state becomes
(6)
x
a
(
t
)
=
[
x
(
t
)
p
(
t
)
]
T
.The system dynamics, which is augmented by the unknown parameters, can be expressed as(7)
x
˙
(
t
)
=
f
~
c
(
x
(
t
)
,
p
(
t
)
,
u
(
t
)
)
+
w
s
(
t
)
,
p
˙
(
t
)
=
w
p
(
t
)
,
where w
s
(
t
) is the system process noise and w
p
(
t
) is the parameter process noise. Equation (7) can be combined to get the augmented dynamics
(8)
x
˙
a
(
t
)
=
f
c
(
x
a
(
t
)
,
u
(
t
)
)
+
w
(
t
)
,
(9)
y
(
t
)
=
h
(
x
a
(
t
)
,
u
(
t
)
)
+
v
(
t
)
=
[
1
0
0
0
0
0
]
x
a
(
t
)
+
v
(
t
)
,
where the function f
c represents the continuous augmented dynamics, h
c is the continuous measurement function, and w
(
t
)
=
[
w
s
(
t
)
w
p
(
t
)
]
T expresses system and parameter process noise. The process noise and measurement noise have the properties of white noise; in other words, they have a sequentially uncorrelated Gaussian distribution with a zero mean:
(10)
w
(
t
)
~
N
(
0
,
Q
f
t
)
,
v
(
t
)
~
N
(
0
,
R
f
)
,
where Q
f
t is a process noise covariance matrix and R
f is a measurement noise covariance matrix.Note that, by augmenting the dynamics of the mass-spring-damper system by the parameters in (8) and (9), we essentially turn the otherwise linear state-space model into a nonlinear formulation, since the parameters will be used in multiplicative operations inside the model [24, 25].Even though the real system dynamics is continuous in nature, in practice the output can be only observed at discrete sampling timest
=
k
T, k
=
1,2
,
3
,
…, where T is the sampling interval. Therefore, we will consider a numerically simulated nonlinear system [24] described and propagated by
(11)
x
a
(
k
+
1
)
=
f
(
x
a
(
k
)
,
u
(
k
)
)
+
w
(
k
)
,
y
(
k
)
=
h
(
x
a
(
k
)
,
u
(
k
)
)
+
v
(
k
)
=
[
1
0
0
0
0
0
]
x
a
(
k
)
+
v
(
k
)
.


## 2.2. Hybrid Extended Kalman Filter

The EKF algorithm is commonly used to estimate the parameters of nonlinear dynamic systems and the variant summarized here is known in literature as the continuous-discrete or hybrid EKF [24]. The EKF [24, 25] is initialized with the initial estimate of the augmented state
(12)
x
^
a
0
+
=
E
[
x
a
0
]
and the covariance matrix of the initial state error estimate
(13)
P
0
+
=
E
[
(
x
a
0
-
x
^
a
0
+
)
(
x
a
0
-
x
^
a
0
+
)
T
]
.To obtain the a priori state estimate (denoted by the- subscript) in the continuous part of the filter, the state estimate of the augmented dynamic system (8) is simulatively propagated from time t
=
(
k
-
1
) one step ahead to the next time instant t
=
(
k
)
(14)
x
^
a
(
k
)
-
=
f
(
x
^
a
(
k
-
1
)
+
,
u
(
k
)
)
,using the system dynamics of the mass-spring-damper augmented by the change in parametersp
(
t
)
=
[
k
(
t
)
b
(
t
)
c
(
t
)
F
e
(
t
)
]
T
(15)
q
˙
(
t
)
=
x
˙
1
(
t
)
=
x
2
(
t
)
q
¨
(
t
)
=
x
˙
2
(
t
)
=
-
1
m
x
3
(
t
)
x
1
(
t
)
-
1
m
x
4
(
t
)
x
2
(
t
)
+
1
m
x
5
(
t
)
u
(
t
)
+
1
m
x
6
(
t
)
k
˙
(
t
)
=
x
˙
3
(
t
)
=
0
b
˙
(
t
)
=
x
˙
4
(
t
)
=
0
c
˙
(
t
)
=
x
˙
5
(
t
)
=
0
F
˙
e
(
t
)
=
x
˙
6
(
t
)
=
0
.The time update of the a priori covariance matrix estimate [24, 25] is given by the equation
(16)
P
˙
(
t
)
=
Z
(
x
^
a
(
t
)
)
P
(
t
)
+
P
(
t
)
Z
T
(
x
^
a
(
t
)
)
+
Q
f
t
,
where the Jacobian of the state equation (15) is expressed by(17)
Z
(
x
^
a
)
=
∂
f
c
(
x
a
)
∂
x
a
|
x
a
=
x
^
a
=
[
0
1
0
0
0
0
-
1
m
x
3
(
t
)
-
1
m
x
4
(
t
)
-
1
m
x
1
(
t
)
-
1
m
x
2
(
t
)
1
m
u
(
t
)
1
m
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
]
x
a
=
x
^
a
.Concluding the continuous part of the hybrid EKF, the covariance matrix estimate is obtained by simulative propagation of (16) from time t
=
(
k
-
1
) to the sample at t
=
(
k
)
(18)
P
(
k
)
-
=
g
(
P
(
k
-
1
)
+
,
Z
(
x
^
a
(
k
-
1
)
+
)
)
.The Kalman gain matrix at sample (k) is
(19)
K
f
(
k
)
=
P
(
k
)
-
L
T
[
L
P
(
k
)
-
L
T
+
M
R
f
M
T
]
-
1
.The augmented state estimate is updated to obtain the a posteriori estimate (denoted by the+ subscript) along with the a posteriori update of the error covariance matrix
(20)
x
^
a
(
k
)
+
=
x
^
a
(
k
)
-
+
K
f
(
k
)
[
y
(
k
)
-
h
(
x
^
a
(
k
)
-
,
u
(
k
)
)
]
,
(21)
P
(
k
)
+
=
[
I
-
K
f
(
k
)
L
]
P
(
k
)
-
[
I
-
K
f
(
k
)
L
]
T
+
K
f
(
k
)
M
R
f
M
T
K
f
(
k
)
T
,
using the Jacobians of the measurement equation h
(
x
a
(
t
)
) with respect to x
a
(
t
) and v
a
(
t
)
(22)
L
(
k
)
=
∂
h
(
x
a
(
t
)
)
∂
x
a
(
t
)
|
x
a
(
t
)
=
x
^
a
(
t
)
-
=
[
1
0
0
0
0
0
]
M
(
k
)
=
∂
h
(
x
a
(
t
)
)
∂
v
a
(
t
)
|
x
a
(
t
)
=
x
^
a
(
t
)
-
=
1
.Bothx
^
a
(
k
)
+ and P
(
k
)
+ are stored and used in the next step of the a priori computation.


## 2.3. Model Updates and Discretization

The continuous model in (2) defined by Φ
(
k
) and Γ
(
k
) is updated at discrete times (
k
) based on the latest available estimates of the parameters k
(
k
), b
(
k
), and c
(
k
) by the EKF introduced in the previous section. The linear time-invariant MPC formulation assumed in this work expects a discrete-time model given by
(23)
x
(
k
+
1
)
=
A
(
k
)
x
(
k
)
+
B
(
k
)
u
(
k
)
,
(24)
y
(
k
)
=
C
x
(
k
)
.MatrixA
(
k
)
∈
ℝ
2
×
2 is the dynamics matrix, B
(
k
)
∈
ℝ
2
×
1 is the input matrix, and C
=
[
1
0
] outputs the position y
(
k
)
=
x
1
(
k
)
=
q
(
k
) (m). Using a discrete sampling time T (s), the continuous-time state and input matrices Φ
(
k
) and Γ
(
k
) are discretized [44] using
(25)
A
(
k
)
=
I
+
Φ
(
k
)
T
Ψ
(
k
)
,
(26)
B
(
k
)
=
Ψ
(
k
)
T
Γ
(
k
)
,
where the term Ψ
(
k
) is expressed by the infinite series
(27)
Ψ
(
k
)
=
I
+
Φ
(
k
)
T
2
!
+
Φ
(
k
)
T
3
!
+
Φ
(
k
)
T
4
!
+
⋯
,
which in practice can be approximated using successive iterations j of [44]
(28)
Ψ
(
k
)
≈
I
+
Φ
(
k
)
T
2
(
I
+
Φ
(
k
)
T
3
(
⋯
Φ
(
k
)
T
j
-
1
(
I
+
Φ
(
k
)
T
j
)
)
⋯
)
.


## 2.4. Model Predictive Control

The online MPC optimization problem involves minimizing the cost function subject to constraints on the voltage input. The particular MPC method considered here is known as the infinite horizon constrained dual-mode MPC algorithm [33, 34]. After choosing a horizon of n steps, we may predict the evolution of the future states—compactly denoted by x
→
(
k
)—based on the estimate of the current state x
0
(
k
)
=
x
^
0
(
k
)
+ (excluding parameters p
(
k
)) and the sequence of future inputs u
→
(
k
), at time (
k
) by the recursive substitution of states using (23) [34]. We may express this in a compact notation by [33, 34]
(29)
x
→
(
k
)
=
M
(
k
)
x
0
(
k
)
+
N
(
k
)
u
→
(
k
)
,
where matrices M
(
k
) and N
(
k
) for step (
k
) are [33](30)
M
(
k
)
=
[
A
(
k
)
0
A
(
k
)
1
…
A
(
k
)
n
-
2
A
(
k
)
n
-
1
A
(
k
)
n
]
T
(31)
N
(
k
)
=
[
0
0
…
0
0
0
B
(
k
)
0
…
0
0
0
A
(
k
)
B
(
k
)
B
(
k
)
…
0
0
0
⋮
⋮
⋱
⋮
⋮
⋮
A
(
k
)
n
-
1
B
(
k
)
A
(
k
)
n
-
2
B
(
k
)
…
A
(
k
)
2
B
(
k
)
A
(
k
)
B
(
k
)
B
(
k
)
]
.Let us define the linear, quadratic cost function, expressing the contribution of future statesx
→
(
k
) and future inputs u
→
(
k
) into the numerical control quality indicator J(-). This cost function uses the dual-mode paradigm, utilizing free inputs for the first n steps and a fixed feedback controller afterwards [45]:
(32)
J
(
k
)
=
∑
i
=
0
n
-
1
‍
(
x
(
k
+
i
)
T
Q
x
(
k
+
i
)
+
u
(
k
+
i
)
T
R
u
(
k
+
i
)
)
+
x
(
k
+
n
)
T
P
f
(
k
)
x
(
k
+
n
)
,
where Q
=
Q
T
≥
0 is the user determined penalization matrix for states and R
=
R
T
≥
0 is the penalization matrix for inputs. Furthermore, P
f
(
k
) is the solution of the unconstrained, infinite horizon quadratic regulation problem [33, 34, 45] at sample (
k
). The recursive feasibility of constraints to provide stability guarantees [46, 47] is not implemented in this work, so that the feasibility of the real-time implementation is ensured given the relatively short sampling periods of the vibration control problem [8, 37].Unlike in the case of MPC with a nominal model, here the terminal weightP
f
(
k
) must be recomputed online as the solution of the discrete-time algebraic Riccati equation (DARE) [48]:
(33)
A
(
k
)
T
P
f
(
k
)
A
(
k
)
-
P
f
(
k
)
-
A
(
k
)
T
P
f
(
k
)
B
(
k
)
(
R
+
B
(
k
)
T
P
f
(
k
)
B
(
k
)
)
-
1
B
(
k
)
P
f
(
k
)
A
(
k
)
+
Q
=
0
.For the numerical solution of this classical control problem, amongst others, iterative cost reduction algorithms [49], Schur methods [50], or eigenvalue decomposition methods [48, 51] may be applied, of which we consider the latter procedure based on the work of Pappas et al. [48]. Let us define G
(
k
) as
(34)
G
(
k
)
=
B
(
k
)
R
-
1
B
(
k
)
T
.Consider the generalized eigenvalue problem as(35)
Y
(
k
)
V
(
k
)
=
U
(
k
)
X
(
k
)
V
(
k
)
,
where U
(
k
) is the diagonal matrix containing generalized eigenvalues and the full matrix V
(
k
) corresponding to the generalized principal vectors of (35). We define X
(
k
) and Y
(
k
) as
(36)
X
(
k
)
=
[
I
G
(
k
)
0
A
(
k
)
T
]
Y
(
k
)
=
[
A
(
k
)
0
-
Q
I
]
.After solving (35), the diagonal matrix U
(
k
) can be used to identify the indices of stable eigenvalues. The eigenvectors of V
(
k
) that are corresponding to the identified stable eigenvalues can be inserted in a matrix W
(
k
) used as the basis of the stable eigenspace, which is partitioned into two submatrices
(37)
W
(
k
)
=
[
W
1
(
k
)
W
2
(
k
)
]
.Finally, the solution of the discrete-time algebraic Riccati equation in (33) can be computed as [48]
(38)
P
f
(
k
)
=
W
2
(
k
)
W
1
(
k
)
-
1
.Since the adaptive MPC algorithm considered here does not implement a recursive constraint checking horizon needed to guarantee stability on an infinite horizon [46, 47], it does not explicitly need the LQ gain itself, only the solution of the DARE. Nevertheless, this work uses the unconstrained equivalent LQ gain as a certain diagnostic indicator of the algorithm behavior. The unconstrained equivalent LQ gain can be computed by evaluating
(39)
K
(
k
)
=
(
R
+
B
(
k
)
T
P
f
(
k
)
B
(
k
)
)
-
1
B
(
k
)
T
P
f
(
k
)
A
(
k
)
.The model predictive control algorithm is formulated based on the minimization of the dual-mode cost in (32). To find the solution of the MPC problem at a given sample, one must execute the minimization [33, 34]
(40)
min
⁡
u
k
J
(
u
k
,
x
k
)
=
∑
i
=
0
n
-
1
‍
(
x
(
k
+
i
)
T
Q
x
(
k
+
i
)
+
u
(
k
+
i
)
T
R
u
(
k
+
i
)
)
∑
i
=
0
n
-
1
‍
+
x
(
k
+
n
)
T
P
f
(
k
)
x
(
k
+
n
)
,
subject to the following constraints:
(41)
u
_
≤
u
(
k
+
i
)
≤
u
¯
,
i
=
0,1
,
2
,
…
,
n
-
1
(42)
x
0
(
k
)
=
x
(
k
)
(43)
x
(
k
+
1
+
i
)
=
A
(
k
)
x
(
k
+
i
)
+
B
(
k
)
u
(
k
+
i
)
,
i
≥
0
(44)
y
(
k
+
i
)
=
C
x
(
k
+
i
)
,
i
≥
0
(45)
u
(
k
+
i
)
=
K
(
k
)
x
(
k
+
i
)
,
i
≥
n
,
where (41) expresses the predefined process constraints on the voltage levels, (42) is the observed dynamic state containing position and velocity estimates, (43) and (44) are the discrete linear state-space model updated by the parameters, and (45) is the constraint which is implicitly included in the terminal cost and, due to the receding horizon implementation, it is never realized.In order to make this optimization procedure suitable for numerical solvers, the cost and the constraints are transformed. The cost function in (32) can be compactly denoted [33, 34] by
(46)
J
(
k
)
=
u
→
(
k
)
T
H
(
k
)
u
→
(
k
)
+
2
x
0
(
k
)
T
G
(
k
)
T
u
→
(
k
)
+
x
0
(
k
)
T
F
(
k
)
x
0
(
k
)
.Since the partx
0
(
k
)
T
F
(
k
)
x
0
(
k
) does not depend on the optimization variable u
→
(
k
) and contributes only a fixed value to the cost at each iteration, we may omit it from the optimization procedure. For the variable model structure, H
(
k
) and G
(
k
) are evaluated online at each step (
k
) by [33, 34]
(47)
H
(
k
)
=
∑
i
=
0
n
-
1
‍
N
i
(
k
)
T
Q
N
i
(
k
)
+
N
n
(
k
)
T
P
f
(
k
)
N
n
(
k
)
+
ℛ
,
G
(
k
)
=
∑
i
=
0
n
-
1
‍
N
i
(
k
)
T
Q
M
i
(
k
)
+
N
n
(
k
)
T
P
f
(
k
)
M
n
(
k
)
,
where the index i denotes the ith block row, respectively, and n denotes the last block row of N
(
k
) and M
(
k
) identified, discretized, and computed for time (
k
). Matrix ℛ is the block matrix with the input penalty R on its main diagonal. Using this formulation, the constraints defined by (43)–(45) are included in the cost function itself.The upper and lower constraints in (41) which are enforced on all future input variables u
_
≤
u
→
(
k
)
<
u
¯ until the end of the horizon can be compactly expressed [34] by
(48)
[
I
-
I
]
u
→
(
k
)
≤
[
1
u
¯
-
1
u
_
]
.In this case, the constraints define limits on the allowable voltage potential supplied to the actuators.


## 2.5.  Adaptive-Predictive EKF-MPC Vibration Control Strategy

The resulting adaptive-predictive EKF-MPC control strategy can be summarized in the following algorithm [40].Algorithm 1.
Perform the following set of operations at each sampling instant(
k
): (1)
propagate the statex
^
a
(
k
)
- in (14) and covariance matrix P
˙
(
t
) in (16) in simulation to obtain the a priori estimate;
(2)
sample the actual deflectiony
(
k
) filtered by the low-pass and running mean filter;
(3)
compute the Kalman gainK
f
(
k
) using (19);
(4)
based on the measurement sample and gain, update the statex
^
a
(
k
)
+ through (20) to get a posteriori state estimates; then update the covariance matrix P
(
k
)
+ in (21);
(5)
based on the new parametersp
(
k
), reassemble the continuous system model Φ
(
k
) and Γ
(
k
) by (4);
(6)
approximateΨ
(
k
) through (28) to discretize the system matrices A
(
k
) and B
(
k
) using (25);
(7)
solve the discrete-time algebraic Riccati equation in (33) for the terminal weighting P
f
(
k
) matrix using (34)–(38);
(8)
use the discretized model to compute the state prediction matricesM
(
k
) and N
(
k
) by (30);
(9)
use the state prediction matrices to compute the cost prediction matricesH
(
k
) and G
(
k
) through (46);
(10)
minimize the quadratic MPC cost functionJ
(
k
) in (46) subject to the following input constraints: u
¯
≤
u
(
k
+
i
)
≤
u
_, i
=
0,1
,
2
,
…
,
n
-
1;
(11)
apply the first element of the vector of optimal control movesu
→
(
k
) to the controlled system;
(12)
repeat the procedure from(
1
).


## 3. Laboratory Hardware and Experimental Settings

The experimental hardware (Figure1(a)) consists of a cantilever beam fixed at one end and free at the other. Though it is very simple in its design, this lightly damped cantilever beam used here as the demonstration example may represent a wide variety of real-life structures [8, 38] from the standpoint of dynamic control, including wing surfaces [39], helicopter rotor blades [52], robotic manipulators in space [53], solar panels in space [54], and antenna systems [55].The beam is made of commercially available pure aluminum in the dimensions of550
×
40
×
3mm. The actuating elements are MIDÉ QP16n piezoelectric transducers made of PZT5A piezoceramic material (Figure 1(b)), housed in a film case with prefabricated electrical terminals. Actuators are connected counterphase to a MIDÉ EL-1225 amplifier and receive an analog control signal. The sensor is a Keyence LK-G82 laser triangulation system, connected to the Keyence LK-G3001V filtering and processing unit, providing an analog voltage signal to the controller (Figure 1(c)).The external disturbance resembles transient vibration effects common in aerospace constructions [39] that are often introduced to cantilever-like structures using release tests [56–58]. Repeated impacts are delivered using the stinger mechanism, providing bursts of force at the free end of the beam using a linear motor controlled by a digital input signal (Figure 1(c)).The EKF-MPC vibration control algorithm is developed in the MATLAB/Simulink environment, where all online algorithm stages are made compatible with the real-time capability of the software. Custom blocks such as the EKF algorithm, model reassembly, discretization, prediction, and others are developed in the MATLAB m-file scripting language and recompiled into C language for the Real-Time Workshop target using the Embedded MATLAB Editor.The online Simulink block scheme of the controller algorithm is implemented on a personal computer running the MATLAB xPC Target rapid control software prototyping platform, connected to hardware components with a National Instruments PCI-6030E measurement card. The developed control system is loaded onto the Target platform via the TCP/IP protocol. The overall schematic representation of the experimental system is illustrated in Figure2.Figure 2
Schematic representation of the experimental system.The simplified block scheme of the proposed adaptive vibration control algorithm is featured in Figure3, where the essential blocks are emphasized. The position measurement is sampled in the block named“PCI-6030E AD,” and then it is scaled according to the linear amplification of the laser system and converted to meter units. The incoming signal is low- pass filtered and then the running mean is removed. This processed position measurement is used in the block named“EKF” to estimate the dynamic states and parameters of the system using (5)–(22). The estimated parameters are used to reassemble the continuous state-space model of the system in (2)–(4) by“vibSS,” and then the system matrices are discretized using the procedure in (23)–(28) in the block named“c2d.” The prediction matrices are updated according to the new model in“mpc” based on (29)–(31), along with the current LQ gain and terminal weight based on (33)–(39). The cost function defined by (46) is minimized using quadratic programming in the block named“QP Solver,” which requires the Hessian, gradient, and a state estimate; the rest of the inputs such as the constraints are fixed. The first input is scaled down to the amplifier constant and sent to the output named“PCI-6030E DA/DO.” The rest of the blocks are responsible for data logging and monitoring, timing and driving the stinger mechanism with the periodic disturbance signal fed to the digital output.Figure 3
Simplified block scheme of the proposed adaptive vibration control algorithm.To evaluate the control strategy presented in Section2, four experimental scenarios were considered. In all of these scenarios, during the first 30 s, the initial identification procedure was performed using the fixed LQ gain based on the initial parameter estimates and the settings of the MPC algorithm, excited by disturbance bursts spaced 4 s apart. This was followed by switching to the adaptive MPC control scheme and an 8 s pause; then, the disturbances spaced at 4 s apart continued. Figure 4 illustrates the experiment design with the different experimental settings. Experiment 1 (grey) assumed the original beam without added mass and the controllers switched off (open-loop control). The rest of the experiments used the EKF-MPC algorithm to generate the control inputs to the piezoceramic actuator. Experiment 2 (blue) assumed a small weight added to the end of the beam for the whole experiment, Experiment 3 (green) assumed that the weight was added manually during the pause after ~30 s, and finally Experiment 4 (red) assumed no added weight, just the original beam. The summary of experimental scenarios and settings is given in Table 1.Table 1
Summary of experimental scenarios and settings.
Exp.
Control
Weight
Changet
Color
1
✗
✗
N/A
Grey
2
✓
✓
0 s
Blue
3
✓
✓
~30 s
Green
4
✓
✗
N/A
RedFigure 4
Experiment design: measured position output and actuator input.Global sampling times were setT
=
0.01 s for all experiments, while the simulation step for the continuous part of the EKF algorithm was d
t
=
T
/
500
=
20 μs. The EKF was set with a fixed mass of m
=
0.178 g, which is the real measured mass of the beam section starting from the clamp. An offline grey-box identification procedure [59] may result in a different weight parameter; however, it is important to realize that the dynamic weight may not correspond to the physical weight and is highly dependent on the settings of the EKF.For each experimental scenario, the initial augmented state estimate was(49)
x
^
a
0
+
=
[
0
0
450
0.1
1
E
-
3
0
]
T
,
and the initial error covariance estimate matrix was set to
(50)
P
0
+
=
diag
⁡
[
1
E
-
5
1
E
-
2
1
E
+
4
1
E
-
2
1
E
-
5
1
E
-
1
]
.Measurement noise covariance for the single position output wasR
f
=
R
f
=
1
E
-
9, while the process noise covariance matrix was set to (51)
Q
f
t
=
diag
⁡
[
1
E
-
10
1
E
-
4
1
E
-
2
1
E
-
8
1
E
-
12
E
-
2
]
.Though the measurement noise covariance may seem to be too low, this is a physically meaningful setting; as the position varies in the range of millimeters and the precise optical measurement system implies only very slight inaccuracies. To prevent the estimation of negative parameters, an ad hoc clipping strategy [32, 40, 43] was used in the EKF.The continuous model was discretized with an approximation ofΨ with j
=
10 iterations. The MPC algorithm was set up using an input penalty R
=
R
=
1
E
-
10, state penalty Q
=
C
T
C, input constraints -
100
≤
u
→
(
k
)
<
100 V (see Figure 4(b)), no state constraints, and n
=
10 steps prediction horizon. The constrained quadratic programming problem was solved online using the qpOASES active-set sequential quadratic programming solver [60–62], which is specifically designed for solving optimization problems occurring in MPC. The measurement chain included a low-pass FIR filter with the passband edge at 30 Hz and stopband edge at 50 Hz and a windowed running mean filter with a 0.1 Hz fundamental frequency.


## 4. Results and Discussion

Figure5 shows a detailed view of the measured position x
1
(
k
)
=
q
(
k
) (Figure 5(a)), estimated velocity x
2
(
k
)
=
q
˙
(
k
) (Figure 5(b)), input voltage u
(
k
) (Figure 5(c)), and disturbance force F
e
(
k
) (Figure 5(d)) for a selected time window of 15 s. The initial identification procedure ends at the 30 s, just after the second disturbance, and the adaptive EKF-MPC algorithm continues to take over the control of the beam.Figure 5
Detail of the estimated system states, control input, and disturbance force.From the measured positions shown in Figure5(a), it is clear that in all experiments with the EKF-MPC (Experiments 2
–
4) the vibration levels are attenuated very effectively; in fact, settling times are shorter by approximately a factor of ten when compared to the open-loop case (Experiment1). The settling time for closed-loop control is approximately t
s
≅
1.9 s, while the open-loop system does not settle until the next impact is delivered from the stinger (t
s
≫
4 s, where t
s is 5% of the position of the largest amplitude at impact). The difference between the considered experimental scenarios is more subtle but still visible on both Figures 5(a) and 5(b), where the lightest beam with the nominal weight produces the smallest vibration amplitudes and speeds (Experiment 4), while the perturbed physical parameters caused by a weight (Experiments 2,3) result in slightly larger vibration amplitudes and speeds.In spite of the added weight and increased deformation amplitudes, the EKF-MPC algorithm compensates for the change in the dynamics by using more agile inputs, which is demonstrated in the longer lasting evolution of input voltages having greater amplitudes outside the constrained region. Note the intensity of the inputs in Figure5(c), where this effect is especially evident in the case of Experiment 3 (before and after the mass has been added) and the difference between the input behavior in the case of Experiment 2 versus Experiment 4. The larger mass in Experiment 2 results in more aggressive inputs than in Experiment 4, while the effects of the mass increase in Experiment 3 are demonstrated after a short adaptation period.These results suggest that the proposed EKF-MPC algorithm enables the system to maintain the time required to settle comparable in all closed-loop control cases (Experiments2
–
4). These matching settling times are a consequence of increased input activity on the piezoceramic elements. In addition to the adaptation feature, the constrained infinite horizon model predictive algorithm respects the process constraints imposed on the inputs throughout the experiments (Figure 5(c)). Figure 5(d) shows the estimated force disturbance. The shock-like disturbance force F
e
(
k
) is somewhat distorted by the low-pass and running mean filtering on the measured inputs; however, as one would expect, it remains consistent throughout all experimental scenarios (Experiments 1 and 4).Figure6 shows the estimated parameters and the equivalent LQ gain for the whole duration of the experiments. The estimated spring constant k
(
k
) is shown in Figure 6(a), damping constant b
(
k
) in Figure 6(b), and the force conversion constant c
(
k
) in Figure 6(c). The unconstrained equivalent LQ gain is shown for diagnostic purposes only in Figure 6(d) in order to illustrate the adaptive behavior of the EKF-MPC control algorithm.Estimated model parameters and the equivalent unconstrained LQ gain.
(a)
Stiffness
(b)
Damping
(c)
Force conversion
(d)
Equivalent unconstrained LQ gainAll the parameters simultaneously demonstrate a change to actively compensate for the mass added to the system, since the mass itself is not an extra identified but a fixed parameter. The dynamic model used in the MPC algorithm represents an equivalent point-mass-damper system instead of the real cantilever beam; therefore, the variations in the mass translate to a hypothetical change in other equivalent model parameters. These are to be understood merely as analogous changes in stiffness, damping, or actuator efficiency.Overall, the equivalent stiffness, force conversion constant, and theK
1,1 component of the LQ gain drop with increased mass, while damping increases. The sudden spikes visible in all parameters are due to the mismatch between the impulse-like disturbance of the stinger mechanism and the formulation of the EKF algorithm, which is expecting only centered Gaussian noise as disturbance [24–26]. Of the three identified parameters, stiffness was the most stable and consistent one, while the online identification of the damping parameter proved to be somewhat challenging.In case the mass is kept constant during the entire duration of the test (Experiments1, 2, and 4), all the parameters adapt to the changes in the system dynamics during the initial identification phase and then are kept on these levels even after the EKF-MPC algorithm takes over. When the mass is added to the beam just before 40 s (Experiment 3), the adaptation process starts only after this time. The parameter variations are translated as the model updates in the MPC algorithm.In practical applications, especially with joint state and parameter estimation, the excitation of the dynamic system may be absent at certain operation times. The proposed method works well with information-rich processes, however, it has not been tested with nonpersistent excitation signals. The EKF is known to degrade in performance or even become unstable in the absence of excitation as a result of ill-conditioned numerical operations [63].


## 5. Conclusions

An adaptive active vibration control algorithm based on the extended Kalman filter for real-time parameter estimation and the infinite horizon dual-mode constrained model predictive algorithm has been proposed for fixed-free cantilevers in this paper. The EKF-MPC algorithm has been implemented and experimentally tested in real time on the laboratory system featuring an aluminum beam with piezoceramic actuation and position feedback.Four experimental scenarios were used to evaluate the proposed EKF-MPC adaptive control algorithm. The closed-loop settling times were reduced by a factor of ten in all cases, while the controller promptly reacted to changes in the beam dynamics. The mass increase, emulated by the added weight, induced change in the estimated parameters of the equivalent point-mass-damper system. This mass increase resulted in the harder-to-control system requiring increased control inputs. In fact, according to the experimental results presented here, the parameter changes translated to the MPC algorithm resulted in more aggressive control moves.The proposed EKF-MPC algorithm unifies the key features of both methods, introducing adaptivity and constraint handling to the vibration control of flexible cantilever beam-like structures.


### 5.1. Future Work

Certain aspects of the unmodeled dynamics, namely, the effect of the outside disturbance, may cause the identification algorithm to diverge (not shown in the experiments here). These drawbacks will be addressed in an upcoming work [64] by the introduction of spectrum shaping filters [25, 26] into the augmented models and possibly by selective error covariance matrix update disabling during disturbances.The performance and stability of the proposed algorithm have not been tested in operation modes without significant outside excitation. The nonpersistent excitation may result in an ill-posed numerical problem and stability problems in the online EKF component. A possible solution to this problem may be the use of the moving horizon observer instead of the EKF, as it is suggested here. Instead of producing estimates based on the last know measurement and a recursive procedure, the MHO uses a moving window of past measurements. While this alone may increase the efficiency of the estimation procedure faced with nonpersistent data, the MHO variants enhanced by regularization mechanisms to cope with this situation have been suggested [63]. The price of eliminating these problems using the MHO is a heavy increase in computation cost, as the MHO requires the use of online nonlinear constrained optimization procedures. The use of MHO to estimate the states and parameters of vibration dynamics has been demonstrated in simulation [32] and using offline measurement data for a nanopositioner mechanism [65]. Real-time use of the MHO for vibration mechanics and adaptive vibration control is currently under development, but an online experimental implementation of the estimator alone has been recently proposed [66].Another relevant matter requiring attention in upcoming research is the question of the stability of the MPC component of the proposed algorithm. Stability guarantees may be given for linear systems and nominal models by enforcing recursive constraints beyond the control horizon [46, 47]; however, the stability issue of MPC in face of uncertainties—also known as stochastic MPC—is a more complex issue [67]. Recent advancements in the field of uncertain and robust MPC formulations can address this problem effectively [68–70].An important aspect of self-reliance in an adaptive structural control scheme is process monitoring and fault diagnostics; thus, in addition to the adaptive features discussed in this paper, online system diagnostics can be also easily included in the formulation, since adaptive control requires some form of online system identification. In addition to the advanced model-based fault diagnostics methods [23, 71], recent advances in data-driven fault tolerant control [72–74] can be an attractive way to increase the self-reliance of active vibration control.


## 5.1. Future Work

Certain aspects of the unmodeled dynamics, namely, the effect of the outside disturbance, may cause the identification algorithm to diverge (not shown in the experiments here). These drawbacks will be addressed in an upcoming work [64] by the introduction of spectrum shaping filters [25, 26] into the augmented models and possibly by selective error covariance matrix update disabling during disturbances.The performance and stability of the proposed algorithm have not been tested in operation modes without significant outside excitation. The nonpersistent excitation may result in an ill-posed numerical problem and stability problems in the online EKF component. A possible solution to this problem may be the use of the moving horizon observer instead of the EKF, as it is suggested here. Instead of producing estimates based on the last know measurement and a recursive procedure, the MHO uses a moving window of past measurements. While this alone may increase the efficiency of the estimation procedure faced with nonpersistent data, the MHO variants enhanced by regularization mechanisms to cope with this situation have been suggested [63]. The price of eliminating these problems using the MHO is a heavy increase in computation cost, as the MHO requires the use of online nonlinear constrained optimization procedures. The use of MHO to estimate the states and parameters of vibration dynamics has been demonstrated in simulation [32] and using offline measurement data for a nanopositioner mechanism [65]. Real-time use of the MHO for vibration mechanics and adaptive vibration control is currently under development, but an online experimental implementation of the estimator alone has been recently proposed [66].Another relevant matter requiring attention in upcoming research is the question of the stability of the MPC component of the proposed algorithm. Stability guarantees may be given for linear systems and nominal models by enforcing recursive constraints beyond the control horizon [46, 47]; however, the stability issue of MPC in face of uncertainties—also known as stochastic MPC—is a more complex issue [67]. Recent advancements in the field of uncertain and robust MPC formulations can address this problem effectively [68–70].An important aspect of self-reliance in an adaptive structural control scheme is process monitoring and fault diagnostics; thus, in addition to the adaptive features discussed in this paper, online system diagnostics can be also easily included in the formulation, since adaptive control requires some form of online system identification. In addition to the advanced model-based fault diagnostics methods [23, 71], recent advances in data-driven fault tolerant control [72–74] can be an attractive way to increase the self-reliance of active vibration control.



---

*Source: 741765-2014-05-07.xml*