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

Authors: Tain-Sou Tsay
Journal: Journal of Control Science and Engineering
 (2011)

---

## 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 one2×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.


---

## Body

## 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 [1–10]. It has been shown that for multivariable systems, over arbitrary ranges of amplitudes (Ai), 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 Ai, ω, 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 single-output (SISO) nonlinear feedback control systems [11–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 aren+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–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–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 Section2, 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 Figure1. The relation between plant transfer function matrix G(s) and nonlinearities N(a⃗) is (1)y⃗=G(s)N(a⃗)(r⃗-a⃗),
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(2)ai=Aisin(ωt+θi),(i=1,2,…,n),
where Ai are amplitudes of ai, ω 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(3)[G(s)N(a⃗)+I]a⇀|s=jω=0⃗,
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 (Ai, 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 (4)∑j=1n[∑k=1ngik(s)nkj(aj)]aj=-ai,(5)∑j=1n[∑k=1ngik(s)nkj(aj)]Ajej(ω+θj)=-Aiej(ω+θi),(6)∑j=1n[∑k=1ngik(s)nkj(aj)]Ajejθj=-Aiejθi,
where gij(s) is the (i,j)th element of G(s) and  nkj(aj) is (k,j)th element of N(a⃗). Equation (6) represents ith harmonic-balance equation. Let a1 is the reference signal; that is, θ1=0, then the another n-1 simultaneous equations are derived by (6) for finding solutions (i.e., Ai, i=1,2,…,n, and ω).Figure 1
A general block diagram of nonlinear multivariable feedback control systems.Note that nonlinearities in the off-diagonal and on-diagonal terms are dependent for they have same input signal. For instance, nonlinearities (ni1(a1), j=1,2…,n) are dependent for they have same input a1. Nonlinearities in ith feedback loop, outputs of (gji(jω),j=1,2,…,n) and nii(ai) are dependent also. Therefore, nonlinearities in the diagonal will be discussed in this paper only.For illustration, assume that a2×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(7)A1ejθ1N1(a1)g11(jω)+A2ejθ2N2(a2)g12(jω)=-A1ejθ1,(8)A1ejθ1N1(a1)g21(jω)+A2ejθ2N2(a2)g22(jω)=-A2ejθ2,
respectively. Assume that the input of N1 is the reference input (i.e., θ1=0), (7) gives(9)ejθ2=-A1[1+N1(a1)g11(jω)]A2N2(a2)g12(jω),(10)|ejθ2|=(A1A2)|1+N1(a1)g11(jω)N2(a2)g12(jω)|≡Mθ2=1.
Similarly, (8) gives(11)ejθ2=-A1N1(a1)g21(jω)A2[1+N2(a2)g22(jω)],(12)|ejθ2|=(A1A2)|N1(a1)g21(jω)1+N2(a2)g22(jω)|≡Mθ2=1.
Equating (9) and (11) gives(13)1+N1(a1)g11(jω)+N2(a2)g22(jω)+N1(a1)N2(a2)[g11(jω)g22(jω)-g12(jω)g21(jω)]=0.
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(14)1+N1(a1)g11(s)+N2(a2)g22(s)+N1(a1)N2(a2)[g11(s)g22(s)-g12(s)g21(s)]=0,
in s-domain also. Multiplying least common multiplier (LCM) of denominators of g11(jω), g22(jω), and detG(jω) to (13), and taking real and imaginary parts of it gives the two following equations for limit-cycle evaluation:(15)B1(ω)+N1(a1)C1(ω)+N2(a2)D1(ω)+N1(a1)N2(a2)E1(ω)=0,(16)B2(ω)+N1(a1)C2(ω)+N2(a2)D2(ω)+N1(a1)N2(a2)E2(ω)=0,
where Bi(ω), Ci(ω), Di(ω), and Ei(ω) are polynomials of ω. They will be illustrated by a simple numerical example. Equation (15) gives(17)N2(a2)=-B1(ω)+N1(a1)C1(ω)D1(ω)+N1(a1)E1(ω),(18)N1(a1)=-B1(ω)+N2(a2)D1(ω)C1(ω)+N2(a2)E1(ω),
alternatively. Equation (16) gives(19)N2(a2)=-B2(ω)+N1(a1)C2(ω)D2(ω)+N1(a1)E2(ω),(20)N1(a1)=-B2(ω)+N2(a2)D2(ω)C2(ω)+N2(a2)E2(ω),
alternatively. Equating (17) and (19) gives(21)[C2(ω)E1(ω)-C1(ω)E2(ω)]N1(a1)2+[C2(ω)D1(ω)+B2(ω)E1(ω)-C1(ω)D2(ω)-B1(ω)E2(ω)]N1(a1)+[B2(ω)D1(ω)-B1(ω)D2(ω)]=0.
Equating (18) and (20) gives(22)[D2(ω)E1(ω)-D1(ω)E2(ω)]N2(a2)2+[C1(ω)D2(ω)-C2(ω)D1(ω)+B2(ω)E1(ω)-B1(ω)E2(ω)]N2(a2)+[B2(ω)C1(ω)-B1(ω)C2(ω)]=0.
For specified value of frequency ω, the value of N1(a1) can be found by solving (21); the corresponding value of N2(a2) can be found by (17) or (19).  Similarly, N1(a1) can be found form  (18) or (20) for N2(a2) is found by (22). For a number of suitable values of ω, real solutions of N1(a1) and N2(a2) can be plotted in a N1(a1) versus N2(a2) plane, that is, parameter plane [18]. Note that N1(a1) and N2(a2) represents equivalent gains of nonlinearities.Figure 2
Block diagram of a2×2 nonlinear multivariable feedback control system.A useful equivalent gain expression of nonlinearity is the sinusoidal-input describing function (SIDF)(23)Ni(ai)=Fo+∑n=1∞(Pn+jRn)=Nir(ai)+jN1i(ai),
where(24)Fo=1Ai∫02πY(t)d(ωt),Pn=1Ai∫02πY(t)cos⁡(nωt)d(ωt),Rn=1Ai∫02πY(t)sin(nωt)d(ωt),
and Y(t) is the time function of nonlinearity with respect to input signal Aisinωt. Equation (23) is a function of amplitude Ai of sinusoidal input only. Assume the nonlinearity is symmetric, then the DC component Fo is equal to zero. In general, fundamental components P1 and R1 are used to describe the nonlinearity [7–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 R1=0 and called “doubled-valued nonlinearity” for R1≠0. Inverting N1(a1) and N2(a2) get A1 and A2, then solutions of  (17)–(21) can be plotted in A1 versus A2 plane also.The above statement will be illustrated by a simple numerical example. Consider a2×2 plant with the transfer function matrix [18](25)G(s)=K[1s(s+1)20.3s(s+1)2-0.2s(s+1)1s(s+1)2],
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 (26)Ni(ai)=4MπAi(1-d2Ai2)1/2,Ai≥d,i=1,2,
where M=1 and d=1. The characteristic equation of the closed-loop system in s-domain is(27)s6+4s5+6s4+4s3+s2+KN1(a1)(s3+2s2+s)+KN2(a2)(s3+2s2+s)+K2N1(a1)N2(a2)(0.006s+1.06)=0.
The real and imaginary parts of (27) for s=jω are(28)-ω6+6ω4-ω2+KN1(a1)(-2ω2)+KN2(a2)(-2ω2)+K2N1(a1)N2(a2)(1.06)=0,(29)4ω5-4ω3+KN1(a1)(-ω3+ω)+KN2(a2)(-ω3+ω)+K2N1(a1)N2(a2)(0.06ω)=0.
For K=2, and a number of frequency ω, simultaneous solutions (N1(a1),N2(a2),ω) 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 (N1(a1),N2(a2),ω) satisfy (28) and (29). But, only one set of solution (N1(a1),N2(a2),ω) 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.Figure 3
Root loci of limit cycles in the parameter plane.The criteria for checking existence of a limit cycle will be explained by use of the illustrating example discussed above, and applied to several2×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 N1(a1), N2(a2), 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 ofNi(ai) are less than the maximal gain Ni(ai)max of nonlinearities Ni. Equation (26) gives the ranges of Ni(ai) are between 0 and 0.6366. Now, possible solutions of limit cycle are reduced on the segment of the root-loci between points Q2 and Q3 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 amplitudeAi increase (decrease). In other words, the system becomes stable (unstable) when the amplitude Ai 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 fors=+σ±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− 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:(30)∂σ∂Ai=(∂σ∂Ni(ai))(∂Ni(ai)∂Ai)<0,i=1,2.
The derivatives ∂σ/∂Ai≤0 represent amplitude Ai of a limit cycle is increased by unknown disturbance, then real parts (σ) of characteristic roots of (27) become positive. It implies that magnitude Ai will be converged by system damp. In another way, magnitude Ai of a limit cycle is decreased by unknown disturbance, then real part (σ) of characteristic roots of (27) becomes negative. It implies that amplitudeAi will be diverged. Note that ∂Ni(ai)/∂Ai of (19) can be evaluated as (31)∂Ni(ai)∂Ai=4MπAi2[-(1-d2Ai2)1/2+d2Ai2(1-d2Ai2)-1/2].
Criteria 1–3 give possible solutions of a stable limit cycle are at segment of the locus between Q2 and Q3; that is, they give ranges of frequency ω and Ni(ai). 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,
(32)θ2{9}-θ2{11}=0,
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,
(33)Mθ2{10}-Mθ2{12}=0.
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–10], then modeling errors of the found Ai by inverting describing functions of Ni(ai) make magnitudes of  (10) and (12) are not equal to unities exactly. Therefore, Mθ2{10}=Mθ2{12}=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 N1(a1) are real and equal to each other. This condition gives
(34)[C2(ω)D1(ω)+B2(ω)E1(ω)-C1(ω)D2(ω)-B1(ω)E2(ω)]2-4[C2(ω)E1(ω)-C1(ω)E2(ω)]×[B2(ω)D1(ω)-B1(ω)D2(ω)]=0.
Similar equation can be derived for N2(a2) with (22). Figure 3 shows the maximal frequency  ωmax of the found upper root locus is 1.38824 rad/s at point Q0(2.248,2.266); and the minimal frequency ωmin of the lower root locus is 0.78881 rad/s at Point Q1(0.5719,0.5691). Q0 is a impossible solution for it violates Criteria 2 and 3. Q1 is the unique solution satisfies Criteria 2–5 and (34). Therefore, the unique solution is found.From the root loci shown in Figure3, (34) can be described by a graphical rule also. It is(35)∂Ni(ai)∂ω=0.
Equation (35) represents the departure point ωmin (point Q1(0.5719,0.5691) in Figure 3) of the root locus with respect to the frequency ω, or the approaching point ωmax (point Q0(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. Table1 gives calculated results of point Q1(0.5719,0.5691). Two sets of (A1,A2) satisfy found N1(a1) and N2(a2). First set of (A1,A2) = (1.9039,1.8885) is the desired solutions. Second set of (A1,A2) = (1.175,1.1785) is impossible for its ∂N1(a1)/∂A1 and  ∂N2(a2)/∂A2 violate Criterion 3. Calculated results for Q3 are given in Table 1 also for illustrating it is an unstable limit cycle. Note that (Ai) 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 Sk to left and right side of  (10) for |e-jθ2|=1, then (10) becomes(36)Sk(A1A2)|1+N1(a1)g11(jω)N2(a2)g12(jω)|≡SkMθ2=1.
An approximate formulation for Sk is (37)Sk≈1+(1-Mθ2)/21-(1-Mθ2)/2=1.5-0.5Mθ20.5+0.5Mθ2.
The error of SkMθ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 (Aim) of (Ai) are (38)A1m=A1[1+(1-Mθ2)2]=A1(1.5-0.5Mθ2),A2m=A2[1-(1-Mθ2)2]=A2(0.5+0.5Mθ2).
Using  (38), the modified values are A1m=1.9706 and A2m=1.8229. Figure 4 shows simulation verification result of the considered system in which gives A1=1.9858, A2=1.8549, ω=0.7883 rad/s, and θ2=-69.97°. They give that calculated results (A1m=1.9706 and A2m=1.8229) corrected by  (38)  give accurate prediction of the stable limit cycle.Table 1
Calculated results of a stable (PointQ1) and a unstable limit cycle (Point Q3).
PointN1(a1)N2(a2)ωA1A2∂N1(a1)/∂A1∂N2(a2)/∂A2θ2{9}θ2{11}Mθ2Q10.57190.56910.78881.9041.889−0.185-0.1849−70.39°−70.39°0.931.1751.179+0.788+0.762−70.39°−70.39°0.92Q30.63660.53130.80021.4141.1360+1.149−88.49°−88.49°1.121.4142.1120−0.179−88.49°−88.49°0.61Figure 4
Time responses of the illustrating example.If the loop gainK 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 N1(a1)max, and N2(a2)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.Table 2
The gainsK for just having a limit cycle.
MethodsGainKProposed method1.7915Aizerman conjecture1.79Hirsch plot1.25Mee plot1.50Digital simulation1.787Figure 5
Root locus analyses forK=1.7915.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,R1≠0 in (23). Then, the characteristic (14) is rewritten as (39)1+[N1r(a1)+jN1i(a1)]g11(s)+[N2r(a2)+jN2i(a2)]g22(s)+[N1r(a1)+jN1i(a1)][N2r(a2)+jN2i(a2)]×[g11(s)g22(s)-g12(s)g21(s)]=0,
in s-domain, and the real and imaginary parts of (39) with s=jω are (40)B1(ω)+[N1r(a1)C1(ω)-N1i(a1)C2(ω)]+[N2r(a2)D1(ω)-N2i(a2)D2(ω)]+{[N1r(a1)N2r(a2)-N1i(a1)N2i(a2)]E1(ω)-[N1r(a1)N2i(a2)+N1i(a1)N2r(a2)]E2(ω)}=0,(41)B2(ω)+[N1r(a1)C2(ω)+N1i(a1)C1(ω)]+[N2r(a2)D2(ω)+N2i(a2)D1(ω)]+{[N1r(a1)N2r(a2)-N1i(a1)N2i(a2)]E2(ω)+[N1r(a1)N2i(a2)+N1i(a1)N2r(a2)]E1(ω)}=0.
Then, parameter analyses in the N1(a1) versus N2(a2) plane are replaced by those of in the A1 versus A2 plane. For example, nonlinearities are replaced by two double-valued nonlinearities shown in Figure 6. The describing functions are (42)Ni(ai)=Nir(ai)+jNii(ai)=2MπAi[(1-d2Ai2)1/2+(1-p2Ai2)1/2]-j2M(p-d)πAi2,Ai≥p,
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 checking criteria gives Point Q4(2.164.2.163) represents the unique solution (ωmin=0.7476 rad/sec, A1=2.164, A2=2.163, θ=-71.59°, and Mθ2=0.9172). The corrected amplitudes by (38) with Mθ2=0.9172 are A1m=2.2536 and A2m=2.0735. The simulation verification result is shown in Figure 8, in which gives (ω=0.7506 rad/sec, A1=2.2592, A2=2.0923 and θ=-67.82°). It can be seen that calculated results give accurate prediction for double-valued nonlinearities also.Figure 6
The double-valued nonlinearity.Figure 7
Root loci of limit cycles with two doubled-valued nonlinearities.Figure 8
Time responses with two double-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 following2×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]
(43)G(s)=[12.8e-s16.7s+118.9e-3s21s+16.6e-7s10.9s+119.4e-3s14.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 (Q6,Q8) and two ωmin (Q5,Q7) points of root loci. They represent possible solutions of the stable limit cycle. But only the Q5 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.Table 3
Calculated and simulated results of Example1.
Osci. Freq. (rad/s)Channel no. 1Channel no. 2θ2 (deg)Mθ2Calculation0.4875Nic(ai)0.45410.2929−53.30.95Aic1.09612.1390Simulation0.4836Ais1.06072.2454−54.4Figure 9
Nonlinearities of Example1.Figure 10
Root loci analyses of limit cycles of Example1.Figure 11
Time responses of Example1.Example 2.
Consider a3×3 multivariable process [24] given by
(44)G(s)=K[-2e-s10s+11.5e-ss+1e-ss+11.5e-ss+1e-ss+1-2e-s10s+1e-ss+1-2e-s10s+11.5e-ss+1],
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
(45)[A1ejθ1A2ejθ2A3ejθ3]=[N1(a1)g11(s)N2(a2)g12(s)N3(a3)g13(s)N1(a1)g21(s)N2(a2)g22(s)N3(a3)g23(s)N1(a1)g31(s)N2(a2)g32(s)N3(a3)g33(s)]×[-A1ejθ1-A2ejθ2-A3ejθ3],(46)[1+N1(a1)g11(s)N2(a2)g12(s)N3(a3)g13(s)N1(a1)g21(s)1+N2(a2)g22(s)N3(a3)g23(s)N1(a1)g31(s)N2(a2)g32(s)1+N3(a3)g33(s)]×[A1ejθ1A2ejθ2A3ejθ3]=0⃗,
where gij(s) is the (i,j)th element of G(s). Equation (46) gives the characteristic equation of the system as
(47)1+N3(a3)g33(s)+{g11(s)+N3(a3)[g11(s)g33(s)-g13(s)g31(s)]}N1(a1)+{g22(s)+N3(a3)[g22(s)g33(s)-g32(s)g23(s)]}N2(a2)+{g11(s)g22(s)-g12(s)g21(s)+N3(a3)Dg(s)}×N1(a1)N2(a2)=0,
where Dg(s) represents the determinant of the transfer function matrix G(s). For a specified value of N3(a3), the characteristic equation is function of N1(a1), N2(a2), 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-N3(a3) loci between N3(a3)max=1.273 and N3(a3)=0.00 for K=1.0. Each constant-N3(a3) locus shows the maximal frequency ωmax which is given in Table 4. Intersecting points between the dash-dot line and constant-N3(a3) loci give ωmax of constant-N3(a3) loci. Figure 12(b) shows ωmax locus versus N3(a3). It gives the maximal frequency with respect to N3(a3) is ωmax=2.061 rad/s. The corresponding values of Ni(ai) are the point Q9(0.9968,0.9968) shown in Figure 13. It is the unique solution of the stable limit cycle. The found Ai are (A1,A2,A3)=(1.151,1.151,1.151). They are found by inverting the describing functions. Figure 13 shows simulation verification results in which gives (A1,A2,A3)=(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.Table 4
The maximal frequency pointωmax of each constant-N3(a3) locus of Example 2.
N3(a3)1.2731.201.151.100.99680.900.700.300.100.00ωmax (rad/s)2.0402.0492.0542.0582.0612.0592.0542.0442.0412.039(a) Root loci analyses of limit cycles of Example2. (b) ωmax versus N3(a3) of Example 2.
(a)(b)Figure 13
Time responses of Example2.Example 3.
Consider a3×3 multivariable feedback control system with the transfer function matrix [25]
(48)G(s)=110[119e-5s21.7s+140e-5s337s+12.1e-5s10s+177e-5s50s+176.7e-3s28s+15e-5s10s+193e-5s50s+1-36.7e-5s166s+1103.3e-4s25s+1].
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 N1(a1) versus N2(a2) plane for specified values of N3(a3). The ωmax-locus shows connections of each ωmax point of constant-N3(a3) locus. The maximal value of the ωmax-locus shown in Figure 15(b) gives ωmax=0.3593 rad/s; that is, point Q10. The point Q10 represents existence of a stable limit cycle, that is, ω=0.3593 rad/s, N1(a1)=0.6578, N2(a2)=1.6919, and N3(a3)=0.910. The corresponding amplitudes are A1=1.835, A2=0.7735, and A3=1.2215. They are found by inverting the describing functions. Figure 16 shows simulation verifications give A1=1.965, A2=0.8357, A3=1.249, and ω=0.3541 rad/s. It shows that calculated results give accurate prediction of a stable limit cycle.Figure 14
Nonlinearities of Example3.(a) Root loci analyses of limit cycles of Example3.  (b) ωmax versus N3(a3) of Example 3.
(a)(b)Figure 16
Time responses of Example3.


## 4. Discussions

From analyses and simulated results of Examples1 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.


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



---

*Source: 169848-2011-05-30.xml*