Alexandre Dolgui
email: dolgui@emse.fr

Anatol Pashkevich
email: anatol.pashkevich@emn.fr

Prof Alexandre Dolgui

Manipulator motion planning for high-speed robotic laser cutting

Keywords: PATH PLANNING, OPTIMIZATION, SIMULATION PATH PLANNING, OPTIMIZATION Robot path planning, Laser beam cutting, Multiple-criteria optimisation Robot path planning, Laser beam cutting, Multiple-criteria optimisation

Recent advances in laser technology, and especially the essential increase of the cutting speed, motivate amending the existing robot path methods, which do not allow the complete utilisation of the actuator capabilities as well as neglect some particularities in the mechanical design of the wrist of the manipulator arm. This research addresses the optimisation of the 6-axes robot motions for continuous contour tracking while considering the redundancy caused by the tool axial symmetry. The particular contribution of the paper is in the area of multi-objective path planning using graph-based search space representation. In contrast to previous works, the developed optimisation technique is based on the dynamic programming and explicitly incorporates the verification of the velocity/acceleration constrains. This allows the designer to define interactively their importance with respect to the path-smoothness objectives. In addition, this optimisation technique takes into account the capacity of some manipulator wrist axes for unlimited rotation in order to produce more economical motions. The efficiency of the developed algorithms has been carefully investigated via computer simulation. The presented results are implemented in a commercial software package and verified for real-life applications in the automotive industry.

Laser cutting processes have become increasingly important in a wide range of industrial applications due to their many advantages over other cutting methods. The key features of this technology include high speed, narrow kerf width, small thermal distortions and very low roughness on the edges that allow non-contact processing of complex 3D shapes without additional finishing operations. To ensure high productivity and flexibility, the laser is usually manipulated by either a CNC machine or a 5/6-axis industrial robot, which is programmed off-line using available commercial CAD/CAM systems [START_REF] Nof | Handbook of Industrial Robotics[END_REF][START_REF] Gropp | Laser beam cutting[END_REF][START_REF] Defaux | Welding and cutting: Robots take the lead in 3D[END_REF][START_REF] Ion | Laser processing of engineering materials: principles, procedure and industrial applications[END_REF].

However, recent advances in the laser cutting technologies motivate enhancements of existing programming methods and encourage developments of dedicated optimisation algorithms that are in the focus of this paper.

The most important aspect of this problem is related to essential increase of the technologically admissible limit for the cutting speed. At present, depending on the laser configuration, and also on a type and thickness of the processing material, the cutting speed can reach more than 50 m/min. This is comparable with the kinematic capabilities of modern industrial robots, which are used to manipulate the laser beam [START_REF] Steen | Laser material processing[END_REF][START_REF] Rooks | Laser processing of plastics[END_REF][START_REF] Schlueter | Advances in industrial high power lasers[END_REF]. For instance, for a typical 6-axis anthropomorphic robot, the forearm axis speeds range from 200 to 300 deg/sec and the wrist axes may be actuated with speed from 300 to 500 deg/sec. However, the maximum Cartesian velocity of the tool is usually lower than 1 m/s and significantly varies over the workspace, imposing critical constrains for the implementation of the high-speed laser cutting process. Hence, in this area, the problem of optimal motion planning for a multi-axe manipulator becomes crucial.

However, as follows from our experience, direct application of existing algorithms may produce unfeasible results with non-smooth velocity/acceleration profiles (in spite of a quite satisfactory performance for the lower desired cutting speed). For example, as reported in [START_REF] Ghany | Using a Nd:YAG laser and six axes robot to cut zinc-coated steel[END_REF], some existing path planning methods may produce cutting head vibrations at high speeds. their analysis. In Section 6, an industrial implementation is reported. Finally, Section 7 summarizes the main contributions and suggests some prospective research directions.

Related works

Since the debut of robotic laser cutting, most of the related research focused on robot off-line programming (Backes et al., 1995;[START_REF] Kaierle | Design and manufacturing tools for laser beam processing[END_REF]. The main reason is that, for this technology, the conventional "manual teaching" method proved to be very tedious and time-consuming, requiring a temporary exclusion of the robot from the manufacturing process and performing the operatorcontrolled movement of the tool along the cutting path, which has to be properly marked beforehand. In contrast, off-line programming can generate the control code by means of computer graphics, away from the factory floor. Thus, the down time may be significantly reduced, enabling very small batch sizes to become economically feasible [START_REF] Sendler | Offline-programmiertes Laserschneiden im Automobilbau[END_REF][START_REF] Wittenberg | Developments in offline programming: An overview[END_REF][START_REF] Mitsi | Off-line programming of an industrial robot for manufacturing[END_REF].

At the moment, there are a number of commercial robotic off-line programming/simulation systems on the market. Some of the more commonly used are em-Workplace (RobCAD) from Tecnomatix Technologies (www.tecnomatix.com), IGRIP from Delmia (http://www.delmia.com), CimStation from Silma (www.silma.com),and Workspace (www.workspace5.com). They offer 3D graphical simulation environments with visual programming capabilities and a wide range of convenient design tools, such as robot selection, robot placement, motion simulation, cycle time calculation, and collision-free path planning. Some of these systems include process-specific modules for laser cutting. In particular, the "Robcad/Cut&Seal" module enables the semiautomatic creation of cutting contours and their conversion into linear, circular and spline motions. A similar software tool is also available from Alma (www.alma.fr) --"act/cut3D", which possesses some additional convenient features, including defining the head orientation for bevelled cutting, imposing "technological" tolerances, and interactive modification of the tool orientation for each point of the programme. However there still exists a considerable gap between the capabilities of the commercial robotic CAD systems and the requirements of a particular application. Up to now, the robot programs for many 3D cutting applications (especially for high-speed cutting) are constructed interactively and then are verified using simulation/visualization tools incorporated in the corresponding CAD system. The ultimate goal is the generation of reliable programs automatically from designs and drawings, similar to CNC-machine programming methods [START_REF] Sun | A modified analytical model for optimization of NC-tool cutting path[END_REF]Tsai, 1994, Bohez et al., 2000;[START_REF] Xu | Tool-path generation for five-axis machining of free-form surfaces based on accessibility analysis[END_REF][START_REF] Kim | Integrated machining tool path planning using feature free spaces[END_REF][START_REF] Anotaipaiboon | Tool path generation for five-axis NC machining using adaptive space-filling curves[END_REF]Wang andXie, 2005, Bi and[START_REF] Bi | Automated robotic programming for products with changes[END_REF].

One of the first contributions in the off-line programming for the laser cutting robots was done by M. Geiger and his co-workers (University of Erlangen-Nuremberg, Germany). They developed technology-oriented approaches, which simultaneously consider the part geometry, process parameters and some properties of the machine tool [START_REF] Geiger | Laser beam cutting with solid state lasers and industrial robots[END_REF]Geiger and Kolléra, 1994;[START_REF] Bauer | Offline-programming for 3D-laser systems[END_REF]. However, their techniques may be applied to the non-redundant kinematic structures only and based on the five-axis robots.

For manipulators with six degrees of freedom (which are redundant with respect to 3D cutting requirements), the motion planning problem was firstly addressed by Abe, Shibata and Tanie [START_REF] Abe | Motion planning for 3D cutting by a manipulator with 6 degrees of freedom -optimization by genetic algorithm[END_REF][START_REF] Shibata | Motion planning by genetic algorithm for a redundant manipulator using a model of criteria of skilled operators[END_REF]. These authors proposed a genetic algorithm (GA) that optimises the cutting tool orientation using the evaluation function extracted from the experience of skilled operators. Informally, the operator preferences were defined as follows: "Motion of three basic (1st to 3rd) links with low response should be least so as to minimize movement of centre axis of the end-effector." Formally, the local optimisation criteria was defined as the weighted squared sum of all axis velocities (from 1st to 6th) together with the velocity of the reference point for the 5th axis.

Furthermore, for the overall optimisation with a global criterion, all local criteria were summed over the path. As follows from the presented results, the authors succeeded in the generation of manipulator motions for a relatively slow cutting speed (2m/min). However, there are a number of open questions with the application of this technique to current real-life problems (choosing of weighting factors, taking into account manipulator velocity/acceleration constraints and joint limits, as well as the convergence of algorithms for high-dimensional input data exported from CAD, etc.). Some recent techniques that employ the general end-effector constraints concept [START_REF] Yao | Path planning with general end-effector constraints[END_REF]) also suffer from this drawback.

The above mentioned approach has been enhanced in our previous paper [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF], which presents a detailed study of multiple optimisation objectives and their incorporation in the global criterion. Compared to other works, our motion planning algorithm possesses essential advantages, since it is based on a graph representation of the search space and on the dynamic programming. The latter ensures a much higher computation speed than the GA, which is very important for industrial applications. This algorithm was successfully implemented on a factory floor and later verified in the automotive industry. Nevertheless, recent advances in laser technology and the significant increase of the technological limit on the cutting speed inspire further improvement of our method.

In the area of robotic motion planning, it is also worth mentioning some other works that are not directly related to cutting, but contain some of theoretical background used in this paper. In robotic literature, a lot of research focuses on the time-optimal (or energy-optimal) motion planning with obstacle avoidance for (i) point-to-point and (ii) specified-path formulations. These researches originated from pioneering papers of [START_REF] Bobrow | Time-optimal control of robotic manipulators along specified paths[END_REF], Shin and Mckay (1985), [START_REF] Pfeiffer | Concept for manipulator trajectory planning[END_REF][START_REF] Pfeiffer | Concept for manipulator trajectory planning[END_REF][START_REF] Shiller | On computing the global time-optimal motions of robotic manipulators in the presence of obstacles[END_REF][START_REF] Shiller | On computing the global time-optimal motions of robotic manipulators in the presence of obstacles[END_REF]. In most of these publications, redundancy is used to optimise secondary objectives, such as manipulability measures or distance to singularities.

Comprehensive reviews on this topic can be found in [START_REF] Chettibi | Synthesis of dynamic motions for robotic manipulators with geometric path constraints[END_REF][START_REF] Antonelli | SmartMove4: An industrial implementation of trajectory planning for robots[END_REF].

Another related research area is the robotic motion sequencing (scheduling) where most work concentrates on specific non-trivial cases of the travelling-salesman problem (TSP). For robotic applications, the scheduling problem with a minimum-time objective was first addressed by [START_REF] Rubinovitz | Task level off-line programming system for robotic arc welding-an overview[END_REF]. Later, the TSP-based method was enhanced by [START_REF] Kim | Scheduling for an arc-welding robot considering heatcaused distortion[END_REF], who suggested a number of heuristics and a problem-oriented genetic algorithm for robotic welding. For other robotic applications, task scheduling concentrates mainly on mechanical assembly and processing, as well as grasping-type tasks [START_REF] Yuan | An integration of robot programming and sequence planning[END_REF]Gu 1999, Ben-Arieh et al. 2003). Some recent results are presented in [START_REF] Bagchi | A review of TSP based approaches for flowshop scheduling[END_REF][START_REF] Dolgui | Cluster-level operations planning for the out-ofposition robotic arc-welding[END_REF]. It is obvious that the TSP-based formulations can not be directly used for robotic cutting, where both the processing contour (i.e. the motion sequence) and the tool speed are specified. However, there are some useful modifications of discrete dynamic programming [START_REF] Jouaneh | Trajectory planning for coordinated motion of a robot and a positioning table --II: Optimal trajectory specification[END_REF][START_REF] Lee | Dynamic programming approach to near minimum-time trajectory planning for two robots[END_REF], which are used in this paper.

