# Recursive Reduced-Order Algorithm for Singularly Perturbed Cross Grammian Algebraic Sylvester Equation

Authors: Intessar Al-Iedani; Zoran Gajic
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A new recursive algorithm is developed for solving the algebraic Sylvester equation that defines the cross Grammian of singularly perturbed linear systems. The cross Grammian matrix provides aggregate information about controllability and observability of a linear system. The solution is obtained in terms of reduced-order algebraic Sylvester equations that correspond to slow and fast subsystems of a singularly perturbed system. The rate of convergence of the proposed algorithm isOε, where ε is a small singular perturbation parameter that indicates separation of slow and fast state variables. Several real physical system examples are solved to demonstrate efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

Singularly perturbed systems have multiple time scales corresponding to fast and slow state space variables. For a system with two time scales, the slow time scale is related to the eigenvalues that are close to the imaginary axis and that represent the slow state space variables (slow modes) of the system, while the fast time scale is related to those that are far from the imaginary axis and that represent the fast state space variables (fast modes) of the system. Many algorithms exist in the literature for solving diverse problems related to analysis and control of singularly perturbed linear systems. Fixed point recursive numerical methods were first proposed in [1] and used in [2–4] to solve the closed and open loop optimal control problems. Those methods led thereafter to the Hamiltonian approach, which solves the linear-quadratic optimal control and filtering problem by decomposing the algebraic Riccati equations into pure-slow and pure-fast reduced-order algebraic Riccati equations [5]. The exact decomposition into pure-slow and pure-fast subsystems led to the use of parallel algorithms [6, 7] to solve the algebraic Riccati equation of the linear-quadratic optimal control problem. Moreover, some iterative methods were also used to solve this problem (see, e.g., [8] and the references therein). Most of the previous studies consider solving the algebraic Riccati equation, as it represents the most important equation of the optimal control and filtering problems.The system under investigation in this paper must be asymptotically stable, controllable, and observable. The test for controllability and observability of the system is usually done separately using the controllability and observability Grammians. In many applications, the reduced-order system is welcome to lower computational complexity. Model order reduction retains only state space variables that are both strongly controllable and strongly observable. This requires investigating the behavior of state space variables and balancing the controllability and observability Grammians, such that they are diagonal and identical. It has been shown in [9] that studying the controllability and observability of the system, separately, can be misleading; a method that directly assesses the combination of the two properties is preferred. Therefore, the cross Grammian matrix was defined in [10] as an alternative approach to the existing controllability and observability Grammian matrices. Unlike the controllability and observability Grammians, the cross Grammian contains information about both controllability and observability of the system.In this paper, a new recursive algorithm is proposed to solve the algebraic Sylvester equation of linear singularly perturbed systems whose solution defines the cross Grammian matrix. The algorithm is obtained in terms ofreduced-order algebraic Sylvester equations corresponding to slow and fast subsystems. The solutions of full-order algebraic Sylvester equations for finding the cross Grammian matrix were considered in [11, 12].The remainder of the paper is organized as follows. Section2 reviews the controllability, observability, and the cross Grammian matrices. The proposed recursive algorithm is then described in Section 3. In Section 4, several case studies are considered to demonstrate the performance of the proposed algorithm. Then, the conclusions follow in Section 5.


## 2. The Cross Grammian Matrix

Consider a linear dynamic system(1)dxdt=Axt+But,yt=Cxt+Dut,where x(t)∈Rn are state variables, u(t)∈Rm are control inputs, and y(t)∈Rp are measured outputs. Assume that system (1) is asymptotically stable, controllable, and observable. Controllability and observability of the system can be measured using the controllability and the observability Grammians defined, respectively, as (2)WC=∫0∞eAtBBTeATt,(3)WO=∫0∞eATtCTCeAt.For asymptotically stable, controllable, and observable systems, Grammians (2) and (3) are positive definite and represented the solutions of the algebraic Lyapunov equations(4)AWC+WCAT=-BBT,ATWO+WOA=-CTC.Assuming system (1) is square, that is, the number of inputs equals the number of outputs m=p, the cross Grammian matrix was defined in [10] for single-input single-output (SISO) systems as(5)WX=∫0∞eAtBCeAtand represented by the solution to the algebraic Sylvester equation(6)AWX+WXA+BC=0.In this context, the Sylvester algebraic equation (6) has a unique solution if and only if A and -A have distinct eigenvalues [13]. Several numerical solutions for the Sylvester equation were proposed in the literature; see, for example, [14–16] and the references therein. The definition in (5) was extended in [17–19] to include multi-input multioutput (MIMO) systems.Furthermore, for MIMO symmetric systems, the relation between controllability and observability, on the one hand, and the cross Grammian, on the other hand, is given by [20](7)WX2=WCWO.


