Smooth Curve Fitting of Mobile Robot Trajectories Using Differential Evolution

Mobile robots have recently attracted the attention and applicability in field areas ubiquitously. Within the context of autonomous navigation, path planning is relevant for comfortability, safety, execution time and energy savings. In this paper, we propose an approach to suggest smooth paths from observed robot trajectories by optimizing fitting and smoothness criteria using Differential Evolution with distinct modes of initialization, selection pressure, exploration and exploitation. Our rigorous computational experiments using a relevant set of real-world robot trajectories from the Boe-Bot mobile robot architecture show the feasibility and efficiency of our approach in computing smooth curves, suggesting the superior performance of the greedy initialization scheme based on the triangular convex hull of the robot trajectory, and Differential Evolution with exploitative and parameter adaptation schemes such as Rank-Based Differential Evolution (RBDE), Adaptive Differential Evolution with External Archive (JADE) and Strategy Adaptation Differential Evolution (SADE). Our obtained results offer the building blocks to further advance towards developing data-driven curve fitting and path planning algorithms, which may find use in several real-world applications in Robotics and Operations Research.


I. INTRODUCTION
In recent years, mobile robots have extended their application scope in field areas, including agriculture, forestry and manufacturing. In line of the above, to realize efficient navigation, path planning with smoothness considerations is relevant for comfortability, safety, minimum execution time, and energy savings.
Path planning has been extensively studied in the literature. Often, in the context of autonomous robot navigation, collision-free trajectories are generated considering safety and optimality. In line of these achievements, well-known methods such as Rapidly-exploring Random Tree (RRT) and Probabilistic Roadmaps (PRM) guarantee probabilistic completeness, while RRT* guarantees asymptotic optimality. However, challenges in scalability often arise, since approaches based on RRT and PRM have to extend over large maps to find optimal trajectories. To avoid covering the entire map, often approaches based on optimization using heuristic sampling and gradient of the cost function (length) are used.
The associate editor coordinating the review of this manuscript and approving it for publication was Ruofei Ma .
Examples of these approaches include CHOMP [1], STOMP [2], and TrajOpt [3]. However, these methods are sensitive to initial conditions (initial trajectory). Also, according to the presence or absence of knowledge of the environment, it is well-known the division between global and local methods in path planning. Also, the incremental methods such as Rapidly Exploring Random Trees (RRT), and their extensions, are well-known in the literature for its exploratory features. Furthermore, path planning algorithms using graph-theory and search heuristics are also well-known, e.g. Dijkstra and A*, visibility graphs combined with sampling heuristics [4]- [8].
Also, researchers have tackled the path planning problem by using geometric and graph-based approaches [9]- [13]. In line of these schemes, it has been shown that path planning in triangulated space is highly accurate and efficient. For a map having polygonal obstacles with n vertices, it is possible to use the Delaunay Triangulation of the free space to render a connected graph with O(n) nodes, each of which represents the triangles conforming the free-space. Then, path planning can be efficiently achieved by graph search methods on the adjacency graph of the triangulation. If one uses the Dijkstrabased search [9], time complexity is bounded by O(n.log(n)). VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Complementarily, it is possible to use A* search [10], [11] along with the funneling algorithm [12], [13] to achieve a quasi-linear time complexity.
On the other hand, in order to consider safety and comfortability in mobile robot driving, the inclusion of smoothness factors in path planning has attracted the attention of the community [4], [14], [15]. In terms of safety, path planning methods usually maximize the distance between the vehicle and obstacles [16], and in terms of comfortability, path planning methods usually consider the change and continuity of curvature of the path. For instance, in [17] paths using Bézier curves are smoothed, in which theta* RRT is able to render feasible any-angle paths by using geometric information of the workspace and satisfying the non-holonomic constraints. Also, in [18], the Dijkstra algorithm is used to compute shortest paths, and then smooth trajectories are generated by cubic Bézier curves. In [19], clothoid curves are used to achieve paths which are shorter and have lower curvature derivatives. In [20], various types of Bézier curves are used to generate smooth road trajectories. Furthermore, in [21], vehicle and road constraints are considered in order to generate collisionfree paths. In [14], Support Vector Machines are used to compute the critical points in obstacle-prone environments, and the control points of Bézier curves are optimized to compute smooth paths. In [22], Bézier curves were smoothed by Genetic Algorithm with diversity factors. In [23], the latticebased motion planners were proposed to compute feasible paths, in which the ACADO Toolkit has the role of path smoothing and a Model Predictive Controller (MPC) has the role of optimizing control outputs to follow the trajectory. In [24], the global path planning was developed based on Visibility Graphs and Bézier curves were used to minimize the maximum curvature in order to generate a smooth path.
In line with the above, in [25], the path planning for corridors with curvature discontinuities is proposed by composing half-S-shaped paths of the Four Parameter Logistic (4PL) family achieving zero-end curvature and collision-avoidance with passages. In [26], quintic Bézier curves compute the smooth sequence of curves of A*-based paths by solving the constrained problem and globally-convergent movingasymptotes (MMA) algorithm. In [27], the generation of smooth paths by gradient-informed post-smoothing algorithm (GRIPS) is proposed, in which locally placed vertices are optimized by considering the safe distance between the robot and the obstacles.
In order to generate smooth trajectories, the conventional methods basically have focused on the knot placement of polynomial curves to compute smooth paths and then perform trajectory following to ensure compliance with the generated path. In this paper, we propose an approach to generate smooth trajectories from given robot trajectory points. This problem is related to the well-known curve fitting and reconstruction problem, that is fitting data points to curves, which has attracted the attention of the non-linear optimization community due to the gradient-free and ability to deal with multimodal optimization approaches, e.g. Particle Swarm Optimization [28], Artificial Immune Systems [29], Lasso [30], Estimation of Distribution Algorithms [31], Invasive weed optimization [32], Differential Evolution [33], and Hybrids [34].
Although the above-mentioned schemes have rendered competitive curve fitting metaheuristics, it is unclear whether initialization, selection pressure, parameter adaptation, exploration or exploitation play key roles for the competitive performance to fit and compute fair curves to robot trajectory data. Thus, in this paper, to tackle the abovementioned line of inquiry, our focus and contribution are as follows: • We propose an approach to generate smooth paths from observed robot trajectories by optimizing the fitting and the smoothness performance with Differential Evolution with distinct modes of initialization, parameter adaptation, selection pressure, exploration and exploitation.
In particular we used the following five relevant classes: -DERAND: DE/rand/1/bin strategy -RBDE: Rank-Based Differential Evolution -JADE: Adaptive Differential Evolution with External Archive -SADE: Strategy Adaptation Differential Evolution -DESIM: Differential Evolution with Similarity Based Mutation • The exhaustive experiments using all feasible scenarios of torque control of Boe-Bot mobile robot show (1) the feasibility and efficiency of our approach to generate smooth curves from given trajectory data, (2) the improvement of the convergence rate by a greedy initialization scheme considering the triangular convex hull of trajectory data, and (3) the competitive convergence and curvature profiles being attainable by the class of Differential Evolution algorithms with exploitative and parameter adaptation schemes: RBDE, JADE and SADE.
• Our proposed approach offers insights on the performance of relevant Differential Evolution algorithms generating smooth paths for mobile robots complying with fitting and smoothness criteria, thus being potential to enable comfortability, energy-efficiency, and computationally efficient controllers. In the following sections, we describe the basic idea of our proposed approach, and discuss our computational experiments.

