Limit Cycle Predictions of Nonlinear Multivariable Feedback Control Systems with Large Transportation Lags

A practical method is developed for limit-cycle predictions in the nonlinear multivariable feedback control systems with large transportation lags. All nonlinear elements considered are linear independent. It needs only to check maximal or minimal frequency points of root loci of equivalent gains for finding a stable limit-cycle. This reduces the computation effort dramatically. The information for stable limit-cycle checking can be shown in the parameter plane also. Sinusoidal input describing functions with fundamental components are used to find equivalent gains of nonlinearities. The proposed method is illustrated by a simple numerical example and applied to one 2 × 2 and two 3 × 3 complicated nonlinear multivariable feedback control systems. Considered systems have large transportation lags. Digital simulation verifications give calculated results provide accurate limit cycle predictions of considered systems. Comparisons are made also with other methods in the current literature.


Introduction
In current literature, for nonlinear multivariable systems the Nyquist, inverse Nyquist and numerical optimization, methods are usually used to predict the existence of limit cycles.These methods are based upon the graphical or numerical solutions of the linearized harmonic-balance equations [1][2][3][4][5][6][7][8][9][10].It has been shown that for multivariable systems, over arbitrary ranges of amplitudes (A i ), frequency (ω) and phases (θ i ), an infinite number of possible solutions may exist.Gray has proposed a sequential computational procedure to seek the solutions for only specified ranges of discrete values of A i , ω, and θ i , these specified ranges are determined by use of the Nyquist or inverse Nyquist plots [4,5].Although the aforementioned methods are powerful, large computational efforts are usually expected.
In general, real and imaginary parts of the characteristic equation are used as two simultaneous equations to find the solution of the limit cycle for single-input singleoutput (SISO) nonlinear feedback control systems [11][12][13][14][15][16][17].Therefore, single nonlinearity in the system can be solved easily to find two parameters, that is, oscillation amplitude (A) and frequency (ω) of a limit cycle.The accuracy of calculation is dependent on the accuracy of equivalent gain of the nonlinearity.The complexity of computation for multiple nonlinearities in the system is dependent on the connections of nonlinearities.Two nonlinearities subjected to the same input cannot be considered independent of each other.The relation between two nonlinearities separated by a linear transfer function can be found by evaluating magnitude and phase of the linear transfer function.Therefore, they are dependent of each other also.The complexity of computation for dependent nonlinearities is the same as the single nonlinearity.
However, nonlinearities in multivariable feedback systems are usually independent of each other.Therefore, infinite number of solutions of limit cycles satisfy the characteristic equation for phase shifts are not in the characteristic equation and the number of parameters to be found is always greater than two.The number of parameters to be found are n + 1 for an n × n multivariable feedback control system with n nonlinearities in the diagonal terms, that is, one for oscillating frequency and n for amplitudes of inputs of n nonlinearities.It needs another n − 1 simultaneously equations.The n harmonic-balance equations include phase shifts and input amplitudes of nonlinearities will be used.
The proposed method is based on the parameter-plane analyses method [18] with the characteristic equation.The nonlinearities are replaced by sinusoidal-input describing functions (SIDFs) with fundamental components [8][9][10]; that is, quasilinear gains.An infinite number of possible limit cycles found by real and imaginary parts of the characteristic equation and shown by root loci in the parameter plane first.Then, six criteria developed from the characteristic equation and harmonic-balance equations are used to find the unique solutions [18][19][20][21][22].It is deduced to check maximal frequency (ω max ) or minimal frequency (ω min ) points of root loci only for finding a stable limit cycle.The information for six criteria are (1) stable and unstable region separated by the root-loci; (2) maximal values of SIDFs of nonlinearities; (3) phase shifts of inputs of nonlinearities; (4) derivatives of equivalent gains.The proposed method is applied to nonlinear multivariable feedback systems with large transportation lags.The transportation lag is periodic function of the frequency (ω).Therefore, checking ω max or ω min points in root loci can reduce the computation effort dramatically.
This merit of the work rather than the previous work [18] is the six criteria are applied to check maximal frequency (ω max ) or minimal frequency (ω min ) points of root loci only for finding a stable limit cycle.It can reduce the computation effort dramatically.The extreme property with respect to frequency can be found by numerical equations to determine the stability, maximal and minimal frequencies rather than by the graphical method.Therefore, it has the potential to be applied to higher dimensional systems.
This paper is organized as follows: (1) in Section 2, the basic approach of the proposed method is discussed and illustrated by a 2 × 2 numerical example.Six criteria for finding unique solution are developed; (2) in Section 3, the proposed method are applied to one 2 × 2 and two 3 × 3 complicated nonlinear multivariable feedback control systems.Calculated results are verified by digital simulation verifications.It will be seen that calculated results provide accurate limit-cycle predictions of considered systems with large transportation lags.Comparisons are made also with other methods in the current literature.

