Simplified modelling and backstepping control of the long arm agricultural rover

This paper presents the development of the simplified modelling and control of a long arm system for an agricultural rover, which also extends the modelling methodology from the previous work. The methodology initially assumes a flexible model and, through the use of the integral-based parameter identification method, the identified parameters are then correlated to an energy function to allow a construction of the friction induced nonlinear vibration model. To also capture the effect of the time delay, a delay model was also considered in the form of a second order delay differential equation. Both families of models were applied to identify and characterise a specialised long arm system. The nonlinear model was found to give significant improvement over the standard linear model in data fitting, which was further enhanced by the addition of the time delay consideration. A backstepping controller was also designed for both model families. Results show that the delay model expends less control efforts than the lesser non-delay model.


Introduction
The advent of robotics technology in recent decades has fuelled rapid growth in agricultural robotics, not only to meet the increasing demands for alternatives to human labour in agricultural production due to the difficulty of finding and retaining workers [1], but also to satisfy environmental and food safety needs [2]. Such growth has garnered significant research in recent years [3,4]. Different types of agricultural robots have been developed from the days of the Gerrish tractor robot in 1984 [5]. These robots operate in a wide range of agriculture processes, including harvesting [6][7][8], weed control [9][10][11] and spraying [12][13][14]. It is obvious that the use of a robot arm is essential to reaching the required targets.
Agricultural processes such as harvesting and spraying on tropical fruits such as rambutan (Nephelium lappaceum) and durian (Durio zibethinus), however, necessitate a mobile robotic platform with elongated arms that appends to at least four metres in height. As the arms themselves move, significant vibrations are felt at the tip of the arms which must be mathematically modelled and controlled, an issue also prevalent on the computer numerical control (CNC) machine tools in industrial robots [15].
Oscillatory vibrations usually assume a linear damping model, in which the responses are readily decomposed into various modal frequencies. The common algorithm that follows this line of approach is Prony's method [16,17]. This approach is similar to the wellknown concept of the Fourier transform, except that the exponential decay term is added to the trigonometric basis function. Variants of Prony's method include the use of the total least squares technique instead of the ordinary least squares [18] and the matrix pencil method [19]. Other approaches have also included the Kalman-based estimations [20][21][22], the distributed frequency domain optimisation [23] and the second order generalised integrator [24,25].
More advanced types of vibration modelling also include friction induced models, which can be separated into two main types. The first type views the vibration from a tribological viewpoint, where the variation of the friction coefficient changes against the relative velocity between the structures in contact. Such variation could be described by the Stribeck model, polynomial or even exponential functions [26,27]. The second type gives emphasis to the structural aspects in the process of vibration generations. In this regard, the prediction of the FIV has been conducted through sensitivity approaches and probabilistic models [28,29], hybrid meta models [30] and fuzzy approaches [31]. More recent approaches to estimate the FIV include the use of observer designs to estimate the states [32][33][34].
The surveyed models and methods generally tend to presuppose a definite complicated model structure first, then fit a complex model to the data. Should this initial model assumption fail to capture the responses, then more complex observers are imbued. The approach taken to analyse the vibration of the agricultural mobile rover in this work extends the concept of the work done by the first author in 2015 [35]. In other words, Sect. 2 of this paper presents the integral-based identification methodology for the vibration model without including the delay effects. The same formulation and concept is then used for the vibration model with the delay being included. A corresponding integral-based identification method for the linear model is also presented for the purpose of comparison. Theoretical analyses are also given on the Lyapunov stabilities of the delay as well as the non-delay models. In Sect. 3, the developed model and identification methods are applied to specialised rover arm system, where more complex dynamics are uncovered based on the identified damping and stiffness values. A validation procedure of the developed model is also given as proof of concept. A backstepping controller is designed in Sect. 3.6 based on the developed model of Sect. 3.2. This paper is then concluded in Sect. 4.

