# Limit Cycle Prediction Based on Evolutionary  Multiobjective Formulation

Authors: M. Katebi; H. Tawfik; S. D. Katebi
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

This paper is concerned with an evolutionary search for limit cycle operation in a class of nonlinear systems. In the first part, single input single output (SISO) systems are investigated and sinusoidal input describing function (SIDF) is extended to those cases where the key assumption in its derivation is violated. Describing function matrix (DMF) is employed to take into account the effects of higher harmonic signals and enhance the accuracy of predicting limit cycle operation.
In the second part, SIDF is extended to the class of nonlinear multiinput multioutput (MIMO) systems containing separable nonlinear elements of any general form. In both cases linearized harmonic balance equations are derived and the search for a limit cycle is formulated as a multiobjective problem. Multiobjective genetic algorithm (MOGA) is utilized to search the space of parameters of theoretically possible limit cycle operations. Case studies are presented to demonstrate the effectiveness of the proposed approach.


---

## Body

## 1. Introduction

The theory of linear dynamic systems
is now well understood and is widely applied to many fields of engineering such
as robotics, processes control, ship stirring to name a few. However, nonlinear
systems have received less attention, the reason being the diversity and
complexity of these systems. With the advent of fast and powerful digital
computers, research for a more precise and accurate analysis of nonlinear
systems has grown considerably [1–3]. One such method which traditionally has
been applied is the replacement of nonlinear behavior with a quasi-linear gain
called describing function (DF) [4].Describing function theory and
techniques represent a powerful mathematical approach for analyzing and designing
nonlinear systems. The main motivation for DF techniques is the need to understand
the behavior of nonlinear systems, which in turn is based on the simple fact
that every system is nonlinear, except in very limited operating regimes. The
basic philosophy of DF is to replace each nonlinear element with a quasi-linear
descriptor or describing function. The functional form of such a descriptor is
governed by several factors, the type of input signal, (which is assumed in
advance), and the approximation criterion, for example, minimization of mean squared
error. This technique is dealt with very thoroughly in a number of texts for
the case of nonlinear systems with a single nonlinearity [4]. One category of
DFs that has been particularly successful is the sinusoidal input describing function
(SIDF).The fundamental ideas and use of the
SIDF approach can best be introduced by overviewing the most common
application, limit cycle analysis for a system with a single nonlinearity. A
limit cycle (LC) is a periodic signal,xLC(t+T)=x(t) for all t and for some T  (the period), such
that perturbed solutions either approach xLC  (a stable limit cycle) or diverge from it (an
unstable one). Recently several analytical techniques have been proposed in the
framework of the Bifurcation theory for investigating different aspects of
limit cycle operation  [5, 6]. Nonlinear techniques for modeling the periodic
signal in general and limit cycles in particular have also been suggested [7]. However, an approach to LC analysis that has gained widespread acceptance is
the graphical frequency-domain SIDF method [8, 9].In this paper the describing function
method is first extended to the case of single input single output (SISO)
systems, where the condition for elementary SIDF application is not strictly satisfied
[10]. When the linear system element does not attenuate the super-harmonic components around the feedback loop, the input
to the nonlinear element will not be a pure sinusoid and will be a distorted
waveform containing higher harmonics. The first part of this paper uses the
(DFM), to account for higher harmonics and widen the scope of applications of
the SIDF method [11].The second part of the paper further
extends the SIDF techniques to a class of multiloop nonlinear systems in which
the nonlinear elements are separable from the linear part. In both cases,
emphasis is placed on the multiobjective formulation of predicting the limit
cycle operation and the subsequent solution of the harmonically linearised
system equations by the multiobjective genetic algorithms (MOGA).


## 2. Harmonic Analysis

The
harmonic balance equation for the autonomous feedback configuration of Figure1, in which the nonlinearity exists as a separable element in an otherwise
linear system, is given as (2.1)1+N(a,ω)G(jω)=0, where N(a,ω) is the SIDF representing the nonlinear element
and G(jω) is the frequency transfer function of the
linear part. In a simple harmonic analysis, some form of the solution of (2.1) is sought. However, the valid application of the SIDF requires that the
input signal to the nonlinear element be essentially sinusoidal in form. This
is a condition which imposes an overall low-pass frequency characteristic on
the linear system elements such that the super-harmonic signals are attenuated
around the feedback loop. Based on the assumption that the input to the
nonlinear element is a pure sinusoidal, the SIDF is derived as (2.2)N(a,ω)=(b1+jc1)ejωtaejωt=1a(b1+jc1), where b1 and c1 are the first harmonics in the Fourier series
of  N,
and a is the amplitude of the input sinusoidal. In general the SIDF is a function of both input amplitude
and input frequency, however, for single valued, sector bounded nonlinearities
which are   the subject of this work, the describing function   is real and only a
function of input amplitude a,
therefore for brevity N(a) is used to denote the SIDF.Figure 1
Feed back configuration
for nonlinear systems.


## 3. Extension of Higher
Harmonic Analysis to  SISO Systems

