STABILITY AND DWELL TIME ANALYSIS OF
SWITCHED TIME DELAY SYSTEMS
a thesis
submitted to the department of electrical and
electronics engineering
and the institute of engineering and sciences
of bilkent university
in partial fulfillment of the requirements
for the degree of
master of science
By
Osman Siraceddin Tapkan
September 2007
I certify that I have read this thesis and that in my opinion it is fully adequate, in scope and in quality, as a thesis for the degree of Master of Science.
Prof. Dr. Hitay ¨OZBAY(Supervisor)
I certify that I have read this thesis and that in my opinion it is fully adequate, in scope and in quality, as a thesis for the degree of Master of Science.
Prof. Dr. ¨Omer MORG ¨UL
I certify that I have read this thesis and that in my opinion it is fully adequate, in scope and in quality, as a thesis for the degree of Master of Science.
Assoc. Prof. Dr. Mehmet ¨Onder EFE
Approved for the Institute of Engineering and Sciences:
Prof. Dr. Mehmet Baray
ABSTRACT
STABILITY AND DWELL TIME ANALYSIS OF
SWITCHED TIME DELAY SYSTEMS
Osman Siraceddin Tapkan
M.S. in Electrical and Electronics Engineering
Supervisor: Prof. Dr. Hitay ¨
OZBAY
September 2007
In this thesis we deal with stability analysis of switched feedback system with time delays. We assume that, at any given time for each “candidate” system a controller is designed and a fixed feedback system is obtained until the next switching instant. We investigate the conservativeness of an LMI-based stability test for the time delay systems. This test is used for the dwell time analysis. After obtaining the limitations of this test, we find the exact bounds of allowable parameters appearing in the LMI-based test, in order to optimize the dwell time. For this purpose we consider simple first order systems and higher order systems separately. We also consider the LQR-based switched feedback controllers with time delays and investigate the effects of weighting matrices Q and R on the dwell time.
Keywords: Switched Time-Delay Systems, Dwell Time, Stability Analysis,
¨
OZET
ANAHTARLAMALI ZAMAN GEC˙IKMEL˙I S˙ISTEMLER˙IN
KARARLILIK VE DURMA ZAMANI ANAL˙IZ˙I
Osman Siraceddin Tapkan
Elektrik ve Elektronik M¨uhendisli¯gi B¨ol¨um¨u Y¨uksek Lisans
Tez Y¨oneticisi:
Prof. Dr. Hitay ¨
OZBAY
Eyl¨ul 2007
Bu tez kapsamında zaman gecikmeli anahtarlamalı geribeslemeli kontrol sis-temlerinin kararlılık analizinden bahsedilmi¸stir. Bilinen herhangi bir zamanda, her “aday” sistem i¸cin bir kontrol birimi tasarlandı˘gı farz edilmi¸s ve bir son-raki anahtarlama anına kadar de˘gi¸smez bir geribeslemeli sistem elde edilmi¸stir. Zaman gecikmeli sistemlerin kararlılı˘gını test eden LMI tabanlı bir testin ko-runumlulu˘gu incelenmi¸stir. Bu test, durma zamanı analizi i¸cin kullanılmaktadır. Bu testin sınırlamaları elde edildikten sonra, durma zamanını eniyile¸stirmek i¸cin testte ge¸cen serbest bırakılabilir parametrelerin kesin sınırları bulunmu¸stur. Bu ama¸cla basit tek dereceli sistemler ve daha y¨uksek dereceli sistemler ayrı ayrı dikkate alınmı¸stır. Aynı zamanda, zaman gecikmeli LQR tabanlı anahtarlamalı geribeslemeli sistemler dikkate alınarak a˘gırlıklandırma matrisleri Q ve R’nin durma zamanı ¨uzerindeki etkileri incelenmi¸stir.
Anahtar Kelimeler: Zaman Gecikmeli Anahtarlamalı Sistemler, Durma Zamanı,
ACKNOWLEDGMENTS
I would like to express my gratitude to Prof. Dr. Hitay ¨Ozbay for his unlim-ited supervision, support and guidance throughout my graduate studies.
I would also like to thank Prof. Dr. ¨Omer Morg¨ul and Assoc. Prof. Dr. Mehmet ¨Onder Efe for accepting to be on my thesis committee and providing useful advice.
I would like to thank H¨useyin Yavuz, C¸ a˘gatayhan C¸ olako˘glu, Meysun Avcı ¨
Ozg¨un and Melih G¨unay from ASELSAN Inc. for their unlimited support.
At last, but not the least, I would like to thank my family for their supports and endless love through my life.
Contents
1 INTRODUCTION 1
1.1 Literature Review . . . 1
1.2 Problem Statement . . . 2 1.3 Contribution and Organization . . . 3
2 PROBLEM DEFINITION AND PRELIMINARY RESULTS 5
2.1 Controller Model . . . 6
2.2 Observer Model (Dual Model for Controller) . . . 7 2.3 Maximum Allowable Delay in LQR Design . . . 8
2.4 Maximum Allowable Delay for Observer Model Determined from the Small Gain Theorem . . . 11
3 STABILITY ANALYSIS FOR DELAY SYSTEMS 14
3.1 Feasibility Analysis of the LMI-Based Test in (3.3) . . . 16
3.2 Conservatism Analysis of the LMI-based Test For First Order Sys-tems . . . 28
4 DWELL TIME ANALYSIS 32
4.1 Dwell Time Analysis for First Order Systems . . . 33
4.2 Minimum Dwell Time For nth Order Systems . . . 36
5 CONCLUSIONS 40
APPENDIX 42
A MATLAB CODE - SECOND ORDER SYSTEM 42
List of Figures
2.1 Plant Model for Feedback Control System . . . 7
2.2 Allowable Delay for the choice of qr . . . 9
2.3 Placement of the real part of the eigenvalues for the choice of q r . 10 2.4 Conservative Analysis . . . 12
2.5 Conservatism Analysis for Allowable Delay found using the Small Gain Theorem . . . 13
3.1 Solution set for the inequality given in (3.7) . . . 17
3.2 Solution set for the inequality given in (3.14) . . . 20
3.3 Solution set for the inequality given in (3.15) . . . 20
3.4 Solution set for the inequality given in (3.18) . . . 22
3.5 Solution set for the inequality given in (3.25) . . . 24
3.6 Solution set for the inequality given in (3.26) . . . 25
3.7 Solution set for the inequality given in (3.29) . . . 27
3.9 Conservativeness Analysis of LMI Based Test for First Order
Sys-tems . . . 31
4.1 Allowable Range for √k2− 1 . . . 34
4.2 Minimum µ versus qr for minimum dwell time . . . 35
4.3 A versus¯ qr for minimum dwell time . . . 36
Chapter 1
INTRODUCTION
This thesis deals with the stability analysis of switched time delay systems. In [21] by using an LMI-based stability approach for each candidate system a dwell time is obtained for stability of the switched system. That is, switched system is stable under arbitrary switching between these stable candidate systems provided that the smallest time interval between these switching instants is greater than a certain dwell time computed in [21]. We also investigate for simple systems how conservative the LMI-based test of [2], and how we can minimize the dwell time.
1.1
Literature Review
The analysis of time-delay systems have attracted attention especially in the last decade [2], [6], [9], [13], [16]. Delays appear in many engineering applications such as information network systems, process control, guidance and navigation. In the literature on time delay systems, stability is analyzed in two ways as delay-independent and delay dependent stability. Most popular approaches for the stability analysis of the delay systems are Lyapunov-like methods based either on the Razumikhin or the Krasovskii technique [2], [6], [9], [13], [16]. A numerical
analysis also provided using the bifurcation theory and a toolbox for MATLAB ”DDE-BIFTOOL” is presented in [3], [4], [5], [17]. Other numerical techniques also available [8], [14], [19], [20].
Switched systems are hybrid systems consisting of a family of continuous-time “candidate” systems and discrete-time logic, i.e. switching signal. Switching systems are used to improve the transient response and to achieve the stability when it is hard to achieve with a single system. Switching control has a various applications areas, e.g. mechanical systems, automotive industry and air traffic control. The stability of the candidate systems, does not always means the stability of the switched system [10]. In [7] it is shown that switching among stable systems results in a stable switched system, provided that the switching is slow on the average. Average dwell time introduced for degree of slowness for the switching process. In [11], it is stated that existence of common Lyapunov functions for each candidate system ensures the arbitrary switching between the candidate systems and a gradient algorithm is supplied to find common Lyapunov functions. Because it is usually hard to find common Lyapunov functions for each candidate systems, piecewise continuous Lyapunov functions are introduced in [15] and [22].
In this thesis, switched time delay systems are investigated regarding the im-provement of the transient response of the modeled system. In [18] switched time delay system is investigated using an extension of common Lyapunov approach, whereas in [21] piecewise Lyapunov Razumikhin functionals are used along with the notion of minimum dwell time.
1.2
Problem Statement
In this thesis we deal with stability analysis of switched feedback systems with time delays. We assume that at any given time for each “candidate” system
a controller is designed and a fixed feedback system is obtained until the next switching instant. On each fixed time intervals between switching times, system is assumed to be in the fixed form:
˙x = Ax(t) + ¯Ax(t − τ ) + Bu(t) (1.1) where τ > 0, u(t) is the input, x(t) is the state variable, A, ¯A, B are appropriate
size matrices. For this system, stability analysis is done using an LMI-based test form [2],[13]. Then using [21] we investigate the smallest dwell time which guarantees stability of the switching system.
1.3
Contribution and Organization
Our contributions can be summarized as follows:
• We investigate the conservativeness of an LMI-based stability test for the
time delay systems. This test is used in [21] for dwell time analysis. There-fore we come up with the limitations of the dwell time analysis. For this purpose we consider simple first order systems to illustrate the level of conservatism.
• We find the exact bounds of allowable parameters appearing in the above
mentioned LMI-based test, for a stable switched time delay system in order to optimize the dwell time.
• We also consider the LQR-based switched feedback controllers with time
delays and investigate the effects of weighting matrices Q and R on the dwell time.
In Chapter 2, we first express feedback control problem with delays in terms of state feedback and state estimate models. We give preliminary results in
this chapter. In Chapter 3 stability conditions of the LMI-based test mentioned above is investigated. The conservatism analysis for this test is provided. In Chapter 4, the results from previous chapters are processed to find a minimum dwell time for a first order control system and second order state estimation system. Furthermore, effects of the weighting matrices Q and R, on the minimum dwell time and ¯A are investigated.
Chapter 2
PROBLEM DEFINITION AND
PRELIMINARY RESULTS
We assume that the switched system consists of ` models. Between switching time instants the system is in the form:
˙xσ(t) = Aσ(t)xσ(t)(t) + Bσ(t)uσ(t)(t)
yσ(t)(t) = Cσ(t)xσ(t)(t − τσ(t)) + Dσ(t)wσ(t)(t) (2.1)
Between each consequent switching time instants ti and ti+1, switching signal
selects one of ` models. Switching signal is described as follows. S : i = σ(t) ∈ {1, 2, ..., `}
The switching signal causes an arbitrary selection between ”candidate” sys-tems. The switching signal design for control purposes is out of this thesis’ scope. In other words we discuss what happens under arbitrary switching, which is de-termined externally or internally but out of our control. In particular, we will be interested in finding a dwell time for stability under arbitrary switching.
Assumptions:
1. There are ` candidate models in the form (2.1) where Ai, Bi, Ci and Di
are fixed matrices, ui is the control input and wi is the noise.
2. The delay, τi > 0, is also assumed to be fixed and known.
3. Feedback system is formed by
ui(t) = −Kiyi(t) + vi(t)
where vi(t) is the disturbance input. Then we can write this system as
˙x(t) = Aix(t) + ¯Aix(t − τ ) + Bivi(t) − KiDiwi(t) (2.2)
where ¯Ai = −BiKiCi.
4. In the scalar case, A > 0 and ¯A < 0, i.e. uncontrolled system is unstable,
and we analyze the effect of Ki stability of each candidate systems.
We will return to this model later. In the rest of this chapter we assume that we have only one system and drop the subscripts.
2.1
Controller Model
Let us consider the simple first order plant model with transfer function for the plant shown in Figure 2.1:
Figure 2.1: Plant Model for Feedback Control System
τ > 0, A > 0.
Writing state-space realization for (2.3) in the form (2.2) gives ˙x(t) = Ax(t) + Bu(t)
u(t) = −Kx(t − τ ) + v(t) (2.4) The closed loop controlled state equation is
˙x(t) = Ax(t) − BKx(t − h) + Bv(t) (2.5) Doing the transformation −BK → ¯A, we can express (2.5) in the form of (2.2)
2.2
Observer Model (Dual Model for
Con-troller)
The state-space model for the typical state estimation problem with delay is in the form:
y(t) = Cx(t − τ ) + Dw(t) (2.6) The observer equation is
˙ˆx(t) = Aˆx(t) + L(y(t) − C ˆx(t − τ)) (2.7)
where L is the Kalman gain matrix. Let the error function be
e(t)= x(t) − ˆ4 x(t)
then
˙e(t) = Ax(t) + Bv(t) − Aˆx(t) − L {C [x(t − τ ) − ˆx(t − τ ) + Dw(t)]} (2.8) ˙e(t) = Ae(t) − LCe(t − τ ) + Bv(t) − LDw(t) (2.9)
Again with the transformation −LC → ¯A, problem can be expressed in the
form of (2.2).
2.3
Maximum Allowable Delay in LQR Design
In order to obtain a stable system, roots of the equation
det(sI − (A + ¯Ae−τ s)) = 0 (2.10)
should be in C−, where ¯A = −LC for observer model and ¯A = −BK for the
controller model.
For the numerical analysis given below we chose the observer design for the standard constant velocity vehicle model where
A = 0 1 0 0 , and C = I2×2
Design of L (or K) for given A and τ = 0
In LQR design define the cost function: J(u) = Z ∞ 0 ¡ xTQx + uTRu¢dt (2.11) Where u(t) = −Kx(t) and K minimizes J.
In order to find L, dual of the LQR problem is used as AT
observer → A,
CT
observer → B, and KT → L.
After we find the observer gain L, we end up with a stable closed-loop system for
τ = 0. The next step is to find the largest allowable delay so that the feedback
system is stable. In other words, we need to find the the minimum de-stabilizing delay.
Let Q = qI2×2 and R = rI2×2. The MATLAB program DDE-BIFTOOL is used to obtain the minimum de-stabilizing delay numerically. The allowable delay and the eigenvalues of the closed loop system changes with the choice of r and
q. This change is illustrated for the observer model in Figure 2.2 and Figure 2.3
10−2 10−1 100 101 100
q/r
Maximum Allowable Feedback Delay
10−2 10−1 100 101 −100
q/r
Real Part of Eigenvalues of Feedbak System
Figure 2.3: Placement of the real part of the eigenvalues for the choice of q r
Observing Figure 2.2 and Figure 2.3 we come up with the following results:
1. As qr increases, eigenvalues of the closed loop system come closer to the imaginary axis
2. Large qr results in a faster transient response however the system becomes more aggressive and less robust to the delay. This conclusion is ensured with the Figure 2.2, as it is seen with increasing ratio of q
r allowable delay
decreases significantly.
If the weighting matrix R in (2.11) (here R is a scalar), is increased then the controller gain K is decreased. With a small K, the it is harder to make the
closed loop system stable. This means the feedback control system can tolerate the smaller delays as given in Figure 2.2.
2.4
Maximum Allowable Delay for Observer
Model Determined from the Small Gain
Theorem
In order to have a stable closed loop system with delay τ , (2.10) must be satis-fied. We can express (2.10) as
det(sI − (A − LCe−τ s)) = 0
⇒
det¡sI − (A − LC) − LC(1 − e−τ s)¢ = 0
⇒
det (sI − (A − LC)) det¡I − C(sI − (A − LC))−1L(1 − e−τ s)¢= 0
The eigenvalues coming from the first part of the above equation are on the left-half plane because we have a stable closed loop system without delay as described earlier. Thus, now we are interested in the eigenvalues of the system shown in Figure 2.4 whose characteristic equation is
det¡I − C (sI − (A − LC))−1L¡1 − e−τ s¢¢= 0 (2.12)
Let
G(s) = C (sI − (A − LC))−1L¡1 − e−τ s¢
According to the Small Gain Theorem closed loop system is stable if
Figure 2.4: Conservative Analysis If we rewrite G(s) as G(s) = C (sI − (A − LC))−1Ls µ 1 − e−τ s s ¶ Since ° ° ° ° µ 1 − e−τ s s ¶° ° ° ° ∞ ≤ τ
to guarantee stability using the Small Gain Theorem we need ° °sC (sI − (A − LC))−1L°° ∞ < 1 τ which is equivalent to τ <°°sC (sI − (A − LC))−1L°°−1∞ . (2.13)
Therefore, maximum allowable delay found from this analysis is the quantity on the right hand side of 2.13. Figure 2.5 shows the conservativeness of (2.13) with respect to the allowable delay found by using DDE-BIFTOOL toolbox.
10−2 10−1 100 101 100
q/r
Allowable Delay
° ° °sC(sI − (A − LC)) −1 L°° ° −1 ∞ Exact BoundFigure 2.5: Conservatism Analysis for Allowable Delay found using the Small Gain Theorem
Chapter 3
STABILITY ANALYSIS FOR
DELAY SYSTEMS
Let us begin with a review of some basic concepts from the Linear Algebra.
Minor: The i × j minor of an n × n matrix, X, denoted |Mij|, is the determinant
of the (n−1)×(n−1) matrix obtained by deleting the ithrow and the jthcolumn
of X.
Leading Principal Minor: The kth order principal leading minor of n × n matrix
X, denoted by |Mk|, is the determinant of the first k rows and columns of X
Theorem: n × n symmetric matrix X is negative definite if and only if
(−1)k|Mk| > 0, k ∈ {1, 2, ..., n}
Now consider the stability test used in [21] taken from [2] ˙x(t) = Ax(t) + ¯Ax(t − τ )
x(θ) = φ(θ), ∀θ ∈ [−τ, 0]
(3.1)
The triplet Σ := (A, ¯A, τ ) ∈ Rn×n× Rn×n× R+is asymptotically stable dependent
of delay if the following lemma holds:
Ω P ¯AM
MTA¯TP −R
where Ω = τ−1£(A + ¯A)TP + P (A + ¯A)¤+ p(α + β)P , M = h A A¯ i , R = diag(αP, βP ),
and α > 0, β > 0 and p > 1 are scalars
P ∈ Rn×n is symmetric positive definite matrix.
If we assume that α, β and p are to be fixed, as α = α∗, β = β∗ and p = p∗,
then (3.2) becomes an LMI whose decision variable is the matrix P . For dwell time analysis given in [21], we need to find feasible solution set for (P, α, β, p) satisfying (3.2). In order to find a feasible set, random and linear searches are done assuming fixed values for α, β and p, searching for positive definite P matrix using LMI-toolbox ([12]) developed for MATLAB. This tests, especially for nth
order systems where n > 1, show us it is very difficult to find a feasible set for (3.2). This lead us to the need for analysis of the conservativeness of the test given in (3.2).
In this section we test the conservatism of the test given in (3.2) on a simple first order system where A and ¯A are scalars.
Thus, applying the theorem to (3.2) using the first order controller model described in Section 1.1, we obtain the following matrix inequality (because all of the variables are scalar and P multiplies each non-zero element, we can erase P from each element):
τ−1£2(A + ¯A)¤+ p (α + β) A ¯A A¯2 A ¯A −α 0 ¯ A2 0 −β < 0 (3.3) As before assume that A > 0 and ¯A < 0
3.1
Feasibility Analysis of the LMI-Based Test
in (3.3)
Preliminaries:
Consider the second order polynomial, with coefficients a, b and c,
P (x) = ax2+ bx + c 1. c
a is the multiplication of the roots of P (x) = 0
2. −b
a is the summation of the roots of P (x) = 0
3. If the discriminant (∆ = b2− 4ac) is negative and a > 0, then the polyno-mial is always positive for all x
4. If the discriminant is negative and a < 0, then the polynomial is always negative for all x
First Leading Principal Minor
According to the theorem the following inequalities must hold first:
2τ−1(A + ¯A) + p(α + β) < 0
or
p(α + β) < −2τ−1(A + ¯A) (3.4)
According to (3.4), because p, α, β and τ are positive, (A+ ¯A) should be negative.
This means:
¯
¯ ¯A¯¯ > A (3.5)
Second Leading Principal Minor
Now checking the second leading principal minor, we should have
Rewriting (3.6) as treating α as the variable of the polynomial,
pα2+£2τ−1(A + ¯A) + pβ¤α + (A ¯A)2 < 0 (3.7) Let us investigate the properties of the inequality given in (3.7):
• Because the coefficient of the α2 term is positive, the discriminant of the polynomial in (3.7) should be greater than 0 (if the discriminant is negative, then the polynomial takes positive values for all α). Thus, the polynomial has two real roots, namely α1 and α2.
• Both the constant term and the coefficient of the term α2 of the polynomial in (3.7) is positive, so is the multiplication of the roots of the polynomial. This means, the roots α1 and α2 are either both negative or both positive.
• By definition, α is positive and according to (3.7) solution set of α lies
between the roots α1 and α2 on the real axis (see Figure 3.1).
Using the above information, we conclude that both of the roots must be positive. To obtain positive roots for the polynomial, the coefficient of the term
α should be negative. It is easily verified that the term
£
2τ−1(A + A) + pβ¤
is negative using (3.4).
Using the inequality given in (3.7) we obtain an upper bound and lower bound for α, as it is depicted in Figure 3.1. Let αlower1 = α1 and αupper1 = α2.
Let us go back to the requirement that the discriminant of (3.7) should be positive. We come to the inequality below:
£ 2τ−1(A + ¯A) + pβ¤2− 4(A ¯A)2p > 0 (3.8) We can express (3.8) as ¡ −£2τ−1(A + ¯A) + pβ¤+ 2A ¯A√p¢× ¡ −£2τ−1(A + ¯A) + pβ¤− 2A ¯A√p¢> 0 (3.9)
Because the second multiplier is always positive, the first one should also be positive ⇒ ¡ −£2τ−1(A + A) + pβ¤+ 2A ¯A√p¢ > 0 ⇒ β < 2A ¯A √ p − 2τ−1(A + ¯A) p . (3.10)
This inequality defines an upper bound for β, namely βupper1. Because β is positive by definition, this bound should also be positive.
2A ¯A√p − 2τ−1(A + ¯A) > 0
⇒
√
p < τ−1(A + ¯A)
A ¯A . (3.11)
By definition p > 1, then the upper bound found in the previous step should be also greater than 1. This gives us the following bound on τ :
(A + ¯A)
A ¯A > τ. (3.12) Third Leading Principal Minor
Let us now check the Third Leading Principal Minor:
£
2τ−1(A + ¯A) + p(α + β)¤αβ −£− ¯A4α − β(A ¯A)2¤ < 0 (3.13)
Rewriting (3.13) as treating β as the variable of the polynomial
pαβ2+£2τ−1α(A + ¯A) + pα2+ (A ¯A)2¤β + α ¯A4 < 0 (3.14)
Let us investigate the properties of the inequality given in (3.14):
• Because the coefficient of the β2 term is positive, the discriminant of the polynomial in (3.14) should be greater than 0 (if the discriminant is nega-tive, then the polynomial takes positive values for all β). Thus, the poly-nomial has two real roots, namely β1 and β2.
• Both the constant term and the coefficient of the term β2 of the polynomial in (3.14) is positive, so is the multiplication of the roots of the polynomial. This means, the roots β1 and β2 are either both negative or both positive.
• By definition, β is positive and according to (3.14) solution set of β lies
between the roots β1 and β2 on the real axis (see Figure 3.2).
Using the above information, we conclude that both of the roots must be positive. To obtain positive roots for the polynomial, the coefficient of the term β should be negative.
Figure 3.2: Solution set for the inequality given in (3.14)
Using the inequality given in (3.14) we obtain an extra upper bound and a lower bound for β, as it is depicted in Figure 3.2. Let βlower1 = β1 and
βupper2 = β2.
The condition for (3.14) to have positive roots is:
pα2+ 2τ−1α(A + ¯A) + (A ¯A)2 < 0 (3.15)
• Because the coefficient of the α2 term is positive, the discriminant of the polynomial in (3.15) should be greater than 0 (if the discriminant is nega-tive, then the polynomial takes positive values for all α). Thus, the poly-nomial has two real roots, namely α3 and α4.
• It is easily verified that both of the roots of the polynomial in (3.15) are
positive.
Figure 3.3: Solution set for the inequality given in (3.15)
Using the inequality given in (3.15) we obtain an extra upper bound and an extra lower bound for α, as it is depicted in Figure 3.3. Let αlower2 = α3 and
αupper2 = α4.
Let us go back to the requirement that the discriminant of the polynomial in (3.15) should be positive. 4τ−2(A + ¯A)2− 4p(A ¯A)2 > 0 ⇒ p < (A + ¯A) 2 τ2(A ¯A)2 (3.16)
Notice that we conclude with the same upper bound for p in (3.11)
Let us go back to the requirement that the discriminant of (3.14) should be positive, we come with the inequality below:
£ 2τ−1α(A + ¯A) + pα2 + (A ¯A)2¤2− 4α2p ¯A4 > 0 ⇒ ¡ −£2ατ−1(A + ¯A) + pα2+ (A ¯A)2¤− 2α ¯A2√p¢× ¡ −£2ατ−1(A + ¯A) + pα2+ (A ¯A)2¤+ 2α ¯A2√p¢> 0 (3.17) Because the second multiplier is always positive, the first one should also be positive
⇒
pα2+£2 ¯A2√p + 2τ−1(A + ¯A)¤α + (A ¯A)2 < 0 (3.18)
Let us investigate the properties of the inequality given in (3.18):
• Because the coefficient of the α2 term is positive, the discriminant of the polynomial in (3.18) should be greater than 0 (if the discriminant is nega-tive, then the polynomial takes positive values for all α). Thus, the poly-nomial has two real roots, namely α5 and α6.
• Both the constant term and the coefficient of the term α2 of the polynomial in (3.18) is positive, so is the multiplication of the roots of the polynomial. This means, the roots α5 and α6 are either both negative or both positive.
• By definition, α is positive and according to (3.18) solution set of α lies
between the roots α5 and α6 on the real axis (see Figure 3.4).
Using the above information, we conclude that both of the roots must be positive. To obtain positive roots for the polynomial, the coefficient of the term α should be negative.
Figure 3.4: Solution set for the inequality given in (3.18)
Using the inequality given in (3.18) we obtain an extra upper bound and an extra lower bound for α, as it is depicted in Figure 3.4. Let αlower3 = α5 and
αupper3 = α6.
The condition for (3.18) to have positive roots is:
2τ−1(A + ¯A) + 2√p ¯A2 < 0
⇒
√
p < −(A + ¯A)
τ ¯A2 . (3.19)
By definition p > 1, then the upper bound found in the previous step should be also greater than 1. This gives us the following bound on τ :
(A + ¯A)
¯
A2 > τ. (3.20) Let us go back to the requirement that the discriminant of (3.18) should be positive, we come to the inequality below:
£ 2 ¯A2√p + 2τ−1(A + ¯A)¤2− 4p(A ¯A)2 > 0 ⇒ ¡ −£2 ¯A2√p + 2τ−1(A + ¯A)¤− 2A ¯A√p¢× ¡ −£2 ¯A2√p + 2τ−1(A + ¯A)¤+ 2A ¯A√p¢ (3.21)
Because the first multiplier is always positive, the second one should also be positive ⇒ −2 ¯A2√p − 2τ−1(A + ¯A) + 2A ¯A√p > 0 (3.22) ⇒ √ p < −(A + ¯A) τ ( ¯A2− A ¯A) (3.23)
By definition p > 1, then the upper bound found in the previous step should be also greater than 1. This gives us the following bound on τ :
−(A + ¯A)
( ¯A2 − A ¯A) > τ (3.24)
Rewriting (3.13) treating α as the variable of the polynomial, we find new bounds for β, p and α:
pβα2+£2τ−1β(A + ¯A) + pβ2+ ( ¯A)4¤α + β(A ¯A)2 < 0 (3.25)
• Because the coefficient of the α2 term is positive, the discriminant of the polynomial in (3.25) should be greater than 0 (if the discriminant is nega-tive, then the polynomial takes positive values for all α). Thus, the poly-nomial has two real roots, namely α7 and α8.
• Both the constant term and the coefficient of the term α2 of the polynomial in (3.25) are positive, so is the multiplication of the roots of the polynomial. This means, the roots α7 and α8 are either both negative or both positive.
• By definition, α is positive and according to (3.25) solution set of α lies
between the roots α7 and α8 on the real axis (see Figure 3.5).
Using the above information, we conclude that both of the roots must be positive. To obtain positive roots for the polynomial, the coefficient of the term α should be negative.
Figure 3.5: Solution set for the inequality given in (3.25)
Using the inequality given in (3.25) we obtain an extra upper bound and an extra lower bound for α, as it is depicted in Figure 3.5. Let αlower4 = α7 and
αupper4 = α8.
The condition for (3.25) to have positive roots is:
• Because the coefficient of the β2 term is positive, the discriminant of the polynomial in (3.26) should be greater than 0 (if the discriminant is nega-tive, then the polynomial takes positive values for all β). Thus, the poly-nomial has two real roots, namely β3 and β4.
• It is easily verified that both of the roots of the polynomial in (3.26) are
positive.
Figure 3.6: Solution set for the inequality given in (3.26)
Using the inequality given in (3.26) we obtain an extra upper bound and an extra lower bound for β, as it is depicted in Figure 3.6. Let βlower2 = β3 and
βupper3 = β4.
Let us go back to the requirement that the discriminant of (3.26).
4τ−2(A + ¯A)2− 4p ¯A4 > 0 ⇒
p < (A + ¯A)
2
τ2A¯4 (3.27)
Notice that we conclude with the same upper bound for p in (3.19).
Let us go back to the requirement that the discriminant of (3.25) should be positive, we come with the inequality below:
£
⇒
¡
−£2βτ−1(A + ¯A) + pβ2+ ¯A4¤− 2β(A ¯A)√p¢× ¡
−£2βτ−1(A + ¯A) + pβ2+ ¯A4¤+ 2β(A ¯A)√p¢> 0 (3.28) Because the first multiplier is always positive, the second one should also be positive.
⇒
pβ2+£2(A ¯A)√p + 2τ−1(A + ¯A)¤β + ¯A4 < 0 (3.29)
Let us investigate the properties of the inequality given in (3.29):
• Because the coefficient of the β2 term is positive, the discriminant of the polynomial in (3.29) should be greater than 0 (if the discriminant is nega-tive, then the polynomial takes positive values for all β). Thus, the poly-nomial has two real roots, namely β5 and β6.
• Both the constant term and the coefficient of the term β2 of the polynomial in (3.29) are positive, so is the multiplication of the roots of the polynomial. This means, the roots β5 and β6 are either both negative or both positive.
• By definition, β is positive and according to (3.29) solution set of β lies
between the roots β5 and β6 on the real axis (see Figure 3.7).
Using the above information, we conclude that both of the roots must be posi-tive. To obtain positive roots for the polynomial, the coefficient of the term β should be negative.
Using the inequality given in (3.29) we obtain an extra upper bound and an extra lower bound for β, as it is depicted in Figure 3.7. Let βlower3 = β5 and
Figure 3.7: Solution set for the inequality given in (3.29) 2τ−1(A + ¯A) + 2√p(A ¯A) < 0 ⇒ √ p < (A + ¯A) τ (A ¯A) (3.30)
Notice that we conclude with the same upper bound for p in (3.11).
Let us go back to the requirement that the discriminant of (3.29) should be pos-itive, we come to the inequality below:
£ 2(A ¯A)√p + 2τ−1(A + ¯A)¤2− 4p ¯A4 > 0 ⇒ ¡ −£2(A ¯A)√p + 2τ−1(A + ¯A)¤− 2 ¯A2√p¢× ¡ −£2(A ¯A)√p + 2τ−1(A + ¯A)¤+ 2 ¯A2√p¢ (3.31) Because the second multiplier is always positive, the first one should also be positive ⇒ −2A ¯A√p − 2τ−1(A + ¯A) − 2 ¯A2√p > 0 (3.32) ⇒ √ p < −(A + ¯A) τ ( ¯A2+ A ¯A) (3.33)
By definition p > 1, then the upper bound found in the previous step should be also greater than 1. This gives us the following bound on τ :
−(A + ¯A)
( ¯A2+ A ¯A) > τ (3.34)
3.2
Conservatism Analysis of the LMI-based
Test For First Order Systems
Let’s define the delay system as˙x(t) = Ax(t) + ¯Ax(t − τ ) (3.35) where τ is the delay introduced to the system.
Assume that; A > 0 and ¯A < 0
Let’s define − ¯A = kA, k > 1 according to the stability condition given in(3.5).
We can represent (3.35) in Laplace domain as
s − A + kAe−τ s = 0
⇒
1 + kAe
−τ s
s − A = 0 (3.36)
Applying Nyquist Criteria, we require a diagram similar to Figure 3.8 (i.e. the point (−1 + j0) should be encircled once in the counter clock-wise direction).
To achieve this, at the cross-over frequency ωc, the following phase condition
should be met: −π < −τ ωc− ³ π − tan−1³ωc A ´´
Figure 3.8: Nyquist diagram for stable first order delay system ⇒ tan−1³ωc A ´ > τ ωc (3.37) where ¯ ¯ ¯ ¯jω−kA c− A ¯ ¯ ¯ ¯ = 1 ⇒ ωc= √ k2A2− A2 = A√k2− 1 In order to define a bound on τ A, (3.37) can be expressed as
tan−1¡√k2− 1¢
√
k2− 1 > τ A , (3.38) which is shown in Figure 3.9 as the exact bound.
Using (3.12) we found a similar bound on τ A. If we express (3.12) as 1 + A¯ A ¯ A A > τ A,
and use the same definition for k, where k > 1, we find a bound on τ A as
k − 1
k > τ A. (3.39)
Using (3.20) we found another bound on τ A with the same definition for k:
k − 1
k2 > τ A (3.40) Using 3.24 we find the last bound on τ A:
k − 1 k2(1 + 1
k)
> τ A (3.41)
Finally, using 3.34 we find the last bound on τ A:
k − 1 k2(1 − 1
k)
> τ A (3.42)
Among all of the four bounds given in (3.39), (3.40), (3.41) and (3.42), the one in (3.41) is most conservative one. It is shown in Figure 3.9 as the conservative bound. The level of conservativeness of the LMI-based test given in (3.1) can be viewed in Figure 3.9.
In particular, Figure 3.9 shows that, for example when τ A = 0.2 we cannot find a solution using (3.2) (yet for this case there exists an ¯A = −kA with
1 < √k2− 1 < 7.2 such that the feedback system is stable. For the values of
10−1 100 101 102 103 10−3 10−2 10−1 100 √k2 −1 M a x im u m Al lo w a b le τ A Exact Bound Conservative Bound
Chapter 4
DWELL TIME ANALYSIS
Let us begin with the results on dwell time obtained in [21]. The switched delay system consists of ` triplets as Σi := (Ai, ¯Ai, τi), where i ∈ {1, 2, ..., `}. The
switched time delay system is asymptotically stable if all triplets are asymptoti-cally stable. The following definitions are provided:
Si := −{Pi(Ai+ ¯Ai) + (Ai + ¯Ai)TPi +τiα−1PiA¯iAiPi−1ATi A¯Ti Pi +τiβ−1Pi( ¯Ai)2P−1( ¯AiT)2Pi+ τipi(αi+ βi)Pi} (4.1) κi := σmin[Pi] ¯ κi := σmax[Pi] wi := σmin[Si] λ := maxi κκ¯ii µ := maxiwκ¯ii
Then the dwell time τD is defined as
where
T∗ = λµbλ − 1 ¯
p − 1 + 1c (4.3)
and ¯p := mini{pi}, τmax = maxiτi
In (4.1), the scalars α, β, p and the matrix P are chosen to satisfy the test given in (3.2).
Lemma[21]: Switched time delay system is stable under arbitrary switching if the difference between consecutive switching time instants is strictly grater than the dwell time τD.
4.1
Dwell Time Analysis for First Order
Sys-tems
Advantages of the first order system :
1. Simple model
2. Dwell time analysis is reduced to analyze µ parameter
3. Because the problem is reduced to finding an optimum (in this case the minimum) µ, analysis can be done for each candidate system separately, the one with maximum dwell time will dominate the overall dwell time of the system.
Because P will be scalar with a first order system, λ = 1. Design parameters are only included in T∗we can focus on this parameter. With the fact that λ = 1,
T∗ is reduced to
10−1 100 101 102 103 10−3 10−2 10−1 100
√
k
2− 1
M a x im u m A ll o w a b le τ A Exact Bound Conservative BoundAllowable interval for √k2
− 1 for given τ A
Figure 4.1: Allowable Range for √k2 − 1
⇒ µ = ¯¯ 1 ¯2(A + A) + τ α−1(AA)2+ τ β−1A4+ τ p(α + β) ¯ ¯ ¯
According to (3.2) S is always negative. Thus
T∗ = −1
2(A + A) + τ α−1(AA)2+ τ β−1A4+ τ p(α + β) (4.4)
For the numerical analysis A is chosen as 0.1. For each fixed (τ A) there exists a range√k2− 1, i.e. ¯A = −kA, feedback system is stable as shown in Figure 4.1.
Now assume that ¯A is obtained from LQR design. For different values of q r
we obtain different ¯A. For each ¯A, maximum allowable delay (τmax) is obtained
using Figure 4.1. Then we select τ = τmax
10 . The boundaries for α, β, p and ¯A are used to find the minimum dwell time with the help of MATLAB Optimization
Toolbox. The results are given with Figure 4.2.
0 2 4 6 8 10 0 0.5 1 1.5
q/r
µ
Figure 4.2: Minimum µ versus qr for minimum dwell time
In Figure 4.2, we observe that with increasing ratio of qr the minimum dwell time is decreased. Recall that in Figure 2.2, the increasing ratio of qr causes increased maximum allowable delay due to the increased robustness of the system to the delay type disturbance. Here, we can assume that the decreased dwell time shows us the degree of robustness of the switched time delay system.
10−1 100 101 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0
q/r
¯ A
Optimum curve for dwell time minimization problem
Optimum curve for LQR design
Figure 4.3: ¯A versus qr for minimum dwell time
In Figure 4.3, we observe that the optimum value of ¯A to minimize the dwell
time is different from ¯A found initially to minimize the cost function of LQR
design in (2.11).
4.2
Minimum Dwell Time For n
thOrder
Sys-tems
In Chapter 3 we illustrated how conservative the LMI-test given in (3.2) for some choice of τ A. Furthermore the conservativeness analysis requires complicated calculations even for first order systems. Thus, we regard some extra assumptions in order to analyze the dwell time characteristics for stable switched time delay systems.
Let us assume that P = γIn×n and p is a constant, say p∗. Rewriting 3.2 results in: X : τ−1γ£(A + ¯A) + (A + ¯A)T¤+ pγ (α + β) γ ¯AA γ ¯A2 γATA¯T −γαI n×n 0n×n γ( ¯AT)2 0 n×n −γβIn×n < 0 (4.5)
With this selection of P and p∗,
• λ = 1. As in the scalar case, the problem of finding minimum dwell time
for given system is reduced to find the minimum µ problem.
• 3.2 becomes an LMI.
• Because S in (4.1) is positive definite matrix, singular values σk, k ∈
{1, 2, ..., n} of S are equal to the eigenvalues λk, µ is independent of γ
• Because γ multiplies each non-zero term in (4.7), γ can be canceled out
and the overall problem becomes independent of choice of γ.
Now, the problem of ”minimizing dwell time” can be expressed as ”maximize the smallest singular value (in this case the minimum eigenvalue) of the matrix
S”. Considering the variables α and β as the decision variables, we can state the
problem as,
maximize z
subject to
X < 0 (4.6)
α > 0, β > 0
We can express the nonlinear constraint S − zI in terms of LMI using Schur Complement property ([1]) as:
τ−1γ£(A + ¯A) + (A + ¯A)T − zI¤+ pγ (α + β) γ ¯AA γ ¯A2 γATA¯T −γαI n×n 0n×n γ( ¯AT)2 0 n×n −γβIn×n < 0 (4.7)
For the numerical analysis, we chose the observer design for the standard constant velocity vehicle model where
A = 0 1 0 0 , and C = I2×2. ¯
A is obtained using dual of LQR problem along with the transformation −LC →
¯
A. Using DDE-BIFTOOL maximum allowable delay (τmax) is computed. Tests
for finding minimum dwell time, showed that for only smaller delays than the al-lowable maximum delay, a feasible solution for the problem given in (4.7) exists. This reminds us the conservativeness of the test in (3.2).
In Figure 4.4, change of the parameter µ with respect to the parameter p is depicted. For values greater than approximately 3.7, the solution becomes in-feasible which means we have upper bound for p as in first order system. Note that, as p decreases, µ also decreases.
1 1.5 2 2.5 3 3.5 4 0 2 4 6 8 10 12 14 16 18 20
p
µ
Chapter 5
CONCLUSIONS
In this thesis, we investigate the stability requirements of switched time delay systems. The stability analysis is basically done in terms of dwell-time.
First, we analyze linear delay systems. For this analysis, we used first or-der feedback control systems and second oror-der constant velocity vehicle models. At first step we determine the conservativeness of our techniques on the stabil-ity with respect to the maximum allowable delays found using DDE-BIFTOOL toolbox.
Second, we investigated the limitations of the LMI-based stability test given in (3.2) for a first order feedback control system. We found the bounds for the stability in terms of α, β and p. Using the upper bounds found for p, we built the relationship with the stability analysis of first order delay system using Nyquist criteria and outlined the conservativeness of the LMI-based test. We found that for small τ A values the test is quite conservative.
In the final chapter we gave numerical results for dwell time using the results of the previous chapters for first order controller model and second order constant velocity vehicle observer model. Due to the lack of stability analysis for nth order
systems where n > 1, we could give the minimum dwell time analysis with respect to the scalar p. The tests showed us that for maximum delay found by DDE-BIFTOOL the minimization problem is infeasible as in first order system. For a feasible problem, the maximum allowable delay was decreased, which led us to interrogate the conservativeness of the test given in (3.2) as an open problem.
For first order systems, dwell time minimization problem was solved as as-suming ¯A as a variable just like α, β and p. In Figure 4.3, it is shown that the
optimum ¯A for minimum dwell time is different from the value of ¯A which
min-imizes the cost function of LQR design given in (2.11). This bring us to a trade off between optimum minimum dwell time problem and LQR problem, which is again an open problem to be tackled down. The same situation is valid for nth
APPENDIX A
MATLAB CODE - SECOND
ORDER SYSTEM
search.m
%This script searches for the maximum allowable delay %iteratively using bisection method
%Searches until the maximum real part of the roots %becomes closest possible to the imaginary axis on the %left half plane
lqrDes; %First construct the matrices from LQR design mx=3; mn=0; nm=0.5; while mx-mn > 0.05 roots_=dri(A,L,C,nm); if roots_<0 mn=nm; else mx=nm;
end nm=(mx+mn)/2; end nm=nm*1000; nm=ceil(nm); nm=nm/1000; roots_=dri(A,L,C,nm); if roots_<0 nm=nm+0.001; while roots_<0 roots_=dri(A,L,C,nm); nm=nm+0.001; end roots_=nm-0.002 else nm=nm-0.001; while roots_>=0 roots_=dri(A,L,C,nm); nm=nm-0.001; end roots_=nm+0.001 end lqrDes.m
%Finds L (or K matrix) using LQR design A=[ 0 , 1; ... 0 , 0]; C=eye(2); R=r*eye(2); Q=q*eye(2); sys_=ss(A’,C’,zeros(2),0);
L=lqr(sys_,Q,R); L=L’;
dri.m
%For given system parameters and delay(tau), returns %the maximum real part of the roots of the delayed system
function [roots_]=dri(A,L,C,tau) LC=L*C; stst.kind=’stst’; stst.parameter=[A(1,1) A(1,2) LC(1,1) LC(1,2)...
A(2,1) A(2,2) LC(2,1) LC(2,2) tau]; stst.x=[0 0]’; method=df_mthod(’stst’); [stst,success]=p_correc(stst,[],[],method.point); stst.x; stst.stability=p_stabil(stst,method.stability); roots_=max(real(stst.stability.l0)); figure(1); clf; p_splot(stst); sys_init.m
%Initialize the delayed system in order to use the
%DDE_BIFTOOL toolbox, declare the name and the dimensions function [name,dim]=sys_init()
name=’max_delay’; dim=2;
%path for the DDE_BIFTOOL toolbox files, i.e. .../ddebiftool path(path,’C:\Documents and Settings\...’);
sys_rhs.m
%The right hand side of the delayed system
%PAR contains the parameters including delay, XX contains the present and the %past states (here the states are the error driven from
%observer and state equations) function f=sys_rhs(xx,par)
% PAR: [ A11 A12 LC11 LC12 A21 A22 LC21 LC22 tau ] % XX : [ e1(t) e1(t-tau) ; e2(t) e2(t-tau) ]
f(1,1)= par(1) * xx(1,1) + par(2) * xx(2,1) ... - par(3) * xx(1,2)-par(4) * xx(2,2); f(2,1)= par(5) * xx(1,1) + par(6) * xx(2,1) ... -par(7) * xx(1,2) - par(8) * xx(2,2); return; sys_deri.m
%Defines the first order partial derivatives wrt parameters function J=sys_deri(xx,par,nx,np,v)
% PAR: [ A11 A12 LC11 LC12 A21 A22 LC21 LC22 tau ] % XX : [ e1(t) e1(t-tau) ; e2(t) e2(t-tau) ]
J=[];
if length(nx)==1 & length(np)==0 & isempty(v) % first order derivatives wrt state variables if nx==0 % derivative wrt x(t)
J(1,1)=par(1); J(1,2)=par(2); J(2,1)=par(5);
J(2,2)=par(6);
elseif nx==1 % derivative wrt x(t-tau1) J(1,1)=-par(3); J(1,2)=-par(4); J(2,1)=-par(7); J(2,2)=-par(8); end; end; if isempty(J) err=[nx np size(v)]
error(’SYS_DERI: requested derivative could not be computed!’); end;
return;
sys_tau.m
%Declares the order of the delay term in parameters vector function tau=sys_tau()
% PAR: [ A11 A12 LC11 LC12 A21 A22 LC21 LC22 tau ] tau=[9];
return;
minimize_mu.m
%This script minimize the mu parameter, assuming alpha
%and beta parameter are the decision variables, by maximizing
%minimum singular value of the S matrix. Searches minimum mu parameter %for different p values
P=gamma*eye(2); counter=1;
for p=1.001:0.001:10
[alpha,beta,z]=minc(A,Abar,p,tau); S = (A+Abar)’*P + P*(A+Abar) + tau ...
* ( 1/alpha * P * Abar * A * inv(P) * A’ * Abar’ * P... + 1/beta * P * Abar^2 * inv(P) * (Abar’)^2 * P + p ... * (alpha+beta) * P ); S=-S; svd1=svd(P); kappa=min(svd1); kappa_bar=max(svd1); w=min(svd(S)); mu=kappa_bar/w; result(counter,1)=p; result(counter,2)=mu; counter=counter+1; end minc.m
%Used by the driver script minimize_mu.m script
%Maximize the objective function for the constraints outlined %in Chapter 4 of the thesis.
function [alpha,beta,z]=minc(A,Abar,p,tau); setlmis([]);
alpha=lmivar(1,[1 1]); beta=lmivar(1,[1 1]);
z=lmivar(1,[1 1]); %Parameter to be maximized
lmiterm([1 1 1 alpha],p*eye(2),eye(2)); lmiterm([1 1 1 beta],p*eye(2),eye(2)); lmiterm([1 1 1 z],(1/tau)*eye(2),eye(2)); lmiterm([1 2 1 0],-A’*Abar’); lmiterm([1 3 1 0], -Abar’*Abar’); lmiterm([1 2 2 alpha],-eye(2),eye(2)); lmiterm([1 3 2 0], zeros(2,2)); lmiterm([1 3 3 beta],-eye(2),eye(2)); lmiterm([-2,1,1,alpha],1,1); %0<alpha lmiterm([-3,1,1,beta],1,1); %0<beta lmiterm([4 1 1 0], (1/tau)*(A+Abar)’+(1/tau)*(A+Abar)); lmiterm([4 1 1 alpha],p*eye(2),eye(2)); lmiterm([4 1 1 beta],p*eye(2),eye(2)); lmiterm([4 2 1 0], A’*Abar’); lmiterm([4 3 1 0], Abar’*Abar’); lmiterm([4 2 2 alpha],-eye(2),eye(2)); lmiterm([4 3 2 0], zeros(2,2)); lmiterm([4 3 3 beta],-eye(2),eye(2)); lmis = getlmis; c=[0;0;-1]; [copt,xopt]=mincx(lmis,c); alpha = dec2mat(lmis,xopt,alpha); beta = dec2mat(lmis,xopt,beta); z = dec2mat(lmis,xopt,z);
APPENDIX B
MATLAB CODE - FIRST
ORDER SYSTEM
minimize_dwellTime.m
%This script is used to investigate the effect of Q and %R in LQR design. Finds initial Abar value and fixed delay, %then minimizes the dwell time
global A_ global tau counter=1; for k=0.1:0.1:10 counter if counter==5 counter; end lqr_des search Abar=-L*C; A_=A; tau=delay/10; driver_min result(counter,1)=k;
result(counter,2)=-L; result(counter,3:6)=sol_’; result(counter,7)=fval; result(counter,8)=tau; counter=counter+1; end lqr_des.m
%LQR design for first order controller A_=0.1; C=1; R=1; Q=k; L=lqr(A_’,C’,Q,R); sys_rhs.m
%The right hand side of the delayed system
%PAR contains the parameters including delay, XX contains the present and the %past states (here the states are the error driven from
%observer and state equations) function f=sys_rhs(xx,par) % PAR: [ A11 LC11 tau ] % XX : [ e1(t) e1(t-tau) ]
f(1,1)= par(1) * xx(1,1) - par(2) * xx(1,2); return;
sys_deri.m
%First order partial derivatives are defined function J=sys_deri(xx,par,nx,np,v)
% PAR: [ A11 LC11 tau ] % XX : [ e1(t) e1(t-tau) ] J=[];
if length(nx)==1 & length(np)==0 & isempty(v) % first order derivatives wrt state variables if nx==0 % derivative wrt x(t)
J(1,1)=par(1);
elseif nx==1 % derivative wrt x(t-tau1) J(1,1)=-par(2);
end; end;
if isempty(J)
err=[nx np size(v)]
error(’SYS_DERI: requested derivative could not be computed!’); end;
return;
driver_min.m
%Minimize the mu parameter using the constraints for %alpha, beta, p and Abar defined in nonlcon1.m
global A_ global tau A_=A; clear A
’MaxIter’,1000000,’TolCon’,0.00001); bounds_Abar Abar_ust=min([Abar1,Abar3,Abar5,Abar6]); Abar_alt=max([Abar2,Abar4,Abar7]); Abar_init=(Abar_alt+Abar_ust)/2; Abar=Abar_init; p_bounds p_ust=min([p1,p2,p3,p4]); p_alt=1; p_init=(p_ust+p_alt)/2; p=p_init; bounds_beta beta_ust=min([beta1,beta3,beta5,beta7]); beta_alt=max(beta2,beta4); beta_init=(beta_alt+beta_ust)/2; beta=beta_init; bounds_alpha alpha_ust= min([alpha1,alpha3,alpha5,alpha7,alpha9]); alpha_alt=max([alpha2,alpha4,alpha6,alpha8]); alpha_init=(alpha_ust+alpha_alt)/2;
x0=[p_init alpha_init beta_init Abar_init];
[sol_,fval,exitflag,output] = fmincon(@myfun,x0,... [],[],[],[],[],[],@nonlcon1,options); p=sol_(1); alpha=sol_(2); beta=sol_(3); Abar=sol_(4);
%Checks the feasibility of the found parameters treating %P as the decision variable.
[P,flag,tmin,lhs1,rhs1]=findDelay_func_delayDepend... (A_,Abar,p,alpha,beta,tau);
myfun.m %Objective function function [mu]=myfun(x) global A_ global tau %x:[p;alpha;beta;Abar] p=x(1); alpha=x(2); beta=x(3); Abar=x(4); mu=-(2*(A_+Abar)+tau/alpha*(Abar*A_)^2+tau/beta*(Abar)^4... + tau*p*(alpha+beta)); mu=1/mu; nonlcon1.m
%nonlinear constraint function function [c,ce] = nonlcon1(x) ce=[]; global A_ global tau %x:[p;alpha;beta;Abar] p=x(1); beta=x(3); alpha=x(2); Abar=x(4);
p_bounds p_ust=min([p1,p2,p3,p4]); p_alt=1; p_step=(p_ust-p_alt)/100; p_ust=p_ust-p_step; p_alt=p_alt+p_step; bounds_beta beta_ust=min([beta1,beta3,beta5,beta7,beta8]); beta_alt=max([beta2,beta4,beta9]); beta_step=(beta_ust-beta_alt)/100; beta_ust=beta_ust-beta_step; beta_alt=beta_alt+beta_step; bounds_alpha alpha_ust=min([alpha1,alpha3,alpha5,alpha7,alpha9]); alpha_alt=max([alpha2,alpha4,alpha6,alpha8]); alpha_step=(alpha_ust-alpha_alt)/100; alpha_ust=alpha_ust-alpha_step; alpha_alt=alpha_alt+alpha_step; bounds_Abar Abar_ust=min([Abar1,Abar3,Abar5,Abar6]); Abar_alt=max([Abar2,Abar4,Abar7]); Abar_step=(Abar_ust-Abar_alt)/100; Abar_ust=Abar_ust-Abar_step; Abar_alt=Abar_alt+Abar_step; %Constraints stored in c vector c(1)=p_alt-x(1); c(2)=-p_ust+x(1); c(3)=alpha_alt-x(2); c(4)=-alpha_ust+x(2); c(5)=beta_alt-x(3); c(6)=-beta_ust+x(3);
c(7)=Abar_alt-x(4); c(8)=-Abar_ust+x(4); p_bounds.m p1=(Abar+A_)^2/tau^2/(Abar*A_)^2; p2=(Abar+A_)^2/tau^2/(Abar)^4 ; p3=(Abar+A_)^2/(tau^2*(Abar^2-A_*Abar)); p4=(Abar+A_)^2/(tau^2*(Abar^2+A_*Abar)); bounds_alpha.m disc1=((2*(A_+Abar)/tau+p*beta))^2- 4*(A_*Abar)^2*p; alpha1=(-((2*(A_+Abar)/tau+p*beta))+sqrt(disc1))/(2*p); alpha2=(-((2*(A_+Abar)/tau+p*beta))-sqrt(disc1))/(2*p); disc2=(2*(A_+Abar)/tau)^2- 4*(A_*Abar)^2*p; alpha3=(-((2*(A_+Abar))/tau)+sqrt(disc2))/(2*p); alpha4=(-((2*(A_+Abar))/tau)-sqrt(disc2))/(2*p); disc3=(2*(A_+Abar)*beta/tau+p*beta^2+(Abar)^4)^2 -4... * (Abar*A_)^2* beta^2*p; alpha5=(-(2*(A_+Abar)*beta/tau+p*beta^2+(Abar)^4)... +sqrt(disc3))/(2* p*beta); alpha6=(-(2*(A_+Abar)*beta/tau+p*beta^2+(Abar)^4)... -sqrt(disc3))/(2* p*beta); disc4=(
2*Abar^2*sqrt(p)+2/tau*(A_+Abar))^2 - 4* p*(A_*Abar)^2;
alpha7= (-(2*Abar^2*sqrt(p)+2/tau*(A_+Abar) )+sqrt(disc4)) / (2*p); alpha8= (-(2*Abar^2*sqrt(p)+2/tau*(A_+Abar) )-sqrt(disc4)) / (2*p); alpha9=-2*(Abar+A_)/tau/p-beta;
bounds_beta.m disc_beta1=(2*(A_+Abar)/tau)^2- 4* Abar^4*p; beta1=(-((2*(A_+Abar))/tau)+sqrt(disc_beta1))/(2*p); beta2=(-((2*(A_+Abar))/tau)-sqrt(disc_beta1))/(2*p); disc_beta2=(2*(A_*Abar)*sqrt(p)+2/tau*(A_+Abar))^2-... 4 * p*Abar^4; beta3=(-(2*(A_*Abar)*sqrt(p)+2/tau*(A_+Abar))+... sqrt(disc_beta2))/(2*p); beta4=(-(2*(A_*Abar)*sqrt(p)+2/tau*(A_+Abar))-... sqrt(disc_beta2))/(2* p); beta5=(2*A_*Abar*sqrt(p)-2/tau*(A_+Abar))/p; beta7=-2*(A_+Abar)/tau/p; if exist(’alpha’,’var’)==1 disc_beta3=(2*(A_+Abar)*alpha/tau+p*alpha^2+(Abar*A_)^2)^2... -4* (Abar)^4 * alpha^2*p; beta8= (-(2*(A_+Abar)*alpha/tau+p*alpha^2+(Abar*A_)^2)+... sqrt(disc_beta3))/ (2 * p*alpha); beta9= (-(2*(A_+Abar)*alpha/tau+p*alpha^2+(Abar*A_)^2)-... sqrt(disc_beta3))/ (2 * p*alpha); end bounds_Abar.m delta_Abar1=(1-A_*tau)^2-4*A_*tau; Abar1=(-(1-A_*tau)+sqrt(delta_Abar1))/(2*tau); Abar2=(-(1-A_*tau)-sqrt(delta_Abar1))/(2*tau); delta_Abar2=(1+A_*tau)^2-4*A_*tau; Abar3=(-(1+A_*tau)+sqrt(delta_Abar2))/(2*tau);
Abar4=(-(1+A_*tau)-sqrt(delta_Abar2))/(2*tau); Abar5=-A_/(1-tau*A_); delta_Abar3=1-4*tau*A_;
Abar6=(-1+sqrt(delta_Abar3))/(2*tau); Abar7=(-1-sqrt(delta_Abar3))/(2*tau);
findDelay_func_delayDepend.m
%Checks the feasibility of the given parameters, treating P %as the decision variable
function [P,flag,tmin,lhs1,rhs1]=findDelay_func_delayDepend... (A_,Abar,p,alpha,beta,tau);
setlmis([]);
P=lmivar(1,[1 1]);
lmiterm([1 1 1 P], (A+Abar)’, (1/tau),’s’); lmiterm([1 1 1 P],p*(alpha+beta),1); lmiterm([1 2 1 P], A’*Abar’, 1); lmiterm([1 3 1 P], Abar’*Abar’, 1); lmiterm([1 2 2 P], (-1)*alpha, 1); lmiterm([1 3 2 0], 0); lmiterm([1 3 3 P], (-1)*beta,1); lmiterm([-2,1,1,P],1,1); %0<P lmis=getlmis; [tmin,xfeas] = feasp(lmis,[0,0,-1,0,0]); if tmin<0 %feasible solution exist
P = dec2mat(lmis,xfeas,P); evals = evallmi(lmis,xfeas); [lhs1,rhs1] = showlmi(evals,1); else flag=0; P=0;
Bibliography
[1] Boyd S., L. E. Ghaoui, E. Feron, V. Balakrishnan Linear Matrix Inequalities in System and Control Theory, SIAM, Phidelphia 1994.
[2] Dugard, L., E. I. Verriest Eds., Stability and Control of Time Delay Systems, Springer, London, New York 1998.
[3] Engelborghs K. T., T. Luzyanina, D. Roose, “Numerical Bifurcation Analysis of Delay Differantial Equations using DDE-BIFTOOL”, ACM Transactions on Math-ematical Software, vol.28, pp.1–21, 2002.
[4] Engelborghs K. T., T. Luzyanina, G. Samaey “DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations”, Report TW 330, Katholieke Univ. Leuven, 2001
[5] Engelborghs K. T., D. Roose, “On stabillity of LMS methods and charecteristic roots of delay differential equations”, SIAM J. Numer. Anal., vol.40, pp.629–650, 2002.
[6] Gu K., V. L. Kharitonov, J. Chen Stability of Time-Delay Systems, Birk¨auser, 2003.
[7] Hespanha J., A. Morse “Stability of switched systems with average dwell-time”, Proceedings of the 38. Conference on Decision and Control, Phoenix, USA, 1999 [8] Jarlebring E., “Computing the stability region in delay-space of a TDS using
polynomial eigenproblems”, Proc. of IFAC Workshop on Time-Delay Systems, L’Aquila, Italy, 2006.
[9] Kharitonov V. L. “Robust stability analysis of time delay systems: A Survey”, Annual Reviews in Control, vol.23, pp.185–196, 1999.
[10] Liberzon D., A. S. Morse “Basic problems in stability and design of switched systems,”, IEEE Control Systems Magazine, vol.19, pp.59–70, 1999.
[11] Liberzon D., R. Tempo “Gradient algorithms for finding common Lyapunov func-tions”, Prodings ofthe 42nd IEEE Conference on Decision and Control, Hawaii, USA, 2003.
[12] Nemirovski A., P. Gahinet, “The Projective Method for Solving Linear Matrix Inequalities”, Proc. Amer. Contr. Conf., pp. 840–844, Baltimore, Maryland, USA, 1994.
[13] Niculescu S. I. Delay Effects on Stability, A Robust Control Approach, Springer-Verlag, London 2001.
[14] Olga¸c N., R. Sipahi, “A practical method for analysing the stability of neutral type LTI-time delayed systems”, Automatica, vol.40, pp.847–853, 2004.
[15] Peleties P., R. DeCarlo “Asymptotic stability of m-switched systems using lya-punov functions”, Proceedings of the 31. Conference on Decision and Control, Tuc-son, USA, 1992.
[16] Richard J. P., “Time-delay systems: an overview of some recent advances and open problems”, Automatica, vol.39, pp.1667–1694, 2003.
[17] Roose D., K. T. Engelborghs, T. Luzyanina, W. Michiels, “Software for Sta-bility and Bifurcation Analysis of Delay Differentatial Equations to Stabilization ”, Lecture Notes in Computational Science and Engineering, vol.38, pp.167–181, 2004.
[18] Tong Q. Y., G. F. Yan, G. Z. Zhao, “Stability analysis of hybrid systems with time-varing delayed perturbation via single lyapunov function”, Proceedings of the Second International Conference on Machine Learning and Cybernetics, 2003.
[19] Vyhlidal T., P. Zitek, “Mapping the spectrum of a retarded time delay systems utilizing root distribution features ”, Proc. of IFAC Workshop on Time Delay Systems, LAquila, Italy, 2006.
[20] Vyhlidal T., P. Zitek, “Quasipolynomial mapping based rootfinder for analy-sis of time delay systems ”, Proc. of IFAC Workshop on Time Delay Systems, Rocquencourt, France, 2003.
[21] Yan P., H. ¨Ozbay, “Stability Analysis of Switched Time-Delay Systems,” Proc. of IFAC World Congress, Prague, Czech Rep., 2005.
[22] Ye H., N. Michel, L. Hou, “Stability theory for hybrid dynamical systems”, IEEE Transactions on Automatic Control, vol.43, pp.461–474, 1998.
