Influence of Control Valve Delay and Dead Zone on the Stability of a Simple Hydraulic Positioning System

Abstract

This paper deals with the PI control of a highly simplified dynamic model of a hydraulic cylinder. It is assumed that the hydraulic fluid is incompressible and that the pump provides constant flow rates, which results in the possibility of velocity control. Two types of anomalies are taken into account: (a) the time delay due to the controller computations and the internal pressure dynamics and (b) the dead zone of the controller valve. This results in a nonlinear system are described by a piecewise linear discontinuous map. Nonlinear behavior of the system is explored and the practically globally stable parameter domains are identified.

1. Introduction

Hydraulic systems are widely used in heavy-duty industrial applications, where the exertion of high forces with large stiffness is needed in a robust way. Although there is a considerable effort on developing advanced control strategies (see e.g., [1–4]), PID control still remains the most popular choice. However, it is well known that strong nonlinearities are present in these systems, such as pressure-flow rate relationship, dead zone of the control valves (see e.g., [1]), dry friction [3] or impact dynamics [5]. The discrete sampling time of the closed-loop control introduces additional complexity together with the response lag due to internal (mostly pressure) dynamics. Thus designing and tuning a PID controller of a hydraulic system is a highly challenging task mostly because the conventional ways are based on linear system theory. Moreover, some of the above-mentioned nonlinearities (e.g., dead zone or impact dynamics) cannot be coped with using linearization techniques.

The mathematical modelling of these systems often leads to equations with nonsmooth or even discontinuous right-hand side. Fortunately the progress in the theory of non-smooth dynamical systems (see e.g., [6] for an overview) provides a toolbox, albeit it is still far from being general. This is especially true for systems of higher dimensions (three, four, etc.) with several regions of different dynamics. For example, [7] gives a general theory on the existence of periodic and dense orbits for a bilinear one-dimensional map with a slight extension towards two-dimensional maps with delay and backlash. In [8], the authors study the border collision bifurcation in 𝑛-dimensional maps with two regions. Chaotic oscillations are also identified in these systems [9]. Based on numerical simulations, [10] describes an example on the effect of delay and backlash together.

This paper studies a highly simplified model of a hydraulic positioning system, which, despite its simplicity (from the engineering point of view) and linearity, poses interesting mathematical problems.

2. Mathematical Model

The subject of our investigation is a digitally controlled hydraulic system that consists of a differential hydraulic cylinder, a proportional directional valve, a linear potentiometer as position transducer, a gear pump, and a PC. The PC provides the PI (proportional-integral) controller. It receives the signal from the position transducer, calculates the error signal, and drives the hydraulic valve; see Figure 1. A typical characteristic of a directional proportional valve is shown in Figure 2.

Figure 1 

The hydraulic positioning system. 1: hydraulic cylinder, 2: position transducer, 3: PC on which the PI controller is implemented, 4: proportional directional valve, 5: electromotor, 6: gear pump, 7: pressure limiting valve.

Figure 2 

Dead zone in a typical flow rate/control signal characteristics of a directional proportional hydraulic valve. In the interval $(-u_{\min \hspace{0.17em}},u_{\min \hspace{0.17em}})$ there is no fluid flow.

The mass of the piston is neglected in this study, and the Newtonian dynamics of the system is further simplified by not considering frictional forces at the sealing of the piston rod. Clearly, the latter one has an essential influence on the nonlinear dynamics of the system due to small positive or even negative damping values. When we carry out the investigation with zero damping, we analyze the critical case which already presents an intricate dynamics due to the modeled delay, dead zone, and sampling.

The continuous physical process is sampled in time intervals 𝑡_s (sampling time), thus the position of the piston rod 𝑥(𝑡) is discretized in time as 𝑥_𝑛=𝑥(𝑛𝑡_𝑠). This position is fed into the PC which computes the error signal ℎ. The time needed for this computation is denoted by 𝑡_c. Due to the internal dynamics of the hydraulic system-(notably pressure dynamics), the variation of the velocity of the piston rod follows a second (or even higher) order lag system that is approximated by another delay denoted by 𝑡_id, with subscript referring to internal dynamics. Thus, the overall delay between the previous sampling instant and its effect is 𝑡_𝑑=𝑡_𝑐+𝑡_id. Note that the actual values and ratio of 𝑡_𝑐 and 𝑡_id are irrelevant. Since the Newtonian dynamics is neglected, the velocity of the piston is piecewise constant, and it can be discretized in the following way: 𝑣_𝑛(𝑡)≡𝑣((𝑛−1)𝑡_𝑠+𝑡_𝑑), 𝑡∈[(𝑛−1)𝑡_𝑠+𝑡_𝑑,𝑛𝑡_𝑠+𝑡_𝑑), accordingly, subscripts 𝑛 refer to different time instants for positions and velocities as it is represented in the controlling scheme of Figure 3.