## 3. A Recursive Algorithm for Finding Cross Grammians for Singularly Perturbed Linear Systems

The singularly perturbed structure can be obtained by partitioning the system matrices in (1) as follows [6, 21]:(8)A=A1A2A3εA4ε,B=B1B2ε,C=C1C2,where ε is a small positive singular perturbation parameter. A, B, and C are constant matrices of appropriate dimensions. Based on the singular perturbation theory [6, 21], a singularly perturbed linear system in the explicit state variable standard form is given by(9)dx1tdt=A1x1t+A2x2t+B1ut,εdx2tdt=A3x1t+A4x2t+B2ut,yt=C1x1t+C2x2t+Dut,where x1t∈Rn1 are the slow state variables and x2t∈Rn2 are the fast state variables. Assuming that A4 is nonsingular, the eigenvalues of matrix A consist of two disjoint groups: one corresponds to the slow subsystem λs(A) and the other corresponds to the fast subsystem λf(A). If the two subsystems have a mixture of slow and fast eigenvalues, then a technique has to be applied to convert the system into its standard singularly perturbed form defined in (9). We will give examples on this case in Sections 4.2 and 4.3.The nature of the cross Grammian matrixWX defined in (6) corresponding to the system singularly perturbed form defined in (9) is(10)WX=W1εW2W3W4.Using (8) and (10) in (6), we get the partitioned form of the algebraic Sylvester equation as follows:(11)A1W1+A2W3+W1A1+W2A3+B1C1=0,εA1W2+A2W4+W1A2+W2A4+B1C2=0,A3W1+A4W3+εW3A1+W4A3+B2C1=0,εA3W2+A4W4+εW3A2+W4A4+B2C2=0.Setting ε=0, we get the following approximate algebraic equations: (12)A1W10+A2W30+W10A1+W20A3+B1C1=0,A2W40+W10A2+W20A4+B1C2=0,A3W10+A4W30+W40A3+B2C1=0,A4W40+W40A4+B2C2=0.The solution of equations (12) is given in terms of the following reduced-order algebraic Sylvester equations corresponding to the slow and fast subsystems:(13)A4W40+W40A4+B2C2=0,A0W10+W10A0+G0=0.In addition, we have from (12)(14)W20=-B1C2+A2W40+W10A2A4-1,W30=-A4-1B2C1+A3W10+W40A3,where (15)A0=A1-A2A4-1A3,(16)G0=-A2A4-1B2C1+W40A3-A2W40+B1C2A4-1A3+B1C1.To find a unique solution of (13), we impose the following assumption.Assumption 1.
MatricesA0 and A4 are asymptotically stable. In consequence, unique solutions of (13)-(14) exist.Defining the approximation error as(17)W1=W10+εE1,W2=W20+εE2,W3=W30+εE3,W4=W40+εE4and subtracting (12) from (11), we get the following error equations, after some algebra:(18)A4E4+E4A4=-A3W20+εE2-W30+εE3A2,A3E1+A4E3+E4A3=-W30+εE3A1,A2E4+E1A2+E2A4=-A1W20+εE2,A1E1+A2E3+E1A1+E2A3=0.From the first equation in (18), we can observe that the unknown errors E2 and E3 are multiplied by a small parameter ε. A similar situation is in the second and the third equations of (18). Therefore, we propose the following algorithm for solving error equations (18).


### 3.1. The Proposed Algorithm

