Finite-Horizon Kinetic Energy Optimization of a Redundant Space Manipulator

Featured Application: The proposed optimal inverse kinematics method is suitable for space robotics applications, e.g., space servicing and active debris removal (ADR). It is also suitable for industrial robotics applications, e.g., in the planar operations of SCARA robots. Abstract: The minimization of energy consumption is of the utmost importance in space robotics. For redundant manipulators tracking a desired end-effector trajectory, most of the proposed solutions are based on locally optimal inverse kinematics methods. On the one hand, these methods are suitable for real-time implementation; nevertheless, on the other hand, they often provide solutions quite far from the globally optimal one and, moreover, are prone to singularities. In this paper, a novel inverse kinematics method for redundant manipulators is presented, which overcomes the above mentioned issues and is suitable for real-time implementation. The proposed method is based on the optimization of the kinetic energy integral on a limited subset of future end-effector path points, making the manipulator joints to move in the direction of minimum kinetic energy. The proposed method is tested by simulation of a three degrees of freedom (DOF) planar manipulator in a number of test cases, and its performance is compared to the classical pseudoinverse solution and to a global optimal method. The proposed method outperforms the pseudoinverse-based one and proves to be able to avoid singularities. Furthermore, it provides a solution very close to the global optimal one with a much lower computational time, which is compatible for real-time implementation.