Another potential approach to the considered problem is based on the standard redundancy resolution techniques, which usually is founded on the generalised inverse of the manipulator Jacobian or the priority task decomposition [START_REF] Whitney | Resolved motion rate control of manipulators and human prostheses[END_REF][START_REF] Klein | Review of pseudoinverse control for use with kinematically redundant manipulators[END_REF][START_REF] Yoshikawa | Basic optimization methods of redundant manipulators[END_REF]Doty et al., 1996). Here, the basic idea is to introduce additional objectives/constraints to obtain a well-defined problem and then solve it by conventional methods. However, in robotic cutting, applying the kinematic-based redundancy resolution methods may engender chaotic joint motion and high accelerations, which leads to erratic behaviour with unpredictable arm configurations [START_REF] Duarte | Motion chaos in the pseudoinverse control of redundant robots[END_REF]. In addition, these methods are applied independently (locally) to different points of the path and, consequently, they often fail to get a global optimal solution [START_REF] Hwang | Global motion planner for curve-tracing robots[END_REF].

Summarising the analysis of the related works, most of existing path planning techniques operate with a single objective, such as the path length or travel time. There are few papers which deals directly with multi-objective path planning, and they are mainly for mobile robot navigation [START_REF] Fujimura | Path planning with multiple objectives[END_REF] or surface manufacturing [START_REF] Chen | Multi-objective Optimal Robot Path Planning in Manufacturing[END_REF]. To our knowledge, only our previous paper [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF] addresses this issue to robotic cutting. Finally, current technological advances motivate further enhancements of corresponding path planning techniques which are in the focus of the paper.

Problem statement

Generally, a basic problem in automated robot programming for laser cutting is to find optimal manipulator motions from given CAD models of the (i) processing part, (ii) cutting tool, (iii) robot and (iv) work cell environment in order to satisfy specified design objectives and constraints [START_REF] Chedmail | La CAO en Robotique : outils et methodologies[END_REF]. Based on these data, the trajectory planner generates an optimal tool path in both the Cartesian and joint coordinate space. This path is further validated by a simulation module.

Finally, it is translated into a robot control code. The core for the program generation tool is a set of optimisation routines that are addressed in this section.

Cutting task model

Let us assume that the desired Cartesian path, along which the cutting tool is to be moved, is imported from a CAD system and described by two vector functions as follows:

{ } n , , , k T t k t t t , t K 2 1 0 ; ] [0, ; 1 ) ( : ) ( ) ( = = = = n n p C ( 1
)
where t is a scalar argument (time); p(t) R 3 defines the Cartesian coordinates of the tool tip, and n(t) R 3 is the unit vector of the tool axis direction, which must be normal to the part surface (figure 1). These data can be directly extracted from the graphical model of the part, by defining the processing contour as an "augmented line". for example. The path is assumed to be closed and time-uniformly sampled into the sequence of nodes } 0,... { ; ,

n k k k = n p
, where the first and the last coincide (

n p p = 0 ; n n n = 0
), and the time-interval length is t .

[Insert figure 1 about here]

To describe the tool spatial location, let us also define the Cartesian displacement along the path

; 1 k k k p p p = + 1 2 1 0 = n , , , k K and introduce a unit direction vector || || k k k p p a =
, which is tangent to the part surface and to the direction of the points for the tool motion. For the last node, let us define this vector as 0 a a = n . Then, assuming that the vectors a k and n k are mutually orthogonal, each node may be associated with the Cartesian frame in which the x-axis is directed along the path, the z-axis is directed along the cutting tool, and the y-axis is computed in such way that these three axes form a right-handed coordinate frame. The corresponding homogenous transformation matrix is

4 4 1 0 0 0 × × = k k k k k k p n a n a H
(2)

where "×" denotes the vector product. It should be noted that, generally, the CAD-imported vectors n k and the computed vectors n k may be slightly non-orthogonal. In this case, it is necessary to adjust {a k } using the classical Gram-Schmidt orthogonalisation:

|| || ; ) ( k k k k k T k k k a a a n a n a a .
The introduced frame sequence } 0,...

{ n k k = H
is used as a pivot for defining the complete pose of the robotic tool, which is usually determined by six independent parameters (three Cartesian coordinates and three Euler angles). However, since the cutting tool is axially symmetric, the frames H k can be rotated around corresponding z k -axes without any influence on the technological process. This one-dimensional redundancy leads to an infinite set of admissible tool locations described by the matrix product

] , (- , ) ( ) ( = k k k z k k H R L , ( 3
)
where k is an arbitrary scalar parameter and R z ( ) is the standard z-axis rotation matrix.

Another source of redundancy is related to the manipulator posture µ (or the configuration index), which is required for the unique mapping from the task space to the joint coordinate space.

So, in total, the robotic task is described by a sequence of locations (3), while the design parameters are represented by the sequence } 0,...

{ n k k , k = µ
. At this step, the design problem can be formulated in the terms of non-linear programming; however this straightforward approach is not prudent because of high search space in this problem.

Constraints

While planning the robot motions for laser cutting and other curve-tracking applications, a designer must take into account three types of constraints: task, robot kinematic and collision constraints [START_REF] Hwang | Global motion planner for curve-tracing robots[END_REF]. Task constrains are expressed in the terms of the tool position/orientation required to perform the assigned task and, for the considered problem, they are described by expressions (3). Robot kinematics constraints are caused by the manipulator geometry and expressed as workspace limits as well as limits on the range the manipulator joint coordinates.

And collision constraints arise from the need to avoid collisions between the robot links, workpiece and workcell components.

To define the kinematic constrains more precisely, let us briefly review some results from robot kinematics. For the anthropomorphic serial manipulators that are usually employed with laser cutting, the direct kinematic model, which defines the tool location L corresponding to the joint coordinate vector q , may be always expressed in a closed form, as

base i i i i tool q T T T q L = = 6 1 1 ) ( ) ( , (4)
where i-1

T i is the transformation matrix from the (i-1)th to the ith link, q i is the corresponding joint coordinate, the matrix T tool defines the tool tip location with respect to the manipulator mounting flange, and the matrix T base defines the robot base location in the world coordinate system.

Particular expressions for the homogenous matrices i-1 T i for various manipulators can be found in common reference books [START_REF] Ceccarelli | Fundamentals of mechanics of robotic manipulation[END_REF][START_REF] Spong | Robot modeling and control[END_REF], usually they are expressed using the Denavit-Hartenberg (DH) notation:

4 4 1 1 0 0 0 ) os( ) in( 0 ) in( ) in( ) ( ) ( ) os( ) in( ) os( ) in( ) in( ) ( ) in( ) ( ) ( × = i i i i i i i i i i i i i i i i i i i i d c s q s a s q cos cos q c q s q c a s q s cos q s q cos q T F o r P e e r R e v i e w O n l y
where a i , d i , i are the DH-parameters associated with ith link and ith joint (the offset distance, link length, and twist angle, respectively).

However, the inverse transformation (from the task coordinates space to the manipulator joint space) is generally non-unique and may not be expressed in a closed form. In practice, robot designers prefer special types of manipulator architecture, which admit closed-form solutions. Otherwise, there exist numerical techniques to deal with this problem [START_REF] Manocha | Efficient inverse kinematics for general 6R manipulators[END_REF][START_REF] Pashkevich | Real-time inverse kinematics for robots with offset and reduced wrist[END_REF][START_REF] Husty | A new and efficient algorithm for the inverse kinematics of a general serial 6R manipulator[END_REF]. Within our study, and independent of the solution type (analytical or numerical), the inverse kinematical transformation will be presented in an abstract form as

) ( 1 µ , c L f q =
, where µ is the configuration index that determines the manipulator posture, and the set contains all combinations of admissible postures (shoulder right/left, elbow up/down, wrist plus/minus). For most of 6-axis industrial manipulators, the number of different configurations is equal to 8, but generally this can rise up to 16 [START_REF] Manseur | Robot manipulator with 16 real inverse kinematic solution sets[END_REF].

It should be stressed that typical robot controllers do not allow changing the manipulator configuration while mowing between successive nodes and the use of the Cartesian space on-line interpolation (LIN and CIR commands). But for cutting applications, this specific kinematic constraint can be released by temporarily switching to the joint-space interpolation for these specific path segments (PTP command). So, in contrast to our previous work [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF]), here we do not use separate manifolds for each value of µ. This also causes some modifications at the final stage, when the optimal path is translated into a robot control code.

Using the inverse kinematic transformation, the sequence of the admissible tool locations ( 5) can be mapped onto the joint coordinate space as follows:

{ } n k , C k k k k z c k k k Q K 1 ], , (- ) , ) ( ( ) , ( 1 = = = µ µ µ H R f Q (5)
considering both the workspace dimension constraints (i.e. inverse kinematics existence) and the joint limits

max i i min i q q q < < , I i , 6} {1, K = q I
. For computational convenience, the constraint violation case is denoted as

= ) ( k k

Q

(this notation is easily implemented in data structures From a practical point of view, the desired path planning algorithm should produce "smooth motions at reasonable speeds and at reasonable accelerations". Hence, many of aforementioned authors impose constraints on the joint torques. However, addressing the torque constraints requires a rather accurate dynamic model, which is problematic in real-life industrial projects. A practical solution consists of transforming torque constraints into acceleration ones using the manipulator specification data provided by the manufacturer.

Within the frames of the discretised path presentation (1), the joint velocity/acceleration constraints may be expressed via the finite-difference approximation as:

) ( 1 v i v k , i k , i q q q < , ( 6
) ) ( 2 1 2 a i a k , i k , i k , i q q q q < + , ( 7
) where t q q max i v i = & ) ( ; 2 ) ( t q q max i a i = & & ; the notations

Design objectives

In the qualitative terms, the desired manipulator motion should be as smooth as possible while satisfying the contour-tracking and actuator-dependent constraints ( 7) -( 9). This means that the qualitative performance measures should be based on some type of the "smoothness" and "economy" measures applied to all manipulator joints (Edan and Nof, 1997;[START_REF] Shoval | Analysis of robot motion performance and implications to economy principles[END_REF].

For the considered problem, which is based on the discrete path presentation, the degree of smoothness can be evaluated by the following criteria:

• total displacement ( ) % = = n k k k k , i k k k , i s i q q J 1 1 1 1 ) , ( - ) , ( ) , ( µ µ µ , (8)
• maximum increment (speed)

( ) ) , ( - ) , ( ) , ( 1 1 1 = k k k , i k k k , i k v i q q max J µ µ µ , (9)
• coordinate range

( ) [ ] [ ] ) , ( ) , ( ) , ( k k k , i k k k k , i k i q min q max J µ µ = µ , ( 10
)
where i is the joint number, and , µ are the vectors composed of the design parameters 0 , 1 , … n and µ 0 , µ 1 , …µ n respectively.

A geometrical explanation of these performance measures is the following. The displacement criterion evaluates the total amount of joint motions (without regard to the direction of the motion).

The increment characterizes the coordinate variation during the time-interval t and is proportional to the speed. Finally, the range evaluates the extent of the coordinate variation during the motion.

Intuitively, the minimization of each of these criteria should lead to a smoother generated path.

However, as follows from related studies, the objectives ( 8) -( 10) may be in competition with each other, i.e. minimizing one of them may degrade the performance in others. Besides, these assessments are computed for each joint coordinate. Thus, the resulting performance measures form a vector. The designer must choose one of the techniques that are usually used for multi-criteria problems (Steuer, 1985;[START_REF] Cheng | Multiple-goal priority considerations of redundant manipulators[END_REF][START_REF] Ehrgott | Multicriteria Optimization[END_REF]. In this paper, instead of giving preference to a particular criterion or optimisation technique, all the following options are admitted: (i) defining priority of partial objectives or the primary objective; (ii) applying minimax technique, i.e. the worst-case optimisation; (iii) assigning weights to combine multiple criteria in a linear function. Independent of the chosen option, the corresponding vector-optimisation engine must include the scalar-optimisation routines that are developed in the following section.

Path planning algorithm

Search space discretisation

A practical method to obtain an optimal (or feasible) solution under the constraints described in sub-section 3.2 is the sampling of the allowable domain for the redundant parameter . Such an approach transforms the continues search space into an acyclic directed graph, with the vertices uniquely representing the matrix of the tool location L and the vector of joint coordinates Q.

Let us assume that the parameter [-, ] is sampled with the step = 2 /m 0 , m 0 Z and, for each discrete value of and each admissible configuration µ , the locations (5) are tested for kinematic and collision constraints. Then, the admissible locations, which satisfy the constraints, are included in the search space in such manner that each Cartesian-path node {p k , n k } produces several elements of the data structure {L, Q}:

U k m j j , k j , k k k 1 , } { = ( ) * + , - . Q L n p (11)
where m k is the number of the successful locations for the kth node. It should be noted that the "path-smoothness" constraints (6), ( 7) are not tested at this stage yet, since they are associated with several successive locations. The test result can be presented in the plane t×( , µ), where the abscissa (index k) corresponds to the time t, and the ordinate (index j) corresponds to the combination of the redundant parameter and the configuration index µ (see figure 2).

[Insert figure 2

U U k m j j , k j , k n k V 1 0 = = ( ) * + , - = Q L (12)
and edges

U U U 1 1 1 1 1 1 = = = ( ) * + , - ( ) * + , - = k m l l , k l , k j , k j , k k m j n k , E Q L Q L , (13)
which connect only successive tool locations from the entire set (5). Hence, the considered pathplanning task is reduced to the following network optimisation problem.

Design Problem. For given set of vertexes V and set of edges E, find the "best" path of length n

/ ( / ) * / + / , - . / ( / ) * / + / , - . / ( / ) * / + / , - = 0 n j , n n j , n j , j , j , j , n j j j Q L Q L Q L K K 1 1 1 1 0 0 0 0 1 0 ) , , ( (14)
with initial and final states

U 0 1 0 0 0 m j j , j , V = ( ) * + , - Q L and U n m j j , n j , n n V 1 = ( ) * + , - Q L ,
which minimizes a specified performance index.

It should be stressed that, in this formulation, both the initial and final states are not unique, but the problem can be transformed to the classical problem by adding the virtual start and end nodes.

Besides, the initial and the final states may not coincide with each other, as it follows from the technical nature of the problem. It is also worth mentioning that all nodes are solvable for the inverse kinematics and admissible for the collision tests (otherwise they are excluded from the graph (V,E) at the stage of the graph generation).

To solve this network optimisation problem, it is also necessary to define the distance matrix, which assigns certain lengths to the graph edges. Since in the Cartesian space the tool path is sampled uniformly, the distance should be based on a joint space metric. While defining this metric, it is necessary to take into account that the actuator capacities (i.e. maximum speed, acceleration, etc.) are different for different manipulator axes. So, the displacement components q i corresponding to the joint displacement vector, which are ) , ( 6

1 q q = K q
, should be weighted.

The most appropriate way of assigning the weights, is using the maximum axis speed max i q &

specified by the manufacturer. In this case 1 -

) ( t q w max i i = &
and the product 1 q w i shows the degree of actuator capacity utilisation with respect to speed. Finally, combining the evaluations for all axes, the distance in the joint space may be defined using the Euclidean, Manhattan or Chebychev metrics [START_REF] Deza | Dictionary of distances[END_REF] that are usually used because of both the computational convenience and their clear physical meaning [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF]:

% = = 6 1 2 ) ( ) ( i i i E q w q 2 ; % = = 6 1 | | ) ( i i i M q w q 2 ; | | max ) ( i i i C q w = q 2 (15)
While applying the above expressions, it is necessary to consider that for cutting applications, the 4 th and 6 th axes may allow unlimited rotation. In this case, before computing the distance 2(.), the differences i q must be minimised in accordance with the following expression:

{ } 6} {4, 2 + = i , p q min q i Z p i ( 16
)
where p is an integer number. For example, for the axes 1, 2, 3, 5, the movement from ° = 120

) (k i q to °+ = + 120 ) 1 (k i q corresponds to the joint coordinate displacement °+ = 240 i q
. However, for the axes 4 and 6 it is equal to

° = 120 i q
and the target joint coordinate is must be changed to

° = + 240 ) 1 (k i q
. In practice, there is an open question: which metrics of (.) (.), (.),

C M E 2 2 2
, should be chosen? So, it is proposed to give some freedom to the designer and to include the choice of the corresponding option in the set of parameters for the optimisation algorithm.

Generation of optimal path

Since there is no combinatorial optimisation technique, which is able to solve this multiobjective problem directly, the vector performance measures ( 8) -( 10) should be converted into an aggregate scalar criterion. At this step, a relevant metric is applied to the sequence of the path segments whose length is computed using one of the expressions (15). Hence, the desired scalar criteria may be based on the non-weighted Manhattan/Chebychev norm. The weights w i , which are used for calculating the segment distance 2(.), are altered between optimisation runs to produce a set of Pareto-optimal solutions.

Corresponding expressions for the objective functions may be written as follows:

( )

% = = n k k k a s j k j k J 1 1 ) ( ) , 1 ( - ) , ( ) , ( Q Q µ 2 (17) (
)

) , 1 ( - ) , ( ) , ( 1 ) ( = k k a k j k j k max J Q Q µ 2 (18) where } , , { (.), C M E a a 2
is the distance function, the notation j k defines the values of the redundant parameters , µ at the kth node and, for further convenience, the vectors of the joint coordinates Q k,j are denoted as Q(k,j).

For the additive objective ) (s J , an optimal path can be found by means of dynamic programming. This is a standard shortest-path problem in an acyclic directed graph. Nevertheless due to the graph specific properties (node clustering), it is not reasonable to apply known algorithms (Dijkstra's , Bellman-Ford, etc.). This problem is also closely related to the generalised travelling salesman problem (GTSP) employed in the motion and process planning for which there are some efficient optimisation routines, see for example [START_REF] Ben-Arieh | Process planning for rotational parts using the generalized travelling salesman problem[END_REF][START_REF] Dolgui | Cluster-level operations planning for the out-ofposition robotic arc-welding[END_REF].

To describe the proposed algorithm in terms of graph theory, let us introduce the following

notations. Let ) ( 1 0 E , S S S G n K U U =
be an acyclic directed graph with the set of vertices

k n k S V 0 = = U
partitioned into n+1 mutually exclusive and exhaustive subsets (clusters) ). This is in agreement with usual engineering practices.

} ..., 1 { k kj k m , , j v S = = of sizes k k m S = | | .
For such a formulation, a vertex may be coded as pair (k, j), an edge is identified by quadruplet (k 1 , j 1 , k 2 , j 2 ), and the solution is defined by the array of j-indices {J opt (k), k=0,1,…n} of the sequential visiting vertices. Also, following this notation, the cluster sizes are defined in the array {J max (k), k=0,1,…,n} and the cost of the edges are denoted as 2 (k 1 , j 1 , k 2 , j 2 ).

A basic expression for the developed algorithm is derived in the following way. Let us consider a reduced order sub-problem with k clusters S 0 , S 1 , … S k-1 , and let

. . . k S S S K that links a vertex 1 1 k j , k S v to the nearest vertex 0 S u
, assuming that all the clusters are visited exactly once. Then, using the dynamic programming, the optimal solutions for the sub-problem with k+1 clusters S 0 , S 1 , … S k can be found by choosing the best edge ) , (

j , k v u , 1 k S u
that links the vertex An outline of the path planning procedure based on this technique is presented below. The procedure consists of five basic steps where the first two, (1) and (2), implement the recursion (19).

k j , k S v with the cluster S k-1 : { } n , , k , v v d min d j , k p , k p , k p j , k K 2 1 ) , ( 1 1 = + =
These steps deal with creating matrices of optimal distances M dist (k,j) and the corresponding matrix M ind (k, j) of the optimal j-indices of the preceding clusters. In steps (3) and ( 4), the minimum value in the nth column of the distance matrix M dist (.) is computed, in order to select the best vertex in the 5), the best vertices of the preceding clusters are iteratively extracted from the matrix M ind (.) to generate the optimal path described by the array of the optimal indexes J opt (.). The procedure also uses the array J max (.) containing the upper range for the j-index and the vertex distance function 2(.) defined previously. The notation d 6 defines the infinity (or the largest number for computer implementation).

Procedure: Path planning (dynamic programming)

(1) For j = 1 to J max (0) do Set M dist (0, j) := 0; M ind (0, j) := 0 (2) For k = 1 to n do For j = 1 to J max (k) do A distinct feature of the developed procedure, which takes into account the problem-specific requirements, is contained in sub-steps (2b7) and (2b ) that verify the path-smoothness constraints (6) and ( 7). Sub-step (2b7) incorporates the function ) (k,j,p f v , which evaluates the velocity constraints (6) for all manipulator axes while the manipulator moves from the location L k-1,p to the location L k,j (a non-zero value of this function corresponds to violating one of the velocity constraints). Similarly, at the sub-step (2b ), the function ) ( l , k,j,p f a verifies the acceleration constraints (7) for the location sequence L k-2,l , L k-1,p , L k,j , where the location L k-2,l is assumed to be . It is obvious that this is a penalty-function approach. The definition of the path-segment of the double length strictly corresponds to the optimality principle incorporated in dynamic programming.

(a) Set d min := d 6 (b) For p = 1 to J max (k-1) do ( ) Set d cur := M dist (k-1, p) + 2(k, j, k-1, p) (7) If k > 0 & f v (k, j, p) ! 0 Set d cur := d 6 ( ) If k > 1 & f a ( k, j, p, M ind (k-1, p) ) ! 0 Set d cur := d 6 (
The computational complexity of this algorithm is quite acceptable for practical applications. It is evaluated as O(n m 2 ), where n is the number of clusters and m is the maximum cluster size. This estimation also includes the calculation of the cost components. This calculation is executed for the neighbouring clusters only.

A similar algorithm can be applied for the minimax design objective ) ( J . The only modification needed deals with the sub-step (2b ), which must be re-written in accordance with the following expression:

{ } n , , k , v v d max min d j , k p , k p , k p j , k K 2 1 ) , ( ; 1 1 = = 5 5 2 (20)
that is also based on the dynamic programming principle. Combining two of the path generation options (objectives ) (s J and ) ( J ) and altering the distance metrics weights w i together with the constraint weights v and a (see expressions ( 6) and ( 7) ), the designer can generate a collection of the candidate solutions to be included in the Pareto-optimal set.

Generation of control programs

As assumed in previous sections, the sequence {n k , p k }, describing the tool path in the Cartesian space, is extracted from the CAD workpiece model by uniformly sampling the cutting contour. In performing such a conversion, the designer must specify an appropriate distance S between nodes, which should be small enough to ensure desired accuracy of the contour approximation. However, during off-line programming, there exists a lower bound on the time-sampling step, which is determined by the control unit parameters. Let us examine this problem in detail.

In industrial robots, the control code describes the sequence of the path segments (with a specified interpolation type within each: linear/circular for Cartesian space and linear for joint space). Some advanced robot controllers implement also spline interpolation. It is obvious, that for tracking the cutting contour, the Cartesian space interpolation is used the most, while the joint space interpolation is utilised only for travelling through a singularity.

To reduce the control program size, see for example [START_REF] Korb | Data reduction for manipulator path planning[END_REF], the intermediary nodes within the line/circular segments are not specified. Instead, they are generated on-line using the trapezoid velocity profiles that typically include three sections (acceleration, uniform motion, and deceleration). Their duration depends on the desired displacement and the velocity/acceleration constraints imposed on both each joint variable and on the Cartesian coordinates. Moreover, in the continuous motion mode (i.e. without stopping at the end of the segment), the successive segments are joined in such way that the acceleration section of the succeeding segment coincides with the deceleration section of the preceding one. As a result, the velocity is maintained at the same level for all three time sections.

However, for short segments, the trapezium reduces to a triangle. In addition, the travel time for each section is sampled using the controller "clock time" and is rounded up to a greater value.

Therefore, for very short distances, the basis of the triangle (time) is fixed while its height (velocity) is adjusted, to ensure the square be equal to a specified displacement. Thus, the joining of very short segments yields unexpected velocity reduction, which is not admissible for the considered technological equipment.

To eliminate this effect, the distance between the successive nodes must be not less than

min min v S 9 =
, where v is the cutting speed, and 9 min is the minimal duration of the acceleration/deceleration sections. In practice, the value of min S is high enough to be considered during robot programming. For example, for typical industrial case studies (cutting speed 10 m/min; clock time 0.005 ms; at least 4 time-samples per acceleration section) the minimal distance between the nodes is equal to 3.3 mm. At the same time, in the model (1) the cutting contour is usually sampled with a step less then 1 mm.

To comply with this constraint on the inter-node distance, it is necessary to aggregate short segments into straight and circular portions, which are suitable for on-line implementations. This and the middle index

k m is recomputed as )/2 ( e s m k k k + =
, where . denotes the integer part. Then, each current path segment is verified for its capability to the line and circular approximation:

(1) if the angle between the vectors ) ( ) (

s m k - k p p and ) ( ) ( m e k - k p p is less than max ,
then the segment admits a straight-line approximation, and the index k e is increased.

(2) if the radius of the circle defined by the points ) (

s k p , ) ( m k p and ) ( e k p is less than R max ,
then the segment admits a circular approximation, and the index k e is increased.

(3) otherwise, the segment enlargement is stopped and the last accepted triple ) ( , etc. However, before applying this algorithm, the node sequence is divided into the sections with the same configuration index µ , and the subsequent sections are connected by the arcs using joint-space interpolation (command PTP). The length of the connecting segment must be also greater or equal to min S , so this segment is successively expanded in both directions.

For each segment, the motion arguments (tool locations in the Cartesian or joint space) are defined by the node L(k e ) for the LIN/PTP commands, and by the pair {L(k e ), L(k m )} for the CIR command. It should be noted, that at this stage there is no need to modify the joint angles for the 4th and 6th axes, which admit unlimited rotations (this procedure is executed by the controller on-line).

As a result, for the specified tolerances max and R max , the optimum path with the uniform short sampling step S is aggregated into longer sections, which are described by the robot control

Efficiency of the developed technique

Case-study: robot motion planning for 2D cutting

To demonstrate the efficiency of the proposed technique, first let us apply it to the planar cutting performed by a Scara robot, which possesses a 1 d.o.f. redundancy with respect to this task. The geometric/kinematic parameters of the manipulator are as follows: the link lengths are l 1 = l 2 = 1000 mm; the tool offset (horizontal) is d = 250 mm; the maximum speed of the axes involved in the planar motion is 120, 120, and 320 deg/sec, respectively. The cutting contour is a square 800 mm × 800 mm with rounded angles (radius 100 mm). Its centre is located at point (1000 mm, 1000 mm) and is surrounded by and an obstacle described in detail in [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF].

After the sampling, the contour is presented as a set of 60 uniformly distributed nodes (this corresponds to the cutting speed 30 m/min, sampling distance 50 mm, and sampling time 0.1 sec).

The direct kinematic model of this manipulator is described by the equations

3 2 1 2 1 2 1 1 2 1 2 1 1 ) ( ) ( ) ( ) ( ) ( ) (
q q q sin d q q sin l q sin l y cos d q q cos l q cos l x

+ + = + + + = + + + = : : : , (21)
where q 1 , q 2 , q 3 are the joint coordinates, and x, y, : are the tool location parameters (two Cartesian coordinates and an orientation angle). It should be noted that one of the joint coordinates (vertical translation of Scara manipulator) is not included in the model, since the z-axis motion is not redundant for the considered planar task.

To derive the inverse kinematic model, let us define ) ( 0 :

cos d x x = ; ) ( 0 : sin d y y =
and sequentially solve the obtained equations for q 2 , q 1 and q 3 : where µ = sgn(q 2 ) is the configuration index.

2 1 2 2 2 1 2 0 2 0 2 2 l l l l y x cos a q + = µ ; 2 2 0 2 2 1 0 2 2 0 2 2 1 0 1 ) ( ) ( 2 q sin l x q cos l l y q sin l y q cos l l x tan a q + + + = (22) 2 1 3 q q q = : (23)
Using the inverse model and altering the tool orientation : with the step of 10 o , a set of 1385 feasible tool locations {L k,j } and the corresponding set of joint coordinates {Q k,j } have been generated. To detect collisions, both the manipulator and the obstacle were described by a set of line segments. Then, all pairs from these sets were examined for intersections. The obtained manifold of the admissible motions is presented in figure 3, where the forbidden regions (with detected collisions) are shaded. Within the frames of this geometrical interpretation, the problem is to find the best closed path on the manifold surface, which encloses this figure and avoids the forbidden region. Another interpretation of the problem is given in figure 4, where the admissible and forbidden tool locations are shown in the set of the ( , t)-planes and the goal is to find the best route from the left to the right vertical lines. It should be noted that, in spite of the apparent simplicity of the last interpretation, the optimisation problem is essentially complicated by the non-Euclidian distances.

[Insert figures 3, 4 about here]

For the considered robotic task, a set of solutions was obtained. These solutions are presented in tables 1 & 2 and differ by both the optimisation criteria and the weights: w i , v and a . For each partial criterion, there are given the best values of single-objective optimisations (in brackets). As the tables show, for the relaxed velocity/acceleration constraints, tuning of the weights w i can barely produce acceptable results. Most of the solutions have a tendency for oscillations in the orientation axis, which is undesirable in practice.

In contrast, the explicit imposing of constraints on the velocity and acceleration (and adjusting corresponding weights v , a ) simplifies obtaining smooth solutions with the balanced values of the partial criteria. This result is illustrated in figures 5 and 6, which contain the trajectories generated for the weights 1 -

) ( t q w max i i = &
. They also demonstrate rather high oscillation of q 3 for the unconstrained case (see figure 5; joint velocity plots) and very smooth trajectories for the constrained optimisation (similar plots in figure 6). Hence, the proposed approach gives the possibility to eliminate oscillations in the orientation axes.

[Insert Tables 1, 2 about here]

[Insert figures 5, 6 about here]

The final comparison of the obtained results can be done using a multi-axes plot which presents two typical solutions with the values of nine partial design objectives

3 2 1 J J , J v v

K

(see figure 7).

The first solution (a) corresponds to the best result of a known technique, which corresponds to the very tedious adjusting of the weights w i . The second solution (b) was generated straightforwardly, using the developed method with the standard values 1 -

) ( t q w max i i = &
and with a very fast and simple tuning of the algorithm parameters v and a . As follows from the comparison, the developed technique reduces essentially oscillations of the orientation axis (see axes J 3s , J 3v and J 3 in figure 7, where relevant performance measures are lower by the factors 2…3), while maintaining almost the same values for the rest of the objectives.

[Insert figure 7 about here]

Comparison with other techniques

The main advantage of the proposed technique compared to the previous work [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF] is the ability to generate a smooth path without the tedious adjusting of weights w i . This useful feature is achieved by modifying the standard dynamic programming routine: including an implicit verification of the velocity/acceleration constraints while creating the backtracking matrix.

From this perspective, the new technique is robust with respect to weights assigned by the designer.

This advantage has been confirmed by a number of industrial projects where the developed technique has been successfully applied.

Another advantage, which was demonstrated by several 3D cutting case-studies, relates to the algorithm's ability to generate a robot path with large rotations of the manipulator orientation axes (more than 360°). This allows the exploitation of recent advances in the kinematics of cutting The developed technique has also advantages compared to the smooth path generation method recently proposed by [START_REF] Gasparetto | A new method for smooth trajectory planning of robot manipulators[END_REF]. In their paper, the objective function contains a term proportional to the integral of the squared jerk (defined as the derivative of the acceleration)

along the trajectory, while the second term is proportional to the total execution time. However, this formulation does not take into account the particularities of the considered application problem (where the execution time is fixed) and may produce local oscillations in the neighbourhood of the singularities caused by the integral-type performance measure used.

Hence, the developed technique yields a significant reduction of the robot programming time while respecting simultaneously both the smoothness objectives and the actuators velocity/acceleration constraints. In the future, this approach will be incorporated in more general design algorithms: (i) optimal path generation for several cutting contours, and also (ii) automatic robot placement in the cutting manufacturing cell. Another prospective application is the redundancy resolution for multi-manipulator robotic systems with path/velocity/acceleration constraints defined in Cartesian space. This corresponds to current industrial needs, since some recent designs for cutting cells include both an industrial robot (cutting tool manipulator) and a positioning table (workpiece manipulator).

Industrial implementation

The developed algorithms have been successfully implemented on the manufacturing floor, in ROBOMAX CAD package, which is a powerful system for robotic computer-aided design and offline programming. It is already used in automotive industry and has been successfully applied for the design of a number of manufacturing lines. The package is completely integrated with Autodesk TM products and includes some auxiliary tools, which enable the user to design a robotic cell and generate a control program taking into account the particularities of the employed technology.

With respect to laser cutting, the Robomax/Laser subsystems enable to design a workcell layout and optimise robot motion using multi-objective optimisation techniques. At the beginning, using standard routines of the Autodesk Mechanical Desktop (AMD), a mathematical description of the cutting contour is presented in the form of the "augmented line". The designer may define either the desired distance between vertices or their total number. In addition, he/she can estimate the minimum number of vertices required keeping the accuracy within tolerances. Relevant software tools allow also to aggregate separate segments in a common cutting contour and perform the unification of sampling distances.

The main design procedure consists of three iteratively repeated steps. The first is the selection of a proper manufacturing environment and location of the manufacturing task within the robot workspace. Here, the designer can create and optimise the layout of the cell interactively, using functions incorporated in the system. The required components (robotic manipulators, workpiece positioners, clamping devices, protective fences, columns, etc.) may be selected from the library. To optimise the position of the workpiece, a set of performance measures is used that considers such important features as the distance to the obstacles, closeness to joint limits, required range of joint motions.

The second step deals with path planning using the algorithms described in the previous sections. The goal is to utilize the kinematics redundancy of the robotic manipulator with respect to the cutting operations. To balance competing objectives, the designer can interactively assign priorities or weights, as well as activate/ deactivate constraints and criteria. In particular, it is possible to choose between the displacement/increment minimisation, define the distance metrics, and to change weights after visual analysis of the obtained path. In addition, there is an interactive editor allowing the adjustment of some nodes in accordance with other preferences, which are not presented in the formal criteria used in the model.

In the third step, the obtained solution is verified using realistic 3D simulations of the manipulator motions in the manufacturing environment. The ROBOMAX interactive debugger enables the user to evaluate both the total path and its separate segments, corresponding to elementary motions between the nodes. Also, it is possible to inspect joint, velocity and acceleration profiles for all manipulator axes. If the designer is not satisfied with the current solution, he/she can return to the first or the second step to change the design parameters to get a new solution.

After finishing this optimisation cycle, the generated path is converted into a control code in the robot programming language. At this stage, it is also possible to optimise the program, in order to eliminate some nodes (that belong to straight-line or circular segments), or to add additional commands required for control of technological parameters. At the end, the created program is uploaded to the robotic cell controller.

A recent application of the ROBOMAX/Laser is the off-line programming of a robotic cutting station for a minivan production based on KUKA robots and relevant positioning and clamping devices. The station is used for small-batch manufacturing, which requires frequent reprogramming.

As follows from our industrial experience, the developed technique provides adequate capabilities for generating robot control codes and crucially reduces the programming and changeover time.

Conclusion

The optimal path generation for the redundant manipulators is an important issue for off-line programming of the laser cutting robotic systems. Recent advances in laser technology, and especially essential increase of the cutting speeds, motivate amending the existing methods. At present, these do not completely utilise the actuator capabilities as well as neglect some particularities in the mechanical design of the manipulator wrist.

This work contributes to a multi-objective optimisation of the manipulator motions for the continuous contour tracking with an axis-symmetric technological tool. In contrast to previous work, this optimisation technique explicitly incorporates the verification of the velocity/acceleration constrains, enabling the designer to interactively define their importance with respect to the pathsmoothness objectives. In addition, the proposed approach takes into account the capacity of some wrist axes for unlimited rotation in order to produce more efficient motions.

The developed path planning algorithm is based on dynamic programming with some extensions, which incorporate the constraint checking for each segment of a candidate solution and the penalty assignment if some constraints are violated. To generate a smooth motion, each joint trajectory is evaluated by a set of performance indices such as the coordinate deviation, maximum increment, and total displacement. The search space is converted into a directed acyclic graph and the problem is re-formulated in terms of combinatorial optimisation theory. To evaluate the distance between the successive tool locations, three types of metrics are used, each of which has the weight assigning option for the corresponding joint axis. During the optimisation, the vector objective is converted into a scalar one using the sum or minimax approaches. Then the weights are altered to generate a set of Pareto-optimal solutions. The proposed algorithm has been implemented in an industrial software package and verified on the manufacturing floor in the automotive industry. Its application for off-line programming of robotic laser cutting stations yielded significant reduction of the process planning and changeover time.

Future research will focus on the development of technology-oriented algorithms for optimal robot path planning for several cutting contours. This problem arises from a number of industrial applications where the workpiece contains numerous contours to be cut by a single robot or by several robots. The corresponding path-planning technique requires optimal scheduling of the cutting contours, load balancing among robots, collision avoidance during simultaneous operation, etc. Another prospective research direction is in the automation of the robot placement for cutting manufacturing cells, which requires optimisation of the robot location with respect to the workpiece and also the integration of this technique with the capabilities of the relevant robotic CAD systems. Manipulator motion planning for high-speed robotic laser cutting (Dolgui, Pashkevich) 37

F

Joint coordinates Joint velocities q 1 q 2 q 3 q 1 . q 2 . q 3 .

Joint coordinates

Joint velocities q 1 q 2 q 3 q 1 . q 2 . q 3 .

  algorithms). Concerning the joint limits, it worth mentioning that the latest laser-cutting manipulators allow unlimited rotation of the 4 th and 6 th axes, so in this locations are inevitably excluded from a feasible set. The collision detection functions are standard routines of industrial robotic CAD packages, together with the direct/inverse kinematics of the robotic manipulators.

  ith joint velocity and acceleration; and v and a are the scaling factors to be adjusted by the designer. Hence, the complete set of constraints arising from the technical nature of the problem are summarised in inequality ! ) ( k kQand in the expressions (6), (7), which ensure the path existence and its admissible curvature in the joint coordinate space.

  desired solution for n+1 clusters S 0 , S 1 , … S n can be obtained sequentially, starting from k = 0 with j the k-index up to n, computing

8 )

 8 If d cur < d min then Set d min := d cur ; j opt := p (c) Set M dist (k, j) := d min ; M ind (k,j) := j opt ; (3) Set d min := inf; (4) For j := 1 to J max (n) do (a) Set d cur := M dist (n, j) (b) If d cur < d min then Set d min := d cur ; j opt := j (5) For k = 0 to n do Set J opt (n-k) := j opt ; j opt := M ind (n-k, j opt );

  using our heuristic algorithm, which is based on the sequential extracting of maximum-length sub-sequences with a fixed initial point that can be approximated by a straight or circular line within given tolerances.At each run, the algorithm deals with three nodes: p(k s ), p(k m ) and p(k e ) , where k s , k m and k e are the indices of the start, middle and end nodes of the current segment, respectively. Within the search loop, the index k e is successively increased from the initial value 2 + s k

  of the linear or circular type (commands LIN and CIR).This loop is then repeated for the rest of the nodes, assigning new values for the

  , CIR and PTP , with required location arguments) and may be reproduced on-line by a robot controller with the desired Cartesian velocity.

  as the endless rotation of wrist joints. Besides, because the specific generation of the search space (see figure2) and relevant distance computing for the angles, the feasible solutions may successfully pass through singularities. This allows the changing of the configurations indices while in motion.

  Figure 1. Definition of the task frames

Figure 5 .

 5 Figure 5. Optimisation without the velocity/acceleration constraints

Figure 6 .Figure 7 .

 67 Figure 6. Optimisation with the velocity/acceleration constraints

  Under such assumptions, the feasible search space can be represented by a directed graph with the vertices

	F o r
	P
	e
	e r
	R e v
	i e
	w
	O n l
	y
	about here]

  each cluster. It should be noted that the above formulation does not admit cluster intersection (i.e.

	S k	I	S l	=	,	4	k	!	l
					F o r
									P
									e
										e r
										R e v
										i e
										w
										O n l
										y
										Here	v kj =	( j , k	)	denotes the jth vertex of the kth cluster.
	It is assumed that each edge	v u ) , (	E	is evaluated by some travelling distance	2	) ( v , u	. The clusters
	must be visited in a strictly predefined order	S 0	.	S 1	.	K	.	n S	, passing via exactly one vertex in

Table 1 .

 1 Performance measures for the design problem

	o r
	P
	e
	e r
	R e v
	i e
	w
	O n l
	y

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

Acknowledgements. The second author acknowledges the Ecole des Mines de Saint Etienne where he was visiting professor. The authors thank also Chris Yukna for his help in English.


Satchin Soorkia

Frederic Le Quere
email: lequere@univ-mlv.fr

Céline Léonard

Detlev

Soorkia

F Le Quéré

D Figgen

Dr Frédéric Le Quere

Ab initio study of the spin-orbit coupling between the A 1 Σ + u

Keywords: spin orbit coupling, Wave packet propagation, Dipole moment

come

Ab initio study of the spin-orbit coupling between the A1Σ +u and b 3Πu electronic states of Na2

Electronic structure calculations have reached such a level of precision that a totally ab initio treatment can be envisioned to simulate experiments involving several electronic states that may be coupled. The increasing use of laser pulses in chemical physics [1,2,3], in our case to induce transition between electronic states, requires a complete wave packet dynamical treatment of the system in order to fully describe the simulated process. The couplings between electronic states, such as spin-orbit interaction, are often very weak but play an important role in a lot of experiments [4].

Extracting informations on these couplings via analysis of experimental data generally requires a huge amount of work [5,6] while they can be readily obtained from theory.

In this paper, we present such an ab initio "theoretical experiment" that will try to reproduce experimental data from Rutz et al. [4] involving three electronic states of Na 2 . The ground state, X 1 Σ + g , is coupled to the A 1 Σ + u state by a laser pulse interacting with the transition dipole moment between these two states, and the A 1 Σ + u state is coupled to the b 3 Π u state by spin-orbit interaction. The full description of the experiment requires a precise determination of the electronic states, of the transition dipole moment and of the spin-orbit coupling, computed from electronic structure theory and a time dependent simulation of the evolution of the nuclear wave packet from which the observables of interest will be computed. Na 2 has been the subject of numerous experimental and theoretical studies: the electronic potential energy functions have been derived from ab initio studies by several authors [7,8] and Rydberg-Klein-Rees (RKR) type potentials have been built from experimental data [9,10,11]. Transition dipole moments are also available from theoretical [12] and experimental [13] sources. The interaction between the A 1 Σ + u and b 3 Π u states due to the spin-orbit coupling has been deeply studied. Stolyarov et al. [14] observed this interaction in g factors. Dulieu and Julienne [15] predicted and Effantin et al. [5] observed perturbations in the spectra of these coupled states, however we did not find in the literature any coordinate dependent spin-orbit coupling function.

All these references give accurate information concerning our "theoretical experiment" but we chose to compute all the needed electronic structure data in order to ensure the global coherence of our calculations and to determine the complete spin-orbit coupling function. However, we compared our results to all available previous calculation.

For the nuclear motion study, the time dependent wave packet is obtained from the Schrödinger Where the time dependent hamiltonian, Ĥ is the sum

Ĥ = Ĥe + ĤSO + T -ˆ µ. E(t)
In the second section of this article, we will describe the electronic structure calculation and derive representations for the scalar-relativistic electronic hamiltonian, Ĥe , the spin-orbit interaction term ĤSO and the transition dipole moment operator ˆ µ. In the third section we will compute bound vibrational states of these electronic states. They will be used in the last section where the full time dependent Schrödinger equation will be solved to simulate the experiment.

Electronic structure

The electronic energies were computed using the MOLPRO package [16] for internuclear distances in the range of 3 < r < 20 bohr. The 10 core electrons of each Na atom are replaced by the effective core potential (ECP) developed by Fuentealba et al. [17]. The valence electron is represented by 8 s, 6 p and 1 d even tempered basis functions optimized with the ECP. The core-valence correlation is added by the corresponding core polarisation potential (CPP). The energy of each electronic state (X 1 Σ + g , A 1 Σ + u , b 3 Π u ) and the transition dipole moment for the transition X 1 Σ + g -A 1 Σ + u are determined using the internally contracted MRCI [18,19] where the reference comes from a MCSCF [20,21] calculation with a molecular active space formed by all the valence plus the 4s and 4p orbitals of the two Na atoms. In this MCSCF step, the two electronic components of b 3 Π u were averaged together, as well as the X 1 Σ + g and A 1 Σ + u electronic states.

The spin-orbit coupling terms are the matrix elements of the spin-orbit part of the pseudopotentials, ĤSO , evaluated in the basis set formed by the MRCI electronic wavefunctions of A 1 Σ + u and b 3 Π u . As references, the MRCI calculations used the MCSCF wavefunctions determined with an active space built with the valence orbitals. For these calculations, the same ECP as above is used without the CPP part but with inclusion of an accompanying spin-orbit potential that was adjusted.

Since the one-component ECP of Fuentealba et al. was adjusted to experimental energy data, the terms of the spin-orbit ECP were also adjusted to experimental spin-orbit splittings. We used the spin-orbit splittings of the 3p and 3d configuration [START_REF] Ch | Atomic Energy Levels[END_REF] as reference data for the adjustment of the p-and d-term, respectively; both adjustments were done seperately. The spin-orbit terms were adjusted in two-component valence average-level calculations with a formally non-relativistic Hamiltonian (including the ECP) for these configurations; this was done in numerical calculations using the GRASP code [START_REF] Dyall | GRASP: atomic numerical MCDHF package[END_REF]. By this method we generate two-component pseudopotentials V l,l+1/2 and V l,l-1/2 from which the spin-orbit potential V SO can easily be obtained:

V SO = 2 l=1 l 2l + 1 (V l,l+1/2 -V l,l-1/2 )P l l • ŝP l
P l is the projetion operator on the Hilbert subspace of angular symmetry l:

P l = l m l =-l |lm l lm l |
In accordance with the one-component ECP that consists of simple Gaussians, the terms of the spin-orbit potentials ∆V l = V l,l+1/2 -V l,l-1/2 were also chosen to be simple Gaussians; moreover, we fixed the exponents of the spin-orbit terms to the exponents of the respective one-component term and solely optimized the coefficients of the Gaussian functions. Table 1 shows the exponents, B and coefficients, β of both the one-component and the spin-orbit potential.

============================== Table 1 near here

==============================

We consider the following spin-orbit matrix elements

LS x = 1 Σ + u , M S = 0| ĤSO | 3 Π uy , M S = 1 LS y = 1 Σ + u , M S = 0| ĤSO | 3 Π ux , M S = 1 LS z = 3 Π ux , M S = 1| ĤSO | 3 Π uy , M S = 1 with ĤSO = 2 λ=1 V SO λ
where λ denotes one of the two Na atoms. LS x and LS z are pure imaginary terms and LS x = iLS y , the z axis being the internuclear axis.

The electronic energies, the transition dipole moment between X 1 Σ + g and A 1 Σ + u , and the spin-orbit coupling elements LS y and -iLS z were fitted by 14th order polynomial expansions of morse-type coordinate,

Q = 1 -exp[-0.3(r -r e )]
the 0.3 exponent being found suitable for all the functions. r e is the equilibrium internuclear distance of the given electronic state. For the coupling functions, r e is chosen to be the equilibrium internuclear distance of the electronic ground state X 1 Σ + .

The coefficients of the potential fits are given in Table 2. Note that the first coefficient gives the 2 a "shifted" coefficient that was used in all the bound states and dynamical calculations in order to fit the correct experimental transition energies (see below in Table 3). Figure 1 shows the energy of the three electronic states as a function of internuclear distance. The crossing of the A 1 Σ + u and b 3 Π u states occurs at r = 7.02 bohr.

==============================

Table 2 and Figure 1 near here ============================== Equilibrium geometries, r e , of the three electronic functions are given in Table 3, with other spectroscopic data. In order to compare our results to previous RKR potentials, we also derived the harmonic wavenumbers, ω e , and the anharmonic constants, ω e x e [START_REF] Papousek | Molecular vibrational-rotational spectra[END_REF], the equilibrium rotational constants, B e , the dissociation energies D e and the energy differences between the ground state and each excited state equilibrium positions, T e .

For comparison, non-relativistic CCSD(T) [START_REF] Hampel | [END_REF]26] calculations on the electronic ground state correlating the sub-shells 2s, 2p and 3s of each Na atoms have been carried out using the cc-pCVQZ [27] basis set. The corresponding values of r e = 5.835 bohr and ω e = 158.3 cm -1 are in nice agreement with the ECP + CPP / MRCI results which validates the accuracy of the present ab initio ECP + CPP / MRCI potentials. For these potentials, Table 3 shows that the equilibrium distances are always smaller by 0.01-0.02 bohr than the experimental values, the harmonic wavenumbers larger by 1-2 cm -1 and the energies D e are always smaller by 80-200 cm -1 . The discrepancies are stronger for the spectroscopic constants of the A 1 Σ + u state due the flatness of the potential. The equilibrium transition energies, T e , are 13 cm -1 below the experimental value for b 3 Π u and 78 cm -1 for A 1 Σ + u . The error in T e value is different for each state because the correlation energy treatment is different for single and triplet electronic states in MRCI. However, our results are in better agreement with experimental data than previous accurate ab initio data obtained by Magnier et al. [8] using configuration interaction calculations and pseudopotentials. The energy difference between the X 1 Σ + g and b 3 Π u electronic states at dissociation gives 16872 cm -1 which is in good agreement with the experimental value of 16956 cm -1 determined by Moore [START_REF] Ch | Atomic Energy Levels[END_REF] for the energy difference between X 2 S and A 2 P electronic states of Na. The coefficients for the fits of the coupling functions are given in Table 4. The transition dipole moment between X 1 Σ + g and A 1 Σ + u is plotted for r between 3 and 20 bohr in Figure 2 and is compared with other calculations [12,13]. The three curves present a similar behavior for values of r up to r e and our function is very close to recent ab initio data from Ahmed et al. [13] who calculated the transition dipole moment with the method used by Magnier et al. [8] for the electronic states. For larger values of r the theoretical data of Konowalow et al. [12] obtained at the MCSCF level does not tend towards the experimental atomic limit [12] of -3.56 ± 0.11 a.u. in contrary to the other ab initio functions. The transition dipole moment displays an almost linear behavior in the Franck Condon region of the ground state, around 6 bohr, but this approximation will not remain valid if one chooses to study transitions from an excited vibrational state of the ground electronic state that reaches broader part of the configuration space. ============================== Figure 2 and Table 4 near here ============================== The spin-orbit coupling elements between A 1 Σ + u and b 3 Π u are shown in Figure 3. Their variation in the crossing region of these two electronic states is very small and a linear approximation might again give reasonable results although we used the computed couplings in the calculations of the next sections. Up to now, no geometry dependent coupling was published. However, Dulieu and Julienne [15] used a variable function for the fine structure coupling 2LS y which varies from √ 2/3∆E F S = 8.11 cm -1 for r → ∞ to 5.47 cm -1 for r = 8.5 bohr, where ∆E F S = 17.1963 cm -1

==============================

is the atomic fine structure splitting of the first 2 P state of Na [START_REF] Ch | Atomic Energy Levels[END_REF]. Our ab initio calculations give ∆E F S = 17.21 cm -1 and 2LS y = 8.1 cm -1 at r → ∞, 2LS y = 5.48 cm -1 at r = 8.5 bohr and 2LS y = 6.01 cm -1 at r = 7.02 bohr (the crossing point of the potential electronic curves of A 1 Σ + u and b 3 Π u ). These values are in remarkable agreement due to the fact that the parameters of the pseudopotential (see Table 1) were fitted in orderto reproduce as closely as possible the atomic spin-orbit splitting. These values can also be compared to experimental values extracted from vibrational analysis giving constant couplings of 5.91±0.14 cm -1 , 5.769±0.006 cm -1 and 5.97±0.15 cm -1 [5,14,28]. In the basis of 1 Σ + u , M S = 0 ; 3 Π ux , M S = 1, 0, -1 ; 3 Π uy , M S = 1, 0, -1 , the sum of the scalarrelativistic electronic Hamiltonian, Ĥel and the spin-orbit coupling operator, ĤSO is associated to the matrix :

==============================

H el + H SO =           V ( 1 Σ + u ) LS y 0 LS y LS x 0 -LS x LS y V ( 3 Π u ) 0 0 LS z 0 0 0 0 V ( 3 Π u ) 0 0 0 0 LS y 0 0 V ( 3 Π u ) 0 0 -LS z -LS x -LS z 0 0 V ( 3 Π u ) 0 0 0 0 0 0 0 V ( 3 Π u ) 0 LS x 0 0 LS z 0 0 V ( 3 Π u )          
The eigenfunctions of Lz , labeled |M L , M S ,

|0, 0 = | 1 Σ + u , M S = 0 | ± 1, j = 1 √ 2 | 3 Π ux , M S = j ± i| 3 Π uy , M S = j with j = 1, 0 or -1, are introduced. In the basis {|0, 0 , |1, -1 , |-1, 1 , |-1, 0 , |1, 0 , |-1, -1 , |1, 1 },
and using the fact that LS x = iLS y , the above matrix becomes:

H ′ el +H ′ SO =           V ( 1 Σ + u ) √ 2LS y √ 2LS y 0 0 0 0 √ 2LS y V ( 3 Π u ) -iLS z 0 0 0 0 0 √ 2LS y 0 V ( 3 Π u ) -iLS z 0 0 0 0 0 0 0 V ( 3 Π u ) 0 0 0 0 0 0 0 V ( 3 Π u ) 0 0 0 0 0 0 0 V ( 3 Π u ) + iLS z 0 0 0 0 0 0 0 V ( 3 Π u ) + iLS z          
This matrix is now block-diagonal and highlights the fact that the b 3 Π u state splits into 6 components labeled Ω = 2, 1, 0, 0, -

, 0 , 1 √ 2 [|1, -1 + | -1, 1 ] , 1 √ 2 [|1, -1 -| -1, 1 ]}: H ′′ el + H ′′ SO =   V ( 1 Σ + u ) 2LS y 0 2LS y V ( 3 Π u ) -iLS z 0 0 0 V ( 3 Π u ) -iLS z  
This matrix confirms that the only spin electronic component of b 3 Π u which interacts with A 1 Σ + u via spin-orbit coupling is Ω = 0 + corresponding to This interaction creates an avoided crossing between the two states. In order to reduce the number of spin-electronic states involved in the dynamics presented in the next sections, the preceeding matrix is reduced into a two dimensional space containing only the Ω = 0 + component of b 3 Π u and the A 1 Σ + u state. The reduced matrix in the basis of 1 Σ + u (Ω = 0 + ), 3 Π u (Ω = 0 + ) is:

| 3 Π u (Ω = 0 + ) = 1 2 -| 3 Π ux , M S = 1 -| 3 Π ux , M S = -1 + i| 3 Π uy , M S = 1 -i| 3 Π uy , M S = -1
H el + H SO = V ( 1 Σ + u ) 2LS y 2LS y V ( 3 Π u ) -iLS z
This diabatic representation of two coupled states was used in the following wave packet study.

Vibrational states determination

We computed bound vibrational states on the ground 1 Σ + g , the A 1 Σ + u , the b 3 Π u and the b 3 Π u (Ω = 0 + ) states using a Lanczos diagonalisation scheme [START_REF] Cullum | Lanczos algorithm for large symmetric eigenvalue computation[END_REF] on a grid of 512 points between 2.5 and 20 bohr.

Experimental energies of the A 1 Σ + u bound states were compiled by Gerber and Möller [10] and used to create a RKR potential for this electronic state, and a RKR potential for the b 3 Π u state was established by Whang et al. [11]. Table 5 shows the energies determined from our potentials and the RKR energies. The RKR potentials were made from experimental data associated to highly excited rovibrational states, and spin-orbit deperturbated in order to get uncoupled A 1 Σ + u and b 3 Π u states. There is a nice agreement between experimental and theoretical values although the difference tends to slowly increase up to 22 cm -1 at v = 20 for the singlet and up to 11 cm -1 at v = 20 for the triplet showing that our potentials are slightly less anharmonic than the RKR potentials. The energy splitting, A v , in the normal spin-orbit multiplet of a b 3 Π u vibrational state before interaction with the A 1 Σ + u state has been derived by Shimizu and Shimizu [6] in a perturbative expansion as a function of v :

A v = A e + α A (v + 1 2 )
where A e = 7.08 ± 0.06 cm -1 and α A = -0.0156 ± 0.0034 cm -1 .

In our case A v can be evaluated from the difference between the numbers in the columns associated with b 3 Π u and b 3 Π u (Ω = 0 + ) of Table 5 since the b 3 Π u is associated with the unperturbed terms b 3 Π u (Ω = ±1). A linear regression of this difference as a function of (v + 1/2) gave values of A e = 7.13 cm -1 and α A = -0.0171 cm -1 in agreement with Shimizu and Shimizu [6] numbers. This nice agreement shows that our LS z function is of good quality in a large part of configuration space, spanned by highly excited vibrational states. The spin-orbit coupling between A 1 Σ + u and b 3 Π u has been experimentally studied with great accuracy by Effantin et al. [5] who obtained state to state spin-orbit interaction parameters, ξ vv ′ , defined by

==============================

ξ vv ′ = 1 Σ + u , v|2LS y | 3 Π u (Ω = 0 + ), v ′ .
Using eigenstates obtained in the previous section and the spin-orbit coupling function of Table 4, we computed ξ vv ′ for v values from v = 0 to v = 17 and v ′ values from v ′ = 0 to v ′ = 50. The largest value is ξ 0,5 = 2.67 cm -1 but it is interesting to note that ξ vv ′ remains non negligible for a broad set of v and v ′ numbers and that the general coupling scheme does not seem to be reduced to the coupling of a particular A 1 Σ + u vibrational state to a particular b 3 Π u vibrational state. For example ξ 0,v ′ > 10 -2 cm -1 for v ′ values up to v ′ = 16 and ξ 17,v ′ > 10 -2 cm -1 for all v ′ < 40. Table 6 shows a selection of ξ vv ′ compared to those obtained by Effantin et

Time dependent analysis

In order to further check the validity of our electronic calculations, we tried to simulate the time resolved experiment of Rutz et al. [4]. A pump laser pulse excites the molecule from the ground electronic state of Na 2 to the A 1 Σ + u state and a probe pulse identical to the pump ionizes the excited molecule to the ground state of Na + 2 via a two photon process involving the 1 Π g state. The ion signal is measured as a function of the delay between the pump and probe pulses and an oscillatory pattern is obtained, reflecting the oscillating movement of the wave packet on the A 1 Σ + u surface. An analysis of this signal gives information on the position of the vibrational energy levels of the system, and the influence of the spin-orbit coupling on the level positions is detected when the pulse carrier frequency is set to λ = 642 nm but is not seen for λ = 620 nm. These values of excitation energies are associated with resonances of the ionization step when the wave packet is at the outer or inner classical turning point of the A 1 Σ + u state, respectively. We used the Short Iterative Lanczos (SIL) [START_REF] Park | [END_REF] algorithm to propagate the wave packet on the three coupled states. The grid is the same as the one used in bound states determination above and the initial wave packet is the ground vibrational state of the ground electronic state, a reasonable guess since the vibrational temperature in the experiment is between 10 K and 50 K. The pump pulse is the same one as in the experiment and we carried out the evolution for times up to 3 ps. The experimental pump-probe pulse delays of ref [4] varies up to 60 ps, but the few oscillations that we describe with our propagation are sufficient to discuss most of the features observed. The probing part of the experiment was not simulated but the wave packet evolving on the A 1 Σ + u contains all the relevant information.

The theoretical oscillation period was computed from the oscillations of the expectation value of the bond length as a function of time. This value is equal to the measured oscillation period, especially during the first oscillations, when the wave packet spreading has not taken place. We used a pulse identical to the experimental one with a sech 2 envelope (FWHM= 110 fs, peak power 0.5 GW cm -2 ) and a given wavelength, λ. The complete pulse has the form

E(t) = A sech 2 [α(t -t 0 )] cos(2πt/λ)
Table 7 summarizes the pulse parameters used in our calculations and the theoretical and experimental oscillation periods. A good agreement is obtained, indicating that the first step of the process, i.e. the transition from the ground electronic state to the A 1 Σ + u state is well described in our simulation.

==============================

Table 7 near here ============================== The next step is the analysis of the populations of the vibrational levels of the A 1 Σ + u and b 3 Π u (Ω = 0 + ) states, done by projecting the time dependent wave packet |ψ(t) on the previously computed

eigenstates |ψ v |c v (t)| 2 = | ψ(t)|ψ v | 2
In the λ = 620 nm experiment, the most populated states are in the v = 12 to v = 14 range in the A 1 Σ + u electronic state and no effect of the spin-orbit coupling was detected. This is in agreement with our calculations (see figure 4). These levels are populated during the pulse duration (roughly 250 fs) and the population of all the triplet vibrational states are below 5 • 10 -4 , even after 3 ps when the wave packet has reached the coupling region several times in its oscillating movement.

In the λ = 642 nm experiment, the most populated states are in the v = 7 to v = 9 range and small deviations in the position of the lines were attributed to the spin-orbit coupling with triplet vibrational states around v ′ = 14. The plot of vibrational states population of the singlet (figure 5)

shows an important increase of the population of v = 7 to v = 9 states during the pulse excitation The first order perturbation term of two states coupled via spin-orbit interaction is

P vv ′ = ξ 2 vv ′ E v -E v ′ .
From Table 5, we can see that the energy difference between v = 13 and v ′ = 18 is less than 3 cm -1

and that the energy difference between v = 8 and v ′ = 14 is less than 6 cm -1 , indicating in both case a possibility of resonance. However, the ξ vv ′ coupling element from Table 6 is an order of magnitude stronger in the case of the v = 8 v ′ = 14 interaction which explains the larger population of this state in the λ = 642 nm experiment.

============================== Figure 4 and figure 5 near here

==============================

The largest experimental energy shift [4] of a vibrational level on the A 1 Σ + u due to spin-orbit coupling was observed for v = 8. The value of the shift was evaluated at δ 8 = 0.8 cm -1 . To obtain this shift, we must add the contribution of all the v ′ states to the perturbation of a given v state.

δ v = v ′ P vv ′
As noted above, states which are not in close energetic resonance may still contribute to the global shift of the level because the ξ vv ′ coupling terms remain non negligible for a large set of v and v ′ values. We obtain a global shift of δ 8 = 0.3 cm -1 . A probable reason of the difference between this value and the experimental one is molecular rotation that was not taken into account in our dynamical approach. Experimental data shows that rovibrational states up to J=25 are involved in the process and the summation in the preceeding equation should be on rovibrational states rather than only vibrational states.

Conclusion

The ab initio potential, transition dipole moment and spin-orbit coupling for the electronic states were computed with high accuracy and form a coherent package of data that can be used to simulate dynamical processes involving these states. We were able to reproduce experimental results from Rutz et al [4] concerning the oscillating period of the wave packet on the A 1 Σ + u state for two different excitation frequencies from the ground electronic state and the pattern of coupling between [13] and Konowalow et al. [12]. 3. The reference geometry, r e for each fit comes from Table 3.

Coefficients X 1 Σ + g A 1 Σ + u b 3 Π u C 0 0.
Σ + g , shifted A 1 Σ + u , shifted b 3 Π u and b 3 Π u (Ω = 0 + ) (obtained from shifted b 3 Π u ) electronic states of Na 2 . Level v X 1 Σ + g A 1 Σ + u A 1 Σ + u [10] b 3 Π u b 3 Π u [11] b 3 Π u (Ω = 0 + ) 0 79

  of the functions and we included in Table

Figure 3

 3 Figure 3 near here ==============================

  -2 due to internal spin-orbit coupling. The |Ω| = 2 components are associated with | -1, -1 , |1, 1 and the |Ω| = 1 with | -1, 0 , |1, 0 respectively. The spin electronic components of the b 3 Π u which interacts with A 1 Σ + u via spin-orbit coupling are: |1, -1 and | -1, 1 corresponding to Ω = 0. If the linear combinations of these components are used to express Ĥel + ĤSO , the following matrix is obtained in the basis {|0

  of the triplet state show an important increase in the population of v ′ = 14 state.

12 Page

 12 Figure captions:

Figure 1 :

 1 Figure 1: Variation of the potential energy of the X 1 Σ + g , A 1 Σ + u and b 3 Π u (Ω = 0 + ) electronic states of Na 2 with the internuclear distance.

Figure 2 :

 2 Figure 2: Variation of the transition dipole moment between the X 1 Σ + g and A 1 Σ + u electronic states of Na 2 with the internuclear distance r, compared to previous ab-initio data from by Ahmed et al.

Figure 3 :

 3 Figure 3: Variation of the spin-orbit coupling elements between the A 1 Σ + u and b 3 Π u electronic states of Na 2 with the internuclear distance.

Figure 4 :

 4 Figure 4: Population of the most populated vibrational levels on a) A 1 Σ + u and b) b 3 Π u after excitation at λ = 620 nm.

