# Frequency Weighted Model Order Reduction Technique and Error Bounds for Discrete Time Systems

Authors: Muhammad Imran; Abdul Ghafoor; Victor Sreeram
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Model reduction is a process of approximating higher order original models by comparatively lower order models with reasonable accuracy in order to provide ease in design, modeling and simulation for large complex systems. Generally, model reduction techniques approximate the higher order systems
for whole frequency range. However, certain applications (like controller reduction) require frequency weighted approximation, which introduce the concept of using frequency weights in model reduction techniques. Limitations of some existing frequency weighted model reduction techniques include lack of stability of reduced order models (for two sided weighting case) and frequency response error bounds. A new frequency weighted technique for balanced model reduction for discrete time systems is proposed. The proposed technique guarantees stable reduced order models even for the case when two sided weightings are present. Efficient technique for frequency weighted Gramians is also proposed. Results are compared with other existing frequency weighted model reduction techniques for discrete time systems. Moreover, the proposed technique yields frequency response error bounds.


---

## Body

## 1. Introduction

Model reduction has played a significant role in modern control system design and caught a lot of attention in the last few decades [1–5]. It is desirable that reduced order model preserves the fundamental properties of original system like stability, passivity, and so forth. Moreover, the approximation error between original and reduced order system is required to be small. Balanced truncation [6] is prominent model reduction scheme, which not only ensures stability of reduced order systems but also provides frequency responsea priori error bounds. Other schemes like Hankel norm approximation [7], Pade approximation [8], Krylov technique [9, 10], linear matrix inequality (LMI) technique [11, 12], and so forth, are also useful for model reduction. Hankel norm approximation has complex implementation and does not preserve steady state [13]. Pade and Krylov approximation sometimes provide unstable reduced order models and there exist no global error bounds [14, 15]. LMI technique is based on mathematical iterative methods (i.e., convex optimization and bisection algorithm); therefore, much computational power is consumed [16]. However, LMI technique has been applied on various model reduction problems for different types of systems (including time delay [17], discrete state delay [18, 19], switched hybrid [20], and nonlinear stochastic [21]). Applications of model reduction are not only restricted to control engineering but also find utility in medical [22] and text summarization [23] areas, and so forth. Most model reduction algorithms tend to minimize the reduction error over the whole frequency range; however there are certain situations (like controller reduction [24]), wherein the approximation error is more critical over a certain band of frequencies.Enns [25] has enhanced the balanced truncation [6] to incorporate frequency weights for model reduction of higher order systems. This technique can work with input, output, and two sided weightings. However, in two sided weighting case, this technique may give unstable reduced order models [26]. To circumvent instability issue in the presence of two sided weightings, many techniques appear in the literature (including Lin and Chiu [27], Wang et al. [28], Varga and Anderson [29], Ghafoor and Sreeram [30], etc.).In Lin and Chiu [27] technique weightings are strictly proper. It was improved to include more general proper weightings in [26]. Varga and Anderson [29] technique yields proper model for strictly proper systems also. Ghafoor and Sreeram [30] technique is a parameterized technique. Wang et al.'s technique is relatively useful, since other techniques [27, 29] are not applicable for controller order reduction. Moreover, Wang et al.'s technique yields easily computable expression fora priori error bounds.Various partial fraction based techniques appear in the literature [31–35] that works for continuous as well as discrete time systems. Unfortunately, most of these techniques yield large frequency response error [30] as compared to Enns technique. However, [35] incorporates free parameters to reduce the approximation error.Most frequently, frequency weighted model reduction problem is treated in continuous time; however, there are few papers (like [36–38]) which deal with this problem explicitly in discrete time. Sahlan et al.'s [36] technique (discrete time version of [31]) is a modified version of Lin and Chiu [27] technique. Campbell et al.'s [37] technique is discrete time counterpart of Wang et al.'s [28] technique, which not only provides stability for two sided case but also gives easily computable error bounds. Campbell et al.'s [37] technique involves taking absolute values of eigenvalues. This may introduce a large change to the system and hence a larger error if it contains some negative eigenvalues. Varga and Anderson's [38] tend to minimize the distance between Enns and Campbell et al.'s techniques Gramians by eliminating negative eigenvalues.In this work, a technique is developed for frequency weighted balanced model order reduction for discrete time systems. The large change in eigenvalues is avoided by yielding similar effect on all eigenvalues. The stability is guaranteed even for double sided weighting. The proposed technique provides comparable frequency response error and yields easily computablea priori error bounds. Numerical examples are given to show the usefulness and comparison of proposed technique with the existing frequency weighted balanced reduction techniques.


## 2. Preliminaries

