# Efficient Techniques for Solving the Periodic Projected Lyapunov Equations and Model Reduction of Periodic Systems

Authors: Mohammad-Sahadet Hossain; M. Monir Uddin
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We have presented the efficient techniques for the solutions of large-scale sparse projected periodic discrete-time Lyapunov equations in lifted form. These types of problems arise in model reduction and state feedback problems of periodic descriptor systems. Two most popular techniques to solve such Lyapunov equations iteratively are the low-rank alternating direction implicit (LR-ADI) method and the low-rank Smith method. The main contribution of this paper is to update the LR-ADI method by exploiting the ideas of the adaptive shift parameters computation and the efficient handling of complex shift parameters. These approaches efficiently reduce the computational cost with respect to time and memory. We also apply these iterative Lyapunov solvers in balanced truncation model reduction of periodic discrete-time descriptor systems. We illustrate numerical results to show the performance and accuracy of the proposed methods.


---

## Body

## 1. Introduction

Periodic systems and control theory have received a lot of attention in many areas of science and engineering, especially in the areas where the periodic control is deserved, such as spacecraft attitude control, control of industrial processes and mechanical systems, helicopter design, and nonlinear circuit design [1–6]. In this paper, we consider the linear periodic discrete-time descriptor systems with time-varying dimensions in the form(1)Ekxk+1=Akxk+Bkuk,yk=Ckxk,k∈Z,where Ek∈Rμk+1×nk+1, Ak∈Rμk+1×nk, Bk∈Rμk+1×mk, Ck∈Rqk×nk are periodic with a period K≥1, ∑k=0K-1μk=∑k=0K-1‍nk=n, ∑k=0K-1‍mk=m, and ∑k=0K-1‍qk=q. If all the Ek matrices are nonsingular, then (1) can be transformed into a periodic standard system.Model reduction means to replace the given mathematical model (1) by a much smaller model which preserves certain crucial properties of the original model, such as stability and passivity, and it is then amenable for simulation and analysis. Analysis and reduced-order modeling of system (1) may require solving large-scale sparse projected periodic discrete-time Lyapunov equations. Model reduction procedures of such periodic systems using direct solution approaches have been considered in [6–10]. Recently, attention has been devoted to iterative solutions of large-scale sparse Lyapunov equations using the alternating directions implicit (ADI) method [11, 12], the Smith method [12–14], and Krylov subspace methods [15, 16]. Generalization of all these methods has also been considered in [17, 18].However, low-rank iterative methods and their application in model reduction of such periodic descriptor systems have been considered in [19, 20]. The drawbacks of the ADI method proposed that there are the suboptimal shifts computations and the slow convergence of the iterative solutions to the exact solutions. In this paper, we proposed the adaptive shift parameters [21] computation instead of suboptimal parameters and explored the efficient handing of the complex shift [22] parameters. These two proposed schemes reduce the computational cost and fasten the convergence of the iterative solutions to exact solutions. Numerical results are discussed to show the efficiency of the techniques.The paper outline is as follows. In Section2, we briefly review the analogous time-invariant representation of the periodic systems, known aslifted representation of periodic systems. We discuss the cyclic lifted representation of periodic descriptor system (1). We reviewed the causal-noncausal decomposition, the periodic projected Lyapunov equations [8], and their corresponding lifted forms [7]. We also review some important concepts of the discrete-time descriptor system in linear time-invariant (LTI) form and the balanced truncation (BT) model reduction in LTI case. We also discuss the basics of the generalized ADI method in this section. The main idea behind this is that the concepts in the LTI structures will help more precisely to understand their corresponding periodic interpretations. In Section 3, we discuss the efficient solution of the causal lifted Lyapunov equations using the adaptive shifts based LR-ADI method. We also reviewed the Smith method from [20] for noncausal lifted Lyapunov equations in Section 4. We consider a balanced truncation model reduction method for periodic descriptor systems in Section 5. Section 6 contains numerical results that illustrate the efficiencies of the described iterative methods in model reduction of periodic descriptor system. Section 7 deficits the conclusion and remarks.


## 2. Preliminaries

In this section, we review the analogous time-invariant representation of the periodic systems, known aslifted representation. We also introduce some notations and briefly review the fundamental concept of balanced truncation model reduction for LTI discrete-time systems. We, afterward, generalize these concepts and results in the periodic setting.


### 2.1. Cyclic Lifted Representation of Periodic Systems

Analysis and modeling of periodic discrete-systems are often described by an analogous time-invariant representation of the periodic systems, known aslifted representation [1, 2, 6]. Here we consider thecyclic lifted representation [4] of periodic descriptor system (1), given by(2)EXk+1=AXk+BUk,Yk=CXk,where(3)E=diagE0,E1,…,EK-1,B=diagB0,B1,…,BK-1,A=0⋯0A0A10⋱⋮0AK-10,C=0⋯0C0C10⋱⋮0CK-10.The descriptor vector, system input, and output of (2) are related to those of (1) via (4)Xk=x1⋮xK-1x0,Uk=u0u1⋮uK-1,Yk=y0y1⋮yK-1,respectively. We can define the important concepts and dynamics of the periodic systems, such as regularity and asymptotic stability, by using their corresponding lifted forms.


### 2.2. The PPDLEs and Their Lifted Forms

