Fractional-order lead networks to avoid limit cycle in control loops with dead zone and plant servo system

2.1 Formalization and properties of the controller

Analytic formulas have made the design of a classical, integer-order lead network easy, mainly because the mathematical expressions relate the controller parameters to the frequency of the maximum phase lead. In the frequency domain, this relation is important because the maximum lead has to fall exactly at the gain crossover frequency of the compensated system. Unfortunately, it is not easy to find similar relationships for general multiple-stage, shifted in frequency, fractional-order lead networks. However, if the analysis is restricted to two-stage SFLNs, the next arguments show that the frequency ( u m 12 ) of the maximum lead is the real, positive root of a linear, irreducible, cubic equation. Moreover, the coefficient values of the Cardano’s equation depend on the SFLN parameters. To prove the above statement, consider the transfer functions of the first and second stages that are, respectively, H 1 ( s ) = 1 + τ s 1 + τ Δ s ν 1 and H 2 ( s ) = 1 + τ F s 1 + τ F Δ s ν 2 , where ν 1 = ν 2 = ν is set to obtain equal but shifted stages. Here, τ and τ Δ are the numerator and denominator time constants (in seconds) of the first stage, F represents the shift of the second stage with respect to the first, and the noninteger (fractional) exponent ν , with 0 < ν < 1 , is the so-called fractional order. Then, the frequency response of the two-stage SFLN is given as follows:

where the variable u = ω τ is the dimensionless angular frequency and ω is the angular frequency (in radians/seconds). The minimum value of Δ in a classical lead network can be fixed to 0.05 because it is limited by the physical construction of the compensator and by the limitations of the maximum phase lead and high-frequency gain that are provided by the network [28]. Namely, if a maximum phase lead greater than 6 0 ° is needed, two cascaded lead networks with moderate values of Δ are advised. Moreover, high-frequency noise signals are amplified by 1 ∕ Δ , while a unitary amplification applies to low-frequency control signals. Then, Δ > 0.07 is often recommended [29]. The usual choice is Δ ∈ ( 0.05 , 0.25 ) [28,29]. Then, this article makes the common choice Δ = 0.1 , which is recommended for lead networks [29], for both stages of the proposed SFLN.

Finally, the parameter F , with 0 < F < 1 , shifts the phase diagram of H 2 ( j u F ) with reference to the position of the phase diagram of H 1 ( j u ) on the u -axis [30–32]. Clearly, the fractional-order lead networks, H 1 ( j u ) and H 2 ( j u F ) , provide their maximum phase leads at the frequencies u m 1 = 1 ∕ Δ and u m 2 = 1 ∕ ( F Δ ) , respectively. These values coincide with the frequencies at which the maximum phase lead occurs in a network with the same base but with an integer-order exponent ν = 1 . They are also the geometric mean of the two corner frequencies, i.e., the middle point, in the logarithmic scale, between the zero and pole corner frequency. So, according to De Moivre’s theorem, the phase of H 1 ( j u ) at u = u m 1 has the value Φ m 1 = ∠ H 1 ( j u m 1 ) given as follows:

Similarly, the second stage provides the maximum value Φ m 2 = ∠ H 2 ( j u m 2 F ) at the frequency u = u m 2 = 1 ∕ ( F Δ ) , with

Moreover, if two generic frequencies u 1 and u 2 , with u 1 < u 2 , satisfy the relationship u 1 = u 2 F , then similar arguments easily show that ∠ H 1 ( j u 1 ) = ∠ H 1 ( j u 2 F ) = ∠ H 2 ( j u 2 ) .

Now the frequency u m 12 , where H 12 ( j u ) provides the maximum phase lead, say Φ m 12 , must be determined. Namely, the pair ( u m 12 , Φ m 12 ) plays a major role in an SFLN design finalized to both avoiding the limit cycle and maintaining the stabilization mode. To determine u m 12 , the following equation must be solved:

Since ∠ H 1 ( j u ) = ν [ tan − 1 ( u ) − tan − 1 ( u Δ ) ] and ∠ H 2 ( j u F ) = ν [ tan − 1 ( u F ) − tan − 1 ( u F Δ ) ] , we obtain:

