Application of Polynomial Transition Curves for Trajectory Planning on the Headlands

: This paper presents a method of polynomial transition curve application for making agricultural aggregate movement paths during headland turn drives as well as within the ﬁeld. Four types of agricultural aggregate paths in ﬁve di ﬀ erent variant designs are discussed. Each path is composed of only two curves, making the so-called transition bi-curve. The curvature described by the linear function as well as the third, ﬁfth, seventh, and ninth degree polynomials was designated. Moreover, a trajectory planning algorithm in which the movement proceeds along two transition curves composing the so-called bi-curve was proposed. The simulation was carried out applying the MATLAB program in which the 4th order Runge–Kutta method was used. The results were presented by means of ﬁgures showing the proposed paths and kinematic quantity courses in the displacement function. The obtained trajectories were compared regarding the size and kinematic quantities. The trajectories, whose curvature is described by the 3 ◦ polynomial, were found to possess the smallest absolute values of maximal acceleration and jerk and to lack jerk discontinuity. The proposed solutions can be applied for planning trajectory of not only agriculture machines and aggregates but also autonomous vehicles or mobile robots.


Introduction
Effective and ecological agriculture incorporates numerous projects consisting of modification of those already existing and implementation of new cultivation processes including essential agricultural works (sowing, fertilization, cultivation, and harvest) as well as rational activities such as protection of natural environment and wild animals [1,2]. The mechanical harvest directly connected with working and non-working movements of cultivable and subsidiary aggregates is one of the most energy consuming processes. Reduction of the distance and time of machine drive along the field results in a decrease of production costs and greenhouse gas emissions into the environment [3][4][5].
Most papers available in the literature study minimization of non-working distance on the headlands. The other criteria of the vehicle movement optimization are non-working time of headland turn maneuvers and labor costs. Bochtis and Vougioukas applied the algorithm allowing to choose the most favorable trajectory of working and non-working drives with the total non-working drive length reduction by 50% [6]. Bochtis et al. used the algorithmically computed sequence of field-work tracks (B-patterns). Their simulations showed the decrease in the non-working distance from 3.73% to 58.65% and increase in the area capacity up to 19.23% when compared with the conventional fieldwork patterns [7]. However, as a result of using an autonomous mechanization system (AMS) in the orchard, the length of the non-working run was reduced to 32.4% and the time was reduced to 40.2% [8]. The curvature of the planar curve is defined as: where: Δ -the angle between the tangents to the curve on the arc ends, Δl -the arc length. The (l)  function is defined as: (l) (l)dl     (2) whereas the coordinates of any point of the transition curve are calculated from: x(l) cos (l)dl   The first two conditions, based on which polynomial coefficients are determined for all analyzed curvatures, are as follows: For the linear curvature these conditions are sufficient for determination of integration constants and the κ(l) function form. For the curvature κ(l) defined by the degree of polynomial i ≥ 3 the other conditions assume the form: The curvature of the planar curve is defined as: (1) where: ∆θ-the angle between the tangents to the curve on the arc ends, ∆l-the arc length. The θ(l) function is defined as: whereas the coordinates of any point of the transition curve are calculated from: The curvature of the planar curve is defined as: where: Δ -the angle between the tangents to the curve on the arc ends, Δl -the arc length. The (l)  function is defined as: whereas the coordinates of any point of the transition curve are calculated from:  For the curvature κ(l) described by the i-th degree polynomial the equation assumes the form: The first two conditions, based on which polynomial coefficients are determined for all analyzed curvatures, are as follows: For the linear curvature these conditions are sufficient for determination of integration constants and the κ(l) function form. For the curvature κ(l) defined by the degree of polynomial i ≥ 3 the other conditions assume the form: For the curvature κ(l) described by the i-th degree polynomial the equation assumes the form: The first two conditions, based on which polynomial coefficients are determined for all analyzed curvatures, are as follows: For the linear curvature these conditions are sufficient for determination of integration constants and the κ(l) function form. For the curvature κ(l) defined by the degree of polynomial i ≥ 3 the other conditions assume the form: where: k = 1, 2, . . . , i−1 2 . The total number of conditions is always larger by one than the degree of the polynomial describing the curvature.
The paper is confined to the transition curves whose curvature is described by the linear function and the third, fifth, seventh, and ninth degree polynomials. Table 1 presents the polynomials coefficients in the simplified equation describing the function of curvature κ(l): a k l k (8) where: κ 1 -the curvature at the initial point, κ 2 -the curvature at the terminal point, l-the position of a given point along the curve.
Using formula (2) one can determine the function of tangent to the curve angle and record it also in the simplified form: θ(l) = a 0 + a 1 l + (κ 2 − κ 1 ) 10 k=2 a k l k , The coefficients of the above equation are listed in Table 2.
Formulas (8) and (9) as well as the values of coefficients listed in Tables 1 and 2 refer to the single transition curve of the length L. The angle formed by the tangent to the curve with the axis x at the initial point is designated θ b .