Consider periodic discrete-time descriptor system (1), where the set of periodic matrix pairs {(Ek,Ak)}k=0K-1 is periodic stable [20, 23]. Then the periodic Gramians of discrete-time periodic descriptor systems satisfy a class of projected discrete-time Lyapunov equations [7, 8].Theorem 1 (see [23]).
(1) The periodic projected discrete-time Lyapunov equations with special right hand sides(5)AkGkcrAkT-EkGk+1crEkT=-PlkBkBkTPlkT,Gkcr=PrkGkcrPrkT,(6)AkGkncrAkT-EkGk+1ncrEkT=QlkBkBkTQlkT,Gkncr=QrkGkncrQrkT,with k=0,1,…,K-1, have the unique symmetric, positive semidefinite periodic solutions Gkcr and Gkncr, respectively.
(2) The periodic projected discrete-time Lyapunov equations with special right hand sides(7)AkTGk+1coAk-Ek-1TGkcoEk-1=-PrkTCkTCkPrk,Gkco=Plk-1TGkcoPlk-1,(8)AkTGk+1ncoAk-Ek-1TGkncoEk-1=QrkTCkTCkQrk,Gknco=Qlk-1TGkncoQlk-1,with k=0,1,…,K-1, have the unique symmetric, positive semidefinite periodic solutions Gkco and Gknco, respectively.The solutions{Gkcr}k=0K-1 of (5) and {Gkncr}k=0K-1 of (6), respectively, are calledthe causal and noncausal reachability Gramians of system (1). Similarly, the solutions {Gkco⁡}k=0K-1 and {Gknco}k=0K-1 of (7) and (8), respectively, are calledthe causal and noncausal observability Gramians of (1).Theorem 2 (see [7]).
Consider a periodic descriptor system (1) and the corresponding lifted form (2), where {(Ek,Ak)}k=0K-1 is periodic stable. Then the matrices Gcr=diag(G1cr,…,GK-1cr,G0cr) and Gncr=diag(G1ncr,…,GK-1ncr) solve the lifted projected discrete-time algebraic Lyapunov equations (LPDALEs):(9a)AGcrAT-EGcrET=-PlBBTPlT,Gcr=PrGcrPrT,(9b)AGncrAT-EGncrET=QlBBTQlT,Gncr=QrGncrQrT,respectively, where E, A, and B are as in (3), and(10)Pl=diagPl0,Pl1,…,PlK-1,Ql=I-Pl,Pr=diagPr1,…,PrK-1,Pr0,Qr=I-Pr.The matricesGcr and Gncr are the solutions of the lifted causal and noncausal reachability Gramians, respectively. A similar result can also be stated for the lifted causal and noncausal observability Gramians [7]. The solution techniques for such Lyapunov equations are shown in [20] using low-rank ADI methods. In the following section we mainly focus on the solution techniques of (9a) more efficiently.


### 2.3. Balanced Truncation for LTI Discrete-Time Systems

In order to understand the BT model reduction method for the time-varying discrete-time systems, we briefly review the idea for the time-invariant case. Consider a LTI discrete-time system:(11)E¯x¯k+1=A¯x¯kt+B¯uk,y¯k=C¯x¯k,in which E¯∈Rnc×nc is nonsingular, and A¯∈Rnc×nc, B¯∈Rnc×mc, and C¯∈Rqc×nc. Assume that the original system is asymptotically stable; that is, all finite eigenvalues of the system pencil lie inside the unite circle; then the controllability Gramian P∈Rnc×nc and the observability Gramian Q∈Rnc×nc can be computed by solving the discrete-time Lyapunov equations (see, e.g., [24]):(12)A¯PA¯T-E¯PE¯T=-B¯B¯T,(13)A¯TQA¯-E¯TQE¯=-C¯TC¯.Equivalently, to compute the controllability and observability Gramians we can solve the continuous-time Lyapunov equations(14)A_PE_T+E_PA_T=-B_B_T,(15)A_TQE_+E_TQA_=-C_TC_,where λA_-E_≕λ(A¯-E¯)-(A¯+E¯) is the Cayley-transformed pencil [25], and B_=2B¯, and C_=2C¯. Note that the solutions of (12) and (13) are identical to the solutions of (14) and (15), respectively [20, 25]. It is true that the Cayley transformation can destroy the system’s sparsity pattern where it needs to be restored. But, using proper permutations and efficient computation, one can recover this difficulty.Since both the Gramians are symmetric and positive semidefinite, they have symmetric Gramians factors(16)P=RRT,Q=LLT.The main drawback of BT based model order reduction (MOR) techniques is the computation of the solutions of the Lyapunov equations that are used in computing the Gramian factors. During the recent decades several approaches have been developed that allow exploiting the fact that often all coefficient matrices are highly sparse and the numbers of inputs and outputs are very small compared to the number of DoFs. Then usually the GramiansP and Q can be expressed in low-rank factored forms; that is, R and L are thin rectangular matrices. Therefore instead of computing the full Gramian factors it is quite reasonable to compute the thin rectangular matrices R and L.


### 2.4. Low-Rank ADI for LTI Discrete-Time Systems