Introduction
In space robotics missions, both for space servicing or debris capture, the minimization of the robot energy consumption and the minimization of the reactions transferred to the base spacecraft are of crucial importance [1]. Both of these targets are related to a longer system lifetime in orbit: while energy minimization allows longer operations by definition, reactions minimization allows to keep the base spacecraft orientation unchanged, which in turn means fuel savings due to the fact that attitude corrections are not necessary. In the case of trajectory tracking of redundant manipulators, inverse kinematics methods are usually adopted, which however only provide local minima along the end-effector trajectory. Thus, as trajectories increase in length, they fail to follow the lowest energy path in the robot joints space, and tend to be singularity prone.
Differential kinematic inversion has been extensively used for the solution of the realtime motion of redundant manipulators, while minimizing a certain cost function [2][3][4][5][6][7]. It has been used in plenty of variations, as it has several interesting properties: it has fast execution times, can be used to prioritize specific tasks over others, and it can be used to optimize a variety of different cost functions while driving the robot end-effector through a predefined path. However, this comes at a price: it only provides local optima, which can be quite far from the global one [2], and it is subject to the problem of singularities. In field applications [8], like space or underwater robotics, avoiding singularities and using as little energy and power as possible is important as robots usually run on limited energy sources, such as batteries or solar panels. In such applications, local inverse kinematics frameworks are not ideal, and, on the other hand, global optimization methods are too computationally demanding. Minimization of energy consumption of robotic systems has been subject of intense research, since it helps reducing operational costs, and allows us to use robots even when resources are limited. Focusing on robot manipulators, many attempts to minimize kinetic energy through inverse kinematics have been focused on minimizing the joint velocities norms. To the best of the authors' knowledge, the first systematic review of this topic has been proposed by Klein and Huang [9], which is a milestone in the research on the use of the pseudoinverse of the Jacobian matrix for robot control. The paper outlines how the pseudoinverse generally keeps the required instantaneous power low, since it minimizes joint velocities, which has been considered as equivalent to approximately minimize the kinetic energy. It also states that the pseudoinverse method, since it minimizes joint velocities, keeps the manipulator away from singularities, which was however proven incorrect by Bailleul et al. [10]. Many of the subsequent works developed their algorithms starting from the concept of weighted pseudoinverse, as presented by Whitney [11]. In particular, several researchers have followed the direction of pointwise minimizing the joint torques as a means to minimize the actuator energy consumption, associating such an algorithm with positive results regarding the minimization of the kinetic energy integral along the whole trajectory. Extensive analysis on this approach is presented by Hollerbach and Suh [12], who however found out that a dynamics-based inverse kinematics resolution method can lead to instabilities on longer trajectories. A more strictly energy-focused approach was followed by Nedungadi et al. [13], who proposed a method for joint torques minimization that holds some global properties. However, it does not allow to minimize kinetic energy when initial conditions on joint displacements are imposed and, although less instable than the previous work, still presents stability issues. The first research providing actuator models for energy consumption minimization was presented by Vukobratovic et al. [14], who introduced the dynamic models of hydraulic and electrical actuators and showed that a major advantage of energy minimization compared to joint velocities minimization is that it tends to move joints with lower inertia more than those with higher one, thus reducing energy consumption. As pointed out by Nenchev [2] in his review, these early studies made it very clear that a trade-off between local optimization performance and global stability is generally unavoidable: pointwise minimization must always be relaxed in order not to get stuck in singularities. Several works have been published on stabilization of joint torques minimization methods: most of them are aimed at minimizing joint torques without considering any robot payload (worthy of notice are [15] by Maciejewski and [16] by Hu et al.); nevertheless, the minimization of torques in the presence of external loads has also been addressed (Woolfrey et al. [7]). However, the solutions they propose are only marginally useful for the methods that minimize kinetic energy. Other recent developments of pseudoinverse-based methods are aimed at addressing specific shortcomings of existing formulations; for example Kelemen et al. [6] presented a solution to minimize computational times by setting a variable priority value for each end-effector task. Alongside pseudoinverse-based methods, other inverse kinematics methods based on Quadratic Programming have been proposed. These are considered promising since they have been proven capable of solving inverse kinematics problems with equality and inequality constraints [17]. These methods, however, have been proposed for joint torques minimization [18] rather than kinetic energy minimization.
Alongside local inverse kinematics methods, offline methods based on optimal control have been pursued. These techniques are slower, but allow for better optima of the cost function, although at the price of much longer computation times. Their cost function is usually a sum of the end-effector position error and a control cost (usually kinetic energy or actuators torque) computed on a specific set of points. The set may correspond to the whole trajectory or to a finite number of points after the last completed integration step.
In [19], Gregory et al. presented an optimal control method for energy optimization along a trajectory subject to holonomic constraints based on Euler-Lagrange equations. However, this method has only proven efficient for non-redundant manipulators. A different optimal control approach, based on Pontryagin maximum principle, is proposed in [20] by Halevi et al. for electromechanical linear actuators. In this work, a mathematical model of the actuators to be included in the cost function is presented, but the adaptability of this approach to manipulators with revolute joints is doubtful. Several other works address the issue of minimizing the energy required to perform the whole robot trajectory using offline methods. Among those, an early attempt was presented by von Stryk et al. [21], featuring an optimal control based minimization of the squared joint velocities along the trajectory, as a means to reduce the kinetic energy in the operations of a non-redundant manipulator. In [22], Saramago et al. investigated an offline path planning method for non-redundant manipulators using a multi-objective cost function involving both trajectory time and mechanical energy. In [23], Field et al. presented a solution of the energy minimization problem based on dynamic programming. Their intuition does not serve the purpose of finding the global optimum, since each one of the iterations only explores a subspace of the search space; however, it manages to avoid mediocre local minima, and keeps the execution time reasonably low avoiding the curse of dimensionality usually involved with dynamic programming. Nurmi et al. [24] recently presented an energy optimal solution for hydraulically powered redundant manipulators based on dynamic programming. Their approach, which has been developed for hydraulic actuators, performs sensibly better than general energy minimizing algorithms in the considered application, proving that actuator models might be beneficial for specific problems. However, the related burden of computational complexity means that this method cannot be used in real time. Recently, Tringali et al. [25] proposed a globally optimal inverse kinematics method for redundant manipulators, which allows for the optimization of different integral cost functions, such as kinetic energy and joint torques norm, can provide solutions with a variety of constraints, both linear and nonlinear, and is suitable for multi-objective optimization. Finally, Faroni et al. [26] presented a fast online predictive method to solve the task-priority differential inverse kinematics of redundant manipulators under kinematic constraints.
In this paper, a novel inverse kinematics method for redundant manipulators is presented, which is able to provide a solution near to the global optimal one, is able to avoid singularities, and is suitable for real-time implementation. The proposed method is based on the optimization of the kinetic energy integral on a limited subset of future endeffector path points, making the manipulator joints to move in the direction of minimum kinetic energy. Rather than optimizing the cost function at every time step, the cost function is updated at larger intervals, and the appropriate choice of the update interval and the number of future exploration points allows a limited computational cost. The proposed method is tested by simulation for a three degrees of freedom (DOF) planar manipulator in a number of test cases, and its performance is compared to the classical pseudoinverse solution and to the one computed with a global optimal method. The paper is organized as follows: Section 2 presents and discusses the proposed optimization method and its main parameters; then, Section 3 presents the simulation setup and results and, in particular, a comparison between the simulation results obtained with the proposed method, the classical pseudoinverse solution, and a global optimal method. Finally, Section 4 concludes the paper.

