Analysis of multiple oscillations in a sampling reduced-order sliding mode control in frequency domain

1. Introduction

Sliding mode (SM) control is one of the most common-used robust control algorithms in control system and nonlinear mechanics due to its prominent robust performance [1-3] and low cost [4]. It is also the important nonlinear algorithm in the neural network and machine learning [5, 6]. However, the undesired phenomenon of oscillations in SM system known as “chattering” or “limit cycle” hinders its applications in engineering since the mechanical wear and resonances caused by oscillations. Hence, previous studies [7-9] paid much attention on suppression or elimination of the oscillations in SM controller. Although the high order SM controller, such as the second-order SM [10] and the adaptive PID SM controller [11], would offer oscillation elimination because of the finite-time convergence, the popular high order SM algorithm [12] still exhibits oscillations in the presence of unmodeled dynamics such as parasitic dynamics, delays and hysteresis. Hence, the oscillation-free SM algorithms merely exist in the ideal system without unmodeled dynamic [13]. Hence, it is worthwhile to provide a methodology to assess and adjust the parameters of those hard-to-remove oscillations for practicing control engineers.

Compared with the Lyapunov method and Poincare maps, the describing function (DF) method provides a relatively simple and efficient solution for the oscillation in the nonlinear control algorithm [14, 15]. For the SM algorithm, the DF method and Routh Criterion were employed to predict the stability of oscillations in the uncertain SM system [16]. Its robustness could be tuned by the intersection angle between the Nyquist plot of the plant and the negative reciprocal DF of controller. Besides, the influences of the delay, hysteresis and saturation on the chattering of SM controller were analyzed by the DF method and stability-equation method in frequency domain [17]. In addition, the oscillation in a second-order SM controller could be set by parameter tuning under the guidance of the DF method [12]. The traditional DF method is effective for these uncoupled nonlinearities in [12, 16, 17], which only depends on the amplitude of oscillations. However, the SM algorithm may depend on both the amplitude and frequency of oscillations with the increase of complexity of controllers, and it is hard for the traditional DF method to analyze the coupled nonlinearity in SM controllers.

On the other hand, the sliding surface is closely related to the stability of SM system and oscillation parameters. The most popular full-order sliding surface is a linear combination of all the state variables of the plant [18]. However, it is unpractical to obtain the accurate measurements of high-order state variables in engineering applications. In this case, the output feedback SM controller is a feasible choice, such as the twisting second-order SM [12, 19], when the relative degree is 2 or 3. However, it is challenging for the output feedback SM controller to provide sufficient robustness [20] because of low switching frequency of the signum function when the relative degree is ≥4. Therefore, the reduced-order sliding surface [21, 22] which is the linear combination of partial state variables is a convenient method to decrease the relative degree. Due to the independence of the sliding surface and the high order state variables, the fulfillment of the reachability and stability may become complicated.

Most of the researches assume that the limit cycles in the SM system are always unique and globally stable. Nevertheless, the nonlinear system usually exhibits the local stability. A nonlinear absorber system in [23] demonstrated a complex dynamic of the bi-stable oscillation. The system could approach to the high-frequency oscillation from the low-frequency one after the perturbation of an excitation force. Similarly, the simulation in a PD fuzzy controller [24] illustrated the attraction basin of a limit cycle. When the state portrait started in the attraction basin, it converged to the origin; when the state portrait started out of the attraction basin, it diverged to infinite. A similar phenomenon that two limit cycles existed in a switching power converter was analyzed in [25]. The DF method is an effective tool for analyzing the complex dynamic behavior in the SM system. The existence of multiple limit cycles in a piecewise gain function was proved based on the DF method [26]. The system could switch from a small limit cycle to a large one when a sufficiently large step signal applied to the input. A generalized DF approach was proposed to estimate multiple limit cycles in a SM system [27, 28], and the relation among these limit cycles were revealed intuitively. Under the influence of the sampling interval, multiple limit cycles were found to coexist in a reduced-order SM system [22]. However, the stability of multiple limit cycles in [22, 26-28] was analyzed by using the graphical DF approach which is only fit for the minimum phase system. Hence, the analysis of the stability of limit cycles needs further perfection.