II. PROPOSED APPROACH
Basically, we tackle the problem of fitting and computing fair curves to prescribed robot trajectory data. Smooth curves are computed by minimizing a cost function comprising the fitting error to trajectory inputs and the smoothness defined by the curvature profile of a curve. To show an example of our proposed scheme, Fig. 1 shows the main elements involved in our approach. Here, the curve r(u) is defined by the control points p i , i = 0, 1, . . . , n, each of which is computed by minimizing the fitting error and the smoothness component of the curve. The fitting error aims at minimizing the shortest distances from prescribed input coordinates q j , j = 1, 2, . . . , m defined by the robot trajectory to the curve r(u), as shown by Fig. 1-(b). The smoothness of the curve is defined by aiming at finding the monotone variation in the curvature profile of r(u), as portrayed by the green-colored region of Fig. 1 showing the magnitude of the radius of curvature along the curve r(u) for u ∈ [0, 1].
In the following section, we briefly describe the key components and dynamics of our proposed approach.

A. BASIC CONCEPT
Let a Bézier curve of degree n be defined by where p i is the i-th control point of the Bézier curve in the domain R 2 , and B n i (u) is the i-th Berstein polynomial of degree n for i ∈ {0, 1, 2, . . . , n}. For simplicity, and without loss of generality, we use the degree n = 3, in which Bernstein polynomials are represented by Basically, the curve r(u) lies within the convex hull of the control points p i , i = 0, 1, . . . , n. Also, each point lying in r(u) is computed by the weighted-average of the control points p i , with weights defined by Bernstein polynomials, and fulfilling the condition n i=0 B n i (u) = 1. Then, smooth paths are computed by minimizing the following cost function min where p i is the i-th control point of the Bézier curve r(u), depicted by nodes connected by lines in pink color on Fig. 1-(a), E is the component related to the fitting error of the curve r(u) to given input trajectory coordinates q j , j = 1, 2, . . . , m, H is the smoothness component (fairness metric) of the curve r(u), and λ is the user-defined parameter to balance the fitting error and the smoothness factors. Concretely speaking, the following relation computes the fitting error: where • m is the number of input trajectory points, • a j is the vector pointing from c j to q j , • q j ∈ R 2 is the j-th input trajectory coordinate obtained from the robot's trajectory measurements, • c j ∈ R 2 is the nearest point in the curve r(u) to the input coordinate q j , which is obtained by the following: where u ∈ [0, 1], and · denotes the Euclidean norm. Furthermore, the smoothness component of the curve is computed by where κ is the curvature of the curve r(u) and s is the arclength of the curve r(u).
In the above expression By the chain rule, r ×r (r · r +r ·r )+ r × r (r × r ) To give a glimpse of the smoothness factor, Fig. 1 also shows an example of the curvature profile of a Bézier curve as well as the main components to compute the error in curve fitting. Here, the radius of curvature ρ = 1/κ is depicted by green color in Fig. 1. The reader may note that the errors in curve fitting E as well as the smoothness component H (fairness metric) are contradicting terms, determining the optimal ratio λ is out of the scope of this paper.
Furthermore, we use the integration of d 2 κ ds 2 2 due to its relevance to compute curvature profiles with linear nature. The reader may note that using other fairness metrics, such as the minimum energy curvature (MEC) and the minimum variation curvature (MVC), is also possible. Also, once curves minimizing the criterion in Eq. 8 are computed, it becomes possible to generate compounded smooth paths with C 0 continuity by joining the curvature profiles of Bézier paths and by matching the tangent angles α and β between two consecutive curvature profiles. However, compounded paths will require further pruning to satisfy optimality of the cost function F (and smoothness component H ) at the connecting node points of two path elements. The study of compounded fair curves is out of the scope of this paper, and left in our future agenda.

