A new method for model reduction of linear systems is presented, based on Chebyshev rational functions, using the Harmony Search (HS) algorithm. First, the full order system is expanded and then a set of parameters in a fixed structure are determined, whose values define the reduced order system. The values are obtained by minimizing the errors between the first coefficients of the Chebyshev rational function expansion of full and reduced systems, using the HS algorithm. To assure stability, the Routh criterion is used as constraints in the optimization problem. To present the ability of the proposed method, three test systems are reduced. The results obtained are compared with other existing techniques. The results obtained show the accuracy and efficiency of the proposed method.
Chebyshev rational functions ; Harmony search algorithm ; Order reduction ; Orthogonal functions ; Routh array
Reduced order models are highly desirable for engineers in system analysis, synthesis and simulation of complicated highorder systems, since the analysis and design of such systems is not an easy task. Various methods for model reduction are reported in the literature in the time and frequency domains. Model reduction was started by Davison in 1966 [1] and followed by Chidambara, who suggested several modifications to Davison’s approach [2] , [3] and [4] . After that, different approaches were proposed using dominant eigenvalue retention [1] and [5] , Routh approximation [6] , Hurwitz polynomial approximation [7] and [8] , the stability equation method [9] and [10] , moments matching [11] , [12] , [13] and [14] , the continued fraction method [15] , [16] and [17] , Pade approximation [18] and etc.
The issue of optimality in model reduction was considered by Wilson [19] and [20] , who suggested an optimization approach based on minimization of the integral squared impulse response error between full and reducedorder models. This attempt was continued by other researchers through other approaches [21] , [22] , [23] and [24] .
In 1981 [25] , the controllability and observability of the states were considered in model reduction by Moore. The suggested approach suffered from steady state errors but the stability of the reduced model was assured if the original system was also stable [26] . Furthermore, the concept of , , and were used for model reduction in [27] , [28] , [29] and [30] .
In recent decades, evolutionary techniques, such as Particle Swarm Optimization (PSO) and Genetic Algorithm (GA), have been used for order reduction of systems [31] , [32] and [33] . In these approaches, the reduced order model parameters are achieved by minimizing a fitness function, which is often Integral Square Error (ISE), Integral Absolute Error (IAE), norm or norm [34] , [35] and [36] .
This paper proposes an alternative method for order reduction based on Chebyshev rational functions, using the HS algorithm. The full order system is expanded and then the first coefficients of Chebyshev rational function expansion are obtained. A desire fixed structure for the reduced order model is considered and a set of parameters are defined, whose values determine the reduced order system. These unknown parameters are determined using the harmony search algorithm by minimizing the errors between the first coefficients of Chebyshev rational function expansion of full and reduced systems. To assure stability, the Routh criterion is applied, as used in [37] , where it states optimization problems as constraints, which, subsequently, are converted to constrained optimization problems. To show the accuracy of the proposed method, three systems are reduced by the proposed method and compared with those available in the literature.
To make a proper background, Chebyshev rational functions and the harmony search are briefly explained in Sections 2 and 3 , respectively. The proposed method is explained in Section 4 . The ability of the proposed approach is shown in Section 5 , and, finally, the paper is concluded in Section 6 .
The Chebyshev rational functions are a sequence of functions that are both orthogonal and rational [38] . A rational Chebyshev function of degree is presented as:

( 1) 
where is a Chebyshev polynomial and defined by a recursive formula as:

( 2) 
Using the orthogonality relationship, an arbitrary function, , can be expanded as:

( 3) 
where is given by the following equations:

( 4) 
in which, equals 2 for and equals 1 for . Also, W is given as:

( 5) 
Therefore, by considering the first terms of Eq. (3) , a good approximant of is obtained.
The HS is based on a natural musical process, which searches for a perfect state of harmony. In general, the HS algorithm works as follows [39] and [40] :
Step 1. Initialization: Define the objective function and decision variables, and input the system parameters and the boundaries of the decision variables. The optimization problem can be defined as:
Minimize subject to where and are the lower and upper bounds for decision variables.
The HS algorithm parameters are also specified in this step. They are the Harmony Memory Size (HMS ) or the number of solution vectors in harmony memory, Harmony Memory Considering Rate (HMCR ), distance bandwidth (bw ), Pitch Adjusting Rate (PAR ), and the number of improvisations (K ), or stopping criterion. K is the same as the total number of function evaluations.
Step 2. Initialize the harmony memory (HM). The harmony memory is a memory location where all the solution vectors (sets of decision variables) are stored. The initial harmony memory is randomly generated in the region . This is done based on the following equation:

( 6) 
where is a random from a uniform distribution of .
Step 3. Improvise a new harmony from the harmony memory. Generating a new harmony, , is called improvisation, which is based on three rules: memory consideration, pitch adjustment and random selection. First of all, a uniform random number, r , is generated in the range . If r is less than HMCR, the decision variable, , is generated by the memory consideration; otherwise, is obtained by a random selection. Then, each decision variable, , will undergo a pitch adjustment with a probability of PAR if it is produced by the memory consideration. The pitch adjustment rule is given as follows:

( 7) 
Step 4. Update harmony memory. After generating a new harmony vector, , the harmony memory will be updated. If the fitness of the improvised harmony vector, , is better than that of the worst harmony, the worst harmony in the HM will be replaced with and become a new member of the HM.
Step 5. Repeat Steps 3–4 until the stopping criterion (maximum number of improvisations K ) is met.
Consider a stable SingleInput SingleOutput (SISO) system described by the transfer function of order n as follows:

( 8) 
where and are constants.
The objective is to obtain a reduced model of order r , where r is smaller than n , such that the principal and important specifications of the full order system are retained in the reduced order model. This reduced order system is presented as follows:

( 9) 
where and are unknown constants.
To obtain the reduced model by the proposed method, firstly, the full order system is expanded based on Chebyshev rational functions. Then, the first coefficients of the Chebyshev rational function expansion of the original system are obtained and shown by . Then, a desired fixed structure is considered for the reduced order model, as defined in Eq. (9) , where and are unknown parameters of the reduced order model that are obtained by HS. The goal of the optimization is to find the best parameters for . Therefore, each harmony is a dimensional vector, in which is . Each harmony is a solution to and for each solution (harmony), the Chebyshev rational function expansions are obtained. Each harmony is evaluated by minimizing the following fitness function:

( 10) 
in which are the coefficients of the Chebyshev rational function expansions of the reduced order system. The algorithm searches for the best harmony until the termination criteria are met. At this stage, the best parameters are given as parameters of the reduced order model.
Furthermore, the reduced model must be stable if the original system is stable. Therefore, the Routh criterion is applied to assure stability. For specifying the stability conditions, following [41] , the denominator polynomial of the reduced order model in Eq. (9) can be shown as below:

( 11) 
which is constructed by taking the coefficients of the first two rows of the Routh array, in which elements of its first column have the following entries:

( 12) 
where k is equal to 1 for even r , and k is equal to 0 for odd r .
Comparing the entries of the array in Eq. (12) with and those of the second row with gives Eq. (13) :

( 13) 
Substituting the above relations in the reduced order model denominators, Eq. (11) is achieved. Therefore, the necessary and sufficient condition for all roots of the reduced system to be strictly in the lefthalf plane is:

( 14) 
and, subsequently:

( 15) 
Thus, to have an optimum stable reduced system, the reduced order model’s parameters are determined by minimizing the following fitness function:

( 16) 
Therefore, the reduced order model is achieved, such that the first coefficients of the Chebyshev rational function expansions of the full order system are equal to the first coefficients of the Chebyshev rational function expansions of the reduced order model.
The proposed method can be summarized in the following steps:
Step 1: The Chebyshev rational functions of the full order system in Eq. (8) are obtained.
Step 2: A desire fixed structure is considered for the reduced order model, as defined in Eq. (9) , where and are unknown parameters of the reduced order model that are obtained in the next step.
Step 3: To obtain the unknown parameters, HS is applied. The goal of the optimization is to find the best parameters for . Therefore, each harmony is a dimensional vector, in which is . Each harmony is a solution to and, for each solution (harmony), the Chebyshev rational functions are obtained. Each harmony is evaluated using the objective function defined by Eq. (16) , searching for the harmony associated with until the termination criteria are met. At this stage, the best parameters are given as parameters of the reduced order model.
To assess the efficiency of the proposed approach, it has been applied to three test systems, where a stepbystep procedure is given for the first test system.
Test system 1: The first system to be reduced is a system of order 6 given by Mukherjee in [42] , where he uses a response matching technique to obtain the reduced system. is given in Box I .

( 17) 
The reducedorder model can be obtained by the following steps:
Step 1: Based on Section 2 , the Chebyshev rational function of the full order system in Eq. (17) is obtained as:

( 18) 
where the fifth order ( ) Chebyshev rational function of Eq. (17) is considered to present .
Step 2: The full order of the system represented in Eq. (17) is going to be reduced to a thirdorder system with the following transfer function:

( 19) 
where and are unknown parameters of the reduced order model.
Step 3: HS is applied to obtain the unknown parameters in Eq. (19) . Since, the goal of the optimization is to find the best parameters for , each harmony is a dimensional vector in which . The HMS is selected to be 6, and the HMCR and evaluation number are set to be 0.9 and 1000, respectively. Each harmony is a solution to and, for each solution (harmony), the Chebyshev rational functions are obtained. Each harmony is evaluated using the objective function defined by Eq. (16) , searching for the best , until the termination criteria are met. At this stage, the best parameters are given for the reduced order model, where the following reduced order model is obtained:

( 20) 
The Chebyshev rational functions of the obtained reduced order model are as:

( 21) 
Comparing Eqs. (18) and (21) shows that the best approximant of is achieved. The step response of the full order system and that of the system with a thirdorder reduced model is shown in Figure 1 . This figure shows that the obtained reduced order model is an adequate loworder model that retains the characteristics of a full order model. Also, to show the efficiency of the proposed method, the step and frequency responses of the obtained reduced model are compared with those available in the literature. Figure 1 and Figure 2 show a comparison of the results obtained with the proposed method by Mukherjee [42] , Optimal Hankel norm approximation (HSV) [43] and Balanced Truncation (BT) [43] , respectively. It should be noted that a reduced order model is called a balanced truncation of the full order system, when a system is balanced and, then, the states corresponding to small Hankel singular values are discarded. Also, Optimal Hankel norm approximation finds a reduced order model, (of degree r ), such that the Hankel norm of the approximation error is minimized.

Figure 1. The step response of full order and reduced order model by the proposed method and other methods for test system 1.


Figure 2. The fequency response of full order and reduced order model by the proposed method and other methods for test system 1.

Figure 1 and Figure 2 show that the achieved results from the proposed method and the proposed method by Mukherjee are very similar to the original system, with respect to HSV and BT methods.
Furthermore, the specifications of the obtained system by the proposed method, such as steady state value, rise time, settling time and maximum overshoot, are compared with those obtained by Mukherjee [42] , HSV and BT, shown in Table 1 . Also, Integral Square Error (ISE) and the norm of the error between the step responses of full order and reduced order models ( ) are given in Table 1 . It is clearly seen that the specifications of the reduced order model that is achieved by the proposed method and the one by Mukherjee are close to the specifications of the original system.
Model  Steady state value  Over shoot (%)  Rise time (s)  Settling time (s)  ISE  Infinitynorm of error 

