Research Article

INTERNATIONAL JOURNAL OF ADVANCED ROBOTIC SYSTEMS

A method for the automatic generation of inverse kinematic maps in modular robotic systems

International Journal of Advanced Robotic Systems September-October 2016: 1-15 © The Author(s) 2016 DOI: 10.1177/1729881416662790 arx.sagepub.com

1 2 Mohammad Mayyas and Rochelle Mellish

Abstract

Flexible manufacturing based on rapidly reconfigurable robotic systems will enable factories to meet time-sensitive and fast-changing industrial demands. However, with the rise of modular systems there is also the need to quickly and easily determine which configuration is optimal for performing a certain task. In this article, we present a path-based ad hoc technique for determining the inverse and forward kinematics map based on relative joint space variable to reduce the computational complexity. The proposed technique is nonsingular and suits kinematic analysis and optimization of robotic systems with undetermined configuration, and it can be extended to solve generalized inverse kinematic of robotics system involving large number of joint variable.

Keywords

Flexible manufacturing, modular robot, inverse kinematics, automation

Date received: 18 February 2016; accepted: 14 June 2016

Topic: Robot Manipulation and Control Topic Editor: Andrey V Savkin

Introduction

In the ever-changing world of technology, electronic products are continually becoming smaller and smaller, implying the need for the micromanufacturing of system components. The microfactory is an industrial robotic setup that allows for the quick changing of a factory's layout to meet the changing needs of industry.1-3 It is popular to use modular components for the construction of robots, as modular components tend to be independent agents, making them amenable to relocation within the robotic infrastructure.4 While downsizing for the purpose of micromanufacturing and modularizing joints for quick reconfigurability, designers must keep in mind the main concerns of the micromanufacturing industry, that is, the research of minimal actuations factory parts and fabrication accuracy.

The fabrication process in an automated manufacturing setting is normally a set of commands given by a central processing unit to a robotic manipulator to follow a certain trajectory. Therefore, the path that a manipulator must take in Cartesian space is often known before a specific

robotic manipulator is designed. At this stage, the problem of fabrication comes down to selecting a configuration of prismatic and rotary robotic links that will minimize path error and also minimize the amount of space and time needed to complete the task. The most common way to obtain the joint angles needed for a robotic system is to solve kinematic problem. The approach can be applied when both the path and the robotic manipulator are known.

The determination of configuration using inverse kinematic is a method widely used in biomechanics for the analysis of human motion. For example, inverse kinematic

1 Department of Engineering Technologies, Bowling Green State University, Bowling Green, OH, USA

2 Department of Aeronautical and Astronautical Engineering, Purdue University, West Lafayette, IN, USA

Corresponding author:

Mohammad Mayyas, Department of Engineering Technologies, Bowling Green State University, 219 Tech. Building, Bowling Green, OH 43403, USA. Email: mmayyas@bgsu.edu