Mathematical Formulation
A kinematically redundant manipulator is a manipulator with more DOF than controlled end-effector variables. The tracking problem of a redundant manipulator is the problem of following a reference end-effector trajectory x ref , expressed as a function of time t in the Cartesian space. This requires to find the joint variables q(t) such that: (1) in which x is the actual end-effector trajectory. The redundancy implies that the problem expressed by (1) may have an infinite number of solutions in the joints space. The method hereby presented is based on the general solution computed with the Moore-Penrose pseudoinverse: in which J is the manipulator Jacobian, J + = J T J −1 J T is its Moore-Penrose pseudoinverse, .
x is the desired end-effector velocity, I is the identity matrix, and . q 0 is an arbitrary vector. The first part of (2) is a least-squares solution based on the minimization of the joint velocities norm and is related to the primary task (end-effector trajectory tracking), while the second part is related to a secondary task that can be imposed by exploiting the nullspace operator [27] I − J + J . The arbitrary vector . q 0 can be chosen so that it results in the minimization of an integral cost function along the specified end-effector trajectory, whereas the precision tracking is provided by the pseudoinverse solution (first term in (2)). A general expression of an integral cost function subject to the kinematic constraint related to the tracking of the desired end-effector trajectory is: where C q, q, t is a cost function to be integrated along the manipulator trajectory, t 0 and t f in are the initial and final times. Final time t f in does not need to be the time when the manipulator motion is completed: it is possible to define a time horizon h , much shorter than the whole motion time, so that the integral cost is optimized in a time window from t 0 to t 0 + h , and then updated during the motion (using each time an updated t 0 ), until t 0 + h corresponds to t f in . This method is used here to optimize the kinetic energy integral of a serial manipulator, leading to the following discrete formulation:  (3) is dropped in favour of a closer time The integral cost function is used to obtain a velocity vector . q prediction pointing to the direction in which the kinetic energy integral is minimum for the chosen trajectory. This velocity vector is used as an extra (secondary) task added to the manipulator motion through the Jacobian null space as per (2).
The steps of the proposed method, called Predictive Minimization of Kinetic Energy (PMKE) method, are presented in Table 1. Step nr. Operation

1
A horizon h is selected (h is a positive integer number).

2
A prediction step size ∆T is selected. It is a multiple of the pseudoinverse algorithm step size ∆t and it is generally much larger.

3
An update interval size I is selected. It is a multiple of the pseudoinverse step size ∆t, but it is generally smaller than ∆T.

4
Expression (4) is optimized with sampling times t = t 0 , t 0 + ∆T, . . . , t 0 + h∆T. Any deterministic optimisation algorithm can be used, and the Matlab fmincon Sequential Quadratic Programming (SQP) function has been used in this work. This provides a prediction of the optimal joint variables for a future window that is defined here as prediction window T p = h * ∆T.

