Abstract

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.

1. 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 [110]. It has been shown that for multivariable systems, over arbitrary ranges of amplitudes (𝐴𝑖), frequency (𝜔) and phases (𝜃𝑖), 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 𝐴𝑖, 𝜔, and 𝜃𝑖, 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 single-output (SISO) nonlinear feedback control systems [1117]. Therefore, single nonlinearity in the system can be solved easily to find two parameters, that is, oscillation amplitude (𝐴) 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 𝑛+1 for an 𝑛×𝑛 multivariable feedback control system with 𝑛 nonlinearities in the diagonal terms, that is, one for oscillating frequency and 𝑛 for amplitudes of inputs of 𝑛 nonlinearities. It needs another 𝑛1 simultaneously equations. The 𝑛 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 [810]; 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 [1822]. 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.

2. The Basic Approach

Consider the nonlinear multivariable feedback system shown in Figure 1. The relation between plant transfer function matrix 𝐺(𝑠) and nonlinearities 𝑁(𝑎) is 𝑦=𝐺(𝑠)𝑁𝑎𝑟𝑎,(1) where 𝐺(𝑠) is the transfer matrix of the linear elements, 𝑁(𝑎) is the transfer matrix of equivalent gains of nonlinear elements, 𝑟 is the reference input vector, and 𝑎 is a column vector of sinusoidal inputs to these nonlinear elements, such that𝑎𝑖=𝐴𝑖sin𝜔𝑡+𝜃𝑖,(𝑖=1,2,,𝑛),(2) where 𝐴𝑖 are amplitudes of 𝑎𝑖, 𝜔 is the oscillating frequency, 𝜃𝑖 are phase angles with respect to a reference input, and 𝑛 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[𝐺(𝑠)𝑁𝑎+𝐼]𝑎|||𝑠=𝑗𝜔=0,(3) for zero reference inputs 𝑟. The determinant det[𝐺(𝑠)𝑁(𝑎)+𝐼] is the characteristic equation of the considered system. It is independent of phase angle 𝜃𝑖 and can be decomposed into two equations by taking real and imaginary parts for 𝑠=𝑗𝜔. The solutions need to be found for the considered nonlinear feedback control system are (𝐴𝑖, 𝑖=1,2,,𝑛) and oscillating frequency 𝜔 of the limit-cycle. The number of parameters 𝑛+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[𝐺(𝑠)𝑁(𝑎)+𝐼]. It needs another 𝑛1 simultaneously equations. For zero inputs, (1) can be rewritten as 𝑛𝑗=1𝑛𝑘=1𝑔𝑖𝑘(𝑠)𝑛𝑘𝑗𝑎𝑗𝑎𝑗=𝑎𝑖,(4)𝑛𝑗=1𝑛𝑘=1𝑔𝑖𝑘(𝑠)𝑛𝑘𝑗𝑎𝑗𝐴𝑗𝑒𝑗(𝜔+𝜃𝑗)=𝐴𝑖𝑒𝑗(𝜔+𝜃𝑖),(5)𝑛𝑗=1𝑛𝑘=1𝑔𝑖𝑘(𝑠)𝑛𝑘𝑗𝑎𝑗𝐴𝑗𝑒𝑗𝜃𝑗=𝐴𝑖𝑒𝑗𝜃𝑖,(6) where 𝑔𝑖𝑗(𝑠) is the (𝑖,𝑗)th element of 𝐺(𝑠) and 𝑛𝑘𝑗(𝑎𝑗) is (𝑘,𝑗)th element of 𝑁(𝑎). Equation (6) represents 𝑖th harmonic-balance equation. Let 𝑎1 is the reference signal; that is, 𝜃1=0, then the another 𝑛1 simultaneous equations are derived by (6) for finding solutions (i.e., 𝐴𝑖, 𝑖=1,2,,𝑛, and 𝜔).

