Finding a Feasible Initial Solution for Flatness-Based Multi-Link Manipulator Motion Planning under State and Control Constraints

TL;DR

This paper introduces a method to generate feasible initial solutions for flatness-based multi-link manipulator motion planning problems with constraints, improving the reliability of solving non-convex optimization tasks.

Contribution

The paper proposes a novel approach to initialize feasible solutions by analyzing final time effects and scaling trajectories, enhancing motion planning under constraints.

Findings

Successfully generated feasible initial solutions for constrained manipulators

Integrated the method with inverse dynamics algorithms for practical testing

Demonstrated effectiveness through numerical experiments

Abstract

In this paper, we present a method to initialize at a feasible point and unfailingly solve a non-convex optimization problem in which a set-point motion is planned for a multi-link manipulator under state and control constraints. We construct an initial feasible solution by analyzing the final time effect for feasibility problems of flatness based motion planning problems. More specifically, we first find a feasible time-optimal trajectory under state constraints without a control constraint by solving a linear programming problem. Then, we find a feasible trajectory under control constraints by scaling the trajectory. To evaluate the practical applicability of the proposed method, we did numerical experiments to solve a multi-link manipulator motion planning problem by combining the method with recursive inverse dynamics algorithms.

Equations

$$\mathbf{y}(t)$$

$$\mathbf{x}(t)$$

$$\mathbf{u}(t)$$

$$\mathbf{\tau }=M(\mathbf{q})\ddot{\mathbf{q}}+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+\mathbf{G}(\mathbf{q})$$

$$\underset{\mathbf{q}(\cdot ),t_f}{\min }J(\mathbf{q}(\cdot ),t_f)$$

$$s.t.$$

$$\mathbf{q}(0)={\mathbf{q}}_0,\mathbf{q}(t_f)={\mathbf{q}}_f$$

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

$$\underline{\mathbf{q}}\le \mathbf{q}(t)\le \overline{\mathbf{q}}$$

$$\underline{\dot{\mathbf{q}}}\le \dot{\mathbf{q}}(t)\le \overline{\dot{\mathbf{q}}}$$

$$\underline{\mathbf{\tau }}\le \mathbf{\tau }(\mathbf{q}(t),\dot{\mathbf{q}}(t),\ddot{\mathbf{q}}(t))\le \overline{\mathbf{\tau }}(t\in [0,t_f])$$

$$where$$

$$\mathbf{\tau }(\mathbf{q},\dot{\mathbf{q}},\ddot{\mathbf{q}})=M(\mathbf{q})\ddot{\mathbf{q}}+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+\mathbf{G}(\mathbf{q})$$

$$\mathbf{q}(t)$$

$$\dot{\mathbf{q}}(t)$$

$$\ddot{\mathbf{q}}(t)$$

$$\tilde{\mathbf{\tau }}(s)$$

$$s$$

$$\dot{s}$$

$$\ddot{s}$$

$$\tilde{\mathbf{q}}(s)$$

$${\tilde{\mathbf{q}}}^{\prime }(s)$$

$${\tilde{\mathbf{q}}}^{\prime \prime }(s)$$

$$\tilde{\mathbf{\tau }}(s)$$

$$\underset{A\in {\mathbb{R}}^{n\times m},t_f\in \mathbb{R}}{\min }J(A,t_f)$$

$$s.t.$$

$$A\mathbf{b}(0)={\mathbf{q}}_0,A\mathbf{b}(1)={\mathbf{q}}_f$$

$$A{\mathbf{b}}^{\prime }(0)={\dot{\mathbf{q}}}_0{t}_f,A{\mathbf{b}}^{\prime }(1)={\dot{\mathbf{q}}}_f{t}_f$$

$$\underline{\mathbf{q}}\le A\mathbf{b}(s)\le \overline{\mathbf{q}}$$

$$t_f\underline{\dot{\mathbf{q}}}\le A{\mathbf{b}}^{\prime }(s)\le t_f\overline{\dot{\mathbf{q}}}$$

