全並聯式線性驅動平台機器人工作空間多餘自由度之最佳動力軌跡規劃

(1)

計畫類別：個別型計畫

計畫編號： NSC94-2212-E-011-007-

執行期間： 94 年 08 月 01 日至 95 年 07 月 31 日執行單位：國立臺灣科技大學機械工程系

計畫主持人：王勵群

計畫參與人員：溫家俊

報告類型：精簡報告

處理方式：本計畫可公開查詢

中華民國 95 年 10 月 31 日

(2)

UNCORRECTED

PROOF

1 2 Optimal dynamic trajectory planning

3 for linearly actuated platform type parallel manipulators

4 having task space redundant degree of freedom

5 Ka-Tjun Oen

¹

, Li-Chun T. Wang

^*

6

Department of Mechanical Engineering, National Taiwan University of Science and Technology, 43, Keelung Road, Taipei 106, Taiwan, ROC

7

Received 3 October 2005; received in revised form 29 May 2006; accepted 31 May 2006

8

9

Abstract

10

Two types of dynamic trajectory planning problems associated with task space redundant degree of freedom of parallel

11

manipulators are investigated. The ﬁrst problem involves synthesizing the point-to-point trajectory of the moving platform

12

such that the distribution of the cutting force of the tool bit along a speciﬁed machining contour is optimized. The second

13

problem is to determine the maximum constant cutting force along the contour while maintaining an optimal distribution

14

of the actuator forces. The numerical algorithms developed for solving these problems are based on the method of approx-

15

imate programming and take into account all of the dynamic and kinematic constraints imposed on the system.

16

2006 Published by Elsevier Ltd.

17

Keywords: Parallel manipulator; Dynamic trajectory planning; Task space; Redundant degree of freedom

18 19 1. Introduction

20 The positioning accuracy and structural rigidity of fully parallel, linearly actuated platform (LAP) type 21 manipulators are superior to conventional serial manipulators, and hence by attaching a high speed spindle 22 on the moving platform, they can be adopted as multi-axis machine tools to carry out heavy load, high 23 precision three-dimensional machining tasks. However, when using a six-degree-of-freedom (DOF) LAP to 24 perform a ﬁve-axis machining process such as contour milling and grinding, one needs only to specify the posi- 25 tion of the tool bit mounted on the spindle and the direction of the spindle axis along the given machining 26 contour, therefore the rotation of the moving platform about the spindle axis becomes a redundant DOF.

27 From a dynamic trajectory planning point of view, this redundant DOF is useful for adjusting the conﬁgura- 28 tion of the driving mechanism along the given trajectory to optimize certain dynamic performance indices of 29 the robot.

0094-114X/$ - see front matter 2006 Published by Elsevier Ltd.

doi:10.1016/j.mechmachtheory.2006.05.006

* Corresponding author. Tel.: +886 2 2737 6446; fax: +886 2 2737 6460.

E-mail address:[email protected](Li-Chun T. Wang).

1 Currently with the Northern Taiwan Institute of Science and Technology.

Mechanism and Machine Theory xxx (2006) xxx–xxx

www.elsevier.com/locate/mechmt

Mechanism and

Machine Theory

(3)

UNCORRECTED

PROOF

30 In this paper, two optimal dynamic trajectory planning problems associated with the task space redundant 31 DOF of fully parallel LAPs are investigated. The ﬁrst problem is to synthesize the optimal point-to-point rota- 32 tional trajectory of the moving platform about the spindle axis such that the cutting force is maximized along 33 the speciﬁed machining contour. The second problem involves synthesizing the optimal rotary trajectory of the 34 moving platform about the spindle axis to not only give the maximum constant cutting force but also optimize 35 the distribution of the actuator forces to avoid structural damage due to sudden jumps of the control forces.

36 Optimal dynamic trajectory planning problems for serial manipulators have been extensively studied by 37 many researchers. The problem of synthesizing time-optimal point-to-point robot motions was first investi- 38 gated by Kahn and Roth [1]. Later, Bobrow et al. [2] and Shin and McKey [3] independently developed 39 phase-plane switching point techniques to solve the minimum-time robot control problem with geometric path 40 and actuator constraints. This method has been subsequently improved by Shiller and Lu [4] and Shiller [5] to 41 handle dynamic singularities. A method for computing the time-optimal motions of non-redundant serial 42 manipulators in the presence of obstacles has also been presented by Shiller and Dubowsky [6]. Although 43 the basic concept behind these methods may be extended to parallel robots as well as systems of multiple coop- 44 erating manipulators without much difficulty [7,8], the application is nevertheless limited to minimum-time 45 trajectory planning problems and is not suitable for implementation in the forenamed maximum-force prob- 46 lems, for in which the profiles of the linear velocity and acceleration of the tool tip are also specified along the 47 given motion contour and hence the execution time can not be altered. In addition, no consideration has been 48 addressed to the issue of redundant DOF in these methods.

49 On the other hand, most of the works developed for trajectory planning problems of redundant serial 50 manipulators were concentrated on the use of the extra joint space DOF to synthesize feasible trajectories 51 for the purpose of, for example, Cartesian space obstacle avoidance [9], configuration singularity eliminations 52 [10], or to compensate for joint and torque limitations [11]. More recently, a genetic algorithm based approach 53 for synthesizing point-to-point motion trajectory to minimize vibration and execution time of flexible redun- 54 dant robots has also been reported [12]. However, these methods are also not suitable for handling the prob- 55 lems investigated in this paper. Not only because they are restricted to the simple open-loop kinematic 56 structure of serial robots, but also because the redundant DOF of the moving platform is specified in the task 57 (Cartesian) space rather than in the joint space.

58 The method presented in this paper for solving the two forenamed trajectory planning problems of the 59 LAPs takes into account the non-linear force–speed characteristics of the actuators and the kinematic con- 60 straints of the mechanism such as link interferences, actuator stroke limitations, and the passive joint bounds 61 as inequality constraints. In addition, the dynamic equations of the robot together with that of the tool spindle 62 are formulated as time varying equality constraints. Due to the highly non-linear nature of the constraints, the 63 optimal solutions to these problems can not be obtained analytically and numerical optimization techniques 64 must be adopted. The solution procedures developed here are based on the method of approximate program- 65 ming [13], which is frequently used to solve non-linear programming problems and has been successfully 66 implemented by Wang and Ravani for synthesizing the point-to-point load-optimal trajectory of serial manip- 67 ulators [14].

