# Combining a Genetic Algorithm and Simulated Annealing to Design a Fixed-Order MixedH2/H∞ Deconvolution Filter with Missing Observations

Authors: Jui-Chung Hung
Journal: Journal of Control Science and Engineering
 (2008)

---

## Abstract

We introduce a new combination approach to a fixed-order mixedH2/H∞ deconvolution filter with missing observations. The missing observations model is based on a probabilistic structure with the probability of the occurrence of missing data modeled as the unknown prior. The aim of the mixed H2/H∞ criterion is to achieve H2 optimal reconstruction and subject the H∞ norm constraint to the transfer function from the channel input to the filter error. For simplicity of implementation, the fixed-order model is interesting for engineers in signal processing and in practical applications. In this situation, the deconvolution filter design becomes a complicated nonlinear estimation problem. In this paper, we combine a genetic algorithm (GA) and simulated annealing (SA) to treat the signal reconstruction problem with missing observations. Finally, a numerical example is presented to illustrate the design procedure and confirm the robustness performance of the proposed method.


---

## Body

## 1. Introduction

A deconvolution filter is used to remove the distortion of a channel and suppress the influence of additive noise in the transmission path. This problem has a wide range of applications in the field of engineering, such as equalization, seismology, and image restoration. When the signal is not corrupted by noise, the channel model is minimum phase, and none of the received data is missing; the inverse system is the simplest deconvolution filter. However, most of the deconvolution filters are not as easy as mentioned above since the signal is usually corrupted by noise; the channel may be nonminimun phase, and the received data may be missing. Therefore, various types of optimal deconvolution filters have been developed by minimizing different criteria [1–3].The deconvolution problems generally have one of two frameworks: one isH2, the other is H∞. Most of the deconvolution problems are solved via the H2 optimal method in the time [4], and frequency domains [5]. When the channels are perturbed, the H2 optimal filter is inadequate. A robust deconvolution filter design based on H∞ theory has received great attention for its robustness against system uncertainties [5–7]. However, the H∞ robustness design cannot achieve optimal performance. The proposed mixed H2/H∞ optimal deconvolution filter design minimizes H2 reconstruction performance subject to a robustness requirement based on the H∞ norm to attenuate the performance degradation due to the system's uncertainty [8]. This problem can be interpreted as a problem of optimal reconstruction filter design subject to a robustness constraint against the deterioration due to parameter variation in the channel.In practical situations, the received data is not consecutive and contains missing observations [9]. The missing observations may be uniform or randomly occur and are caused by a variety of reasons, for example, intermittent sensor measurement failures, loss of collected data, jammed data, or data coming from a high-noise environment. If not properly taken into account, the missing observations can seriously deteriorate the quality of the optimal deconvolution filter. The design problem for a mixed H2/H∞ deconvolution filter with missing observations is more difficult than the case with no missing observations. In general, the pattern of misses can be quite arbitrary. In this study, the missing observations are modeled as random Bernoulli pattern, in which each measurement has a fixed probability of being missing, and the misses are independent.Furthermore, if the order of signal models and the channel dynamics are higher, they may lead to complexities in the deconvolution filter structure. In practical applications, fixed-order mixedH2/H∞ optimal deconvolution filter designs are appealing for their simplicity of implementation and operation time savings. Thus, the fixed-order mixed H2/H∞ deconvolution filter [5] is used; conventional optimization techniques cannot be employed to obtain a closed-form solution of the mixed H2/H∞ optimal deconvolution filter with missing observations. In this study, the fixed-order mixed H2/H∞ deconvolution filter design procedure is divided into two steps. The first step is based on a GA to search the neighborhood of the coefficients of the optimal deconvolution filter. In the second step, SA is developed to search the coefficients of the optimal deconvolution filter estimation.The GA was introduced as an optimal search algorithm [10–14]. It is a parallel global search technique that emulates natural genetic operators such as reproduction, crossover, and mutation. At each generation, it explores different areas of the parameter space, and then directs the search to the region where there is a high probability of finding improved performance. Because the genetic algorithm simultaneously evaluates many points in parameter space, it is more likely to converge toward the global solution. In addition, it need not assume that the search space is differentiable or continuous, and it can also iterate several times on each piece of received data.The SA algorithm [15–19] is a numerical simulation method based on the dynamics of crystallization. A solid may be heated until its constituents can move freely and it melts. Then, the melt is allowed to cool very slowly until it solidifies in a certain arrangement. The SA algorithm makes use of this idea. Essentially, the SA strategy is a serial search method, unlike the GA, which is a parallel search method. S. Geman and D. Geman [19] proved that if the cooling schedule is slow enough, the desired global extreme can always be found. Comparing GA and SA, we can find that GA exhibits fast initial convergence, but its performance deteriorates as it approaches the desired global extreme. Interestingly, SA shows a complementary convergence pattern, in addition to high accuracy. We combine the selected features from GA and SA to achieve weak dependence on initial parameters, parallel search strategy, fast convergence, and high accuracy [20, 21]. The GA/SA starts the search procedure as a pure-GA and ends as a pure-SA. The transition from GA to SA occurs when the fittest individual remains the same for Lg generations. Hence, GA/SA is very suitable to treat the global optimization problem of the nonlinear parameter estimation with corrupting noise and missing observations. Finally, a numerical example is given to illustrate the design procedure of the proposed design method and to confirm the robustness performance of the fixed-order mixed H2/H∞ deconvolution filter with missing observations.


## 2. Fixed-Order MixedH2/H∞ Deconvolution Filter under Missing Observations