Modified integral-based identification for the long arm rover without delay
Consider the normalised second order differential equation: where: θ ≡ θ (t) represents angular movement, c is the damping constant, k is the spring constant, u is the input.
To designate some degree of flexibility for the model of damping as well as stiffness, these variables can be made time-varying. The simplest such time-varying function is the piecewise constant function defined as follows: = c n , T n-1 < t < T end and = k n , T n-1 < t < T end .
The system of Equation (1) can now be rewritten as follows: where the functions c(t) and k(t) are defined in Equations (2) and (3). Furthermore, the following quantities are also defined: For future references, let the measurement times t (j) i be defined as follows: i ≡ Measurement instant t i of section j, i = 1, . . . , N and j = 1, . . . , n.
Defining also the operator: Applying the operator of Equation (8) onto Equation (4) yields Here, the initial conditions are defined as follows: The integral reconstruction model for the vibration system can now be written as follows:  Figure 1 shows a possible scenario at the joins in the neighbourhood of t = T i , i = 1, . . . , n.
In practice, the angle measurements θ meas will be obtained from an encoder, whose additive quantisation noise implies that the value of θ meas may not equal the θ model at the joins.
This phenomenon introduces possible discontinuities at the joins. A simple method of resolving this discontinuity is to proceed to identify the sections by piece and resolve the discontinuity at the end of every section. The integral reconstructor for the first section is written as follows: Equations (15)- (17) are solved by linear least squares subjected to the constraints The result yields the unknown parameters that are the elements of vector p 0 . Substituting the elements of p 0 into Equation (13) yields an integral reconstructor model for the angle θ (t) of the first section from the data segmentation.
To ensure that the C 0 continuity is ensured at join k = i -1, the initial condition θ i can be computed thus: The value of b, which was identified for the first section, is assumed to be constant for all sections. This assumption is based on the fact that the arm is balanced and thus no sudden change of inertia is possible across the different time sections. The knowledge of the initial conditions for the ith section now leaves only three unknowns to be identified in the model of Equation (11). In this respect, setting θ (t) ≡ θ model,i (t) and input u(t) ≡ u data,i (t) for the time instants t ∈ {t (i) 1 , . . . , t (i) N } yields a system of N equations in three unknowns: . . .
Equations (20)- (22) can now be solved by linear least squares subjected to the constraints The result of such an operation will now provide the values for the damping c i and stiffness k i along with the integral reconstructor model for the angular movement θ model (t). Figure 2 shows the algorithm for identifying the time-varying non-delay model.

Integral-based identification for the long arm rover with delay
As a comparison, consider now the normalised second order delay differential equation given by where τ represents the delay, and the variable θ and the associated parameters c and k retain their meaning from the non-delay case in Sect. 2.1. Inserting also the time-varying models for the damping and stiffness as provided by Equations (2) and (3) gives Note that the time segmentation intervals for Equation (25) are again given by Equations (5) and (6). The time delay functions θ (tτ ) and their derivatives are usually difficult to model. However, under the assumption that the time delay τ is small, it is possible to approximate θ (tτ ) and their derivatives by Taylor approximations: Substituting Equation (26a)-(26c) into Equation (25), we obtain: where Applying the operator of Equation (8) to Equation (27) gives where the initial conditions are The discontinuities at every section joins are resolved using the method presented in Sect. 2.1. In this light the integral reconstructor model for the first section is written as follows: Substituting the values of θ ≡ θ data (t) and u(t) ≡ u applied (t) for the times of t ∈ {t (1) 0 , . . . , t (1) N } will give a matrix equation where θ data (t (1) j ) denotes the angular data at t = t j of the first section. The system identification algorithm begins with the solving of Equations (35)-(37) by linear least squares subjected to the conditions The result of this process yields the unknown parameters that belong to the elements of p 0 , whose vector is then substituted into Equation (34) to obtain the integral reconstructor model for the angle θ (t) for the first section of the segmentation. The initial conditions for the beginning of the ith segmentation are evaluated thus: Again note that the variable b is assumed to be constant throughout all sections. The knowledge of the initial conditions implies that only six parameters are now required to be identified in the model of Equation (30). Setting θ (t) ≡ θ model,i (t) and input u(t) ≡ u data,i (t) for the time instants t ∈ {t (i) 1 , . . . , t (i) N } yields the matrix equation . . .
Equations (40)-(42) can now be solved by linear least squares subjected to the constraints The result of such an operation will now provide the values for the parameters a i,1 , a i,2 and a i,3 along with the integral reconstructor model for the angular movement θ model (t). Figure 3 shows the algorithm for identifying the time delay model. Note that Equations (43b)-(43d) place the Lipschitz constraints on the derivatives of a i,1 , a i,2 and a i, 3 to make sure that they are bounded. Were these constraints not placed on the derivatives, it would be possible to choose very small such that the modelled response y model (t) gets very close Algorithm 2: algorithm for identifying the time-varying a 1 (t), a 2 (t) and a 3 (t) functions of Equation (25) to y true (t), yet the values of a i,1 , a i,2 and a i,3 do not resemble the true function. In fact, previous work by Wongvanich et al. [35] has shown that the identified parameters significantly oscillate without bound about the true values, yet the modelled response matches very well with the true data.