B. DIFFERENTIAL EVOLUTION
The cost function depicted by Eq. 8 is of non-linear landscape, and computing its gradient is of non-trivial nature. Thus, we use we use the class of gradient-free optimization algorithms considering features of balance between exploration and exploitation. In this paper, we use the relevant classes of Differential Evolution (DE) [35] algorithms due to the flexibility and the versatility to realize diversity of exploration and exploitation during search. In line of the above, we used five relevant classes of Differential Evolution as optimization heuristics, each of which denotes distinct modes of selection pressure and balance between exploration and exploitation. Considering other nature-inspired heuristics is straightforward, yet including a large number of such heuristics is out of the scope of this paper.
Solutions in a minimization problem are sampled by: where t denotes the iteration index, • is the element-wise Hadamard product, v t is the mutant vector at iteration t, b t is a D-dimensional binary vector, r t,j is a random number with uniform distribution U [0, 1], krand is a random integer uniformly distributed in U [1, D], and CR is the cross-over probability. Due to handling Bézier curves of degree n = 3 and due to control points laying in the plane (p i ∈ R 2 , i ∈ {0, 1, 2, 3} in the formulation of Eq. 1 and Eq. 8), the (individual) vector x t encodes the x − y coordinates of the control points of the Bézier curve, thus D = 2(n + 1) = 8. Generally speaking, Differential Evolution considers the evolution of the set P of individuals, with NP = |P| denoting the population size, in which three relevant processes are performed: • initialization to generate the population P, • mutation to generate the vector v t , • crossover to generate the trial vector u t (Eq. 17), • selection to generate the vector x t+1 (Eq. 16). As for initialization of the population, we used two different procedures, as follows: • Random Initialization. Control points p i ∈ R 2 are set arbitrarily within the convex hull defined by the x − y coordinates of − q m and q m , for a constant ≥ 1. This procedure is in line with the conventional initialization schemes in Differential Evolution, in which solutions are sampled from the hypercube defined by the lower and upper bound of the search space. Since the coordinate q m represents the last input coordinate from the observed robot trajectory, new coordinates sampled within − q m and q m ensure the feasibility of the search space.
• Greedy Initialization. Here, the first (last) control point is set equal to the first (last) input trajectory; and p 2 and p 3 are computed from interpolating the vertices of the triangular convex hull of the coordinates q j (j ∈ [m]), as portrayed by Fig. 2. Concretely speaking, this procedure is realized as follows: where q max is the farthest coordinate from the segment p 0 p 3 , ≥ 1 is a user-defined constant to estimate the triangular convex hull of the input coordinates q 1 , q 2 , . . . , q m , and d 1 , d 2 ∈ [0, 1] are interpolation constants. This procedure enables to define the control points of the Bézier curve such that the coordinates defined by the robot trajectory lie within the convex hull of the polygon defined by p 0 , p 1 , p 2 and p 3 . After initialization, procedures related to mutation, crossover and selection are performed until the number of maximum number of evaluations of the cost function is attained. In line with our goal of evaluating distinct modes of selection pressure, exploration and exploitation, we describe the relevant heuristics used in our study in the following sections.