In cases where the requirements of
the overall low-pass characteristics of the linear system elements are met, the
conventional graphical technique in the frequency domain has been shown to give
acceptable results [9]. This elementary SIDF analysis and the associate
graphical solution are well documented [4]. However, in those cases where the
amplitude of  higher harmonics of the
Fourier series at the output of the nonlinear element are significant and are
not suppressed around the feedback loop, the key assumption in derivation of the
SIDF is violated [10]. The reason being that the input to the nonlinear element
will be a distorted waveform containing the fundamental as well as other
harmonics. One method to remedy this circumstance is to find error bounds for
the elementary SIDF. Based on Tsypkin's method, a strategy with which to find
systems with a low-pass linear part for which the describing function technique
erroneously predicts limit cycles is outlined in [10]. A more attractive method that widens the scope
of the SIDF applications, as well as achieving a more accurate limit cycle
prediction is to account for additional harmonics. An appropriate technique is
the so called “describing function matrix” (DFM), which is a mathematical
formulation that allows an arbitrary number of harmonics to be taken into account
in the SIDF analysis [11, 12].Consider a nonlinear elementN which gives output y=N(x) for an input x. It is required to examine the
behavior of N under an input which is exactly periodic. Let
this be(3.1)x=Re[∑r=0∞arejrωt], where a¯=(a0,a1,a2…)T  is the vector of amplitudes of the harmonic components
at the input to the nonlinear element and for a finite number of harmonics (r=m),a¯m=(a0,a1,a2,…,am)T, T
denotes transposition of a¯. Let C be an analogous vector of complex Fourier coefficients of the nonlinearity
output y, then an infinite matrix N (the describing function matrix) may be
defined as follows:The 0th column is given byNk0=C(a0)/a0,The 1st column is given byNk1=[(Ck(a1)−Ck(a0))/a1], where, for example, the (3.1),  (2.1)
element represents the amplitude of the 3rd output harmonic divided by that of the 1st  input harmonic. In general, the jth  column is given by, (3.2)Nkj=[Ck(aj)−Ck(aj−1)aj], where k signifies the row number and
is the index of harmonic components. If the nonlinear element is not frequency
dependent and no bias level is present, then the 0th column can be ignored, and in this case the
first element of the first column is the normal SIDF. Also note that if the
nonlinear element is symmetrical and odd valued, then even rows (k=2,4,6,…)
can be ignored.Now define the linear partGm=Diag(g(0jω),g(jω),g(j2ω),…,g(jmω))  and representing the nonlinear element with
DFM for m harmonic components, the condition for existence of oscillation is that
the following simultaneous harmonic balance equation  (3.3)(3.3)(I+NmGm)a¯m=0 has a solution   ψm=[ωam]. The reason for G  being a diagonal matrix is that different harmonics do not interact in
passing through the linear element.Consider the nonlinear SISO system
shown in Figure1 in which G(jω)andN(a)  are scalars. Assume a fundamental plus a third
harmonic to be present at the input to the nonlinear element. Further, assume
that the nonlinear element is not frequency dependent, but is symmetrical and
odd valued function. The describing function matrix for an input of the following
form(3.4)x(t)=a1sinωt+a3sin(3ωt+ϕ)
can be constructed as (3.5)N≜[c1(a1)a1c1(a1,a3)−c1(a1)a3c3(a1)a1c3(a1,a3)−c3(a3)a3], where ϕ is the phase shift between the two harmonics
and c1(a1),c3(a3) are respectively, the first and the third coefficients
of Fourier components at the nonlinear output, due to the component a1sin(ωt) only, and c1(a1,a3),c3(a1,a3) are the corresponding coefficients due to the
sum of the two signals at the input to the nonlinear element.In
the case of only two harmonics such asx(t) given by  (3.4), the harmonic
balance equation  (3.3)  is reduced to the following two simultaneous equations (3.6)(I+N2G2), where N2 is the second order (first and the third
harmonic) DFM and G2 is a 2×2 diagonal matrix (3.7)(1001)+(n11(a1)n12(a1,a3)n21(a1)n22(a1,a3))(g(jω)00g(j3ω))(a1a3∠ϕ)=0, equation (3.7)  may be written as (3.8)(1+n11(a1)g(jω))a1+n12(a1,a3)g(j3ω)a3∠ϕ=0,n21(a1)g(jω)a1+(1+n22(a1,a3)g(j3ω))a3∠ϕ=0, where  n11,n12,n21,n22 represent the elements of matrix N as defined by (3.5).Rearranging the second
equation in the equation set (3.8)(3.9)a3∠ϕ=−|n21(a1)g(jω)1+n22(a1,a3)g(j3ω)|a1. A search has to be carried out in the
space of  (a1,ϕ,ω,a3   in that order) to find those values of a3 and ϕ which will satisfy (3.9). If such
values are found, the result will be substituted into the first equation of equation
set (3.8), to test for this equation which governs the fundamental oscillation.Substituting (3.9) into the first equation of equation
set (3.8) and rearranging yields (3.10)n11(a1)g(jω)a1+G1=−1,  where G1=−n12(a1,a3)g(j3ω)n21(a1)g(jω)/(1+n22(a1,a3)g(j3ω)),(3.10) represents the overall harmonic balance
equation for both the fundamental and the third harmonic components. Hence, a
limit cycle with parameters (a1,a3,ω,ϕ)  exists when the following two equations are
simultaneously satisfied; (3.11)|n11(a1)g(jω)a1+G1|=1,∠(n11(a1)g(jω)a1+G1)=180±2kπ,k=0,1,…. Alternatively, (3.8) may be rearranged
as follows and be directly used as the objective (fitness function) in the subsequent
MOGA search as the satisfaction of these objectives implies the solution of the
two simultaneous harmonic balance equations (3.8)(3.12)Obj1=|(1+n11(a1)g(jω))a1+n12(a1,a3)g(j3ω)a3exp(jϕ)|≤ε,Obj2=|n21(a1)g(jω)a1+(1+n22(a1,a3)g(j3ω))a3exp(jϕ)|≤ε.  For the exact solution of the simultaneous equations
(3.12), values of Obj1andObj2 should reach zero but ε is a small positive number (as near to zero as
possible). An intelligent
search based on MOGA with the above defined objective (fitness) function is
carried out over the space of (a1,a3,ω,ϕ) to facilitate an efficient and an accurate
prediction of limit cycle in the presence of higher harmonics.