Algorithm of Trajectory Planning
The kinematic quantities were analyzed for four paths presented in Figure 3. All paths are composed of two transition curves BT and TE constituting the so-called transition bi-curve. The total length of the transition bi-curve is 2L. The curvature at the initial point of trajectory B is equal to κ = 0, at the tangency point T, half of the total path length being equal to the value of circle curvature κ = 1/R, Agriculture 2020, 10, 144 5 of 16 but at the terminal point E, similarly to the initial point it equals κ=0. The curvature at the point T is characterized by the largest value. It can be calculated using the minimal radius of the agricultural machine turn R min . point T is characterized by the largest value. It can be calculated using the minimal radius of the agricultural machine turn Rmin.
Individual paths differ in the angle α between the tangents in the initial and terminal points. The paths for which the angle α is π/2, π, 3π/2, and 2π were analyzed whereby it is possible to form a path for any angle α. In the case of the path called "half of lambda" (Figure 3a), the angle between the tangents is α = π/2 (angle of tangent to the curve at the point T is θT = π/4). For the path presented in Figure 3b called "half of chi" the angle between the tangents α = π but θT = π/2. The path in Figure  3c is called "rho". The angle between the tangents for this path is α = 3π/2, but θT = 3π/4. The path in Figure 3d is called "omega". In the case of this trajectory the angle between the tangents α = 2π but the tangent to the curve angle θT = π in the midlength of total motion trajectory. Due to the fact that velocities, accelerations, and jerks are determined as the first, second, and third displacement derivatives in relation to time, the simulation was made in the time domain. Assuming the constant maximal velocity of motion along the path vmax and introducing the substitution l=vmaxτ and dl=vmaxdτ into relation (3), rectangular coordinates of any point of the single transition curve in the function of time t are obtained from: The coordinates of centre S(xS, yS) and the radius of circle R in the global coordinates system as well as the angle α are the initial data.
The preliminary calculations:  the increase in the coordinates x and y of the single transition curve is calculated from relations (10) and (11) by inserting the adequate values of arguments Individual paths differ in the angle α between the tangents in the initial and terminal points. The paths for which the angle α is π/2, π, 3π/2, and 2π were analyzed whereby it is possible to form a path for any angle α. In the case of the path called "half of lambda" (Figure 3a), the angle between the tangents is α = π/2 (angle of tangent to the curve at the point T is θ T = π/4). For the path presented in Figure 3b called "half of chi" the angle between the tangents α = π but θ T = π/2. The path in Figure 3c is called "rho". The angle between the tangents for this path is α = 3π/2, but θ T = 3π/4. The path in Figure 3d is called "omega". In the case of this trajectory the angle between the tangents α = 2π but the tangent to the curve angle θ T = π in the midlength of total motion trajectory.
Due to the fact that velocities, accelerations, and jerks are determined as the first, second, and third displacement derivatives in relation to time, the simulation was made in the time domain. Assuming the constant maximal velocity of motion along the path v max and introducing the substitution l = v max τ and dl = v max dτ into relation (3), rectangular coordinates of any point of the single transition curve in the function of time t are obtained from: The coordinates of centre S(x S , y S ) and the radius of circle R in the global coordinates system as well as the angle α are the initial data.
The preliminary calculations: • the angle tangent to the curve θ T at the point T • the length of single transition curve L L = 2θ T R (13) • the time of motion along the single transition curve Agriculture 2020, 10, 144 6 of 16 • the increase in the coordinates x and y of the single transition curve is calculated from relations (10) and (11) by inserting the adequate values of arguments • the coordinates of the initial point B in the global system of coordinates The equations of the displacement in the global system of coordinates for the transition bi-curve are determined from the relations: As mentioned earlier, velocities, accelerations, and jerks are obtained as the first time, second time, and third time derivatives of displacement. The simulation followed the choice of values of turn radius and working velocity. The assumed value of turn radius is within that commonly used for agricultural tractors (Mc CorimckCX60L-R min = 6 m; New Holland TN65S-R min = 6.84 m; New Holland TN55D-R min = 7.26 m; Renault Palés 210-R min = 7.5 m; John Deere 5115M-R min = 9.2 m). Simulation was undertaken by means of the program MATLAB using the fixed time step (10 −2 s) and the 4th order Runge-Kutta (ode4 Runge-Kutta) method. There were applied the following input data: the coordinates of centre S-x S = 20 m, y S = 20 m, the constant maximal velocity of movement along the path-v max = π m/s, and the minimal turn radius-R = 7.5 m.