The main contribution of this paper is the prediction of multiple oscillations in the sampling reduced-order SM system by an improved DF method. The characteristic equation of the nonlinear SM system is established and its pole locations are used as the criterion to determine the stability of oscillations. The criterion is valid even the poles located in the right hand $s$-plane (RHP), hence, the improved DF method is applicable to the minimum and non-minimum phase system. In addition, the parameter tuning method is also introduced to eliminate the undesired oscillation guided by the existence condition of oscillations. In this way, a locally stable system can be tuned as a globally stable one, and the performance of the system can be improved by sufficient gain. The simulations show the existence of three limit cycles which confirms the predictions, and demonstrate the benefit of the parameter tuning method.

The paper is organized as follows. Firstly, the reduced-order sampling SM system is considered. Then the frequency-domain model of the controller is calculated and the improved DF method is derived. In the following section, the complex dynamics including multiple oscillations and their start-oscillation conditions are estimated and simulated. Finally, the coefficient tuning method for eliminating undesired oscillations is presented.

2. Design of the control system

Consider a plant given by the following equation:

where $x\in R^n$, $u\in R^1$, $A\in R^{n\times n}$, $B\in R^{n\times 1}$; $\sigma$ is the SM surface and $C\in R^{1\times n}$ is the vector of the sliding surface coefficients, and the control $u$ generated by the SM controller is defined as:

where gain $M$ is a positive constant; $k$ is a sampling counter; $T_e$ is a sampling period.

Assume that the high-order state-variables are immeasurable, which means that the coefficients $c_i$ corresponding to the immeasurable state-variables equal to zero, then a reduced-order sliding surface could be designed to reduce the complexity of the system. However, the controller may lead to complex dynamic and stability. The following sections will focus on the prediction and adjustment of the multiple oscillations in the reduced-order SM controller.

3. Model of the SM control system in frequency domain

The system can be divided into linear and nonlinear elements. The output of the linear element $\sigma \left(t\right)$ is the input to the nonlinear element, and the output of the nonlinear element $u\left(t\right)$ is also the input to the linear element. The linear element is the composition of the plant$\dot{x}=Ax+Bu$ and the sliding surface $\sigma =Cx\text{.}$ Hence, the transfer function of the linear element is $W\left(s\right)=C{\left(sI-A\right)}^{-1}B\text{,}$ which can also be written as the ration of two polynomials $W\left(s\right)=P\left(s\right)/Q\left(s\right)$.

To model the nonlinear element $u\left(t\right)=f\left[-\sigma \left(t\right)\right]$, its DF should be calculated. Firstly, suppose that a sinusoidal signal $-\sigma \left(t\right)=A\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega t\right)$ ($A>$0) is applied to the input of $M$ sign $\left[-\sigma \left(k{T}_e\right)\right]$, where $A$ and $\omega$ denote the amplitude and angular frequency (For simplicity, $\omega$ is called the frequency in this paper), respectively. And then, the waveforms of $-\sigma \left(t\right)$ and $u\left(t\right)$ are plotted in Fig. 1. The signal $-\sigma \left(t\right)$ is the excitation of $M$ sign $\left[-\sigma \left(k{T}_e\right)\right]$ and $u\left(t\right)$ is its response. Finally, the nonlinear element $f\left[-\sigma \left(t\right)\right]$ can be approximated as a describing function as shown in Eq. (2):

where $a_1\left(\omega \right)$ and $b_1\left(\omega \right)$ are fundamental harmonic Fourier Series which can be obtained by the integral in a period of $\sigma \left(t\right)$:

According to the waveforms in Fig. 1, $u\left(t\right)$ lags behind $-\sigma \left(t\right)$ with uncertain interval $\mathrm{\Delta }{T}_1\in \left[0,T_e\right]$ and $\mathrm{\Delta }{T}_2\in \left[0,T_e\right]$. For simplicity, let $\mathrm{\Delta }{T}_1=\mathrm{\Delta }{T}_2=T_e$ and thus the fundamental harmonic Fourier Series can be written as:

Compared with the uncoupled nonlinearities $N\left(A\right)$ which only depends on the amplitude, Eq. (2) is a more complex coupled nonlinearity $N\left(A,\omega \right)$ which is closely related to the frequency and amplitude. Hence, it is difficult for the traditional DF method to estimate the periodic motion in the coupled nonlinearity. For a simpler analysis procedure, the following section will discuss an improved DF method.

Fig. 1Waveforms of the sinusoidal excitation and response

4. Improved DF method

To estimate the oscillation parameters, the system is presented as a basic feedback structure which consists of the nonlinear element $N\left(A,\omega \right)$ described by Eq. (2) and the linear element as follows:

Assume that the sliding surface $\sigma \left(t\right)$ moves on the oscillation with frequency ${\omega }_n$ and amplitude $A_n$; hence, $\sigma \left(t\right)=A_n\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_n\right)t$, and the negative of $\sigma \left(t\right)$ can be written as $-\sigma \left(t\right)=A_n\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_{n}t-\pi \right)$. After pass through the sampling term and switching function, the signal $u\left(t\right)$ is changed in both the magnitude and phase as follows:

where $\left|N\left(A_n,{\omega }_n\right)\right|$ and ${\theta }_N\left({\omega }_n\right)$ are the magnitude and phase angle of $N\left(A_n,{\omega }_n\right)$, respectively.

When $u\left(t\right)$ goes through the linear plant $W\left(s\right)$, it becomes:

where $\left|W\left(j{\omega }_n\right)\right|$ and ${\theta }_W\left({\omega }_n\right)$ are the magnitude and phase angle of $W\left(s\right)$, respectively, at the frequency ${\omega }_n$. Substitute:

into Eq. (6a), $\sigma \left(t\right)$ can be rewritten as:

Hence, the amplitude $A_n$ and frequency ${\omega }_n$ can be estimated according the following equations:

However, Eq. (7a) cannot be solved straightforward for a high-order plant. Therefore, the interative operations or function plot can be used to estimate ${\omega }_n$. It proposes that the phase-frequency characteristic $\phi \left(\omega \right)={\theta }_W\left(\omega \right)-\omega T_e$ can be plotted in the frequency $0\le \omega \le \pi /T_e$ to find the ${\omega }_n$ such that $\phi \left({\omega }_n\right)=\left(-2k+1\right)\pi$, ($k=$1,…, $n$)$.$

The analysis above is helpful for to determine the possible oscillations. However, the stability of the oscillations needs to be further determined. To analyze the stability of the system, the transfer function of the closed-loop system can be written as:

The characteristic equation $D\left(s\right)=1+W\left(s\right)N\left(A,\omega \right)=Q\left(s\right)+P\left(s\right)N\left(A,\omega \right)$ determines its stability. Replacing the term $j{A}_1$ by $\left(A_1/\omega \right)s$, then $D\left(s\right)$ can be represented as:

where $a_n^{\text{'}}$ is the function of $A$ and $\omega$. According to the small amplitude perturbation analysis method, we get the following stability criterion.

Criterion I: The closed-loop system has a stable oscillation with parameters $A_n$ and ${\omega }_n$ if there exists an arbitrarily small positive number $∆A$ such that the characteristic equation $D\left(s\right)=Q\left(s\right)+P\left(s\right)N\left(A_n+∆A,{\omega }_n\right)$ is stable and $D\left(s\right)=Q\left(s\right)+P\left(s\right)N\left(A_n-∆A,{\omega }_n\right)$ is unstable; otherwise the oscillation is unstable.

Based on the criterion, the stability of an oscillation can be predicted via the following two steps.

Step 1: If there exists an arbitrarily small positive number $∆A$ such that all the roots of $D\left(s\right)=Q\left(s\right)+P\left(s\right)N\left(A_n+∆A,{\omega }_n\right)$ are in the left hand $s$-plane (LHP) and any of the roots of $D\left(s\right)=Q\left(s\right)+P\left(s\right)N\left(A_n-∆A,{\omega }_n\right)$ is in the RHP, then the oscillation $\left(A_n,{\omega }_n\right)$ is stable.

Step 2: If there exists an arbitrarily small positive number $∆A$ such that any of the roots of $D\left(s\right)=Q\left(s\right)+P\left(s\right)N\left(A_n+∆A,{\omega }_n\right)$ are in the RHP or all the roots of $D\left(s\right)=Q\left(s\right)+P\left(s\right)N\left(A_n-∆A,{\omega }_n\right)$ is in the LHP, then the oscillation $\left(A_n,{\omega }_n\right)$ is unstable.

The improved DF method could intuitively predict multiple oscillations by the curve $\phi \left(\omega \right)$. It also gives us a guidance to adjust the oscillation amplitude and frequency by analyzing the variation of $\phi \left(\omega \right)$ when coefficients are tuned. Compared with the method in [27, 28], the improved DF method is valid for the minimum and non-minimum phase system. In addition, the roots of $D\left(s\right)$ can be calculated by the computer software.

5. Analysis of multiple oscillations

The main task of this section is to evaluate the accuracy of the improved DF method and to demonstrate the behavior of multiple oscillations. Firstly, consider a fourth-order canonical form plant:

where $a_1=$ 2, $a_2=$ 5, $a_3=$ 12, $a_4=$ 7. Assume that high-order state-variable $x_4$ is immeasurable, hence, the sliding surface coefficient is designed as a reduced-order one $C=\left[c_1,c_2,c_3,c_4\right]=$[4, 1, 0.5, 0], and the controller is set as $T_e=\text{10}$ ms and $M=\text{1}$. Hence, the transfer function of the linear element can be written as:

According to Eq. (2) and (3), the describing function of the SM controller can be written as:

Next, Eqs. (11) and (12) will be substituted into the improved DF method to estimate the periodic motions.

5.1. Prediction of multiple oscillations

In accordance with Eq. (7a), the oscillation frequency ${\omega }_n$ satisfies the following equation:

$\mathrm{t}{\mathrm{g}}^{-1}\frac{{\omega }_n}{4-0.5{{\omega }_n}^2}+\mathrm{t}{\mathrm{g}}^{-1}\frac{7{{\omega }_n}^3-5{\omega }_n}{{{\omega }_n}^4-12{{\omega }_n}^2+2}-0.01{\omega }_n=\left(-2k+1\right)\pi .$

By plotting the phase-frequency characteristic:

13

$\phi \left(\omega \right)=\mathrm{t}{\mathrm{g}}^{-1}\frac{\omega }{4-0.5{\omega }^2}+\mathrm{t}{\mathrm{g}}^{-1}\frac{7{\omega }^3-5\omega }{{\omega }^4-12{\omega }^2+2}-0.01\omega ,$

As shown in Fig. 2, we can get ${\omega }_n$ directly according to the condition $\phi \left({\omega }_n\right)=\left(-2k+1\right)\pi$. In Fig. 2, three oscillation frequencies are estimated as ${\omega }_{n1}=\text{1.06}$ rad/s, ${\omega }_{n2}=\text{2.54}$ rad/s and ${\omega }_{n3}=\text{22.12}$ rad/s. Hence, three oscillations may exist in the system. By substituting the frequencies ${\omega }_{n1}$, ${\omega }_{n2}$ and ${\omega }_{n3}$ into Eq. (7b), the corresponding amplitudes are calculated as $A_{n1}=$ 4.30×10^-1, $A_{n2}=$ 3.14×10^-2 and $A_{n3}=$ 1.25×10^-3. Finally, the stability of the three oscillations needs to be determined. The characteristic equation $D\left(s\right)$ can be written as:

14

$D\left(s\right)=Q\left(s\right)+P\left(s\right)N(A_n,{\omega }_n)$
$=s^4+a_4{s}^3+a_3{s}^2+a_{2}s+a_1+\frac{4M}{A_n\pi }\left(\mathrm{c}\mathrm{o}\mathrm{s}\right({\omega }_n{T}_e)-j\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_n{T}_e\left)\right)(c_3{s}^2+c_{2}s+c_1)$
$=a_4^{\text{'}}{s}^4+a_3^{\text{'}}{s}^3+a_2^{\text{'}}{s}^2+a_1^{\text{'}}s+a_0^{\text{'}}.$

where:

$a_4^{\text{'}}=1,a_3^{\text{'}}=a_4-\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_n{T}_e\right)\frac{c_3}{{\omega }_n}\frac{4M}{\pi A_n},$
$a_2^{\text{'}}=a_3-\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_n{T}_e\right)\frac{c_2}{{\omega }_n}\frac{4M}{\pi A_n}+\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_n{T}_e\right)c_3\frac{4M}{\pi A_n},$
$a_1^{\text{'}}=a_2-\mathrm{s}\mathrm{i}\mathrm{n}\left({\omega }_n{T}_e\right)\frac{c_1}{{\omega }_n}\frac{4M}{\pi A_n}+\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_n{T}_e\right)c_2\frac{4M}{\pi A_n},a_0^{\text{'}}=a_1+\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_n{T}_e\right)c_1\frac{4M}{\pi A_n}.$

Fig. 2Diagram of φω with three oscillations

Table 1Variation of the roots of polynomial Ds with the small perturbation on oscillation amplitudes

Oscillation	Characteristic equations $D\left(s\right)$	Root 1	Root 2	Root 3	Root 4
$\omega ={\omega }_{n1}=$ 1.06 rad/s $A=A_{n1}+∆A=$0.4343	$D\left(s\right)=s^4+6.9853s^3$ $+13.4365s^2+7.8143s$ $+13.7262$	–3.4916 + j0.3393	–3.4916 – j0.3393	–0.0010 + j1.0561	–0.0010 – j1.0561
$\omega ={\omega }_{n1}=$ 1.06 rad/s $A=A_{n1}-\mathrm{\Delta }A=$ 0.4257	$D\left(s\right)=s^4+6.9850s^3$ $+13.4655s^2+7.8711s$ $+13.9631$	–3.4937 + j0.3811	–3.4937 – j0.3811	0.0011 + j1.0633	0.0011 – j1.0633
$\omega ={\omega }_{n2}=$ 2.54 rad/s $A=A_{n2}+\mathrm{\Delta }A=$ 0.0317	$D\left(s\right)=s^4+6.7993s^3$ $+31.6659s^2+43.5289s$ $+162.5384$	–3.4038 + j3.7029	–3.4038 - j3.7029	0.0042 + j2.5348	0.0042 – j2.5348
$\omega ={\omega }_{n2}=$ 2.54 rad/s $A=A_{n2}-\mathrm{\Delta }A=$ 0.0311	$D\left(s\right)=s^4+6.7952s^3$ $+32.0632s^2+44.3072s$ $+165.7816$	–3.3937 + j3.7409	–3.3937 – j3.7409	–0.0039 + j2.5492	–0.0039 – j2.5492
$\omega ={\omega }_{n3}=$22.12 rad/s $A=A_{n3}+\mathrm{\Delta }A=$ 0.0013	$D\left(s\right)=s^4+2s^3+494s^2$ $+948.9s+3937.7$	–0.0229 + j22.0401	–0.0229 – j22.0401	–0.9763 + j2.6745	–0.9763 – j2.6745
$\omega ={\omega }_{n3}=$22.12 rad/s $A=A_{n3}-\mathrm{\Delta }A=$ 0.0012	$D\left(s\right)=s^4+1.9s^3$ $+503.7s^2+968s$ $+4017.2$	0.0281 + j22.2644	0.0281 – j22.2644	–0.9768 + j2.6739	–0.9768 – j2.6739

For simplicity, $∆A$ is considered as a hundredth of $A_n$. By using the numerical software, the roots of six characteristic equations are listed in Table 1 in accordance with Step 1 and Step 2. Table 1 shows that oscillation 1 $\left(A_{n1},{\omega }_{n1}\right)=$ (4.30×10^-1, 1.06 rad/s) is stable because four roots are located in LHP when $A=A_{n1}+∆A$, and root 3 and root 4 are located in RHP when $A=A_{n1}-∆A$. It also shows that oscillation 2 $\left(A_{n2},{\omega }_{n2}\right)=$(3.14×10^-2, 2.54 rad/s) is unstable because root 3 and root 4 are located in RHP when $A=A_{n2}+∆A$. Similarly, oscillation 3 $\left(A_{n3},{\omega }_{n3}\right)=$ (1.25×10^-3, 22.12 rad/s) is stable. The improved DF method predicts that the system has two stable oscillations and one unstable oscillation.

5.2. Simulation of multiple oscillations

Simulation I is carried out to verify the three oscillations and their stability by using the numerical software. The state space in Eq. (10) is discretized by $G=I+A\times T_e$ and $H=B\times T_e$. The controller updates its output $u\left(t\right)$ once in every $T_e$. The limit cycles and frequency spectrums can be plotted by using the data $\sigma \left(t\right)$ in the steady-state process. The time series, phase portraits and frequency spectrums of $\sigma$ are plotted in Fig. 3(a), (b) and (c), in which the oscillation amplitudes and frequencies can be measured directly.

Fig. 3a) Three types of motions in time domain, b) the phase portraits of the steady state of three oscillations, c) spectrum analysis of the steady state of three oscillations, d) stability of oscillation 2

