A Robust Computational Technique for Model Order Reduction of Two-Time-Scale Discrete Systems via Genetic Algorithms

A robust computational technique for model order reduction (MOR) of multi-time-scale discrete systems (single input single output (SISO) and multi-input multioutput (MIMO)) is presented in this paper. This work is motivated by the singular perturbation of multi-time-scale systems where some specific dynamics may not have significant influence on the overall system behavior. The new approach is proposed using genetic algorithms (GA) with the advantage of obtaining a reduced order model, maintaining the exact dominant dynamics in the reduced order, and minimizing the steady state error. The reduction process is performed by obtaining an upper triangular transformed matrix of the system state matrix defined in state space representation along with the elements of B, C, and D matrices. The GA computational procedure is based on maximizing the fitness function corresponding to the response deviation between the full and reduced order models. The proposed computational intelligence MOR method is compared to recently published work on MOR techniques where simulation results show the potential and advantages of the new approach.


Introduction
Model order reduction (MOR) of multi-time-scale systems has been an important subject area in control engineering for many years [1,2]. In many industrial control systems, simple controllers are preferable. However, derivation of the mathematical model often leads to detailed description of a complex model in the form of high order differential equations [2]. Due to this point of view along with other different design objectives, model order reduction has been an active research area in the control society since the 1960s where a large number of model order reduction methods have been introduced in literature for single input single output (SISO) as well as MIMO type systems. The reduction operation is to search for the possibility of finding some lower order equations of the same type that may be considered to adequately reflect the dominant characteristics of the original system. The objective of simplification is to obtain a low order model of the existing high order model such that both are equivalent in terms of system response and being close to each other in some physical representation means. Model reduction problems have attracted much attention in recent years; for example, the model reduction problem has been investigated using artificial neural networks [3], genetic algorithms [4], and invasive weed optimization [5]. It was also used in nonlinear systems [6], gain scheduling [7], linear time-varying systems [8,9], and linear parameter-varying systems [10].
To obtain a model of lower order, a significant number of methods have been proposed in recent and earlier years, some for continuous time systems [3][4][5] and some for discretetime systems [1,[11][12][13][14][15]. Some methods, such as model order reduction by matching Markov parameters [16], were introduced to ensure stability of the reduced order model. A popular technique for obtaining reduced order models is the Krylov subspace [17]; however, stability of the reduced model is not guaranteed. Another important group of reduction algorithm is the eigenvalue preservation technique [3][4][5]11] where important eigenvalues of the system are retained to find suitable lower order models.

Computational Intelligence and Neuroscience
A numerous number of MOR methods are available for continuous systems, but very few have been devoted to the discrete-time systems MOR. The discrete-time system MOR may be performed in two different ways. The first one is performed based on transforming a continuous time model into another form using different types of transformation as seen in [6,18]. In this group of MOR, the reduction process is completely performed in the continuous time form. The discrete reduced order model is then obtained by the corresponding inverse transformation of the continuous time reduced model. The second method for obtaining a discrete reduced order model, which is known as a direct method [14], is deriving the discrete reduced order model directly without using any type of transformation. Some of these methods perform the MOR using canonical expansion of -transfer function and stable optimal methods [13,19], power decomposition and system identification [20], and multipoint step response matching [21]. New optimization techniques, particle swarm optimization [22], and artificial neural networks [11] have also been introduced for MOR of discrete-time systems.
GA-based MOR, on the other hand, has received some of the researchers' attention as well. Recently, Ponda et al. [23] employed a particular swarm optimization technique to obtain a reduced order model of SISO large scale linear systems. Their technique is based on the integral square error (ISE). Vishwakarma and Prasad [24] proposed a mixed method for reducing the order of large-scale linear systems. They have synthesized the denominator of the reduced order transfer function using modified pole clustering while the coefficients of the numerator elements are computed using GA. Parmar et al. [25] presented a technique for model order reduction using GA for SISO linear time systems. They have focused on obtaining a reduced order model that maintains stability and retains the steady state value. In spite of the methods available in literature, each method has advantages and disadvantages when tried on a particular system. In addition to that, no approach always gives the best results for all systems. It is important to mention that GAs have also been used for model system identification, where order and parameters are set to be determined, as we have investigated in [26,27]. In this paper, however, and as motivated by the singular perturbation method which has the characterization of multi-time-scale systems, the GA procedure is performed with the advantages of retaining the exact dominant dynamics in the designed model, obtaining a new robust model with a lower order, and maintaining a minimum steady state response error.
The work in this paper is organized as follows: Section 2 presents problem formulation of the discrete full and reduced order models. In Section 3, the genetic algorithm approach for MOR of multi-time-scale discrete systems is presented. Illustrative examples utilizing the new approach along with simulation comparative results of different MOR techniques are presented in Section 4. Section 5 presents an overall conclusion of the proposed MOR method.