68 It should be pointed out that several methods for formulating the dynamic equations of parallel manipu-

69 lators have been developed in the past, based on the Newton–Euler formulation, the Lagrangian equations,

70 and the generalized D’Alembert’s principle [15]. However, the generalized coordinates used in these formula-

71 tions are either the Cartesian coordinates of the individual links or the joint variables which describe the rel-

72 ative motions between the links. Both require elaborate numerical coordinate partitioning algorithms to

73 separate the independent and the dependent generalized coordinates [16,17]. In addition, it is diﬃcult to rep-

74 resent the given task space trajectory of the tool bit in terms of the independent generalized coordinates used

75 in these formulations. The dynamic formulation presented in this paper is novel in the sense that it not only

76 takes the multiple closed kinematic loops of the robot into account, but also uses the six Cartesian coordinates

77 of the moving platform as the independent generalized coordinates of the system. With this approach, the

78 kinematically admissible virtual displacements and rotations of the individual links and the driving joints

79 can all be easily transformed into the task space by means of transformation Jacobian matrices, and hence

80 a concise parametric form of the dynamic equations in terms of the task space redundant DOF can be

81 obtained in a more straightforward manner.

(4)

UNCORRECTED

PROOF

82 2. Parametric form of the dynamic equations

83 The kinematic structure of fully parallel LAPs can be generally classified into two types based on the 84 arrangements of the linear actuators. As illustrated in Fig. 1(a) and (b), the linear actuators of the Stewart- 85 type LAP are attached to the limbs, and those of the Hexaglide-type [18] LAP are attached to the sliders along 86 the rails fixed on the base platform. Although the procedure for kinematic analysis of the two types of robots 87 is different [19], their dynamic equations can nevertheless be derived by using the same approach as described 88 below.

89 From the dynamics point of view, the LAP together with the spindle attached to the moving platform can 90 be considered as a system of interconnected rigid bodies, and the generalized principle of D’Alembert states 91 that for such systems the total virtual work done by the eﬀective forces along any kinematically admissible 92 virtual displacement and rotation is zero [17], or

93 dW

a

þ dW

e

¼ 0; ð1Þ

95 95

96 where dW

_a

is the total virtual work done by the gravity forces, inertia forces, and inertia moments, and dW

_e

is 97 that done by the externally applied forces and moments.

98 As shown in Fig. 2, the origin of the coordinate system U–V–W attached to the moving platform is denoted 99 as point P, which is located at the center of mass of the moving platform. The spindle is rigidly attached to the 100 moving platform and its axis is coincident with the W axis of the moving frame. The position vector of the 101 center of mass of the spindle is given by

P

_s

¼ P þ qw; ð2Þ

103 103

X Y Z V U

W

P

V U

Z Y O X

P

O

a

b

Fig. 1. Kinematic structure of fully parallel linearly actuated platform manipulators: (a) Stewart platform, (b) Hexaglide.

(5)

UNCORRECTED

PROOF

104 where q is the distance between the two mass centers and w is a unit vector which denotes the positive direction 105 of the W axis. Consequently, the virtual displacement of the center of mass of the spindle and its virtual rota- 106 tion can be, respectively, written as

107 dP

s

¼ dP þ dW qw ð3Þ

109 109 110 and 111

dW

s

¼ dW; ð4Þ

113 113

114 where dP is the virtual displacement of the mass center of the moving platform, and dW is a 3 · 1 vector rep- 115 resenting the virtual rotation of the moving platform [20]. Eqs. (3) and (4) can be combined together and writ- 116 ten into matrix notation as

d b X

s

¼ ½J

s

d b X; ð5Þ

118 118

119 where d b X

_s

¼ ½dP

^T_s

dW

^T_s

^T

and d b X ¼ ½dP

^T

dW

^T

, and the superscript T denotes the transpose of a vector or a 120 matrix; [J

s

] is a 6 · 6 transformation Jacobian matrix given by

; 122

122 123 in which [I

3

] and [0

3

], respectively, denote the 3 · 3 identity matrix and the null matrix, and ½~ w is a 3 · 3 skew 124 symmetric matrix in the form of [21]

½~ w ¼

0 w

z

w

_y

w

_z

0 w

x

w

y

w

_x

0 2

6 4

3 7 5;

126 126

127 where w

x

, w

y

, and w

z

are the components of w, expressed with respect to the ﬁxed X–Y–Z frame. Accordingly, 128 the total virtual work done by the inertia and gravity forces and the inertia moment of the spindle can be writ- 129 ten as

V U

W

ρ P

X Y Z

O Ps

S

P

φs

Fig. 2. Deﬁnition of coordinate systems.

(6)

UNCORRECTED

PROOF

130 dW

^s_a

¼ d b X

^T

½J

s

^T

bF

s

; ð6Þ

132 132

133 where b F

_s

is a 6 · 1 composite vector given by 134

bF

s

¼ m

_s

ða

s

gÞ

½I

s

a

s

þ x

s

½I

s

x

s

136 ð7Þ 136

137 in which g is the gravity acceleration vector, m

s

and [I

s

] are, respectively, the mass and the inertia matrix of the 138 spindle evaluated with respect to the ﬁxed coordinate system; a

s

, x

s

, and a

s

, respectively, denotes the linear 139 acceleration of the mass center, the angular velocity, and the angular acceleration of the spindle.

140 By using the same approach, it is not diﬃcult to show that the virtual works done by the inertia and gravity 141 forces and inertia moment of the moving platform and those of the links can also be written into a form sim- 142 ilar to that of Eq. (6) as

143 dW

^p_a

¼ d b X

^T

bF

p

ð8Þ

145 145 146 and 147

dW

^l_a

¼ d b X

^T

X

¹²

i¼1

½J

i

^T

bF

i

; ð9Þ

149 149

150 where b F

p

and b F

i

are 6 · 1 composite vectors containing the inertia and gravity forces and inertia moment of 151 the moving platform and that of each link, and are of the same form as Eq. (7), and [J

i

] (i = 1–12) are the 152 transformation Jacobian matrices of the links, the derivations of these matrices are given in Appendix A of 153 this paper.