$$\underline{\mathbf{\tau }}\le \tilde{\mathbf{\tau }}(s;A,t_f)\le \overline{\mathbf{\tau }}(s\in [0,1])$$

$$t_f\ge 0$$

$$where$$

$$\tilde{\mathbf{\tau }}(s;A,t_f)=[M(\tilde{\mathbf{q}}){\tilde{\mathbf{q}}}^{\prime \prime }+C(\tilde{\mathbf{q}},{\tilde{\mathbf{q}}}^{\prime }){\tilde{\mathbf{q}}}^{\prime }]t_f^{-2}+\mathbf{G}(\tilde{\mathbf{q}})$$

$$\{A^{LP},t_f^{LP}\}=\underset{A\in {\mathbb{R}}^{n\times m},t_f\in \mathbb{R}}{\mathrm{argmin}}{t}_f$$

$$s.t.$$

$$A\mathbf{b}(0)={\mathbf{q}}_0,A\mathbf{b}(1)={\mathbf{q}}_f$$

$$A{\mathbf{b}}^{\prime }(0)={\dot{\mathbf{q}}}_0{t}_f,A{\mathbf{b}}^{\prime }(1)={\dot{\mathbf{q}}}_f{t}_f$$

$$\underline{\mathbf{q}}\le A\mathbf{b}(s_i)\le \overline{\mathbf{q}}$$

$$t_f\underline{\dot{\mathbf{q}}}\le A{\mathbf{b}}^{\prime }(s_i)\le t_f\overline{\dot{\mathbf{q}}}(i\in \{1,2,...,N\})$$

$$t_f\ge 0$$

$$\underset{t_f\to \infty }{\lim }\tilde{\mathbf{\tau }}(s;A,t_f)=\mathbf{G}(A\mathbf{b}(s))$$