One popular iterative method for computing the low-rankR and L is the generalized low-rank Cholesky factor-alternating directions implicit (GLR-ADI) iteration [11, 26, 27]. In the following we briefly introduce the GLR-ADI methods to solve continuous-time algebraic Lyapunov equations (14) and (15). For model reduction of system (1) we have to solve such a kind of Lyapunov equations.The GLR-ADI iterations successively compute the columns of the two Gramian factors following the simple steps [28]:(17)Z1=V1=2Reμ1F+μ1N-1G,Vi=Vi-1-μi+μ¯i-1F+μiN-1NVi-1,Zi=Zi-12ReμiVi,i=2:imax.In this iteration either(N,F,G)=(E_,A_,B_) is applied to generates Z=RT, or (N,F,G)=(E_T,A_T,C_T) is used for Z=LT.A set of optimalADI shift parameters or simplyshift parameters{μi}i=1J⊂C- are necessary for fast convergence of the GLR-ADI iterations. We will discuss the shift parameter selection criterion later in this section. In these iterations we also see that if the maximum number of iterations imax is greater than the number of shifts J, then the shift parameters are used in a cyclic way.In the GLR-ADI iterations all of the input matricesN, F, and G are real, due to the complex shift parameters in each iteration step, the updated Zi store complex data, which increases the overall complexity and memory requirements of the method. Moreover, in the balancing based methods using these complex Gramian factors yields complex reduced systems by performing some complex arithmetic operations. This problem is resolved in [22]. In this regard, the important assumption is that the selected ADI shift parameters should be proper.Definition 3.
The ADI shift parameters{μi}i=1J⊂C- are called proper if μi and μi+1 are complex conjugates of each other when one of them is complex.In [22] it is shown that at the (i+1)-st iteration of the GLR-ADI iterations, Vi+1 can be computed by(18)Vi+1=V¯i+2δIm⁡Vi,where δ=Re(μi)/Im(μi). This identity states that, for μi+1=μ¯i, Vi+1 can be computed explicitly from Vi, which releases us from solving a shifted linear system with F+μ¯iN. This idea also results in the following theorem.Theorem 4.
Let us assume a set of proper ADI shift parameters. For a pair of complex conjugate shifts{μi,μi+1≔μ¯i}, the two subsequent block iterates Vi and Vi+1 of GLR-ADI iterations satisfy [22](19)Vi,Vi+1=-2Re⁡μiRe⁡Vi+δIm⁡Vi,-2Re⁡μiδ2+1Im⁡Vi.This theorem reveals that for a pair of complex conjugate shifts at any iteration step in the GLR-ADI iterationZi can be updated by(20)Zi+1=Zi-1,-2Re⁡μiRe⁡Vi+δIm⁡Vi,-2Re⁡μiδ2+1Im⁡Vi.A version of the GLR-ADI algorithm is summarized in [29, Algorithm 5], which computes low-rank real Gramian factors. Additionally, the GLR-ADI iterations can be stopped whenever the norm of theADI-residual(21)FZi=FiZiTNT+NZiZiTFT+GGTbecomes very small. But computing F(Zi) in Frobenius-norm or 2-norm in each iteration step is an expensive task, since in each iteration the resulting residual (21) is an nc×nc matrix. Recently, in [30] the authors show that, in the ith iteration, the ADI-residual in (21) can be represented as (22)FZi=WiWiH,with(23)Wi=∏j=1iF-μ¯jNF+μjN-1G,and the Vi in the GLR-ADI iterations can be expressed as [30](24)Vi=F+μiN-1Wi-1.The authors in [30] also show that the ADI-residual factor defined in (23) can be computed by the expression(25)Wi=Wi-1-2Re⁡μiNVi.In the case of real setting when we consider μi+1≔μ¯i, one must compute(26)Wi+1=Wi-1-4Re⁡μiNRe⁡Vi+δIm⁡Vi,where δ is defined in (18). The rank of Wi is at most mc, that is, the number of columns of G. Therefore, the computation of the Frobenius-norm or 2-norm of WiWiT=WiTWi in each iteration is extremely cheap. Applying these strategies (computation of real Gramian factors and low-rank residual based stopping techniques), the updated GLR-ADI will be briefly discussed in Section 3.


## 2.1. Cyclic Lifted Representation of Periodic Systems

Analysis and modeling of periodic discrete-systems are often described by an analogous time-invariant representation of the periodic systems, known aslifted representation [1, 2, 6]. Here we consider thecyclic lifted representation [4] of periodic descriptor system (1), given by(2)EXk+1=AXk+BUk,Yk=CXk,where(3)E=diagE0,E1,…,EK-1,B=diagB0,B1,…,BK-1,A=0⋯0A0A10⋱⋮0AK-10,C=0⋯0C0C10⋱⋮0CK-10.The descriptor vector, system input, and output of (2) are related to those of (1) via (4)Xk=x1⋮xK-1x0,Uk=u0u1⋮uK-1,Yk=y0y1⋮yK-1,respectively. We can define the important concepts and dynamics of the periodic systems, such as regularity and asymptotic stability, by using their corresponding lifted forms.


## 2.2. The PPDLEs and Their Lifted Forms

Consider periodic discrete-time descriptor system (1), where the set of periodic matrix pairs {(Ek,Ak)}k=0K-1 is periodic stable [20, 23]. Then the periodic Gramians of discrete-time periodic descriptor systems satisfy a class of projected discrete-time Lyapunov equations [7, 8].Theorem 1 (see [23]).
(1) The periodic projected discrete-time Lyapunov equations with special right hand sides(5)AkGkcrAkT-EkGk+1crEkT=-PlkBkBkTPlkT,Gkcr=PrkGkcrPrkT,(6)AkGkncrAkT-EkGk+1ncrEkT=QlkBkBkTQlkT,Gkncr=QrkGkncrQrkT,with k=0,1,…,K-1, have the unique symmetric, positive semidefinite periodic solutions Gkcr and Gkncr, respectively.
(2) The periodic projected discrete-time Lyapunov equations with special right hand sides(7)AkTGk+1coAk-Ek-1TGkcoEk-1=-PrkTCkTCkPrk,Gkco=Plk-1TGkcoPlk-1,(8)AkTGk+1ncoAk-Ek-1TGkncoEk-1=QrkTCkTCkQrk,Gknco=Qlk-1TGkncoQlk-1,with k=0,1,…,K-1, have the unique symmetric, positive semidefinite periodic solutions Gkco and Gknco, respectively.The solutions{Gkcr}k=0K-1 of (5) and {Gkncr}k=0K-1 of (6), respectively, are calledthe causal and noncausal reachability Gramians of system (1). Similarly, the solutions {Gkco⁡}k=0K-1 and {Gknco}k=0K-1 of (7) and (8), respectively, are calledthe causal and noncausal observability Gramians of (1).Theorem 2 (see [7]).
Consider a periodic descriptor system (1) and the corresponding lifted form (2), where {(Ek,Ak)}k=0K-1 is periodic stable. Then the matrices Gcr=diag(G1cr,…,GK-1cr,G0cr) and Gncr=diag(G1ncr,…,GK-1ncr) solve the lifted projected discrete-time algebraic Lyapunov equations (LPDALEs):(9a)AGcrAT-EGcrET=-PlBBTPlT,Gcr=PrGcrPrT,(9b)AGncrAT-EGncrET=QlBBTQlT,Gncr=QrGncrQrT,respectively, where E, A, and B are as in (3), and(10)Pl=diagPl0,Pl1,…,PlK-1,Ql=I-Pl,Pr=diagPr1,…,PrK-1,Pr0,Qr=I-Pr.The matricesGcr and Gncr are the solutions of the lifted causal and noncausal reachability Gramians, respectively. A similar result can also be stated for the lifted causal and noncausal observability Gramians [7]. The solution techniques for such Lyapunov equations are shown in [20] using low-rank ADI methods. In the following section we mainly focus on the solution techniques of (9a) more efficiently.


## 2.3. Balanced Truncation for LTI Discrete-Time Systems