After simple but tedious calculations, Eq. (5) and the positions α = F Δ , β = ( 1 + Δ 2 ) , and γ = ( 1 − F + F 2 ) lead to:

So putting u 2 = q in Eq. (6), an obvious choice of symbols leads to

with

The coefficients of Eq. (8) only depend on F and Δ and can be easily determined. Since c 1 , c 2 , and c 3 are real, Eq. (7) has at least one real root. Moreover, by Descartes’ Rule of signs, the number of positive real roots of Eq. (7) is either equal to the sign changes in the sequence { 1 , c 1 , c 2 , c 3 } or is less than this number by a positive even integer [33]. The values F ∈ ( 0.1 , 1 ) and Δ = 0.1 used by this article (see next section), give c 1 > 0 , c 2 < 0 , and c 3 < 0 , so that the number of sign variations is 1. Hence, Eq. (7) has exactly one real positive root (say q 1 ). The value of q 1 can be obtained as a numerical solution of Eq. (7) or as a result of Cardano’s analytical procedure [34] and leads to the frequency of the maximum phase lead. In this article, q 1 is obtained as a numerical solution of Eq. (7). Hence, u m 12 = q 1 .

Now comments on further characteristics of H 12 ( j u ) are in order. The following analysis tries to fix the fractional order ν on the basis of preliminary design decisions. Thus, to provide the phase lead the control loop requires for compensation, the value of ν is first selected. For example, in the case of the next section, a good first choice is ν = 0.7 . Moreover, the parameters F and Δ influence u m 12 and Φ m 12 = ∠ H 12 ( j u m 12 ) , which are the main frequency-domain measures assessing performance. Hence, if Δ is fixed (e.g., Δ = 0.1 ), each value given to F will identify a member Φ 12 ( u , F ) = ∠ H 12 ( j u ) = f F ( u ) of a function family. As shown in Figure 2, each phase-frequency plot is marked by an F -value. Moreover, the lower values of F make the phase curves flatter around their maximum, Φ m 12 , whereas the more the value of F rises, the more Φ m 12 increases and u m 12 decreases (Figure 2 and Table 1).

Figure 2

Phase and amplitude Bode diagrams of H 12 ( j u ) as a function of u and for different values of F , with ν = 0.7 and Δ = 0.1 .

In the second relation, each value assigned to Δ leads to a one-to-one function u m 12 ( F , Δ ) = f Δ ( F ) and to a curve of Figure 3. Lower values of Δ lead to higher values of u m 12 corresponding to assigned F .

Figure 3

Value of the dimensionless angular frequency u m 12 ( F , Δ ) for the maximum phase lead as function of F and with ν = 0.7 .

Figures 2 and 3 with formulas (6), (7), and (8) are the basis of an efficient grapho-analytical design procedure of H 12 ( j u ) . Of course, different pairs of the above figures can be obtained for different values of the fractional order ν .

2.2 Realization of the controller

The model of the SFLN is irrational and infinite-dimensional. However, a finite-dimensional, rational transfer function realization can be obtained in three steps. The first one introduces the rational transfer function approximating the irrational derivative operator x ν :

where x ∈ C (complex plane) and 0 < ν < 1 . Note that the coefficients a N − i and b N − i depend on ν and can be obtained by one of the existing methods (see [35] and references therein). However, this article proposes the following simple, effective closed-formulas for a N − i and b N − i :

where

define the Pochhammer functions, with ( ν + N + 1 ) ( 0 ) = ( ν − N ) ( 0 ) = 1 (see [36,37] for more details). In the second step, the transformation

is introduced in Eq. (9). Then, some elementary operations lead to:

Furthermore, a routine-based algebra yields

where A N − k ( B N − k ) depends on a N − i ( b N − i ) according to the following linear combination:

with

and where m 1 = max { 0 , k + i − N } and m 2 = min { i , k } because of the restriction imposed by the binomial coefficients: 0 ≤ m ≤ i and 0 ≤ k − m ≤ N − i [30].

In the third step, the rational transfer function of the second stage H 2 ( s ) can be immediately obtained by replacing s with F s :

In conclusion, H 12 ( s ) ≈ G 1 ( s ) G 2 ( s ) .