Figure 5 :

 5 Figure 5: Population of the most populated vibrational levels on a) A 1 Σ + u and b) b 3 Π u after excitation at λ = 642 nm.

Figure 1 20

 1 Figure 1

Figure 2 21

 2 Figure 2

Table 3 near here ============================== Page 5 of 24 URL: http://mc.manuscriptcentral.com/tandf/tmph

 3

	Molecular Physics

Table 5

 5

	near here
	==============================
	8

  al. The agreement is again very nice, showing the good quality of the LS y component of the spin orbit term, and of the calculated

	wavefunctions.
	==============================
	Table 6 near here
	==============================

Table 1 :

 1 Coefficients B and exponents β of the one-component ( V SDF ,[17]) and spin-orbit ( V SO )

	pseudopotential for Na.
	l	B	β
	V SDF
	0 10.839000 1.378000
	1	2.303000 0.663900
	2 -1.777000 0.924900
	V SO
	1	0.015867 0.663900
	2 -0.001224 0.924900

Table 2 :

 2 Polynomial expansion coefficients of the electronic energies (in hartree) of the X 1 Σ + g , A 1 Σ + u , b 3 Π u electronic states of Na 2 . The C ′ 0 value shifts the surfaces in agreement with experimental transition energies of Table

Table 3 :

 3 Equilibrium geometries (/bohr) and spectroscopic constants (/cm -1 ) of the X 1 Σ + g , A 1 Σ + u , b 3 Π u electronic states of Na 2 .

	State		r e	ω e	ω e x e	B e	D e	T e
	X 1 Σ + g	RKR a Theo. A b	5.818631* 159.17732 0.760159 0.1546855 6022.6 5.83 156.8 5892
		Theo. B b	5.77	159.3			5725
		ECP+CPP/MRCI	5.7995	159.893	0.7487	0.1557	5950
	A 1 Σ + u	RKR c Theo. A b	6.87474 6.86	117.2703 117.5	0.3534801 0.1108142	8310 14680.682 8284 14575
		Theo. B b	6.85	119.5			8118	14581
		ECP+CPP/MRCI	6.8593	119.147	0.3810	0.1113	8068	14602
	b 3 Π u	RKR d	5.87075	154.209	0.47682	0.15195	9475.3 13520.946
		Theo. A b	5.87	154.6			9411	13447
		Theo. B b	5.81	153.0			9396	13304
		ECP+CPP/MRCI	5.8561	154.946	0.4722	0.1527	9316	13507