154 The external forces applied to the system are the driving forces of the linear actuators and the cutting force 155 exerted on the tip of the tool bit mounted on the spindle. It is well known that for fully parallel LAPs, the 156 kinematically admissible virtual displacement of the driving joints is related to the virtual displacement and 157 rotation of the moving platform by [15]

dQ ¼ ½J

p

d b X; ð10Þ

159 159

160 where dQ is a 6 · 1 vector containing the virtual displacement of the joint variables, and [J

p

] is the 6 · 6 overall 161 Jacobian matrix of the robot. Therefore, the virtual work done by the driving forces can be written as 162

dW

^q_e

¼ d b X

^T

½J

p

^T

s; ð11Þ

164 164

165 where s is a 6 · 1 vector containing the driving forces provided by the linear actuators.

166 On the other hand, as shown in Fig. 3, the external force exerted on the tip of the tool bit during a machin- 167 ing process can be written as [22]

f ¼ f ðu

t

þ k

r

u

_n

Þ; ð12Þ

169 169

170 where f is the magnitude of the cutting force and u

t

is a unit vector along the tangential direction of the cutting 171 path; k

r

is a constant coeﬃcient, the value of which is dependent on various factors such as the tool shape, chip 172 thickness, feed rate, spindle speed, etc., and u

n

= w · u

t

is a unit vector which denotes the direction normal to 173 the cutting path. Consequently, the virtual work done by the cutting force can be written as

174 dW

^f_e

¼ d b X

^T

bF

f

; ð13Þ

176 176

177 where b F

f

is a 6 · 1 composite vector given by 178

bF

f

¼ f u

_t

þ k

r

u

_n

‘

_t

ðu

n

k

r

u

t

Þ þ r

t

w

; ð14Þ

180 180

181 in which ‘

t

is the distance from the center of mass of the moving platform to the tip of the tool bit, and r

t

is the 182 radius of the circular cross-section of the tool bit.

183 Substituting Eqs. (6), (8), (9), (11), and (13) into Eq. (1) yields

(7)

UNCORRECTED

PROOF

184 d b X

^T

½J

s

^T

bF

s

þ b F

_p

þ X

¹²

i¼1

½J

i

^T

bF

i

þ ½J

p

^T

s þ b F

_f

( )

¼ 0: ð15Þ

186 186

187 Since d b X contains only the kinematically admissible virtual displacement of the mass center and the virtual 188 rotation of the moving platform, which can be arbitrarily assigned for a six DOF LAP, therefore, Eq. (15) 189 can be satisﬁed if and only if

½J

s

^T

bF

s

þ b F

p

þ X

¹²

i¼1

½J

i

^T

bF

i

þ ½J

p

^T

s þ b F

f

¼ 0;

191 191 192 or 193

s ¼ ½J

¹_p

^T

bF

p

þ ½J

s

^T

bF

s

þ X

¹²

i¼1

½J

i

^T

bF

i

þ b F

_f

!

: ð16Þ

195 195

196 Eq. (16) is the basic form of the governor dynamic equation of the system, which is expressed in the task space 197 in terms of the Cartesian coordinates of the moving platform. If the system is used to perform a ﬁve-axis 198 machining task along a speciﬁed Cartesian contour, then the position and the linear velocity and acceleration 199 of the mass center of the moving platform are known, and can be expressed as functions of the execution time.

200 In addition, the orientation and the angular velocity and acceleration of the moving platform at any time in- 201 stance can be, respectively, written as

202 ½R

p

ð/Þ ¼ ½Rðw; /Þ½R

⁰_p

; ð17Þ

x

p

¼ x

⁰_p

þ _/w; ð18Þ

a

p

¼ a

⁰_p

þ x

⁰_p

_/w þ € /w; ð19Þ

204 204

205 where [R(w, /)] is a 3 · 3 spatial rotation matrix [20], in which / is the rotation angle of the moving platform 206 about the w axis, and ½R

⁰_p

¼ ½u j v j w is a 3 · 3 orthogonal matrix which represents the given orientation of 207 the moving platform. x

⁰_p

and a

⁰_p

are, respectively, the speciﬁed angular velocity and acceleration of the moving 208 platform, and _ / and € / are, respectively, the angular speed and angular acceleration of the moving platform 209 about the W axis.

U

V

n t

tw−r u

t

P W

f ut

Work piece rt

Feed

n

fk ur

f

Fig. 3. Vector representation of the cutting force.

(8)

UNCORRECTED

PROOF

210 Noting that the linear acceleration of the mass center of the tool bit can be written as (see Fig. 2) 211

a

s

¼ a

p

þ a

s

qw þ x

s

ðx

s

qwÞ; ð20Þ

213 213

214 where a

p

is the speciﬁed linear acceleration of the mass center of the moving platform, and 215

x

s

¼ x

p

þ _/

s

w; ð21Þ

a

s

¼ a

p

þ x

p

_/

s

w; ð22Þ

217 217

218 are the angular velocity and the angular acceleration of the tool bit, in which _ /

_s

is the constant rotation speed 219 of the spindle. Besides, according to the rotation axis theorem, the inertia matrix of the tool bit with respect to 220 the ﬁxed coordinate system at any instance is given by

221 ½I

s

¼ ½R

p

ð/Þ½

^s

I

_s

½R

p

ð/Þ

^T

; ð23Þ

223 223

224 where [

^s

I

_s

] is the constant inertia matrix of the tool bit evaluated with respect to the moving coordinate system.

225 Substituting Eqs. (17)–(19) into (20)–(23), and then substituting the resultant equations into (7) yields 226

bF

s

¼ ðA

^s₄

ð/Þ€ / þ A

^s₃

ð/Þ _/

²

þ A

^s₂

ð/Þ _/ þ A

^s₀

ð/ÞÞ; ð24Þ 228

228 229 where

A

^s₀

ð/Þ ¼ m

_s

ða

p

þ a

⁰_p

qw þ x

⁰_p

ðx

⁰_p

qwÞ gÞ

½I

s

ða

⁰_p

þ _/

s

x

⁰_p