## 4. Extension of Harmonic Analysis to
MIMO System

Extension of describing function
techniques to multiloop nonlinear systems is not new and follows the
development of the frequency domain multivariable linear theory [13, 14]. A
numerical technique with a graphical interpretation for quantifying the limit
cycle parameters in multivariable systems in the frequency domain is reported
in [15].  The extension of the graphical
techniques to multiloop systems with coupled multivalued nonlinear elements
is reported in [16]. Another graphical method based on the phasor diagram which
particularly gives accurate limit cycle prediction for relay systems is
developed in [8]. A computer aided design (CAD)  tool for limit cycle prediction
aimed at educational purposes is reported in [9]. Another novel numerical
technique based on deriving the least damped eigenvalue to the imaginary axis
for nonlinear systems with multiple nonlinearities is suggested in [17]. Based
on defining an appropriate error function, the authors use both eigenvalue and
eigenvectors to formulate a generalized Newton-Raphson method to solve for the
state variable amplitude in a minimum norm sense [17]. Most of the above
mentioned techniques are essentially dependent on the graphical displays in the
frequency domain. While this is useful for getting insight into the subsequent
compensator design, an efficient technique for accurately quantifying the limit
cycle parameters is still highly desirable.In this paper, a multiobjective
formulation is presented to search numerically for limit cycles in a class of multiloop
nonlinear systems. The approach is computationally efficient and is based on multiobjective genetic
algorithms (MOGA).


## 5. Limit Cycle
Prediction in Nonlinear Multivariable Systems

