Integrated Trajectory Planning and Sloshing Suppression for Three-Dimensional Motion of Liquid Container Transfer Robot Arm

Abstract

For liquid transfer system in three-dimensional space, the use of multijoint robot arm provides much flexibility. To realize quick point-to-point motion with minimal sloshing in such system, we propose an integrated framework of trajectory planning and sloshing suppression. The robot motion is decomposed into translational motion of the robot wrist and rotational motion of the robot hand to ensure the upright orientation of the liquid container. The trajectory planning for the translational motion is based on cubic spline optimization with free via points that produces smooth trajectory in joint space while it still allows obstacle avoidance in task space. Input shaping technique is applied in the task space to suppress the motion induced sloshing, which is modeled as spherical pendulum with moving support. It has been found through simulations and experiments that the proposed approach is effective in generating quick motion with low amount of sloshing.

Wisnu Aribowo 1, 2 and Takahito Yamashita 1 and Kazuhiko Terashima 1

Academic Editor:Rene V. Mayorga

1, Department of Mechanical Engineering, Toyohashi University of Technology, Toyohashi 441-8580, Japan
2, Department of Industrial Engineering, Bandung Institute of Technology, Bandung 40132, Indonesia

Received 17 November 2014; Accepted 14 April 2015; 18 May 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In this research, we address a case where a robot arm is used to move a liquid filled container in three-dimensional space. Some example application cases for this setting are transfer of molten metal in casting industries and liquid containers handling by service robots. The use of articulated robot arm with high degree of freedom for these purposes provides more dexterity and flexibility when compared with the more dedicated type of liquid transfer devices. There are two main performance characteristics to consider in this system: quick motion and minimal sloshing. Those two criteria are conflicting, as quick motion tends to generate much sloshing of the contained liquid. Sloshing may cause the liquid to spill out of the container, and in the case of molten metal too much sloshing may affect the physical properties of the liquid. Moreover, other task, for example, pouring, cannot be performed when residual sloshing still occurs, thereby potentially reducing the system productivity. Thus the challenge in this problem is to generate trajectories with concurrent consideration of motion time minimization and sloshing suppression for point to point motion in a three-dimensional working space where static obstacles may exist.

The motion speed of robot arm is mainly constrained by its kinematic and dynamic limit of the joints. One method to generate trajectories for such system is by posing the optimal control problem as a nonlinear optimization problem where the system states are discretized in some way. To reduce the large solution space and its associated computation time, the trajectory is often parameterized into a known structure, for example, polynomial function. Among the polynomial splines, cubic spline is the lowest order polynomial trajectory that guarantees continuation in acceleration and velocity as well as limited jerk. Higher order polynomials would result in smoother profile but longer motion time, while lower order ones would result in unlimited jerk.

Numerous cubic spline algorithms and optimization schemes have been proposed in the literature [1-5]. Most of them try to optimize either the motion time or the location of knots, depending on the application needs. Relatively few addressed the optimization of time and location of knots concurrently [5]. In the latter approach, the location of knots is not fixed, which is suitable for trajectory planning problem with obstacle avoidance. In addition, it often results in a quicker motion because there is less restriction for the optimization algorithm to reach optimal velocity profile.

The smooth cubic spline trajectories induce less vibration when compared to nonsmooth or lower order trajectories. However, in cases where minimal sloshing is required, smooth trajectory alone is not enough, and it needs to be supplemented with explicit sloshing control strategy. There have been numerous researches dealing with sloshing suppression in liquid container transfer system. Feedback control based approaches present good and relatively robust sloshing control ability [6, 7]. However, in some practical cases, sensor measurements cannot always be reliable or even are difficult to perform, for example, in the case of controlling sloshing of high-temperature molten metal. In these cases, feedforward approaches [8-14] are more useful. Based on sloshing models, the behavior of the system can be predicted, and thus the trajectory or motion path can be designed accordingly. Examples include generalized bang-bang control [9], infinite impulse response (IIR) filter [8], input shaping filter [10-13], and Hybrid Shape Approach filters [14]. Among them, the input shaping approach is particularly interesting because of its simplicity and effectiveness in practice. It works by decomposing the command input signal into two or more, using only the knowledge of natural frequency and damping ratio of the system [15, 16]. Over time, various type of input shapers with different characteristics have been proposed, and they have been successfully applied in practice for many kind of vibration problems.

With respect to the motion path, traditionally sloshing is analyzed for one-dimensional straight motion, where the sloshing is often modeled as a simple pendulum. The motion planning and sloshing suppression is relatively simple and straightforward for this kind of problem [9, 13]. Researches that address suppression of sloshing or pendulum sway in higher dimensional space are relatively few. Tzamtzi et al. [17] modeled the sloshing as simple pendulum in 2D planar vertical motion. Williams et al. [18] used a dynamic programming approach to solve the sway-free motion of suspended spherical pendulum in 2D planar horizontal motion. Yano and Terashima [14] proposed HSA approach to control sloshing in three degree of freedom Cartesian liquid container transfer system.