5
The optimization at point 4 outputs a set of h joint configurations, for the h + 1 sampling times t = t 0 , t 0 + ∆T, . . . , t 0 + h∆T. The set of these vectors is from now on referred as q prediction . A specific vector in the set is referred with its sampling time as a subscript, with a notation such as q h∆T for the sampling time t 0 + h∆T (q 0 for t 0 ).

6
A 4th order spline interpolation from q 0 to q h∆T is computed, which intersect all the points in the set q prediction .
Endpoint conditions of the interpolation are set so that . q and ..
q at q 0 are always set equal to the ones reached by the manipulator at time t 0 , and are set to zero at q h∆T . This ensures that the resulting spline has not large derivatives in the interval of interest. The resulting spline is from now on referred as q prediction = f(t).

7
For each step i in the update interval I, the joint velocities are computed as per (2), using the error between the previous joint variables q i−1 and the computed value of q prediction divided per ∆t as arbitrary vector: Once time I has passed, steps 4 to 7 are repeated for the new update interval, using the last time step of previous interval as t 0 . 9 Steps 4-8 are repeated until the amount of time before the end of the trajectory is T p . 10 The steps left at this point are computed with the last computed q prediction .

Parameters Description
The main parameters involved in the PMKE method are the inverse kinematics discrete time step ∆t, the prediction step size ∆T, the horizon h, and the update interval I. This section explains their meaning and the value that has been chosen for the simulations performed for the purposes of this work.
Inverse kinematics discrete time step ∆t: this parameter is the time step used to numerically perform the integration of the inverse kinematics equation (2). It has been set to 0.01 s for the implementation hereby presented.
Prediction step size ∆T and horizon h: these altogether determine the prediction window T p , which is related to the ability of the method to take into account future (up to T p after the current time t 0 ) issues that could undermine the quality of the solution. The two parameters must be balanced considering that, for example, for a fixed T p , having a higher h and a smaller ∆T allows for better resolution of the prediction, making it less likely for the method to pursue local improvements over global ones. However, it also implies longer computational times. For the 3-DOF manipulator model used in the simulations (see Section 3), good results are reached with h = 2. Thus, this value is used to avoid the computational complexity related to higher h values. This means that ∆T = T p /2, and leaves the question open about how to choose T p . Increasing it does not change the computational complexity, which is only influenced by h, however a prediction too far in the future might feature a reduced influence on the current path point. Furthermore, it is hard to establish a proper T p because a value that is good for certain trajectories does not necessarily fit other ones: some of them might benefit from a longer prediction window, while some others might benefit from a shorter one, which makes it difficult to develop a tuning strategy. The value of the prediction window is set to 20% of the total simulation time. This means that, for trajectories with simulation time of 10 s such as the ones used in the following section, the method will perform the optimization of the cost function up to 2 s, corresponding to 200 inverse kinematics time steps. This means that ∆T is equal to 1 s. Update interval I: the update interval is the variable that decides how often the optimization of the cost function is updated. It should be ideally set as big as the time step ∆t, but this would slow down the method. For the simulations presented in this paper, an update interval of 0.5 s is chosen to show that this is sufficient to obtain good results. However, much shorter I intervals are also possible if computational time is not a concern. The computational time required for the PMKE method is discussed in more detail in Section 3.2.
The values of the main parameters of the PMKE method used in this work are summarized in Table 2.