Provided that certain limitations are
placed on the form of the linear system elements, the extension of the of
harmonic linearization to multivariable systems is conceptually
straightforward. The equation governing limit cycle operation in the autonomous
multivariable nonlinear feedback system of Figure1 can be expressed as (5.1)(T(a,jω)+I)a=0, where T(a,jω)=N(a,jω)G(jω) and N(a,ω) is the matrix of sinusoidal input describing
function corresponding to the nonlinear elements of N,
and a is the column vector of the sinusoid at the
inputs to these elements. The use of a single sinusoidal describing function
analysis implies the following assumptions:(a)
Each element ofG(jω)  acts as a low pass system so that higher
harmonic signal components are effectively suppressed.(b)
If a limit-cycle is present, then all loops will oscillate at the same
frequency. Experience indicates that this is a reasonable assumption,
particularly if the nonlinear elements are similar and if the dominant linear
elements have approximately the same frequency characteristics.Equation (5.1) will have a nontrivial solution
only if (5.2)det[N(a,ω)G(jω)+I]=0.  Thus for
no limit cycle to exist, no eigenvalue of  N(a,ω)G(jω) can equal (−1, j0).The conventional graphical frequency domain method for the
solution of (5.2) for single input single output systems with a single
nonlinear element may be extended to MIMO systems employing Nyquist or the
Inverse Nyquist Array.  Invoking the Inverse
Nyquist Array method [13], the possibility of limit cycle existence may be
examined by studying the following inequalities: (5.3)|nkk(ak)+g^kk(jω)|>∑j≠k|njk(ak)+g^jk(jω)|,(5.4)|nkk(ak)+g^kk(jω)|>∑j≠k|njk(ak)|+∑j≠k|g^jk(jω)|, where  njk(a) and g^jk(jω) are the jkth   elements of the matrices, N(a) for single valued nonlinear elements and G−1(jω),
respectively. Inequality (5.3)  implies that no limit cycle operation is
theoretically possible if Gershgorin bands associated with each diagonal
element of (N(a)+G−1(jω)) do not encompass the origin in the complex
frequency domain. Inequality (5.4)  implies the same result, provided that the
bands traced out by the Gershgorin discs on the loci of g^ii(jω) and nii(a) do not intercept for every i. The former representation is computationally more demanding, but gives less
conservative results because of the weakening step between inequalities (5.3)
and (5.4). The Gershgorin discs set bounds to the eigenvalue locations and any
method of limit cycle prediction based on band intersection should have a
tendency to be conservative.An alternative approach lies in calculating eigenvalue of
the harmonically linearised return ratio system equation. For a general system
which contains both on and off diagonal nonlinear elements, the computational
effort involved in determining the eigenvalue becomes almost formidable as, at
any frequency, the return ratio matrix is a function of the signal amplitude at
the input to the nonlinear elements [16].A numerically based technique called the sequential loop
balance method has also been devised [15]. This method is based on (5.1)
for which T(a,jω) has elements of the form tij(aj,jω),
and a is a column vector at the input to the
nonlinear elements such that aj=Ajexp(jϕi),j=1,2,…,n,  for all j,i=1,2,…m over arbitrary ranges A,ω, and ϕ,
a possible infinite number of solutions may exist for (5.1). For a
specified value of frequency and a specified range of discrete values of the
reference signal a1=A1ej0,A1>A2,A3,…,An  a finite number of p sets of sinusoids (a1k,a2k,…,ank),k=1,2,…,p  which will satisfy the condition for harmonic
signal balance in the nth  system loop as given by  (5.5) are found; (5.5)(tnn+1)an+∑j=1n−1tnjaj=0. Next, the (n−1)th loop is considered but using only those finite
sets of solution sinusoids derived from (5.5). If there are q solution
sets of sinusoids (q≤p)  which satisfy both equations, the (n−2)th loop equation is next examined using only the
q sets of solution sinusoids and the process is repeated in a sequential manner
until loop 1 is reached. Clearly this method is also computationally demanding
due to an indirected search over the specified ranges of parameters. Further,
it does not search for the phase difference between the oscillating loops, and as discrete data are used in the computation, there is the possibility of
solution sets which exist over the data intervals.Considering the2×2 autonomous system shown in
Figure 2, the set of equations governing this model is given by (5.6)(5.6)(1+n11(a1)g11(jω))a1+n21(a2)g21(jω)a2ejϕ=0,(n12(a1)g12(jω)a1+1+n22(a2)g22(jω))a2ejϕ=0. In this case there are two equations and four variables
(unknown). The solution of these equations is sought over a range of specific
values of a1,a2,ω and φ,
where a1 and a2 are amplitudes of limit cycles in loop 1 and 2,
respectively, ω is frequency of oscillation for both loops and ϕ is the phase-shift between the loops.Figure 2
A two-input two-output nonlinear system.


## 6. Multiobjective Genetic Algorithms

Genetic algorithms
(GAs) are search procedures based on the evolutionary process in nature. They
differ from other approaches in that they use probabilistic and not
deterministic criteria for progressing the search. The idea is that GA operates
on a population of individuals, each individual representing a potential solution
of the problem, and applies the principle of survival of the fittest to the
population, so that the individuals evolve towards better solutions of the
problem. Each individual is given a chromosoidal representation, which
corresponds to the genotype of an individual in nature. Three operations can be
performed on individuals in the population,selection, cross-over and mutation. These correspond to the selection of individuals
in nature for breeding, where the fitter members of a population breed and
pass-on their genetic material. The cross-over corresponds to the combination
of genes by mating, and mutation to genetic mutation in nature. The selection
is biased so that the “fittest" individuals are more likely to be selected for
cross-over, the fitness being a function of the criteria which is being
minimized. By means of these operations, the population will evolve towards a
solution. Most GAs have been used for single objective optimization problems
[18], although several multiobjective schemes have been proposed [19, 20]. Applications of multiobjective evolutionary schemes to control systems
analysis and design have been widespread [21, 22]. A formulation called the
multiobjective genetic algorithm (MOGA) maintains the genuine multiobjective
nature of the problem, and is essentially the scheme adapted here [23]. Further
details for the MOGA can be found in [20]. The design philosophy of the MOGA
differs from other methods in that a set of simultaneous solutions are sought and
the designer then selects the best solution from the set. The idea behind the
MOGA is to develop a population of Pareto-optimal or near Pareto-optimal
solutions. The aim is to find a set of solutions which are nondominated and satisfy
a set of inequalities. An individual j with a set of objective functions φj=(φ1j,…,φnj)  is said to be nondominated if for a
population of N  individuals there are no other individuals k=1,N;k≠j
such
that (6.1)(a)φik≤φij∀i=1,2,…,n,(b)φik<φijfor at least onei.   The MOGA is set
into a multiobjective context by means of the fitness function. Individuals
are ranked on the basis of the number of other individuals they are dominated
by for the unsatisfied inequalities. Each individual is then assigned fitness
according to their rank. The mechanism is described in detail in [20]. To
summarize, the MOGA problem could be stated as:Find a set    ofM  admissible points  Pj,j=1,…,Msuch that(6.2)φij≤εi(j=1,2,…,m),(i=1,2,…,n).And such thatφj
(j=1,…,M) are nondominated.Genetic algorithms
are naturally parallel and hence lend themselves well to multiobjective settings. They also work well on nonsmooth objective functions. Thus MOGA can be used to
search the existence of any possible limit cycle operation in nonlinear MIMO
systems.