In this paper, we propose the integration of the trajectory planning of seven degrees of freedom liquid container transfer robot arm and the sloshing suppression of the contained liquid. We use the Mitsubishi PA 10-7C robot arm as the experiment device. The container is allowed to rotate about vertical axis while the pitch and roll angles are regulated to maintain upright orientation of the liquid container. Such system setup is different from that used in other related papers (e.g., [14, 17, 18]); in that translational motion in three-dimensional space plus orientation of the container needs to be controlled. The output of the proposed method is joint trajectories, to be used as a reference of the robot motion, thus obviating the need to change the closed-loop controller. Consequently, the proposed solution is easier to be applied to different kind of systems.

Because the robot arm has spherical wrist, the motion can be conveniently decoupled into translational motion of the wrist locus and rotational motion of robot hand. This decomposition is useful in our application case, because for sloshing suppression we want to take full control of the container orientation. For the trajectory generation of the wrist translational motion, we define the problem as a nonlinear optimization model where each joint trajectory is parameterized as cubic spline. The locations of knots are part of decision variables, thus not fixed, to allow obstacle avoidance in task space and quick motion. It is an extension to a more basic cubic spline optimization reported in our previous paper [5] in two points: motion in three-dimensional space is addressed by using multijoint robot arm, and obstacle avoidance constraints are added. As for the sloshing suppression, we use a new strategy of applying input shaping. The decoupled translational and rotational motions are each shaped by using its respective suitable input shaper and then combined as the final trajectory. We use the spherical pendulum model to simulate the response of the trajectory. In addition, a numerical sloshing simulation model is built for the purpose of examining the sloshing behavior in response to vertical motion of the liquid container in 3D space. Finally, experiments are carried out for a few example cases.

2. Sloshing Model and Suppression

2.1. Equivalent Sloshing Model