In order to understand the BT model reduction method for the time-varying discrete-time systems, we briefly review the idea for the time-invariant case. Consider a LTI discrete-time system:(11)E¯x¯k+1=A¯x¯kt+B¯uk,y¯k=C¯x¯k,in which E¯∈Rnc×nc is nonsingular, and A¯∈Rnc×nc, B¯∈Rnc×mc, and C¯∈Rqc×nc. Assume that the original system is asymptotically stable; that is, all finite eigenvalues of the system pencil lie inside the unite circle; then the controllability Gramian P∈Rnc×nc and the observability Gramian Q∈Rnc×nc can be computed by solving the discrete-time Lyapunov equations (see, e.g., [24]):(12)A¯PA¯T-E¯PE¯T=-B¯B¯T,(13)A¯TQA¯-E¯TQE¯=-C¯TC¯.Equivalently, to compute the controllability and observability Gramians we can solve the continuous-time Lyapunov equations(14)A_PE_T+E_PA_T=-B_B_T,(15)A_TQE_+E_TQA_=-C_TC_,where λA_-E_≕λ(A¯-E¯)-(A¯+E¯) is the Cayley-transformed pencil [25], and B_=2B¯, and C_=2C¯. Note that the solutions of (12) and (13) are identical to the solutions of (14) and (15), respectively [20, 25]. It is true that the Cayley transformation can destroy the system’s sparsity pattern where it needs to be restored. But, using proper permutations and efficient computation, one can recover this difficulty.Since both the Gramians are symmetric and positive semidefinite, they have symmetric Gramians factors(16)P=RRT,Q=LLT.The main drawback of BT based model order reduction (MOR) techniques is the computation of the solutions of the Lyapunov equations that are used in computing the Gramian factors. During the recent decades several approaches have been developed that allow exploiting the fact that often all coefficient matrices are highly sparse and the numbers of inputs and outputs are very small compared to the number of DoFs. Then usually the GramiansP and Q can be expressed in low-rank factored forms; that is, R and L are thin rectangular matrices. Therefore instead of computing the full Gramian factors it is quite reasonable to compute the thin rectangular matrices R and L.


## 2.4. Low-Rank ADI for LTI Discrete-Time Systems

One popular iterative method for computing the low-rankR and L is the generalized low-rank Cholesky factor-alternating directions implicit (GLR-ADI) iteration [11, 26, 27]. In the following we briefly introduce the GLR-ADI methods to solve continuous-time algebraic Lyapunov equations (14) and (15). For model reduction of system (1) we have to solve such a kind of Lyapunov equations.The GLR-ADI iterations successively compute the columns of the two Gramian factors following the simple steps [28]:(17)Z1=V1=2Reμ1F+μ1N-1G,Vi=Vi-1-μi+μ¯i-1F+μiN-1NVi-1,Zi=Zi-12ReμiVi,i=2:imax.In this iteration either(N,F,G)=(E_,A_,B_) is applied to generates Z=RT, or (N,F,G)=(E_T,A_T,C_T) is used for Z=LT.A set of optimalADI shift parameters or simplyshift parameters{μi}i=1J⊂C- are necessary for fast convergence of the GLR-ADI iterations. We will discuss the shift parameter selection criterion later in this section. In these iterations we also see that if the maximum number of iterations imax is greater than the number of shifts J, then the shift parameters are used in a cyclic way.In the GLR-ADI iterations all of the input matricesN, F, and G are real, due to the complex shift parameters in each iteration step, the updated Zi store complex data, which increases the overall complexity and memory requirements of the method. Moreover, in the balancing based methods using these complex Gramian factors yields complex reduced systems by performing some complex arithmetic operations. This problem is resolved in [22]. In this regard, the important assumption is that the selected ADI shift parameters should be proper.Definition 3.
The ADI shift parameters{μi}i=1J⊂C- are called proper if μi and μi+1 are complex conjugates of each other when one of them is complex.In [22] it is shown that at the (i+1)-st iteration of the GLR-ADI iterations, Vi+1 can be computed by(18)Vi+1=V¯i+2δIm⁡Vi,where δ=Re(μi)/Im(μi). This identity states that, for μi+1=μ¯i, Vi+1 can be computed explicitly from Vi, which releases us from solving a shifted linear system with F+μ¯iN. This idea also results in the following theorem.Theorem 4.
Let us assume a set of proper ADI shift parameters. For a pair of complex conjugate shifts{μi,μi+1≔μ¯i}, the two subsequent block iterates Vi and Vi+1 of GLR-ADI iterations satisfy [22](19)Vi,Vi+1=-2Re⁡μiRe⁡Vi+δIm⁡Vi,-2Re⁡μiδ2+1Im⁡Vi.This theorem reveals that for a pair of complex conjugate shifts at any iteration step in the GLR-ADI iterationZi can be updated by(20)Zi+1=Zi-1,-2Re⁡μiRe⁡Vi+δIm⁡Vi,-2Re⁡μiδ2+1Im⁡Vi.A version of the GLR-ADI algorithm is summarized in [29, Algorithm 5], which computes low-rank real Gramian factors. Additionally, the GLR-ADI iterations can be stopped whenever the norm of theADI-residual(21)FZi=FiZiTNT+NZiZiTFT+GGTbecomes very small. But computing F(Zi) in Frobenius-norm or 2-norm in each iteration step is an expensive task, since in each iteration the resulting residual (21) is an nc×nc matrix. Recently, in [30] the authors show that, in the ith iteration, the ADI-residual in (21) can be represented as (22)FZi=WiWiH,with(23)Wi=∏j=1iF-μ¯jNF+μjN-1G,and the Vi in the GLR-ADI iterations can be expressed as [30](24)Vi=F+μiN-1Wi-1.The authors in [30] also show that the ADI-residual factor defined in (23) can be computed by the expression(25)Wi=Wi-1-2Re⁡μiNVi.In the case of real setting when we consider μi+1≔μ¯i, one must compute(26)Wi+1=Wi-1-4Re⁡μiNRe⁡Vi+δIm⁡Vi,where δ is defined in (18). The rank of Wi is at most mc, that is, the number of columns of G. Therefore, the computation of the Frobenius-norm or 2-norm of WiWiT=WiTWi in each iteration is extremely cheap. Applying these strategies (computation of real Gramian factors and low-rank residual based stopping techniques), the updated GLR-ADI will be briefly discussed in Section 3.


## 3. Efficient Solutions of Causal LPDALEs Using Low-Rank ADI Methods

