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

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.


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

Mathematical Problems in Engineering
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.

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.
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 t s sampling time , thus the position of the piston rod x t is discretized in time as x n x nt s .This position is fed into the PC which computes the error signal h.The time needed for this computation is denoted 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 id represents the internal dynamics of the hydraulic system see text for details , t d t c t id is the overall delay, and t s denotes sampling time.by t 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 t id , with subscript referring to internal dynamics.Thus, the overall delay between the previous sampling instant and its effect is t d t c t id .Note that the actual values and ratio of t c and t 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: v n t ≡ v n − 1 t s t d , t ∈ n − 1 t s t d , n t s t d , accordingly, subscripts n refer to different time instants for positions and velocities as it is represented in the controlling scheme of Figure 3.
Assume the cases when we have 0 ≤ t d ≤ t s .Integrating the piecewise constant velocities, we arrive at a relationship between two neighbouring sampled piston rod positions, which can be expressed as 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 h n is calculated in the following form: where P is the proportional gain I is the integral one, and is the discrete integral of the position function.The piston rod velocities v n 1 and v n are calculated from h n and h n−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 −α.
We introduce the dimensionless variables by means of 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,

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 t d > t s , 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 ≤ t d ≤ t s .

Stability Analysis of the Linear System
Eliminating the valve dead zone δ 0 , a linear valve characteristic means regarding to the piston rod velocity: 3.1 from equation 2.1 : Introducing z n , that consists of the actual and the previous piston positions, and integral values: with matrix formalism: The stability of the system depends on the eigenvalues of matrix A, all of the absolute values of the eigenvalues have to be less than 1:

3.6
The characteristic polynomial of matrix A is The polynomial has one root that equals to zero; therefore; it can be divided by μ / μ 4 0. Since the stability criteria of polynomials are determining the coefficients of polynomial so Mathematical Problems in Engineering that all of the roots should be on the left side of the complex plain, Moebius transformation has been applied: As a consequence of the transformation, the case of the absolute values of the eigenvalues of matrix A less than one is equal to the case of the roots of the transformed polynomial on the left side of the complex plain: The transformed characteristic polynomial is and the coefficients are

3.11
According to Routh-Hurwitz stability criterion, all of the polynomial coefficients 3.11 and the determinant of matrix H 2 3.13 should be positive: where

3.14
Considering 0 < t d < t s , P > 0 and I > 0, the necessary condition is:

Dynamics of the Piecewise Linear System
According to equation 2.1 , the upcoming value of x n 1 depends on the actual v n 1 and the previous v n values of velocities.Considering Equation 2.5 , x n 1 is the piecewise linear function of two previous error signals h n and h n−1 .Since each past value can fall into 3 cases, our system is described by 9 scalar equations: 4.1 The piecewise system can be written in a compact form, where F i is a linear operator, that shortens the matrix formalism: In 4.2 , the elements of A i matrices and b i vectors i 1, . . ., 9 can be calculated as the coefficients of scalar equations listed in 4.1 .See Appendix A.1 for details.The appropriate F i is selected according to the previous two error signals, shown in Table 1.
Extracting h n according to 2.2 , the 3 intervals of h n result in 3 domains in the x-y plane: Figure 6 represents the x-y plane with the dead zone in the middle region.In the dead zone, x is constant, meaning that the piston rod is stopped, since the proportional directional valve is closed.Fixed points are on the y axis in the interval of −δ/I, δ/I .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 y has no importance if we managed to reach the desired x 0 position.

Periodic Orbits
Numerical simulations were carried out with sampling time t s 0.1 s, time delay t d 0.04 s, and dimensionless dead zone δ 12.By solving the system of algebraic equations shown in 4.8 , starting the system from initial condition z 2 x 2 x 1 y 2 y 1 T 1 1 1 1 T , one finds that at the values of integral gain I cr 24.982 and proportional gain P cr 3.5714, a periodic orbit exists, and the invariant set is stable for P cr < P stable, and unstable for P < P 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.
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 9 represents a periodic orbit in a general case.The numbering of the points starts from the first step outside the dead zone region; k is the number of steps in the upper region; j 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 k j 2.
Initial conditions are arbitrarily chosen: For a symmetric half-orbit, we have The first step is from point number 2 to number 3. In this case, both h 2 > δ and h 1 > δ, therefore operator F 1 is used to calculate point number 3, just as on the next k − 1 steps.As we reach the dead zone, according to Table 1, F 4 is applied on the border, and then F 5 j − 1 times.The last step until the half-cycle is made with operator F 8 .Therefore, the half-cycle is formed in general: where Using 4.8 , one can generate periodic orbits as follows.We fix t d , t s , and δ, furthermore, the initial condition x 1 , y 1 and the "shape" of the orbit with k and j, then by solving the second and fourth components of 4.8 P and I can be calculated, with which x 2 and y 2 can be easily determined.In Figure 10, three optional periodic orbits are shown.