Figures 2 and 4 show the amplitude and phase plots of an SFLN with some values assigned to ν and Δ . Figure 4 shows a fourth-order ( N = 4 ) approximation by a dash-dotted line. Note that the discrepancy between the irrational network H 12 ( j ω ) and its approximation G 1 ( j ω ) G 2 ( j ω ) in the Bode phase diagram is not important. Namely, it is negligible for u < u m 12 and gives no problem for u = u m 12 , where the maximum phase Φ m 12 is achieved and used for design purpose (Section 3). The phase difference only consists in few degrees for u > u m 12 , but it does not significantly affect the design that requires specifications at u m 12 and a certain flatness in the immediate range around u m 12 .

Figure 4

Phase and amplitude Bode diagrams of H 12 ( j u ) (continuous line) and G 1 ( j u ) G 2 ( j u ) (dashed line).

As for the phase plot concerns, curve C ( F = 0.7 ) in Figures 2 and 4 shows a sufficiently flat phase behavior in a sufficiently large frequency range both for the irrational controller H 12 ( s ) and for its rational transfer function realization. The phase diagram is not really “flat” (constant phase) but its variation is limited, especially with reference to the corresponding phase change provided by integer-order networks.

3.1 The system with limit cycle

This section proposes the positional control of a common type of plant, i.e., a servomotor, as case study (Figure 5). The transfer functions of the controller H 12 ( s ) and servomotor G p ( s ) are the constituent linear parts of the control loop. Preliminarily, however, Figure 5(a) is considered without H 12 ( s ) to point out the presence of limit cycles in the autonomous system (with r = 0 ) including a memoryless time-invariant nonlinearity. So let

with K p = 10 , T 1 = 0.15 s, and T 2 = 0.95 s. Without H 12 ( s ) , at the phase and gain crossover frequencies ( ω p c P = 2.65 rad/s and ω g c P = 3.01 rad/s), respectively, the values of the gain and phase margins of G p ( j ω ) are poor (i.e., G M P = 0.77 and P M P = − 5.0 2 ° ).

Figure 5

Block schemes describing the FOC of a servomotor with a dead zone. (a) The uncompensated system including the DZNL with μ = 1 and a = 0.8 and the plant transfer function G p ( s ) and (b) the controlled system: H 12 ( s ) is the controller with two, shifted in frequency, fractional-order lead networks, DZNL is the dead zone nonlinearity, and P is the plant.

Moreover, a DZNL in the loop is described as follows:

where μ = 1 and a = 0.8 are here used, w is the input to the dead zone, z is the output from the dead zone, and Figure 5(a) explains the meaning of the symbols μ (slope) and a (half of dead zone).

Now the describing function approach can be successfully applied because, by replacing s → j ω in Eq. (19), G p ( j ω ) shows the necessary higher harmonics filtering capacity. Hence, the harmonic balance equation yields

where W is the amplitude of the sinusoidal input w ( t ) = W sin ( ω t ) applied to the nonlinear element and

is the describing function of the DZNL element [38,39]. Since the dead zone is a single-valued nonlinearity, N ( W ) and its negative inverse V ( W ) are real and their plots lie on the real axis (Figure 6). Moreover, the intersection point P 1 of the curve G p ( j ω ) and the curve V ( W ) corresponds to a limit cycle that is defined by the pair ( W c , ω c ) , with N ( W c ) = ( T 1 + T 2 ) ∕ ( K p T 1 T 2 ) and ω c = 1 ∕ T 1 T 2 . The amplitude W c corresponds to P 1 on the V ( W ) curve, and the same point P 1 on the plot of G p ( j ω ) identifies the limit cycle frequency, ω c (Figure 6). In this case, it holds ( W c , ω c ) = ( 4.42 , 2.65 ) . The curve V ( W ) plays the role of a “critical curve” that is comparable to the critical point ( − 1 + j 0 ) in the simplified Nyquist stability criterion [38], according to Eq. (21) and the simplified scheme in Figure 5(a). Then, if no control is applied, the time-domain control system response exhibits a limit cycle oscillation, as shown in Figure 7. Finally, since the plot of V ( W ) meets the curve of G p ( j ω ) from left to right of the cross-point P 1 , then a practical rule states that the limit cycle is unstable [38,40]. As remarked in [39], the importance of the Nyquist plot and Nyquist criterion in control theory made the graphical implementation of the describing function method a popular tool available to control engineers, when facing nonlinear control systems [39].