Start withE20=0 and E30=0 and recursively evaluate(19)A4E4i+1+E4i+1A4=-A3W20+εE2i-W30+εE3iA2,A0E1i+1+E1i+1A0=A2A4-1W30+εE3iA1+A1W20+εE2iA4-1A3+A2A4-1E4i+1A3+A2E4i+1A4-1A3,E2i+1=-A1W20+εE2i+A2E4i+1+E1i+1A2A4-1,E3i+1=-A4-1W30+εE3iA1+A3E1i+1+E4i+1A3,fori=0,1,2,….Theorem 2.
Assuming that matricesA0 and A4 are asymptotically stable, algorithm (19) converges to the exact solution of (18) with a rate of convergence Oε; that is,(20)Eji+1-Eji=Oε,Eji-Ej=Oεi,for j=1,2,3,4 and i=1,2,….Therefore, the exact solutionWX can be obtained with an accuracy of Oεi after performing i iterations on the proposed algorithm (19) as follows:(21)Wji=Wj+εEji,forj=1,2,3,4,i=1,2,….Proof.
Using Assumption1, that is, A0 and A4 are asymptotically stable, it can be shown that (19) represents a contraction mapping [22]; that is,(22)Eji-Ej=Oεi,forj=1,2,3,4,i=1,2,….
Formula (22) will be also valid if (23)E2i+1=-A1W20+εE2i+A2E4i+E1iA2A4-1,E3i+1=-A4-1W30+εE3iA1+A3E1i+E4iA3,fori=0,1,2,…,E10=0,E40=0.Formula (22) implies that algorithm (19) is convergent. Using Ej∞, for j=1,2,3,4, in (19) and comparing it to (18), it can be seen that algorithm (19) converges to the unique solution of (18).


## 3.1. The Proposed Algorithm

Start withE20=0 and E30=0 and recursively evaluate(19)A4E4i+1+E4i+1A4=-A3W20+εE2i-W30+εE3iA2,A0E1i+1+E1i+1A0=A2A4-1W30+εE3iA1+A1W20+εE2iA4-1A3+A2A4-1E4i+1A3+A2E4i+1A4-1A3,E2i+1=-A1W20+εE2i+A2E4i+1+E1i+1A2A4-1,E3i+1=-A4-1W30+εE3iA1+A3E1i+1+E4i+1A3,fori=0,1,2,….Theorem 2.
Assuming that matricesA0 and A4 are asymptotically stable, algorithm (19) converges to the exact solution of (18) with a rate of convergence Oε; that is,(20)Eji+1-Eji=Oε,Eji-Ej=Oεi,for j=1,2,3,4 and i=1,2,….Therefore, the exact solutionWX can be obtained with an accuracy of Oεi after performing i iterations on the proposed algorithm (19) as follows:(21)Wji=Wj+εEji,forj=1,2,3,4,i=1,2,….Proof.
Using Assumption1, that is, A0 and A4 are asymptotically stable, it can be shown that (19) represents a contraction mapping [22]; that is,(22)Eji-Ej=Oεi,forj=1,2,3,4,i=1,2,….
Formula (22) will be also valid if (23)E2i+1=-A1W20+εE2i+A2E4i+E1iA2A4-1,E3i+1=-A4-1W30+εE3iA1+A3E1i+E4iA3,fori=0,1,2,…,E10=0,E40=0.Formula (22) implies that algorithm (19) is convergent. Using Ej∞, for j=1,2,3,4, in (19) and comparing it to (18), it can be seen that algorithm (19) converges to the unique solution of (18).


## 4. Case Studies

Three case studies are considered to demonstrate the proposed algorithm: a fourth-order aircraft example whose mathematical model is in the explicit singularly perturbed form defined in (9) in which with accuracy of Oε the slow eigenvalues are all contained in the approximate slow subsystem represented by A0 and all fast eigenvalues are contained in the approximate fast subsystem represented by A4; a fifth-order chemical plant model given in implicit singularly perturbed form (it has two slow and three fast eigenvalues, but the state variables have to be reordered to achieve explicit singularly perturbed form defined in (9)); a tenth-order hydrogen gas reformer used to provide hydrogen to a fuel cell from hydrogen rich fuels (natural gas, methanol).


### 4.1. L-1011 Aircraft

