Luận án On the state estimation problem for some classes of dynamical systems and its application

n ] . (2.32) Accordingly, we present an effective algorithm to transform a general n-order time-delay system (n > 3) with single output into the observable form (2.15)-(2.16). Algorithm 2.1 Step 1: Obtain matrices Xn and Yn according to (2.22). Check if condition (2.23) is satisfied or not. If so, obtain χn where χn = YnX + n , where X + n denotes the Moore-Penrose inverse of Xn. Step 2: Substitute βk (k = 1, 2, . . . , n − 1) and α` (` = 1, 2, . . . , n − 2) into (2.26)-(2.27) and obtain Zn and Tn. Check if condition (2.24) is satisfied or not. If so, obtain ζn = TnZ + n , where Z + n denotes the Moore-Penrose inverse of Zn. Step 3: From (2.9)-(2.11), obtain matrices Mi and Ni (i = 1, 2, . . . , n) and hence the state transformation (2.8). Finally, obtain a transformed system according to (2.15)-(2.16). 27 Remark 2.1.1. Once a transformed system as described by (2.15)-(2.16), we can easily apply any Luenberger-typed state observers design method (see, for example, [80]) to design a state observer to estimate z(t) since the matrix pair (A¯, C¯) is now observable. After a satisfactory state observer z(t) has been designed, we can use the method of backward state transformations reported in [47]. 2.2 Application to the GI model 2.2.1 State transformation for the GI model In this section, we will apply the results obtained in the previous section to the GI model (2.5)-(2.7). By following the steps (Step 1-Step 3) of Algorithm 2.1, we obtain M1 = [ 1 0 0 0 ] , N1 = 01,4, M2 = 01,4, N2 = [ 0 a1 0 0 ] , M3 = 01,4, N3 = [ 0 −a1a2 0 a1a3 ] , M4 = [ a2a4a6 2 0 0 0 ] , N4 = [ 0 a1a 2 2 a1a3a5 −a1a3(a2 + a4) ] . Hence, we obtain the following state transformations z1(t) = ξ(t), z2(t) = a1x2(t− τ), z3(t) = −a1a2x2(t− τ) + a1a3x4(t− τ), z4(t) = a2a4a6 2 ξ(t) + a1a 2 2x2(t− τ) + a1a3a5x3(t− τ) −a1a3(a2 + a6)x4(t− τ) and a transformed system of the forms (2.15)-(2.16), where A¯ =  0 1 0 0 0 0 1 0 0 0 0 1 0 γ1 γ2 γ3  , B¯ =  0 0 0 a1a3a5  , C¯ = [ 1 0 0 0 ] , γ1 = −a2a4a6 2 , γ2 = −(a2a6 + a2a4 + a4a6), γ3 = −(a2 + a4 + a6), B¯1 = 04,1, Γ =  0 0 −a2a4a62 a2a4a6(a2+a4+a6) 2  , Γ3 =  1 0 0 0 0 0 0 0 0 0 0 0 a2a4a6 2 0 0 0  , Γ1 = 04,1, Γ4 =  0 0 0 0 0 a1 0 0 0 −a1a2 0 a1a3 0 a1a 2 2 a1a3a5 −a1a3(a2 + a6)  . 28 2.2.2 State observer design for the GI model Since the matrix pair (A¯, C¯) is observable, it is easy to design a state observer to estimate any linear function of the state vector z(t). Let h(t) =  h1(t)h2(t) h3(t)  = Fz(t) =  0 1 0 00 0 1 0 0 0 0 1  z(t) be a vector that is required to be estimated. To reconstruct the state function, h(t), we consider a functional observer of order 3 as follows: hˆ(t) = ω(t) + Ey¯(t), (2.33) ω˙(t) = Nω(t) + Jy¯(t) +Hu(t− τ) +Lµ¯(y¯(t), y¯(t− τ), y¯(t− τ − τg)), t > τmax, (2.34) where ω(t) ∈ R3, hˆ(t) ∈ R3 is the estimate of h(t), E, N , J , H and L are observer parameters to be determined. Let us define the following error vectors (t) and e(t) as (t) = ω(t)− Lz(t), (2.35) e(t) = hˆ(t)− Fz(t). (2.36) Based on [80, Theorem 3.1], hˆ(t) converges asymptotically to Fz(t) if the following conditions are satisfied N is Hurwitz, (2.37) NL+ JC¯ − LA¯ = 0, (2.38) H − LB¯ = 0, (2.39) F − EC¯ − L = 0. (2.40) Accordingly, for the given matrices A¯, C¯, and B¯ as above, we can easily solve (2.37)-(2.40) to obtain the following matrices: N = A¯22 + L1A¯12, L1 is chosen such that N is Hurwitz, E = −L1, J = −NL1, H = LB¯, where L = [ L1 L2 ] , A¯12 = [ 1 0 0 ] , A¯22 =  0 1 00 0 1 γ1 γ2 γ3  and L2 =  1 0 00 1 0 0 0 1 . Note that the matrix pair (A¯22, A¯12) is observable and thus L1 can be easily found to ensure that N is stable with any prescribed eigenvalues. Upon hˆ(t) has been obtained, then based the method of backward state transformations (Case 2) reported in [47], we obtain xˆ2(t− τ) = 1 a1 hˆ1(t), (2.41) xˆ3(t− τ) = 1 a1a3a5 [a2a6hˆ1(t) + (a2 + a6)hˆ2(t) +hˆ3(t)− a2a4a6 2 y¯(t)], (2.42) xˆ4(t− τ) = 1 a1a3 [a2hˆ1(t) + hˆ2(t)]. (2.43) 29 2.2.3 Simulation results In order to obtain simulation results, we consider the nonlinear time-delay GI model (2.1)-(2.3) with a set of parameters, the initial conditions and the input u(t) are as follows: a1 = 3.11 × 10−5, a2 = 1.211 × 10−2, a3 = 10.25×55 , a4 = a5 = a6 = 155 , a7 = 1, τ = 3min, g(y(t), y(t − τg)) = 3 187 1.573 0.25 ( x1(t−τg) 9 )3.205 1+ ( x1(t−τg) 9 )3.205 0 0 T , τg = 4min, x1(θ) = 10.66, x2(θ) = 49.29, x3(θ) = 0, x4(θ) = 0 for all θ ∈ [−4, 0], ω1(ζ) = 20e−0.07t, ω2(ζ) = 2e−t, ω3(ζ) = 5e−t for all t > 0, ζ ∈ [−8, 0] and u(t) = { sin t+ 50, 0 6 t 6 100, sin t+ 3, 100 < t 6 180. Let us now apply the reduced-order state observer (2.33)-(2.34) for this example. The eigenvalues of matrix N are chosen as, say, λ1 = −0.05, λ2 = −0.07, λ3 = −0.08, hence by using LMI Toolbox in Matlab, we obtain L =  −0.1515 1 0 0−0.0050 0 1 0 0.0001 0 0 1 , J =  0.0180−0.0008 0 , E =  0.15150.0050 −0.0001  and H =  00 0.4112× 10−5 . Figure 2.1 shows the responses of x2(t) and its delayed-estimation, i.e., xˆ2(t − 3), while, Figure 2.2 shows the responses of x2(t−3) and its estimation, i.e., xˆ2(t−3). It is clear from Figure 2.2 that the designed observer able to track the delayed version of the state vector, as expected. 30 0 20 40 60 80 100 120 140 160 180 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Time(min) Th ird −o rd er s ta te o bs er ve r x2(t) xˆ2(t − 3) Figure 2.1: Responses of x2(t) and xˆ2(t− 3) 31 0 20 40 60 80 100 120 140 160 180 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Time(min) Th ird −o rd er s ta te o bs er ve r xˆ2(t − 3) x2(t − 3) Figure 2.2: Responses of xˆ2(t− 3) and x2(t− 3) 32 Conclusion of Chapter 2 In this chapter, we have proposed a novel procedure for designing a state observer of a general nonlinear time-delay Glucose-Insulin model. The reported result is significant as the two-stage design process transforms a nonlinear time-delay model into a new observable form which allows a third-order delayed state observer to be easily designed. Simulation results have been given to illustrate the effectiveness of our results. The main advantage of the proposed method in this chapter is that the state transformations can be used to transform nonlinear time-delay systems into new coordinates where all the time-delay terms in the system description are associated with the output and input only. Therefore, in the new coordinate system, a Luenberger-type state observer can be readily designed. However, the advantage of this method is that there are some cases instead of an instantaneous state observer; we only obtain a delayed or a mixed instantaneous and delayed state observer of nonlinear time-delay systems. 33 Chapter 3 A new observer for interconnected time-delay systems and its applications to fault detection problem Faults of control systems may very often lead to a drastic reduction of system performance or a loss of stability, which even cause damages to the physical system [8], [56], [74]. Therefore, increasing attention has recently been devoted to the problem of fault detection of control systems (see, for example, [2], [10], [25], [58], [32], [19], [20]). To date and to the best of knowledge, the problem of designing functional observers to detect actuator faults of interconnected time-delay systems has not yet received adequate attention. In [78], the authors presented an approach to the design of fault detection scheme to detect actuator faults for a class of interconnected time-delay systems consisting of N subsystems with unknown input vector di(t) ∈ Rnid and actuator fault vector fi(t) ∈ Rnif entering from the i-th subsystem input defined as follows: x˙i(t) = Aiixi(t) + K∑ j=1(j 6=i) Aijxj(t− τji) +Adiixi(t− τii) +Biui(t) +Didi(t) + Eifi(t), t > 0, (3.1) xi(θ) = φi(θ), θ ∈ [−τmax, 0], (3.2) yi(t) = Cixi(t), (3.3) where xi(t) ∈ Rni , xj(t) ∈ Rnj , ui(t) ∈ Rmi , and yi(t) ∈ Rpi are the local state vector, the remote state vector, the control input vector, and the measured output vector of the i-th subsystem. Each φi(θ) (i = 1, 2, . . . ,K) is an initial function. Aii ∈ Rni×ni , Adii ∈ Rni×ni , Aij ∈ Rni×nj , Bi ∈ Rni×mi , Di ∈ Rni×nid , Ei ∈ Rni×nif and Ci ∈ Rpi×ni are real known system matrices. Without loss of generality, 34 it is assumed that rankCi = pi. In (3.2), τmax is defined as τmax = max 1≤i,j≤K,i6=j {τii, τji}. di(t) ∈ Rnid is the unknown disturbance and fi(t) ∈ Rnif is an unpredictable fault signal. To detect the faults fi(t) in the i-st subsystem, the authors in [78] proposed the following residual generator: ri(t) = Tiωi(t) + Fiyi(t), (3.4) ω˙i(t) = Niωi(t) +Ndiωi(t− τii) +Giyi(t) +Gdiyi(t− τii) +Hiui(t) + K∑ j=1(j 6=i) Gijyj(t− τji), t > 0, (3.5) ωi(θ) = ψi(θ), θ ∈ [−τmax, 0], (3.6) where ri(t) ∈ R is the local residual generator that is generated based on ωi(t) and the local output vector yi(t), Ti ∈ R1×qi and Ei ∈ R1×pi . Here ωi(t) ∈ Rqi , 1 6 qi 6 ni is the functional observer state vector, Ni ∈ Rqi×qi , Ndi ∈ Rqi×qi , Hi ∈ Rqi×mi , Gi ∈ Rqi×pi , Gdi ∈ Rqi×pi , Gij ∈ Rqi×pj are observer parameters that need to be determined such that ωi(t) asymptotically converges to a linear function of the local state vector Lixi(t) when there are no faults in the system, and Li ∈ Rqi×ni is a functional matrix to be determined for the purpose of fault detection. Nevertheless, there are instances where such a residual generator does not exist. To support this statement, let us consider the following motivated example: x˙1(t) = A11x1(t) +Ad11x1(t− τ11) +A12x2(t− τ21) +B1u1(t) +D1d1(t) + E1f1(t), t > 0, (3.7) x1(θ) = φ1(θ), θ ∈ [−τmax, 0], (3.8) y1(t) = C1x1(t), (3.9) x˙2(t) = A22x2(t) +Ad22x2(t− τ22) +A21x1(t− τ12) +B2u2(t) +D2d2(t) + E2f2(t), t > 0, (3.10) x2(θ) = φ2(θ), θ ∈ [−τmax, 0], (3.11) y2(t) = C2x2(t), (3.12) where τmax = max{τ11, τ22, τ12, τ21}, and x1(t) =  x11(t)x12(t) x13(t)  , x2(t) =  x21(t)x22(t) x23(t)  , A11 =  0.1 0.2 0.10.4 0.1 0 0.6 0.1 0.01  , A12 =  0.3 0.4 0−0.1 −0.2 0 0.2 0 0.1  , A21 =  0.1 0 00 0.2 0 0.6 0 0.2  , A22 =  0.1 0.4 −0.1−0.4 0.1 0 0.7 0.3 0.01  , Ad11 =  0.2 0.3 00.1 0.2 0 0.4 −0.5 0.1  , Ad22 =  0.2 0.1 00.1 −0.2 0 0.5 0.2 0.1  , B1 =  12 3  , B2 =  43 2  , 35 D1 =  12 3  , D2 =  −1−2 −3  , E1 =  00 1  , E2 =  00 2  , C1 = [ 1 0 0 0 1 0 ] , C2 = [ 1 0 0 0 1 0 ] . For 1-st subsystem, there exists a residual generator r1(t) that is constructed from the following reduced-order unknown input functional observer: r1(t) = T1ω1(t) + F1y1(t), (3.13) ω˙1(t) = N1ω1(t) +Nd1ω1(t− τ11) +G1y1(t) +Gd1y1(t− τ11) +H1u1(t) +G12y2(t− τ21), t > 0, (3.14) ω1(θ) = ψ1(θ), θ ∈ [−τmax, 0], (3.15) provided that all steps of Design Algorithm 1 in [78] are implemented. In the following, we will show that the method in [78] cannot apply to this example. Indeed, according to Step 1 of Design Algorithm 1 in [78], condition (13) is found to be satisfied since ν1 = rank(Ψ1) = 4 < p2+n1 = 5, where Ψ1 = [ C2 0 −A12 D1 ] . Accordingly, the nontrivial solution, [ G12 L1 ] 6= 0, to (11) exist. From (14), X1 is obtained, where X1 = [ 0.2748 0.3925 0.7851 −0.3925 0 ] . According to Step 2 of Design Algorithm 1 in [78], we search for the minimum order q1, where 1 6 q1 6 nullityΨ1 = 2. First, we assign q1 = 1 and matrix L1, extracted from X1 [ 0p2×n1 In1 ] , can be L1 = [ 0.7851 −0.3925 0 ] . Note that condition (23) does not hold since rank [ Ω¯1 Φ1 ] = 5 6= 4 = rank(Ω¯1). Similarly, it is easy to verify that the conditions in Step 2 for subsystem 2 are also not satisfied. Hence, the design method [78] does not work for this example. While, we will present in this chapter that, our method still works for this example. So far, to the best of our knowledge, no other observer design method is developed in the literature that enables the design of two stand-alone distributed functional observers for detecting faults for this example. Therefore, our main aim of this chapter is to propose a new observer design method. For this, we will develop the state transformation in [50] to obtain the general one which allows to solve easily the problem of fault detection for this kind of systems. For m,n ∈ N, n > 1 and arbitrary matrices M ∈ Rn×n and N ∈ R1×n, the following notations will be used throughout this chapter: • nζ = n+ 2p− 1, MT denotes the transpose of M, 0m,n denotes the m× n zero matrix, In denotes the identity matrix of dimension n× n; and • N = [ [N ]O [N ]U ] , where [N ]O ∈ R1×p and [N ]U ∈ R1×n are sub-matrices of N . 36 3.1 New state transformation Without loss of generality, we assume that matrix Ci takes the canonical form, where Ci =  C1i C2i ... Cpii  = [ Ipi 0pi,ni−pi ] . (3.16) Let us define the following new state transformation vector ζi(t) =  ζi1(t) ζi2(t) ... ζiniζ (t)  ∈ Rniζ , for i = 1, 2, . . . ,K, ζi(t) = [ T i T dii ] [ xi(t) xi(t− τii) ] + K∑ `=1(` 6=i) T `ix`(t− τ`i), t > τmax, (3.17) T i =  T i1 ... T iniζ  ∈ Rniζ×ni , T dii =  T dii1 ... T diiniζ  ∈ Rniζ×ni , T `i =  T `i1 ... T `iniζ  ∈ Rniζ×ni , (3.18) and the row matrices T iσ, T dii σ and T `i σ (σ = 1, 2, . . . , niζ ) are obtained according to the following T iσ =  Cσi, σ = 1, 2, . . . , pi, T iσ−piAii − ∑pi s=1 α s σT i s , σ = pi + 1, . . . , 2pi, T iσ−1Aii − ∑pi `=1 α ` σT i ` , σ = 2pi + 1, . . . , niζ , (3.19) T idσ =  01×ni , σ = 1, 2, . . . , pi, T iσ−piAdii − ∑pi s=1 β s σiT i s , σ = pi + 1, . . . , 2pi, T iσ−1Adii + T i d(σ−1)Aii − ∑pi `=1 β ` σiT i ` , σ = 2pi + 1, . . . , niζ , (3.20) T `iσ =  01×ni , σ = 1, 2, . . . , pi, T iσ−piAi` − ∑pi s=1 β s σ`T `i s , σ = pi + 1, . . . , 2pi, T iσ−1Ai` + T `i σ−1A`` − ∑pi s=1 β s σ`T `i s , σ = 2pi + 1, . . . , niζ , (3.21) where αsσ, β s σi, β s σ` (σ = pi + 1, . . . , niζ , s = 1, . . . , pi) are scalars to be determined later. We now denote the following matrices Hi1k = T i kAii, H i2 k = T i kAdii + T i dkAii, H i3 k = T i dkAdii , H i4 k = T i dkAi`, (3.22) Hi5k = T i kAi` + T `i k A``, H i6 k = T `i k Ad`` , H i7 k = T `i k A`i, H if k = T i dkEi, (3.23) Hidk = T i dkDi, H `f k = T `i k Ei, H `d k = T `i k Di. (3.24) 37 Remark 3.1.1. In [50], the authors proposed a state transformation method (state transformation (11) in [50]) to transform a class of single systems with a time delay (system (5)-(7) in [50]) into an observable canonical form. Note that if Aij = 0, ∀j = 1, . . . ,K, j 6= i and T `i = 0, ∀` = 1, . . . ,K, ` 6= i, then the system (3.1)-(3.3) and the state transform (3.17) in this chapter are reduced to the system (5)-(7) and the state transformation (11) in [50], respectively. Hence, the state transformation considered in [50] can be regarded as a special case of the one in this chapter. Remark 3.1.2. In contrast to the state transformations in [81], the state transformation (3.17) in this chapter utilizes multiple delayed states x`(t− τ`i), ` = 1, . . . ,K, ` 6= i and has a more general structure. This trade-off feature overcomes some drawbacks in [81] and enables distributed functional observers to be designed for a wider class of time-delay interconnected systems. The following theorem presents conditions for the existence of the state transformation vector (3.16) to transform each subsystem of system (3.1)-(3.3) into a novel canonical form. Theorem 3.1.1. For some scalars αji , β j i1, β j i2 (i = pi+1, pi+2, . . . , niζ , j = 1, 2, . . . , pi), if the following equations hold [ Hi1k ] U = 01×ni , k = pi + 1, . . . , 2pi − 1, k = niζ (3.25)[ Hi2k ] U = 01×ni , k = pi + 1, . . . , 2pi − 1, k = niζ , (3.26)[ Hi3k ] U = 01×ni , k = pi + 1, . . . , niζ , (3.27)[ Hi4k ] U = 01×n` , k = pi + 1, . . . , niζ , (3.28)[ Hi5k ] U = 01×n` , k = pi + 1, . . . , 2pi − 1, k = niζ , (3.29)[ Hi6k ] U = 01×n` , k = pi + 1, . . . , niζ , (3.30)[ Hi7k ] U = 01×ni , k = pi + 1, . . . , niζ , (3.31) Hifk = 01×nif , k = pi + 1, pi + 2, . . . , niζ , (3.32) Hidk = 01×niξ , k = pi + 1, pi + 2, . . . , niζ , (3.33) H`fk = 01×nif , k = pi + 1, pi + 2, . . . , niζ , (3.34) H`dk = 01×nid , k = pi + 1, pi + 2, . . . , niζ , (3.35) where Hi1k , H i2 k , H i3 k , H i4 k , H i5 k , H i6 k , H i7 k , H id k , H `f k and H `d k are as defined in (3.22)-(3.24), then the transformation vector (3.17) transforms each subsystem of system (3.1)-(3.3) into the following observable form ζ˙i(t) = A¯iζi(t) + B¯iui(t) + B¯ 1 i ui(t− τii) + K∑ `=1(` 6=i) B¯`iu`(t− τ`i) +D¯idi(t) + E¯ifi(t) + Γiyi(t) + Γ 1 i yi(t− τii) + Γ2i yi(t− 2τii) + K∑ `=1(` 6=i) Γ1`iy`(t− τ`i) + K∑ `=1( 6`=i) Γ2`iy`(t− τ`` − τ`i) 38 +K∑ `=1(` 6=i) Γ3`iy`(t− τii − τ`i) + K∑ `=1( 6`=i) Γ¯`siys(t− τs` − τ`i), t > 2τmax, (3.36) yi(t) = C¯iζi(t), (3.37) where A¯i =  0pi,pi Ipi 0pi,(niζ−2pi) 0pi−1,pi 0pi−1,pi 0pi−1,(niζ−2pi) 0niζ−2pi,pi 0niζ−2pi,pi Iniζ−2pi 01,pi 01,pi 01,niζ−2pi  , C¯i = [ Ipi 0pi,(niζ−pi) ] , and B¯i ∈ Rniζ×mi , B¯1i ∈ Rniζ×mi , B¯2i ∈ Rniζ×mi , B¯`i ∈ Rniζ×m2 , Γi ∈ Rniζ×pi , Γ1i ∈ Rniζ×pi , Γ2i ∈ Rniζ×pi , Γ¯`si ∈ Rniζ×pi , Γ1`i ∈ Rniζ×p` , Γ2`i ∈ Rniζ×pi , Γ3`i ∈ Rniζ×p` are defined as below B¯i =  T i1Bi T i2Bi ... T iniζ Bi  , B¯1i =  T id1Bi T id2Bi ... T iniζ Bi  , B¯`i =  T `i1 B` T `i1 B` ... T `iniζ B`  , D¯i =  T i1Di T i2Di ... T iniζ Di  , E¯i =  T i1Ei T i2Ei ... T iniζ Ei  , Γ2`i =  0pi,1 0pi,1 . . . 0pi,1[ Hi6pi+1 ] O1 [ Hi6pi+1 ] O2 . . . [ Hi6pi+1 ] Op` ... ... . . . ...[ Hi6niζ ] O1 [ Hi6niζ ] O2 . . . [ Hi6niζ ] Op`  , Γ3`i =  0pi,1 0pi,1 . . . 0pi,1[ Hi4pi+1 ] O1 [ Hi4pi+1 ] O2 . . . [ Hi4pi+1 ] Op` ... ... . . . ...[ Hi4niζ ] O1 [ Hi4niζ ] O2 . . . [ Hi4niζ ] Op`  , 39 Γi =  α1pi+1 α 2 pi+1 . . . α pi pi+1 ... ... . . . ... α12pi α 2 2pi . . . α pi 2pi[ Hi1pi+1 ] O1 [ Hi1pi+1 ] O2 . . . [ Hi1pi+1 ] Opi ... ... . . . ...[ Hi12pi−1 ] O1 [ Hi12pi−1 ] O2 . . . [ Hi12pi−1 ] Opi α12pi+1 α 2 2pi+1 . . . α pi 2pi+1 ... ... . . . ... α1niζ α2niζ . . . αpiniζ[ Hi1niζ ] O1 [ Hi1niζ ] O2 . . . [ Hi1niζ ] Opi  , Γ1i =  β1(pi+1)i β 2 (pi+1)i . . . βpi(pi+1)i ... ... . . . ... β1(2pi)i β 2 (2pi)i . . . βpi(2pi)i[ Hi2pi+1 ] O1 [ Hi2pi+1 ] O2 . . . [ Hi2pi+1 ] Opi ... ... . . . ...[ Hi22pi−1 ] O1 [ Hi22pi−1 ] O2 . . . [ Hi22pi−1 ] Opi β1(2pi+1)i β 2 (2pi+1)i . . . βpi(2pi+1)i ... ... . . . ... β1niζ i β2niζ i . . . βpiniζ i[ Hi2niζ ] O1 [ Hi2niζ ] O2 . . . [ Hi2niζ ] Opi  , Γ2i =  0pi,1 0pi,1 . . . 0pi,1[ Hi3pi+1 ] O1 [ Hi3pi+1 ] O2 . . . [ Hi3pi+1 ] Opi ... ... . . . ...[ Hi32pi ] O1 [ Hi32pi ] O2 . . . [ Hi32pi ] Opi[ Hi32p−1 ] O1 [ Hi32p−1 ] O2 . . . [ Hi32p−1 ] Opi[ Hi32pi+1 ] O1 [ Hi32pi+1 ] O2 . . . [ Hi32pi+1 ] Opi ... ... . . . ...[ Hi3niζ ] O1 [ Hi3niζ ] O2 . . . [ Hi3niζ ] Opi  , Γ¯`si =  0pi,1 0pi,1 . . . 0pi,1[ Hi7pi+1 ] O1 [ Hi7pi+1 ] O2 . . . [ Hi7pi+1 ] Opi ... ... . . . ...[ Hi7niζ ] O1 [ Hi7niζ ] O2 . . . [ Hi7niζ ] Opi  , 40 Γ1`i =  β1(pi+1)` β 2 (pi+1)` . . . βp`(pi+1)` ... ... . . . ... β1(2pi)` β 2 (2pi)` . . . βp`(2pi)`[ Hi5pi+1 ] O1 [ Hi5pi+1 ] O2 . . . [ Hi5pi+1 ] Op` ... ... . . . ...[ Hi52pi−1 ] O1 [ Hi52pi−1 ] O2 . . . [ Hi52pi−1 ] Op`[ Hi52pi ] O1 [ Hi52pi ] O2 . . . [ Hi52pi ] Op`[ Hi52pi+1 ] O1 [ Hi52pi+1 ] O2 . . . [ Hi52pi+1 ] Op` ... ... . . . ...[ Hi5niζ ] O1 [ Hi5niζ ] O2 . . . [ Hi5niζ ] Op`  . Proof. By taking the derivatives of (3.17), we obtain ζ˙i(t) = [ T i T dii ] [ x˙i(t) x˙i(t− τii) ] + K∑ `=1(` 6=i) T `ix˙`(t− τ`i), t > τmax. (3.38) For σ = 1, 2, . . . , pi, from (3.19)-(3.21), we have ζ˙iσ(t) = T i σx˙i(t) = T i σ(Aiixi(t) +Adiixi(t− τii) + K∑ `=1( 6`=i) Ai`x`(t− τ`i) +Biui(t) +Didi(t) + Eifi(t)) = T iσAiixi(t) + T i σAdiixi(t− τii) + K∑ `=1( 6`=i) T iσAi`x`(t− τ`i) +T iσBiui(t) + T i σDidi(t) + T i σEifi(t) = T i(σ+pi)xi(t) + T i d(σ+pi) xi(t− τii) + K∑ `=1(` 6=i) T `i(σ+pi)x`(t− τ`i) + pi∑ k=1 αkσ+piy i k(t) + pi∑ k=1 βk(σ+pi)iy i k(t− τii) + K∑ `=1(` 6=i) p∑` r=1 βr(σ+pi)`y ` r(t− τ`i) + T iσBiui(t) + T iσDidi(t) + T iσEifi(t) = ζiσ+pi(t) + pi∑ k=1 αkσ+piy i k(t) + pi∑ k=1 βk(σ+pi)iy i k(t− τii) + K∑ `=1( 6`=i) p∑` r=1 βr(σ+pi)`y ` r(t− τ`i) + T iσBiui(t) + T iσDidi(t) + T iσEifi(t). (3.39) For σ = pi + 1, . . . , 2pi − 1, from (3.19)-(3.21), we obtain ζ˙iσ(t) = T i σAiixi(t) + (T i σAdii + T i dσAii)xi(t− τii) + T idσAdiixi(t− 2τii) 41 +T iσ K∑ `=1( 6`=i) Ai`x`(t− τ`i) + T idσ K∑ `=1(` 6=i) Ai`x`(t− τii − τ`i) + K∑ `=1( 6`=i) T `iσ A``x`(t− τ`i) + K∑ `=1(` 6=i) T `iσ Ad``x`(t− τ`` − τ`i) + K∑ `=1( 6`=i) T `iσ K∑ s=1(s6=`) A`sxs(t− τs` − τ`i) + T iσBiui(t) + T iσDidi(t) +T iσEifi(t) + T i dσBiui(t− τii) + K∑ `=1(` 6=i) T `iσ B`u`(t− τ`i) = [ Hi1σ ] O yi(t) + [ Hi2σ ] O yi(t− τii) + [ K∑ `=1( 6`=i) K∑ s=1(s6=`) Hi7σ ] O ys(t− τs` − τ`i) + [ K∑ `=1( 6`=i) Hi5σ ] O y`(t− τ`i) + [ Hi3σ ] O yi(t− 2τii) + [ K∑ `=1( 6`=i) Hi4σ ] O y`(t− τii − τ`i) + [ K∑ `=1(` 6=i) Hi6σ ] O y`(t− τ`` − τ`i) + T iσBiui(t) + T iσDidi(t) + T iσEifi(t) +T idσBiui(t− τii) + K∑ `=1(` 6=i) T `iσ B`u`(t− τ`i). (3.40) For σ = 2pi, from (3.19)-(3.21), we obtain ζ˙iσ(t) = T i σ(Aiixi(t) +Adiixi(t− τii) + K∑ `=1( 6`=i) Ai`x`(t− τ`i) +Biui(t) +Didi(t) + Eifi(t)) +T idσ(Aiixi(t− τii) +Adiixi(t− 2τii) + K∑ `=1( 6`=i) Ai`x`(t− τii − τ`i) +Biui(t− τii) +Didi(t− τii) + Eifi(t− τii)) + K∑ `=1(` 6=i) T `iσ (A``x`(t− τ`i) +Ad``x`(t− τ`` − τ`i) + K∑ s=1(s6=`) A`sxs(t− τs` − τ`i) +B`u`(t− τ`i) +D`d`(t− τ`i)) + E`f`(t− τ`i)) = T iσAiixi(t) + (T i σAdii + T i dσAii)xi(t− τii) +T idσAdiixi(t− 2τii) + T iσ K∑ s=1(s6=`) Ai`x`(t− τ`i) + T idσ K∑ s=1(s6=`) Ai`x`(t− τii − τ`i) + K∑ s=1(s6=`) T `iσ A``x`(t− τ`i) + K∑ s=1(s6=`) T `iσ Ad``x`(t− τ`` − τ`i) 42 +K∑ s=1(s6=`) T `iσ K∑ s=1(s 6=`) A`sxs(t− τs` − τ`i) + T iσBiui(t) +T iσDidi(t) + T i σEifi(t) + T i dσBiui(t− τii) + K∑ s=1(s6=`) T `iσ B`u`(t− τ`i) = ζi2pi+1(t) +  α12pi+1 ... αpi2pi+1  T yi(t) +  β1(2pi+1)i ... βpi(2pi+1)i  T yi(t− τii) + K∑ `=1(` 6=i)  β1(2pi+1)` ... βp`(2pi+1)`  T y`(t− τ`i) + [ Hi3σ ] O yi(t− 2τii) + [ K∑ `=1(` 6=i) K∑ s=1(s6=`) Hi7σ ] O ys(t− τs` − τ`i) + [ K∑ `=1( 6`=i) Hi4σ ] O y`(t− τii − τ`i) + [ K∑ `=1(` 6=i) Hi6σ ] O y`(t− τ`` − τ`i) + T iσBiui(t) + T iσDidi(t) + T iσEifi(t) +T idσBiui(t− τii) + K∑ `=1( 6`=i) T `iσ B`u`(t− τ`i). (3.41) For σ = 2pi + 1, . . . , niζ − 1, from (3.38), we obtain ζ˙iσ(t) = T i σ(Aiixi(t) +Adiixi(t− τii) + K∑ `=1(` 6=i) Ai`x`(t− τ`i) +Biui(t) +Didi(t) + Eifi(t)) + T i dσ(Aiixi(t− τii) +Adiixi(t− 2τii) + K∑ `=1( 6`=i) Ai`x`(t− τii − τ`i) +Biui(t− τii) +Didi(t− τii) + Eifi(t− τii)) + K∑ `=1( 6`=i) T `iσ (A``x`(t− τ`i) +Ad``x`(t− τ`` − τ`i) + K∑ s=1(s 6=`) A`sxs(t− τs` − τ`i) +B`u`(t− τ`i) +D`d`(t− τ`i) + E`f`(t− τ`i)) = T iσAiixi(t) + (T i σAdii + T i dσAii)xi(t− τii) + K∑ `=1( 6`=i) (T iσAi` + T `i σ A``)x`(t− τ`i) + [ K∑ `=1( 6`=i) Hi6σ ] O y`(t− τ`` − τ`i) + [ K∑ `=1( 6`=i) K∑ s=1(s6=`) Hi7σ ] O ys(t− τs` − τ`i) + [ K∑ `=1( 6`=i) Hi4σ ] O y`(t− τii − τ`i) + T iσBiui(t) 43 +T iσDidi(t) + T i σEifi(t) + T i dσBiui(t− τii) + K∑ `=1( 6`=i) T `iσ B`u`(t− τ`i). (3.42) From (3.19)-(3.21), we have for σ = 2pi, 2pi + 1, . . . , niζ − 1 and ` = 1, . . . , N, ` 6= i, T iσAiixi(t) = T i (σ+1)xi(t) +  α1σ+1 ... αpiσ+1  T yi(t), (3.43) (T iσAdii + T i dσAii)xi(t− τii) = T id(σ+1)xi(t− τii) +  β1(σ+1)i ... βpi(σ+1)i  T yi(t− τii), (3.44) (T iσAi` + T `i σ A``)x`(t− τ`i) = T `i(σ+1)x`(t− τ`i) +  β1(σ+1)` ... βp`(σ+1)`  T