Figure 3 

Scheme of the control. Instants signed with green circles represent the sampling; red circles show the moments of velocity actuations. For $h_n$, see equation (2.2); $t_c$ stands for computational time, $t_{\text{id}}$ represents the internal dynamics of the hydraulic system (see text for details), $t_d=t_c+t_{\text{id}}$ is the overall delay, and $t_s$ denotes sampling time.

Assume the cases when we have 0≤𝑡_𝑑≤𝑡_𝑠. Integrating the piecewise constant velocities, we arrive at a relationship between two neighbouring sampled piston rod positions, which can be expressed as 𝑥_𝑛+1=𝑥_𝑛+𝑣_𝑛𝑡_𝑑+𝑣_𝑛+1(𝑡_𝑠−𝑡_𝑑). (2.1) Although this expression is similar to the Euler-discretization of the governing equations, this discrete form is the exact solution of the real physical system controlled digitally. Considering the proportional-integral controller, the error signal ℎ_𝑛 is calculated in the following form: ℎ_𝑛=𝑃𝑥_𝑛+𝐼𝑦_𝑛, (2.2) where 𝑃 is the proportional gain 𝐼 is the integral one, and 𝑦_𝑛=𝑦_𝑛−1+𝑡_s𝑥_𝑛−1 (2.3) is the discrete integral of the position function. The piston rod velocities 𝑣_𝑛+1 and 𝑣_𝑛 are calculated from ℎ_𝑛 and ℎ_𝑛−1, respectively, according to the simplified characteristics of the proportional directional valve. Figure 4 represents this reduced, saturation-free characteristics. The interval of closure is (−Δ,Δ) and the slope is characterized by −𝛼.

Figure 4 

Simplified, saturation-free valve characteristics. In the interval $(-\Delta ,\Delta )$ there is no fluid flow.

We introduce the dimensionless variables by means of 𝑥𝑥=𝐿,𝑣𝑣=𝐿,𝑦𝑦=𝐿,𝛿=𝛼Δ𝐿,ℎ_𝑛=𝛼ℎ_𝑛𝐿,𝑃=𝛼𝑃,𝐼̂=𝛼𝐼 (2.4) and by abuse of the notation we drop the hats immediately. According to Figure 4, the velocity is a piecewise function of the error signal. With dimensionless quantities, 𝑣_𝑛=𝛿−ℎ_𝑛−1,0−𝛿−ℎ_𝑛−1,ℎ_𝑛−1>𝛿,if|ℎ_𝑛−1ℎ|≤𝛿,_𝑛−1<−𝛿. (2.5)

In the subsequent sections, we are going to construct a 4-dimensional linear mapping for the backlash-free system. In the presence of backlash, we derive a piecewise linear mapping which is compiled from 9 linear maps of dimension 4. If one also investigates the case 𝑡_𝑑>𝑡_𝑠, then similar linear and piecewise linear mappings can be constructed, but their dimensions increase extremely. In order to represent the method in a compact mathematical form, we restrict the description to the basic case 0≤𝑡_𝑑≤𝑡_𝑠.

3. Stability Analysis of the Linear System