Figure 6

Nyquist plot of G p ( j ω ) (red solid line); negative inverse of the describing function of the DZNL (blue solid line); point P 1 ( W c , ω c ) : limit cycle position; Nyquist plot of H 12 ( j ω ) G p ( j ω ) (black solid line): no limit cycle occurs.

Figure 7

Output oscillation (limit cycle) when no control is applied.

3.2 Controller design

The robust limit cycle avoidance is committed both to the Gain Margin (GM) and Phase Margin (PM) of G H P ( j ω ) = H 12 ( j ω ) G p ( j ω ) . Namely, one or another between the PM and GM is a suitable design metric because G H P ( j ω ) is minimum-phase with a unique relation between the amplitude and phase. Hence, the controller H 12 ( j ω ) must shape G H P ( j ω ) to make the compensated Nyquist plot be at a suitable distance from P 1 ( W c , ω c ) .

However, some difficulties arise. Namely, combining H 12 ( j ω ) in cascade with G p ( j ω ) shifts the gain crossover frequency of G H P ( j ω ) , say ω g c H P , on the right of the gain crossover frequency ω g c P of G p ( j ω ) . Hence, H 12 ( j ω ) must provide additional leads to compensate for the larger phase lag introduced by G p ( j ω ) at ω = ω g c H P . A well-established practice suggests an increment of ε = 2 0 ° ∼ 2 5 ° to the specified desired margin (i.e., to P M s = 4 0 ° ∼ 5 0 ° ) [29]. The choice of the controller parameters of H 12 ( j ω ) is in six steps.

1. First step: setting Δ . To begin with, the value of Δ is chosen for H 12 ( j ω ) . In common, integer-order, phase lead networks, the value Δ = 0.1 reaches a good compromise between the phase request and an acceptable sensitivity to high-frequency noise [40]. Hence, the same value is chosen for H 12 ( j ω ) , as it was shown and anticipated in Section 2.

2. Second step: setting ν and F to obtain a specified phase margin. The values of the parameters F and ν must guarantee both a phase margin of the system G H P ( j ω ) = H 12 ( j ω ) G p ( j ω ) , namely P M H P = P M s ≥ 4 0 ° ∼ 5 0 ° , and a nearly flat-phase interval. Hence, to compensate a large phase lag at the (still unknown) gain crossover frequency ω g c H P of the open-loop transfer function G H P ( j ω ) , the controller H 12 ( j ω ) should provide a maximum phase lead Φ m 12 = P M H P − P M P + ε = 7 0 ° ∼ 7 5 ° . As for ν , since the fractional order deeply affects the phase leads provided by H 12 ( j ω ) (Section 2.2), the required values can suitably be obtained with ν ∈ ( 0.7 , 0.8 ) . According to Figure 2, which is drawn for ν = 0.7 , the curves A-D, corresponding to some values of F in the range F ∈ ( 0.3 , 0.9 ) , guarantee the required leads with nearly flat phase patterns around the frequency of the maximum phase lead. Since flatter curves correspond to lower leads, flatness and phase lead values must be balanced. Hence, after some attempts, F = 0.7 and ν = 0.7 are chosen so as to obtain both the maximum lead Φ m 12 = 76 . 3 ° in correspondence of u m 12 = 3.78 and a flat phase in a given frequency interval (Figure 3 and Table 1).

3. Third step: computing the shift of the gain crossover frequency. The time constant τ of the first stage H 1 ( j ω ) = 1 + j ω τ 1 + j ω τ Δ of H 12 ( j ω ) is chosen so that the maximum lead of the controller will contribute wholly to the phase margin. Hence, let τ ω g c H P = u m 12 and use the parameter h > 1 to express ω g c H P in terms of ω g c P as follows: ω g c H P = h ω g c P . Then, τ can be written as:

   (23) τ ( h ) = u m 12 h ω g c P .

   To obtain out of the circular difficulty of determining the (unknown) time constant τ in terms of the (unknown) parameter h , a simple numerical code is used to find the element h = h ¯ giving both the smallest value E ( h ¯ ) = min h ∣ 1 − G H P ( j h ω g c P ) ∣ ≅ 0 and its position in the vector E ( h ) . The result gives E ( h ¯ ) = 0.003 for h ¯ = 2.21 and ω g c H P = h ¯ ω g c P = 6.67 rad/s.