a)

b)

c)

d)

Firstly, three sets of initial conditions $\left(x_1,x_2,x_3,x_4\right)=$(0.1, 0.008, 0,0), $\left(x_1,x_2,x_3,x_4\right)=$(0.01, 0.008, 0,0), $\left(x_1,x_2,x_3,x_4\right)=$(0.01, 0.006, 0,0) are selected to generate the three kinds of distinctly different motions in time domain (as shown in Fig. 3(a) and in phase plane $\sigma -\dot{\sigma }$ (as shown in Fig. 3(b). After the transient processes, the three motions approach their respective periodic oscillations. The measured amplitudes of the three oscillations are $A_{n1}=$ 4.12×10^-1, $A_{n2}=$ 3.18×10^-2) and $A_{n3}=$ 1.32×10^-3. The frequency spectrum analysis is also made as shown in Fig. 3(c), by which we obtain spectrum components for each oscillation and measure the frequencies of the three oscillations: ${\omega }_{n1}=$1.09 rad/s, ${\omega }_{n2}=$2.49 rad/s and ${\omega }_{n3}=$ 21.09 rad/s. Figs. 3(a), (b) and (c) demonstrate that the system exists three oscillations and it convergences to which oscillation depends on the initial conditions.

Next, the stability of the three oscillations needs to be verified. Obviously, oscillations 1 and 3 are locally stable. But the stability of oscillation 2 is not clear according to the above simulations. Fig. 3(d) shows the stability of oscillation 2 when disturbance is applied to the controller. Affected by external disturbances, the system is driven to escape oscillation 2 and approaches to oscillation 1 or 3.

Table 2Comparison of parameters of multiple oscillations estimated by DF method with the results of simulation I

Plant and sliding mode surface coefficients	Oscillations	Stability	Amplitude $A_n^{}$		Frequency ${\omega }_n$ (rad/s)
			Theory predict	Simulation I result	Theory predict	Simulation I result
$[a_1,a_2,a_3,a_4]=$ [2, 5, 12, 7] $\left[c_1,c_2,c_3,c_4\right]=$ [4, 1, 0.5, 0]	1	Locally stable	4.30×10-1	4.12×10-1	1.06	1.09
	2	Unstable	3.14×10-2	3.18×10-2	2.54	2.49
	3	Locally stable	1.25×10-3	1.32×10-3	22.12	21.09
$M=$ 1; $T_e=$ 10 ms

Hence, it shows that oscillation 2 is unstable. Finally, the estimated and simulated parameters of the three oscillations are listed in Table 2, and the theoretical results and the simulations are in good agreement. However, the errors between the theory predictions and simulations may mostly originate from the fundamental harmonic approximation of the DF method. In fact, the third and fifth harmonic which is approximately assumed to have been completely suppressed by the low-pass performance of $W\left(s\right)$ may affect the phase angle of the closed loop. Hence, the accuracy of predictions depends on the intensity of residual high-order harmonics in $\sigma \left(t\right)$.

Fig. 4a) Attraction zones for each oscillation in x1-x2 plane, b) the relation between the attraction zones and the oscillation 1 and 3