wÞ þ ðx

⁰_p

þ _/

s

wÞ ½I

s

ðx

⁰_p

þ _/

s

wÞ

" #

;

A

^s₂

ð/Þ ¼ 0

½I

s

ðx

⁰_p

wÞ þ x

⁰_p

½I

s

w þ w ½I

s

x

⁰_p

" #

;

A

^s₃

ð/Þ ¼ 0 w ½I

s

w

231 231 232 and

A

^s₄

ð/Þ ¼ 0

½I

s

w

234 : 234

235 By using the same approach, it can be shown that the inertia and gravity forces of the moving platform and 236 that of the links, b F

_p

and b F

_i

(i = 1–12), can also be derived into the same form as Eq. (24) [23]. Consequently, 237 by substituting Eqs. (14) and (24), and the corresponding expressions of b F

_p

and b F

_i

into Eq. (16), after rear- 238 ranging and combining terms, it can be transformed into a parametric form as

239 s ¼ A

4

ð/Þ€ / þ A

3

ð/Þ _/

²

þ A

2

ð/Þ _/ þ A

1

ð/Þf þ A

0

ð/Þ; ð25Þ 241

241 242 where the coeﬃcients A

i

(/) (i = 0–4) are 6 · 1 vectors and their elements are non-linear functions of / only.

243 3. Formulation of the trajectory planning program

244 Noting that Eq. (25) represents a system of second-order non-linear diﬀerential equations governing the 245 dynamic relationships between the control force s, the cutting force f, and the rotation angle /. As for a 246 ﬁve-axis machining task the rotation of the moving platform about the spindle axis becomes a redundant 247 DOF, therefore / can be considered as a free parameter to be used for optimizing the dynamic performance 248 of the system. In this regard, two optimal dynamic trajectory planning problems are studied in this paper, as 249 stated below:

250 P1: Synthesize the optimal trajectory of / and the corresponding control forces in-between the two end 251 points of the speciﬁed machining contour, such that the cutting force is maximized along the contour.

252 P2: Find the optimal trajectory of / and the maximum allowable constant cutting force along the machining

253 contour while keeping the control forces as close to the mean value of their upper and lower bounds as

254 possible.

(9)

UNCORRECTED

PROOF

255 256 In addition to the governing dynamic equations, it should be noted that both of the forenamed opti- 257 mization problems are also subjected to the limitations of the driving force and the rotational angle given 258 by

259 s

^l

ðq; _qÞ 6 s 6 s

^u

ðq; _qÞ ð26Þ

261 261 262 and 263

/

^l

6 / 6 /

^u

; ð27Þ

265 265

266 where s

^u

ðq; _qÞ and s

^l

ðq; _qÞ are the upper and lower bounds of the driving force, which are in general non-linear 267 functions speciﬁed by the force–speed characteristics of the linear actuators. /

^u

and /

^l

are, respectively, the 268 upper and lower bounds of the rotation angle, which are dependent on the kinematic conditions of link inter- 269 ferences, passive joint constraints and the stroke limitations of the linear actuators [24,25].

270 4. Linearization of the parametric dynamic equation

271 The solution procedures developed in this paper for solving the two optimal control problems formulated 272 in the last section are based on the method of approximate programming (MAP) [13]. It involves first con- 273 structing an initial discretized trajectory of the rotation angle and then linearizing the dynamic equation about 274 each of the discrete points, such that the original non-linear programming problem can be replaced by a sim- 275 plified linear or quadratic programming problem. The discretized trajectory is then updated iteratively based 276 on the solution to these simplified problems until the approximation error of the linearization converges to a 277 specified degree of precision.

278 By using the ﬁrst-order approximation, the discretized form of the motion trajectory of the rotation angle / 279 at the k + 1th iteration can be written as

280 /

^kþ1_j

¼ /

^kþ1_j1

þ h _/

^kþ1_j1

; for j ¼ 2 to m þ 1 ð28Þ

/ _

^kþ1_j

¼ _/

^kþ1_j1

þ h€ /

^kþ1_j1

; for j ¼ 2 to m ð29Þ

282 282

283 where m is the total number of the discretized time segments, and h denotes the time interval of each segment.

284 By repeatedly substituting the diﬀerence equations of the succeeding segments into the preceding ones and not- 285 ing that /

₁

¼ _/

1

¼ 0, Eqs. (28) and (29) can be written into the expanded form as

286 /

^kþ1_j

¼ h

²

X

^j2

i¼1

ðj i 1Þ€ /

^kþ1_i

; ð30Þ

/ _

^kþ1_j

¼ h X

^j1

i¼1

/ €

^kþ1_i

: ð31Þ

288 288

289 At each set point of the discretized trajectory, expanding Eq. (25) into Taylor series with respect to the points 290 of the previous (kth) trajectory and neglecting higher order terms yields

291 s

^kþ1_j

¼ B

j

/ €

^kþ1_j

þ C

j

/ _

^kþ1_j

þ D

j

/

^kþ1_j

þ E

j

f

_j^kþ1

þ K

j

ð32Þ 293

293 294 for j = 1 to m + 1, where B

_j

¼ os

o€ / ¼ A

4

ð/

^k_j

Þ;

C

_j

¼ os

o _/ ¼ 2A

3

ð/

^k_j

Þ _/

^k_j

þ A

2

ð/

^k_j

Þ;

D

_j

¼ os o/ ; E

_j

¼ os

of ¼ A

1

ð/

^k_j

Þ 296

296

(10)

UNCORRECTED

PROOF

297 and

K

_j

¼ s

^k_j

B

j

/ €

^k_j

C

j

/ _

^k_j

D

j

/

^k_j

E

j

f

_j^k

: 299

299 300 Substituting Eqs. (30) and (31) into (32) for the ﬁrst m set points gives 301

s

^kþ1_j

¼ B

j

/ €

^kþ1_j

þ hC

j

/ €

^kþ1_j1

þ X

^j2

i¼1

a

j;i

/ €

^kþ1_i

þ E

j

f

_j^kþ1

þ K

j

ð33Þ

303 303

304 for j = 1 to m, where a

j,i

= hC

j

+ (j i 1)h

²

D

j

. In addition, by noting that both the angular speed _ /

