Automatica 48 (2012) 583–594
Contents lists available at SciVerse ScienceDirect
Automatica journal homepage: www.elsevier.com/locate/automatica
A control-theoretic study on iterative solutions to nonlinear equations for applications in embedded systems✩ Ying Yang a,1 , Steven X. Ding b a
State Key Lab for Turbulence and Complex Systems, Department of Mechanics and Aerospace Engineering, College of Engineering, Peking University, Beijing 100871, PR China
b
Institute for Automatic Control and Complex Systems, University of Duisburg-Essen, 47057 Duisburg, Germany
article
info
Article history: Received 22 November 2010 Received in revised form 5 August 2011 Accepted 18 August 2011 Available online 3 March 2012 Keywords: Numerical analysis Stability theory Fixed point iteration Newton’s methods Observer design LMI technique
abstract In this paper, the fixed point iteration and Newton’s methods for iteratively solving nonlinear equations are studied in the control theoretical framework. This work is motivated by the ever increasing demands for integrating iterative solutions of nonlinear functions into embedded control systems. The use of the well-established control theoretical methods for our application purpose is inspired by the recent controltheoretical study on numerical analysis. Our study consists of two parts. In the first part, the existing fixed point iteration and Newton’s methods are analysed using the stability theory for the sector-bounded Lure’s systems. The second part is devoted to the modified iteration methods and the integration of sensor signals into the iterative computations. The major results achieved in our study are, besides some academic examples, applied to the iterative computation of the air path model embedded in the engine control systems. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction Iterative computation is one of the standard techniques for solving nonlinear equations (Quarteroni, Sacco, & Saleri, 2000; Stoer & Bulirsch, 2002). It is a powerful mathematical tool widely used in engineering applications (Eich-Soellner & Führer, 1998; Hoffmann, Marx, & Vogt, 2005). Thanks to the rapid development of microprocessor technology, more and more iterative solutions of nonlinear equations are implemented on electrical control units (ECUs) in real-time embedded systems. For instance, for the real-time control and on-board-diagnosis (OBD) of an internal combustion engine numerous iterative computation blocks are integrated into the ECU (Kiencke & Nielsen, 2005). Nonlinear equations like x = ϕ(x)
and f (x) = 0
(1)
✩ This work is supported by the National Natural Science Foundation of China under Grant 60874011, 61174052, 90916003. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Yoshito Ohta under the direction of Editor Roberto Tempo. E-mail addresses:
[email protected] (Y. Yang),
[email protected] (S.X. Ding). 1 Tel.: +86 10 62751815; fax: +86 10 62764044. She has completed this study
during her sabbatical year in Duisburg. 0005-1098/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2012.01.007
are most typical forms met in engineering applications (EichSoellner & Führer, 1998; Hoffmann et al., 2005). The so-called fixed point iteration described by x(k + 1) = ϕ(x(k)) and Newton’s methods with the general iterative form x(k + 1) = x(k) − Ψ (x(k))f (x(k)) are standard algorithms for solving (1) iteratively, where k stands for the iterative number and Ψ (x(k)) is some matrix (Quarteroni et al., 2000; Stoer & Bulirsch, 2002). Under certain conditions, the iteration will converge to the solution of the equations x∗ , i.e. lim x(k) = x∗ ,
k→∞
x∗ = ϕ(x∗ ) or f (x∗ ) = 0.
For the real-time applications in embedded control systems, the nonlinear equations in (1) build typically sub-models embedded in a complex functional block and thus are often a function of system inputs and parameters (see also the example in Section 5). In this context, the nonlinear equations in (1) can be extended to x = ϕ(x, p) and
f (x, p) = 0
(2)
with p being a vector representing these (time-varying) system inputs and parameters. The iterative computation will then be triggered by each update of p. For applications in embedded control systems, such computation often demands for high realtime ability. Although ECUs of the new generation are becoming
584
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
more powerful, the low cost requirement on the one hand and ever increasing demands for the high system performance on the other hand call for more attention to the real-time implementation of iterative algorithms. Our study presented in this paper is mainly motivated by the real-time implementation of control and OBD algorithms on the ECU for the internal combustion engine control (Weinhold, 2007). In this industrial application, we have been confronted with the following two problems: (a) the convergence rate of the existing iterative solutions for some nonlinear equations does not satisfy the real-time requirement (b) the quantisation errors due to the use of look-up tables instead of analytical functions may considerably affect the computation performance. As a result, poor control performance and, in the worst case, instability can be observed. As investigating solutions for these problems, we notice the recent efforts on applying the modern control theoretical methods to numerical analysis (Bhaya & Kaszkurewicz, 2006, 2007; Kashima & Yamamoto, 2007; Schaerer & Kaszkurewicz, 2001; Söderlind, 2002). In particular, encouraged by the result on the Newton’s method reported in Kashima and Yamamoto (2007), we focus our study on (a) the convergence conditions of the fixed point iteration and Newton’s methods for nonlinear functions satisfying the so-called sector conditions (Khalil, 2002) (b) the robustness issue with respect to computation errors e.g. caused by the quantisation. In this study, vector-valued nonlinear equations are considered, and the well-established nonlinear control theory (Khalil, 2002) and the LMI (linear matrix inequality) technique (Boyd, Ghaoui, & Feron, 1994) are applied as the main analysis and design tool. This effort allows us to express the convergence conditions in terms of some LMIs, which can then be checked using some standard software. Driven by the demands for high system performance and reliability, a trend can be observed in the area of embedded control systems that the number of the integrated sensors is continuously increasing. The intuitive idea of our further study on solving nonlinear equations is to integrate those sensor signals, which are available in the embedded system, into the iterative computation. This work has been strongly motivated by the moderate and, in some cases, poor performance delivered by the standard iterative methods as they were applied in the engine control and OBD algorithms running on the ECU. It is reasonable to expect that the additional information provided by the sensors will improve the convergence rate and enhance the robustness. For our purpose, an observer-like iterative method is proposed, which modifies the existing methods and allows us to integrate the available sensor signals so that both the convergence rate and robustness are improved. From the control-theoretical viewpoint, this work is a combination of our study on the existing iterative solutions and the well-established observer design technique. By testing our new approach to the practical real-time implementation problems, we notice that the sensor noises may have a strong influence on the computation performance. Considering that e.g. in automotive systems the sensors embedded in the control loops and in OBD are often low-cost products, the last topic in our study is dedicated to the robustness analysis with respect to the measurement noises. The paper is organised as follows. In Section 2, the needed preliminaries in numerical analysis and stability theory are briefly presented, based on which the main problems to be addressed in this paper are formulated. Section 3 is devoted to the controltheoretical study on the existing fixed point iterative solution and Newton’s methods with a focus on the convergence and robustness issues. In Section 4, we shall first propose a unified approach for modifying the existing iterative algorithms and integrating sensor signals aiming at improving the convergence rate. It is followed by an analysis of the influence of the measurement noises on
the computation performance. The last part in this section deals with the optimisation of the iterative schemes with respect to the convergence rate and robustness. To illustrate the major results, some (academic) examples are included in each section, and an application example from the real-time implementation of the air path model on the ECU is presented in Section 5. Notation. The notation adopted throughout this paper is fairly standard. Rn denotes the n-dimensional Euclidean space and Rn×m the set of all n × m real matrices. The superscript ‘‘T ’’ stands for the transpose of a matrix. ‘‘I’’ and ‘‘0’’ denote the identity and zero matrix with appropriate dimension, respectively. ‘‘P > 0 (≥ 0)’’ means matrix P is positive definite (semi-definite). σ¯ (·), σ (·) denote the maximum and minimum singular value of a matrix respectively. ∥ · ∥ stands for the Euclidean norm. For vector x ∈ Rn , √ ∥x∥ = xT x. We use E(·) to denote the expectation of a statistical or stochastic variable. 2. Preliminaries, basic ideas and problem formulation 2.1. Fixed point iteration and Newton’s methods Let the function ϕ : Rn → Rn have a fixed point x∗ : ϕ(x∗ ) = x∗ ∈ Rn . The fixed point iteration algorithm is given e.g. by Stoer and Bulirsch (2002) x(k + 1) = ϕ(x(k)) ∈ Rn .
(3)
The following definition and theorem are standard results in numerical analysis (Quarteroni et al., 2000; Stoer & Bulirsch, 2002). Definition 1. ϕ : D ⊆ Rn → Rn is called contractive on a set Do ⊂ D if there exists a constant α < 1 such that for all x, y ∈ Do
∥ϕ(x) − ϕ(y)∥ ≤ α∥x − y∥.
(4)
Theorem 1 (Eich-Soellner and Führer (1998), Banach’s Fixed Point Theorem). Let ϕ : D ⊆ Rn → Rn be contractive in a closed set Do ⊂ D and suppose ϕ (Do ) ⊂ Do . Then ϕ has a unique fixed point x∗ ∈ Do . Moreover, for any x(0) = x0 ∈ Do , the iteration (3) converges to x∗ . The distance to the solution is bounded by
∥x∗ − x(k)∥ ≤
αk ∥x(1) − x(0)∥. 1−α
(5)
It is worth mentioning that under the same conditions given in the above theorem it is easy to prove that
∥x∗ − x(k)∥ ≤ α k ∥x∗ − x(0)∥. n
(6) n
∗
Let the function f : R → R have a unique solution x ∈ Rn : f (x∗ ) = 0. The standard Newton’s method is an iterative algorithm described by x(k + 1) = x(k) − (Df (x(k)))−1 f (x(k))
(7)
with Df (x(k)) ∈ R as the Jacobian matrix of f (x(k)). The following theorem describes the major properties of Newton’s method (Stoer & Bulirsch, 2002). n×n
Theorem 2. Given f : D ⊆ Rn → Rn and the convex set Do ⊆ D, let f be differentiable for all x ∈ Do and continuous for all x ∈ D. For xo ∈ Do let positive constants r , α, β, γ be given with the following properties: Sr (xo ) = {x|∥x − xo ∥ < r } ⊆ Do , h = αβγ /2 < 1,
r =
α , 1−h
and let f (x) have the following properties:
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
585
(a) ∥Df (x) − Df (y)∥ ≤ γ ∥x − y∥ for all x, y ∈ Do , (b) (Df (x))−1 exists and satisfies ∥(Df (x))−1 ∥ ≤ β for all x ∈ Do , (c) ∥(Df (x(0)))−1 f (x(0))∥ ≤ α. Then x(k), k = 0, 1, . . . , is well defined and satisfies x(k) ∈ Sr (xo ), limk→∞ x(k) = x∗ , f (x∗ ) = 0, and
∀k ≥ 0,
∥x(k) − x∗ ∥ ≤ α
k h2 −1 k
1 − h2
.
(8)
It is evident that computing (Df (x(k)))−1 in each iteration needs considerable computational effort and has to be avoided for a realtime application. The so-called simplified Newton method (Stoer & Bulirsch, 2002) is a mostly used solution in this case, which is given by x(k + 1) = x(k) − (Df (x(0)))−1 f (x(k)),
(9)
where (Df (x(0)))−1 is a constant matrix. 2.2. Sector-bounded nonlinearities and Lure’s systems In their pioneering study, Kashima and Yamamoto (2007) have investigated the convergence condition of Newton’s method for scalar nonlinear equations satisfying the so-called sector conditions. It is demonstrated that replacing the Lipschitz (contractive) condition given in Theorem 1 by the sector boundedness may relax the convergence conditions. Motivated by this result, we shall study the convergence conditions of the fixed point iteration and Newton’s methods for vector-valued nonlinear equations with sector-bounded nonlinearities. The following definition of the sector condition is given by Khalil (2002), which is widely used in the study on system passivity, nonlinear control (Khalil, 2002), observer (Arcak & Kokotovic, 2001) and filter design (Wang, Liu, & Liu, 2008). Definition 2. Function ψ : Rn → Rn is said to belong to the sector [K1 , K2 ] on D with K = K T = K2 − K1 > 0, K1 , K2 ∈ Rn×n , and denoted by ψ ∈ [K1 , K2 ] if
(ψ(x) − K1 x)T (ψ(x) − K2 x) ≤ 0 for all x ∈ D.
(10)
The sector boundedness on nonlinearities is quite general that also includes the widely used Lipschitz conditions as a special case (Khalil, 2002). One of the applications of the sector condition is the stability study on the so-called Lure’s system, which is in general described by y(z ) = G(z )u(z ) ∈ Rn ,
u(k) = −ψ(y(k)) ∈ Rn
(11)
Fig. 1. A loop transformation.
nonlinear equations satisfying the sector condition. Inspired by the work in Kashima and Yamamoto (2007), where scalar nonlinear equations are dealt with, we shall address this topic in the stability framework using the LMI technique. The major objective of this study is to derive some new convergence conditions, which may relax the existing ones, and to build the basis for investigating possible modifications of the existing methods, the integration of the sensor signal into the iterative computation and the robustness issues. As previously mentioned, in real-time embedded systems with low-cost ECUs, there may exist computation errors at each iterate. These can be caused by the quantisation error εx (k + 1) = x(k + 1) − xo (k + 1), where xo (k + 1) represents the ideal computation being free of the round-off error, or by the approximation error of the nonlinear function εϕ (k) = ϕ(x(k)) − ϕo (x(k)) with ϕo (x(k)) denoting the original function and ϕ(x(k)) its approximation by e.g. look-up tables. In practice, such computation errors may remarkably affect the computation performance. Similar to the study in Kashima and Yamamoto (2007), we model these errors in terms of an additional term ε(k) ∈ Rn in the iteration as follows
x(k + 1) = ϕ(x(k)) + ε(k), fixed point iteration, x(k + 1) = x(k) − ϕ(x(k)) + ε(k), Newton’s method with ϕ(x(k)) = (Df (x(k)))−1 f (x(k)).
(14)
We assume that ε(k) is unknown but bounded by sup ∥ε(k)∥ ≤ δε (> 0).
(15)
k
It is our task to find the maximum influence of ε(k) on limk ∥x(k) − x∗ ∥ and to minimise it. Suppose that the sensors, being available in the embedded control system and used for our purpose of improving the iterative computation, is modelled by
where G(z ) is the n × n transfer function matrix of an LTI (linear time invariant) system and ψ(y(k)) a nonlinear mapping satisfying the sector condition ψ ∈ [K1 , K2 ]. It is well-known that by means of a loop transformation as sketched in Fig. 1 (Khalil, 2002), we have
where y represents the vector of the measurements, η the measurement noise. It is assumed that
ψ˜ T (y(k))y(k) ≥ 0 with
(12)
Eη = 0,
˜ (z )˜u(z ), y(z ) = G
(13)
matrix C ∈ Rm×n is known. In real applications, a sensor model may be nonlinear. In that case, a linearised model can be used. Also, advanced nonlinear observer methods (Arcak & Kokotovic, 2001) are helpful. Recall that in embedded systems iterations are often triggered by each update of the system inputs and parameters. According to (2), the solution x∗ is in fact a function of p. As a result, (16) can be rewritten into
˜ y(k)), u˜ (k) = −ψ( ˜G(z ) = K (I + G(z )K1 )−1 G(z ) + I , ˜ y(k)) = ψ(ς (k)) − K1 ς (k), ψ(
ς (k) = K −1 (y(k) − u˜ (k)).
System (13) satisfying (12) is called passive and the nonlinearity
˜ y(k)) is said belong to [0, ∞] and denoted by ψ˜ ∈ [0, ∞] ψ( (Khalil, 2002).
y = Cx∗ + η,
y ∈ R m , m < n,
E(ηT η) = E∥η∥2 ≤ δη2 ,
(16)
δη > 0,
(17)
2.3. Basic ideas and problem formulation
y = Cx∗ (p) + η.
Our first task is to study the convergence conditions of the fixed point iteration and Newton’s methods for vector-valued
For instance, the sampling time of today’s ECU for the modern internal combustion engine control system is 10 ms, which means
(18)
586
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
e.g. the angle of the throttle may change every 10 ms. It is required, in turn, that all those iterative computations of the nonlinear functions parameterised by the angle of the throttle (see also the example in Section 5) should be completed within a couple of milliseconds and before a new sampling. Thus, although the sensor signals are time sequences, y is constant during the whole iteration running between two sampling instants. Also, x∗ (p) is a function of the system parameters, but remains constant between two sampling instants. In the context of iterative computations of nonlinear functions, which depend on system inputs and parameters and are implemented in an embedded system, the integration of sensor signals modelled by (18) will enhance the real-time ability of the iteration schemes. The basic idea of integrating the sensor signals into the iterative computation is to add an additional term r (k) = L(y − Cx(k)) ∈ Rn in the iterative computations. Such an iterative computation system is called observer-like iteration due to its similar structure with an observer. Inspired by this structure, we will propose a unified form for the iterative computation, in which the sensor signals and further possible modifications can be included in a generalised form. Our task consists in (a) analysing the convergence rate and robustness of the unified iteration solution and (b) developing an approach to optimise the iterative computation by finding matrix L ∈ Rn×m and the possible modifications. 3. Study on the existing iterative solutions In this section, we will study the fixed point iteration and Newton’s methods, and discuss about the issues related to the needed computation efforts and robustness. 3.1. On the fixed point iteration
e(z ) = (zI − I )−1 u(z ),
u(k) = −ψ(e(k)).
By the above described loop transformation, we have
˜G : e(k + 1) = Ae(k) + u˜ (k), A := (I − K1 ) ˜ y(k)), y(k) = Ke(k) + u˜ (k), u˜ (k) = −ψ( ˜ y(k)) = ψ(ς (k)) − K1 ς (k), ς (k) = K −1 (y(k) − u˜ (k)) ψ(
(24)
˜ y) ∈ [0, ∞) on Do . Next, define V (k) = eT (k)Pe(k), P > with ψ( 0, 1V (k) = V (k + 1) − V (k). It holds 1V = eT (AT PA − P )e + ψ˜ T (y)P ϕ( ˜ y) − 2ψ˜ T (y)PAe.
(25)
Note that in the above equation and also in the remaining of this proof, we drop variable k for simplifying the notation, if no confusion is caused. Since (22) is, by Schur complement (Boyd et al., 1994), equivalent to AT PA − P < −ν I + (PA − K )T (2I − P )−1 (PA − K ),
(26)
it turns out
1V < −eT (PA − K )T (2I − P )−1 (PA − K )e ˜ y). −ν eT e + ψ˜ T (y)P ϕ( ˜ y) − 2eT AT P ψ(
˜ T (y)Ke and 2ψ˜ T (y)ψ( ˜ y) to and from the Adding and subtracting 2ψ right side of the above inequality lead to 1V < −ν eT e − eT (PA − K )T (2I − P )−1 (PA − K )e ˜ y) − ψ˜ T (y)(2I − P )ψ( ˜ y) −2eT (PA − K )T ψ(
˜ y)) = −ν eT e − z T z − 2ψ˜ T (y)y, −2ψ˜ T (y)(Ke − ψ(
For our purpose, we rewrite (3) into x(k + 1) = x(k) − (x(k) − ϕ(x(k))).
(19)
Let e(k) = x(k) − x be the iteration error. It holds ∗
e(k + 1) = e(k) − ψ(e(k))
(20)
with a new function ψ(e(k)) ∈ Rn defined by
ψ(e(k)) = e(k) − (ϕ(x(k)) − ϕ(x∗ )) = e(k) − (ϕ( ˆ e(k)) − ϕ(x∗ )),
Proof. The proof consists of three steps. In the first step, a loop transformation, as shown in Fig. 1, is carried out. To this end, write (20) into the form given in (11)
ϕ( ˆ e(k)) = ϕ(x(k)). (21)
Note that ψ(0) = ϕ(x∗ ) − ϕ(x∗ ) = 0. We now present and prove the first result of this work, whose proof plays an essential role in the subsequent study.
˜ y). Since where z = (2I − P )−1/2 (PA − K )e + (2I − P )1/2 ψ( ˜ y) ∈ [0, +∞) means ψ˜ T (y)y ≥ 0 for all e ∈ Do (Khalil, 2002), ψ( we have 1V < −ν eT e < 0
(27)
for e(k) ̸= 0. Following the Lyapunov stability theory (Khalil, 2002), system (24) and equivalently system (20) are asymptotically stable at origin, i.e. limk→∞ e(k) = 0 ⇐⇒ limk→∞ x(k) = x∗ . In the third step, the inequality (23) will be proven. It follows from (27) that V (k + 1) < V (k) − ν eT (k)e(k) ≤ σ 2 V (k). Note that P > ν I ⇒ 0 < σ¯ (νP ) < 1. Hence,
ν 0<σ = 1− σ¯ (P )
1/2
< 1.
(28)
Theorem 3. Consider the iterative algorithm (3) and the associated Eq. (20). Suppose that ϕ : Rn → Rn has a unique fixed point x∗ ∈ D ⊆ Rn and ψ ∈ [K1 , K2 ] on Do , where Do = {x − x∗ , x ∈ D}. If there exist ν > 0 and P = P T > 0 such that the following linear matrix inequality holds
As a result, we have
then for an arbitrary (start) point with x(0) ∈ D the iteration converges to the fixed point x∗ . Moreover, it holds
Theorem 3 provides us with an alternative convergence condition different from the contractive condition (4). Moreover, it demonstrates that the convergence rate of the iterative algorithm under the sector condition (10) is determined by σ and analogous to (6). In the following corollary and example we are going to check if the sector condition (10) relaxes the contractive condition (4).
∥e(k)∥ < τ σ k ∥e(0)∥, k = 1, . . . , 1/2 where 0 < σ = 1 − σ¯ (νP ) < 1, τ = σσ¯ ((PP )) .
Corollary 1. Suppose that ϕ : Rn → Rn has a unique fixed point x∗ ∈ D ⊆ Rn and the contractive condition (4) holds for all x, y ⊆ D. Then the convergence conditions given in Theorem 3 are satisfied.
νI − P −K P − PK1
−K −2I P
P − K1T P < 0, P −P
(22)
(23)
V (k) < σ 2k V (0) H⇒ σ 1/2 (P )∥e(k)∥ < σ¯ 1/2 (P )σ k ∥e(0)∥
⇐⇒ ∥e(k)∥ < τ σ k ∥e(0)∥. This ends the proof of Theorem 3.
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
587
Proof. It follows from the contractive condition (4) that
∀x + x∗ ∈ D,
∥ϕ(x + x∗ ) − ϕ(x∗ )∥ ≤ α∥x∥.
On the other hand,
∥ϕ(x + x∗ ) − ϕ(x∗ )∥ = ∥ϕ(x + x∗ ) − ϕ(x∗ ) − x + x∥ = ∥x − (ϕ(x + x∗ ) − ϕ(x∗ )) − x∥ = ∥ψ(x) − x∥ H⇒ ∥ψ(x) − x∥2 − α 2 ∥x∥2 ≤ 0 H⇒ (ψ(x) − (1 − α)x)T (ψ(x) − (1 + α)x) ≤ 0,
(29)
i.e. ψ(x) ∈ [K1 , K2 ] on Do = {x − x , x ∈ D} with K1 = (1 − α)I , K2 = (1 + α)I and so K = 2α I . To check if the LMI condition ∗
(22) is satisfied, we study (26), which is, by Schur complement lemma, equivalent to (22) AT PA − P + ν I PA − K
AT P − K P − 2I
2 α P − P + νI ⇐⇒ α(P − 2I )
<0
α(P − 2I ) <0 P − 2I
⇐⇒ (2α 2 + ν)I < P < 2I .
Fig. 2. Iterative computation results in Example 1.
(30)
Since P > 0, ν > 0 satisfying (30) exist, we can conclude that both the sector condition ψ ∈ [K1 , K2 ] and the LMI (22) are satisfied. The corollary is thus proven. Corollary 1 illustrates that the result with the sector condition given in Theorem 3 includes the one satisfying the contractive condition (4). The following example adopted from Kashima and Yamamoto (2007) will show that the convergence conditions given in Theorem 3 are less strict than the contractive condition (4). Example 1. Consider x = ϕ(x),
ϕ(x) =
2x3 + x2 − 1 3x2 + 2x + 1
− 4e + 2
(ψ(e) − ε e)(ψ(e) − (2 − ε)e) ≤ 0,
∀e ∈ R
satisfies the sector condition
i.e. ψ ∈ [ε, 2 − ε] on R, where 0 < ε ≪ 1 is a constant enough close to zero. With K1 = ε, K = 2 − 2ε , it is straightforward to check that for β > ε and 0 < β − ε ≪ 1 being enough close to zero P =
1−ε
=1+
1−β 1−ε
<2
and a small enough ν > 0 solve the LMI (22). Thus, it can be concluded that starting from any point in R, the iteration error will converge to zero. Fig. 2 shows the results with three different initial starts: x(0) = 0, x(0) = 1, x(0) = −2. All three iterations converge to −1.
νI − P −K P − PK1
−K −2I P
P − K1T P < 0, P −P
(32)
Next, we demonstrate the relations between the convergence conditions given in Theorems 2 and 4.
Proof. In Stoer and Bulirsch (2002), Section 5.3, it is proven that conditions (a) and (b) in Theorem 2 lead to
∥D−1 f (y)(f (x) − f (y) − Df (y)(x − y))∥ ≤
Consider the Newton’s method given in (7). Let
βγ 2
∥x − y∥2
(33)
for all x, y ∈ Do , where β, γ > 0. Let 0 < α˜ < 1 such that D∗o ⊆ Do be a neighbourhood of x∗ D∗o =
x| ∥x − x∗ ∥ ≤
α˜ βγ ,ω = ω 2
.
Now, setting x = x∗ , y = e + x∗ ∈ D∗o in (33) yields
∥D−1 f (e + x∗ )f (e + x∗ ) − e∥ ≤
3.2. On Newton’s methods
βγ 2
∥e∥2 ≤ α∥ ˜ e∥.
(34)
Note that (34) can be rewritten as
e(k) = x(k) − x∗ ,
∥ψ(e) − e∥ ≤ α∥ ˜ e∥,
ψ(e(k)) = (Df (x(k))) f (x(k)) = (Df (e(k) + x∗ ))−1 f (e(k) + x∗ ) −1
H⇒ e(k + 1) = e(k) − ψ(e(k)).
Corollary 2. Consider the iterative algorithm (7) and the associated equation (31). Suppose that the conditions in Theorem 2 are satisfied. Then there exist a neighbourhood of x∗ , D∗o ⊆ Do , and associated with it D∗ = {x − x∗ , x ∈ D∗o }, ν > 0 and P = P T > 0 so that ψ ∈ [K1 , K2 ] on D∗ and the LMI (32) is satisfied.
e3 − 2e2 + 2e 3e2
2−ε−β
Theorem 4. Consider the iterative algorithm (7) and the associated equation (31). Suppose that f : Rn → Rn has a unique solution x∗ ∈ D ⊆ Rn and ψ ∈ [K1 , K2 ] on Do , where Do = {x − x∗ , x ∈ D}. If there exist ν > 0 and P = P T > 0 such that the following linear matrix inequality holds
then for an arbitrary (start) point x(0) with x(0) ∈ D the iteration converges to the solution x∗ .
,
for which we have solution x∗ = −1. It is clear that the contractive condition (4) is not satisfied e.g. by checking |ϕ(1) − ϕ(0)| = 1 + 13 > 1. On the other hand,
ψ(e) = e − (ϕ(x) − ϕ(−1)) =
(31) is identical with (20). Thus, we have the following theorem whose proof is identical with the one of Theorem 3 and therefore omitted.
(31)
e ∈ D∗ ,
which is of the same form like (29). Thus, using the results given in the proof of Corollary 1 we can conclude that ψ ∈ [K1 , K2 ] on D∗ and the LMI (32) is satisfied.
588
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
Remember that the real-time ability is the major concern in our study on the application of the iterative solutions to the embedded systems. Among the modified Newton’s methods (Quarteroni et al., 2000), the one given in (9) is,widely used in the embedded system applications. Based on Theorem 4, the following convergence conditions can be derived.
for P > 0, ν > 0, K = K T > 0 and Q , where B¯ = I for the fixed point and Newton’s methods and B¯ = (Df (x(0)))−1 for the modified Newton’s method (9). (b) setting
Corollary 3. Consider the iterative algorithm (9). Suppose that f :
and K2 = K + K1 . The LMI (37) and Eq. (38) follow directly from (22) and (35) and can be solved using e.g. the above-mentioned MATLAB LMI toolbox.
Rn → Rn has a unique solution x∗ ∈ D ⊆ Rn and f˜ ∈ [K1 , K2 ]
on Do , where f˜ (x) = f (x + x∗ ), Do = {x − x∗ , x ∈ D}. Let B = (Df (x(0)))−1 . If there exist ν > 0, P = P T > 0 such that the following linear matrix inequality holds
νI − P
−K
P − K1T BT P
−K P − PBK1
−2I
BT P −P
PB
< 0,
(35)
then for x(0) ∈ D the iteration converges to x∗ .
(36)
˜G : e(k + 1) = Ae(k) + Bu˜ (k), A = (I − BK1 ) ˜ y(k)), y(k) = Ke(k) + u˜ (k), u˜ (k) = −ψ( ˜ y(k)) = f˜ (ς (k)) − K1 ς (k), ς (k) = K −1 (y(k) − u˜ (k)). ψ( The remaining proof of LMI (35) can be done along with the lines in the proof of Theorem 3 and thus omitted. 3.3. Computation issues In our previous study, no modification has been made on the iterative algorithms and thus, in the real-time application, no additional on-line computation is needed. On the other hand, different (off-line) computational efforts for checking the proposed convergence conditions are expected. Below, we will analyse the needed (off-line) computational costs briefly and propose an algorithm for a systematic checking procedure. It follows from Theorems 3 and 4 and Corollary 3 that the computations introduced above for checking the convergence conditions consist of two steps: (a) checking the sector condition (10) (b) solving an LMI. When a scalar nonlinear function is addressed, the sector condition (10) can be checked analytically or graphically as shown in Example 1. The latter is preferred in engineering applications, since the computation range is, in general, physically well defined and considerably limited. In case of dealing with a n-dimensional vector of nonlinear functions, Khalil (2002) suggests to set K1 , K2 to be n-dimensional diagonal matrices if the nonlinear functions are decoupled. Thus, the problem is reduced to check n scalar sector conditions. In some cases, choosing K1 , K2 suitably may result in decoupled nonlinear functions, see e.g. Examples 2 and 3 in this paper. In our study, the LMIs given in Theorems 3 and 4 and Corollary 3 are solved using the MATLAB LMI toolbox (Gahinet, Nemirovski, Laub, and Chilali, 1995). It should be pointed out that the solvability of the LMIs given in (22), (32) and (35) depends on K1 and K2 (K ). For instance, in Example 1, although ψ ∈ [0, ∞] on R, the LMI (22) is solvable only for ψ ∈ [ε, 2 − ε] on R. In other words, it is often necessary to determine suitable K1 and K2 (K ). For this purpose, we propose the following algorithm: (a) solving
−K −2I P B¯
P − QT B¯ T P < 0 −P
(37)
(38)
3.4. Convergence error In this subsection, we analyse the convergence error caused by the computation error ε(k) in each iterate. Let e(k) = x(k) − x∗ . It follows from (14) that for both the fixed point iteration and Newton’s method
ψ(e(k)) =
A loop transformation shown in Fig. 1 results in
νI − P −K P −Q
P −1 Q , for the iterations (3), (7) Df (x(0))P −1 Q , for the iteration (9)
e(k + 1) = e(k) − ψ(e(k)) + ε(k)
Proof. Let e(k) = x(k) − x∗ . It holds e(k + 1) = e(k) − Bf (x(k)) = e(k) − Bf˜ (e(k)).
K1 =
(39)
ˆ e(k)) − ϕ(x∗ ) , e(k) − ϕ(
for the fixed point iteration
ˆ e(k)), for the Newton’s method with ϕ( ϕ( ˆ e(k)) = ϕ(x(k) + x∗ ).
We suppose that ψ ∈ [K1 , K2 ] on D, and the convergence conditions in Theorems 3 and 4 are satisfied respectively for the fixed point iteration and Newton’ method, which ensure x(k) converges to x∗ if ε(k) = 0. Next, we analyse how large the convergence error will become for the bounded computation error ε(k). Let V (k) = eT (k)Pe(k), P > 0. It holds
1V = V (k + 1) − V (k) = (e − ψ(e))T P (e − ψ(e)) − eT Pe + ε T P ε + 2(e − ψ(e))T P ε H⇒ 1V ≤ (e − ψ(e))T P (e − ψ(e)) − eT Pe + ε T P ε + 2σ¯ (P )∥e − ψ(e)∥∥ε∥. Note that ψ ∈ [K1 , K2 ] is equivalent to
(ψ(e) − Me)T (ψ(e) − Me) − eT N T Ne ≤ 0 M = (K1 + K2 )/2, N = (K2 − K1 ) /2 H⇒ ∥e − ψ(e)∥ ≤ ∥ψ(e) − Me∥ + ∥(M − I )e∥ ≤ (eT N T Ne)1/2 + ∥(M − I )e∥. Recalling that under convergence conditions given in Theorem 3 or Theorem 4 we have V (k) + (e − ψ(e))T P (e − ψ(e)) − eT Pe < σ 2 V (k) with σ = (1 − σ¯ (νP ) )1/2 < 1, it yields V (k + 1) < σ 2 V (k) + γ1 δε
V (k) + σ¯ (P )δε2 ,
γ1 = 2σ¯ (P )σ¯ (P −1/2 )(σ¯ (N ) + σ¯ (M − I )).
(40)
Let δ > 0 satisfy
ξ = σ 2 + γ1 /δ + σ¯ (P )/δ 2 ≤ 1.
(41)
As a result, if V (k) ≥ δ δε , we have V (k + 1) < ξ V (k). That means (a) if x(0) is selected so that V (0) ≥ δ 2 δε2 , 2 2
lim V (k) ≤ δ 2 δε2 H⇒ lim ∥e(k)∥ ≤ σ¯ (P −1/2 )δδε , k
k
(42)
(b) if x(0) is selected so that V (0) < δ 2 δε2 , V (k) ≤ δ 2 δε2 H⇒ ∥e(k)∥ ≤ σ¯ (P −1/2 )δδε . In summary, we have proven the following theorem.
(43)
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
Theorem 5. Consider (39) with bounded ε(k) whose boundedness is given in (15). Suppose that the convergence conditions in Theorems 3 and 4 are satisfied. Then, lim ∥e(k)∥ ≤ ϑδε , k
ϑ = σ¯ (P −1/2 )δ > 0
(44)
with δ satisfying (41). (44) gives an estimation of the upper bound of the convergence error. Recall that σ in (41) is a key parameter which decides the convergence rate of the iteration (see also Theorem 3). If σ in (41) is very close to 1, δ and so ϑ will become very large to ensure (41). In other words, if the error-free iteration (i.e. ε(k) = 0) converges (very) slowly to its solution x∗ , the convergence error caused by ε(k) may become (very) large.
which is, to simplify the notation, unified written as e(k + 1) = AL e(k) − Bψ(e(k)), A = I − LC , ψ(e(k)) = f (e(k) + x∗ ) L
We now study improving the real-time ability of the iterative solutions by modifying the algorithms and embedding the sensor signals into the iterative computation. 4.1. A unified form of the observer-like solutions
(47)
for the modified Newton’s iteration (45) ψ(e(k)) = x∗ − ϕ(e(k) + x∗ ) for the modified fixed point iteration (46).
AL = I − B − LC ,
We now introduce a theorem that provides us with the sufficient conditions for the convergence of (45) and (46). Theorem 6. Consider the iterative algorithm (45) or (46). Suppose that they have a unique solution x∗ ∈ D ⊆ Rn and ψ ∈ [K1 , K2 ] on Do , Do = {x − x∗ , x ∈ D}. If there exist ν > 0, P = P T > 0 and Q , R such that the following linear matrix inequality holds
νI − P −K P − Q K˜ − RC
4. A unified observer-like iterative solution
589
P − K˜ T Q T − C T RT < 0, QT −P
−K −2I Q
(48)
where K˜ = K1 for the modified Newton’s method (45) and K˜ = I + K1 for the modified fixed point iteration (46), then for B = P −1 Q , L = P −1 R and x(0) ∈ D the iteration converges to x∗ with
1/2 σ¯ (P ) ν τ= , σ = 1− . σ (P ) σ¯ (P )
For iteratively solving nonlinear equations f (x) = 0 and x = ϕ(x), we propose the following two algorithms: (a) modified Newton’s method:
∥e(k)∥ < τ σ ∥e(0)∥,
x(k + 1) = x(k) − Bf (x(k)) + L(y − Cx(k)),
The proof of this theorem can be done along with the lines in the proof of Theorem 3 and is thus omitted. In order to demonstrate the application of the above result, we give the following example which shows the convergence of the observer-like Newton’s method.
(45)
(b) modified fixed point iteration: x(k + 1) = x(k) − B(x(k) − ϕ(x(k))) + L(y − Cx(k)), n×n
(46)
k
(49)
n×m
where B ∈ R and L ∈ R are design parameter matrices to be determined. The introduction of the term −Bf (x(k)) in (45) is a natural extension of (9), while adding −B(x(k)−ϕ(x(k))) in (46) is inspired by the study on the fixed point iteration in Section 3.1. It is worth noting that these two additional terms act as a feedback of the estimation error of the nonlinear equations in (1) at each iterate:
−Bf (x(k)) = B(0 − f (x(k))), −B(x(k) − ϕ(x(k))) = B(0 − (x(k) − ϕ(x(k)))). It is known from the control theory that suitably selecting B may improve the convergence rate. Both (45) and (46) are of the observer structure (Chen, 1984), whose core is the feedback of signal r (k) = L(y − Cx(k)) under use of the sensor model y = Cx∗ + η. It is well known that by a suitable selection of matrix L, which is also called observer matrix, feeding back r (k) will significantly improve the convergence rate of the iteration. On the other hand, a most important pre-condition for achieving nice convergence performance is the observability which depends on the output matrix C and the system matrix. The introduction of matrix B, which is a part of the system matrix, also serves for this purpose. It is evident that limk→∞ x(k) in (45) or (46) respectively solves f (x) = 0 or x = ϕ(x) if x(k) converges to x∗ .
Example 2. Consider the iteration algorithm (45) with
x31 + 2x21 + 2x1 + 1
2x2 + 2x + 1 1 1 . x32 + x1 x22 + x2 + x1
f (x) =
x22 + 1 It holds ψ = f (e + x∗ ) ∈ [K1 , K2 ] on R2 with ∗
x =
−1 1
,
0 K1 = −1
0 , 0
4 K2 = −1
0 . 3
Using the MATLAB LMI toolbox, it can be found out that (48) is infeasible for C = 0. That means we cannot find B to ensure the convergence of the iteration (45). Now, assume that we have an embedded sensor signal y = Cx∗ with C = [2 1]. It can be verified that (48) is feasible with 0.0906 B= −0.1795
−0.1951 , 0.3818
0.3236 L= . 0.3483
Fig. 3 shows the computation results of x(k) starting from the initial value x(0) = [1 − 1]T . 4.3. Convergence error caused by measurement noise
4.2. Convergence study For our purpose, we first assume η = 0 and let e(k) = x(k)− x∗ . It follows respectively from (45) and (46) that
e(k + 1) = (I − LC )e(k) − Bf (e(k) + x∗ )
for the modified Newton’s iteration (45),
∗ ∗ e(k + 1) = (I − B − LC )e(k) − B(x − ϕ(e(k) + x ))
for the modified fixed point iteration (46),
Considering that the influence of ε(k) on the convergence error limk ∥e(k)∥ can be estimated along with the discussion in Section 3.4, we shall, in this subsection, focus on the analysis of the influence of the measurement noise η on the convergence performance for both iteration algorithm (45) and (46). To this end, consider the unified form of the iteration errors e(k + 1) = AL e(k) − Bψ(e(k)) + Lη.
(50)
590
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
we have finally V (k + 1) < ϖ 2 V (k) + ληT η
H⇒ V (k + 1) < ϖ 2k V (0) + lim (eT (k)e(k)) <
k→∞
λ(1 − ϖ 2k ) T η η, 1 − ϖ2
λ ηT η σ (P )(1 − ϖ 2 )
H⇒ lim E(eT (k)e(k)) < k→∞
λ δ2 . σ (P )(1 − ϖ 2 ) η
(56)
The main result from the above discussion can be summarised in the following theorem. Theorem 7. Consider (50) with measurement noise η whose boundedness is given in (17). Suppose that the convergence conditions in Theorem 6 are satisfied. Then, Fig. 3. Computation results of x(k) with observer-like Newton’s method.
lim E∥e(k)∥2 <
k→∞
λ δ2 σ (P )(1 − ϖ 2 ) η
(57)
e(k) and ψ(e(k)) are correlated with η and thus
with ϖ , λ satisfying (54)–(55).
E(η e(k)) ̸= E(η )E(e(k)) = 0,
We would like to remark that ϖ will (very) close to 1 if σ in (49) (see Theorem 6) is (very) close to 1. Moreover, in this case ζ will become (very) small, which leads to a (very) large λ. In other words, if the noise-free iteration (i.e. η = 0) converges (very) slowly to its solution x∗ , the convergence error caused by η will become (very) large. Also, due to the term λ = σ¯ (LT ( ζ1 P + I )L) a large σ¯ (LT PL) will cause large convergence error.
T
T
E(ηT ψ(e(k))) ̸= E(ηT )E(ψ(e(k))) = 0. As a result, the standard methods known in the filtering theory cannot be directly used. Let V (k) = eT (k)Pe(k), P > 0. It turns out
1V (k) = V (k + 1) − V (k) = (AL e − Bψ(e) + Lη)T P (AL e − Bψ(e) + Lη) − eT Pe = (AL e − Bψ(e))T P (AL e − Bψ(e)) + 2(AL e − Bψ(e))T PLη + ηT LT Lη − eT Pe.
4.4. Design of the observer-like solution
Since for any ζ > 0, 2(AL e − Bψ(e))T PLη
≤
1
ζ
ηT LT PLη + ζ (AL e − Bψ(e))T P (AL e − Bψ(e))
(51)
Wang, Xie, and de Souza (1992), it holds
1V ≤ (ζ + 1)(AL e − Bψ(e))T P (AL e − Bψ(e)) 1 − eT Pe + ηT LT P + I Lη. ζ
(52)
Recall that if the convergence conditions given in Theorem 6 are satisfied, then we have
(AL e − Bψ(e))T P (AL e − Bψ(e)) − eT Pe < −ν ∥e∥2 ⇒ (AL e − Bψ(e))T P (AL e − Bψ(e)) < (σ¯ (P ) − ν)∥e∥2 .
min
ν,σ¯ (P )
As a result, setting ζ (> 0) such that ζ (σ¯ (P ) − ν) − ν < 0 leads to
+ ηT LT
1
ζ
P +I
Lη.
(53)
ν − ζ (σ¯ (P ) − ν) ϖ = 1− σ¯ (P ) 1 T λ = σ¯ L P +I L , ζ
1/2
1−
⇐⇒ min (σ¯ (P ) − ν),
ν >0 σ¯ (P )
σ¯ (P ) − ν > 0.
ν,σ¯ (P )
min ς subject to σ¯ (P ) < ς + ν,
ν,σ¯ (P )
Note that for ν I − P < 0, ζ (σ¯ (P ) − ν) − ν > −σ¯ (P ) results in ν−ζ (σ¯ (P )−ν) 1>1− > 0. Let σ¯ (P )
ν 1− , σ¯ (P )
< 1,
(58)
Let σ¯ (P ) − ν < ς . Then, we are able to write optimisation problem (58) as
V (k + 1) < V (k) + (ζ (σ¯ (P ) − ν) − ν)∥e∥2
Remember that for the real-time applications in embedded systems the performance of an iterative computation will be evaluated in terms of (a) the convergence rate (b) the robustness against computation errors and (c) in case of an integration of sensor signals, the robustness against noises. In this context, we propose an approach for the design of observer-like iterative solutions (45) and (46). For our purpose, we first define a cost function. It follows from (49) in Theorem 6 that the convergence rate is determined by σ = (1 − σ¯ (νP ) )1/2 , and the smaller σ is, the higher the convergence rate becomes. Moreover, the discussions at the end of Sections 3.4 and 4.2 show that enhancing the convergence rate by reducing σ will also reduce the convergence error. These facts motivate us to optimise the (observer-like) iteration by minimising σ , which can then be equivalently formulated as
(54)
(55)
min ς subject to ν,P
−(ς + ν)I P
ς > 0 ⇐⇒ P
−(ς + ν)I
< 0,
(59)
ς > 0.
(60)
Recall that σ¯ (LT PL) should be bounded in order to limit the influence of the measurement noise. To this end, we introduce σ¯ (LT PL) ≤ χ , for some given χ > 0. Together with the convergence conditions given in Theorem 6, we finally formulate our optimisation problem as follows.
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
591
Theorem 8. Consider iterative computation (45) or (46) whose error system is unified described by (50). Suppose the conditions given in Theorem 6 are satisfied. If there exist ν > 0, P = P T > 0 and Q , R such that the following linear matrix inequality optimisation problem has a solution: min ς subject to
νI − P −K P − Q K˜ − RC
−(ς + ν)I
R
RT −P
T
T
T
Q
P − K˜ Q − C R < 0, QT −P
P
−K −2I
−(ς + ν)I
P
−χ I
(61) T
< 0,
(62)
ς > 0,
(63)
< 0,
(64)
for given χ > 0, where in (62) K˜ =
Fig. 4. Computation of the modified fixed point iteration with and without feedback of sensor data.
K1 , for the mod. Newton’s method (45) I + K1 , for the mod. fixed point iteration (46)
then for any x(0) ∈ D the iteration converges to x∗ with the minimum convergence rate ς and the optimal solution for L, B is given by B = P −1 Q , L = P −1 R. Example 3. Consider the algorithm (3) and its modified form (46) with
− ϕ(x) =
x31 + 4x21 + 3x1 + 2
3x21 + 2x1 + 1
,
x1 + 0.5 sin(x2 + 1)
Fig. 5. Schematic description of the air path in the SI engines.
ignition) engines, which is integrated into the embedded engine control systems (Guzzella & Onder, 2010; Weinhold, 2007).
for which we have the solution x∗ = [−1 − 1]T . We assume that a sensor signal y = Cx∗ is available with C = [10 5]. Moreover, we set χ = 1. Since
ϕˆ ∈ [K1 , K2 ] on R2 ,
K1 =
0 −1
0 2 , K2 = −0.5 −1
0 , 0.5
we get the largest convergence rate with σmin = 0.5876 by solving (61)–(64) using the MATLAB LMI toolbox. For a comparison, we consider the algorithm (3), i.e. L = 0, B = 0. It turns out σmin = 0.865, which indicates a lower convergence rate. It is clear from Fig. 4 that the algorithm (46) with B, L delivers a better convergence performance than the algorithm (3) regarding to the same initial value x(0) = [1 2]T . 4.5. On the computation issue In comparison with the standard fixed point and Newton’s iteration algorithms (3) and (9), the realisation of the observerlike forms (45) and (46) needs the following additional on-line computations: L(y − Cx(k)) in (45) and (I − B)x(k) + L(y − Cx(k)) as well as Bϕ in (46). These computations will not significantly increase the on-line computation cost. For the needed off-line computations, the check of the sector boundedness requires the same computational efforts as mentioned in Section 3.3. In Examples 2, 3, K 1, K 2 are chosen in such a way: the off-diagonal items are chosen to make the nonlinear function decoupled, the diagonal items are determined graphically as the lower and upper bounds of the sector which the nonlinearity belongs to. The LMI (48) and the optimisation problem in Theorem 8 are standard ones (Scherer, Gahinet, & Chilali, 1997) and can be solved using the MATLAB LMI toolbox (Gahinet et al., 1995). 5. An application example In this section, we describe the application of our major results to the real-time computation of the air path model of the SI (spark
5.1. A brief system description and problem formulation Fig. 5 sketches the air path of the SI engine with
˙ LF = 3600 m
100pU (pU − pv DK ) cfr TU
,
0.5 pS ˙ DK = 3600A(αDK )100pvDK √ m Ψ pv DK RTv DK 0.5 · 100pS 1 ˙ Zyl = 3600VZyl m n · η(n). RTS 60
,
Variable
Meaning
Unit
˙ LF m ˙ DK m ˙ Zyl m pU pv DK pS n
air flow through the air filter air flow entering the intake manifold air flow entering the cylinder ambient air pressure = 1013 pressure in manifold behind air filter pressure in the intake manifold engine speed angle of the throttle
kg/h kg/h kg/h hPa hPa hPa rpm
αDK
o
The coefficients and physical constants are listed in Appendix.
Ψ ( p pS ), A(αDK ), η(n) represent three nonlinear functions, for v DK which only look-up tables shown in Fig. 6 are available. αDK , n may change strongly depending on the driver. In the steady state, we have
√
A(αDK ) RTS
pS = √ pv DK · Ψ Tv DK n · η(n) pv DK = pU −
25cfr TU pU
pS
pv DK
A(αDK )pv DK
√
RTv DK
Ψ
,
(65) pS pv DK
2
.
(66)
592
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594 1.5 1 0.5 0 –0.5
0
5
10
15
20
25
30
35
40
35
40
45
50
Iteration numbr (k) 1100 1000 900 800 700 0
5
10
15
20
25
30
45
50
Iteration numbr (k)
Fig. 7. The (unstable) computation result using the standard fixed point iteration.
Fig. 6. Look-up tables for the three nonlinear functions.
For the control and OBD purposes, pS , pv DK are iteratively computed using (65)–(66). This is done at each sampling instant and depending on n, αDK . Our major task consists in a reliable and real-time applicable iterative computation of pS , pv DK in the whole operating range defined by αDK , n. In addition, the computation should be as robust as possible against uncertainties in the lookp up tables for η(n), A(αDK ), Ψ ( p S ). v DK
5.2. The iterative computations and the results In order to reduce the coupling between the computations of pS , pv DK , we rewrite (65) into
√
pS pv DK
A(αDK ) RTS
= √
Tv DK n · η(n)
Ψ
pS pv DK
,
(67)
pS
and iteratively compute p instead of pS . For our purpose, the v DK standard fixed point iteration (13) has been first applied, which leads to
µ(k + 1) = c1 Ψ (µ(k)),
µ(k) =
pS (k) pv DK (k)
,
(68)
pv DK (k + 1) = pU − c2 (pv DK (k)Ψ (µ(k)))2 ,
(69)
pS (k) = µ(k)pv DK (k),
(70)
√ A(αDK ) RTS c1 = √ , Tv DK n · η(n)
c2 =
25cfr TU A2 (αDK ) pU
RTv DK
.
We have noticed that in some operating ranges, the iterations do not converge. Fig. 7 shows a representative (unstable) iteration at the operating point n = 2000, αDK = 8. The initial condition is pv DK (0) = 900, µ(0) = 700/900. This result has motivated us to check the sector conditions given in Theorem 3 for the fixed point iteration. In Fig. 8, the values of the nonlinear functions
ψ1 (e1 ) = e1 − c1 (αDK , n)(Ψ (e1 + µ∗ ) − Ψ (µ∗ )) = µ − c1 (αDK , n)Ψ (µ), ψ2 (e2 )∗ 2 2 (e2 + pvDK ) Ψ (e1 + µ∗ ) = e2 − pU + c2 (αDK ) −(p∗vDK Ψ (µ∗ ))2 = pvDK − pU + c2 (αDK )((pvDK Ψ (µ))2 ) e1 = µ∗ − µ, e2 = pv DK − p∗v DK are sketched for the (physically) realistic values of µ, pv DK , representatively at the operating points αDK = 8, n = 1500, 2000
Fig. 8. ψ1 (e1 ), ψ2 (e2 ) in the operating ranges.
respectively. It is evident that ψ1 ∈ [0, 4], ψ2 ∈ [0, 2], and thus the known sufficient sector condition [0, 2] is not satisfied for ψ1 . Motivated by this observation, we have decided to apply the modified fixed point iteration given in (46), first without integrating the sensor data (i.e. L = 0):
µ(k + 1) µ(k) = pv DK (k + 1) pv DK (k) µ(k) − c1 (αDK , n)Ψ (µ(k)) −B 2 . pv DK (k) − pU + c2 (αDK )(pv DK (k)Ψ (µ(k)))
Applying Theorem 6 to our system results in B = diag(0.39, 0.30). As shown in Fig. 9, the modified fixed point iteration with the above B delivers a reliable and fast iterative computation of pv DK , µ(pS ) over the operating range. Having ensured the stable iterative computation of pv DK , µ(pS ), we have further focused on enhancing the computation robustness against uncertainties in η(n), A(αDK ) by means of the observerlike iterative solution studied in Section 4. To this end, the measurement pS delivered by an embedded pressure sensor is integrated into the iteration. Considering that µ = pS /pv DK is a nonlinear function, which can be approximated by a standard linearisation around pv DK (k), pS 1(k), 1(k) = µ − pS /pvDK (k), µ ≈ pS /pvDK (k) − 2 pv DK (k) we have used y = pS /pv DK (k) as an approximated measurement of µ in order to simplify the design and analysis. On this assumption, applying Theorem 6 results in B = diag(0.32, 0.30), L = [0.19 0]T . As a result, the iterative computations are
µ(k + 1) = µ(k) − 0.32(µ(k) − c1 (αDK , n)Ψ (µ(k))) + 0.19(ps /pvDK (k) − µ(k)),
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594
593
signals into the iterative computations. This work is strongly motivated by the ever increasing demands for integrating reliable and fast converging iterative solutions into the ECUs in the embedded control systems. To this end, the major focus of our study is on the improvement of the convergence rate, and on enhancing the computation robustness against uncertainties in the nonlinear equations as well as the measurement noises. We have derived sufficient conditions for the convergence of the so-called observer-like iterative algorithms as well as the estimations of the convergence rate and the boundedness of the convergence errors. To illustrate the theoretical results, three academic examples are included in the paper. Our final result is a successful application to the iterative computation of the air path model embedded in the engine control systems. Fig. 9. Iterative computations of µ, pv DK with different initial conditions. The operating point in this example: αDK = 8, n = 2000. The real solution: p∗v DK = 887.5, p∗S = 718.5, µ∗ = 718.5/887.5 = 0.81.
Acknowledgements The authors thank Mr. de Moll for his help in the application study described in Section 5. The authors are also grateful to the reviewers for the constructive and valuable comments and suggestions. Appendix. Coefficients in the air path model
Fig. 10. Robustness comparison. The operating point in this example: αDK = 8, n = 2000. The real solution: p∗v DK = 887.5, p∗S = 718.5, µ∗ = 718.5/887.5 = 0.81 The final values of the iteration without measurement: pvDS = 882.9, µ = 0.720. The final values of the iteration with measurement: pv DS = 888.0, µ = 0.743.
pv DK (k + 1) = pv DK (k) − 0.3(pv DK (k)
− pU + c2 (αDK )(pvDK (k)Ψ (µ(k)))2 ). In our study, the values of η(n), A(αDK ) given by the look-up tables have been manipulated up to ±20% to simulate the natural uncertainties in η(n), A(αDK ). It can be observed that also in this case all iterations are stable and converge to the values, which may be different from the real solutions. In that case, the iterations with the integrated measurement (sensor signal pS ) deliver values with smaller convergence errors, although a linearised sensor model has been applied. In Fig. 10, an example of a such iteration is presented. 6. Conclusions In this paper, we have studied the iterative solutions for nonlinear equations and their applications to embedded systems in the control theoretical framework. The first part of this study deals with the analysis of the existing fixed point iteration and Newton’s methods using the stability theory for the sectorbounded Lure’s systems. As a result of this study, sufficient conditions for the convergence of the iterations have been derived, which are expressed in terms of the sector-boundedness and some LMIs. Moreover, estimations of the convergence rate and the boundedness of the convergence errors caused by uncertainties in the nonlinear equations or by the quantisation are provided. The second part of our study addresses the modified iteration methods with a focus on (a) the feedback of the estimation error of the nonlinear equations (b) the integration of sensor
Parameter
Value
Unit
Cfr TU = TS = Tv DK R Vzyl A(αDK ) η(n)
500,000,000 273 + 20 284 0.0014 Look-up table Look-up table
1/(K m2 s2 ) K m2 /(K s2 ) m3 m2 –
Look-up table
–
Ψ
pS pv DK
References Arcak, M., & Kokotovic, P. (2001). Nonlinear observers: a circle criterion design and robustness analysis. Automatica, 37, 1923–1930. Bhaya, A., & Kaszkurewicz, E. (2006). Control perspectives on numerical algorithms and matrix problems. Society for Industrial and Applied Mathematics. Bhaya, A., & Kaszkurewicz, E. (2007). A control-theoretic approach to the design of zero finding numerical methods. IEEE Transactions on Automatic Control, 52, 1014–1026. Boyd, S., Ghaoui, L., & Feron, E. (1994). Linear matrix inequalities in systems and control theory. SIAM. Chen, C. T. (1984). Linear system theory and design. Winston: Holt Rinehart. Eich-Soellner, E., & Führer, C. (1998). Numerical methods in multibody dynamics. B. G. Teubner Stuttgart. Gahinet, P., Nemirovski, A., Laub, A. J., & Chilali, M. (1995). LMI control toolbox user’s guide. The MathWorks, Inc. Guzzella, L., & Onder, C. H. (2010). Introduction to modeling and control of internal combustion engine systems (2nd ed.). Springer. Hoffmann, A., Marx, B., & Vogt, W. (2005). Mathematik für ingenieure. Pearson Studium. Kashima, K., & Yamamoto, Y. (2007). System theory for numerical analysis. Automatica, 43, 1156–1164. Khalil, H. K. (2002). Nonlinear systems. Prentice-Hall. Kiencke, U., & Nielsen, L. (2005). Automotive control systems. Springer-Verlag. Quarteroni, A., Sacco, R., & Saleri, F. (2000). Numerical mathematics. Springer-Verlag. Schaerer, C. E., & Kaszkurewicz, E. (2001). The shooting method for the solution of ordinary differentila equations: a control-theoretical perspective. International Journal of Systems Science, 32, 1047–1053. Scherer, C., Gahinet, P., & Chilali, M. (1997). Multiobjective output-feedback control via LMI optimization. IEEE Transactions on Automatic Control, 42(7), 896–911. Söderlind, G. (2002). Automatic control and adaptive time-stepping. Numerical Algorithms, 31, 281–310. Stoer, J., & Bulirsch, R. (2002). Introduction to numerical analysis. Springer-Verlag. Wang, Z., Liu, Y., & Liu, X. (2008). H-infinity filtering for uncertain stochstic timedelay systems with sector-bounded nonlinearities. Automatica, 44, 1268–1277. Wang, Y., Xie, L., & de Souza, C. E. (1992). Robust control of a class of uncertain nonlinear systems. Systems and Control Letters, 19, 139–149. Weinhold, N. (2007). Einbettung modellgestuetzter fehlerdiagnose in regelungssysteme und Deren Anwenung Fuer die on-board diagnose in Fahrzeugen. Ph.D. thesis. University of Duisburg-Essen.
594
Y. Yang, S.X. Ding / Automatica 48 (2012) 583–594 Ying Yang received her Ph.D. degree in control theory from Peking University, China in 2002. From January 2003 to November 2004, she worked as a Postdoctoral Researcher at Peking University. Since 2005, she has been an associate professor at the Department of Mechanics and Aerospace Engineering, College of Engineering, Peking University. Her research interests include robust and optimal control, nonlinear systems control, numerical analysis, fault detection and fault tolerant systems.
Steven Ding received Ph.D. degree in electrical engineering from the Gerhard-Mercator University of Duisburg, Germany, in 1992. From 1992 to 1994, he was an R&D engineer at Rheinmetall GmbH. From 1995 to 2001, he was a professor of control engineering at the University of Applied Science Lausitz in Senftenberg, Germany, and served as vice president of this university during 1998–2000. He is currently a full professor of control engineering and the head of the Institute for Automatic Control and Complex Systems (AKS) at the University of Duisburg-Essen, Germany. His research interests are model-based and data-driven fault diagnosis, fault tolerant systems, real-time control, and their application in industry with a focus on automotive systems and chemical processes.