## 7.  Applications

### 7.1. SISO System

The
feedback configuration of the single input single output system considered in
this example is shown in Figure3.Figure 3
SISO nonlinear feedback system.A Fourier
analysis of the nonlinearity output shows that the amplitude of the third
harmonic component is significant. Therefore, the elementary describing
function analysis does not give accurate limit cycle prediction. Based on the
calculation of the required Fourier coefficients at the output of the nonlinear
element for the input of the form given by (3.4), the describing
function matrix as defined by  (3.5) is derived as, (7.1)DFM=[1.91a1−1.12a30.64a10.89a3]a1 and a3 denote the amplitude of the first and third
harmonics at the input to the nonlinear element. The harmonic   balance equation
is formulated in the presence of the third harmonic using (3.12) with ε=0.01. The range of limit cycle parameters are specified as (7.2)0.1≤ω≤1.5,3≤a1≤8.0,0.05≤a3≤0.5,0.1≤ϕ≤2π. The
population is initialized randomly and the MOGA program is run to search the
space of these parameters for the simultaneous satisfaction of objectives
(3.12). The parameters of MOGA are shown in Table 1  and after
18 generations the algorithm converged with the value of ε=0.012. The first three solutions are ranked   based on the concept of non dominance and are
listed in Table 2. The results obtained from simulation are also shown in bold
in Table 2   and are comparable with those predicted by MOGA. The phase space of
the predicted limit cycle is shown in Figure 4.Table 1
MOGA parameters for the SISO system.
Number of generationPopulation sizeSelection methodRepresentation methodMultiobjective methodMutation probabilityCross over probability50100TournamentBinary 8 bit for each parameterPareto0.020.7Table 2
Results of limit cycle prediction with MOGA for
SISO system.
Best rankedω  rad./sec.a1a3φ  radiansObjective 1Objective  210.686.8530.1035.410.010.01520.686.8540.1035.410.0120.00630.686.8580.1035.410.0090.08Simulation (actual)0.696.90.15.21Figure 4
Phase space of limit cycle for SISO example.The computation
time required depends largely on the number of generations, population size and
on the number of bits specified for each parameter. For this example with the
above specifications, 50 generations took 6.2 seconds on a Toshiba
Centrino 1.6 with 512 M Ram and 2 M cashe.


### 7.2. MIMO Systems