The basic ADI method for linear systems has been proposed in [31] and then used in [11, 12, 17, 32–34] for solving (projected) continuous-time and discrete-time Lyapunov equations. Consider the LPDALE (9a), where the matrix E=diag(E0,…,EK-1) is singular. For such equations both the ADI and Smith iterations fail to converge [20]. This problem is circumvented by considering ageneralized Cayley transformation [25]:(27)CE,A=λA-E-A+E,which transforms the LPDALE (9a) to an equivalent projected continuous-time algebraic Lyapunov equation (PCALE):(28)EGcrAT+AGcrET=-2PlBBTPlT,Gcr=PrGcrPrT.Here λE-A=λ(A-E)-(A+E) is the Cayley-transformed pencil. In this transformation the finite stable eigenvalues of λE-A are mapped to the finite stable eigenvalues of λE-A, and the infinite eigenvalue of λE-A is mapped to λ=1. More details of that transformation can be found in [20, 25]. The solutions of PCALE (28) and (9a) are unique and they are solvable for every BBT [20].


### 3.1. Efficient Computation of Low-Rank Gramian Factors

We present the updated version of low-rank ADI methods for Lyapunov equation (28), which is also the solution of (9a). An updated algorithm for the low-rank ADI method of generalized state space systems can be found in [29]. Here the idea is generalized for the periodic descriptor system.Applying ADI method the solutionGcr of Lyapunov equation (28) can be obtained by performing the following iteration steps, beginning with G0cr=0 and i=1:(29a)A+μiEiGi-1/2cr=-Gi-1crAT-μiE-B~B~T,(29b)A+μiEμ-iGi-1/2cr=-Gi-1/2crAT-μiE-B~B~T,where μi∈C-, with i=1,2,…,imax, is a shift parameter and B~=2PlB. Steps (29a) and (29b) can be written in a single step(30)Gicr=AμiGi-1crAμi∗-2Re⁡μiB~μiB~μi∗,where Aμi=(A-μ-iE)(A+μiE)-1 and B~μi=(A+μiE)-1B~. The error bound in iteration i is given by (see, e.g., [35])(31)ΔGicr=Gcr-Gicr=AμiGcr-Gi-1crAμi∗.Since G0cr=0, this leads to(32)ΔGicr=Gcr-Gicr=∏l=1iAμlGcr∏l=1iAμl∗.In ith iteration the residual can be computed by(33)Resi=AGicrET+EGicrAT+B~B~T=AGicr-GcrET+EGicr-GcrAT=∏l=1iAμl-AGcrET+EGcrAT∏l=1iAμl∗=∏l=1iAμlB~B~T∏l=1iAμl∗=WiWi∗,where Wi=∏l=1iAμlB~.We observe that the computedGicr in (30) can be decomposed as(34)Gicr=ZiZi∗,where Zi∈Rn×kl and kl≪n. Since the number of inputs is very small compared to the dimension of the system, the right hand side of (9a) is a low-rank matrix. Therefore often this decomposition is possible, and Zi can be computed by following iteration steps:(35a)Z1=2Reμ1A+μ1E-1B~,(35b)Zi=AμiZi-12ReμiA+μiE-1B~,with i=2,3,…,imax, and imax is maximum number of iterations. In each iteration in (35b) the number of columns of Zi is increased by the number of columns Zi-1. Therefore after few steps the number of columns of Zi grows dramatically. In order to resolve this drawback a further optimization occurs, where the number of columns to modify Zi is kept constant by using the clever reordering of the ADI-parameters. This updated version of GLR-ADI [36] is given below:(36a)Z1=V1=2Reμ1A+μ1E-1B~,(36b)Vi=Vi-1-μi+μ¯i-1A+μiE-1EVi-1,(36c)Zi=Zi-12ReμiVi,i=2,3,…,imax.Unfortunately, the GLR-ADI iteration does not preserve the block diagonal structure at every iteration step in Algorithm1. But Gicr=ZiZiT has almost block diagonal form since it approximates the solution of the LPDALE (9a) given by Gcr=diag(G1cr,…,GK-1cr,G0cr), where Gkcr are the periodic solutions of (5). We then consider the block partitioning (37)Zi=R1iT,…,RK-1,iT,R0iTT,with Rki∈Rnk×ri, and the causal reachability Gramian Gkcr can be approximated by the matrices RkiRkiT for k=0,1,…,K-1.Algorithm 1:Updated GLR-ADI (UGLR-ADI) for LPDALEs.
Input:  A,E,B~. Output:  Z=Zi, such that Gicr≈ZiZiT. (1)W0=B~,  Z0=,  i=1; (2)Vi=A+μiE-1Wi-1; (3)if  Imμi=0  then (4)Zi=Zi-1-2μiVi; (5)Wi=Wi-1-2μiEVi.; (6)else (7)γ=-2Reμi,  δ=Reμi/Imμi; (8)Zi+1=Zi-12γReVi+δImVi2γδ2+1ImVi; (9)Wi+1=Wi-1+2γEReVi+δImVi; (10)i=i+1; (11)i=i+1;We compute low-rank Gramian factorZ, so that ZZT≈Gcr by following the steps (36a)–(36c). We summarize the resulting iteration using column compression in Algorithm 1. It should be noted that in Algorithm 1 the computation of the matrices E and A can be avoided by rewriting the GLR-ADI iteration in terms of the original matrices E and A [20]. We should also note that we do not compute the inverses of the matrices at step (2) of Algorithm 1 explicitly. We represent the matrices E and A with the original sparse matrices with E and A and use the sparse inversion technique in step (2) [20].The iteration can be stopped in theith step if the residual norm, that is, Resi, where Resi is defined in (33), is sufficiently small. If we can compute Wi as a low-rank factor then computing Res is very cheap. Here, we discuss how to compute the low-rank residual-factor from each iteration in the GLRCF-ADI algorithm. Following the procedure as shown in Section 2.4, iterate Vi in (36b) can also be written as(38)Vi=A+μiE-1A-μ¯i-1EVi-1=A+μiE-1A-μ¯i-1EA+μi-1E-1A-μ¯i-2EVi-2=A+μiE-1∏l=1i-1AμlB~=A+μiE-1Wi-1,using33.Also,(39)Wi=A-μ¯iEA+μ¯iE-1∏l=1i-1AμlB~=A-μ¯iE-1Vi=A-μ¯iEA+μ¯iEWi-1=I-2ReμiEE+μiA-1Wi-1=Wi-1-2ReμiEVi,Wi+1=Wi-2Reμ¯iEVi+1=Wi-1-2ReμiEVi-2Reμ¯iEReVi-jImVi+2δImVi=Wi-1-2ReμiEReVi-jImVi-2Reμ¯iEReVi-jImVi+2δImVi=Wi-1-4ReμiEReVi+δImVi.