Creative Commons CC-BY: This article is distributed under the terms of the Creative Commons Attribution 3.0 License (http://www.creativecommons.org/licenses/by/3.0/) which permits any use, reproduction and distribution of the work without further permission provided the original work is attributed as specified on the SAGE and Open Access pages (https://us.sagepub.com/en-us/nam/ open-access-at-sage).

analysis may be used to reconstruct anatomical parameters during some motor tasks such as handwriting or upper limb motor strategies. Vimercati et al.5 observed that participants with Down syndrome tend to draw faster but with less accuracy, or that the motor sequence adopted in a pointing and tapping task may be also affected by the pathology.6 Sale et al.7 adopted an end effector-based robot to control the kinematics of the lower limb in patients with Parkinson's disease in order to improve their mobility. Another example is the work of Ancillao et al.8 who developed a protocol to reconstruct the kinematics of the upper limb and upper body during a drawing/handwriting task.

At the qualitative level, previous work in determining configurations for robotic and kinematic tasks has taken an agent-based approach, where each link is an agent and the objective is to coordinate the behaviors of all agents. For example, in 1996, Carnegie Mellon University's Robotics Institute began a program called the Agile Assembly Architecture, whose purpose was to create an Internet interface through which users could input which robotic components they desired to use; a decentralized cooperative framework enabled each component to be modular, and thus allowed the system to be easily reconfigurable, while maintaining high operability.9-12 In similar work by Chen,13 the tools for developing a reconfigurable robotic workcell were considered (with a workcell being a specific system component or robot). A model for automatic robot generation was presented, which was modular and geometry independent. Optimization of the robotic workcell was specified, and a communication and control framework presented for coordination among the system workcells, with the end goal being a CAD-like program that would enable the appropriate robotic configuration to be determined for a specific task.13

At the quantitative level, Barraquand et al.14 presented a robot path-planning method to handle high degrees of freedom (DOFs). In their method, a graph of the local minima of a potential function was constructed incrementally during the search for the goal configuration. Their method relied on numerical methods for a solution. Optimization from a slightly different approach was carried out by Sny-man and Berner15 who optimized over robotic parameters (link lengths, masses, etc.) rather than positions.

Other work includes the development of inverse kinematic mappings for ''hyper-redundant'' robots. In the study by Chir-ikjian,16 the author presents a ''backbone curve'' technique to compare manipulator configurations. Goldenberg et al.17 present a technique based on the Newton-Raphson numerical method to iteratively find an inverse kinematic mapping of a kinematically redundant robot. On the other hand, one work that closely addresses accuracy concerns is that by Park and Bobrow,18 in which the authors develop methods for optimizing robotic configurations via the ''rigid body guidance problem.'' Fault tolerance, that is, trajectory planning that takes into account the failure of a joint, was proposed by Paredis and Khosla.19 The authors discussed that added redundancy in the form of multiple actuators and joints can aid the robot in

completing a task. Their most important point was that certain joint configurations make it difficult to use the manipulator.

In recent years, neural network-based technique has found various applications for inverse kinematics.20 The work by Li et al.21 proposes a strategy with hierarchical tree generation on topology and nonlinear neural dynamics to control multiple manipulator with only partial command. The method is globally stable and the optimality is guaranteed, which makes it scalable to large number and types of manipulators including redundant DOF system. To solve redundant manipulator under physical constraints, Zhang et al.22 proposed an online approach based on quadratic programming problem formulation for manipulator torque optimization using various levels of acceleration and velocity schemes.

Another approach to solving the inverse kinematics problem for redundant robots was given by Tabandeh et al. ,23 who used a genetic algorithm to solve the inverse kinematics problem for a redundant robot by way of a niching method. They also used a numerical inverse kinematics solver to attain arbitrary accuracies ofjoint angles. Work by Chen2 and Chen and Yang,24 addresses and solves the inverse kinematics problem for redundant serial and tree-type modular robots. Assuming that each link has only 1-DOF, is symmetric to all other links, the Jacobian is used to solve differential kinematics equations for the joint angles using the Newton-Raphson method. The work by Pac et al.25 uses interval analysis in inverse kinematics problem to provide a range ofjoint space variables for a desired end effector position. Wu et al.26 describe a rooted tree in module-level approach to generate kinematics of modular robots based on task-oriented concepts.

Meredith and Maddock27 describe four different types of inverse kinematic mappings, one of which is called the analytical method; analytical methods are described as those that ''reduce the inverse kinematics problem to a mathematical equation that need only be evaluated in a single step to produce a result.''27 One of the drawbacks of analytical methods is that they are often not scalable enough to handle large chains of integrators.27 In many robotics applications, particularly in the case of a modular robot, it is economically efficient to determine the robotic configuration that will be the best for a certain task from an existing pool of components. For a system with many links, determining an inverse kinematic mapping for each possible configuration is time consuming, and limiting to the time-sensitive efficiency of many factories. Thus, an inverse kinematic mapping that can be applied to many different configurations of robotic components is needed. Given a path, a global definition of the translation and rotation necessary to move from one point to another must be reiterated in terms of the joint variables and joint limits.

Many of the concerns for modular robotic systems have been addressed in the literature, particularly in developing globalized inverse kinematic maps whose coordinates can be mapped to various robotic structures. In addition, many continuous inverse kinematic problems based on linearization method have been addressed in literature to compute a

Figure 1. Projections of the center of mass on a hemisphere, with t = 0.1 The black squares represent the paths points p and q , and the lines between points p and q represent the line segments l,. The dotted lines are the end effector vectors e,, and the black circles are the center of mass.

configuration of a manipulator capable of producing a prescribed trajectory in the task space.28-30 However, much of this work relies on the Homotopy technique or Newton-Raphson method, whose success depends on how close the initial guess of the solution is to the actual solution.31 For initial guesses that are close to the actual solution, the run time of the algorithm is low, and Newton-Raphson techniques are ideal. For large joint displacements over small periods of time, in which smaller sampling intervals must be considered to maintain reasonable initial guesses for Newton-Raphson, computational complexity will increase.

Research aim

The first goal of this work is to present an algorithm that estimates the amount of rotation and translation needed to traverse a planar path using the path's geometry. The second goal of this work is to outline a mapping scheme from a simplified robot to an actual 6-DOF manipulator, and then compare the manipulator's performance for different configurations.

The contribution of this article is to present a framework in which the joint angles may be retrieved for motion along planar paths in Cartesian coordinates. This method is different from previously derived inverse kinematics methods in that no knowledge of the robot's structure is needed to determine the joint angles. The method used is a non-gradient-based optimization problem that focuses on changes in the joint angles rather than the actual joint angles of the robots at each point during the task. Since these coordinates may be mapped to the workspace of any robotic configuration, this framework is ideal for reconfigurable robots. The contribution is also to compare the performance of a more

complicated robotic arm to a simplified system, enabling the fast determination of manipulator efficiency along a path.

This article is organized as follows: we present the materials and methods in the ''Materials and methods'' section. We introduce inverse kinematic mapping to find joint angles along a path. Then, we introduce a simple robot with four links, which is used to implement an existing inverse kinematics algorithm. Then, we map the joint angles from the inverse kinematics map to a physical robot by introducing constraint equations. We present the results and discussion in the ''Results and discussion'' section. In this section, we discuss simulation results highlighting the use of the inverse kinematic map and the results of mapping joint angles to a physical robot. Finally, in the ''Conclusion'' section, we offer concluding remarks and discuss future directions for research.

Materials and methods

Path-based determination of the inverse kinematic map

The inverse kinematic mapping of a robot—serial, parallel, and mixed—is based on a method that involves approximating the curve by sampling points along the curve, and using the joint center of gravity and the vectors between two consecutive path points to obtain translation and rotation of the joint. For consecutive movement increments, the amount of rotation and translation is divided among the joints in proportion with the total amount of movement the joint is capable of. Suppose we have a path P 2 R3 parameterized by time, so that P = P(t). At each time step t, the path is specified by a pointp(t) 2 R3,

where t = 1,2,......, n and n is the total number of time-

parameterized data points that comprise path P as illustrated in Figure 1. Let two consecutive points be given byp and q,where

p = p(t) andq = p(t + 1). The homogeneous transformation of the end effector from point p to point q is given by

Tp = q (1)

where p = |p, 1]T, q = [q, 1]T, and where T 2 SE(3) is the homogeneous transformation

et = U x lt

If we let d to be the distance from the path to the center of mass of the final link, then et is given by

j et d Ñ

Theorem 1. Recall that the translation of the mth link's center of mass is given by r, and let us denote the change in the center of mass's location as dg(t) = g(t + 8t)— g(t). Then

lim r- dg(t) = 0

R = R(O) 2 SO(3) is a rotational transformation and r = [rx, ry, rz] is the translation of the center of mass.

Notice that to solve this equation, two quantities must be specified: the translational joint angles, given by r, and the rotational joint angles, given by O. With only knowledge of the path, finding this information reduces to a geometric problem, detailed below.

Translation of the end effector

Suppose that between every consecutive pair of points p and q there is a straight line lt connecting pointp to point q. This line segment has a length |lt |, and an orientation in the base frame O = O(dt, pt, ), where d is the angle of rotation about x, p is the angle of rotation about y, and ^ is the rotation about z. We assume that the end effector of the manipulator is defined in coordinate frame with a vector in R3 denoted by et, is always orthogonal to the line segment lt and axis of rotation for the path U, that is, it lies in the plane of rotation. Therefore, the end effector et originates at the center of mass of the final link and ends at the path point q. To find the line segment et, we assume that the axis of rotation for the path, is known; in fact, it is the axis orthogonal to the plane of the path. With knowledge of the path's axis of rotation (assumed to be parallel to one of the global axes X, y, or z), we can form a 3-axis body coordinate system in which the axis of rotation is the first axis, and lt is the second axis. We would like for the final axis to point inward toward the center of rotation, thus serving as the end effector et and completing the 3-axis coordinate system. Therefore, the body frame's coordinate system is given by b = b(u, lt, et). We have that

Proof. In appendix 1. The results of this proof indicate that using a large number of data points enhances the accuracy of the inverse kinematic mapping. On the other hand, the number of data points is proportional to the amount of time taken for the robot to complete the path. If we let the time needed to compute a single data point be tc, the desired time limit be T, and the time increment be dt, then the proportionality relationship is given by

tpath dftc

Thus, although smaller time increments increase accuracy, they also increase the run time of the robot along the path. The engineer must balance throughput with accuracy by tweaking such parameters. □

Rotation of the end effector