In this
section, two2×2 nonlinear systems are presented. In the first
example the elements of the nonlinear matrix are similar and consist of four
ideal relays.In both examples the numerical
solution of    (5.6) is formulated as a 2-objective problem, and the
satisfaction   of these objectives implies the solution of the simultaneous
harmonic balance (5.6)(7.3)Obj1=|(1+n11(a1)g11(jω))a1+n21(a2)g21(jω)a2eϕ|≤ε,Obj2=|n12(a1)g12(jω)a1+(1+n22(a2)g22(jω))a2ejϕ|≤ε,where Obj1 and Obj2 are the values of the right-hand side of
(5.6)
which ideally should tend towards zero for specific values of a1,a2,ω, and ϕ. Upper and   lower bounds are specified for ω,a1,a2, and ϕ,
then the real generational MOGA with specified selection method, cross over,
mutation rate and population size  is
called to search over the  parameter
space of ω,a1,a2, and ϕ. The required numerical accuracy may be achieved by specifying the number of
genes (binary bits) for each individual in accordance with the value of ε. If conditions for inequalities (7.3) exist, then  MOGA converges to the correct values of limit
cycle parameters after a number of generations. If finer and more accurate
parameter values are required, the bounds on the parameters ω,a1,a2,  and ϕ may be tightened and the number of genes
increased on the subsequent run of the MOGA program. In order to avoid the
local minima and premature convergence the mutation probability rate may be
taken higher at the beginning and decreased exponentially toward the end of the
run. One advantage of limit cycle prediction based on MOGA is that multiple
solutions may be distinguished by using the Niching mechanism [19].Due to the inherent approximation in
using the SIDF, and also the nature of multiobjective formulation, it may not
be possible to reach the exact minimum which is zero, as required by the
equation set (7.3). Therefore MOGA may converge to a set of Pareto optimal
solutions and in the following examples the three best ranked are given as the
final solutions.Example 7.1.
The nonlinear control system with two inputs and two outputs has
the configuration of Figure2
with four similar nonlinear elements (n11=n12=n21=n22) of the form shown in Figure 5.
The linear system matrix is(7.4)G=[0.5s(s2+1.5s+0.5)−0.15s2+s+10.8s2+3s+21s(s2+1.6s+0.8)].  In order to make the MOGA search more
realistic, upper and lower bounds are specified for all parameters of possible
limit cycle operation, namely, the frequency and amplitude for each loop. For
this example, lower and upper bounds for ω,a1anda2 are specified as 0.1 and 1.5 for all three
parameters and 0.1 to 2π radians for ϕ and ε=0.1. Parameters of MOGA were also specified as in Table 3. The search is monitored
interactively and if within the specified range of the parameters feasible
solutions exist, MOGA converges to solutions which are ranked according to
their fitness values measured by the concept of non dominance, that is, Pareto
optimal. For a population size of 50, it was observed that initially a large
number of solutions exist. However, as the search progressed most members of the
population tend to be concentrated around a narrow range of parameters. After 25
generations the algorithm converged and the results for the first three solutions
are given in Table 4 and the convergence curve is shown in Figure 6. In
Figure 6 the cost axis is a measure of fitness and the doted curve represent
the sum of the average of the two objectives (average  (Obj1)  + Average (Obj2))
and the solid curve represent the sum of the best values of the two objectives
(best (Obj1) + best (Obj2)). With the parameters given in Table 3 and the
above specified ranges for parameters, the MOGA algorithm converged within 25
generations, consuming 5.2 seconds of CUP time on a Toshiba
Centrino 1.6 with 512 M Ram and 2 M cashe. For the purpose of having an appropriate axes for the
convergence curve the program was allowed to execute all 50 generations which
took 11 seconds on the same machine.
In order to validate the results
given by Table4, the nonlinear system was simulated. The results of
simulation (actual LC) compare well with the results obtained by MOGA. This is
also shown in bold in Table 4 and in Figure 7. The reason for the very
accurate limit cycle prediction in this example is two folds, one is that all
nonlinear elements have the same characteristics and the other is that the
linear system elements have a low pass frequency characteristics.Table 3
MOGA parameters for MIMO Example7.1.
Number of
generationPopulation
sizeSelection
methodRepresentation
methodMultiobjective methodMutation probabilityCross over probability5050TournamentBinary 8 bit for each parameterPareto0.020.7Table 4
Results of limit cycle prediction with MOGA for Example7.1.
Best rankedωa1a2φObj 1Obj 210.863960.6177851.422032.174770.03179520.015062320.8638650.6177871.414852.156230.03368970.0061906430.8458350.6196581.440052.486680.009174390.0861531Simulation
(actual)0.8630.6201.4062.168Figure 5
Ideal relay-nonlinear elements for MIMO Example7.1.Figure 6
Convergence curve for MIMO Example7.1.Figure 7
Simulation results for MIMO Example7.1.Example 7.2.
The nonlinear system for this example
is of the same configuration as in Figure2, and is shown in Figure 8. In
this case the nonlinear matrix consists of different nonlinear behavior such as
saturation, ideal relay, and dead zone.
The same MOGA parameters were chosen as in Table3 and the
following ranges were specified for the parameters (7.5)0.1≤ω≤1.0,0.5≤a1≤1.5,1.0≤a2≤2.0,0.1≤ϕ≤2π.
For this example, the algorithm
converged within 30 generations, consuming 31.2 seconds of CUP time on the
same machine. For all 50 generations 52 seconds of the CPU time is consumed on the same
machine. The convergence curve is shown in Figure 9.
Contrary to Example7.1 the results
of simulation for this example, do not agree well with the results predicted by
MOGA. This is shown (bold) in Table 5 and in Figure 10. As seen in Figure 10,
the oscillation in loop 2 is distorted due to the higher harmonic components. This is caused by the frequency characteristics of element g21 of the linear part, which do not strictly
conform to the filter hypothesis, thus reducing the describing function
accuracy and hence the reason for the small discrepancy between simulation and
MOGA results. The development of similar   analysis as for SISO systems will improve
the SIDF accuracy  significantly.Table 5
Results of limit-cycle prediction with MOGA for
MIMO Example7.2.
ωa1a2φObj 1Obj 210.5177281.185281.732963.67720.02291420.061144620.5176381.180251.733963.677180.01810930.065582830.5301391.162281.682923.676890.04106110.0383831Simulation
(actual)0.54001.0731.5433.500Figure 8
Two dimensional nonlinear systems for Example7.2.Figure 9
Convergence curve for MIMO Example7.2.Figure 10
Simulation results for MIMO Example7.2.


## 7.1. SISO System