Analyses
This section gives the theoretical analyses for the vibration model, both without the delay and with the delay.

Model without delay
Lemma 1 Consider the homogeneous second order differential equation The required Lyapunov function can be written as follows: and Proof The proof of this lemma is similar to the one given in [36].
Hence the resulting upper bound Q is max sup

Model with delay Lemma 3 Consider the third order homogeneous system
r ≡ r(t) and 0 < r(t) < r m .
The required Lyapunov function is written as follows: where the states are chosen as follows: Proof The proof of this lemma is similar to the one given in [37].

Theorem 4
Consider the third order differential equation Proof The following Lyapunov function is written by applying Lemma 3: To keep the Lyapunov function of Equation (55) finite, first choose M 1 and M 2 so that the coefficients of (z 1 -z 2 2 F 1,true ) 2 and z 2 2 are finite: The resulting upper bound M is thus: Having established the global asymptotic stability for the system of Equation (27), it is now possible to establish the convergence of our integral reconstructor model. In this respect, we propose the following theorem.
Theorem 5 Consider the following third order nonlinear differential equation: where Define also the following functions: If the parameter a k i,j,model , i = 1, 2, 3, are functions satisfying the Lipschitz condition, that is, The limit of θ model,n will approach the true angular function θ true . In addition, Proof The integral reconstructor model for the system of Equation (56) is as follows: We will firstly consider the case where t = 0. In this case, there existN and δ > 0 such Since the constituent functions of F i,model,k , i = 1, 2, and K i,model,k are a i,model,k , i = 1, 2, 3, which are Lipschitz functions, there will exist a time t ∈ [0, dt * ] regardless of k such that F i,model,k and K i,model,k will intersect with the true functions F i,true and K true respectively. Therefore, And if θ (0) ≤ 0, then it is possible to choosedt * < dt * such that the value of θ true (t) ≤ 0, t ∈ [0, dt * ]. Hence the value of F i,model,k -F i,true (t) and K model,k -K true or F true -F i,model,k and K true -K i,model,k cannot change sign in that time. Thus, the error between the integral reconstruction function and the true value is as follows: Equation (61) contradicts the assumption that the limit of θ model,k will approach θ true .
For the case of t > 0, we also prove by contradiction. Suppose now that there exists a time t 0 > 0, which is the smallest time such that F 1,model,k does not approach F 1,true (0), F 2,model,k does not approach F 2,true (0) and K model,k does not approach K true (0). Hence, Using the same concept as for the case of t = 0, since the functions a i,j,model , i = 1, 2, 3 are Lipschitz function, there exists dt * > 0 such that The above statement implies that it is possible to find the smallest t < t 0 such that the limit of F 1,model,k will not approach F 1,true (t), the limit of F 2,model,k will not approach F 2,true (t) and K model,k will not approach K 2,true (t). This contradicts the assumption that t 0 is the smallest such value. Statement is thus proved for the case of t > 0.