Eliminating the valve dead zone (𝛿=0), a linear valve characteristic means regarding to the piston rod velocity: 𝑣_𝑛=−ℎ_𝑛−1, (3.1) from equation (2.1): 𝑥_𝑛+1=𝑥_𝑛−ℎ_𝑛−1𝑡_𝑑−ℎ_𝑛(𝑡_𝑠−𝑡_𝑑). (3.2) Introducing 𝐳_𝑛, that consists of the actual and the previous piston positions, and integral values: 𝐳_𝑛=𝑥_𝑛𝑥_𝑛−1𝑦_𝑛𝑦_𝑛−1^𝑇, (3.3) with matrix formalism: 𝐳_𝑛+1=𝐀𝐳_𝑛, (3.4)𝑥_𝑛+1𝑥_𝑛𝑦_𝑛+1𝑦_𝑛=1−𝑃(𝑡_𝑠−𝑡_𝑑)−𝑃𝑡_𝑑−𝐼(𝑡_𝑠−𝑡_𝑑)−𝐼𝑡_𝑑𝑡1000_𝑠0100010_𝐀𝑥_𝑛𝑥_𝑛−1𝑦_𝑛𝑦_𝑛−1. (3.5) The stability of the system depends on the eigenvalues of matrix 𝐀, all of the absolute values of the eigenvalues have to be less than 1: det(−𝜇𝐈+𝐀)=0,|𝜇_i|<1,𝑖=1,…,4. (3.6) The characteristic polynomial of matrix 𝐀 is 𝜇⁴+(𝑃(𝑡_s−𝑡_𝑑)−2)𝜇³+(1+2𝑃𝑡_𝑑−(𝑃+𝐼𝑡_𝑑)𝑡_s+𝐼𝑡²_s)𝜇²+𝑡_𝑑(𝐼𝑡_s−𝑃)𝜇=0. (3.7) The polynomial has one root that equals to zero; therefore; it can be divided by 𝜇≠𝜇₄=0. Since the stability criteria of polynomials are determining the coefficients of polynomial so that all of the roots should be on the left side of the complex plain, Moebius transformation has been applied: 𝜇=1+𝜂1−𝜂. (3.8) As a consequence of the transformation, the case of the absolute values of the eigenvalues of matrix 𝐀 less than one is equal to the case of the roots of the transformed polynomial on the left side of the complex plain: |𝜇_𝑖|<1,𝑖=1,…,4⇔Re(𝜂_𝑖)<0,𝑖=1,…,3. (3.9) The transformed characteristic polynomial is 𝑎₃𝜂³+𝑎₂𝜂²+𝑎₁𝜂+𝑎₀=0 (3.10) and the coefficients are 𝑎₃=4−(𝑡_𝑠−2𝑡_𝑑)(2𝑃−𝐼𝑡_𝑠𝑎),₂=4(1−𝑃𝑡_𝑑)+𝐼𝑡_𝑠(4𝑡_𝑑−𝑡_s𝑎),₁=(2𝑃−𝐼(2𝑡_𝑑+𝑡_𝑠))𝑡_𝑠,𝑎₀=𝐼𝑡²_𝑠. (3.11) According to Routh-Hurwitz stability criterion, all of the polynomial coefficients (3.11) and the determinant of matrix 𝐇₂ (3.13) should be positive: 𝑎_𝑖>0,𝑖=0,…,3,det(𝐇₂)>0, (3.12) where 𝐇₂=𝑎₂𝑎₀𝑎₃𝑎₁. (3.13) Extracting the determinant, 𝑎₁𝑎₂−𝑎₀𝑎₃>0. (3.14) Considering 0<𝑡_𝑑<𝑡_𝑠, 𝑃>0 and 𝐼>0, the necessary condition is: 0<𝑡_𝑑(𝑃(𝑡_𝑑+𝑡_𝑠)−2𝐼𝑡_𝑠𝑡²_𝑑−1)−𝑡_𝑠+(𝑃𝑡_𝑑−1)²((𝑡_𝑠−𝑡_𝑑)²−4𝑡_𝑑𝑡_𝑠𝑃𝑡_𝑑−1), (3.15) in addition, one of conditions (3.16)–(3.18) should be also satisfied for the stability: 4𝑡_𝑑=𝑡_𝑠,𝑃<4(𝑡_𝑠−3𝑡_𝑑)(𝑡_𝑠−2𝑡_𝑑)², (3.16) or 4𝑡_𝑑<𝑡_𝑠,𝑃<4(𝑡_𝑠−3𝑡_𝑑)(𝑡_𝑠−2𝑡_𝑑)²,𝐼<4+2𝑃(2𝑡_𝑑−𝑡_𝑠)𝑡_𝑠(𝑡_𝑠−2𝑡_𝑑), (3.17) or 4𝑡_𝑑>𝑡_𝑠1,𝑃<𝑡_𝑑. (3.18) In the above formulas, the continuous extension is to be used, when 𝑡_𝑑→𝑡_𝑠/2. Figure 5 shows parametric stability charts of the linear system at sampling time 𝑡_𝑠=0.1s and for various time delays in the range 𝑡_𝑑∈[0,𝑡_𝑠]. It is easy to prove that the stability boundaries are straight lines when 𝑡_𝑑=0. On the left stability boundary, |𝜇₁|<1 and |𝜇_2,3|=1. It can also be shown that only this type of stability boundary exists when 𝑡_𝑑≥𝑡_𝑠/4. The rightmost stability boundary for 𝑡_𝑑<𝑡_𝑠/4 is always a straight line, where |𝜇₁|<1, |𝜇₂|<1, and |𝜇₃|=1. As the system is overdetermined, 𝜇₄≡0.