Simulation Setup
A 3-DOF planar manipulator model (with three revolute joints) is considered, which corresponds to a real robot prototype for space applications that has been used in previous works [3][4][5]25]. The two end-effector Cartesian coordinates are controlled, x and y, thus leaving one degree of kinematic redundancy. Geometrical and inertial characteristics of the manipulator, which is depicted in Figure 1, are presented in Table 3.
thermore, it is hard to establish a proper because a value that is good for certain trajectories does not necessarily fit other ones: some of them might benefit from a longer prediction window, while some others might benefit from a shorter one, which makes it difficult to develop a tuning strategy. The value of the prediction window is set to 20% of the total simulation time. This means that, for trajectories with simulation time of 10 s such as the ones used in the following section, the method will perform the optimization of the cost function up to 2 s, corresponding to 200 inverse kinematics time steps. This means that is equal to 1 s. Update interval : the update interval is the variable that decides how often the optimization of the cost function is updated. It should be ideally set as big as the time step , but this would slow down the method. For the simulations presented in this paper, an update interval of 0.5 s is chosen to show that this is sufficient to obtain good results. However, much shorter intervals are also possible if computational time is not a concern. The computational time required for the PMKE method is discussed in more detail in Section 3.2.
The values of the main parameters of the PMKE method used in this work are summarized in Table 2.

Simulation Setup
A 3-DOF planar manipulator model (with three revolute joints) is considered, which corresponds to a real robot prototype for space applications that has been used in previous works [3][4][5]25]. The two end-effector Cartesian coordinates are controlled, and , thus leaving one degree of kinematic redundancy. Geometrical and inertial characteristics of the manipulator, which is depicted in Figure 1, are presented in Table 3.  A simulator of the manipulator inverse kinematics has been developed using the Matlab programming language, except for the code used for steps 4-5 of the PMKE method described in Section 2.1, which corresponds to the minimization of the kinetic energy integral: the routine for this optimization has been first written in Matlab exploiting the fmincon function, and then translated into C, compiled, and interfaced to the inverse kinematics simulator.
In order to show the effectiveness of the presented method, four trajectories have been simulated with both the PMKE and a traditional method, based on the Moore-Penrose pseudoinverse and performing the minimization of the joint velocity vectors through local least-squares optimization. This method is called here Least Squares Velocities (LSV) minimization method. Moreover, a further comparison is made with the results computed with a global optimal method [25] recently published by the authors.
Simulations 1-2 feature two short end-effector trajectories aimed at demonstrating the PMKE method general performance. Simulation 3 features a trajectory in which the LSV method encounters a singularity. Finally, simulation 4 features a trajectory whose length causes a significant degradation of the performance of the LSV method. All of the proposed trajectories start from the same end-effector position x init = [0.468; 0] m and robot initial configuration q init = [0; 0.327; −0.754] rad, and run for a total time T = 10 s. The considered trajectories and their characteristic dimensions are presented in Table 4. The end-effector velocity profiles have been chosen in order to have sufficiently smooth accelerations (continuous and with continuous derivative), as shown in Figure 2. The equation of the derivative of the curvilinear abscissa of the end-effector used is: where T is the total motion time, L the overall length of the considered trajectory, and β = 2π T .
Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 16 Figure 2. End-effector curvilinear abscissa and its derivatives. Table 5 compares the kinetic energy integral for all the considered methods and trajectories, both in absolute and relative terms. The PMKE method presented in this paper   Table 5 compares the kinetic energy integral for all the considered methods and trajectories, both in absolute and relative terms. The PMKE method presented in this paper is used as a reference and highlighted in grey. PMKE computational times have also been included, as measured on a machine featuring an Intel i7 ninth generation exa-core processor with 32 GB RAM, SSD mass storage, and using Windows 10 and Matlab version 2019b.