Identification with the pure linear model
The analyses given in the previous section means that the use of constraints on the identified parameters are possible. This conjecture is due to the fact that the nonlinear least squares problem has duly been converted to a corresponding linear least squares problem by integral reconstruction, that is guaranteed to yield a unique solution. Two possible constraints exist for the models considered in this work, one for the non-delay model and the other for the delay model.

Non-delay model
For the non-delay model, the simplest such constraint is simply the equality constraints c 1 = c 2 = · · · = c n and k 1 = k 2 = · · · = k n . In other words, the damping and stiffness parameters are assumed constant across the entire data set. This constraint also represents a full linear model without delay. Setting θ (t) = θ model (t) and u(t) = u applied (t) for the times of t ∈ {t 0 , . . . , t N } gives a matrix equation which is written as follows: Equations (64)-(66) can duly be solved by linear least squares to yield the linear model without delay.

Delay model
Applying the equality constraints to the delay case means that a 1,1 = a 1,2 = · · · = a 1,n , a 2,1 = a 2,2 = · · · = a 2,n and a 3,1 = a 3,2 = · · · = a 3,n . Setting θ (t) = θ model (t) and u(t) = u applied (t) for the times of t ∈ {t 0 , . . . , t N }, as was done for the non-delay case, now gives a matrix equation which is written as follows: where M lin,d = 1 (tt 0 ) (tt 0 ) 2 -I 1,2,3 u appl , In a similar fashion to the non-delay model, Equations (67)-(70) are again solved by linear least squares. The result from this solving process yields the unknown parameters as well as the model for the delay case. Figure 4 depicts the algorithm for identifying the linear model, both for the non-delay and the delay cases. Figure 5 shows the schematic of the specialised robot arm system used for orchard spraying. The arm contains two links where the first link is denoted by L 1 and the second link A separate data acquisition system is developed where an inertial measurement unit (IMU) is fixated on the centre of gravity of both arms. The IMU used has an on-board accelerometre as well as a gyroscope, and it is also connected to a microcontroller via an I2C protocol, where the microcontroller of the acquisition system is also connected to the workstation via a wireless transmitter. Once a PWM command is sent from a workstation, the acquisition system then transmits the angular velocity ω, as well as the angular acceleration α signals for the roll, pitch and yaw axes to the workstation.

Data preparation and preprocessing
The data used are the angular velocity data for the roll, pitch and yaw axes for both links, where the command is given by a step function from 25 degrees and ends at 35 degrees. The command u(t) in symbols is given as follows: where H(t) is the heaviside step function and t s = 8.5 s is the time of the step change. Figure 6 shows the raw data resulting from the command of Equation (71) for L 1 . Figure 7 depicts the raw data resulting from the command of Equation (71) for L 2 . Since these data will have to be integrated with respect to time to yield the angular position for the three axes, incurring integration drift in the process, a compensation mechanism must be in place. In this respect consider the following model for the roll, pitch and yaw axes of both links: Figure 7 Raw angular velocity data for roll, pitch and yaw axes of link L 2 θ L 2 ,roll (t) = a L 2 ,roll + b L 2 ,roll t.
As an example, to find out the parameters of Equations (72a), setting θ L 1 ,yaw ≡ θ L 1 ,yaw,meas (t) for t ∈ {t 0 , . . . , t N } yields the matrix equation Solving Equations (74) by linear least squares yields the parameters of vector p which can be substituted into Equation (72a) to obtain the model for θ L 1 ,yaw . The model is then subtracted from the numerically integrated ω L 1 ,yaw data to achieve the required data for parameter identification. This procedure can be reiterated for θ L 1 ,roll and θ L 1 ,pitch and θ L 2 ,roll , θ L 2 ,pitch and θ L 2 ,yaw . The result of this preprocessing is shown in Figs. 8-9. It is apparent from Figs. 8-9 that a step command at joint J 2 from 25 to 35 induces significant vibration which is most visible on the roll and yaw axes of link L 2 . The trend for both signals is an exponential decay, which can be fitted by the methods presented in Sect. 2. The pitch vibration for both links, however, is infinitesimal; its signal is consumed by the accelerometer noise and is therefore unusable. Hence our result presentations and discussion will focus on the application of the methods presented in Sect. 2 to Link L 2 .