For any given vector in R3, there are 3 rotational DOFs that may be used to characterize its orientation; they are rotation about the x-, y-, and z-axes. The overall pose of the end effector is often a coupling of all 3 rotational DOFs, making the individual determination of all 3 rotational DOFs elusive at best. In this section, we present a procedural and algorithmic approach to determine the rotation of the end effector in the plane.

Since we have approximated the path with line segments between path-points, we may estimate the amount of rotation that the manipulator has undergone using the dot product between two consecutive end effector orientations et+i and et. Letting a represent the rotational angle needed to define the orientation of the end effector, we have

a= arccos

/(ft+bft) \let+illet I

We assume that et acts through point q; thus, defining the center of mass to be the variable g, the location of the center of mass in R3 at point p is given by gt = lt + et. The change in location of the center of mass is what gives us the incremental translation of the manipulator. Hence, we have the translation as

This gives us the amount of rotation in the plane formed by the two vectors et+1 and et.

We measure the accuracy of the robot along the path by comparing the "fit" of the path generated from the inverse kinematics mapping—which is defined by a polygon that is passing through a set of vertices—to the original path, that is, we measure the Euclidean distance between points on the generated path and the corresponding points on the original path. Since error is cumulative, we compare performance at different path parcings by examining the average error per time increment. The error at time t is given by

gt = gt+i - gt

;(t)=p - T(0M0)-T(t)p(t)

-E i"(t)i

(n - k)!

8z. The second rotation in sequence, which we can denote 6x, is about the local x-axis a2 = [1,0,0]T. Subproblem 2 requires finding an intermediate point h between p and q that can be reached by either of the transformations32

T2 (0x)p = h Ti(-ez)q = h

Figure 2. The simple robot in the (1234) configuration shown in the top view and side view. Links 1, 2, and 3 are the x-, y-, and z-directional linear joints, respectively. Link 4 is the 2-DOF rotary joint. The white dot on link 4 is representative of an end effector.

where T(0)r(0) is the initial location of the end effector. Therefore, for a total of n data points, the average error per unit time along the path is

We find the vector, z, that connects h to a predetermined point on the axis of rotation of the first rotation; here, we've chosen that point to be r = [0, 0, 0, 1]T. Thus, z = h — r. We also can define vectors between pointp and r and point q and r as w0 = p — r, and r0 = q — r. Letting a3 = a1 x a2, and then using a linear combination of the triad a1, a2, and a3 gives z as

z = aâi + ftâ2 + ga3 The coefficients a, ft, and g are found by

(âi, ci2)ctT2wo - âivo

a =-2-

((âi,â2)) - 1

Global joint angles from a simplified manipulator

In this section, we demonstrate the use of joint angles from an inverse kinematic map using a simple robot. We consider a robotic manipulator that consists of three linear joints: one that traverses the x direction, one that traverses the y direction, and one that traverses the z direction. It also consists of one rotary joint that is assumed to rotate about the local x- or z-axis. Such kind of robots essentially boil down to the inverse kinematics of a spherical dyad. The four links are shown in Figure 2. There are 64 configurations possible from these components, that is