_mþ1

and 305 angular acceleration € /

_mþ1

should be zero at the last set point, it can be shown that

306 s

^kþ1_mþ1

¼ X

^m1

i¼1

a

_mþ1;i

/ €

^kþ1_i

þ E

_mþ1

f

_mþ1^kþ1

þ K

_mþ1

; ð34Þ

308 308

309 where a

m+1,i

= h

²

(m i)D

m+1

and K

_mþ1

¼ s

^k_mþ1

D

_mþ1

/

^k_mþ1

E

_mþ1

f

_mþ1^k

. Eqs. (33) and (34) are the basic form 310 of the k + 1th linearized dynamic equations of the system about the kth discretized trajectory, and they can be 311 combined together and written into matrix form as

312 T ¼ ½M

1

M

₂

X

₁

X

₂

þ T

^h

; ð35Þ

314 314

315 where T ¼ ½s

^kþ1₁ ^T

s

^kþ1₂ ^T

s

^kþ1_mþ1^T

^T

is a 6(m + 1) · 1 vector containing the control forces of the actuators of the 316 k + 1th iteration and T

^h

¼ ½K

^T₁

K

^T₂

K

^T_mþ1

^T

is a 6(m + 1) · 1 constant vector; X

1

¼ ½€ /

^kþ1₁

/ €

^kþ1₂

€ /

^kþ1_m

^T

and 317 X

2

¼ ½f

₁^kþ1

f

₂^kþ1

f

_mþ1^kþ1

^T

are, respectively, m · 1 and (m + 1) · 1 vectors containing the angular acceleration of 318 the moving platform and the cutting forces. The sub matrices M

1

and M

2

are, respectively, 6(m + 1) · m and 319 6(m + 1) · (m + 1) matrices given by

M

₁

¼

B

₁

0 0 0

hC

2

B

2

0 0

.. . .. .

.. . a

m;1

a

m;2

a

m;3

B

m

a

_mþ1;1

a

_mþ1;2

a

_mþ1;3

0 2

6 6 6 6 6 6 6 4

3 7 7 7 7 7 7 7 5

and M

₂

¼

E

₁

0 0 0 0 E

2

0 0 0 0 E

₃

0 .. . .. .

.. . .. .

.. . 0 0 0 E

_mþ1

2 6 6 6 6 6 6 6 4

3 7 7 7 7 7 7 7 5 :

321 321

322 It should also be noted that at each set point of the discretized trajectory, the non-linear constraints 323 imposed on the actuator forces (Eq. (26)) could be written as

324 s

^l_j

ðq

^kþ1_j

; _q

^kþ1_j

Þ 6 s

^kþ1_j

6 s

^u_j

ðq

^kþ1_j

; _q

^kþ1_j

Þ ð36Þ 326

326 327 for j = 1 to m + 1, in which q

^kþ1_j

and _q

^kþ1_j

are the displacement and velocity vectors of the driving joints eval- 328 uated at the jth set point of the k + 1th trajectory, and hence are unknowns. Although it is possible to trans- 329 form these variables into functions of /

^kþ1_j

and _ /

^kþ1_j

based on the inverse kinematics transformations, this 330 approach would considerably increase the problem complexity. By noting that the diﬀerence between two suc- 331 cessive discretized trajectories would be asymptotically eliminated when the procedure converges, a simpler 332 way to handle these constraints is to approximate Eq. (36) by

333 s

^l_j

ðq

^k_j

; _q

^k_j

Þ 6 s

^kþ1_j

6 s

^u_j

ðq

^k_j

; _q

^k_j

Þ ð37Þ

335 335

336 for j = 1 to m + 1, where q

^k_j

and _q

^k_j

are the joint displacement and velocity of the jth point at the kth iteration, 337 and which can be evaluated by solving the inverse kinematics problems once the kth discretized trajectory is 338 established. Consequently, the bounds on the control forces can be treated as constant bounds for the k + 1th 339 trajectory.

340 By substituting Eq. (35) into Eq. (37) and combining terms, one obtains 341

T

^l

6 ½M

1

M

2

X

1

X

₂

6 T

^u

; ð38Þ

343

(11)

UNCORRECTED

PROOF

344 where

T

^l

¼

s

^l₁

ðq

^k₁

; _q

^k₁

Þ s

^l₂

ðq

^k₂

; _q

^k₂

Þ

.. . s

^l_mþ1

ðq

^k_mþ1

; _q

^k_mþ1

Þ 2

6 6 6 6 6 4

3 7 7 7 7 7 5

T

^h

and T

^u

¼

s

^u₁

ðq

^k₁

; _q

^k₁

Þ s

^u₂

ðq

^k₂

; _q

^k₂

Þ

.. . s

^u_mþ1

ðq

^k_mþ1

; _q

^k_mþ1

Þ 2

6 6 6 6 6 4

3 7 7 7 7 7 5

T

^h

: 346

346 347 On the other hand, the side constraints imposed on the rotation angle (Eq. (27)) at each set point of the 348 discretized trajectory can be written as

349 /

^l_j

6 /

_j

6 /

^u_j

ð39Þ

351 351

352 for j = 3 to m + 1, where /

^u_j

and /

^l_j

are, respectively, the upper and lower bounds of the rotation angle at the 353 jth set point, and which can be explicitly evaluated in advance based on the kinematic constraints [23]. Substi- 354 tuting Eq. (30) into Eq. (39) and written into matrix form gives

355 U

^l

6 ½M

3

X

1

6 U

^u

; ð40Þ

357 357

358 where U

^l

¼ 1

h

²

/

^l₃

/

^l₄

/

^l_mþ1

^T

; U

^u

¼ 1

h

²

/

^u₃

/

^u₄

/

^u_mþ1

T

360 360 361 and

½M

3

¼

1 0 0 0 0

2 1 0 0 0

3 2 1 0 0

.. .

.. . .. .

.. . .. . m 1 m 2 m 3 1 0 2

6 6 6 6 6 6 4

3 7 7 7 7 7 7 5 :

363 363

364 5. Computational algorithm for problem P1

365 Based on the trajectory discretization and noting that the magnitude of the cutting force should always be a 366 non-negative value at each of the set points, the objective function of problem P1 presented in Section 3 for 367 each iteration of the MAP can be formulated as to maximize