Identification with non-delay model
Consider the application of Algorithm 1 on the roll and yaw data of L 2 , which is the most visible. The data was taken from t = 8.5 s to t = 12 s, with a sampling period of T s = 0.02 s. The value of the time interval t used is taken to be 0.1 s, which gives N = 35 values for c(t) and k(t). The values are then plotted in Figs. 10-11 for the roll and yaw axes, respectively.
It is apparent from Figs. 10 and 11 that both c roll and c yaw exhibit an increasing trend. The stiffness k roll shows an exponentially decreasing trend as time progresses, while k yaw does not increase beyond k yaw = 15. These trends show that it is possible to correlate c roll , c yaw , k roll and k yaw to energy used in order to accurately model the responses.

Damping and stiffness as a function of energy
A change in damping and stiffness with respect to time as was seen in Figs. 10 and 11 suggests that there must also be a change of kinetic energy used by the system, possibly to overcome the stiction in the gears as well as mechanical backlashes. To model these changes in energy, suppose that the damping and stiffnesses are modelled in terms of the kinetic energy-like function To find the relationships of the form given in Equations (76) and (77), the dampings and stiffnesses are then plotted against the average velocity-squared function defined as fol-  Figures 12 and 13 plot the damping and stiffnesses against the average kinetic energy-like functions. It is seen from these figures that as the value of the kinetic energy-like function approaches zero, the damping c roll quickly increases with respect to the energy-like function. This slope becomes significantly flatter as the value of the energy-like approaches about 20. The value of k roll exhibits a diode-like phenomenon, that is, the k roll value is close to zero initially, and it exponentially increases once the kinetic energy-like function approaches 40. This phenomenon is again seen in the function of c yaw , where its value   is close to zero initially and increases exponentially with a decreasing slope as its value reaches about 0.7. The values of k yaw are close to zero initially and quickly approach the value of 15 once the kinetic energy-like function reaches 0.2. Each of the four parameters exhibits a nonlinear phenomenon, which could be mathematically described by a nonlinear piecewise function The parameters of Equations (79)-(82) were identified through nonlinear regression analyses. Figures 14-15 compare the c roll , k roll , c yaw and k yaw data against the fitted relations, showing a close fit between the parameters and the models.
Once the nonlinear relationships between c roll , k roll , c yaw and k yaw have been modelled, they can now be substituted into the original differential equation of Equation (1) to construct a nonlinear model for the long arm rover. To measure the performance of the matches, a mean absolute error metric is used. Figures 16-17 plot the true angular movement data for the yaw and roll axes against the resimulations of Equation (1), with c roll , k roll , c yaw and k yaw given in Equations (79)-(82). It is apparent from both figures that the model provides an accurate match to the data. Specifically, the percentage match of the roll data is 97.8%, while the corresponding match for the yaw data was 97.5%.