The Basic Approach
Consider the nonlinear multivariable feedback system shown in Figure 1.The relation between plant transfer function matrix G(s) and nonlinearities where G(s) is the transfer matrix of the linear elements, N ( − → a ) is the transfer matrix of equivalent gains of nonlinear elements, − → r is the reference input vector, and − → a is a column vector of sinusoidal inputs to these nonlinear elements, such that where A i are amplitudes of a i , ω is the oscillating frequency, θ i are phase angles with respect to a reference input, and n is the dimension of the considered multivariable feedback system.The n linearized harmonic-balance equations governing the existence of limit cycles can be expressed as for zero reference inputs − → r .The determinant det[G(s)N ( − → a ) + I] is the characteristic equation of the considered system.It is independent of phase angle θ i and can be decomposed into two equations by taking real and imaginary parts for s = jω.The solutions need to be found for the considered nonlinear feedback control system are (A i , i = 1, 2, . . ., n) and oscillating frequency ω of the limit-cycle.
The number of parameters n + 1 to be found is larger than that of two decomposed characteristic equations.It implies that there are an infinite number of solutions satisfy the characteristic equation; that is, det[G(s)N ( − → a ) + I].It needs another n − 1 simultaneously equations.For zero inputs, (1) can be rewritten as where g i j (s) is the (i, j)th element of G(s) and n k j (a j ) is (k, j)th element of N ( − → a ).Equation ( 6) represents ith harmonic-balance equation.Let a 1 is the reference signal; that is, θ 1 = 0, then the another n − 1 simultaneous equations are derived by (6) for finding solutions (i.e., A i , i = 1, 2, . . ., n, and ω).Note that nonlinearities in the off-diagonal and ondiagonal terms are dependent for they have same input signal.For instance, nonlinearities (n i1 (a 1 ), j = 1, 2 . . ., n) are dependent for they have same input a 1 .Nonlinearities in ith feedback loop, outputs of (g ji ( jω), j = 1, 2, . . ., n) and n ii (a i ) are dependent also.Therefore, nonlinearities in the diagonal will be discussed in this paper only.
For illustration, assume that a 2 × 2 nonlinear multivariable feedback system with two single-valued nonlinearities in the diagonal terms is considered.The block diagram is shown in Figure 2.For s = jω, harmonic-balance equations of channel 1 and channel 2 are A 1 e jθ1 N 1 (a 1 )g 11 jω + A 2 e jθ2 N 2 (a 2 )g 12 jω = −A 1 e jθ1 , ( A 1 e jθ1 N 1 (a 1 )g 21 jω + A 2 e jθ2 N 2 (a 2 ) respectively.Assume that the input of N 1 is the reference input (i.e., θ 1 = 0), (7) gives Similarly, (8) gives Equating ( 9) and ( 11) gives Equation ( 13) is the characteristic equation of the considered system in ω.It is independent on the phase angle θ 2 .Equation (13) can be expressed as 1 + N 1 (a 1 )g 11 (s) + N 2 (a 2 )g 22 (s) in s-domain also.Multiplying least common multiplier (LCM) of denominators of g 11 ( jω), g 22 ( jω), and det G( jω) to ( 13), and taking real and imaginary parts of it gives the two following equations for limit-cycle evaluation: where , and E i (ω) are polynomials of ω.They will be illustrated by a simple numerical example.Equation (15) gives alternatively.Equation (16) gives alternatively.Equating ( 17) and ( 19) gives Equating ( 18) and (20) gives For specified value of frequency ω, the value of N 1 (a 1 ) can be found by solving (21); the corresponding value of N 2 (a 2 ) can be found by (17) or (19).Similarly, N 1 (a 1 ) can be found form (18) or (20) for N 2 (a 2 ) is found by (22).For a number of suitable values of ω, real solutions of N 1 (a 1 ) and N 2 (a 2 ) can be plotted in a N 1 (a 1 ) versus N 2 (a 2 ) plane, that is, parameter plane [18].Note that N 1 (a 1 ) and N 2 (a 2 ) represents equivalent gains of nonlinearities.
A useful equivalent gain expression of nonlinearity is the sinusoidal-input describing function (SIDF) where and Y (t) is the time function of nonlinearity with respect to input signal A i sin ωt.Equation ( 23) is a function of amplitude A i of sinusoidal input only.Assume the nonlinearity is symmetric, then the DC component F o is equal to zero.In general, fundamental components P 1 and R 1 are used to describe the nonlinearity [7][8][9][10].Therefore, there is a modeling error between describing function and the real nonlinear element.It affects the accuracy of limit-cycle prediction.The modeling error can be reduced by taking extra more high order harmonic components of (23).It implies that limitation of applying describing function for nonlinearity is dependent on extra harmonic components used [21].In this paper, fundamental components are used.
The modeling error will be corrected by correction formula developed by (10) and (12).Note that the nonlinearity is called "single-valued nonlinearity" for R 1 = 0 and called "doubled-valued nonlinearity" for R 1 / = 0. Inverting N 1 (a 1 ) and N 2 (a 2 ) get A 1 and A 2 , then solutions of ( 17)-( 21) can be plotted in A 1 versus A 2 plane also.
The above statement will be illustrated by a simple numerical example.Consider a 2 × 2 plant with the transfer function matrix [18] where K is the loop gains.Nonlinearities are two identical on-off relays with dead-zones having unit switching level (d) and unit height (M).Six criteria will be developed and illustrated by this numerical example, systematically.Describing functions with fundamental components of nonlinearities are where M = 1 and d = 1.The characteristic equation of the closed-loop system in s-domain is The real and imaginary parts of (27) for s = jω are For K = 2, and a number of frequency ω, simultaneous solutions (N 1 (a 1 ), N 2 (a 2 ), ω) of ( 28) and ( 29) are calculated and represented by two root loci (solid line), as shown in Figure 3.The root-loci show there are an infinite sets of possible solutions (N 1 (a 1 ), N 2 (a 2 ), ω) satisfy ( 28) and (29).But, only one set of solution (N 1 (a 1 ), N 2 (a 2 ), ω) satisfies for the considered system, that is, stable limit cycle.Other solutions are called as "unstable limit cycle".Therefore, criteria for checking the existence of a stable limit cycle must be developed.This is the motivation of the paper.The criteria for checking existence of a limit cycle will be explained by use of the illustrating example discussed above, and applied to several 2 × 2 and 3 × 3 complicated numerical examples.By use of Figure 3, six criteria of the system having a stable limit cycle are developed and explained as follows.
Criterion 1.Every point on the root loci evaluated by ( 21), (28), and (29), as shown in Figure 3, represents a set of N 1 (a 1 ), N 2 (a 2 ), and ω, which can satisfy the condition of having a limit cycle.Note that infinite possible solutions are found.
Criterion 2. A limit cycle may exist only if the values of N i (a i ) are less than the maximal gain N i (a i ) max of nonlinearities N i .Equation (26) gives the ranges of N i (a i ) are between 0 and 0.6366.Now, possible solutions of limit cycle are reduced on the segment of the root-loci between points Q 2 and Q 3 only.Criterion 3. If the root loci separate the stable and unstable regions, then a stable limit cycle may exist at the root loci.The reason is that the system will become stable (unstable) when amplitude A i increase (decrease).In other words, the system becomes stable (unstable) when the amplitude A i increase (decrease), a stable limit cycle may exist on the stability boundary, that is, on the root loci.
The stable and unstable regions are identified by the root loci found for s = +σ ± jω and s = −σ ± jω with (27).The value σ is a small positive value-0 + represents the solution found with s = +σ ± jω, and small negative value-0 − represents the solution found with s = −σ ± jω. Figure 3 shows small positive/negative values-0 + /0 − rootlocus classifies the stability of the system in the parameter plane.
The descriptions of a stable limit cycle can be expressed mathematically by the following equation:  The derivatives ∂σ/∂A i ≤ 0 represent amplitude A i of a limit cycle is increased by unknown disturbance, then real parts (σ) of characteristic roots of (27) become positive.It implies that magnitude A i will be converged by system damp.In another way, magnitude A i of a limit cycle is decreased by unknown disturbance, then real part (σ) of characteristic roots of (27) becomes negative.It implies that amplitudeA i will be diverged.Note that ∂N i (a i )/∂A i of ( 19) can be evaluated as (31) Criteria 1-3 give possible solutions of a stable limit cycle are at segment of the locus between Q 2 and Q 3 ; that is, they give ranges of frequency ω and N i (a i ).But it still has an infinite number of solutions.Criterion 4. A stable limit cycle exists only for phase angles found by ( 9) and ( 11) that are equal to each other; that is, where θ 2 {9} and θ 2 {11} represent phase angles found by ( 9) and (11), respectively.This criterion will reduce the number of possible solutions of limit cycles.Criterion 5. A stable limit cycle exists only for magnitudes found by (10) and (12) are equal; that is, Note that ( 9) and (11) give magnitudes of them are equal to unities; that is, represented by ( 10) and ( 12).However, if nonlinearities are described by the sinusoidal input describing function with fundamental components [7][8][9][10], then modeling errors of the found A i by inverting describing functions of N i (a i ) make magnitudes of ( 10) and ( 12) are not equal to unities exactly.Therefore, 1 cannot be used as criterion to find the solution except that exact description of nonlinearity is used.Naturally, M θ2 equals to unity is expected, and it is dependent on the accuracy of the nonlinearity described by SIDF.Note that a rule of thumb for expects value of M θ2 greater than 0.80 is used in this paper.Two correction equations will be developed to correct the mathematical errors of describing functions with fundamental components.Criteria 4 and 5 reduced the number of possible solutions.The next criterion will be developed for finding unique solution.
Criterion 6.The unique solution of a stable limit cycle is at the unique frequency point of the root locus; that is, the solutions of (21) for N 1 (a 1 ) are real and equal to each other.This condition gives Similar equation can be derived for N 2 (a 2 ) with (22). Figure 3 shows the maximal frequency ω max of the found upper root locus is 1.38824 rad/s at point Q 0 (2.248, 2.266); and the minimal frequency ω min of the lower root locus is 0.78881 rad/s at Point Q 1 (0.5719, 0.5691).Q 0 is a impossible solution for it violates Criteria 2 and 3. Q 1 is the unique solution satisfies Criteria 2-5 and (34).Therefore, the unique solution is found.
From the root loci shown in Figure 3, (34) can be described by a graphical rule also.It is Equation ( 35) represents the departure point ω min (point Q 1 (0.5719, 0.5691) in Figure 3) of the root locus with respect  to the frequency ω, or the approaching point ω max (point Q 0 (2.248, 2.266) in Figure 3) of root locus with respect to the frequency ω.
If the solution satisfies all six criteria for a stable limit cycle, then a stable limit cycle will exist.Table 1 gives calculated results of point Q 1 (0.5719, 0.5691).Two sets of (A 1 , A 2 ) satisfy found N 1 (a 1 ) and N 2 (a 2 ).First set of (A 1 , A 2 ) = (1.9039,1.8885) is the desired solutions.Second set of (A 1 , A 2 ) = (1.175,1.1785) is impossible for its ∂N 1 (a 1 )/∂A 1 and ∂N 2 (a 2 )/∂A 2 violate Criterion 3. Calculated results for Q 3 are given in Table 1 also for illustrating it is an unstable limit cycle.Note that (A i ) are found from (26), that is, describing function of the relay with dead band; therefore, M θ2 found by (10) or (12) are usually not equal to unities for mathematical errors of the nonlinearities.By multiplying a scaling factor S k to left and right side of (10) An approximate formulation for S k is The error of S k M θ2 − 1 is less than 0.5% for 0.9 < M θ2 < 1.1 (1.2% for 0.85 < M θ2 < 1.15).Equations ( 36) and (37) give the modified values (A im ) of (A i ) are Using (38), the modified values are A 1m = 1.9706 and A 2m = 1.8229.Figure 4 shows simulation verification result of the considered system in which gives A 1 = 1.9858,A 2 = 1.8549, ω = 0.7883 rad/s, and θ 2 = −69.97• .They give that calculated results (A 1m = 1.9706 and A 2m = 1.8229) corrected by (38) give accurate prediction of the stable limit cycle.
If the loop gain K is an adjustable parameter, then the minimal value of K just having a stable limit cycle can be found by the same evaluating procedures and criteria.The found value is 1.7915.The root locus for K = 1.7915 is shown in Figure 5.It implies that there will have no intersection between root locus and constant N 1 (a 1 ) max , and N 2 (a 2 ) max lines.The system is asymptotically stable for K is less than 1.7915.Therefore, the proposed method can be used for designing nonlinear multivariable feedback control  systems also, that is, not only for analyses.The comparisons with other methods [6] for minimal K are given in Table 2.
The procedures for finding a stable limit cycle have been developed and illustrated by a simple example with two "single-valued nonlinearities".If the nonlinearities are "two double-valued nonlinearities"; that is, R 1 / = 0 in (23).Then, the characteristic ( 14) is rewritten as 1 + N 1r (a 1 ) + jN 1i (a 1 ) g 11 (s) + N 2r (a 2 ) + jN 2i (a 2 ) g 22 (s) + N 1r (a 1 ) + jN 1i (a 1 ) N 2r (a 2 ) + jN 2i (a 2 ) × g 11 (s)g 22 (s) − g 12 (s)g 21 (s) = 0, (39) in s-domain, and the real and imaginary parts of (39) with s = jω are Then, parameter analyses in the N 1 (a 1 ) versus N 2 (a 2 ) plane are replaced by those of in the A 1 versus A 2 plane.For example, nonlinearities are replaced by two double-valued nonlinearities shown in Figure 6.The describing functions are where M = 1, d = 0.9 and p = 1.1.Figure 7 shows the root-loci of the system with new two-valued nonlinearities.Six criteria for finding a stable limit cycle have been developed for nonlinear multivariable feedback control systems with single-and double-valued nonlinearities.The same analyzing and design procedures with six checking criteria will be applied to following 2 × 2 and 3 × 3 complicated nonlinear multivariable feedback systems.Note that six criteria are deduced to check the ω max or ω min point of root loci which satisfies Criteria 2 to 5.This reduce the computing effort dramatically.

Numerical Examples
Example 1.Consider a nonlinear multivariable system with transfer function matrix [23] 12.8e −s 16.7s + 1 18.9e −3s 21s + 1 6.6e −7s 10.9s + 1 19.4e −3s 14.4s + 1 Two nonlinearities are shown in Figure 9. Equation (43) gives the considered system is a large transportation lag system.Similar to the procedure stated in Section 2, the found root loci are shown in Figure 10.There are two ω max (Q 6 , Q 8 ) and two ω min (Q 5 , Q 7 ) points of root loci.They represent possible solutions of the stable limit cycle.But only the Q 5 is the solution for it satisfies Criterions 2 to 5. The simulation verification is shown in Figure 11.Comparison of the calculated and simulated results is given in Table 3.It can be seen that calculated results give accurate prediction of the considered system.Note that the transportation lag is a periodic function of frequency ω.Therefore, Figure 10 gives four maximal ad minimal frequency points of root loci.The illustrating example stated in Section 2 gives one maximal and one minimal frequency points (Figure 2) only for it is a 2 × 2 system and has no transportation lag.Example 1 gives the proposed method give an effect way to find the N 1 (a 1 ) exact solution.The developed criteria for 2 × 2 systems are extended to following 3 × 3 systems.
Example 2. Consider a 3 × 3 multivariable process [24] given by where K is the loop gain.There are three relay nonlinearities in the diagonal terms.The deadband (d) and magnitude (M) of each nonlinearity are 0.5 and 1.0, respectively.The describing functions of them are given in (26).The harmonic-balance equations of the system are given by N 1 (a 1 )g 11 (s) N 2 (a 2 )g 12 (s) N 3 (a 3 )g 13 (s) N 1 (a 1 )g 21 (s) N 2 (a 2 )g 22 (s) N 3 (a 3 )g 23 (s) N 1 (a 1 )g 31 (s) N 2 (a 2 )g 32 (s) N 3 (a 3 )g 33 (s) 5 where D g (s) represents the determinant of the transfer function matrix G(s).For a specified value of N 3 (a 3 ), the characteristic equation is function of N 1 (a 1 ), N 2 (a 2 ), and ω only.Therefore, the same analyzing procedures for 2 × 2 nonlinear multivariable systems described by ( 17)-( 22) and six criteria can be applied.Figure 12(a) shows parameter analyses of several constant-N 3 (a 3 ) loci between N 3 (a 3 ) max = 1.273 and N 3 (a 3 ) = 0.00 for K = 1.0.Each constant-N 3 (a 3 ) locus shows the maximal frequency ω max which is given in Table 4. Intersecting points between the dashdot line and constant-N 3 (a 3 ) loci give ω max of constant-N 3 (a 3 ) loci. Figure 12 There are three nonlinearities in the diagonal.Figure 14 shows nonlinearities.Q 10 represents existence of a stable limit cycle, that is, ω = 0.3593 rad/s, N 1 (a 1 ) = 0.6578, N 2 (a 2 ) = 1.6919, and N 3 (a 3 ) = 0.910.The corresponding amplitudes are A 1 = 1.835,A 2 = 0.7735, and A 3 = 1.2215.They are found by inverting the describing functions.Figure 16 shows simulation verifications give A 1 = 1.965,A 2 = 0.8357, A 3 = 1.249, and ω = 0.3541 rad/s.It shows that calculated results give accurate prediction of a stable limit cycle.

Discussions
From analyses and simulated results of Examples 1 to 3, procedures for finding the stable limit cycle can be simplified by only checking one or two points of root loci whether they satisfy six criteria or not.Those points are minimal or maximal frequency point (ω min , ω max ) of the root loci of the characteristic equation.Such that the proposed method reduces the computational efforts largely.The simplified procedures are given below: (1) plotting root loci from the characteristic equation in the parameter plane (or space); (2) finding the maximal and minimal frequency points (ω in , ω max ) of root loci; (3) checking points (ω in , ω max ) satisfy Criteria 2 to 5 or not.

Conclusions
In this paper, a practical method for limit-cycle predictions in nonlinear multivariable feedback control system has been presented and found to be much simpler than other methods given in the current literature.It has been shown that calculated results give accurate predictions for consider nonlinear multivariable feedback control systems.

Figure 3 :
Figure 3: Root loci of limit cycles in the parameter plane.

Figure 10 :Figure 11 :
Figure 10: Root loci analyses of limit cycles of Example 1.

Table 1 :
Calculated results of a stable (Point Q 1 ) and a unstable limit cycle (Point Q 3 ).

Table 2 :
The gains K for just having a limit cycle.

Table 3 :
Calculated and simulated results of Example 1.

Table 4 :
The maximal frequency point ω max of each constant-N 3 (a 3 ) locus of Example 2.
(b)shows ω max locus versus N 3 (a 3 ).It gives the maximal frequency with respect to N 3 (a 3 ) is ω max = 2.061 rad/s.The corresponding values of N i (a i ) are the point Q 9 (0.9968, 0.9968) shown in Figure13.It is the unique solution of the stable limit cycle.The found A i are (A 1 , A 2 , A 3 ) = (1.151,1.151,1.151).They are found by inverting the describing functions.Figure13shows simulation verification results in which gives (A 1 , A 2 , A 3 ) = (1.264,1.264, 1.264) and ω = 2.058 rad/s.It can be seen that calculated results are quite closed to simulated results for this nonlinear 3 × 3 multivariable feedback control system.