Note that nonlinearities in the off-diagonal and on-diagonal terms are dependent for they have same input signal. For instance, nonlinearities (𝑛𝑖1(𝑎1), 𝑗=1,2,𝑛) are dependent for they have same input 𝑎1. Nonlinearities in 𝑖th feedback loop, outputs of (𝑔𝑗𝑖(𝑗𝜔),𝑗=1,2,,𝑛) and 𝑛𝑖𝑖(𝑎𝑖) 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 𝑠=𝑗𝜔, harmonic-balance equations of channel 1 and channel 2 are𝐴1𝑒𝑗𝜃1𝑁1𝑎1𝑔11(𝑗𝜔)+𝐴2𝑒𝑗𝜃2𝑁2𝑎2𝑔12(𝑗𝜔)=𝐴1𝑒𝑗𝜃1,𝐴(7)1𝑒𝑗𝜃1𝑁1𝑎1𝑔21(𝑗𝜔)+𝐴2𝑒𝑗𝜃2𝑁2𝑎2𝑔22(𝑗𝜔)=𝐴2𝑒𝑗𝜃2,(8) respectively. Assume that the input of 𝑁1 is the reference input (i.e., 𝜃1=0), (7) gives𝑒𝑗𝜃2𝐴=11+𝑁1𝑎1𝑔11(𝑗𝜔)𝐴2𝑁2𝑎2𝑔12(𝑗𝜔),(9)||𝑒𝑗𝜃2||=𝐴1𝐴2||||1+𝑁1𝑎1𝑔11(𝑗𝜔)𝑁2𝑎2𝑔12(||||𝑗𝜔)𝑀𝜃2=1.(10) Similarly, (8) gives𝑒𝑗𝜃2𝐴=1𝑁1𝑎1𝑔21(𝑗𝜔)𝐴21+𝑁2𝑎2𝑔22(𝑗𝜔),(11)||𝑒𝑗𝜃2||=𝐴1𝐴2||||𝑁1𝑎1𝑔21(𝑗𝜔)1+𝑁2𝑎2𝑔22(||||𝑗𝜔)𝑀𝜃2=1.(12) Equating (9) and (11) gives1+𝑁1𝑎1𝑔11(𝑗𝜔)+𝑁2𝑎2𝑔22(𝑗𝜔)+𝑁1𝑎1𝑁2𝑎2𝑔11(𝑗𝜔)𝑔22(𝑗𝜔)𝑔12(𝑗𝜔)𝑔21(𝑗𝜔)=0.(13) Equation (13) is the characteristic equation of the considered system in 𝜔. It is independent on the phase angle 𝜃2. Equation (13) can be expressed as1+𝑁1𝑎1𝑔11(𝑠)+𝑁2𝑎2𝑔22(𝑠)+𝑁1𝑎1𝑁2𝑎2𝑔11(𝑠)𝑔22(𝑠)𝑔12(𝑠)𝑔21(𝑠)=0,(14) in 𝑠-domain also. Multiplying least common multiplier (LCM) of denominators of 𝑔11(𝑗𝜔), 𝑔22(𝑗𝜔), and det𝐺(𝑗𝜔) to (13), and taking real and imaginary parts of it gives the two following equations for limit-cycle evaluation:𝐵1(𝜔)+𝑁1𝑎1𝐶1(𝜔)+𝑁2𝑎2𝐷1(𝜔)+𝑁1𝑎1𝑁2𝑎2𝐸1(𝐵𝜔)=0,(15)2(𝜔)+𝑁1𝑎1𝐶2(𝜔)+𝑁2𝑎2𝐷2(𝜔)+𝑁1𝑎1𝑁2𝑎2𝐸2(𝜔)=0,(16) where 𝐵𝑖(𝜔), 𝐶𝑖(𝜔), 𝐷𝑖(𝜔), and 𝐸𝑖(𝜔) are polynomials of 𝜔. They will be illustrated by a simple numerical example. Equation (15) gives𝑁2𝑎2𝐵=1(𝜔)+𝑁1𝑎1𝐶1(𝜔)𝐷1(𝜔)+𝑁1𝑎1𝐸1(𝑁𝜔),(17)1𝑎1𝐵=1(𝜔)+𝑁2𝑎2𝐷1(𝜔)𝐶1(𝜔)+𝑁2𝑎2𝐸1(𝜔),(18) alternatively. Equation (16) gives𝑁2𝑎2𝐵=2(𝜔)+𝑁1𝑎1𝐶2(𝜔)𝐷2(𝜔)+𝑁1𝑎1𝐸2(𝑁𝜔),(19)1𝑎1𝐵=2(𝜔)+𝑁2𝑎2𝐷2(𝜔)𝐶2(𝜔)+𝑁2𝑎2𝐸2(𝜔),(20) alternatively. Equating (17) and (19) gives𝐶2(𝜔)𝐸1(𝜔)𝐶1(𝜔)𝐸2𝑁(𝜔)1𝑎12+𝐶2(𝜔)𝐷1(𝜔)+𝐵2(𝜔)𝐸1(𝜔)𝐶1(𝜔)𝐷2(𝜔)𝐵1(𝜔)𝐸2𝑁(𝜔)1𝑎1+𝐵2(𝜔)𝐷1(𝜔)𝐵1(𝜔)𝐷2(𝜔)=0.(21) Equating (18) and (20) gives𝐷2(𝜔)𝐸1(𝜔)𝐷1(𝜔)𝐸2𝑁(𝜔)2𝑎22+𝐶1(𝜔)𝐷2(𝜔)𝐶2(𝜔)𝐷1(𝜔)+𝐵2(𝜔)𝐸1(𝜔)𝐵1(𝜔)𝐸2𝑁(𝜔)2𝑎2+𝐵2(𝜔)𝐶1(𝜔)𝐵1(𝜔)𝐶2(𝜔)=0.(22) For specified value of frequency 𝜔, the value of 𝑁1(𝑎1) can be found by solving (21); the corresponding value of 𝑁2(𝑎2) can be found by (17) or (19). Similarly, 𝑁1(𝑎1) can be found form (18) or (20) for 𝑁2(𝑎2) is found by (22). For a number of suitable values of 𝜔, real solutions of 𝑁1(𝑎1) and 𝑁2(𝑎2) can be plotted in a 𝑁1(𝑎1) versus 𝑁2(𝑎2) plane, that is, parameter plane [18]. Note that 𝑁1(𝑎1) and 𝑁2(𝑎2) represents equivalent gains of nonlinearities.