Identification with the delay model
Consider now the application of Algorithm 2 on the same angular roll and yaw data of Fig. 9. Here the data is taken from t = 8.5 s to t = 12 s. The value of t used is chosen to be 0.1167 s, yielding N = 30 values for c(t) and k(t). These values are again plotted in Figs. 18 and 19 for the roll and yaw axes respectively. It is seen from Figs. 18 and 19 that the value of a 2,yaw (t) is small when t < 0.5 seconds and saturates at a 2,yaw (t) = 15. The functions a 1,roll (t) and a 1,yaw (t) take on very small values, but exponentially increase as time progresses. The value of a 3,roll (t), however, stays at 15 for about 2.5 seconds, and since then exponentially decreases to zero. These figures suggest that it is possible to explain these phenomena of the changes in a 1 , a 2 , a 3 parameters through correlating with the kinetic-like function as was done with the non-delay case. Specifically, we again suppose that the a 1 , a 2 and a 3 functions are defined in terms In this light, Figs. 20-21 plot the a 1,yaw , a 2,yaw , a 3,yaw , a 1,roll , a 2,roll and a 3,roll responses as a function of the average kinetic energy. It is apparent from these figures that both a 1,yaw and a 1,roll functions are close to zero for energy levels less than 0.15, and they increase at an exponential rate thereafter. Both functions appear to saturate at a higher energy point. The functions a 2,yaw and a 2,roll and a 3,yaw and a 3,roll all follow hyperbolic tangent trends, where the parameters take on a very small value for low kinetic energy and exponentially increase once a threshold value is reached, before saturating at a higher energy value.
The six functions clearly exhibit nonlinear phenomena, which are again explained by nonlinear piecewise functions of the form Figure 19 The identified coefficients for the roll axis a 3,roll = α a3,roll tanh β a3,roll v 2 + γ a3,roll .
The parameters of Equations (85)-(90) were again identified through nonlinear regression analyses. Figures 22-23 compare the a i,yaw and a i,roll , i = 1, 2, 3, data against the fitted models, illustrating a close fit between the two. The nonlinear models for the row and yaw axes can now be written as follows: θ roll (t) + a 1,roll (kE)θ roll (t) + a 2,roll (kE)θ roll (t) + a 3,roll (kE)θ roll (t) = bu(t).
(91) θ yaw (t) + a 1,yaw (kE)θ yaw (t) + a 2,yaw (kE)θ yaw (t) + a 3,yaw (kE)θ yaw (t) = bu(t). (91) and (92) against the angular movement data. It is apparent that the fit was again very close. Specifically the mean absolute error of the fit was 0.5% for the yaw angle and 0.4% for the roll angle. Hence the nonlinear time delay model yields a more accurate description of the long arm rover than the non-delay counterparts.

System identification of the linear model
Consider the application of Algorithm 3 of Fig. 4 to the roll and yaw data of Fig. 11 with the use of the non-delay model. The identified parameters are c r = 8.48 and k r = 30.2. The corresponding identified parameters for the yaw data are c y = 4.99 and k y = 19.87. Figure 26 plots the comparisons for the response data against their identified models. It is seen that the non-delay linear model captures the descent of both angles accurately, but could not capture the oscillation that occurs after t = 0.5 seconds.

Figure 20
The identified coefficients for the roll axis Figure 27 compares the responses between the true measured data against both linear models. It is seen that the delay model did indeed improve the fit of the non-delay model, as evidenced by an improved value of the mean squared error of 0.0028 for the yaw angle response, compared to 0.0036 for the non-delay model. Both models, however, could not capture the oscillations occurring after t = 1.5 seconds.