In this section we review some of the existing frequency weighted model reduction techniques for discrete time systems which include Enns [25], Campbell et al.’s [37], and Varga and Anderson's [38].Consider a stable full order original system with transfer functionH(z)=C(zI-A)-1B+D, a stable input weighting system with input transfer function Vi(z)=Ci(zI-Ai)-1Bi+Di, and a stable output weighting system with output transfer function Wo(z)=Co(zI-Ao)-1Bo+Do; the augmented systems are given by
(1)H(z)Vi(z)=Cdi(zI-Adi)-1Bdi+Ddi,Wo(z)H(z)=Cdo(zI-Ado)-1Bdo+Ddo,
where
(2)[AdiBdiCdiDdi]=[ABCiBDi0AiBiCDCiDDi],[AdoBdoCdoDdo]=[AoBoCBoD0ABCoDoCDoD].Let the Gramians(3)Pdi=[PenP12P12TPV],Qdo=[QWQ12TQ12Qen]
satisfy the following Lyapunov equations:
(4)AdiPdiAdiT-Pdi+BdiBdiT=0,(5)AdoTQdoAdo-Qdo+CdoTCdo=0.Expanding the (1,1) and (2,2) blocks of above equations, we get(6)APenAT-Pen+Xen=0,ATQenA-Qen+Yen=0,
where
(7)Xen=AP12CiTBT+BCiP12TAT+BCiPVCiTBT+BDiDiTBT,Yen=CTBoTQ12TA+ATQ12BoC+CTBoTQWBoC+CTDoTDoC.


### 2.1. Enns Technique [25]

LetT be contragredient matrix obtained as
(8)TTQenT=T-1PenT-T=diag⁡{σ1,σ2,σ3,…,σn},
where σj≥σj+1,j=1,2,3,…,n-1 and σl>σl+1. By transforming and then partitioning the original system, we have
(9)A^=T-1AT=[A11A12A21A22],B^=T-1B=[B1B2],C^=CT=[C1C2],D^=D,
where A11∈Rl×l. The reduced order system is obtained as follows: Hl(z)=C1(zI-A11)-1B1+D.Remark 1.
It is not guaranteed to ensure thatXen≥0 and Yen≥0; the reduced order models obtained using Enns technique may not remain stable for both sided weightings.


### 2.2. Campbell et al.’s Technique [37]

Campbell et al.’s [37] (a discrete time version of [28]) guarantees the positive semidefiniteness of symmetric matrices Xen and Yen to ensure stability. Let the new controllability PCS and observability QCS Gramians, respectively, be calculated by solving the following Lyapunov equations:
(10)APCSAT-PCS+BCSBCST=0,ATQCSA+QCS+CCSTCCS=0.
which are used to obtain contragredient matrix T as
(11)TTQCST=T-1PCST-T=Σ,
where Σ=diag⁡{σ1,σ2,σ3,…,σn} and σj≥σj+1,j=1,2,3,…,n-1 and σl>σl+1. The fictitious input BCS and output CCS matrices shown in the above Lyapunov equations are defined as BCS=UCS|SCS|1/2 and CCS=|RCS|1/2VCST, respectively. Since the expressions UCS,  SCS,  VCS, and RCS are calculated by orthogonal eigen decomposition Xen=UCSSCSUCST  and  Yen=VCSRCSVCST,  where  SCS=diag⁡(s1,s2,…,sn),  RCS=diag⁡(r1,r2,…,rn),  |s1|≥|s2|≥⋯≥|sn|≥0 and |r1|≥|r2|≥⋯≥|rn|≥0. The reduced order systems are calculated by transforming and then partitioning the original system.Remark 2.
The stability of reduced order models in the presence of both input and output weightings is guaranteed and the following error bound holds [37]
(12)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LCS∥∞∥KCSVi(z)∥∞∑j=l+1n‍σj,
where
(13)LCS=CVCSdiag⁡(|r1|-1/2,|r2|-1/2,…,|rli|-1/2,0,…,0),KCS=diag⁡(|s1|-1/2,|s2|-1/2,…,|sko|-1/2,0,…,0)UCSTB.li=rank⁡[Xen]  and  ko=rank⁡[Yen].


### 2.3. Varga and Anderson's Technique [38]

Note that the Gramians satisfyPen≤PCS and Qen≤QCS. For minimizing the distances between the Gramians, Pen-PCS and Qen-QCS, Varga and Anderson proposed the following technique.Let new controllability and observability GramiansPVd and QVd, respectively, be calculated as the solutions to Lyapunov equations
(14)APVdAT-PVd+BVdBVdT=0,ATQVdA-QVd+CVdTCVd=0
which are used obtain contragredient matrix T as
(15)TTQVdT=T-1PVdT-T=Σ,
where Σ=diag⁡{σ1,σ2,σ3,…,σn}, and σj≥σj+1,  j=1,2,…,n-1, σl>σl+1. The new fictitious input BVd and output CVd matrices in the above Lyapunov equations are defined as BVd=UVd1SVd11/2 and CVd=RVd11/2VVd1T, respectively. The terms UVd1,  SVd1,  VVd1  , and RVd1 are calculated from the orthogonal eigen decomposition of symmetric matrices
(16)Xen=[UVd1UVd2][SVd100SVd2][UVd1TUVd2T],Yen=[VVd1VVd2][RVd100RVd2][VVd1TVVd2T],
where [SVd100SVd2]=diag⁡{s1,s2,…,sn},  [RVd100RVd2]=diag⁡{r1,r2,r3,…,rn},  s1≥s2≥s3≥⋯≥sn,  r1≥r2≥r3≥⋯≥rn,  SVd1=diag⁡{s1,s2,s3,…,sk},  RVd1=diag⁡{r1,r2,r3,…,rk},  s1≥s2≥s3≥⋯≥sk≥0,  r1≥r2≥r3≥⋯≥rk≥0.The reduced order systems are calculated by transforming and then partitioning the original system.Remark 3.
The stability of the reduced system is guaranteed and the following error bound holds [38]
(17)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LVd∥∞∥KVdVi(z)∥∞∑j=l+1n‍σj,
where LVd=CVVd1RVd1-1/2 and KVd=SVd1-1/2UVd1TB.