4. Fourth step: setting time constants of the two stages. The outcome h ¯ of the third step gives τ = u m 12 ∕ ω g c H P = 0.57 s and τ 1 = F τ = 0.40 s.

5. Fifth step: checking the phase margin. The resulting irrational network H 12 ( j ω ) leads to ∠ G H P ( j ω g c H P ) = − 132 ° and to the phase margin P M H P = 4 8 ° . By using the fourth-order, rational transfer function realization G 1 ( j ω ) G 2 ( j ω ) , the phase margin P M G P = 45.1 ° is obtained. A better approximation to P M H P can be easily obtained with rational transfer functions of higher order N (Section 2.2).

6. Sixth step: determination of the phase crossover frequency and gain margin. To determine the gain margin, the vector ∣ ∠ H 12 ( j ω ) G p ( j ω ) + π ∣ is defined and its component with minimum value is determined. Based on the previously determined parameters for H 12 ( j ω ) , the minimum returns ∠ H 12 ( j ω p c H P ) G p ( j ω p c H P ) = − 3.15 radians at the phase crossover frequency ω p c H P . The correspondent ∣ H 12 ( j ω p c H P ) G p ( j ω p c H P ) ∣ = 0.20 leads to the margin G M H P = 14 dB.

The gain margin obtained with the controller transfer function H 12 ( j ω ) establishes a useful robustness measure, which shows the efficiency of the design approach. It guarantees a suitable distance between the intersection of the Nyquist plot of the compensated linear part with the negative real axis (square point in Figure 8) and the plot of the negative inverse of the dead zone describing function. Namely, Figure 8 indicates the design robustness to plant parametric variations. Part (a) considers a ± 20 % change in the plant gain, and part (b) considers a ± 20 % change in the lower time constant T 1 . These changes determine an uncertainty region whose limits are specified by the dashed red curves around the Nyquist plot of G p ( j ω ) . Then, the resulting region of the compensated system has limits given by the dashed black curves around the Nyquist plot of H 12 ( j ω ) G p ( j ω ) . The first square point on the nominal Nyquist plot (solid line) located between the limit curves intersects the unit circle, thus providing the gain crossover frequency and the phase margin. The second square point on the nominal plot corresponds to the phase crossover frequency and the gain margin. Then, as the figures show, the design is quite robust with respect to parametric variations.

Figure 8

Evaluation of robustness to plant parametric uncertainties. (a) Sensitivity to ± 20 % variation in the plant gain and (b) sensitivity to ± 20 % variation in the plant lower time constant.

3.2.1 Integer-order network

An integer-order lead network can be designed by a classical procedure [28,29]. The proposed fractional-order control design draws inspiration from this classical method. In the integer-order case, since the required phase lead ( 7 5 ° ) in the example is higher than the maximum achievable by a lead compensator (about 6 5 ° ), a double integer-order lead network, i.e., two equal lead networks with the same parameters are designed. Then, the controller is represented by H ( s ) = 1 + τ a s 1 + τ a α s 2 . The classical design procedure gives α = 0.24 and τ a = 0.36 s. Then, a phase margin P M = 4 5 ° is obtained at ω g c = 5.69 rad/s , and the gain margin is G M = 11.2 dB at ω p c = 13.1 rad/s . Although the frequency-domain specifications are met and the Nyquist plot does not intersect the plot of the negative inverse of the describing function, the time-domain control error is not satisfactory because the steady-state error is not zero, as shown in Figure 9, where a reference step input r = 4.5 is considered to generalize the case of r = 0 . We may conclude that the integer-order controller is not suitable, thus showing the superiority of the FOC.

Figure 9

Control error with noninteger-order and integer-order networks.