Figure 5 

Linear stability chart at sampling time $t_s=0.1\hspace{0.17em}\text{s}$ and time delay $t_d=0\text{,}$$t_d=t_s/7\text{,}$$t_d=t_s/5\text{,}$$t_d=t_s/4\text{,}$$t_d=t_s/3\text{,}$$t_d\to t_s/2\text{,}$$t_d=t_s.$

4. Dynamics of the Piecewise Linear System

According to equation (2.1), the upcoming value of 𝑥_𝑛+1 depends on the actual 𝑣_𝑛+1 and the previous 𝑣_𝑛 values of velocities. Considering Equation (2.5), 𝑥_𝑛+1 is the piecewise linear function of two previous error signals ℎ_𝑛 and ℎ_𝑛−1. Since each past value can fall into 3 cases, our system is described by 9 scalar equations: 𝑥_𝑛+1=𝑥_𝑛+(𝛿−ℎ_𝑛−1)𝑡_𝑑+(𝛿−ℎ_𝑛)(𝑡_𝑠−𝑡_𝑑𝑥),_𝑛+(𝛿−ℎ_𝑛)(𝑡_𝑠−𝑡_𝑑)𝑥_𝑛+(−𝛿−ℎ_𝑛−1)𝑡_𝑑+(𝛿−ℎ_𝑛)(𝑡_𝑠−𝑡_𝑑ℎ),if^𝑛−1>𝛿,|ℎ_𝑛−1|≤𝛿,ℎ_𝑛ℎ>𝛿,_𝑛−1𝑥<−𝛿,_𝑛+(𝛿−ℎ_𝑛−1)𝑡_𝑑,𝑥_𝑛𝑥_𝑛+(−𝛿−ℎ_𝑛−1)𝑡_𝑑ℎif^𝑛−1>𝛿,|ℎ_𝑛−1|≤𝛿,|ℎ_𝑛ℎ|≤𝛿,_𝑛−1𝑥<−𝛿,_𝑛+(𝛿−ℎ_𝑛−1)𝑡_𝑑+(−𝛿−ℎ_𝑛)(𝑡_𝑠−𝑡_𝑑𝑥),_𝑛+(−𝛿−ℎ_𝑛)(𝑡_𝑠−𝑡_𝑑)𝑥_𝑛+(−𝛿−ℎ_𝑛−1)𝑡_𝑑+(−𝛿−ℎ_𝑛)(𝑡_𝑠−𝑡_𝑑ℎ),if^𝑛−1>𝛿,|ℎ_𝑛−1|≤𝛿,ℎ_𝑛ℎ<−𝛿,_𝑛−1<−𝛿. (4.1) The piecewise system can be written in a compact form, where 𝐹_𝑖 is a linear operator, that shortens the matrix formalism: 𝐳_𝑛+1=𝐀_i𝐳_𝑛+𝐛_𝑖=𝐹_𝑖(𝐳_𝑛),𝑖=1,…,9. (4.2) In (4.2), the elements of 𝐀_𝑖 matrices and 𝐛_𝑖 vectors (𝑖=1,…,9) can be calculated as the coefficients of scalar equations listed in (4.1). See Appendix (A.1) for details. The appropriate 𝐹_𝑖 is selected according to the previous two error signals, shown in Table 1.

Table 1 