Stability of Periodic Orbits
We can reduce equation 4.8 to one single operation, introducing A and b: z k j 2 A z 2 b, 4.9  where

4.10
The stability of this reduced dynamical system shown in 4.9 depends on the eigenvalues of matrix A. Since δ does not appear explicitly in A, the stability boundary of the invariant set including the desired x 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 A 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 P, I 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 .
Figure 12 shows the behavior of three periodic or dense orbits for three different values of P 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 x 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 P .In other words, for each combination of the matrices A i corresponding to any kinds of orbits, the eigenvalues of A 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 x-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.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.

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 P, I 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 P and I 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 A i matrices and b i vectors.

Figure 2 :
Figure 2: Dead zone in a typical flow rate/control signal characteristics of a directional proportional hydraulic valve.In the interval −u min , u min there is no fluid flow. x

Figure 6 :
Figure 6: Three domains of the x-y plane, in the middle region the error signal is less then the threshold of the valve, there is no piston motion.

Figure 8 :
Figure 8: Three simulation results, both with initial conditions z 2 1 1 1 1T , t s 0.1 s, t d 0.04 s, I cr 24.982 and P cr 3.5714.Periodic orbit occurs at δ 12 black , dense orbits with increased period length at different dead zone sizes, δ 20 green and δ 4 red .

Figure 9 :
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.

Figure 10 :
Figure 10: Three optional periodic orbits, both orbits step through x n 1 and y n 1, and have t s 0.1 s sampling time, t d 0.04 s time delay, and δ 12 dead zone width.Black colored orbit is produced with P cr 3.5714 and I 24.982, red one is with parameters P cr 3.1859, I 22.968, green one is with parameters P cr 4.6586 and I 32.259.

Figure 11 :
Figure 11: Left: linear stability chart at t s 0.1 s, t d 0.04 s.Black colored point shows P cr 3.5714 and I 24.982; blue one shows parameters P cr 3.1859, I 22.968; red one shows parameters P cr 4.6586 and I 32.259, both with δ 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, I 24.982, b P 3.5714, I 24.982, c P 3.9, I 24.982.

Figure 12 :
Figure 12: Simulation results with δ 12 corresponding to the parameters a , b , and c in Figure 11.

1 − 1 − 4 , b 6 −b 4 ,A 7 A 1 , b 7 −b 3 ,A 8 A 2 , b 8 −b 2 ,A 9 A 1 , b 9 −b 1 ,A. 1
P t s − t d −Pt d −I t s − t d −It d P t s − t d 0 −I t s − t d t s /4, t d t s /3, t d → t s /2, t d t s .inaddition,one of conditions 3.16 -3.18 should be also satisfied for the stability:In the above formulas, the continuous extension is to be used, when t d → t s /2.Figure5shows parametric stability charts of the linear system at sampling time t s 0.1 s and for various time delays in the range t d ∈ 0, t s .It is easy to prove that the stability boundaries are straight lines when t d 0. On the left stability boundary, |μ 1 | < 1 and |μ 2,3 | 1.It can also be shown that only this type of stability boundary exists when t d ≥ t s /4.The rightmost stability boundary for t d < t s /4 is always a straight line, where |μ 1 | < 1, |μ 2 | < 1, and |μ 3 | 1.As the system is overdetermined, μ 4 ≡ 0.
Figure 5: Linear stability chart at sampling time t s 0.1 s and time delay t d 0, t d t s /7, t d t s /5, t d

Table 1 :
Linear operator selection based on the error signals.

1 Figure 7 :
Three simulation results, both with initial conditions z 2 1 1 1 1 T , t s 0.1 s, t d 0.04 s, δ 12,and I cr 24.982.Periodic orbit occurs at P cr 3.5714 black ; the control is stable at P 4.1 green , unstable at P 3.3 red .
1, and have t s 0.1 s sampling time, t d 0.04 s time delay, and δ 12 dead zone width.Black colored orbit is produced with P cr 3.5714 and I 24.982, red one is with parameters P cr 3.1859, I 22.968, green one is with parameters P cr 4.6586 and I 32.259.
1 s, t d 0.04 s.Black colored point shows P cr 3.5714 and I 24.982; blue one shows parameters P cr 3.1859, I 22.968; red one shows parameters P cr 4.6586 and I 32.259, both with δ 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, I 24.982, b P 3.5714, I 24.982, c P 3.9, I 24.982.