Validation
An important step in the modelling and identification is the process of validation. To subject the proposed algorithms to the validation test, the algorithms are separately applied to the step responses with the inputs being the following functions: where t s,40 = 5.6 s and t s,60 = 4.7 s are the step change times for each of the data sets. The identified parameters for the linear constant damping model without delay are as follows: For the purpose of comparison, we define the following models for the linear and nonlinear categories. For the linear models without delay, the models are: For the nonlinear models with delay, the models are: NLMD yaw,j ≡ Model of Equation (27) with a 1 = a 1,yaw,j v 2 , NLMD roll,j ≡ Model of Equation (27) with a 1 = a 1,roll,j v 2 , To compare the linear models without delay against the nonlinear models, the models LM yaw,40 and LM roll,40 , LM yaw,60 and LM roll,60 , as well as their nonlinear counterparts, are tested against the 25-35 degrees step change data for the yaw and roll axes, respectively. Similarly, the 40-60 degrees step change data is used as validation data for the models LM yaw,25 , LM roll,25 , LM yaw,60 and LMroll, 60, as well as their nonlinear counterparts. Lastly, the 60-80 degrees step change data is used as validation data for the models LM yaw,25 , LM roll,25 , LM yaw,40 and LM roll,40 and their nonlinear counterparts. Table 1 depicts the mean  Model match between the identified nonlinear model with delay against the measured roll angle response absolute error for the validating data against the models for the roll responses. Table 2 compares the mean absolute error for the validating data against the models for the yaw responses, respectively. It is apparent from these tables that the nonlinear models gave significantly less errors for both the yaw and roll responses, for all the validating data.
To compare the linear models with delay, the models LMD yaw,40 and LMD roll,40 , LMD yaw,60 and LMD roll,60 , LMD yaw,80 and LMD roll,80 , and their nonlinear counterparts,  are again tested against the 25-35 degrees step change data, as was done with the nondelay case. The 40-60 degrees step change data is used as validation data for the models LMD yaw,40 , LMD roll,40 , LMD yaw,60 , LMDroll, 60, LMD yaw,80 and LMD roll,80 , as well as their nonlinear versions. The 60-80 degrees step change data is used as validation data for the     Table 3 depicts the mean absolute error for the validating data against the models for the roll responses, while Table 4 gives the mean absolute errors for the validation data against the yaw responses. It is again apparent that the nonlinear models gave significantly less errors than the corresponding linear models. Note also that the these errors are also less than the corresponding errors for the non-delay case, thereby completing the validation for the proposed models. Note also that the concept undertaken in this paper, for both families of models, is to initially assume a flexible model structure, and through the use of the system identification, unveil yet more complicated relationships between the underlying physical quantities. This concept is different to the ones normally seen in the literature, where a complex structure of vibration induced models, including finite elements and statistical distributions, must firstly be assumed.

Backstepping controller design
Once the models for the long arm rover have been attained, a controller can now be designed for both the non-delay model and the model with delay. To simplify the controller process and as proof of concept, the roll and yaw angles are controlled separately. Since both models could be transformed into a strict feedback form, the backstepping controller is then suitable. This section thus details the control design and presents some results.

Non-delay model
Consider the system of Equation (1) with the understanding that the variable θ is taken to mean either θ roll or θ yaw , the parameters c ≡ c roll (kE) or c yaw (kE) and k ≡ k roll (kE) or k yaw (kE). The states are thus: The state space equations are theṅ In order to stabilise the first equation, we can then define two new state variables where the variable α 1 represents a virtual controller. Define a Lyapunov candidate function Differentiating Equation (102) along the trajectory of Equation (101a)-(101b) yieldṡ The variable z 2 inside the parenthesis of Equation (103) will be designed in the next stage of the controller to be zero. Hence, designing the variable α 1 = -c 1 z 1 , where c 1 > 0, and with the assumption that z 2 = 0, we will then haveV 1 = -c 1 z 2 1 < 0, which means x 1 will now be globally exponentially stable.
The derivative of z 2 with respect to time iṡ The new Lyapunov function V 2 is now defined as follows: Differentiating Equation (105) along the trajectory of Equation (101a)-(101b) yieldṡ Designing u(t) as u(t) = 1 b (cx 2 + kx 1z 1c 1 α 1 -(c 1 + c 2 )z 2 ) will then make the derivative of the Lyapunov function equal toV 2 = -c 1 z 2 1c 2 z 2 2 < 0, ensuring the global asymptotically stability of the closed loop system.