Linear operator selection based on the error signals.

Extracting ℎ_𝑛 according to (2.2), the 3 intervals of ℎ_𝑛 result in 3 domains in the 𝑥-𝑦 plane: 𝑃𝑥_𝑛+𝐼𝑦_𝑛<−𝛿,−𝛿≤𝑃𝑥_𝑛+𝐼𝑦_𝑛≤𝛿,𝛿<𝑃𝑥_𝑛+𝐼𝑦_𝑛, (4.3) bordered by 2 lines: 𝑃𝑦=−𝐼𝛿𝑥±𝐼. (4.4) Figure 6 represents the 𝑥-𝑦 plane with the dead zone in the middle region. In the dead zone, 𝑥 is constant, meaning that the piston rod is stopped, since the proportional directional valve is closed. Fixed points are on the 𝑦 axis in the interval of [−𝛿/𝐼,𝛿/𝐼]. This invariant set corresponds to the trivial solution of the backlash-free linear system. From practical viewpoint, we are interested in the stability of the invariant set, since the actual value of 𝑦 has no importance if we managed to reach the desired 𝑥=0 position.

Figure 6 

Three domains of the $x\text{-}y$ plane, in the middle region the error signal is less then the threshold of the valve, there is no piston motion.

4.1. Periodic Orbits

Numerical simulations were carried out with sampling time 𝑡_𝑠=0.1s, time delay 𝑡_𝑑=0.04s, and dimensionless dead zone 𝛿=12. By solving the system of algebraic equations shown in (4.8), starting the system from initial condition 𝐳₂=(𝑥₂𝑥₁𝑦₂𝑦₁)^𝑇=(1111)^𝑇, one finds that at the values of integral gain 𝐼_cr=24.982 and proportional gain 𝑃_cr=3.5714, a periodic orbit exists, and the invariant set is stable for 𝑃_cr<𝑃 stable, and unstable for 𝑃<𝑃_cr. This behavior is similar to a degenerate Hopf bifurcation, being neither sub- nor supercritical. Figure 7 represents the stable, unstable runs and one periodic orbit from those, which exist in the critical case.

Figure 7 

Three simulation results, both with initial conditions ${\mathbf{z}}_2={(1\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}1)}^T\text{,}$$t_s=0.1\hspace{0.17em}\text{s}\text{,}$$t_{\text{d}}=0.04\hspace{0.17em}\text{s}\text{,}$$\delta =12$,and $I_{\text{cr}}=24.982$. Periodic orbit occurs at $P_{\text{cr}}=3.5714$ (black); the control is stable at $P=4.1$ (green), unstable at $P=3.3$ (red).

According to the simulation results shown in Figure 8 for different modified dead zone sizes, the periodic orbit remains either periodic or it becomes a quasiperiodic dense orbit.

Figure 8 

Three simulation results, both with initial conditions ${\mathbf{z}}_2={(1\hspace{0.17em}1\hspace{0.17em}1\hspace{0.17em}1)}^T\text{,}$$t_s=0.1\hspace{0.17em}\text{s}\text{,}$$t_d=0.04\hspace{0.17em}\text{s}\text{,}$$I_{\text{cr}}=24.982$ and $P_{\text{cr}}=3.5714$. Periodic orbit occurs at $\delta =12$ (black), dense orbits with increased period length at different dead zone sizes, $\delta =20$ (green) and $\delta =4$ (red).

Figure 9 represents a periodic orbit in a general case. The numbering of the points starts from the first step outside the dead zone region; 𝑘 is the number of steps in the upper region; 𝑗 is half of the number of the steps inside the dead zone. Due to being a symmetric system, the examination of the half of the periodic orbit is satisfactory, with the end point subscript 𝑘+𝑗+2.

Figure 9 

General periodic orbit, $k$ is half of the number of steps outside the middle region, $j$ is half of the number of steps inside the middle region.

