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

Authors: Othman M. K. Alsmadi; Zaer S. Abo-Hammour
Journal: Computational Intelligence and Neuroscience
 (2015)

---

## Abstract

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 ofB, 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.


---

## Body

## 1. 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–5] and some for discrete-time systems [1, 11–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–5, 11] where important eigenvalues of the system are retained to find suitable lower order models.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 z-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: Section2 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.


## 2. 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(1)
y
k
+
a
1
y
k
-
1
+
a
2
y
k
-
2
+
⋯
+
a
n
y
k
-
n
=
b
0
u
k
+
b
1
u
k
-
1
+
⋯
+
b
n
-
u
k
-
n
-
,where u
(
k
) is the input and y
(
k
) is the output of the system at the kth sampling instant. Equation (1) can be written in the form of a pulse transfer function as(2)
G
(
z
)
=
Y
(
z
)
U
(
z
)
=
b
0
z
n
-
+
b
1
z
n
-
-
1
+
⋯
+
b
n
-
z
n
+
a
1
z
n
-
1
+
⋯
+
a
nwith n
-
≤
n. The characteristic polynomial contains the system dominant and nondominant poles (distinct, repeated, or complex) where their number, n, is referred to as the model order. The corresponding desired reduced rth order model is given by(3)
G
r
z
=
Y
z
U
z
=
b
0
z
r
-
+
b
1
z
r
-
-
1
+
⋯
+
b
r
-
z
r
+
a
1
z
r
-
1
+
⋯
+
a
r
,where some of the coefficients a
i
(
i
=
1,2
,
…
,
r
) and b
i
(
i
=
0,1
,
2
,
…
,
r
-
) may be zeros as long as r
-
≤
r. For the MIMO systems, consider the following nth order discrete-time system:(4)
x
k
+
1
=
A
x
k
+
B
u
k
,
(5)
y
k
=
C
x
k
+
D
u
k
,where k is the time index, x
∈
R
n is the state vector, u
∈
R
p and y
∈
R
m are the input and output vectors, respectively, and A
∈
R
n
×
n, B
∈
R
n
×
p, C
∈
R
m
×
n, D
∈
R
m
×
p are matrices of appropriate dimensions with n, p, and m being the system order, number of inputs, and number of outputs, respectively. The corresponding desired reduced rth order model is obtained as follows: (6)
x
r
k
+
1
=
A
r
x
r
k
+
B
r
u
k
,
(7)
y
r
k
=
C
r
x
r
k
+
D
r
u
k
,where x
r
(
k
) is r-state vector, y
r
(
k
) is the reduced order model output, and A
r, B
r, C
r, and D
r are matrices with appropriate dimensions.


## 3. 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 matrixA
r 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 A
r, in (6), is designed to have the following decoupling format:(8)
A
r
=
λ
1
a
12
a
13
a
14
⋯
a
1
r
0
λ
2
a
23
a
24
⋯
⋱
⋱
⋮
·
0
λ
b
-
⋮
·
0
σ
1
α
1
·
0
-
α
1
σ
1
0
0
⋱
a
(
r
-
2
)
r
0
σ
p
-
α
p
-
0
·
·
·
0
-
α
p
-
σ
p
-
,where the original system dominant poles (real and/or complex) are preserved in the diagonal, seen as λ
i, i
=
1,2
,
…
,
b
- (real) and σ
i
±
α
i, i
=
1,2
,
…
,
p
- (complex). Notice that, for this reduced order model, r
=
(
b
-
+
2
p
-
)
<
n. To insure that the dominant poles are preserved in the reduced order model and for further order reduction, the following condition is satisfied:(9)
λ
dominant
:
=
λ
1
>
λ
2
>
⋯
>
λ
b
-
,
λ
>
σ
,
σ
1
±
α
1
>
σ
2
±
α
2
>
⋯
>
σ
p
-
±
α
p
.Take into account that if λ
i
<
σ
i
±
α
i, 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 a
i
j (i and j
=
1
,
…
,
r) are set to zero.The GA will determine the parameters of theB
r and C
r (and D
r if necessary) in (6) and (7). Hence, the total number of elements that the GA will need to find is given by (10)
n
p
=
n
r
·
m
+
p
·
n
r
.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 (n
p), 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 n
p individuals. Given a population size (n
pop), a matrix consisting of n
pop rows is formed with each row containing one set of solutions for the unknown parameter values. This would result in a matrix containing n
pop × n
p 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.


### 3.1. Initialization