Problem Formulation
In this paper, MOR is investigated for discrete LTI systems of both SISO and MIMO type models. For SISO systems, a transfer function model is used, while the state space representation is used for MIMO systems.
For the SISO systems, consider the discrete-time system described by where ( ) is the input and ( ) is the output of the system at the th sampling instant. Equation (1) can be written in the form of a pulse transfer function as with ≤ . The characteristic polynomial contains the system dominant and nondominant poles (distinct, repeated, or complex) where their number, , is referred to as the model order. The corresponding desired reduced th order model is given by where some of the coefficients ( = 1, 2, . . . , ) and ( = 0, 1, 2, . . . , ) may be zeros as long as ≤ . For the MIMO systems, consider the following th order discretetime system: where is the time index, ∈ R is the state vector, ∈ R and ∈ R are the input and output vectors, respectively, and ∈ R × , ∈ R × , ∈ R × , ∈ R × are matrices of appropriate dimensions with , , and being the system order, number of inputs, and number of outputs, respectively. The corresponding desired reduced th order model is obtained as follows: where ( ) is -state vector, ( ) is the reduced order model output, and , , , and are matrices with appropriate dimensions.

Genetic Algorithms with MOR
GAs are based on principles inspired from the genetic and evolution mechanisms observed in natural systems. Their basic principle is the maintenance of a population of solutions to the problem that evolves in time. They are based on the triangle of genetic reproduction, evaluation, and selection [12].
Genetic reproduction is performed by means of two basic genetic operators: crossover and mutation. Evaluation is performed by means of the fitness function that depends on the specific problem. Selection is the mechanism that selects parent individuals with probability proportional to their relative fitness.
In this paper, using the computational intelligence of GA, we will obtain the reduced order model based on only the dominant dynamics of the system. For dynamic decoupling, the GA will set the reduced order model state matrix in the modal form where all of the selected eigenvalues (system dynamics), real and/or complex, are placed on the diagonal. Thus, the reduced order model state matrix , in (6), is designed to have the following decoupling format: where the original system dominant poles (real and/or complex) are preserved in the diagonal, seen as , = 1, 2, . . . , (real) and ± , = 1, 2, . . . , (complex). Notice that, for this reduced order model, = ( + 2 ) < . To insure that the dominant poles are preserved in the reduced order model and for further order reduction, the following condition is satisfied: Take into account that if | | < | ± |, then (9) is to be redefined accordingly if necessary. For simplicity, the modal form is chosen, which implies that all elements seen in (8) as ( and = 1, . . . , ) are set to zero. The GA will determine the parameters of the and (and if necessary) in (6) and (7). Hence, the total number of elements that the GA will need to find is given by It is to be noted that all of the parameters that the GA will have to find are restricted to be real values. Now, based on the number of unknown parameters ( ), the GA creates a population of individuals, where each parameter is basically an individual in this population. The population consists of different "sets" of solutions. Each solution set is called a chromosome, which contains individuals. Given a population size ( pop ), a matrix consisting of pop rows is formed with each row containing one set of solutions for the unknown parameter values. This would result in a matrix containing pop × elements. Each element in this matrix contains a value pertaining to one unknown parameter, where each row presents one set of solutions.
The genetic algorithm used in this work will operate as follows.