a)

b)

The simulation indicates that the system converges to which oscillation depends on the initial conditions located in which region. In another word, each stable limit cycle has its own attraction region. However, predicting the attraction basin is difficult. The attraction region in the fuzzy system was plotted using simulations [24]. For simplicity, we only consider the influence of two state-variables $x_1$ and $x_2$ on the dynamic of the system. Equally spaced initial conditions are picked in the region $R:\left\{-0.02\le x_1\le 0.02,-0.04\le x_2\le 0.04,x_3=0,x_4=0\right\}$, and the steady-state processes corresponding to each point are simulated. We classified the processes as oscillation 1 or 2 or 3 based on their spectrum components. Subsequently, the region R can be divided into three attraction zones as shown in Fig. 4(a). The domain outside of the blue curve is the attraction zone of oscillation 1, and the domain in the blue curve is the attraction zone of oscillation 3. In addition, the attraction zone of oscillation 2 is very small which is represented as the only two red points.

The phenomenon of multiple oscillations can be interpreted as the inconsistent stability of the system, which is the significant difference compared with the linear system. The equivalent gain of the switching function $M$ sign $\left[-\sigma \left(t\right)\right]$ can be regarded as $k_n=2M/\pi A_n\left[20\right]$. $k_n$ changes along with the oscillation amplitude. Suppose that the system starts from point A as shown in Fig. 4(b), then the system is unstable and its trajectory is called the unstable spiral. In the motion, the gain $k_n$ continues to reduce when $A$ is less than $A_{n1}$ until $A$ equals $A_{n1}.$ Similarly, when the system starts from point B as shown in Fig. 4(b), then the gain $k_n$ continues to increase when $A>A_{n3}$ until it converges to oscillation 3. Hence, the blue curve in Fig. 4(b) is similar with the watershed.

