A Computationally Efficient Algorithm to Find Time-Optimal Trajectory of Redundantly Actuated Robots Moving on a Specified Path

TL;DR

This paper introduces a computationally efficient algorithm for determining the time-optimal trajectory of redundantly actuated robots on a specified path, addressing unresolved issues in solution curve calculation.

Contribution

It proposes a new method to efficiently establish solution curves and calculate maximum or minimum accelerations based on actuator saturation patterns.

Findings

Reduced computational cost in trajectory planning

Effective calculation of maximum velocity curves

Validated with numerical examples

Abstract

A time-optimal problem for redundantly actuated robots moving on a specified path is a challenging problem. Although the problem is well explored and there are proposed solutions based on phase plane analysis, there are still several unresolved issues regarding calculation of solution curves. In this paper, we explore the characteristics of the maximum velocity curve (MVC) and propose an efficient algorithm to establish the solution curve. Then we propose a straightforward method to calculate the maximum or minimum possible acceleration on the path based on the pattern of saturated actuators, which substantially reduces the computational cost. Two numerical examples are provided to illustrate the issues and the solutions.

Tables

Table 1. Table 1: Physical characteristics of the System

length (m)		mass (kg)
$l_0$	0.2	$m_0$	1
$l_1$	0.5	$m_1$	1
$l_2$	0.6	$m_2$	1
$l_3$	0.3	$m_3$	0.3
$b_0$	1.4
limit torques (N.m)
${𝝉}_{limit}=\pm {[35,25,10,35,25,10]}^T$

Table 2. Table 2: Evaluating the efficiency of proposed time-optimal algorithm by comparison with reference [ 22 ]

index	time consumption (sec)
example	I		II
method	ref [22]	proposed	ref [22]	proposed
MVC	5.4	1.2	696	594
solution curve	12.6	5.6	1608	402
time-optimal path	18	6.8	2304	996
reduction of time consumption (%)	62.2		56.8

Equations

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}+\mathbf{h}(\mathbf{q},\dot{\mathbf{q}})+\mathbf{g}(\mathbf{q})=\mathbf{B}(\mathbf{q})\mathbf{\tau }$$

$$x_i=f_i(s)i=1,\dots ,r$$

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}+\mathbf{h}(\mathbf{q},\dot{\mathbf{q}})+\mathbf{g}(\mathbf{q})=\mathbf{B}(\mathbf{q})\mathbf{\tau }$$

$$\mathbf{x}=\mathbf{f}(s)$$

$${\mathbf{\tau }}_{min}\le \mathbf{\tau }\le {\mathbf{\tau }}_{max}$$

$$\mathbf{q}(t_i)={\mathbf{q}}_i,\dot{\mathbf{q}}(t_i)={\dot{\mathbf{q}}}_i$$

$$\mathbf{q}(t_f)={\mathbf{q}}_f,\dot{\mathbf{q}}(t_f)={\dot{\mathbf{q}}}_f.$$

$$\mathbf{c}(s)\ddot{s}+\mathbf{d}(s){\dot{s}}^2+\mathbf{e}(s)=\mathbf{B}(s)\mathbf{\tau }$$

$${\mathbf{\tau }}_{min}\le \mathbf{\tau }\le {\mathbf{\tau }}_{max}$$

$$s(t_i)=0,\dot{s}(t_i)={\dot{s}}_i$$

$$s(t_f)=1,\dot{s}(t_f)={\dot{s}}_f.$$

$$\mathbf{c}(s)\ddot{s}+\mathbf{d}(s){\dot{s}}^2+\mathbf{e}(s)=\mathbf{B}(s)\mathbf{\tau }$$

$${\mathbf{\tau }}_{min}\le \mathbf{\tau }\le {\mathbf{\tau }}_{max}.$$

$$\mathbf{c}(s)\ddot{s}+\mathbf{d}(s){\dot{s}}^2+\mathbf{e}(s)=\left[\begin{array}{cc}{\mathbf{B}}_{ns}& {\mathbf{B}}_s\end{array}\right]\left[\begin{array}{c}{\mathbf{\tau }}_{ns}\\ {\mathbf{\tau }}_s\end{array}\right]$$

$$\left[\begin{array}{cc}\mathbf{c}(s)& -{\mathbf{B}}_{ns}\end{array}\right]\left[\begin{array}{c}\ddot{s}\\ {\tau }_{ns}\end{array}\right]=-\mathbf{d}(s){\dot{s}}^2-\mathbf{e}(s)+{\mathbf{B}}_s{\mathbf{\tau }}_s.$$

$${\overline{c}}_i(s^{†})=0.$$

$$\mathbf{c}(s)\ddot{s}+\mathbf{d}(s){\dot{s}}^2+\mathbf{e}(s)=\left[\begin{array}{cc}{\mathbf{B}}_a& {\mathbf{B}}_b\end{array}\right]\left[\begin{array}{c}{\mathbf{\tau }}_a\\ {\mathbf{\tau }}_b\end{array}\right].$$

$$\overline{\mathbf{c}}(s)\ddot{s}+\overline{\mathbf{d}}(s){\dot{s}}^2+\overline{\mathbf{e}}(s)={\mathbf{\tau }}_a+\overline{\mathbf{B}}{\mathbf{\tau }}_b={\mathbf{\tau }}^{\prime }$$