### 3.2. Selection of the ADI Shift Parameters

The convergence rate of the low-rank ADI algorithms discussed above depends on a set of ADI shift parameters. The ADI shift parameters were originally introduced by Wachspress in [37] to solve Lyapunov equations using the ADI methods. The author shows that a set of optimal ADI shift parameters {μi}i=1J for system (2) can be computed by solving the so-calledADI min-max problem [37, 38]:(40)μ1,μ2,…,μi=arg⁡minμ1,μ2,…,μi∈C-∖-1⁡maxλl∈SpE-1A∖1⁡∏j=1iμ-j-λlμj+λl,where Sp(E-1A) is the spectrum of the matrix E-1A.For a large-scale system, determining the entire spectrum of(A,E) is almost impossible. Therefore, in the literature several techniques are proposed; see, for example, [26, 36, 39, 40] to solve the min-max problem (40) using a much smaller part of the spectrum. Currently, one such commonly used technique is Penzl’sheuristic approach introduced in [26], where k+ Ritz values (see, e.g., [41]) and k- (k+,k-≪n) reciprocal Ritz values with respect to E-1A and A-1E, respectively, are employed. A complete procedure of the heuristic approach can be found in [20]. However, for large-scale descriptor systems, computing the Ritz values with respect to E-1A is a challenging task since E is then not invertible.Another ADI shift selection criterion that we emphasize here is theadaptive approach which is introduced in [21]. The performances of this adaptive shift selection approach for a large-scale system are also investigated in [29, 42, 43]. There the authors show that the approach is superior to the heuristic approach, especially for the (generalized) descriptor systems. To initialize the shifts, we first compute the eigenvalues of the projected pencil (41)λUTEU-UTAU,λ∈C,where U is formed by the basis of range of any nc×r (r≪nc) random matrix. Then select some desired number of optimal shifts by solving the ADI min-max problem like in the heuristic procedure. Then whenever all the initial shifts have been used, the pencil is projected to the span of the current Vi (see, e.g., Algorithm  1 of [21]) and computes some optimal shifts from the eigenvalues of the projected pencil and uses them as the set of next shift parameters. This approach is repeated while the algorithm has not converged to the given tolerance.


## 3.1. Efficient Computation of Low-Rank Gramian Factors

We present the updated version of low-rank ADI methods for Lyapunov equation (28), which is also the solution of (9a). An updated algorithm for the low-rank ADI method of generalized state space systems can be found in [29]. Here the idea is generalized for the periodic descriptor system.Applying ADI method the solutionGcr of Lyapunov equation (28) can be obtained by performing the following iteration steps, beginning with G0cr=0 and i=1:(29a)A+μiEiGi-1/2cr=-Gi-1crAT-μiE-B~B~T,(29b)A+μiEμ-iGi-1/2cr=-Gi-1/2crAT-μiE-B~B~T,where μi∈C-, with i=1,2,…,imax, is a shift parameter and B~=2PlB. Steps (29a) and (29b) can be written in a single step(30)Gicr=AμiGi-1crAμi∗-2Re⁡μiB~μiB~μi∗,where Aμi=(A-μ-iE)(A+μiE)-1 and B~μi=(A+μiE)-1B~. The error bound in iteration i is given by (see, e.g., [35])(31)ΔGicr=Gcr-Gicr=AμiGcr-Gi-1crAμi∗.Since G0cr=0, this leads to(32)ΔGicr=Gcr-Gicr=∏l=1iAμlGcr∏l=1iAμl∗.In ith iteration the residual can be computed by(33)Resi=AGicrET+EGicrAT+B~B~T=AGicr-GcrET+EGicr-GcrAT=∏l=1iAμl-AGcrET+EGcrAT∏l=1iAμl∗=∏l=1iAμlB~B~T∏l=1iAμl∗=WiWi∗,where Wi=∏l=1iAμlB~.We observe that the computedGicr in (30) can be decomposed as(34)Gicr=ZiZi∗,where Zi∈Rn×kl and kl≪n. Since the number of inputs is very small compared to the dimension of the system, the right hand side of (9a) is a low-rank matrix. Therefore often this decomposition is possible, and Zi can be computed by following iteration steps:(35a)Z1=2Reμ1A+μ1E-1B~,(35b)Zi=AμiZi-12ReμiA+μiE-1B~,with i=2,3,…,imax, and imax is maximum number of iterations. In each iteration in (35b) the number of columns of Zi is increased by the number of columns Zi-1. Therefore after few steps the number of columns of Zi grows dramatically. In order to resolve this drawback a further optimization occurs, where the number of columns to modify Zi is kept constant by using the clever reordering of the ADI-parameters. This updated version of GLR-ADI [36] is given below:(36a)Z1=V1=2Reμ1A+μ1E-1B~,(36b)Vi=Vi-1-μi+μ¯i-1A+μiE-1EVi-1,(36c)Zi=Zi-12ReμiVi,i=2,3,…,imax.Unfortunately, the GLR-ADI iteration does not preserve the block diagonal structure at every iteration step in Algorithm1. But Gicr=ZiZiT has almost block diagonal form since it approximates the solution of the LPDALE (9a) given by Gcr=diag(G1cr,…,GK-1cr,G0cr), where Gkcr are the periodic solutions of (5). We then consider the block partitioning (37)Zi=R1iT,…,RK-1,iT,R0iTT,with Rki∈Rnk×ri, and the causal reachability Gramian Gkcr can be approximated by the matrices RkiRkiT for k=0,1,…,K-1.Algorithm 1:Updated GLR-ADI (UGLR-ADI) for LPDALEs.
Input:  A,E,B~. Output:  Z=Zi, such that Gicr≈ZiZiT. (1)W0=B~,  Z0=,  i=1; (2)Vi=A+μiE-1Wi-1; (3)if  Imμi=0  then (4)Zi=Zi-1-2μiVi; (5)Wi=Wi-1-2μiEVi.; (6)else (7)γ=-2Reμi,  δ=Reμi/Imμi; (8)Zi+1=Zi-12γReVi+δImVi2γδ2+1ImVi; (9)Wi+1=Wi-1+2γEReVi+δImVi; (10)i=i+1; (11)i=i+1;We compute low-rank Gramian factorZ, so that ZZT≈Gcr by following the steps (36a)–(36c). We summarize the resulting iteration using column compression in Algorithm 1. It should be noted that in Algorithm 1 the computation of the matrices E and A can be avoided by rewriting the GLR-ADI iteration in terms of the original matrices E and A [20]. We should also note that we do not compute the inverses of the matrices at step (2) of Algorithm 1 explicitly. We represent the matrices E and A with the original sparse matrices with E and A and use the sparse inversion technique in step (2) [20].The iteration can be stopped in theith step if the residual norm, that is, Resi, where Resi is defined in (33), is sufficiently small. If we can compute Wi as a low-rank factor then computing Res is very cheap. Here, we discuss how to compute the low-rank residual-factor from each iteration in the GLRCF-ADI algorithm. Following the procedure as shown in Section 2.4, iterate Vi in (36b) can also be written as(38)Vi=A+μiE-1A-μ¯i-1EVi-1=A+μiE-1A-μ¯i-1EA+μi-1E-1A-μ¯i-2EVi-2=A+μiE-1∏l=1i-1AμlB~=A+μiE-1Wi-1,using33.Also,(39)Wi=A-μ¯iEA+μ¯iE-1∏l=1i-1AμlB~=A-μ¯iE-1Vi=A-μ¯iEA+μ¯iEWi-1=I-2ReμiEE+μiA-1Wi-1=Wi-1-2ReμiEVi,Wi+1=Wi-2Reμ¯iEVi+1=Wi-1-2ReμiEVi-2Reμ¯iEReVi-jImVi+2δImVi=Wi-1-2ReμiEReVi-jImVi-2Reμ¯iEReVi-jImVi+2δImVi=Wi-1-4ReμiEReVi+δImVi.