## 2.1. Enns Technique [25]

LetT be contragredient matrix obtained as
(8)TTQenT=T-1PenT-T=diag⁡{σ1,σ2,σ3,…,σn},
where σj≥σj+1,j=1,2,3,…,n-1 and σl>σl+1. By transforming and then partitioning the original system, we have
(9)A^=T-1AT=[A11A12A21A22],B^=T-1B=[B1B2],C^=CT=[C1C2],D^=D,
where A11∈Rl×l. The reduced order system is obtained as follows: Hl(z)=C1(zI-A11)-1B1+D.Remark 1.
It is not guaranteed to ensure thatXen≥0 and Yen≥0; the reduced order models obtained using Enns technique may not remain stable for both sided weightings.


## 2.2. Campbell et al.’s Technique [37]

Campbell et al.’s [37] (a discrete time version of [28]) guarantees the positive semidefiniteness of symmetric matrices Xen and Yen to ensure stability. Let the new controllability PCS and observability QCS Gramians, respectively, be calculated by solving the following Lyapunov equations:
(10)APCSAT-PCS+BCSBCST=0,ATQCSA+QCS+CCSTCCS=0.
which are used to obtain contragredient matrix T as
(11)TTQCST=T-1PCST-T=Σ,
where Σ=diag⁡{σ1,σ2,σ3,…,σn} and σj≥σj+1,j=1,2,3,…,n-1 and σl>σl+1. The fictitious input BCS and output CCS matrices shown in the above Lyapunov equations are defined as BCS=UCS|SCS|1/2 and CCS=|RCS|1/2VCST, respectively. Since the expressions UCS,  SCS,  VCS, and RCS are calculated by orthogonal eigen decomposition Xen=UCSSCSUCST  and  Yen=VCSRCSVCST,  where  SCS=diag⁡(s1,s2,…,sn),  RCS=diag⁡(r1,r2,…,rn),  |s1|≥|s2|≥⋯≥|sn|≥0 and |r1|≥|r2|≥⋯≥|rn|≥0. The reduced order systems are calculated by transforming and then partitioning the original system.Remark 2.
The stability of reduced order models in the presence of both input and output weightings is guaranteed and the following error bound holds [37]
(12)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LCS∥∞∥KCSVi(z)∥∞∑j=l+1n‍σj,
where
(13)LCS=CVCSdiag⁡(|r1|-1/2,|r2|-1/2,…,|rli|-1/2,0,…,0),KCS=diag⁡(|s1|-1/2,|s2|-1/2,…,|sko|-1/2,0,…,0)UCSTB.li=rank⁡[Xen]  and  ko=rank⁡[Yen].


## 2.3. Varga and Anderson's Technique [38]

Note that the Gramians satisfyPen≤PCS and Qen≤QCS. For minimizing the distances between the Gramians, Pen-PCS and Qen-QCS, Varga and Anderson proposed the following technique.Let new controllability and observability GramiansPVd and QVd, respectively, be calculated as the solutions to Lyapunov equations
(14)APVdAT-PVd+BVdBVdT=0,ATQVdA-QVd+CVdTCVd=0
which are used obtain contragredient matrix T as
(15)TTQVdT=T-1PVdT-T=Σ,
where Σ=diag⁡{σ1,σ2,σ3,…,σn}, and σj≥σj+1,  j=1,2,…,n-1, σl>σl+1. The new fictitious input BVd and output CVd matrices in the above Lyapunov equations are defined as BVd=UVd1SVd11/2 and CVd=RVd11/2VVd1T, respectively. The terms UVd1,  SVd1,  VVd1  , and RVd1 are calculated from the orthogonal eigen decomposition of symmetric matrices
(16)Xen=[UVd1UVd2][SVd100SVd2][UVd1TUVd2T],Yen=[VVd1VVd2][RVd100RVd2][VVd1TVVd2T],
where [SVd100SVd2]=diag⁡{s1,s2,…,sn},  [RVd100RVd2]=diag⁡{r1,r2,r3,…,rn},  s1≥s2≥s3≥⋯≥sn,  r1≥r2≥r3≥⋯≥rn,  SVd1=diag⁡{s1,s2,s3,…,sk},  RVd1=diag⁡{r1,r2,r3,…,rk},  s1≥s2≥s3≥⋯≥sk≥0,  r1≥r2≥r3≥⋯≥rk≥0.The reduced order systems are calculated by transforming and then partitioning the original system.Remark 3.
The stability of the reduced system is guaranteed and the following error bound holds [38]
(17)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LVd∥∞∥KVdVi(z)∥∞∑j=l+1n‍σj,
where LVd=CVVd1RVd1-1/2 and KVd=SVd1-1/2UVd1TB.