The
feedback configuration of the single input single output system considered in
this example is shown in Figure3.Figure 3
SISO nonlinear feedback system.A Fourier
analysis of the nonlinearity output shows that the amplitude of the third
harmonic component is significant. Therefore, the elementary describing
function analysis does not give accurate limit cycle prediction. Based on the
calculation of the required Fourier coefficients at the output of the nonlinear
element for the input of the form given by (3.4), the describing
function matrix as defined by  (3.5) is derived as, (7.1)DFM=[1.91a1−1.12a30.64a10.89a3]a1 and a3 denote the amplitude of the first and third
harmonics at the input to the nonlinear element. The harmonic   balance equation
is formulated in the presence of the third harmonic using (3.12) with ε=0.01. The range of limit cycle parameters are specified as (7.2)0.1≤ω≤1.5,3≤a1≤8.0,0.05≤a3≤0.5,0.1≤ϕ≤2π. The
population is initialized randomly and the MOGA program is run to search the
space of these parameters for the simultaneous satisfaction of objectives
(3.12). The parameters of MOGA are shown in Table 1  and after
18 generations the algorithm converged with the value of ε=0.012. The first three solutions are ranked   based on the concept of non dominance and are
listed in Table 2. The results obtained from simulation are also shown in bold
in Table 2   and are comparable with those predicted by MOGA. The phase space of
the predicted limit cycle is shown in Figure 4.Table 1
MOGA parameters for the SISO system.
Number of generationPopulation sizeSelection methodRepresentation methodMultiobjective methodMutation probabilityCross over probability50100TournamentBinary 8 bit for each parameterPareto0.020.7Table 2
Results of limit cycle prediction with MOGA for
SISO system.
Best rankedω  rad./sec.a1a3φ  radiansObjective 1Objective  210.686.8530.1035.410.010.01520.686.8540.1035.410.0120.00630.686.8580.1035.410.0090.08Simulation (actual)0.696.90.15.21Figure 4
Phase space of limit cycle for SISO example.The computation
time required depends largely on the number of generations, population size and
on the number of bits specified for each parameter. For this example with the
above specifications, 50 generations took 6.2 seconds on a Toshiba
Centrino 1.6 with 512 M Ram and 2 M cashe.


## 7.2. MIMO Systems