J

_P1

¼ X

^mþ1

j¼1

f

j

; ð41Þ

369 369

370 subjected to the inequality constraints (38) and (40), and f

j

P 0 for j = 1 to m + 1. As the objective function 371 and all the constraint equations are linear, this problem can be solved by using standard linear programming 372 algorithms such as the simplex method [26]. Consequently, the solution to the original non-linear optimization 373 problem can be obtained by using the following iterative algorithm:

374 Algorithm 1.

375 Step 1. Initialize the trajectory of / as a fourth-order polynomial by using the method presented in Appendix 376 B. Select a moderate number m and discretize the trajectory into m equally spaced segments, and then 377 evaluate the bounds U

^l

and U

^u

. Specify a convergence criterion e 1 and initialize the iteration coun- 378 ter k = 1, and then go to the next step.

379 Step 2. Compute the coefﬁcients matrices [M

1

], [M

2

], [M

3

] and vectors T

^h

, T

^l

, T

^u

of the constraint equations

380 (38) and (40).

(12)

UNCORRECTED

PROOF

381 Step 3. Solve the linear programming problem to evaluate the state variables X

1

and X

2

. Update the discret- 382 ized trajectory by using Eqs. (28) and (29), and then compute the control force vector T from Eq. (25).

383 Step 4. If maxfj€ /

^k_j

€ /

^kþ1_j

j; j ¼ 1 mg 6 e then the process is converged, stop the iteration and then output 384 the results; otherwise set k = k + 1 and then go to step 2.

385 386

387 6. Computational algorithm for problem P2

388 Based on the concept of MAP, the multi-objective optimization problem P2 formulated in Section 3 can be 389 solved by using an iterative algorithm which consists of two loops. The inner loop involves of ﬁnding the 390 optimal distribution of the control forces for any given value of constant cutting force, and the outer loop 391 performs a one-dimensional bisector search to ﬁnd the maximum allowable value of the constant cutting 392 force.

393 Noting that by specifying a constant value to the cutting force f, the linearized dynamic equation (Eq. (32)) 394 would be reduced to

s

^kþ1_j

¼ B

j

/ €

^kþ1_j

þ C

j

/ _

^kþ1_j

þ D

j

/

^kþ1_j

þ K

j

ð42Þ

396 396

397 for j = 1 to m + 1, and hence Eqs. (35) and (38) can be, respectively, simpliﬁed to 398

T ¼ ½M

1

X

1

þ T

^h

ð43Þ

400 400 401 and 402

T

^l

6 ½M

1

X

1

6 T

^u

: ð44Þ

404 404

405 Consequently, the optimization problem for each iteration step of the inner loop can be formulated as to 406 minimize

407 J

_P2

¼ ðT UÞ

^T

ðT UÞ; ð45Þ

409 409

410 where U ¼ 1=2½ðs

^l₁

þ s

^u₁

Þ ðs

^l_mþ1

þ s

^u_mþ1

Þ

^T

, subject to the equality constraint equation (43) and the inequality 411 constraint equations (40) and (44). By substituting Eq. (43) into (45), the objective function can be expanded 412 into the quadratic form as

413 J

_P2

¼ X

^T₁

½HX

1

þ 2G

^T

X

₁

þ C; ð46Þ

415 415

416 where [H] = [M

1

]

^T

[M

1

], G = [M

1

]

^T

(T

^h

U), and C = (T

^h

U)

^T

(T

^h

U). Since the state variable vector X

1

is 417 subjected to linear inequality constraint given by Eq. (44) only, the problem can be solved by using standard 418 quadratic programming algorithms [27]. The complete procedure for solving problem P2 is summarized 419 below:

420 Algorithm 2.

421 Step 1. Initialize the trajectory of / and evaluate the bounds U

^l

and U

^u

as in step 1 of Algorithm 1. Set the 422 initial value of the constant cutting force to f = 0, and specify an increment Df and the convergence 423 criterion e. Set both the iteration counter of the inner loop k and that of the outer loop s to 1, and then 424 go to the next step.

425 Step 2. Compute the coefﬁcients matrices [M

1

], [H], [M

3

] and vectors, G, U, T

^h

, T

^l

of Eqs. (40), (44) and (46).

426 Step 3. Solve the quadratic programming problem to evaluate the state variable vector X

1

. If no feasible solu- 427 tion exists, then the magnitude of the cutting force must be reduced, go to Step 7; otherwise go to the

428 next step.

429 Step 4. Update the discretized trajectory by using Eqs. (28) and (29), and then compute the control force vec- 430 tor T from Eq. (25).

431 Step 5. If maxfj€ /

^k_j

€ /

^kþ1_j

j; j ¼ 1 mg 6 e, then the inner loop is converged, set k = 1 and save the current

432 value of the constant cutting force as f

_s

, and then go to the next step; otherwise set k = k + 1 and then

433 go to step 2.

(13)

UNCORRECTED

PROOF

434 Step 6. If s > 1 and jf

s

f

s1

j 6 e, then the outer loop is converged, terminate the iteration, and then output 435 the result; otherwise if s = 1 then set f = f + Df, else set s = s + 1 and f = f + (1/2)

^s

Df and then go to

436 step 2.

437 Step 7. Set s = s + 1 and f = f (1/2)

^s

Df, and then go to step 2.

438

Table 1

Dimension parameters of the Stewart platform

Limb number i Bi(cm) ^UVWri(cm)

X Y Z U V W

1 1.743 19.923 0.000 8.191 5.735 0.000

2 1.743 19.923 0.000 8.191 5.735 0.000

3 18.126 8.452 0.000 9.063 4.226 0.000

4 16.383 11.471 0.000 0.871 9.961 0.000

5 16.383 11.471 0.000 0.871 9.961 0.000

6 18.126 8.452 0.000 9.063 4.226 0.000

Table 2

Dimension parameters of the Hexaglide

Limb number i B⁰_i ðcmÞ ^UVWri(cm) ‘i(cm) q^d_i ðcmÞ

X Y Z U V W

1 0.000 15.000 0.000 6.428 7.660 0.000 60.0 60.605

2 0.000 15.000 0.000 9.397 3.420 0.000 60.0 60.797