Here, we consider the lateral axis equations of the rigid body model of L-1011 aircraft at cruise condition [23]. The state variables are the bank angle, roll rate, yaw rate, and sideslip angle, which are represented in the state vector xt=x1x2x3x4T in the same order. The input vector consists of two variables, the rudder deflection δr and the aileron deflection δa, and is given as ut=[δrδa]T. The system matrices are given as(24)A=01000-1.890.39-5.550-0.034-2.982.430.034-0.0011-0.99-0.21,B=000.36-1.6-0.95-0.0320.030,C=10000100.The eigenvalues of the matrix A are −0.1016, −1.4811±0.6292i, and −2.0162. The system is asymptotically stable (all eigenvalues are in the left half plane), controllable, and observable. Moreover, there is only one slow mode with eigenvalue −0.0899, and there are three fast modes with eigenvalues −1.4891±0.7686i and −2.1017. The singular perturbation parameter ε=0.07, which is the ratio between the fastest slow eigenvalue and the slowest fast eigenvalue. Solving the algebraic Sylvester’s equation (6), the cross Grammian matrix can be obtained as follows: (25)WX=-3.77168467243-2.25320146563-4.6028545113112.68652326564-0.43134179103-0.62940797540-0.278099839141.332572206696-0.143139895690.00048741612-0.01060057667-0.021220541030.209670365570.06408641105-0.038614510040.00250243909.Using the proposed algorithm, the initial cross Grammian matrixW(0) (first-order approximate solution) is obtained as follows: (26)W0=-3.98315817315-2.25882783809-4.4622998610912.32553354098-0.41906814039-0.42418654922-0.214432523620.883962065453-0.14908955443-0.002849079280.00235249833-0.005866460180.2023309755780.0015539806820.000836843461-0.00419856136.Comparing the exact solutionWX to the first-order approximate solution of the cross Grammian matrix by calculating the error norm, we get(27)WX-W02=0.624920763816596.Then, the cross Grammian matrix is calculated using the proposed recursive algorithm. The error norm at each iteration is shown in Table 1. By taking the error norm, it can be seen that the algorithm converges rapidly to the exact solution.Table 1
Error norm values for each iteration for L-1011 aircraft system.
Number of iterations,i W X - W X i 2 2 7.503328625983908 e - 04 3 2.611125417476592 e - 05 4 9.085062505535487 e - 07 5 3.161097999142058 e - 08 6 1.099918110580769 e - 09 7 3.830668858621788 e - 11 8 1.368874658220410 e - 12 9 9.598118726733237 e - 14 10 6.474407682370149 e - 14


### 4.2. Chemical Plant

In this section, the linearized chemical plant considered in [24] is chosen to explain the behavior of the proposed algorithm when the linear singularly perturbed system is not in the explicit standard form (9). The system matrices are given as follows:(28)A=-0.10940.06280001.306-2.1360.98070001.595-3.1491.547000.03552.632-4.2571.85500.002700.1636-0.1625,B=000.063800.0838-0.13960.1004-0.2060.0063-0.0128,C=1000000001.The eigenvalues of matrix A are −5.9822, −2.8408, −0.8954, −0.0774, and −0.0141, which indicates that this system has two slow modes (eigenvalues). The singular perturbation parameter ε is chosen as the ratio of the fastest slow eigenvalue to the slowest fast eigenvalue and is equal to ε=0.086=0.0774/0.8954. Introducing the following permutation matrix that exchanges the second row of matrix A, with its fifth row, that is,(29)P=1000000001001000001001000,the explicit singularly perturbed form of the system matrices can be obtained as follows:(30)ASP=PAP,BSP=PB,CSP=CP.Thereby, they are calculated as(31)ASP=-0.10940000.06280-0.162500.16360.002700-3.1491.5471.595001.85502.632-4.2570.03551.306000.98070-2.1320,BSP=000.0063-0.01280.0838-0.13960.1004-0.2060.06380,CSP=1000001000.Matrices ASP and A have the same eigenvalues and the same number of slow and fast modes. However, the slow and fast eigenvalues are now correctly separated into two disjoint groups. The slow approximate subsystem (see (15)) has the eigenvalues −0.0788 and −0.0160, and the fast approximate subsystem has the fast eigenvalues −5.9521, 2.7925, and -0.7933 (eigenvalues of A4).Solving the algebraic Sylvester’s equation (6), the cross Grammian matrix can be obtained as follows:(32)WX=-0.07530866-0.11561658-0.00985227-0.00806596-0.00985946-0.23194836-0.49463876-0.03897379-0.03313413-0.03676743-0.07341430-0.42776115-0.02785663-0.02552098-0.02384797-0.12954567-0.52436133-0.03510447-0.03163354-0.03102749-0.05734203-0.26232141-0.01916962-0.01692499-0.01682689.Calculating the initial cross Grammian matrix W(0) using the proposed algorithm and comparing the result to the exact solution of the cross Grammian, using the error norm, we get the following results:(33)W0=-0.07894318-0.13640657-0.01124158-0.00932742-0.0110635-0.25878609-0.56732132-0.04537161-0.03829073-0.04292241-0.05767992-0.50330195000-0.12521917-0.60940965000-0.04496552-0.31507279000,WX-W02=0.160161520985517.The cross Grammian matrix is then calculated using the proposed recursive algorithm. The error norm for each iteration is shown in Table 2.Table 2
Error norm values for each iteration for the chemical plant.
Number of iterations,i W X - W X i 2 2 6.549205248602 e - 3 4 3.305462486403183 e - 04 6 1.644553914568401 e - 05 8 7.961048437784984 e - 07 10 3.791658707506841 e - 08 12 1.790737974576206 e - 09 14 8.421250159940589 e - 11 16 3.944852245039569 e - 12 18 1.781054537942190 e - 13 20 1.746372245754009 e - 14It can be seen from Table2 that the proposed algorithm converges to the exact solution according to the convergence result stated in Theorem 2.