Initial conditions are arbitrarily chosen: 𝐳₂=𝑥₂𝑥₁𝑦₂𝑦₁^𝑇. (4.5) For a symmetric half-orbit, we have 𝐳_𝑘+𝑗+2=(−𝑥₂−𝑥₁−𝑦₂−𝑦₁)^𝑇=−𝐳₂. (4.6) The first step is from point number 2 to number 3. In this case, both ℎ₂>𝛿 and ℎ₁>𝛿, therefore operator 𝐹₁ is used to calculate point number 3, just as on the next 𝑘−1 steps. As we reach the dead zone, according to Table 1, 𝐹₄ is applied on the border, and then 𝐹₅ 𝑗−1 times. The last step until the half-cycle is made with operator 𝐹₈. Therefore, the half-cycle is formed in general: 𝐳_𝑘+𝑗+2=𝐹₈(𝐹₅(⋯(𝐹₅(𝐹₄(𝐹₁(⋯(𝐹₁(𝐳₂))))))))−−−→^def.𝐹₈(𝐹₅^(𝑗−1)(𝐹₄(𝐹₁^(𝑘−1)(𝐳₂)))), (4.7) where 𝐹_𝑖^(𝑘) means applying 𝐹_𝑖 operator 𝑘 times. Extracting (4.7), 𝐳_𝑘+𝑗+2=−𝐳₂=𝐀₈(𝐀₅(⋯(𝐀₅(𝐀₄(𝐀₁(⋯(𝐀₁𝐳₂+𝐛₁)+⋯)+𝐛₁)+𝐛₄)+𝐛₅)+⋯)+𝐛₅)+𝐛₈. (4.8) Using (4.8), one can generate periodic orbits as follows. We fix 𝑡_𝑑, 𝑡_𝑠, and 𝛿, furthermore, the initial condition 𝑥₁, 𝑦₁ and the “shape” of the orbit with 𝑘 and 𝑗, then by solving the second and fourth components of (4.8) 𝑃 and 𝐼 can be calculated, with which 𝑥₂ and 𝑦₂ can be easily determined. In Figure 10, three optional periodic orbits are shown.

Figure 10 

Three optional periodic orbits, both orbits step through $x_n=1$ and $y_n=1$, and have $t_{\text{s}}=0.1\hspace{0.17em}\text{s}$ sampling time, $t_d=0.04\hspace{0.17em}\text{s}$ time delay, and $\delta =12$ dead zone width. Black colored orbit is produced with $P_{\text{cr}}=3.5714$ and $I=24.982$, red one is with parameters $P_{\text{cr}}=3.1859\text{,}$$I=22.968$, green one is with parameters $P_{\text{cr}}=4.6586$ and $I=32.259$.

4.2. Stability of Periodic Orbits

We can reduce equation (4.8) to one single operation, introducing 𝐀 and 𝐛: 𝐳_𝑘+𝑗+2=𝐀𝐳₂+𝐛, (4.9) where 𝐀=𝐀₈𝐀₅^𝑗−1𝐀₄𝐀₁^𝑘−1,𝐛=𝐀₈𝐀₅^𝑗−1𝐀₄𝐀₁^𝑘−2𝐛₁+𝐀₈𝐀₅^𝑗−1𝐀₄𝐀₁^𝑘−3𝐛₁+⋯+𝐀₈𝐀₅^𝑗−1𝐀₄𝐛₁+𝐀₈𝐀₅^𝑗−1𝐛₄+𝐀₈𝐀₅^𝑗−2𝐛₅+⋯+𝐀₈𝐀₅𝐛₅+𝐀₈𝐛₅+𝐛₈. (4.10) The stability of this reduced dynamical system shown in (4.9) depends on the eigenvalues of matrix 𝐀. Since 𝛿 does not appear explicitly in 𝐀, the stability boundary of the invariant set including the desired 𝑥=0 position in the system with dead zone coincides with the stability boundary of the linear system with linear valve characteristics derived in Section 3. However, the structure of 𝐀 could change if 𝛿 is large enough, and this way it can still affect the stability regions.

In the left panel of Figure 11, three pairs of those control parameters (𝑃,𝐼) are denoted at the limit of linear stability, where the periodic orbits exist as presented in Figure 10. Note that their numerical values were calculated with the previously described method (solving (4.8)).

(a)

(b)

Figure 11 