Initialization. An initial population comprising
individuals is randomly generated. The GA type used in this paper is a binary genetic algorithm, where each value in the solution set consists of a number of bits (genes). The number of bits used to encode each numerical value depends on three variables, lower parameter bound ( ), upper parameter bound ( ), and accuracy ( ). Hence, the number of bits used is defined as follows [30]: Given that each parameter value will consist of bits number of bits, each solution set will consist of , = ⋅ bits bits. This represents one row of the entire pop × , matrix.
The GA starts by randomly initializing a binary matrix with pop rows and , columns. Each row (set of solutions) is made up of multiple values decoded into binary and placed next to each other as illustrated in Figure 1 for one chromosome. This chromosome consists of three parameters (individuals) with each individual being made of three genes. These genes (bits) can be later decoded back into decimal values, which in return provide the desired parameters' values.
The initial population matrix consists of random numbers within the lower and the upper bounds of the parameters. The population will be split into rows (chromosomes), each constituting one solution set. These solution sets are each taken to have their fitness evaluated, as seen next.

Evaluation.
The fitness, a nonnegative measure of quality, is used as a measure to reflect the degree of goodness of the individual and is calculated for each individual in the population. This measure of quality is calculated based on minimizing the following related cost function: which is the deviation between the full and reduced order models' responses given in (5) and (7). The fitness for each chromosome (solution set) is then evaluated as follows [31]: where 4 Computational Intelligence and Neuroscience  where MSE is the mean-square-error and is the number of elements in the output vector(s). The higher this fitness is, the closer the reduced model is to mimic the original model.

Selection.
In the selection process, chromosomes are chosen from the current population to enter a mating pool devoted to the creation of new children (offspring) for the next generation. The chance of a given individual to be selected to mate is proportional to its relative fitness That is, the larger the fitness value of a chromosome, the higher the probability of the chromosome to contribute one or more children in the next generation. First, the population fitness and associated chromosomes are ranked from highest fitness to lowest fitness. Then, only the best are selected to continue, while the rest are left. The selection rate ( ) is a rank based ratio, which is the percentage of the best rank individual that should move to the mating pool where pairs of the mating pool are selected for crossover process. Deciding how many chromosomes to keep is somewhat arbitrary. Letting only a few chromosomes survive to the next generation limits the available genes in the offspring (children). Keeping too many chromosomes gives bad performers a chance to contribute their traits to the next generation. In the proposed algorithm, an initial selection rate = 0.60 was used and the top 60% fitness chromosomes were placed in a mating pool.
There are several methods for choosing the chromosome pairs to be mated from the mating pool. In the proposed algorithm, random pairing was chosen, which uses a uniform random number generator to select chromosomes that enter the mating pool from kept chromosomes. The algorithm randomly chooses chromosomes from the mating pool for mating, while making sure that all pairs are unique. After the chromosome pairs are chosen, the next operation will be crossover.

3.4.
Crossover. Crossover provides the means by which valuable information is shared among the population. In the crossover operation, a pair of chromosomes is mated to produce two offspring in the process that inherit genes from their parents. To perform crossover, the chromosomes need to be in their gene format (i.e., binary representation). The two parent chromosomes are crossed over at random points to segment each chromosome into a number of parts. These parts of each chromosome pair are swapped between each other to produce two new chromosomes, the offspring. This is shown in Figure 2.
Crossover is not applied to all pairs of chromosomes selected for mating. For every parent chromosome pair, a crossover rate decides whether or not crossover occurs. If crossover is not applied, children are produced simply by duplicating the parents. This gives each chromosome a chance of passing on its genes without the disruption of crossover. In the proposed algorithm, the crossover rate is set to 0.80, meaning that 80% of parent pairs produce offspring. The number of crossover points for the parent chromosomes is chosen as a random integer for each generation. As a crossover result, a new children population will be performed with the same size as the parent population.

Mutation.
Mutation is often introduced to guard against premature convergence. Generally, over a period of several generations, the gene pool tends to become more and more homogeneous. Therefore, further mutation is introduced to the offspring to guarantee that the offspring will have new qualities while retaining a similar overall structure. Mutation is applied to each child individually after crossover, where it randomly flips any bit (gene) with a small mutation probability (between 0.1 and 0.001). Mutation provides a small amount of random search to guard against premature convergence. In the proposed methodology, the mutation rate is set to 0.0125.

Replacement. Replacement operation takes place once the crossover of all parent chromosome pairs is performed.
The parent population is totally or partially replaced by the children population depending on the replacement scheme used. This completes the life cycle of the population. At this stage, the population is ready to enter the next generation and undergo a new round genetic operation. The decision in which offspring replaces which population individual is made based on the fitness evaluation. Fitness is evaluated for all the resulting offspring as well as for all of the original population individuals and a replacement factor decides how many of the offspring with the highest fitness values are to replace the main population individuals. In the proposed approach, a replacement factor of 0.70 is used. This process continues until reaching the end of the population size or best desired fitness.