An initial population comprisingN
p 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 (p
l), upper parameter bound (p
u), and accuracy (ϕ). Hence, the number of bits used is defined as follows [30]:(11)
n
bits
=
lo
g
10
p
u
-
p
l
/
ϕ
lo
g
10
2
.Given that each parameter value will consist of n
bits number of bits, each solution set will consist of n
b
,
t
=
n
p
·
n
bits bits. This represents one row of the entire n
pop × n
b
,
t matrix.The GA starts by randomly initializing a binary matrix withn
pop rows and n
b
,
t 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.Figure 1
A set of individuals in a GA chromosome.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.


### 3.2. 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:(12)
e
(
k
)
=
y
(
k
)
-
y
r
(
k
)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]:(13)
Fitness
=
1
1
+
MSE
×
100
%
,where (14)
MSE
=
1
N
∑
i
=
1
N
y
i
-
y
r
i
2
,where MSE is the mean-square-error and N 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.


### 3.3. 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 (X
r) 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 X
r = 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 Figure2.Figure 2
Chromosomes random crossover operation.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.


### 3.5. 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.


### 3.6. 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.


### 3.7. 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 theB
r, C
r, and D
r matrices to produce the reduced order model.


## 3.1. Initialization

An initial population comprisingN
p 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 (p
l), upper parameter bound (p
u), and accuracy (ϕ). Hence, the number of bits used is defined as follows [30]:(11)
n
bits
=
lo
g
10
p
u
-
p
l
/
ϕ
lo
g
10
2
.Given that each parameter value will consist of n
bits number of bits, each solution set will consist of n
b
,
t
=
n
p
·
n
bits bits. This represents one row of the entire n
pop × n
b
,
t matrix.The GA starts by randomly initializing a binary matrix withn
pop rows and n
b
,
t 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.Figure 1
A set of individuals in a GA chromosome.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.


## 3.2. 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:(12)
e
(
k
)
=
y
(
k
)
-
y
r
(
k
)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]:(13)
Fitness
=
1
1
+
MSE
×
100
%
,where (14)
MSE
=
1
N
∑
i
=
1
N
y
i
-
y
r
i
2
,where MSE is the mean-square-error and N 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.


## 3.3. 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 (X
r) 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 X
r = 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 Figure2.Figure 2
Chromosomes random crossover operation.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.


## 3.5. 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.


## 3.6. 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.


## 3.7. 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 theB
r, C
r, and D
r matrices to produce the reduced order model.


## 4. Illustrative Examples