$$\{\begin{array}{l}x(s)=0.6(s-0.5)+0.7\\ y(s)=-3.46s^5+8.66s^4-5.77s^3+0.58s+0.7\\ \gamma (s)=-0.45(s-0.5)\end{array}$$

$$\{\begin{array}{l}x(s)=0.2cos(2\pi s)+0.7\\ y(s)=0.2sin(2\pi s)+0.62\\ \gamma (s)=0.7s\end{array}$$

$$\{\begin{array}{l}x(s)=0.2cos(2\pi s)+0.77\\ y(s)=0.15sin(2\pi s)+0.122\\ z(s)=0.1cos(2\pi s)-0.1\\ \alpha (s)=\frac{\pi }{8}sin(2\pi s)\\ \beta (s)=\frac{\pi }{3}sin(2\pi s)+\frac{\pi }{2}\\ \gamma (s)=\frac{\pi }{12}sin(2\pi s)+\frac{\pi }{6}\end{array}$$

Full text

††volume: xx††issue: x

{summary}

A time-optimal problem for redundantly actuated robots moving on a specified path is a challenging problem. Although the problem is well explored and there are proposed solutions based on phase plane analysis, there are still several unresolved issues regarding calculation of solution curves. In this paper, we explore the characteristics of the maximum velocity curve (MVC) and propose an efficient algorithm to establish the solution curve. Then we propose a straightforward method to calculate the maximum or minimum possible acceleration on the path based on the pattern of saturated actuators, which substantially reduces the computational cost. Two numerical examples are provided to illustrate the issues and the solutions.

1 Introduction

Time-optimal motion planning remains an ongoing challenge for redundantly actuated robots, such as cooperative multi-manipulator systems, or humanoid robots during double support phase. While several approaches have been proposed to tackle this issue, large computational cost has yet remained as a major obstacle. Here, we propose a computationally efficient algorithm to address the time-optimal planning problem for robots with actuation redundancy.

Several studies have examined the problem of minimum time motion for serial and parallel manipulators while operating on a specified path. In general, there are two approaches to deal with such problem, including, the Pontryagin maximum principle, and convex optimization. In the former approach, a system that moves in a specified path is regarded as having only one degree of freedom (DOF). Accordingly, it is proven that a bang-bang solution, in terms of the acceleration of the tip along the path, $\ddot{s}$, accounts for the problem. This means that the solution curve in the $s-\dot{s}$ plane could be obtained by successively integrating the maximum and minimum values of $\ddot{s}$ [1]. This method, mainly known as the numerical integration method [2], has been widely used in previous studies for time-optimal solutions in both normally actuated [1, 3, 4, 5, 6, 7, 8] and over-actuated [9, 10, 11, 12, 13, 14, 15, 16, 17] robots. Guaranteed convergence of solution, low computational cost, and providing the solution directly in the trajectory space are of main advantages to this approach [22].

In the convex optimization approach, the time-optimal path tracking problem is transformed into a convex optimal control problem by discretization the path $s$ into $N$ segments [18, 19, 20]. Recently, this method has been used to find the minimum time motion of humanoid robots in multi-contact tasks [21]. The convex optimization method is indeed a more generic approach, as it not only applies to finding time-optimal solutions, but can also be used to minimize any other index function. Moreover, there are available commercial packages that exclusively deal with the convex optimization problem, making it a more accessible approach.

In this paper, we mainly focus on the direct integration method for the time-optimal solution of redundantly actuated robots, while moving on a specified path. We aim to try and further reduce the computational costs associated with this method. One should note that parallel manipulators are an important example of such systems for which the proposed method can be used.

In the numerical integration method, the solution process consists of two parts: finding the extremum acceleration in each step of solution, and extracting the switching points. A primary approach to finding the extremum acceleration was first proposed in [1] for non-redundant serial manipulators. The admissible range of $\ddot{s}$ is determined at each ($s$, $\dot{s}$) under actuation constraints. A more rigorous approach to find the maximum or minimum $\ddot{s}$ was proposed in [3] by taking advantage of geometrical concepts. In the latter method, inequality constraints on actuators are converted into a convex polygon in ${\dot{s}}^{2}-\ddot{s}$ plane. The solution provided in this study, however, was applied in serial robots with no redundant actuators. Recently, this concept has been extended to over-actuated robots based on polytope projection technique [22]. Alternatively, the linear programming (LP) method [9, 10, 11, 12] has also been used to determine extremum acceleration for redundantly actuated systems. LP, however, is a search based algorithm, which makes it numerically inefficient, and leads to large computational costs in the solution process.

The first issue addressed in this paper is how to reduce the computational cost of calculating maximum/minimum acceleration along the path. To this end, we propose a method based on the pattern of saturated actuators. According to this method, the saturation pattern of actuators is preserved in intervals of $s$. It, hence, reduces the challenges of solving LP [9, 10, 11, 12] or forming polygon [22] at each step of integration, to a much simpler problem of finding the pattern of saturated actuators at some specific $s$. It then follows by solving a linear set of equations of motion at each step.

The second significant problem is to find the switching points in cases where more than one switching point exists. Such problem occurs when the velocity of the manipulator along the path exceeds a maximum feasible value. In such cases, one needs to establish the maximum velocity curve (MVC) and find possible switching points on this curve, known as critical points [3]. One approach to find MVC is to search for the maximum feasible $\dot{s}$ at each point of the path. To this end, one needs to find $\dot{s}$ for which $\ddot{s}_{max}=\ddot{s}_{min}$. This technique has been previously applied for both serial and parallel manipulators in [1] and [9, 10, 11, 12], respectively. However, it is evident that this method is numerically inefficient due to successive calculation of $\ddot{s}_{max}$ and $\ddot{s}_{min}$. To reduce the computational cost, a method for direct calculation of MVC was suggested in [3] for non-redundant serial manipulators. Recently, the same technique has also been used to determine MVC for over-actuated robots [22]. Nevertheless, the construction of MVC for whole domain of $s$ is still costly and numerically expensive.

Among the previous studies, three different strategies have been proposed to determine the critical points. The first strategy that theoretically made an important advancement is based on the shooting method [1]. Here, critical points are determined by looking for a solution curve which comes in contact with the MVC without crossing it. However, from a practical point of view this approach was very difficult to apply and numerically inefficient. Another strategy which was introduced in [3] determines critical points by comparing the slope of MVC with the slope of solution curve for all points located on MVC. In both of the above-mentioned strategies, the construction of the MVC is necessary, which adds to the computational demands of the problem. This is while the MVC curve is not, per se, a solution curve, and nor does it play an important role in finding the switching points. What is crucial for finding the solution, however, is the lower boundary of the trapped area, as well as the locked area [4], which are both located under the MVC. This lower boundary, which itself is a solution trajectory, is in fact the switching curve. It seems, therefore, that it would be more efficient if one tries to construct the switching curve directly while skipping the construction of MVC as long as it is possible.

In this direction, a straightforward approach has been proposed in [14] to find jump points on the MVC (a subset of critical points that are also known as zero-inertia points [23, 3, 2]), to establish the switching curve based on that. Such advancement could substantially reduce the computational cost. In some cases, however, such as finding the time-optimal solution for a humanoid robot during double support phase [15], there remains a part of MVC whose construction is inevitable. Such a case motivated this study to further examine this issue and try to enhance the method proposed in [14]. In this method, we introduce a new strategy to find the critical points by constructing MVC for minimum domain of $s$. This strategy is based on direct construction of solution curve, taking advantage of analytical calculation of zero-inertia points, and to find the domain of $s$, if any, for witch no value of $\dot{s}$ on solution curve is found during previous step. Construction of MVC is only necessary for such value of $s$. We developed a new method to convert actuator constraints into a convex polygon in ${\dot{s}}^{2}-\ddot{s}$ plane for over-actuated robots. The proposed method which directly constructs the polygon reduces the computational cost compared with the method proposed in [22], which is based on iteration.

To summarize, in this study we tried to improve computational efficiency of solutions based on integration method for minimum time trajectory planning. The proposed approach is used for cooperative multi-manipulator system (CMMS) with redundant actuators, moving on a specified path. Three main contributions are made which are described in section 3. In section 3.1, an efficient method is proposed for computation of $\ddot{s}_{max}$ and $\ddot{s}_{min}$, which is to be found in each step of integration. In section 3.2, we discuss which part of MVC could be skipped by direct construction of solution curve, and which part is necessary to be constructed. The final contribution in this study is explained in section 3.3, where we describe a novel method to directly establish the convex polygon of constraints for redundantly actuated robot as an essential tool to construct MVC. This also helps to find pattern of saturated actuators. Our proposed method is applied to solving two examples presented in section 5, and the elapsed computation time is compared with one of presented methods in [22]. The results of this comparison are given in section 6.

2 Time Optimal Problem

An important set of redundantly actuated robots is CMMS which is considered to solved as an example in this paper. Assume that the CMMS shown in Fig. 1 is supposed to move a payload from an initial point to a final point on a specified path in minimum time subject to the actuator’s limits. To state the time-optimal problem let us consider a CMMS composed of $v$ similar, rigid non-redundant serial manipulators working on a single object. The motion of the payload in task space is defined by $r$ coordinates; $\boldsymbol{x}=[x_{1},\dotso,x_{r}]^{T}\in\mathbb{R}^{r}$ and the motion of the system is defined with $n$ variables, $q_{1},\dotso,q_{n}$ where $\boldsymbol{q}=[\theta_{11},\dotso,\theta_{1r},\dotso,\theta_{v1},\dotso,\theta_{vr}]^{T}\in\mathbb{R}^{vr}$ and $\theta_{ij}$ denotes the $j$th joint coordinates of the $i$th manipulator. Due to the kinematic chain configuration, the system is subject to $p$ constraints. The reduced form of dynamic equations of such system can be written as

[TABLE]

where $\mathbf{M}(\boldsymbol{q})\in\mathbb{R}^{n-p\times n}$ is inertia matrix, $\boldsymbol{h}(\boldsymbol{q},\dot{\boldsymbol{q}})\in\mathbb{R}^{n-p}$ represents coriolis and centrifugal terms, $\boldsymbol{g}(\boldsymbol{q})\in\mathbb{R}^{n-p}$ indicates gravitational effects, $\mathbf{B}(\boldsymbol{q})\in\mathbb{R}^{n-p\times m}$ is actuators coefficient matrix and $\boldsymbol{\tau}\in\mathbb{R}^{m}$ is input torques [14].

A desired path can be defined by

[TABLE]

in which $x_{i}$ is one of the task space coordinates of manipulators and $s$ is an independent non-dimensional parameter that uniquely describes the position of the system along the path. Also, the value of $r$ is to be equal to $n-p$ due to the fact that each manipulator is kinematically non-redundant.

We may state the time-optimal problem as:

$Problem\ I:$ Find the optimum path $\boldsymbol{q}^{*}(t)$ to minimize the elapsed time $\Gamma=\int_{t_{i}}^{t_{f}}dt$ subject to

[TABLE]

with initial and final condition of

[TABLE]

To solve $problem\ I$ we restate it in terms of non-dimensional path parameter $s$. To this end, one must substitute for $\boldsymbol{q}$, $\dot{\boldsymbol{q}}$ and $\ddot{\boldsymbol{q}}$ in terms of $s$, $\dot{s}$ and $\ddot{s}$ in (1) to get

[TABLE]

in which all ${\boldsymbol{c}}$, ${\boldsymbol{d}}$ and ${\boldsymbol{e}}$ are members of $\mathbb{R}^{n-p}$ [14].

The relation (4) represents $n-p$ equations with two states $[s,\dot{s}]$. Any motion of the system which follows the prescribed path must satisfy all above $n-p$ equations. Now, the optimization problem can be rewritten as:

$Problem\ II:$ Find the optimum path $s^{*}(t)$ to minimize the elapsed time $\Gamma=\int_{t_{i}}^{t_{f}}dt$ subject to

[TABLE]

with initial and final condition of

[TABLE]

As we stated earlier in introduction, finding the solution of $problem\ II$ is equivalent to solve two major problems of finding extremum acceleration and extracting switching points. In the next section, after a detailed review of some background, we propose an efficient method to calculate the extremum acceleration based on pattern of saturated actuators.

3 Optimization Method

3.1 Saturation Pattern of Actuators

The form of $problem\ II$ resembles that of the time-optimal for non-redundant serial manipulators with the main difference being in matrix $\mathbf{B}$. In case of systems with no redundant actuators, this matrix is invertible and could be reduced into a unity matrix; which means that only one actuator does appear in each equation. However, for CMMS which has redundancy in actuators, $\mathbf{{B}}$ is a non-square matrix, which means that there are always more than one actuator torques in each equation of $problem\ II$. That is why the provided methods for serial manipulators in [1] and [3] cannot be used to calculate ${\ddot{s}}_{max}$ or ${\ddot{s}}_{min}$ for CMMS. This issue was independently addressed in [9] and [10], where the problem of finding extremum acceleration was introduced as:

$Problem\ III:$ Find $\boldsymbol{\tau}$ to maximize or minimize $\ddot{s}$ subject to

[TABLE]

This problem could be solved using a linear programming (LP) strategy. One should note that this problem must be solved for all values of $s$ and $\dot{s}$ on the solution curve, which, clearly, requires a huge amount of computation.

It is demonstrated in [13] that for a robotic system with $n$ generalized coordinates, $p$ holonomic constraints and $m$ actuators, at least $m-n+p+1$ actuators are in the saturated state. Therefore, solving $problem\ III$ results in finding $m-n+p+1$ saturated actuators which leads to $\ddot{s}$ being maximized or minimized. The solution gives the value of $\ddot{s}$ and $\boldsymbol{\tau}$. Here, we introduce a method based on the saturation pattern of actuators, which reduces the solution of $problem\ III$ into the solution of a set of linear equations of motion.

Considering the fact that when $\ddot{s}$ accepts its extremum, at least $m-n+p+1$ actuators are saturated [13], one can easily compute non-saturated actuators and $\ddot{s}$ from (4). To this end, we will have

[TABLE]

in which $\boldsymbol{\tau}_{ns}\in\mathbb{R}^{n-p-1}$ and $\boldsymbol{\tau}_{s}\in\mathbb{R}^{m-n+p+1}$.

For every point along the solution trajectory in $s-\dot{s}$ plane, one can always find a neighborhood in which the saturation pattern of actuators does not change. Therefore, we may use a set of saturated actuators at one point to calculate $\ddot{s}$ and $\boldsymbol{\tau}_{ns}$ for other points in such a neighborhood. So the problem of determining extremum acceleration, $problem\ III$, reduces to

$Problem\ IV:$ Assuming $\boldsymbol{\tau}_{s}$ to be known, find $\boldsymbol{\tau}_{ns}$ and $\ddot{s}$ from

[TABLE]

Clearly, the solution of (8) which is a linear set is much easier and involves less computation compared to (6), which is a linear programming problem. To solve the problem with hand, first, the saturated actuators set which leads to maximizing or minimizing $\ddot{s}$ is determined for one point by constructing the polygon of constraints in ${\dot{s}}^{2}-\ddot{s}$ plane as will be explained in section 3.3, or by solving LP problem, i.e. $problem\ III$, for that point. The next point in $s-\dot{s}$ plane is found by integration of equations of motion using this value of $\ddot{s}$ and $\boldsymbol{\tau}$. For the next step, we do not need to solve the LP problem. Instead, because $\boldsymbol{\tau}_{s}$ remains the same, we can calculate $\ddot{s}$ and $\boldsymbol{\tau}_{ns}$ from (8).

If none of the computed values for $\boldsymbol{\tau}_{ns}$ exceeds the actuator’s bound, the solution is acceptable; and therefore, we may proceed to the next step. However, if at least one of the values of $\boldsymbol{\tau}_{ns}$ exceeds its bound, it means that the pattern of saturation is changed and we need to establish the polygon to find real set of saturated actuators. Fig. 2 shows the flow diagram of the proposed algorithm which we would call SPA hereafter.

By applying this innovative method, while the $problem\ IV$ is solved in each integration step, determining the pattern of saturated actuators which needs more computational effort, is done only in some limited points, where the pattern of saturation changes. This method, hence, significantly reduces the amount of computation.

3.2 Critical Points

An important step to find the solution is determining the points to switch between the minimum and maximum acceleration. For this purpose, the first stage is calculating the critical points which are candidates for switching points on the MVC. The concept of critical points was first introduced in [3]. Points on the MVC are categorized as sink or source depending on whether the slope of solution curve in phase plane, $\ddot{s}/\dot{s}$, is greater or smaller than the slope of the MVC, $k^{\prime}$. Any solution starting from a source continues in feasible region, whereas any solution starting from a sink immediately leaves the feasible region. By critical points, we mean the points for which the nature of MVC is neither a sink nor a source.

In this section, we propose a new strategy to find the critical points. This method is an improvement of the method presented in [14]. This method involves two main parts. In the first part, we analytically determine zero-inertia points, i.e. the points which satisfy the following equation [14]

[TABLE]

Then we construct the solution curve using these zero-inertia points. In the second part, we construct the portion of MVC has no points on the solution curve with the same value of $s$. This will be discussed in more details in section 4 while describing the solution algorithm. One should note that not all $s^{\dagger}$ which satisfy (9) are feasible zero-inertia point. To check if $s^{\dagger}$ is a feasible point, it is just enough to check whether $\dot{s}^{\dagger}$ associated with $s^{\dagger}$ violates actuator constraints or not.

One should also note that not all feasible zero-inertia points located on MVC are candidates for switching points. Depending on the change of the nature of boundary on two sides of a zero-inertia point, it might be categorized as one of the four types of sink to source, source to sink, sink to sink and source to source. Out of these, only sink to source zero-inertia points are critical points. To realize the type of a zero-inertia point without constructing the MVC, one needs to establish the solution from such points in backward and forward directions with $\ddot{s}_{min}$ and $\ddot{s}_{max}$, respectively. It is done by integrating equations of motion backward and forward starting from that point. If the zero-inertia point is of the sink–sink type, for the forward direction, the actuator’s bounds are immediately violated, which means that the forward path instantly enters the non-feasible region (NFR). Also, if the point is of the source–source type, the actuator’s limits are immediately violated for the backward direction. Similarly, if the point is of the source–sink type, for both the forward and backward directions actuator’s bounds are violated immediately. However, solutions starting from a sink–source point remain in the feasible area.

In the developed strategy, in order to determine the switching points, only some parts of the MVC are constructed. It is more efficient than the other studies, e.g. [22], which establish MVC in the whole domain of $s$.

3.3 Polygon of Constraints

As mentioned in the previous section, there are points on the path for which we need to find $\ddot{s}_{max}$ or $\ddot{s}_{min}$ and the pattern of saturated actuators. In addition to that, one might also need to find the maximum possible $\dot{s}$ at each $s$ to establish MVC. One efficient method to solve both of the above mentioned problems is to establish the polygon of constraints which was first reported in [3] for serial manipulators.

To construct the polygon for CMMS systems, we attempt to extend the method proposed in [3]. For this, consider relation (4) which contains $n-p$ equations. Choosing $n-p$ independent actuators $\boldsymbol{\tau}_{a}\in\mathbb{R}^{n-p}$, one might rewrite relation (4) as

[TABLE]

Here by independent actuator, we mean a set of them which are sufficient to statistically balance the system at each $s$, regardless of the choice of $\boldsymbol{\tau}_{b}$. In above equation, $\boldsymbol{\tau}_{b}\in\mathbb{R}^{m-n+p}$ depicts the remaining actuators. Premultiplying (10) by ${{\mathbf{B}}_{a}}^{-1}$ leads to

[TABLE]

in which ${\overline{\boldsymbol{c}}}={{\mathbf{B}}_{a}}^{-1}{\boldsymbol{c}}\in\mathbb{R}^{n-p}$, ${\overline{\boldsymbol{d}}}={{\mathbf{B}}_{a}}^{-1}{\boldsymbol{d}}\in\mathbb{R}^{n-p}$, ${\overline{\boldsymbol{e}}}={{\mathbf{B}}_{a}}^{-1}{\boldsymbol{e}}\in\mathbb{R}^{n-p}$ and ${\overline{\mathbf{B}}}={{\mathbf{B}}_{a}}^{-1}{{\mathbf{B}}}_{b}\in\mathbb{R}^{n-p\times m-n+p}$. One should note that existence of ${{\mathbf{B}}_{a}}^{-1}$ is guaranteed by the virtue of independency of actuators $\boldsymbol{\tau}_{a}$.

In any equation of (11), all ${\boldsymbol{\tau}}_{b}$ and one of ${{\tau}_{a}}_{i}$, which equal to $m-(n-p)+1$ actuators appear. By selecting the appropriate values of $\boldsymbol{\tau}_{b}$ and $\tau_{{a}_{i}}$ from their lower or upper bound, one could obtain the minimum and maximum values for ${\tau}^{\prime}_{i}$. Therefore, each equation of (11) shows two straight lines on the plane $\dot{s}^{2}-\ddot{s}$. The area out of these two lines indicates a region in which for any $(\dot{s}^{2},\ddot{s})$ at least one of the actuator’s limits is violated so this point is non-feasible. The feasible area is restricted to points inside this polygon. Changing the choice of $\boldsymbol{\tau}_{b}$, other polygons are also generated. The common area of all these polygons is considered as the feasible area of CMMS. This area constitutes a unique polygon which works in the same way as the one established for serial manipulators [3].

The polygon constructed for CMMS is the common area restricted by $C(m,m-(n-p))\times n-p$ pairs of lines. Here, $C$ depicts the combination function. For any given $(\dot{s}^{2},\ddot{s})$ inside the polygon, there is one feasible $\boldsymbol{\tau}$ that satisfies the actuators saturation limits. On the edge of the polygon, $m-(n-p)+1$ actuators are in saturation conditions. $m-(n-p)$ are those of $\boldsymbol{\tau}_{b}$ and the other one is the corresponding ${{\tau}_{a}}_{i}$. At any vertex of the polygon, at least $m-(n-p)+2$ actuators are in saturation conditions. The set of saturated actuators on each edge is the common ones at two vertices located at the ends of the edge. Therefore, by knowing the saturated actuators at vertices of the polygon, the pattern of saturation on each edge is determined.

In this paper, we present a new method to establish the polygon for CMMS by finding its vertices. By applying this method, in addition to the position of vertices in ${\dot{s}}^{2}-\ddot{s}$ plane, the set of saturated actuators at each vertex is also obtained. Our proposed strategy to construct the polygon is as follows:

• $i.$

Choose a state of $\boldsymbol{\tau}_{b}$ and extract relation (11) for it.

• $ii.$

For each equation of (11), establish the corresponding pair of lines on $\dot{s}^{2}-\ddot{s}$ plane.

• $iii.$

Calculate all intersection points of each two lines which have the same value of $\boldsymbol{\tau}_{b}$.

• $iv.$

For each intersection point, compute the value of $\boldsymbol{\tau}_{a}$ from (11) and check whether $\boldsymbol{\tau}_{a}$ violates its saturation limits. If yes, this point could not be a vertex of polygon and should be ignored.

• $v.$

Continue step $i$ to $iv$ until all possible states of $\boldsymbol{\tau}_{b}$ are selected. After this, all vertices of the polygon are obtained.

The developed method is a straightforward and efficient strategy to find the polygon. It is different from the approach proposed in [22] that is based on the polytope projection technique. This technique iteratively attempts to find polygon by expanding the edges which leads to a trade-off between accuracy and computational cost.

4 Solution Algorithm

In this section, we provide the final algorithm to solve the time-optimal problem of redundantly actuated robots moving on a specified path. The algorithm utilizes strategies introduced in previous sections to find critical points and determine the pattern of saturated actuators. In phase plane, the higher solution trajectory has the shorter traveling time. By considering this fact, we propose the following algorithm to find the time-optimal trajectory:

$Step\ 1:$ Using SPA strategy, start from initial point, $D_{i}(0,\dot{s}_{i})$, in forward direction with $\ddot{s}_{max}$. Similarly, start from final point, $D_{f}(1,\dot{s}_{f})$, in backward direction with $\ddot{s}_{min}$. Continue these two procedures until either the line of $s=0$ or $s=1$ is crossed or actuator’s bounds are violated, i.e. solution trajectories enter the NFR (points $b_{i}$ and $b_{f}$ in Fig. 3). If one or both of the curves intersect the line of $s=0$ or $s=1$, solution trajectories cross each other at feasible region (FR). Therefore, the algorithm ends and the solution curve is obtained with one switching point. Otherwise, zero-inertia points should be determined, like the case shown in Fig. 3.

$Step\ 2:$ Apply the strategy proposed in section 3.2 to calculate a set of zero-inertia points with the type of sink–source which are located on MVC between $b_{i}$ and $b_{f}$, e.g. $C_{2}$ and $C_{3}$ in Fig. 3. If there are no such points skip to $Step\ 6$.

$Step\ 3:$ Run SPA strategy by starting from the lowest point of the set, $C_{2}$, in backward and forward directions with $\ddot{s}_{min}$ and $\ddot{s}_{max}$, respectively. This process is continued until either the solution trajectories generated in previous steps is crossed ($S_{1}$) or no acceptable solution fulfills the actuator’s limits (point $b_{2}$ on MVC).

$Step\ 4:$ Ignore all other points of the set, $C_{3}$, which are located above the generated solution curve inside the inadmissible area, i.e. trapped and locked area introduced in [4].

$Step\ 5:$ Repeat $Steps\ 3$ and $4$ from the lowest remaining points of the set. Continue the same process for other points of the set until either the generated solution curve covers the whole domain of $s$ from [math] to $1$ or there remains no zero-inertia point outside the trapped and locked area, like the case shown in Fig. 3. If both of the initial and final solution trajectories which are established in $Step\ 1$ are crossed with the generated curve, the algorithm ends and the solution curve is obtained. Otherwise, some parts of MVC should be constructed, like the case shown in Fig. 3.

$Step\ 6:$ Construct portions of MVC which are not covered by the previously generated solution curve, i.e. dashed curve $b_{2}b_{f}$. It is done by establishing polygons for such intervals using the algorithm proposed in Section 3.3.

$Step\ 7:$ By calculating $\ddot{s}/\dot{s}-k^{\prime}$, find a set of smooth critical points where the characteristics of the MVC constructed in $Step\ 6$ changes from sink to source, e.g. points $C_{4}$ and $C_{5}$.

$Step\ 8:$ Repeat $Step\ 5$ for the set of smooth critical points to obtain the last portions of the solution curve, i.e. $B_{3}$ and $B_{4}$ trajectories. After doing this, lowest portions of solution curves $B_{1}$ to $B_{n}$ constitute the solution curve, e.g. curve $D_{i}S_{1}C_{2}S_{3}C_{4}S_{5}D_{f}$ in Fig. 3.

A switching point located on MVC could either be a smooth critical point or zero-inertia point. Using the proposed algorithm, depending on the variety of switching points, three general cases may happen which we call them indirect, direct and semi-direct procedures. When there is no zero-inertia points on MVC, the procedure is indirect. In this case, all of the switching points are to be determined by constructing the portion of MVC which is not covered by initial and final solution curves generated in $Step\ 1$. The direct procedure happens when all of the switching points are of the zero-inertia type. In this case, the solution curve is calculated without the need to construct any part of MVC. Semi-direct procedure occurs when switching points are the combination of smooth critical points and zero-inertia points. In this case, the solution is obtained by constructing some portions of MVC which are not covered by the generated solution from zero-inertia points in $Step\ 5$. In the next section, to clarify the developed solution algorithm, two examples are given to represent two of these cases.

5 Numerical Example

In this section, two numerical examples are presented. We consider two different CMMS systems which one of them is planar and the other is three dimensional. The first example is an example which leads to a semi-direct procedure. The second example illustrates a case which leads to a completely indirect procedure to solve the time-optimal problem for 3D CMMS that clearly shows the effectiveness and efficiency of the proposed algorithm.

5.1 Example I: Semi-Direct Procedure

Fig. 4 shows the schematic of a CMMS composed of two planar manipulators handling a payload. Each manipulator has 3 DOFs which rigidly grasped a payload such that no slipping or rotation is possible at contact points. For this system, the values of $n$, $p$ and $m$ are 6,3 and 6, respectively. The physical properties of the system are listed in Table 1. The system is assumed to move the payload on a prescribed path defined as

[TABLE]

with $\dot{s}_{i}=4$ and $\dot{s}_{f}=4$.

The procedure begins by running SPA strategy in the forward direction by starting from the initial point, $D_{i}$(0, 4), which results branch $B_{i}$, as shown in Fig. 5a. The solution is continued by performing the same strategy in the backward direction by starting from the final point, $D_{f}$(1, 4), which makes branch $B_{f}$. Actuator’s bounds are violated at $b_{i}$(0.1434, 5.7960) and $b_{f}$(0.9301, 5.2179) that means the solution enters NFR at these points. By following the solution algorithm, the only zero-inertia point is calculated as $C_{2}$(0.8526, 4.1395). Then, SPA strategy is applied in the backward and forward direction by starting from $C_{2}$ that makes branches $B_{1}$ and $B_{2}$. When the branch $B_{1}$ reaches to $b_{1}$(0.3758, 6.7744) the actuator’s limits are broken. Also, branch $B_{2}$ crosses the curve $B_{f}$ at $S_{5}$(0.9630, 4.7920).

The solution curve is unknown between $b_{i}$ and $b_{1}$. Therefore, to complete the solution, a portion of MVC which is not covered by the generated solution curve is constructed. This part of MVC is shown with dashed curve in Fig. 5a. It should be noted that portions of MVC shown in Fig. 5a by dotted line are not necessary for calculating of the solution curve and are drawn only for better understanding.

Values of $\ddot{s}/\dot{s}$ and slope of MVC are compared to categorize points on MVC as sink or source. The behavior of the MVC changes smoothly from sink to source at $s_{c_{1}}=0.2672$. To obtain the solution curve, we execute SPA strategy by starting from $C_{1}$, forwardly with $\ddot{s}_{max}$ and backwardly with $\ddot{s}_{min}$ respectively, up to the point where generated solution curve in previous steps is crossed in $S_{1}$(0.0945, 5.2431) and $S_{3}$(0.4345, 6.3232). The solution curve has been shown in Fig. 5b. Also, the minimum time calculated for moving on this trajectory is $0.206$ sec.

The solution curve includes five switching points of $S_{1}$(0.0945, 5.2431), $S_{2}$(0.2672, 5.4703), $S_{3}$(0.4345, 6.3232), $S_{4}$(0.8526, 4.1395) and $S_{5}$(0.9630, 4.7920). Points $S_{1}$, $S_{3}$ and $S_{5}$ are located out of MVC and the acceleration is switched from the maximum value to its minimum. However, two other switching points, $S_{2}$ and $S_{4}$, are placed on MVC and the acceleration is switched from the minimum value to its maximum. At both critical points, the acceleration changes continuously from the minimum value to its maximum. However, the acceleration is discontinuous at the switching points which are located out of MVC. Fig. 6 shows time-discrete snapshots of robot motion. The angular position and velocity of the joints are shown in Fig. 7. Also, Fig. 8 shows the actuator torques of each robot arms while moving on the prescribed path with the minimum time solution. In terms of continuity, the actuator torques has the same behavior of the acceleration at the switching points. Also, the torque limitation has been fulfilled and on each point of the path, at least four actuators are always saturated.

The solution depends on the system and payload path whether it is obtained by a direct, semi-direct or indirect procedure. For instance, it can be easily check that the solution could be obtained directly, i.e. no portion of MVC is needed to construct, for payload path defined as

[TABLE]

with $\dot{s}_{i}=0$ and $\dot{s}_{f}=0$.

5.2 Example II: Indirect Procedure

The CMMS system considered at this example is composed of two PUMA 560 manipulators which are handling a payload, as shown in Fig. 9. Each manipulator has 6 DOFs which rigidly grasped a payload such that no slipping or rotation is possible at contact points. For this system, the values of $n$, $p$ and $m$ are 12, 6 and 12, respectively. The physical characteristics and Denavit-Hartenberg parameters of PUMA 560 are given in [24]. The system is assumed to move the payload on a prescribed path defined as

[TABLE]

with $\dot{s}_{i}=1$ and $\dot{s}_{f}=1$.

The procedure begins by running SPA strategy in the forward and backward direction from $D_{i}$(0, 1) and $D_{f}$(1, 1) respectively which results in branches $B_{i}$ and $B_{f}$, as shown in Fig. 10a. Actuator’s bounds have been violated at $b_{i}$(0.0564, 1.4694) and $b_{f}$(0.9488, 1.6534) that means the solution enters to NFR at these points. In the range of $b_{i}b_{f}$, there is no zero-inertia point of the sink-source type. Thereby to generate the solution, a portion of MVC placed between $b_{i}$ and $b_{f}$ is constructed. This part of MVC is shown with dashed curve in Fig. 10a.

The values of $\ddot{s}/\dot{s}$ and slope of MVC are compared to determine the smooth critical points. The behavior of the MVC changes smoothly from sink to source at points $s_{c_{1}}=0.0854$, $s_{c_{2}}=0.5001$ and $s_{c_{3}}=0.8922$. The solution curves are established by performing SPA strategy in the backward and forward direction by starting from points $C_{1}$, $C_{2}$ and $C_{3}$ respectively. The generated curves cover whole domain of $s$ from 0 to 1 and the solution algorithm ends. The solution curve consists of seven switching points $S_{1}$(0.0412, 1.3573), $S_{2}$(0.0854, 1.2833), $S_{3}$(0.2365, 1.7672), $S_{4}$(0.5001, 1.4047), $S_{5}$(0.7319, 1.8045), $S_{6}$(0.8922, 1.4868) and $S_{7}$(0.9720, 1.3734) as shown in Fig. 10b. The minimum time calculated, also, amounts to $0.674$ sec. The angular position and velocity of the joints are shown in Fig. 11 and Fig. 12, respectively. The actuator torques of each robot arms for the minimum time solution is shown in Fig. 13. The torque limitation has been fulfilled on each point of the path and at least seven actuators are always located in their saturation bound.

The importance of using an efficient algorithm to solve the minimum time problem is more clear in $Example\ II$. It can be investigated in two different aspects. First, the CMMS system which is composed of two PUMA 560 manipulators is more complicated than planar robot considered in the previous example. Second, for the path prescribed in (14), there is no zero-inertia point, therefore, to determine the critical points, a large portion of MVC should be constructed.

6 Efficiency Evaluation

In this section, the efficiency of the proposed algorithm is compared with one of the most recent method. We implemented the presented method in [22] for the numerical examples solved in section 5. Both methods are run on Matlab/Simulink 9.0 enabling parallel computing on a PC (Intel i7-5930 CPU and 24 GB Memory). The forward and backward integrations performed with 1 $ms$ as step time. The MVC is also constructed by considering 1000 points along the path. The computation time for calculating two major parts of the time-optimal path, i.e. MVC and solution curve, are given in Table 2 for two methods as well as total computation time. The results show that the computation time are reduced 62.2% and 56.8% for examples I and II, respectively. The time consumed to calculate the MVC for the examples is reduced which is due to skip of the whole construction of MVC by implementing our algorithm. Using the SPA algorithm, the computation time for calculating solution curve is substantially reduced by looking for saturated actuators in only a limited number of points, instead of determining the saturated actuators in each point which is performed in [22]. For example II, the equations of motion become more complicated which results a large computation cost to find the time-optimal path. By implementing the proposed method, the computation time for applying the method presented in [22], which was 2304 $sec$, is drastically reduced to 996 $sec$. This example clearly shows the efficiency of our algorithm by saving of 1308 $sec$.

We also implemented our method for the model of biped robot reported by Liu et al. [25]. They used parametric optimization method to find periodic steady state trajectory for a biped walking at a specified speed of 1.8 $km/h$, and tried to re-optimize the resulted trajectory by using dynamic differential programming method. As proposed in our previous study [15], we formulated the minimum time problem with stability and non-slip conditions along with actuator’s limits expressed as some inequality constraints. In addition, certain kinematic constraints are considered in terms of hip joint position that ensure an acceptable walking pattern. By solving this problem with the method presented in this paper, the maximum achievable speed is 5.6 $km/h$.

7 Conclusion

The problem of finding the time-optimal trajectory of a cooperative multi-manipulator system with redundant actuators moving on a specified path is studied in this paper. The focus of this study is to reduce the computational cost of the solutions obtained based on phase plane analysis. Two important issues which are sources of large computational costs were considered and resolved. The first problem is how to compute minimum or maximum acceleration along the path. Using the presented method, instead of determining the saturated actuators in each point which is computationally expensive, we can look for saturated actuators only in a limited number of points. The second problem is the calculation of the switching points which is the most difficult part of the time-optimal solution. This paper provides a comprehensive study on this subject and examines all possible situations. Several examples are provided to show different possibilities which may be encountered, and how they are solved by means of the proposed method.