### 4.3. Natural Gas Hydrogen Reformer

In this section, we investigate the behavior of the proposed algorithm in case of higher order singularly perturbed systems. The linearized 10th-order mathematical model of the gas hydrogen reformer introduced and studied in [25–27] is chosen. The system matrices are given as follows [25]:(34)A=-0.07400000-3.531.0748000-1.468-25.3000002.558213.91100-15600000033.586000-124.5212.630112.69112.69000000-3.3330000000000-32.4332.30432.3040000000331.8-344-34109.9042000221.9700-253.2-254.9032.526002.03540001.83091.214-0.358-3.3040.018808.16420005.60435.39940-13.61,B=00000.1200000000000.18340000T,C=100000000000.994-0.0880000000.The eigenvalues of the system matrixA are as follows: −660.68, −157.9, −89.137, −12.175, −3.33, -2.77±0.547i, −1.468, −0.358, and −0.0861. All real parts of those eigenvalues lie in the left part of the complex plane; hence, the system is asymptotically stable. Moreover, the system has multiple time scales (slow and fast) since there are three eigenvalues located very close to the imaginary axis while the other seven eigenvalues are located far from that axis. The singular perturbation parameter ε is chosen as the ratio of the fastest slow eigenvalue to the slowest fast eigenvalue and is equal to ε=0.52=1.468/2.77.What are supposed to be slow mode eigenvalues, obtained via (15), −1.468, −0.0838, and −128.4495, and what are supposed to be the fast mode eigenvalues (eigenvalues of A4), −0.358, −660.68, −89.144, −13.954, −2.829±0.886, and −3.333, do not display the slow-fast time scale separation. Similar to the chemical plant example in the previous section, we can notice that the slow and fast eigenvalues are not well separated into two disjoint groups. However, since the gas former system is of higher order, it will be a little cumbersome to come up with a permutation matrix that would convert the system into its explicit singularly perturbed form. Therefore, the algorithm presented in [21, 28] and used by [27] to study the slow and fast dynamics of the gas reformer system is also considered here to convert the system from an implicit singularly perturbed form to the explicit singularly perturbed form. The algorithm in [21, 28] is based on introducing a similarity transformation T that transforms the general linear system in implicit singularly perturbed form into the explicit singularly perturbed form defined in (9). The similarity transformation T is given by [27](35)T=4.85572.1797-0.35351.7829113.74-52.659-5.14681000000000010-0.37940.7980-0.0771005.26670.51480010010000000000100000000001000000000010000000000100000000001000000000001.The new transformed system matrices are defined as(36)ASP=TAT-1,BSP=TB,CSP=CT-1.Thereby, they are calculated as(37)ASP=-0.31925.5761-3.138-0.3530.568136.3094-0.259-0.0505.51563.13790-0.35802.03540001.8311.2140-3.304-0.08662.0413-1.2315-0.1250.15449.84731.92520.1901-0.31951.2314000-124.50000033.5860000-124.5212.630112.69112.69000000-3.3330000000000-32.4332.30432.3040000221.9700331.8-344-3419.9042002.03540221.9700-253.2-254.932.5260.018808.16420-0.0057-0.36290.2135.6255.3962-13.61,BSP=13.648800000.120000-9.657700.96590000.1834000T,CSP=0.16970-0.46360.0243-0.302-19.30311.37881.112-0.16970.46360.080201.02660.0195-0.143-9.1238-1.1825-0.115-0.0802-1.0266.The matricesASP and A have the same eigenvalues, since they are preserved under the similarity transformation. The slow mode eigenvalues are −1.468, −0.358, and −0.08552, while the fast mode eigenvalues are −660.682, −157.89, 89.137, −12.174, −3.333, and −2.7697±0.60087. They are clearly separated now into two disjoint groups.Using our proposed algorithm to calculate the initial cross Grammian matrixW(0) and compare the result to the exact solution of the algebraic Sylvester equation, we get the error norm as(38)WX-W02=53.489147248362102.The cross Grammian matrix is then calculated using the proposed recursive algorithm. The error norm for each iteration is shown in Table3.Table 3
Error norm values for each iteration for the gas reformer system.
Number of iterations,i W X - W X i 2 2 56.670747861070758 5 13.260463185751986 10 0.386489799910691 12 50.958037951316 e - 02 15 5.949048480951 e - 03 20 6.505261253303921 e - 04 23 6.630466299301564 e - 05 25 7.912793405030453 e - 06 30 7.210694103862186 e - 07 35 3.050169378906093 e - 08 40 8.638446874677160 e - 10