6. Elimination of undesired oscillation

According to the analysis above, the system may switch between oscillations 1 and 3 with the variation of external conditions. We hope to eliminate oscillation 1 due to its small equivalent gain which may lead to poor precision and weak robustness. Furthermore, the system could become the globally stable one if oscillation 1 can be eliminated.

6.1. Coefficient tuning method

The fact that $\phi \left(\omega \right)$ crosses the line $\phi =-\pi$ three times leads to multiple oscillations. $\phi \left(\omega \right)$ can be rewritten as:

15

$\phi \left(\omega \right)=\angle P\left(j\omega \right)+\angle \frac{1}{Q\left(j\omega \right)}+\angle N\left(A,\omega \right),$

where:

$\angle P\left(j\omega \right)=\mathrm{t}{\mathrm{g}}^{-1}\left(\frac{c_2\omega -c_4{\omega }^3}{c_1-c_3{\omega }^2}\right),\angle \frac{1}{Q\left(j\omega \right)}=\mathrm{t}{\mathrm{g}}^{-1}\frac{a_4{\omega }^3-a_2\omega }{{\omega }^4-a_3{\omega }^2+a_1},\angle N\left(A,\omega \right)=\omega T_e.$

Assuming that the plant (as shown in Eq. (10) and sampling period are invariant, tuning the sliding surface is a feasible method to eliminate the undesired oscillation. In the following, we will discuss the influence of the sliding surface on $\phi \left(\omega \right)$, and eliminate oscillation 1 by tuning the coefficients $C=\left[c_1,c_2,c_3,c_4\right]$. As shown in Fig. 2, if the minimum value of $\phi \left(\omega \right)$ in the low-frequency band can be increased to more than $–\pi$, then there is only one intersection corresponding to oscillation 3. Increasing $c_3$ or $c_2$, or decreasing $c_1$ or $c_4$ will make $\phi \left(\omega \right)$ rise faster in the low-frequency band. Hence, such tuning method could achieve the purpose. For simplicity, we only consider $c_3$ as an adjustable coefficient and also set $c_1=$ 4, $c_2=$ 1, $c_4=$ 0, $a_1=$ 2, $a_2=$ 5, $a_3=$ 12, $a_4=$ 7 and $T_e=$ 10 ms.

As shown in Fig. 5, the minimum value of $\phi \left(\omega \right)$ in the low-frequency band increases with the variation of $c_3$ from 0.5 to 2. When $c_3=$ 0.5 and $c_3=$ 0.8, there exists three intersections between $\phi \left(\omega \right)$ and the straight line $\phi =–\pi$. When $c_3\ge$1.1, there remains only one intersection in the high-frequency band, and the system becomes a globally stable one. Table 3 also shows the relation among the parameters of oscillations, coefficient $c_3$ and zero placement. It is found that oscillation 1 is eliminated and there is little variation in the frequency of oscillation 3 when $c_3$ increases.

From the perspective of the phase-frequency characteristic of the sliding surface, we find the relations between $\angle P\left(j\omega \right)$ and the number of oscillations. According to the zero placement, the damping coefficient is increased and the natural frequency is reduced with the increase of $c_3$. Hence, $\angle P\left(j\omega \right)$ could provide greater increase rate for $\phi \left(\omega \right)$ at the lower frequency band with the increase of $c_3$. Therefore, $\phi \left(\omega \right)$ can be uplifted in the low-frequency band. The undesired oscillation 1 can be eliminated when the minimum value of $\phi \left(\omega \right)$ is greater than $–\pi$.

Fig. 5Variation of curves φω with c3

Table 3The relation among the coefficient c3, parameters of oscillations and zero placement

$c_3$	${\omega }_{n1}$(rad/s)	$A_{n1}$	${\omega }_{n3}$(rad/s)	$A_{n3}$	zero placement
0.5	1.06	4.30×10-1	22.12	1.25×10-3	–1.0000 ± j2.6458
0.8	1.11	3.45×10-1	23.80	1.70×10-3	–0.6250 ± j2.1469
1.1	–	–	24.52	2.30×10-3	–0.4545 ± j1.8520
1.4	–	–	24.88	2.80×10-3	–0.3571 ± j1.6521
1.7	–	–	25.20	3.30×10-3	–0.2941 ± j1.5055
2.0	–	–	25.36	3.90×10-3	–0.2500 ± j1.3919

6.2. Verification for tuning method

Simulation II is carried out to verify the stability of the tuned system when the initial condition is set as $\left(x_1,x_2,x_3,x_4\right)=$(0.1, 0.7, 0, 0) and $c_3$ is set as a number of different values. Fig. 6 shows the system responses in time series when $c_3$ equals to 0.5, 0.8, 1.1, 1.4, 1.7, 2.0, 3.0 and 5.0. When $c_3=$ 0.5 and $c_3=$ 0.8, the system convergences to oscillation 1 which is a large and slow periodic motion. When $c_3\ge$ 1.1, the system convergences to oscillation 3 since oscillation 1 has been eliminated, and oscillation 3 becomes a globally stable oscillation.

Fig. 6a) System responses when c3 equals 0.5, 0.8, 1.1 and 1.4, b) system responses when c3 equals 1.7, 2.0, 3.0 and 5.0