Left: linear stability chart at $t_{\text{s}}=0.1\hspace{0.17em}\text{s}\text{,}$$t_{\text{d}}=0.04\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{s}$. Black colored point shows $P_{\text{cr}}=3.5714$ and $I=24.982$; blue one shows parameters $P_{\text{cr}}=3.1859$,$\mathrm{\hspace{0.17em}\hspace{0.17em}}I=22.968$; red one shows parameters $P_{\text{cr}}=4.6586$ and $I=32.259$, both with $\delta =12$. Right: one point is picked on the boundary of linear stability, one inside the stable, one in the unstable region: (a) $P=3.3\text{,}$$I=24.982$, (b) $P=3.5714\text{,}$$I=24.982$, (c) $P=3.9\text{,}$$I=24.982$.

Figure 12 shows the behavior of three periodic or dense orbits for three different values of 𝑃 represented in the right panel of Figure 11. As it was shown above, the stability of the linear system is preserved by the invariant set including 𝑥=0 in the system with dead zone. This means that all orbits will spiral outwards or inwards corresponding to the unstable or stable linear behavior independently whether the orbits are periodic or dense at the critical values of 𝑃. In other words, for each combination of the matrices 𝐀_𝑖 corresponding to any kinds of orbits, the eigenvalues of 𝐀 behave similarly in terms of stability. In the phase space this means that the trajectories inside the dead zone are purely vertical and do not change the value of 𝑥—apart from the steps when the trajectory enters and leaves the dead zone. Roughly speaking, the dead zone only “cuts” and “extracts” an already existing orbit. This also explains why the orbits with dead zone behave similarly to the backlash-free case.

(a)

(b)

(c)

Figure 12 

Simulation results with $\delta =12$ corresponding to the parameters (a), (b), and (c) in Figure 11.

Further investigation is needed to study how these dynamical properties will change with slight perturbation caused by the damping. It is likely that the above structure of periodic and dense orbits will not survive, but some of them may exist either in the linearly stable or unstable domain, depending on whether the damping is slightly positive or negative.

5. Conclusions

In this paper, the PI control of a hydraulic positioning system with cylinder was studied with an emphasis on the interaction of digital sampling, time delay due to finite computational time and internal dynamics, and backlash due to valve characteristics affecting the global dynamics of the controlled system. The stability boundary of the backlash-free system was computed analytically and represented on the (𝑃,𝐼) control parameter plane. Then, it was shown that the dynamics is described by a piecewise linear system with 9 possible states caused by the presence of backlash. It was shown that both periodic and dense orbits are present in the system when the parameters are tuned to the stability boundary of the backlash-free system.

An analytical method was presented which, for a given sequence of switchings, initial conditions, sampling and delay time, and dead zone width, allows the computation of the corresponding 𝑃 and 𝐼 parameters ensuring the existence of a periodic orbit. Moreover, it was shown that as the linear coefficient matrices of the piecewise linear system are independent of the dead zone width, so does the stability of the periodic or dense orbits. In other words, the stability boundary of the linear system provides a practical stability margin for the system with backlash, too.

Appendix

List of 𝐀_𝑖 matrices and 𝐛_𝑖 vectors. 𝐀₁=1−𝑃(𝑡_𝑠−𝑡_𝑑)−𝑃𝑡_𝑑−𝐼(𝑡_𝑠−𝑡_𝑑)−𝐼𝑡_𝑑𝑡1000_𝑠0100010,𝐛₁=𝛿𝑡_𝑠000^𝑇,𝐀₂=1−𝑃(𝑡_𝑠−𝑡_𝑑)0−𝐼(𝑡_𝑠−𝑡_𝑑𝑡)01000_𝑠0100010,𝐛₂=𝛿(𝑡_𝑠−𝑡_𝑑)000^𝑇,𝐀₃=𝐀₁,𝐛₃=𝛿(𝑡_𝑠−2𝑡_𝑑)000^𝑇,𝐀₄=1−𝑃𝑡_𝑑0−𝐼𝑡_𝑑𝑡1000_𝑠0100010,𝐛₄=𝛿𝑡_𝑑000^𝑇,𝐀₅=𝑡10001000_𝑠0100010,𝐛₅𝐀=𝟎,₆=𝐀₄,𝐛₆=−𝐛₄,𝐀₇=𝐀₁,𝐛₇=−𝐛₃,𝐀₈=𝐀₂,𝐛₈=−𝐛₂,𝐀₉=𝐀₁,𝐛₉=−𝐛₁, (A.1)