1) DE/rand/1/bin (DERAND)
This algorithm is one of the well-known strategies able to handle nonlinear cost functions [35], offering exploration abilities during search. Here, the mutant vector v t are sampled as follows: where r 1 , r 2 , r 3 ∈ U [1, NP], with r 1 = r 2 = r 3 = k, NP is the population size, and F d is the scaling factor. In the above, exploration is tunable by the factor F d .

2) RANK-BASED DIFFERENTIAL EVOLUTION (RBDE)
Here, the mutant vector v t is computed by [36]: where x 1 t , x 2 t and x 3 t are individual vectors selected from the population by using the Whitley Distribution. That is, for i w = {1, 2, 3}: where P S c is the c-th individual in the population sorted by fitness from best to worst, θ is a user-defined bias term, r is a random number uniformly distributed in U [0, 1], and · denotes the floor function.
The above-mentioned sampling scheme focuses on the selective pressure to allow the improved exploitation by using a rank-based ordering of the population. The range of θ ∈ (1, 3] is found to be favorable to improve the convergence in non-separable problems [36].

3) ADAPTIVE DIFFERENTIAL EVOLUTION WITH EXTERNAL ARCHIVE (JADE)
This algorithm uses the successful historical references to update the learning parameters CR and F a [37]. The mutation is based on DE/current-to-pbest/1 with an archive of inferior solutions, as follows: where F x t is the mutation scaling factor associated to vector x t , x pbest is a random individual from the 100p% best of the population P for a constantp ∈ (0, 1], r 1 is an (integer) index chosen randomly from [1, NP],x r 2 t is a random individual from P ∪ A (r 2 ∈ [1, NP + |A|], r 1 = r 2 = k), in which: • A denotes the set of inferior solutions, • A is empty during initialization, • The vector x t is added to the set A if the selection of u t succeeds in (Eq. 16), and • Elements of the set A are deleted arbitrarily if |A| > NP. Also, the scaling factor and crossover rate are updated by the following rules: where N(·, ·) is a random number with normal distribution, Cauchy is a random number with Cauchy distribution (truncated at [0, 1]), and CR x t is the crossover rate associated to vector x t . During initialization, µ CR = µ F = 0.5; σ CR = σ F = 0.1. Then, after selection (Eq. 16), the following is computed: where c ∈ [0, 1] is a user-defined constant for averaging, S CR and S F are sets of successful parameters corresponding to CR x t and F x t , respectively, S CR is the arithmetic mean of S CR , and S F is the Lehmer mean of S F as follows:

4) STRATEGY ADAPTATION DIFFERENTIAL EVOLUTION (SADE)
SADE allows a pool of strategies learn selection probabilities over a number of generations, wherein probabilities are updated in line with successful mutation [38]. Thus, rather than implementing one mutation strategy, SADE uses a set of strategies. Here, v g,t denotes the vector generated by the gth strategy at iteration t. We use the original pool of strategies as those used in [38], that is the DE/rand/1/bin strategy defined by Eq. 25, and four additional strategies, as follows: DE/current-to-rand/1/bin (with CR = 1) where r 1 , r 2 , r 3 , r 4 , r 5 are mutually exclusive random integers in the uniform range [1,NP]. Furthermore, the probability to select strategies, the crossover probability CR and the scaling factor F d are computed adaptively following the same learning rules proposed in SADE [38]. VOLUME 8, 2020

5) DIFFERENTIAL EVOLUTION WITH SIMILARITY BASED MUTATION (DESIM)
DESIM considers the similarity-based mutation strategy, and was shown to outperform explorative and exploitative variants of DE [39]. The mutation is as follows: where x si is the si-th vector from the population sorted by similarity in descending order with respect to the best (fittest) individual. Here, similarity is computed in terms of the Euclidean distance. Also, in the above, si is an integer uniformly sampled from U [si l , si u ], is the maximum number of function evaluations, ω is the number of function evaluations at iteration t, δ is the intensification parameter (δ = {5, 10, 15} are found to be favorable [39]). Furthermore, the parameters F d and CR are computed by the adaptation principles from JADE. Initialization of the population is realized by the Opposition Based Learning.

III. COMPUTATIONAL EXPERIMENTS
In order to evaluate the performance and the feasibility of our proposed approach, we performed a set of experiments in which our goal is to evaluate the convergence ability when minimizing the cost function, and the degree of the smoothness and curvature profiles of the obtained paths.