a)

b)

In addition, the comparison of these responses in Fig. 6 shows that sufficient large $c_3$ can also improve system performances in overshoot and settling time. The results of the two performances are given in Table 4.

Table 4Overshoot and settling time

$c_3$	1.1	1.4	1.7	2.0	3.0	5.0
Overshoot	3.935	3.830	3.725	3.622	3.284	2.649
Settling time (s)	42.1	29.73	24.69	22.01	18.91	16.58

7. Conclusions

This paper has presented a frequency-domain approach to analyze the stability, oscillation frequency and amplitude, and attraction zone of the reduced-order SM system by using an improved DF method. The method allows for the predictions of oscillation parameters in the coupled nonlinearities. The stability of oscillations can be determined by the pole placement; hence, the method is suitable for the minimum and non-minimum phase system. The analysis focuses on the multiple oscillations in a fourth-order system in which there exists one unstable and two locally stable oscillations. More specifically, the reduced-order SM controller and zero-order holder are the major cause of multiple periodic motions. Here, it should be pointed out that the system converges to which oscillations depends on the location of its initial condition, and a method to eliminate the undesirable oscillation by tuning coefficients is presented. The simulations in time domain support the prediction of multiple oscillations and their parameters. Although the inaccuracy of the estimation of oscillations cannot be neglected due to the inherent approximation of the DF approach, this study provides guidance for the design of the SM system.