## 3.2. Selection of the ADI Shift Parameters

The convergence rate of the low-rank ADI algorithms discussed above depends on a set of ADI shift parameters. The ADI shift parameters were originally introduced by Wachspress in [37] to solve Lyapunov equations using the ADI methods. The author shows that a set of optimal ADI shift parameters {μi}i=1J for system (2) can be computed by solving the so-calledADI min-max problem [37, 38]:(40)μ1,μ2,…,μi=arg⁡minμ1,μ2,…,μi∈C-∖-1⁡maxλl∈SpE-1A∖1⁡∏j=1iμ-j-λlμj+λl,where Sp(E-1A) is the spectrum of the matrix E-1A.For a large-scale system, determining the entire spectrum of(A,E) is almost impossible. Therefore, in the literature several techniques are proposed; see, for example, [26, 36, 39, 40] to solve the min-max problem (40) using a much smaller part of the spectrum. Currently, one such commonly used technique is Penzl’sheuristic approach introduced in [26], where k+ Ritz values (see, e.g., [41]) and k- (k+,k-≪n) reciprocal Ritz values with respect to E-1A and A-1E, respectively, are employed. A complete procedure of the heuristic approach can be found in [20]. However, for large-scale descriptor systems, computing the Ritz values with respect to E-1A is a challenging task since E is then not invertible.Another ADI shift selection criterion that we emphasize here is theadaptive approach which is introduced in [21]. The performances of this adaptive shift selection approach for a large-scale system are also investigated in [29, 42, 43]. There the authors show that the approach is superior to the heuristic approach, especially for the (generalized) descriptor systems. To initialize the shifts, we first compute the eigenvalues of the projected pencil (41)λUTEU-UTAU,λ∈C,where U is formed by the basis of range of any nc×r (r≪nc) random matrix. Then select some desired number of optimal shifts by solving the ADI min-max problem like in the heuristic procedure. Then whenever all the initial shifts have been used, the pencil is projected to the span of the current Vi (see, e.g., Algorithm  1 of [21]) and computes some optimal shifts from the eigenvalues of the projected pencil and uses them as the set of next shift parameters. This approach is repeated while the algorithm has not converged to the given tolerance.


## 4. Smith Method for Noncausal LPDALEs

Consider the LPDALE (9b). For nonsingular A, this equation is equivalent to the LPDALE:(42)Gncr-A-1EGncrA-1ET=QrA-1BBTA-TQrT,Gncr=QrGncrQrT.Note that ifη is the index of the system [20], then the approximate solution of the noncausal LPDALE (9b) can be obtained by the Smith iterations [14]:(43)Gηncr=∑j=0η-1A-1EjQrA-1BBTA-TQrTA-1ETj=Gncr.Therefore, the Cholesky factor Xη of the solution Gncr=XηXηT of (42) and also of the LPDALE (9b) takes the form(44)Xη=QrA-1B,A-1EQrA-1B,…,A-1Eη-1QrA-1B.Thus, the solution of (9b) can be obtained with few computations for systems of low index. The solution of the noncausal LPDALE (9b) is briefly discussed in [20].Similar to the causal case, we consider the block partitioning(45)XηXηT=Gncr=diagG1ncr,…,GK-1ncr,G0ncr,such that the noncausal reachability Gramians of system (1) can be computed in factored form Gkncr=X~kX~kT for k=0,1,…,K-1.Remark 5.
The causal and noncausal observability GramiansGkco⁡ and Gknco of system (1) and the corresponding the low-rank Cholesky factors Lk and Y~k, where Gkco⁡≈LkLkT and Gknco=Y~kY~kT, can also be determined from the corresponding LPDALEs that are dual to the LPDALEs (9a) and (9b); see [7] for details.


## 5. Model Reduction of Periodic Systems