## 4.1. L-1011 Aircraft

Here, we consider the lateral axis equations of the rigid body model of L-1011 aircraft at cruise condition [23]. The state variables are the bank angle, roll rate, yaw rate, and sideslip angle, which are represented in the state vector xt=x1x2x3x4T in the same order. The input vector consists of two variables, the rudder deflection δr and the aileron deflection δa, and is given as ut=[δrδa]T. The system matrices are given as(24)A=01000-1.890.39-5.550-0.034-2.982.430.034-0.0011-0.99-0.21,B=000.36-1.6-0.95-0.0320.030,C=10000100.The eigenvalues of the matrix A are −0.1016, −1.4811±0.6292i, and −2.0162. The system is asymptotically stable (all eigenvalues are in the left half plane), controllable, and observable. Moreover, there is only one slow mode with eigenvalue −0.0899, and there are three fast modes with eigenvalues −1.4891±0.7686i and −2.1017. The singular perturbation parameter ε=0.07, which is the ratio between the fastest slow eigenvalue and the slowest fast eigenvalue. Solving the algebraic Sylvester’s equation (6), the cross Grammian matrix can be obtained as follows: (25)WX=-3.77168467243-2.25320146563-4.6028545113112.68652326564-0.43134179103-0.62940797540-0.278099839141.332572206696-0.143139895690.00048741612-0.01060057667-0.021220541030.209670365570.06408641105-0.038614510040.00250243909.Using the proposed algorithm, the initial cross Grammian matrixW(0) (first-order approximate solution) is obtained as follows: (26)W0=-3.98315817315-2.25882783809-4.4622998610912.32553354098-0.41906814039-0.42418654922-0.214432523620.883962065453-0.14908955443-0.002849079280.00235249833-0.005866460180.2023309755780.0015539806820.000836843461-0.00419856136.Comparing the exact solutionWX to the first-order approximate solution of the cross Grammian matrix by calculating the error norm, we get(27)WX-W02=0.624920763816596.Then, the cross Grammian matrix is calculated using the proposed recursive algorithm. The error norm at each iteration is shown in Table 1. By taking the error norm, it can be seen that the algorithm converges rapidly to the exact solution.Table 1
Error norm values for each iteration for L-1011 aircraft system.
Number of iterations,i W X - W X i 2 2 7.503328625983908 e - 04 3 2.611125417476592 e - 05 4 9.085062505535487 e - 07 5 3.161097999142058 e - 08 6 1.099918110580769 e - 09 7 3.830668858621788 e - 11 8 1.368874658220410 e - 12 9 9.598118726733237 e - 14 10 6.474407682370149 e - 14


## 4.2. Chemical Plant