3 0.000 0.000 0.000 9.848 1.736 0.000 40.0 36.915

4 0.000 15.000 0.000 1.736 9.848 0.000 60.0 49.498

5 0.000 15.000 0.000 3.420 9.397 0.000 60.0 58.004

6 0.000 0.000 0.000 7.660 6.428 0.000 40.0 40.267

Y X

Z

W

V U

k

(0,0,15cm) 10 cm

O

P

) φ(t )

θ(t

P

Fig. 4. Deﬁnition of the desired machining contour.

(14)

UNCORRECTED

PROOF

439 7. Numerical examples

440 The LAPs used for the examples are the Stewart platform and the Hexaglide as shown in Fig. 1. The geo- 441 metric dimensions of these robots are given in Tables 1 and 2. The speciﬁed machining contour is a circle with 442 10 cm radius parallel to the x–y plane of the ﬁxed frame and centered at (0, 0, 15 cm), as shown in Fig. 4. The

0 2 4 6 8 10

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Time (Sec.) Rolling Path

Bounds φ^(rad)

Fig. 5. Optimal trajectory of /—Stewart platform (problem P1).

0 2 4 6 8 10

-1.5 -1 -0.5 0 0.5 1 1.5

Bounds (rad)

φ

Fig. 6. Optimal trajectory of /—Hexaglide (problem P1).

(15)

UNCORRECTED

PROOF

443 total execution time of the trajectory is speciﬁed to 10 s, and the orientation matrix and the angular velocity 444 and acceleration of the moving platform are, respectively, speciﬁed by

½R

⁰_p

¼ ½Rðk; hðtÞÞ

0:832 0 0:555

0:0 1 0:0

0:555 0 0:832 2

6 4

3 7 5;

x

⁰_p

ðtÞ ¼ _hðtÞk 446

446 447 and

0 2 4 6 8 10

-250 -200 -150 -100 -50 0 50 100 150 200 250

Time (Sec.)

Driving Force (N.)

Upper Bound

Lower Bound

τ

1

τ

2

τ

3

0 2 4 6 8 10

-250 -200 -150 -100 -50 0 50 100 150 200 250

Time (Sec.)

Driving Force (N.)

Upper Bound

Lower Bound

τ

4

τ

₅

τ

⁶

Fig. 7. Optimized control force proﬁles—Stewart platform (problem P1).

(16)

UNCORRECTED

PROOF

a

⁰_p

ðtÞ ¼ € hðtÞk;

449 449

450 where

hðtÞ ¼

pð1:2t

³

1:6t

⁴

þ 0:6t

⁵

Þ; 0 6 t < 1 s;

0:2pt; 1 s 6 t 6 9 s;

2p pð1:2ð10 tÞ

³

1:6ð10 tÞ

⁴

þ 0:6ð10 tÞ

⁵

Þ; 9 s < t 6 10 s 8 >

<

> : 452

452 453 and k = [0 0 1]

^T

. As shown in the ﬁgure, the direction of the tool bit is coincident with the W axis of the 454 moving frame and the spindle has a constant angular speed of 5 · 10

⁴

rpm. The torque-speed characteristic 455 equation of the linear actuators of the Stewart platform are given by 200 5 _q

i

6 s

i

6 200 5 _q

i

(for

0 2 4 6 8 10

-500 -400 -300 -200 -100 0 100 200 300 400 500

Time (Sec.)

Driving Force (N.)

Upper Bound

Lower Bound

τ

1

τ

3

τ

2

0 2 4 6 8 10

-500 -400 -300 -200 -100 0 100 200 300 400 500

Time (Sec.)

Driving Force (N.)

Lower Bound Upper Bound

τ

4

τ

5

τ

6

Fig. 8. Optimized control force proﬁles—Hexaglide (problem P1).

(17)

UNCORRECTED

PROOF

456 i = 1–6), and that of the Hexaglide are 400 10 _q

i

6 s

_i

6 400 10 _q

i

(for i = 1–6). The convergence criterion 457 for the algorithms is speciﬁed as e = 10

⁴

.

458 For problem P1, Algorithm 1 converged after 6 and 10 iterations, respectively, for the Stewart Platform and 459 the Hexaglide. The computational results are shown in Figs. 5–10. Figs. 5 and 6 are the optimal trajectory of / 460 of the two robots, and Figs. 7 and 8 are the plots of the corresponding proﬁles of the control forces. The initial 461 and the optimized distributions of the cutting force with respect to the speciﬁed machining contour are plotted 462 in Figs. 9 and 10. It can be seen that the optimization algorithm successfully increased the minimum value of 463 the cutting force of the Stewart platform from 153 N to 165 N and that of the Hexaglide from 30 N to 48 N.

0 2 4 6 8 10

150 160 170 180 190 200 210 220 230 240

Time (Sec.)

Cutting Force (N.)

Optimum

Initial

*

165*

153

Fig. 9. Comparison of the initial and the optimized cutting force distributions of the Stewart platform.

0 2 4 6 8 10

0 50 100 150 200 250 300 350 400 450

Time (Sec.)

Cutting Force (N.)

Optimum

Initial *

*30 48

Fig. 10. Comparison of the initial and the optimized cutting force distributions of the Hexaglide.

(18)

UNCORRECTED

PROOF

464 For problem P2, the maximum constant cutting force returned from Algorithm 2 for the Stewart platform 465 is 250 N and that for the Hexaglide is 175 N. The corresponding optimal trajectories of / are, respectively, 466 shown in Figs. 11 and 12. The optimal distributions of the control forces for the two robots are plotted in

0 2 4 6 8 10

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

Bounds (rad)

φ

Fig. 11. Optimal trajectory of /—Stewart platform (problem P2).

0 2 4 6 8 10

-1.5 -1 -0.5 0 0.5 1 1.5

Time (Sec.)

Rolling Path Bounds (rad)

φ

Fig. 12. Optimal trajectory of /—Hexaglide (problem P2).

(19)

UNCORRECTED

PROOF

467 Figs. 13 and 14, by comparing to Figs. 7 and 8, it can be seen that the bang–bang eﬀect has been eﬀectively 468 eliminated.

469 Since the problems dealt with in this paper have not been investigated previously and currently no hardware 470 facilities are available in the authors’ institution, therefore it is difficult to obtain experimental data to compare 471 with the computational results. Nevertheless, the accuracy of the presented dynamic model and the optimiza- 472 tion algorithms can be justified based on the principle of work and energy. For the closed circular trajectory 473 specified in the examples, the total change of the potential energy and that of the kinetic energy of the system is 474 zero, as the initial and the final states are the same, therefore the sum of the work done by the actuator forces 475 and that done by the cutting force should be zero. Fig. 15 shows the total work evaluated with respect to the 476 optimal constant cutting force trajectories of the two robots synthesized on the basis of different number of 477 time segments. It can be seen that the inaccuracy of the algorithm is primarily due to the discretization error

τ

3

τ

2

τ

1

0 2 4 6 8 10

-250 -200 -150 -100 -50 0 50 100 150 200 250

Time (Sec.)

0 2 4 6 8 10

Time (Sec.)

Lower Bound

Lower Bound Upper Bound Upper Bound

Driving Force (N.)

-250 -200 -150 -100 -50 0 50 100 150 200 250

Driving Force (N.)

τ

4

τ

5

τ

6

Fig. 13. Optimized control force proﬁle—Stewart platform (problem P2).

(20)

UNCORRECTED

PROOF

478 of the difference equations, and which can be effectively reduced by increasing the sampling points. Further- 479 more, although no attempt has been made in this paper to prove the convergence of the MAP algorithm the- 480 oretically, it is found from numerical experiments that the robustness of the algorithm is dependent mostly on 481 the feasibility of the guessed initial trajectory. The method presented in Appendix B ensures the initial trajec- 482 tory is within the bounds of the kinematic constraints and no convergence difficulty has been encountered for 483 any of the problems tested. It should, however, be noted that if there are any singular points (i.e. where the 484 determinate of the overall transformation Jacobian matrix [J

p

] is zero) occurred within the initial approxima- 485 tion or any of the subsequent linearized trajectories, then the iterative procedure can not be continued and the 486 trajectory should be altered or perturbed by enforcing additional kinematic constraints at the neighborhood of 487 these points.

0 2 4 6 8 10

-500 -400 -300 -200 -100 0 100 200 300 400 500

Time (Sec.)

Driving Force (N.)

1

2

3

Upper Bound

Lower Bound

0 2 4 6 8 10

-500 -400 -300 -200 -100 0 100 200 300 400 500

Time (Sec.)

DrivingForce(N.)

4

5 6

Upper Bound

Lower Bound

τ τ

τ

τ τ

Fig. 14. Optimized control force proﬁles—Hexaglide (problem P2).

(21)

UNCORRECTED

PROOF

488 It should also be pointed out that the dynamic model of the system established in this study is based on the 489 assumptions of rigid links and perfect joints. In real applications, the response of the system will be affected by 490 link flexibilities and joint backlashes. In addition, due to the iterative nature of the computational algorithms, 491 the computing time is dependent not only on the CPU speed but also on the number of iterations, and which 492 in turn is dependent on the geometric complexity of the specified motion contour. Therefore the presented 493 methods may not suitable for real-time applications, yet the optimal driving forces evaluated off-line can be 494 served as feed-forward terms in an appropriate on-line control strategy.

495 8. Conclusions

496 In this paper, a systematic approach for formulating the dynamic equations of fully parallel, linear actuated 497 platform type manipulators is presented. The generalized coordinates used in this approach are the Cartesian 498 coordinates of the moving platform, and hence the resulting equations of motion can be easily transformed 499 into a parametric form in terms of the task space redundant DOF of the robot.

500 Based on this formulation, two optimal dynamic trajectory planning problems associated with the task 501 space redundant DOF have been studied. The ﬁrst problem is of theoretical interest as it would provide 502 knowledge of the optimal distribution of the cutting force along any speciﬁed machining contour. The solution 503 to the second problem is useful in a practical sense, as it would not only maintain a maximum constant cutting 504 force for the entire machining process but would also eliminate possible discontinuities and sudden jumps of 505 control forces to avoid structural damage.

506 The numerical algorithms developed for solving these problems are based on the method of approximate 507 programming, and taking into account the non-linear bounds imposed on the control forces and the rotation 508 angle of the moving platform. The eﬀectiveness of these algorithms is demonstrated by two numerical 509 examples.

510 Acknowledgement

511 The ﬁnancial support of this work by the National Science Council of the Republic of China under Grant 512 94-2212-E-011-007 is gratefully acknowledged.

0 50 100 150 200 250 300 350

0 2 4 6 8 10 12

Number of Time Segment

Total Work (J.)

Stewart Hexaglide

Fig. 15. Comparison of computational accuracy versus sampling rate.

(22)

UNCORRECTED

PROOF

513 Appendix A

514 As shown in Fig. 16(a) and (b), each limb of the Stewart platform and the Hexaglide consist of two links.

515 For the Stewart platform, the virtual displacement of the mass center of two links can be, respectively, written 516 as

517 dP

1i

¼ dW

i

q

_1i

l

_i

ðA:1Þ

519 519 520 and 521

dP

2i

¼ dP

i

dW

i

q

_2i

l

_i

; ðA:2Þ

523 523

524 where l

i

is a unit vector along the direction of the limb and dW

i

is the virtual rotation of the limb. d P

i

is the 525 virtual displacement of the center point of the spherical joint p

i

attached to the moving platform, and which is 526 related to the virtual rotation and virtual displacement of the center of mass of the moving platform by 527

dP

_i

¼ dP þ dW r

i

: ðA:3Þ

529 529

530 On the other hand, since the absolute position of joint p

_i

can also be expressed as

q_i

i

ρ2

bi

Link 1

Xi

Yi

Zi

Lin k 2

O

r

i

l

i i

P

e_i B_i P_1i

P_2i P_i

P p_i

ri

bi

Xi

Yi

Zi

pi

Link 2

Link1

X Z

Z

Y O Y

i

ρ

2

i

ρ

1

q

i

l

i

P

P_1i P_2i

P_i P

B_i

a

b

Fig. 16. Schematic diagram and geometrical parameters of the limbs: (a) Stewart platform, (b) Hexaglide.