A useful equivalent gain expression of nonlinearity is the sinusoidal-input describing function (SIDF)𝑁𝑖𝑎𝑖=𝐹𝑜+𝑛=1𝑃𝑛+𝑗𝑅𝑛=𝑁𝑖𝑟𝑎𝑖+𝑗𝑁1𝑖𝑎𝑖,(23) where𝐹𝑜=1𝐴𝑖02𝜋𝑃𝑌(𝑡)𝑑(𝜔𝑡),𝑛=1𝐴𝑖02𝜋𝑅𝑌(𝑡)cos(𝑛𝜔𝑡)𝑑(𝜔𝑡),𝑛=1𝐴𝑖02𝜋𝑌(𝑡)sin(𝑛𝜔𝑡)𝑑(𝜔𝑡),(24) and 𝑌(𝑡) is the time function of nonlinearity with respect to input signal 𝐴𝑖sin𝜔𝑡. Equation (23) is a function of amplitude 𝐴𝑖 of sinusoidal input only. Assume the nonlinearity is symmetric, then the DC component 𝐹𝑜 is equal to zero. In general, fundamental components 𝑃1 and 𝑅1 are used to describe the nonlinearity [710]. 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 𝑅1=0 and called “doubled-valued nonlinearity” for 𝑅10. Inverting 𝑁1(𝑎1) and 𝑁2(𝑎2) get 𝐴1 and 𝐴2, then solutions of (17)–(21) can be plotted in 𝐴1 versus 𝐴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]1𝐺(𝑠)=𝐾𝑠(𝑠+1)20.3𝑠(𝑠+1)20.21𝑠(𝑠+1)𝑠(𝑠+1)2,(25) where 𝐾 is the loop gains. Nonlinearities are two identical on-off relays with dead-zones having unit switching level (𝑑) and unit height (𝑀). Six criteria will be developed and illustrated by this numerical example, systematically. Describing functions with fundamental components of nonlinearities are 𝑁𝑖𝑎𝑖=4𝑀𝜋𝐴𝑖𝑑12𝐴2𝑖1/2,𝐴𝑖𝑑,𝑖=1,2,(26) where 𝑀=1 and 𝑑=1. The characteristic equation of the closed-loop system in 𝑠-domain is𝑠6+4𝑠5+6𝑠4+4𝑠3+𝑠2+𝐾𝑁1𝑎1𝑠3+2𝑠2+𝑠+𝐾𝑁2𝑎2𝑠3+2𝑠2+𝑠+𝐾2𝑁1𝑎1𝑁2𝑎2(0.006𝑠+1.06)=0.(27) The real and imaginary parts of (27) for 𝑠=𝑗𝜔 are𝜔6+6𝜔4𝜔2+𝐾𝑁1𝑎12𝜔2+𝐾𝑁2𝑎22𝜔2+𝐾2𝑁1𝑎1𝑁2𝑎2(1.06)=0,(28)4𝜔54𝜔3+𝐾𝑁1𝑎1𝜔3+𝜔+𝐾𝑁2𝑎2𝜔3+𝜔+𝐾2𝑁1𝑎1𝑁2𝑎2(0.06𝜔)=0.(29) For 𝐾=2, and a number of frequency 𝜔, simultaneous solutions (𝑁1(𝑎1),𝑁2(𝑎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 (𝑁1(𝑎1),𝑁2(𝑎2),𝜔) satisfy (28) and (29). But, only one set of solution (𝑁1(𝑎1),𝑁2(𝑎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 𝑁1(𝑎1), 𝑁2(𝑎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 𝑁𝑖(𝑎𝑖) are less than the maximal gain 𝑁𝑖(𝑎𝑖)max of nonlinearities 𝑁𝑖. Equation (26) gives the ranges of 𝑁𝑖(𝑎𝑖) are between 0 and 0.6366. Now, possible solutions of limit cycle are reduced on the segment of the root-loci between points 𝑄2 and 𝑄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 𝐴𝑖 increase (decrease). In other words, the system becomes stable (unstable) when the amplitude 𝐴𝑖 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 𝑠=+𝜎±𝑗𝜔 and 𝑠=𝜎±𝑗𝜔 with (27). The value 𝜎 is a small positive value-0+ represents the solution found with 𝑠=+𝜎±𝑗𝜔, and small negative value-0 represents the solution found with 𝑠=𝜎±𝑗𝜔. Figure 3 shows small positive/negative values-0+/0 root-locus 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:𝜕𝜎𝜕𝐴𝑖=𝜕𝜎𝜕𝑁𝑖𝑎𝑖𝜕𝑁𝑖𝑎𝑖𝜕𝐴𝑖<0,𝑖=1,2.(30) The derivatives 𝜕𝜎/𝜕𝐴𝑖0 represent amplitude 𝐴𝑖 of a limit cycle is increased by unknown disturbance, then real parts (𝜎) of characteristic roots of (27) become positive. It implies that magnitude 𝐴𝑖 will be converged by system damp. In another way, magnitude 𝐴𝑖 of a limit cycle is decreased by unknown disturbance, then real part (𝜎) of characteristic roots of (27) becomes negative. It implies that amplitude𝐴𝑖 will be diverged. Note that 𝜕𝑁𝑖(𝑎𝑖)/𝜕𝐴𝑖 of (19) can be evaluated as 𝜕𝑁𝑖𝑎𝑖𝜕𝐴𝑖=4𝑀𝜋𝐴𝑖2𝑑12𝐴𝑖21/2+𝑑2𝐴𝑖2𝑑12𝐴𝑖21/2.(31) Criteria 13 give possible solutions of a stable limit cycle are at segment of the locus between 𝑄2 and 𝑄3; that is, they give ranges of frequency 𝜔 and 𝑁𝑖(𝑎𝑖). 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, 𝜃2{9}𝜃2{11}=0,(32) 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, 𝑀𝜃2{10}𝑀𝜃2{12}=0.(33) 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 [710], then modeling errors of the found 𝐴𝑖 by inverting describing functions of 𝑁𝑖(𝑎𝑖) make magnitudes of (10) and (12) are not equal to unities exactly. Therefore, 𝑀𝜃2{10}=𝑀𝜃2{12}=1 cannot be used as criterion to find the solution except that exact description of nonlinearity is used. Naturally, 𝑀𝜃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 𝑀𝜃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 𝑁1(𝑎1) are real and equal to each other. This condition gives 𝐶2(𝜔)𝐷1(𝜔)+𝐵2(𝜔)𝐸1(𝜔)𝐶1(𝜔)𝐷2(𝜔)𝐵1(𝜔)𝐸2(𝜔)2𝐶42(𝜔)𝐸1(𝜔)𝐶1(𝜔)𝐸2×𝐵(𝜔)2(𝜔)𝐷1(𝜔)𝐵1(𝜔)𝐷2(𝜔)=0.(34) Similar equation can be derived for 𝑁2(𝑎2) with (22). Figure 3 shows the maximal frequency 𝜔max of the found upper root locus is 1.38824 rad/s at point 𝑄0(2.248,2.266); and the minimal frequency 𝜔min of the lower root locus is 0.78881 rad/s at Point 𝑄1(0.5719,0.5691). 𝑄0 is a impossible solution for it violates Criteria 2 and 3. 𝑄1 is the unique solution satisfies Criteria 25 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𝜕𝑁𝑖𝑎𝑖𝜕𝜔=0.(35) Equation (35) represents the departure point 𝜔min (point 𝑄1(0.5719,0.5691) in Figure 3) of the root locus with respect to the frequency 𝜔, or the approaching point 𝜔max (point 𝑄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 𝑄1(0.5719,0.5691). Two sets of (𝐴1,𝐴2) satisfy found 𝑁1(𝑎1) and 𝑁2(𝑎2). First set of (𝐴1,𝐴2) = (1.9039,1.8885) is the desired solutions. Second set of (𝐴1,𝐴2) = (1.175,1.1785) is impossible for its 𝜕𝑁1(𝑎1)/𝜕𝐴1 and 𝜕𝑁2(𝑎2)/𝜕𝐴2 violate Criterion 3. Calculated results for 𝑄3 are given in Table 1 also for illustrating it is an unstable limit cycle. Note that (𝐴𝑖) are found from (26), that is, describing function of the relay with dead band; therefore, 𝑀𝜃2 found by (10) or (12) are usually not equal to unities for mathematical errors of the nonlinearities. By multiplying a scaling factor 𝑆𝑘 to left and right side of (10) for |𝑒𝑗𝜃2|=1, then (10) becomes𝑆𝑘𝐴1𝐴2||||1+𝑁1𝑎1𝑔11(𝑗𝜔)𝑁2𝑎2𝑔12(||||𝑗𝜔)𝑆𝑘𝑀𝜃2=1.(36) An approximate formulation for 𝑆𝑘 is 𝑆𝑘1+1𝑀𝜃2/211𝑀𝜃2=/21.50.5𝑀𝜃20.5+0.5𝑀𝜃2.(37) The error of 𝑆𝑘𝑀𝜃21 is less than 0.5% for 0.9<𝑀𝜃2<1.1 (1.2% for 0.85<𝑀𝜃2<1.15). Equations (36) and (37) give the modified values (𝐴𝑖𝑚) of (𝐴𝑖) are 𝐴1𝑚=𝐴11+1𝑀𝜃22=𝐴11.50.5𝑀𝜃2,𝐴2𝑚=𝐴211𝑀𝜃22=𝐴20.5+0.5𝑀𝜃2.(38) Using (38), the modified values are 𝐴1𝑚=1.9706 and 𝐴2𝑚=1.8229. Figure 4 shows simulation verification result of the considered system in which gives 𝐴1=1.9858, 𝐴2=1.8549, 𝜔=0.7883 rad/s, and 𝜃2=69.97. They give that calculated results (𝐴1𝑚=1.9706 and 𝐴2𝑚=1.8229) corrected by (38) give accurate prediction of the stable limit cycle.

If the loop gain 𝐾 is an adjustable parameter, then the minimal value of 𝐾 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 𝐾=1.7915 is shown in Figure 5. It implies that there will have no intersection between root locus and constant 𝑁1(𝑎1)max, and 𝑁2(𝑎2)max lines. The system is asymptotically stable for 𝐾 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 𝐾 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, 𝑅10 in (23). Then, the characteristic (14) is rewritten as 𝑁1+1𝑟𝑎1+𝑗𝑁1𝑖𝑎1𝑔11+𝑁(𝑠)2𝑟𝑎2+𝑗𝑁2𝑖𝑎2𝑔22(+𝑁𝑠)1𝑟𝑎1+𝑗𝑁1𝑖𝑎1𝑁2𝑟𝑎2+𝑗𝑁2𝑖𝑎2×𝑔11(𝑠)𝑔22(𝑠)𝑔12(𝑠)𝑔21(𝑠)=0,(39) in 𝑠-domain, and the real and imaginary parts of (39) with 𝑠=𝑗𝜔 are 𝐵1𝑁(𝜔)+1𝑟𝑎1𝐶1(𝜔)𝑁1𝑖𝑎1𝐶2+𝑁(𝜔)2𝑟𝑎2𝐷1(𝜔)𝑁2𝑖𝑎2𝐷2(+𝑁𝜔)1𝑟𝑎1𝑁2𝑟𝑎2𝑁1𝑖𝑎1𝑁2𝑖𝑎2𝐸1𝑁(𝜔)1𝑟𝑎1𝑁2𝑖𝑎2+𝑁1𝑖𝑎1𝑁2𝑟𝑎2𝐸2(𝐵𝜔)=0,(40)2𝑁(𝜔)+1𝑟𝑎1𝐶2(𝜔)+𝑁1𝑖𝑎1𝐶1+𝑁(𝜔)2𝑟𝑎2𝐷2(𝜔)+𝑁2𝑖𝑎2𝐷1(+𝑁𝜔)1𝑟𝑎1𝑁2𝑟𝑎2𝑁1𝑖𝑎1𝑁2𝑖𝑎2𝐸2+𝑁(𝜔)1𝑟𝑎1𝑁2𝑖𝑎2+𝑁1𝑖𝑎1𝑁2𝑟𝑎2𝐸1(𝜔)=0.(41) Then, parameter analyses in the 𝑁1(𝑎1) versus 𝑁2(𝑎2) plane are replaced by those of in the 𝐴1 versus 𝐴2 plane. For example, nonlinearities are replaced by two double-valued nonlinearities shown in Figure 6. The describing functions are 𝑁𝑖𝑎𝑖=𝑁𝑖𝑟𝑎𝑖+𝑗𝑁𝑖𝑖𝑎𝑖=2𝑀𝜋𝐴𝑖𝑑12𝐴𝑖21/2+𝑝12𝐴𝑖21/2𝑗2𝑀(𝑝𝑑)𝜋𝐴𝑖2,𝐴𝑖𝑝,(42) where 𝑀=1, 𝑑=0.9 and 𝑝=1.1. Figure 7 shows the root-loci of the system with new two-valued nonlinearities. Six checking criteria gives Point 𝑄4(2.164.2.163) represents the unique solution (𝜔min=0.7476 rad/sec, 𝐴1=2.164, 𝐴2=2.163, 𝜃=71.59, and 𝑀𝜃2=0.9172). The corrected amplitudes by (38) with 𝑀𝜃2=0.9172 are 𝐴1𝑚=2.2536 and 𝐴2𝑚=2.0735. The simulation verification result is shown in Figure 8, in which gives (𝜔=0.7506 rad/sec, 𝐴1=2.2592, 𝐴2=2.0923 and 𝜃=67.82). It can be seen that calculated results give accurate prediction for double-valued nonlinearities also.

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.

3. Numerical Examples

Example 1. Consider a nonlinear multivariable system with transfer function matrix [23] 𝐺(𝑠)=12.8𝑒𝑠16.7𝑠+118.9𝑒3𝑠21𝑠+16.6𝑒7𝑠10.9𝑠+119.4𝑒3𝑠14.4𝑠+1.(43) 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 (𝑄6,𝑄8) and two 𝜔min (𝑄5,𝑄7) points of root loci. They represent possible solutions of the stable limit cycle. But only the 𝑄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 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 𝐺(𝑠)=𝐾2𝑒𝑠10𝑠+11.5𝑒𝑠𝑒𝑠+1𝑠𝑠+11.5𝑒𝑠𝑒𝑠+1𝑠𝑠+12𝑒𝑠𝑒10𝑠+1𝑠𝑠+12𝑒𝑠10𝑠+11.5𝑒𝑠𝑠+1,(44) where 𝐾 is the loop gain. There are three relay nonlinearities in the diagonal terms. The deadband (𝑑) and magnitude (𝑀) 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 𝐴1𝑒𝑗𝜃1𝐴2𝑒𝑗𝜃2𝐴3𝑒𝑗𝜃3=𝑁1𝑎1𝑔11(𝑠)𝑁2𝑎2𝑔12(𝑠)𝑁3𝑎3𝑔13𝑁(𝑠)1𝑎1𝑔21(𝑠)𝑁2𝑎2𝑔22(𝑠)𝑁3𝑎3𝑔23𝑁(𝑠)1𝑎1𝑔31(𝑠)𝑁2𝑎2𝑔32(𝑠)𝑁3𝑎3𝑔33(×𝑠)𝐴1𝑒𝑗𝜃1𝐴2𝑒𝑗𝜃2𝐴3𝑒𝑗𝜃3,(45)1+𝑁1𝑎1𝑔11(𝑠)𝑁2𝑎2𝑔12(𝑠)𝑁3𝑎3𝑔13𝑁(𝑠)1𝑎1𝑔21(𝑠)1+𝑁2𝑎2𝑔22(𝑠)𝑁3𝑎3𝑔23𝑁(𝑠)1𝑎1𝑔31(𝑠)𝑁2𝑎2𝑔32(𝑠)1+𝑁3𝑎3𝑔33(×𝐴𝑠)1𝑒𝑗𝜃1𝐴2𝑒𝑗𝜃2𝐴3𝑒𝑗𝜃3=0,(46) where 𝑔𝑖𝑗(𝑠) is the (𝑖,𝑗)th element of 𝐺(𝑠). Equation (46) gives the characteristic equation of the system as 1+𝑁3𝑎3𝑔33+𝑔(𝑠)11(𝑠)+𝑁3𝑎3𝑔11(𝑠)𝑔33(𝑠)𝑔13(𝑠)𝑔31(𝑁𝑠)1𝑎1+𝑔22(𝑠)+𝑁3𝑎3𝑔22(𝑠)𝑔33(𝑠)𝑔32(𝑠)𝑔23𝑁(𝑠)2𝑎2+𝑔11(𝑠)𝑔22(𝑠)𝑔12(𝑠)𝑔21(𝑠)+𝑁3𝑎3𝐷𝑔(𝑠)×𝑁1𝑎1𝑁2𝑎2=0,(47) where 𝐷𝑔(𝑠) represents the determinant of the transfer function matrix 𝐺(𝑠). For a specified value of 𝑁3(𝑎3), the characteristic equation is function of 𝑁1(𝑎1), 𝑁2(𝑎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-𝑁3(𝑎3) loci between 𝑁3(𝑎3)max=1.273 and 𝑁3(𝑎3)=0.00 for 𝐾=1.0. Each constant-𝑁3(𝑎3) locus shows the maximal frequency 𝜔max which is given in Table 4. Intersecting points between the dash-dot line and constant-𝑁3(𝑎3) loci give 𝜔max of constant-𝑁3(𝑎3) loci. Figure 12(b) shows 𝜔max locus versus 𝑁3(𝑎3). It gives the maximal frequency with respect to 𝑁3(𝑎3) is 𝜔max=2.061 rad/s. The corresponding values of 𝑁𝑖(𝑎𝑖) are the point 𝑄9(0.9968,0.9968) shown in Figure 13. It is the unique solution of the stable limit cycle. The found 𝐴𝑖 are (𝐴1,𝐴2,𝐴3)=(1.151,1.151,1.151). They are found by inverting the describing functions. Figure 13 shows simulation verification results in which gives (𝐴1,𝐴2,𝐴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.

Example 3. Consider a 3×3 multivariable feedback control system with the transfer function matrix [25] 𝐺1(𝑠)=10119𝑒5𝑠21.7𝑠+140𝑒5𝑠337𝑠+12.1𝑒5𝑠10𝑠+177𝑒5𝑠50𝑠+176.7𝑒3𝑠28𝑠+15𝑒5𝑠10𝑠+193𝑒5𝑠50𝑠+136.7𝑒5𝑠166𝑠+1103.3𝑒4𝑠25𝑠+1.(48) There are three nonlinearities in the diagonal. Figure 14 shows nonlinearities. Figure 15(a) shows root loci of possible solutions of limit cycles in the 𝑁1(𝑎1) versus 𝑁2(𝑎2) plane for specified values of 𝑁3(𝑎3). The 𝜔max-locus shows connections of each 𝜔max point of constant-𝑁3(𝑎3) locus. The maximal value of the 𝜔max-locus shown in Figure 15(b) gives 𝜔max=0.3593 rad/s; that is, point 𝑄10. The point 𝑄10 represents existence of a stable limit cycle, that is, 𝜔=0.3593 rad/s, 𝑁1(𝑎1)=0.6578, 𝑁2(𝑎2)=1.6919, and 𝑁3(𝑎3)=0.910. The corresponding amplitudes are 𝐴1=1.835, 𝐴2=0.7735, and 𝐴3=1.2215. They are found by inverting the describing functions. Figure 16 shows simulation verifications give 𝐴1=1.965, 𝐴2=0.8357, 𝐴3=1.249, and 𝜔=0.3541 rad/s. It shows that calculated results give accurate prediction of a stable limit cycle.

4. 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 (𝜔𝑖𝑛,𝜔max) of root loci;(3)checking points (𝜔𝑖𝑛,𝜔max) satisfy Criteria 2 to 5 or not.

5. 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.