## 3. Main Results

In Campbell et al.’s [37] technique, the symmetric matrices Xen and Yen are guaranteed positive/semipositive definite by taking the square root of absolute values of the eigenvalues obtained by eigen decomposition of symmetric matrices Xen and Yen. This may lead to a large change in some eigenvalues and may not affect other eigenvalues. Although in Varga and Anderson's [38] technique, this large change was slightly improved by eliminating negative eigenvalues, but the problem persists with the other eigenvalues. In the following, a new technique is proposed in which a similar effect on all eigenvalues of indefinite matrices Xen and Yen guarantees stability, error bound, and improved frequency response error.


### 3.1. Proposed Technique

Let a new controllabilityPIG and observability QIG Gramians, respectively, be calculated by solving the following Lyapunov equations:
(18)APIGAT-PIG+BIGBIGT=0,(19)ATQIGA-QIG+CIGTCIG=0.The matricesBIG and CIG are the new fictitious input and output matrices, respectively, and are defined as
(20)BIG={UIG(SIG-snI)1/2forsn<0UIGSIG1/2forsn≥0,CIG={(RIG-rnI)1/2VIGTforrn<0RIG1/2VIGTforrn≥0.The termsUIG,  SIG,  VIG, and RIG are calculated by the orthogonal eigen decomposition of symmetric matrices Xen=UIGSIGUIGT and Yen=VIGRIGVIGT, where SIG=diag⁡(s1,s2,s3,…,sn), RIG=diag⁡(r1,r2,r3,…,rn), s1≥s2≥⋯≥sn, and r1≥r2≥⋯≥rn.Let a contragradient transformation matrixT be obtained as
(21)TTQIGT=T-1PIGT-T=Σ,
where Σ=diag⁡{σ1,σ2,σ3,…,σn}, and σj≥σj+1,  j=1,2,…,n-1, σl>σl+1. The reduced order system is calculated by transforming and partitioning the original system.Remark 4.
SinceXen≤BIGBIGT≥0, Yen≤CIGTCIG≥0, PIG>0 and QIG>0, therefore, the realization (A,BIGS,CIGS) is minimal. Moreover, the reduced order models are stable.Theorem 5.
The following error bound for the proposed technique holds if the rank conditionsrank⁡[BIGB]=rank⁡[BIG] and rank⁡[CIGC]=rank⁡[CIG] (which follows from [2, 28, 29, 38]) are satisfied
(22)(i)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LIG∥∞∥KIGVi(z)∥∞∑j=l+1n‍σj(ii)∥(H(z)-Hl(z))Vi(z)∥∞≤2∥KIGVi(z)∥∞∑j=l+1n‍σj(iii)∥Wo(z)(H(z)-Hl(z))∥∞≤2∥Wo(z)LIG∥∞∑j=l+1n‍σj,
where
(23)LIG={CV(RIG-rnI)-1/2forrn<0CVRIG-1/2forrn≥0,KIG={(SIG-snI)-1/2UTBforsn<0SIG-1/2UTBforsn≥0.Proof.
We show proof of (i) (whereas (ii) and (iii) are special cases of (i)). Since the rank conditionsrank⁡[BIGB]=rank⁡[BIG] and rank⁡[CIGC]=rank⁡[CIG] are satisfied, the relationships B=BIGKIG and C=LIGCIG hold. By partitioning BIG=[BIG1BIG2],  CIG=[CIG1CIG2] and then substituting B1=BIG1KIG, C1=LIGCIG1, respectively, yields
(24)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞=∥Wo(z)(C(zI-A)-1B-C1(zI-A11)-1B1)Vi(z)∥∞=∥Wo(z)(-LIGCIG1(zI-A11)-1BIG1KIGLIGCIG(zI-A)-1BIGKIGmmmmmmn-LIGCIG1(zI-A11)-1BIG1KIG)Vi(z)∥∞=∥Wo(z)LIG(-CIG1(zI-A11)-1BIG1CIG(zI-A)-1BIGmmmmmmmmn-CIG1(zI-A11)-1BIG1)KIGVi(z)∥∞=∥Wo(z)LIG∥∞×∥(CIG(zI-A)-1BIG-CIG1(zI-A11)-1BIG1)∥∞×∥KIGVi(z)∥∞.
If {A11,BIG1,CIG1} is reduced order model obtained by partitioning a balanced realization {A,BIG,CIG}, we have [7, 39]
(25)∥(CIG(zI-A)-1BIG-CIG1(zI-A11)-1BIG1)∥∞≤2∑j=l+1n‍σj.
Therefore,
(26)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LIG∥∞∥KIGVi(z)∥∞∑j=l+1n‍σj.Remark 6.
For the case when symmetric matricesXen≥0 and Yen≥0, then Pen=PCS=PVd=PIG and Qen=QCS=QVd=QIG. However, when matrices Xen and Yen are indefinite, then Pen<PIG and Qen<QIG. Moreover, frequency weighted Hankel singular values satisfy (λj[PenQen])1/2≤(λj[PIGQIG])1/2.Remark 7.
For the case when inputVi(z) weights are co-inner and output Wo(z) weights are inner [40], then P=Pen=PCS=PVd=PIG and Q=Qen=QCS=QVd=QIG, where P and Q are unweighted Gramians defined as
(27)APAT-P+BBT=0,ATQA-Q+CTC=0.Remark 8.
For the case when the symmetric matricesXen≥0 and Yen≥0, the reduced order models obtained using Enns [25], Campbell et al.'s [37], Varga and Anderson's [38], and proposed technique are same.


### 3.2. Computational Aspects

The frequency weighted balanced truncation model reduction techniques balance the original system and then truncate the balanced realization to get the desired reduced order system. The balancing procedure involves computation of transformation matrix from controllability and observability Gramians. Sometimes these matrices become numerically low rank especially in large scale systems (possibly) due to rapid decay of their eigenvalues [3]. Due to this reason, balancing procedure becomes inefficient. Accuracy enhancing techniques for different frequency weighted model reduction techniques appear in [4, 29].For unweighted case, Hammarling's technique [41] is used to obtained Cholesky factors of Gramian matrices from original system realization without actually computing controllability and observability Gramian matrices, respectively.In frequency weighted techniques, Cholesky factors of the Gramian matrices are obtained from the augmented system realizations. LetS¯ and R¯ be the Cholesky factors of the augmented system Gramians matrices Pbi and Qbo of (4) and (5), respectively,
(28)Pbi=S¯S¯T=[S11S120S22][S11T0S12TS22T]=[S11S11T+S12S12TS12S22TS22S12TS22S22T]=[PenP12P12TPV],Qbo=R¯TR¯=[R11T0R12TR22T][R11R120R22]=[R11TR11R11TR12R12TR11R22TR22+R12TR12]=[QWQ12TQ12Qen].By making use of Cholesky factorsS¯ and R¯ calculated above, the Cholesky factors corresponding to Gramians in frequency weighted model reduction techniques like Enns [25], Campbell et al.'s [37], Varga and Anderson's [38], and proposed technique can be obtained as follows:(
1) Enns Technique. The Cholesky factors Sen=[S11S12] and Ren=[R12R22] satisfy [29]
(29)Pen=SenSenT=S11S11T+S12S12T=[S11S12][S11TS22T],Qen=RenTRen=R22TR22+R12TR12=[R22TR12T][R22R12].(
2
) Campbell et al.'s Technique. The Cholesky factors S¯CS and R¯CS satisfy PCS=S¯CSS¯CST and QCS=R¯CSTR¯CS [4].(
3
) Varga and Anderson's Technique. The Cholesky factors S¯Vd and R¯Vd satisfy PVd=S¯VdS¯VdT and QVd=R¯VdTR¯Vd [4].(
4
) Proposed Technique. The Cholesky factors S¯IG and R¯IG satisfy PIG=S¯IGS¯IGT and QIG=R¯IGTR¯IG.In the following we establish a relationship between Cholesky factors of Gramian matrices used in Enns and proposed techniques. Equations (18) and (19) can be expressed as
(30)A(Pen+Pad)AT-(Pen+Pad)+(Xen-snI)=0,forsn<0,APenAT-Pen+Xen=0,forsn≥0,AT(Qen+Qad)A-(Qen+Qad)+(Yen-rnI)=0,forrn<0,ATQenA-Qen+Yen=0,forrn≥0,APadAT-Pad-snI=0,forsn<0,ATQadA-Qad-rnI=0,forrn<0.
Since
(31)XIG=U(S-snI)1/2(S-snI)1/2UT=Xen-snI,forsn<0,XIG=U(S)1/2(S)1/2UT=Xen,forsn≥0,YIG=VT(R-rnI)1/2(R-rnI)1/2V=Yen-rnI,forrn<0YIG=VT(R)1/2(R)1/2V=Yen,forrn≥0.
By using Hammarling technique to calculate Cholesky factors of Gramians Pad and Qad from realization {A,-snI,-rnI,D}, we can write Pad=S¯adS¯adT and Qad=R¯adTR¯ad. Therefore, frequency weighted controllability PIG (18) and observability QIG (19) Gramians can be expressed as
(32)PIG=S¯IGS¯IGT=Pen+Pad=S11S11T+S12S12T+SadSadT=[S11S12Sad][S11TS12TSadT],(33)QIG=R¯IGTR¯IG=Qen+Qad=R22TR22+R12TR12+RadTRad=[R22TR12TRadT][R22R12Rad].Remark 9.
Note that, Cholesky factors for Enns and proposed technique are computed directly from augmented system realization using Hammarling technique without calculating augmented system realization Gramian matricesPbi and Qbo.


## 3.1. Proposed Technique

Let a new controllabilityPIG and observability QIG Gramians, respectively, be calculated by solving the following Lyapunov equations:
(18)APIGAT-PIG+BIGBIGT=0,(19)ATQIGA-QIG+CIGTCIG=0.The matricesBIG and CIG are the new fictitious input and output matrices, respectively, and are defined as
(20)BIG={UIG(SIG-snI)1/2forsn<0UIGSIG1/2forsn≥0,CIG={(RIG-rnI)1/2VIGTforrn<0RIG1/2VIGTforrn≥0.The termsUIG,  SIG,  VIG, and RIG are calculated by the orthogonal eigen decomposition of symmetric matrices Xen=UIGSIGUIGT and Yen=VIGRIGVIGT, where SIG=diag⁡(s1,s2,s3,…,sn), RIG=diag⁡(r1,r2,r3,…,rn), s1≥s2≥⋯≥sn, and r1≥r2≥⋯≥rn.Let a contragradient transformation matrixT be obtained as
(21)TTQIGT=T-1PIGT-T=Σ,
where Σ=diag⁡{σ1,σ2,σ3,…,σn}, and σj≥σj+1,  j=1,2,…,n-1, σl>σl+1. The reduced order system is calculated by transforming and partitioning the original system.Remark 4.
SinceXen≤BIGBIGT≥0, Yen≤CIGTCIG≥0, PIG>0 and QIG>0, therefore, the realization (A,BIGS,CIGS) is minimal. Moreover, the reduced order models are stable.Theorem 5.
The following error bound for the proposed technique holds if the rank conditionsrank⁡[BIGB]=rank⁡[BIG] and rank⁡[CIGC]=rank⁡[CIG] (which follows from [2, 28, 29, 38]) are satisfied
(22)(i)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LIG∥∞∥KIGVi(z)∥∞∑j=l+1n‍σj(ii)∥(H(z)-Hl(z))Vi(z)∥∞≤2∥KIGVi(z)∥∞∑j=l+1n‍σj(iii)∥Wo(z)(H(z)-Hl(z))∥∞≤2∥Wo(z)LIG∥∞∑j=l+1n‍σj,
where
(23)LIG={CV(RIG-rnI)-1/2forrn<0CVRIG-1/2forrn≥0,KIG={(SIG-snI)-1/2UTBforsn<0SIG-1/2UTBforsn≥0.Proof.
We show proof of (i) (whereas (ii) and (iii) are special cases of (i)). Since the rank conditionsrank⁡[BIGB]=rank⁡[BIG] and rank⁡[CIGC]=rank⁡[CIG] are satisfied, the relationships B=BIGKIG and C=LIGCIG hold. By partitioning BIG=[BIG1BIG2],  CIG=[CIG1CIG2] and then substituting B1=BIG1KIG, C1=LIGCIG1, respectively, yields
(24)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞=∥Wo(z)(C(zI-A)-1B-C1(zI-A11)-1B1)Vi(z)∥∞=∥Wo(z)(-LIGCIG1(zI-A11)-1BIG1KIGLIGCIG(zI-A)-1BIGKIGmmmmmmn-LIGCIG1(zI-A11)-1BIG1KIG)Vi(z)∥∞=∥Wo(z)LIG(-CIG1(zI-A11)-1BIG1CIG(zI-A)-1BIGmmmmmmmmn-CIG1(zI-A11)-1BIG1)KIGVi(z)∥∞=∥Wo(z)LIG∥∞×∥(CIG(zI-A)-1BIG-CIG1(zI-A11)-1BIG1)∥∞×∥KIGVi(z)∥∞.
If {A11,BIG1,CIG1} is reduced order model obtained by partitioning a balanced realization {A,BIG,CIG}, we have [7, 39]
(25)∥(CIG(zI-A)-1BIG-CIG1(zI-A11)-1BIG1)∥∞≤2∑j=l+1n‍σj.
Therefore,
(26)∥Wo(z)(H(z)-Hl(z))Vi(z)∥∞≤2∥Wo(z)LIG∥∞∥KIGVi(z)∥∞∑j=l+1n‍σj.Remark 6.
For the case when symmetric matricesXen≥0 and Yen≥0, then Pen=PCS=PVd=PIG and Qen=QCS=QVd=QIG. However, when matrices Xen and Yen are indefinite, then Pen<PIG and Qen<QIG. Moreover, frequency weighted Hankel singular values satisfy (λj[PenQen])1/2≤(λj[PIGQIG])1/2.Remark 7.
For the case when inputVi(z) weights are co-inner and output Wo(z) weights are inner [40], then P=Pen=PCS=PVd=PIG and Q=Qen=QCS=QVd=QIG, where P and Q are unweighted Gramians defined as
(27)APAT-P+BBT=0,ATQA-Q+CTC=0.Remark 8.
For the case when the symmetric matricesXen≥0 and Yen≥0, the reduced order models obtained using Enns [25], Campbell et al.'s [37], Varga and Anderson's [38], and proposed technique are same.


## 3.2. Computational Aspects

The frequency weighted balanced truncation model reduction techniques balance the original system and then truncate the balanced realization to get the desired reduced order system. The balancing procedure involves computation of transformation matrix from controllability and observability Gramians. Sometimes these matrices become numerically low rank especially in large scale systems (possibly) due to rapid decay of their eigenvalues [3]. Due to this reason, balancing procedure becomes inefficient. Accuracy enhancing techniques for different frequency weighted model reduction techniques appear in [4, 29].For unweighted case, Hammarling's technique [41] is used to obtained Cholesky factors of Gramian matrices from original system realization without actually computing controllability and observability Gramian matrices, respectively.In frequency weighted techniques, Cholesky factors of the Gramian matrices are obtained from the augmented system realizations. LetS¯ and R¯ be the Cholesky factors of the augmented system Gramians matrices Pbi and Qbo of (4) and (5), respectively,
(28)Pbi=S¯S¯T=[S11S120S22][S11T0S12TS22T]=[S11S11T+S12S12TS12S22TS22S12TS22S22T]=[PenP12P12TPV],Qbo=R¯TR¯=[R11T0R12TR22T][R11R120R22]=[R11TR11R11TR12R12TR11R22TR22+R12TR12]=[QWQ12TQ12Qen].By making use of Cholesky factorsS¯ and R¯ calculated above, the Cholesky factors corresponding to Gramians in frequency weighted model reduction techniques like Enns [25], Campbell et al.'s [37], Varga and Anderson's [38], and proposed technique can be obtained as follows:(
1) Enns Technique. The Cholesky factors Sen=[S11S12] and Ren=[R12R22] satisfy [29]
(29)Pen=SenSenT=S11S11T+S12S12T=[S11S12][S11TS22T],Qen=RenTRen=R22TR22+R12TR12=[R22TR12T][R22R12].(
2
) Campbell et al.'s Technique. The Cholesky factors S¯CS and R¯CS satisfy PCS=S¯CSS¯CST and QCS=R¯CSTR¯CS [4].(
3
) Varga and Anderson's Technique. The Cholesky factors S¯Vd and R¯Vd satisfy PVd=S¯VdS¯VdT and QVd=R¯VdTR¯Vd [4].(
4
) Proposed Technique. The Cholesky factors S¯IG and R¯IG satisfy PIG=S¯IGS¯IGT and QIG=R¯IGTR¯IG.In the following we establish a relationship between Cholesky factors of Gramian matrices used in Enns and proposed techniques. Equations (18) and (19) can be expressed as
(30)A(Pen+Pad)AT-(Pen+Pad)+(Xen-snI)=0,forsn<0,APenAT-Pen+Xen=0,forsn≥0,AT(Qen+Qad)A-(Qen+Qad)+(Yen-rnI)=0,forrn<0,ATQenA-Qen+Yen=0,forrn≥0,APadAT-Pad-snI=0,forsn<0,ATQadA-Qad-rnI=0,forrn<0.
Since
(31)XIG=U(S-snI)1/2(S-snI)1/2UT=Xen-snI,forsn<0,XIG=U(S)1/2(S)1/2UT=Xen,forsn≥0,YIG=VT(R-rnI)1/2(R-rnI)1/2V=Yen-rnI,forrn<0YIG=VT(R)1/2(R)1/2V=Yen,forrn≥0.
By using Hammarling technique to calculate Cholesky factors of Gramians Pad and Qad from realization {A,-snI,-rnI,D}, we can write Pad=S¯adS¯adT and Qad=R¯adTR¯ad. Therefore, frequency weighted controllability PIG (18) and observability QIG (19) Gramians can be expressed as
(32)PIG=S¯IGS¯IGT=Pen+Pad=S11S11T+S12S12T+SadSadT=[S11S12Sad][S11TS12TSadT],(33)QIG=R¯IGTR¯IG=Qen+Qad=R22TR22+R12TR12+RadTRad=[R22TR12TRadT][R22R12Rad].Remark 9.
Note that, Cholesky factors for Enns and proposed technique are computed directly from augmented system realization using Hammarling technique without calculating augmented system realization Gramian matricesPbi and Qbo.


## 4. Illustrative Examples

In this section, using numerical illustrative examples we show the usefulness of the proposed technique in comparison with existing frequency weighted balanced model reduction techniques for discrete time systems. Note that, proposed work deals with frequency weighted model reduction problem for discrete time systems, therefore, comparison is done with existing frequency weighted balanced model reduction techniques only.Example 1.
Consider (example C appeared in [26]) a 4th order stable discrete time system
(34)H(z)=z3z4+1.1z3-0.01z2-0.275z-0.06
with the following weightings
(35)Vi(z)=Wo(z)=z+0.9z+0.1.
The first order reduced model obtained by Enns [25] technique is unstable while reduced order model obtained by Campbell et al.'s, Varga and Anderson's,and proposed techniques is stable yielding frequency response errors 112.9338, 100.8739, and 94.116, respectively. Note that, proposed technique provides stability with relatively lower error when compared to other techniques.Example 2.
Consider a 6th order stable low pass digital elliptic filter with0.2 dB of peak-to-peak ripple and a minimum stopband attenuation of 20 dB represented by
(36)H(z)=(0.1054z6-0.1944z5+0.1187z4-0.1187z2+0.1944z-0.1054)×(z6-2.9621z5+4.8325z4-4.9819z3mmi+3.5245z2-1.5262z+0.3657)-1
with the following input and output weightings, respectively,
(37)Vi(z)=z3+3.0081z2+1.9944z+1.0325z3+0.2z2+0.75z+0.2Wo(z)=z3+2.97z2+2.9403z+0.9703z3+1.1619z2+0.6959z+0.1378.Table1 gives the comparison of error and error bounds for reduced order systems obtained by Enns, Campbell et al.'s, Varga and Anderson's, and proposed techniques for the input and two sided weighting cases. Note that, the proposed technique mostly yields lower error as compared to other techniques.Table 1
Frequency weighted errors and error bounds comparison for reduced order models.
Weighting
Order
Enns technique [25]
Campbell et al.’s technique [37]
Varga and Anderson’s technique [38]
Proposed technique
Error
Error bound
Error
Error bound
Error
Error bound
Two sided
1
21.1254
20.5953
1634.2
21.6491
725.7718
15.6065
2549.0
2
31.9647
32.8319
978.34
32.8863
433.8123
18.4571
1623.9
3
35.0441
32.3860
590.48
33.9063
255.7017
26.1274
998.92
4
28.7611
31.4710
236.41
30.4331
102.7518
30.7929
427.13
5
12.7538
29.5760
117.81
12.7660
50.4647
25.6547
203.36
Input
1
7.0257
7.1275
145.811
7.2356
92.6748
7.0140
242.5983
2
10.4643
10.7354
87.9789
10.7694
55.8603
10.6714
149.7017
3
11.2055
10.3852
53.0816
10.8477
32.9346
9.9857
90.8656
4
8.9654
10.0342
21.8067
9.6079
13.5877
8.4277
41.1182
5
2.4435
3.1720
10.4718
2.8446
6.3761
3.2061
16.5445Example 3.
Consider a 4th order stable discrete time system [42]
(38)H(z)=10-3(3.315z3-4.9695z2+2.1668z-0.24002)z4-3.7035z3+5.1957z2-3.2718z+0.77986
with the following input weighting:
(39)Vi(z)=z2-0.1z-0.05z2-0.9z+0.75.
Table2 gives the comparison of error and error bounds for reduced order systems obtained by Enns, Campbell et al.'s Varga and Anderson's and proposed techniques for the input weighting case. Note that, the proposed technique compares well and yields relatively lower error as compared to other techniques.Table 2
Frequency weighted errors and error bounds comparison for reduced order models.
Weighting
Order
Enns technique [25]
Campbell et al.’s technique [37]
Varga and Anderson’s technique [38]
Proposed technique
Error
Error bound
Error
Error bound
Error
Error bound
Input
1
0.0216
0.0245
0.6024
0.0241
0.3321
0.0240
0.6783
2
0.0021
0.0027
0.2582
0.0026
0.1403
0.0025
0.2853
3
0.0015
0.0025
0.0453
0.0023
0.0247
0.0021
0.0508Example 4.
Consider a 4th order stable low pass digital Chebychev type 1 filter with 0.1 dB of peak-to-peak ripples in the passband represented by(40)H(z)=0.49z4-0.9799z2+0.49z4-0.2893z3-0.6629z2+0.0246z+0.2904
with the following output weighting:
(41)Wo(z)=z-0.2z2-0.4z+0.5.Table3 gives the comparison of error and error bounds for reduced order systems obtained by Enns, Campbell et al.'s Varga and Anderson's and proposed techniques for output weighting. Note that, the proposed technique compares well and yields relatively lower error as compared to other techniques.Table 3
Frequency weighted errors and error bounds comparison for reduced order models.
Weighting
Order
Enns technique [25]
Campbell et al.’s technique [37]
Varga and Anderson’s technique [38]
Proposed technique
Error
Error bound
Error
Error bound
Error
Error bound
Output
1
1.7905
1.7898
8.2727
1.7900
5.6369
1.7867
6.2509
2
0.8967
0.8126
2.9832
0.8520
1.9243
0.7296
2.5083
3
0.5098
0.4868
1.4488
0.4979
0.9219
0.4740
0.9543


## 5. Conclusion

A new frequency weighted technique for model reduction of discrete time systems is explored. The reduced order models obtained in the presence of input, output, and two sided weightings are stable. A comparison with existing schemes shows that proposed technique provides comparable results (mostly produces lower error) for reduced order models.



---

*Source: 498453-2014-03-18.xml*