The aim of this brief is to present a finite-step algorithm for the positive realization of a rational

transfer function H(z). In comparision with previously described algorithms we emphasize that we do

not make an a priori assumption on (but, instead, include a finite step procedure for checking) the non-

negativity of the impulse response sequence of H(z). For primitive transfer functions a new method for

reducing the pole order of the dominant pole is also proposed

Matolcsi, Máté

Nagy, Béla

English

B. Nagy

M. Matolcsi

Crossref

Algorithm for positive realization of transfer functions

Repository of the Academy's Library

1Algorithm for positive realization of transferfunctionsBe´la Nagy, Ma´te´ MatolcsiBe´la Nagy is Professor at the Mathematics Department of Technical University, Budapest. E-mail:bnagy@math.bme.huMa´te´ Matolcsi is a research assistant in Applied Mathematics at the Alfred Renyi Institute, Budapest. E-mail:matomate@renyi.hu27. November, 2002 DRAFT2AbstractThe aim of this brief is to present a finite-step algorithm for the positive realization of a rationaltransfer function H(z). In comparision with previously described algorithms we emphasize that we donot make an a priori assumption on (but, instead, include a finite step procedure for checking) the non-negativity of the impulse response sequence of H(z). For primitive transfer functions a new method forreducing the pole order of the dominant pole is also proposed.KeywordsPositive linear systems, impulse response sequence, positive realization.I. IntroductionThe positive realization problem of a given rational transfer function H(z) of a discretetime-invariant linear system is to find a triple A ∈ RN×N+ , b ∈ RN+ , c ∈ RN+ (with non-negative entries) such that H(z) = cT (zI −A)−1b holds. The nonnegativity restriction onthe entries of A, b, c reflect physical constraints in applications. Such positive systems ap-pear for example in modelling bio-systems, chemical reaction systems, and socio-economicsystems, as described in detail in the monograph by Farina and Rinaldi [5]. A recentapplication of positive systems in the construction of CRN’s (Charge Routing Networks)was given by Benvenuti, Farina and Anderson [2].In [9] Ohta, Maeda and Kodama reduced the problem of positive realizability of H(z) tofinding an appropriate convex polyhedral cone in the room sandwiched by the reachabilityand observability cones in the state space of an arbitrary minimal realization of H(z).However, the problem of constructing such a polyhedral cone turned out to be highlynontrivial (a characterization of all such cones is still lacking). In [1] Anderson, Deistler,Farina and Benvenuti proved that such a cone is always possible to construct if H(z) is aprimitive transfer function with nonnegative impulse response. Finally, the case of non-primitive transfer functions H(z) with nonegative impulse response was settled by Farinain [4] by the method of downsampling the impulse response of H(z) (see also Kitano andMaeda [8], and Fo¨rster and Nagy [6]). A common feature of the results mentioned above isthat the impulse response of H(z) is assumed to be nonnegative. However, it has not beenshown so far how one can check, in finite steps, the nonnegativity of the (clearly infinite)impulse response sequence corresponding to H(z) (this open problem was raised in [1]).27. November, 2002 DRAFT3In the course of this brief we aim to supplement the theory of positive realizability bytackling this problem in Section 2. In Section 3 we propose a new method of constructinga positive realization of a primitive transfer function with multiple dominant pole. InSection 4 we illustrate our results by an example.II. The nonnegativity of the impulse response sequenceLetH(z) =p1zn−1 + ...+ pnzn + q1zn−1 + ... + qn=r∑j=1nj∑i=1c(i)j(z − λj)i (1)be a strictly proper rational transfer function, where λ1 denotes the nonnegative pole ofH(z) with greatest modulus. Note that the coefficients pj and qj are assumed to be real,but the poles λj (j 6= 1) can be complex. We will use the notation hk (k = 1, 2, 3...) forthe impulse response sequence of H(z), i.e. H(z) =∑∞k=1 hkz−k. A minimal realizationof H(z) will be denoted by (g, F, h). We will describe the structure of our algorithm inseveral steps. The first two steps are standard in the theory of positive realizations (cf.[1] and [4]), but we include them for completeness and convenience.Step 1. It is well-known (cf. [1]) that a necessary condition for the existence of a positiverealization of H(z) is that λ1 be a dominant pole of H(z), i.e. the modulus of no poleexceed λ1. Therefore, if there is no nonnegative pole of H(z), or λ1 is not dominant, wecan conclude that H(z) does not have positive realizations. Assume that λ1 is dominant.If λ1 = 0, then the realization problem is trivial, and if λ1 > 0, then it is also well-knownthat we may (and will) assume without loss of generality that λ1 = 1 (cf. [1]).Step 2. If H(z) is not primitive (i.e. H(z) has dominant poles other than λ1 = 1),then a necessary condition for the existence of a positive realization is that the dominantpoles of H(z) be cyclic (see [1]), i.e. there exist p ∈ N such that all the dominant poles ofH(z) satisfy the equation zp = 1. If the dominant poles are not cyclic, then we concludethat there is no positive realization of H(z). If the dominant poles are cyclic with indexp (the smallest of the p’s above), then the necessary and sufficient condition for positiverealizability of H(z) is that all the “downsampled” transfer functions H(j)(z) := g(zI −F p)−1F jh (for j = 0, 1, ..., p − 1) be positively realizable (see [4] and [8]). Notice thatH(z) =∑p−1j=0 zp−1−jH(j)(zp).27. November, 2002 DRAFT4If some of the functions H(j)(z) are not primitive, then we apply the downsampling stepagain to these functions (provided that they are cyclic), and, following the algorithm of[4] (cf. also [8]), we arrive (after a finite number of steps) at a decomposition of the formH(z) =∑szβsH(s)(zαs) (2)where 0 ≤ βs < αs, and all the functions H(s)(z) are either primitive or not cyclic. If forany s the function H(s)(z) is not cyclic, then we conclude that H(z) does not have positiverealizations (cf. [4]). Assume therefore that all H(s)(z) are primitive. In this case, H(z) ispositive realizable if and only if the impulse response of H(z) is nonnegative (cf. [4]).Step 3. In this step we give an upper estimate on the finite number of the terms of theimpulse response sequence hk of H(z) whose nonnegativity we need to check in order toconclude that the whole impulse response sequence is nonnegative.Instead of checking the impulse response of H(z) directly, we take the decompositionH(z) =∑s zβsH(s)(zαs) of Step 2, and check the impulse response of each H(s)(z). If theimpulse response of each H(s)(z) is non-negative, then clearly so is the impulse responseof H(z). The advantage of this method is that all the functions H(s)(z) are primitive.For the sake of simplicity we will still use the notation H(z) instead of H(s)(z), but wewill assume that H(z) is primitive.We shall use several times the following observation: if a transfer function H(z) has theform H(z) =∑nk=1ek(z−s)k , then the impulse response sequence ishm =n∑k=1(m− 1k − 1)sm−kek m = 1, 2, 3, . . . , (3)(with the convention that for α < β we define(αβ):= 0, and(00):= 1 and 00 := 1).This can be proved by using for each k = 1, 2, . . . n and |z| > |s| the formula 1(z−s)k =(1z)k[1+ sz+( sz)2+( sz)3+ . . . ]k. Since we have∑nk=1ek(z−s)k = H(z) =∑∞m=1hmzm, comparingcoefficients yields the stated formula.Recall now that for our transfer function H(z) the only dominant pole is λ1 = 1.If c(n1)1 < 0, then it is clear (see the explicit formula for hk,1 and the estimates on hkbelow) that for large k we have hk < 0, hence there exists no positive realization of H(z).Therefore we can assume that c(n1)1 > 0, and, without loss of generality, that c(n1)1 = 1.27. November, 2002 DRAFT5We use the notation h(i)k,j for the impulse response sequence corresponding to the functionc(i)j(z−λj)i , and we let hk,j :=∑nji=1 h(i)k,j. Hence we shall have hk =∑rj=1 hk,j. The idea behindthe forthcoming calculations is that hk,1 turns out to be ’dominant’ in the long termbehaviour of hk, and h(n1)k,1 will be ’dominant’ in the long term behaviour of hk,1.First we find an index N0 such that hk,1 ≥ 1 for all k ≥ N0. If n1 = 1, then hk,1 = 1for all k ≥ 1, therefore we can take N0 = 1. Assume n1 > 1. Since we have λ1 = 1 andc(n1)1 = 1, it follows from (3) thathk,1 =(k−1n1−1)+ c(n1−1)1(k−1n1−2)+ ...+ c(2)1(k−11)+ c(1)1 =1 + [(k−1n1−1)+ c(n1−1)1(k−1n1−2)+ ... + c(2)1(k−11)+ c(1)1 − 1].Let C := max{|c(n1−1)1 |, . . . , |c(2)1 |, |c(1)1 |+1}, and assume that k ≥ (n1C+1)(n1−1) =: N0.(Note that N0 ≥ 3n1 − 3 ≥ 2n1 − 1, so the finite sequence of the binomial coefficients(k−1n1−1),(k−1n1−2), ...,(k−10)is strictly monotonically decreasing.) For k this large we havek−n1+1n1−1 ≥ n1C and this means that(k−1n1−1)/(k−1n1−2) ≥ n1C. Hence for any 2 ≤ j ≤ n1 wehave(k−1n1−1)/(k−1n1−j) ≥ n1C. This means thathk,1 = 1 + [(k−1n1−1)+ c(n1−1)1(k−1n1−2)+ ...+ c(2)1(k−11)+ c(1)1 − 1] ≥1 +(k−1n1−1)− C∑n1−1i=1 (k−1i−1) ≥ 1,as desired.Next we find an index M0 such that∑rj=2 |hk,j| ≤ 1 for all k > M0. It follows from (3)that hk,j =∑nji=1 h(i)k,j =∑nji=1 c(i)j λk−ij(k−1i−1). Therefore|hk,j| ≤∑nji=1 |h(i)k,j| ≤∑nji=1 |c(i)j ||λj|k−i(k−1i−1).Now, there are altogether N1 :=∑rj=2 nj coefficients of type h(i)k,j so it is enough toensure that the modulus of each of them is not greater than 1/N1. That is, we want|c(i)j ||λj|k−i(k−1i−1) ≤ 1N1to hold. If λj = 0, then this is obviously true for k ≥ i+1. Assumeλj 6= 0. To simplify forthcoming calculations we use the notation ρ := max{|λj| : j =2, 3, . . . , r}, γ := max{|c(i)j | : j = 2, 3, . . . , r; i = 1, 2, . . . , nj} and η := max{nj : j =2, 3, . . . , r}. The desired inequality |c(i)j ||λj|k−i(k−1i−1) ≤ 1N1is implied by γρk−i(k−1i−1) ≤ 1N1,which is equivalent to ρk/2ρk/2(k−1i−1) ≤ ρiN1γ. It is easy to check that for fixed i the value ofρk/2(k−1i−1)is monotonically decreasing (in k) for k ≥ i−11−ρ1/2 + 1 =: N(i).27. November, 2002 DRAFT6We use the notation C(i) := ρ(N(i)/2)(N(i)−1i−1). If k ≥ N (i), then it is sufficient that(ρ1/2)k ≤ ρiN1γC(i)=: K(i),therefore we can takeM0 := max{η,N (i), logρ1/2K(i), i = 1, 2, . . . η}(Note that we include η because we must not forget about the possible pole at 0.)This means that if k > max{N0,M0}, then hk ≥ hk,1 −∑rj=2 |hk,j| ≥ 0, therefore itis definitely enough to check the nonnegativity of the first K0 := max{N0,M0} terms ofthe impulse response sequence. We remark that the calculations above show that it issufficient to know an upper bound on the values of ρ and γ instead of the exact values.This means that it is enough to determine the approximate locations of the poles, and theapproximate values of the partial fraction coefficients, and Step 3 can already be applied.III. A nonnegative realizationStep 4. In order to construct a positive realization of H(z) it is sufficient to find apositive realization for each H(s)(z) in the decomposition (2), and then apply the methodof [4] or [8]. As all H(s)(z) are primitive, we could find a positive realization of H(s)(z)by applying the results of [1] directly. However, in this Step we propose a method forreducing the pole order of the dominant pole, and apply the construction of [1] only whenthe dominant pole is simple. This step seems to simplify the construction of [1] in the casewhen H(s)(z) has a multiple dominant pole.Again, for the sake of simplicity we will use the notation H(z) instead of H(s)(z), andwe will assume that H(z) is primitive with nonnegative impulse response.Denote in this Step the transfer function corresponding to the “shifted” impulse responsesequence (hk, hk+1, ...) by Hk(z), i.e.Hk(z) :=∑∞j=1 hj+k−1z−j .(Note that H(z) = H1(z).) First we make the following observation:If Hk(z) has a nonnegative realization in N dimensions (for some k > 1), then so doesHk−1(z) in N + 1 dimensions. (For an easy proof see [7].)In order to construct a nonnegative realization of H(z) we take the following guidelines:27. November, 2002 DRAFT7We find a positive integer k1 so that we can construct a nonnegative realization ofHk1(z)in some dimensions N , and then use the obsevation above to construct a nonnegativerealization of H(z) in dimensions N + k1 − 1. The index k1 will be chosen so thatHk1(z) = [∑n1m=1em(z−1)m ] + [d1z−1 +∑rj=2∑nji=1e(i)j(z−λj)i ] =: f(3)(z) + f (4)(z)holds, where em ≥ 0 for all 1 ≤ m ≤ n1 and the whole impulse response sequence of f (4)(z)is nonnegative. Here f (3)(z) has a trivial positive realization in n1 dimensions, and f(4)(z)is primitive with a simple dominant pole and nonnegative impulse response, so that theconstruction of [1], Theorem 4.1 can be applied. We remark that the index k1 is notuniquely determined. In order to minimize the dimension of the realization it is importantto determine the optimal value of k1. This may be easy to do for a particular transferfunction H(z), but a general formula for the optimal value of k1 does not seem possibleto find. The proof below shows only that such an index k1 always exists.Assume that n1 > 1. Notice thatHk(z) =∑∞j=1 hk+j−1z−j = zk−1H(z)−∑k−1r=1 hrzk−r−1.In particular, the formula H2(z) = zH(z) − h1 = (z − 1)H(z) +H(z)− h1 shows that inthe partial fraction decomposition of H2(z) the part corresponding to λ1 = 1 is given by1(z−1)n1 +∑n1−1j=1c(j)1 +c(j+1)1(z−1)jFrom this it follows by mathematical induction with respect to k that the partial fractionpart corresponding to λ1 = 1 in Hk(z) is given by1(z − 1)n1 +n1−1∑j=1∑n1i=j(k−1i−j)c(i)1(z − 1)j . (4)If k ≥ N0 (as in Step 3), then(k−1j)/(k−1j−1) ≥ n1C for all j = 1, 2, ..., n1 − 1. Therefore,we see as in Step 3 that for k ≥ N0 all the numerators in (4) are not less than 1. Therefore,HN0(z) = [∑n1m=2dm(z−1)m ] + [d1z−1 +∑rj=2∑nji=1d(i)j(z−λj)i ] =: f(1)(z) + f (2)(z).where dm ≥ 1 for all 1 ≤ m ≤ n1 (but the first few terms of the impulse responsesequence of f (2)(z) may be negative!). We can now apply the estimates of Step 3 to thefunction f (2)(z) and obtain a number K˜0 such that the impulse response of the function27. November, 2002 DRAFT8f (2)(z) becomes nonnegative if k ≥ K˜0. This means that we can take k1 = N0 + K˜0, andthe desired decomposition Hk1(z) = f(3)(z) + f (4)(z) holds. Now, the construction of [1]Theorem 4.1 applies to f (4)(z), and we can combine the positive realizations of f (3)(z) andf (4)(z) to get a positive realization of Hk1(z). Then we obtain a positive realization ofH(z) by applying the observation above.We remark that an upper bound on the dimension of the positive realization constructedin [1], Theorem 4.1 has not yet been presented. Such dimension estimate is possible toprove, but the proof is fairly long and mathematically involved. It will be presented in aforthcoming publication.IV. ExampleWe illustrate the steps of the algorithm by the following example. LetH(z) =1+z− 14z2z3−1 +1/3(z−1)2 +2/3z−1/2 +2/3z+1/2− 1z+0.9.Then p = 3, and the downsampled functions areH(0)(z) =−1/4z−1 +1(z−1)2 +2/3z−1/8 +2/3z+1/8− 1z+0.729H(1)(z) =4/3z−1 +1(z−1)2 +1/3z−1/8 − 1/3z+1/8 + 0.9z+0.729H(2)(z) =5/3z−1 +1(z−1)2 +1/6z−1/8 +1/6z+1/8− 0.81z+0.729.Next we check the nonnegativity of the impulse response of H(z) by checking eachH(s)(z) (s = 0, 1, 2). For s = 0 we have n1 = 2, C = 5/4 and N0 = 4. Further, followingthe definitions, we have N1 = 3, ρ = 0.729, γ = 1, η = 1, N(1) = 1, C(1) =√0.729,K(1) =√0.7293, and logρ1/2K(1) = 7.95. Therefore it is sufficient to check the first K0 = 8terms of the impulse response sequence ofH(0)(z). Following similar calculations we deducethat it is enough to check the nonnegativity of the first 8 terms of the impulse responsesequence of H(1)(z), and the first 7 terms in the case of H(2)(z).Next, we construct a positive realization for each H(s)(z) (s = 0, 1, 2). In the caseof H(0)(z) we apply Step 4 with k1 = 1. The ’shifted’ transfer function is given byH2(z) = [1(z−1)2 ] + [3/4z−1 +1/12z−1/8 − 1/12z+1/8 + 0.729z+0.729 ] = f (3)(z) + f (4)(z)Now, a positive realization for f (3)(z) and f (4)(z) is given by the triplets (c3, A3, b3) and(c4, A4, b4), respectively, where27. November, 2002 DRAFT9c3 =(1 0), A3 = 1 10 1 , b3 = 01,c4 =(1 0 1.972 0.291412), A4 =0 1 0 01/64 0 0 00 0 0 0.7290 0 1 0.271, b4 =01/480.750.Now we apply the observation of Step 4 (cf. [7]), and give a positive realization ofH(0)(z) byc0 =(1/12, c3, c4), A0 =0 0 0b3 A3 0b4 0 A4 , b0 =100The cases of H(1)(z) and H(2)(z) are ’trivial’ because we can take k1 = 0 in Step 4, andwriteH(1)(z) = [1(z−1)2 ] + [4/3z−1 +1/3z−1/8 − 1/3z+1/8 + 0.9z+0.729 ]H(2)(z) = [1(z−1)2 ] + [5/3z−1 +1/6z−1/8 +1/6z+1/8− 0.81z+0.729].In both cases positive realizations of dimension 6 can be given as follows:c1 =(1 0 1 0 1.675 0.507925), A1 =1 1 0 0 0 00 1 0 0 0 00 0 0 1 0 00 0 1/64 0 0 00 0 0 0 0 0.7290 0 0 0 1 0.271, b1 =0101/124/30c2 =(1 0 1 0 0.514 1.354294), A2 =1 1 0 0 0 00 1 0 0 0 00 0 0 1 0 00 0 1/64 0 0 00 0 0 0 0 0.7290 0 0 0 1 0.271, b2 =011/305/30A positive realization of the original H(z) is then possible to construct as in [4] or [8].27. November, 2002 DRAFT10V. ConclusionsIn this brief we provided a general finite step procedure for checking the nonnegativity ofthe impulse response sequence of H(z), which answers an open problem raised in [1]. Forprimitive transfer functions a new method of positive realization was proposed by reducingthe pole order of the dominant pole.References[1] B.D.O. Anderson, M. Deistler, L. Farina, and L. Benvenuti, “Nonnegative realization of a linear system withnonnegative impulse response,” IEEE Trans. Circuits Syst. I, vol. 43, no.2, pp 134-142, Feb. 1996.[2] L. Benvenuti, L. Farina, and B.D.O. Anderson, “Filtering through a combination of positive filters,” IEEETrans. Circuits Syst. I, vol. 46, no.12, pp 1431-1440, Dec. 1999.[3] L. Benvenuti and L. Farina, “An example of how positivity may force realizations of ’large’ dimensions,” Syst.Contr. Lett., vol. 36, no.4, pp 261-266, 1999.[4] L. Farina, “On the existence of a positive realization,” Syst. Contr. Lett., vol. 28, no.4, pp 219-226, 1996.[5] L. Farina and S. Rinaldi, Positive linear systems: Theory and applications. New York: Wiley, 2000.[6] K.-H. Foerster, B. Nagy, “Nonnegative realizations of matrix transfer functions,” Linear Algebra Appl., vol.311, no.1-3, pp 107-129, 2000.[7] C. Hadjicostis, “Bounds on the size of minimal nonnegative realizations for discrete-time LTI systems,” Syst.Contr. Lett., vol. 37, no.1, pp 39-43, 1999.[8] T. Kitano and H. Maeda, “Positive realization of discrete-time system by geometric approach,” IEEE Trans.Circuits Syst. I , vol. 45, no.3, pp 308-311, Mar. 1998.[9] Y. Ohta, H. Maeda, and S. Kodama, “Reachability, observability and realizability of continuous-time positivesystems,” SIAM J.Control Optim., vol. 22, no.2, pp. 171-180, 1984.27. November, 2002 DRAFT

An example of how positivity may force realizations of ’large’ dimensions,”

Bounds on the size of minimal nonnegative realizations for discrete-time LTI systems,”

Filtering through a combination of positive ﬁlters,”

Nonnegative realization of a linear system with nonnegative impulse response,”

Nonnegative realizations of matrix transfer functions,”

On the existence of a positive realization,” Syst.

Positive linear systems: Theory and applications.

Positive realization of discrete-time system by geometric approach,”

Reachability, observability and realizability of continuous-time positive systems,”

http://real.mtak.hu/7970/1/algorithm%20for.pdf

Algorithm for positive realization of transfer functions

Abstract

Similar works

Full text

Available Versions

Crossref

Repository of the Academy's Library