In this section, the linearized chemical plant considered in [24] is chosen to explain the behavior of the proposed algorithm when the linear singularly perturbed system is not in the explicit standard form (9). The system matrices are given as follows:(28)A=-0.10940.06280001.306-2.1360.98070001.595-3.1491.547000.03552.632-4.2571.85500.002700.1636-0.1625,B=000.063800.0838-0.13960.1004-0.2060.0063-0.0128,C=1000000001.The eigenvalues of matrix A are −5.9822, −2.8408, −0.8954, −0.0774, and −0.0141, which indicates that this system has two slow modes (eigenvalues). The singular perturbation parameter ε is chosen as the ratio of the fastest slow eigenvalue to the slowest fast eigenvalue and is equal to ε=0.086=0.0774/0.8954. Introducing the following permutation matrix that exchanges the second row of matrix A, with its fifth row, that is,(29)P=1000000001001000001001000,the explicit singularly perturbed form of the system matrices can be obtained as follows:(30)ASP=PAP,BSP=PB,CSP=CP.Thereby, they are calculated as(31)ASP=-0.10940000.06280-0.162500.16360.002700-3.1491.5471.595001.85502.632-4.2570.03551.306000.98070-2.1320,BSP=000.0063-0.01280.0838-0.13960.1004-0.2060.06380,CSP=1000001000.Matrices ASP and A have the same eigenvalues and the same number of slow and fast modes. However, the slow and fast eigenvalues are now correctly separated into two disjoint groups. The slow approximate subsystem (see (15)) has the eigenvalues −0.0788 and −0.0160, and the fast approximate subsystem has the fast eigenvalues −5.9521, 2.7925, and -0.7933 (eigenvalues of A4).Solving the algebraic Sylvester’s equation (6), the cross Grammian matrix can be obtained as follows:(32)WX=-0.07530866-0.11561658-0.00985227-0.00806596-0.00985946-0.23194836-0.49463876-0.03897379-0.03313413-0.03676743-0.07341430-0.42776115-0.02785663-0.02552098-0.02384797-0.12954567-0.52436133-0.03510447-0.03163354-0.03102749-0.05734203-0.26232141-0.01916962-0.01692499-0.01682689.Calculating the initial cross Grammian matrix W(0) using the proposed algorithm and comparing the result to the exact solution of the cross Grammian, using the error norm, we get the following results:(33)W0=-0.07894318-0.13640657-0.01124158-0.00932742-0.0110635-0.25878609-0.56732132-0.04537161-0.03829073-0.04292241-0.05767992-0.50330195000-0.12521917-0.60940965000-0.04496552-0.31507279000,WX-W02=0.160161520985517.The cross Grammian matrix is then calculated using the proposed recursive algorithm. The error norm for each iteration is shown in Table 2.Table 2
Error norm values for each iteration for the chemical plant.
Number of iterations,i W X - W X i 2 2 6.549205248602 e - 3 4 3.305462486403183 e - 04 6 1.644553914568401 e - 05 8 7.961048437784984 e - 07 10 3.791658707506841 e - 08 12 1.790737974576206 e - 09 14 8.421250159940589 e - 11 16 3.944852245039569 e - 12 18 1.781054537942190 e - 13 20 1.746372245754009 e - 14It can be seen from Table2 that the proposed algorithm converges to the exact solution according to the convergence result stated in Theorem 2.


## 4.3. Natural Gas Hydrogen Reformer

