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.
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, for all and for some (the period), such that perturbed solutions either approach (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 Figure 1, in which the nonlinearity exists as a separable element in an otherwise linear system, is given as where is the SIDF representing the nonlinear element and 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 where and are the first harmonics in the Fourier series of , and 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 , therefore for brevity is used to denote the SIDF.
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 element which gives output for an input . It is required to examine the behavior of under an input which is exactly periodic. Let this be where is the vector of amplitudes of the harmonic components at the input to the nonlinear element and for a finite number of harmonics , denotes transposition of . Let be an analogous vector of complex Fourier coefficients of the nonlinearity output , then an infinite matrix (the describing function matrix) may be defined as follows:
The 0th column is given by The 1st column is given by , 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 th column is given by, where 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 () can be ignored.
Now define the linear part and representing the nonlinear element with DFM for harmonic components, the condition for existence of oscillation is that the following simultaneous harmonic balance equation (3.3) has a solution . The reason for being a diagonal matrix is that different harmonics do not interact in passing through the linear element.
Consider the nonlinear SISO system shown in Figure 1 in which 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 can be constructed as where is the phase shift between the two harmonics and are respectively, the first and the third coefficients of Fourier components at the nonlinear output, due to the component only, and 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 as given by (3.4), the harmonic balance equation (3.3) is reduced to the following two simultaneous equations where is the second order (first and the third harmonic) DFM and is a diagonal matrix equation (3.7) may be written as where represent the elements of matrix as defined by (3.5).
Rearranging the second equation in the equation set (3.8) A search has to be carried out in the space of in that order) to find those values of 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 where (3.10) represents the overall harmonic balance equation for both the fundamental and the third harmonic components. Hence, a limit cycle with parameters exists when the following two equations are simultaneously satisfied; 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) For the exact solution of the simultaneous equations (3.12), values of 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 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 Figure 1 can be expressed as where and is the matrix of sinusoidal input describing function corresponding to the nonlinear elements of , and 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 of 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 Thus for no limit cycle to exist, no eigenvalue of can equal (−1, ).
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: where and are the th elements of the matrices, for single valued nonlinear elements and , respectively. Inequality (5.3) implies that no limit cycle operation is theoretically possible if Gershgorin bands associated with each diagonal element of 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 and do not intercept for every . 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 has elements of the form , and is a column vector at the input to the nonlinear elements such that , for all over arbitrary ranges , 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 a finite number of sets of sinusoids which will satisfy the condition for harmonic signal balance in the th system loop as given by (5.5) are found; Next, the th loop is considered but using only those finite sets of solution sinusoids derived from (5.5). If there are solution sets of sinusoids () which satisfy both equations, the th loop equation is next examined using only the 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 the autonomous system shown in Figure 2, the set of equations governing this model is given by (5.6) 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 and , where and 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.
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 is said to be nondominated if for a population of N individuals there are no other individuals such that 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 of admissible points such thatAnd such that () 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 Figure 3.
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, and 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 . The range of limit cycle parameters are specified as 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 . 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.
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, two 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)where and are the values of the right-hand side of (5.6) which ideally should tend towards zero for specific values of , and . Upper and lower bounds are specified for , 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 , 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 , 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 Figure 2
with four similar nonlinear elements of the form shown in Figure 5.
The linear system matrix is 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 are specified as 0.1 and 1.5 for all three
parameters and 0.1 to radians for and . 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 () + Average ())
and the solid curve represent the sum of the best values of the two objectives
(best () + best ()). 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 Table 4, 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.
Example 7.2. The nonlinear system for this example
is of the same configuration as in Figure 2, 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 Table 3 and the
following ranges were specified for the parameters
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 Example 7.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 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.
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.