Relative stability of linear systems

RELATIVE STABILITY OF LINEAR SYSTEMS Jorge S´a Esteves ∗,1Paula Rocha ∗,2 ∗University of Aveiro, Portugal Abstract: In this paper we analyze the evolution of the relative error in the state of linear systems subject to perturbations. For this purpose we introduce the property of relative stability. The relevance of this property is illustrated by means of a case study, namely the computation of the Erlang-B function, which is well-known in the context of Queueing Theory. Keywords: Linear time-varying systems, stability, relative error, Erlang-B function 1. INTRODUCTION The stability properties of a dynamical system reflect, in a certain sense, the sensitivity of the system to small perturbations in the initial conditions. For a system with trajectories x(., x(0)), where x(0) denotes the initial condition, the sensitivity to an error e(0) is measured in terms of the absolute error kx(., x(0)+e(0))−x(., x(0))k. This issue has been widely studied and the obtained results both for linear and nonlinear systems are nowadays classical in the theory of dynamical and control systems, (Kailath, 1980), (Hirsch et al., 2004). However, for numerical purposes it makes often more sense consider the relative error, rather than the absolute one. In this paper we consider linear time-varying state space systems and analyze the evolution of the relative error when the initial condition is subject 1The research is supported by the Center for Research on Optimization and Control (CEOC), University of Aveiro, Portugal, from the “Funda¸c˜ao para a Ciˆencia e a Tecnologia” (FCT), co-financed by the European Community Fund FEDER/POCI 2010. 2The research is partially supported by the Unidade de Investiga¸c˜ao Matem´atica e Aplica¸c˜oes (UIMA), University of Aveiro, Portugal, through the Programa Operacional “Ciˆencia e Tecnologia e Inova¸c˜ao”(POCTI) of the Funda¸c˜ao para a Ciˆencia e a Tecnologia (FCT), cofinanced by the European Union fund FEDER. to perturbations. Unlike what happens for the classical stability property (to which we shall often refer as absolute stability) we allow the input to be nonzero. If, for a given input u, the relative error remains bounded, we say that the system is relatively stable with respect to the input u. In case, additionally, the relative error tends to zero in time we say that the system is relatively asymptotically stable. Although our purpose is not to present a detailed analysis of this property, we obtain some preliminary results concerning scalar systems and higher dimensional systems with sufficient excitation (a property of the input that we shall clarify later on). It turns out that scalar linear systems are always relatively stable, but not relatively asymptotically stable, with respect to the zero input. On the other hand, higher dimensional systems that are absolutely asymptotically stable are also relatively asymptotically stable with respect to sufficiently exciting inputs. An interesting situation even in case a system is only relatively asymptotically stable, is when the relative error is strongly attenuated in time with respect to the initial one. This situation is illustrated with the case study of the computation of the Erlang-B function in the context of Teletraffic Theory. 2. RELATIVE STABILITY Let Σ be a system described by the time-varying state space equations ½x(k+ 1) = A(k)x(k) + u(k) y(k) = x(k)(1) where, for k= 0,1, . . .,A(k)∈Rn×n, and u(k), x(k), y(k)∈Rn. As is well-known, the state trajectory generated by an input u(.) and an initial condition x(0) can be written as x(k, x(0), u) = xfree(k, x(0)) + xforced(k, u), where the free response xfree(., x(0)) and the forced response xforced(., u) are respectively given by xfree(k, x(0)) = A(k−1) · · · A(0)x(0) (2) and xforced(k, u) = k−2 X j=0 A(k−1) · · · A(j+ 1)B(j)u(j) +B(k−1)u(k−1). (3) Suppose now that the initial condition x(0) is affected by an error e(0). Then, the corresponding solution x(., ˆx(0), u), with ˆx(0) = x(0) + e(0), becomes x(k, ˆx(0), u) = xfree(k, ˆx(0)) + xforced(k, u) and hence the error e(.). =x(., ˆx(0), u)−x(., x(0), u) is given by e(k) = xfree(k, ˆx(0)) −xfree(k, x(0)) =xfree(k, ˆx(0) −x(0)) =A(k−1) · · · A(0)e(0). Under the assumption that the state vector x(k) never becomes null, we may define the relative error ²(.) by means of ²(k) = e(k) kx(k)k. Thus, ²(k) = A(k−1) · · · A(0)e(0) kA(k−1) · · · A(0)x(0) + xforced(k, u)k. (4) Definition 1. •The system Σ described by (1) is said to be relatively stable with respect to the input u, or to have a relatively stable response to the input u, if there exists a positive real number Lsuch that, for all x(0), e(0) ∈Rn, k²(k)k< Lk²(0)k. •The system Σ described by (1) is said to be relatively asymptotically stable with respect to the input u, or to have a relatively asymptotically stable response to the input u, if, in addition to the previous property of relative stability it satisfies the condition that, for all x(0), e(0) ∈Rn, limk→∞ ²(k) = 0. 2 Note that when the input is zero, i.e., u(k) = 0, k = 0,1, . . ., then ²(k) = A(k−1) · · · A(0)e(0) kA(k−1) · · · A(0)x(0)k. In the particular case of scalar systems, i.e, for n= 1, this implies that k²(k)k=ke(0)k kx(0)k=k²(0)k, showing that the relative error remains constant. Hence, the free response of a scalar system is relatively stable (but not relatively asymptotically stable) no matter what its the (absolute) stability properties are. However the situation is different for the vector case, as shown in the following example. Example 1. Let Σ be a system described by (1), with A(k) = ·λ10 0λ2¸, for all k. For the initial condition x(0) = [0 1]Tand the initial error e(0) = [1 0]T, we have k²(k)k=¯¯¯¯ λ1 λ2¯¯¯¯ k . Clearly, if |λ1|>|λ2|, then k²(k)k → ∞. On the other hand, for the initial condition x(0) = [1 0]Tand the initial error e(0) = [0 1]T, we have k²(k)k=¯¯¯¯ λ2 λ1¯¯¯¯ k and therefore if |λ2|>|λ1|, then we also have that k²(k)k → ∞. So, the system will be relatively unstable if |λ1| 6= |λ2|. Moreover, it is not difficult to show that if |λ1|=|λ2|the system is relatively stable. 2 A detailed analysis of the relative stability of the free response of a system will be carried out elsewhere. We now turn to the study of the relative stability of the response to a non-zero input, which is more directly connected to the case study presented in Section 3. For this purpose, rewrite (4) as ²(k) = A(k−1) · · · A(0)²(0) kA(k−1) · · · A(0) x(0) kx(0)k+xforced(k,u) kx(0)kk. If the system is (absolutely) asymptotically stable, then lim k→∞ A(k−1) ···A(0)²(0) = 0 ∈Rn and lim k→∞ A(k−1) · · · A(0) x(0) kx(0)k= 0 ∈Rn. Thus, ²(k) will tend to zero if xforced(k, u) does not go to zero as kgoes to infinity. The fact that the forced response to the input udoes not go to zero can be interpreted, in a certain sense, as the ability that this input has to excite the system. In this way, speaking in very loose terms, we may say that an (absolutely) asymptotically stable system has relatively asymptotically stable responses to inputs that are sufficiently exciting. Let us consider once more the scalar case (n= 1). Now (4) takes the form |²(k)|=|a(k−1) · · · a(0)e(0)| |a(k−1) · · · a(0)x(0) + xforced(k, u)| =|²(0)| |1 + xforced(k,u) a(k−1)···a(0)x(0) | by (3) z}|{ =|²(0)| |1 + Pk−1 j=0 u(j) a(j)···a(0)x(0) | (5) where the low case a(j) are used instead of the capitals A(j). A particularly interesting situation for our case study is when u(j) = 1, j= 0,1, . . ., i.e., when the input is a step. In this case it follows from (5) that |²(k)|=1 |1 + S(k)||²(0)|,(6) with S(k). =Pk−1 j=0 1 a(j)···a(0)x(0) . Once more, it is easily seen that if the system is (absolutely) asymptotically stable, then xfree(j, x(0)) →0 as j→ ∞ and consequently S(k)→ ∞. This implies that |²(k)| → 0, meaning that the system is relatively asymptotically stable with respect to the step input. If, on the contrary, xfree(j, x(0)) has a quick growth, so that S(k) converges to a value S∈R, then |²(k)| → 1 |1 + S||²(0)|. If S(k) takes on large values, the relative error will be attenuated with respect to its initial value after ktime units. As we shall see in the next section, this fact can be used in order to increase the computational efficiency in certain procedures. 3. CASE STUDY — COMPUTATION OF THE ERLANG-B FUNCTION The Erlang B and C formulas are true probability classics. Indeed, much of the theory was developed by A. K. Erlang (Erlang, 1909; Erlang, 1917) and his colleagues prior to 1925 (Brockmeyer et al., 1948). The subject has been extensively studied and applied by telecommunications engineers and mathematicians ever since. A nice introductory account, including some of the telecommunications subtleties, is provided by (Cooper, 1981). The Erlang B (or loss) formula gives the (steadystate) blocking probability in the Erlang loss model, i.e., in the M/M/s/0 model. This model has shomogeneous servers working in parallel and no extra waiting space. Customers arriving when all sservers are busy are blocked (lost) without affecting future arrivals; e.g., there are no customer retrials. This model has a Poisson arrival process and IID (independent and identically distributed) service times (which are also independent of the arrival process) with an exponential distribution having finite mean. (The two M’s in M/M/s/0 are for Markov, referring to the “lack-of-memory” property of the exponential distribution. Both the interarrival times and the service times have exponential distributions.) Following convention, let the arrival rate be denoted by λand let the mean service time be denoted by 1/µ. Thus, the (individual) service rate is µ. Since at most scustomers can be in the system at any time, the stochastic process representing the number of busy servers as a function of time has a proper steady-state distribution for all (positive) values of the parameters λand µ. The Erlang loss model has an insensitivity property implying that the blocking probability is independent of the service-time distribution beyond its mean. The steady-state distribution of the number of busy servers also does not depend on the units we use to measure time. Thus the blocking probability depends on the arrival rate λand the service rate µonly through their ratio, the offered load, a. =λ/µ. (7) As indicated above, the Erlang B formula gives the steady-state blocking probability of a typical arrival, and is given by (see for example (Cooper, 1981, pp. 5 and 79)): Ba(s). =as/s! Ps j=0 aj/j!, s ∈N0, a ∈R+.(8) The numerical studies regarding this formula are usually based on its analytical continuation, ascribed to R. Fortet (Sysky, 1960, pag.602) (Jagerman, 1974): Ba(s). =µaZ+∞ 0 e−az (1 + z)sdz¶−1 ,(9) which is valid for traffic offered a∈R+and s∈R+ 0servers. An important related quantity (because it tends to be easier to analyze) is the reciprocal, or the inverse probability of blocking: Ia(s). =1 Ba(s).(10) A major known result, due to D. L. Jagerman (Jagerman, 1984), is the following recursion obtained by partial integration of (9): Ia(s+ 1) = s+ 1 aIa(s)+1, s ∈R+ 0,(11) Since Ba(0) = Ia(0) = 1 for all a∈R+,Ia(s) may be calculated by recursion (11) for any positive integer s. The analysis of the propagation of the error associated with the initial value, throughout the proposed iterative procedure, is a key question for evaluating the effectiveness of the method and will be dealt with in this section. Some related numerical methods (including for calculating derivatives of the Erlang B function) have been presented in previous works (Esteves et al., 1995; Esteves et al., 1997). It will be assumed that even Ia(s) which is an exact value for s= 0 (Ia(0) = 1) may have an associated error. This will enable to obtain perfectly general results, namely for situations where the initial values are not established for s= 0. This happens, for example, if we want to analyze the effect of perturbations at a generic point s∈R+ 0. Moreover, it will be assumed that the rounding errors inherent to the recursive calculation (11) are disregardable. As a result of the error associated with the initial value the value calculated for Ia(s) will be denoted by e Ia(s) such that e(s) is the corresponding error: e Ia(s) = Ia(s) + e(s).(12) The relative error is ²(s): e Ia(s) = [ 1 + ²(s) ] Ia(s).(13) After calculating an approximation e Ia(s) with relative error ²(s) it is possible to compute an approximation f Ba(s) of Ba(s): f Ba(s) = 1 e Ia(s) f Ba(s) = 1 [1 + ²(s)]Ia(s) f Ba(s) = 1 1 + ²(s)Ba(s) Using the binomial expansion (1 + v)−1= 1 −v+ v2−v3+v4+··· (convergent if |v|<1), it may be written for small ²(s) (good approximations): f Ba(s) = £1−²(s) + [²(s)]2−[²(s)]3+· · · ¤Ba(s) f Ba(s)≈[1 −²(s)] Ba(s).(14) From (14), it follows that the absolute value of relative error of f Ba(s) is also nearly equal to |²(s)|. Thence, it may be said that the calculation of approximations of the values of the Erlang B function is made with an accuracy of the same order of the calculations of its reverse Ia(s). This fact justifies that in the following our only objective is to analyze the propagation of the relative error ²(.) of the inverse probability Ia(.) along the recursive procedure (11). We shall concentrate on the computation of Ia(s) for integer values of sand replace the variable sby the variable kused in the previous section to denote the discrete time instants. With this change of notation we obtain: Ia(k+ 1) = k+ 1 aIa(k)+1, k = 0,1,2, . . . , with initial condition Ia(0) = 1. Note that this is precisely a state space system of the form (1), with state Ia(instead of x), input u(.)≡1, and scalar A(k) = k+1 a. Thus, applying (6), and taking into account that 1 + S(k) is in our case always positive, we can write, after some simple computations: ²(k) = 1 1 + S(k)²(0), with S(k) = k−1 X j=0 aj+1 (j+ 1)!. This leads to the following result. Theorem 1. If the recursive relation (11) is used to calculate approximations to Ia(s), s= 1,2,3, . . ., then: lim s→∞ ²(s) = e−a²(0) . Thus, although the system is not relatively asymptotically stable, it is relatively stable and the relative error is attenuated in time. In fact, for aand ssufficiently high, Ba(s) may be calculated with great accuracy even if the initial value of the iteration is a very rough estimate. The following example will illustrate this in a perhaps surprising manner. If recursion (11) is used to calculate B100(100) by starting, as usual, with the initial value I100(0) = 1, one obtains (using double precision arithmetic): [I100(100)]−1= 0.075 700 452 710 860 97 (15) Considering now absurd initial values, such as 1027 or −1027, the result obtained after 100 recursive steps is exactly the same as (15). In the sequel we study how the relative error is propagated throughout a finite number of successive steps of the recursion, and show that the attenuation of the relative error does not only happen from s= 0 to s=k, but also from an arbitrary step sto step s+l. This study is the starting point for establishing a method that allows a fast computation with good accuracy. In particular, good bounding techniques are needed to obtain an efficient process of calculating adequate bounds for the relative error. The proposed approach is based on the following two preparatory lemmas. Lemma 1. If |²(s)| ≤ ξand l≥1, then |²(s+l)| ≤ ξ l Y j=1 s+j a. Proof: Ba(s) is a decreasing function in s, thus Ia(s+1) > Ia(s) for all x∈R+ 0.The result follows from the fact that ²(s+ 1) = (s+1)Ia(s) aIa(s+1) ²(s). 2 Since it is trivial that |²(s+l)|< ξ, Lemma 1 will be used only in cases such that s+l≤a. In other words, we are only able to quantify the decay of the relative error in this situation. The following inequality allows an important simplification in the expressions for the bound of the relative error. Note that the given expression is specially simple, since it is a quadratic function of l. Lemma 2. If a∈R+,s∈R+ 0and l≥n, then: ln   l Y j=n s+j a ≤ − 2 (a−s)−l−n 2a(l−n+1). Proof: Applying the arithmetic-geometric mean inequality, we obtain: P= l Y j=n s+j a≤ 1 a(l−n+ 1) l X j=n (s+j)  l−n+1 . Since Pl j=n(s+j) is the sum of (l−n+ 1) terms of an arithmetic progression, it follows that: ln P≤(l−n+ 1) ln µ2s+n+l 2a¶. The result follows applying the known inequality ln z≤z−1, z > 0.2 For easily estimating the number of correct digits of the approximation ˜ Ia(s+l) it is important to establish an efficient process for calculating a bound for ln |²(s+l)|. The next lemma suggests a simple analytic expression, which gives an adequate bound for the intended purpose. Lemma 3. If |²(s)| ≤ ξand l≥1, s≥0, then ln |²(s+l)| ≤ − 2(a−s)−l−1 2al+ ln ξ . Proof: Applying Lemma 1, and Lemma 2 the result is easily proved. 2 Note now that in order to compute Ia(s∗), with s∗=s∗ 0+land s∗ 0having the same fractional part as s∗, one may proceed as follows. Take the (crude) approximation Ia(s∗ 0)≈0. Then ²(s∗ 0) = 0−Ia(s∗ 0) Ia(s∗ 0)=−1. In other words, zero is an approximation of any positive quantity with exactly 100% of error, that is |²(s∗ 0)|=ξ= 1. Recall that, from Lemma 3, after lrecursive steps we have calculated an approximation of Ia(s+l) with precision |²(s+l)|such that: ln |²(s+l)|≤−2(a−s)−l−1 2al . (16) Thus if a >> s this approximation is very accurate. Additionally, we can guarantee that any perturbation introduced in the values of the recursive calculations decreases rapidly. This procedure obviously decreases the computational burden, since it reduces the number of iterations, maintaining a good accuracy. 4. CONCLUSION We have introduced the concept of relative stability to study the evolution of the relative error of the state trajectory of a linear system when the initial condition is subject to perturbations. This issue has been analyzed for the recursive computation of the Erlang-B function, in the context of Teletraffic Theory. It turns out that in this case the relative error is strongly attenuated in a small number of steps, enabling the use of an efficient computational method. REFERENCES Brockmeyer, E., H. L. Halstrom and A. Jensen (1948). The Life and Works of A. K. Erlang. Danish Academy of Technical Sciences. Copenhagen. Cooper, R. (1981). Introduction to Queueing Theory. North Holland. Erlang, A. K. (1909). The theory of probabilities and telephone conversations. (Danish) Nyt Tidsskrift for Matematik [English Translation in the P. O. Elect. Engrs. Journal 10, 1917– 1918, 197–203] 20, 33–39. Erlang, A. K. (1917). Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges. (Danish) Elektroteknikeren [English Translation in the P. O. Elect. Engrs. Journal 10, 1917–1918, 189–197] 13, 5–28. Esteves, Jorge S´a, J. Craveirinha and D. Cardoso (1995). Computing Erlang-B function derivatives in the number of servers — a generalized recursion. ORSA Communications in Statistics, Stochastic Models 11(2), 311–331. Esteves, Jorge S´a, J. Craveirinha and D. Cardoso (1997). A reduced recursion for computing Erlang-B function derivatives. In: Proceedings of the 15th International Teletraffic Congress, Washington D.C., USA (V. Ramaswami and P. E. Wirth (AT&T), Eds.). number 2 Teletraffic Contributions for the Information Age In: Teletraffic Science and Engineering. Elsevier Science B. V.. Boston. pp. 1315–1326. Hirsch, M. W., S. Smale and R. L. Devaney (2004). Differential Equations, Dynamical Systems & An Introduction to Chaos. Elsevier Science, Academic Press. Jagerman, D. L. (1974). Some properties of the Erlang loss function. The Bell System Technical Journal 53(3), 525–551. Jagerman, D. L. (1984). Methods in traffic calculations. The Bell System Technical Journal 63(7), 1283–1310. Kailath, T. (1980). Linear Systems. Prentice-Hall. Sysky, R. (1960). Congestion Theory in Telephone Systems. Oliver and Boyd.