Full order  0.1  0  2.44  4.34  –  – 
ChebyshevHS  0.1  0  2.45  4.51  1.59×10−6  0.0067 
Mukherjee  0.1  0  2.45  4.5  1.34×10−6  0.0074 
BT  0.114  0  5.83  11.4  0.0379  0.0141 
HSV  0.114  0  3.63  7.07  0.0375  0.0139 
Also, the plot of is given for the full order and reduced systems in Figure 3 . This figure illustrates that the obtained error by the proposed method in this paper and the proposed method by Mukherjee is very similar, and less than HSV and BT methods.

Figure 3. The plot of for the full order and reduced systems by the proposed method and other methods for test system 1.

Test system 2: In [44] , a procedure is presented to obtain a reduced order system by Routh–Pade approximation using the Luus–Jaakola algorithm. To compare the proposed method in this paper with the Luus–Jaakola algorithm, the system given in [44] is adopted, which is a thirdorder system:

( 22) 
Based on the explanations given for test system 1, the obtained reduced system by the algorithm is as follows:

( 23) 
The step response of the original system and the obtained reduced model is shown in Figure 4 . In this figure, the responses of the system with secondorder primary reduced models obtained by other methods are also included for comparison. Also, the plot of is given for the full order and reduced systems in Figure 5 .

Figure 4. The step response of full order and reduced order model by the proposed method and other methods for test system 2.


Figure 5. The plot of for the full order and reduced systems by the proposed method and other methods for test system 2.

Once again, maximum overshoot, rise time, settling time, steady state value, ISE and norm of the error ( ) are given in Table 2 . It is clearly seen that the specifications of the reduced order model achieved by the proposed method are close to the specifications of the original system, and better than other methods.
Steady state value  Overshoot (%)  Rise time (s)  Settling time (s)  ISE  Infinitynorm of error  

Original model  1  86.5  0.129  6.74  –  – 
Chebyshev HS  0.99  87.9  0.127  2.35  0.0286  0.1139 
Luus–Jaakola  1  66.1  0.13  1.71  0.1404  0.3425 
BT  0.836  123  0.103  3.15  0.3802  0.1635 
HSV  0.836  115  0.118  3.44  0.4043  0.1635 
Test system 3: The third system to be reduced is a system given in [45] by Tao, where a procedure is presented to obtain a reduced system. is given in Box II .

( 24) 
Based on the explanations given for test system 1, the obtained reduced system by the algorithm is the one which is shown in Box III .

( 25) 
The comparison of the proposed method, the one proposed by Tao in [45] , HSV and BT methods are shown in Figure 6 and Figure 7 and Table 3 , which illustrate that the achieved results from the proposed method and the proposed method by Tao are very similar to the original system, with respect to HSV and BT methods.

Figure 6. Step response of full order and reduced order model by the proposed method and other methods for test system 3.


Figure 7. The plot of for the full order and reduced systems by the proposed method and other methods for test system 3.

Steady state value  Over shoot (%)  Rise time (s)  Settling time (s)  ISE  Infinitynorm of error  

Original system  −0.00281  −0.00814  0.0185  7.56  –  – 
ChebyshevHs  −0.00281  −0.0078  0.0164  9.13  5.8043×10−7  9.8248×10−4 
Tao method  −0.00281  −0.00738  0.0236  9.13  5.8235×10−7  7.5718×10−4 
BT  −0.00183  −0.00743  0.0146  7.31  1.4339×10−5  0.0011 
HSV  −0.00218  −0.00805  0.0719  7.18  5.9557×10−6  0.0010 
This paper introduces a new method for order reduction, using orthogonal polynomials, through Chebyshev rational functions and the harmony search algorithm. To get an optimum stable reduced system, Routh array is applied as constraints in the formulation. The proposed method is compared with some order reduction techniques, where the results obtained show that the proposed approach has high accuracy, and which results in an adequate loworder model that retains the characteristics of the full order model. The obtained results confirm that the proposed method can be used as an alternative method for order reduction, and the ability of other orthogonal functions can be investigated for the order reduction of the system.
Published on 06/10/16
Licence: Other
Are you one of the authors of this document?