Simulation Results
It can be noticed that the LSV method is outperformed by the PMKE one for all the trajectories: non-singular trajectories feature relative differences ranging from 30.1% to 47.4%, while the trajectory that is singular using the LSV method (trajectory 3) is not singular using the PMKE method, resulting in a very high improvement of kinetic energy integral (relative difference of 187.36%). The absolute value of the difference between PMKE solutions and global optimal solutions, on the other hand, is very small (0.13%) for short trajectories and grows with the length of the trajectory (up to 23.6%). This can be explained by observing that, the longer the distance covered by the end-effector, the higher the chance of an energy consuming motion taking place beyond the PMKE prediction window. Kinetic energy integral profiles are shown in detail in Figure 3.  It has also to be noticed that the execution times of the PMKE method (see Table 5) allow for a real-time kinematic inversion, as it is discussed below. The core of the PMKE method performs two operations: the optimization (steps 4-5 of the method as described in Section 2.1), and the generation of the coefficients of the 4th order splines used for interpolation (step 6). A hybrid C-Matlab implementation has been deployed to perform these operations. As a result, the optimization of steps 4-5 is performed in C in an average time of ~0.0011 s, while the computation of the spline coefficients of step 6 is performed in Matlab in an average time of ~0.0010 s. The overall average computational time for the It has also to be noticed that the execution times of the PMKE method (see Table 5) allow for a real-time kinematic inversion, as it is discussed below. The core of the PMKE method performs two operations: the optimization (steps 4-5 of the method as described in Section 2.1), and the generation of the coefficients of the 4th order splines used for interpolation (step 6). A hybrid C-Matlab implementation has been deployed to perform these operations. As a result, the optimization of steps 4-5 is performed in C in an average time of~0.0011 s, while the computation of the spline coefficients of step 6 is performed in Matlab in an average time of~0.0010 s. The overall average computational time for the two operations is~0.0021 s, which is compatible with real-time implementation, since the considered time step ∆t is equal to 0.01 s. Even faster computation can be obtained by implementing both of these operations in C/C++ and optimizing the code.
In order to make a more detailed comparison of the results obtained with the PMKE and LSV methods, Figures 4-7 show (top to bottom) the stroboscopic views of the robot motion, joint positions, and joint velocities for both the PMKE (left) and LSV (right) methods. In each stroboscopic plot, the joint 1 (which allows the rotation between link 1 and the base frame) is in the origin and, in the joint positions and velocities plots, joint 1 is shown in blue, joint 2 in red, and joint 3 in yellow.   PMKE results for trajectories 1 and 2 feature a behavior that has been previously noticed when using kinetic energy saving algorithms [5,13]: the joints that are closer to the end-effector feature higher velocities compared to those closer to the base of the manipulator. This behavior is explained by considering that, when a joint is rotated, the same rotation is also transferred to the successive joints; thus, the moment of inertia of a link is   PMKE results for trajectories 1 and 2 feature a behavior that has been previously noticed when using kinetic energy saving algorithms [5,13]: the joints that are closer to the end-effector feature higher velocities compared to those closer to the base of the manipulator. This behavior is explained by considering that, when a joint is rotated, the same PMKE results for trajectories 1 and 2 feature a behavior that has been previously noticed when using kinetic energy saving algorithms [5,13]: the joints that are closer to the end-effector feature higher velocities compared to those closer to the base of the manipulator. This behavior is explained by considering that, when a joint is rotated, the same rotation is also transferred to the successive joints; thus, the moment of inertia of a link is constituted by the sum of its own moment of inertia and an additional term for each of the successive joints. While the LSV method simply minimizes the joints velocity norm, the PMKE method, similarly to other methods that minimize kinetic energy, takes into account the fact that the motion of the joints closer to the base are more energetically expensive. Stroboscopic and joint velocities plots in Figures 4 and 5 confirm that the LSV method features a joint motion more uniform than that of the PMKE method. In fact, for trajectory 1, the first joint features a maximum velocity of 0.063 rad/s in the LSV solution and 0.017 rad/s in the PMKE solution, while the third joint reaches up to 0.197 rad/s in the LSV solution and 0.203 rad/s in the PMKE solution. While the maximum velocity of the third joint is around three times the velocity of the first one in the LSV solution, the third joint is more than 10 times faster than the first one in the PMKE solution. Similarly, for trajectory 2, the first joint has a peak velocity of 0.196 rad/s in the LSV solution and a peak velocity of 0.095 rad/s in the PMKE solution, while the third joint velocity reaches 0.529 rad/s in the LSV solution and 0.613 rad/s in the PMKE solution.
This behavior is typical of local kinetic energy minimization algorithms. However, PMKE is able to depart from it when required to obtain improved optimization results. A very clear example of this case is provided by trajectory 3, for which the simulation results are presented in Figure 6. This trajectory is singular using the LSV method, while, on the other hand, the PMKE method is able to avoid the singularity. This remarkable behavior is allowed by the fact that the PMKE method, although generally striving to keep the first joint velocity low, moves it (and the other ones) away from singularities when required, if they fall within its prediction window. The stroboscopic plots in Figure 6 show that, while LSV (right) hits a singularity at t = 6.     This behavior is typical of local kinetic energy minimization algorithms. However, PMKE is able to depart from it when required to obtain improved optimization results. A very clear example of this case is provided by trajectory 3, for which the simulation results are presented in Figure 6. This trajectory is singular using the LSV method, while, on the other hand, the PMKE method is able to avoid the singularity. This remarkable behavior is allowed by the fact that the PMKE method, although generally striving to keep the first joint velocity low, moves it (and the other ones) away from singularities when required, Finally, the PKME-LSV comparison for trajectory 4 is presented in Figure 7. This is the longest trajectory that has been simulated, and it would be impossible to cover its whole length while also limiting the motion of the first joint. In this case PMKE and LSV feature similar joint velocity profiles, and PMKE features a consistently higher first joint maximum velocity (sitting at 0.780 rad/s, against 0.512 rad/s in the LSV solution). In this case, the PMKE method provides a solution with a lower kinetic energy integral by making the robot to assume configurations that feature lower coefficients of the inertia matrix, compared to those obtained with the LSV method.
In order to verify that the coefficients of the inertia matrix are lower for the PMKE solution, a possible indicator is the maximum eigenvalue of the manipulator inertia matrix, which is known as its spectral radius. In fact, the spectral radius of a square matrix A is defined as: ρ(A) = max {|λ 1 |, . . . , |λ n |} (6) where λ 1 , . . . , λ n are the eigenvalues of the matrix A. Figure 8 shows the evolution of the spectral radius of the inertia matrix of the robot manipulator along trajectory 4 for both the LSV (dotted line) and PMKE (dashed line) solutions. if they fall within its prediction window. The stroboscopic plots in Figure 6 show that, while LSV (right) hits a singularity at = 6. Finally, the PKME-LSV comparison for trajectory 4 is presented in Figure 7. This is the longest trajectory that has been simulated, and it would be impossible to cover its whole length while also limiting the motion of the first joint. In this case PMKE and LSV feature similar joint velocity profiles, and PMKE features a consistently higher first joint maximum velocity (sitting at 0.780 rad/s, against 0.512 rad/s in the LSV solution). In this case, the PMKE method provides a solution with a lower kinetic energy integral by making the robot to assume configurations that feature lower coefficients of the inertia matrix, compared to those obtained with the LSV method.
In order to verify that the coefficients of the inertia matrix are lower for the PMKE solution, a possible indicator is the maximum eigenvalue of the manipulator inertia matrix, which is known as its spectral radius. In fact, the spectral radius of a square matrix A is defined as: where , … , are the eigenvalues of the matrix A. Figure 8 shows the evolution of the spectral radius of the inertia matrix of the robot manipulator along trajectory 4 for both the LSV (dotted line) and PMKE (dashed line) solutions. The figure shows that the maximum eigenvalue of the inertia matrix for the PMKE trajectory is always lower than the maximum eigenvalue for the LSV one, indicating that the matrix coefficients are generally lower. Particularly, the maximum eigenvalue for PMKE is as low as 0.0556 for = 6.33 s, while for LSV its minimum value is 0.0998 for = 5.71 s. The figure shows that the maximum eigenvalue of the inertia matrix for the PMKE trajectory is always lower than the maximum eigenvalue for the LSV one, indicating that the matrix coefficients are generally lower. Particularly, the maximum eigenvalue for PMKE is as low as 0.0556 for t = 6.33 s, while for LSV its minimum value is 0.0998 for t = 5.71 s.
While the comparison with the LSV method can help to understand how the PMKE method stands with respect to another local method, it is also worth noticing that the solutions obtained with the PMKE method for the considered trajectories are very similar to the globally optimal ones. Figures 9-12 show the PMKE solutions on the left and the global optimal solutions on the right, both stroboscopic plots and joint variables.      For trajectories 1 and 2, the similarity can be easily noticed in both the stroboscopic and joint variables plots. The longer trajectories (3 and 4) feature some difference in terms of joint positions; however, the solutions are still very similar as it can be noticed in the stroboscopic plots. A complete overlap is not possible with the limited prediction window of the PMKE method: the longer the trajectory is, the higher the chance that local effects will influence the final energetic cost. However, the PMKE method is able to compute the best joint velocities to perform kinetic energy minimization with the chosen prediction window, and these are very similar to the ones of the global optimal solution.