Results
The simulation results are presented in the form of courses of kinematic parameters with regard to displacement l. Figure 4 presents the "half of lambda" (HoL) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
Agriculture 2020, 10, x FOR PEER REVIEW 6 of 15  the coordinates of the initial point B in the global system of coordinates The equations of the displacement in the global system of coordinates for the transition bi-curve are determined from the relations: As mentioned earlier, velocities, accelerations, and jerks are obtained as the first time, second time, and third time derivatives of displacement. The simulation followed the choice of values of turn radius and working velocity. The assumed value of turn radius is within that commonly used for agricultural tractors (Mc CorimckCX60L-Rmin =6 m; New Holland TN65S-Rmin = 6.84 m; New Holland TN55D-Rmin = 7.26 m; Renault Palés 210-Rmin = 7.5 m; John Deere 5115M-Rmin = 9.2 m). Simulation was undertaken by means of the program MATLAB using the fixed time step (10 -2 s) and the 4th order Runge-Kutta (ode4 Runge-Kutta) method. There were applied the following input data: the coordinates of centre S-xS = 20 m, yS = 20 m, the constant maximal velocity of movement along the path-vmax = π m/s, and the minimal turn radius-R = 7.5 m.

Results
The simulation results are presented in the form of courses of kinematic parameters with regard to displacement l. Figure 4 presents the "half of lambda" (HoL) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials. Figure 4.The trajectories "half of lambda" (HoL) for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials. Figure 5 presents the courses of components x and y of HoL trajectory kinematic quantities for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials and Figure 6 shows the courses of resultant acceleration and jerk. The trajectories "half of lambda" (HoL) for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
Agriculture 2020, 10, 144 7 of 16 Figure 5 presents the courses of components x and y of HoL trajectory kinematic quantities for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials and Figure 6 shows the courses of resultant acceleration and jerk.    Figure 7 presents the "half of chi" (HoC) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials. Figure 7.The "half of chi" (HoC) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the HoC trajectory kinematic quantities for the curvatures described by the linear function and the third, fifth, seventh, and ninth degree polynomials are presented in Figure 8 and those of the resultant acceleration and jerk in Figure 9.    Figure 7 presents the "half of chi" (HoC) motion trajectories for the curvatures described by the inear function as well as the third, fifth, seventh, and ninth degree polynomials. Figure 7.The "half of chi" (HoC) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the HoC trajectory kinematic quantities for the curvatures escribed by the linear function and the third, fifth, seventh, and ninth degree polynomials are resented in Figure 8 and those of the resultant acceleration and jerk in Figure 9.
(a) (b) Figure 7. The "half of chi" (HoC) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the HoC trajectory kinematic quantities for the curvatures described by the linear function and the third, fifth, seventh, and ninth degree polynomials are presented in Figure 8 and those of the resultant acceleration and jerk in Figure 9.
The "rho" (RHO) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials are presented in Figure 10. Figure 7.The "half of chi" (HoC) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the HoC trajectory kinematic quantities for the curvatures described by the linear function and the third, fifth, seventh, and ninth degree polynomials are presented in Figure 8 and those of the resultant acceleration and jerk in Figure 9.   The "rho" (RHO) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials are presented in Figure 10. iculture 2020, 10, x FOR PEER REVIEW 10 o Figure 10.The "rho" (RHO) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the kinematic quantities of the RHO trajectory for t vatures described by the linear function as well as the third, fifth, seventh, and ninth deg ynomials are presented in Figure 11 but those of the resultant acceleration and jerk in Figure 1  The "rho" (RHO) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the kinematic quantities of the RHO trajectory for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials are presented in Figure 11 but those of the resultant acceleration and jerk in Figure 12. Figure 13 presents the "omega" (OMG) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of the x and y components of the OMG trajectory kinematic quantities for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials are given in Figure 14 but those of the resultant acceleration and jerk in Figure 15.
Analyzing the courses of acceleration and jerk it was found that the linear curvature causes discontinuity in the course of jerk components (Figures 5g,h, 8g, 11g,h and 14h). The exceptions are only the jerks along the axis y for the path HoC (8h) and those along the axis x for the path OMG (Figure 14g). On the other hand, considering the courses of both acceleration and jerk components, the largest absolute values of their local extremes were found using the curvature described by the polynomial 9 • . These values decrease with the decreasing degree of the polynomial describing the curvature. Thus, the trajectory, whose curvature is designated by the polynomial 3 • , is the most favorable taking the motion kinematics into account. Figure 10.The "rho" (RHO) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of x and y components of the kinematic quantities of the RHO trajectory for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials are presented in Figure 11 but those of the resultant acceleration and jerk in Figure 12.  (g) (h) Figure 11.The course of kinematic quantity components for the RHO trajectory: (a) displacement along the axis x; (b) displacement along the axis y; (c) velocity along the axis x; (d) velocity along the axis y; (e) acceleration along the axis x; (f) acceleration along the axis y; (g) jerk along the axis x; (h) jerk along the axis y.
(a) (b)  Figure 13 presents the "omega" (OMG) motion trajectories for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials. Figure 13.The "omega" (OMG) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
The courses of the x and y components of the OMG trajectory kinematic quantities for the curvatures described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials are given in Figure 14 but those of the resultant acceleration and jerk in Figure 15.  ure 13 presents the "omega" (OMG) motion trajectories for the curvatures described nction as well as the third, fifth, seventh, and ninth degree polynomials. ure 13.The "omega" (OMG) motion trajectories for the bi-curves whose curvature is described linear function as well as the third, fifth, seventh, and ninth degree polynomials. e courses of the x and y components of the OMG trajectory kinematic quantities res described by the linear function as well as the third, fifth, seventh, and ninth ials are given in Figure 14 but those of the resultant acceleration and jerk in Figure 1 Figure 13. The "omega" (OMG) motion trajectories for the bi-curves whose curvature is described by the linear function as well as the third, fifth, seventh, and ninth degree polynomials.
Application of the curvature described by the polynomials of increasingly higher degrees results in the step "clinching" of path shape which is particularly evident for the paths HoC, RHO, and OMG whereby the open path HoC becomes "narrower" and the closed ones RHO and OMG "tighter". The greatest differences were found for the trajectories whose curvature is described by the 1-5 degree polynomials and the smallest ones for the 7 and 9 degree polynomials (Figures 4, 7, 10 and 13). In the case of the curvature described by the polynomial 9 • one can find lengthening of the initial and terminal parts of the segment whose deviation from the straight line is insignificant (deviation from the rectilinear segment is 5 mm for about 7.5 m). This fact is of significant importance as it enables inclusion of these segments into the working area.  Analyzing the courses of acceleration and jerk it was found that the linear curvature causes discontinuity in the course of jerk components (Figures 5g,h, 8g, 11g,h, and 14h). The exceptions are only the jerks along the axis y for the path HoC ( Figure 8h) and those along the axis x for the path OMG (Figure 14g). On the other hand, considering the courses of both acceleration and jerk components, the largest absolute values of their local extremes were found using the curvature described by the polynomial 9°. These values decrease with the decreasing degree of the polynomial describing the curvature. Thus, the trajectory, whose curvature is designated by the polynomial 3°, is the most favorable taking the motion kinematics into account.
Application of the curvature described by the polynomials of increasingly higher degrees results in the step "clinching" of path shape which is particularly evident for the paths HoC, RHO, and OMG whereby the open path HoC becomes "narrower" and the closed ones RHO and OMG "tighter". The greatest differences were found for the trajectories whose curvature is described by the 1-5 degree polynomials and the smallest ones for the 7 and 9 degree polynomials (Figures 4, 7, 10, and 13). In the case of the curvature described by the polynomial 9° one can find lengthening of the initial and terminal parts of the segment whose deviation from the straight line is insignificant (deviation from the rectilinear segment is 5 mm for about 7.5 m). This fact is of significant importance as it enables inclusion of these segments into the working area.