For periodic descriptor system (1), we construct a reduced-order model of dimension r of the form(46)E~kx~k+1=A~kx~k+B~kuk,y~k=C~kx~k,k∈Z,where E~k∈Rγk+1×rk+1, A~k∈Rγk+1×rk, B~k∈Rγk+1×mk, C~k∈Rqk×rk are K-periodic matrices, ∑k=0K-1‍γk=∑k=0K-1‍rk=r, and r≪n. Apart from having a much smaller state space dimension, the reduced-order model also preserves important physical properties of the original system such as regularity and stability and the small approximation error. Balanced truncation for discrete-time periodic descriptor system has been considered in [7, 20]. Some important concepts related to balanced truncation of discrete-time periodic descriptor systems has been presented in [8].Consider{(Ek,Ak)}k=0K-1 is periodic stable, and the Cholesky factors of the causal and noncausal Gramians satisfy(47)Gkcr=RkRkT,Gkco⁡=LkLkT,Gkncr=R^kR^kT,Gknco=L^kL^kT.Then thecausal andnoncausal Hankel singular valuesσk,j and θk,j of periodic descriptor system (1) are defined as(48)σk,j=ζjLkTEk-1Rk,θk,j=ζjL^k+1TAkR^k,respectively, where ζj(·) denotes the singular values of the corresponding product matrices. For a balanced system, one truncates the states related to the small causal Hankel singular values because they do not change system properties essentially. Unfortunately, truncation of small nonzero noncausal Hankel singular values may lead to unstable reduction [7].Hence, we compute the singular value decompositions of the product matrices(49)LkTEk-1Rk=Uk,1,Uk,2Σk,1Σk,2Vk,1,Vk,2T,Y~k+1TAkX~k=Uk,3ΘkVk,3T,where [Uk,1,Uk,2], [Vk,1,Vk,2], Uk,3, and Vk,3 are orthogonal:(50)Σk,1=diagσk,1,…,σk,rkf,Σk,2=diagσk,rkf+1,…,σk,nkf,with σk,1≥⋯≥σk,rkf>σk,rkf+1≥⋯≥σk,nkf, and Θk=diag(θk,1,…,θk,rk∞) is nonsingular for k=0,1,…,K-1.We then compute reduced-order system (46) as(51)E~k=Sk,rTEkTk+1,r,A~k=Sk,rTAkTk,r,B~k=Sk,rTBk,C~k=CkTk,r,where the projection matrices Sk,r and Tk,r have the form(52)Sk,r=Lk+1Uk+1,1Σk+1,1-1/2,Y~k+1Uk,3Θk-1/2∈Rμk+1×γk+1,Tk,r=RkVk,1Σk,1-1/2,X~kVk,3Θk-1/2∈Rnk×rk,with rk=rkf+rk∞ and γk+1=rk+1f+rk∞.LetH(z)=C(zE-A)-1B be the transfer function of lifted system (2) and H~(z) be the transfer function of the corresponding reduced-order lifted system formed from the reduced-order matrices in (51). Then the H∞-norm error bound for the reduced-order system is given by (53)H-H~H∞≤2tracediagΣ0,2,…,ΣK-1,2,where Σk,2,k=0,1,…,K-1, contains the truncated causal Hankel singular values from (49).


## 6. Numerical Experiments

In this section, we present some numerical results. We consider an artificial periodic discrete-time descriptor system which is reformed from its original model in [8, Example 1]. The reformulation of the model is described in [20]. Here all the computations are carried out using MATLAB® R2012b on an Intel CORE i5 processor with a 2.90-GHz clock and 16 GB of RAM.To show the comparisons of GLR-ADI (i.e., Algorithm  2 in [20]) and our UGLR-ADI (i.e., Algorithm 1), we compute the low-rank Gramian factors by solving the LPDALEs applying both the algorithms. We also investigate the performances of both theheuristic andadaptive shifts to implement these algorithms. For implementing the GLR-ADI iteration we chose 40 optimal heuristic shifts out of 50 large and 40 small magnitude Ritz values. To select 40 heuristic shifts it takes 3.98 seconds of CPU time. In case of the adaptive shifts (for implementing UGLR-ADI method), in each cycle 40 proper shift parameters are selected following the procedure discussed above. We require only 0.03 seconds of CPU time to compute this number of adaptive shift parameters. Note that for computing the initial shifts first we project the pencil (A-λE) onto the column space of an nc×100 random matrix.Figures1 and 2, respectively, show the rate of convergence and the computational time of the GLR-ADI and UGLR-ADI methods. In both cases, we see that UGLR-ADI method performs better than the GLR-ADI method. We compute the reduced-order model of the periodic descriptor system by using the low-rank Gramian factors obtained from UGLR-ADI method. The reformed periodic descriptor system has μk=nk=404, mk=2, and qk=3 for k=0,1,…,9. The set of periodic matrix pairs {(Ek,Ak)}k=0K-1 is periodic stable with nkf=400 and nk∞=4 for k=0,1,…,9. The original lifted system has order n=4040. The computed reduced-order model has subsystems of orders (9,8,9,9,9,9,9,9,9,9). The comparisons of the frequency responses and the error bounds of the original model and the reduced model are shown in Figure 3. Figure 3(a) shows the frequency responses of the original and reduced models match nicely. The estimated error is below the error bound which is reflected in Figure 3(b).Figure 1
Comparisons of the old algorithm (GLR-ADI) and new algorithm (UGLR-ADI) in terms of convergence rate in computing the low-rank Gramian factors.
(a)
Controllability Gramian factor (b)
Observability Gramian factorFigure 2
Comparisons of the old algorithm (GLR-ADI) and new algorithm (UGLR-ADI) in terms of required time in computing the low-rank Gramian factors.
(a)
Controllability Gramian factor (b)
Observability Gramian factorFigure 3
Frequency responses and error bounds of the original and the reduced-order model.
(a)
Frequency responses (b)
Absolute error and error bound


## 7. Conclusions

In this article we have suggested the adaptive shifts based ADI method and reviewed the Smith method for computing the low-rank factors of the periodic Gramians for periodic discrete-time descriptor systems. Later, these factors are used in a balancing based model reduction technique to find a reduced-order model for the periodic discrete-time descriptor system. The proposed techniques save memory and the computational time to generate a reduced-order model. The capability and efficiency of the methods have been reflected in the numerical results.



---

*Source: 4362641-2017-02-14.xml*