In this
section, two2×2 nonlinear systems are presented. In the first
example the elements of the nonlinear matrix are similar and consist of four
ideal relays.In both examples the numerical
solution of    (5.6) is formulated as a 2-objective problem, and the
satisfaction   of these objectives implies the solution of the simultaneous
harmonic balance (5.6)(7.3)Obj1=|(1+n11(a1)g11(jω))a1+n21(a2)g21(jω)a2eϕ|≤ε,Obj2=|n12(a1)g12(jω)a1+(1+n22(a2)g22(jω))a2ejϕ|≤ε,where Obj1 and Obj2 are the values of the right-hand side of
(5.6)
which ideally should tend towards zero for specific values of a1,a2,ω, and ϕ. Upper and   lower bounds are specified for ω,a1,a2, and ϕ,
then the real generational MOGA with specified selection method, cross over,
mutation rate and population size  is
called to search over the  parameter
space of ω,a1,a2, and ϕ. The required numerical accuracy may be achieved by specifying the number of
genes (binary bits) for each individual in accordance with the value of ε. If conditions for inequalities (7.3) exist, then  MOGA converges to the correct values of limit
cycle parameters after a number of generations. If finer and more accurate
parameter values are required, the bounds on the parameters ω,a1,a2,  and ϕ may be tightened and the number of genes
increased on the subsequent run of the MOGA program. In order to avoid the
local minima and premature convergence the mutation probability rate may be
taken higher at the beginning and decreased exponentially toward the end of the
run. One advantage of limit cycle prediction based on MOGA is that multiple
solutions may be distinguished by using the Niching mechanism [19].Due to the inherent approximation in
using the SIDF, and also the nature of multiobjective formulation, it may not
be possible to reach the exact minimum which is zero, as required by the
equation set (7.3). Therefore MOGA may converge to a set of Pareto optimal
solutions and in the following examples the three best ranked are given as the
final solutions.Example 7.1.
The nonlinear control system with two inputs and two outputs has
the configuration of Figure2
with four similar nonlinear elements (n11=n12=n21=n22) of the form shown in Figure 5.
The linear system matrix is(7.4)G=[0.5s(s2+1.5s+0.5)−0.15s2+s+10.8s2+3s+21s(s2+1.6s+0.8)].  In order to make the MOGA search more
realistic, upper and lower bounds are specified for all parameters of possible
limit cycle operation, namely, the frequency and amplitude for each loop. For
this example, lower and upper bounds for ω,a1anda2 are specified as 0.1 and 1.5 for all three
parameters and 0.1 to 2π radians for ϕ and ε=0.1. Parameters of MOGA were also specified as in Table 3. The search is monitored
interactively and if within the specified range of the parameters feasible
solutions exist, MOGA converges to solutions which are ranked according to
their fitness values measured by the concept of non dominance, that is, Pareto
optimal. For a population size of 50, it was observed that initially a large
number of solutions exist. However, as the search progressed most members of the
population tend to be concentrated around a narrow range of parameters. After 25
generations the algorithm converged and the results for the first three solutions
are given in Table 4 and the convergence curve is shown in Figure 6. In
Figure 6 the cost axis is a measure of fitness and the doted curve represent
the sum of the average of the two objectives (average  (Obj1)  + Average (Obj2))
and the solid curve represent the sum of the best values of the two objectives
(best (Obj1) + best (Obj2)). With the parameters given in Table 3 and the
above specified ranges for parameters, the MOGA algorithm converged within 25
generations, consuming 5.2 seconds of CUP time on a Toshiba
Centrino 1.6 with 512 M Ram and 2 M cashe. For the purpose of having an appropriate axes for the
convergence curve the program was allowed to execute all 50 generations which
took 11 seconds on the same machine.
In order to validate the results
given by Table4, the nonlinear system was simulated. The results of
simulation (actual LC) compare well with the results obtained by MOGA. This is
also shown in bold in Table 4 and in Figure 7. The reason for the very
accurate limit cycle prediction in this example is two folds, one is that all
nonlinear elements have the same characteristics and the other is that the
linear system elements have a low pass frequency characteristics.Table 3
MOGA parameters for MIMO Example7.1.
Number of
generationPopulation
sizeSelection
methodRepresentation
methodMultiobjective methodMutation probabilityCross over probability5050TournamentBinary 8 bit for each parameterPareto0.020.7Table 4
Results of limit cycle prediction with MOGA for Example7.1.
Best rankedωa1a2φObj 1Obj 210.863960.6177851.422032.174770.03179520.015062320.8638650.6177871.414852.156230.03368970.0061906430.8458350.6196581.440052.486680.009174390.0861531Simulation
(actual)0.8630.6201.4062.168Figure 5
Ideal relay-nonlinear elements for MIMO Example7.1.Figure 6
Convergence curve for MIMO Example7.1.Figure 7
Simulation results for MIMO Example7.1.Example 7.2.
The nonlinear system for this example
is of the same configuration as in Figure2, and is shown in Figure 8. In
this case the nonlinear matrix consists of different nonlinear behavior such as
saturation, ideal relay, and dead zone.
The same MOGA parameters were chosen as in Table3 and the
following ranges were specified for the parameters (7.5)0.1≤ω≤1.0,0.5≤a1≤1.5,1.0≤a2≤2.0,0.1≤ϕ≤2π.
For this example, the algorithm
converged within 30 generations, consuming 31.2 seconds of CUP time on the
same machine. For all 50 generations 52 seconds of the CPU time is consumed on the same
machine. The convergence curve is shown in Figure 9.
Contrary to Example7.1 the results
of simulation for this example, do not agree well with the results predicted by
MOGA. This is shown (bold) in Table 5 and in Figure 10. As seen in Figure 10,
the oscillation in loop 2 is distorted due to the higher harmonic components. This is caused by the frequency characteristics of element g21 of the linear part, which do not strictly
conform to the filter hypothesis, thus reducing the describing function
accuracy and hence the reason for the small discrepancy between simulation and
MOGA results. The development of similar   analysis as for SISO systems will improve
the SIDF accuracy  significantly.Table 5
Results of limit-cycle prediction with MOGA for
MIMO Example7.2.
ωa1a2φObj 1Obj 210.5177281.185281.732963.67720.02291420.061144620.5176381.180251.733963.677180.01810930.065582830.5301391.162281.682923.676890.04106110.0383831Simulation
(actual)0.54001.0731.5433.500Figure 8
Two dimensional nonlinear systems for Example7.2.Figure 9
Convergence curve for MIMO Example7.2.Figure 10
Simulation results for MIMO Example7.2.


## 8. Conclusion

The single sinusoidal input describing
function is extended by means of describing function matrix to account for
higher harmonic   signal components. In the analysis of single loop nonlinear
systems, the describing function matrix provides a rigorous method, enhances
limit cycle   prediction and widens the scope of the valid application of the SIDF
to those cases for which the assumption of the filter hypothesis is not
strictly met. It is shown how, by formulating limit cycle prediction in the
presence of higher harmonics as a multiobjective problem and using the MOGA
procedure, a solution to the DF matrix equation can be obtained efficiently in
numerical form. This is useful in illustrating both how a multiplicity of
possible solutions may arise, and for emphasizing the effects of higher
harmonic signal content.Next, the single sinusoidal input describing
function technique is extended to multiloop nonlinear systems. For the class
of separable nonlinear elements of any general form, the linearized harmonic
balance equations are derived. The search for limit cycle is formulated as a multiobjective
problem. The multiobjective genetic algorithm is employed to obtain quantitative
values for parameters of possible limit cycle operation. The MOGA search space
is the space of the limit cycle parameters such as amplitudes, frequency and
phase difference between the interacting loops. The algorithm is
computationally efficient and provides very accurate results when the conditions
for the SIDF application are satisfied.In both cases of SISO and MIMO systems,
emphasis is placed on the formulation of limit cycle parameters as a multiobjective problem and subsequent solution by the evolutionary MOGA method. Finally,
examples of use are given to demonstrate the effectiveness of the proposed
approach.



---

*Source: 816707-2009-03-23.xml*