Discussion
Four types of paths of agricultural aggregate movement in five different designs in Figures 6, 7, 10, and13 can constitute the chosen paths or fragments of paths during the working movements (HoL, RHO) or headland turns (HoC, RHO). On the other hand, the paths HoL and OMG can be used during the agricultural and auxiliary machine drives on the driveways for transport aggregates e.g., trans-shipping trailers, containers with spray fluids, etc.
A typical field turn at the angle 90° for both forward and backward drives, HoL (Figures 3 and  4) can be regarded as a fragment of the cultivation aggregate trajectory in the field and thus path curvature adjustment can ensure a smaller area in the field corner.
Depending on the sequence of drives in the field for the turn at the angle 180°, the path HoC (Figures 3b and 7), Ω-turn [19,[23][24][25], ∏-turn [24], T-turn [19,20,24], or the path HoL (Figures 3a and  4) with the addition of e.g., rectilinear fragments can be applied for obtaining a required distance between the parallel working drives in the field [19]. The effects of such trajectories planning are discussed in [6] where different designs of aggregate movement are presented taking into account the sequence of field drives. Based on the above, it can be stated that with the minimization of both working and non-working movement path numbers, determination of the sequence of trajectories and turns is of significant importance. Its value is directly affected by tractor turn radius, kind of agricultural work, as well as kind and width of cultivation aggregate, practical width of headland turns, field shape, and drive orientation-along or across the field.

Discussion
Four types of paths of agricultural aggregate movement in five different designs in Figures 6, 7, 10 and 13 can constitute the chosen paths or fragments of paths during the working movements (HoL, RHO) or headland turns (HoC, RHO). On the other hand, the paths HoL and OMG can be used during the agricultural and auxiliary machine drives on the driveways for transport aggregates e.g., trans-shipping trailers, containers with spray fluids, etc.
A typical field turn at the angle 90 • for both forward and backward drives, HoL (Figures 3 and 4) can be regarded as a fragment of the cultivation aggregate trajectory in the field and thus path curvature adjustment can ensure a smaller area in the field corner.
Depending on the sequence of drives in the field for the turn at the angle 180 • , the path HoC (Figures 3b and 7), Ω-turn [19,[23][24][25], -turn [24], T-turn [19,20,24], or the path HoL (Figures 3a and 4) with the addition of e.g., rectilinear fragments can be applied for obtaining a required distance between the parallel working drives in the field [19]. The effects of such trajectories planning are discussed in [6] where different designs of aggregate movement are presented taking into account the sequence of field drives. Based on the above, it can be stated that with the minimization of both working and non-working movement path numbers, determination of the sequence of trajectories and turns is of significant importance. Its value is directly affected by tractor turn radius, kind of agricultural work, as well as kind and width of cultivation aggregate, practical width of headland turns, field shape, and drive orientation-along or across the field.
The proposed RHO track (Figures 3c and 10) can be one of alternative trajectories enabling change of movement direction without stoppage of cultivation aggregate and making a maneuver not leaving a free area in the field corner. As for the drives in the perpendicular system, only four maneuvers of RHO connected with the field border can be made on the headland turns. Depending on the type of the agricultural work and aggregate width, the other movement paths can be applied also with the function of avoiding obstacles in the field area [26]. This results in the reduction of working distance and can be exploited e.g., in cereal harvest or grass mowing. The turn direction is essential and depends on the kind of agricultural work and the type of machine. The so-called turn agility of asymmetric aggregate in the form of tractor and mowing machines depends on determination of direction and minimal turn radius. Thus, the application of the proposed path RHO with the turn to the right can reduce the non-working distance significantly [27,28].
The loop maneuver OMG (Figures 3d and 13) is the trajectory of machine movement on the headland turns which can be used during the agricultural and auxiliary machines drives on the driveways for the transport aggregates e.g., trans-shipping trailers, containers with spray fluids, etc.
As follows from the simulation, the proposed movement trajectories enable smooth transition from the rectilinear path (κ = 0) into the curvilinear path (κ = 1/R) and the other way around without the step change of normal acceleration component. Taking into consideration discontinuity of kinematic quantities course and local extremes values, the trajectory whose curvature is described by the polynomial 3 • proves to be the most favorable. For the path, where the polynomial 9 • was used, the headland turns can be reduced by including the initial and terminal segments into the working area.
Author Contributions: M.B.: conceptualization, data curation, formal analysis, visualization, writing-original draft. P.K.: data curation, formal analysis, investigation, methodology, project administration, resources. K.G.: funding acquisition, project administration, supervision, validation, writing-original draft, writing-review & editing. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.