(âi, â2)â[vo

((âi, ai2))

2 |wo| - a2 - ft2 - 2aft(âi, â2) g =-^-T"[2-

|âi X â2|

Using z = hi — r, we can solve for the point h. To solve for the joint angles, we follow the steps outlined in Murray32 by solving subproblem 1 for dx and 0z. To find 0z, if we define

w = wo

âiâi wo

âiâivo

configurations possible from the four links, and we compare them based on average error and computational cost, where n is the original set of joints. For a given path we measure the error in the manipulator's performance under the inverse kinematic map based on the equation (10).

We obtain the joint angles from the paths using the well-known Paden-Kahan subproblems produced by Murray et al.32 Since there are essentially only two revolute joints, we determine joint angles using Paden-Kahan subproblem 2, featured. The joint angle trajectories are determined from the path by applying subproblem 2 repeatedly to consecutive points along the path. Again, let the first point sampled be denoted by p, and the second point by q. We wish to find the amount of rotation and translation that has occurred between both points during the time increment dt. Starting with rotation, we determine the axis of revolution for the first rotation in sequence to be about the z-axis, whose unit vector is given by ai = [0,0, 1]T; we denote this angle by

An identical procedure may be used to find To find the amount of translation that has occurred between two points, we first use homogeneous transformations using the two rotational angles on the point p, and then subtract the result from the point q. Thus,

vv = vvo

then the rotation may be found as

, /âT(w x v) 80z= tan-1' 1 ( v

cos(ftz sin(0z

- sin(6>z)

cos($z)

o cos(6>x o sin( ex)

- sin($x) cos(0x)

so that Ti =

Ö3X1 1

and T2

R2 0i 3

03X1 1

we have that the translation is given by d = q — T1 T2p.An example using this method is given in the ''Results'' section.

Mapping joint angles to a physical system

The variables that we wish to find at each time step are the changes in the robot's joint angles. The robot's movement in 3-space can then be defined as f = f (8x, Sy, 8z, 80), where 8x, 8y, 8z, and 80 are the joint variables associated with each of the manipulator's links; each variable can be a multidimensional vector depending on the redundancy of links exhibiting such motion. We define the global position variables obtained from the inverse kinematic mapping to be 8x giobai, 8y giobai, and 8z giobai, while the global rotational displacement is given by 8g global. Constraints between the physical robot and the three components of Cartesian space are given by equating the global position variables to the projections of the robot's motion onto the three ordinal directions. That is, for the inertial unit vectors x, y, and z

8x global = fR ' x 89 global = fR • 9 8z global = fR • Z

We also represent constraints on the revolute motion of the robot as

global

where m is the number of revolute variables affected by the constraint, g global is the constraint condition, 80,- is the i th revolute variable, and is a weight.

To determine the optimal solution angles based on minimum path error between time t and t + 1, we propose constrained optimization, with the Lagrangian33

L = /r(8x, 8y, 8z, 80) + Ai (8xglobal -/r • X)

+ A2(8y global - fR • >0+ A3(8z global - fR •

+ A4(8g global -J2

By taking partial derivatives of the Lagrangian with respect to each variable and setting the resulting equations equal to zero, we have

d8x dL

dL @8y dL d\~2

d8z dL dÄ3

dL d 89 dL d\A

With this framework, we are able to solve for each of the desired joint angles of the target robot.

Figure 3. 6-DOF robot. Joints that enable rotation about the local z-axis are denoted by the cylindrical prisms, while joints that enable rotation about the local x-axis are denoted by circles. 6-DOF: 6 degrees of freedom.

Example using a 6-DOF manipulator

Inverse kinematics using the simplified robotic structure gives us global definitions of the amount of rotation about the x- and y-axes, and translation in all three directions. The next step is to map these globally referenced joint variables to the 6-DOF robotic manipulator. Unlike the simplified representation, only 6-DOF robot has rotary links.

To accomplish the mapping, we note that we have two rotational DOFs centered at the same center of rotation; during our comparison process, we sampled different configurations of the simplified robot, and moved the center of rotation to various positions within the robot's structure. We were able to compare which rotation center produced the most accurate following of the three-dimensional (3-D) path. As can be seen from Figure 3, the manipulator has six joints like the center of rotation for the simplified robot, each successive pair of joints can rotate in about two directions. Thus, in the mapping process, we can choose which two joints we would like to be the center of rotation and assign the two rotational joint angles to that pair. For the remaining joint angles, a combination of the unused rotational DOFs would need to be used to produce the needed translational

motion. Let 80,- (i = 1,......6) represent the i th joint

angle of the 6-DOF manipulator, and let the optimal joint angles be given by 80* = [80*, 80*, 80*, S04*, S05*, 80*]T. Choosing 80* to be fixed and the manipulator's 5th and 6th rotational joints, 805* and 806*, respectively, to be the center of rotation, we may remove these angles from the vector 80. For the remaining rotational DOFs, we have the following constraint relationships

8xglobal = 12 cos(892) + l3 cos(803) + l5 cos(805) 8y global = I2 sin(802) + I3 sin(803) + I5 sin(805 )

We now use constrained optimization to determine the remaining joint. The Lagrangian is given by

L = 80 + A1 ^8xglobaj — /2 cos(802) — /3 cos(803) — /5 cos(805)^

+ A^8yglobal — /2 sin(802)— /3 sin(803)— /5 sin(805) j + A3(0B — 806 — 804 — 801) + A4(0a — 802 — 803 — 805)

where 0A and 0B are the two rotational DOFs retrieved from the all DOF robot, and A1, A2, A3, and A4 are Lagrange multipliers. The partial derivatives of the Lagrangian are given by

9L 9L 9L

980! 9806 9804 9L

= 1 - A3 = 0

9802 9L

9803 9L

9805 9L 98A1 9L 98A2 9L 98A3 9L 98 A4

= 1 + A1 /2 sin(802)- A2/2 cos(802)- A4 = 0 = 1 + A1 /3 sin(803)- A2/3 cos(803)- A4 = 0 = 1 + A1 /5 sin(805)- A2I5 cos(805)- A4 = 0 = 8x — /2 cos(802)- /3 cos(803)- /5 cos(805) = 0 = 8y — /2 sin(802) — /3 sin(803) — /5 sin(805) = 0 = 0B — 806 — 804 — 801 = 0 = 0A — 802 — 803 — 805 = 0

Suppose we choose joints 5 and 6 to be the center of rotation for the manipulator, and since joint 1 is at the base of the manipulator, we fix it so that it has zero rotation at all times. We can then reduce the size of the problem to one in which we have only three unknown joint angles to solve for. Thus, along with the Lagrange multipliers, the optimization variables we want to find are given by 80 = [802,803, 804, A1; A2, A3, A4]t, and the new reduced system to be optimized is

/1 = 1 + A1/2 sin(802) — A2/2 cos(802) — A4 /2 = 1 + A1/3 sin(803) — A2/3 cos(803) — A4

/3 = 1 — A3

/4 = 8xglobal — /2 cos(802)— /3 cos(803) — /5 cos(805) /5 = 8yglobal — /2 sin(802)— /3 sin(803)— /5 sin(805) /6 = 0B — 806 — 804 — 801 /7 = 0A — 802 — 803 — 805

The procedure for finding the core set of partial derivatives when the center of rotation is either links 3 and 4 or 1 and 2 is a similar process of removing those equations for the joint variables that are already known.

Equation (29) must be solved at every time step, and they have no closed form solution. Performing a Newton-Raphson technique at every time step could be computationally time consuming when path is discretized into very small step size. Thus, since we are taking small increments of movement and time, we linearize the equation (29) about the manipulator's current pose, and solve the linearized system. Let "•" denote the dot product, A = [A1, A2], and K(80,-) = [cos(80,-), sin(80,-)]. With the Jacobian matrix

/2A • K(802) 0 0

/2 sin(802) —/2 cos(802) 0

/3A • K(803) 0

/3 sin(803) —/3 cos(803) 0

/2 sin(802) /3 sin(803) 0

/5 sin(805) —/5 cos(805) 0 0

—/2 cos(802) —/3 cos(803) 0 0 0 0 0

—1 —1 0 0 0 0 0

The linearized system is given by J80* = 0. Thus, we survey the null space of equation (30) to find our optimal solution 80* = [802*, 80*, 804*, A1, A2, A3, A4]T at each time step. To show that the joint angles obtained from the null space of the linearized system are suitable for equation (31), we have the following theorem.

Theorem 2. As the size of the state 80 * ! 0, the linear approximation of the partial derivatives given in equation (31) approaches the nonlinear partial derivatives of equation (30).

Proof. In Appendix 2. For sufficiently small 80, we can use the null space of the Jacobian matrix J to find the joint angles 80*. This raises the question of what is "sufficiently small'' for the constructed framework, and how do we measure whether or not the joint angles commanded by the optimization fall within limits for which the linearization is valid? These questions lead to the next theorem. □

Theorem 3. The Jacobian matrix J given by equation (31) is evaluated at the current value of the state 80. That is

dfiL d 80

J 80 = J (80 *)80 *

Therefore, the linearization is valid for arbitrarily large 80 *.

Proof. In appendix 3. Upon finding the optimal joint angles, we substitute them into the forward kinematic map for the 6-DOF robot. If we let the matrix describing local rotations about the y-axis be given by

Ry(80i)

cos(80i) 0 sin(80i ) 0 1 0 - sin(80i) 0 cos(80i )

and the matrix describing local rotations about the z-axis be given by

Rz(80i) =

cos(80,) - sin(80,) 0 sin(80i) cos(80i) 0 0 0 1

and furthermore, let the center of mass of the i th link be given by di, the respective homogeneous transformations for rotary links about the y-axis and about the z-axis are given by

Ty (80,) =

and Tz (80i) =

the simple robot. In particular, we have included the initial location of the simple robot's end effector as T simple(0)p(0) and subtracted it from q(t) to ''realign'' the path traced by the physical robot with the commanded path given by the equation. Likewise, the error measured by the 6-DOF manipulator is given by

"6 dof (t) = q(t) — T6dof(0M0) — T6DOF(t)p(t) (38)

Let £ simple = Tsimple(0M0) and £6DOF =T6 DOF (°M0).

With these substitutions, and dropping the (t) from the defined variables for simplicity, we may solve equation (37) for the position p as

e(q - Cs

' simple )

simple (q — £ simple

Substituting equation (39) into equation (38) for p(t) gives

"6 DOF = q - T6 DOF Tsimple (q - C simple

£ simple ) C6 DOF

Thus, if p =[x,y, z, 1]T is the current location of the robotic manipulator, the transformation from the manipulator's frame (m) to the robot's base frame (b) is given by

Thm (0) = Ty(801*)Tz(802*)Tz(803*)Ty (80l)Tz(805s)Ty(80l)p

Since the system is linearized, using small time increments would be the most advantageous for increased accuracy along the path. □

Theorem 4. An optimal solution always exists for the system J 80,* = 0.

Proof. Since any system Ax = 0 has at least the solution x = 0, the system (35) always has a solution. Therefore, the inverse kinematic mapping based on the relative joint space method has no singularities.

Transforming error of the simple robot to the 6-DOF manipulator

What if, given the average error per unit time of the simplified robot, we would like to predict the average error of the 6-DOF manipulator? Based on equation (13), we define the error measured by the simple robot to be

" simple (t) = q(t) — Tsimple(0)p(0) — Tsimple (t)p(t) (37)

where p(t) and q(t) are the beginning and end points of the travel segment, Tsimple(t) is the homogeneous transformation from the end effector's frame to the base frame of

This transformation gives a mapping at time t from the error measured by the simple robot to the error measured by the 6-DOF manipulator. Since the transformations Tsimple and T6DOF are in the group SE(3), they are guaranteed to be invertible for all time t.

Cost function

Different configurations have dissimilar path-following attributes that may make certain configurations more desirable than others. For instance, the more a linear joint moves along its DOFs, the more actuation is needed to produce the motion. Also, the heavier a joint, the more stress on the actuators to overcome the joint's inertia. Finally, the most important factor—how well the robot, using the inverse kinematic map, follows the path, or the accuracy of the robot along the path— must be considered in a cost function. In this section, we focus on formulating a cost function that encapsulates all of these concepts in order to aid us in identifying which configuration has the most efficient behavior along a given path.

The cost function is composed of the joint variables, the inertia penalty, and the path error. We choose an energylike cost function so that it will be amenable to control design. To obtain the cost, we use

C = 2 (XTMxX + y Myy + zTMzz + Ip + Ex2 + E2 + El)

where x, y, and z are the n x 1 vectors containing the x-y-z variables for each joint of the robot, Mx, My, Mz are positive semi-definite weighting matrices, where we have chosen to set each equal to

Mx = My = Mz =

m1 0 0 0

0 m2 0 0

0 0 0

0 0 0 m

and is the inertia penalty. Let 0 be the n x 1 vector containing the angular displacement for each of the n joints, let 11,12, ..., In be the inertia of each joint about its own center of mass, and let

Ir,i(t) Ir,2(t)

miE di,k (t)2

E^2,k (t)2

Ir,n(t) = m^Edn,k (t)2

where di*, ¿2. k, ... , dn ,k are the time-varying distances of the joints from each of the K centers of rotation. Note that only joints between the rotary joint and the path will be affected by the rotary inertia. Placing

Figure 4. Tracing of a curved planar sinusoidal trajectory using the proposed inverse kinematic mapping. The time increments of 8t = 0.1 and 8t = 0.001 give average path errors of e = 0.1609 and e = 0.0016, respectively. The error cost function is linearly proportional to the time increment, and both have the same order of magnitude.

all of these components together, we obtain the inertia penalty as

'h o o o o o o

o 12 o o o E,2 (t) o o

0(t) h0(t)T 0(t)

o o o o o o

o o o In. o o o 1r,n (t) -

Ex, Ey, and Ez are the x-, y-, and z-projections of the path error, given by

Ex = (E, x) (45)

Ey = (E, y) Ez = (E, £)

where E is the difference between the actual pathp actuai and the path generated from the inverse kinematics joint variables p IK , given by

E = p actual - PI

Results and discussion

In this section, we present the results of simulations using the proposed methodology for determining two-dimensional (2D) path inverse kinematic maps, and for the case study of multiple configurations of a simplified robot whose joint angles have been mapped to a 6-DOF robotic arm. We implemented the algorithms in MATLAB (T. M. MATLAB 8.5, Inc., Natick, Massachusetts, United States).6

2-D path generation along a planar sinusoid—error cost function

The planar path that we consider is given by

x(t) = cos(t) y(t) = sin(t) z(t) = o

where o < t < in steps of size 8t. The simplified robot is in the (1234) configuration in Figure (3). Following the analysis in section ''Materials and methods,'' we assumed that the lengths of the three linear links are Lx = Ly = Lz = o.i m.

The results are shown in Figure 4, and confirm that smaller time increments—which correspond to small 80 increments in Theorem 3—produce higher accuracy between the generated curve and the curve produced using the output of the inverse kinematic map.

2-D path generation along a planar sinusoid— Generalized cost function

To test the inverse kinematic mapping, we use a three stages ''joints'' robot. Each slot is stacked vertically (along the y-direction). We have two linear joints, the first with one redundant DOF (x) and the other with two DOFs (x-z), and a rotary joint. We make the simplifying assumption that the rotary joint is a disk, and the two linear joints are rectangular plates as shown in Figure 5. We will use this system to test the inverse kinematic mapping along a sinusoid. The joint dimensions, movement capabilities, and slot information are shown in Table 1.

For this exercise, we study the robot in the 1-2-3 configuration, whose forward kinematic mapping is given by

19 21'

cos(6>3) 0 - sin(6>3) Su + S21 + S31 + «i + «2 + ^7: cos(03)+ —

400 80

tO (0) =

S12 + S22 + S32

sin(03) 0 cos(03) S13 + S23 + S33 + d2 + 4UU sin(03) 0 0 0 1

Figure 5. Physical representation of each stage "joint".

We can already see that for the selection of joint types and their slot orientations, the path of the end effector will always be planar, regardless of configuration.

The general equation for the sinusoid that we use is

z = p start + A sin(

where p start is the insertion point in the x-z plane of the sinusoid function, A is the amplitude of the sinusoidal wave, U is the frequency, and B is the wavelength. Using the position in TeO(0), we were able to identify values of xt and zt at t = 0 (given by ^1jX and ^1jZ, respectively) and solve forp start. With the parameters A = 0.01, B = 0.5, and U = 2-fl-, we determined p start to be given by the equation

p start = ,P1,z - A sin[ B_P1,x

differences in cost can be seen for the six configurations that use only two of the joints and the three configurations that use only one of the joints. The reason for this pattern in the six configurations that use all three joints may be that regardless of the joint sequence, the same net movement is undertaken by the robot. This seems to indicate that the order in which parts are assembled does not matter, and could be due to the commutative nature of many of the joints involved and the fact that the approximate inverses always have a "trivial" unique solution. Increases in the inertia penalty due to configurations that place the rotary joint in either the first slot or the second slot are leveraged by smaller translations made by the joints. The ''sweeping motion'' of the rotary joint, particularly when the rotary j oint is further in position from the path, eliminates much of the needed translation.

The costs for two-joint configurations is largely dependent on the joints being used. For instance, configuration

1-3 has a much lower cost than configuration 2-3 because

2-3 cannot traverse the entire sinusoid, and therefore accrues much more error than configuration 1-3. Configurations 1-2 and 2-1 actually have fairly low costs, even without the help of the rotary joint, because they are independently able to traverse the entire sinusoidal curve.

3-D path generation along a curved sinusoid

To test the efficiency of the various configurations of the simple robot, we consider using a curved sinusoid, given by the parametric representation

where ^1jZ is the z-component of the insertion point, and is the x-component.

Table 2 and Figure 6 compare the accuracy of the inverse kinematic map along the sinusoidal curve for increasing numbers of data points.

We use the 1-norm for comparison because using a single number allows us to quickly see the net effect of error over time. Taking the 1-norm of the difference between the desired and actual paths confirms that increasing the number of data points increases accuracy of the estimated path. Using the automatic-inverse kinematic map generation algorithm, we were able to compare n! configurations for the three-joint robot. The results are shown in Table 3.

Under the chosen cost function, all of the six configurations that use all three joints have similar costs, while marked

x(t) = cos(t) y(t) = sin(t) z(t) = sin(15t)

where ^/3 < t < 2^/3, with increments of 8t = 0.001. We tested both the simple robot and the 6-DOF manipulator along the path. Starting with the simple robot of four groups of configurations shown in Figure 7, we chose to make it both massless and dimensionless, that is, all motion is concentrated at a single moving point. All 64 configurations were tested. Configuration sets that have the same cardinality (4,3,2, or 1 link(s)) also have the same average error, because of the small increments in movement taken. When infinitesimally small changes in the joint angles are used in majority of simplified manipulator, the end effector's motion is independent of the sequence of links. The findings in Table 4 indicate that differences in average error between configurations depends on

Table 1. Part and slot information.

Part no. Range Symbol End effector Dimensions Mass

1 x = 0- -1.0 m al 0.125 m 1:1.395 m; w:0 .225 nl; h:0.l00 m 44.l kg

2 x = 0- -0.15 m; z = 0-0.15 m a2; d2 0.1375 m 1:0.275 m; w:0.l30 m; h:0.04 m 4.6 kg

3 0 = 0- -2- rad 03 0.0475 m h:0.055 m; r:0.0475 m 2.0 kg

Slot no. Rel. Loc. (m) Symbol Rel. Rot. Relative to:

1 [0, 0, 1.950] [Sll, Sl2, Sl3] I3 x 3 Origin

2 [0, 0.15, 0] [S2l, S22, S23] I3 x 3 l

3 [0, 0.15, 0] [S3l, S32, S33] I3 x 3 2

Table 2. Comparison of accuracy with 10, 20, 40 and 100 data points.

Number of data points l-Norm of position error

l0 0.8609

20 0.3687

40 0.l268

l00 0.0927

error of the simple robot at time t, the end point q of a path segment at time t, and the initial positions £ simple and £6DOF of each robot's end effector. To generate Figure 8, the result of equation (40) was subtracted from the original path's coordinates to obtain the actual path followed by the 6-DOF robot. As expected, since we use all four links of the simple robot, the error is relatively small, where there is no difference observed between the plotted curves in the graph.

the number of links used. The same idea holds true for computational efficiency. The number of computational operations depends primarily on the implementation of the algorithm, and this number is independent of permutation when configurations of the same cardinality are used. Between cardinalities, the difference lies in the length of the chain of exponentials used to compute the forward kinematic map. For the multiplication of two matrices A and B of size b x n and n x p, matrix multiplication alone has bp (2n — 1) operations (often referred to as ''flops'').

Depending on how the forward kinematics algorithms are implemented, matrix computations can be the most costly. Thus, the programmer must weigh the cost of production, judging between accuracy and computational efficiency.

The patterns in accuracy and computational efficiency of the simple robot can be extended to the 6-DOF manipulator when considered as separate cases. However, generating a correspondence relationship between the simple robot's four links and the links of the 6-DOF manipulator is elusive at best, since several of the 6-DOF robot's joint angles contribute to multiple global DOFs. To see an example of this, suppose that we set 8x global to zero to represent the removal of the simple robot's x-directional link. The unevaluated Jacobian matrix equation (31) will still be the same as if we had not. The only way to alter the Jacobian would be to set 8x global and 8y global to zero (removing the x-and y-directional links of the simple robot). Doing this would leave the equation 0B — 806 — 804 — 80i = 0, and there are an infinite number of solutions to this equation.

Mapping error from the simple robot to the 6-DOF manipulator is more reliable. Using equation (40) and only considering the case when all four links of the simple robot and all six links of the 6-DOF robot are used, the results of the error mapping are shown in Figure 7. The inputs were the forward kinematic maps of the simple and 6-DOF robots, the

Conclusion

In this work, we have described a computationally efficient joint space algorithm that determines the configuration of a modular robotic system for a given path. We presented a framework in which the joint angles may be retrieved for motion along planar paths in Cartesian coordinates, where we reconstruct robot configurations using the change in the joint angles rather than the absolute values. To accomplish this, we first determine the inverse kinematic mapping based solely on the path traveled; then, applying the inverse kinematic map to all configurations of the robot, we compare each configuration based on accuracy in following the path, and minimal computation required. The developed algorithm allows users to input information about the joints of a robotic platform and receive an automatically generated forward kinematic mapping for any configuration of the joints, as well as the time-dependent joint variables along a given path. The generalized forward kinematic mapping can handle a set of distinct robotic joints as well as duplicate joints and still produce the appropriate mapping by the utilization of numerical part markers that the algorithm uses to query information about the joints considered. Using this inverse kinematic mapping has proven to render accurate reproductions of a path with bounded curvature and without singularity when the measured joint variables were placed through the forward kinematic mapping.

We also successfully map joint angles obtained from an existing 3-D inverse kinematics method to a physical 6DOF manipulator and set up the framework for comparing multiple configurations of this robot. Although such simple manipulator may have closed form solution, but our method is suitable for solving inverse kinematic of robotics system involving large number of joint variable, example may include self-assembly of molecular scale robots. This

Figure 6. (a) 40 data points and (b) 100 data points.

Table 3. Comparison of 1-norm of costs for all configurations when n — 3.

Permutation 1-Norm of cost

(3-i-2)>(3-2-i)>(i-3-2)>(i-2-3)>(2-i-3)>(2-3-i) 546.3976

(1-3), (3-1) 550.959!

(2-3^(3-2) i2.3874

(2-i)> (i-2) 556.i739

(3) I6.42I6

(1) 537.3598

(2) 12.4337

new comparison metric will enable designers and micro-machinists to build manufacturing platforms with prior knowledge of the configurations that produce the optimal output for the task, hence speeding up the rapid prototyping process significantly.

In future work, we will extend the two dimensional inverse kinematics scheme to three dimensions in order to completely specify a mechanism for the fast determination of optimal configurations in modular robotic platforms. We will also develop a cost function to more fully characterize the tradeoff between computational efficiency and accuracy.

Figure 7. Schematic of simple robot with four groups of configurations.(a) Group of one rotational joint, (b) group of two rotational joints, (c) group of three rotational joints, and (d) group of four rotational joints.

Table 4. Four groups of configurations for the simple robot based on the number of links. The time increment is St = 0.001.

Average error (mm),

Assigned number Configurations (in order) computational operations

1,2,3,4,5,6,7,8,9,10,11,12,13; (I432),(I423),(I342),(I324), (I234),(I243), (4I32),(4I23), 0.0000I4344, 47I

14,15,16,17,18,19,20,21, (43I2),(432I),(423I),(42I3), (34I2),(342I),(3I42),(3I24),

22,23 ;24 (32I4),(324I), (243I),(24I3), (234I),(23I4),(2I34),(2I43)

25,26,27,28,29,30,31,32, (34I),(3I4),(43I),(4I3),(I43), (I34),(24I),(2I4),(42I),(4I2), 0.002I54348, 407

33,34,35,36,35,37,38, (I42),(I24),(23I),(2I3),(32I), (3I2),(I32),(I23),(234),(243),

39,40,41,42,43,44,45,46,47,48 (324),(342),(432),(423)

49,50,51,52,53,54,55,56,57,58,59,60 (4I),(I4),(3I),(I3),(2I),(I2), (34),(43),(24),(42),(23),(32) 0.00994306I, 343

61,62,63,64 (I),(4),(3),(2) 0.009929444, 279

Figure 8. Comparison of path following given the original path p(t), the path generated by the 6-DOF manipulator's inverse kinematics map, the path generated by p(t) — e6DOF(t), and the path generated by p(t) — e simpie(t). 6-DOF: 6 degrees of freedom.

Acknowledgement

The authors would like to thank Dr Woo Ho Lee and Dr John Shin from the Robotics Division of the University of Texas at Arlington Research Institute for the valuable feedback. The authors declare that there is not conflict of interest that may have affected the outcomes of this work.

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) disclosedreceipt of the following financial support for the research, authorship, and/or publication of this article: This work was partially funded under a grant from Office of Naval Research.

References

1. Mishima N, Nakano S, Ashida K, et al. Development of a downsized modular manufacturing system and its efficiency analysis. In: 42nd C/RP Conference on ManM/actMring Systems SUstainable Development o/ManM/actMring Systems, Grenoble, France, 24-26 June 2009.

2. Chen I-M, Yang G and Kang I-G. Numerical inverse kinematics for modular reconfigurable robots. JoMrnal o/Robotic Systems 1999; 16: 213-225.

3. Chen I-M. Modular robots. In: Nee AYC (ed) Handbook o/ manM/actMring engineering and technology. London: Springer, 2015, pp. 2129-2168.

4. Horbach S, Ackermann J, Muuller E, et al. Building blocks for adaptable factory systems. Robot CompMt/ntegrManM/2011; 27: 735-740.

5. Vimercati S, Galli M, Stella G, et al. Clumsiness in fine motor tasks: evidence from the quantitative drawing evaluation of children with Down Syndrome. J/ntellDisabilRes 2015; 59: 248-256.

6. Vimercati SL, Galli M, Rigoldi C, et al. Motor strategies and motor programs during an arm tapping task in adults with Down syndrome. Exp Brain Res 2013; 225: 333-338.

7. Sale P, De Pandis MF, Sova I, et al. Robot-assisted walking training for individuals with Parkinson's disease: a pilot randomized controlled trial. BMC Neurol 2013; 13: 1.

8. Ancillao A, Galli M, Vimercati SL, et al. An optoelectronic based approach for handwriting capture. Comput Meth Prog Biomed 2013; 111: 357-365.

9. Gowdy J and Rizzi A. Programming in the architecture for agile assembly. In: Robotics and Automation, 1999. Proceedings. 1999 IEEE International Conference on, Detroit, MI, 1999, vol. 4. pp. 3103-3108. doi: 10.1109/ROBOT.1999.774070.

10. Rizzi AA and Hollis RL. Opportunities for increased intelligence and autonomy in robotic systems for manufacturing. In: Shirai Y and Hirosw S (eds) Robotics research. London: Springer, 1998, pp. 141-151.

11. Kume S and Rizzi A. A high-performance network infrastructure and protocols for distributed automation. In: Robotics and automation, 2001. Proceedings 2001 ICRA. IEEE international conference on, 2001, pp. 3121-3126.

12. Hollis R, Rizzi A, Brown H, et al. Toward a second-generation minifactory for precision assembly. In: International advances robotics program workshop on microrobots, micromachines and microsystems, Moscow, Russia, 2003.

13. Chen I-M. Rapid response manufacturing through a rapidly reconfigurable robotic workcell. Robot ComputlntegrManu/ 2001; 17: 199-213.

14. Barraquand J, Langlois B and Latombe J-C. Numerical potential field techniques for robot path planning. IEEE Trans Syst Man Cybernet 1992; 22: 224-241.

15. Snyman J and Berner D. The design of a planar robotic manipulator for optimum performance of prescribed tasks. Struct Optim 1999; 18: 95-106.

16. Chirikjian GS. Hyper-redundant manipulator dynamics: a continuum approximation. Adv Robot 1994; 9: 217-243.

17. Goldenberg A, Benhabib B and Fenton RG. A complete generalized solution to the inverse kinematics of robots. IEEE J Robot Autom 1985; 1: 14-20.

18. Park FC and Bobrow JE. Geometric optimization algorithms for robot kinematic design. J Robot Syst 1995; 12: 453-463.

19. Paredis CJ and Khosla PK. Fault tolerant task execution through global trajectory planning. Reliab Eng Syst Saf 1996; 53: 225-235.

20. Li S, He J, Li Y, et al. Distributed Recurrent Neural Networks for Cooperative Control of Manipulators: A Game-Theoretic Perspective. IEEE Trans Neural Netw Learn Syst 2016; PP(99): 1-12.

21. Li S, Cui H, Li Y, et al. Decentralized control of collaborative redundant manipulators with partial command coverage via locally connected recurrent neural networks. Neural Comput Appl 2013; 23: 1051-1060.

22. Zhang Y, Ge SS and Lee TH. A unified quadratic-programming-based dynamical system approach to joint torque optimization of physically constrained redundant manipulators. IEEE Trans Syst Man Cybernet B Cybernet 2004; 34: 2126-2132.

23. Tabandeh S, Melek WW and Clark CM. An adaptive niching genetic algorithm approach for generating multiple solutions of serial manipulator inverse kinematics with applications to modular robots. Robotica 2010; 28: 493-507.

24. Chen I-Ming and Yang Guilin. Inverse kinematics for modular reconfigurable robots. Robotics and Automation, 1998. Proceedings. 1998IEEE International Conference on,Leuven, 1998. vol. 2, pp. 1647-1652. doi: 10.1109/ROBOT.1998.677390.

25. Pac MR, Rakotondrabe M, Khadraoui S, et al. Guaranteed Manipulator Precision via Interval Analysis of Inverse Kinematics. In: Proc. of ASME 2013 International Design Engineering Technical Conferences, Portland, Oregon, August 2013.

26. Wu Wenqiang, Guan Yisheng, Li Huaizhu, et al. Task-oriented inverse kinematics of modular reconfigurable robots. In: 2013 IEEE/ASME International Conference on Advanced Intelligent Mechatronics, Wollongong, Australia, 9-2 July 2013, pp. 1187-1192. IEEE.

27. Meredith M and Maddock S. Real-time inverse kinematics: the return of the Jacobian. Technical Report No. CS-04-06, Sheffiel: Department of Computer Science, University of Sheffield, 2004.

28. Sentis L. Compliant control of whole-body multi-contact behaviors in humanoid robots. In: Harada Kensuke, Yoshida Eiichi and Yukoi Kazuhito (eds) Motion planning for humanoid robots. London: Springer, 2010, pp. 29-66. ISBN 978-184996-220-9.

29. Chang P-H. A closed-form solution for inverse kinematics of robot manipulators with redundancy. IEEE J Robot Autom 1987; 3(5): 393-403.

30. Siciliano B. A closed-loop inverse kinematic scheme for online joint-based robot control. 1990; Robotica 8: 231-243.

31. Deuflhard P. Newton methods for nonlinear problems: affine invariance and adaptive algorithms. Vol. 35. Berlin: Springer Science & Business Media, 2011.

32. Murray RM, Li Z, Sastry SS, et al. A mathematical introduction to robotic manipulation. Boca Raton: CRC press, 1994.

33. Rockafellar RT. Lagrange multipliers and optimality. SIAM Rev 1993; 35(2): 183-238.

Appendix 1

Al. Proof of theorem 1

Proof. Let et = \et 1]T. Then we have

r - Sg(t) = (q(t)- T(e)p(t)) - (e(t)+ et+At - p(t)— et) = -T(e)jr(t)+ p(t)- et+At + et

= p(t)(I - T(e)) + et - et+At

=p(t)(yI - T(e)) - l(t) =¿(t)(I - T(e)) - (e(t)-p(t)) =p(t)(yI - t(e))- (p(t + so-zKt))

WeknowthatasSt ! 0,p(t + St) !p(t),whichimplies that p(t + St)-p(t) ! 0. To show that p(t)(j - T(e)) vanishes, recall that

'het, et+At)\

arccos

\et\\et+At 1/

As St —► 0

{et, et+At) ; {et, et) ; Iet|2 = 1

Ietj| et+At1 |et ||et1 |et |2 ! S0 = arccos(1) = 0

! T (Se) = T (0) = I

! p(t)(I - T(se)) = p(t)(I -1) = 0

This completes the proof. A2. Proof of theorem 2

Proof. Let dfi for i = 1,......, 7 represent the Taylor

series expansion of the system (34), where we expand

each equation 1,......, 7 separately. If we consider

equilibrium to be 0 = 0 (that is, when the linearization is theoretically equivalent to the nonlinear representation), then performing the Taylor series expansion on equation (30), and choosing the scenario in which 0,* = 0 to be the reference condition, we have

df1 = l2 sin02 + S02/2A1 cos02 - SA2/2 cos02

+ S02A2/2 sin02 - SA4 = 0 df2 = SA1/3 sin03 + A1/3S03 cos03 - SA2/3 cos03

+ A2l3S03 sin03 - SA4 = 0 df3 = -SA3 = 0

df4 = l2S02 sin02 + l3S03 sin03 + Z5S05 sin05 = 0 f = -l2S02 cos02 - l3S03 cos03 - l5S05 cos05 = 0 @f6 = -S04 = 0 df7 = -S02 - S03 = 0

Plugging in 0i* = 0 for the appropriate 0i gives

f = S02/2A1 - SA4 = 0

@f2 = A1/3S03 - SA2/3 - SA4 = 0

df3 = -SA3 = 0

df4 = l5S05 sin05 = 0 (2.2)

@f = -S02/2 - l3S03 - l5S05cos05 = 0 @f6 = -S04 = 0

df7 = -S02 - S03 = 0

For the linear system, we may use the Jacobian matrix given in equation (31), to form the differential equation

dS0 dt

Multiplying out the terms of this equation gives d/,L

d/,£ dt

d/4,L dt

d/5,L dt

d/6,L dt

■ = cos02 + 802/2A2 sin02 + 8A1/2 sin02

—8A2/2 cos02 — 8A4 = 0 = 803/3A1 cos03 + 803/3A2 sin03 + 8A1/3 cos03

—8A2/2 cos02 — 8A4 = 0

= —8A3 = 0

= 802l2 sin02 + 803l3 sin03 + 8A1l5 sin05 = 0 = —802/2 cos02 — 803/3 cos03 — 8A1/5 cos05 = 0 = — 804 = 0 = — 802 — 803 = 0

Since Ai and 05 are known, any term pre-multiplied by them goes to zero.

If we let @f represent the i th partial derivative of equation (30), where i = 1,......, 7, and d/ijL/dt represent

the i th differential equation of equation (38), then we have that

3/(80)-/ ! 0as80,- ! 0 (2.6)

Assuming that 80 * = 0 is our reference condition, this is

t'80=080 = 0

Again, plugging in 0* = 0 for the appropriate 0,- gives

df1,L dt

However, this is precisely equation (2.7). This completes the proof.

A3. Proof of theorem 3

Proof. Considering that the first-order linear approximation of a system near the reference state 80 * is given by

d/4,L dt

d/5,L dt

d/6,z dt

d/?,L dt

802/2A1 — 8A4 = 0 /1(80)" /1(80 T

A1/3803 — 8A2/3 — 8A4 = 0 ../7(80). /7(80 *).

• 90 (3.1)

—8A3 = 0 l5805 sin05 = 0 —802l2 —13803 —804 = 0 —802 — 803 = 0

We can readily see that for our framework, 80 = 80 *. Hence we are left with the zeroith-order term for (2.5) fi,......f ]T .That is

2l2 — l3803 — l58A1 cos05 = 0 71(80)" /1(80 y

4=0 /7(80). /7(80 y

This completes the proof.