a Ref.

[9]

, b Ref.

[8]

c Ref.

[10]

, d Ref.

[11]

.

* calculated from B e .

Table 4 :

 4 Polynomial expansion coefficients of the transition dipole moment between the X 1 Σ + g , A 1 Σ + u electronic states and of the spin-orbit coupling between the A 1 Σ + u and b 3 Π u electronic states of Na 2 .

	Coefficients Dipole moment (a.u.) LS y (/10 -5 ) (a.u.) -iLS z (/10 -5 ) (a.u.)
	C 0	-3.662650	1.635586	-2.286497
	C 1	-1.005127	-1.012301	2.009225
	C 2	-0.685894	0.350559	-0.908788
	C 3	-0.172353	0.748145	-0.006000
	C 4	1.694909	-0.649386	-0.344029
	C 5	4.232320	-4.163085	-2.639413
	C 6	-3.950725	3.259819	1.915624
	C 7	-12.639680	16.838510	9.289735
	C 8	11.693785	-4.351840	-3.316573
	C 9	29.699469	-31.730642	-18.524042
	C 10	-8.743500	-1.191377	-0.961422
	C 11	-33.603118	29.706639	16.776839
	C 12	-4.957655	7.978938	5.141925
	C 13	13.275650	-10.805054	-5.429050
	C 14	5.234212	-4.774218	-2.496550
	RMS error	0.00084	9.92×10 -9	4.11×10 -9

Table 5 :

 5 Vibrational states energies (/cm -1 ) of the X 1

Table 6 :

 6 State to state spin-orbit interaction parameters, ξ vv ′ (/cm -1 ).

	.52 14741.89 14739.24 13597.70 13597.92	13590.57
	1

Table 7 :

 7 Pulse parameters and oscillation period of the wave packet.

	fs)

λ (/nm)

A (/a.u.) α (/(fs -1 )) t 0 (/fs) T calc (/fs) T exp (/URL: http://mc.manuscriptcentral.com/tandf/tmph