Delay model
The generic form of the system to be controlled under the delay model is Again it is understood that a 1 ≡ a 1,yaw (KE) or a 1,roll (KE), a 2 ≡ a 2,yaw (KE) or a 2,roll (KE), a 3 ≡ a 3,yaw (KE) or a 3,roll (KE). The goal of the control is to design a controller u(t) that brings θ (t) to zero with global asymptotic stability. The state variables for the system of Equation (107) are The system of Equation (107) can be rewritten as follows: In order to stabilise the first equation, we can then define two new state variables as follows: where the variable α 1 in Equation (108) is the virtual controller. Defining the Lyapunov candidate function in much the same way as was defined in Equation (102) and differentiating the Lyapunov function along the trajectory of the system yieldṡ Designing α 1 = -c 1 z 1 will thus makė Since we are designing for z 2 to eventually reach zero, the value ofV 1 will be less than zero, and thus the top equation is stabilised. Furthermore, We next design another virtual controller to stabilise the second equation. In this light the third new state variable is defined as follows: The second Lyapunov candidate function is defined as follows: The derivative of Equation (111) with respect to time along the trajectory of the system iṡ The value ofż 2 iṡ Substituting the value ofż 2 in Equation (112) yieldṡ Designing α 2 to be results inV 2 that is guaranteed to be negative definite, meaning z 2 will be stabilised. Furthermore, α 1 = -c 1ż1 -(c 1 + c 2 )ż 2 = -c 1 (z 2 + α 1 ) + (c 1 + c 2 ) x 3 + c 1 (z 2 + α 1 ) .

Controller response results
For the non-delay model, Fig. 28 shows the results of the backstepping control responses when c 1 = 5 and c 2 = 5. It is seen that the yaw angle increases from zero to about 0.02 before settling at zero by t = 3 seconds. The roll angle gradually decays and finally settles to zero at around one second. Figure 29 shows the case of using c 1 = 20, c 2 = 100. The yaw angle in this case decreases to zero within 1 second, while the roll angle decreases from about 2 degrees to zero in about 1 seconds, all the while staying at the zero equilibrium.
For the use of the delay model, Fig. 30 shows the case of using c 1 = 10, c 2 = 5, c 3 = 5. In this case the yaw angle deviates about 0.008 degrees before returning to zero in three seconds. The roll angle decreases from about 4 degrees to zero within one second and stays at the zero equilibrium thereafter. Similar results were obtained for the case of using c 1 = 10, c 2 = 10, c 3 = 5 which is shown in Fig. 31. In this case, however, the deviation of the yaw angle decreases by 0.002 degrees to 0.006 degrees before returning to the zero equilibrium.
Overall, it is seen that using the delay model uses less controlling efforts to keep the yaw angle stabilised at about the zero equilibrium, while similar efforts must be exerted to bring the roll angle down to the zero equilibrium. These results imply that, even though the backstepping controller does a good job of handling the nonlinearities and modelling errors, the controller itself, like many such nonlinear controllers, is designed based on the principle of Lyapunov stability, which does not take into account the controlling effort penalties. In this respect the use of a better model would lessen the controller efforts required in the design of the nonlinear controllers.

Conclusion
This work has developed and extended a modelling methodology to characterise nonlinear damping and stiffness for second order systems from previous research. The methodology initially assumes a flexible model and, through the use of the integral-based parameter identification method, identifies a constant piecewise damping and stiffness parameters. The identified parameters are then correlated to the square of the velocity, effectively an energy-like function, to allow a construction of a nonlinear vibration model. To capture the effect of the time delay, a time delay model in the form of a second order delay differential equation was also considered. Theoretical analyses were also investigated into the global stabilities of both model families, namely the non-delay model and the delay model.
For both model families, coefficients were firstly identified across the entire data with equally spaced intervals. These coefficients were then correlated to the energy-like function to yield a nonlinear piecewise hyperbolic function, which effectively reveals the structural aspects of the vibration, namely that the energy levels also matter in the vibration generations. Both families of models were applied to a specialised robot arm system designed for orchard spraying. The nonlinear model gave significant improvements over the standard linear model in data fitting, which were further enhanced by the addition of the time delay consideration. This concept is different from the friction induced models normally seen in the literature, where initially a complicated structure must be assumed to model the system. To demonstrate the model in action with control, a backstepping controller was designed for both the non-delay model and the delay model. It was found that overall the use of the delay model expends less controlling efforts to keep both the yaw and roll angle about the zero equilibrium. Even though the backstepping controller, being a nonlinear controller, is very strong in handling the nonlinearities and modelling errors, the use of a better model lessens the required control efforts.