Inverse kinematics of asymmetric octahedral variable geometry truss manipulator with obstacle avoidance through inexact interior point optimization

Abstract

This article studies the inverse kinematics for asymmetric octahedral variable geometry truss manipulator with obstacle avoidance. The inverse kinematics problem is cast as a nonconvex optimization that having quadratic objective function subject to quadratic constraints. This article uses an inexact interior point optimization to solve it, which is developed on the basis of the imprecise algorithm Ipopt. According to the particularity of our actual optimization problem, each iteration undergoes specific modifications so as to minimize the memory consumption as well as computation time. Utilizing the sparse and binary characteristics of the coefficient matrix, respectively, the algorithm allocates the computation to the finite sparse matrix vector multiplication and changes the storage form, which greatly reduces the memory space. Based on the unique rules of inverse kinematics, the iteration direction of the algorithm becomes more clear. With the aid of mechanical constraints inherent in the manipulator, the algorithm omits the feasibility recovery part that embedded in the solver Ipopt. All these make us save the operation time greatly while utilizing Ipopt algorithm. To demonstrate the effectiveness of the proposed approach, the scheme was applied to obstacle avoidance inverse kinematics of variable geometry truss manipulator with three modules.

Keywords

Variable geometry truss manipulator inverse kinematics obstacle avoidance interior point optimization

Introduction

The variable geometry truss manipulator (VGTM) is a type of hyper-redundant robotic mechanism. It is a statically determinate truss-type structure whose length of the actuators can be changed. The robotic arm studied in this article consists of a number of asymmetrical octahedral modules. Each module uses three SANMOTION F2 2-phase stepper motors as the power drive. By adjusting the length of the linear actuator, the robot arm can perform the specified operation. Due to its advantage of low weight and high stiffness, it is anticipated to be applied in the fields of future space station, submarine exploration, and detection.¹

As early as 1980s, the concept of VGT evolved from the adaptive structure concept.² Since then, many scholars have devoted themselves to the study of kinematics,^3,4 analysis of dynamics and control,⁵ design of conceptual and actual structures of VGTM,^6,7 workspace analysis,^8,9 and so on. Concerning nonsymmetric VGTM, the number of solutions to its inverse kinematics is infinite. In other words, for a preset end effector position, the configuration of VGTM may change with the end effector position remaining constant. Consequently, it is essential to choose an effective selection criteria.

In general, there are three main methods of inverse kinematics for manipulator: the geometric method,¹⁰ the numerical approach,^11,12 and the analytical method.^13,14 Samer Yahya et al. proposed a new geometric approach to solve the inverse kinematics of hyper-redundant manipulators.¹⁰ This method is accurate and fast and can determine all possible solutions. But it needs to decompose the spatial geometric parameters of the manipulator into planar ones. If a closed solution is to be obtained, the manipulator must satisfy a certain geometric structure, and the closed-form solution of this class of the robotic arm cannot be used for other manipulators with different geometries. Therefore, the geometric approach is suitable for solving inverse kinematics of robots with simple joint structure model. For a robotic arm with complex structure, the numerical approach can also obtain one solution from the infinitely multiple inverse solutions of the manipulator through iterative calculation.¹¹ However, the long computation time makes it unsuitable for online use, meanwhile, the algorithm usually does not have global convergence, and it is prone to singularity.¹² The analytical method has the advantages of accurate calculation results and can find all solutions. However, it requires a large amount of algebra and matrix operations, and the derivation process is complicated.¹³ Moreover, analytical solutions are difficult to obtain due to the multiparameters, nonlinearity and coupling of the solution and the complex nonlinear constrained algebraic equations involved in the inverse kinematics of the robot.¹⁴

For the high redundancy of VGTM, many studies have attempted to decompose it, which captured the characteristics of macroscopic geometry of VGTM through the use of “reference spline” or “backbone curve”. Zanganeh and Angeles¹⁵ put forward a solution method based on spline. For a prescribed optimization criterion, redundancy can be resolved by looking for the optimum shape of the skeleton curve. However, this approach requires a priori assumption that, along any one of the candidate curves, there exists a fixed coordinate system distribution of intermediate platforms. Chirikjian and Burdick^16,17 proposed two approaches to determining an optimal configuration of VGTM. In the first proposed method of modal hyper-redundancy resolution, the shape function of the backbone curve was restrained to only one set of modes. After that, they put forward the other method on basis of a continuum model for kinematics. Since this method did not provide a direct way of satisfying the constraints associated with the pose of the end effector, the modes needed to be selected carefully. In addition, the complex nature of this approach makes it rather costly to solve the inverse kinematics in real time for hyper-redundant manipulators in three-dimensional space. Fahimi et al.¹⁸ extended the modal approach. By changing the original mode using a new form of function, they applied the expanded method to solve the inverse kinematics of spatial hyper-redundant manipulators. However, this method is only applicable to the case that the hyper-redundant mechanical arm is a serial chain.

In recent years, on the other hand, the research on VGTM mainly focused on two major types of structural forms. One type is manipulators driven by binary actuating mechanisms with two steady states,^19,20 and the other type is those whose actuators are continuous but not installed on moving platforms.^21
–23 For the first type of VGTM, Bayram and Ozgoren^19,20 studied the concept design together with the direct kinematics along with position control, respectively, through the use of inverse kinematics of a binary hyper-redundant mechanical arm in space. For VGTMs with the second structure, Iwatsuki et al.²¹ analyzed the kinematics of a hyper-redundant mechanical arm and proposed a motion control method for it. Mintenbeck and Estanna²² investigated the configuration design, physical modeling and motion control of a hyper-constrained 3-revolute prismatic spherical (RPS) parallel robot. Chibani et al.²³ presented a deterministic optimizing process to generate optimal reference kinematic configurations for a spatial hyper-redundant parallel robot with multiple rigid modules. Aviles et al.²⁴ described two procedures to identify the optimum damper locations in VGTMs, and then they developed a novel method of analyzing the inverse dynamics of a VGTM containing elastic elements.²⁵ Rost et al.²⁶ presented the joint design and hydraulic drives for a VGTM consisting of three-degree-of-freedom octahedral modules.

However, for the VGTM constituted by asymmetric octahedral modules that to be studied, it is almost impossible to create a fixed coordinate system for each module because its actuators are mounted on a mobile platform. Therefore, it is extremely difficult to formulate the forward kinematics with explicit equations, which makes the methods used in the above two structures not applicable for nonsymmetric octahedral VGTMs.

In the nearly 30 years since the interior point algorithm has been invented, it has become the kernel of many existing nonlinear optimization solvers.²⁷ Due to its superlinear convergence as well as possible combination with inexact Newton’s method and powerful linear solver, the algorithm can even solve general nonconvex nonlinear optimization problems in an ideal time.²⁸ In addition, because the dominant computing cost in the interior point algorithm is to solve the linear equation,²⁹ its parallel computation and compression algorithms also play a significant role in quickly solving practical problems. However, few articles focus on applying the interior point method to inverse kinematics optimization. Based on the framework of the interior point algorithm, Vigerske Stefan and Andreas Wächter developed Ipopt as a solver for solving large-scale sparse problems.³⁰ Before solving the optimization problem through the Ipopt solver, it is necessary to provide the solver sparseness to speed up the computing speed. While a large number of large-scale sparse matrices exist in the optimized model established in this article, which is just right for solving our actual problems.