$$b_{i,1}(s)=\{\begin{array}{ll}1& \text{if si≤s<si+1}\\ 0& \text{otherwise}\end{array}$$

$$b_{i,k}(s)=$$

$$\frac{s-s_i}{s_{i+k-1}-s_i}{b}_{i,k-1}(s)+\frac{s_{i+k}-s}{s_{i+k}-s_{i+1}}{b}_{i+1,k-1}(s)$$

$$\begin{gathered} b_{i,k}^{\prime }(s)=(k-1)\times \\ \left[\frac{1}{s_{i+k-1}-s_i}{b}_{i,k-1}(s)-\frac{1}{s_{i+k}-s_{i+1}}{b}_{i+1,k-1}(s)\right] \end{gathered}$$

$$\begin{gathered} \underline{\mathbf{q}}{1}_m^{\top }\underset{i,j}{\le }A\underset{i,j}{\le }\overline{\mathbf{q}}{1}_m^{\top }\Rightarrow \\ \underline{\mathbf{q}}\le A\mathbf{b}(s)\le \overline{\mathbf{q}}\forall s\in [0,1] \end{gathered}$$

$$\begin{gathered} \frac{1}{k-1}{t}_f\underline{\dot{\mathbf{q}}}({\mathbf{v}}_{k+1:k+m-1}-{\mathbf{v}}_{2:m}{)}^{\top } \\ \underset{i,j}{\le }A\left(\left(\begin{array}{c}{I}_{m-1\times m-1}\\ 0_{m-1}^{\top }\end{array}\right)-\left(\begin{array}{c}{0}_{m-1}^{\top }\\ I_{m-1\times m-1}\end{array}\right)\right) \\ \underset{i,j}{\le }\frac{1}{k-1}{t}_f\overline{\dot{\mathbf{q}}}({\mathbf{v}}_{k+1:k+m-1}-{\mathbf{v}}_{2:m}{)}^{\top } \\ \Rightarrow \underline{\dot{\mathbf{q}}}\le A\dot{\mathbf{b}}(s)\le \overline{\dot{\mathbf{q}}}\forall s\in [0,1] \end{gathered}$$

Full text

1 Introduction

To make a highly automated system, motion planning is more important than following a reference trajectory given by humans. To plan the motion of robotic manipulators, the motion can be ineffective and sometimes physically impossible if we inexactly deal with nonlinear dynamics of the manipulators and constraints like joint angle and torque limit. In this paper, we tackle a continuous-time set-point motion planning problem of multi-link manipulators under state and control constraints by using nonlinear dynamics exactly.

A typical method for motion planning that consider dynamics under some constraints discretizes state equations by using small time grids and solves their constrained optimization problems (e.g. Multiple Shooting Methods). However, such methods suffer from huge numbers of variables and numerical integration errors when the methods are used to try to solve long time horizon problems. We adopt flatness based motion planning methods because they can formulate long time horizon problems that have few parameters while considering nonlinear dynamics precisely and avoiding numerical integration errors.

In planning using a flatness based method under control constraints, its optimization problem generally becomes non-convex due to nonlinear dynamics. Therefore, if we solve the problem with a general nonlinear programming (NLP) algorithm like sequential quadratic programming (SQP) or interior point methods, the behavior of the algorithm varies depending on the initialization. Thus, there is no guarantee that a feasible solution can be found even if the problem has feasible regions. [7] [2]. (In our numerical experiments of manipulator planning problems, NLP algorithms often failed to find a feasible solution when we tried to use several initial points.)

We aim to reveal details of the optimization problem structure of the flatness based motion planning algorithm under state and control constraints and construct a method that can unfailingly find feasible solutions of it.

We derive a method to initialize a solution algorithm with a feasible initial point by parameterizing a trajectory with a linear combination of basis functions and by time transformation. Our contributions are summarized as follows.

• •

By combining time scaling transformation with parameterization of a linear combination of an arbitrary basis function, we convert the planning problem to analyze the effect of final time on the optimization problem structure.

• •

We describe a time-optimal feasible solution under state constraints that uses linear programming(LP) and use it to construct a method to find feasible initial points under control constraints.

• •

We present details of practical implementations to solve multi-link manipulator planning problems with numerical simulation by combining our initialization method and recursive inverse dynamics algorithms.

2 Related Work

A method that deals with a state equation constraint by locating sample times in the control horizon and calculating the states by numerical integration is called a “shooting method” [1]. There are many types of research that extend this family of methods and reduce their computational cost. We focus on a method that does not use numerical integration of a state equation because a state equation of a Multi-Link robot manipulator is difficult to numerically integrate both very accurately and quickly.

The differential flatness concept was proposed by Fliess et al.[5]. Flatness has been studied from the viewpoint of motion planning problems [12]. Louembet et al.[10] developed a convex programming algorithm for a motion planning problem under semi-infinite inequality constraints by approximating a feasible region to a convex polytope. Van Looke et al.[18] developed a method to approximate polynomial semi-infinite inequality constraints by using the convex hull property of B-splines [3]. Furthermore, they applied the methods to a flexible robot arm planning problem [17]. We follow this line of research further and tackle more practical problem settings on the basis of more complex dynamics.

In robotics, Martin and Bobrow [11], Wang et al.[19], and Lee et al.[9] have solved motion planning problems of multi-joint robot manipulators by parameterizing joint trajectories with B-splines, though they did not use the flatness concept explicitly. The theoretical basis of these methods becomes clearer from the viewpoint of flatness with our formulation, and their formulation can be basically understood to be the same as our formulation of flatness-based motion planning.

Their methods are based on the Lie group formulation of multi-link robot manipulator modeling and effective calculation of the inverse dynamics of a manipulator. The theory of the Lie group formulation of robot dynamics is described in detail elsewhere [13], [14]. A recursive algorithm for the Lie group formulation of robot dynamics was developed [14] and then analyzed more deeply and extended [15] [16]. Our method references Kim and Polland’s implementation [8].

Wang et al. [19] solved a problem by formulating control constraints as costs, which are soft constraints. We deal with control constraints as hard constraints to avoid the weight of control constraints, which vary depending on cost weight tuning.

3 Problem Formulation

In this section, we provide a basic formulation of a flatness based motion planning of robotic manipulators. In this paper, all vectors are column vectors and ${\top}$ means vector or matrix transposition.

3.1 Flatness Based Motion Planning of Robot Manipulator

A nonlinear system $\dot{\bm{x}}(t)=\bm{f}(\bm{x}(t),\bm{u}(t))$ is differentially flat or flat if there exists one-to-one mapping from a tuple of state and control $(\bm{x}(t),\bm{u}(t))$ to flat outputs $\bm{y}(t)$ that satisfy the below equations.

Definition 1** (Differential Flatness [5]).**

[TABLE]

The equation of the motion of a mechanical system, like that of a robot manipulator, is generally formulated as follows [13]. (The friction term is omitted for the sake of simplicity).

[TABLE]

where $\bm{q}(t)\in{\mathbb{R}}^{n}$ are joint trajectories, $M(\bm{q})\in{\mathbb{R}}^{n\times n}$ is a generalized inertia matrix, $C(\bm{q},\dot{\bm{q}})\in{\mathbb{R}}^{n\times n}$ is a Coriolis matrix, $\bm{G}(\bm{q})\in{\mathbb{R}}^{n}$ is a potential term, and $\bm{\tau}\in{\mathbb{R}}^{n}$ are torque control inputs. Robot manipulators are flat systems if we regard joint trajectories $\bm{q}(t)$ as flat outputs. In other words, $\bm{q}(t),\dot{\bm{q}}(t),\ddot{\bm{q}}(t)$ and $\bm{\tau}(t)$ can be derived only by the $\bm{q}(t)$ trajectory. Therefore, we can plan both state and control trajectories of a manipulator by solving a function optimization problem only of a joint trajectory $\bm{q}(t)$ that is second order differentiable. In this paper, we focus on the following a set-point problem under linear state and control constraints.

Problem 1** (Flatness Based Robot Motion Planning Problem).**

[TABLE]

$\underline{\bm{q}},\overline{\bm{q}},\underline{\dot{\bm{q}}},\overline{\dot{\bm{q}}}\in{\mathbb{R}}^{n},\quad\underline{\bm{\tau}},\overline{\bm{\tau}}\in{\mathbb{R}}^{n}$ are lower or upper bounds constants with respect to joint angles, joint angular velocities, and input torques of each joint. The cost function $J$ is specified in accordance with task settings. $t_{f}\in{\mathbb{R}_{+}}$ is the final time of the control horizon, We formulate both free and fixed final time problems in a unified way by assuming $t_{f}$ of Problem 1 as a variable or a constant.

3.2 Time Scaling Transformation

To analyze our optimization problem structure of both free and fixed final time and separate parameterization of a trajectory and final time, we normalize control time horizon $t\in[0,t_{f}]$ to $s\in[0,1]$ by time scaling transformation. Then our trajectories of joint angle, angular velocity, angular acceleration, and torque are formalized as follows.

[TABLE]

In this paper, we write the time normalized function of variable $s$ as a tilde over a character like “$\;\tilde{}\;$” and “ ${}^{\prime}\;$” for derivatives of $s$ .

3.3 Trajectory Parameterization with Linear Combination of Basis Functions

We parameterize $\tilde{q}_{i}(s)$ as a linear combination of $m$ basis functions $b_{j}(s)$ as $\tilde{q}_{i}(s)=\sum_{j=1}^{m}a_{ij}b_{j}(s)$ and $s(t)$ as linear scaling of $t$ to represent a large class of trajectory and to parameterize simply. Then, we can derive the following trajectory formulation.

[TABLE]

where $\bm{a}_{j}\in{\mathbb{R}}^{n}$, $A=[\bm{a}_{1},\bm{a}_{2},\cdots,\bm{a}_{m}]\in{\mathbb{R}}^{n\times m}$, and $\bm{b}(s)\in{\mathbb{R}}^{m}$ . Time scaled and trajectory parameterized Problem 1 is formulated as follows.

Problem 2** (Time Scaled Motion Planning Problem).**

[TABLE]

The main challenge in this paper is this formulation of a problem. Problem 2 has the following difficulty.

• •

Inequality constraints are semi-infinite constraints on a continuous interval of $s$.

• •

Torque $\tilde{\bm{\tau}}(s;A,t_{f})$ constraint is generally nonlinear and non-convex.

We can relax the former problem by evaluating finite sample points on $s_{i}\in[0,1]$ . The latter non-convex problem is the fundamental difficulty of our optimization problem, and in our numerical experiments, we often failed to solve the problem by directory applying a standard nonlinear programming(NLP) algorithm like SQP and interior point methods.

Problem 2 has no obvious feasible points and if we solve the problem by using NLP with a infeasible initial point, the algorithm may stick at an infeasible point near it, and there is no guarantee to find a feasible point even if the problem has a feasible region.

4 Initialization Method with a Feasible Initial Point

To unfailingly find a solution for our non-convex Problem 2 in feasible settings, we think that one reliable solution is initializing with a feasible point. To construct the feasible point, first we solve a convex sub-problem that satisfies state constraints and then modify its solution to satisfy control constraints. Our method is composed of three steps and overviewed as Algorithm 1. We next describe details of the algorithm and how the algorithm can find a feasible point.

4.1 Time-Optimal Feasible Point under State Constraints

We first ignore torque constraints and solve the time-optimal trajectory under state constraints that are relaxed to evaluate finite sample points in control intervals. This sub-problem can be effectively solved with linear programming(LP).

Problem 3** (Time-Optimal Motion Planning Problem under State Constraints (Linear Programming)).**

[TABLE]

4.2 Feasible Initial Point under State and Control Constraints

Next, we construct a feasible initial point under control constraints without violating state constraints using $\{A^{LP},t_{f}^{LP}\}$. We use the following property.

Property 1** (Relaxed Torque Constraints Check).**

If $\dot{\bm{q}}_{0}$ and $\dot{\bm{q}}_{f}=0$ and whole trajectory of the solution of Problem 3: $\tilde{\bm{q}}(s)=A^{LP}\tilde{\bm{b}}(s)$ satisfies $\underline{\bm{\tau}}<\bm{G}(A^{LP}\bm{b}(s))<\overline{\bm{\tau}}$ , there exists $t_{f}>=t_{f}^{LP}$ that satisfies torque constraints.

Proof.

Using torque expression in (14)

[TABLE]

Therefore, if $t_{f}$ is large, the torque constraints are $\underline{\bm{\tau}}\leq\epsilon+\bm{G}(A\bm{b}(s))\leq\overline{\bm{\tau}}$ for $|\epsilon|\ll 1$ . This inequality is satisfied if $\underline{\bm{\tau}}<\bm{G}(A^{LP}\bm{b}(s))<\overline{\bm{\tau}}$ is always satisfied for $\forall{s}\in[0,1]$ . If we make $t_{f}$ larger, the range of constraints $t_{f}\underline{\dot{\bm{q}}}\leq A\bm{b}^{\prime}(s_{i})\leq t_{f}\overline{\dot{\bm{q}}}$ expands and is not violated. ∎

Practically, we find suitable $t_{f}^{feas}$ by a line search in the second step of the algorithm. In the last step, we apply a NLP algorithm using this feasible initial solution $\{A^{LP},t_{f}^{feas}\}$ and improve the solution.

Remark 1** (About Assumption in Property 1).**

Assumption $\dot{\bm{q}}_{0}$ and $\dot{\bm{q}}_{f}=0$ is satisfied if the motion is from a stop state to a stop state, which is the case in a practical situation. Assumption $\underline{\bm{\tau}}<\bm{G}(A^{LP}\bm{b}(s))<\overline{\bm{\tau}}$ is always satisfied in practical task settings for an industrial manipulator. Thus, we think these assumptions cover a large class of practical applications. A general case is an interesting problem but is future work.

We summarize Algorithm 1 and call the solutions of each step in the algorithm as follows.

• •

solution LP: $\bm{q}(t)^{LP}=A^{LP}\bm{b}((t_{f}^{LP})^{-1}t)$

a time-optimal solution that satisfies state constraints solved with LP

• •

**solution feas:111Solution feas is made by expanding solution LP, thus variable A is not changed. ** $\bm{q}(t)^{feas}=A^{LP}\bm{b}((t_{f}^{feas})^{-1}t)$

a solution that satisfies both state and control constraints

• •

solution opt: $\bm{q}(t)^{opt}=A^{opt}\bm{b}((t_{f}^{opt})^{-1}t)$

a solution solution feas improved with NLP by minimizing cost in accordance with task settings

We have described the main procedure of our algorithm. However, to implement our algorithm for practical multi-link manipulators, inequality constraints need to be efficiently checked. We next describe the techniques to do this.

5 Practical Implementation Techniques to Calculate Torque of Multi-Link Manipulators

To check torque constraints at normalized time $s_{i}$, we have to calculate $\tilde{\bm{\tau}}(\tilde{\bm{q}}(s_{i}),\tilde{\bm{q}}^{\prime}(s_{i}),\tilde{\bm{q}}^{\prime\prime}(s_{i}),t_{f})$ by using an inverse dynamics algorithm. The problem here is the very high computational cost to calculate inverse dynamics for a multi-link robot manipulator such as a 6-degrees of freedom (DoF) one because a motion equation of the manipulator is extremely complex and almost impossible to write down in closed form explicitly. To calculate torque efficiently, we use a recursive inverse dynamics algorithm formulated by a Lie group.

5.1 Lie Group Formulation of Recursive Inverse Dynamics

An equation of motion of a multi-link robot manipulator obtains a high-level, coordinate-free view when it is formulated in a Lie group [13], [14]. We use a Lie group formulation of robot dynamics because it is useful to analyze and improve from a geometric viewpoint. An order of the computational cost of calculating inverse dynamics by using a closed-form equation of motion is $O(n^{4})$, where $n$ is the number of joints, which is very large [4][15]. We use a recursive calculation algorithm for this, which is known to have a low computational cost at order $O(n)$ .

In the recursive algorithm, an equation of motion is divided into recursive equations by introducing intermediate variables. First, the algorithm computes states of joints from a base of a manipulator to the end-effector by forward recursion. Then, torque values for the states are calculated from the end-effector to the base by backward recursion.

In this paper, the algorithm is implemented with reference to the work of Kim and Pollard [8] and modified with time scaling as Algorithm 2.

$\bm{\xi}=(\bm{\omega},\bm{v})\in{\mathbb{R}}^{6},\bm{V}_{i}\in{\mathbb{R}}^{6}$ are generalized velocities (or twist) and $\bm{\xi}$ is generalized velocity of joint frame at a reference frame, $\bm{F}_{i}\in{\mathbb{R}}^{6}$ are generalized forces, $J_{i}\in{\mathbb{R}}^{6\times 6}$ is generalized inertia matrix, $\bm{g}\in{\mathbb{R}}^{3}$ is a gravitational constant vector. $\hat{\cdot}:{\mathbb{R}}^{6}\rightarrow se(3)$ is defined as $\hat{\bm{\xi}}=\begin{bmatrix}\hat{\bm{\omega}}&\bm{v}\\ \bm{0}^{\top}&0\\ \end{bmatrix}$ where $\hat{\bm{\omega}}=\begin{bmatrix}0&-{\omega}_{3}&{\omega}_{2}\\ {\omega}_{3}&0&-{\omega}_{1}\\ -{\omega}_{2}&{\omega}_{1}&0\end{bmatrix}$. For $se(3)$, an adjoint mapping $Ad:SE(3)\times se(3)\rightarrow se(3)$ is defined as $Ad_{G}(\hat{\bm{\xi}})=G\hat{\bm{\xi}}G^{-1}$ , and a Lie bracket $ad:se(3)\times se(3)\rightarrow se(3)$ is defined as $ad_{\hat{\bm{\xi}}_{1}}(\hat{\bm{\xi}}_{2})=\hat{\bm{\xi}}_{1}\hat{\bm{\xi}}_{2}-\hat{\bm{\xi}}_{2}\hat{\bm{\xi}}_{1}$ . If $\hat{\bm{\xi}}\in se(3)$ is represented as a vector $\bm{\xi}\in{\mathbb{R}}^{6}$ , $Ad$ and $ad$ are represented in matrix form with using $G=\begin{bmatrix}R&\bm{p}\\ \bm{0}^{\top}&1\\ \end{bmatrix}$ as $Ad_{G}(\bm{\xi})=\begin{bmatrix}R&0_{3\times 3}\\ \hat{\bm{p}}R&R\\ \end{bmatrix}\bm{\xi}$, $ad_{\bm{\xi}_{1}}(\bm{\xi}_{2})=\begin{bmatrix}\hat{\bm{\omega}}_{1}&0_{3\times 3}\\ \hat{\bm{v}}_{1}&\hat{\bm{\omega}}_{1}\\ \end{bmatrix}\bm{\xi}_{2}$ . See details in the work of Lee et al. [9] and Park et al. [14].

5.2 Efficient Inequality Constraint Checking Method for B-Spline Basis

In our method, we can choose an arbitrary basis function such as a polynomial, orthogonal polynomial, wavelet, spline, etc. If we adopt the B-Spline basis function, we can check semi-infinite linear inequality constraints without finite dense sampling. We describe the efficient inequality checking method applicable for B-Spline. In our proposed method, the B-spline order $k\in\mathbb{N}_{+}$ and knot vector $\bm{v}=(0=s_{1},s_{2},\cdots,s_{m+k}=1)^{\top}\in{\mathbb{R}}^{m+k}$ are given by users. The B-spline basis of order $k$ is as follows.

Definition 2** (B-Spline Basis of Order $k$).**

[TABLE]

The joint angular velocity is $\tilde{\dot{\bm{q}}}(s)=A\bm{b}^{\prime}(s)t_{f}^{-1}$, where a derivative of B-spline basis $b^{\prime}_{i}(s)$ is as follows [3].

[TABLE]

Since a B-spline has a convex hull property [3][18], if the convex hull is in a feasible region, the feasible region always includes the B-Spline’s trajectory. Thus, the (finite number of) inequality constraints below is a sufficient condition for semi-infinite constraints (15)(16)[18].

Property 2** (Inequality of $\tilde{\bm{q}}(s)=A\bm{b}(s)$ Linear Constraints for B-Spline Parameterization).**

[TABLE]

where $\bm{1}_{m}=(\underbrace{1,\cdots,1}_{m})^{\top}\in{\mathbb{R}}^{m}$ .

Property 3** (Inequality of $\tilde{\bm{q}}^{\prime}(s)t_{f}=A\bm{b}^{\prime}(s)$ Linear Constraints for B-Spline Parameterization).**

[TABLE]

where $\underset{i,j}{\leq}$ is elementwise inequality of a matrix, $\bm{v}_{i:j}$ means a vector part of indices $\{i,\cdots,j\}$ 222Matlab like notation, $I_{m-1\times m-1}\in{\mathbb{R}}^{m-1\times m-1}$ is an identity matrix, and $\bm{0}_{m-1}=(\underbrace{0,\cdots,0}_{m-1})^{\top}\in{\mathbb{R}}^{m-1}$ . Semi-infinite inequality for the whole continuous range $[0,1]$ can be checked by using the above linear inequality of $A$ . Note that the above property is a sufficient condition for checking (15)(16). Thus, this method could produce conservative constraints [18].

6 Numerical Experiments

We conducted numerical experiment to assess whether our algorithm can be calculated correctly with practical computational time and generate an appropriate solution. First, we evaluated basic behavior by applying the algorithm to a 2-link manipulator. We next applied the algorithm to a practical 6-DoF manipulator that has many links to assess the scalability of our algorithm. The algorithm is implemented by Matlab and optimization problems are solved by fmincon with SQP. The Derivative of torque constraints, which is a very complicated expression, is calculated by using an automatic differentiation [6] module we made. We checked correctness of inverse dynamics implementation and visualized result motions by using simscape multibody.

6.1 In the case of a 2-Link Manipulator

Initial states, final states, and the planned motion result are shown in Fig. 2. Problem settings are as follows.

• •

Cost function $J(A,t_{f})=t_{f}$, thus the problem is a time-optimal set-point problem.

• •

State and control constraints are

[TABLE]

(torque limit satisfies the assumption of enough force for the gravity term)

• •

Basis function is a polynomial basis, and order is 10.

• •

The number of sample points $N$ for checking inequality is 24.

The planned trajectory details are shown in Fig.1. The implemented algorithm worked correctly, and observed behaviors were as expected.

• •

Solution LP satisfies all state constraints but violates torque constraints.

• •

Solution feas satisfies all state and torque constraints and is a feasible but not optimal solution.

• •

Solution opt satisfies all state and torque constraints and lowers cost $t_{f}$.

We actually made torque constraints slightly stricter in step 2 of the algorithm because we observed some cases in which if solution feas is very close to the torque constraints, the NLP algorithm can barely improve the solution.

6.2 In the case of a 6-DoF Manipulator

Initial states, final states, and the planned motion result are shown in Fig. 4. The problem settings are as follows.

• •

Cost function $J(A,t_{f})=t_{f}$

• •

The manipulator has 6-links and 6-DoFs as shown in Fig. 3

• •

Basis function is B-spline, order is $k=9$, and the number of knots is $24$

Both trajectories of the states and controls of the result successfully satisfied the inequality constraints for them. We show some torque trajectories that move largely in Fig.5.

7 Discussion

We found that our method can correctly construct feasible initial points and prevent failure to find a feasible solution with NLP. The difference between our method and using multiple heuristic initial points is its reliability: our method can find a feasible solution for our non-convex problem without fail.

However, there are several problems that our method does not solve.

The first problem is that improved solution opt is local minima, and how much global optimum and solution opt differ is not clear. An important research direction is to find solutions nearer to the global optimum. For example, one idea is apply a tight convex relaxation to torque expression.

The second unsolved problem is which basis function is best for our problem. We tried polynomial basis, B-spline basis, and orthogonal polynomial basis (e.g. Chebyshev polynomial, Legendre polynomial). In our experiments, polynomial and B-Spline bases generate better results, but their mechanisms are not clear. An interesting direction is to use another type of basis such as wavelet spline that has both short and wide region base.

The third problem is that our algorithm correctly works under the assumption in Property 1, but we do not know whether it works correctly under harder constraints. Thus, the assumption should be removed and the algorithm should be applied under harder constraints. In our preliminary experiment under harder torque constraints, we observed an energy storing motion that first uses back swing to store motion and next achieves a final state that cannot be reached directly. Researching smart motions under harder constraints is an interesting direction.

8 Conclusion and Future Work

We presented a method for unfailingly finding a feasible solution for a flatness-based robot manipulator motion planning problem under state and control constraints when it is a non-convex optimization problem. Our method can reliably find a feasible solution under practical assumptions even though common methods that initialize multiple infeasible initial points and use general nonlinear programming (NLP) often fail and cannot reliably find a feasible solution.

Our method first finds a solution that satisfies state constraints, transforms it to satisfy control constraints, and lowers its cost. We found that our method correctly works for a multi-link robot manipulator that has two links or six-degrees of freedom (DoFs). Future directions are to find the best basis function, find a solution nearer to the global optimum, and construct a reliable method under constraints harder than those assumed in this paper.

9 Acknowledgments

The authors gratefully acknowledge the contributions of the reviewers.