Consider a discrete deconvolution system with missing observations, as shown in Figure1. The received signal ym(k) is(1)ym=g(k)y(k)+v(k)=g(k)[x(k)+r(k)]+v(k)=g(k)[H(z-1)C(z-1)ω(k)+N(z-1)n(k)]+v(k), where z-1 is the inverse of z (i.e., unit delay), v(k), ω(k), and n(k) are assumed to be zero-mean white Gaussian noises with the following covariances:(2)E[v2(k)]=σv2,E[ω2(k)]=σw2,E[n2(k)]=σn2,E[ω(k)v(k)]=0,E[ω(k)n(k)]=0,E[n(k)v(k)]=0, where E denotes the expectation operator, σv2, σω2, and σn2 are assumed to be positive scales, and the system is assumed to have reached a statistical steady state.Figure 1
Deconvolution filter system model with irregular missing observations.Theg(k) is the model of missing observations:(3)g(k)={0,ify(k)ismissing,1,otherwise.Thus,ym can be regarded as the measurements of y(k) with missing observations. The sequence g(k) is assumed to be asymptotically stationary and independent of y(k). Furthermore, they are mutually independent. The probability of a missing measurement is [9](4)E[g(k)]=Pr[g(k)=1]=p,E[g(k)g(t)]={p,ifk=t,p2,ifk≠t, where p is a fixed probability, independent of time. We let u(k) denote the estimate of u(k). From Figure 1, the estimation error is given by(5)e(k)=u(k)-u^(k)=[1-g(k)D(z-1)H(z-1)]×C(z-1)ω(k)-g(k)D(z-1)×N(z-1)n(k)-D(z-1)v(k),where D(z-1) is the fixed-order deconvolution filter. The objective of this design is to find a fixed-order IIR deconvolution filter: (6)D(z-1)=b0+b1z-1+⋯+biz-i+⋯+bnz-n1+a1z-1+⋯+aiz-i+⋯+anz-n.The power spectral density ofe(k) is given by(7)Φ(z-1)=p[1-D(z-1)H(z-1)]C(z-1)σω2C*(z-1)×[1-D(z-1)H(z-1)]*+p⋅D(z-1)N(z-1)×σn2N*(z-1)D*(z-1)+D(z-1)σv2D*(z-1), where superscript * denotes the complex conjugate and Φ(z-1) is the power spectral density of e(k).In the fixed-orderH2 optimal deconvolution filter design case, a stable filter D(z-1) will be specified to minimize the mean square error (MMSE) [4]. The following expression for the performance index is derived according to Parseval's Theorem [22]:(8)minD(z-1),pE{e2(k)}=minD(z-1),p12πi∮|z|=1‍Φ(z-1)dzz, where Φ(z-1) is defined in (7), ∮‍ denotes integration around the unit circle. We assume the Laurent expansion of Φ(z-1)/z about the singular point, z-1=zj for j=1,2,…,n. By the Residue theorem, the H2 optimal deconvolution problem in (8) becomes the following minimization problem:(9)minD(z-1),pI2=min∑j=1nRezj‍s(Φ(z-1)z).In practical deconvolution systems, the noise variancesσv2, σω2, and σn2 may vary and the observations may be missing. Eliminating the performance degradation due to noise uncertainties and channel perturbation to guarantee the reconstruction performance is an important topic in practical signal deconvolution problems. The H∞ robustness design can guarantee the worst-case effect of these noise uncertainties or channel perturbation on the reconstruction performance to be less than a prescribed level [7], that is, if the H∞ norm of error spectrum Φ(z-1) is less than ε, the sensitivity of reconstruction error e(k) to the noise uncertainties or channel perturbation must be less than ε from the energy point of view. In this study, in order to take advantage of both H2 optimal reconstruction and less sensitivity to noise uncertainties and channel perturbation, a fixed-order deconvolution filter D(z-1) is specified to achieve the H2 minimization of reconstruction error in (8) and at the same time to satisfy the following H∞ robustness requirement,(10)I∞=∥Φ(z-1)∥∞=:supw∈[0,π]|Φ(e-jw)|<ε, where ε is a positive value. In other words, the optimal reconstruction in (8) and H∞ robustness in (10) should be satisfied simultaneously.Searching the coefficients ofD(z-1) via a genetic algorithm to solve the mixed H2/H∞ deconvolution filter design problem is the key point in our design. This nonlinear complicated design problem will be treated by GA/SA in the next section.


## 3. Deconvolution Filter Design

### 3.1. GA for Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

The genetic algorithm is composed of three operations: (1) reproduction, (2) crossover, and (3) mutation [10–14]. These operations are implemented by performing the basic tasks of copying strings, exchanging portion of strings, and changing the state of bit from 1's to 0's or from 0's to 1's. These operations ensure that the optimum members of the population survive and their information contents are preserved and combined to generate better offspring, which then improves the next generation's performance. We describe the genetic algorithm in the next subsection [14].


#### 3.1.1. Fitness and Cost Function

In this study, the cost function (or energy) is defined as follows:(11)K=E{e2(k)}.Our objective is to search(a1⋯an,b0⋯bn,p) to satisfy the constraint I∞<ε and then to achieve the minimization of (8). In a genetics-based design procedure, a chromosome generates a cost function K and returns a value. The value of the cost is then mapped into a fitness value Fit so as to fit into the genetic algorithm. The fitness value is a reward based on the performance of the possible solution represented by the string, or it can be thought of as how well a fixed-order deconvolution filter can be tuned according to the string to minimize the cost function K. In this paper, we use windowing mapping [14]. Windowing mapping does not allow the fitness difference between good and bad chromosomes to become smaller then the decay of chromosomes that cause the performance of the algorithm to decrease. In this paper, the relationship between K and Fit is expressed as the following [5]:(12)Fit=Fitworst+Fitbest-FitworstKbest-Kworst(K-Kworst), where Kbest and Kworst are the smallest and the largest values evaluated in the generation, and Fitbest and Fitworst are the corresponding fitness values [5, 14].The following three operations are employed in the genetic algorithm to search for the global optimal solution (i.e., the best fitness) in (12) without becoming trapped at local minima [10–14].


#### 3.1.2. Reproduction

Reproduction is based on the principle of survival of the fitness. The fitnessFiti of the ith string is assigned to each individual string in the population where a higher Fiti implies a better fitness. Strings with a large fitness have a greater number of copies in the new generation, for instance, in the roulette wheel selection, the ith string with high fitness value Fi is given a proportionately high probability of reproduction, Ri. The distribution of Ri can be presented as follows: (13)Ri=Fiti∑Fiti‍.Once the strings are reproduced or copied for possible use in the next generation in a mating pool, they wait for the action of the other operators: crossover and mutation.


#### 3.1.3. Crossover

If the chromosomes are only reproduced, they search for the optimum existing individual and do not create any new individuals. Crossover provides a mechanism for strings to mix and match the desirable qualities through a random process. After reproduction, simple crossover proceeds in three steps. First, two newly-reproduced strings are selected from the mating pool, as produced by reproduction. Second, a position with the two strings is uniformly selected at random. The third step involves exchanging all characters by following the crossing site [14].


#### 3.1.4. Mutation

Reproduction and crossover are the primary functions that govern the search power of genetic algorithms. The third operation, mutation, enhances the ability of genetic algorithms to search for the optimal solution. Mutation is the occasional alternation of a value at a particular string position and an insurance policy against the permanent loss of any simple bit.


### 3.2. SA for Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

#### 3.2.1. Introduction to Simulated Annealing

Metropolis [15] introduced a placeMonte Carlo approach to simulate the evolution to thermal equilibrium of a solid for a given temperature. Kirkpatrick et al. [16] found the analogy between minimizing the cost function of a combinatorial optimization problem and the slow cooling of a solid until it reaches its low-energy ground state. They found that the optimization process can be realized by applying the metropolis algorithm, by substituting cost for energy, configuration for state and viewing temperature as a control parameter. We can also view it as an enhanced version of the technique of local optimization, or iterative improvement, in which an initial solution is repeatedly improved by making small local alterations until no alteration can yield a better solution. SA is an optimization technique that simulates the physical process of heating up a solid and then cooling it down slowly until it crystallizes [17, 18]. It was implemented using a stepwise-exponential decrease of temperature. As the temperature decreases, the atomic energies decrease. A crystal with regular structure is obtained at the state where the system has minimum energy. When the cooling is carried out very quickly, which is known as rapid quench, extensive irregularity and defects are seen in the crystal structure. The system does not reach the minimum energy state and ends in a crystalline state, which has a higher energy. Given a temperature, the probability distribution of system energies can be described by the Boltzmann probability [17],(14)ps(K)∝exp(-KBT), where B is the Boltzmann's constant, T is the temperature and ps(K) is the probability that the system is in a state with energy K.Note thatexp(-K/BT) is a number in the interval (0,1) when B and T are both positive, and so it can be interpreted as a probability that depends on B and T. When T is very high, ps(K) is closed to 1 for all energy states according to equation. It can also be seen that there exists a small probability that the system might have high energy even at low temperatures. Therefore, the statistical distribution of energies allows the system to escape from a local energy minimum.SA can be viewed as an algorithm that generates a sequence of Markov chains for a sequence of decreasing temperature values [17, 18]. At each temperature, the generation process is repeated until the probability distribution of the system states approaches the Boltzmann distribution. If the temperature is decreased slowly enough, the Boltzmann distribution tends to converge to global minimal states. The analysis of SA based on infinite length Markov chains has been carried out in the literature [18]. However, in any implementation of the algorithm the Markov chain is of finite length. Therefore, asymptotic convergence can only be approximated. Due to these approximations, the SA is no longer guaranteed to find a global minimum with probability 1.


#### 3.2.2. SA-Base Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

Simulated annealing consists of three components: (1) generation of neighbor candidate solution by perturbing the current solution, (2) accept the solution based on Boltzmann probability, and (3) iterative procedure.The neighbor candidate solution of the first component is [16, 17]  (15)(a1,…,an,b0,…,bn,p)′=(a1,…,an,b0,…,bn,p)+Norm(0,δa2), where Norm(0,δa2) is a random Gaussian number with zero mean and variance of δa2. The δa2 is formulated as [17] (16)δa2={Cdδa2,ifr(k)<0.2,Ciδa2,ifr(k)>0.2,δa2,ifr(k)=0.2, where r(k) is the success ratio of the changes made during the last k iterations, Cd,Ci are near 1, and Cd<1,Ci>1.In the second part, a new solution is accepted as the current solution when its cost is cheaper than that of the current solution. If a new solution is of higher cost (K is lower), it is accepted with a probability of acceptance ps(K), given by(17)ps(K)=11+exp(-K/BT).In the third part, the last accepted candidate solution becomes the initial solution for the next iteration. The temperature of the next iteration is reduced according to a cooling schedule. The temperature updating (decreasing) rule employed is the following [17]:(18)T=Ct×T, where Ct approaches 1 and Ct<1.From the above analysis, the SA-based deconvolution filter design is divided into the following steps.(1)
Initialize the annealing schedule parameters(T,Cd,Ct,Ci,k,σa2), randomly choose an initial state (a1,…,an,b0,…,bn,p) to satisfy the constraint I∞<ε and compute energy K in (12).(2)
Generate the new-state by(19)(a1,…,an,b0,…,bn,p)′=(a1,…,an,b0,…,bn,p)+Norm(0,δa2).DenotethenewenergyK′ as(20)K′=K((a1,…,an,b0,…,bn,p)′).(3)
Computeps(K′)=1/(1+exp(-K′/BT)); if ps(K')> random [0,1] or new-state energy (K′)< original-state energy (K) then the current-state = new-state.(4)
If the temperatureT is the same in k iterations, then decrease the temperature as T=Ct×T, else T=T.(5)
Repeat step 2 to step 4 until the system is frozen. The flowchart of the SA procedure is presented in Figure2.Figure 2
Flow chart of design procedure.


### 3.3. GA/SA for Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

In order to associate the initial fast convergence and weak dependence on initial parameters of GA with the high accuracy of SA, we combined both methods. The combined GA/SA always starts the search procedure as a pure GA and ends as a pure SA. The transition from GA to SA occurs when the following condition is satisfied. The fittest individual remains the same forLg generations [20]. This condition is satisfied whenever the algorithm converges to an intermediate solution. The solution thus far constitutes a good initial guess to SA. The SA's initial and final temperature, as well the step length can adjust as a function of the fittest individual's energy.Based on the above analysis, the design procedure of GA/SA parameter estimation with noises and missing observations is divided into the following steps.(1)
Given the received dataym(k), the order n of the deconvolution filter D(z-1), and robustness constraint ε and generate a random population (a1,…,an,b0,…,bn,p) of Γ chromosomes.(2)
Check theH∞ robustness constraint I∞(a1,…,al,b0,…,bn,p)<ε. If the robustness constraint is not satisfied, then renew the chromosomes.(3)
Compute theH2 performance I2(a1⋯al,b0⋯bn,p) from (8).(4)
Compute the corresponding fitness valueFit from (12).(5)
Use the GA operators (reproduction, crossover, and mutation) to produce chromosomes of next generation.(6)
Repeat the procedure from step 2 to step 5 until the fittest individual remains the same forLg generations.(7)
Given the best chromosome,(T,Cd,Ct,Ci,k,σa2) are the initial SA values.(8)
Generate the new state(a1,…,an,b0,…,bn,p) from (19).(9)
Check theH∞ robustness constraint I∞(a1⋯al,b0⋯bn,p)<ε. If the robustness constraint is not satisfied, then renew the new state.(10)
Compute the new-state energyK' from (20).(11)
Letps(K')=1/(1+exp(-K'/BT)); if ps(K')> random [0,1) or new-state energy (K′)< original-state energy (K) then the current-state = new-state.(12)
If the temperatureT is the same in k iterations, then decrease the temperature as T=CtT, and compute the δa2 from (16). Then, repeat the procedure from step 8 to step 12 until system is frozen.


## 3.1. GA for Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

The genetic algorithm is composed of three operations: (1) reproduction, (2) crossover, and (3) mutation [10–14]. These operations are implemented by performing the basic tasks of copying strings, exchanging portion of strings, and changing the state of bit from 1's to 0's or from 0's to 1's. These operations ensure that the optimum members of the population survive and their information contents are preserved and combined to generate better offspring, which then improves the next generation's performance. We describe the genetic algorithm in the next subsection [14].


### 3.1.1. Fitness and Cost Function

In this study, the cost function (or energy) is defined as follows:(11)K=E{e2(k)}.Our objective is to search(a1⋯an,b0⋯bn,p) to satisfy the constraint I∞<ε and then to achieve the minimization of (8). In a genetics-based design procedure, a chromosome generates a cost function K and returns a value. The value of the cost is then mapped into a fitness value Fit so as to fit into the genetic algorithm. The fitness value is a reward based on the performance of the possible solution represented by the string, or it can be thought of as how well a fixed-order deconvolution filter can be tuned according to the string to minimize the cost function K. In this paper, we use windowing mapping [14]. Windowing mapping does not allow the fitness difference between good and bad chromosomes to become smaller then the decay of chromosomes that cause the performance of the algorithm to decrease. In this paper, the relationship between K and Fit is expressed as the following [5]:(12)Fit=Fitworst+Fitbest-FitworstKbest-Kworst(K-Kworst), where Kbest and Kworst are the smallest and the largest values evaluated in the generation, and Fitbest and Fitworst are the corresponding fitness values [5, 14].The following three operations are employed in the genetic algorithm to search for the global optimal solution (i.e., the best fitness) in (12) without becoming trapped at local minima [10–14].


### 3.1.2. Reproduction

Reproduction is based on the principle of survival of the fitness. The fitnessFiti of the ith string is assigned to each individual string in the population where a higher Fiti implies a better fitness. Strings with a large fitness have a greater number of copies in the new generation, for instance, in the roulette wheel selection, the ith string with high fitness value Fi is given a proportionately high probability of reproduction, Ri. The distribution of Ri can be presented as follows: (13)Ri=Fiti∑Fiti‍.Once the strings are reproduced or copied for possible use in the next generation in a mating pool, they wait for the action of the other operators: crossover and mutation.


### 3.1.3. Crossover

If the chromosomes are only reproduced, they search for the optimum existing individual and do not create any new individuals. Crossover provides a mechanism for strings to mix and match the desirable qualities through a random process. After reproduction, simple crossover proceeds in three steps. First, two newly-reproduced strings are selected from the mating pool, as produced by reproduction. Second, a position with the two strings is uniformly selected at random. The third step involves exchanging all characters by following the crossing site [14].


### 3.1.4. Mutation

Reproduction and crossover are the primary functions that govern the search power of genetic algorithms. The third operation, mutation, enhances the ability of genetic algorithms to search for the optimal solution. Mutation is the occasional alternation of a value at a particular string position and an insurance policy against the permanent loss of any simple bit.


## 3.1.1. Fitness and Cost Function

In this study, the cost function (or energy) is defined as follows:(11)K=E{e2(k)}.Our objective is to search(a1⋯an,b0⋯bn,p) to satisfy the constraint I∞<ε and then to achieve the minimization of (8). In a genetics-based design procedure, a chromosome generates a cost function K and returns a value. The value of the cost is then mapped into a fitness value Fit so as to fit into the genetic algorithm. The fitness value is a reward based on the performance of the possible solution represented by the string, or it can be thought of as how well a fixed-order deconvolution filter can be tuned according to the string to minimize the cost function K. In this paper, we use windowing mapping [14]. Windowing mapping does not allow the fitness difference between good and bad chromosomes to become smaller then the decay of chromosomes that cause the performance of the algorithm to decrease. In this paper, the relationship between K and Fit is expressed as the following [5]:(12)Fit=Fitworst+Fitbest-FitworstKbest-Kworst(K-Kworst), where Kbest and Kworst are the smallest and the largest values evaluated in the generation, and Fitbest and Fitworst are the corresponding fitness values [5, 14].The following three operations are employed in the genetic algorithm to search for the global optimal solution (i.e., the best fitness) in (12) without becoming trapped at local minima [10–14].


## 3.1.2. Reproduction

Reproduction is based on the principle of survival of the fitness. The fitnessFiti of the ith string is assigned to each individual string in the population where a higher Fiti implies a better fitness. Strings with a large fitness have a greater number of copies in the new generation, for instance, in the roulette wheel selection, the ith string with high fitness value Fi is given a proportionately high probability of reproduction, Ri. The distribution of Ri can be presented as follows: (13)Ri=Fiti∑Fiti‍.Once the strings are reproduced or copied for possible use in the next generation in a mating pool, they wait for the action of the other operators: crossover and mutation.


## 3.1.3. Crossover

If the chromosomes are only reproduced, they search for the optimum existing individual and do not create any new individuals. Crossover provides a mechanism for strings to mix and match the desirable qualities through a random process. After reproduction, simple crossover proceeds in three steps. First, two newly-reproduced strings are selected from the mating pool, as produced by reproduction. Second, a position with the two strings is uniformly selected at random. The third step involves exchanging all characters by following the crossing site [14].


## 3.1.4. Mutation

Reproduction and crossover are the primary functions that govern the search power of genetic algorithms. The third operation, mutation, enhances the ability of genetic algorithms to search for the optimal solution. Mutation is the occasional alternation of a value at a particular string position and an insurance policy against the permanent loss of any simple bit.


## 3.2. SA for Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

### 3.2.1. Introduction to Simulated Annealing

Metropolis [15] introduced a placeMonte Carlo approach to simulate the evolution to thermal equilibrium of a solid for a given temperature. Kirkpatrick et al. [16] found the analogy between minimizing the cost function of a combinatorial optimization problem and the slow cooling of a solid until it reaches its low-energy ground state. They found that the optimization process can be realized by applying the metropolis algorithm, by substituting cost for energy, configuration for state and viewing temperature as a control parameter. We can also view it as an enhanced version of the technique of local optimization, or iterative improvement, in which an initial solution is repeatedly improved by making small local alterations until no alteration can yield a better solution. SA is an optimization technique that simulates the physical process of heating up a solid and then cooling it down slowly until it crystallizes [17, 18]. It was implemented using a stepwise-exponential decrease of temperature. As the temperature decreases, the atomic energies decrease. A crystal with regular structure is obtained at the state where the system has minimum energy. When the cooling is carried out very quickly, which is known as rapid quench, extensive irregularity and defects are seen in the crystal structure. The system does not reach the minimum energy state and ends in a crystalline state, which has a higher energy. Given a temperature, the probability distribution of system energies can be described by the Boltzmann probability [17],(14)ps(K)∝exp(-KBT), where B is the Boltzmann's constant, T is the temperature and ps(K) is the probability that the system is in a state with energy K.Note thatexp(-K/BT) is a number in the interval (0,1) when B and T are both positive, and so it can be interpreted as a probability that depends on B and T. When T is very high, ps(K) is closed to 1 for all energy states according to equation. It can also be seen that there exists a small probability that the system might have high energy even at low temperatures. Therefore, the statistical distribution of energies allows the system to escape from a local energy minimum.SA can be viewed as an algorithm that generates a sequence of Markov chains for a sequence of decreasing temperature values [17, 18]. At each temperature, the generation process is repeated until the probability distribution of the system states approaches the Boltzmann distribution. If the temperature is decreased slowly enough, the Boltzmann distribution tends to converge to global minimal states. The analysis of SA based on infinite length Markov chains has been carried out in the literature [18]. However, in any implementation of the algorithm the Markov chain is of finite length. Therefore, asymptotic convergence can only be approximated. Due to these approximations, the SA is no longer guaranteed to find a global minimum with probability 1.


### 3.2.2. SA-Base Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

Simulated annealing consists of three components: (1) generation of neighbor candidate solution by perturbing the current solution, (2) accept the solution based on Boltzmann probability, and (3) iterative procedure.The neighbor candidate solution of the first component is [16, 17]  (15)(a1,…,an,b0,…,bn,p)′=(a1,…,an,b0,…,bn,p)+Norm(0,δa2), where Norm(0,δa2) is a random Gaussian number with zero mean and variance of δa2. The δa2 is formulated as [17] (16)δa2={Cdδa2,ifr(k)<0.2,Ciδa2,ifr(k)>0.2,δa2,ifr(k)=0.2, where r(k) is the success ratio of the changes made during the last k iterations, Cd,Ci are near 1, and Cd<1,Ci>1.In the second part, a new solution is accepted as the current solution when its cost is cheaper than that of the current solution. If a new solution is of higher cost (K is lower), it is accepted with a probability of acceptance ps(K), given by(17)ps(K)=11+exp(-K/BT).In the third part, the last accepted candidate solution becomes the initial solution for the next iteration. The temperature of the next iteration is reduced according to a cooling schedule. The temperature updating (decreasing) rule employed is the following [17]:(18)T=Ct×T, where Ct approaches 1 and Ct<1.From the above analysis, the SA-based deconvolution filter design is divided into the following steps.(1)
Initialize the annealing schedule parameters(T,Cd,Ct,Ci,k,σa2), randomly choose an initial state (a1,…,an,b0,…,bn,p) to satisfy the constraint I∞<ε and compute energy K in (12).(2)
Generate the new-state by(19)(a1,…,an,b0,…,bn,p)′=(a1,…,an,b0,…,bn,p)+Norm(0,δa2).DenotethenewenergyK′ as(20)K′=K((a1,…,an,b0,…,bn,p)′).(3)
Computeps(K′)=1/(1+exp(-K′/BT)); if ps(K')> random [0,1] or new-state energy (K′)< original-state energy (K) then the current-state = new-state.(4)
If the temperatureT is the same in k iterations, then decrease the temperature as T=Ct×T, else T=T.(5)
Repeat step 2 to step 4 until the system is frozen. The flowchart of the SA procedure is presented in Figure2.Figure 2
Flow chart of design procedure.


## 3.2.1. Introduction to Simulated Annealing

Metropolis [15] introduced a placeMonte Carlo approach to simulate the evolution to thermal equilibrium of a solid for a given temperature. Kirkpatrick et al. [16] found the analogy between minimizing the cost function of a combinatorial optimization problem and the slow cooling of a solid until it reaches its low-energy ground state. They found that the optimization process can be realized by applying the metropolis algorithm, by substituting cost for energy, configuration for state and viewing temperature as a control parameter. We can also view it as an enhanced version of the technique of local optimization, or iterative improvement, in which an initial solution is repeatedly improved by making small local alterations until no alteration can yield a better solution. SA is an optimization technique that simulates the physical process of heating up a solid and then cooling it down slowly until it crystallizes [17, 18]. It was implemented using a stepwise-exponential decrease of temperature. As the temperature decreases, the atomic energies decrease. A crystal with regular structure is obtained at the state where the system has minimum energy. When the cooling is carried out very quickly, which is known as rapid quench, extensive irregularity and defects are seen in the crystal structure. The system does not reach the minimum energy state and ends in a crystalline state, which has a higher energy. Given a temperature, the probability distribution of system energies can be described by the Boltzmann probability [17],(14)ps(K)∝exp(-KBT), where B is the Boltzmann's constant, T is the temperature and ps(K) is the probability that the system is in a state with energy K.Note thatexp(-K/BT) is a number in the interval (0,1) when B and T are both positive, and so it can be interpreted as a probability that depends on B and T. When T is very high, ps(K) is closed to 1 for all energy states according to equation. It can also be seen that there exists a small probability that the system might have high energy even at low temperatures. Therefore, the statistical distribution of energies allows the system to escape from a local energy minimum.SA can be viewed as an algorithm that generates a sequence of Markov chains for a sequence of decreasing temperature values [17, 18]. At each temperature, the generation process is repeated until the probability distribution of the system states approaches the Boltzmann distribution. If the temperature is decreased slowly enough, the Boltzmann distribution tends to converge to global minimal states. The analysis of SA based on infinite length Markov chains has been carried out in the literature [18]. However, in any implementation of the algorithm the Markov chain is of finite length. Therefore, asymptotic convergence can only be approximated. Due to these approximations, the SA is no longer guaranteed to find a global minimum with probability 1.


## 3.2.2. SA-Base Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

Simulated annealing consists of three components: (1) generation of neighbor candidate solution by perturbing the current solution, (2) accept the solution based on Boltzmann probability, and (3) iterative procedure.The neighbor candidate solution of the first component is [16, 17]  (15)(a1,…,an,b0,…,bn,p)′=(a1,…,an,b0,…,bn,p)+Norm(0,δa2), where Norm(0,δa2) is a random Gaussian number with zero mean and variance of δa2. The δa2 is formulated as [17] (16)δa2={Cdδa2,ifr(k)<0.2,Ciδa2,ifr(k)>0.2,δa2,ifr(k)=0.2, where r(k) is the success ratio of the changes made during the last k iterations, Cd,Ci are near 1, and Cd<1,Ci>1.In the second part, a new solution is accepted as the current solution when its cost is cheaper than that of the current solution. If a new solution is of higher cost (K is lower), it is accepted with a probability of acceptance ps(K), given by(17)ps(K)=11+exp(-K/BT).In the third part, the last accepted candidate solution becomes the initial solution for the next iteration. The temperature of the next iteration is reduced according to a cooling schedule. The temperature updating (decreasing) rule employed is the following [17]:(18)T=Ct×T, where Ct approaches 1 and Ct<1.From the above analysis, the SA-based deconvolution filter design is divided into the following steps.(1)
Initialize the annealing schedule parameters(T,Cd,Ct,Ci,k,σa2), randomly choose an initial state (a1,…,an,b0,…,bn,p) to satisfy the constraint I∞<ε and compute energy K in (12).(2)
Generate the new-state by(19)(a1,…,an,b0,…,bn,p)′=(a1,…,an,b0,…,bn,p)+Norm(0,δa2).DenotethenewenergyK′ as(20)K′=K((a1,…,an,b0,…,bn,p)′).(3)
Computeps(K′)=1/(1+exp(-K′/BT)); if ps(K')> random [0,1] or new-state energy (K′)< original-state energy (K) then the current-state = new-state.(4)
If the temperatureT is the same in k iterations, then decrease the temperature as T=Ct×T, else T=T.(5)
Repeat step 2 to step 4 until the system is frozen. The flowchart of the SA procedure is presented in Figure2.Figure 2
Flow chart of design procedure.


## 3.3. GA/SA for Fixed-Order MixedH2/H∞ Optimal Deconvolution Filter Design

In order to associate the initial fast convergence and weak dependence on initial parameters of GA with the high accuracy of SA, we combined both methods. The combined GA/SA always starts the search procedure as a pure GA and ends as a pure SA. The transition from GA to SA occurs when the following condition is satisfied. The fittest individual remains the same forLg generations [20]. This condition is satisfied whenever the algorithm converges to an intermediate solution. The solution thus far constitutes a good initial guess to SA. The SA's initial and final temperature, as well the step length can adjust as a function of the fittest individual's energy.Based on the above analysis, the design procedure of GA/SA parameter estimation with noises and missing observations is divided into the following steps.(1)
Given the received dataym(k), the order n of the deconvolution filter D(z-1), and robustness constraint ε and generate a random population (a1,…,an,b0,…,bn,p) of Γ chromosomes.(2)
Check theH∞ robustness constraint I∞(a1,…,al,b0,…,bn,p)<ε. If the robustness constraint is not satisfied, then renew the chromosomes.(3)
Compute theH2 performance I2(a1⋯al,b0⋯bn,p) from (8).(4)
Compute the corresponding fitness valueFit from (12).(5)
Use the GA operators (reproduction, crossover, and mutation) to produce chromosomes of next generation.(6)
Repeat the procedure from step 2 to step 5 until the fittest individual remains the same forLg generations.(7)
Given the best chromosome,(T,Cd,Ct,Ci,k,σa2) are the initial SA values.(8)
Generate the new state(a1,…,an,b0,…,bn,p) from (19).(9)
Check theH∞ robustness constraint I∞(a1⋯al,b0⋯bn,p)<ε. If the robustness constraint is not satisfied, then renew the new state.(10)
Compute the new-state energyK' from (20).(11)
Letps(K')=1/(1+exp(-K'/BT)); if ps(K')> random [0,1) or new-state energy (K′)< original-state energy (K) then the current-state = new-state.(12)
If the temperatureT is the same in k iterations, then decrease the temperature as T=CtT, and compute the δa2 from (16). Then, repeat the procedure from step 8 to step 12 until system is frozen.


## 4. Design Example

A numerical example is given to demonstrate the feasibility of applying GA/SA to the design of fixed-order mixedH2/H∞ optimal deconvolution filter with missing observations. The parameters of GA/SA are set as follows [4, 14, 21]:(21)Γ=200,pc=0.71,pm=0.01,Lg=7,T=2,k=100,Ct=0.99,Cd=0.95,Ci=1.05,δa2=1, where pc is crossover probability and pm is the mutation probability. The number of evaluations is set to 10 000. In the case of the genetic-based estimation, the population size is Γ=200, and the generation number is equal to 50. The simulated results are obtained by averaging 50 independent Monte Carlo(MC) runs.Consider a nonminimum phase deconvolution system in Figure1 with the received signal corrupted by colored noise such that(22)ym(k)=g(k)[H(z-1)u(k)+N(z-1)n(k)]+v(k),u(k)=C(z-1)ω(k), where the missing sequence g(k) is a random process using Bernoulli modulations given in (3), with p=0.9,0.8,0.7. The robustness constraint, signal model, channel model, and noise model are described as follows:(23)ε=0.6,C(z−1)=1+0.8z−11+0.5z−1,H(z−1)=0.9+1.45z−1+0.5z−2(1+0.1z−1)(1-0.1z-1),N(z−1)=1+1.7z−1+0.72z−2(1+0.1z−1)(1-0.1z−1).The driving signalω(k) with disturbance noise n(k) and measurement noise v(k) are also assumed to be independent, stationary, and white with zero mean and variances σω2=1, σn2=0.1, σv2=0.1, respectively. The proposed algorithm is used to design two kinds of fixed-order deconvolution filters. Those results are summarized in Table 1. The second-order IIR(2) filters generated by the proposed algorithm are labeled D1, D2, and D3 with p=0.9, p=0.8, and p=0.7, respectively. The tenth order FIR(10) filters by proposed algorithm are labeled D2, D4, and D6 with p=0.9, p=0.8, and p=0.7, respectively. D7 is a third-order IIR(3) filter with no missing data using the MMSE criterion [4]. The second-order (2) filter generated by the pure SA are labeled D8, D9, D10 with p=0.9, p=0.8, and p=0.7, respectively. The error variances of various criteria and robustness performance are tabulated in Tables 2, 3. Inspecting Tables 2, 3, the mean square error of the proposed method is smaller than that of the full order optimal deconvolution method [4] under different missing observations. Obviously, the reconstruction performance of the conventional optimal deconvolution filter deteriorates because of the missing observations. The reconstruction performance is improved significantly if the missing probability is considered in the design procedure in the case of missing observations. From Table 2, 3, it is evident that the robustness is greater than that of the MMSE deconvolution filter when the channel suffers from small perturbations, so the proposed fixed-order mixed H2/H∞ is more appealing for practical applications. The convergence of the cost function K of the second-order mixed H2/H∞ optimal filter via GA/SA is shown in Figures 3, 5 with different probabilities of missing data. Note that the cost functions of the GA-based estimation method have exponential and rapid convergence at the beginning of generation, but its performance improves very slowly as it approaches the desired global extreme. The proposed GA/SA-based estimation algorithm always starts the search procedure as a pure-GA parameter estimation and ends as a pure-SA parameter estimation. From Figures 3, 5, the shift of algorithm occurs at the 15th generation in the case of p=0.9, at the 23th generation in the case of p=0.8, and at the 21th generation in the case of p=0.7, respectively. From Figures 3, 5, it is seen that the convergence of the GA/SA-based algorithm is rapid initially and the performance is more accurate than that of the GA-based algorithm.Table 1
Filter style.
MethodMissing probabilityFilter styleDeconvolution filterProposed methodp=0.9IIR(2)D1=0.4446−0.0194z−1+0.0348z−21+0.7627z−1+0.268z−2FIR(10)D2=0.3718−0.2672z−1+0.1542z−2−0.0365z−3+0.0388z−4+0.0449z−5−0.0003z−6−0.0409z−7+0.0653z−8−0.0366z−9+0.0824z−10p=0.8IIR(2)D3=0.3483−0.0108z−1+0.0609z−21+0.3667z−1+0.383z−2FIR(10)D4=0.2271−0.1437z−1+0.0027z−2+0.0588z−3+0.0203z−4−0.0693z−5+0.0308z−6+0.0422z−7−0.0297z−8+0.0719z−9−0.0161z−10p=0.7IIR(2)D5=0.3115−0.01874z−1+0.0255z−21+0.2714z−1+0.091z−2FIR(10)D6=0.2686−0.0925z−1−0.098z−2+0.1106z−3−0.0337z−4−0.0081z−5+0.0976z−6+0.0208z−7−0.001z−8+0.0289z−9+0.04z−10MMSE [4]Not missingIIR(3)D7=1.11+0.7737z−1−0.0111z−2−0.77z−31.1+2.42z−1+1.727z−2+0.39z−3SA with initial =1−1z−1+0.25z−21−1z−1+0.25z−2p=0.9IIR(2)D8=0.4233+0.0269z−1+0.0393z−21+0.787z−1+0.0393z−2p=0.8IIR(2)D9=0.2370−0.4954z−1+0.1406z−21−0.7008z−1−0.1758z−2p=0.7IIR(2)D10=0.5189−0.5675z−1−0.0823z−21−0.5675z−1−0.0823z−2Table 2
Performance for various filters and the robustness property of deconvolution filters under various channel perturbations.
Missing probabilityDeconvolution filterOriginal channelH=(0.9+1.08z−1)(1+0.55z−1)(1+0.1z−1)(1−0.1z−1)H=(0.9+1.08z−1)(1+0.55z−1)(1+0.15z−1)(1−0.2z−1)p=0.9D10.62970.70540.713D20.69080.75210.7617D723.928929.41831.6052D80.68280.70720.718p=0.8D30.7830.79880.8289D40.86030.85690.8738D741.255246.345649.3608D91.37981.46521.5395p=0.7D50.85440.86130.873D60.94010.94850.9561D755.240156.769963.4307D100.97081.0541.0936Table 3
Performance for various missing probability and the robustness property of deconvolution filters under various channel perturbations.
Deconvolution filterp=0.9p=0.8p=0.7D10.62970.80070.9099D20.69080.80480.9327D30.71230.7830.8558D40.80050.86030.9178D50.72810.78780.8544D60.82800.87350.9401D723.928941.255255.2401D80.68280.81280.9448D91.37951.37981.4026D100.84010.91420.9708The cost function of the proposed GA/SA estimation algorithm withp=0.9 and an IIR(2) filter.
(a)(b)The cost function of the proposed GA/SA estimation algorithm withp=0.8 and an IIR(2) filter.
(a)(b)The cost function of the proposed GA/SA estimation algorithm withp=0.7 and an IIR(2) filter.
(a)(b)


## 5. Conclusion

In this paper, design methods of fixed-order mixedH2/H∞ optimal deconvolution filters with missing observations have been introduced via GA/SA. This deconvolution filter design method takes advantage of H2 optimal reconstruction performance and H∞ robustness against channel variation and noise uncertainties. Simulation results indicate that the proposed method behaves well even in cases with 10% to 30% missing observations. It is obvious that the reconstruction performance is improved significantly if the missing probability is considered in the proposed deconvolution filter design procedure in the case of missing data. The GA/SA design algorithm rapidly converges and the reconstruction performance is acceptable even if the order of the deconvolution filter is lower. The proposed design methods are suitable for lower order IIR and FIR deconvolution filter designs In addition, the proposed design method is easy to implement and saves operation time. Moreover, it is useful for practical application in the signal reconstruction problems with high-order channels and signal models and missing observations.



---

*Source: 530803-2009-03-22.xml*