In this section, we consider a few discrete system examples, which have been investigated by different researchers in recent years. The examples are given as SISO and MIMO linear time invariant systems. Results of investigations are compared with the existing ones to show the potential of the new approach.Example 1.
As a first example, we consider the following 7th order SISO system investigated by Telescu et al. [32] given as(15)
G
z
=
2.0434
z
7
-
4.98255
z
6
+
6.57
z
5
-
5.8189
z
4
+
3.636
z
3
-
0.00088
z
2
-
1.4105
z
+
0.2997
·
z
7
-
2.4
6
6
+
3.433
z
5
-
3.333
z
4
+
2.546
z
3
-
1.584
z
2
+
0.7478
z
-
0.252
-
1with system dynamics (poles): [0.8913, 0.6843
±
0.5820
i, 0.2988
±
0.7574
i, -
0.1987
±
0.6993
i]. Telescu et al. [32] used the Laguerre functions for their proposed method and obtained a 5th order reduced model:(16)
G
r
(
z
)
=
2.043
z
5
-
3.057
z
4
+
2.195
z
3
-
1.545
z
2
+
0.8617
z
z
5
-
1.518
z
4
+
1.270
z
3
-
1.032
z
2
+
0.7539
z
-
0.3156with poles given as [0.8320, −0.2318
±
0.7612
i, 0.5748
±
0.5183
i], which are unrelated to the original system dynamics. On the other hand, our first advantage of the proposed technique is that a lower dimension (3rd order) than the Telescu reduced model was obtained:(17)
G
r
(
z
)
=
1.51511
z
3
-
3.12377
z
2
+
2.40413
z
-
0.63117
z
3
-
2.2598
z
2
+
2.0266
z
-
0.71915with poles given as [0.8913, 0.6843
±
0.5820
i], 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.Figure 3
Impulse responses for the full and reduced order models.Our fourth advantage is clearly seen when simulating the full and reduced models to a step input. As seen in Figure4, 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.Figure 4
Step responses for the full and reduced order models.Example 2.
In this example, we consider the following 8th order SISO system investigated by Yadav et al. [29]:(18)
G
z
=
0.1625
z
7
+
0.125
z
6
-
0.0025
z
5
+
0.00525
z
4
-
0.02263
z
3
-
0.00088
z
2
+
0.003
z
-
0.000413
·
z
8
-
0.6307
z
7
-
0.4185
z
6
+
0.078
z
5
-
0.057
z
4
+
0.1935
z
3
+
0.09825
z
2
-
0.0165
z
+
0.00225
-
1with system dynamics: [0.8797
±
0.2442
i, −0.0542
±
0.6558
i, −0.5875
±
0.0959
i, 0.0773
±
0.1078
i]. Using the proposed method, the dominant dynamics [0.8797
±
0.2442
i] are preserved in the reduced order model which is given by(19)
G
r
z
=
0.021129
z
2
+
0.12259
z
-
0.063666
z
2
-
1.7594
z
+
0.83347
.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.Table 1
Model order reduction method comparison.
System
Dominant dynamics
Original 8th order
0.8797 ± 0.2442i, −0.5875 ± 0.0959i
0.0773 ± 0.1078i, −0.0542 ± 0.6558i
Proposed 2nd order
0.8797 ± 0.2442i
Artificial bee colony [24]  2nd order
0.8778 ± 0.2432i
Differential evolution [28] 2nd order
0.8768 ± 0.2484i
Genetic algorithm [29] 2nd order
0.8854 ± 0.2398iAs a second comparison, the full and reduced order models were simulated for a step input with results as seen in Figure5. 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.Figure 5
Step responses for the 8th and reduced 2nd order models.Figure 6
Error comparison of the reduced order models with the full order.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:(20)
x
k
+
1
=
0.5034
0.1768
-
0.2340
-
0.1406
0.5814
0.0096
0.5498
-
0.0362
-
0.6744
2.2496
0.0337
0.2546
0.0984
-
0.4051
1.3599
-
0.2709
0.1470
0.3249
0.0484
0.6356
-
0.0909
0.0491
0.1075
-
0.1019
0.5681
x
(
k
)
+
0.3306
0.1700
0.8951
0.3442
0.5487
0.2143
0.8748
0.8821
0.5217
0.4479
u
k
,
y
k
=
3.0622
-
0.9986
-
0.7126
6.4339
-
10.4291
3.0396
-
0.9913
-
0.7073
5.2369
-
8.4887
x
k
+
0
0
0
0
u
k
.This system is given with the dynamics [0.5614, 0.5008, 0.3389, 0.2295, 0.1375]. Li performed the Coprime Factor MOR and obtained the following 4th order reduced model:(21)
x
r
k
+
1
=
0.3830
0.0865
-
0.3771
0.7548
-
0.1633
0.6255
-
0.0157
0.1556
-
0.6731
0.5150
0.3777
-
0.2194
-
0.3692
0.2832
0.0094
0.2388
x
r
k
+
1.0572
0.0349
-
0.2642
-
0.5105
-
0.6643
-
1.5790
0.9864
0.0412
u
k
,
y
k
=
-
0.2745
-
0.0669
-
0.6423
-
0.2425
-
0.2725
-
0.0664
-
0.5228
-
0.1973
x
r
k
+
0
0
0
0
u
(
k
)with poles [0.5567, 0.5004, 0.3385, 0.2294], which are close to the full order model, but not the exact values. Using the proposed approach, the following 4th order reduced model was obtained:(22)
x
r
k
+
1
=
0.5614
0
0
0
0
0.5008
0
0
0
0
0.3389
0
0
0
0
0.2295
x
r
(
k
)
+
0.4285
0.2628
0.4157
0.2551
0.3405
0.2348
0.3111
0.2671
u
k
,
y
k
=
29.7314
-
20.0311
-
37.9382
27.0723
26.5158
-
20.4912
-
28.5003
21.6793
x
r
k
+
1
0
-
3
0.0061
-
0.1580
0.0046
-
0.1191
u
(
k
)with poles [0.5614, 0.5008, 0.3389, 0.2295], which are the exact dominant dynamics of the original model. Hence, the superiority of the proposed method is clearly seen. Next, we investigate for quality performance of the new method compared with the Li approached reduced order 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.Figure 7
Full and reduced order model’s responses for mixed input signals.Figure 8
Error of the reduced order models for the system first output response.Figure 9
Error of the reduced order models for the system second output response.In Figures8 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.


## 5. 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 theB, C, and D 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.



---

*Source: 615079-2015-03-08.xml*