Conclusions
In this paper, a novel inverse kinematics method for redundant manipulators, called Predictive Minimization of Kinetic Energy (PMKE) method, is presented, which is able to provide a solution near to the global optimal one, is able to avoid singularities, and is  For trajectories 1 and 2, the similarity can be easily noticed in both the stroboscopic and joint variables plots. The longer trajectories (3 and 4) feature some difference in terms of joint positions; however, the solutions are still very similar as it can be noticed in the stroboscopic plots. A complete overlap is not possible with the limited prediction window of the PMKE method: the longer the trajectory is, the higher the chance that local effects will influence the final energetic cost. However, the PMKE method is able to compute the best joint velocities to perform kinetic energy minimization with the chosen prediction window, and these are very similar to the ones of the global optimal solution.

Conclusions
In this paper, a novel inverse kinematics method for redundant manipulators, called Predictive Minimization of Kinetic Energy (PMKE) method, is presented, which is able to provide a solution near to the global optimal one, is able to avoid singularities, and is For trajectories 1 and 2, the similarity can be easily noticed in both the stroboscopic and joint variables plots. The longer trajectories (3 and 4) feature some difference in terms of joint positions; however, the solutions are still very similar as it can be noticed in the stroboscopic plots. A complete overlap is not possible with the limited prediction window of the PMKE method: the longer the trajectory is, the higher the chance that local effects will influence the final energetic cost. However, the PMKE method is able to compute the best joint velocities to perform kinetic energy minimization with the chosen prediction window, and these are very similar to the ones of the global optimal solution.

Conclusions
In this paper, a novel inverse kinematics method for redundant manipulators, called Predictive Minimization of Kinetic Energy (PMKE) method, is presented, which is able to provide a solution near to the global optimal one, is able to avoid singularities, and is suitable for real-time implementation. The PMKE method has been tested by simulation for a three degrees of freedom (DOF) planar manipulator in a number of test cases, and its performance has been compared to the classical pseudoinverse solution and to the one computed with a global optimal method. The simulations showed that the proposed method outperforms the pseudoinverse-based one and is suitable to avoid singularities. Furthermore, it provides a solution very close to the global optimal one with a much lower computational time, which is compatible for real-time implementation. The optimization of different cost functions (e.g., torques norm integral), the simulation of the performance in 3D and higher-DOF manipulators, a more optimized C/C++ implementation of the method to allow for faster computation, and the implementation in a real robot prototype will be part of future work.
Author Contributions: Conceptualization, S.C. and A.T.; methodology, S.C. and A.T.; software, A.T.; investigation, S.C. and A.T.; writing-original draft preparation, S.C. and A.T.; writing-review and editing, S.C. and A.T.; supervision, S.C. All authors have read and agreed to the published version of the manuscript.