The present disclosure relates generally to automatic control of machines and, more particularly, to full coverage path planning over three-dimensional terrain.
Automation of machines, such as agricultural machines, can provide benefits that can improve crop yield and increase efficiency of crop processing. For example, automation of fertilizer spreaders can improve the speed and accuracy of fertilizer application. Automation of a machine requires data to control operation of the machine, see for example “Automatic guidance of a farm tractor relying on a single CP-DGPS” by Thuilot, B., Cariou, C., Martinet, P., and Berducat, M., in Autonomous Robots, 2002, no. 1, p. 53-71. Determining a path that a machine should follow for a particular operation requires planning that takes into account multiple factors such as machine capabilities. For example, a machine using an Ackermann steering mechanism with a constrained maximum steering angle can lead to unprocessed areas of a field due to geometric infeasibility of strict parallel paths. Machine and path planning limitations can also result in an undesirable amount of overlap of paths. What is needed is a method for path planning for a machine that takes into account multiple factors that affect efficient operation of the machine.
The present disclosure pertains to path planning for a machine to traverse an area (e.g., a field). In one embodiment, a spline trajectory based on a plurality of control points of a first path is calculated. A subset of the plurality of control points having an equal step is selected. The equal step can be based on a Euclidian metric (see, for example, “The Theory of Matrices”, 2nd edition, by Lancaster, P and Tismenetsky, M. Academic Press 1985). A direction of the normal to the spline trajectory for each of the selected points is determined. Control points within the subset that are a solution to a second order cone programming class optimization problem along each normal to the spline trajectory are searched for and the spline trajectory is extended to a border of the field to create a second path adjacent to the first path based on the control points. The optimization problem can minimize the average curvature at junction points of elementary sections of the spline trajectory and/or can minimize the average width overlap of adjacent paths. The optimization problem comprises an upper bound of a curvature constraint that can be based on a maximum steering angle of a machine to follow a planned path. The first path can be identified by an azimuth and a point in the area. The first path can be a section of an area boundary of the field. An apparatus and computer readable medium for path planning are also described.
In one embodiment, the Boustrophedon method is used to process field 102 (i.e., traversing field 102 in a back-and-forth manner). Field 102 is divided into parallel swaths each of which are further processed. A pre-selected section of a field boundary is used as the starting path (i.e., a seed path for use by an algorithm) to minimize the number of turning points. There is a large amount of literature on this subject. For example, see “On complete coverage path planning algorithms for non-holonomic mobile robots: Survey and Challenges” by Khan, I. Noreen, and Z. Habib in Journal of Information Science and Engineering, 2017, no. 1, pp. 101-121; “Region filling operations with random obstacle avoidance for mobile robots” by Z. L. Cao, Y. Huang, and E. L. Hall in Journal of Robotic systems, 1988, no. 2, pp. 87-102; “Driving angle and track sequence optimization for operational path planning using genetic algorithms” by Hameed, I. A., Bochtis, D., and Sorensen, C. in Applied Engineering in Agriculture, 2011, no. 6, pp. 1077-1086; “Coverage path planning algorithms for agricultural field machines,” by Oksanen, T., and Visala, A. in Journal of field robotics, 2009, no. 8, pp. 651-668; “Side-to-side 3d coverage path planning approach for agricultural robots to minimize skip/overlap areas between swaths” by Hameed, I. A., la Cour-Harbo, A., and Osen, O. L. in Robotics and Autonomous Systems, 2016, vol. 76, pp. 36-45; “Coverage path planning on three-dimensional terrain for arable farming” by Jin, J., and Tang, L. in Journal of field robotics, 2011, no. 3, pp. 424-440; “Automated generation of guidance lines for operational field planning” by Hameed, I., Bochtis, D., Sørensen, C., and Nøremark, M. in Biosystems engineering, 2010, no. 4, pp. 294-306; “Farmwork path planning for field coverage with minimum overlapping” by Fabre, S., Soures, P., Taix, M., and Cordesses, L. in 8th International Conference on Emerging Technologies and Factory Automation (ETFA), 2001, vol. 2, pp. 691-694; “Simulation study on coverage path planning of autonomous tasks in hilly farmland based on energy consumption model” by Shen, M., Wang, S., Wang, S., and Su, Y. in Mathematical Problems in Engineering, Hindawi, vol. 2020; and “Intelligent coverage path planning for agricultural robots and autonomous machines on three-dimensional terrain” by Hameed, I. A. in Journal of Intelligent and Robotic Systems, 2014, no. 3-4, pp. 965-983. Different heuristic algorithms are proposed to cover the field by equidistant paths, though there is a lack of strict problem settings and computationally efficient algorithms. The path is approximated using basis splines (“B-splines”), for which the first and second derivatives are continuous. The spline trajectory is represented by a set of elementary sections, each of which is constructed from four control points and lies in their convex hull.
The constrained maximum steering angle of machine 100 due to the Ackermann steering mechanism leads to a constraint on the maximum allowable curvature of the trajectory when planning paths. This can lead to unprocessed areas of the field due to geometric infeasibility of strict parallel paths.
In one embodiment, a path planning method consists of exclusion of unprocessed areas due to the constrained overlap of adjacent tracks (for example, 5-20% of the initial width). The exact value of the minimum overlap percentage can be determined iteratively. This will be the first value for which the constructing a parallel path results in the feasible solution). The method is reduced to solving an optimization problem of the Second Order Cone Programming (SOCP) class, in which the weighted sum of the average curvature at the junction points of the elementary sections of the splines and/or the average width of the overlap of adjacent tracks is minimized, while the values of the overlap and the curvature are constrained. As a result, due to the overlap, the paths are straightened, and the curvature constraint is satisfied.
Paths adjacent to the starting path are determined using an algorithm. In one embodiment, an adjacent path scheduling algorithm is as follows: A spline trajectory is constructed from control points of an original path, and then points are selected on the trajectory with an equal step in a Euclidean metric. For these points, the directions of the normal to the spline path are determined. Along the calculated normal directions, a search is made for points that are a solution to an optimization problem. The trajectories are extended to the border of a field (e.g., field 102). The obtained solution to the optimization problem is used as a set of control points for the next path. In one embodiment, the upper bound of the curvature constraint is taken into account when solving the optimization problem. Curvature constraints are expressed in the form of the second order cone constraints which casts the optimization problem into the SOCP class. This makes the problem numerically tractable as SOCP is convex, see “Lectures on Modern Convex Optimization: Analysis, Algorithms, and Engineering Applications” by Ben-Tal, A. and Nemirovski, A. in MOS-SIAM Series on Optimization, 2001, for additional information about convex optimization. There are freeware (for example, CVX, as described on the CVX Research website) and commercial packages (for example, GUROBI, as described on the Gurobi Optimization website) available for efficient solution of convex optimization problems.
In one embodiment, the terrain of field 102 is modeled as follows. A local coordinate system East North Up (“ENU”) is used, with the origin of coordinates being referenced to one of the points of field 102, the axis z (U) is directed upwards, and axes x and y correspond to east and north.
A digital elevation model (“DEM”) using the ENU coordinate system is defined by a set of control points in three-dimensional space, which are obtained by measurements, and the function z(x, y), which allows determination of the height of a point not contained in the measurement set by its coordinates in the X-Y plane. For construction of z(x, y) 2D B-spline interpolation is used. The spline coefficients are calculated once based on a set of control points, and then are used whenever determination of a height of a point needs to be determined.
To approximate the path, uniform cubic B-splines are used, which are first constructed in the X-Y plane, and then projected onto the surface using a function. Each path consists of elementary B-splines, and the first and second derivatives of the path are continuous in the spline joining currents. Each elementary B-spline is constructed using four control points ri−1, ri, ri+1, ri+2 ∈ R2 the following way:
where t is a spline parameter, 0≤t≤1. Thus, the target trajectory is given by a set of control points r1, r2, . . . , rn, where each four consecutive points defines an elementary segment of the trajectory. The set of control points of the trajectory must be supplemented with additional points r−1, r0, rn+1, rn+2, for example:
r
−1=3r1−2r2,r0=2r1−r2,rn+1=2rn−rn−1,rn+2=3rn−2rn−1
In this case, the set of control points is complemented so that the last and first four points lie on a straight line. The start and end of the spline path will match the first and last original control points (r1 and rn). This is due to the fact that the convex hull of the first and last four control points degenerates into a segment.
Calculation of the derivative with respect to the spline parameter determines the direction of the tangent to it
Unlike natural parametrization for B-splines, the length of the first derivative vector may differ from one, but the direction also coincides with the tangent. Accordingly, the direction of the normal to the projection of the trajectory (in the X-Y plane) is determined by the rotation of the tangent vector r(i)′(t) to ±π/2.
Second derivative with respect to the parameter is:
r
(i)
(t)=b0n(t)ri−1+b1n(t)ri+b2n(t)ri+1b3n(t)ri+2
b
0
n(t)=1−t,b1n(t)=3t−2,
b
2
n(t)=1−3t,b3n(t)=t.
Due to the spline being parameterized not by the length, but by the spline parameter t, 0≤t≤1 the angle between vectors r(i)′(t) and r(i)
Curvature of the path at a point r(i)(t) ∈R2 is defined as
using a cross product, the norm of vectors is Euclidean, φ(t)−acute (or right) angle between vectors r(i)′(t) and r(i)
Since the target vehicle has the Ackermann steering mechanism, the constrained angle of steering wheels, assumes the constraint on the curvature of the path: k(i)(t)≤kmax∀i=1, . . . ,n.
In one embodiment, a starting path in the form of a straight line can be determined by a point in field 102 and an azimuth. In this embodiment, further adjacent paths are constructed on the left and right sides of the initial one. In order to optimize the number of turns, a section of the field boundary can be used as the starting path. In this scenario, only one adjacent path is planned further, by only one side of the initial one.
In one embodiment, the parallel path planning problem is formulated as follows. Given some initial path, at approximately equidistant points are set on it r1, r2, . . . , rn ∈R3 (e.g., the path was resampled). These points are projections of the spline path in the X-Y plane onto an uneven surface. In one embodiment, it is required in three-dimensional space to find the control points of a neighboring path for which the curvature constraint would be satisfied. In this case, the width of the track cannot exceed the specified distance
N
i
z
=z(xi+Nix,yi+Niy)−z(xi,yi).
In one embodiment, the three-dimensional vectors Ni=(Nix,Niy,Niz) will also be considered. The following simplifying assumption can be taken: the piece of surface in the ∥Ni∥− neighborhood around the current point is sufficiently precise treated as a plane. Hereinafter, it will be assumed that all norms of vectors are Euclidean. The control point of the new spline will be ri+diNi, where di>0. The swath width constraints are:
α
Referring to
Referring to
In one embodiment, the optimization problem of finding a neighboring path is determined as follows. For simplicity, constraints on the curvature of a flat trajectory will be applied, i.e. for B-splines in the X-Y plane. Proof for the following statement is provided below.
Statement 1. The maximum value of the norm of the second derivative vector with respect to the spline parameter ∥r(i)
r
(i)
(t)=(1−t)ri−1+(3t−2)ri+(1−3t)ri+1+(t)ri+2=(ri−1−2ri+ri+1)+t(−ri−1−2ri+ri+1+ri+2)
Dependency ∥r(i)
Curvature of the elementary segment is
The curvature constraint must be satisfied at all points of the new path. Its upper bound {circumflex over (k)}(i)(t)≥k(i)(t) is
Further, assume that ∥r(i)′(t)∥ is slowly varying with an approximately equidistant location of control points. This allows ∥r(i)′(t)∥ to be consider approximately constant on each elementary spline. Based on this, the constraint on the estimate of the curvature is required only at the connection points of the elementary sections of the splines:
For uniformly distributed control points ds=∥ri−ri+1∥=∥ri−ri+2∥
∥ri+1−ri−1∥≤2ds
∥r(i)′(0)∥2=∥−0.5ri−1+0.5ri+1∥≤ds
The condition below follows from the inequality:
In one embodiment, if in a problem whose objective function minimizes Σiki2 includes conditions:
then straightening of trajectories is enforced, i.e. ∥r(i)
In one embodiment, second order cone programming (SOCP) path planning with path straightening is formulated as follows:
i=1, 2, . . . ,n. where β is a weighting coefficient that determines the rate of straightening of trajectories when passing from row to row. If β=0 no straightening is assumed, only the averaged overlap is intended to be minimized.
In one embodiment, a more accurate estimate can be calculated by accounting for the actual value ∥r(i)′(0)∥2 Let di=1−qi, i=1, 2, . . . ,n.
All terms but (Ni−1xqi−1−Ni+1xqi+1)2 and (Ni−1yqi−1−Ni+1yqi+1)2 are linearly dependent on their variables and allow for formulation of SOCP.
The methods, calculations, and operations described herein can be performed on a computer. A high-level block diagram of such a computer is illustrated in
The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the inventive concept disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the inventive concept and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the inventive concept. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the inventive concept.
Filing Document | Filing Date | Country | Kind |
---|---|---|---|
PCT/RU2021/000410 | 9/23/2021 | WO |