Termination.
The GA is terminated when some convergence criterion is met. The termination condition could be considered as specified fitness value, reaching maximum number of generations, or a set progress limit. As one generation has gone by, depending on the genetic algorithm's termination condition, the algorithm could stop any time and identifies the chromosome with the highest fitness value as the optimal solution set. On the other hand, it may repeat the entire process from the selection procedure to continue another generation of the genetic algorithm. The termination condition for the proposed algorithm is reaching a fitness value of 99.999% and an average fitness value of at least 99.9% in the entire population set. When the algorithm terminates, the highest fitness chromosome is distributed across the , , and matrices to produce the reduced order model.

Illustrative Examples
In this section, we consider a with poles given as [0.8913, 0.6843 ± 0.5820 ], which are exactly the dominant dynamics of the full order system as seen above, and this is our second advantage. The result was obtained for a population of 1000 and 300 generations with a fitness of 99.992. To investigate the system behavior, the full and reduced order models were excited by an impulse input (as performed by Telescu) with results of simulation shown in Figure 3, which shows our third advantage as the propose 3rd order reduced model response converges in about one second, which is much faster than the Telescu's response.
Our fourth advantage is clearly seen when simulating the full and reduced models to a step input. As seen in Figure 4, Telescu et al. [32] 5th order model has a huge steady state error while the proposed 3rd order model's error can barely be seen at the steady state which shows the robustness of the proposed method.
For result evaluation, the proposed reduced order model was compared with recent research for MOR, that is, artificial bee colony [28], differential evolution optimization algorithm, and real coded genetic algorithm [33]. Hence, we first compare for dominant dynamics of the reduced order models with results obtained as in Table 1. As can be seen, the proposed method provides the dominant dynamics of the full order model retained exactly in the reduced order, while in the other methods, the reduced order dynamics are close to the full order model, but not as close as our proposed method's results.
As a second comparison, the full and reduced order models were simulated for a step input with results as seen in Figure 5. Observing the results in Figure 5, we can see that the 2nd order reduced model of all methods seems to be close to the original full order model. However, when we take a closer look, we see the differences very clearly as present in Figure 6. In this figure, it is seen that the proposed approach provides the least error among all. The error is seen a little at the beginning and then, at about one second, the error becomes very close to zero.
Example 3. In this example, we consider the 5th order MIMO discrete system investigated by Li [2] as given by the following state space model:   Simulating the full and reduced order models for a mixed type signal, the two system output responses are presented in Figure 7. As seen in this figure, the two reduced order models seem to perform with same quality. Therefore, in order to observe the differences between the proposed and existed methods, we present the response accuracy error, as seen in Figures 8 and 9 for the two output responses. In Figures 8 and 9, notice that the Coprime Factor method response errors are much far from zero when compared with the proposed method response errors. The response errors of the proposed GA approach can barely be seen as they are almost zeros, which clearly shows the superiority of the proposed approach.
As seen in the three previous examples, the advantages of the new approach are present in all three of them. That is, the reduction process provides new reduced order models, maintains the exact dominant dynamics of the original model in the reduced order model, and provides reduced order models with the least response errors as compared with other methods. In addition to that, an advantage that the new approach is applicable to SISO and MIMO discrete type systems.  Figure 9: Error of the reduced order models for the system second output response.

Conclusion
In this paper, a robust computational intelligence approach is proposed using GA for MOR of discrete MIMO type systems with the advantage of dominant dynamic preservation. The reduction process is performed based on transforming the system state matrix (in a state space model) while decoupling the multi-time-scale dynamics. The dominant dynamic preservation is performed by retaining the dominant poles of the original system as a subset in the reduced order model utilizing the transformed system state matrix. Once the reduced order model state matrix is designed, the GA MOR new technique will search for the , , and matrices of the reduced order model. Different examples are presented with a comparison of recently published work, such as differential evolution, artificial bee colony, and the Coprime Factor MOR methods. Comparison results show the robustness of the proposed method, as it outperforms such techniques, which points out the potential of the new approach.