In this article, on basis of the inherent constraints existing in the VGTM structure, first, the inverse kinematics optimization model is established as a general nonconvex quadratic constrained quadratic programming (QCQP) problem. And then, according to the unique characteristics of our actual optimization model, the Ipopt algorithm is adjusted correspondingly to solve this specific optimization model. The sparsity and binarity of the coefficient matrix in our actual optimization model greatly reduce the storage space of the algorithm. The fact that the inverse kinematics follows certain observable rules makes the iteration direction clearer. The fixed mechanical structure constraints existing in VGTM can make the algorithm omit the feasibility recovery part. All of these have saved us a great amount of computation time in the solving process. To demonstrate the effectiveness of the presented approach, the adjusted algorithm was applied to obstacle avoidance inverse kinematics of an asymmetrical octahedral VGTM with three modules.

The remainder of this article is arranged as follows: “Configuration of VGTM” section introduces the kinematic constraints that exist in VGTM. In “Inverse kinematics” section, firstly, the inverse kinematics optimization of VGTM with n modules is formulized as a generalized nonconvex QCQP problem. And then based on inexact interior point optimization, the inverse kinematics optimization gets solved. After that, the numerical simulation results are provided in “Numerical simulations for inverse kinematics” Section. Finally, the conclusions are given in “Conclusion” section.

Configuration of VGTM

This section presents the kinematics restriction exists in an asymmetric octahedral VGTM based on its specific geometric construction. It will be divided into two parts and introduced separately: geometric construction and kinematical constraints of a single VGTM module.

Geometric construction of a single module

The geometric construction of the i th VGTM module is displayed in Figure 1. There are three independent linear actuators at the bottom and top of the octahedral module respectively. As long as the length of any actuator is changed, the platform on top of the module will move relative to that at the bottom. Let $Q_{i - 1, j}$ and $Q_{i, j} (i = 1, 2, \dots, n, j = 1, 2, 3)$ denote the node positions of this VGTM module. Referring to Figure 1, the geometric relationship contains six inequality constraints and six equality constraints.

Figure 1.

Geometric construction of the i th VGTM module. VGTM: variable geometry truss manipulator.

Kinematical constraints in a single VGTM module

Assume that the length of fixed rod on the side of VGTM module is $L_{d}$ . The maximum and minimum length of the actuator are $L_{max}$ and $L_{min}$ respectively. According to the parameters described in Figure 1, the six equality and six inequality constraints describing the VGTM module structure can be expressed as follows

\begin{array}{l} {(Q_{i - 1, 1} - Q_{i,1})}^{T} (Q_{i - 1, 1} - Q_{i,1}) = L_{d}^{2} \\ {(Q_{i - 1, 1} - Q_{i,2})}^{T} (Q_{i - 1, 1} - Q_{i,2}) = L_{d}^{2} \\ {(Q_{i - 1, 2} - Q_{i,2})}^{T} (Q_{i - 1, 2} - Q_{i,2}) = L_{d}^{2} \\ {(Q_{i - 1, 2} - Q_{i,3})}^{T} (Q_{i - 1, 2} - Q_{i,3}) = L_{d}^{2} \\ {(Q_{i - 1, 3} - Q_{i,1})}^{T} (Q_{i - 1, 3} - Q_{i,1}) = L_{d}^{2} \\ {(Q_{i - 1, 3} - Q_{i,3})}^{T} (Q_{i - 1, 3} - Q_{i,3}) = L_{d}^{2} \end{array}

\begin{array}{l} L_{min}^{2} \leq {(Q_{i,1} - Q_{i,2})}^{T} (Q_{i,1} - Q_{i,2}) \leq L_{max}^{2} \\ L_{min}^{2} \leq {(Q_{i,2} - Q_{i,3})}^{T} (Q_{i,2} - Q_{i,3}) \leq L_{max}^{2} \\ L_{min}^{2} \leq {(Q_{i,1} - Q_{i,3})}^{T} (Q_{i,1} - Q_{i,3}) \leq L_{max}^{2} \end{array}

After determining the location of each node in the VGTM module, the position of the top platform center of the i th VGTM module $Q_{i}$ can be described as

Q_{i} = (Q_{i,1} + Q_{i,2} + Q_{i,3}) / 3

When the nodes of all VGTM modules were identified, a configuration of VGTM is confirmed. Then the inverse kinematics of VGTM is to locate each node of the VGTM to make sure that the end of the VGTM reached the specified location.

Inverse kinematics

In this section, firstly, the inverse kinematics of VGTM is formulized as a minimization problem and simplified to a QCQP problem in general form. Then, according to the particularity of the VGTM structure, the Ipopt algorithm is adjusted to solve the optimization problem by making the most of the salient features of the optimized model.

Establishment of inverse kinematic optimization model

The nonsymmetrical octahedral VGTM studied in this article is shown in Figure 2. Suppose the node positions are represented by $Q_{i, j} (i = 1, 2, ..., n, j = 1, 2, 3)$ . The centre of the end of the manipulator is expressed as $Q_{n}$ . Then, the optimization model of inverse kinematics can be formulated as

The spatial distance between the desired position $Q_{t}$ and the terminal location $Q_{n}$ of VGTM must be minimized.

The structural constraints of VGTM have to be satisfied.

The VGTM will not collide with any of the obstacles. And assume that the sphere center and radius of the k th spherical obstacle are $O_{k}$ and $L_{k} (k = 1, 2, \dots, m)$ separately. In order to ensure the VGTM does not collide with any obstacle, it is assumed that each node of the manipulators is far enough away from the obstacle, and the distance between the obstacle and the midpoint of each fixed rod or actuator is also far enough (see Figure 3).

Figure 2.

Universal structure of the VGTM. VGTM: variable geometry truss manipulator.

Then the inverse kinematics problem can be described as the following nonconvex quadratic programming problem with quadratic constraints

\begin{array}{l} min_{Q_{i}} F (Q_{i}) = ∥ Q_{n} Q_{t} ∥^{2} \\ s.t. {(Q_{i - 1,1} - Q_{i,1})}^{T} (Q_{i - 1,1} - Q_{i,1}) = L_{d}^{2} \\ {(Q_{i - 1,1} - Q_{i,2})}^{T} (Q_{i - 1,1} - Q_{i,2}) = L_{d}^{2} \\ {(Q_{i - 1,2} - Q_{i,2})}^{T} (Q_{i - 1,2} - Q_{i,2}) = L_{d}^{2} \\ {(Q_{i - 1,2} - Q_{i,3})}^{T} (Q_{i - 1,2} - Q_{i,3}) = L_{d}^{2} \\ {(Q_{i - 1,3} - Q_{i,1})}^{T} (Q_{i - 1,3} - Q_{i,1}) = L_{d}^{2} \\ {(Q_{i - 1,3} - Q_{i,3})}^{T} (Q_{i - 1,3} - Q_{i,3}) = L_{d}^{2} \\ L_{min}^{2} \leq {(Q_{i,1} - Q_{i,2})}^{T} (Q_{i,1} - Q_{i,2}) \leq L_{max}^{2} \\ L_{min}^{2} \leq {(Q_{i,2} - Q_{i,3})}^{T} (Q_{i,2} - Q_{i,3}) \leq L_{max}^{2} \\ L_{min}^{2} \leq {(Q_{i,1} - Q_{i,3})}^{T} (Q_{i,1} - Q_{i,3}) \leq L_{max}^{2} \\ {(Q_{i, j} - O_{k})}^{T} (Q_{i, j} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i,1} + Q_{i,2}}{2} - O_{k})}^{T} (\frac{Q_{i,1} + Q_{i,2}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i,2} + Q_{i,3}}{2} - O_{k})}^{T} (\frac{Q_{i,2} + Q_{i,3}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i,1} + Q_{i,3}}{2} - O_{k})}^{T} (\frac{Q_{i,1} + Q_{i,3}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i - 1,1} + Q_{i,1}}{2} - O_{k})}^{T} (\frac{Q_{i - 1,1} + Q_{i,1}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i - 1,1} + Q_{i,2}}{2} - O_{k})}^{T} (\frac{Q_{i - 1,1} + Q_{i,2}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i - 1,2} + Q_{i,2}}{2} - O_{k})}^{T} (\frac{Q_{i - 1,2} + Q_{i,2}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i - 1,2} + Q_{i,3}}{2} - O_{k})}^{T} (\frac{Q_{i - 1,2} + Q_{i,3}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i - 1,3} + Q_{i,1}}{2} - O_{k})}^{T} (\frac{Q_{i - 1,3} + Q_{i,1}}{2} - O_{k}) \geq L_{k}^{2} \\ {(\frac{Q_{i - 1,3} + Q_{i,3}}{2} - O_{k})}^{T} (\frac{Q_{i - 1,3} + Q_{i,3}}{2} - O_{k}) \geq L_{k}^{2} \end{array}