A. SETTINGS
The robot trajectories are realized by the Boe-Bot hardware [40] with a marker-based localization through an RGB camera @30FPS located on top. Our algorithms were implemented in Matlab 2018a, in which our computing environment consisted of i7-4930K @3.40GHz. The evaluated algorithms are labeled as follows: • DERAND: DE/rand/1/bin strategy • RBDE: Rank-Based Differential Evolution • JADE: Adaptive Differential Evolution with External Archive • SADE: Strategy Adaptation Differential Evolution • DESIM: Differential Evolution with Similarity Based Mutation Our main rationale in using the above described heuristics is due to our aim in evaluating whether initialization, selection pressure, parameter adaptation, exploration and exploitation play key roles in the competitive performance to fit and fair curves to robot trajectory data. Basically, our approach extends the conventional nonlinear least squares curve fitting problem to ordered data [41] by allowing the inclusion of a functional fairness to obtain monotonically varying robot trajectories. Although it is possible to extend the local algorithms for curve fairing [42], the conjugate gradient-based methods [43], and the nonlinear constrained optimization approaches such as Sequential Quadratic Programming [44], their performance is limited by the ability to approximate a local minimum, the differentiability of the cost function, the sequential nature of sampling solutions and the computational cost involved in computing gradients with respect to all degrees of freedom of the curve. Thus, our focus of interest span metaheuristics due to the nonlinearity of the cost function [41], [43], the inherent parallelization scheme to sample solutions through a population-based approach, and the potential to scape local optima, thus enabling to compute fit and fair curves in the global perspective. The evaluation of gradient-based and local heuristics is out of the scope of this paper, whose study and integration to metaheuristic approaches is left for future work in our agenda.
The key motivations of using the above parameters are as follows: • Crossover probability with CR = 0.5 enables to consider equal importance to historical search directions.
• Small population size NP = 10 and number of evaluations up to 2000 and 10000 are used in order to evaluate the (competitive) performance of the gradientfree classes of Differential Evolution algorithms under tight evaluation budgets. This setting allows to evaluate the feasibility on the fast convergence of Differential Evolution during distinct modes of selection pressure.
• Also, as for initialization of the population, in order to allow the relevant initialization schemes with reasonable geometry of the search space, we use η = 2 and interpolation constants d 1 = d 2 = 2/3. The optimal selection of initialization parameters is out of the scope of the paper. The parameter for smoothness preference in Eq. 8 was set at λ = 0.01, which showed the reasonable results, after a number of trials, in balancing the fitting error and the smoothness of curves. Fine-tuning of the above-described parameters is out of the scope of this paper.
Furthermore, the mobile robot (Boe-Bot hardware [40]) is able to generate trajectories as byproduct of navigation. Basically, by fixing the ratio of Pulse Width Modulation (PWM) on the servo motors, the expected trajectory is circumference; yet due to inherent noises in current, the mechanical configuration of the gearboxes and the interactions between the wheel and the floor, the realized trajectories are rather noisy. Thus, in order to collect real-world robot trajectories, we employed the following setup: • Boe-Bot [40] uses discretized combinations of torque through Pulse-Width Modulation ratio of the servo motors of the robot.
• Pulse-Width Modulations are discretized by using a microcontroller (Arduino UNO) attached to the robot hardware.
• Trajectories are tracked by an RGB camera located on the ceiling, and markers located on top of the robot. For simplicity and without loss of generality, we collected trajectories by using feasible values of Pulse-Width Modulations in one side of the servo motor. In order to show the kind of collected trajectories, Fig. 3 shows the set of coordinates and trajectories for all discretized Pulse-Width Modulations of the torque in the left side of the servo motor (inducing a curved motion to the right). Overall, 28 trajectories were collected, in which the origin was set at locations close to (0, 0) in Fig. 3, and the end locations are positioned at the right side of the plot. The reader may note that the accurate positioning of the robot at (0, 0) is irrelevant since curve modeling considers the relative coordinates of the point q 1 (initial point in the curve r(u). We observed that the number of points m (Eq. 9) were variable in each trajectory, and the following range was obtained m ∈ [311, 766], which is mainly due to the inherent nature of the noisy and curved paths.

B. CONVERGENCE PERFORMANCE
To evaluate the convergence ability of our proposed approach, Fig. 4 and Fig. 5 show the convergence of the minimization of the cost function (Eq. 8) under = 2000 and = 10000 function evaluations, respectively. In these figures, the x-axis denotes the number of function evaluations, and the y-axis denotes the value of the cost function F. Both Fig. 4 and Fig. 5 report the mean of the convergence of each of the 28 trajectories cases over 20 independent runs.
• Furthermore, in order to portray the effect of the greedy initialization on the convergence performance within the first = 2000 function evaluations, Fig. 4 shows the convergence of the cost function when minimizing Eq. 8 for 28 trajectories. Here, greedy initialization is depicted by blue color, whereas the convergence under random initialization is portrayed in red (DERAND), black (RBDE), green (SADE), purple (JADE) and orange (DESIM) colors. • Also, in order to show the convergence behaviour under large number of function evaluations, Fig. 5 shows the convergence under = 10000 function evaluations. VOLUME 8, 2020 By observing Fig. 4 and Fig. 5, we note the following facts: • The number of iterations required for convergence is around 2000 in the best scenarios, in which RBDE, JADE and SADE show the best convergence speed.
• Reasonable convergence fulfilling the termination criteria occurs in all studied cases when minimizing the FIGURE 6. Test to evaluate the statistical significance of using the greedy initialization scheme. Wilcoxon tests were performed at 5% significance level over 20 independent runs and overall robot trajectories. cost function over all collected trajectories and classes of Differential Evolution algorithms.
• As shown by Fig. 4, the improved convergence performance is observed within the stages of the search process when the population uses the greedy initialization scheme (Eq. 20 -Eq. 24). The above-mentioned observations show the amenability of the Differential Evolution classes to find solutions within the basins of the cost function F, suggesting the ability to compute smooth curve configurations within small number of function evaluations.
In order to evaluate the statistical significance of the improvement of the greedy initialization scheme, Fig. 6 shows the Wilcoxon rank test at 5% significance level over the converged cost function F under = 2000 function evaluations. Here, the x-axis shows the robot trajectory scenario, and the y-axis shows the Differential Evolution class. As shown by Fig. 6, dark blue colors denote that the greedy initialization scheme (Eq. 20 -Eq. 24) is significantly better than the random initialization at 5% significance level. Also, the greedy initialization improves the performance significantly in RBDE (21 out of 28 cases), SADE (23 out of 28 cases) and JADE (26 out of 28 cases). However, no significant improvement was observed in DERAND and DESIM. We believe that the greedy initialization is advantageous for better convergence in RBDE, SADE and JADE due to the exploitation abilities of these schemes compared to DERAND and DESIM. Likewise, the exploration mechanism of DERAND and DESIM in early generations enables to counter-balance the lack of greedy solutions.
Furthermore, to show the performance comparison among the Differential Evolution algorithms, Fig. 7 shows the statistical comparison of the converged cost function after = 2000 function evaluations. Here, Wilcoxon tests were performed at the 5% significance level over 20 independent runs and overall robot trajectories. By observing the results in Fig. 7 we observe the following facts: • All studied Differential Evolution algorithms achieve similar performance when using random initialization.
• SADE was shown to outperform DESIM in 21% of the trajectory cases, yet found to achieve similar results in 79% of the trajectory cases.
• In greedy initialization, RBDE, SADE and JADE attained the best results, outperforming DERAND and DESIM. Also, DERAND and DESIM were found to attain similar results in 86% of the trajectory cases.
• When using greedy initialization, SADE and JADE achieved similar results in all trajectory cases. However, JADE was shown to outperform RBDE and SADE in 64% and 50% of the trajectory cases, respectively.
The above-mentioned facts confirm the feasibility and the efficiency for competitive convergence. Also, the above results show the advantageous properties of the greedy initialization to improve the convergence of Differential Evolution algorithms with exploitative and parameter adaptation mechanism such as RBDE, SADE and JADE.

C. CURVATURE PROFILES
In order to show the kind of smooth curves which our proposed approach is able to compute, Fig. 8 shows examples of curves and curvature profiles rendered as a result of the minimization problem in Eq. 8. In these figures, the following elements are portrayed: • A selected number of points in the robot trajectory are denoted by nodes colored in red.
• The rendered smooth curves approximating the robot trajectories are shown in dark color.
• The control points p i of the curve r(u) are denoted by vertices and lines colored in purple.
• The radius of curvature are shown as curved shapes in blue, green and orange color, respectively. By observing Fig. 8, we can note that smooth curves reasonably fit the collected robot trajectories. Also, in order to show the properties of the curvature profiles and to portray the (variation of) curvature in all cases of smooth robot trajectory scenarios, Fig. 9 -Fig. 11 show the profiles of the attained smooth curves. Here, the x-axis denotes the arc-length of the smooth curve, and the y-axis denotes the value of the radius of curvature ρ = 1 κ . By observing Fig. 9 -Fig. 11, we note the following facts: • Under small number of evaluations ( = 2000), RBDE, SADE and JADE attain a linear-like variation of radius VOLUME 8, 2020    of curvature over all studied cases of smooth robot trajectories.
• In particular, the linear-like behaviour is attined when the greedy initialization scheme is used. The above suggests a monotonically varying curvature, which is due to the presence of the smoothness factor H (Eq. 12).
• Under large number of evaluations ( = 2000), all algorithms attain the linear-like variation of the radius of curvature over all studied cases. The above observations confirm the feasibility of generating curves which not only are able to approximate robot trajectories reasonably, but also are able to show linear-like variation of curvature. These features are useful to suggest alternative robot trajectories which comply with linear variation of torque, and whose control laws are computationally efficient.
Furthermore, in order to exemplify the ability of compounding trajectories generated by smooth curves within the context of mobile robot navigation (Boe-Bot robot architecture), Fig. 12 shows examples of compounded trajectories by a plural number of fitted smooth curves. Here, curves are compounded considering C 0 continuity and by matching the gradient of the curveṙ8(u) at the connecting node points. Fig. 12 shows the compounded path segments in green color and the prescribed input coordinates from the robot trajectory in blue-colored nodes. In scenario 1, 5 arbitrary path segments were used to generate an arbitrary compounded trajectory, whereas in scenario 2, 6 arbitrary path segments were used. By observing the profile in Fig. 12, we can note that it becomes possible to compute compounded smooth paths by using the fitted smooth trajectories. However, as expressed in section II, the compounded paths will require further pruning to satisfy the minimal cost function F at the node joints. Also, the compounded paths are unable to satisfy collisionfree navigation. Our results offer the scheme to generate alternative navigation trajectories satisfying not only the closeness to real-world achieved trajectories, which are byproduct of the interaction of the robot dynamics and the environment, but also showing the linear-varying curvature profile, which is relevant for comfortability and safety in navigation.
In future work, we aim at deriving efficient methods for collision-free path planning of smooth paths, and integrating with the proposed aesthetic B-spline Curves [45]- [48]. Also, in our agenda is the development of compact trajectories in the plane and its application in resource distribution and interaction of Multi-Agent Sytems. Our obtained results offer the building blocks to further advance towards developing data-driven path planning algorithms for field areas, which may find use in several real-world applications in Robotics, Operations Research and Multi-Agent Systems.

IV. CONCLUSION
In this paper, we have proposed an approach to generate smooth paths from observed robot trajectories by optimizing criteria for fitting and smoothness using Differential Evolution under distinct modes of initialization of the population, selection pressure, exploration and exploitation during sampling. In particular, we the DE/rand/1/bin strategy (DERAND), the Rank-Based Differential Evolution (RBDE), the Adaptive Differential Evolution with External Archive (JADE), the Strategy Adaptation Differential Evolution (SADE), and the Differential Evolution with Similarity Based Mutation (DESIM).
Our rigorous experiments using all feasible cases of pulse width modulation on one side of the (servo) motor of the Boet-Bot mobile robot architecture showed the feasibility of our approach to generate smooth curves efficiently, in which fast convergence to the basins of the search space occurs with the greedy initialization scheme and Differential Evolution with exploitative and parameter adaptation schemes, such as RBDE, SADE and JADE.
Our proposed approach is useful to suggest alternative trajectories for mobile robots complying with linear variation of the radius of curvature, and whose control realization would be computationally efficient, offering the data-driven planning mechanisms for comfortability in riding.