In this section, we investigate the behavior of the proposed algorithm in case of higher order singularly perturbed systems. The linearized 10th-order mathematical model of the gas hydrogen reformer introduced and studied in [25–27] is chosen. The system matrices are given as follows [25]:(34)A=-0.07400000-3.531.0748000-1.468-25.3000002.558213.91100-15600000033.586000-124.5212.630112.69112.69000000-3.3330000000000-32.4332.30432.3040000000331.8-344-34109.9042000221.9700-253.2-254.9032.526002.03540001.83091.214-0.358-3.3040.018808.16420005.60435.39940-13.61,B=00000.1200000000000.18340000T,C=100000000000.994-0.0880000000.The eigenvalues of the system matrixA are as follows: −660.68, −157.9, −89.137, −12.175, −3.33, -2.77±0.547i, −1.468, −0.358, and −0.0861. All real parts of those eigenvalues lie in the left part of the complex plane; hence, the system is asymptotically stable. Moreover, the system has multiple time scales (slow and fast) since there are three eigenvalues located very close to the imaginary axis while the other seven eigenvalues are located far from that axis. The singular perturbation parameter ε is chosen as the ratio of the fastest slow eigenvalue to the slowest fast eigenvalue and is equal to ε=0.52=1.468/2.77.What are supposed to be slow mode eigenvalues, obtained via (15), −1.468, −0.0838, and −128.4495, and what are supposed to be the fast mode eigenvalues (eigenvalues of A4), −0.358, −660.68, −89.144, −13.954, −2.829±0.886, and −3.333, do not display the slow-fast time scale separation. Similar to the chemical plant example in the previous section, we can notice that the slow and fast eigenvalues are not well separated into two disjoint groups. However, since the gas former system is of higher order, it will be a little cumbersome to come up with a permutation matrix that would convert the system into its explicit singularly perturbed form. Therefore, the algorithm presented in [21, 28] and used by [27] to study the slow and fast dynamics of the gas reformer system is also considered here to convert the system from an implicit singularly perturbed form to the explicit singularly perturbed form. The algorithm in [21, 28] is based on introducing a similarity transformation T that transforms the general linear system in implicit singularly perturbed form into the explicit singularly perturbed form defined in (9). The similarity transformation T is given by [27](35)T=4.85572.1797-0.35351.7829113.74-52.659-5.14681000000000010-0.37940.7980-0.0771005.26670.51480010010000000000100000000001000000000010000000000100000000001000000000001.The new transformed system matrices are defined as(36)ASP=TAT-1,BSP=TB,CSP=CT-1.Thereby, they are calculated as(37)ASP=-0.31925.5761-3.138-0.3530.568136.3094-0.259-0.0505.51563.13790-0.35802.03540001.8311.2140-3.304-0.08662.0413-1.2315-0.1250.15449.84731.92520.1901-0.31951.2314000-124.50000033.5860000-124.5212.630112.69112.69000000-3.3330000000000-32.4332.30432.3040000221.9700331.8-344-3419.9042002.03540221.9700-253.2-254.932.5260.018808.16420-0.0057-0.36290.2135.6255.3962-13.61,BSP=13.648800000.120000-9.657700.96590000.1834000T,CSP=0.16970-0.46360.0243-0.302-19.30311.37881.112-0.16970.46360.080201.02660.0195-0.143-9.1238-1.1825-0.115-0.0802-1.0266.The matricesASP and A have the same eigenvalues, since they are preserved under the similarity transformation. The slow mode eigenvalues are −1.468, −0.358, and −0.08552, while the fast mode eigenvalues are −660.682, −157.89, 89.137, −12.174, −3.333, and −2.7697±0.60087. They are clearly separated now into two disjoint groups.Using our proposed algorithm to calculate the initial cross Grammian matrixW(0) and compare the result to the exact solution of the algebraic Sylvester equation, we get the error norm as(38)WX-W02=53.489147248362102.The cross Grammian matrix is then calculated using the proposed recursive algorithm. The error norm for each iteration is shown in Table3.Table 3
Error norm values for each iteration for the gas reformer system.
Number of iterations,i W X - W X i 2 2 56.670747861070758 5 13.260463185751986 10 0.386489799910691 12 50.958037951316 e - 02 15 5.949048480951 e - 03 20 6.505261253303921 e - 04 23 6.630466299301564 e - 05 25 7.912793405030453 e - 06 30 7.210694103862186 e - 07 35 3.050169378906093 e - 08 40 8.638446874677160 e - 10


## 5. Conclusions

The algorithm was developed to solve the algebraic Sylvester equation whose solution defines the cross Grammian of singularly perturbed linear systems. The algorithm is very efficient, defined in terms of reduced-order subproblems corresponding to slow and fast subsystems, and converges rapidly to the required solution. The efficacy of the algorithm is demonstrated on three real physical examples. The algorithm can be directly applied to singularly perturbed systems in the explicit standard forms. A similarity transformation needs to be applied to singularly perturbed systems in implicit forms to convert them into their explicit form before the proposed algorithm can be applied to this class of systems.



---

*Source: 2452746-2016-12-12.xml*