Discrete-time extremum-seeking for Wiener-Hammerstein plants
Automatica 50 (2014) 2998–3008
Contents lists available at ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
Discrete-time extremum-seeking for Wiener–Hammerstein plants✩ Rohan C. Shekhar 1 , William H. Moase, Chris Manzie Department of Mechanical Engineering, The University of Melbourne, Victoria 3010, Australia
article
info
Article history: Received 24 September 2013 Received in revised form 23 March 2014 Accepted 16 July 2014 Available online 27 October 2014 Keywords: Extremum-seeking Adaptive control Optimisation Engine control
abstract This paper presents a new formulation of extremum-seeking control for discrete-time Wiener– Hammerstein plants. A novel method of analysis using Linear Parameter-Varying (LPV) system theory demonstrates semi-global stability of the control scheme. Assuming only limited plant knowledge, the stability result ensures convergence of the plant output in steady state to a point in an arbitrarily small neighbourhood of the extremum, for appropriately chosen controller parameters. The behaviour of the control scheme is analysed on a simple simulated system, prior to being implemented on an internal combustion engine. Experiments demonstrate how the scheme is able to maximise engine output torque in the presence of an uncertain fuel composition by modifying the spark timing. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Extremum-Seeking (ES) is a closed-loop control paradigm that seeks to optimise the steady-state output of a given nonlinear plant, where a plant model is unavailable or only approximately known. Assuming that the input-to-output map at steady state contains an extremum (maximum or minimum) in the region of interest, ES manipulates the plant input so that the output is contained within a small neighbourhood of this extremum at equilibrium. For a small cross-section of different ES schemes, see Chichka, Speyer, and Park (1999), Choi, Krstić, Ariyur, and Lee (2002), Drakunov, Ozguner, Dix, and Ashrafi (1995), Guay and Zhang (2003), Moase, Manzie, and Brear (2010) and Zhang and Ordonez (2007). A comprehensive review of the evolution of ES schemes throughout history can be found in Tan, Moase, Manzie, Nešić, and Mareels (2010). Most popular ES schemes emulate some form of gradient descent, using a dither signal to estimate local derivative information and an adaptation gain to regulate the speed at which the plant input approaches the steady-state optimising value (Ariyur & Krstić, 2003). The first analyses of local stability were performed for continuous-time ES schemes with nonlinear dynamics (Krstić
✩ This work was supported under the Australian Research Council’s Discovery Projects funding scheme (Project DP120101830). The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Raul Ordóñez under the direction of Editor Miroslav Krstic. E-mail addresses: [email protected] (R.C. Shekhar), [email protected] (W.H. Moase), [email protected] (C. Manzie). 1 Tel.: +61 3 83447714; fax: +61 3 8344 4290.
http://dx.doi.org/10.1016/j.automatica.2014.10.030 0005-1098/© 2014 Elsevier Ltd. All rights reserved.
& Wang, 2000). Non-local, or semi-global stability results were subsequently derived in Tan, Nešić, and Mareels (2006). Analysis of these ES schemes generally involves singular perturbation analysis and averaging theory; however, recent approaches have also used Lie bracketing techniques (Dürr, Stanković, & Johansson, 2011). Grey-box ES schemes using limited plant knowledge have also recently been shown to guarantee semi-global stability (Nešić, Mohammadi, & Manzie, 2013; Sharafi, Moase, Shekhar, & Manzie, 2013). Semi-globally stable ES for continuous-time Hammerstein and Wiener–Hammerstein plants was considered in Moase and Manzie (2012a,b) respectively. Such plants consist of a static nonlinear element with either output LTI dynamics (Hammerstein) or both input and output LTI dynamics (Wiener–Hammerstein). These plant models arise in the black-box identification of nonlinear dynamic systems (Nešić, 1999) and have proven useful for representing biological systems (Hunter & Korenberg, 1986), diesel engines (Perez, Blasco, García-Nieto, & Sanchis, 2006) and chemical processes (Eskinat, Johnson, & Luyben, 1991). Various discrete-time ES schemes have also been presented in the literature. Specifically, Teel and Popović (2001) outline an approach that waits for a sufficient period of time until system dynamics have settled before applying subsequent control actions, effectively treating the plant as a static nonlinearity. Discrete dynamics are considered in Choi et al. (2002), which demonstrates local stability for a Wiener–Hammerstein plant with a quadratic nonlinearity. In contrast to these methods, a dither-free sampleddata approach is considered in Khong, Manzie, Nešić, and Tan (2012). It should be noted that all of these approaches rely on explicit time-scale separation arguments so that the nonlinear
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
system approximates a static map; this can limit convergence speed. It is possible to somewhat relax time-scale separation requirements for ES on Wiener–Hammerstein plants in continuous time, as demonstrated by Moase and Manzie (2012a). In discrete time, Shekhar, Moase, and Manzie (2013) laid the foundation for achieving the same goal, using a new method of stability analysis for ES on a static nonlinear plant. Specifically, Linear Difference Inclusion (LDI) theory and Lyapunov arguments are used to show semiglobal stability for sufficiently small adaptation gains, assuming a strongly-convex plant map. This paper builds upon the work in Shekhar et al. (2013), presenting a new formulation of ES for discrete-time Wiener– Hammerstein plants with nonlinear maps satisfying the same convexity requirement. Extending the techniques developed in Shekhar et al. (2013), Linear Parameter-Varying (LPV) system theory is applied to prove semi-global stability, which improves upon the local stability results of Choi et al. (2002). The only prior knowledge required is the sign of the system’s dynamic response evaluated at the dither frequency, which can be experimentally determined for real systems; this is less restrictive than the knowledge requirement for grey-box schemes. No explicit timescale separation between the plant dynamics and gradient estimation is required, potentially improving convergence speed as demonstrated in simulation. The scheme presented in this paper also exhibits convergence of the plant output to a point in an explicitly characterised neighbourhood of the extremum, in contrast to most existing ES approaches that result in limit-cycle equilibria. An experimental example is used to illustrate the performance of the ES scheme on a spark-ignition engine. In particular, the ES scheme is used to maximise the mean torque output at steady state by modifying the spark advance angle. Different ES schemes have already been shown to improve efficiency through modification of engine tuning parameters (Mohammadi, Manzie, & Nešić, 2012; Popovic, Jankovic, Magner, & Teel, 2003), making this an ideal candidate application to demonstrate the efficacy of the scheme developed in this paper. The scheme is shown to improve torque output, motivating its use for closed-loop tuning of alternative fuel engines where the fuel composition may vary. 1.1. Nomenclature The sets of natural numbers, integers and real numbers are denoted as N, Z and R respectively. The set of all non-negative integers is denoted as Z≥0 . The column vector [aT , bT , c T ]T is written as [a; b; c ]. In represents the n × n identity matrix. The operator ⊕ represents the direct sum of matrices. The interior of a set is indicated by int(·). For some P ≻ 0 and a square matrix A of the same dimension, define the discrete Lyapunov operator by
R(P , A) := P − AT PA. The operator λi (M ) denotes the ith eigenvalue of a matrix M and vi (M ) denotes the corresponding eigenvector, where the eigenvalues are enumerated in some specified order. λ(M ) and λ(M ) denote the maximum and minimum eigenvalues of a symmetric matrix M respectively. The 2-norm of a matrix or vector A is written as ∥A∥. The ⋆ symbol is used in symmetric block matrices to denote the transpose of the diagonally opposite block. To illustrate this, given compatible matrices A = AT , C = C T and B,
A B
⋆ C
:=
A B
BT . C
2. Problem formulation This paper is concerned with optimising the steady-state output of the Wiener–Hammerstein plant shown in Fig. 1. The plant is
2999
Fig. 1. Discrete-time Wiener–Hammerstein plant.
described by a static nonlinear map f : R → R, with linear input dynamics G(z ) and output dynamics H (z ). The dynamics are assumed to have minimal state-space realisations defined by the equations xG (k + 1) = AG xG (k) + BG u(k) p(k) = CG xG (k) + DG u(k) xH (k + 1) = AH xH (k) + BH q(k) y(k) = CH xH (k) + DH q(k), for states xG and xH of the input and output dynamics respectively. Whilst explicit forms for the mapping and LTI systems are unknown, they are assumed to satisfy the following properties: Assumption 1. There exist positive constants α, β such that α ≤ f ′′ (p) ≤ β . This implies that the nonlinear map f (·) is twice differentiable and strongly convex. It also guarantees the existence of a point p∗ ∈ R such that f ′ (p∗ ) = 0. Assumption 2. The matrices AG and AH have spectral radii strictly less than unity. Assumption 3. Without loss of generality, the LTI dynamics satisfy G(1) = H (1) = 1, so that any steady-state gain is considered as part of f (·). Remark 4. Assumption 1 guarantees that p∗ is the unique minimiser of f (·). Assumption 2 ensures that G(·) and H (·) are stable. Remark 5. Assumption 1 is less restrictive than assuming a quadratic steady-state map, which is used extensively in the literature for local stability analysis (e.g. Choi et al., 2002). Indeed, it will be shown that α can be arbitrarily small and β can be arbitrarily large whilst still ensuring closed-loop stability. 2.1. Controller design To find the steady-state input that minimises f (·), the control scheme in Fig. 2 will be used. The input u¯ (k) represents an estimate for the minimiser of f (·) at time step k. A simple dither signal a(−1)k is added to this estimate and applied to the Wiener–Hammerstein plant, where a > 0 is the amplitude. The resulting output is then ‘‘demodulated’’ by the signal (−1)k , before being passed through a two-step FIR averaging filter and scaled by an adaptation gain κ > 0 having sign σ ∈ {−1, 1}. This quantity, which is approximately proportional to f ′ (¯u(k)), is then subtracted from u¯ (k), essentially performing discrete-time integration. The sign σ is used to correct for any phase shift induced by G(·) and H (·). The closed-loop dynamics of Fig. 2 are described by u(k) = u¯ (k) + a(−1)k
(1a)
xG (k + 1) = AG xG (k) + BG u(k)
(1b)
p(k) = CG xG (k) + DG u(k)
(1c)
q(k) = f (p(k))
(1d)
xH (k + 1) = AH x(k) + BH q(k)
(1e)
y(k) = CH x(k) + DH q(k)
(1f)
v(k) = y(k)(−1)
(1g)
k
σκ
(v(k − 1) + v(k)), (1h) 2 where the initial conditions u¯ (0), xG (0), xH (0) and v(−1) are arbitrary.
u¯ (k + 1) = u¯ (k) −
3000
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Fig. 2. ES scheme for Wiener–Hammerstein plant.
Fig. 3. Loop transformation.
Remark 6. This ES scheme is one of many approaches to find the extremum of a Wiener–Hammerstein plant. More traditional schemes would use an integrator after the demodulator, but it will be seen that the averaging filter provides desirable properties. The analysis presented in this paper provides a general methodology that could be easily applied to other discrete-time ES schemes with different kinds of filters. Remark 7. The ES scheme and analysis to follow can be trivially modified to find the steady-state maximum when f (·) is a strongly concave function, by replacing κ with −κ .
˘ (z ) = G(z ). G
(8)
The modified dither signal is then given by a˜ = G(−1)a.
(9)
Note that this may be negative depending on any phase inversion introduced by G(·). For the output dynamics, (1e) and (1f) can be written in the form xH (k + 1)(−1)k+1 = −AH xH (k)(−1)k − BH q(k)(−1)k y(k)(−1)k = CH xH (k)(−1)k + DH q(k)(−1)k = v(k).
3. Analysis of stability Analysis of (1) will be carried out in four steps. Firstly, a loop transformation will be applied to Fig. 2 in order to shift the dynamics into the gradient estimator portion of the closed loop (Section 3.1). Following this, Fourier analysis will be utilised in order to identify an equilibrium solution (Section 3.2). In terms of this solution, various state transformations will then be applied to write the dynamics in terms of an error system having an equilibrium at the origin (Section 3.3). Finally, the error system dynamics will be cast into the form of an LPV system (Section 3.4) that will be analysed with Lyapunov techniques (Section 3.5). To simplify notation, the input-to-state transfer functions
ΓG (z ) := (zI − AG ) BG
(2)
ΓH (z ) := (zI − AH )−1 BH
(3)
−1
which demonstrates that
By defining x˘ H (k) := xH (k)(−1)k
(10)
q˘ (k) := q(k)(−1)k ,
˘ (z ) can be written in state-space form as the system H x˘ H (k + 1) = −AH x˘ H (k) − BH q˘ (k)
(11)
v(k) = CH x˘ H (k) + DH q˘ (k),
(12)
which shows that
˘ (z ) = H (−z ). H
(13)
This completes the transformation. 3.2. Identifying an equilibrium solution
will be used extensively, together with the identities BG − z ΓG (z ) = −AG ΓG (z )
(4)
BH − z ΓH (z ) = −AH ΓH (z ),
(5)
To find an equilibrium solution for u¯ , note that from Fig. 3, q(k) = f (˘p(k) + a˜ (−1)k ).
which can be easily verified by rearranging (2) and (3).
The Discrete-Time Fourier Series (DTFS) can now be used to write
3.1. Loop transformation
q(k) = qˆ 0 (˘p(k)) + qˆ 1 (˘p(k))(−1)k , where the DTFS coefficients as functions of p˘ are given by
Analysis of the closed-loop system (1) can be simplified by swapping the order of the G(·) and the summation block, as well as H (·) and the product block, as shown in Fig. 3. From (1a) and (1b),
qˆ 0 (˘p) :=
xG (k + 1) = AG xG (k) + BG (¯u(k) + a(−1)k ).
qˆ 1 (˘p) :=
Now define the state x˘ G (k) := xG (k) − aΓG (−1)(−1)k , which is the error between xG (k) and its steady-state value to the excitation caused by the dither signal. In terms of this new state, ˘ (·) can be shown to have a state-space realisation G x˘ G (k + 1) = AG x˘ G (k) + BG u¯ (k)
(6)
p˘ (k) = CG x˘ G (k) + DG u¯ (k),
(7)
1 2 1 2
(f (˘p + a˜ ) + f (˘p − a˜ ))
(14a)
(f (˘p + a˜ ) − f (˘p − a˜ )).
(14b)
For a fixed p˘ , this harmonically decomposes q into a mean component and an oscillating component. In terms of these components, q˘ (k) = qˆ 1 (˘p(k)) + qˆ 0 (˘p(k))(−1)k ,
(15)
which has essentially swapped the harmonic components of q. Theorem 8 exploits the harmonic decomposition to identify an equilibrium solution to the closed-loop system (1).
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
forcing signal q˘ (k) (15). Using x¯ H (k, p˘ ) to denote the quasi-steady solution,
Theorem 8. There exists a unique point p ∈ (p − |˜a|, p + |˜a|) 0
∗
3001
∗
x¯ H (k, p˘ ) := ΓH (−1)ˆq1 (˘p) + ΓH (1)ˆq0 (˘p)(−1)k .
such that
This is then used to define the error state
u¯ (k) = p0
(16a)
xG (k) = ΓG (1)p0 + aΓG (−1)(−1)k
x˜ H (k) := x˘ H (k) − x¯ H (k, p˘ (k)).
(16b)
xH (k) = ΓH (1)ˆq0 (p0 )
(16c)
y(k) = qˆ 0 (p0 )
(16d)
From (10) and (16c), it can be seen that the equilibrium value for x˜ H is zero if p˘ (k) = p˘ (k − 1) = p0 . Using (5) and (11), the dynamics of this error state can be written in the form
is an equilibrium solution to (1).
x˜ H (k + 1) = x˘ H (k + 1) − x¯ H (k + 1, p˘ (k + 1))
= −AH x˜ H (k) − ΓH (−1)1qˆ 1 (˘p(k + 1))
Proof. Applying the mean value theorem to (14b) and invoking Assumption 1, a˜ −1 qˆ ′1 (˘p) =
1
1qˆ n (˘p(k)) := qˆ n (˘p(k)) − qˆ n (˘p(k − 1)),
where c ∈ (˘p − |˜a|, p˘ + |˜a|). This shows that a q1 (·) is strictly monotonically increasing. Then, since p∗ is the global minimiser of f (·), it can be seen that
˜ −1 ˆ
a˜ −1 qˆ 1 (p∗ − |˜a|) = a˜
(21)
where
(f ′ (˘p + a˜ ) − f ′ (˘p − a˜ ))
2a˜ = f ′′ (c ) > 0,
−1
+ ΓH (1)1qˆ 0 (˘p(k + 1))(−1)k ,
qˆ 1 (p + |˜a|) = ∗
1 2|˜a| 1 2|˜a|
(f (p∗ ) − f (p∗ − 2|˜a|)) < 0 (f (p + 2|˜a|) − f (p )) > 0. ∗
∗
qˆ 1 (p0 ) = 0.
for n ∈ {0, 1}. Eq. (12) can then be written in the form
v(k) = CH x˜ H (k) + H (−1)ˆq1 (˘p(k)) + qˆ 0 (˘p(k))(−1)k ,
(17)
Furthermore, this point is unique by monotonicity. Having established the existence of a unique p0 , (16) can be verified by substitution into (1), applying Assumption 3 and (17). In particular, it can be seen that the input to the FIR filter is given by
v(k) = qˆ 0 (p0 )(−1)k . Since this is completely removed by the filter, u¯ remains unchanged with value p0 .
(23)
which follows from (13) and Assumption 3. Finally, by defining the translations u˜ := u¯ − p0
(24)
0
(25)
p˜ := p˘ − p
Hence, by the intermediate value theorem, there exists a point p0 ∈ (p∗ − |˜a|, p∗ + |˜a|) such that
(22)
and substituting these into (19), (20), (21), (23), and (1h), the dynamics of the error system are given by x˜ G (k + 1) = AG x˜ G (k) − ΓG (1)1u˜ (k + 1)
(26a)
p˜ (k) = CG x˜ G (k) + u˜ (k)
(26b)
x˜ H (k + 1) = −AH x˜ H (k) − ΓH (−1)1qˆ 1 (˜p(k + 1) + p ) 0
+ ΓH (1)1qˆ 0 (˜p(k + 1) + p0 )(−1)k v(k) = CH x˜ H (k) + H (−1)ˆq1 (˜p(k) + p ) + qˆ 0 (˜p(k) + p0 )(−1)k σκ u˜ (k + 1) = u˜ (k) − (v(k − 1) + v(k)) ,
(26c)
0
2
(26d) (26e)
where 1u˜ (k) := u˜ (k) − u˜ (k − 1). 3.3. Error system formulation For a constant or slowly-varying input u¯ , (2) and (8) allow the ˘ (·) to be written as quasi-steady state of G x¯ G (¯u) := ΓG (1)¯u.
3.4. LPV system formulation
Then define the new state x˜ G (k) := x˘ G (k) − x¯ G (¯u(k)).
(18)
Using (4) and (6), it can be seen that x˜ G (k + 1) = x˘ G (k + 1) − x¯ G (¯u(k + 1))
= AG x˜ G (k) − ΓG (1)1u¯ (k + 1),
The error system will be further transformed into LPV form. In order to accomplish this, define the functions g (˜p) :=
(19)
where 1u¯ (k) := u¯ (k) − u¯ (k − 1). In addition, (7), (8), (18) and Assumption 3 can be combined to give p˘ (k) = CG (˜xG (k) + x¯ G (¯u(k))) + DG u¯ (k)
= CG x˜ G (k) + u¯ (k).
Remark 9. These equations describe the closed-loop dynamics (1) exactly, without any approximation. The quasi-steady solutions have only been used to formulate an error system for facilitating stability analysis.
(20)
Now, observe that the transfer function from q˘ to x˘ H is given by ΓH (−z ), from (3) and (13). Then, if p˘ is fixed or slowly varying with time, the quasi-steady solution for x˘ H can be evaluated as the frequency response of ΓH (−z ) to the approximately periodic
1 a˜
p˜
qˆ 1 (µ + p0 )dµ
(27)
0
h(˜p) := qˆ 0 (˜p + p0 ).
(28)
The following lemmas describe some important properties of these functions. Lemma 10. For all p˜ ∈ R, g (·) satisfies the bounds
α ≤ g ′′ (˜p) ≤ β. It also has a unique minimum at the origin. Proof. See Appendix A.
3002
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Lemma 11. For all p˜ ∈ R, h(·) satisfies the bounds
where
α ≤ h (˜p) ≤ β.
φ(θ (k)) := a˜ H (−1)g ′′ (ξ1 (k)) + h′ (ξ2 (k))(−1)k .
′′
Proof. See Appendix B.
To complete the transformation, (26b), (26d) and (37) are substituted into (26e) to give
σκ
Remark 12. These results demonstrate that g (˜p) and h(˜p) are closely related to f (˜p + p0 ), being strongly convex and satisfying the same second derivative bounds as f (·).
1u˜ (k + 1) =
Lemma 13. The function h(·) satisfies the bounds
which uses the fact that
2
(2a˜ H (−1)ζ (k)(CG x˜ G (k) + u˜ (k))
− w(k) + CH x˜ H (k)),
(40)
from (29) and (33). Combining (26a), (26e) and (34)–(40), the LPV system description is given by
In terms of (25), (27) and (28), the DTFS coefficients (14) are given by qˆ 1 (˘p(k)) = a˜ g ′ (˜p(k))
(29)
qˆ 0 (˘p(k)) = h(˜p(k)).
(30)
To allow the LPV system to be defined, it is useful to apply the mean value theorem to (22), using (26a), (26b), (29) and (30) to write
1qˆ 1 (˘p(k + 1)) = a˜ g ′′ (ξ1 (k))1p˜ (k + 1) 1qˆ 0 (˘p(k + 1)) = h′ (ξ2 (k))1p˜ (k + 1),
X (k + 1) = A(θ (k), κ)X (k),
1p˜ (k + 1) := p˜ (k + 1) − p˜ (k) = CG (AG − I )˜xG (k) + DG 1u˜ (k + 1)
(31)
ξ1 (k), ξ2 (k) ∈ (min{˜u(k), u˜ (k + 1)}, max{˜u(k), u˜ (k + 1)}).
(32)
Also define p˜ (k) ̸= 0
(33)
p˜ (k) = 0.
(41)
where X (k) := [˜u(k); w(k); x˜ H (k); x˜ G (k)] A(θ , κ) := A0 (θ ) −
1 0 A0 (θ ) := 0 0
0 0 0 0
σκ ¯ A(θ )
2 0 −C H −AH 0
0 φ(θ )CG (AG − I ) Φ (θ )CG (AG − I ) −AG
T
2aG(−1)H (−1)ζ (k) DG φ(θ ) −1 ¯A(θ ) := D Φ (θ ) . CH G −ΓG (1) 2aG(−1)H (−1)ζ (k)CG
where
′ g (˜p(k)) ζ (k) := p˜ (k) ′′ g (0)
2
qˆ 1 (˜p(k) + p0 ) = 2a˜ ζ (k)˜p(k),
1
|h′ (˜p)| ≤ |f ′ (˜p + p0 )| + a˜ (β − α). Proof. See Appendix C.
(39)
1
Note that zeros in the matrices may represent blocks of zeros having appropriate dimensions. 3.5. LPV Lyapunov analysis This section will prove the following convergence result, which is the primary contribution of this paper.
It can be easily shown that g ′ (·) is sector bounded, so
α ≤ ζ (k) ≤ β
Theorem 14 (ES Scheme Convergence). As k → ∞,
for all k. In terms of (27), (28) and (33), define the parameter vector
u¯ (k) → p0 ∈ (p∗ − |˜a|, p∗ + |˜a|)
θ(k) := [θ1 (k); θ2 (k); θ3 (k); θ4 (k); θ5 (k)] := a˜ ζ (k); a˜ g ′′ (ξ1 (k)); h′ (ξ2 (k))(−1)k ;
exponentially and y(k) → qˆ 0 (p0 ), for any initial conditions of (1), any positive α, β ∈ Z≥0 , α < β and any a˜ ∈ Z, provided κ is sufficiently small and
a˜ ζ (k)g ′′ (ξ1 (k)); a˜ ζ (k)h′′ (ξ2 (k))(−1)k .
(34)
These definitions can then be used to write (26c) in the form x˜ H (k + 1) = −AH x˜ H (k) + Φ (θ (k))(CG (AG − I )˜xG (k)
+ DG 1u˜ (k + 1)),
(35)
σ G(−1)H (−1) > 0. Proof. The intermediate results required to prove the theorem are detailed in Section 3.5.1, while the proof of the theorem itself is provided in Section 3.5.2.
where
Φ (θ(k)) := −˜ag ′′ (ξ1 (k))ΓH (−1) + h′ (ξ2 (k))(−1)k ΓH (1).
(36)
Now define the new state
w(k) := −v(k − 1) + H (−1)ˆq1 (˜p(k) + p0 ) −ˆq0 (˜p(k) + p0 )(−1)k .
(37)
Using (26d) and (31), it follows that
w(k + 1) = −CH x˜ H (k) + φ(θ (k))(CG (AG − I )˜xG (k) + DG 1u˜ (k + 1)),
(38)
Remark 15. Theorem 14 presents a semi-global stability result, since it states that for any initial condition, there is always a κ small enough to ensure convergence. The requirement σ G(−1)H (−1) > 0 ensures that any inversion of the dither signal induced by G(·) and H (·) is compensated by σ . It also implies that both LTI systems cannot have zeros at z = −1, which would remove the dither signal completely in steady state. Note also that the theorem does not require |˜a| to be small, unlike most other ES schemes in the literature. Hence, a can be chosen based on the desired proximity to the true optimum, but large values may necessitate the use of a smaller κ .
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
3003
3.5.1. Supporting results
Proof. See Appendix F.
Lemma 16. If for some κ0 there exists a matrix P that satisfies the parameter-dependent LMI condition
Since G(·) and H (·) are stable, there exist positive-definite matrices PG and PH such that
P
PA(θ , κ0 )
⋆
P
R(PG , AG ) ≻ 0
≻ 0,
R(PH , AH ) ≻ 0.
for all θ contained in some parameter set Θ , then V (X ) := X PX T
(42)
is a Lyapunov function for the LPV system (41) on a subset of the state space X corresponding to θ ∈ Θ , for κ = κ0 . Proof. This result follows from applying the Schur complement to the Discrete Lyapunov Equation (DLE), P − AT (θ , κ0 )PA(θ , κ0 ) = R(P , A(θ , κ0 )) ≻ 0.
(43)
From elementary Lyapunov theory, if there exist P ≻ 0 and κ0 > 0 satisfying (43) for all θ ∈ Θ , then (42) is a Lyapunov function for the LPV system on X. Lemma 17 (Block Upper-Triangular System Stability). Given a block upper-triangular matrix with square diagonal entries
A :=
A11 0
A12 , A22
if there exist P11 ≻ 0 and P22 ≻ 0 such that
Using these matrices, stability of the LPV system can now be shown for suitably chosen κ . Theorem 20 (Existence of Lyapunov Matrix). If ∥Ψ (θ )∥ is bounded for some θ and σ G(−1)H (−1) > 0, then
R(P , A(θ , κ)) ≻ 0, where P := 1 ⊕ δγ ⊕ δ PH ⊕ PG ,
Proof. The proof proceeds by considering the case κ = 0, and showing that as κ is increased by a small amount, the DLE (47) is satisfied. It utilises the fact that 1 0 A(θ , 0) = A0 (θ ) = 0 0
0
A1 :=
0 φ(θ )CG (AG − I ) Φ (θ )CG (AG − I ) −AG
−CH , −A H
0 0
P1 := γ ⊕ PG ,
R(P , A) ≻ 0,
(44)
for some sufficiently small γ > 0. Lemma 17 then guarantees that R(P1 , A1 ) ≻ 0. Consider now the larger submatrix
for all −1 , AT11 )) ∥A12 ∥−2 . 0 < γ < λ(R(P22 , A22 ))λ(R(P11
A2 (θ ) := (45)
Proof. See Appendix D.
Lemma 18 (Eigenvalue Derivatives). Consider a symmetric matrixvalued function M : R → RN , assumed to be differentiable on some interval D ⊆ R. The derivative of any distinct eigenvalue λi (M (κ)) with respect to κ for κ ∈ int(D ) is given by the expression
= vTi (M (κ))
Proof. See Appendix E.
dκ
A1 0
=
dM (κ)
where (46) has been substituted. Since ∥Ψ (θ )∥ is bounded, it follows that the quantity
∥Ψ (θ )CG (AG − I )∥ is also bounded, so Lemma 17 can be applied again to show that
= δγ ⊕ δ PH ⊕ PG satisfies R(P2 , A2 (θ )) ≻ 0, for a sufficiently small δ > 0. The overall candidate Lyapunov matrix
Lemma 19. If |˜u| ≤ L1 and ∥xG ∥ ≤ L2 , then
a˜ H (−1) ∥Ψ (θ)∥ ≤ −ΓH (−1) β H (1) 1 + β(L1 + ∥CG ∥ L2 ) + a˜ (β − α) , ΓH (−1) 2
P := 1 ⊕ P2
where
M (θ , κ) :=
is a matrix comprising the parameters defined in (36) and (39).
φ(θ )CG (AG − I ) Φ (θ )CG (AG − I ) 0 −AG Ψ (θ )CG (AG − I ) , −AG −CH −A H
P2 := δ P1 ⊕ PH
vi (M (κ)).
Ψ (θ) := [φ(θ ); Φ (θ )]
0 0 0
Furthermore, such a γ always exists.
dκ
0 −CH −A H 0
define
satisfies the DLE
dλi (M (κ))
0 0 0 0
the matrix 0 P22
(48)
for any positive α, β and sufficiently small κ, γ and δ .
R(P22 , A22 ) ≻ 0,
γ P11
(47)
has a block upper-triangular structure. Considering the submatrix
R(P11 , A11 ) ≻ 0
P :=
= 1 ⊕ δγ ⊕ δ PH ⊕ PG
(49)
is now defined in terms of γ , δ, PG and PH . Applying the LMI stability condition of Lemma 16, define the matrix function
(46)
PA(θ , κ)
=
⋆
P
P
PA0 (θ )
P
⋆ P
−
σκ 2
0 P A¯ (θ )
⋆ 0
.
3004
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Satisfaction of (47) is then equivalent to M (θ , κ) being positive definite. Given the block diagonal structure of P, it is useful to permute the blocks of M (·, ·). Note that for any invertible matrix S,
corresponding to the initial states of the original system. Then, define a level set of the quadratic form X T PX passing through this state by the ellipsoid
M (θ, κ) ≻ 0 ⇐⇒ S −1 M (θ , κ)S ≻ 0.
˜ (0)}, E := {X |X T PX ≤ X T (0)PX
This is effectively a similarity transformation, which preserves the eigenvalues of M (θ , κ) for all parameter values. Since P (49) is expressed as a 2 × 2 block matrix, M (·, ·) can be expressed in the form of a 4 × 4 block matrix. Choosing S to be the (orthogonal) permutation matrix that reorders the block columns {1, 2, 3, 4} to {1, 3, 2, 4}, the new matrix
where P is chosen as in (48), for some γ and δ , and P˜ := 1 ⊕ 1 ⊕ PH ⊕ PG ,
˜ (θ, κ) := S T M (θ , κ)S M
which is the matrix P with γ = δ = 1. Any positive value can ˜ but they are set to unity for ease of be assigned to γ and δ in P, analysis. Given the block diagonal structure of P, it is useful to define E in the equivalent form
can be defined. For κ = 0, this matrix has the form
E = {(˜u, w, x˜ H , x˜ G ) | u˜ 2 + δγ w 2 + δ˜xTH PH x˜ H
˜ (θ, 0) = Q1 ⊕ Q2 (θ ) M where
Q1 :=
1 1
1 , 1
Q2 (θ ) :=
P2 P2 A2 (θ )
⋆ P2
.
˜ (0)}. + x˜ TG PG x˜ G ≤ X T (0)PX For any γ , δ ≤ 1, it is clear that X (0) ∈ E . Note also that independent of the values given to γ and δ , the block-diagonal structure of P guarantees that
sup |˜u| =
Since the eigenvalues of a block diagonal matrix are simply the ˜ (θ , 0) is positive semidefieigenvalues of the individual blocks, M nite. It has the same eigenvalues as Q1 , namely 0 and 1, in addition to the eigenvalues of Q2 , which are all strictly positive for sufficiently small γ and δ .
˜ (·, ·) as a 4 × 4 block matrix, the zero Continuing to treat M ˜ (θ , 0)), which is unique from Assumption 2, has eigenvalue λ1 (M the corresponding normalised eigenvector (in block form) ˜ (θ, 0)) = √1 [1; −1; 0; 0]. v1 ( M 2 ˜ (θ , 0)) with respect to Applying Lemma 18, the derivative of λ1 (M
κ at κ = 0 is given by
˜ ˜ (θ , 0)) ∂λ1 (M ˜ (θ , 0)) ∂ M (θ , 0) v1 (M ˜ (θ , 0)) = vT1 (M ∂κ ∂κ = σ G(−1)H (−1)aζ (k). Since α ≤ ζ (·) ≤ β , it is evident that
˜ (θ , 0) ∂M σ G(−1)H (−1)aα ≤ ≤ σ G(−1)H (−1)aβ. ∂κ Furthermore, since σ G(−1)H (−1) > 0, the derivative is positive, independent of the values taken on by the other components of θ . ˜ (θ , κ)) will move into the right half-plane This means that λ1 (M as κ is increased by a small amount from zero. The continuity of ˜ (θ, κ) in κ also ensures that there will be some small κ for which M the remaining positive eigenvalues will still be situated in the right ˜ (θ , κ) and hence M (θ , κ) positive definite half-plane, rendering M for sufficiently small κ . By increasing κ beyond some value κmax however, it will become indefinite; this is certain to occur when κ is large enough to cause A(θ , κ) to have eigenvalues outside the unit disc. It can therefore be concluded that for a sufficiently small γ and δ , the Lyapunov condition (47) will be satisfied for κ < κmax . The original system (1) can now be shown to converge, for a sufficiently small κ that depends on the initial conditions. 3.5.2. Proof of ES scheme convergence (Theorem 14) Consider any initial state of the LPV system X (0) = [˜u(0); w(0); x˜ H (0); x˜ G (0)],
˜ (0) X T (0)PX
X ∈E
sup x˜ G =
˜ (0) X T (0)PX
X ∈E
λ(PG )
.
Hence, for any state within or on the boundary of E , Lemma 19 guarantees that ∥Ψ (θ )∥ < N, for some N < ∞. This bound will also take into account the initial conditions u˜ (−1) and x˜ G (−1), since |h′ (ξ2 (0))| depends on them (see the proof of Lemma 19 in Appendix F). Then, from Lemma 17 and Theorem 20, it is possible to select κ , γ ≤ 1 and δ ≤ 1 sufficiently small without affecting this bound to ensure that R(A(θ , κ), P ) ≻ 0, for all θ such that ∥Ψ (θ )∥ < N; this is true for all X ∈ E . Lemma 16 then implies that the quadratic form V (X ) = X T PX is a Lyapunov Function for all X ∈ E . Hence, E is now a level set of V (·), making it an invariant set for the LPV system (41). Given that X (0) ∈ E , exponential stability of the LPV system is guaranteed, meaning that as k → ∞, u˜ → 0 exponentially. Since u¯ is related to u˜ by a translation (24), it can be concluded that u¯ exponentially converges to the equilibrium solution p0 , where
|p0 − p∗ | < |˜a|. This in turn ensures that y → qˆ 0 (p0 ) (16d). Remark 21. This ES Scheme results in a plant output converging to a constant value. This is a markedly different behaviour to most other dither-based schemes, which exhibit periodic orbits. The FIR filter with its zero at z = −1 is the main reason for this behaviour, as is evident from the nature of the equilibrium solution described in Theorem 8. Remark 22. A practitioner implementing (1) simply has to choose the dither amplitude a based upon the desired proximity to the optimum, σ to compensate for any phase-shift introduced by the dynamics and κ to be sufficiently small for stability. This is a notable improvement on existing ES methods, since Theorem 14 gives convergence bounds on u¯ , namely |˜a| either side of p∗ , rather than simply quantifying it as O(|˜a|). The requirement on κ being small is similar to other schemes in the literature; analysing stability for larger κ would require more knowledge of the system than has been assumed in this paper. Note that even if the bound on ∥Ψ (θ )∥ is arbitrarily large, as determined by the initial conditions of (1), a sufficiently small κ always exists to guarantee stability, since all of the terms containing |h′ (ξ2 )| are multiplied by κ .
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
3005
4. Example Two examples will be presented to demonstrate the efficacy of the ES scheme developed in this paper. The first example shows the behaviour of the controller on a simulated Wiener–Hammerstein plant, for various initial conditions. An experimental implementation of the scheme is then detailed, illustrating its ability to increase the mean maximum brake torque of a gas fuelled spark-ignition engine operating at constant speed, by adjusting the spark advance angle. 4.1. Simulated system Consider the Wiener–Hammerstein plant described by f (p) =
AG CG
1
1
3
(p − 5)2 − sin(p − 5) + (p − 5) 4 2 2 1 1 −2 1 2 BG 1 1 = 12 , 4 2 DG
AH
BH
CH
DH
5 4
=
1 2 1 3
0 1 1 3
0
(50a) Fig. 4. ES scheme input trajectories for different initial conditions.
(50b)
.
The static nonlinear map in this case is asymmetric about the location of the extremum p∗ = 5 (f (5) = 0). Selecting a = 1, κ ∈ {0.2, 0.6} and σ = −1, the ES scheme (1) is applied to this system with several initial values of u¯ . For each of these initial values, xG and xH are initialised at their corresponding steady-state values. Figs. 4–5 illustrate the resulting trajectories of u¯ and y respectively. As can be seen, the equilibrium of p0 ≈ 4.9681 is close to the true optimising value. In this case, the effective dither amplitude is given by a˜ ≈ −0.8824, where a phase shift has been introduced by G(·). For κ = 0.6, the dynamics are seen to initially generate transient behaviour which subsides as the system approaches steady state. Convergence is achieved after approximately 60 steps to qˆ 0 (p0 ) ≈ 0.0013, which is very close to the actual minimum of zero. It is evident that the transient dynamics are not excited by using a smaller κ , leading to smoother trajectories, at the expense of slower convergence. Further insight into convergence characteristics can be seen in Fig. 6, which illustrates the input–output behaviour of (50) in the closed loop, together with the steady-state map (50a). It can be seen that for initial conditions far from the extremum, there is an initial ‘‘fast’’ convergence behaviour exhibited with κ = 0.6, before the trajectories begin to approximately follow the path of steepest descent. The κ = 0.2 case exhibits more conventional steepest descent trajectories.
Fig. 5. ES scheme output trajectories for different initial conditions.
4.2. Experimental results Spark-ignition engines provide a practical platform upon which the ES scheme described in this paper can be experimentally tested. For a given fuel composition, an engine operating at a constant speed condition has a mean steady-state spark advance to maximum brake torque mapping which is approximately quadratic (Mohammadi et al., 2012). In engines using alternative fuels such as natural gas, LPG and ethanol-gasoline blends, the fuel in the tank is subject to compositional variation due to cost, geographical region and availability of different fuels, amongst other factors. It is therefore desirable to vary the spark advance in response to changing composition to ensure that the engine is operating at peak efficiency. The ES scheme (1) is implemented to vary the spark angle on a modified 6-cylinder Ford Falcon BF MY2006 engine coupled to
Fig. 6. Plant output vs. input for different initial conditions.
an eddy-current dynamometer, at an initial operating condition of 800 RPM and 32 N.m torque output. An ATI VISION system is used to modify the spark advance commanded by the ECU, measured in terms of Crank Angle Degrees before Top Dead Centre (CAD BTDC). By considering only mean-value engine behaviour (Guzzella & Onder, 2004), nonlinear dynamics arising from combustion and cyclic variations can be neglected, for an appropriately selected sampling rate. The overall engine and dynamometer system
3006
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Acknowledgements The authors would like to thank the Advanced Centre for Automotive Research and Testing (ACART) for the use of their dynamometer facility, as well as Alireza Mohammadi for his assistance with running the experiments. Appendix A. Proof of Lemma 10 For the first property, twice differentiating (27) and applying the mean value theorem gives 1
(f ′ (˜p + p0 + a˜ ) − f ′ (˜p + p0 − a˜ )) 2a˜ = f ′′ (c1 ),
g ′′ (˜u) =
Fig. 7. Convergence of (top) engine input; and (bottom) measured output towards their respective optima following implementation of the proposed scheme to online spark calibration of a CNG engine.
can then be adequately represented by a Hammerstein model, with output dynamics arising from the torque sensor on the dynamometer. Fig. 7 shows how the output torque of the engine can be optimised by modifying the spark advance angle, for running on pure methane. A sampling rate of 0.2 Hz has been used for updating the control action. From an initial spark advance angle of 26 °, which is not the optimal angle for methane fuel, the results show an improvement of approximately 3 % in mean maximum brake torque from the initial condition by adjusting the spark advance angle to just above 34 °. The effect of the discrete dither signal is clearly seen in the plot of spark angle. Note the noisy torque signal arises from cyclic variations in engine operation. Its effect on the ES scheme is minimised by using the average torque measurement over each sample. Remark 23. Using the alternating dither signal for spark timing modification has no detrimental effect on the system, since the spark plugs are driven electronically and contain no moving parts. If mechanical wear is of concern for other applications, the dither amplitude can be reduced to a suitably small value. In addition, the continuous-time realisation of the dither signal need not necessarily be a square wave obtained by applying a zeroorder hold. Other choices are possible, but will lead to different discretisations of the dynamics.
for some c1 ∈ (˜p + p0 − |˜a|, p˜ + p0 + |˜a|). Hence, g ′′ (·) must satisfy the same bounds as f ′′ (·). Furthermore, since a˜ g ′ (˜p) = qˆ 1 (˜p + p0 ), the second property follows directly from the first property and the equilibrium solution defined in Theorem 8. Appendix B. Proof of Lemma 11 Twice differentiating (28) and applying the intermediate value theorem shows that 1 h′′ (˜p) = (f ′′ (˜p + p0 + a˜ ) + f ′′ (˜p + p0 − a˜ )) 2 = f ′′ (c2 ), for some c2 ∈ (˜p + p0 − |˜a|, p˜ + p0 + |˜a|). The result follows directly as in Lemma 10. Appendix C. Proof of Lemma 13 Differentiating (28) and using a Taylor series expansion for each term with the Lagrangian form for the remainder, it can be shown that a˜ h′ (˜p) = f ′ (˜p + p0 ) + (f ′′ (c1 ) − f ′′ (c2 )), 2 for some c1 and c2 . Applying the triangle inequality and the bounds on f ′′ (·) from Assumption 1 gives the desired result. Appendix D. Proof of Lemma 17 Expanding (44), it is required that
5. Conclusions and future work This paper has presented a new discrete-time extremumseeking algorithm for discrete-time Wiener–Hammerstein plants with stable linear dynamics. It has shown how limited plant knowledge can be used to guarantee semi-global stability for appropriate choices of tuning parameters, assuming the system possesses a strongly-convex steady-state mapping. A practitioner can utilise this theoretical result for tuning the ES controller by choosing the dither amplitude based upon the desired proximity to the true optimum and then selecting the adaptation gain to be sufficiently small. The algorithm has been demonstrated in simulation as well as experimentally on a spark-ignition engine. Future work will analyse Wiener–Hammerstein-like systems with a Multi-Input Single-Output (MISO) nonlinear maps. It will also investigate whether upper bounds on the maximum adaptation gain for stability can be provided, as in Shekhar et al. (2013). An additional avenue of investigation will be the incorporation of input and output constraints into the formulation.
γ R(P11 , A11 ) −γ AT12 P11 A11
⋆ ≻ 0. T R(P22 , A22 ) − γ A12 P11 A12
Using the Schur complement, this is equivalent to the condition
R(P22 , A22 ) − γ AT12 Q(P11 , A11 )A12 ≻ 0,
(D.1)
where
Q(P , A) := P + PAR(P , A)−1 AT P . Applying the matrix inversion lemma and the Schur complement again, it follows that −1 Q(P11 , A11 ) = R(P11 , AT11 )−1 ≻ 0,
so (D.1) can be written in the form −1 R(P22 , A22 ) − γ AT12 R(P11 , AT11 )−1 A12 ≻ 0,
which can be satisfied by choosing γ such that 0
Contents lists available at ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
Discrete-time extremum-seeking for Wiener–Hammerstein plants✩ Rohan C. Shekhar 1 , William H. Moase, Chris Manzie Department of Mechanical Engineering, The University of Melbourne, Victoria 3010, Australia
article
info
Article history: Received 24 September 2013 Received in revised form 23 March 2014 Accepted 16 July 2014 Available online 27 October 2014 Keywords: Extremum-seeking Adaptive control Optimisation Engine control
abstract This paper presents a new formulation of extremum-seeking control for discrete-time Wiener– Hammerstein plants. A novel method of analysis using Linear Parameter-Varying (LPV) system theory demonstrates semi-global stability of the control scheme. Assuming only limited plant knowledge, the stability result ensures convergence of the plant output in steady state to a point in an arbitrarily small neighbourhood of the extremum, for appropriately chosen controller parameters. The behaviour of the control scheme is analysed on a simple simulated system, prior to being implemented on an internal combustion engine. Experiments demonstrate how the scheme is able to maximise engine output torque in the presence of an uncertain fuel composition by modifying the spark timing. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Extremum-Seeking (ES) is a closed-loop control paradigm that seeks to optimise the steady-state output of a given nonlinear plant, where a plant model is unavailable or only approximately known. Assuming that the input-to-output map at steady state contains an extremum (maximum or minimum) in the region of interest, ES manipulates the plant input so that the output is contained within a small neighbourhood of this extremum at equilibrium. For a small cross-section of different ES schemes, see Chichka, Speyer, and Park (1999), Choi, Krstić, Ariyur, and Lee (2002), Drakunov, Ozguner, Dix, and Ashrafi (1995), Guay and Zhang (2003), Moase, Manzie, and Brear (2010) and Zhang and Ordonez (2007). A comprehensive review of the evolution of ES schemes throughout history can be found in Tan, Moase, Manzie, Nešić, and Mareels (2010). Most popular ES schemes emulate some form of gradient descent, using a dither signal to estimate local derivative information and an adaptation gain to regulate the speed at which the plant input approaches the steady-state optimising value (Ariyur & Krstić, 2003). The first analyses of local stability were performed for continuous-time ES schemes with nonlinear dynamics (Krstić
✩ This work was supported under the Australian Research Council’s Discovery Projects funding scheme (Project DP120101830). The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Raul Ordóñez under the direction of Editor Miroslav Krstic. E-mail addresses: [email protected] (R.C. Shekhar), [email protected] (W.H. Moase), [email protected] (C. Manzie). 1 Tel.: +61 3 83447714; fax: +61 3 8344 4290.
http://dx.doi.org/10.1016/j.automatica.2014.10.030 0005-1098/© 2014 Elsevier Ltd. All rights reserved.
& Wang, 2000). Non-local, or semi-global stability results were subsequently derived in Tan, Nešić, and Mareels (2006). Analysis of these ES schemes generally involves singular perturbation analysis and averaging theory; however, recent approaches have also used Lie bracketing techniques (Dürr, Stanković, & Johansson, 2011). Grey-box ES schemes using limited plant knowledge have also recently been shown to guarantee semi-global stability (Nešić, Mohammadi, & Manzie, 2013; Sharafi, Moase, Shekhar, & Manzie, 2013). Semi-globally stable ES for continuous-time Hammerstein and Wiener–Hammerstein plants was considered in Moase and Manzie (2012a,b) respectively. Such plants consist of a static nonlinear element with either output LTI dynamics (Hammerstein) or both input and output LTI dynamics (Wiener–Hammerstein). These plant models arise in the black-box identification of nonlinear dynamic systems (Nešić, 1999) and have proven useful for representing biological systems (Hunter & Korenberg, 1986), diesel engines (Perez, Blasco, García-Nieto, & Sanchis, 2006) and chemical processes (Eskinat, Johnson, & Luyben, 1991). Various discrete-time ES schemes have also been presented in the literature. Specifically, Teel and Popović (2001) outline an approach that waits for a sufficient period of time until system dynamics have settled before applying subsequent control actions, effectively treating the plant as a static nonlinearity. Discrete dynamics are considered in Choi et al. (2002), which demonstrates local stability for a Wiener–Hammerstein plant with a quadratic nonlinearity. In contrast to these methods, a dither-free sampleddata approach is considered in Khong, Manzie, Nešić, and Tan (2012). It should be noted that all of these approaches rely on explicit time-scale separation arguments so that the nonlinear
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
system approximates a static map; this can limit convergence speed. It is possible to somewhat relax time-scale separation requirements for ES on Wiener–Hammerstein plants in continuous time, as demonstrated by Moase and Manzie (2012a). In discrete time, Shekhar, Moase, and Manzie (2013) laid the foundation for achieving the same goal, using a new method of stability analysis for ES on a static nonlinear plant. Specifically, Linear Difference Inclusion (LDI) theory and Lyapunov arguments are used to show semiglobal stability for sufficiently small adaptation gains, assuming a strongly-convex plant map. This paper builds upon the work in Shekhar et al. (2013), presenting a new formulation of ES for discrete-time Wiener– Hammerstein plants with nonlinear maps satisfying the same convexity requirement. Extending the techniques developed in Shekhar et al. (2013), Linear Parameter-Varying (LPV) system theory is applied to prove semi-global stability, which improves upon the local stability results of Choi et al. (2002). The only prior knowledge required is the sign of the system’s dynamic response evaluated at the dither frequency, which can be experimentally determined for real systems; this is less restrictive than the knowledge requirement for grey-box schemes. No explicit timescale separation between the plant dynamics and gradient estimation is required, potentially improving convergence speed as demonstrated in simulation. The scheme presented in this paper also exhibits convergence of the plant output to a point in an explicitly characterised neighbourhood of the extremum, in contrast to most existing ES approaches that result in limit-cycle equilibria. An experimental example is used to illustrate the performance of the ES scheme on a spark-ignition engine. In particular, the ES scheme is used to maximise the mean torque output at steady state by modifying the spark advance angle. Different ES schemes have already been shown to improve efficiency through modification of engine tuning parameters (Mohammadi, Manzie, & Nešić, 2012; Popovic, Jankovic, Magner, & Teel, 2003), making this an ideal candidate application to demonstrate the efficacy of the scheme developed in this paper. The scheme is shown to improve torque output, motivating its use for closed-loop tuning of alternative fuel engines where the fuel composition may vary. 1.1. Nomenclature The sets of natural numbers, integers and real numbers are denoted as N, Z and R respectively. The set of all non-negative integers is denoted as Z≥0 . The column vector [aT , bT , c T ]T is written as [a; b; c ]. In represents the n × n identity matrix. The operator ⊕ represents the direct sum of matrices. The interior of a set is indicated by int(·). For some P ≻ 0 and a square matrix A of the same dimension, define the discrete Lyapunov operator by
R(P , A) := P − AT PA. The operator λi (M ) denotes the ith eigenvalue of a matrix M and vi (M ) denotes the corresponding eigenvector, where the eigenvalues are enumerated in some specified order. λ(M ) and λ(M ) denote the maximum and minimum eigenvalues of a symmetric matrix M respectively. The 2-norm of a matrix or vector A is written as ∥A∥. The ⋆ symbol is used in symmetric block matrices to denote the transpose of the diagonally opposite block. To illustrate this, given compatible matrices A = AT , C = C T and B,
A B
⋆ C
:=
A B
BT . C
2. Problem formulation This paper is concerned with optimising the steady-state output of the Wiener–Hammerstein plant shown in Fig. 1. The plant is
2999
Fig. 1. Discrete-time Wiener–Hammerstein plant.
described by a static nonlinear map f : R → R, with linear input dynamics G(z ) and output dynamics H (z ). The dynamics are assumed to have minimal state-space realisations defined by the equations xG (k + 1) = AG xG (k) + BG u(k) p(k) = CG xG (k) + DG u(k) xH (k + 1) = AH xH (k) + BH q(k) y(k) = CH xH (k) + DH q(k), for states xG and xH of the input and output dynamics respectively. Whilst explicit forms for the mapping and LTI systems are unknown, they are assumed to satisfy the following properties: Assumption 1. There exist positive constants α, β such that α ≤ f ′′ (p) ≤ β . This implies that the nonlinear map f (·) is twice differentiable and strongly convex. It also guarantees the existence of a point p∗ ∈ R such that f ′ (p∗ ) = 0. Assumption 2. The matrices AG and AH have spectral radii strictly less than unity. Assumption 3. Without loss of generality, the LTI dynamics satisfy G(1) = H (1) = 1, so that any steady-state gain is considered as part of f (·). Remark 4. Assumption 1 guarantees that p∗ is the unique minimiser of f (·). Assumption 2 ensures that G(·) and H (·) are stable. Remark 5. Assumption 1 is less restrictive than assuming a quadratic steady-state map, which is used extensively in the literature for local stability analysis (e.g. Choi et al., 2002). Indeed, it will be shown that α can be arbitrarily small and β can be arbitrarily large whilst still ensuring closed-loop stability. 2.1. Controller design To find the steady-state input that minimises f (·), the control scheme in Fig. 2 will be used. The input u¯ (k) represents an estimate for the minimiser of f (·) at time step k. A simple dither signal a(−1)k is added to this estimate and applied to the Wiener–Hammerstein plant, where a > 0 is the amplitude. The resulting output is then ‘‘demodulated’’ by the signal (−1)k , before being passed through a two-step FIR averaging filter and scaled by an adaptation gain κ > 0 having sign σ ∈ {−1, 1}. This quantity, which is approximately proportional to f ′ (¯u(k)), is then subtracted from u¯ (k), essentially performing discrete-time integration. The sign σ is used to correct for any phase shift induced by G(·) and H (·). The closed-loop dynamics of Fig. 2 are described by u(k) = u¯ (k) + a(−1)k
(1a)
xG (k + 1) = AG xG (k) + BG u(k)
(1b)
p(k) = CG xG (k) + DG u(k)
(1c)
q(k) = f (p(k))
(1d)
xH (k + 1) = AH x(k) + BH q(k)
(1e)
y(k) = CH x(k) + DH q(k)
(1f)
v(k) = y(k)(−1)
(1g)
k
σκ
(v(k − 1) + v(k)), (1h) 2 where the initial conditions u¯ (0), xG (0), xH (0) and v(−1) are arbitrary.
u¯ (k + 1) = u¯ (k) −
3000
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Fig. 2. ES scheme for Wiener–Hammerstein plant.
Fig. 3. Loop transformation.
Remark 6. This ES scheme is one of many approaches to find the extremum of a Wiener–Hammerstein plant. More traditional schemes would use an integrator after the demodulator, but it will be seen that the averaging filter provides desirable properties. The analysis presented in this paper provides a general methodology that could be easily applied to other discrete-time ES schemes with different kinds of filters. Remark 7. The ES scheme and analysis to follow can be trivially modified to find the steady-state maximum when f (·) is a strongly concave function, by replacing κ with −κ .
˘ (z ) = G(z ). G
(8)
The modified dither signal is then given by a˜ = G(−1)a.
(9)
Note that this may be negative depending on any phase inversion introduced by G(·). For the output dynamics, (1e) and (1f) can be written in the form xH (k + 1)(−1)k+1 = −AH xH (k)(−1)k − BH q(k)(−1)k y(k)(−1)k = CH xH (k)(−1)k + DH q(k)(−1)k = v(k).
3. Analysis of stability Analysis of (1) will be carried out in four steps. Firstly, a loop transformation will be applied to Fig. 2 in order to shift the dynamics into the gradient estimator portion of the closed loop (Section 3.1). Following this, Fourier analysis will be utilised in order to identify an equilibrium solution (Section 3.2). In terms of this solution, various state transformations will then be applied to write the dynamics in terms of an error system having an equilibrium at the origin (Section 3.3). Finally, the error system dynamics will be cast into the form of an LPV system (Section 3.4) that will be analysed with Lyapunov techniques (Section 3.5). To simplify notation, the input-to-state transfer functions
ΓG (z ) := (zI − AG ) BG
(2)
ΓH (z ) := (zI − AH )−1 BH
(3)
−1
which demonstrates that
By defining x˘ H (k) := xH (k)(−1)k
(10)
q˘ (k) := q(k)(−1)k ,
˘ (z ) can be written in state-space form as the system H x˘ H (k + 1) = −AH x˘ H (k) − BH q˘ (k)
(11)
v(k) = CH x˘ H (k) + DH q˘ (k),
(12)
which shows that
˘ (z ) = H (−z ). H
(13)
This completes the transformation. 3.2. Identifying an equilibrium solution
will be used extensively, together with the identities BG − z ΓG (z ) = −AG ΓG (z )
(4)
BH − z ΓH (z ) = −AH ΓH (z ),
(5)
To find an equilibrium solution for u¯ , note that from Fig. 3, q(k) = f (˘p(k) + a˜ (−1)k ).
which can be easily verified by rearranging (2) and (3).
The Discrete-Time Fourier Series (DTFS) can now be used to write
3.1. Loop transformation
q(k) = qˆ 0 (˘p(k)) + qˆ 1 (˘p(k))(−1)k , where the DTFS coefficients as functions of p˘ are given by
Analysis of the closed-loop system (1) can be simplified by swapping the order of the G(·) and the summation block, as well as H (·) and the product block, as shown in Fig. 3. From (1a) and (1b),
qˆ 0 (˘p) :=
xG (k + 1) = AG xG (k) + BG (¯u(k) + a(−1)k ).
qˆ 1 (˘p) :=
Now define the state x˘ G (k) := xG (k) − aΓG (−1)(−1)k , which is the error between xG (k) and its steady-state value to the excitation caused by the dither signal. In terms of this new state, ˘ (·) can be shown to have a state-space realisation G x˘ G (k + 1) = AG x˘ G (k) + BG u¯ (k)
(6)
p˘ (k) = CG x˘ G (k) + DG u¯ (k),
(7)
1 2 1 2
(f (˘p + a˜ ) + f (˘p − a˜ ))
(14a)
(f (˘p + a˜ ) − f (˘p − a˜ )).
(14b)
For a fixed p˘ , this harmonically decomposes q into a mean component and an oscillating component. In terms of these components, q˘ (k) = qˆ 1 (˘p(k)) + qˆ 0 (˘p(k))(−1)k ,
(15)
which has essentially swapped the harmonic components of q. Theorem 8 exploits the harmonic decomposition to identify an equilibrium solution to the closed-loop system (1).
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
forcing signal q˘ (k) (15). Using x¯ H (k, p˘ ) to denote the quasi-steady solution,
Theorem 8. There exists a unique point p ∈ (p − |˜a|, p + |˜a|) 0
∗
3001
∗
x¯ H (k, p˘ ) := ΓH (−1)ˆq1 (˘p) + ΓH (1)ˆq0 (˘p)(−1)k .
such that
This is then used to define the error state
u¯ (k) = p0
(16a)
xG (k) = ΓG (1)p0 + aΓG (−1)(−1)k
x˜ H (k) := x˘ H (k) − x¯ H (k, p˘ (k)).
(16b)
xH (k) = ΓH (1)ˆq0 (p0 )
(16c)
y(k) = qˆ 0 (p0 )
(16d)
From (10) and (16c), it can be seen that the equilibrium value for x˜ H is zero if p˘ (k) = p˘ (k − 1) = p0 . Using (5) and (11), the dynamics of this error state can be written in the form
is an equilibrium solution to (1).
x˜ H (k + 1) = x˘ H (k + 1) − x¯ H (k + 1, p˘ (k + 1))
= −AH x˜ H (k) − ΓH (−1)1qˆ 1 (˘p(k + 1))
Proof. Applying the mean value theorem to (14b) and invoking Assumption 1, a˜ −1 qˆ ′1 (˘p) =
1
1qˆ n (˘p(k)) := qˆ n (˘p(k)) − qˆ n (˘p(k − 1)),
where c ∈ (˘p − |˜a|, p˘ + |˜a|). This shows that a q1 (·) is strictly monotonically increasing. Then, since p∗ is the global minimiser of f (·), it can be seen that
˜ −1 ˆ
a˜ −1 qˆ 1 (p∗ − |˜a|) = a˜
(21)
where
(f ′ (˘p + a˜ ) − f ′ (˘p − a˜ ))
2a˜ = f ′′ (c ) > 0,
−1
+ ΓH (1)1qˆ 0 (˘p(k + 1))(−1)k ,
qˆ 1 (p + |˜a|) = ∗
1 2|˜a| 1 2|˜a|
(f (p∗ ) − f (p∗ − 2|˜a|)) < 0 (f (p + 2|˜a|) − f (p )) > 0. ∗
∗
qˆ 1 (p0 ) = 0.
for n ∈ {0, 1}. Eq. (12) can then be written in the form
v(k) = CH x˜ H (k) + H (−1)ˆq1 (˘p(k)) + qˆ 0 (˘p(k))(−1)k ,
(17)
Furthermore, this point is unique by monotonicity. Having established the existence of a unique p0 , (16) can be verified by substitution into (1), applying Assumption 3 and (17). In particular, it can be seen that the input to the FIR filter is given by
v(k) = qˆ 0 (p0 )(−1)k . Since this is completely removed by the filter, u¯ remains unchanged with value p0 .
(23)
which follows from (13) and Assumption 3. Finally, by defining the translations u˜ := u¯ − p0
(24)
0
(25)
p˜ := p˘ − p
Hence, by the intermediate value theorem, there exists a point p0 ∈ (p∗ − |˜a|, p∗ + |˜a|) such that
(22)
and substituting these into (19), (20), (21), (23), and (1h), the dynamics of the error system are given by x˜ G (k + 1) = AG x˜ G (k) − ΓG (1)1u˜ (k + 1)
(26a)
p˜ (k) = CG x˜ G (k) + u˜ (k)
(26b)
x˜ H (k + 1) = −AH x˜ H (k) − ΓH (−1)1qˆ 1 (˜p(k + 1) + p ) 0
+ ΓH (1)1qˆ 0 (˜p(k + 1) + p0 )(−1)k v(k) = CH x˜ H (k) + H (−1)ˆq1 (˜p(k) + p ) + qˆ 0 (˜p(k) + p0 )(−1)k σκ u˜ (k + 1) = u˜ (k) − (v(k − 1) + v(k)) ,
(26c)
0
2
(26d) (26e)
where 1u˜ (k) := u˜ (k) − u˜ (k − 1). 3.3. Error system formulation For a constant or slowly-varying input u¯ , (2) and (8) allow the ˘ (·) to be written as quasi-steady state of G x¯ G (¯u) := ΓG (1)¯u.
3.4. LPV system formulation
Then define the new state x˜ G (k) := x˘ G (k) − x¯ G (¯u(k)).
(18)
Using (4) and (6), it can be seen that x˜ G (k + 1) = x˘ G (k + 1) − x¯ G (¯u(k + 1))
= AG x˜ G (k) − ΓG (1)1u¯ (k + 1),
The error system will be further transformed into LPV form. In order to accomplish this, define the functions g (˜p) :=
(19)
where 1u¯ (k) := u¯ (k) − u¯ (k − 1). In addition, (7), (8), (18) and Assumption 3 can be combined to give p˘ (k) = CG (˜xG (k) + x¯ G (¯u(k))) + DG u¯ (k)
= CG x˜ G (k) + u¯ (k).
Remark 9. These equations describe the closed-loop dynamics (1) exactly, without any approximation. The quasi-steady solutions have only been used to formulate an error system for facilitating stability analysis.
(20)
Now, observe that the transfer function from q˘ to x˘ H is given by ΓH (−z ), from (3) and (13). Then, if p˘ is fixed or slowly varying with time, the quasi-steady solution for x˘ H can be evaluated as the frequency response of ΓH (−z ) to the approximately periodic
1 a˜
p˜
qˆ 1 (µ + p0 )dµ
(27)
0
h(˜p) := qˆ 0 (˜p + p0 ).
(28)
The following lemmas describe some important properties of these functions. Lemma 10. For all p˜ ∈ R, g (·) satisfies the bounds
α ≤ g ′′ (˜p) ≤ β. It also has a unique minimum at the origin. Proof. See Appendix A.
3002
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Lemma 11. For all p˜ ∈ R, h(·) satisfies the bounds
where
α ≤ h (˜p) ≤ β.
φ(θ (k)) := a˜ H (−1)g ′′ (ξ1 (k)) + h′ (ξ2 (k))(−1)k .
′′
Proof. See Appendix B.
To complete the transformation, (26b), (26d) and (37) are substituted into (26e) to give
σκ
Remark 12. These results demonstrate that g (˜p) and h(˜p) are closely related to f (˜p + p0 ), being strongly convex and satisfying the same second derivative bounds as f (·).
1u˜ (k + 1) =
Lemma 13. The function h(·) satisfies the bounds
which uses the fact that
2
(2a˜ H (−1)ζ (k)(CG x˜ G (k) + u˜ (k))
− w(k) + CH x˜ H (k)),
(40)
from (29) and (33). Combining (26a), (26e) and (34)–(40), the LPV system description is given by
In terms of (25), (27) and (28), the DTFS coefficients (14) are given by qˆ 1 (˘p(k)) = a˜ g ′ (˜p(k))
(29)
qˆ 0 (˘p(k)) = h(˜p(k)).
(30)
To allow the LPV system to be defined, it is useful to apply the mean value theorem to (22), using (26a), (26b), (29) and (30) to write
1qˆ 1 (˘p(k + 1)) = a˜ g ′′ (ξ1 (k))1p˜ (k + 1) 1qˆ 0 (˘p(k + 1)) = h′ (ξ2 (k))1p˜ (k + 1),
X (k + 1) = A(θ (k), κ)X (k),
1p˜ (k + 1) := p˜ (k + 1) − p˜ (k) = CG (AG − I )˜xG (k) + DG 1u˜ (k + 1)
(31)
ξ1 (k), ξ2 (k) ∈ (min{˜u(k), u˜ (k + 1)}, max{˜u(k), u˜ (k + 1)}).
(32)
Also define p˜ (k) ̸= 0
(33)
p˜ (k) = 0.
(41)
where X (k) := [˜u(k); w(k); x˜ H (k); x˜ G (k)] A(θ , κ) := A0 (θ ) −
1 0 A0 (θ ) := 0 0
0 0 0 0
σκ ¯ A(θ )
2 0 −C H −AH 0
0 φ(θ )CG (AG − I ) Φ (θ )CG (AG − I ) −AG
T
2aG(−1)H (−1)ζ (k) DG φ(θ ) −1 ¯A(θ ) := D Φ (θ ) . CH G −ΓG (1) 2aG(−1)H (−1)ζ (k)CG
where
′ g (˜p(k)) ζ (k) := p˜ (k) ′′ g (0)
2
qˆ 1 (˜p(k) + p0 ) = 2a˜ ζ (k)˜p(k),
1
|h′ (˜p)| ≤ |f ′ (˜p + p0 )| + a˜ (β − α). Proof. See Appendix C.
(39)
1
Note that zeros in the matrices may represent blocks of zeros having appropriate dimensions. 3.5. LPV Lyapunov analysis This section will prove the following convergence result, which is the primary contribution of this paper.
It can be easily shown that g ′ (·) is sector bounded, so
α ≤ ζ (k) ≤ β
Theorem 14 (ES Scheme Convergence). As k → ∞,
for all k. In terms of (27), (28) and (33), define the parameter vector
u¯ (k) → p0 ∈ (p∗ − |˜a|, p∗ + |˜a|)
θ(k) := [θ1 (k); θ2 (k); θ3 (k); θ4 (k); θ5 (k)] := a˜ ζ (k); a˜ g ′′ (ξ1 (k)); h′ (ξ2 (k))(−1)k ;
exponentially and y(k) → qˆ 0 (p0 ), for any initial conditions of (1), any positive α, β ∈ Z≥0 , α < β and any a˜ ∈ Z, provided κ is sufficiently small and
a˜ ζ (k)g ′′ (ξ1 (k)); a˜ ζ (k)h′′ (ξ2 (k))(−1)k .
(34)
These definitions can then be used to write (26c) in the form x˜ H (k + 1) = −AH x˜ H (k) + Φ (θ (k))(CG (AG − I )˜xG (k)
+ DG 1u˜ (k + 1)),
(35)
σ G(−1)H (−1) > 0. Proof. The intermediate results required to prove the theorem are detailed in Section 3.5.1, while the proof of the theorem itself is provided in Section 3.5.2.
where
Φ (θ(k)) := −˜ag ′′ (ξ1 (k))ΓH (−1) + h′ (ξ2 (k))(−1)k ΓH (1).
(36)
Now define the new state
w(k) := −v(k − 1) + H (−1)ˆq1 (˜p(k) + p0 ) −ˆq0 (˜p(k) + p0 )(−1)k .
(37)
Using (26d) and (31), it follows that
w(k + 1) = −CH x˜ H (k) + φ(θ (k))(CG (AG − I )˜xG (k) + DG 1u˜ (k + 1)),
(38)
Remark 15. Theorem 14 presents a semi-global stability result, since it states that for any initial condition, there is always a κ small enough to ensure convergence. The requirement σ G(−1)H (−1) > 0 ensures that any inversion of the dither signal induced by G(·) and H (·) is compensated by σ . It also implies that both LTI systems cannot have zeros at z = −1, which would remove the dither signal completely in steady state. Note also that the theorem does not require |˜a| to be small, unlike most other ES schemes in the literature. Hence, a can be chosen based on the desired proximity to the true optimum, but large values may necessitate the use of a smaller κ .
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
3003
3.5.1. Supporting results
Proof. See Appendix F.
Lemma 16. If for some κ0 there exists a matrix P that satisfies the parameter-dependent LMI condition
Since G(·) and H (·) are stable, there exist positive-definite matrices PG and PH such that
P
PA(θ , κ0 )
⋆
P
R(PG , AG ) ≻ 0
≻ 0,
R(PH , AH ) ≻ 0.
for all θ contained in some parameter set Θ , then V (X ) := X PX T
(42)
is a Lyapunov function for the LPV system (41) on a subset of the state space X corresponding to θ ∈ Θ , for κ = κ0 . Proof. This result follows from applying the Schur complement to the Discrete Lyapunov Equation (DLE), P − AT (θ , κ0 )PA(θ , κ0 ) = R(P , A(θ , κ0 )) ≻ 0.
(43)
From elementary Lyapunov theory, if there exist P ≻ 0 and κ0 > 0 satisfying (43) for all θ ∈ Θ , then (42) is a Lyapunov function for the LPV system on X. Lemma 17 (Block Upper-Triangular System Stability). Given a block upper-triangular matrix with square diagonal entries
A :=
A11 0
A12 , A22
if there exist P11 ≻ 0 and P22 ≻ 0 such that
Using these matrices, stability of the LPV system can now be shown for suitably chosen κ . Theorem 20 (Existence of Lyapunov Matrix). If ∥Ψ (θ )∥ is bounded for some θ and σ G(−1)H (−1) > 0, then
R(P , A(θ , κ)) ≻ 0, where P := 1 ⊕ δγ ⊕ δ PH ⊕ PG ,
Proof. The proof proceeds by considering the case κ = 0, and showing that as κ is increased by a small amount, the DLE (47) is satisfied. It utilises the fact that 1 0 A(θ , 0) = A0 (θ ) = 0 0
0
A1 :=
0 φ(θ )CG (AG − I ) Φ (θ )CG (AG − I ) −AG
−CH , −A H
0 0
P1 := γ ⊕ PG ,
R(P , A) ≻ 0,
(44)
for some sufficiently small γ > 0. Lemma 17 then guarantees that R(P1 , A1 ) ≻ 0. Consider now the larger submatrix
for all −1 , AT11 )) ∥A12 ∥−2 . 0 < γ < λ(R(P22 , A22 ))λ(R(P11
A2 (θ ) := (45)
Proof. See Appendix D.
Lemma 18 (Eigenvalue Derivatives). Consider a symmetric matrixvalued function M : R → RN , assumed to be differentiable on some interval D ⊆ R. The derivative of any distinct eigenvalue λi (M (κ)) with respect to κ for κ ∈ int(D ) is given by the expression
= vTi (M (κ))
Proof. See Appendix E.
dκ
A1 0
=
dM (κ)
where (46) has been substituted. Since ∥Ψ (θ )∥ is bounded, it follows that the quantity
∥Ψ (θ )CG (AG − I )∥ is also bounded, so Lemma 17 can be applied again to show that
= δγ ⊕ δ PH ⊕ PG satisfies R(P2 , A2 (θ )) ≻ 0, for a sufficiently small δ > 0. The overall candidate Lyapunov matrix
Lemma 19. If |˜u| ≤ L1 and ∥xG ∥ ≤ L2 , then
a˜ H (−1) ∥Ψ (θ)∥ ≤ −ΓH (−1) β H (1) 1 + β(L1 + ∥CG ∥ L2 ) + a˜ (β − α) , ΓH (−1) 2
P := 1 ⊕ P2
where
M (θ , κ) :=
is a matrix comprising the parameters defined in (36) and (39).
φ(θ )CG (AG − I ) Φ (θ )CG (AG − I ) 0 −AG Ψ (θ )CG (AG − I ) , −AG −CH −A H
P2 := δ P1 ⊕ PH
vi (M (κ)).
Ψ (θ) := [φ(θ ); Φ (θ )]
0 0 0
Furthermore, such a γ always exists.
dκ
0 −CH −A H 0
define
satisfies the DLE
dλi (M (κ))
0 0 0 0
the matrix 0 P22
(48)
for any positive α, β and sufficiently small κ, γ and δ .
R(P22 , A22 ) ≻ 0,
γ P11
(47)
has a block upper-triangular structure. Considering the submatrix
R(P11 , A11 ) ≻ 0
P :=
= 1 ⊕ δγ ⊕ δ PH ⊕ PG
(49)
is now defined in terms of γ , δ, PG and PH . Applying the LMI stability condition of Lemma 16, define the matrix function
(46)
PA(θ , κ)
=
⋆
P
P
PA0 (θ )
P
⋆ P
−
σκ 2
0 P A¯ (θ )
⋆ 0
.
3004
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Satisfaction of (47) is then equivalent to M (θ , κ) being positive definite. Given the block diagonal structure of P, it is useful to permute the blocks of M (·, ·). Note that for any invertible matrix S,
corresponding to the initial states of the original system. Then, define a level set of the quadratic form X T PX passing through this state by the ellipsoid
M (θ, κ) ≻ 0 ⇐⇒ S −1 M (θ , κ)S ≻ 0.
˜ (0)}, E := {X |X T PX ≤ X T (0)PX
This is effectively a similarity transformation, which preserves the eigenvalues of M (θ , κ) for all parameter values. Since P (49) is expressed as a 2 × 2 block matrix, M (·, ·) can be expressed in the form of a 4 × 4 block matrix. Choosing S to be the (orthogonal) permutation matrix that reorders the block columns {1, 2, 3, 4} to {1, 3, 2, 4}, the new matrix
where P is chosen as in (48), for some γ and δ , and P˜ := 1 ⊕ 1 ⊕ PH ⊕ PG ,
˜ (θ, κ) := S T M (θ , κ)S M
which is the matrix P with γ = δ = 1. Any positive value can ˜ but they are set to unity for ease of be assigned to γ and δ in P, analysis. Given the block diagonal structure of P, it is useful to define E in the equivalent form
can be defined. For κ = 0, this matrix has the form
E = {(˜u, w, x˜ H , x˜ G ) | u˜ 2 + δγ w 2 + δ˜xTH PH x˜ H
˜ (θ, 0) = Q1 ⊕ Q2 (θ ) M where
Q1 :=
1 1
1 , 1
Q2 (θ ) :=
P2 P2 A2 (θ )
⋆ P2
.
˜ (0)}. + x˜ TG PG x˜ G ≤ X T (0)PX For any γ , δ ≤ 1, it is clear that X (0) ∈ E . Note also that independent of the values given to γ and δ , the block-diagonal structure of P guarantees that
sup |˜u| =
Since the eigenvalues of a block diagonal matrix are simply the ˜ (θ , 0) is positive semidefieigenvalues of the individual blocks, M nite. It has the same eigenvalues as Q1 , namely 0 and 1, in addition to the eigenvalues of Q2 , which are all strictly positive for sufficiently small γ and δ .
˜ (·, ·) as a 4 × 4 block matrix, the zero Continuing to treat M ˜ (θ , 0)), which is unique from Assumption 2, has eigenvalue λ1 (M the corresponding normalised eigenvector (in block form) ˜ (θ, 0)) = √1 [1; −1; 0; 0]. v1 ( M 2 ˜ (θ , 0)) with respect to Applying Lemma 18, the derivative of λ1 (M
κ at κ = 0 is given by
˜ ˜ (θ , 0)) ∂λ1 (M ˜ (θ , 0)) ∂ M (θ , 0) v1 (M ˜ (θ , 0)) = vT1 (M ∂κ ∂κ = σ G(−1)H (−1)aζ (k). Since α ≤ ζ (·) ≤ β , it is evident that
˜ (θ , 0) ∂M σ G(−1)H (−1)aα ≤ ≤ σ G(−1)H (−1)aβ. ∂κ Furthermore, since σ G(−1)H (−1) > 0, the derivative is positive, independent of the values taken on by the other components of θ . ˜ (θ , κ)) will move into the right half-plane This means that λ1 (M as κ is increased by a small amount from zero. The continuity of ˜ (θ, κ) in κ also ensures that there will be some small κ for which M the remaining positive eigenvalues will still be situated in the right ˜ (θ , κ) and hence M (θ , κ) positive definite half-plane, rendering M for sufficiently small κ . By increasing κ beyond some value κmax however, it will become indefinite; this is certain to occur when κ is large enough to cause A(θ , κ) to have eigenvalues outside the unit disc. It can therefore be concluded that for a sufficiently small γ and δ , the Lyapunov condition (47) will be satisfied for κ < κmax . The original system (1) can now be shown to converge, for a sufficiently small κ that depends on the initial conditions. 3.5.2. Proof of ES scheme convergence (Theorem 14) Consider any initial state of the LPV system X (0) = [˜u(0); w(0); x˜ H (0); x˜ G (0)],
˜ (0) X T (0)PX
X ∈E
sup x˜ G =
˜ (0) X T (0)PX
X ∈E
λ(PG )
.
Hence, for any state within or on the boundary of E , Lemma 19 guarantees that ∥Ψ (θ )∥ < N, for some N < ∞. This bound will also take into account the initial conditions u˜ (−1) and x˜ G (−1), since |h′ (ξ2 (0))| depends on them (see the proof of Lemma 19 in Appendix F). Then, from Lemma 17 and Theorem 20, it is possible to select κ , γ ≤ 1 and δ ≤ 1 sufficiently small without affecting this bound to ensure that R(A(θ , κ), P ) ≻ 0, for all θ such that ∥Ψ (θ )∥ < N; this is true for all X ∈ E . Lemma 16 then implies that the quadratic form V (X ) = X T PX is a Lyapunov Function for all X ∈ E . Hence, E is now a level set of V (·), making it an invariant set for the LPV system (41). Given that X (0) ∈ E , exponential stability of the LPV system is guaranteed, meaning that as k → ∞, u˜ → 0 exponentially. Since u¯ is related to u˜ by a translation (24), it can be concluded that u¯ exponentially converges to the equilibrium solution p0 , where
|p0 − p∗ | < |˜a|. This in turn ensures that y → qˆ 0 (p0 ) (16d). Remark 21. This ES Scheme results in a plant output converging to a constant value. This is a markedly different behaviour to most other dither-based schemes, which exhibit periodic orbits. The FIR filter with its zero at z = −1 is the main reason for this behaviour, as is evident from the nature of the equilibrium solution described in Theorem 8. Remark 22. A practitioner implementing (1) simply has to choose the dither amplitude a based upon the desired proximity to the optimum, σ to compensate for any phase-shift introduced by the dynamics and κ to be sufficiently small for stability. This is a notable improvement on existing ES methods, since Theorem 14 gives convergence bounds on u¯ , namely |˜a| either side of p∗ , rather than simply quantifying it as O(|˜a|). The requirement on κ being small is similar to other schemes in the literature; analysing stability for larger κ would require more knowledge of the system than has been assumed in this paper. Note that even if the bound on ∥Ψ (θ )∥ is arbitrarily large, as determined by the initial conditions of (1), a sufficiently small κ always exists to guarantee stability, since all of the terms containing |h′ (ξ2 )| are multiplied by κ .
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
3005
4. Example Two examples will be presented to demonstrate the efficacy of the ES scheme developed in this paper. The first example shows the behaviour of the controller on a simulated Wiener–Hammerstein plant, for various initial conditions. An experimental implementation of the scheme is then detailed, illustrating its ability to increase the mean maximum brake torque of a gas fuelled spark-ignition engine operating at constant speed, by adjusting the spark advance angle. 4.1. Simulated system Consider the Wiener–Hammerstein plant described by f (p) =
AG CG
1
1
3
(p − 5)2 − sin(p − 5) + (p − 5) 4 2 2 1 1 −2 1 2 BG 1 1 = 12 , 4 2 DG
AH
BH
CH
DH
5 4
=
1 2 1 3
0 1 1 3
0
(50a) Fig. 4. ES scheme input trajectories for different initial conditions.
(50b)
.
The static nonlinear map in this case is asymmetric about the location of the extremum p∗ = 5 (f (5) = 0). Selecting a = 1, κ ∈ {0.2, 0.6} and σ = −1, the ES scheme (1) is applied to this system with several initial values of u¯ . For each of these initial values, xG and xH are initialised at their corresponding steady-state values. Figs. 4–5 illustrate the resulting trajectories of u¯ and y respectively. As can be seen, the equilibrium of p0 ≈ 4.9681 is close to the true optimising value. In this case, the effective dither amplitude is given by a˜ ≈ −0.8824, where a phase shift has been introduced by G(·). For κ = 0.6, the dynamics are seen to initially generate transient behaviour which subsides as the system approaches steady state. Convergence is achieved after approximately 60 steps to qˆ 0 (p0 ) ≈ 0.0013, which is very close to the actual minimum of zero. It is evident that the transient dynamics are not excited by using a smaller κ , leading to smoother trajectories, at the expense of slower convergence. Further insight into convergence characteristics can be seen in Fig. 6, which illustrates the input–output behaviour of (50) in the closed loop, together with the steady-state map (50a). It can be seen that for initial conditions far from the extremum, there is an initial ‘‘fast’’ convergence behaviour exhibited with κ = 0.6, before the trajectories begin to approximately follow the path of steepest descent. The κ = 0.2 case exhibits more conventional steepest descent trajectories.
Fig. 5. ES scheme output trajectories for different initial conditions.
4.2. Experimental results Spark-ignition engines provide a practical platform upon which the ES scheme described in this paper can be experimentally tested. For a given fuel composition, an engine operating at a constant speed condition has a mean steady-state spark advance to maximum brake torque mapping which is approximately quadratic (Mohammadi et al., 2012). In engines using alternative fuels such as natural gas, LPG and ethanol-gasoline blends, the fuel in the tank is subject to compositional variation due to cost, geographical region and availability of different fuels, amongst other factors. It is therefore desirable to vary the spark advance in response to changing composition to ensure that the engine is operating at peak efficiency. The ES scheme (1) is implemented to vary the spark angle on a modified 6-cylinder Ford Falcon BF MY2006 engine coupled to
Fig. 6. Plant output vs. input for different initial conditions.
an eddy-current dynamometer, at an initial operating condition of 800 RPM and 32 N.m torque output. An ATI VISION system is used to modify the spark advance commanded by the ECU, measured in terms of Crank Angle Degrees before Top Dead Centre (CAD BTDC). By considering only mean-value engine behaviour (Guzzella & Onder, 2004), nonlinear dynamics arising from combustion and cyclic variations can be neglected, for an appropriately selected sampling rate. The overall engine and dynamometer system
3006
R.C. Shekhar et al. / Automatica 50 (2014) 2998–3008
Acknowledgements The authors would like to thank the Advanced Centre for Automotive Research and Testing (ACART) for the use of their dynamometer facility, as well as Alireza Mohammadi for his assistance with running the experiments. Appendix A. Proof of Lemma 10 For the first property, twice differentiating (27) and applying the mean value theorem gives 1
(f ′ (˜p + p0 + a˜ ) − f ′ (˜p + p0 − a˜ )) 2a˜ = f ′′ (c1 ),
g ′′ (˜u) =
Fig. 7. Convergence of (top) engine input; and (bottom) measured output towards their respective optima following implementation of the proposed scheme to online spark calibration of a CNG engine.
can then be adequately represented by a Hammerstein model, with output dynamics arising from the torque sensor on the dynamometer. Fig. 7 shows how the output torque of the engine can be optimised by modifying the spark advance angle, for running on pure methane. A sampling rate of 0.2 Hz has been used for updating the control action. From an initial spark advance angle of 26 °, which is not the optimal angle for methane fuel, the results show an improvement of approximately 3 % in mean maximum brake torque from the initial condition by adjusting the spark advance angle to just above 34 °. The effect of the discrete dither signal is clearly seen in the plot of spark angle. Note the noisy torque signal arises from cyclic variations in engine operation. Its effect on the ES scheme is minimised by using the average torque measurement over each sample. Remark 23. Using the alternating dither signal for spark timing modification has no detrimental effect on the system, since the spark plugs are driven electronically and contain no moving parts. If mechanical wear is of concern for other applications, the dither amplitude can be reduced to a suitably small value. In addition, the continuous-time realisation of the dither signal need not necessarily be a square wave obtained by applying a zeroorder hold. Other choices are possible, but will lead to different discretisations of the dynamics.
for some c1 ∈ (˜p + p0 − |˜a|, p˜ + p0 + |˜a|). Hence, g ′′ (·) must satisfy the same bounds as f ′′ (·). Furthermore, since a˜ g ′ (˜p) = qˆ 1 (˜p + p0 ), the second property follows directly from the first property and the equilibrium solution defined in Theorem 8. Appendix B. Proof of Lemma 11 Twice differentiating (28) and applying the intermediate value theorem shows that 1 h′′ (˜p) = (f ′′ (˜p + p0 + a˜ ) + f ′′ (˜p + p0 − a˜ )) 2 = f ′′ (c2 ), for some c2 ∈ (˜p + p0 − |˜a|, p˜ + p0 + |˜a|). The result follows directly as in Lemma 10. Appendix C. Proof of Lemma 13 Differentiating (28) and using a Taylor series expansion for each term with the Lagrangian form for the remainder, it can be shown that a˜ h′ (˜p) = f ′ (˜p + p0 ) + (f ′′ (c1 ) − f ′′ (c2 )), 2 for some c1 and c2 . Applying the triangle inequality and the bounds on f ′′ (·) from Assumption 1 gives the desired result. Appendix D. Proof of Lemma 17 Expanding (44), it is required that
5. Conclusions and future work This paper has presented a new discrete-time extremumseeking algorithm for discrete-time Wiener–Hammerstein plants with stable linear dynamics. It has shown how limited plant knowledge can be used to guarantee semi-global stability for appropriate choices of tuning parameters, assuming the system possesses a strongly-convex steady-state mapping. A practitioner can utilise this theoretical result for tuning the ES controller by choosing the dither amplitude based upon the desired proximity to the true optimum and then selecting the adaptation gain to be sufficiently small. The algorithm has been demonstrated in simulation as well as experimentally on a spark-ignition engine. Future work will analyse Wiener–Hammerstein-like systems with a Multi-Input Single-Output (MISO) nonlinear maps. It will also investigate whether upper bounds on the maximum adaptation gain for stability can be provided, as in Shekhar et al. (2013). An additional avenue of investigation will be the incorporation of input and output constraints into the formulation.
γ R(P11 , A11 ) −γ AT12 P11 A11
⋆ ≻ 0. T R(P22 , A22 ) − γ A12 P11 A12
Using the Schur complement, this is equivalent to the condition
R(P22 , A22 ) − γ AT12 Q(P11 , A11 )A12 ≻ 0,
(D.1)
where
Q(P , A) := P + PAR(P , A)−1 AT P . Applying the matrix inversion lemma and the Schur complement again, it follows that −1 Q(P11 , A11 ) = R(P11 , AT11 )−1 ≻ 0,
so (D.1) can be written in the form −1 R(P22 , A22 ) − γ AT12 R(P11 , AT11 )−1 A12 ≻ 0,
which can be satisfied by choosing γ such that 0