The theoretical natural frequency of sloshing inside an upright cylindrical container can be derived from Navier-Stokes equations as [19] [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the gravity constant, [figure omitted; refer to PDF] is the container radius, [figure omitted; refer to PDF] is the liquid height, and [figure omitted; refer to PDF] is the root of the derivative of Bessel function of the first kind where [figure omitted; refer to PDF] is the sloshing mode. For the fundamental first mode, [figure omitted; refer to PDF] equals 1.841.

For control analysis, the phenomenon of sloshing is often modeled by its simpler mechanical equivalents, for example, pendulum and mass-spring-damper system. Here, we model the sloshing as a simple pendulum, where the sway of the pendulum corresponds to the elevation of the liquid surface. A simple pendulum system, as depicted in Figure 1(a), consists of a point mass suspended by a massless rigid cable of length [figure omitted; refer to PDF] to a friction-less movable support. For the case where the support moves in straight line, it can be described by the following equation of motion: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the gravity constant and [figure omitted; refer to PDF] is the distance of the support from origin.

Figure 1: (a) Simple pendulum and (b) spherical pendulum.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

In cases where the support moves in two-dimensional horizontal space, the system can be modeled as a spherical pendulum. Figure 1(b) illustrates the spherical pendulum model, where the pendulum bob can freely swing in two degrees of freedom spherical space as response to movement of the support in horizontal planar space. Angle [figure omitted; refer to PDF] is the angle between plane ZY and the mass, and [figure omitted; refer to PDF] is the angle in the perpendicular direction.

The position of the pendulum bob in Cartesian coordinate system is [figure omitted; refer to PDF]

The equations of motion of the spherical pendulum can be derived as follows [18]: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the position of the support. Indeed, when there is no acceleration in Y direction, and initially [figure omitted; refer to PDF] , the above equations of motion will reduce to the simple pendulum case in (2).

2.2. Input Shaping

In the simple pendulum case, the natural frequency and damping ratio can be obtained by linearization of the equation of motion around equilibrium as [figure omitted; refer to PDF]

When excited by an impulse, the sway response of the pendulum is decaying sinusoidal: [figure omitted; refer to PDF]

Such sway can be suppressed by shaping the impulse input such that the sum of the response is zero. It is performed by splitting the original impulse into two or more impulses and delays them such that their responses cancel each other [15, 16]. For example, in the simplest case of splitting into two impulses, the amplitudes of the first and second impulses are [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively, where [figure omitted; refer to PDF] and time between impulse equals one-half damped period of vibration: [figure omitted; refer to PDF]

This set of impulses is called an input shaper, and the two-impulse input shaper above is known as zero vibration (ZV) shaper. A more robust input shaper type can be obtained by convolving two ZV input shapers together, resulting in three-impulse input shaper, known as zero vibration and derivative (ZVD) input shaper. Other than those basic input shapers, there are still many more input shapers developed with different characteristics and capabilities to tackle different type of vibration problems.

In a linear system, any signal can be represented as a collection of scaled and shifted impulses. Therefore the system response to a signal is the linear sum of the responses of each component impulse. It follows that the input shaping can be applied by convolving the system input signal with the desired input shaper, as shown in Figure 2. As a result, the motion time becomes longer, as much as the length of the input shaper. For example, applying ZV input shaper increases the motion time by half damped period of vibration, while with ZVD input shaper the motion time becomes one damped period longer. For general planar horizontal motion, the input signal may be in the form of acceleration or velocity reference of the motion.

Figure 2: Shaping arbitrary signal by convolution.

[figure omitted; refer to PDF]

For systems whose motion is commanded by rotational commands, for example, rotary crane and other rotating mechanical systems, regular input shapers do not work well for suppressing motion induced vibration. Those problems exhibit different characteristics in that the angular command impulse excites sway in radial as well as tangential directions. As the angular velocity increases, the natural frequency bifurcates from the nominal natural frequency, resulting in two sway modes in the system. For this kind of problem, Lawrence and Singhose [20] developed ZV2lin shaper. The input shaper consists of three impulses, where the amplitude and timing are determined by solving an optimization problem which seeks minimization of vibration at the end of impulse train. Those are determined not only by natural frequency and damping ratio, but also by the nominal angular velocity and rotation radius. The input shaper length is approximately equal to one damped period, more or less comparable to the length of ZVD shaper.

3. Trajectory Planning

The liquid transfer system concerned uses the Mitsubishi PA 10-7C robot arm. The frames placement and link dimension of the robot arm are shown in Figure 3. Table 1 lists the limit value of each joint angle and angular velocity. A cylindrical liquid container is attached to an extension at the robot tip. The diameter and height of the liquid container are 0.15 mm and 0.25 m, respectively. The distance between robot tip and the center of the liquid container is 0.12 meter; thus the length of the last link is practically 0.2 meter. For our current application, one of the arm joints, the S3 joint, is locked so that only six joints are actively used. Table 2 lists the Denavit-Hartenberg parameters of the robot arm configuration. The transformation matrix of a frame relative to the preceding frame is as follows: [figure omitted; refer to PDF]

Table 1: Limit values of joint velocities and angles.

Limit	Joint
S1	S2	S3	E1	E2	W1	W2
[figure omitted; refer to PDF] [rad/s]	1	1	2	2	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
[figure omitted; refer to PDF] [degree]	177	91	174	137	255	165	360

Table 2: Denavit-Hartenberg parameters.

Parameter	Joint
S1	S2	S3	E1	E2	W1	W2
[figure omitted; refer to PDF] [rad]	π /2	-π /2	π /2	-π /2	π /2	-π /2	0
[figure omitted; refer to PDF] [m]	0.315	0	0.45	0	0.5	0	0.2
[figure omitted; refer to PDF] [m]	0	0	0	0	0	0	0

Figure 3: (a) Frames placement and (b) dimension of the liquid container system and the base frame.

(a) Frames placement

[figure omitted; refer to PDF]

(b) Dimension of links and liquid container (in meter)

[figure omitted; refer to PDF]

The position and orientation of a link with respect to other link can be obtained by chain product of the transformation matrices: [figure omitted; refer to PDF]

The upper left [figure omitted; refer to PDF] rotation matrix [figure omitted; refer to PDF] describes the link orientation, while the vector [figure omitted; refer to PDF] represents the position.

The links E2, W1, and W2 form a spherical wrist, where the three rotational joint axes intersect at a common point. This robot configuration simplifies the inverse kinematic calculation because it can be decoupled into two simpler subsystems. The first three rotational links S1, S2, and E1 provide the positioning in three-dimensional space, while the spherical wrist adjusts the orientation of end effector.

In the inverse kinematics for position, the position of wrist locus in the task space (i.e. [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] ) is known, and the associated configuration of S1, S2, and E1 joints (i.e. [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] ) are to be sought. The values can be calculated geometrically as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the length of link [figure omitted; refer to PDF] , atan2 is the arc tangent function with two arguments, and [figure omitted; refer to PDF]

Meanwhile, the inverse kinematics for orientation calculates the configuration of E2, W1, and W2 joints (i.e. [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] ) when the desired orientation of end effector relative to the base frame [figure omitted; refer to PDF] is known. First, from the inverse kinematics for position, the rotation matrix describing the orientation of the third link can be obtained [figure omitted; refer to PDF] . Then, the rotation operation of the wrist joints required to realize the end effector posture is calculated as [figure omitted; refer to PDF]

The three Euler angles can be calculated from the matrix, and finally we can obtain the configuration of the wrist joints: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the element of the matrix [figure omitted; refer to PDF] at row [figure omitted; refer to PDF] and column [figure omitted; refer to PDF] . When there are more than one possibilities of angle values, choose the one suitable for the application.

Accordingly, the trajectory planning of the liquid container transfer system is decoupled into planning of the translation of the robot wrist position and rotation of the robot hand. The translation part is designed as joint space cubic spline trajectory, while the hand rotation is designed in task space to maintain upright orientation of the liquid container. Upright orientation means that the container has only one degree of freedom: the rotation along the vertical z -axis of fixed coordinate system. Those translation and rotation trajectories are then combined to make the final quick slosh-free trajectory solution. The trajectory planning steps are shown in Figure 4, which clearly depicts whether each step is done in joint space or task space. Transformation between joint and task space is carried out using direct and inverse kinematics.

Figure 4: Trajectory planning steps for (a) translational motion and (b) hand orientation.

(a) Trajectory planning for translation motion

[figure omitted; refer to PDF]

(b) Trajectory planning for hand orientation

[figure omitted; refer to PDF]

Figure 4(a) shows the overview of planning the translational motion. It starts with specifying the start point and end point, as well as several knots between them in task space, which are then transformed to joint space as the initial solution to the cubic spline optimization. The optimization step encompasses joint space as well as task space. Although the trajectory is generated in joint space, it is possible to define task space constraints in the optimization, for example, obstacle avoidance. The obtained solution is joint space trajectory, which is transformed back to task space, where the command shaping takes place. The solution of the translational motion is obtained as joint trajectories of links S1, S2, and E1.

Figure 4(b) shows the steps of planning the rotational motion of the robot hand. The input is the start and end angle of the container along z -axis relative to XZ vertical plane. The shaped rotation trajectory is then transformed to joint space as the solution of the rotation part of the motion in form of joint trajectories of links E2, W1, and W2. Although rotation trajectory generation is independent of translation trajectory generation, it has to be performed after because the former needs the orientation history of the translational motion as well as the total motion time in order to coordinate the overall motion.

3.1. Translational Motion Planning

In this section, we outline the cubic spline optimization approach for trajectory planning. It builds from the same basic approach as in a previous paper [5] for two-dimensional problem and then adapted to the current application of motion in three-dimensional space and is complemented with input shaping in the task space.

A cubic spline trajectory of a robot joint is comprised of several curve segments, as illustrated in Figure 5, where each segment is represented by a cubic polynomial function of time. The general form of the cubic function in a segment [figure omitted; refer to PDF] is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes the time, which ranges from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] .

Figure 5: Cubic spline trajectory.

[figure omitted; refer to PDF]

The curve segments are connected by knots. Thus for a trajectory that consists of [figure omitted; refer to PDF] curve segments, there are [figure omitted; refer to PDF] prespecified joint angle values [figure omitted; refer to PDF] , [figure omitted; refer to PDF] which consist of a start point, an end point, and [figure omitted; refer to PDF] knots between them. The trajectory is required to pass through all those points. Other requirements are continuity in velocity and acceleration at every knot. The relationship between joint angle [figure omitted; refer to PDF] , joint velocity [figure omitted; refer to PDF] , joint acceleration [figure omitted; refer to PDF] , and time interval [figure omitted; refer to PDF] is as follows: [figure omitted; refer to PDF]

In a joint trajectory with [figure omitted; refer to PDF] segments (or [figure omitted; refer to PDF] points), there are [figure omitted; refer to PDF] such equations. Together with the boundary requirements of acceleration ( [figure omitted; refer to PDF] ), we have a total of [figure omitted; refer to PDF] equations. The robot motion usually also requires that velocities at start and end points equal zero. To accommodate this, two more unknowns are added by inserting two virtual points: one after the start point and another one just before the end point. Furthermore, because [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are already known as zero, they can be removed from the equations. Joint angles [figure omitted; refer to PDF] are defined beforehand by the trajectory designer. The knots may be specified as points in Cartesian task space and later transformed into joint space by inverse kinematics. The above system of linear equations forms a symmetric tridiagonal system, which can be solved efficiently using the tridiagonal matrix algorithm (TDMA) to find acceleration at each knot.

The trajectories of S1, S2, and E1 joints are each designed as piecewise cubic splines. The joint trajectories are coupled to each other by using the same value of time interval variables [figure omitted; refer to PDF] . We formulate the trajectory generation problem as a nonlinear optimization problem. The objective function is combination of motion time and weighted squared acceleration. There are several constraints to consider in the trajectory generation: robot position in the task space, joint velocity, and joint jerk. According to (18), there are three kinds of variables involved: segment time, joint acceleration, and joint angle at knots. We let both the segment time and the joint angles (knots location) as the optimization variables. The optimization problem is thus as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the vector of optimization variables which contains the segment times as well as the joint angles excluding the fixed start point, end point, and the two virtual points: [figure omitted; refer to PDF]

Here, all internal knots are part of the optimization variables instead of predefined points. This provides flexibility for the optimization to generate quick motion as well as obstacle avoidance at the same time. If the planned motion is required to pass certain fixed points on the way, they can be added as additional equality constraints. The first term of the objective function is to minimize the total motion time, which is the main objective. The second term is used to avoid unnecessary motion of robot joints, where [figure omitted; refer to PDF] is the weight parameter. In our subsequent simulations and experiments, the parameter value equals [figure omitted; refer to PDF] .

The first and second constraints are joint angle and velocity constraints, according to the physical joint limits in Table 1. The third is jerk constraint, to prevent excessive value of joint jerk. The maximum jerk of each joint is set equal to 15 m/s³ . The last constraint is for obstacle avoidance, where the obstacle and robot links are modeled as rigid capsules with radius [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively. The function [figure omitted; refer to PDF] is the minimum distance between link [figure omitted; refer to PDF] and static obstacle [figure omitted; refer to PDF] in segment [figure omitted; refer to PDF] . The obstacle and links' axis are discretized into several representative points. The distance between a link and obstacle is defined as the minimum distance between the points in the obstacle and the points in the link.

We use the MATLAB Optimization Toolbox implementation of Sequential Quadratic Programming (SQP) for solving the above constrained nonlinear trajectory planning optimization problem. The SQP is one of the most used methods in solving general constrained nonlinear optimization problems. The method solves the problem by iteratively solving a quadratic subproblem, which is the quadratic approximation of the Lagrange function of the original nonlinear problem. The solution of the quadratic subproblem is then used to form a new set of optimization variables for the next iteration of the original problem. The calculation time differs depending on the computer used for calculation; it is around 15 seconds for a typical problem size in our Core2 Duo computer.

Figure 6 shows an example motion path of a cubic trajectory generated by the above method in perspective view and top view. The PA-10 robot arm is stationed at origin of the world frame. The robot has to move its wrist from start point [figure omitted; refer to PDF] to end point [figure omitted; refer to PDF] in the task space, which correspond to [figure omitted; refer to PDF] rad and [figure omitted; refer to PDF] rad in the joint space. The robot posture shown in the figure is the configuration at start point. A cylindrical obstacle with radius 0.05 m lies vertically in the workspace at [figure omitted; refer to PDF] . For this problem, four initial via points are specified randomly between the start and end points. The line going from the start point and the end point is the generated motion path in the task space, which successfully avoids the obstacle. The small spheres on the line denote the location of output via points. The total motion time is 2.7 seconds. Figure 7 shows the joint trajectories of the example case: the piecewise cubic joint angles and the piecewise quadratic joint velocities. The small circles in the figure denote the respective values at the knots, including at the two virtual knots. From the velocity figure, we can understand that the motion time is limited by joint S1, where the velocity hits its limit most of the time.

Figure 6: Motion path of example case of translation motion.

(a) Perspective view

[figure omitted; refer to PDF]

(b) Top view (X -Y plane)

[figure omitted; refer to PDF]

Figure 7: Trajectory of example case of translational motion: (a) joint angle (b) joint velocity.

(a) Joint angle

[figure omitted; refer to PDF]

(b) Joint velocity

[figure omitted; refer to PDF]

The next step in the translational motion planning is to suppress the sloshing generated by the trajectory obtained above. For analysis, we use the pendulum model in Section 2.1 to represent the sloshing, where the pendulum support is attached at the robot wrist locus and the pendulum sways freely as the robot moves. For the case where the pendulum support can move in a horizontal planar space, we first decompose the input signal to its projection into orthogonal Cartesian axes and then apply input shaping independently to each input signal. By doing this, the spherical pendulum is considered as two independent simple pendulums which act on orthogonal vertical planes such that the sway behavior in one subsystem is not affected by the states in the other subsystem. This simplified approach works well in small angle regions, and indeed linearization of (4) by small angle approximations and removing quadratic and higher terms results in two independent equations of simple pendulum.

We use the ZVD input shaping for sloshing suppression. For that, the velocity trajectory is decomposed into two orthogonal velocities parallel to x -axis and y -axis, respectively. Each velocity trajectory is then convolved individually with ZVD input shaper. The shaped results are then added back and integrated to produce the shaped trajectory in task space. Joint space solution can be obtained by inverse kinematics.

As an example, we simulate a spherical pendulum of length 0.2 meter moving according to the trajectory obtained above. The natural frequency and the vibration period of the pendulum model are 7.0 rad/sec and 0.8973 seconds, respectively. In sloshing problems, the damping is usually very small; therefore here we assume zero damping in the pendulum model. Figure 8 shows the original and ZVD shaped Cartesian velocities in X and Y directions. The motion time is 2.7 seconds and 3.6 seconds in the unshaped and ZVD shaped case, respectively. Figure 9 shows the sway angles in the spherical coordinate system. At the end of motion, sway is much suppressed in the shaped trajectory. Figure 10 compares the projection of the pendulum position at horizontal plane for the two cases: original and shaped trajectories. It is clear from the figure that the shaped trajectory generates significantly less sway than its unshaped counterpart.

Figure 8: Original and ZVD shaped velocity in task space: (a) X direction (b) Y direction.

(a) X direction

[figure omitted; refer to PDF]

(b) Y direction

[figure omitted; refer to PDF]

Figure 9: Pendulum sway angles: (a) without shaping (b) with ZVD shaping.

(a) Without shaping

[figure omitted; refer to PDF]

(b) ZVD input shaping

[figure omitted; refer to PDF]

Figure 10: Projection of pendulum bob position in XY plane: (a) without shaping (b) with ZVD shaping.

(a) Without shaping

[figure omitted; refer to PDF]

(b) ZVD input shaping

[figure omitted; refer to PDF]

3.2. Rotational Motion Planning

The liquid container has to be maintained upright all time. This means that there is only one task space DOF: the rotation along the vertical z -axis. Thus the hand rotational motion can be simply parameterized by the angle along z -axis relative to the XZ plane. Here, robot hand means any parts after the wrist, including the last robot link, the container holder, and the liquid container itself. The length of the robot hand is the distance between the wrist locus and the center of container, which equals 0.2 meter. The trajectory is generated by applying ZV2lin input shaping to a rectangular velocity profile, in which the maximum angular velocity is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the time of the unshaped translation motion and [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the desired start and finish angle value, respectively. This is to coordinate the translational and rotational motions, so that both motions end at approximately the same time.

As an example, we use the same spherical pendulum model to simulate the motion of the pendulum support on a circular trajectory of length [figure omitted; refer to PDF] radian, where the robot hand moves from angle [figure omitted; refer to PDF] radian to 0 (relative to XZ plane). The unshaped translation motion time is 2.7 seconds; thus the maximum angular velocity for the rotational motion equals 0.58 rad/sec. Upon applying ZV2lin input shaper, the rotational motion time becomes 3.47 seconds.

Figure 11 shows the comparison of the projections of the pendulum bob on XY plane of the unshaped trajectory and the trajectory shaped by ZV2lin input shaper. The circular motion starts from the origin and stops at (0.2, 0.2). In the shaped case, the pendulum sway is much suppressed, both the sway while the support is moving and the residual sway after the support stops. This shows the effectiveness of the shaper in suppressing vibration for motions generated by angular command.

Figure 11: Pendulum position in XY plane of circular motion for (a) without shaping (b) with shaping.

(a) Without shaping

[figure omitted; refer to PDF]

(b) ZV2lin shaped

[figure omitted; refer to PDF]

The trajectories of the three spherical wrist links can be obtained using the inverse kinematic relations in (13)-(16), the shaped angle value, and the orientation history of the translation motion part. The rotation matrix [figure omitted; refer to PDF] is the desired orientation of the liquid container: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the angle about z -axis.

The translational and rotational motions are independent as each is handled by different set of robot joints. Because of that, the final trajectory can be obtained simply by combining the joint trajectories of the translational and rotational motions. Once again, we simulate cases where spherical pendulum is attached to the robot tip, where the robot wrist locus moves from location (0.3, -0.8, 0.2) to (0.5, 0.5, 0.2) while at the same time the robot hand rotates with respect to the wrist from angle [figure omitted; refer to PDF] rad to 0 rad. In short, the trajectory is the combination of the simulated trajectories at Sections 3.1 and 3.2. Figure 12 compares the sway angles of the case without shaping and the shaped case. The shaping is done individually: the translation motion is shaped using ZVD input shaper and the rotational motion is shaped using ZV2lin input shaping. The motion time is the same as in Section 3.1: 2.7 seconds and 3.6 seconds for the unshaped and shaped case, respectively. Figure 13 shows the motion path of the robot tip and the projection of the pendulum bob in XY plane. From those two figures, we can understand that the strategy of combining individually shaped trajectories result in much reduced sway of the pendulum. In that way, for different conditions (in this case, translational and rotational motion), we can use suitable input shaping for each condition and then combined the generated trajectories.

Figure 12: Pendulum sway angles: (a) without shaping (b) combined input shaping.

(a) Without shaping

[figure omitted; refer to PDF]

(b) Combined input shaping

[figure omitted; refer to PDF]

Figure 13: Projection of pendulum bob position in XY plane: (a) without shaping (b) combined input shaping.

(a) Without shaping

[figure omitted; refer to PDF]

(b) Combined input shaping

[figure omitted; refer to PDF]

4. Numerical Simulation of Sloshing

In order to analyze the effect of vertical motion to sloshing, we built a numerical simulation of sloshing inside a moving liquid container based on exact distributed parameter computational fluid dynamics (CFD) model, not the simple pendulum model. The liquid container is cylindrical, with 150 mm diameter and 250 mm in height. The depth of liquid is 170 mm. The simulation model is developed in the open source OpenFOAM software package. The liquid oscillation is measured during simulation time with measurement sampling time 0.01 seconds by two virtual probes that record the liquid height near the container wall at the direction of x -axis and y- axis, respectively. By using the measured sloshing data, we identified the dominant fundamental mode of the sloshing at frequency 2.3 Hz.

As explained in Section 3.1, input shaping for translational motion is implemented by decomposing the trajectory into two Cartesian components (x -axis and y -axis) and then applying input shaping to those two trajectory components individually. The shaped trajectories are then combined as the final shaped trajectory. By using the spherical pendulum model, it has been shown that the approach works well in suppressing the pendulum sway. However, when we consider motion in 3D space, the motion component in z -axis also has to be taken into consideration. To examine the effect of acceleration in z -axis toward sloshing, we run a simulation comparing trajectories with and without vertical acceleration. In both trajectories, the liquid container moves horizontally for the first one second to generate some sloshing. Then, in trajectory A, we let the container stay still, while in trajectory B, the liquid container moves upward. Figure 14 shows the comparison of the sloshing generated by both trajectories. From second one, the sloshing frequency differs slightly, where it is a bit higher when the liquid container accelerates upward. Figure 15 shows the comparison in frequency domain. This simulation result agrees with the theory that states natural frequency of sloshing is directly proportional to square root of vertical acceleration (1). In the usual case of lateral movement, it consists of only the constant gravitational acceleration. But when it accelerates upward, the vertical acceleration would be higher than [figure omitted; refer to PDF] , and as a result natural frequency becomes higher. The theoretical natural frequency under constant vertical acceleration of 1 m/s² equals 2.4 Hz, which is approximately equal with the simulation result.

Figure 14: Comparison of the simulated sloshing under vertical acceleration.

[figure omitted; refer to PDF]

Figure 15: Comparison of the sloshing frequency under vertical acceleration.

[figure omitted; refer to PDF]

The consequence of this to the sloshing suppression is that we have to consider the frequency shift in the input shaper calculation. From simulation, we understand that the range of the frequency shift is not wide. With current setup of the robot arm, the vertical acceleration of a typical motion is under ±4 m/s² , which approximately corresponds to maximum frequency shift of ±0.4 Hz. To deal with this small shift, we use zero vibration and damping (ZVD) input shaper, which is a more robust version of the input shaper, instead of the basic ZV shaper.

To examine the sloshing in a vertical dominant motion trajectory, we setup another simulation case where liquid container is moved from location (0.2, 0.6, 0.2) to location (0.3, 0.2, 1.1). In this trajectory, vertical movement is relatively dominant. Thus, we expect the natural frequency to shift along the trajectory. Figure 16 shows the sloshing comparison. Consistent with the previous example case, the shaped trajectories result in less vibration. In addition, the robust ZVD shaped trajectory improves the sloshing suppression. The maximum residual sloshing is 46.22 mm, 4.04 mm, and 1.76 mm for original, ZV shaped, and ZVD shaped trajectories, respectively. From this simulation study, we understand that vertical motion affects the sloshing in that small shift of natural frequency occurs. The small shift justifies the use of the zero vibration and derivative (ZVD) input shaping, which is more robust than the default ZV input shaping, in our proposed solution.

Figure 16: Comparison of sloshing in vertically dominant case.

[figure omitted; refer to PDF]

5. Experiment Results and Discussion

In order to assess the effectiveness of the overall proposed system, including the trajectory planning for translational and rotational motion as well as the sloshing, we setup a few experiment cases. Table 3 lists the start and end position and orientation. The start points and end points are the position of the center of the container liquid. The angles are measured relative to XZ plane. Figure 17 illustrates the motion path of the three experiment cases, where the start and end points of each case are shown by small spheres.

Table 3: The start and end configuration of the experiment cases.

Case	Point	Angle
Start	End	Start	End
1	(0.3, -0.8, 0.2)	(0.8, 0.1, 0.3)	-π /2	0
2	(0.8, 0.1, 0.3)	(0.1, 0.9, 0.2)	0	π /2
3	(0.1, 0.9, 0.2)	(0.5, 0.6, 0.9)	π /2	0

Figure 17: Motion path of the experiment cases.

[figure omitted; refer to PDF]

The generated sloshing, as measured by level sensor is shown in Figure 18. It is measured for ten seconds from the start of motion. The motion time, residual vibration, and vibration reduction of the three cases are shown in Table 4. The motion time is the total time required to move from the start point to end point. The maximum residual sloshing for both unshaped and ZVD shaped trajectories is measured from the ZVD shaped motion time until the end of measurement, so as to make fair comparison.

Table 4: Result of experiments.

Case	Motion time [figure omitted; refer to PDF] s [figure omitted; refer to PDF]	Max residual sloshing [figure omitted; refer to PDF] mm [figure omitted; refer to PDF]	Sloshing reduction [figure omitted; refer to PDF] % [figure omitted; refer to PDF]
Unshaped	ZVD	Unshaped	ZVD
1	1.791	2.197	25.735	6.123	76.207
2	1.783	2.189	22.1	5.929	73.172
3	1.156	1.562	45.997	9.921	78.431

Figure 18: Sloshing measurement of (a) Case 1, (b) Case 2, and (c) Case 3.

(a) Case 1

[figure omitted; refer to PDF]

(b) Case 2

[figure omitted; refer to PDF]

(c) Case 3

[figure omitted; refer to PDF]

In general, we can see that the shaped trajectories generate much less sloshing at expense of longer motion times. The additional motion time is 0.406 seconds for all cases. That amount is equivalent to the length of the ZVD input shaper, which equals one vibration period. The sloshing reduction is more or less the same in all cases, where the largest is as much as 78.431% in Case 3. The motion path of Case 3 is mostly in upward direction. This shows that the solution does not have problem in handling trajectory that contains vertical motion.

6. Conclusion and Further Work

This paper has examined the case of liquid container transfer system using multi joint robot arm where quick motion and minimal sloshing are concurrently addressed. The proposed framework integrates trajectory planning in three-dimensional space and input shaping for sloshing suppression. The trajectory planning algorithm, based on cubic spline optimization, has been built to generate trajectories for translation motion in three-dimensional space with joint kinematic constraints and obstacle avoidance. The sloshing suppression is built into the algorithm, based on principle of input shaping and decoupling of the robot motion into translation and hand rotation subsystems. It has been shown using simulations and experiments that the generated trajectories induce much less sloshing, thus demonstrating the benefit of the suppression strategy. Further work will mainly be directed to improve the current framework by incorporating the dynamics consideration into the algorithm. Another planned improvement is to provide redundancy by utilizing one additional joint of the robot arm for better maneuverability.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] C.-S. Lin, P.-R. Chang, J. Y. S. Luh, "Formulation and optimization of cubic polynomial joint trajectories for industrial robots," IEEE Transactions on Automatic Control , vol. 28, no. 12, pp. 1066-1074, 1983.

[2] A. Piazzi, A. Visioli, "Global minimum-jerk trajectory planning of robot manipulators," IEEE Transactions on Industrial Electronics , vol. 47, no. 1, pp. 140-149, 2000.

[3] T. Chettibi, H. E. Lehtihet, M. Haddad, S. Hanchi, "Minimum cost trajectory planning for industrial robots," European Journal of Mechanics-A/Solids , vol. 23, no. 4, pp. 703-715, 2004.

[4] A. Gasparetto, V. Zanotto, "A technique for time-jerk optimal planning of robot trajectories," Robotics and Computer-Integrated Manufacturing , vol. 24, no. 3, pp. 415-426, 2008.

[5] W. Aribowo, K. Terashima, "Cubic spline trajectory planning and vibration suppression of semiconductor wafer transfer robot arm," International Journal of Automation Technology , vol. 8, no. 2, pp. 265-274, 2014.

[6] K. Yano, K. Terashima, "Robust liquid container transfer control for complete sloshing suppression," IEEE Transactions on Control Systems Technology , vol. 9, no. 3, pp. 483-493, 2001.

[7] M. Reyhanoglu, J. Rubio Hervas, "Nonlinear modeling and control of slosh in liquid container transfer via a PPR robot," Communications in Nonlinear Science and Numerical Simulation , vol. 18, no. 6, pp. 1481-1490, 2013.

[8] J. Feddema, C. Dohrmann, G. Parker, "A comparison of maneuver optimization and input shaping filters for robotically controlled slosh-free motion of an open container of liquid," in Proceedings of the American Control Conference, Albuquerque, NM, USA, 1997.

[9] L. Consolini, A. Costalunga, A. Piazzi, M. Vezzosi, "Minimum-time feedforward control of an open liquid container," in Proceedings of the 39th Annual Conference of the IEEE Industrial Electronics Society (IECON '13), pp. 3592-3597, IEEE, Vienna, Austria, November 2013.

[10] K. Terashima, K. Yano, "Sloshing analysis and suppression control of tilting-type automatic pouring machine," Control Engineering Practice , vol. 9, no. 6, pp. 607-620, 2001.

[11] M. Hamaguchi, T. Taniguchi, "Transfer control and curved path design for cylindrical liquid container," in Proceedings of 15th IFAC World Congress, pp. 1479-1479, Barcelona, Spain, 2002.

[12] A. Aboel-Hassan, M. Arafa, A. Nassef, "Design and optimization of input shapers for liquid slosh suppression," Journal of Sound and Vibration , vol. 320, no. 1-2, pp. 1-15, 2009.

[13] W. Aribowo, T. Yamashita, K. Terashima, H. Kitagawa, "Input shaping control to suppress sloshing on liquid container transfer using multi-joint robot arm," in Proceedings of the 23rd IEEE/RSJ 2010 International Conference on Intelligent Robots and Systems (IROS '10), pp. 3489-3494, Taipei, Taiwan, October 2010.

[14] K. Yano, K. Terashima, "Sloshing suppression control of liquid transfer systems considering a 3-D transfer path," IEEE/ASME Transactions on Mechatronics , vol. 10, no. 1, pp. 8-16, 2005.

[15] N. C. Singer, W. P. Seering, "Preshaping command inputs to reduce system vibration," Journal of Dynamic Systems, Measurement, and Control , vol. 112, no. 1, pp. 76-82, 1990.

[16] T. Singh, S. R. Vadali, "Robust time-delay control," Journal of Dynamic Systems, Measurement and Control , vol. 115, no. 2A, pp. 303-306, 1993.

[17] M. P. Tzamtzi, F. N. Koumboulis, N. D. Kouvakas, "A two stage robot control for liquid transfer," in Proceedings of the 12th IEEE International Conference on Emerging Technologies and Factory Automation, pp. 1324-1333, Patras, Greece, September 2007.

[18] C. Williams, G. Starr, J. Wood, R. Lumia, "Curvilinear transport of suspended payloads," in Proceedings of the IEEE International Conference on Robotics and Automation (ICRA '07), pp. 4537-4543, IEEE, Roma, Italy, April 2007.

[19] R. Ibrahim Liquid Sloshing Dynamics, Theory and Applications , Cambridge University Press, New York, NY, USA, 2005.

[20] J. Lawrence, W. Singhose, "Command shaping slewing motions for tower cranes," Journal of Vibration and Acoustics , vol. 132, no. 1, 2010.