where $Q_{i, j} (i = 1, 2, ..., n, j = 1, 2, 3)$ represent the node positions of the i th VGTM module, $L_{d}$ is the length of the fixed bar in the module, while $O_{k}$ and $L_{k} (k = 1, 2, \dots, m)$ are supposed as the center and radius of the k th sphere obstacle.

Figure 3.

Details for obstacle avoidance.

On the basis of the physical implication of every parameter in the above model, next, the simplification of the optimization model (4) will be detailed.

Suppose the optimized variable $Q^{T} = [Q_{1, 1}^{T}, Q_{1, 2}^{T}, Q_{1, 3}^{T}, \dots, Q_{n,1}^{T}, Q_{n,2}^{T}, Q_{n,3}^{T}]$ and define matrixes $A_{i, j} \in ℝ^{3 \times 9 n}$

\begin{array}{l} ​ \\ A_{i, j} = {\begin{cases} [O, O, \dots, I, \dots, O], i = 1, 2, \dots, n, j = 1, 2, 3 \\ 1, 2, \dots,3 (i - 1) + j, \dots,3 n \\ O^{3 \times 9 n}, i = 0 \end{cases} \end{array}

where O and I are zero matrix and unit matrix of 3 order, respectively. Then formulas $Q_{i, j} = A_{i, j} Q (i = 1, 2, \dots, n, j = 1, 2, 3)$ can be drawn and the square of the distance between $Q_{n}$ and $Q_{t}$ is $∥ Q_{n} Q_{t} ∥^{2} = \frac{1}{9} Q^{T} {(A_{n,1} + A_{n,2} + A_{n,3})}^{T} (A_{n,1} + A_{n,2} + A_{n,3}) Q - \frac{2}{3} Q_{d}^{T} (A_{n,1} + A_{n,2} + A_{n,3}) Q + Q_{d}^{T} Q_{d}$

Analogously, the six equality describing fixed rod length in (1), the six inequality constraints representing the actuator length of the VGTM structure in (2), and obstacle avoidance conditions in (4) can also be expressed as quadratic constraints on the optimization variables Q as follows

f_{i - 1, i}^{s t} = ∥ Q_{i - 1, s} Q_{i, t} ∥^{2} - L d^{2} = Q^{T} {(A_{i - 1, s} - A_{i, t})}^{T} (A_{i - 1, s} - A_{i, t}) Q - L d^{2} = 0

with $s = 1, t = 1, 2; s = 2, t = 2, 3; s = 3, t = 1, 3$ .

g_{i}^{s t} = ∥ Q_{i, s} Q_{i, t} ∥^{2} = Q^{T} {(A_{i, s} - A_{i, t})}^{T} (A_{i, s} - A_{i, t}) Q \in [L_{min}^{2}, L_{max}^{2}]

where $s = 1, t = 2; s = 2, t = 3; s = 1, t = 3$ .

\begin{matrix} h_{i, j} = L_{k}^{2} - ∥ Q_{i, j} O_{k} ∥^{2} \\ = L_{k}^{2} - Q^{T} A_{i, j}^{T} A_{i, j} Q + 2 O_{k}^{T} A_{i, j} Q - O_{k}^{T} O_{k} \leq 0 \end{matrix}

\begin{array}{l} h_{i}^{s t} = L_{k}^{2} - ∥ \frac{Q_{i, s} + Q_{i, t}}{2} O_{k} ∥^{2} \\ = L_{k}^{2} - \frac{1}{4} Q^{T} {(A_{i, s} + A_{i, t})}^{T} (A_{i, s} + A_{i, t}) Q \\ + O_{k}^{T} (A_{i, s} + A_{i, t}) Q - O_{k}^{T} O_{k} \leq 0 \end{array}

here, $s = 1, t = 2; s = 2, t = 3; s = 1, t = 3$ .

\begin{matrix} h_{i - 1, i}^{s t} = L_{k}^{2} - ∥ \frac{Q_{i - 1, s} + Q_{i, t}}{2} O_{k} ∥^{2} \\ = L_{k}^{2} - \frac{1}{4} Q^{T} {(A_{i - 1, s} + A_{i, t})}^{T} (A_{i - 1, s} + A_{i, t}) Q \\ + O_{k}^{T} (A_{i - 1, s} + A_{i, t}) Q - O_{k}^{T} O_{k} \leq 0 \end{matrix}

with $s = 1, t = 1, 2; s = 2, t = 2, 3; s = 3, t = 1, 3$ .

Consequently, the minimization problem (4) finally becomes

\begin{array}{l} min_{Q \in ℝ^{9 n}} F (Q) = \frac{1}{9} Q^{T} {(A_{n,1} + A_{n,2} + A_{n,3})}^{T} (A_{n,1} + A_{n,2} + A_{n,3}) Q \\ - \frac{2}{3} Q_{d}^{T} (A_{n,1} + A_{n,2} + A_{n,3}) Q + Q_{d}^{T} Q_{d} \\ s.t. f_{i - 1, i}^{s t} = 0, g_{i}^{s t} \in [L_{min}^{2}, L_{max}^{2}] \\ h_{i, j} \leq 0, h_{i}^{s t} \leq 0, h_{i - 1, i}^{s t} \leq 0, i = 1, 2, \dots, n, j = 1, 2, 3 \end{array}

Apparently, both the objective function and the constraints of the above minimization problem are quadratic, that is, it is a general nonconvex QCQP problem.

In order to solve the problem more conveniently in the later chapters, the above specific problem will be simplified to the following nonconvex QCQP optimization with more general form

\begin{array}{l} min_{Q \in ℝ^{9 n}} F (Q) \\ s.t. f_{i} (Q) = 0, i = 1, 2, \dots,6 n \\ g_{j} (Q) \geq 0, j = 1, 2, \dots,18 n \end{array}

where the objective function $F (Q)$ is defined the same as in the problem (5), while the sets of equality and inequality constraints are severally defined as below

\begin{array}{l} f (Q) = [f_{0, 1}^{11}, f_{0, 1}^{12}, f_{0, 1}^{22}, f_{0, 1}^{23}, f_{0, 1}^{31}, f_{0, 1}^{33}, \\ f_{1, 2}^{11}, f_{1, 2}^{12}, f_{1, 2}^{22}, f_{1, 2}^{23}, f_{1, 2}^{31}, f_{1, 2}^{33}, \\ \dots \dots, \\ f_{n - 1, n}^{11}, f_{n - 1, n}^{12}, f_{n - 1, n}^{22}, f_{n - 1, n}^{23}, f_{n - 1, n}^{31}, f_{n - 1, n}^{33}] \end{array}

\begin{array}{l} g (Q) = [g_{1}^{12} - L_{min}^{2}, g_{1}^{23} - L_{min}^{2}, g_{1}^{13} - L_{min}^{2}, L_{max}^{2} - g_{1}^{12}, L_{max}^{2} - g_{1}^{23}, L_{max}^{2} - g_{1}^{13}, \\ h_{1,1}, h_{1,2}, h_{1,3}, h_{1}^{12}, h_{1}^{23}, h_{1}^{13}, h_{0,1}^{11}, h_{0,1}^{12}, h_{0,1}^{22}, h_{0,1}^{23}, h_{0,1}^{31}, h_{0,1}^{33}, \\ g_{2}^{12} - L_{min}^{2}, g_{2}^{23} - L_{min}^{2}, g_{2}^{13} - L_{min}^{2}, L_{max}^{2} - g_{2}^{12}, L_{max}^{2} - g_{2}^{23}, L_{max}^{2} - g_{2}^{13}, \\ h_{2,1}, h_{2,2}, h_{2,3}, h_{2}^{12}, h_{2}^{23}, h_{2}^{13}, h_{1,2}^{11}, h_{1,2}^{12}, h_{1,2}^{22}, h_{1,2}^{23}, h_{1,2}^{31}, h_{1,2}^{33}, \\ \dots \dots \dots \dots, \\ g_{n}^{12} - L_{min}^{2}, g_{n}^{23} - L_{min}^{2}, g_{n}^{13} - L_{min}^{2}, L_{max}^{2} - g_{n}^{12}, L_{max}^{2} - g_{n}^{23}, L_{max}^{2} - g_{n}^{13}, \\ h_{n,1}, h_{n,2}, h_{n,3}, h_{n}^{12}, h_{n}^{23}, h_{n}^{13}, h_{n - 1, n}^{11}, h_{n - 1, n}^{12}, h_{n - 1, n}^{22}, h_{n - 1, n}^{23}, h_{n - 1, n}^{31}, h_{n - 1, n}^{33}] \end{array}

Solving inverse kinematics via interior point algorithm with inexact step computation

Typically, a generalized nonconvex QCQP problem described as in (6) is NP-hard in general. The existing methods solve nonconvex problems by linearizing their nonconvex parts or relaxing the nonconvexity via semi-definite relaxation. Nevertheless, when the quadratic constraint matrices are indefinite, these approaches are rarely successful in even just acquiring a feasible solution. Nearly 30 years since the interior point algorithm has been invented, it has become a core part of many existing nonlinear optimization algorithms. If its superlinear convergence is combined with powerful linear solver and inaccurate Newton method, this algorithm may even quickly solve a generalized nonconvex nonlinear optimization problem. Next, we will first introduce the uniqueness of the actual problem established in this article. And then, we will describe how to apply the inexact interior point method to solve this problem in detail.

Uniqueness of the optimization problem

The inverse kinematics optimization problem of the VGTM established in this article is an interesting topic and has the following three distinctive features.

The first is sparsity and binarity. As stated previously, as the number of VGTM modules increases, the size of the matrix Q in the target function and the number of constraints in the optimization problem can become extremely large. Despite these matrices appear to be huge, due to sparsity, the valid values stored in the calculation process are actually quite limited. Vigerske Stefan and Andreas Wächter designed a solver Ipopt to solve large-scale sparse problems.³⁰ As can be seen from the modeling process of the problem, there are many large-scale sparse matrices in the optimization model. While the solver Ipopt needs to be provided with sparsity to speed up the calculation before optimization begins, which is quite suitable for solving our problem. However, another important uniqueness of this practical problem has also drawn our interest in that all matrices Q are not only highly sparse but also contain only two possible valid values of one and minus one besides zero. While the solver Ipopt stores each entry in the matrix as a floating point form, which greatly increases the cost of space even if it is provided the same matrix, because it does not know that all the values to be stored are integers, not to mention binary. Consequently, these sparse matrices are also closely related to binarity, and this property can lead to improvements of the algorithm.

The second is its special way of convergence. In contrast to a pure theoretical problem, the inverse kinematics of the VGTM follows some observable rules, which makes the direction of the iteration clearer. And what is also unique is the way in which the optimization problem converges to the optimal solution, so the following issue is presented. Since the VGTM is a robotic arm and needs to determine the status (such as the length of each actuator) before it performs any operational tasks, the amendment algorithm will not consider the generation of initial points. In addition, the constraint conditions will invariably be within the feasible region that satisfies the initial conditions until it achieves the optimal or converges to the stationary point nearest the optimal one. The desired ideal points are the ones that closest to the local minimum point or the boundary points nearest the initial status. This is because the deformation of the VGTM is limited by its fixed mechanical limitations, and in any case there will never be any actuator beyond its physical limits without destroying the manipulator. Consequently, the feasibility recovery part embedded within the solver Ipopt is needless to our specific problem. And the conditions for implementation can also be adjusted appropriately according to this particular problem. It will also be discussed in the following algorithm section.

Another important uniqueness is that, for VGTM or other general asymmetric VGT types of robots, once its structure is established, the system can be described using an adjacency matrix. And the transformation from the mapping to Hessian matrix and Jacobian matrix is directly available without any manual computations to obtain these important information for the solver. That is to say, even though the system may have variety and asymmetry, as long as its structure is determined, for all VGT-type robots, this procedure can be accomplished automatically. The principle is not hard to understand because it is build on the mechanical coupling between the joints of the VGT structure.

Solving QCQP problem via inexact interior point optimization

In this section, the inaccurate interior point algorithm for the inverse kinematics of the VGTM will be discussed in detail. This algorithm is modified specifically on the basis of the inexact algorithm Ipopt in Curtis et al.,^31,32 in particular, the storage cost and CPU or GPU time in each iteration are minimized.

In order to solve the original optimization problem (6), a barrier parameter μ is first introduced, and then the original problem is transformed into the following sub-problem with the barrier parameter decreasing

\begin{array}{l} ​ & min_{Q, r} F (Q) - μ \sum_{j = 1}^{18 n} ln (r_{j}) \\ s.t. & f_{i} (Q) = 0, i = 1,2, \dots,6 n \\ ​ & g_{j} (Q) = r_{j}, r_{j} > 0, j = 1,2, \dots, 18 n \end{array}

Here the barrier parameter μ is usually expected to be iteratively reduced. On the basis of the further amendment of the damping process of μ in Gould et al.,³³ the update process can require fewer sequences of μ and achieve superlinear convergence speed. In addition, by introducing relaxation variables in (7), inequality constraints can be transformed into equality ones.

Suppose that the Lagrangian multipliers of the equality and inequality constraints are $λ_{E} \in ℝ^{6 n}$ and $λ_{I} \in ℝ^{18 n}$ respectively, then the Lagrangian function of the optimization problem (7) can be expressed as

L (Q, r, λ_{E}, λ_{I}; μ) = F (Q) - μ \sum_{j = 1}^{18 n} ln (r_{j}) + λ_{E}^{T} f (Q) + λ_{I}^{T} (g (Q) - r)

where the objective function $F (Q)$ , the constraint functions $f (Q)$ and $g (Q)$ are all second-order differentiable. Therefore, the first-order optimality conditions (Karush–Kuhn–Tucker, KKT conditions) for the nonlinear optimization problem (7) are

\begin{matrix} \nabla F (Q) + \nabla f (Q) λ_{E} + \nabla g (Q) λ_{I} = 0 \\ - μ R^{- 1} e - λ_{I} = 0 \\ f (Q) = 0 \\ g (Q) - r = 0 \end{matrix}

where $r > 0$ and $R = diag (r)$ . While e is defined as an 18n-dimensional vector consisting of element 1.

When the first order optimality conditions in formula (9) are not satisfied, that is, the original optimization problem is infeasible, we assume that the optimal solution of the algorithm converges toward the stationary point of the feasible problem as follows

min_{Q} \frac{1}{2} ∥ f (Q) ∥^{2} + \frac{1}{2} ∥ g {(Q)}^{-} ∥^{2}

with $g {(Q)}^{-} = max {0, - g (Q)}$ and the piecewise function max is assumed to be element-wise. In addition, since $∥ g {(Q)}^{-} ∥^{2}$ is differentiable and with its derivative $\nabla (∥ g {(Q)}^{-} ∥^{2}) = - 2 \nabla g (Q) g {(Q)}^{-}$ , another expression for $g {(Q)}^{-}$ can be obtained as follows

\begin{array}{l} g {(Q)}^{-} = r - g (Q), - R (g (Q) - r) = 0 \end{array}

The latter equation in formula (11) guarantees that $g_{j} (Q)$ is not restricted by this equation when $r_{j} = 0$ , but only when r_j is not equal to zero can $g_{j} (Q) = r_{j}$ be obtained.

By differentiating the formula (10), it can be concluded that the point that converges at its first-order optimality can also be characterized as the solution of the following nonlinear equations

\begin{array}{l} \nabla f (Q) f (Q) - \nabla g (Q) g {(Q)}^{-} = 0 \end{array}

When $r \geq 0$ and $g (Q) - r \leq 0$ , this equation is equivalent to the following one

\begin{matrix} \nabla f (Q) f (Q) + \nabla g (Q) (g (Q) - r) = 0 \\ - R (g (Q) - r) = 0 \end{matrix}

In order to extend the original and dual iterates Q and r into a single vector, define the following augmented iterates as

q : = [\begin{matrix} Q \\ r \end{matrix}] and λ : = [\begin{matrix} λ_{E} \\ λ_{I} \end{matrix}]

respectively, with $λ > 0$ . While the barrier objective function and constraints are respectively defined as

φ_{μ} (q) = F (Q) - μ \sum_{j = 1}^{18 n} ln (r_{j}) and c (q) = [\begin{matrix} f (Q) \\ g (Q) - r \end{matrix}]

So the Lagrange function in formula (8) becomes

L (q, λ; μ) = φ_{μ} (q) + λ^{T} c (q)

Notice that, with a penalty parameter $η > 0$ , by replacing $c (q)$ with its norm, the new penalty function below will have the same convergent process as the Lagrange equation in (14)

ϕ (q; μ, η) = φ_{μ} (q) + η ∥ c (q) ∥

The corresponding first derivatives of the objective function and constraint are

\nabla φ_{μ} (q) = [\begin{matrix} \nabla F (Q) \\ - μ e \end{matrix}] and \nabla c (q) = [\begin{matrix} \nabla f {(Q)}^{T} & 0 \\ \nabla g {(Q)}^{T} & - R \end{matrix}]

And the Hessian matrix of the Lagrangian at $(q, λ)$ is

\nabla_{Q Q}^{2} L (q, λ; μ) = [\begin{matrix} \nabla_{Q Q}^{2} F & 0 \\ 0 & μ e^{T} R^{- 2} e \end{matrix}] + \sum_{i = 1}^{6 n} λ_{E}^{(i)} [\begin{matrix} \nabla_{Q Q}^{2} f_{i} & 0 \\ 0 & 0 \end{matrix}] + \sum_{j = 1}^{18 n} λ_{I}^{(j)} [\begin{matrix} \nabla_{Q Q}^{2} g_{j} & 0 \\ 0 & 0 \end{matrix}]

where $f_{i}, λ_{E}^{(i)}$ and $g_{j}, λ_{I}^{(j)}$ denote severally the ith and jth constraint conditions and their corresponding dual variables. In the meantime, after the primal and dual iterates are merged, the KKT conditions (9) turn into the following form

\begin{matrix} \nabla φ_{μ} (q) + \nabla c {(q)}^{T} λ = 0 \\ c (q) = 0 \end{matrix}

According to the definition of $c (q)$ in the formula (13), the system of nonlinear equations (12) can be rewritten as the following matrix form

\nabla c {(q)}^{T} c (q) = 0

In the next step, in order to derive the Newton iteration d for the sequence of barrier subproblems (7) from an iterate $(q_{k}, λ_{k})$ , it is necessary to start with solving the following system of linear equations

\nabla G_{μ} (q_{k}, λ_{k}) d = - G_{μ} (q_{k}, λ_{k})

with $G_{μ} (q_{k}, λ_{k}) = [\begin{matrix} \nabla φ_{μ} (q_{k}) + \nabla c {(q_{k})}^{T} λ_{k} \\ c (q_{k}) \end{matrix}]$ .

As described in Akrotirianakis and Rustem,³⁴ the equation (17) can be converted to

[\begin{matrix} \nabla_{Q Q}^{2} L (q_{k}, λ_{k}) & \nabla c {(q_{k})}^{T} \\ \nabla c (q_{k}) & 0 \end{matrix}] [\begin{matrix} d_{k} \\ δ_{k} \end{matrix}] = - [\begin{matrix} \nabla ϕ_{μ} (q_{k}) + \nabla c {(q_{k})}^{T} λ_{k} \\ c (q_{k}) \end{matrix}]

Note: Obtaining the exact solution of equation (18) is computationally costly. Consequently, our objective is to solve the equation (18) by using an iterative linear system solver, under the condition that the imprecise solution is permitted and the global convergence of the algorithm is ensured. For this, we define the dual and constrained residual vector corresponding to the first and second equations in equation (18), respectively, as

\begin{matrix} ρ_{k} (d_{k}, δ_{k}) : = \nabla φ_{μ} (q_{k}) + \nabla_{Q Q}^{2} L (q_{k}, λ_{k}) d_{k} + \nabla c {(q_{k})}^{T} (λ_{k} + δ_{k}) \\ ζ_{k} (d_{k}) : = c (q_{k}) + \nabla c (q_{k}) d_{k} \end{matrix}

If (18) is accurately solved, then $ρ_{k} (d_{k}, δ_{k})$ and $ζ_{k} (d_{k})$ are zero. But in our algorithm, we only require them to be relatively small in norm. That is, suppose that the overall primal-dual relative residual is

Ψ_{k} (d_{k}, δ_{k}) : = ‖ [\begin{matrix} ρ_{k} (d_{k}, δ_{k}) \\ ζ_{k} (d_{k}) \end{matrix}] ‖ / ‖ [\begin{matrix} \nabla φ_{μ} (q_{k}) + \nabla c {(q_{k})}^{T} λ_{k} \\ c (q_{k}) \end{matrix}] ‖

then in order to promote rapid convergence, this residual must be small. Hence, our algorithm aims at calculating the step that the relative residual is lower than the expected threshold. However, we will remove this constraint and select a potentially less precise solution that satisfies our other terminating conditions, on condition that the iterated linear system solver cannot reach this precision after a specified number of iterations. Furthermore, only when $Ψ_{k} (d_{k}, δ_{k})$ is very small, we consider triggering a test to modify the Hessian matrix (to insure a descent), because otherwise the Hessian may be modified unnecessarily according to an imprecise intermediate solution in the iterative solver.

In order to handle the probable rank deficiency appropriately, the linear equation is replaced by an optimization problem. This not only avoids the factorization process but also makes it possible to calculate the inexact step, which still applies to nonconvex programming. Therefore, replacing d_k with the sum of a normal step v_k and a tangential step u_k gives the following formula

d_{k} = v_{k} + u_{k}

Here v_k is a move toward satisfying the linear model of constraints that within a trust region and is determined by solving the following optimization problem

\begin{array}{l} min_{v} \frac{1}{2} ∥ c (q_{k}) + \nabla c (q_{k}) v ∥^{2} \\ s.t. ∥ v ∥ \leq ω ∥ \nabla c {(q_{k})}^{T} c (q_{k}) ∥ \end{array}

Here $ω > 0$ is chosen as a constant. From formula (16), it can be clearly seen that the displacement of a stationary point in v_k is zero. It should be noted that this equation takes the same form as a trust domain problem, and so it can be handled in the same way. Many of the existing algorithms are still applicable, and inexact dogleg method is proved superior to other candidates.³⁵ After v_k has been obtained, substituting its decreasing direction(namely $- \nabla c (q_{k}) v_{k}$ ) for $c (q_{k})$ , then d_k can be calculated by solving the disturbed Newton system as follows

[\begin{matrix} \nabla_{Q Q}^{2} L (q_{k}, λ_{k}) & \nabla c {(q_{k})}^{T} \\ \nabla c (q_{k}) & 0 \end{matrix}] [\begin{matrix} d_{k} \\ δ_{k} \end{matrix}] = - [\begin{matrix} \nabla φ_{μ} (q_{k}) + \nabla c {(q_{k})}^{T} λ_{k} \\ - \nabla c (q_{k}) v_{k} \end{matrix}]

Solving the linear equation above takes up most of the memory costs and CPU time. And the performance of optimization algorithm depends, to a great extent, on the selected linear solver. In general, there are two methods for solving linear equations, one is a direct approach based on Cholesky factorization, and the other is an indirect and iterative approach based on the conjugate gradient (CG) method. Because the cost of memory is a major problem encountered in implementation, we are particularly interested in implementing a most advanced iteration solver, which aims to make the total memory cost minimized while maintaining the scalability of parallelization. Compared to the original technique presented in Curtis et al.,³² such an implementation is expected to reduce memory usage significantly. Exploring begins by first converting it to a QP issue through the extension and insertion of (19).

\nabla_{Q Q}^{2} L (q_{k}, λ_{k}) (v_{k} + u_{k}) + \nabla c {(q_{k})}^{T} δ_{k} = - (\nabla ϕ_{μ} (q_{k}) + \nabla c {(q_{k})}^{T} λ_{k})

\nabla c (q_{k}) u_{k} = 0

Integrating formula (21) with respect to variable u_k, according to equation (22), then the primal linear system of equations can be converted into a QP problem about u as follows

\begin{array}{l} min_{u} {(\nabla φ_{μ} (q_{k}) + \nabla_{Q Q}^{2} L (q_{k}, λ_{k}) v_{k})}^{T} u + \frac{1}{2} u^{T} \nabla_{Q Q}^{2} L (q_{k}, λ_{k}) u \\ s.t. \nabla c (q_{k}) u = 0 \end{array}

In general, on account that any inertia information of the primal-dual matrix cannot be obtained from the iterative linear system solver used in the disturbed Newton system (20), we may not assure that $\nabla_{Q Q}^{2} L (q_{k}, λ_{k})$ must be positive definite in the null space of matrix $\nabla c (q_{k})$ . When the matrix $\nabla_{Q Q}^{2} L (q_{k}, λ_{k})$ is positive definite in the null space of $\nabla c (q_{k})$ , the solution of the perturbed Newton system (20) corresponds to the solution of the QP problem (23). If matrix $\nabla_{Q Q}^{2} L (q_{k}, λ_{k})$ is not adequately positive definite, then it needs to be modified or replaced by a fully positive definite approximate matrix to guarantee that the consequent solution of system (20) is a suitable search direction. In our algorithm implementation, we adopt a frequently used technique to increase all the eigenvalues of the Hessian matrix, that is, adding a multiple of a positive definite diagonal matrix to $\nabla_{Q Q}^{2} L (q_{k}, λ_{k})$ , and modify the Hessian matrix according to this method to make it to be fully positive definite and uniformly bounded after a limited number of amendments. Obviously, the constraint conditions ensure that the positive definiteness (PD) needs to be satisfied only in the null space of $\nabla c (q_{k})$ , which is more relaxed than simply just requiring PD conditions. In addition, the KKT condition for the QP problem (23) above can be stated as

[\begin{matrix} \nabla_{Q Q}^{2} L (q_{k}) & \nabla c {(q_{k})}^{T} \\ \nabla c (q_{k}) & 0 \end{matrix}] [\begin{matrix} u \\ δ \end{matrix}] = - [\begin{matrix} \nabla φ_{μ} (q_{k}) + \nabla c {(q_{k})}^{T} λ_{k} \\ 0 \end{matrix}]

Next, the preconditioned CG approach that suitable for sparse and large-scale problems is adopted to deal with the above problem (24), which is mainly because CG is proved to be the optimal choice when the variable size exceeds $10^{9}$ , and it is also an iterative approach that maintains a high efficiency on large-scale problems but with fairly low storage requirements. In addition, it can also well allocate computational efficiency to matrix-vector multiplication, which makes it possible to utilize a state-of-art sparse matrix vector multiplication approach. And even from the problem description stage of the QCQP problem, it makes an overall memory reduction possible.

After getting the decent search direction, the maximum step length $α_{k}^{max} \in (0, 1]$ will now be determined based on the fraction-to-the-boundary rule, rather than simply using 1 as a step length. That is, $α_{k}^{max}$ is determined according to the following inequality

r_{k} + α_{k}^{max} R_{k} d_{k}^{r} \leq (1 - τ_{1}) r_{k}

Where $τ_{1}$ is a constant between 0 and 1. Then, in the case that the penalty function in formula (15) satisfies the Armijo condition below, the actual step length $α_{k} \in (0, α_{k}^{max})$ can be determined according to the following inequality which gives a fundamental sufficient decline requirement for $α_{k}$

ϕ (q_{k} + α_{k} {\tilde{d}}_{k}; μ_{j}, η_{k}) \leq ϕ (q_{k}; μ_{j}, η_{k}) + τ_{2} α_{k} \nabla ϕ_{k} ({\tilde{d}}_{k}; μ_{j}, η_{k})

Here $τ_{2} \in (0, 1)$ , while $\nabla ϕ_{k} ({\tilde{d}}_{k}; μ_{j}, η_{k})$ is the directional derivative of the penalty function ϕ at variable q_k along ${\tilde{d}}_{k}$ and $η_{k}$ can guarantee $\nabla ϕ_{k} ({\tilde{d}}_{k}; μ_{j}, η_{k}) < 0$ . In addition, since the original search direction d_k defined in formula (18) is not standard, the original iteration will be updated by using the line search coefficient $α_{k}$ as follows

q_{k + 1} \leftarrow q_{k} + α_{k} {\tilde{d}}_{k} with {\tilde{d}}_{k} \leftarrow [\begin{matrix} d_{k}^{Q} \\ R_{k} d_{k}^{r} \end{matrix}]

where $d_{k}^{Q}$ and $d_{k}^{r}$ denote the components of the search direction d_k corresponding to the optimization variables Q and r separately. Here, of particular note is the termination criterion for the algorithm is beyond our consideration and some recent research results in this respect can refer to the study by Curtis et al.³⁶ and Grote et al.²⁸

In conclusion, a flowchart of the proposed inverse kinematics optimization algorithm for VGTM can be described as:

Numerical simulations for inverse kinematics

In this part, the results of proposed inverse kinematics optimization for VGTM composed of three modules are presented. We simulate a VGTM with three modules using a personal computer that equipped with Windows 10, a CPU of i7-4790k and a memory of 16 GB DDR3. The parameters of the VGTM are given. For a VGTM with 3 modules, it has 18 lateral passive rods and 9 actuators. The length of each passive member is $L_{d}$ = 475 mm and the motion range of each actuator length is 500 mm = $L_{min} < L < L_{max}$ = 700 mm, as stated in Table 1. To illustrate the validity of our presented approach, T different terminal positions of VGTM are chosen.

Table 1.

Preset parameters of VGTM for simulation.

Parameters	Values (mm)
S	82
$L_{d}$	475
$L_{min}$	500
$L_{max}$	700

VGTM: variable geometry truss manipulator.

Assuming that T sampling points and n VGTM modules are taken for simulation, then the absolute error can be worked out based on the formula below

\begin{array}{l} E_{a b s} = \frac{1}{T} \sum_{j = 1}^{T} \sqrt{{(P_{n x}^{j} - P_{t x}^{j})}^{2} + {(P_{n y}^{j} - P_{t y}^{j})}^{2} + {(P_{n z}^{j} - P_{t z}^{j})}^{2}} \end{array}

In the light of the conclusion obtained in “Solving inverse kinematics via interior point algorithm with inexact step computation” section, it is assumed that the centre and radius of the spherical obstacle are $O_{m} {= (150, 260, 260)}^{T}$ and $L_{m} = 150$ , respectively. Zhao et al.³⁷ and Yang et al.³⁸ have studied manipulators consisting of double octahedral modules. The limiting position of the end of such a manipulator that composed of three modules is equivalent to that composed of six modules in this article, that is, the expected terminal positions used in Zhao et al.³⁷ and Yang et al.³⁸ are beyond the reachable range of the manipulator consisting of three modules described in this article. Therefore, we did not use benchmark trajectories to test the algorithm. Instead, we randomly selected five terminal positions in the reachable workspace of the three-module manipulator to test the algorithm as displayed in Figures 4 to 8. The data about corresponding actuator length for each VGTM module are listed in Table 2. In addition, the comparison results on positioning errors with the methods used in literatures Zhao et al.³⁷ and Yang et al.³⁸ are listed in Table 3. All these data demonstrate that the optimized terminal locations converge to the expected reference positions.

Figure 4.

Configuration of asymmetrical octahedral VGTM for obstacle avoidance with Qd=[765;350;290]. VGTM: variable geometry truss manipulator.

Figure 5.

Configuration of asymmetrical octahedral VGTM for obstacle avoidance with Qd=[785;350;350]. VGTM: variable geometry truss manipulator.

Figure 6.

Configuration of asymmetrical octahedral VGTM for obstacle avoidance with Qd=[920;350;350]. VGTM: variable geometry truss manipulator.

Figure 7.

Configuration of asymmetrical octahedral VGTM for obstacle avoidance with Qd=[1025;0;0]. VGTM: variable geometry truss manipulator.

Figure 8.

Configuration of asymmetrical octahedral VGTM for obstacle avoidance with Qd=[1100;150;150]. VGTM: variable geometry truss manipulator.

Table 2.

Actuator length of each VGTM module-three modules.

Sample points	Desired position (X, Y, Z/mm)	Actuator length/mm			Actual terminal position (X, Y, Z/mm)	Error (mm)
Sample points	Desired position (X, Y, Z/mm)	Module 1	Module 2	Module 3	Actual terminal position (X, Y, Z/mm)	Error (mm)
1	765	500.0032	500.0000	503.8201	764.5881	0.6236
	350	592.5435	700.0000	644.2914	350.4198
	290	700.0000	667.1490	685.3796	290.2072
2	785	500.0000	500.0000	500.0000	784.7064	0.5484
	350	540.8883	700.0000	604.5915	349.5387
	350	700.0000	628.0208	645.9678	349.9585
3	920	500.0000	500.0000	500.0000	919.2154	1.1441
	350	525.0215	633.4575	500.0000	349.8959
	350	700.0000	500.0000	500.0000	349.1739
4	1025	618.5270	515.9884	636.4290	1025.5819	0.6441
	0	609.8607	500.0000	685.0384	−0.1245
	0	612.6441	500.0197	691.7717	−0.2464
5	1100	500.0000	500.0000	500.0000	1100.1558	0.9981
	150	523.6896	521.2740	519.8418	150.9405
	150	604.7051	500.0000	514.5251	150.2956

Table 3.

Comparison results of positioning errors with methods in Zhao et al.³⁷ and Yang et al.³⁸

Sample points	Desired position	Our method	Method in Zhao et al.³⁷	Method in Yang et al.³⁸	Our errors	Errors in Zhao et al.³⁷	Errors in Yang et al.³⁸
1	765	764.5881	762.0208	757.9678	0.6236	3.1201	7.2455
	350	350.4198	349.1678	348.5159
	290	290.2072	289.5915	289.0816
2	785	784.7064	781.9884	778.2587	0.5484	3.1322	6.8572
	350	349.5387	349.4575	349.1591
	350	349.9585	349.3318	349.0678
3	920	919.2154	916.0197	912.4293	1.1441	4.0550	7.6734
	350	349.8959	349.5467	349.1384
	350	349.1739	349.3716	349.0924
4	1025	1025.5819	1021.9488	1017.4829	0.6441	3.1928	7.6039
	0	−0.1245	−0.5621	−0.6854
	0	−0.2464	−0.7538	−0.9177
5	1100	1100.1558	1096.8418	1092.1525	0.9981	3.1757	7.9231
	150	150.9405	149.7240	149.2740
	150	150.2956	149.8146	149.1846

Based on the data obtained from Tables 2 and 3, the average absolute positioning error of the actual terminal location calculated from the physical platform of the VGTM is 0.7917 mm, while the errors in Zhao et al.³⁷ and Yang et al.³⁸ are respectively 3.3352 mm and 7.4606 mm. Furthermore, as can be seen from Table 4, the average running time of the algorithm proposed in this article is 0.1689 s, while those in Zhao et al.³⁷ and Yang et al.³⁸ are 0.265 s and 0.233 s, respectively. That is, the VGTM can reach a given position with great accuracy in a relatively fast time.

Table 4.

Comparison results of computation time with methods in Zhao et al.³⁷ and Yang et al.³⁸

Sample points	Desired position (X, Y, Z/mm)	$t_{CPU}$ /s
Sample points	Desired position (X, Y, Z/mm)	Our Method	Methods in Zhao et al.³⁷	Methods in Yang et al.³⁸
1	${[765, 350, 290]}^{T}$	0.1382	0.3112	0.2418
2	${[785, 350, 350]}^{T}$	0.1935	0.2457	0.2612
3	${[920, 350, 350]}^{T}$	0.1795	0.2396	0.2176
4	${[1025, 0, 0]}^{T}$	0.1646	0.2784	0.2154
5	${[1100, 150, 150]}^{T}$	0.1687	0.2501	0.2290

Conclusion

The inverse kinematics with obstacles avoidance for an asymmetric octahedral VGTM has been accomplished in our article. Considering the structural constraints inherent in VGTM and obstacle information, a new inverse kinematics optimization model was proposed. The inverse kinematics is cast as a generalized nonconvex QCQP problem. By well utilizing the distinctive characteristics of the actual optimization problem established, the inexact interior point optimization has been adjusted correspondingly, which greatly saved the computation time of the algorithm. Simulation results of obstacle avoidance inverse kinematics for a nonsymmetrical octahedral three-module VGTM validated our conclusion.

Footnotes

Acknowledgements

The authors will be greatly indebted to the editors and each anonymous reviewer for their constructive suggestions that helped to improve this article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The article reported here have the joint support of both the National Natural Science Foundation of China (grant 61175028, 61374161, 61673262) and the Key Project of Natural Science Foundation of Shanghai (16JC1401100).

ORCID iD

Yanchun Zhao

References

Miura

Furuya

Suzuki

. Variable geometry truss and its application to deployable truss and space crane arm. Acta Astronaut 1985; 12: 599–607.

Miura

Furuya

. Adaptive structure concept for future space applications. AIAA J 1988; 26: 995–1002.

Agirrebeitia

Avilés

De Bustos

. A method for the study of position in highly redundant multibody systems in environments with obstacles. IEEE Trans Robot Autom 2002; 18(2): 257–262.

Agirrebeitia

Avilés

Bustos

IFD

. Inverse position problem in highly redundant multibody systems in environments with obstacles. Mech Mach Theory 2003; 38(11): 1215–1235.

Lee

Sanderson

. Dynamic analysis and distributed control of the tetrobot modular reconfigurable robotic system. Auton Robots 2001; 10: 67–82.

Hanahara

Tada

. Variable geometry truss with SMA wire actuators (Basic consideration on kinematical and mechanical characteristics). In: 19th international conference on adaptive structures and technologies, Ascona, Switzerland, 6–9 October 2008, pp. 96–100. Zurich: Swiss Federal Laboratories for Materials Science and Technology (Empa).

Rost

Uhlemann

Modler

. On the joint design and hydraulic actuation of octahedron VGT robot manipulators. In: IEEE/ASME international conference on advanced intelligent mechatronics, Wollongong, Australia, 9–12 July 2011, pp. 92–97. Piscataway, NJ: IEEE/ASME.

Badescu

Mavroidis

. Workspace analysis of reconfigurable robotic arms using parallel platforms as modules. In: Proceedings of ASME international mechanical engineering congress exposition, Washington, DC, 16–21 November 2003, pp. 717–726. American Society of Mechanical Engineers (ASME).

Badescu

Mavroidis

. New performance indices and workspace analysis of reconfigurable hyper-redundant robotic arms. Int J Robot Res 2004; 23(6): 643–659.

10.

Yahya

Moghavvemi

Mohamed

HAF

. Geometrical approach of planar hyper-redundant manipulators: inverse kinematics, path planning and workspace. Simul Model Pract Theory 2011; 19(1): 406–422.

11.

Rocco

Eklund

Sommese

. Algebraic C*-actions and the inverse kinematics of a general 6R manipulator. Appl Math Comput 2010; 216(9): 2512–2524.

12.

Chiddarwar

Babu

. Comparison of RBF and MLP neural networks to solve inverse kinematic problem for 6R serial robot by a fusion approach. Eng Appl Artif Intell 2010; 23(7): 1083–1092.

13.

Tutsoy

Barkana

Colak

. Learning to balance an NAO robot using reinforcement learning with symbolic inverse kinematic. Transact Institute Measure Control 2016; 39(11): 1735–1748.

14.

Önder

. Cpg based RL algorithm learns to control of a humanoid robot leg. Acta Press 2015; 30(2): 178–183.

15.

Zanganeh

Angeles

. The inverse kinematics of hyper-redundant manipulators using splines. In: IEEE international conference on robotics and automation, vol. 3, Nagoya, Japan, 21–27 May 1995, pp.2797–2802. New York: IEEE.

16.

Chirikjian

Burdick

. A hyper-redundant manipulator. IEEE Robot Autom Mag 1994; 1: 22–29.

17.

Chirikjian

Burdick

. Kinematically optimal hyper-redundant manipulator configurations. In: IEEE international conference on robotics and automation, Nice, France, 12–14 May 1992, pp.415–420. New York: IEEE.

18.

Fahimi

Ashrafioun

Nataraj

. An improved inverse kinematic and velocity solution for spatial hyper-redundant robots. IEEE Trans Robot Autom 2002; 18(1): 103–107.

19.

Bayram

Özgören

Bayram

. The conceptual design of a spatial binary hyper-redundant manipulator and its forward kinematics. Proc Instit Mech Eng C J Mech Eng Sci 2012; 226(1): 217–227.

20.

Bayram

Özgören

Bayram

. The position control of a spatial binary hyper-redundant manipulator through its inverse kinematics. Proc Institution Mech Eng Part C: J Mech Eng Sci 2013; 227(2): 359–372.

21.

Iwatsuki

Nishizaka

Morikawa

. Motion control of hyper-redundant manipulator built by serially connecting many parallel mechanism units with a few DOF. Int J Autom Tech 2010; 4(4): 364–371.

22.

Mintenbeck

Estanna

. Design modelling and control of a hyper-redundant 3-RPS parallel mechanism. In: IEEE international conference on robotics and biomimetics, Tianjin, China, 14–18 December 2010, pp. 591–596. New York: IEEE.

23.

Chibani

Mahfoudi

Chettibi

. Generating optimal reference kinematic configurations for hyper-redundant parallel robots. Proc Instit Mech Eng C J Mech Eng Sci 2015; 229(9): 867–882.

24.

Bilbao

Avilés

Aguirrebeitia

. Eigensensitivity-based optimal damper location in variable geometry trusses. AIAA J 2009; 47(3): 576–591.

25.

De Zárate

Aguirrebeitia

Avilés

. A method for the inverse dynamics analysis of variable geometry trusses. In: International conference of numerical analysis and applied mathematics, Rhodes, Greece, 19–25 September 2010, vol. 1281 (1), pp. 1275–1280. American Institute of Physics, publishing.

26.

Rost

Uhlemann

Modler

. On the joint design and hydraulic actuation of octahedron VGT robot manipulators. In: IEEE/ASME international conference on advanced intelligent mechatronics, Budapest, Hungary, 3–7 July 2011, pp. 92–97. IEEE Xplore.

27.

Gondzio

. Interior point methods 25 years later. Eur J Oper Res 2012; 218(3): 587–601.

28.

Grote

Huber

Kourounis

. Inexact interior-point method for PDE-constrained nonlinear optimization. SIAM J Sci Comput 2014; 36(3): A1251–A1276.

29.

Kang

Chiang

Laird

. Nonlinear programming strategies on high-performance computers. In: IEEE conference on decision and control (CDC), Osaka, Japan, 15–18 December 2015, pp. 4612–4620. IEEE.

30.

Stefan

Wächter

Introduction to IPOPT: A tutorial for downloading, installing, and using IPOPT, 2422nd ed. 2013.

31.

Curtis

Nocedal

Wächter

. A matrix-free algorithm for equality constrained optimization problems with rank-defcient Jacobians. SIAM J Optimiz 2009; 20(3): 1224–1249.

32.

Curtis

Schenk

Wächter

. An interior-point algorithm for large-scale nonlinear optimization with inexact step computations. SIAM J Sci Comput 2010; 32(6): 3447–3475.

33.

Gould

NIM

Orban

Sartenaer

. Superlinear convergence of primal-dual interior point algorithms for nonlinear programming. SIAM J Optimiz 2000; 11(4): 974–1002.

34.

Akrotirianakis

Rustem

. Globally convergent interior-point algorithm for nonlinear programming. J Optimiz Theory App 2005; 125(3): 497–521.

35.

Pawlowski

Simonis

Walker

. Inexact newton dogleg methods. SIAM J Numer Anal 2008; 46(4): 2112–2132.

36.

Curtis

Huber

Schenk

. A note on the implementation of an interior-point algorithm for nonlinear optimization with inexact step computations. Math Program 2012; 136(1): 209–227.

37.

Zhao

Yang

. Inverse kinematics for the variable geometry truss manipulator via a Lagrangian dual method. Int J Adv Robot Syst 2016; 13(6): 1–11.

38.

Yang

Jiang

. Fast trajectory planning for VGT manipulator via convex optimization. Int J Adv Robot Syst 2015; 12(9): 1–17.