Keywords: Redundant robotic system, Motion planning, Time-optimal trajectory, Dynamic programming, Automated composite lay-up iii

The thesis deals with the optimal motion planning in redundant robotic systems for automation of the composite lay-up processes. The primary goal is to improve the lay-up workcell productivity by developing a novel methodology of optimizing coordinated motions of the robotic manipulator, workpiece positioner and workspace extension unit, which ensure the shortest processing time and smooth movements of all mechanical components. In contrast to the previous works, the proposed methodology provides high computational efficiency and also takes into account both the technological constraints and the robotic system constraints, which describe capacities of the actuators and are expressed by the maximum allowable velocities and accelerations in the actuated joints. The developed technique is based on conversion of the original continuous problem into a combinatorial one, where all possible configurations of the mechanical components are represented as a directed multi-layer graph and the desired time-optimal motion is generated using dynamic programming principle for searching the shortest path on the graph satisfying the smoothness constraints. It is also proposed an enhancement of this technique by dividing the optimization procedure in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with large discretization step. Then, the same technique is applied in the local search space, which is created with smaller step in the neighborhood of the obtained trajectory. Alternatively, the second stage may implement a straightforward smoothing of the redundant variable profiles. The advantages of the developed methodology are confirmed by industrial implementation on the factory floor that deals with manufacturing of the high-pressure vessel.

GENERAL INTRODUCTION

Motivation. At present, composite materials are increasingly used in engineering practice. Compared to traditional ones, they have good strength-to-weight ratio, durability, flexibility of shaping and corrosion resistance [START_REF] Christensen | Mechanics of composite materials[END_REF][START_REF] Jones | Mechanics of composite materials[END_REF][START_REF] Kaw | Mechanics of composite materials[END_REF][START_REF] Lubin | Handbook of composites[END_REF][START_REF] Peters | Handbook of composites[END_REF]. For these reasons, they are extremely attractive for fabrication of large dimensional parts in aerospace, automotive and marine industries [START_REF] Gay | Composite materials: design and applications[END_REF][START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF].

A conventional way of manufacturing such composite components is based on the manual lay-up process. However, this labor-intensive procedure is quite slow, expensive and does not allow achieving required repeatability [START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF]). An alternative solution implementing such process is usually based on automated lay-up of the composite tape or fiber, which are placed layer-by-layer to ensure full coverage of the relevant surface. Compared to manual techniques, such automated lay-up processes ensure higher productivity, repeatability and better product quality.

The automated lay-up processes can be implemented by using either specifically designed machines or robotic systems, which in this application are usually redundant. In the first case, general CNC machines equipped with specific lay-up head are used. They have no essential limitations on the workpiece size, but usually are quite expensive and require large work floor areas [START_REF] Gallet-Hamlyn | Multiple-use robots for composite part manufacturing[END_REF][START_REF] Rabeneck | Boosting composite production[END_REF]. Compared to the process-dedicated machines, the robotic systems composed of a 6-axis serial robot (mounted on a linear track/gantry), an actuated positioner and a dedicated technological tool are relatively cheap and flexible allowing changing the product type easily [START_REF] Krombholz | Improving the production quality of the advanced automated fiber placement process by means of online path correction[END_REF][START_REF] Shirinzadeh | Fabrication process of open surfaces by robotic fibre placement[END_REF]. For these reasons, the robotic lay-up systems are increasingly used in practice. However, the manufacturing process planning and manipulator motion programming for such robotic systems is not trivial. The main difficulty here arises because the kinematic redundancy of the robotic system, which usually contains two extra degrees of freedom provided by the rotational positioner and the translational unit (linear track or gantry). This redundancy does not allow generating the manipulator trajectories straightforwardly and also complicates robot programming, which is still a challenge in productivity improvement for the automated lay-up process. On the other hand, the redundancy gives user some flexibility in motion planning and optimal utilization of the actuator capabilities allowing generating faster trajectories.

In robotic literature, the issue of the kinematic redundancy resolution has been studied for several decades. A conventional way is based on the pseudo inverse of the kinematic Jacobian [START_REF] Flacco | Discrete-time redundancy resolution at the velocity level with acceleration/torque optimization properties[END_REF][START_REF] Kazerounian | Redundancy resolution of serial manipulators based on robot dynamics[END_REF][START_REF] Wu | Redundancy coordination of multiple robotic devices for welding through genetic algorithm[END_REF][START_REF] Buss | Introduction to inverse kinematics with jacobian transpose, pseudoinverse and damped least squares methods[END_REF]. It provides a unique solution for the differential kinematic equations in the sense of least squares corresponding to the smallest Euclidean norm of the displacement vector in the joint space. Some other approaches [START_REF] Tabarah | Optimal motion coordination of two robots-a polynomial parameterization approach to trajectory resolution[END_REF][START_REF] Gan | Off-line programming techniques for multirobot cooperation system[END_REF][START_REF] Gan | Cooperative path planning and constraints analysis for master-slave industrial robots[END_REF] for non-complex shape tasks are based on coordinating of robot/positioner motions, where the motions of the "master" are assigned firstly and the conjugate trajectories of the "slave" are then determined. For complex shape workpiece taking into account constraints imposed by actuators, several alternative optimization based methods have been developed for the spot-welding application (Gueta et al., 2011a, Gueta et al., 2011b). They are based on converting the original problem into a discrete form, where the robotic manipulator and the positioner joint spaces are discretized and desired trajectory can be represented as a path on the corresponding graph. Slightly different approach was proposed in (Dolgui andPashkevich, 2009, Pashkevich et al., 2004) for the laser cutting applications. It minimizes the oscillation in actuator vi velocities when assuming the tool speed is constant. Hence, the known techniques cannot be straightforwardly applied to the considered lay-up technology where the robot tool speed should be as high as possible. Besides, required functionality is not included in existing commercial software packages (such as Robcad, Robotmaster, CATFiber, etc) allowing generating certain coordinated motions for the redundant robotic systems. For these reasons, the problem of optimal motion planning for robotic lay-up systems becomes very important and corresponds to current needs of the composite product fabrication technology.

Thesis goal and research problems. The thesis focuses on the optimal motion planning in redundant robotic systems allowing improving the productivity of the automated composite lay-up workcell. Special attention is paid to the motion coordination of the robotic manipulator, workpiece positioner and workspace extension unit, which ensures the shortest manufacturing time and smooth movements of all mechanical components. To achieve this goal, the following tasks have to be solved:

Task #1: Analysis of existing systems for composite product manufacturing, detecting their weak points and comparative study of known methods for the motion planning in redundant robotic systems.

Task #2:

Creating kinematic model of a typical redundant robotic system used in the composite layup technology and formalization of the related optimal motion planning problem.

Task #3: Development of the method for the optimal coordinated motion planning in the redundant robotic system composed of the robotic manipulator, workpiece positioner and workspace extension unit, which allows minimizing the total motion time and ensuring smooth movements of all mechanical components. Task #4: Application of the developed algorithms to real industrial problems, development of the robot control programs implementing generated time-optimal trajectories, simulation of the coordinated movements in 3D environment, and experimental validation of the developed techniques on the factory floor.

Thesis organization. To solve the above defined tasks, the thesis is organized as follows:

Chapter 1 presents the state of the art in the field of robot based composite lay-up technology and describes existing techniques of motion planning in robotic systems with kinematic redundancies. The main purpose of this chapter is to detect the weak points of existing methods and justify necessity of developing of new ones, which allow defining the thesis main goal and related research problems.

Chapter 2 focuses on the kinematic modeling of a typical redundant robotic system used in the composite lay-up technology and formalization of the related optimal motion planning problem. The considered problem is presented as finding of the time-optimal trajectory in robotic system joint space vii under specific constraints related to both the robotic system and technological task. These constraints allow taking into account limitations of the robotic system (joint limits, maximum joint velocities/accelerations), practical requirements for the manipulator postures (via singularities and collision constraints) as well as the Cartesian path constraints issued from the composite lay-up process. The main purpose of this chapter is to formalize the considered technological problem and to present it in the form common for the mathematical optimization theory.

Chapter 3 is devoted to the development of a new motion planning method for the redundant robotic systems utilized in the automated composite lay-up process. It is based on the transforming the original problem into a combinatorial one and applying the dynamic programming principle. First, the task graph is created by discretizing the redundant variables and sequentially applying to all task locations the direct kinematics of the positioner (and workspace extension unit) as well as inverse kinematics of the robot. Then, the developed algorithm based on dynamic programming generates the time-optimal motion taking into account all relevant constraints. To reduce the computing time, the time-optimal motion planning is divided in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with large discretization step. Then, the same technique is applied in the local search space, which is created with smaller step in the neighborhood of the obtained trajectory. Alternatively, the second stage may implement a straightforward smoothing of the redundant variable profiles. Efficiency of them and advantages compared to the conventional techniques are confirmed by several case studies.

Chapter 4 deals with industrial implementation of the developed motion planning method on the factory floor (for manufacturing of a high-pressure vessel). Using the proposed method, there were generated time-optimal smooth motions for the manipulator and positioner of the SPIDE-TP robotic platform allowing speeding up the thermoplastic tape lay-up process. Correctness of these motions was verified in 3D simulation environment of CATIA software package. Besides, it was created a program in KRL language implementing these motions and ensuring coordinated control of KUKA KR210 robot and AFPT workpiece positioner. There are presented results of industrial experiments carried out with this program, which confirmed smoothness of the manipulator/positioner movements and essential reduction of the time required for the thermoplastic tape lay-up process.

Finally, Conclusion summarizes the main contributions of the thesis and defines perspectives for future research work.

Main contributions of the thesis. Theoretical contributions of this thesis are in the area of redundancy resolution and time-optimal motion planning for robotic systems composed of the robotic manipulator, workpiece positioner and workspace extension unit. The most essential results can be summarized as follows:

(i)

Formalization of the optimal motion planning problem for typical robotic lay-up platforms.

The proposed approach presents the considered problem as finding of the time-optimal trajectory in the joint space of the robotic system under specific constraints related to the technological task. The latter include both conventional kinematic constraints (joint limits, maximum joint velocities/accelerations), allowable distances to the singularities and obstacles as well as the Cartesian path constraints issued from the composite lay-up process.

viii (ii)

A new method of coordinated time-optimal motion planning for redundant robotic systems ensuring smooth movements of all mechanical components. It is based on transformation of the original optimization problem into a combinatorial one and application of the dynamic programming principle. At the first stage, the redundant variables are discretized and the task graph is created by sequentially applying to all task locations the direct kinematics of the positioner (and workspace extension unit) as well as inverse kinematics of the robot. Then, the developed optimization algorithm based on dynamic programming generates the time-optimal motion taking into account all constraints related to the system kinematics and considered technology. Efficiency of the developed method and its advantages compared to the conventional techniques are confirmed by several case studies.

(iii) Enhancement strategies for the developed motion planning method. To reduce the computing time required for the motion planning, it was proposed to avoid combinatorial search in high-dimensional space. To achieve this target, the time-optimal motion planning is divided in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with relatively large discretization step. Then, the same technique is applied in the local search space, which is created with small discretization step in the neighborhood of the trajectory obtained at the previous stage. As an alternative, the second stage may implement a straightforward smoothing of the redundant variable profiles based on the spline approximation. As follows from relevant study, the proposed enhancement strategies allow essentially reducing the computing time, down to the level acceptable in engineering practice.

(iv) Application of the developed optimal motion planning method to real industrial problems.

Using the thesis results, there were generated time-optimal smooth motions for the manipulator and positioner of the SPIDE-TP robotic platform allowing speeding up the thermoplastic tape lay-up process for manufacturing of the composite high-pressure vessel.

For these motions, it was created a program in KRL language ensuring coordinated control of KUKA KR210 robot and AFPT workpiece positioner. Experimental evaluation of the motions described by this program confirmed smoothness of the manipulator/positioner movements and essential reduction of the time required for the thermoplastic tape lay-up process.

Dissemination of research results. The main results obtained in this thesis have been published in four works and have been presented in three international conferences. Among them, there is a paper in an international journal (Mechanism and Machine Theory) and three papers in international conference proceedings (EUCOMES'2016 -European Conference on Mechanism Science of IFTOMM, ICMIT'2017 -International Conference on Manufacturing and Industrial Technologies, and ICOME'2017 -International Conference on Mechanical Engineering).

The scientific contribution of this thesis have been developed and validated experimentally in the framework of strong cooperation with R&D Department of Polymer and Composite of the technical center CETIM (Centre Technique des Industries Mécaniques), Nantes. 2 Strength plotted against density (yield strength for metals and polymers; compressive strength for ceramics, tear strength for elastomers and tensile strength for composites) [START_REF] Ashby | Materials selection in mechanical design[END_REF] 8 Summary of the related works for motion coordination redundant robotic system ................. Table 9 Software packages for computer aided design of industrial robotic systems ........................... Table 10 Definition of the configuration index for serial 6-axis manipulator. ...................................... Table 11 Kinematic and dynamic constraints describing capabilities of robot KUKA KR210 ............ Table 12 Kinematic constraints of redundant robotic system for the benchmark example ................... Table 13 Comparison of the conventional and developed algorithms using the benchmark example (Robot Motion Time vs. Computing Time) .. This chapter presents the state of the art in the field of robot based composite lay-up technology and describes existing techniques of motion planning in robotic systems with kinematic redundancies. The main purpose of this chapter is to detect the weak points of existing methods and justify necessity of developing of new ones, which allow defining the thesis main goal and related research problems.

FIGURES

COMPOSITES PARTS FABRICATION TECHNIQUES

Composite materials are being increasingly used for primary structures in industrial, aerospace and marine structures because of their outstanding properties and features. Among the existing composites parts fabrication techniques, automated composite lay-up is attractive since it suits the workpiece with arbitrary surfaces. For this process, a robotic system (an industrial robot collaborating with an actuated positioner) replaces the general CNC machines with the advantages of no limitations on the object size and geometry. This section presents the utilizations of the composite materials in industry and the automation of fabrication.

Composite Materials in Manufacturing

Composite materials are combinations of at least two constituent materials with significantly different physical or chemical properties. One of the material is called the reinforcing phase that is usually in the form of fibers (continuous or short), sheets, or particles, and is embedded in the other material called the matrix phase [START_REF] Lubin | Handbook of composites[END_REF][START_REF] Peters | Handbook of composites[END_REF]. Usually, continuous fiber composites are stronger and stiffer than particulate composites [START_REF] Campbell | Structural composite materials[END_REF]. Figure 1 briefly shows how typical continuous fiber composites are made. The two constituents are a reinforcing and a matrix. The reinforcing phase provides the strength and stiffness. In most cases, the reinforcement is harder, stronger, and stiffer than the matrix [START_REF] Campbell | Structural composite materials[END_REF]. The common constitutes of fiber reinforced composites are shown in Table 1. [START_REF] Gay | Composite materials: design and applications[END_REF][START_REF] Christensen | Mechanics of composite materials[END_REF][START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF][START_REF] Jones | Mechanics of composite materials[END_REF][START_REF] Kaw | Mechanics of composite materials[END_REF][START_REF] Lubin | Handbook of composites[END_REF][START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF][START_REF] Peters | Handbook of composites[END_REF][START_REF] Rabeneck | Boosting composite production[END_REF][START_REF] Rao | Carbon composites are becoming competitive and cost effective[END_REF]. Generally, the matrix phase is the main constituent of composite materials (may be up to 70% by volume) mainly responsible for its overall mechanical properties. Accordingly, with the different types of matrix phase, the composite materials may have very different properties, see Table 2. Polymer matrix composites (PMCs) have established themselves as engineering structural materials, which are prominent class of composites compared to other composite materials in commercial applications. The property of high strength (stiffness)-to-weight ratio makes PMCs more attractive than conventional materials in manufacturing, see Figure 2.

Figure 2 Strength plotted against density (yield strength for metals and polymers; compressive strength for ceramics, tear strength for elastomers and tensile strength for composites) [START_REF] Ashby | Materials selection in mechanical design[END_REF] Glass fiber reinforced polymers represent the largest class of PMCs [START_REF] Chawla | Composite materials: science and engineering[END_REF]. An important application is in the use of manufacturing pressure vessels (200kPa), see Figure 3. Heavy steel cylinders usually result in a reduced payload, which are desired to be replaced by much lighter PMC cylinders. Other common applications involve pipes for transportation of water or petroleum. In the offshore industry, they are widely used for risers, stress joints and fluid handling since PMCs offer many advantages over metal because of their high specific strength and stiffness, good durability, low thermal conductivity and good corrosion resistance [START_REF] Pham | A review on design, manufacture and mechanics of composite risers[END_REF][START_REF] Salama | The first offshore field installation for a composite riser joint[END_REF][START_REF] Tamarelle | High-performance composite tubes for offshore applications[END_REF][START_REF] Garoushi | Fiber-Reinforced Composites[END_REF]. Carbon fiber reinforced PMCs are the most important structural composites; especially in the aerospace field [START_REF] Chawla | Composite materials: science and engineering[END_REF][START_REF] Nicolais | Composite materials: a vision for the future[END_REF][START_REF] Deborah | Composite materials: Science and applications[END_REF][START_REF] Pilato | Advanced composite materials[END_REF]. For example, the Boeing 787 and Airbus 350 make great use of composite materials in its airframe and primary structure (an airframe comprising more than 50% carbon fiber reinforced epoxy and other composites) [START_REF] Hale | Boeing 787 from the ground up[END_REF][START_REF] Marsh | Airbus takes on Boeing with reinforced plastic A350 XWB[END_REF], as is shown in Figure 4. By using the composites, it obviously saves the weight (on average of 20% compared to more conventional aluminum designs), fuel and extends range of flying. In addition, such large airplanes, flying at high altitudes, have pressurized cabins. The limit of pressurization depends on the strength of the fuselage material. A fuselage made of PMCs can withstand higher pressure (corresponding to 1800 m altitude) than one made of aluminum (corresponding to 2400 m) [START_REF] Chawla | Composite materials: science and engineering[END_REF]. In wind energy industry, the rotor blade of a wind turbine has the shape of an aerofoil like the wing of an airplane. The material for the skins of the aerofoil needs to be strong, stiff, but light. These requirements lead to fiber reinforced polymer composites as the optimum materials. Goubalt and Mayes [START_REF] Goubalt | Comparative analysis of metal and composite materials for the primary structures of a patrol craft[END_REF] compared composite materials to steel and aluminum to be used for the primary structures of a patrol craft. The corresponding studies find that the structural weight of a patrol craft made of glass-reinforced plastic materials is 10% lighter than an aluminum one and 36% lighter than a steel one of similar size [START_REF] Mäkinen | The development of sandwich structures for naval vessels during 25 years[END_REF][START_REF] Mouritz | Review of advanced composite structures for naval ships and submarines[END_REF][START_REF] Goubalt | Comparative analysis of metal and composite materials for the primary structures of a patrol craft[END_REF]. A prediction is made that the cost of operating a composite craft will be less than for a steel design because of the maintenance reduction (less corrosion) and lower fuel consumption. Additionally, calculated life cycle costs of a composite craft are 7% less than that of a steel one of the same size.

Non-polymer matrix composites here include metal matrix composites and ceramic matrix composites. Metal matrix composites (MMCs) usually use aluminum, titanium or magnesium alloys as continuous matrix, with high modulus of elasticity [START_REF] Pank | Metal-matrix composite processing technologies for aircraft engine applications[END_REF]. In the Hubble telescope, pitch based continuous carbon fiber reinforced aluminum was used for waveguide booms because this composite is very light, has a high elastic modulus and has a low coefficient of thermal expansion. In the US Trident missile, beryllium has been replaced by SiCp/Al composite [START_REF] Chawla | Composite materials: science and engineering[END_REF]. Metal matrix composites can be tailored to have optimal thermal and physical properties to meet the requirements of electronic packaging systems (e.g., cores, substrates, carriers, and housings). Continuous boron fiber reinforced aluminum composites made by diffusion bonding have been used as heat sinks in chip carrier multilayer boards.

Ceramic matrix composites (CMCs) in general have a very attractive property, i.e., resistance to high temperature (800°C1650°C) and corrosive environment. CMCs are more commonly used in heat engines, components requiring resistance to aggressive environments, special electronic/electrical applications, energy conversion, and military systems.

In general, because metal and ceramic matrix composites require very high temperatures and sometimes high pressures for processing, they are normally much more expensive than polymer matrix composites. For this reason, compared to PMCs, the applications of CMCs and MMCs are relatively limited [START_REF] Chung | COMPOSITE MATERIALS: Functional Materials for Modern Technologies[END_REF][START_REF] Campbell | Structural composite materials[END_REF]. In this work, the discussed composite material is common high-performance polymer matrix composites that are widely used in industry.

Automation of Composite Products Manufacturing

In industry, there are several techniques originally developed for making fiber reinforced polymer matrix composites. Since the matrix type affects processing, fabrications for PMCs can be classified into two categories: thermo-set composite processing and thermoplastic composite processing (see Table 3). A thermo-set is a resin with low-viscosity that reacts and cures during processing, forming an intractable solid. A thermoplastic is a high-viscosity resin that is processed by heating it above its melting temperature. Because a thermo-set resin sets up and cures during processing, it cannot be reprocessed by reheating. By comparison, a thermoplastic can be reheated above its melting temperature for additional processing [START_REF] Campbell | Structural composite materials[END_REF][START_REF] Chawla | Composite materials: science and engineering[END_REF][START_REF] Park | Manufacture of Carbon Fiber Composites[END_REF][START_REF] Bratukhin | Composite manufacturing technology[END_REF]. Figure 5 shows the main manufacturing processes for the both.

Lay-up is the most common polymer processing techniques. It can be done either manually (known as hand lay-up) or by automated devices. Fiber reinforcements can be laid onto a mold and the resin (unsaturated polyester is one of the most common) is applied with a brusher or a roller. The schematic is shown in Figure 5a. Hand lay-up is suitable for making a wide variety of composites products from very small to very large, however, it is usually time consuming and labor intensive, which leads to high costs [START_REF] Elkington | Hand layup: understanding the manual process[END_REF]. The labor cost associate with creating composites takes a large portion of the total manufacturing cost [START_REF] Lindbäck | New Automated Composite Manufacturing Process:: Is it possible to find a cost effective manufacturing method with the use of robotic equipment[END_REF][START_REF] Frketic | Automated manufacturing and processing of fiber-reinforced polymer (FRP) composites: An additive review of contemporary and modern techniques for advanced materials manufacturing[END_REF]. A highly trained technician creates around 1kg composite material per hour, but even the best can make a mistake [START_REF] Sloan | ATL & AFP: Defining the megatrends in composite aerostructures[END_REF]. Human error can introduce voids and irregularities into the composite part during production, which negatively affects mechanical properties, making parts unusable for standard operation. For those reasons, manufacturers are transitioning to automated lay-up to improve production volume per mold and cost efficiency [START_REF] Grimshaw | Advanced technology tape laying for affordable manufacturing of large composite structures[END_REF][START_REF] Skinner | Trends, advances and innovations in filament winding[END_REF].

Thermoplastic compression molding

Spray-up is quite similar to hand lay-up in its suitability for making boats, tanks, transportation components, and tub/shower units in a large variety of shapes and sizes. Frequently, resin and fibers (chopped) are fed and sprayed together onto the mold surface through a chopper gun [START_REF] Chawla | Composite materials: science and engineering[END_REF]. The schematic is shown in Figure 5b. This process uses simple, low cost tooling and simple processing. Portable equipment permits on-site fabrication with virtually no part size limitations. Similar to hand lay-up, the production volume per mold is not high, but it also could be improved by automation.

Vacuum bag molding is designed to improve the mechanical properties of multiple layers of fiber reinforcement bonded with a resin after hand lay-ups (resin rich problem). It can be used with wet laminating and pre-impregnated composites molding. In wet laminating case, the fiber reinforcement is saturated using hand lay-up, and then a vacuum bag is mounted on the mold to force out trapped air and excess resin, compact the laminate. In the case of pre-impregnated composites molding, the preimpregnated material is laid up on the mold, a vacuum bag is mounted and the mold is heated or the mold is placed in an autoclave that applies both heat and external pressure, adding to the force of atmospheric pressure. The schematic of this process is shown in Figure 5c. The procedure of preimpregnating, vacuum bagging and autoclave curing is most often used to create advanced composite aircraft and military products [START_REF] Chawla | Composite materials: science and engineering[END_REF][START_REF] Mallick | Fiber-reinforced composites: materials, manufacturing, and design[END_REF].

Compression molding is a process in which the preheated materials are placed over a mold that mounted in a hydraulic or mechanical molding press, then, two halves are closed and pressure is applied. After the hydraulic press releases, the finished piece will be ejected out of the mold. The schematic is shown in Figure 5d. Compression molding enables part design flexibility and features such as inserts, ribs, bosses and attachments. Good surface finishes are obtainable, contributing to lower part finishing cost. Subsequent trimming and machining operations are minimized in compression molding and labor costs are low [START_REF] Mallick | Fiber-reinforced composites: materials, manufacturing, and design[END_REF]. Resin transfer molding (RTM), also called liquid molding, is an intermediate volume molding process. In this method, dry reinforcement material is laid inside the mold, the mold is clamped, and resin is pumped in (through injection ports) under pressure. The schematic is shown in Figure 5e. This process produces complex parts with smooth finishes on all exposed surfaces. The process can be simple or highly automated-and cycle times are speedy. The automotive industry uses RTM to be a cost-effective process for large scale processing, e.g., composite parts of Dodge Viper automobile are made by RTM [START_REF] Chawla | Composite materials: science and engineering[END_REF].

Injection molding is a process in which the resin with short fibers is forced into a heated chamber, and then pushed through a nozzle at the end of the barrel that is pressed against the mold. Once the part inside the mold cools completely, the mold opens, and the part is ejected. The schematic is shown in Figure 5f. An extension of this process is reinforced reaction injection molding (RRIM), in which two (or more) resins are heated separately, combined with short fibers, and then injected into a mold under high pressure and compressed. The types of products that are processed by this technique are wide and varied, ranging from large automotive parts to tiny gears [START_REF] Chawla | Composite materials: science and engineering[END_REF].

Pultrusion forms composites into long, straight shapes like rods or bars. Continuous strands of reinforcement are pulled through a resin bath to saturate them, and then pulled through heated steel molds that sculpt the composites into continuous lengths. The schematic is shown in Figure 5h. It is a continuous process and can be readily automated. Labor costs are low and finished products are very strong. Pultrusion is utilized to make products such as beams, channels, pipes, tubing, fishing rods and golf club shafts (www.compositeslab.com).

Filament winding is an automated process commonly used for manufacturing axial symmetric structures by winding continuous filaments under tension around a rotating mandrel [START_REF] Skinner | Trends, advances and innovations in filament winding[END_REF], as shown in Figure 5g and Figure 6. Filament winding may use either dry fiber passed through a resin bath (wet winding) or pre-impregnated materials (dry winding) for production [START_REF] Abdalla | Design and fabrication of low cost filament winding machine[END_REF]. After the composite part reach the desired thickness, the mandrel then can be cured and removed from the part, leaving the hollow final product. For some products, the mandrel is a permanent part of the finished product as a barrier to protect the composite from the fluid/gas to be stored. Filament winding is well suited to industrial automation. The controlled variables for winding are fiber type, winding angle, tow or bandwidth and thickness of the fiber bundle. The winding angle has an effect on the properties of the final product. A high angle (hoop) will provide circumferential strength, while lower angle patterns (polar or helical) will provide greater longitudinal/axial tensile strength [START_REF] Laval | CADWIND 2006-20 years of filament winding experience[END_REF][START_REF] Fleischer | Joining automotive space frame structures by filament winding[END_REF][START_REF] Hernandez-Moreno | Influence of winding pattern on the mechanical behavior of filament wound composite cylinders under external pressure[END_REF][START_REF] Rousseau | The influence of winding patterns on the damage behaviour of filament-wound pipes[END_REF]. Products currently being produced using this technique range from pipes, oars, bicycle forks, power and transmission poles, pressure vessels to missile casings, aircraft fuselages and lamp posts and yacht masts.

The above presented processes are all used for polymer matrix composite product manufacturing, and their properties are briefly shown in Table 4. It is obvious that, compared to the others, lay-up technique can be applied in a wider range of products, from very small to very large, and it well suits complex and large dimensional components fabrication especially in aerospace fields. However, the production rate and the repeatability of hand lay-up are quite low, which cannot meet the requirements of mass production volumes in practice. For these reasons, to gain the productivity, automation of layup techniques is becoming necessary [START_REF] Rao | Carbon composites are becoming competitive and cost effective[END_REF][START_REF] Mallick | Fiber-reinforced composites: materials, manufacturing, and design[END_REF][START_REF] Olsen | Automated composite tape lay-up using robotic devices[END_REF]. As follows from the literature analysis, at present, there are two main mechanized lay-up techniques used in industry: automated tape laying (ATL) and automated fiber placement (AFP) [START_REF] Frketic | Automated manufacturing and processing of fiber-reinforced polymer (FRP) composites: An additive review of contemporary and modern techniques for advanced materials manufacturing[END_REF]. These two techniques use a computer-aided design and/or computer-aided manufacturing model to build up a specific structure, in a layer by layer process (layer of laminated tapes or tows, resin-impregnated continuous fibers) [START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF][START_REF] Groover | Automation, production systems, and computer-integrated manufacturing[END_REF]. In an automated lay-up process (ATL or AFP), the product geometry is programmed into control system and commonly the molds themselves are made of composite materials.

Automated tape laying (ATL) is a technique mainly used to produce composite structures for large noncomplex (flat or single curvature) parts [START_REF] Beakou | Modelling slit tape buckling during automated prepreg manufacturing: A local approach[END_REF][START_REF] Sloan | ATL & AFP: Defining the megatrends in composite aerostructures[END_REF]. ATL machines lay a resin pre-impregnated tape or continuous fabric strips with widths ranging from 75 to 300 mm onto a flat surface in various orientations [START_REF] Sloan | ATL & AFP: Defining the megatrends in composite aerostructures[END_REF][START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF], as shown in Figure 7. The tool of the lay-up machine usually consists of spools of tape, a winder, winder guides, a compaction roller, tape cutter and a positioner sensor. The first step to produce composite is depositing a starting amount of pre-impregnated tape onto the mold using a soft silicone roller. Then, the machine deposits the tape according to tool paths defining the part geometry. At the end of the lay-up, the tape is cut automatically by rotating pinching blades. Low-rail gantry platforms are also available in medium and large size ranges that can be matched to customer part size and floor space requirements. The machines feature tape-deposit speeds of up to 800 mm/s and a high degree of placement accuracies. This process is realized by computer aided design (CAD). Using the CAD system, the product to be manufactured is developed mathematically onto a surface, which is further broken down into layers to be fabricated by laying tapes side by side. Control software is used to place strips in each layer via a series of numerical control steps to develop the final product shape. Such equipments were initially found in defense related applications, e.g., the wing skin panels of F-22 Raptor fighter jet. With Boeing 787 and Airbus 380 and other advanced aircrafts, automated processing techniques have moved into civilian aircraft construction as well, i.e., make aircraft parts such as wing stringers, spars, skins and elevators, tail skins and horizontal planes, engine cowls, fuselage skins and belly fairings [START_REF] Chawla | Composite materials: science and engineering[END_REF]. Boeing uses the tape laying process wherein strips of carbon fiber are laid on a spinning mandrel by multiple robotic laying heads. Layer upon layer of pre-impregnated strips are laid on the spinning mandrel until desired shape and thickness is obtained. This is followed by curing for about 2 h in an autoclave at around 250 °C [START_REF] Chawla | Composite materials: science and engineering[END_REF]. Figure 8 shows the Boeing 787's unique one-piece composite barrel construction processed with ATL. Automated fiber placement (AFP) is a technique much similar to ATL and it is employed when producing large composite structures with much complex geometries [START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF]. That is because AFP places a number of narrow pre-impregnated fiber tows (2 to 32 tows) that can be steered over sharply curved surfaces whereas wider tapes cannot be so placed without wrinkling some of the fibers [START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF][START_REF] Blom | A theoretical model to study the influence of tow-drop areas on the stiffness and strength of variable-stiffness laminates[END_REF]. This difference is shown in Figure 9. Each tow here is normally driven individually and can be clamped, cut and restarted during processing with low tow tensions. This enables lay-up over complex shapes and tows steering, and is suitable for example in structures such as fuselage sections with window cut-outs, or wing skins with numerous pad-ups and valleys [START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF][START_REF] Izco | High speed tow placement system for complex surfaces with cut/clamp/& restart capabilities at 85 m/min[END_REF][START_REF] Devlieg | High-speed fiber placement on large complex structures[END_REF][START_REF] Denkena | Automated Fiber Placement Head for Manufacturing of Innovative Aerospace Stiffening Structures[END_REF], see Figure 10. Similar to ATL, the tows here can be laid in any orientations and positions so that laminate can be tailored to deliver the strength and stiffness required by the designers at various parts of the structure. Low-rail gantry platforms are also available here in medium and large size ranges that can be matched to customer part size and floor space requirements. Compared to the ATL, the processing speed of AFP is slightly slower (around half of ATL) [START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF], but AFP can place material more effectively over contoured surfaces, it will be key to high-volume production of composite structures with complex shapes [START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF]. As follows from the literature, the above presented automated lay-up techniques (ATL and AFP) are becoming industry standards for producing fiber reinforced polymer matrix composites. AFP improves on ATL by allowing direct lay-up of more complex components. In addition, material wastage rates are reduced and productivity for complex parts is even higher due to the unique operating per tow and steering capabilities. The future for AFP probably lies with machines that are considerably faster. Increases in manufactured part size and complexity, together with the high rates at which the industry needs to fabricate composite parts, have demanded the fiber placement at 0.85m/s and more [START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF].

Industrial Systems for AFP Based Composite Manufacturing

The automated fiber placement process was firstly realized in practice by using specially designed CNC-controlled machines, and some of the AFP system suppliers are companies rooted in CNC machine tools, e.g. Cincinnati, Ingersoll and MTorres. A specific placement head can be driven with at least 6 axes of motion to lay the fiber reinforcements. These machines are usually delivered under different architectures, i.e. gantry type, column type and cantilever type.

A well-known system is the VIPER™ platform made by Cincinnati (USA), as is shown in Figure 11a. With 7 axes of motion, this device particularly suited to flat and highly contoured structures such as cowls, ducts, fuselage panels and barrels, bulkheads, wings, payload adaptors, fan blades, spars, frames and stringers. The latest VIPER 6000 series AFP machines can produce fuselage panels up to 6.3 m in diameter. These VIPER 6000s machines are being used in Boeing for the barrel sections of its B787 fuselages [START_REF] Marsh | Automating aerospace composites production with fibre placement[END_REF][START_REF] Marsh | Airbus takes on Boeing with reinforced plastic A350 XWB[END_REF]. Another system for fiber placement is the TORRESFIBERLAYUP machine developed by MTorres (Spain), see Figure 11b. Similar to VIPER™ platform, this dedicated machine is based on a general CNC machine equipped with a placement head. It is being used for building the Airbus A350 XWB wing front spar. Electroimpact (USA) produced a cell built around a control architecture that permitted the use of two CNC-controlled machines equipped modular AFP heads running simultaneously or independently [START_REF] Flynn | Automated Fiber Placement Machine Developments: Modular Heads, Tool Point Programming and Volumetric Compensation Bring New Flexibility in a Scalable AFP Cell[END_REF][START_REF] Flynn | Production Implementation of Multiple Machine, High Speed Fiber Placement for Large Structures[END_REF]. It is employed for the manufacture of large primary aircraft structures.

Compared to the CNC-controlled fiber placement machines, the robotic lay-up systems are economically attractive and flexible, allowing changing the product type easily [START_REF] Rabeneck | Boosting composite production[END_REF][START_REF] Gallet-Hamlyn | Multiple-use robots for composite part manufacturing[END_REF][START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF]. In comparison to their sizes, they can provide a large working area. These robotic lay-up systems are usually composed of a 6-axis serial robotic manipulator equipped with a specific fiber placement end-effector, an actuated workpiece positioner and a workspace extension unit (linear track/gantry). The actuated positioner can have one or two degrees of freedoms. The product to be manufactured is mounted on the positioner flange that adjusts its posture by rotating the positioner axis. In practice, a one-axis positioner can be sufficient if the product is not too large and its shape is simple. In case of large dimensional product, the linear rail unit is activated to increase the robot workspace.

In the robotic lay-up system developed by Coriolis Composite (France), a placement head is drove by a standard 6-axis poly-articulated robot (located on a rail) for laying and the mold is held by an actuated positioner, see Figure 12. The set composed of 8 axes of motion effectively meet the specifications of complex composite part processing. For feeding the fibers, a creel (or bobbin cabinet) situated at the foot of the robot provides all the necessary functions for unwinding the bobbins at high speed with low tension and enables swift ergonomic loading and unloading of the bobbins. The pipes individually feed each fiber from the creel to the compact light head avoiding all risks of twisting or damage to the fiber, while maintaining a low tension (www.coriolis-composites.com). The Coriolis AFP system is also used to develop structural components in thermo-set and thermoplastic composites as well dry fiber performs in National Aerospace Laboratory of Netherlands. Electroimpact (USA) developed a modular placement head that can be mounted on the robot flange, it avoids using the long tow tube back to the refrigerated "creel house" attached to a major structural element mounted on a linear axis [START_REF] Flynn | Automated Fiber Placement Machine Developments: Modular Heads, Tool Point Programming and Volumetric Compensation Bring New Flexibility in a Scalable AFP Cell[END_REF][START_REF] Rower | Robot Driven Automatic Tapehead for Complex Composite Lay-Ups[END_REF], Rudberg et al., 2011), see Figure 13a. Such architecture is also used by NASA for producing composite parts for the agency's aeronautics and space exploration programs, see Figure 13b. Similar architecture is also launched by Automated Dynamics (USA). Mikrosam (Macedonia) is launching a system that integrates AFP head for complex parts and ATL head for flat molds into a single robotic cell for custom development of structural composites. It allows fast transition from automated fiber placement to automated tape laying and vice versa, by simply changing the head.

Current challenges still exist in design of robot based AFP, although such robotic systems have been widely employed in practice with good flexibility and adaptability. The main challenges can be related to several following terms:  Improvement of the fiber feeding system. There are two representative designs of the fiber feeding system. 1) A "creel house" is situated at the foot of the robot and pipes individually feed each fiber from the creel to the head. It is developed by Coriolis Composite, see Figure 12b. 2) A modular head is made by embedding all the fiber creels into the head. This design is delivered by Electroimpact, see Figure 13a. For the former, feeding fibers from the creel to the head limits the robot movements; and in the latter case, the increased head dimension and weight by modular designing require to use heavy payload robot (relatively slower and less accurate) (Gallet-Hamlyn, 2011). The design of feeding system might be future optimized.  Automation of robot programming. Productivity improvement can be expected from improved robot programming [START_REF] Dirk | The engineering aspects of automated prepreg layup: History, present and future[END_REF]. Since the robotic system usually contains two external axes of motion from rotational positioner and linear rail track, the system is kinematically redundant with respect to the fiber placement task. These two redundancies create some difficulties and complicate the robot programming, but it also gives some space for optimizing the robot and positioner movement to improve the productivity, for example, generating of the fastest trajectory by using the redundant degrees of freedom.

In the thesis, the challenge of productivity improvement via generating time-optimal motion at the stages of manufacturing process planning is mainly concerned and the problem of utilizing the redundant motion axes in best way is solved by generating the fastest lay-up motion for the redundant robotic system.

REDUNDANCY RESOLUTION IN ROBOTIC SYSTEMS

In order to analyze the redundancy problem, let us model the redundant robotic system with a set of links connected by rotational/translational joints. The configurations of joints are described by single scalar angle values. The complete configuration of the robotic system is specified by 8×1 vector [ , , ] T L P R qq  q q where L q describes the configuration coordinate of the linear rail track; P q describes the joint coordinate of the positioner; T R q q q q q q ) , , , , , (

6 5 4 3 2 1
 q describes the angles of the robot joints. A 6dimensional lay-up task location can be described as 6×1 vector

T z y x ] , , , , , [     χ
, where

T z y x ) , , (
is the position coordinate relative to the world frame and

T ) , , (   
is the Z-Y-X Euler angles representing the orientation. Thus, a specific posture of the end-effector can be represented in two ways, i.e.

q from joint space and χ from task space.

The direct kinematic transformation is the mapping from the joint space to the task space. A point in joint space represents a unique location of the end-effector referred to the robot base frame in task space. The task kinematics can be described by

) (q χ g  where ) (
g is the geometric transformation function. The standard Denavit-Hartenberg (DH) technique can be used to build the geometric model of the robotic system [START_REF] Hartenberg | A kinematic notation for lower pair mechanisms based on matrices[END_REF], and the locations of the end-effector are presented as 4×4 homogeneous transformation matrices.

The inverse kinematic transformation is the mapping from the task space to the joint space. Compared to the direct kinematics, the inverse one is not trivial since the solution is not unique and cannot be expressed in a closed form. For this reason, special types of manipulator architecture are used to satisfy the Peiper condition [START_REF] Peiper | The kinematics of manipulators under computer control[END_REF]) that ensures the closed form solution. Otherwise, there are also some numerical techniques to deal with this problem [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]. In this thesis, the serial robot with last three intersecting axes is employed, which corresponds to the majority of industrial applications. In order to obtain a unique solution for the robot, the task kinematics referred to robot base can be described by

) , ( 1 μ g R χ q  
where  is the configuration index that determines the posture of the robot shoulder, elbow and wrist. However, for the whole system with redundancies, it is still impossible to get a unique solution for

) , ( 1 μ g χ q  
because the dimension of q is higher than that of χ .

For automated lay-up processes, the architecture of 6-axis robot plus 1-axis positioner and 1-axis workspace extensioner essentially complicates the preparation stage of robot programming since a location on the product can be reached with infinite numbers of q . To generate the desired motion profile for a given task (lay-up path that is predefined as an augmented line), it is required decomposing the path in task space into robot motion and positioner motion in joint space. In literature, there are several works dealing with the redundancy resolution problems. Some representative approaches are presented as follows.

Redundancy Resolution via Generalized Inverse of the Kinematic Jacobian

The existence of the redundant motion axes leads to the insolvability of the system inverse kinematics, i.e.

) , ( 1 μ g χ q   , that is because the dimension of q is higher than that of χ . In order to generate unique solution for ) , ( 1 μ g χ q   , a redundancy resolution technique based on the generalized inverse of the kinematic Jacobian can be found in literature [START_REF] Kazerounian | Redundancy resolution of serial manipulators based on robot dynamics[END_REF][START_REF] Buss | Introduction to inverse kinematics with jacobian transpose, pseudoinverse and damped least squares methods[END_REF][START_REF] Andres | Implementation and testing of a CAM postprocessor for an industrial redundant workcell with evaluation of several fuzzified Redundancy Resolution Schemes[END_REF][START_REF] Flacco | Discrete-time redundancy resolution at the velocity level with acceleration/torque optimization properties[END_REF][START_REF] Wang | Inverse kinematics and control of a 7-DOF redundant manipulator based on the closed-loop algorithm[END_REF][START_REF] Fahimi | Autonomous robots: modeling, path planning, and control[END_REF][START_REF] Patel | Control of redundant robot manipulators: theory and experiments[END_REF][START_REF] Chiaverini | Redundant Robots[END_REF][START_REF] Siciliano | Kinematic control of redundant robot manipulators: A tutorial[END_REF][START_REF] Fernandez | A generalized method for automatic downhand and wirefeed control of a welding robot and positioner[END_REF][START_REF] Nenchev | Redundancy resolution through local optimization: A review[END_REF]. The most widely known generalized inverse of the kinematic Jacobian used for solving this problem is the Moore-Penrose pseudo-inverse.

The Moore-Penrose pseudo-inverse was proposed as a general way to find the solution to the equation

x A b   where n m n m R R R     A x b ; ;
. In the case of n m  , a unique solution can be obtained with the form

b A x   
, where

1 ) (      T T A A A A
is the Moore-Penrose pseudo-inverse of A (Eldén, 1982, Ben-Israel and[START_REF] Ben-Israel | Generalized inverses: theory and applications[END_REF]. This particular solution has the property of minimizing the Euclidean norm of x , i.e. min ... : [START_REF] Buss | Introduction to inverse kinematics with jacobian transpose, pseudoinverse and damped least squares methods[END_REF].

2 2 2 2 1      n x x x x
To apply this technique, let 8 R  q be the joint coordinates of the considered redundant robotic system and 6 R  χ be the Cartesian coordinates describing the 6-dimensional lay-up task. The task kinematics is given by the direct transformation

) (q χ g 
. By differentiating it with respect to time, the first-order kinematic can be expressed as follows:

dt d dt d q q J χ   ) ( (1)
where 8 6 R J is the Jacobian matrix. It should be noticed that since the system here is redundant, the Jacobian matrix is not square. Then, from the above equations, the rates of displacement q d can be obtained by using the pseudo-inverse of Jacobian as follows

χ q J q d d    ) ( (2)
where 1 ) ( (Whitney, 1969, Tucker and[START_REF] Tucker | Generalized inverses for robotic manipulators[END_REF]. Then, by defining respectively represent the previous and the current system configurations and the corresponding locations in task space, the joint coordinates for each task location can be sequentially computed by using the following equation

     T T J J J J
0 q q q  
) ( 0 0 χ χ J q q      c c (3)
It should be mentioned that the initial 0 0 , χ q have to be pre-determined.

Application Example. To verify the properties of solution from Moore-Penrose pseudo-inverse, let us apply it to a redundant robotic system composed of a two-axis planar robot and one-axis positioner. The geometrical parameters of the systems are given in Figure 14. This system possesses a 1-dof redundancy with respect to the planar task. The ellipsis-type lay-up path was described by a set of task locations, i.e. From the geometry of the robotic system, the position of the end-effector on the workpiece surface is written as the following function of joint parameters

                 ) ( sin ) ( sin sin cos ) ( cos ) ( cos cos sin 2 1 3 2 1 3 1 3 3 3 2 1 2 3 1 1 3 3 q q q L q q L q D q H y q q q L q q L q D q H x (4)
that is the direct kinematics of the robotic system. Then, by differentiating those equations, the Jacobian matrix is written as follows

                                  3 3 1 2 3 2 1 3 1 1 2 3 2 1 2 3 2 1 3 1 3 3 3 2 1 2 3 1 1 3 2 1 2 3 2 1 2 3 1 1 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( q C H Cq D q q q C L q q C L q q q C L q q q C L q q C L q C H q S D
q q q S L q q S L q q q S L q q q S L q q S L J [START_REF]Operating and Programming Instructions for System Integrators[END_REF] where C and S are shorthand for cos and sin.

By applying the equations ( 2) and (3) presented above, a particular solution with the minimization of q d is quickly obtained. However, when implementing the solution in Matlab software, the generated motion leads to collisions (see Figure 15) that are completely not acceptable.

Figure 15 Simulation of the solution from pseudo-inverse for planar redundant robotic system: (a)  (h) From the above application example, it is noticed that the redundancy resolution technique via generalized inverse of Jacobian can quickly generate a solution with the minimization of instantaneous power, since the norm of joint velocities relate to the energy consumption. However, it must be mentioned that the solution from this approach usually leads to unreasonable manipulator behavior [START_REF] Padula | A new pseudoinverse for manipulator collision avoidance[END_REF]Perdereau, 2011, Klein and[START_REF] Klein | Review of pseudoinverse control for use with kinematically redundant manipulators[END_REF], i.e. collisions between system components in the above example. Besides, this pure mathematical technique does not take into account the velocity and acceleration constraints of actuators. For these reasons, the technique based on generalized inverse of Jacobian cannot be straightforwardly applied for our technical problem (automated lay-up process).

Redundancy Resolution via Coordinating of Robot/Positioner Motions

Another idea for solving the inverse kinematics of the redundant system is to geometrically decompose the Cartesian path into a robot motion and a positioner motion. Related research have been carried out to develop robot-positioner coordinated motion planning since the middle of 1990s [START_REF] Holmes | Method and apparatus for controlling manipulator and workpiece positioner[END_REF][START_REF] Jouaneh | A kinematic approach for coordinated motion of a robot and a positioning table[END_REF], Jouaneh et al., 1990a, Jouaneh et al., 1990b[START_REF] Ahmad | Coordinated motion control of multiple robotic devices for welding and redundancy coordination through constrained optimization in Cartesian space[END_REF][START_REF] Alford | Coordinated control of two robot arms[END_REF].

A method of trajectory planning for a robot and a positioning table coordinated motion was proposed by Jouaneh [START_REF] Jouaneh | A kinematic approach for coordinated motion of a robot and a positioning table[END_REF] in welding applications. The planar tool motion was resolved among a 5-axis manipulator and a 2-axis table, by using two geometrically defined parameters a i and b i representing the fraction of the incremental motion in the x and y axis directions, respectively moved by the robot in the i th intervals along the welding path. Then, depending on the given path information, Jouaneh (Jouaneh et al., 1990b, Jouaneh et al., 1990a) extended this approach to two strategies for the motion coordination problem. The first one is driving two devices in opposite direction for simple path shapes, e.g. a straight line or a gentle curve. The second strategy that is used for cornered-paths, resolves the path with sharp corners into smooth ones. By this way, the trajectories for both of the table and the robot can be obtained in an off-line mode.

Another way of coordinated motion planning for redundant robotic systems [START_REF] Tabarah | Optimal motion coordination of two robots-a polynomial parameterization approach to trajectory resolution[END_REF][START_REF] Gan | Off-line programming techniques for multirobot cooperation system[END_REF][START_REF] Gan | Cooperative path planning and constraints analysis for master-slave industrial robots[END_REF] is based on the idea of "master-slave", where the trajectories of the "master" manipulator are assigned firstly, and the corresponding conjugate trajectories of the "slave" are then determined. This idea is simple and computationally efficient, but assigning the master trajectory is not trivial, especially for complex shape objects. Besides, here it is not possible to take into account the actuator constraints in an explicit way.

Application Example. To verify the properties of solution from motion coordination techniques, let us apply it to the same redundant robotic system that is already used in previous example. The geometrical parameters of the systems can be found in Figure 14. The ellipse-type lay-up path was described by a set of task locations, i.e.

T i i i y x ) , (  p where 100 ... , 2 , 1  i .
The kinematics of the robotic system here is separately described by robot and positioner. From the geometry of the robot, the position of the end-effector with respect to the robot base frame (also selected as the world frame) is written as the following function of joint parameters

         ) ( sin sin ) ( cos cos 2 1 2 1 1 2 1 2 1 1 q q L q L y q q L q L x (6)
that is the direct kinematics of the robot. Accordingly, its inverse kinematics is written as follows

             ) 2 ( cos a ) cos , sin ( 2 tan a ) , ( 2 tan a 2 1 2 2 2 1 2 2 2 1 2 2 2 2 1 L L L L y x q L q L q L x y q  (7)
where

1    and 1 2 2 1 2 2 2 1 2 2     L L L L y x
. Also, the perspective of the positioner, the task points can be described as follows

           H q y q x y D q y q x x i i i i i i 3 3 3 3 cos sin sin cos (8)
that is the direct kinematics of the positioner.

     180 180 1 q 1 1 1 20 20          s q s       180 180 2 q 1 2 1 30 30          s q s       180 180 3 q 1 3 1 40 40          s q s 
To generate the coordinated motion, let us assign the positioner motion firstly as determining the robot motion is relatively harder. Here, the positioner actuator is simply driven at the maximum speed (see Table 6) within its pre-determined range, and the robot coordinates with the positioner. The positioner motion is described as

) ( 3 i t q where 100 ... , 2 , 1  i
. Then, by applying the direct kinematics of the positioner and the inverse kinematics of the robot sequentially, corresponding

) ( i t q
can be obtained. The solution implemented in Matlab is shown in Figure 16.

Figure 16 Simulation of the solution from motion coordination for planar redundant robotic system: (a)  (h) From the above application example, it is found that robot and positioner motion generated via geometrically decomposing the Cartesian path shows a reasonable system behavior. Compared to the previous solution from generalized inverse of Jacobian, it avoids the collisions since these physical constrains can be considered when assigning trajectories for the two devices (in this example, if the rotational range of the positioner is limited to a small value, collisions will appear). However, it should be mentioned that for automated lay-up process, it is not always easy to decompose the placement path into a robot motion and a positioner motion since the product shape may be irregular and complex. Besides, this idea does not allow taking into account the velocity and acceleration constraints of actuators. For these reasons, the techniques based on this idea are not suitable for the considered applications.

Optimization Based Techniques for Redundancy Resolution

Besides the previously presented techniques, optimization-based methods are also frequently used for the redundancy resolution problem. This scheme tries to fully exploit the kinematic redundancy of the system in order to improve robot motion performance by using some criteria, e.g. minimizing the positioning error, obtaining smooth solutions in joint space, minimizing the displacements of the actuators, obtaining collision-free joint paths, etc. [START_REF] Za'er | Cartesian path generation of robot manipulators using continuous genetic algorithms[END_REF]. Several related work can be found in literature, and briefly shown in Table 6.

By reviewing the techniques shown above, it can be found that most of such techniques are based on heuristic search algorithms (mostly genetic algorithms, GAs) with the objective of minimizing the positioning errors. In these approaches, the difficulties of solving inverse kinematics of the redundant system are avoided. In order to analyze the optimization-based techniques, let us briefly classify them into two categories: 1) for point-to-point motion generation; 2) for continuous-path generation.

Publication

Method, objective function, constraints [START_REF] Parker | Inverse kinematics of redundant robots using genetic algorithms[END_REF] GA;

Minimize positioning error and the maximum joint displacement;

Point-to-point motion; applied to 4R robot for (x, y, z) control. [START_REF] Nearchou | Application of genetic algorithms to point-to-point motion of redundant manipulators[END_REF][START_REF] Nearchou | Solving the inverse kinematics problem of redundant robots operating in complex environments via a modified genetic algorithm[END_REF][START_REF] Nearchou | A genetic path planning algorithm for redundant articulated robots[END_REF] GA;

Minimize positioning error;

Collision-free constraint;

Point-to-point motion; applied to 7R robot for (x, y) control.

( [START_REF] Za'er | Cartesian path generation of robot manipulators using continuous genetic algorithms[END_REF] Continuous GA;

Minimizing accumulative deviation between the generated and the desired paths;

Continuous-path motion; applied to 3R robot for (x, y) control / 6R robot for (x, y, z, a, b, c). (Ata and[START_REF] Ata | Optimal point-to-point trajectory tracking of redundant manipulators using generalized pattern search[END_REF][START_REF] Ata | Collision-free trajectory planning for manipulators using generalized pattern search[END_REF] Generalized pattern search;

Minimizing positioning error, joint displacement;

Continuous-path motion; applied to 3R robot for (x, y) control; [START_REF] Pires | Manipulator trajectory planning using a MOEA[END_REF], Pires and Machado, 2000a, Pires and Machado, 2000b) GA;

Minimizing positioning error, joint displacement, energy, etc.;

Point-to-point motion; applied to 3R robot for (x, y) control.

(da Graça [START_REF] Da Graça | Trajectory planning of redundant manipulators using genetic algorithms[END_REF] GA with closed-loop pseudo inverse;

Minimizing positioning error, maximum joint displacement, etc.;

Point-to-point motion; applied to 3R robot for (x, y) control. [START_REF] Menasri | A trajectory planning of redundant manipulators based on bilevel optimization[END_REF] Bi-GA;

Minimizing positioning error and maximizing the robot manipulability;

Collision-free constraint;

Point-to-point motion; applied to 3R robot for (x, y) control. [START_REF] Debout | Tool path smoothing of a redundant machine: Application to Automated Fiber Placement[END_REF] Optimizing redundant axis with Levenberg-Marquardt method;

Minimizing tool path length and the curvature variation;

Continuous-path motion; applied to 7R robot for (x, y, z, a, b, c) control. (Gueta et al., 2009a, Gueta et al., 2008b) Graph based search space representation;

Minimizing total travelling time;

Point -to-point motions; applied to 6R robot and 1R table for (x, y, z, a, b, c) control. [START_REF] Dolgui | Manipulator motion planning for high-speed robotic laser cutting[END_REF][START_REF] Dolgui | Enhanced Optimisation Techniques for Off-line Programming of Laser Cutting Robots[END_REF][START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF] Graph based search space representation;

Minimizing joint displacements;

Continuous-path motions; applied to 6R robot for (x, y, z, a, b) control.

To generate point-to-point motions, a representative approach of redundancy resolution was presented by [START_REF] Parker | Inverse kinematics of redundant robots using genetic algorithms[END_REF]. This GA based method minimizes the positioning error and the maximum joint displacement by applying following expression

min } { max      f i f i k q q χ χ (9)
where i χ and f χ represent the initial and final Cartesian position of the end-effector; i q and f q are the initial and final joint coordinate of the robot; factor k scales the contributions from each term. In order to take into account the collision constraint, Nearchou et al. [START_REF] Nearchou | Application of genetic algorithms to point-to-point motion of redundant manipulators[END_REF][START_REF] Nearchou | Solving the inverse kinematics problem of redundant robots operating in complex environments via a modified genetic algorithm[END_REF][START_REF] Nearchou | A genetic path planning algorithm for redundant articulated robots[END_REF] presented a similar GA based approach while the geometries of obstacles are described in constraint function. It allows generating a point-to-point solution subject to a collision-free movement.

To generate joint trajectories for a continuous-path in task space, that is more interested in automated lay-up process, Za'er et al. [START_REF] Za'er | Cartesian path generation of robot manipulators using continuous genetic algorithms[END_REF] proposed a technique based on continuous genetic algorithm taking into account the singularity-free constraints, which minimizes the accumulative deviation between the two paths given by the following formula

      n i N k gc dc i k P i k P E 1 1 ) , ( ) , ( (10)
where n is the number of the task point on the path; N represents the number of the joint axis; P dc is the desired Cartesian location and P gc is the generated Cartesian location. The optimal solution of the problem is obtained when the deviation function, E, approaches zero. [START_REF] Ata | Optimal point-to-point trajectory tracking of redundant manipulators using generalized pattern search[END_REF] presented a multi-objectives optimization technique that is realized by implementing GA and GPS (generalized pattern search) to lead the end-effector along a straight path or a circle path. The main objective is to minimize the sum of the positioning error of the endeffector at each intermediate point, and the sub-objectives are expressed as joint displacement, velocities, etc. The objective function is expressed as follows

    m i i i E w F 1 (11)
where w i is the weighting factor to control the desired configuration which satisfy the constraint   m i i w 1 , and m is the number of the objectives.

Application example. The case study presented above can be also used to test the GA-based technique (see Figure 14), and the geometry of the robotic system is described as Equation ( 4). The problem can be formulated as follows

                                    free on intersecti 100 ... , 2 , 1 180 180 . . ) , , ( ) , , ( ) 1 ( min : find ) ( 3 , 2 , 1 99 1 3 1 ) ( ) ( ) 1 ( ) ( 100 1 ) ( 3 ) ( 2 ) ( 1 ) ( 3 ) ( 2 ) ( 1 i q t s q q q q q g y q q q g x F i i j i j i j i i i i y i i i i x i   (12)
where x g and y g are geometric model of the system, based on Equation (14);

T i i i y x ) , (  p
is the i th point on the Cartesian path, F can be set as the fitness value of genetic algorithm. Then, by applying the standard genetic algorithm, i.e. ga(•) in Matlab, with the setting parameters shown in Table 7 and randomly generated initial solution, a result for 300 decision variables is generated within around 10 hours (Matlab2014b in Win7SP1 with Intel ® i5 @2.67GHz 2.67GHz). Figure 17a visualizes the progress of the algorithm, which is provided in the options of ga(•). It shows that the algorithm halted when reached 3500 generations, with best fitness value F=14.6137. However, such straightforwardly obtained solution is not ideal. As is shown in Figure 17b, the positioner joint coordinate for 100 task points does not vary continuously which brings sharp jumps in the trajectories of the robotic system. This phenomenon is possibly caused by the randomly generated initial values. Another reason is that the motion smoothness is hardly expressed as a constraint during the path searching. Besides, this solution does not provide acceptable positioning accuracy. By analyzing several solutions obtained from GA-based technique, it can be found that, using such GA-based straightforward techniques, only the direct kinematic equation of the robotic system is required, and the difficulties related to the inverse kinematics does not exist here. Also, singular configurations and collision cases can be avoided by adding the constraint function. Obviously, it brings some simplicity for the redundancy resolution. However, these straightforward mathematical methods do not take into account properties of the systems kinematics, which leads to some unreasonable movements. Additionally, it has to be noticed that the result qualities are strongly related to the initial values of the variables. But, in automated lay-up processes, selection of these initial values is non-trivial since the Cartesian path in task space is usually complex. A reference test is done on a simple task, i.e. straight line, with simply estimated initial values, better solutions can be obtained. It proves that this straightforward GA based approach is more suitable for the tasks with non-complex Cartesian path.

Besides the above illustrated techniques, there are some other types of optimization based methods in literature. For example, a method was developed for fiber placement process with a specific designed 7-axis robot [START_REF] Debout | Tool path smoothing of a redundant machine: Application to Automated Fiber Placement[END_REF]. The redundancy resolution is performed by optimizing one axis control value all along the tool path using a least square optimization algorithm. And then, the other 6 axes values can be computed analytically using inverse kinematics. The technique generates a local solution in the objective of fastest tool path after several iterations and evaluations of an objective function which takes into account the variations of several axes. Obviously, compared to the previous optimization-based techniques, this approach ensures exactly positioning, and also contributes to the trajectory smoothness in the configuration space and manufacturing efficiency by minimizing variations of the joint coordinates. But, it does not allow using full capacities of all the actuators in the robotic system. Another efficient strategy is based on converting the problem into a discrete one, where the robotic manipulator and the positioner joint spaces are discretized and the desired trajectory is represented as the shortest path on the corresponding graph. One of such techniques was developed by Gueta (Gueta et al., 2011a, Gueta et al., 2009b, Gueta et al., 2008a, Gueta et al., 2009a, Gueta et al., 2011b, Gueta et al., 2008b[START_REF] Gueta | Compact design of a redundant manipulator system and application to multiple-goal tasks with temporal constraint[END_REF] to generate a set of point-to-point motions for a robot and a rotary table in the objective of minimum travelling time in multi-goal applications. Slightly different approach was proposed in [START_REF] Dolgui | Manipulator motion planning for high-speed robotic laser cutting[END_REF][START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF][START_REF] Dolgui | Enhanced Optimisation Techniques for Off-line Programming of Laser Cutting Robots[END_REF] for laser cutting applications with continuous-path task in Cartesian space. In this case, it minimizes the oscillations in actuator velocities, taking into account the collision constraint when assuming the tool speed is constant.

Table 8 Summary of the related works for motion coordination redundant robotic system

Methods

Properties and applications

Pseudo inverse of the kinematic Jacobian Singular configuration may be generated;

Constraints are hardly applied;

Cannot be straightforwardly used in industry.

Robot/Positioner motion coordination

Cannot use the redundancy optimally;

Can be used in the case of simple path task.

(e.g. some arc-welding applications, etc. )

Optimization-based technique (GA based straightforwardly search)

Unreasonable movements may be generated;

Might be used for pick and place applications;

Might be used in the case of simple path task.

Optimization-based technique (Graph based search space representation)

Can be used for continuous-path task;

(e.g. laser cutting, arc-welding, etc.)

Can be used for multi-goal task;

(e.g. spot-welding, etc.)

By reviewing all the techniques illustrated above, a small summary for the redundancy resolution techniques is presented in Table 8. It shows that only the graph based search space representation can be applied on the automated lay-up processes since it allows generating trajectories for tasks with complex path in task space. However, the aforementioned techniques considering graph based search space representation were respectively designed for multi-goal tasks and continuous-path tasks assuming constant tool speed cases. They do not allow taking into account the kinematic constraints of the robotic system completely, i.e. velocity and acceleration constraints of all the actuators. In addition, they cannot meet the demands of high-speed automated lay-up processes since they do not generate time-optimal motions for a continuous-path task.

To our knowledge, there are no techniques directly addressing the problem of the time-optimal motion planning for robotic lay-up system. It will be very significant to develop a comprehensive methodology to achieve time-optimal smooth motion generation in redundant robotic system for the robotic automated lay-up applications, which is in the focus of the thesis.

ROBOTIC CELL DESIGN AND PROGRAMMING FOR COMPOSITE PRODUCT MANUFACTURING

In order to implement the robot-based automated lay-up process on the factory floor, CAD/CAM software is usually used for both designing a composite product, creation of the robotic system and programming movements of relevant equipment (robot, positioner and linear unit). The CAD software creates models and assemblies that are used by the CAM portion to generate tool paths that drive the robotic system to turn the designs into physical parts. This section provides some details on the manufacturing process preparation as well as robotic system design and programming using modern software tools. It also highlights some difficulties arising here, which are resolved in the following chapters of the thesis.

Manufacturing Process Planning for Robotic Lay-up Applications

Computer aided manufacturing process planning is the systematic determination of manufacturing methods and operations details by which parts can be produced economically and efficiently from raw materials to finished products [START_REF] Leondes | Computer-Aided Design, Engineering, and Manufacturing: Systems Techniques and Applications[END_REF]. According to (ElMaraghy and Nassehi, 2014), the main focuses of this procedure are the optimal selection of equipment, proper tuning of the process parameters and creation of numerical control codes. Usually, the manufacturing process planning involves a series of key steps that can be itemized as follows [START_REF] Alting | Computer Aided Process Planning: the state-of-the-art survey[END_REF]Zhang, 1989, Bagge, 2014): (i) interpretation of the product design data; (ii) selection of the process; (iii) selection of equipments and fixtures; (iv) generation of movements ensuring desired process; (v) creation of process sheets, including numerical control programs.

For the robotic lay-up processes (see Figure 18), the manufacturing process preparation stage includes the workpiece 3D modeling, generation of the required path for the technological tool and the path presentation in the form of sequence of the robot end-effector locations. In this work, it is assumed that the CAD model of the manufacturing task is known and was already created using dedicated computer aided design system. In particular, the desired lay-up path ensuring covering of the workpiece by fiber reinforcements is generated as a 3D augmented curve on the surface of the product. This curve defines both positions and orientations of the technological tool while it is moving along the path. For further convenience, this augmented line is discretized and is presented as a sequence of location vectors (or homogeneous matrices) describing the postures of the technological tool with respect to the product. The density of the discretization is defined by the user, as well as the type of the sampling (regular or irregular).

Figure 18 Manufacturing process planning for robotic lay-up processes

The next stage deals with the design of robotic workcell. It includes selection of the robotic system components (manipulator, positioner, workspace extension unit, technological tool, etc.) as well as their space arrangement. The components selection essentially depends on the product size and its geometry. In particular, an AFP head guiding multiple fibers is more suitable for complex shapes, while an ATL head is more appropriate for lay-up over flat surfaces (with high curvature radius). To manipulate the technological tool, which must be properly positioned and oriented with respect to the product, serial industrial robots with 6 degrees of freedom are usually used. Their selection is mainly based on the workspace size that must suit product geometry. Besides, the robot payload is also taken into account to ensure capacity of the tool handling. To implement the desired lay-up path that usually includes a number of quasi-circular segments, a rotational positioner is commonly employed. Its selection depends on the product size/weight. Besides, for some large dimensional products, a robot workspace extension unit (linear track, etc.) is used, which is selected to be in agreement with the product size. After selecting the workcell components, the designer decides on their mutual location that must ensure accessibility of the lay-up path avoiding collisions between the technological tool, manipulator, positioner and the workpiece. The latter procedure is also often referred to as the robotic cell layout design.

At the third stage, the geometric models of the workcell components and the given lay-up task are integrated. This operation requires direct/inverse kinematic models of all workcell components and relevant geometric presentation of the lay-up task, which must be included in a basic relation describing the closed kinematic loop "manipulator-tool-workpiece-positioner-manipulator". It is clear that this kinematic relation includes some redundant variables and it must be treated as the principle constraint while planning the motions in the considered robotic systems. More details concerning the integration of the task and workcell models are given in Chapter 2.

At the fourth stage, planning of robot and positioner motions, the movements of the workcell components (robotic manipulator, positioner and the workspace extension unit) are determined taking into account the design objectives and technological/physical constraints. The latter include the lay-up path specification, physical limits of the actuators, geometrical limits on the actuating variables as well as some additional limitations describing minimum distances to singularities and collisions. It is obvious that the considered motion planning problem cannot be solved in a unique way because of kinematic redundancy, which gives some space for optimization of workcell components movements. This optimization is the main issue studied in this work; relevant contributions are presented in Chapter 3.

At the fifth stage, robotic cell simulation, the obtained optimal motion is carefully examined using industrial CAD/CAM packages (such as Robcad, DELMIA, etc.). This allows designer to verify accurately the manipulator/positioner joint limits, the collision constraints and also examine the lay-up path visually. If some problems are detected during the simulation, the designer can modify some settings in the motion planning algorithm and repeat the previous stages.

Finally, at the sixth stage, the verified optimal motion is described using the language of the robotic system controller. Relevant post-processors are usually integrated in the robotic simulation packages allowing the user creating the control program, where the motions of the manipulator, positioner and linear unit are presented as the sequence of linear, circular or spline-based segments.

It should be mentioned that the above described procedure of the robotic lay-up process planning is iterative. In particular, some steps can be repeated several times before arriving to the suitable result. Besides, there are a number of commercial software packages that can be applied to some of the process planning stages. Their brief review and applicability analysis is presented in the next subsection.

Software Packages for Computer Aided Design of Industrial Robotic System

At present, there are a number of software packages on the market that help users to design robotic cells for some specific applications such as arc or spot welding, water-jet or laser cutting, painting, spraying, etc. They are available either as standalone applications or as a set of modules embedded in the universal CAD/CAM systems. Usually they offer the users a 3D interactive graphical simulation environment and some convenient tools for robotic cell layout design, manipulator motion planning, workcell simulation and offline programming. The most common of these packages are summarized in Table 9. In early works devoted to this subject, computer aided design of robotic systems was mainly related to the workcell layout [START_REF] Zhang | Challenges and key technologies in robotic cell layout design and optimization[END_REF]. The primary attention to this subject was caused by its high influence on the system productivity, since an ill-placed robot risks inefficient operation and even failure. To make a proper layout design with the aid of graphical simulation software, there are some critical issues highlighted by [START_REF] Lueth | Automated planning of robot workcell layouts[END_REF] who mentioned that the well placed robot must reach all task points without collisions and ensure collision-free movement in their neighborhood. The simplest strategy to find optimum robot location is based on straightforward testing a 2D-grid of possible robot base positions [START_REF] Barral | Optimization tool for robot placement[END_REF]. This test is usually performed either automatically or in an interactive mode, using relevant interface. More sophisticated modern software packages also allows user to generate collision-free movements, optimize the robot cycling time, simulate the manipulator motions, create and debug the robot control program, etc. Let us give some details in the most common software packages for computer aided design of the robotic systems.

One of the most popular commercially available packages, Robcad (www.siemens.com/tecnomatix), is a robotic cell design, simulation and programming system, which provides a platform to handle the task throughout the various steps of the process planning. Robcad has a comprehensive library of standard robots, technological tools, machines and equipments from KUKA, ABB, FANUC, YASKAWA, etc. Its layout design tool allows displaying graphically the robot envelope and indicating whether the task points are reachable or not (showing them in red or green). Using such environment, the user can easily modify the robot placement in an interactive mode to ensure implementation of the given task. In some cases, a suitable robot location can be found automatically. In addition, Robcad is able to generate some typical motions taking into account specific features of the robot controller. With the RRS (realistic robot simulation) module, it offers extremely accurate cycling time calculation and online collision detection. At the final stage, Robcad OLP (offline programming) module transforms the obtained motion into the program codes suited to the relevant control system (over 200 types from KUKA KRC, ABB IRC, YASKAWA YRC, etc.).

It should be also mentioned that Robcad allows data exchange with the most mainstream CAD systems by supporting such data formats as JT, IGES, DXF, STL and STEP. This interoperability allows the manufacturing task created in a dedicated CAD environment to be easily imported. Besides, Robcad includes several process-specific modules such as Robcad Spot, Robcad Arc, Robcad Paint and Robcad Cut, which concentrate on the robot-based automation of spot/arc welding, painting, sand blasting, shot peening, flaming, thermal spraying, laser/water-jet/plasma cutting, sealing and gluing. For example, for the spraying process, Robcad Paint enables generating the robot path, verifying access to all areas, determining coverage parameters and thickness, creating and adjusting process triggers, simulating and downloading the optimized program to the shop floor. However, for the composite layup application studied in this work, there is no specific tool integrated in Robcad.

Another software package, Robotmaster (www.robotmaster.com), provides similar functionalities such as the workcell layout design, automated robot motion optimization, workspace simulation and control code generation. Its application areas include numerous metal processing technologies such as machining, cutting, trimming, polishing, de-burring and also painting, spraying, surface treatment, etc. For example, for the robot-based machining, the user can define the cutter shape, adjust the cutter diameter as well as the depth and the number of cuts. This adjustment is simplified with the modifiable screens that define the interaction, terminology and control setting. At the motion generation stage, the cutter orientations and manipulator configurations can be managed automatically to minimize the wrist rotation, maximize the robot reach and ensure the optimal milling trajectories. At the final stage, the desired robot motion is simulated and programmed automatically. In addition, Robotmaster contains an extensive set of programming and simulation tools for external axes such as rails and rotaries, which allows user to design workcells with robots mounted on linear tracks/gantries or with workpieces mounted on actuated positioners, where up to eight axes are controlled simultaneously. However, Robotmaster cannot address the entire composite layup process considered in this work, and it does not allow generating the desired optimal motions for the technologies studied here.

Among the robotic CAD systems, it is also worth mentioning DELMIA IGRIP (www.3ds.com) from Dassault. It deals with robotic arc welding and allows generating automatically a manipulator end-effector path based on the welding seam geometry and kinematics of the welding positioner, which ensures workpiece optimal orientation with respect to gravity. Its Arc Weld Macro Programming (AMP) module allows the user to generate robot trajectories for complex welding seams, validate them in simulation environment, create the control program, download it into the robot controller and debug the program on the shop floor. Other similar software packages widely used in industry are Workspace5, KUKA.SIM and ABB RobotStudio. The latter two developed by the primary robot manufacturers provide very realistic simulation and programming environment while using their own equipment, but their functionalities for different technological processes are rather limited [START_REF] Hasenjaeger | Programming and simulating automated fiber placement (AFP) CNC machines[END_REF].

For the technology considered in this work, the robot-based composite laying, there is very small number of software systems supporting design and offline programming of robotic cells. One of the most known is CATFiber integrated in CATIA/DELMIA (www.coriolis-composites.com). It allows user to design composite parts in CATIA and simulate automated fiber placement cells with DELMIA offline programming module. Within this package, it is possible to define the desired fiber propagation pattern (parallel, geodesic, etc.) and to set the required fiber placement parameters (steering radius, angular deviation, minimum tape length and fibers overlapping). Similar to other robotic CAD systems, DELMIA includes modules for the workcell components selection (robot, linear track/gantry, technological tool and workpiece positioner from ABB or KUKA), robotic cell layout design, workcell motion planning, simulation and offline programming. However, handling redundant degrees of freedom caused by linear unit and positioner creates some difficulties and requires essential effort of the user who improves the trajectory interactively. For non-CATIA users, Coriolis Composite also provides a package CADFiber directly interfaced with SIEMENSNX. Another software package devoted to fiber placement and tape laying processes, FPM/FPS suite in conjunction with SolidWorks, was developed by Automated Dynamics but the latter has very limited functionality for the kinematic redundancy resolution and optimal motion planning for the layup process.

Hence, there still exists a gap between capabilities of the commercial robotic software packages and requirements of the particular technology. There is no comprehensive robotic CAD system on the market that implement the complete process planning procedures for the robot-based lay-up processes, despite some common functionalities such as workcell layout design, collision check and cycle time evaluation that are available in all above mentioned software. In our case, the robotic system contains one or two redundant degrees of freedom, which makes the motion planning very complicated by using commercial robotic CAD systems. In particular, the existing robotic software packages are not able to generate automatically optimal manipulator/positioner motions desired in the automated layup process. So, in practice, the redundancy resolution is usually performed in an interactive mode [START_REF] Pan | Recent progress on programming methods for industrial robots[END_REF]. For these reasons, this thesis is aimed at developing optimal motion planning methods for robot-based layup processes, which should allow generating time-optimal movements of the robot, linear track and positioner using acceptable computing efforts.

THESIS GOAL AND RESEARCH PROBLEMS

In automated manufacturing of polymer matrix composite products utilization of robotic system is economically attractive since they allow easily changing the product type. Typical robotic systems for such applications are kinematically redundant. Usually they are composed of a 6-axis manipulator, a 1-axis positioner and a 1-axis workspace extension unit (a linear track or gantry), which yields two extra degrees of freedom with respect to the considered technological task. This redundancy provides some convenience for the robotic cell programming but it also creates some difficulties in the motion planning for the mechanical components. On the other hand, the redundancy gives some space for increasing the workcell productivity by optimal coordination of their movements and using full capacities of the system actuators.

In the considered technology, a manipulator usually implements the lay-up process that requires smooth continuous motion of the fiber (tape) laying head along the specified path with the maximum achievable speed. In literature, there are known several motion planning techniques for similar redundant systems but they assume that the manipulator end-effector speed is constant. Hence, they cannot be applied directly for the considered processes where the laying head speed can vary in some degree and minimization of the processing time is required. Besides, existing commercial CAD/CAM software packages cannot solve perfectly the problem of the optimal motion planning for the lay-up process. To our knowledge, at present there are no techniques or software products, which directly solve the problem of the time-optimal motion planning for robotic systems in the lay-up applications.

For the above mentioned reasons, the thesis focuses on the optimal motion planning in redundant robotic systems allowing improving productivity of the automated composite lay-up workcell. Special attention is paid to the motion coordination of the robotic manipulator, workpiece positioner and workspace extension unit, which ensures the shortest manufacturing time and smooth movements of all mechanical components. To achieve this goal, the following tasks have to be solved:

Task #1: Analysis of existing systems for composite product manufacturing, detecting their weak points and comparative study of known methods for the motion planning in redundant robotic systems.

Task #2:

Modeling of typical redundant robotic system used in the composite lay-up technology and formalization of the related optimal motion planning problems.

Task #3: Development of the optimization algorithms for coordinated motion planning in the redundant robotic system composed of the robotic manipulator, workpiece positioner and workspace extension unit, which allow minimizing the total motion time and ensuring smooth movements of all mechanical components.

Task #4:

Application of the developed algorithms to real industrial problems, development of the robot control programs implemented generated time-optimal trajectories, simulation of the coordinated movements in 3D environment and experimental validation of the developed techniques on the factory floor.

To solve these tasks, the remainder of the work is organized as follows. Chapter 2 focuses on the kinematic modeling of a typical redundant robotic system used in the composite lay-up technology and formalization of the related optimal motion planning problem. Chapter 3 is devoted to the development of a new motion planning method for the redundant robotic systems utilized in the considered process. Chapter 4 deals with industrial implementation of the developed method on the factory floor (for manufacturing of a high-pressure vessel). This chapter focuses on the kinematic modeling of a typical redundant robotic system used in the composite lay-up technology and formalization of the related optimal motion planning problem. The considered problem is presented as finding of the time-optimal trajectory in robotic system joint space under specific constraints related to both the robotic system and technological task. These constraints allow taking into account limitations of the robotic system (joint limits, maximum joint velocities/accelerations), practical requirements for the manipulator postures (via singularities and collision constraints) as well as the Cartesian path constraints issued from the composite lay-up process. The main purpose of this chapter is to formalize the considered technological problem and to present it in the form common for the mathematical optimization theory.

CHAPTER 2 ROBOTIC LAY-UP SYSTEM MODEL AND MOTION GENERATION PROBLEM FORMALIZATION

LAY-UP TASK DESCRIPTION AND ITS CREATION IN CAD SYSTEM

Robotic offline programming for automated lay-up process starts from a 3D CAD model of the workpiece that is to be covered with fiber reinforcements. There are various types of CAD files, but most offline programming software packages are capable of converting other types of CAD files to a compatible one [START_REF] Pan | Recent progress on programming methods for industrial robots[END_REF]. For example, in Figure 19, a 3D model is shown in FPM© that is developed by AUTOMATED DYNAMICS ® . This software works in conjunction with SolidWorks. The user can develop the part surface in SolidWorks directly or import the surface geometry as an IGES, STEP, or other CAD translation. Such universal CAD/CAM designs allow users to create the lay-up curve data files and extract position as well as orientation tags directly. The layup task can be also described in CATIA, which is also able to operate with surfaces and augmented lines defining both positions and orientations of the technological tool for the considered task.

Representation of the Lay-up Task in the Form of the Sequence of 4×4 matrices

Let us assume that the desired lay-up curve with respect to workpiece frame, along which the placement head is to be moved, is imported from the composite part CAD system and described by a 3D-augmented line (as the typical task shown in Figure 19). This augmented curve is discretized in n segments and can be expressed as follows

  ; ... , 2 , 1 , n i i i   a p C (13)
where

T i i i i z y x ) , , (  p
defines the Cartesian coordinates of each sampled task point on the lay-up curve and the unit vector

T zi yi xi i a a a ) , , ( ˆ a
defines the normal direction of each task location outside the workpiece surface.

To describe the spatial location of each task point, let us define the displacement along the Cartesian path

1 ... , 2 , 1 ; 1       n i i i i p p p (14)
and derive the unit vector of roller axis by following equation

1 ... , 2 , 1 ; ˆˆ       n i i i i i i p a p a s (15)
Then, the traveling direction of the placement head on each task location is described as the following vector product

1 ... , 2 , 1 ; ˆ    n i i i i a s n (16)
Thus, each task node can be associated with a Cartesian frame in which the frame origin is located at the path point i p ; X-axis is directed along the path and coincides with the path tangent; the Z-axis is along norm direction; the Y-axis perpendicular to the axes X and Z so that these three axes form a right-handed coordinate frame (see Figure 20), i.e.

1 ... , 2 , 1 ; ) (   n i F i task .
The discretized task locations are described by homogeneous transformation matrix from workpiece frame

W F to ) (i task F 1 ... , 2 , 1 ; 1 0 0 0 ˆ4 4 ) (             n i i i i i i task W p a s n T (17)
Using the above definition, the lay-up task can be presented as a set of locations that should be visited sequentially by the robot end-effector. Since each location is already described as a 4×4 homogenous transformation matrix, the considered task is formalized as follows 18), the task locations (also known as the tool locations) are described by the way of 4×4 homogenous transformation matrices. However, it should be mentioned that, to control the industrial manipulator in practice, an equivalent expression of the end-effector location is more frequently used, i.e. 6×1 vector

n i n task W i task W task W task W ... , 2 , 1 ; ... ... ) ( ) ( ) 2 ( ) 1 (       T T T T (18)
T z y x ) , , , , , (   
. These vectors include three position coordinates and three orientation angles. It can be directly computed from the known homogeneous transformation matrices [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF] It has to be stressed that the task is presented uniquely by the sequences of frames and 4×4 matrices, whereas the corresponding 6×1 vectors might be non-unique. To compute the angle values from the matrices, the scalar equations can be extracted. For example, the angle β i can be obtained from the following expressions

3 3 ) , , (  i i i a s n
n i i i i i i i i task W ... , 2 , 1 ; 1 ) ( 1 0 0 0 ˆ4 4 4 4 ) (                      0 p φ R p a s n T (19)
                      1 0 0 0 cos cos
i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i task w z y x                              r T (20) and ) , ( atan2 2 1 , 2 2 1 , 1 1 , 3 r r r i      (21)
In industrial robotics, the positive solution is usually used here, for which [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF]. In the case of

0 cos  i  . And then, i  and i  can be easily obtained if 0 cos  i  (
0 cos  i 
, known as the singular cases, it leads to infinite number of solutions. However, in practice, two possible conventions are usually recommended to ensure unique solution [START_REF] Nof | Handbook of industrial robotics[END_REF][START_REF] Niku | Introduction to robotics: analysis, systems, applications[END_REF][START_REF] Groover | Automation, production systems, and computer-integrated manufacturing[END_REF].

 Utilizing previously-used i  and i  ;  Using the present i i    and the previous 1 1    i i   . For example, 2 ) ( ) ( : ; 2 ) ( ) ( : 1 1 1 1 1 1                 i i i i i i i i i i i i             (22)
By this way, an equivalent expression of Equation ( 18) is obtained as a unique sequence of 6×1 vectors

i i i i i i ... , 2 , 1 ; 2 2 2 2 2 2 1 1 1 1 1 1                                                                                                     (23)
Based on this one-to-one transformation between two representations, both the joint trajectory and Cartesian trajectory can be obtained during the motion generation for automated lay-up application.

ROBOTIC LAY-UP SYSTEM COMPONENTS AND THEIR MODELS

The considered robotic system includes two principle mechanisms: an industrial robot (a technological tool manipulator) and a positioner (a workpiece manipulator). It is assumed that their kinematics is known and spatial locations of the output frames depend on the vectors of actuated joint coordinates. Let us derive the kinematic models of these two mechanisms and then aggregate them with the model of the technological task.

Kinematic Model of Robotic Manipulator

The kinematic model of robotic manipulator defines relations between the actuated variables (joint coordinates) and spatial location of the technological tool (its position and orientation). In the frame of this work, the links of the robot are modeled as rigid bodies and the joints are assumed to provide pure rotation or translation.

In order to derive the kinematic model of an industrial robot, a universal DH-technique proposed by Denavit-Hartenberg [START_REF] Hartenberg | A kinematic notation for lower pair mechanisms based on matrices[END_REF]) can be applied. However, for industrial robots with standard architectures considered in this work, a more convenient description can be used in order to simplify equations and reduce the computing time required for the optimal motion generation algorithm. To obtain the desired model, let us define a number of frames attached to each link as shown in Figure 21:

Rbase F
Robot base frame. This frame is attached to the robot base and does not move. It can be considered the reference frame where the positions of all other link frames are described. In

X when 0 1  R q . 2 F
Frame 2 is attached to the second link. The origin locates on the rotational axis of 2 R q ; the Xdirection coincides with link 2 L ; and the Y-direction is the same to that in

1 F . 3 F
Frame 3 is attached to the third link. The origin locates on the rotational axis of 3 R q ; the Xdirection is parallel to link 3 L ; and the Y-direction is the same to that in 2 F .

4

F Frame 4 origin is located at the wrist center since the axes 4, 5, 6 all intersect at the wrist center point and are mutually orthogonal. X-direction is parallel to link 3 L , and Y as well as Z have the same directions to that in 3 F when

0 4  R q . 5 F
Frame 5 origin is also located at the wrist center. Its Y-direction is along the rotational axis 5. When

0 5  R q , it completely coincides to 4 F . 6 F
Frame 6 is assigned to the wrist center point too. The X-direction is along the rotational axis 6. When

0 6  R q
, it completely coincides to 5 F .

Rtool F

Robot tool frame has its origin at the tip of the tool (tool center point). The orientation is selected in such a way that its Z axis is identical to the tool axis direction and points out of the tool; its Y axis is parallel to In robot modeling, there are two basic types of models that are usually referred to as direct and inverse ones. The robot direct kinematic model describes the spatial location of mounting flange frame 6 F with respect to Rbase F as a function of the actuated joint coordinates. For the serial architecture presented above, the robot direct kinematic model can be expressed as a product of the 4×4 homogenous transformation matrices [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF])

Rtool R R R R R R Rbase R R q q q q q q g T T T T T T T q 6 6 6 5 5 5 4 4 4 3 3 3 2 2 2 1 1 1 ) ( ) ( ) ( ) ( ) ( ) ( ) (        (25)
that depends on the joint variables

6 2 1 , , , R R R q q q 
. It should be mentioned that it is assumed here that the robot base frame is located at the intersection of the axis 1 Z and the floor level, but it can be easily relocated in practice using the relevant commands of robot programming languages (e.g. $ROBROOT in KRL language). Besides, it is also assumed that the last frame is located to the TCP point (so called "tool center point") and orientation of the tool axes

tool tool tool Z Y X , ,
is determined with the accordance of the technological requirements. This definition corresponds to general type of the 4x4 homogeneous matrix

                    1 0 0 0 cos cos
Z C B C B B Y C A C B A C A C B A B A X C A C B A C A C B A B A T (26)
that depends on six parameters 21) is usually included in the tool parameterization in order to avoid useless computations. In practice, these six parameters can be easily modified using dedicated commands of robot programming languages (e.g. $TOOL in KRL language).

For the remaining matrices, relevant expressions can be obtained in a conventional way, using the Denavit-Hartenberg technique where they are presented as a composition of elementary rotations and translations

) ( ) ( ) ( ) ( ) ( Ri 1 1 1 i Z Z i X i X Ri i i d q a q D R D R T        
depending on both the joint variables Ri q and the parameters

i i i d a  , ,
describing the links/joints geometry. However, in this work, with the slight different definition of the frames, it is possible to achieve some simplification and to avoid unnecessary matrix multiplications. In particular, for the first joint that provides rotation about the vertical axis 1 Z , the matrix

) ( 1 1 R Rbase q T
can be expressed as

                1 0 0 0 1 0 0 0 0 cos sin 0 0 sin cos 0 1 1 1 1 1 L q q q q R R R R Rbase T
where the parameter 0 L takes into account the shift between the robot base frame Rbase F and the frame of the first link 1 F . Similarly, the second and the third joints can be described using the rotations around the horizontal axes 2 Y and 3 Y , which yields the following expressions for the matrices

) ( 2 2 1 R q T and ) ( 3 3 2 R q T                 1 0 0 0 0 cos 0 sin 0 0 1 0 sin 0 cos 2 2 1 2 2 2 1 R R R R q q L q q T                 1 0 0 0 0 cos 0 sin 0 0 1 0 sin 0 cos 3 3 2 3 3 3 2 R R R R q q L q q T
where the parameters For the remaining joint angles providing the tool orientation, the desired matrices can be derived in similar way. In particular, the matrix

) ( 4 4 3 R q T
can be expressed using rotation about the axis

5 R R R R q q q q T
It is clear that the above expressions do not include translational parameters, because the origins of the frames 4, 5, 6 are located at the same points (so-called "robot wrist center"). It should be mentioned that the parameter 4 L defining the distance between the mounting flange and the wrist center is usually included in the matrix describing the tool geometry, i.e. the (1,4) element of Rtool T 6 is assumed to be equal to 4 L X t  . The latter also allows us to reduce computational efforts.

For further convenience, the matrix computations included in the direct kinematic model (24) were executed analytically which allowed us to present the product

) ( ) ( ) ( 6 6 5 2 2 1 1 1 R R R Rbase q q q T T T   in the following way          1 0 0 0 ) ( 6 p a s n q T R Rbase (27) where T R R R R R R R q q q q q q ) , , , , , ( 6 5 4 3 2 1  q
is the vector of the robot joint variables and

                
C S C S C S C S S C S C S C S C C S S S C C C n                       ) (
C C C S S S C S C S S S C C C C C S S S C S C C S S S C S S C C C S S C C S S C S s                       
) ( ) ) ( ( ) ( ) ) ( ( ) ( C S S C C C S S C C S S S S C C S C C S S S C C S C C S S S C S S C C S S C S C S a                              2 2 23 23 3 0 2 2 23 23 3 1 1 2 2 23 23 3 1 1 ) ( ) ( S L C d S L L C L S d C L L S C L S d C L L C p .
Here, usual robotic notations are used allowing achieving compact presentation: To generate desired robot motions, the developed algorithm (see Chapter 3) extensively utilizes the inverse kinematic transformations allowing us to compute the vector of the actuated joint variables R q corresponding to the desired robot flange location described by a given homogeneous 4×4 matrix 6 T Rbase . In robotics literature, this transformation is usually called the robot inverse kinematic model, and it will be denoted below as

Ri i Ri i q S q C sin ; cos   and ) ( sin ); ( cos 3 2 23 3 2 23 R R R R q q S q q C     .
(.) 1  R g
. Compared to the direct transformation, the inverse one is not trivial since the solution is usually not unique and can be expressed in a closed form in some cases only. For this reason, special types of manipulator architecture are used in industry to satisfy the Pieper condition [START_REF] Peiper | The kinematics of manipulators under computer control[END_REF]) that ensures the closed-form solution. Otherwise, there are also some numerical techniques to deal with this problem [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][START_REF] Manocha | Efficient inverse kinematics for general 6R manipulators[END_REF] that may produce up to 16 different solutions for the inverse kinematics. In this work, the serial robot with last three intersecting axes is employed, which corresponds to the majority of industrial applications. For this architecture, for most of the allowable robot end-effector locations, there are 8 possible solutions as is shown in Figure 22. In industrial robot programming, these solutions are distinguished using so-called "configuration index". From geometric point of view, the manipulator configurations corresponding to the same tool location differ by the shoulder posture (forward/backward), the elbow posture (up/down) and the wrist posture (up/down). Algebraically, these configurations will be described by the binary vector

) , , ( 3 2 1     μ
with the components 1  that will be used below for selection of particular solutions of relevant trigonometrical equations (see Table 10). To take into account the multiplicity of the inverse kinematic solutions, the corresponding inverse function will be denoted below as

) , ( 1 μ T q Rtool Rbase R R g   (28)
where the 4×4 homogeneous matrix

Rtool

Rbase

T

defines the desired end-effector location; the configuration vector

) 1 , 1 , 1 (     μ
describes the desired manipulator posture and R q is the vector of the actuated joint coordinates. It is worth of mentioning that the vector μ is an output of the robot direct kinematic model, which is determined by the coordinates of

5 3 1 , , R R R q q q
with the rules derived below.

Table 10 Definition of the configuration index for serial 6-axis manipulator.

Value

1  2  3  1
The x-value of the intersection of the wrist axes, relative to F 1 , is positive.

(shoulder forward)   3 q (elbow up) 0 5  q (wrist up) -1
The x-value of the intersection of the wrist axes, relative to F 1 , is negative.

(shoulder backward)   3 q (elbow down) 0 5  q (wrist down)
 depends on the size of the offset between axis 3 and axis 4;

0   when 0  d
To obtain the inverse kinematic solution for the considered manipulator (satisfying the Pieper condition), two classic approaches can be applied: algebraic and geometric. The first of them is based on decomposing the original six-dimensional task into several relatively simple plane geometry problems that may be easily treated analytically [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF]. The second approach is based on solving a set of trigonometric equations derived directly from Equation (25) describing the direct kinematics [START_REF] Craig | Introduction to robotics: mechanics and control[END_REF]. There are some special techniques here allowing separating the unknown variables (left-and right-hand side matrix multiplications, etc.). The most efficient of them is based on the matrix expression

1 6 6 ) ( ) (    Rtool Rtool Rbase R Rbase T T q T (29)
which is equivalent to relocating the target frame Rtool Rbase T to the robot wrist center, where the frames

6 5 4 , , F F F
are located. The latter allows us to separate the variables

3 2 1 , , R R R q q q and 6 5 4 , , R R R q q q
. By this way, to simplify the computation in this work, a method combining the geometric approach and algebraic approach together is presented as follows.

Computing the angle for axis #1. The joint variable 1 R q can be easily computed from algebraic equations straightforwardly derived from last column of the matrix expression ( 29) . The above system of trigonometric equations can be transformed into

) ( ) ( 2 2 23 23 3 1 1 2 2 23 23 3 1 1 C L S d C L L S p C L S d C L L C p y x               (30)
  y x p S p C   1 1 ; (31) where ) ( 2 2 23 23 3 1 C L S d C L L        
. This allows us to present the desired expression for 1 R q in the form

) , ( atan2 ) , ( atan2 1 1 1 1 1 x y R p p C S q       (32)
where

) ( 2 2 23 23 3 1 1 C L S d C L L sign         (33)
is the configuration index defining the manipulator posture with respect to the axis #1 (shoulder forward/backward). It is clear that 1  must be an input variable for the inverse kinematics and an output of the direct kinematics.

It should be mentioned that the above equations cannot be applied directly if 0   , i.e. when the wrist center point is located on axis#1. In this case, the angle 1 R q cannot be determined unambiguously and thus can take any value. In practice, it is named as "shoulder position singularity" (see Figure 23a). To avoid chaotic manipulator motions in the neighborhood of this singularity, industrial robot controllers may assign to 1 R q either default or previous value. For example, in KRL language of KUKA, user can chose one of these two options by setting a special system variable: $ SINGUL_POS[1] that is equal to 0 for the default setting and is equal to 1 for the previous value. Computing angles for axes #2 and #3. After obtaining the actuated angle for the axis #1, it is possible to solve equations for the wrist center point coordinates

2 2 23 23 3 0 2 2 23 23 3 1 1 2 2 23 23 3 1 1 ) ( ) ( S L C d S L L p C L S d C L L S p C L S d C L L C p z y x                      (34)
with respect to 2 R q and 3 R q assuming that 1 C and 1 S are already known. However, to simplify the final expressions, it is more convenient to apply here the geometric approach.

In order to find the desired angles, let us consider the manipulator projection to the vertical plane defined by the axis that allows us to compute the axis #3 angle (this triangle is highlighted in pink in Figure 24). Applying the law of cosines, one can get the following expression for the adjacent angle  at the vertex connecting the edges 2 L and

2 3 2 L d                  2 2 3 2 2 2 2 2 3 2 1 2 0 2 ) ( ) ( ) ( arccos d L L L d L L L p z   (35)
As follows from the figure, the desired angle 3 R q can be easily obtained from the computed value  and the constant angle

) arctan( 3 L d  
describing geometry of the third link, i.e.

      ) ( 3 R q that yields the expression                          3 2 2 3 2 2 1 2 0 2 2 2 2 3 3 arctan 2 ) ( ) ( ) ( arccos L d d L L L L p L d L q z R  (36)
It should be mentioned that here the value of d is negative, but it may be positive for some industrial robots. In case of 0  d , 3 R q should be rewritten as

      ) ( 3 R q
. Besides, compared to the standard DH convention, the positive directions of 3 R q is opposite compared to the usual one, which leads to slightly different inverse kinematic expressions (this definition of 3 R q corresponds to KUKA robots used for the implementation of the algorithms developed in this work). It is also worth mentioning that the above expression for 3 R q is valid for a particular manipulator configuration that usually is denoted as "ELBOW_UP". To take into account multiplicity of the configurations corresponding to the same location of the wrist center point (see Figure 29), the final expression for angle 3 R q can be presented as

                          3 2 2 3 2 2 1 2 0 2 2 2 2 3 2 3 arctan 2 ) ( ) ( ) ( arccos L d d L L L L p L d L q z R   (37)
where

1 2   
is the second configuration index that defines the manipulator elbow configuration (ELBOW UP/DOWN). Similar to 1  , the configuration index 2  must be computed in the direct kinematics. Analyzing Figure 25, it is easy to prove that the second configuration index is expressed as q , let us consider the second triangle that is highlighted in green in Figure 26, which allows us to compute the auxiliary angle  . As follows from the figure, it can be computed from the expression

) ( 3 2     R q sign (38)
) , ( atan2
1 0 L L p z      (39)
that may provide either positive or negative value. Besides, considering the third triangle (highlighted in blue in Figure 26) and assuming that 3 R q is already known, one can compute the auxiliary angle  using the expression

 

) cos( ), sin( atan2

3 2 3 2 2 3 2 3 2            R R q L d L q L d (40)
Then, using the angles  and  , the desired actuated angle can be expressed as    . However, it is worth mentioning that here the positive direction of 2 R q is opposite to the usual one, which leads to slightly different inverse kinematic expressions (this definition of 2 R q corresponds to KUKA robots). Thus,

) ( 2      R q
here and the final expression for

2 R q is   ) , ( atan2 , atan2 1 0 3 2 3 2 2 3 2 3 2 2 L L p C L d L S L d q z R              (41)
where It should be mentioned that the above equations can be applied only on the condition that the inverse function (.) arccos in equation ( 37) exists and provides real value of  , i.e. if

) sin( 3 3     R q S and ) cos( 3 3     R q C .
1 2 ) ( ) ( ) ( 2 2 3 2 2 2 2 2 3 2 1 2 0         d L L L d L L L p z  (42)
Otherwise, the inverse kinematic transformation should provide an error message "Too Far" or "Too Close", because the desired wrist center point is out of the manipulator working area. It is also worth mentioning that if the above expression is equal to 1  , the robot is at the border of its workspace (see Figure 23b). In this case, the wrist center point is located on the line connecting the centers of the axis #2 and axis #3, so independence of the configuration index . In industrial robotics, the above mentioned phenomenon is categorized as "elbow singularity" and expected to be avoided by assigning the target points inside of the manipulator workspace (so-called reachable locations).

Computing angles for axes #4, #5 and #6. After computing the actuated angles for the axes #1, #2 and #3, it is possible to solve equations for the remaining angles 4 R q , 5 R q and 6 R q . To derive relevant equations, the direct kinematic model ( 25) can be rewritten as

    ) ( ) ( ) ( ) ( ) ( ) ( 6 6 5 5 5 4 4 4 3 1 6 1 3 3 2 2 2 1 1 1 R R R Rtool Rtool Rbase R R R Rbase q q q q q q T T T T T T T T          (43)
where the left hand-side is already known and the right hand-side contains the unknowns to be found. This equation can be also expressed in the scalar form using the following notations

                         1 0 0 0 0 0 0 ) ( ) ( ) ( 1 6 1 3 3 2 2 2 1 1 1 z z z y y y x x x Rtool Rtool Rbase R R R Rbase a s n a s n a s n q q q T T T T T (44) and                       1 0 0 0 0 0 0 ) ,
C C C S S S C C C S S C C C S S C S C S C C S S C S S S C q q q R R R T (45)
For the angle 5 R q , let us consider three scalar equations provided by the first columns of the above matrices

5 4 5 4 5 ; ; S C n S S n C n z y x     (46)
which allow us to find explicitly cosine and sine of 5 R q 2 2 5 5 sin ; cos

z y R x R n n q n q     (47)
and easily compute the desired angle using the function (.) atan2 . However, it is necessary to take into account that two symmetric solutions are possible here, which geometrically correspond to different manipulator wrist configuration (WRIST UP/DOWN) and may be described algebraically by the third configuration index

) ( 5 3 R q sign   (48)
Using this notation, the expression for the angle 5 R q can be presented as

  x z y R n n n q , atan2 2 2 3 5     (49)
where

1 3   
is provided either by the user or by the direct kinematics. It is worth mentioning that there is no computational problems (singularities) related to the angle 5 R q because the situation

  0 , 0 atan2
is not possible due to the identity equation

1 2 2 2    z y x n n n .
For the angles 4 R q and 6 R q , it is easy to write similar scalar equations provided by the first lines and columns of the matrices (45) and (46) 6 R q may be found using expressions

    x x R z y R a s q n n q        3 3 6 3 3 4 , atan2 , atan2     (52)
where the sign of 6 R q . Relevant manipulator configurations correspond to so-called the "wrist singularity" for which

  0 5 R q or    5 R q
(it is clear that the second case is not possible in practice, see Figure 23c). For the feasible case when

  0 5 R q , the matrix (45) is reduced to                      1 0 0 0 0 0 0 0 0 0 0 1 ) ,
C C S S S C C S C S S C S S C C q q R R T (53)
and yields the following scalar equations

) sin( ) cos( 6 4 6 4 6 4 6 4 6 4 6 4 R R z R R y q q S C C S s q q S S C C s         (54)
that provide infinite number of solutions for . In practice, the ambiguity related to the wrist singularity is solved by assigning either the default or previous value to one of the angles, 4 R q or 6 R q . For example, in KRL language of KUKA, user can chose one of these two options by setting a special system variable: $ SINGUL_POS[3] that is equal to 0 for the default setting and is equal to 1 for the previous value. Another solution is sharing the required rotation between axis #4 and axis #6 similar to equation ( 22) in the previous section.

Therefore, in the following sections the considered manipulator can be described using either its direct kinematic transformation

) ( ] [ , R R Rtool Rbase g q μ T 
or inverse kinematic transformation

) ( ] [ , 1 μ T q Rtool Rbase R R g   where ) 1 , 1 , 1 (     μ
the vector of the configuration indices and both of the functions (.) R g and

(.)

1  R g are expressed in closed-form.

Kinematic Model of Actuated Positioner

An actuated positioner is a machine employed to adjust the workpiece posture with respect to the robot in order to ensure accessibility of the given task locations or satisfy some technological requirements. For example, in arc-welding application, the positioner allows the welding tool to be oriented almost vertically while keeping the welding seam nearly horizontally. In automated lay-up process studied in this work, the actuated positioners are also utilized to improve access to the task locations, but they provide some additional opportunities for increasing the lay-up process speed by combining motions of the robot and positioner.

Basically, an actuated positioner is a workpiece hold structure that is driven by the robot controller. The simplest architecture contains a head and a tailstock with a servo drive (so-called one-axis positioner). There are varieties of one-axis positioners that differ in sizes and geometries adopted to the work to be done (see Figure 27). Some smaller types are suitable for light objects while others are capable of handling workpiece weighing many metric tons. For example, KUKA positioner presented in Figure 31 are capable to manipulate objects weighting from 250 kg to 4 tons. The two-axis positioners are usually utilized if it is necessary to orientate the component with respect to the vertical axis. They are frequently used in arc-welding applications since one-axis positioners are not capable of providing the full weld orientation with respect to the gravity. Some examples of two-axis positoners from KUKA are presented in Figure 28. For large-scale workpieces, such as engine cowls, airplane fuselage, etc., the robotic manipulator is mounted on the linear tracks or gantries providing additional translational axis and enlarging the robot workspace. Some examples of these equipments from KUKA are presented in Figure 29. It is clear that physically these mechanical components should be included in the model of the robotic manipulator leading to increase of the actuated axis number up to seven (and obviously causing kinematic redundancy). However, in the frame of this work, the linear tracks and gantries are modeled separately, similar to the workpiece positioners, allowing considering the robot as a non-redundant mechanism with finite number of the inverse kinematic solutions. For the technological process studied in this work (automated lay-up), the workpieces usually are shaped as solids of revolution. Typical examples include pressure vessels, engine cowls, airplane components, etc. For this reason, a one-axis positioner with horizontal rotation is usually sufficient for this application. Here, the positioner adjusts the orientation of the product in order to achieve the desired tool orientation and also ensures the workpiece rotation required for the lay-up. In some specific cases where the product is long or large dimensional, the robotic manipulator is mounted on the linear track/gantry while the workpiece is manipulated using a one-axis positioner. The two-axis positoners are not common for the considered technology.

Similar to the robot model presented above, the kinematic model of the positioner defines relations between the actuated variables (joint coordinates) and spatial location of the workpiece mounted on flange (its position and orientation). In the frame of this work, the components of the positioner are modeled as rigid bodies and the joints are assumed to provide pure rotation or translation. To obtain the desired model, let us define a number of frames as shown in Figure 30: is h . The X-axis of this frame has the same direction to that of Pbase F and their Z directions coincide when

0 1  P q . PF F
Positioner flange frame. The origin of this frame is located at the positioner mounting flange center. Its Y-Z plane describes the flange surface and X-direction is along the rotational axis in the model of one-axis positioner (see Figure 30(a)). In the model of two-axis positioner, the X-Y plane describes the flange surface and Z-direction is along the rotational axis (see Figure 30(b)).

Lbase

F

Linear unit base frame. This frame is attached to the linear track/gantry base and does not move (see Figure 30(c)). It completely coincides with the robot base frame when

0  L q .
Usually, for computational convenience, the world frame is also located at this position. as a function of the actuated coordinates. For the one-axis positioner presented in Figure 34a, the direct kinematic model can be expressed as a product of the following 4×4 homogenous transformation matrices

W PF P PF Pbase P P q q g T T   ) ( ) ( (56)
Z C B C B B Y C A C B A C A C B A B A X C A C B A C A C B A B A T (57)
For the two-axis positioner, the direct kinematic model is expressed as the product of three 4×4 homogenous transformation matrices

W PF P PF M P M Pbase P P q q g T T T q    ) ( ) ( ) ( 2 1 (58)
where . Similarly to the above case, W PF T is the constant matrix describing the workpiece location with respect to the mounting flange.

                1 0 0 0 cos sin 0 0 sin cos 0 0 0 0 1 ) ( 1 1 1 1 1 h q q q q q P P P P P M PBase T                 1 0 0 0 1 0 0 0 0 cos sin 0 0 sin cos ) (
For the linear track or gantry, the direct kinematic model is expressed similar to the one-axis positioner case, i.e. as a product of two 4×4 homogenous transformation matrices

Rbase LF L LF Lbase L L q q g T T   ) ( ) ( (59)
where

               1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 ) ( L L Rbase Lbase q q T
and Rbase LF

T

is the constant matrix describing the robot base location with respect to the mounting flange of the linear unit.

It is clear that the inverse kinematic model of the positioner cannot be presented in systematic way because of lack of degrees of freedom. However, the motion planning and optimization technique developed below does not need the inverse transformation for the positioner or linear track/gantry (only inverse kinematics of the robotic manipulator is used extensively). Nevertheless, it is worth mentioning that for some other applications, such as arc welding etc., the positioner inverse kinematics may be solved in a reduced form assuming that it is necessary to ensure the desired workpiece orientation with respect to the gravity only [START_REF] Pashkevich | Kinematic aspects of a robotpositioner system in an arc welding application[END_REF].

Integrated Kinematic Model of the Robotic Lay-up System

Using the robot and positioner kinematic models presented above, the lay-up task locations can be described in two ways, with respect to the positioner base frame , an integrated model of the lay-up task can be derived. The latter is considered as a principle constraint for the optimization problem studied in the following sections. This constraint describes coupled motions of the positioner and the robot, which guarantees that all the task points are visited by the robot end-effector (it is clear that the solution is not unique due to the redundant degrees of freedom). It should be also mentioned that the location of the world frame is arbitrary in general cases and depends on the user references. In this work, for computational convenience, it is assumed that the world frame coincides with the base frame of the linear unit

) ( ) ( R R L L Lbase World Rtool World g q g q T T    (60)
Further, to express the task locations, it is necessary to take into account the frame alignment conditions:

 The X-axes of the tool and task frames coincide, i.e. that are described by the following homogeneous matrix

                 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 task Rtool T (61)
The latter allows us to obtain the final expression for the given task locations as functions of actuated coordinates of the robot and the linear unit

n i g q g task Rtool i R R i L L Lbase World i Rtask World ... , 2 , 1 ; ) ( ) ( ) ( ) ( ) (      T q T T (62)
Similarly, from the positioner perspective, the workpiece location with respect to the world frame can be computed as

W PF PF Pbase Pbase World W World T T T T   
, which after substitution the direct kinematic equations for the positioner After equating these two presentations ( 62) and ( 63), the desired constraints associated with the technological task may be presented in the following form n i q g g q g i task

W i P P Pbase World task Rtool i R R i L L Lbase World ... , 2 , 1 ; ) ( ) ( ) ( ) ( ) ( ) ( ) (        T T T q T (64)
that may be also rewritten as

    n i q g q g g task Rtool i task W i P P Pbase World i L L Lbase World i R R ... , 2 , 1 ; ) ( ) ( ) ( 1 ) ( ) ( 1 ) ( ) (          T T T T q (65)
This equation describes the relations between the actuated coordinates of the robot and the coordinates of the positioner and linear track, which ensure implementation of the given path. In the following section, for the computational convenience, the constraints (65) are presented in slightly different form

      n i q g q g g task Rtool i task W i P P Pbase World i L L Lbase World R i R ... , 2 , 1 ; , ) ( ) ( 1 ) ( ) ( 1 ) ( 1 ) (           μ T T T T q (66)
that is based on the robot inverse kinematic transformation (.)

1  R g
, which is not unique and depends on the manipulator configuration index μ .

It should be noted that the obtained equations ( 65) may be also treated as integrated kinematic model of the lay-up task, which produce at each task point 6 independent scalar constraints related to the position and orientation. It is clear that among 16 scalar relations produced by straightforward equating of 4×4 homogeneous matrices, there are 4 identity equations coming from the last line and 9 dependencies caused by properties of the 3×3 orthogonal sub-matrices describing rotations. On the other side, the obtained model includes 8 variables (6 actuated coordinates of the robot and 2 actuated coordinates of the positioner and the linear track). This redundancy causes multiplicity of the robot/positioner motions implementing the given technological task. Hence, it provides us some space for the motion optimization in the considered robotic workcell.

MOTION GENERATION IN ROBOTIC LAY-UP SYSTEM

Motion generation in robotic system is defined as the planning of robot/positioner motions while taking into account practical constraints [START_REF] Brock | Generating robot motion: The integration of planning and execution[END_REF]. It includes generation of the time profile for Cartesian position/orientation or the time profiles for all actuated coordinates, together with relevant velocities and accelerations. General scheme for this process is presented in Figure 32. It deals with creating inputs for the robot controller allowing executing the desired smooth movements. In robotic lay-up system, the desired path is pre-defined in the workpiece coordinate system as the sequence of 4×4 homogenous matrices (18). The motion generator must produce a sequence of successive joint configurations

n i q q i L i P i R ... , 2 , 1 ); , , ( ) ( ) ( ) (  q
for the robot, the positioner and the linear unit, which ensures that the robot end-effector strictly follows the desired path. It is clear that because of kinematic redundancy, there is no unique trajectory in the joint space corresponding to the given task. For this reason, there arises the problem of optimal redundancy resolution that is addressed in this work.

Formalization of Motion Planning Problem for Robotic Lay-up System

To utilize the redundancy in the best way for the considered lay-up application, it is reasonable to partition the desired motion among the robotic manipulator, the positioner and the linear unit in such a way that the technological tool passes the given path smoothly and as fast as possible. It is obvious that for technological reasons the lay-up speed should be limited, but in practice this limit is usually much higher than the maximum velocity of the end-effector relative to the workpiece that can be achieved in typical industrial robotic cells. For this reason, the principle objective in the motion planning below will be the total processing time (or the tool travelling time along the given path). The secondary objective is the smoothness of the time profiles for all actuated coordinates. It is worth mentioning that such approach to the redundancy resolution essentially differs from the conventional ones where the processing speed is given and the motion planning is targeted at minimizing the total joint displacement or the coordinate ranges associated with the desired path [START_REF] Dolgui | Manipulator motion planning for high-speed robotic laser cutting[END_REF][START_REF] Dolgui | Enhanced Optimisation Techniques for Off-line Programming of Laser Cutting Robots[END_REF][START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF].

To present this problem in a more formal way, let us introduce the functions . Using this notation, the problem can be presented as minimization of the total travelling time

) (t R q , ) (
) ( ), ( ), ( min t q t q t L P R T q  (67)
over the set of continuous functions

) (t R q , ) (t q P , ) (t q L
that satisfy the Cartesian path constraints imposed by equation ( 64) at all considered time instants

n i t q g t g t q g i task W i P P Pbase World task Rtool i R R i L L Lbase World ... , 2 , 1 ; )) ( ( )) ( ( )) ( ( ) (        T T T q T (68)
In addition, some specific constraints describing kinematic and dynamic capacities of the robot, positioner and the linear unit must be also taken into account. Usually, they can be extracted from the manufacturer specifications and presented as the set of the following inequalities

] , [ ) ( max min R R i R t q q q  , ] , [ ) ( max min P P i P q q t q  , ] , [ ) ( max min L L i L q q t q  (69)
describing mechanical constraints (so-called "joint limits") and the actuator capability presented by the maximum allowable speed in the actuating joints

] , [ ) ( max min R R i R t q q q     , ] , [ ) ( max min P P i P q q t q     , ] , [ ) ( max min L L i L q q t q     (70)
and maximum allowable acceleration

] , [ ) ( max min R R i R t q q q        , ] , [ ) ( max min P P i P q q t q        , ] , [ ) ( max min L L i L q q t q        (71)
To give an idea of typical values of the above constraints, Table 11 presents the joint limits, velocities and accelerations for robot KUKA KR210 R3100 that was used at the implementation stage of this work. It should be noted that the maximum accelerations are not usually provided in a straightforward way but they are expressed via the acceleration time (required to achieve maximum speed starting from zero). For the above mentioned robot, the acceleration time is fixed in the controller and it is equal to 0.25 sec. For the actuated positioner and linear track used in our experiments, the maximum velocities are equal to 142°/s and 1.96m/s respectively, while the acceleration time is 0.5 sec. Hence, the considered motion planning problem can be presented as finding eight smooth (bounded with their derivatives) functions

) (t R q , ) (t q P and ) (t q L on ] , 0 [ T t 
describing motions of the robot manipulator, positioner and linear unit. These functions must not only suit some boundary conditions at the initial and final points 0 1  t and T t n  , but must also satisfy the Cartesian path constraints (68) at each intermediate time instant

1 3 2 ... , ,  n t t t .
From general point of view, this problem can be categorized as the optimization in function space with free end-time and specific constraints at the boundary and intermediate points. In literature, there are several techniques for the problems of such type. They include classical calculus of variations (with Euler-Lagrange equation), Pontryagin's maximum principle, and Bellman's dynamic programming (continuous-time version) [START_REF] Sasane | Optimization in Function Spaces[END_REF][START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]. Classical calculus of variations is a traditional approach for maximizing or minimizing numerical objectives over the function space [START_REF] Gelfand | Calculus of variations[END_REF]. However, this method assumes that the optimization is performed in an open space, i.e. the unknown functions as well as their derivatives are unconstrained. For this reason, it cannot be applied because of the numerous inequality constraints imposed on the unknown functions (joint limits, velocity and acceleration limits).

The second technique, the maximum principle, is widely used in optimal control theory because of its capability of maximizing numerical objectives over the closed function space [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]. The considered problem can be rather easily converted into the relevant form by treating the second derivatives of

) (t R q , ) (t q P and ) (t q L
as the control inputs. The latter corresponds to the control system with double saturated integrators, which sequentially compute velocities and joint variables from the accelerations. For this system, the maximum principle allows rather efficiently solve the time-optimal control problem if the additional constraints are applied at the initial and final points only (classical boundary conditions). However in our case, there are a number of specific path constraints applied at the intermediate points that do not allow easily solving corresponding differential equations and maximizing the Hamiltonian function.

The third technique, Bellman's dynamic programming [START_REF] Bellman | Dynamic programming[END_REF], relies on the principle that if the optimal trajectory is divided into the sub-arcs, then any one of them is also optimal for the corresponding sub-problem. In the case of the continuous-time, this idea leads to the so called Hamilton-Jacobi-Bellman equation in partial derivatives that it is not easy to solve for our technical problem taking into account all above mentioned constraints (joint limits, velocity and acceleration limits, and path constraints). On the other hand, discrete-time formulation with sampling at the points n t t t ... , , 2 1 looks attractive for numerical solution and will be explored in the following chapters. To author's knowledge, at present there is no known technique that can be straightforwardly applied to the problem formulated in this section.

Additional Constraints from the Actual Robotic System

In practice, in addition to the above mentioned constraints (joint limits, velocity and acceleration limits, and path constraints), there are some other important limitations to be taken into account while planning optimal motions in considered robotic cell. First of all, it is necessary to avoid collisions between the manipulator and other workcell components. Besides, it is reasonable to avoid movements in the neighborhood of the manipulator singular configurations, where it is difficult to control the robot motions precisely. Let us present relevant techniques allowing excluding some potential manipulator configurations taking into account these limitations.

In literature, there are several schemes addressing collision detection in path planning of industrial robots. A common one is based on approximating the actual robot geometry by certain simple structures that completely enclose the robot. For example, in some works [START_REF] Chang | Collision-free motion planning for two articulated robot arms using minimum distance functions[END_REF][START_REF] Ennen | Efficient collision avoidance for industrial manipulators with overlapping workspaces[END_REF][START_REF] Liu | Algorithmic safety measures for intelligent industrial co-robots[END_REF], each robot link is replaced by a simple cylinder with the radius equal to the maximum distance between the cylinder axis and the link surface. A collision is detected if the minimal distance between the cylinder axes is smaller than the sum of the related radiuses. This idea is illustrated in Figure 33 where two components are described by the cylinders with the axes

where

(.) dist is the function for calculating minimum Euclidean distances between the cylinder axes. It is clear that this technique should be applied to detect collisions between the robot and positioner, as well as between the technological tool and workpiece, positioner and the robot. However, it is not suitable for detection of the robot self-collisions. In the above expression (72), the function (.) dist can be implemented using the robust technique proposed in Eberly's work [START_REF] Eberly | Robust computation of distance between line segments[END_REF]. To present it, let us parameterize the cylinder axes as

) ( ) ( 1 1 1 n n n n P P P s        and ) ( ) ( 1 2 2 m m m m P P P s        where ] 1 , 0 [ 2 , 1  
. Then the distance between two points on the axes is expressed as follows

) ( ) ( ) ( ) , ( 1 2 1 1 2 1 m m n n m n P P P P P P s                , ( 73
)
which allows us to present the distance between the cylinder axes as the solution of the following quadratic optimization problem over the closed set

] 1 , 0 [ ] 1 , 0 [ ) , ( 2 1     2 1 , 2 2 1 2 1 2 2 2 2 2 1 2 min 2 2 2                        
] 1 , 0 [ ] 1 , 0 [ ) , ( 2 1    
or at an interior point where the gradient is equal to zero, i.e.

0 2 2 2 0 2 2 2 1 2 2 2 2 2 2 1 1 2                        
                 (76) provided that 0 ) ( 2 2 2     m n m n e e e e
. It is clear that the condition

] 1 , 0 [ ] 1 , 0 [ ) , ( 2 1    
must be verified to be sure that the obtained point can be accepted for computing the desired distance.

Otherwise, the minimum of

) , ( 2 1 2   s  occurs on the boundary of ] 1 , 0 [ ] 1 , 0 [  .
Since here the level curves of the objective function are ellipses, the minimum value of 2 s  corresponds to the situation when the relevant ellipse just touches the square edge. Hence, four possible cases

) 0 , ( 1  , ) 1 , ( 1  , ) , 0 (
                    (77)
that should be verified while finding the global minimum on the square

] 1 , 0 [ ] 1 , 0 [  .
In addition, the square corners

) 0 , 0 ( , ) 1 , 0 ( , ) 0 , 1 ( and ) 1 , 1 (
should be also evaluated, which finally leads to 9 separate cases to be considered in this optimization problem.

It should be also noted that equations ( 76) can be applied only when 75) are linear dependent and produce the following set of critical points to be verified and select the minimum of 4 corresponding distances. Therefore, the above presented expressions allow us to compute the distance between the cylinder axes required for the collision detection in the robotic cell.

    1 , 0 ) ( ˆ) ( 1 , 0 ˆ2 2 2 1 2 1 2 1                    m n m m
In order to apply this collision detection technique to the considered task, a proper strategy for selecting pairs of cylinders is also important. Here, because of the existence of the joint limits, it is impossible to have collisions inside the robot (between links) or inside the positioner. Besides, for this task, the collision check between the workpiece and the end-effector is meaningless. So, only the following pairs should be tested: 1) robot forearm/workpiece; 2) robot forearm/positioner shaft; 3) robot forearm/tool part B; 4) tool part A/positioner shaft; 5) tool part B/workpiece; 6) tool part B/positioner shaft. The latter corresponds to the robot description presented in Figure 34 that is based on the approximations listed below:  Robot forearm: approximated by cylinder 1 with center axis Using these notations, the desired collision-free condition may be presented as follows However, it should be mentioned that the above set of inequalities is only a sufficient but not necessary condition for no-collision configurations of the considered robotic cell. Nevertheless, it is rather simple and allows easily reducing the set of admissible candidates in the considered motion planning problem. It is clear that in practice it is necessary to verify the remaining configurations carefully using more precise approximation of the workcell components. To detect possible collisions taking into account real geometry of robot, technological tool, workpiece and positioner, some standard functions provided by commercial CAD/CAM software can be used. In CAD/CAM environments, all components of the robotic cell can be created by 3D modelers and there are several schemes to represent these objects. The simplest ones are the wireframe representation that operates with edges and vertices only. Another scheme is the surface representation where a solid is described by its boundaries (using planar faces, NURBS surfaces, triangular meshes etc.). The most advanced technique widely used in modern CAD system is usually called solid modeling and allows presenting all the object geometry, topology and physical properties in the form of a data structure. It is based on the combination of different methods, such as parameterized primitive instancing, spatial occupancy enumeration, sweeping etc. For example, solids may be defined via finite number of regularized Boolean operations on the primitives or by means of primitive sweeping along a certain space trajectory. However, in this work we will use a simple triangular mesh representation allowing essentially speed-up the collision test required for the motion planning for the considered robotic system. Modern commercial CAD systems usually include integrated functions providing transformations of internal representation of the robotic cell into the set of triangle primitives describing the surfaces of workpiece as well as the components of the robot, technological tool and positioner (see Figure 35). Hence, the desired collision test between workcell components is reduced to numerous checks of triangle-to-triangle intersections.

                                   5 3 2 5 1 5 2 3 1 3 4 3 2 4 1 4 2 3 1 3 5 2 2 5 1 5 2 2 1 2 5 1 2 5 1 5 2 1 1 1 4 1 2 4 1 4 2 1 1 1 3 1 2 3 1 3 2 1 1 1 , , , , , , r r dist r r dist r r dist r r dist r r dist r r dist C C C C C C C C C C C C C C C C C C C C C C
In literature, there are several techniques for checking triangle-to-triangle intersections, and the most common of them was developed by Möller [START_REF] Möller | A Fast Triangle-Triangle Intersection Test[END_REF]. To present this technique, let us describe two triangles T 1 and T 2 by their vertices

1 1 V , 1 2 V , 1 3 V and 2 1 V , 2 2 V , 2 3
V respectively, and denote the planes containing the triangles as S 1 and S 2 . Corresponding equations of these planes may be written in the following form

. 2 , 1 ; 0 1      i i i i V n P n (80) where ) ( ) ( 1 3 1 2 i i i i i V V V V n    
and P denotes the point belonging to the considered plane. Using this notation, one can compute the signed distances from the vertices

1 1 V , 1 2 V , 1 3 V of the first triangle to the plane S 2 : . 3 , 2 , 1 ; 2 1 2 1 2 2      i d i S i V n V n (81)
It is clear that if three computed values 2 S i d are either all positive or all negative, then the triangle T 1 is located on one side of S 2 , which guarantees that there is no intersection between the triangles T 1 and T 2 . Similar test should be applied to the triangle T 2 and the plane S 1. This early test applied at the first step allows us avoiding a lot of computations related to the detailed analysis presented below. ). In this case, the problem is reduced to a simple 2D triangle-triangle overlap test. First, the edges of T 1 should be checked for intersection with the edges of T 2 (or vice versa) that can be easily implemented using the above presented technique computing the distance between the line segments (73). It is clear that an alternative technique can be also applied that is based on 2D computer geometry where the line segments intersection is detected if both segment ends lie on the opposite sides of the corresponding infinite line [START_REF] Ericson | Real-time collision detection[END_REF]. For instance, this approach yields the following inequalities for detection intersection of the edges

1 2 1 1 V V and 2 2 2 1 V V         0 ) ( ) ( 0 ) ( ) ( 2 2 1 1 2 2 2 1 1 1 1 1 1 2 2 1 1 1 2 1              s V V s V V s V V s V V (82)
where

1 2 1 1 ) ( v v v s    , 2 2 1 2 ) ( v v v s    , 1 2 1 1 1 V V v   and 2 2 2 1 2 V V v  
. In case that there is no intersections between the triangle edges, it is necessary to check whether the triangle T 1 is located inside of the triangle T 2 or not (and vice versa). In order to do this, all vertices of T 1 can be checked if they all lie in the interior of T 2 . For example, a condition for the vertex 1 i V to be inside of the triangle T 2 is expressed as

0 0 3 2 2 1      w w w w (83) where ) ( ) ( 1 2 2 1 2 1 1 i i V V V V w     , ) ( ) ( 1 2 3 1 2 2 2 i i V V V V w     and ) ( ) ( 1 2 1 1 2 3 3 i i V V V V w     . If all three 1 1 V , 1 2 V , 1 3
V are inside the T 2 , T 1 is inside of T 2 , i.e. the triangles intersect (similar test should be applied to

2 1 V , 2 2 V , 2 3 V ).
The third step deals with the case when some of the distances 2 S i d have opposite signs. Geometrically, it means that both planes S 1 and S 2 intersect corresponding triangles T 2 and T 1 (Figure 40). It is clear that such intersections are the segments of the common line of the planes S 1 and S 2 , which is directed along the vector 2 1 n n  . If these two segments overlap, the triangles intersect as well. To derive relevant test, let us describe this common line by a parametric equation

      t t
), ( 2 and the common P (see Figure 36) can be computed as

) ( 1 1 1 1 1 1 1 2 2 1 1 1 1 S S S S S S p p d d d p t      (84)
where

) ( ) ( 0 1 2 1 P V n n 1     i S i p . After computing 4 1 ... t t
, the problem is reduced to a simple 1D segmentsegment overlap test. It yields the following inequalities for the intersection detection

        2 1 4 3 4 3 2 1 , min , max , min , max t t t t t t t t    (85)
Hence, the three simple tests presented above allow us to detect intersections between the triangles belonging to the surfaces of the workcell components. However, creation of the relevant sets of the triangles is not a trivial task. In some commercial CAD/CAM software packages, for example CATIA used in this work, the collision detection function can be straightforwardly invoked in "DMU Kinematics" workbench. Using the interface of the "Check Clash" module, the user can define the type of collisions (clearance, contact or clash) to be tested between any selected pair of workcell components (see Figure 37). Then, after activating the sensors on the "Kinematics Simulation" panel, the relevant collision data can be recorded into a ".xls" or ".txt" file. This procedure is implemented by checking all triangle pairs belonging to the selected objects. However, as follows from our experience, using a very fine triangular mesh precisely describing the mechanical components leads to a very high computational time. For example, for the application example considered in Chapter 4, it took more than 3 hours to test collisions for 10 5 manipulator configurations if the triangulation accuracy was set to 0.2 mm (number of the triangles is over 5.3×10 5 , processor Intel ® i5 @ 1.70GHz 2.40GHz, 8G memory). On the other side, setting triangulation accuracy to 4.0 mm yielded about 3.5×10 5 triangles and allowed reducing the computing time down to 2 hours. Hence, the collision test is one of the most timeconsuming procedures in the motion planning process studied in this work and it deserves some special efforts of speeding up. , which is equal to zero if there is no collision detected. This allows us to express the non-collision requirement for the desired motions in the robotic system in a compact form

  ] , 0 [ ; 0 ) ( ), ( , ) ( T t t q t q t cols L P R    q (86)
where L P R q q , , q are the joint variables describing current configurations of the manipulator, positioner and linear unit.

In practice, in addition to collision avoidance, the optimal manipulator motions should be rather far from singular configurations where robot loses some of its degrees of freedom and its kinematic control is difficult. In particular, in the neighborhood of a singular location, an enormous change of the joint coordinates is required for a requested normal change of the robot posture. Generally, the manipulator singular configurations can be identified by studying its Jacobian matrix that becomes rank-deficient. For the considered problem of the optimal motion planning, it is necessary to introduce some numerical indices allowing evaluating how close is the manipulator to a singular configuration.

The simplest quantitative measure of this type (so-called manipulability measure) was proposed by Yoshikawa [START_REF] Yoshikawa | Manipulability of Robotic Mechanisms[END_REF] and is expressed as follows

  ) ( ) ( det R T R w q J q J   (87)
where

) ( R q J
is the manipulator kinematic Jacobian corresponding to the joint variables R q . It is clear that 0  w if the configuration is singular. Geometrically, this measure is proportional to the volume of the manipulability ellipsoid that represents an ability of manipulator to move its end-effector in various directions. However, this index can be hardly used for the considered problem because it provides false evaluation of the manipulator capabilities in the case when the ellipsoid is stretched, i.e. some axes are much longer than other ones but the volume is high enough. Kinematically, it means that the end-effector can easily move at higher speed in the direction of the ellipsoid major axis, while only lower speed is achievable in direction of the minor axes [START_REF] Huo | The joint-limits and singularity avoidance in robotic welding[END_REF].

For this reason, it is logical to utilize another index that can be induced from the manipulability ellipsoid, which is the Jacobian condition number [START_REF] Salisbury | Articulated Hands:Force Control and Kinematic Issues[END_REF]. This index is defined as the ratio of the maximum and minimum of the ellipsoid axes and can be expressed via the singular values i  of the Jacobian matrix

  i i R cond   min max ) (  q J (88)
Corresponding singular values

0  i 
are obtained from the singular value decomposition (SVD) of matrix

) ( R q J as follows   6 2 1 ,... , diag ; ) (        S V S U q J T R (89)
where U and V are 6 6 orthogonal matrices, and S is a 6 6 diagonal matrix of the singular values. It can be easily proved that if the manipulator is in a singular state, then at least one of its singular value is equal to zero, i.e.

0  i 
. The corresponding column of U , i u is referred to as the singular direction. In this state, the motion along the singular directions is not possible. So, the condition number can be treated as an index of the ellipsoid directional uniformity. When its value is close to 1, the ellipsoid is close to a sphere and the manipulator configuration is evaluated as well-conditioned. If the minor axis of the ellipsoid is close to zero (i.e.

0 min  i 
), the condition number becomes high and this state is evaluated as ill-conditioned. So, the condition number can be used in our research as the distance criterion to avoid manipulator's singularities in the motion planning. Mathematically, relevant constraint can be presented in the form of the following inequality

  ] , 0 [ ; )) ( ( max T t C t cond i R    q J (90)
where max C is the user defined value that is about 5 for the application example studied in this work. It should be mentioned that in this work, in order to take into account the non-homogeneity of

  6 2 1 ,... ,   
, the definition of the condition number was slightly modified. In particular, the ratio ( 88)

was computed separately for the translational  

3 2 1 , ,    and rotational   6 5 4 , ,   
components and the maximum of them was used as the condition number.

Summarizing this section, it is necessary to propose a practical strategy of verifications of the above collision and singularity constraints. To minimize overall computing time required for both the data preparation and motion planning itself, first it is reasonable to verify the configuration candidates for the singularities (see Figure 38). As it was mentioned above, it is a rather fast procedure that allows us reducing the size of the initial set. Further, it is logical to apply the "rough" collision detection test based on the cylinder-cylinder intersections that also does not require high computational expenses. And finally, the "fine" collision test should be applied that is based on the triangle-triangle intersections, which provides precise evaluations of the possible configurations analyzed at the motion planning stage. It is clear that such sequence of tests, where many configurations are already eliminated before applying the time-consuming "fine" collision test, allows us avoiding some unnecessary computations.

Figure 38 Strategy of verifications on the additional constraints Hence, after integrating the additional constraints presented above, the considered optimization problem related to the motion planning can be finally formulated as follows

                                                  n i T t t where t q t q t cols C t cond t t t t q g t g t q g t. s. T that such n i t q t q t find n i L i P i R i R L P R i L P R L P R L P R i L P R L P R L P R i L P R L P R i task W i P P Pbase World task Rtool i R R i L L Lbase World i L i P i R q q J q q q q q q q q q T T T q T q          (91)
and it aims at finding 8 continuous functions, i.e.

) (t R q , ) (t q P and ) (t q L
describing the robot manipulator, positioner and linear unit motions. Solution of this problem will be presented in the next section.

SUMMARY

This chapter is devoted to the kinematic modeling of a typical redundant robotic system used in the composite lay-up technology and formalization of the related optimal motion planning problem. The main contribution is presenting the considered technological problem in a pure mathematical way, as finding of the time-optimal trajectory in the joint space of the robotic system under specific constraints related to the composite lay-up task. The latter include both conventional kinematic constraints (joint limits, maximum joint velocities/accelerations), allowable distances to the singularities and obstacles as well as the Cartesian path constraints issued from the composite lay-up process.

In more details, the contribution of Chapter 2 can be summarized as follows:

(i)

Formal description of the technological task and its representation in the form of the sequence of homogenous transformation matrices (or the location vector sequence) describing desired position and orientation of the robot end-effector for the considered lay-up process.

(ii) Integrated kinematic model of the redundant robotic lay-up system incorporating description both the workcell mechanical components (robotic manipulator, workpiece positioner, workspace extensioner, technological tool) and the desired technological task.

(iii) Formalization of the motion planning problem for robotic lay-up system and presenting it in the form of optimization in the function space under the specific constraints describing the technological task, capabilities of the robotic manipulator (joint limits, maximum velocities and accelerations in joints) and the requirements for the manipulator postures expressed via the singularity and collision constraints.

For the optimization problem derived in this chapter, no technique exists in literature that can be directly applied because this problem is highly non-linear and includes very specific constraints arising from the corresponding technology. For these reasons, a new method dedicated to motion planning in the redundant robotic system for automated lay-up process will be developed in the following chapter. This chapter is devoted to the development of a new motion planning method for the redundant robotic systems utilized in the automated composite lay-up process. It is based on the transforming the original problem into a combinatorial one and applying the dynamic programming principle. First, the task graph is created by discretizing the redundant variables and sequentially applying to all task locations the direct kinematics of the positioner (and workspace extension unit) as well as inverse kinematics of the robot. Then, the developed algorithm based on dynamic programming generates the time-optimal motion taking into account all relevant constraints. To reduce the computing time, the time-optimal motion planning is divided in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with large discretization step. Then, the same technique is applied in the local search space, which is created with smaller step in the neighborhood of the obtained trajectory. Alternatively, the second stage may implement a straightforward smoothing of the redundant variable profiles. Efficiency of them and advantages compared to the conventional techniques are confirmed by several case studies.

CHAPTER 3 MOTION PLANNING FOR ROBOTIC LAY-UP SYSTEM

PROBLEM TRANSFORMATION INTO DISCRETE FORM

To solve the motion planning problem (91) presented in Chapter 2, which deals with optimization in function space with specific constraints at the boundary and intermediate points, it is reasonable to transform it from the continuous form into a discrete one. This transformation allows us to apply some combinatorial optimization techniques in order to find the time-optimal trajectories describing the desired motions.

For the considered robotic system containing a 6-axis robot, 1-axis positioner and a 1-axis linear track/gantry (optional), there are one or two redundant variables with respect to the given lay-up task that is defined as a sequence of the robot end-effector locations in 6-dimensional space. Let us consider these cases separately.

Search Space Discretization for One Redundant Variable

It is clear that for the one redundant variable case, any of the workcell joint coordinates can be treated as the redundant one, but it is more convenient to consider the positioner joint angle P q as the redundant variable. This allows us, after solving the positioner direct kinematics for any given P q , to apply straightforwardly the manipulator inverse kinematic function

) , ( 1 μ T Rtool Rbase R g 
and find corresponding joint coordinate vectors R q for the robot (which is not unique because of multiplicity of possible configurations defined by the parameter μ , see Section 2.2.1). This approach permits to take into account explicitly the equality constraints presented in (68) and substantially reduce the search space dimension.

To present the problem in a discrete way, let us sample the allowable domain of the redundant variable ] , [ max min P P P q q q  with the step

P q    m k k q q q P P k P ... , 2 , 1 ) 1 ( min ) (       (92) where 1 / ) ( min max     P P P q q q m
.Then, applying sequentially the positioner direct kinematics and the manipulator inverse kinematics in accordance with ( 56) and ( 28), one can get a set of possible configuration states for the robotic system. For the full mapping from the task space to joint space, let us also take into account the configuration index μ that corresponds to the manipulator posture, and specifies the shoulder, elbow and wrist configurations. Using these notations and assumptions, the set of the manipulator configurations corresponding to the given task and sampled redundant variable

m k q k P ... , 2 , 1 , ) (
 can be presented as follows:

        n i m k t q g g t task Rtool i task W i k P P Pbase World Rbase World R i k R ... , 2 , 1 ; ... , 2 , 1 , ) ) ( ( ) ( -1 ) ( ) ( 1 1 ) (          μ T T T T q (93)
where i t specifies the unknown time instant corresponding to the task location denote the positioner direct kinematics and the manipulator inverse kinematics respectively. Therefore, for each task location we can generate a number of configuration states

) , ( i k task L     n i m k n i i k task i task W ... , 2 , 1 ; ... , 2 , 1 ... , 2 , 1 ) , ( ) (     L T (94) where   ) ( , ) ( ) ( ) ( ) , ( i k P i k R i k task t q t q L 
will be further referred to as the task location cell. This procedure is illustrated in Figure 39. Taking into account that the task locations are strictly ordered in time, the original sequence of 18) may be converted into the directed graph presented in Figure 40. It should be noted that some of the configurations generated using ( 93) must be excluded from the graph because of the collision constraints violation or exceeding of the actuator joint limits (69). Besides, from an engineering point of view, it is prudent to avoid some configurations that are very close to the manipulator singular postures. These cases correspond to the "inadmissible nodes" in the graph (see Figure 40), which are not connected to any of the neighbors. It is clear that due to time-irreversibility, the allowable connections between the graph nodes are limited to the subsequent configuration states

) (i task W T described in (
) 1 , ( ) , (      i k task i k task L L
and the edge weights correspond to the minimum travelling time between relevant configurations that are restricted by the maximum actuator velocities and accelerations expressed by the constraints ( 70) and (71).

To generate the above graph, it can be applied specially developed algorithm presented below. It is described in generic programming language, but it has been tested using the MATLAB and C++ environment. The algorithm input is a sequence of 4×4 location matrices

) (i Task corresponding to ) (i task W T .
The upper and lower limits of the redundant variable are defined as max P q and min P q respectively. The discretization density m determines the number of the discrete values of the redundant variable. The algorithm operates with the positioner direct kinematic function by implementing formula (92), and m×n matrix

} ... , 2 , 1 ; ... , 2 , 1 ) ( { P m k n i i k, q   is obtained.
At the second step, the functions

; ) ( { P k i i k, q  
are computed. It should be mentioned that in the algorithm description the configuration index μ is treated as a constant input parameter, while in practice the step ( 2) should be re-executed for each possible μ . The sub-step (2a) checks the solution of the robot inverse kinematics for each

} ; ) ( { R k i i k,   q
and applies the collision test function for each configuration state

} ; ) ( ), , ( { R k i k,i i k q P   q
. Finally, the task graph is generated with nodes

  ) ( ), , ( ) ( R k,i i k q i k, L P q 
. If a collision is detected or the function

(.) 1  R g
returns the "null" denoting that the robot inverse kinematics is not solvable (because of the joint limits violation, etc.), the "null" value is assigned to the current node (it is marked as "inadmissible" one on Figure 40). The "null" value is assigned if the robot inverse kinematic solution does not exist or the configuration is too close to one of the singularities. Using this presentation of the search space, the considered problem can be transformed to the shortest path searching on the directed graph shown in Figure 40, where each column corresponds to the same task location respectively. With this notation, the desired solution can be represented as the sequence of the nodes

} { } { } { ) ( ) 2 ( ) 1 ( 2 1 ,n k task , k task , k task n L L L     (95)
corresponding to the robot and positioner configuration states  

) ( ), ( ) ( ) ( i k P i k R t q t q
, which ensures the shortest travelling time under the velocity/acceleration constraints of the actuator.

In accordance with the actuator constraints, the distance between subsequent graph nodes can be evaluated as the displacement time for the slowest joint

                 max ) ( 1 , ) ( , 6 ,... 0 ) 1 , ( ) , ( 1 1 max , j k i j k i j j i k task i k task q q q dist i i i i  L L (96)
where 0  j corresponds to the positioner joint variable (1) Redundant variable discretization

For n i to 1  do ); 1 /( )) ( ) ( ( : ) ( min max     m i q i q i q P P P For m k to 1  do ); ( ) 1 ( ) ( : ) , ( min i q k i q i k q P P P      (2) Location matrix creation For n i to 1  do For m k to 1  do (a) ) ( )) , ( ( : i Task i k q g T P P P   ; ) ( : P R T Trans T  ; ); , ( : ) , ( 1 μ q R R R T g i k   (b) If ) ) ) , ( ( ( ) 1 )) , ( ), , ( ( ( ) ) , ( ( max C i k cond i k i k q cols null i k R R P R    q q q ; : ) , ( null i k L  else ; } ) , ( ), , ( { : ) , ( i k i k q i k L R P q    ) 1 , ( ) , ( 1 1 1 ,       i k task i k task n i total i i dist T L L (97) that depends on the indices n k k k ... , , 2 1
. To justify the above expressions, let us consider in details the motion commands implemented in typical robot controllers. Usually, they are based on the trapezoidal velocity profiles and coordinated actuation of all joints, as shown in Figure 41. For the point-to-point motions, these profiles are generated to ensure the minimum time movement between two configurations using full capacities of actuators. Corresponding parameters of trapezoid i T and i  are computed using expressions max max max ;

i i i i i i q q q q T         (98)
where i q  is the required joint displacement, max i q  and max i q   are the maximum joint velocity and acceleration respectively. If the displacement is small enough (it is detected using the inequality

i i T   )
, the trapezoid is reduced to a triangle that is also described by the same parameters i T and i  but they are calculated as

i i i i i T q q      ; max   (99)
It is clear that in general case, the minimum displacement times may be different for different manipulator axes. So, to synchronize the movements of all joints, the robot controller adapts the desired motion to the slowest actuator, i.e. it selects the longest trapezoid and scales the remaining ones to ensure equal displacement time for all of them. This procedure is illustrated in Figure 41 where the second trapezoid is selected as the reference and the others are scaled ensuring the required joint displacement. , the synchronized motion of all joint axes is expressed as follows

N i q q t u q t q start i et t i start i i ... , 2 , 1 ) ( ) ( ) ( arg       (100) where ] 1 , 0 [ ) (  t u
is the normalized timing law function defined on the interval

] , 0 [    T t , such that 0 ) 0 (  u , 1 ) (   T u and 0 ) ( ) 0 (     T u u  
. In the trapezoid case, this function is generated as

                ] , [ , 2 ) ( ) ( ) ( ] , [ ), ( ) ( ] , 0 [ , 2 ) ( 2 2      T T t T t a T t v T u T t t v u t at t u (101)
where

T v 1  and ) ( 1  T a  . In triangle case when   T
, the above expression should be written as follows

          ] 2 , [ , 2 ) ( ) ( ) ( ] , 0 [ , 2 ) ( 2 2       t t a t v u t at t u (102)
Therefore, for a single point-to-point displacement with zero velocities at the beginning and the end (so-called start-stop mode), the robot controller generates trajectories with the motion time   T . However for the considered task that includes multiple point-to-point segments, it is reasonable to utilize another functionality provided by most of the robot controllers that is usually called the "continuous motion" mode (or "approximated positioning" in KUKA robot controllers). It is based on the superposition of the velocity profiles for the consecutive segments where the acceleration phase of the current segment is combined with the deceleration phase of the previous one, as shown in Figure 42. In this case, the motion time for each segment is equal to T which justify expressions ( 96) and ( 97) that are used above in the objective function to be minimized in the motion planning for the considered lay-up task. It should be also mentioned that in practice, it is reasonable to avoid triangular profiles that appear if the sampling of the lay-up path is too fine. For such profiles, the motion time should be computed using equations (99) that are not in agreement with ( 96) and ( 97) but here it is not possible to use full capacities of the actuators (such as the maximum speed), so the total motion time is obviously over the minimum value limited by acceleration/velocity constraints.

Figure 42 Superposition of the velocity profiles for continuous motions

It should be also noted that the applied method of edge weights computing for the directed graph presented in Figure 40 takes into account the velocity constraints, but the acceleration constraints must be examined for each considered path separately. This test can be implemented by applying a formula, which is based on the second order approximation of the corresponding functions

) (t R q and ) (t q P on the time interval ] , [ 1 1    i i t t t
. To derive this formula, let us consider a scalar function ) (t q defined at three sequential time instants

1 1 , ,   i i i t t t by its values 1 1 , ,   i i i q q q
and approximate its acceleration by the second order polynomial

2 ) ( ½ ) ( ) ( i i i i i t t f t t f f t q        
. Then, after solution of equations

1 1 ) (    i i q t q , i i q t q  ) ( and 1 1 ) (    i i q t q
with respect to the coefficients

  ) ( ) ( ) ( 2 1 1 1 1 1 i i i i i i i i i i i t t t t t q q t q q f                   where 1     i i i t t t and i i i t t t      1 1
. The latter allows us to present the acceleration constraints on the desired trajectory of the considered robotic system in the following form

      max 1 1 ) ( 1 , ) ( , 1 ) ( , ) ( 1 , 1 1 2 j i i i i k i j k i j i k i j k i j i q t t t t q q t q q t i i i i                     (103)
and where i t  and 1   i t correspond to the edge weights of the above presented graph which were defined as the traveling times between successive task locations, i.e.

 

) 1 ( ) ( 1 ,     ,i k task ,i k task i i i dist t L L and   ) 1 ( ) ( 1 1 ,      ,i k task ,i k task i i i dist t L L .

Search Space Discretization with Two Redundant Variables

For the two redundant variables case, where the linear track/gantry is activated together with the positioner, the above presented methodology can be generalized in the following way. Similarly to the previous case, let us treat the robot joint angles R q as non-redundant variables and the positioner joint angle ] , [ max min P P P q q q  as well as the linear track/gantry coordinate ] , [ max min L L L q q q  as the redundant ones. Their sampling with the steps P q  and L q  respectively produces the following set of the redundant coordinates

    1 2 ( ) ( ) min min 1 1 1 2 2 2
( 1) 1, 2,... ( 1) 1, 2,...

k k P P P L L L q q q k k m q q q k k m              ( 104
)
where

1 / ) ( min max 1     P P P q q q m and 1 / ) ( min max 2     L L L q q q m
. Then, applying sequentially the direct kinematics for the positioner and linear track/gantry as well as the manipulator inverse kinematics in accordance with ( 56), ( 59) and ( 28), one can get a set of possible configuration states for the robotic system. This set of the manipulator configurations corresponding to the given task and sampled redundant variables can be expressed as follows

                               2 1 2 1 -1 ) ( ) ( 1 ) ( 1 ) , ( , 1 , 1 ) , ( ; ... , 2 , 1 , ) ) ( ( )) ( ( ) ( 1 2 2 1 m m k k n i t q g t q g g t task Rtool i task W i k P P Pbase World i k L L Lbase World R i k k R μ T T T T q (105)
where the functions denote the direct kinematics of the linear track/gantry, direct kinematics of the positioner and the inverse kinematics of the manipulator respectively. Therefore, for each task location we can generate a number of configuration states

) , , ( 2 1 i k k task L     n i m k m k n i i k k task i task W ... , 2 , 1 ; ... , 2 , 1 ; ... , 2 , 1 ... , 2 , 1 2 2 1 1 ) , , ( ) ( 2 1      L T ( 106
)
where the arrays

  ) ( ), ( , ) ( ) ( ) ( ) , ( ) , , ( 2 1 2 1 2 1 i k L i k P i k k R i k k task t q t q t q L 
will be further referred to as the task location cells. This procedure is illustrated in Figure 42. Using similar idea as above, the original sequence of the task locations

) (i task W
T may be converted into the 3D directed graph shown in Figure 43. In contrast to the one redundant variable case where each task location corresponds to a column (see Figure 40), here each Cartesian task location is presented by a set of nodes located on the plane

    2 1 2 1 , 1 , 1 ) , ( m m k k  
. In this graph, the desired motion of the workcell components can be presented as a shortest path

} { } { } { ) , ( ) 2 , ( ) 1 , , ( 2 1 2 2 2 1 1 2 1 1 n k , k task k , k task k k task n n L L L     (107)
that sequentially connects the nodes belonging to the neighbor planes (exactly one node for each plane). For this optimization problem, the objective function, the edge weights, and all constraints are similar to ones corresponding to the case of one redundant variable.

To generate the 3D graph, it can be applied specially developed algorithm (Algorithm 1b) presented below. Its input includes is a sequence of 4×4 location matrices

) (i
Task , the discretization densities 1 m and 2 m , and the upper and lower limits of the redundant variables denoted as max P q , min P q and max L q , min L q respectively. The algorithm transforms the task locations by implementing (99), and m 1 ×m 2 ×n matrix

) (i Task
} ... , 2 , 1 ; ... , 2 , 1 ; ... , 2 , 1 ) , ( ), , ( { 2 2 1 1 2 1 L 2 1 P m k m k n i i k , k q i k , k q   
is obtained. Then, at the second step,

; ; ) , ( ), , ( { 2 1 2 1 L 2 1 P k k i i k , k q i k , k q   
are computed. After checking with the inverse kinematic solvability, collision test and the distance to singularities, the task graph is finally generated with nodes

  ) , ( ), , ( , ) , ( ) , ( 2
1 R 2 1 2 1 2 1 i k , k i k , k q i k , k q i k , k L L P q  .
Figure 44 Directed graph describing the motion planning problem with two redundant variables (1) Positioner and linear unit coordinates discretization

2 1 1 2 1 n i m k m k i k , k L    Notations: L q , P q , R q -Robot
For n i to 1  do ); 1 /( )) ( ) ( ( : ) ( min max     n i q i q i q P P P ); 1 /( )) ( ) ( ( : ) ( min max     n i q i q i q L L L For 1 1 to 1 m k  do For 2 2 to 1 m k  do ); ( ) 1 ( ) ( : ) , ( 1 min 2 1 i q k i q i k , k q P P P      ); ( ) 1 ( ) ( : ) , ( 2 min 2 1 i q k i q i k , k q L L L      (2) Location matrix creation For n i to 1  do For 1 1 to 1 m k  do For 2 2 to 1 m k  do (a) ) ( ) ) , ( ( : 2 1 i Task i k , k q g P P P task   T ; P task L R task i k , k q Trans T T   ) ) , ( ( : 2 1 ; ); , ( : ) , ( 1 2 1 μ T q R task R R g i k , k   (b) If ) ) ) , ( ( ( ) 1 ) ) , ( , ) , ( , ) , ( ( ( )
) , ( ( max 2 1 2 1 2 1 2 1 2 1 C i k , k cond i k , k i k , k q i k , k q cols null i k , k R R L P R    q q q ; : ) , ( 2 1 null i k , k L  else ; )} , ( , ) , ( , ) , ( { : ) , ( 2 1 2 1 2 1 2 1 i k , k i k , k q i k , k q i k , k L R L P q 

MOTION PLANNING METHOD BASED ON COMBINATORIAL OPTIMIZATION

Using the graph generation approach presented in Section 3.1, the original optimization problem in function space ( 91) is converted into a discrete one, which is aimed at searching of the shortest path connecting the subsets of the nodes corresponding to the given task locations. These subsets are grouped in columns for the one redundant variable and in planes for the two redundant variables. The path should visit exactly one node in the subsets for each task location.

Straightforward Approaches and Related Difficulties

The classical shortest path problem is to find a path between 2 vertices in an undirected/directed graph such that the total sum of the edges weights is the minimum. It is clear that the generated graph for the considered problem is a directed one because of the time irreversibility. In the simplest case, when all edge weights were equal, this problem could be solved easily using the Breadth-first search. However, for the considered problem the equal-weight assumption is not valid because the robot travelling times between subsequent task locations are usually different. So, the considered problem can be categorized as finding the shortest path in a directed graph with non-equal positive weights.

In literature, there are several approaches that address the general shortest path problem in a directed graph. One of the most common of them is the Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF]. It can be used for finding the shortest paths from a single departure vertex to a single destination vertex. To apply the Dijkstra's algorithm to the considered engineering problem, it is necessary to create two additional virtual vertices to the generated graph at the beginning and the end as shown in Figure 43 where the weights of the new edges are assumed to be equal to each other (Gueta et al., 2009b, Gueta et al., 2008b[START_REF] Dolgui | Enhanced Optimisation Techniques for Off-line Programming of Laser Cutting Robots[END_REF]. Using such technique, the desired solution is presented as the shortest path starting from the virtual beginning vertex to the virtual end one. It is clear that for the considered graph the shortest path should pass through exactly one vertex in each column. Using the Dijkstra's algorithm, it is possible to compute the shortest distance (robot travelling time here) between the initial vertex and any other vertices in a graph. To implement this algorithm, it is initially assumed that a tentative distance value is assigned to every vertex (zero for the initial one and to infinity for the others). First, the initial vertex is considered as the current one, while others are marked as "unvisited". At the following steps, all neighbors of the current vertex are analyzed and their tentative distances are calculated by summing the current vertex distance and the weight of the edge connecting the current vertex with the corresponding neighbor. Then, the smaller value is selected and considered further as the distance of the current vertex, which is marked as a "visited" one. This procedure is repeated until all the vertices have been visited. Generally, the idea of the Dijkstra's algorithm is to pick the unvisited vertex with the lowest distance, calculate the distance through it to each unvisited neighbor and to exclude longer distances when making an update [START_REF] Cormen | Introduction to algorithms[END_REF].

It should be mentioned that the standard Dijkstra's algorithm uses specific structure of the input data, usually in the form "a sparse matrix" (see Matlab function "graphshortestpath" for instance). Relevant transformation of the original data describing our engineering problem into the format required for the Dijkstra's algorithm requires non-negligible amount of computing time. For example, even for the simplest workcell consisting of a 2-axis manipulator and a 1-axis positioner considered in this work as a benchmark example, it takes several hours to transform the data of a 100×360 two dimensional directed graph into the desired format (processor Intel ® i5 1.70 GHz 2.40 GHz, 8G memory). In contrast, the shortest path computation for this graph using the Dijkstra's algorithm requires several seconds only.

Another known technique that solves the shortest path problem is the Bellman Ford's algorithm [START_REF] Bellman | On a routing problem[END_REF][START_REF] Ford Jr | Network flow theory[END_REF]. Similar to the algorithm of Dijkstra, it is also based on the relaxation principle where a tentative distance value is gradually replaced by the minimum of its old value and the length of a newly found path until eventually reaching the optimum solution. However, Dijkstra's algorithm uses a priority queue to select the closest vertex that has not been "visited", and relaxes all of its outgoing edges. By contrast, the Bellman-Ford's algorithm does this relaxation on all the possible edges. In each of these repetitions, the number of vertices with optimum distances grows, and eventually all vertices will have their right distances. However, it brings the difficulties similar to that when applying Dijkstra's algorithm, i.e. too much computational effort for data transformation and preparation.

For the considered engineering problem, motion planning for redundant robotic systems, straightforward application of the known techniques has essential limitations. The first of them has been already mentioned above and is related to the high computational efforts for input data transformation into required structure. The second one comes out when the designer tries to take into account the acceleration constraints expressed by inequalities (71), which are quite important here for technological reasons. In fact, the conventional shortest path algorithms focus on simple minimization of the edge weight sum, while ignoring the smooth motion requirement that comes from robotics and is outside of the classical directed graph with weighted edges.

On the other hand, the data structure describing the discretized motion planning problem includes additional information. In particular, each vertex of the corresponding graph represents the task location containing all joint coordinates of the robotic manipulator, positioner and linear track/gantry that ensure desired positioning of the robot end-effector. These vertex attributes allow us to evaluate the accelerations of all robotic system components for any considered trajectory, which gives some perspectives for enhancement of the shortest path techniques on the directed graphs adapted to the engineering problem studied in this thesis.

Motion Generation using Dynamic Programming Technique

As follows from the previous section, it is not reasonable to apply conventional shortest path algorithms (Dijkstra's, Bellman-Ford's, etc.) to find optimal motions for the considered robotic system because of high computational expenses and neglecting some technical constraints. On the other side, the directed graph related to the engineering problem studied here has very specific matrix structure where the vertices are naturally grouped in planes or columns and the edges connect the neighboring planes/columns only. Besides, the desired optimal path must include a single vertex from each plane/column starting from the first and ending at the last one. These specific properties of the directed graph motivate us to develop a problem-oriented optimization technique presented below.

Because of inherited matrix structure of the data describing robotic lay-up process and an additive objective, it is reasonable to apply here dynamic programming. Generally, it is a method for solving a complex problem by breaking it down into a set of sub-problems of smaller sizes, solving each of those sub-problems only once, and storing their results [START_REF] Bellman | The theory of dynamic programming[END_REF]. Next time the same subproblem appears, instead of computing its solution again, one simply indexes the previously computed solution, thereby saving computation time at the expense of a slight expenditure in storage space. A dynamic programming algorithm will examine the previously solved sub-problems and will combine their solutions to give the best solution for the given problem.

To apply the dynamic programming, the problem must possess two key properties [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]. The first one is called optimal substructure that means the solution to an optimization problem can be obtained by the combination of optimal problems to its sub-problems, and these optimal substructures are able to be formulated in a recursive form. Here, the desired solution of the considered shortest path problem is already represented as the sum of sub-solutions described in (97), which allows us to write it in a type of recursion in future. Secondly, the sub-problems should be overlapping, which means any recursive algorithm solving the problem should solve the same subproblems repeatedly, rather than generating new sub-problems. In the considered problem, the subproblems are always the paths from the initial plane/column to the current one, and it is impossible to generating new sub-problems during the recursion.

To present the developed algorithm based on dynamic programming for the case of one redundant variable, let us denote i k d , as the length of the shortest path connecting one of the initial vertices

} , { 1 ) 1 1 ( k , k task  L to the current vertex } { ) ( i k, task

L

. Then, taking into account the additivity of the objective (97), the shortest path for the nodes corresponding to the next task point

} , { ) 1 ( k i k, task   L
can be found by combining the optimal solutions for the previous column

} , { ) ( k i , k task    L
and the distances between the nodes with the indices i and

1  i
.The latter corresponds to the formula

    ) ( ) 1 , ( , 1 , , min i , k task i k task i k k i k dist d d        L L (108)
that is applied sequentially starting from the second task point, i.e. for

     180 180 P q s s / 50 / 50      P q  2 2 / 200 / 200 s s      P q       180 180 R1 q s s / 10 / 10      R1 q  2 2 / 40 / 40 s s      R1 q       180 180 R2 q s s / 30 / 30      R2 q  2 2 / 120 / 120 s s      R2 q 
As follows from a simulation study, the developed algorithm is able to find the desired time-optimal motions but the computing time is not negligible. In particular, it takes about 10 minutes to find the optimal solution for the simple benchmark problem. It should be stressed that in contrast to the conventional techniques, this algorithm can generate desired trajectories taking into account the acceleration constraints. Nevertheless, since the inequalities (103) are checked for each candidate path during searching, the computation time may be still too high for real life lay-up applications where 6-

m j j r n k k   min arg ; (4) Backtracking For 2 to n i  Set ) ), ( ( : ) 1 ( 0 0 i i k P i k  
axis robots are usually used and the number of the task points is essentially larger. This motivates us to improve the developed algorithm further.

ENHANCEMENT OF THE MOTION PLANNING ALGORITHM

As follows from the simulation study, the motion planning algorithm presented in the previous section should be enhanced to be applied to real industrial problems. In particular, to obtain the desired time-optimal smooth motions the discretization should be very fine, which causes quite high computational efforts. For this reason, this section presents two strategies allowing finding reasonable balance between the algorithm running time and smoothness of the obtained trajectories. In both of them, the basic idea is to apply the rough discretization first in order to generate an initial solution, which is locally improved further using either smaller discretization step or approximation of the redundant variable profiles by polynomial functions.

First Strategy: Progressive reduction of the discretization step

To reduce the computing time for real industrial problems, let us apply the following heuristic technique: solve the problem applying several times the same optimization routine (Algorithm 2) in different search spaces. It is reasonable to start with a rather big discretization step P q  (stage I), and further improve the initial solution in its neighborhood using relatively small P q  (stage II). In this technique, in spite of the repetition of the basic optimization routine, both stages operate with the search space of smaller size compared to the straightforward optimization for small P q  in the fullrange space, which allows us essentially decreasing the amount of the computations. The details of this strategy are presented in Figure 46. The stage I, named here as the global rough search, is based on the full-range discretization of the redundant variable with relatively big increment step (compared to that in the fine search algorithm presented before). By applying the dynamic programming principle expressed by ( 108), an initial solution can be obtained quite rapidly (it takes about one minute for the benchmark example that deals with the two-axis robot and one-axis positioner, 100 task locations). It should be mentioned that the motion obtained at this stage may be not satisfactory from the engineering point of view since the discretization step is large. So, it is reasonable to optimize the obtained solution further.

m k n i ,i k L DP D n i i k           (d) Extract optimal trajectory For n i to 1  do ); ), ( ( ) ( 0 i i k L i OptTraj    q

Notations:

S P q , S R q -Smoothed positioner and robot joint coordinates

 ; order }, 2 1 . ) ( { }, 2 1 . ) ( { } ... , 2 , 1 ) ( { n ... , , i t i InitTraj n ... , , i q i InitTraj AP n i i q P S P     (b) Recompute the robot trajectory using smoothed redundant variable profile For n i to 1  do   ) ( ) ( : i Task i q g T S P P P   ; ) ( : P R T Trans T  ; ); , ( : ) ( 1 μ q R R s R T g i  
Table 13 Comparison of the conventional and developed algorithms using the benchmark example In order to compare the developed algorithms, the positioner joint coordinate P q was treated as the redundant variable and firstly was discretized in the interval

] 180 , 180 [    with the step    0 . 3 P q
. Then, by sequentially applying the positioner direct kinematic and robot inverse kinematic transformations, a column of 121 candidate configuration states was obtained for each of 100 task points. To ensure that no collision and no singular configuration happen, each candidate configuration state was verified. For this operation, the robot links and the path segments were presented as short straight lines, and the intersections of the robot links and path segments were checked. If the intersection was detected, the current configuration state was marked as "inadmissible" one. Besides, similar mark was used for the vertices for which the robot inverse kinematics is insolvable or the manipulator configuration is too close to the singularity. After this procedure, the original sequence of the task points i p was converted into the task graph composed of 12100 vertices containing the joint coordinates of the robot and positioner. The obtained graph was used for the stage I of the optimization algorithm. The computing time for the graph generation was 43 seconds. Similarly, there were generated the task graphs for the smaller discretization steps

   0 . 1 P q and    5 . 0 P q
. In the first case, the task graph was composed of 36100 vertices and it took around 2 minutes to generate it. For the second case, the task was converted into the graph composed of 72100 vertices, it took about 5 minutes.

Using the obtained graphs, a conventional technique based on the Dijkstra's shortest path was applied first to find the time-optimal trajectory. In the case of

   0 . 3 P q
, it took 3.2 minutes to go through the graph generation, data transformation and path planning in order to find a solution with the robot motion time T motion = 7.63 sec. By reducing the discretization step down to

   0 . 1 P q
, a better solution with T motion = 6.15 sec was obtained but it required more than 20 hours of computations. Further, for the discretization step of

   5 . 0 P q
, the conventional technique became extremely slow and did not allow us to complete the optimization in acceptable time. In fact, a rough estimation of T comp. based on polynomial approximation yields the computing time more than 30 days, which is not acceptable in practice. Nevertheless, relevant time-optimal solution with T motion = 6.00 sec for this case was found using the developed one-stage DP based algorithm with released acceleration constraints, which obviously must give the same result as the Dijkstra's shortest path. Hence, the conventional technique, which does not take into account the particularities of the task graph, cannot be applied to the considered engineering problem because of excessive computational expenses. Besides, the trajectory smoothness requirement is also not respected here.

It should be noted that the above estimation of the computing time is not in good agreement with the classical formula

) log ( V V E O  
describing the complexity of the Dijkstra's algorithm via the number of the graph edges E and vertices V [START_REF] Fredman | Fibonacci heaps and their uses in improved network optimization algorithms[END_REF]. This issue is caused by a number of additional computations included in the path planning algorithm that deal with the data transformation and computing the distances between the nodes, which are in fact the travelling times between the manipulator and the positioner configurations corresponding to the considered vertices.

In contrast to the conventional Dijkstra's shortest path technique, the developed one-stage DP algorithm is essentially faster and allows user to take into account the acceleration constraints that ensure smoothness of the robot and positioner motion. In particular, the computing time for the discretization

       5 . 0 , 0 . 1 , 0 . 3 P q
is equal to 2.4, 15.5 and 59.9 minutes respectively. It is understandable that the optimal motion time should be higher here compared to the previous column where the acceleration constraints are omitted. While reducing the discretization step, the optimal motion time gradually decreases from 9.44 to 8.20 and 7.51 seconds. It is worth mentioning that further reduction of the discretization step is not reasonable because it yields almost the same optimal motion time of the robotic system (the minimum value of T motion obtained in computational experiments for smaller P q  is equal to 7.49 sec). Besides, for    5 . 0 P q this algorithm generates smooth motions that satisfies the engineering requirements. Nevertheless, even for this simple benchmark example, the computing time is still too high if the discretization step is small enough, which limits the practical applicability of the developed one-stage DP algorithm.

The first version of the enhanced two-stage DP algorithm which combines the initial search with the rough discretization and the local optimization with the fine discretization provided a slightly higher value of the optimal motion time which is equal to 7.63 sec. The latter is understandable because the local search space is limited here and it does not obligatory include the global optimal solution. Nevertheless, such small increase of the motion time is not very essential for practice while the profit from the reduction of computational expenses is very high. In particular, the total computing time is about 5.2 minutes only, for executing the first stage with . It should be also mentioned that the smoothness of the obtained trajectory satisfies the engineering requirements.

The second version of the enhanced two-stage DP algorithm combining the rough search, local search and the trajectory smoothing yielded higher value of the motion time 8.47 sec but it is faster than the previous one. In particular, for this technique the total computing time is about 3.8 min and the trajectory smoothness is very good. It is clear that in practice, the user can repeat the local optimization several times and then apply the smoothing, which allows to achieve the same motion time as above. Generally, the first and second versions of the enhanced two-stage DP algorithms are complimentary and may be applied in different variations depending on the particular engineering problem.

Advantages of the proposed motion planning techniques are illustrated by Figures 49-52 that present the trajectories generated using the enhanced versions of the developed motion planning algorithms. The first group of them, which includes Figures 49 and50, is related to the version #1 and demonstrates benefit of the local optimization. It also shows improvement of the optimization results with reduction of the discretization step, which gradually decreases from 3.0° to 0.5°. As follows from Figure 49, the first stage that is running with relatively large discretization step

   0 . 3 P q
yields quite smooth motions but the robotic system motion time is not achieved its minimal value. In contrast, Figure 50 showing results of the local optimization for

   5 . 0 P q
with the range  30 confirms benefit of the proposed approach, which for this example allowed reducing the motion time by 20%. However, the trajectory smoothness drops down here. So further reduction of P q  is needed but it obviously requires additional computing time. The second group of illustrations, which includes Figures 51 and52, is related to the version #2 and demonstrates the gain from the smoothing of the redundant variable profile. Here, the local optimization was stopped with the larger discretization step

   0 . 1 P q
that does not provide the satisfactory smoothness, see Figure 51. Further, the function

) (t q P
describing the positioner motion is smoothed by applying the cubic spline approximation without changing the motion time. The robotic manipulator motions are also adjusted using the kinematic equations. It allows us to obtain very smooth trajectories, see Figure 52, which are obviously suboptimal with respect to the motion time. It is clear that the similar smoothing technique can be also integrated into the first version allowing achieving smaller motion time if it is necessary from the engineering point of view. It is also worth mentioning that in some cases the smoothing can cause some small violations of the velocity/acceleration constraints but it was not observed for the considered benchmark example. Therefore, the presented simulation results confirm the efficiency of the developed enhanced algorithms for the motion planning in the redundant robotic system. They can be applied in several ways, combining stages of global/local optimization and smoothing of the redundant variable profiles. While generating the robotic system motions, the user can find reasonable balance between the computing expenses and achieved level of the key design objectives, which are the total motion time and the trajectory smoothness. Another conclusion following from the simulation is related to the selection of the discretization step, which is an issue to be investigated in more details. As follows from the presented case study, a smaller discretization step may provide faster and smoother motions but relevant computation cost is extremely high. On the other side, excessive discretization does not yield significant improvement of the performance index (the robotic system motion time). For instance, for the considered example, reducing the discretization step from

   5 . 0 P q to    3 . 0 P q
gives almost negligible improvement of the motion time, by 0.4% only (tested using the one-stage DP algorithm). At the same time, it requires 3.9 hours of computation instead of 55.4 minutes for

   5 . 0 P q
. Hence, it is reasonable to provide the user with some simple rules for selecting reasonable discretization step in order to avoid too excessive computations.

PARAMETERS TUNING FOR MOTION PLANNNING ALGORITHM

As follows from the previous section, the developed algorithms have some essential advantages such as high computational efficiency and capability to take into account both of the velocity and acceleration constraints. But proper selection of discretization step for the redundant variable is not trivial. An inappropriate discretization step may lead either to a bad solution or to high computational effort. For this reason, some recommendations for selecting the discretization step are presented below to tune the optimization algorithm.

Influence of Discretization Step on the Generated Trajectories

As it was shown before, a smaller discretization step is better for the considered problem and it allows generating faster and smoother motions because of better approximation of the original continuous problem. On the other hand, too excessive discretization may not yield significant improvement of the primary objective, which is the robotic system motion time here, while taking too much computing time. To investigate the influence of discretization step on the generated trajectories, a slightly different benchmark example is considered here dealing with a straight line lay-up task presented in Figure 53. For this task, the desired linear path was uniformly discretized and replaced by 40 discrete points . These points should be visited sequentially by the robot end-effector in minimum time by using in the best way the motion capabilities of both the robot and the positioner. using the above presented one-stage algorithm, which requires more computations compared to the two-stage ones but obviously allows us to obtain the global optimal solution. Relevant results are presented in Table 14 and in Figure 54, which show that decreasing of the discretization step P q  from 2.00° to 0.50° allows achieving the motion time reduction by about 30%. However, further decreasing of P q  down to 0.25° and 0.10° yields negligible influence on the motion time that is the primary optimization objective in the developed algorithms. Hence, for this case study, the optimal discretization step is about 0.5°.

There are also some interesting phenomena that were observed in the optimization results for relatively large discretization step P q  . In particular, in the case of

   2 P q
, the optimization algorithm generates a technically unreasonable solution that does not take advantage of the positioner motion capabilities. For this solution, the desired motion is executed by the robot only while the positioner stands still (locked), i.e.

var const q R P   q ,
. The reason for this is that the discretization step P q  here is so high that time required for the positioner rotation from P q to P P q q   is always larger than the robot alone moving time between the subsequent task points. Another specific phenomenon can be observed for slightly smaller discretization step

   1 P q
, where the algorithm produces a solution with non-smooth intermittent positioner rotation (start-stop motion), i.e. ... ), 2 ( ... ), ( ), ( ... , , P P P P P P P P q q q q q q q q      

. It is clear that both of these solutions are not acceptable in practice because they do not use all advantages of the robotic system redundancy. Hence, the discretization step for the redundant variable P q  should be small enough to avoid both the positioner locking and its start-stop motions. In addition, it was noted that the discretization step reduction does not always yield monotonic decreasing of the robotic system motion time, while it is intuitively follows from the physical nature of the considered motion planning problem. For example, the discretization step reduction from

   2 P q to    1 P q
leads to even worse solution, where the system motion time is about 10% higher (case with acceleration constraints). This non-monotonic phenomenon can be explained by heuristic integration of the acceleration constraints into the optimization algorithm, which may slightly violate the dynamic programming principle. In more details it is explained in Figure 55 where it is presented a portion of the task graph corresponding to 17 th -19 th task locations. For the convenience, the configuration states are denoted by capital letters, and the numbers in each cell indicates the length of the shortest path to the current one. For instance, in the cell F, 1.1334 is the shortest path from previous C to the current one, and the distance between C and F is 0.0667. In the search process, the acceleration constraint is applied on the candidate sub-path composed of three configuration states. For example, here there are three possible ways to J, which are B-E-J, C-F-J and D-G-J, and all of them are checked for the acceleration constraint. Then, the shortest path satisfying the acceleration constraint is recorded in J. In this example, the recorded path is D-G-J with the length 1.6737, since shorter options C-F-J and B-E-J with the lengths 1.2004 and 1.2111 violate the acceleration limits. However, we can notice that another possible path B-F-J with the length 1.2020 was not detected by the algorithm because the recorded previous cell of F is C since C-F is obviously shorter than B-F. So, because of the heuristic integration of acceleration constraints in the DP-based algorithms, some better trajectories may be missed, which explains increasing of the objective function while reducing the discretization step P q  from 2° down to 1°. However, for smaller P q  this phenomenon is not observed. Hence, to apply the developed technique in practice, users need some simple "rules of thumb" that allows setting an initial value of P q  . Then, the optimization algorithm can be applied several times (sequentially reducing P q  ) until the objective function convergence.

Practical Recommendations for Selecting the Discretization Step

To find a reasonable initial value of the discretization step, let us investigate in details robot and positioner motions between two sequential task locations. It is clear that for smooth positioner motions, it is required that corresponding increments of the coordinate P q should include at least several discretization steps. To estimate the maximum value of P q  , let us consider the movement between two adjacent task locations (P i -P i+1 ) and denote corresponding increment of the positioner coordinate as   and the length of the path segment as s  .

Since the workpiece is rotated along the positioner axis, as small segment of processing path can be approximated by a circle arc (see Figure 56 . Moreover, to ensure optimal coordinated motions of the positioner and the robot, it should be satisfied the equality

max max max P R q v r s          
where max R v is the maximum Cartesian speed of the robot end-effector. The latter gives us a rough estimation of the positioner increment

max max max P R q v r s      (110)
which allows setting an initial value of the discretization step

k q P    
where k should be higher than one to avoid the positioner locking or the start-stop motions. For example, for the above benchmark example, this expression gives the increment value . Hence, such simple "rules of thumb" based on the equation ( 110) helps the user to select an initial value of P q  . Then, the optimization algorithm is applied iteratively with sequentially reduced discretization step until the robotic system motion time convergence.

Figure 56 A sample of the lay-up path and its execution by the robot and the positioner Another issue to be discussed here is related to the sampling of the lay-up path (i.e. the task discretization). Since the original task is defined as a continuous Cartesian curve, it seems reasonable to describe it using relatively high number of discrete points. However, as follows from our study, the excessive sampling of the path leads to essential computational efforts while providing just slight benefits for the fiber placement quality, since relevant technology is not extremely sensitive to the positional accuracy of the robot end-effector. In particular, some positional errors can be compensated by the end-effector that provides pressure for the fiber compacting. Besides, this technology requires some overlaying of the successful fiber coats, and the overlaying degree is not very critical here. From our experience, the lay-up path discretization with 100 to 200 nodes is quite sufficient for the considered examples. But this number should be certainly increased for complicated objects that are presented in the Chapter 4.

SUMMARY

This chapter is devoted to the development of a new motion planning method for the redundant robotic systems utilized in the automated composite lay-up processes. The main contributions are in the area of redundancy resolution via optimization techniques. The proposed method is based on the transforming the original problem into a combinatorial one and applying the dynamic programming principle. To reduce the computing time, the time-optimal motion planning is divided in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with large discretization step. Then, the same technique is applied in the local search space, which is created with smaller step in the neighborhood of the obtained trajectory. Alternatively, the second stage may implement a straightforward smoothing of the redundant variable profiles. Efficiency of them and advantages compared to the conventional techniques are confirmed by several case studies.

In more details, contributions of Chapter 3 can be summarized as follows (i)

Transformation of the original continuous optimization problem into a combinatorial one by discretizing the redundant variables (actuated coordinates of the positioner and workspace extension unit) and sequentially applying to all task locations the direct kinematics of the positioner/extensioner as well as inverse kinematics of the robot.

(ii) Development of a new motion planning method which allows minimizing the total motion time and ensuring smooth movements of all mechanical components. This method is based on dynamic programming and, in contrast to the conventional techniques, it allows taking into account the acceleration constraints related to the trajectory smoothness.

(iii) Development of the enhancement strategies for the developed motion planning technique.

To reduce the computing time required for the motion planning, the time-optimal motion planning is divided in two stages combining global and local searches. As an alternative, the second stage may implement a straightforward smoothing of the redundant variable profiles based on the spline approximation. The efficiency of these strategies was confirmed by case studies.

(iv) Development of the practical recommendation for the parameters tuning in the proposed motion planning algorithm, such as selection of the discretization step for the redundant variables and size of the local optimization sub-space.

The motion planning method developed in this chapter allows generating smooth time-optimal trajectories required for the composite lay-up technology. Its advantages are confirmed by an industrial case study presented in the following chapter.

The main results of Chapter 3 have been presented in conferences EUCOMES'2016 and ICMIT'2017, they also have been published in the journal Mechanism and Machine Theory (2017) (Gao et al., 2017c, Gao et al., 2017a, Gao et al., 2017b). This chapter deals with industrial implementation of the developed motion planning method on the factory floor (for manufacturing of a high-pressure vessel). Using the proposed method, there were generated time-optimal smooth motions for the manipulator and positioner of the SPIDE-TP robotic platform allowing speeding up the thermoplastic tape lay-up process. Correctness of these motions was verified in 3D simulation environment of CATIA software package. Besides, it was created a program in KRL language implementing these motions and ensuring coordinated control of KUKA KR210 robot and AFPT workpiece positioner. There are presented results of industrial experiments carried out with this program, which confirmed smoothness of the manipulator/positioner movements and essential reduction of the time required for the thermoplastic tape lay-up process.

MANUFACTURING PROCESS PREPARATION FOR HIGH-PRESSURE VESSEL FABRICATION

High-pressure vessels are containers designed to hold gases or liquids under a pressure substantially higher than the ambient one. Common pressure vessels are made of steel. However, since rolled and forged parts were welded together during fabrication, some mechanical properties could be adversely affected, which leads to increase the thickness to resist the high pressure. So, such all-metal construction is usually very heavy, weighing approximately 1.4kg/L. To save the weight, there are some other manufacturing techniques based on composite materials. For example, the constructions of metal liner with full composite overwrap (generally aluminum with a carbon fiber composite) and allcomposite structure comprising a polymer liner with carbon fiber composite allow to achieve the weight parameter 0.45 kg/L and 0.3 kg/L respectively (www.compositesworld.com). This section presents the manufacturing process and relevant robotic system used for fabrication of a high-pressure vessel covered by thermoplastic fiber at CETIM company (Nantes, France).

Manufacturing Task Description and its Regularization

The high-pressure vessel considered in this chapter is composed of a cylindrical part and two elliptical domes at both ends of the cylinder. The cylinder has the diameter of 168 mm and the length of 400 mm, the maximum pressure is 500 bar. The vessel is covered by several layers of the thermoplastic fiber as shown in Figure 57, which are winded by the robot and consolidated in situ following certain patterns.

Figure 57 The high-pressure vessel with thermoplastic fiber covering, ©CETIM.

To generate the lay-up patterns, special software from Composicad was used that allows user to create the winding path for the variety of products such as pipes, tubes, tanks, vessels, etc. (www.composicad.com). This software includes a number of packages targeted at specific geometry of the workpiece. For the considered pressure vessel, the dedicated package Composicad-Silver was used, which is specially designed for shapes with axial symmetry. This package allows user to specify each layer with the information of lamina type (a helical winding, a circumferential winding, a connector or a transition winding), weight, thickness, etc. Then, the software calculates the minimum number of circuits required for complete coverage. For a number of layers, Composicad also calculates the possible pattern repeats. Finally, the generated pattern is exported to the xls/.txt file that contains the Cartesian positions of the winding path as well as the normal direction vectors, and winding angles. It is also displayed as shown in Figure 58. It is worth mentioning that Composicad provides a functionality of time optimized motion generation while fixing the tool path with respect to the robot base frame (called as the "time optimal fixed trajectory path generation" in Composicad-Silver). Within this functionality, the user can specify the maximum Cartesian speed/acceleration of the technological tool as well as the maximum angular speed/acceleration of the workpiece rotation provided by the positioner. It is clear that, because of constraining the tool path in the global frame, such approach cannot be treated as strictly time-optimal compared to those developed in this work. Besides, the robot kinematics and its motion capabilities expressed by the maximum velocities/accelerations in the actuated coordinates are not taken into account in Composicad. Nevertheless, this quasi-optimal technique, which ignores redundancy of the robotic system, can be acceptable in some practical applications. On the other side, it gives user some preliminary estimation of the fiber-laying time that can be further reduced by applying the proposed motion planning algorithms for the redundant robotic systems.

For the high-pressure vessel considered here, a sample of the laying task ensuring a single circuit placement of helical lamina on the cylinder and two elliptical domes was created using the Composicad-Silver. Relevant data containing description of a discrete 3D-augmented curve were exported into an Excel file. They include Cartesian positions and normal directions of the laying path with respect to the workpiece frame as shown in Figure 59. To describe this path in the desired way, there were applied expressions ( 14)-( 17) allowing presenting the laying task in the form of a sequence of 4×4 homogeneous transformation matrices. Besides, an equivalent representation was obtained using expressions ( 20)-( 22) that yield relevant sequence of 6×1 vectors composed of Cartesian coordinates and Euler angles, which are commonly used in robot control and programming. By analyzing the task data file obtained from Composicad-Silver, it was noticed that the path sampling was not very regular. In particular, the sampling step varies from 1.5 to 21.0 mm and it is higher for cylinder surface and lower for the domes (see Figure 59). In practice, implementation of very short path segments creates certain difficulties for robot controllers that impose strict constraint on the minimum motion time between the trajectory nodes. The latter does not allow achieving the desired speed of the technological tool if the nodes are too close to each other. So, before applying the developed motion planning algorithms, the task data were modified and presented in the form with a regular sampling step. Besides, from analysis of the original data file generated by Composicad-Silver, it was observed that there are rather sharp variations on some Euler angles profiles, especially at the connection of the cylinder and the domes (see Figure 60a). It is clear that this phenomenon will leads to unsmooth robotic system motions if task path is not improved. In order to avoid these, the original task was slightly smoothed using cubic spline before implementing the motion planning algorithms. The final form of the laying task data is presented in Figure 60b.

Selection of Robotic System Components and Workcell Lay-out Design

For manufacturing of the considered high-pressure vessel, the redundant robotic system called the SPIDE-TP platform (CETIM, Nantes) was used. This platform can be utilized for fabricating vessels of different sizes (diameters from 25 mm to 500 mm and lengths up to 3500 mm), including the vessel with the length 400 mm and diameter 168 mm presented in Figure 57. The SPIDE-TP system is composed of a 6-axis robot with 1-axis linear track from KUKA, a 1-axis workpiece positioner and a special technological tool from AFPT (see Figure 61). The technological tool from AFPT used here is specially designed for producing rotational symmetric workpieces (www.afpt.de). It perfectly suits the shape of high-pressure vessel considered in this work. This tool is presented in Figure 62 and is composed of a compaction roller, a thermal camera, a heat source (diode-lasers), a material storage (fiber creel) and a tensioning system inside. The tool mechanics guides the thermoplastic fiber towards a mandrel or mold where the laser is used to heat the fiber to its melting temperature. Then the tape is consolidated under pressure by the compaction roller (adjusted by an air piston-cylinder) and cooled down. During the laying, the fiber temperature is maintained at the requested level using feedback control of the laser power (maximum 4 kW) based on the measurements from the thermal camera. This ensures optimal consolidation and prevents degradation of the thermoplastics. The workpiece positioner (see Figure 61) is also from AFPT company. It has one horizontal rotational axis and the following mandrel parameters: height of 1089.90 mm and length of 4000.00 mm. This positioner rotates the workpiece in the full range 180° simultaneously with the robot actuators, with the maximum velocity 142°/s and the accelerating time 0.5 s. It is clear that such positioner size is excessive for the high-pressure vessel considered in this work (with length of 400 mm and diameter of 168 mm), but it is used here to increase the robotic platform versatility.

The 6-axis robot employed here is KUKA KR210 R3100 ultra, which has a rather high payload of 210 kg and a maximum reach of 3095 mm (www.kuka.com). More detailed information concerning the robot workspace is presented in Figure 63. The manipulator joint limits, velocities and accelerations are given in Table 11. This robot is equipped with the digital controller KR C4, similar to other types of KUKA robots. The robot is installed on the linear track KUKA KL 2000, which allows extending the work envelop up to 4500 mm. This linear unit has maximum payload 2000 kg and provides a translational motion along a horizontal axis with the speed up to 1.96 m/s. The linear track is connected to the robot controller KR C4 coordinating actuation of six robot axes and two external axes (of the linear track and the rotational positioner).

Mutual location of the SPIDE-TP platform components, i.e. the robotic system layout, is predetermined by the platform manufacturer AFPT. Here, the liner track is paralleled the positioner mandrel with the distance 1396.08 mm between their centerlines. The workpiece is attached to the positioner mandrel at the distance of 1077.18 mm to the flange center, which in our geometric model is the distance between the origins of the workpiece frame and the positioner flange frame. Such design ensures the considered vessel to be completely covered by the robot work envelop.

It is clear that the presented robotic system is redundant with respect to the considered technological task, which requires six degrees of freedom only. In fact, in addition to the usual six degrees of freedom of the robotic manipulator, there are two extra degrees of freedom here (rotational and translational ones) that are provided by the positioner and the linear track. For this design, the workpiece rotation allows adjusting its orientation for better reachability of the technological tool while the robot base translation increases the manipulator working range allowing processing large dimensional products.

For considered manufacturing task that deals with rather small workpiece, the robot can easily execute the composite tape laying process without changing its base location. Besides, the minimum travel distance of the linear track is comparable with the vessel length. So, for this product, there is no reason to activate the linear track during the technological process. Nevertheless, the problem of optimal robot placement still exists here while it is limited to the optimization of the linear track coordinate defining the robot base location relative to the workpiece. In this work, this problem was solved in straightforward way, by sampling the linear track coordinate and applying the algorithm developed in Chapter 3 to generate the time-optimal motions for each robot base location. Then, the optimal robot location was found by selecting the smallest motion time of the robotic system. To speed up this procedure, the motion planning Algorithm 2 was applied, and the discretization step for the redundant variable was relatively high (it was equal to 2°). Relevant computational results are presented in Table 15 andFigure 64. As follows from them, the optimal robot location corresponds to the linear track coordinate 3000 mm, which ensures the robotic system motion time 4.11 sec. It is worth mentioning that here the optimal robot placement problem does not have very clear minimum of the objective function, which is only 6.4% less than the biggest value corresponding to the track coordinate 2000 mm. So, in practice while optimizing the robot location, it is reasonable to pay primary attention to accessibility of the task points and collision avoidance. Figure 64 Selecting the optimal robot location on the linear track.

ROBOTIC SYSTEM MOTION PLANNING USING DEVELOPED ALGORITHMS

To generate the time-optimal motions of the considered robotic system, the redundant variable (positioner joint angle) was discretized and it was created the task graph containing all robot/positioner configurations corresponding to the given laying task. Then, the developed DP-based algorithms were applied to this graph to find the desired trajectory.

Generation of the Task Graph for the Thermoplastic Tape Lay-up

To apply the developed algorithms to the time-optimal motion planning, the considered technological task must be presented in the form of the directed graph. Relevant technique is described in Section 3.1, which includes discretization of the redundant variable and creation of the graph vertices that are also referred to as the configuration cells. To generate the desired graph, the redundant variable (the positioner joint angle P q ) was discretized with the step 1°, which produced the initial graph with 57399 vertices that are arranged in 361 rows and 159 columns. For each vertex, to verify the admissibility of the corresponding task point, it was solved the robot inverse kinematics and it was also checked the collisions and the distance to singularities. Application of these constraints allowed reducing the number of the graph vertices down to 24089. More detailed information concerning the task graph generation is presented below. 16 where each line contains the Cartesian coordinates of the corresponding frame origins and Euler angles describing their positions and orientations with respect to the global coordinate system. For computational convenience, all these location vectors were also presented in the form of 4×4 homogeneous transformation matrices. It should be mentioned that in the table the robot base location is defined directly in the global frame (taking into account the linear track coordinate 3000 mm). Besides, mutual location of the task and defines combinations of the angle signs allowing to move axes through angles greater than +180° or less than -180°. More details concerning the configuration indices are given in Table 18 andTable 19. For example, the entry S 'B110' defines the robot configuration with shoulder forward, elbow up and wrist down. Similarly, the entry T 'B001010' means that the angles 2 R q and 4 R q are negative while the others are all positive. In this work, the continuity of the composite laying process did not allow to change the shoulder and elbow postures during the robot motion while the wrist posture might be adjusted. For this reason, only two values of the S index and four values of T index were used: 'B010', 'B110' for the index S and 'B000010', 'B001010', 'B100010' and 'B101010' for the index T.

  0 6 q   0 5 q   0 4 q   0 3 q   0 2 q   0 1 q 1   0 6 q   0 5 q   0 4 q   0 3 q   0 2 q   0 1 q
When generating the task graph, the robot kinematic model was used first, which allowed to check accessibilities of the task points and to evaluate distances to manipulator singular configurations. At this stage, all 57399 vertices were verified and some of them were eliminated from the task graph, if one of the following conditions was satisfied: (i) the inverse kinematic solution did not exist, (ii) the joint coordinates were outside of the limits or (iii) the Jacobian condition number was more than 6. The task graph after applications of kinematic constraints is shown in Figure 66a, it includes 30229 admissible vertices. Then, the difficult configurations were eliminated which allowed reducing the number of the admissible vertices down to 28907, as shown in Figure 66b. Finally, the rough/fine collision constraints were applied using Matlab codes and commercial CAD system CATIA, and the final task graph was generated with 24089 admissible vertices (see Figure 66c).

In order to use collision detection capabilities of CATIA, a detailed 3D model of the SPIDE-TP platform was created and "DMU Kinematics" workbench was activated. Then, using the visual interface of the "Check Clash" module, two types of interferences between the system components were selected: 1) interference between the technological tool and robot/positioner; 2) interference between the robot links and other components of the workcell. Further, by activating the sensors on the "Kinematic Simulation" panel and executing this application, the collisions between the components were detected, if any. Corresponding areas are highlighted in yellow in the CAD graphical window and relevant sensor value (interference flag) becomes non-zero, as shown in Figures 67 and68. Finally, the collision check data for all vertices were recorded to an Excel file describing the task graph that was used further to generate optimal motions of the robot and positioner. It should be mentioned that generation of such graph required quite a lot of computing time. In particular, it took about one minute to generate the graph taking into account the kinematic constraints and to eliminate difficult configurations. However, the rough/fine collision test in CATIA required more than 10 hours. The generated optimal trajectories are presented in Figure 69, which contains the time profiles for all joint coordinates and corresponding velocities estimated using the moving window technique. As follows from detailed analysis, the obtained profiles are quite smooth and satisfy both the joint limits and the velocity/acceleration constraints. It is worth mentioning that at each time interval one of the constraints is active, i.e. either the velocity or acceleration reaches its maximum/minimum value for at least one of the joint coordinate. Besides, it was observed that after the trajectory smoothing, some segments of the velocity profiles were slightly over their limits (for the positioner and axis #4 of the robot). However, such minor violations of the constraints in Figure 69 are usually acceptable in practice and are compensated by the robotic system controller. On the other hand, these violations can be easily eliminated by simple modifications of the trajectories, by means of the motion time extension between the corresponding task points.

Figure 69 Profiles of time-optimal motion for thermoplastic fiber covering of the high-pressure vessel It is worth mentioning that for the considered technological problem, the task graph has rather complicated topology that includes some bottleneck areas. In particular, it can be seen a very narrow gap between the cells corresponding to the task points #89 and #90 (Figure 70). So, there is very limited number of the graph vertices to be included in the optimal path for this area. From engineering point of view, this segment of the laying task is very difficult to be implemented without collisions. In fact, as follows from the simulation of the time-optimal motion in CATIA V5, the technological tool moves very close to the robot forearm during the transition between the task points #89 and #90 (see Figure 70, where the robot wrist was set to be invisible for better view). In addition, it is reasonable to compare the obtained results with ones generated using another technique that was developed for the laser cutting applications [START_REF] Pashkevich | Multiobjective optimization of robot motion for laser cutting applications[END_REF] and assumes constant Cartesian speed of the robotic tool with respect to the workpiece. In the frame of this work, it corresponds to the constant laying speed of the thermoplastic fiber. Using the latter assumption and relevant objective function (maximum weighted joint coordinate increment between the task points), the motion planning problem was solved in different way. The obtained results are presented in Figure 71, which clearly shows advantages of the developed technique that allows reducing the motion time of the robotic system from 8.4 sec down to 4.0 sec.

Figure 71 Laying speeds for the motions generated using the developed and known algorithms It is also worth mentioning that the laying task considered here deals with just one circuit of the laminate. In practice, the laying process includes numerous circuits and robot/positioner motions should be repeated many times, with some small adjustment of the tool path. For this reason, to meet the demand of the industrial partners, it was slightly modified the final step of the Algorithm 2 that deals with selection of the final robot configuration from the candidates belonging to the last column of the task graph matrix. Initially, it was selected to minimize the motion time only. After the modification, the selection of the final configuration was based on the combined criteria, which takes into account both the motion time and the difference between the initial and final robot configurations (by applying the weighted sum method). This modification allowed us to achieve almost identical robot configurations at the circuit beginning/end while it yielded very small increasing of the system motion time, by 0.09 sec only. It should be noted that the time-optimal profiles presented in Figure 69 were generated using the modified version of the Algorithm 2. So for this motion, the initial and final values of the robot joint coordinates are very close to each other.

Figure 72 Configurations of the robotic system for the time-optimal trajectory simulated in CATIA To verify correctness of the obtained trajectory, it was carried out a number of simulations in a 3D visual environment of CATIA. It was used the same 3D model of the SPIDE-TP platform as for the collision test at the task graph generation stage (see Figures 67 and68). To simulate the obtained motion, the "DMU Kinematics" workbench was activated. Then, using the visual interface of the "Simulation with Laws" module, the robotic system model was actuated using the Excel file describing the sequential configurations of the manipulator and positioner. Several pictures captured from this simulation are presented in Figure 72, which shows the robotic system configurations corresponding to the trajectory nodes: #1, #5, #12, #18, #25, #32, #40, #50 and #54. Relevant video is available at the URL https://www.youtube.com/watch?v=ioGFXuV2gDU. As follows from the relevant analysis, the simulated motions of the robot and the positioner are well coordinated and meet the technological requirements of the thermoplastic lay-up process. Besides, the initial and final locations of the technological tool are close to each other (see Figures 72a and72i), which allows in practice easily repeating the same laying motion several times. However, it should be mentioned that this simulation allowed to verify correctness of the obtained laying path only, while evaluation of the laying speed was not possible in this way because 3D animation provided by CATIA does not take into account the time information concerning the examined trajectory.

Hence, the developed algorithm allowed us to generate the desired time-optimal trajectory that should be further implemented by the robotic system controller. Relevant data were presented in the form of the Excel file containing the robot and positioner joint coordinates and corresponding time instances. Using this data, there were also calculated the robot configuration indices and 6×1 location vectors describing the technological tool position and orientation in the global frame, which are required for the robot programming.

OPTIMAL MOTION IMPLEMENTATION IN ROBOTIC SYSTEM

To implement the obtained time-optimal motion to the SPIDE-TP platform, the trajectory was programmed using the KRL language for KUKA robotic system controllers. Then, an experimental evaluation on the factory floor was carried out to verify its applicability.

Motion Implementation using KRL Robot Programming Language

The time-optimal trajectory obtained in the previous section is presented as a sequence of the robot/ positioner joint coordinates and corresponding time instances defining the motion profiles. However, most of modern robot controllers neglect the time information and apply their own built-in algorithms to generate trajectory using only the sequence of the desired end-effector locations, which can be defined by either the joint coordinates or the end-effector position/orientation. Let us give some details concerning typical motion commands implemented in the robot controller and some corresponding analysis (investigation). In the control system KR C4 that is used in SPIDE-TP platform, the configurations of the manipulator and external mechanisms (positioner and liner track) are described by either 8 joint coordinates or 10 numbers defining the end-effector position/orientation, the positioner rotation angle, the linear track displacement as well as the manipulator configuration indices. An example of such descriptions is presented in Figure 73. Here, the first line straightforwardly defines the manipulator joint angles denoted as {A1, A2… A6}, the linear track displacement {E1} and the positioner rotation angle {E2}. In contrast, the second line uses an alternative description of the manipulator posture via the Cartesian coordinates {X, Y, Z} and Euler angles {A, B, C} of the end-effector as well as via the manipulator configuration indices {S} and {T} presented in Tables 18 and19. It is worth mentioning that in our experimental study the linear track coordinate E1 was constant and it was set to its optimal value equal to -3000 mm (as justified in Section 4.1).

To describe motions of the robotic manipulator and external mechanisms, the control system KR C4 uses the programming language KRL. In the frame of this language, there are two possible ways to implement the time-optimal trajectory obtained in Section 4.2. The first of them is based on the conventional robot motion commands PTP/LIN, while the second uses SPLINE blocks. The command PTP implements the fastest linear-path motion in the joint space taking into account the velocity/acceleration constraints for each actuated coordiante. It is clear that in the Cartesian space this command produces a curved path, so it is not reasonable to use it for implementation of the desired motion of the robotic system for the compsite tape laying. The command LIN implements the fastest linear-path motion in the Cartesian space taking into account the velocity/acceleration constraints for each actuated coordiante as well as the Cartesain velocity/acceleration constraints for the robot endeffector. This command produces a linear segment in the Cartesian space, so it can be used for implementation of the desired motion assuming that the piecewise linear interpolation of the composite lay-up path is acceptable from engineering point of view. It should be stressed that, the basic versions of the PTP and LIN commands implement "start-stop" motions that includes three sections (accelerating, constant speed and decelerating) with zero speeds at the beginning and the end. It is clear that such "start-stop" implementation is not accepetable for the considered composite lay-up technology where the end-effector Cartesian speed should be as high as possible and the motion stops at the nodes are not permited. To avoid this diffculty, the programming language KRL allows user to apply so-called approximate positioning at the trajectory end points, which is based on the superposition of deceleration and acceleration sections of the subsequent motion segements. This technique yields non-stop continuous motions for which the end-effector passes the nodes neighbourhood only instead of visiting the nodes exactly (see Figure 74), so it can be hardly accepted for the considered tape laying task. To activate the approximate positioning in KRL language, the robot motion command must include a special suffix, such as C_PTP for the PTP-PTP sequence or C_DIS, C_ORI or C_VEL for the LIN-LIN sequence. The latter define the type of switching criteria that can be based on the translational/ rotational distance to the target node or the velocity level at the beginning of the approximation. The overlapping degree for the approximate positioning is defined via the special system variable $APO.CPTP or $APO.CVEL. For example, the setting $APO.CPTP=100 or $APO.CVEL=100 leads to the complete 100% overlapping of the previous deccelerating and the subsequent accelerating sections (KUKA, 2010).

Durations of the PTP and LIN motions is computed by the KR C4 controller online, using the velocity and acceleration constraints in the joint space (for all axes) and similar constriants in the Cartesian space for the end-effector motion. The default values of these constraints are defined in special tables, but they can be customized by the user by means of the dedicated commands of the KRL language. For instance, the commands $VEL_AXIS [1]=50 and $ACC_AXIS [1]=75 define the maximum velocity and acceleration for the axis#1 on the level of 50% and 75% of their default values correspondingly. Similarly, the commands {$VEL.CP=2, $VEL.ORI 1=400, $VEL.ORI 2=400} define the maximum end-effector translational speed as 2 m/s and the maximum rotational velocities for the orientation angles as 400°/s. Related commands {$ACC.CP, $ACC.ORI 1, $ACC.ORI 2} can be also used to set maximum translational/rotational accelerations in the Cartesian space.

Using the velocity and acceleration constraints, the KR C4 controller generates trapezoidal velocity profiles either in the joint or Cartesian space and applies the superpositioning to avoid the stops at the intermediate task points (see Figure 42). It is also worth mentionning a very important particularity of the motion generation in the robot controller, where a sequence of very short segments cannot be executed with maximum allowable speed. In this case, the velocity profiles have the triangular shapes whose duration is limited by the controller capacity (it is defined in the control software as a constant). So, the resulting average speed for such segments is essentially lower than the expected level, which is defined by the dedicated commands. Hence, it is necessary to avoid very short segments while implementing the time-optimal trajectory via the LIN-LIN sequences (another disadvantage of such implementation is related to the approximate positioning). An alternative way to implenment the desired motion is based on the SPLINE block commands. Their execution employs the quintic polynomial interpolation between the trajectory nodes, instead of the first-order one that is used for LIN or PTP. Besides, compared to the approximated LIN-LIN sequences, here the end-effector passes the given nodes exactly and without stops (see Figure 75). For this reason, the spline based technique is very attractive for implementing of complex paths. In KRL programming language, a spline motion is defined as a set of several individual segments. The segments are grouped together to form the overall motion in a so-called spline block, which is executed by the robot controller as a single motion. The SPLINE block is limited by the SPLINE and ENDSPLINE statements and contains a sequence of commands SPL describing the trajectory nodes. In addition, the desired the motion time for each segment can be defined indirectly by setting the system variables {$VEL.CP, $VEL.ORI 1, $VEL.ORI 2} or in the direct way, using special TIME_BLOCK containing the travelling times between the nodes. An example of such descriptions is presented in Figure 75.

Inside of the SPLINE block it is possible to use the SPL commands both outside and inside of the TIME_BLOCK. If the SPL command is located outside of it, the motion time between the current and target nodes is computed using the maximum velocity/acceleration settings provided by the relevant system varibables, similar to the LIN-LIN sequence. Otherwise, if the SPL command is inside of the TIME_BLOCK, the motion time between the nodes is defined via the statement TIME_BLOCK PART, as a percentage of the total spline motion time that is specified at the TIME_BLOCK END statement. For example, for the KRL program presented in Figure 76 the motion from the node P1 to P54 is executed in 4.7 sec, while for the segment [P1, P2] the motion time is set as 3.0% of the total one, which corresponds to 0.14 sec. It should be also noted that this program also includes the SPL commands before and after the TIME_BLOCK that allow achieving the desired velocity/acceleration values at the spline motion beginning and the end. Hence, as follows from the comparison study presented above and summarized in Table 21, the spline interpolation is the most attractive for implementation of the time-optimal motions generated by the developed algorithms. Its main advantages are the capability to implement continuous motion without stops, exact positioning at each task location as well as possibility to assign the desired motion time for each trajectory segment. Besides, short path segments are also acceptable for the SPLINE blocks available in the KR C4 controller since they do not cause the velocity reduction, in contrast to the LIN-LIN sequences.

Experimental Evaluation of the Implemented Time-Optimal Motions

To evaluate the implemented time-optimal trajectory generated by the developed algorithm, the obtained motion was programmed using the SPLINE block. Relevant data describing the robot and positioner joint coordinates (trajectory nodes) and corresponding time instances were imported from the .xls file into a text editor and presented in the form combining the end-effector location and the joint coordinates for the external axes using the format shown in Figure 73. The latter yielded 159 trajectory nodes denoted in the KRL program as P[1],…P [159], which were used as the arguments in the SPLINE block commands SPL P[1],…SPL P[159]. Besides, there were defined the HOME configuration and the starting point P0 preceding the spline motion. There were also set the system variables defining the maximum allowable speeds/accelerations and added some technological commands for the control of the thermoplastic tape feeding, tensioning and heating (see Figure 77). The obtained text file containing the program in the KRL language was saved with the .src extension that is compatible with KR C4 controller. Finally, the robotic system motion program was downloaded to the controller and executed. It is also worth mentioning that it was discovered after the first experiments that some of the nodes could be eliminated without any infulence on the quality of the thermoplastic tape laying. This allowed us to reduce the number of the trajectory nodes down to 54, which were used in further experimental study.

The experimental study included execution of the above presented KRL program and evaluation of the corresponding motions of the robotic manipulator and positioner. For safety reasons, the robotic system was run with slightly reduced speed, which was achieved by setting the maximum velocities of the actuated axes at the level of 75% of the maximum values. Besides, the workpiece was not mounted on the positioner flange to avoid unexpected collisions. The system motions were registered by a video camera, and the total motion time was estimated using a timer. The experiment shows that even for these settings it took only 7.2 sec to implement the desired laying path, while the Composicad software produced the trajectory with the execution time more than 14.0 sec. It is clear that the motion time for the studied trajectory can be easily reduced down to ~5.4 sec by simple increasing of the maximum velocity settings up to 100%. Nevertheless, it is still higher compared to the 4.0 sec corresponding to the theoretical time-optimal trajectory presented in Section 4.2 (this issue is discussed in details below).

Several photos from this experiment are presented in Figure 78, which shows the robotic system configurations corresponding to the trajectory nodes: #1, #5, #12, #18, #25, #32, #40, #50 and #54. Relevant video is available at the URL https://www.youtube.com/watch?v=9nvdlzpJHFE. As follows from the related analysis, the implemented motions of the robot and the positioner are well coordinated and rather fast. Besides, the initial and final locations of the technological tool (Figures 78a and78i) are close to each other, which allows in practice easily repeating the same laying motion several times. However, there are several segments where the motion smoothness should be impoved if the allowable velocities are increased up to 100% of their maximum values.

; ------------------------------------------------Initialization Section -------------------------------------------------------HOME = {A1 0.00, A2 -120.00, A3 135.00, A4 0.00, A5 75.00, A6 90.00, E1 -3000.00, E2 0.00} ; ---------------------------------------------Trajectory Points Definition ------------------------------------------------- ; -----------------------------------------------------Main Section -----------------------------------------------------------rTapeFeed("on","fly") ; Tape feeding system activated rTapetension("on","fly") ; Tape tensioning system activated rlaserstandby("on","all") ; Laser heating system ready rlaserbeam("on","all") ; Laser beaming rtapeCut("off","all") ; Cutter inactivated Detailed analysis of the non-smooth trajectory segments shows that such imperfect behavior was caused by some drawbacks of the Composicad software, which provided the developed motion planning algorithms with initial data (the Cartesian coordinates and orientation angles of the laying path nodes). In fact, the non-smooth segments are located at the connection of the cylinder and the domes (see Figure 60a) where ComposicaD allowed rather sharp variations of the orientations.

Although the given initial laying path and the generated time-optimal trajectory were slightly smoothed (see Figures 60b and69), the obtained motion is still difficult for implementation using the capabilities of the KR C4 controller.

Figure 78 Configurations of the robotic system for the time-optimal trajectory implemented using SPLINE block with the speed settings

Another important issue to be discussed here is related to the robotic system motion time obtained in the experiments. The problem is that the duration of the implemented motion is essentially higher than the expected theoretical value. In fact, the developed motion planning algorithm yielded the trajectory with the motion time 4.0 sec, while the best value achieved in the experiments was 5.3 sec (using the 100% settings for the allowable speeds). The most probable reason for this motion slowdown is the difference between the nominal values of the maximum velocities/accelerations declared in the technical documents (used as the principle constraints for the motion planning) and their real values specified in the control software. In practice, ordinary users can manipulate with the percentage of the maximum velocity/acceleration for each actuated axis only, while they do not have access to setting their absolute values. Besides, some additional time is required for accelerating/deccelerating when entering/exiting the SPLINE block, within which the velocities/accelerations are maintained at the level corresponding to the time-optimal motion for a single laying circuit.

; -------------------------------------------------------Main Section ---------------------------------------------------------rTapeFeed("on","fly") rTapetension("on","fly") rlaserstandby("on","all") rlaserbeam("on","all") rtapeCut("off","all")

; ………………………………………………………………………………………………………………… P0 = {X 19.62, Y 50.00, Z 102719.62, Y 50.00, Z .28, A 41.11, B 80.86, C 90.10, E1 -3000.00, E2 -138.90, S 2, T 2} PE= {X 13.42, Y 50.00, Z 1027.53, A 60.42, B 80.95, C 106.93, E1 -3000.00, E2 2..53, A 60.42, B 80.95, C 106.93, E1 -3000.00, E2 2. To overcome the above mentioned difficulties, the KRL program of the time-optimal motion was modified using the TIME_BLOCK allowing to specify the motion time for each path segment. Besides, additional commands LIN P0 and SPL PE were included in the main section (before and after the TIME_BLOCK) in order to achieve the desired velocity/acceleration at the spline motion beginning and end (see Figure 79). This modified version of the KRL program will be tested soon by our industrial partner. There are a number of reasons to expect that this modification allows implementing the time-optimal trajectory with better precision and reaching the theoretical level of the motion time. On the other hand, detailed information concerning implementation of the SPLINE and TIME blocks in KR C4 controller is not available in technical literature, so some additional experiments will be conducted in future to reduce distance between the theory and practice.

SUMMARY

This chapter deals with industrial implementation of the developed motion planning method on the factory floor (for manufacturing of a high-pressure vessel). Using the proposed method, there were generated time-optimal smooth motions for the manipulator and positioner of the SPIDE-TP robotic platform allowing speeding up the thermoplastic tape lay-up process. Correctness of these motions was verified in 3D simulation environment of CATIA software package. Besides, it was created a program in KRL language implementing these motions and ensuring coordinated control of KUKA KR210 robot and AFPT workpiece positioner. There are presented results of industrial experiments carried out with this program, which confirmed smoothness of the manipulator/positioner movements and essential reduction of the time required for the thermoplastic tape lay-up process.

In more details, the contribution of Chapter 4 can be summarized as follows (i) Manufacturing process preparation for the given technological task via defining and regularizing the task locations on the lay-up path provided by an industrial partner and optimization of the workcell layout by selecting the best robot base location.

(ii) Generation of the time-optimal coordinated motion for the robotic manipulator and the workpiece positioner that implement the given lay-up task. The obtained time-optimal trajectory allowed essentially reducing the total motion time of the robotic system required for one circuit of the composite laying process.

(iii) Implementation of the generated time-optimal motion. The obtained trajectories were used for creating a motion control program in KRL language, which was tested on the factory floor. The experiment showed that it took only 7.2 sec to implement the desired circuit of the laying path using the time-optimal trajectory, while the commercial Composicad software produced the trajectory with the execution time more than 14.0 sec.

Hence, the industrial experiments confirmed advantages of the developed optimal motion planning method and achieving the main objective of this work, which is targeted at the productivity improvement of the robotic systems for composite lay-up processes. Nevertheless, during these experiments, it was also detected some imperfect behavior of the robot and positioner that gives some perspectives for the future work.

The main results of Chapter 4 have been presented at the conference ICOME'2017.

CONCLUSIONS AND PERSPECTIVES

CONTRIBUTION OF THE THESIS

This thesis is devoted to the optimal motion planning in redundant robotic systems allowing improving the productivity of the automated composite lay-up workcell. Theoretical contributions of the thesis are in the area of redundancy resolution and time-optimal smooth motion generation for robotic systems composed of the robotic manipulator, workpiece positioner and workspace extension unit. The most essential results can be summarized as follows:

(i)

Formalization of the optimal motion planning problem for typical robotic lay-up platforms. The proposed approach presents the considered problem as finding of the timeoptimal trajectory in the joint space of the robotic system under specific constraints related to the technological task. The latter include both conventional kinematic constraints (joint limits, maximum joint velocities/accelerations), allowable distances to the singularities and obstacles as well as the Cartesian path constraints issued from the composite lay-up process.

(ii) A new method of coordinated time-optimal motion planning for redundant robotic systems ensuring smooth movements of all mechanical components. It is based on transformation of the original optimization problem into a combinatorial one and application of the dynamic programming principle. At the first stage, the redundant variables are discretized and the task graph is created by sequentially applying to all task locations the direct kinematics of the positioner (and workspace extension unit) as well as inverse kinematics of the robot. Then, the developed optimization algorithm based on dynamic programming generates the time-optimal motion taking into account all constraints related to the system kinematics and considered technology. Efficiency of the developed method and its advantages compared to the conventional techniques are confirmed by several case studies.

(iii) Enhancement strategies for the developed motion planning method. To reduce the computing time required for the motion planning, it was proposed to avoid combinatorial search in high-dimensional space. To achieve this target, the time-optimal motion planning is divided in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with relatively large discretization step. Then, the same technique is applied in the local search space, which is created with small discretization step in the neighborhood of the trajectory obtained at the previous stage. As an alternative, the second stage may implement a straightforward smoothing of the redundant variable profiles based on the spline approximation. As follows from relevant study, the proposed enhancement strategies allow essentially reducing the computing time, down to the level acceptable in engineering practice.

(iv) Application of the developed optimal motion planning method to real industrial problems.

Using the thesis results, there were generated time-optimal smooth motions for the manipulator and positioner of the SPIDE-TP robotic platform allowing speeding up the thermoplastic tape lay-up process for manufacturing of the composite high-pressure vessel.

For these motions, it was created a program in KRL language ensuring coordinated control of KUKA KR210 robot and AFPT workpiece positioner. Experimental evaluation of the motions described by this program confirmed smoothness of the manipulator/positioner movements and essential reduction of the time required for the thermoplastic tape lay-up process.

The obtained results contribute to the area of redundancy resolution in robotic systems and give the user an efficient method for generation of time-optimal smooth motions whose implementation allow improving productivity of the automated composite lay-up workcell.

LIMITATIONS OF THE OBTAINED RESULTS

In spite of the essential advantages, there are still some limitations that are related both to the developed method of the time-optimal smooth motions planning and their implementation on the factory floor. The most significant of them are presented below.

Limitations of the developed method of optimal motion planning:

(i)
Computing time for optimal motion planning may be too high if the redundancy degree is greater than one. The most important bottleneck is the fine collision test in the CATIA environment that is integrated in the preparation of the task graph. For example, for the industrial case study considered in Chapter 4, it took almost ten hours to test collisions for about 25000 robotic system configurations even setting triangulation accuracy down to 4.0 mm (0.2 mm by default). So, adding even a single additional redundant variable with 100 discrete values leads to increasing of the collision test time up to 1000 hours, which can be hardly accepted in practice. For this reason, the speeding up of the collision test deserves some special efforts.

(ii) In practice, the principle constraints on the velocities and accelerations in the actuated joints of the robot, positioner and workspace extensioner are not known exactly. It was discovered during the industrial experiments that the maximum velocities and accelerations declared in the manual are higher compared to the real values fixed in the KR C4 controller software. The latter leads to over-estimation of the actuator capacities and under-estimation of the motion time obtained using the developed method. For this reason, it is necessary to conduct dedicated experiments with the industrial robots allowing evaluating correct values of the maximum velocities/accelerations used in the motion planning algorithms implemented in the KRL programming language.

Limitations of the implementation on the factory floor:

(iii)

The lay-up task description generated by the commercial software may be imperfect.

For example, for the high-pressure vessel, a single circuit of the laying path generated by Composicad is not perfectly smooth with respect to the tool orientation angles. It includes discontinuities located at the connections of the cylinder and the domes (see Figure 60a). Although the given initial laying path and the generated time-optimal trajectory were slightly corrected (smoothed using the polynomial spline), the obtained motion is still difficult to implement using the KR C4 controller capabilities. As follows from the industrial experiments, the robotic system movements are a little bit edgy in the neighborhood of the connection points. For this reason, the quality of the lay-up task description deserves particular attention.

(iv)

The geometric model of the robotic system from the industrial partner may be imperfect and contain inaccurate parameters. In particular, the TOOL parameters (describing transformation between the robot flange and the tool center point) were not well calibrated.

In our experiments, the latter caused some interference between the technological tool and the workpiece shaft, while no collisions were detected in 3D simulation in the CATIA environment. Besides, the compaction force applied to the thermoplastic by the technological tool was slightly over the desired level at some segments of the laying path. Hence, careful calibration of the geometric model is a very important issue in implementation of the time-optimal trajectories generated using the developed method.

Nevertheless, for the considered application area, the above mentioned limitations are not crucial and the developed method of the optimal motion planning provides essential improvement of the robotic system productivity. On the other hand, these limitations show some research directions for future work.

PERSPECTIVES AND FUTURE WORK

To generalize the obtained results, it is reasonable to continue research in the several directions and to concentrate on the following issues:

(i)
Enhancing the developed method by combining the task graph generation and the optimal path searching on this graph. It will allow reducing the computing time by avoiding exhaustive collision check at the preliminary phase, which is applied to all possible configurations of the considered robotic system satisfying the task constraints. This modification of the proposed motion planning method requires integration of the collision check and verification of the acceleration constraints into dynamic programming based algorithm.

(ii)

Comparing the developed method based on the combinatorial optimization with the continuous nonlinear time-optimal control in the space of the redundant variables and their derivatives or in the space of variables describing the end-effector motion. This approach looks promising but includes numerous difficulties related to the constraints transformation from the robotic system joint space to the considered state space.

(iii) Speeding up the developed method by preliminary segmentation of the given task into several sub-paths corresponding to a single circuit of the lay-up process. This may allow reducing the computing time by sequentially solving a number of low-dimensional combinatorial optimization problems and integrating further the obtained solutions into an aggregated trajectory ensuring implementation of the numerous lay-up circuits. It is worth mentioning that some preliminary work has been already done during industrial

Abstract :

The thesis deals with the optimal motion planning in redundant robotic systems for automation of the composite lay-up processes. The primary goal is to improve the lay-up workcell productivity by developing a novel methodology of optimizing coordinated motions of the robotic manipulator, workpiece positioner and workspace extension unit, which ensure the shortest processing time and smooth movements of all mechanical components. In contrast to the previous works, the proposed methodology provides high computational efficiency and also takes into account both the technological constraints and the robotic system constraints, which describe capacities of the actuators and are expressed by the maximum allowable velocities and accelerations in the actuated joints. The developed technique is based on conversion of the original continuous problem into a combinatorial one, where all possible configurations of the mechanical components are represented as a directed multilayer graph and the desired time-optimal motion is generated using dynamic programming principle for searching the shortest path on the graph satisfying the smoothness constraints. It is also proposed an enhancement of this technique by dividing the optimization procedure in two stages combining global and local searches. At the first stage, the developed algorithm is applied in the global search space generated with large discretization step. Then, the same technique is applied in the local search space, which is created with smaller step in the neighborhood of the obtained trajectory.

The advantages of the developed methodology are confirmed by industrial implementation on the factory floor that deals with manufacturing of the highpressure vessel.

Figure 1

 1 Figure 1 Components of typical fiber reinforced composite material ..................................................... 2 Figure2Strength plotted against density (yield strength for metals and polymers; compressive strength for ceramics, tear strength for elastomers and tensile strength for composites)[START_REF] Ashby | Materials selection in mechanical design[END_REF] ........................................................................................................................................... 3 Figure 3 Composite materials are used to make pressure tanks ©MATERIA (http://www.materiainc.com/) .................................................................................................................................................. 4 Figure 4 Composite materials used throughout the body of the Boeing 787 ©BOEING (https://www.comsol.com) ...................................................................................................................... 4 Figure 5 Common PMCs manufacturing processes in industry .............................................................. 7 Figure 6 Typical filament winding process © CONNOVA (www.connova.com/) ............................... 8 Figure 7 Typical automated tape laying system © MTORRES (www.mtorres.com/) .......................... 10 Figure 8 ATL for Boeing787's unique one-piece composite barrel construction ©BOEING (https://www.boeing.com) ..................................................................................................................... 10 Figure 9 Schematic of automated tape laying and automated fiber placement ..................................... 11 Figure 10 Automated fiber placement process for fuselage section producing with window cut-outs © MTORRES (www.mtorres.com/) .......................................................................................................... 11 Figure 11 (a) Cincinnati VIPER™ platform for automated fiber placement, © Cincinnati VIPER™ (b) TORRESFIBERLAYUP machine for automated fiber placement, © MTorres. ............................. 12 Figure 12 Robot-based AFP system developed by Coriolis Composite © Coriolis Composite (www.coriolis-composite.com) ............................................................................................................. 13 Figure 13 (a) Robot-based AFP system with modular head developed by Electroimpact, ©Electroimpact (b) Robot-based AFP system in the Composites Technology Center of NASA, ©NASA ................................................................................................................................................. 13 Figure 14 Redundant robotic system for application example .............................................................. 16 Figure 15 Simulation of the solution from pseudo-inverse for planar redundant robotic system: (a)  (h) ............................................................................................................................................................... 17 Figure 16 Simulation of the solution from motion coordination for planar redundant robotic system: (a)  (h) ....................................................................................................................................................... 19 Figure 17 (a) visualization of the GA progress for considered application example; (b) Joint coordinates of the positioner from GA based solution. ......................................................................... 23 Figure 18 Manufacturing process planning for robotic lay-up processes.............................................. 25 Figure 19 Typical fiber placement curves defined in 3D CAD model (SolidWorks), with position and normal direction relative to the workpiece frame ©AUTOMATED DYNAMICS ® (www.automateddynamics.com). .......................................................................................................... 33 Figure 20 Definition of task frames with respect to workpiece frame .................................................. 34 Figure 21 Architecture of standard 6-axis industrial robot .................................................................... 37 Figure 22 Different configuration of robotic manipulator for the same location of the end-effector ... 40 Figure 23 Three types of singular configurations in KUKA kinematic system, ©KUKA (KUKA, 2010) ............................................................................................................................................................... 43 Figure 24 Projection of the manipulator links on the vertical plane (Configuration SHOULDER_FORWARD & ELBOW_UP) ........................................................................................ 44 Figure 25 Multiple manipulator configurations for the same wrist center point position ..................... 44 Figure 26 Triangles highlighted for computing the angle of axis #2 .................................................... 45 Figure 27 Typical one-axis positioners from KUKA (www.kuka.com) ............................................... 48 Figure 28 Typical two-axis positioners from KUKA (www.kuka.com) ............................................... 48

Figure 29

 29 Figure 1 Components of typical fiber reinforced composite material ..................................................... 2 Figure2Strength plotted against density (yield strength for metals and polymers; compressive strength for ceramics, tear strength for elastomers and tensile strength for composites)[START_REF] Ashby | Materials selection in mechanical design[END_REF] ........................................................................................................................................... 3 Figure 3 Composite materials are used to make pressure tanks ©MATERIA (http://www.materiainc.com/) .................................................................................................................................................. 4 Figure 4 Composite materials used throughout the body of the Boeing 787 ©BOEING (https://www.comsol.com) ...................................................................................................................... 4 Figure 5 Common PMCs manufacturing processes in industry .............................................................. 7 Figure 6 Typical filament winding process © CONNOVA (www.connova.com/) ............................... 8 Figure 7 Typical automated tape laying system © MTORRES (www.mtorres.com/) .......................... 10 Figure 8 ATL for Boeing787's unique one-piece composite barrel construction ©BOEING (https://www.boeing.com) ..................................................................................................................... 10 Figure 9 Schematic of automated tape laying and automated fiber placement ..................................... 11 Figure 10 Automated fiber placement process for fuselage section producing with window cut-outs © MTORRES (www.mtorres.com/) .......................................................................................................... 11 Figure 11 (a) Cincinnati VIPER™ platform for automated fiber placement, © Cincinnati VIPER™ (b) TORRESFIBERLAYUP machine for automated fiber placement, © MTorres. ............................. 12 Figure 12 Robot-based AFP system developed by Coriolis Composite © Coriolis Composite (www.coriolis-composite.com) ............................................................................................................. 13 Figure 13 (a) Robot-based AFP system with modular head developed by Electroimpact, ©Electroimpact (b) Robot-based AFP system in the Composites Technology Center of NASA, ©NASA ................................................................................................................................................. 13 Figure 14 Redundant robotic system for application example .............................................................. 16 Figure 15 Simulation of the solution from pseudo-inverse for planar redundant robotic system: (a)  (h) ............................................................................................................................................................... 17 Figure 16 Simulation of the solution from motion coordination for planar redundant robotic system: (a)  (h) ....................................................................................................................................................... 19 Figure 17 (a) visualization of the GA progress for considered application example; (b) Joint coordinates of the positioner from GA based solution. ......................................................................... 23 Figure 18 Manufacturing process planning for robotic lay-up processes.............................................. 25 Figure 19 Typical fiber placement curves defined in 3D CAD model (SolidWorks), with position and normal direction relative to the workpiece frame ©AUTOMATED DYNAMICS ® (www.automateddynamics.com). .......................................................................................................... 33 Figure 20 Definition of task frames with respect to workpiece frame .................................................. 34 Figure 21 Architecture of standard 6-axis industrial robot .................................................................... 37 Figure 22 Different configuration of robotic manipulator for the same location of the end-effector ... 40 Figure 23 Three types of singular configurations in KUKA kinematic system, ©KUKA (KUKA, 2010) ............................................................................................................................................................... 43 Figure 24 Projection of the manipulator links on the vertical plane (Configuration SHOULDER_FORWARD & ELBOW_UP) ........................................................................................ 44 Figure 25 Multiple manipulator configurations for the same wrist center point position ..................... 44 Figure 26 Triangles highlighted for computing the angle of axis #2 .................................................... 45 Figure 27 Typical one-axis positioners from KUKA (www.kuka.com) ............................................... 48 Figure 28 Typical two-axis positioners from KUKA (www.kuka.com) ............................................... 48

  ........................................................................................ Table 14 Simulation results for linear lay-up task for different discretization step Δq P ....................... Table 15 Robotic system motion time with respect to different robot base location on the linear track. .............................................................................................................................................................

  PARTS FABRICATION TECHNIQUES ......................................................... 1.1.1 Composite Materials in Manufacturing .................................................................................. 1.1.2 Automation of Composite Products Manufacturing ............................................................... 1.1.3 Industrial Systems for AFP Based Composite Manufacturing ............................................. 1.2 REDUNDANCY RESOLUTION IN ROBOTIC SYSTEMS .................................................... 1.2.1 Redundancy Resolution via Generalized Inverse of the Kinematic Jacobian ...................... 1.2.2 Redundancy Resolution via Coordinating of Robot/Positioner Motions ............................. 1.2.3 Optimization Based Techniques for Redundancy Resolution .............................................. 1.3 ROBOTIC CELL DESIGN AND PROGRAMMING FOR COMPOSITE PRODUCT MANUFACTURING .......................................................................................................................... 1.3.1 Manufacturing Process Planning for Robotic Lay-up Applications ..................................... 1.3.2 Software Packages for Computer Aided Design of Industrial Robotic System ................... 1.4 THESIS GOAL AND RESEARCH PROBLEMS ......................................................................

Figure 1

 1 Figure 1 Components of typical fiber reinforced composite material

Figure 3

 3 Figure 3 Composite materials are used to make pressure tanks ©MATERIA (http://www.materia-inc.com/)

Figure 4

 4 Figure 4 Composite materials used throughout the body of the Boeing 787 ©BOEING (https://www.comsol.com)

Figure 5

 5 Figure 5 Common PMCs manufacturing processes in industry

Figure 7

 7 Figure 7 Typical automated tape laying system © MTORRES (www.mtorres.com/)

Figure 8

 8 Figure 8 ATL for Boeing787's unique one-piece composite barrel construction ©BOEING (https://www.boeing.com)

Figure 9

 9 Figure 9 Schematic of automated tape laying and automated fiber placement

Figure 10

 10 Figure 10 Automated fiber placement process for fuselage section producing with window cut-outs © MTORRES (www.mtorres.com/)

  Figure 11 (a) Cincinnati VIPER™ platform for automated fiber placement, © Cincinnati VIPER™ (b) TORRESFIBERLAYUP machine for automated fiber placement, © MTorres.

Figure 12

 12 Figure 12 Robot-based AFP system developed by Coriolis Composite © Coriolis Composite (www.coriolis-composite.com)

Figure 13

 13 Figure 13 (a) Robot-based AFP system with modular head developed by Electroimpact, ©Electroimpact (b) Robot-based AFP system in the Composites Technology Center of NASA, ©NASA

Figure 14

 14 Figure 14 Redundant robotic system for application example

  Figure 17 (a) visualization of the GA progress for considered application example; (b) Joint coordinates of the positioner from GA based solution.

2. 1

 1 LAY-UP TASK DESCRIPTION AND ITS CREATION IN CAD SYSTEM ......................... 33 2.1.1 Representation of the Lay-up Task in the Form of the Sequence of 4×4 matrices............... 33 2.1.2 Conversion of the Lay-up Task to the Sequence of 6×1 vectors .......................................... 35 2.2 ROBOTIC LAY-UP SYSTEM COMPONENTS AND THEIR MODELS ............................. 36 2.2.1 Kinematic Model of Robotic Manipulator ........................................................................... 36 2.2.2 Kinematic Model of Actuated Positioner ............................................................................. 48 2.2.3 Integrated Kinematic Model of the Robotic Lay-up System ................................................ 51 2.3 MOTION GENERATION IN ROBOTIC LAY-UP SYSTEM ................................................. 53 2.3.1 Formalization of Motion Planning Problem for Robotic Lay-up System ............................ 54 2.3.2 Additional Constraints from the Actual Robotic System ..................................................... 56 2.4 SUMMARY .................................................................................................................................... 65

Figure 19

 19 Figure 19 Typical fiber placement curves defined in 3D CAD model (SolidWorks), with position and normal direction relative to the workpiece frame ©AUTOMATED DYNAMICS ® (www.automateddynamics.com).

Figure 20

 20 Figure 20 Definition of task frames with respect to workpiece frame

F

  . There are many ways to define the orientation angles. In this work, we use the definition in most of the commercial industrial robot, where the angles are defined as Z-Y-X Euler angles.

FL

  Y plane describes the floor and Z-direction is chosen along axis #1.1 Frame 1 is attached to the first link. The distance between origins of in Zdirection. The X-direction is pointed along the link 1 L , and it has the same direction of 0

5 Y

 5 with the same direction. If the tool center point is moved, the tool frame is moved with it.

Figure 21 2 L 2 F origin and 3 F

 21223 Figure 21 Architecture of standard 6-axis industrial robot

2 L 2 X

 22 define the lengths of the manipulator links 1 and 2 along the axes 1 X and . It is worth mentioning that in general case, the link parameters must include additional parameter describing the length in the direction 1 Z but it is reasonable to include it in the above presented parameter 0 L (it allows us to simplify the model and reduce unnecessary computations while keeping the model correctness).

T

  Also it includes translational parameters 3 L and d describing geometry of the link 3. For the remaining axis the desired matrices can be obtained via elementary rotations about the axes 5

Figure 22

 22 Figure 22 Different configuration of robotic manipulator for the same location of the end-effector

  the Cartesian coordinates of the wrist center point relative to the robot base frame Rbase F that can be directly extracted from the last column of the matrix 6

Figure 23

 23 Figure23Three types of singular configurations in KUKA kinematic system, ©KUKA[START_REF] Kuka | [END_REF]

1Z

  and wrist center point as shown in Figure24. For this projection, let us consider the first triangle with the edges

Figure 24

 24 Figure 24 Projection of the manipulator links on the vertical plane (Configuration SHOULDER_FORWARD & ELBOW_UP)

Figure 25

 25 Figure 25 Multiple manipulator configurations for the same wrist center point position

Figure 26

 26 Figure 26 Triangles highlighted for computing the angle of axis #2

2 

 2 , the  angle is equal to either

5 S

 5 is taken into account via the configuration index

  using other pairs of the orientation matrix elements, such as )

Figure 27

 27 Figure 27 Typical one-axis positioners from KUKA (www.kuka.com)

Figure 28

 28 Figure 28 Typical two-axis positioners from KUKA (www.kuka.com)

Figure 29

 29 Figure 29 Typical linear track and gantry from KUKA (www.kuka.com)

PbaseFF

  Positioner base frame. This frame is attached to the positioner base and does not move. It can be considered the reference frame where the positions of all other link frames are described. In Pbase F , X-Y plane describes the floor and Z-axis is vertical and outward the floor. The definitions for one-axis positioner and two-axis positioner are the same, as is shown in Figure30(a)(b). M Positioner intermediate frame. This frame is attached to the first rotational axis of the two-axis positioner, as is shown in Figure 30(b). The distance between origins of M F and Pbase F

Figure 30

 30 Figure 30 Kinematics of typical one-axis positioner, two-axis positioner and linear track

  two presentations with respect to the world frame World F

  in Figure31, the desired tool locations can be expressed in two different ways. From the robot perspective, the homogeneous 4×4 matrix describing the tool location may be computed as

  -axes of the tool and task frames are opposite, i.e.

T

  describes the i th task point in the workpiece frame W F .

Figure 31

 31 Figure 31 Integration of the robot and the positioner kinematic models

Figure 32

 32 Figure 32 Schematic of motion generation in robotic lay-up system

  To ensure that there is no collision between them, the following sufficient condition should be satisfied

Figure 33

 33 Figure 33 Collision detection based on intersection of two cylinders

  differentiable function, the minimum occurs either at the boundary of the square

2 

 2 are computed assuming that only one corresponding derivate is equal to zero. The latter yield the following critical points

Figure 34

 34 Figure 34 Cylinder-based descriptions of the workcell components for collision detection

Figure 35

 35 Figure 35 Triangle-based description of the workcell components for collision detection

P

  is one of the points of this line that can be easily found. Then it is necessary to find four parameters 

Figure 36

 36 Figure 36 Situation of two triangle planes intersecting at a straight line

Figure 37

 37 Figure 37 Structural presentation of the robotic lay-up workcell in CATIA, and component pairs selected for collision detection

3. 1

 1 PROBLEM TRANSFORMATION INTO DISCRETE FORM ............................................... 67 3.1.1 Search Space Discretization for One Redundant Variable ................................................... 67 3.1.2 Search Space Discretization with Two Redundant Variables .............................................. 73 3.2 MOTION PLANNING METHOD BASED ON COMBINATORIAL OPTIMIZATION ..... 77 3.2.1 Straightforward Approaches and Related Difficulties .......................................................... 77 3.2.2 Motion Generation using Dynamic Programming Technique.............................................. 79 3.3 ENHANCEMENT OF THE MOTION PLANNING ALGORITHM ...................................... 82 3.3.1 First Strategy: Progressive reduction of the discretization step............................................ 82 3.3.2 Second Strategy: Smoothing the redundant variable profiles............................................... 84 3.3.3 Performance Evaluation of the Enhanced Algorithms ......................................................... 86 3.4 PARAMETERS TUNING FOR MOTION PLANNNING ALGORITHM ............................. 91 3.4.1 Influence of Discretization Step on the Generated Trajectories ........................................... 91 3.4.2 Practical Recommendations for Selecting the Discretization Step ....................................... 93 3.5 SUMMARY .................................................................................................................................... 95

Figure 39

 39 Figure 39 Transformation of the motion planning problem into discrete form (case of one redundant variable)

  Taskinto the joint space. The procedure is composed of two basic steps. At the first step, the redundant variable is discretized in the interval

}

Figure 40

 40 Figure 40 Directed graph describing the motion planning problem with one redundant variable

  each row represents the same value of the positioner joint coordinate ) ( k P q . In accordance with the physical sense, the initial and final path nodes must belong to the sets

  robot joint coordinates. The latter allows us to present the objective function (total motion time) as the sum of the edge weights Notations: of robot tool and positioner flange in local frames Invoked functions: Robot inverse kinematics

Figure 41

 41 Figure 41 Motion synchronization in a typical robot controller

Figure 43

 43 Figure 43 Transformation of the motion planning problem into discrete form (case of two redundant variables)

  procedure also contains two steps. Firstly, the redundant variables are discretized in

  , and the robot configuration states corresponding to

}

  location, positioner and robot joint coordinates Matrices of robot tool and positioner flange in local frames Invoked functions: Robot inverse kinematics

Figure 45

 45 Figure 45 Directed graph with two virtual vertices for the motion planning problem

.

  At the final step, after selection of the minimum lengths corresponding to the end task point and applying the

Figure 46

 46 Figure 46 Combination of the rough search and the local optimization

  of robot tool and positioner flange in local frames order -Polynomial order Invoked functions: Polynomial approximation Polynomial approximation of the redundant variable profile Set

*

  Obtained using the developed one-stage DP algorithm without acceleration constraints † Estimated using polynomial approximation of the algorithm complexity (leading term) †

Figure 49

 49 Figure 49 Enhanced algorithm #1: trajectories obtained at the Stage I (Δq P = 3°, T motion = 9.44s)

Figure 51

 51 Figure 51 Enhanced algorithm #2: trajectories obtained at the Stage II (Δq P = 1°, T motion = 8.47s)

Figure 53 A

 53 Figure 53 A benchmark task for investigating the influence of Δq P

Figure 55 A

 55 Figure 55 A portion of the task graph corresponding to the task locations 17, 18 and 19

  ), and the task point displacement due to positioner can be approximately expressed as   max rwhere r max is the furthest point distance to the positioner axis.As follows from the physical sense, to avoid undesired intermittent positioner rotations, , since in the Cartesian space the positioner velocity is usually smaller than the velocity of the robot. It can be also rewritten in the form max r s   

  PROCESS PREPARATION FOR HIGH-PRESSURE VESSEL FABRICATION................................................................................................................................... 97 4.1.1 Manufacturing Task Description and its Regularization ...................................................... 97 4.1.2 Selection of Robotic System Components and Workcell Lay-out Design ......................... 100 4.2 ROBOTIC SYSTEM MOTION PLANNING USING DEVELOPED ALGORITHMS ...... 103 4.2.1 Generation of the Task Graph for the Thermoplastic Tape Lay-up ................................... 103 4.2.2 Generation of Time-optimal Coordinated Motions of Robot and Positioner ..................... 108 4.3 OPTIMAL MOTION IMPLEMENTATION IN ROBOTIC SYSTEM ................................ 112 4.3.1 Motion Implementation using KRL Robot Programming Language ................................. 112 4.3.2 Experimental Evaluation of the Implemented Time-Optimal Motions .............................. 116 4.4 SUMMARY .................................................................................................................................. 120

Figure 58 A

 58 Figure58A sample of the laying task for a high-pressure vessel, ©CETIM.

Figure 59

 59 Figure 59 Discrete representation of the lay-up task in the form of the augmented line

Figure 60

 60 Figure 60 Regularization of the lay-up task discrete model (smoothing and replacing an irregular step 1.5 to 21.0 mm by the regular step 8.0 mm)

Figure 61

 61 Figure 61 SPIDE-TP platform for manufacturing high-pressure vessels, ©CETIM.

Figure 62

 62 Figure 62 Technological tool used in SPIDE-TP platform, ©CETIM.

Figure 63

 63 Figure 63 Work envelop of robot KUKA KR210 R3100 ultra, ©KUKA GmbH

Figure 65

 65 Figure 65 Robotic platform SPIDE-TP and arrangement of coordinate frames

Figure 66

 66 Figure66Task graph evolution after applications of kinematic and collision constraints ("blue"-admissible vertices; "white"inadmissible vertices; "red"admissible path).

Figure 67

 67 Figure 67 Collisions between different components of robotic platform SPIDE-TP

Figure 70

 70 Figure 70 Bottleneck areas at the task graph (too close to collisions)

( a )

 a Description via the joint coordinates P = {A1 10, A2 -80.6, A3 -50, A4 0, A5 14.2, A6 0, E1 -3000, E2 200} (b) Description via the end-effector location P = {X 12.3, Y 100.0, Z 50, A 9.2, B 50, C 0, E1 -3000, E2 200, S 'B010', T 'B1010'}

Figure 73

 73 Figure 73 Description of the robotic platform configuration in KR C4 controller

Figure 74

 74 Figure 74 Approximate positioning for PTP-PTP and LIN-LIN motion sequences

Figure 75

 75 Figure 75 Curved path implementation using LIN-LIN sequence and SPLINE block

Figure 76

 76 Figure 76 Description of curved path using spline interpolation

  Figure79The KRL program of the time-optimal motion with time settings at path segments

Titre:

  Planification des mouvements optimaux dans les systèmes robotiques redondants pour un processus d'enroulement filamentaire composite automatisée Mots clés : Système robotique redondant, Planification de mouvements, Trajectoire temps-optimale, Programmation dynamique, Enroulement filamentaire Résumé : La thèse traite de la planification des mouvements optimaux dans les systèmes robotiques redondants pour l'automatisation des processus d'enroulement filamentaire. L'objectif principal est d'améliorer la productivité des cellules de travail en développant une nouvelle méthodologie d'optimisation des mouvements coordonnés du robot manipulateur, du positionneur de pièce et de l'unité d'extension de l'espace de travail. Contrairement aux travaux précédents, la méthodologie proposée offre une grande efficacité de calcul et tient compte à la fois des contraintes technologiques et des contraintes du système robotique, qui décrivent les capacités des actionneurs et s'expriment par les vitesses et accélérations maximales admissibles dans les articulations actionnées. La technique développée est basée sur la conversion du problème continu original en un problème combinatoire, où toutes les configurations possibles des composants mécaniquessont représentées sous la forme d'un graphe multicouche dirigé et le mouvement temporel optimal est généré en utilisant le principe de programmation dynamique. Ce mouvement optimal correspond au plus court chemin sur le graphique satisfaisant les contraintes de lissage. Les avantages de la méthodologie développée sont confirmés par une application industrielle d'enroulement filamentaire pour la fabrication de pièces thermoplastiques au CETIM.Title : Optimal motion planning in redundant robotic systems for automated composite lay-up process Keywords : Redundant robotic system, Motion planning, Time-optimal trajectory, Dynamic programming, Automated composite lay-up

Table 16

 16 Spatial parameters defining the emplacement of the robotic system components ............... Table17Geometric parameters of transformation describing the robot and positioner kinematics ... Table 18 Definition of the configuration index S for 6-axis KUKA robot (Status index) .................. Table 19 Definition of the configuration index T for 6-axis KUKA robot (Turn index) .................... Table 20 Joint limits and maximum velocities/accelerations for the robot and positioner ................. Table 21 Comparison of motion commands of KR C4 controller for curved path implementation ...

Table 1

 1 Constitutes of fiber reinforced composite materials The individual components remain separate and distinct within the finished structure, and the new combined materials have improved properties over the individual materials

	Reinforcing phase	Matrix phase	Interface
	Glass
		Polymer
	Carbon
		Metals	Bonding
	Ceramic
		Ceramics	Surface
	Metallic …	…
	Such new combined materials usually have characteristics different from the individual components.

Table 2

 2 Composite classifications on the basis of matrix phase

		Polymer matrix composite	Ceramic matrix composite	Metal matrix composite
	Reinforcing	Glass, carbon, graphite and aramid	Carbon and aluminum oxide	Graphite, alumina, silicon carbide…
	Main Properties	High strength-to-weight ratio High stiffness-to-weight ratio	Resistance to high temperature and corrosive environment	High modulus of elasticity, ductility, and resistance to elevate temperature
	Applications	Automotive, naval, aeronautical and aerospace…	Jet and automobile engines…	Satellite, missile, space structures…

Table 3

 3 Classification of PMCs processes

		Thermo-set composites processing	Thermoplastic composites processing
	Continuous-fiber composites	Lay-up (manual/automated) Filament winding Resin transfer molding (RTM) Pultrusion	Lay-up (manual/automated) Thermoforming Compression molding
		Spray-up
	Short-fiber	Injection molding	Injection molding
	composites	Compression molding
		Resin transfer molding (RTM)

Table 4

 4 Properties of existing PMCs manufacturing process

	Processes	Properties
	Lay-up	Suitable for a wide variety of products No complicated equipment required (low cost on tooling)
	Spray-up	Only with short fibers
	Vacuum bag molding	An extra process adds cost both in labor and in bagging materials
	Compression molding
	Resin transfer molding	Matched complicated tooling is always required
	Injection molding

Pultrusion

Filament winding

Limited to component shapes (constant cross-section and cylindrical, respectively)

Table 5

 5 Joint limits and velocity limits of the redundant robotic system

	Joint limits	Velocity limits

Table 6

 6 Main optimization-based techniques for redundancy resolution in literature

Table 7

 7

	Parameters setting of ga(•) function in Matlab
	Population Size	200
	Fitness scaling	Rank
	Selection	Stochastic uniform
	Elite count	10
	Mutation	Gaussian (scale =1; shrink =1)
	Crossover	Scattered function; Fraction =0.8
	Migration	Forward (fraction= 0.2, Interval= 20)
	Generation	MaxGenerations = 5000; MaxStallGenerations = 100

Table 9

 9 Software packages for computer aided design of industrial robotic systems

	Software Package	Company	Principle applications
	Robcad	Tecnomatix www.siemens.com/tecnomatix	Welding, cutting, drilling, riveting, painting, spraying, etc.
	DELMIA (IGRIP)	Dassault Systems www.3ds.com	Welding and material handling applications.
	Robotmaster	Hypertherm www.robotmaster.com	Cutting, 3D machining, trimming, welding, polishing, dispensing, de-burring, painting, spraying, etc.
	Workspace5	WAT Solutions www.workspace5.com	Welding; Waterjet cutting.
	CATFiber for CATIA	Coriolis Composite	Fiber placement; Ultrasonic Trimming; Non-
	CADFiber for NX	www.coriolis-software.com	destructive testing.
	FPM/FPS	Automated Dynamics www.automateddynamics.com	Fiber placement; Tape laying
	KUKA.Sim	KUKA www.kuka.com	Welding, palletizing, coating, painting, machining.
	RobotStudio	ABB

http://new.abb.com Application software for cutting, machining, welding and palletizing.

Table 11

 11

		Kinematic and dynamic constraints describing capabilities of robot KUKA KR210
	AXIS	Joint limits max q	Maximum joint velocity max q 	Maximum joint acceleration max q  
		R	R	R
	1	185° to -185°	105°/s	105°s -1 /0.25s
	2	-5° to -140°	101°/s	101°s -1 /0.25s
	3	155° to -120°	107°/s	107°s -1 /0.25s
	4	350° to -350°	136°/s	136°s -1 /0.25s
	5	122.5° to -122.5°	129°/s	129°s -1 /0.25s
	6	350° to -350°	206°/s	206°s -1 /0.25s

Table 12

 12 Kinematic constraints of redundant robotic system for the benchmark example

	Joint limits	Maximum velocity	Maximum acceleration

Table 14

 14 Simulation results for linear lay-up task for different discretization step Δq P Figure 54 Evolution of robotic system motion time with reduction of different discretization step Δq P

		 q		 00 . 2	 q		 00 . 1	 q		 75 . 0	 q		 50 . 0	 q		 25 . 0	 q		 10 . 0
		P			P			P			P			P			P
	T motion	1.897 s	1.836 s	1.540 s	1.298 s	1.290 s	1.290 s
	without a-constr.	(~38s comp.)	(~2min comp.)	(~4min comp.)	(~9min comp.)	(~47min comp.)	(~5h comp.)
	T motion	1.897 s	2.107 s	1.594 s	1.300 s	1.290 s	1.290 s
	with a-constr.	(~67s comp.)	(~4min comp.)	(~8min comp.)	(~17min comp.)	(~1.2h comp.)	(~9h comp.)

Table 15

 15 Robotic system motion time with respect to different robot base location on the linear track.

	Robot location (mm)	2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000
	Motion time (sec)	4.55	4.36	4.27	4.32	4.19	4.11	4.19	4.31	4.41	4.58	4.86

Table 18

 18 Definition of the configuration index S for 6-axis KUKA robot (Status index)

	Value	Bit 2		Bit 1		Bit 0
	0	Wrist down	Elbow down	Shoulder forward
	1	Wrist up	Elbow up	Shoulder backward
	Table 19 Definition of the configuration index T for 6-axis KUKA robot (Turn index)
	Value	Bit 5	Bit 4	Bit 3	Bit 2	Bit 1	Bit 0
	0

Table 21

 21 Comparison of motion commands of KR C4 controller for curved path implementation

	Motion commands	PTP-PTP	LIN-LIN	SPLINE block
	Interpolation space	Joint space	Cartesian space	Cartesian space
	Positioning	Approximate	Approximate	Exact
	Time assignment	Speed/Acceleration settings	Speed/Acceleration settings	Speed/Acceleration & Time settings

  P[1] = {X19.62, Y 17.12, Z 1027.28, A 41.11, B 80.86, C 90.13, E1 -3000.00, E2 -138.89, S 2, T 2} P[2] = {X 33.61, Y 10.09, Z 1028.84, A 16.70, B 78.18, C 47.87, E1 -3000.00, E2 -120.30, S 2, T 2} … P[54] = {X 13.42, Y 23.68, Z 1027.53, A 60.42, B 80.95, C 106.93, E1 -3000.00, E2 2.01, S 2 , T 2}

ACKNOWLEDGEMENTS

Before beginning this manuscript, I would like to express my gratefulness to people who have guided and helped me in my research and who have provided all possible assistance throughout the time of working on the thesis.

Therefore, using sampling of the redundant variables, the original continuous problem is transformed into the discrete form that already incorporates the collision and singularity constraints. It allows us to apply further some combinatorial optimization techniques that should be modified to take into account the remaining velocity/acceleration constraints describing capacities of the actuators. This issue is in the focus of the next section.

Algorithm 1b: Search space generation G 2 (.) Input:

Task location matrices backtracking, one can get the desired optimal path. Therefore, the desired solution for n sets

is obtained by increasing sequentially the i-index, computing . This allows us to present the algorithm basic expression in a similar way

denotes the length of the shortest path connecting one of the initial vertices

In more details, an outline of the developed algorithm for the case of one redundant variable is presented below (see Algorithm 2). Here, the input is the matrix of the locations

, which contains information on the configuration states satisfying the equality constraint (68), the collision constraint (86) and the singularity constraint (90). The algorithm operates with two tables ) and the pointers to the previous locations respectively. The procedure is composed of four basic steps. The step (1) initializes the distance and pointer matrices by defining their first columns. In step [START_REF] Mäkinen | The development of sandwich structures for naval vessels during 25 years[END_REF], the recursive formula (108) is implemented. The computing starts from the second column and tries all possible connections between the nodes in the current column

 and the previous one

. For the admissible configuration states, the acceleration constraints are examined using the expression (103) for each candidate path connecting the nodes with the indices to the previous location in the current path. Then, sub-step (2b) finds the minimum path from the current node

 and records the reference to

into the pointer matrix. In steps (3) and ( 4), the optimal solution is finally obtained and corresponding path is extracted by means of the backtracking.

To demonstrate the efficiency of the proposed methodology, let us apply it to the benchmark example that deals with a redundant robotic system composed of a 2-axis planar manipulator and a 1axis positioner, which is described in details in Subsection 1.2.1 (see Figure 14). This system possesses a 1-dof redundancy with respect to the ellipsis-type planar lay-up task that was sampled in 100 segments uniformly. The segment nodes should be visited sequentially by the end-effector in minimum time by using in the best way the motion capabilities of both the robot and the positioner, which are described by the constraints presented in Table 12. The developed algorithm was implemented and examined using Matlab in Win7SP1 with Intel ® i5 @2.67 GHz 2.67 GHz.

Algorithm 2: DP-based Path planning DP(.) Input:

Matrix of locations

-distance and pointer n m matrices Invoked functions: Distance between nodes

At the stage II (named as the iterative local fine search), a secondary discretization with smaller step is applied in the neighborhood of the initial solution obtained at the stage I. The neighborhood size is defined by the user (see Figure 47). It is clear that this local discretization produces the search space of relatively small size (compared to straightforward approach) allowing us quickly generate an optimal solution that is obviously better than the initial one. For instance, for the above mentioned benchmark example, it takes about two minutes in total (stages I and II) while the one-stage technique (see subsection 3.2.2) requires more than 10 minutes to obtain a similar result. It should be mentioned that the stage II can be repeated several times, sequentially reducing the discretization step.

. At the first stage, the

is generated using the full-range discretization with the density m (step Ia). The initial solution

obtained by applying the path planning function

is generated in the neighborhood of the initial solution (step IIb). After applying again the motion planning function

Second Strategy: Smoothing the redundant variable profiles

An alternative way to speed up the motion generation algorithm is to replace the stage II in the first strategy by simple smoothing the curve

corresponding to the redundant variable without changing the time instants obtained from the stage I (if it is acceptable from technological point of view). This can be achieved by applying the polynomial approximation to the function defined by the nodes

. From our experience, the 3 rd order polynomials are quite satisfactory for this procedure, which generates modified values of the redundant variable by the end-effector, the robot joint coordinates should be also regenerated using the new values of the redundant variable

and sequentially reapplying the positioner direct kinematic and robot inverse kinematic transformations (see subsection 2.2.1 and subsection 2.2.2). This procedure is shown in Figure 48. It is worth mentioning that the smoothing-based technique is very time-efficient because it excludes the optimization at the stage II. However, in most cases studied by the authors the two-stage technique including the local optimization provided better results. An algorithm outline for the second strategy is presented below (see Algorithm 4). It is obvious that the Stage I of this algorithm is exactly the same to that in previous subsection, so the Stage II dealing with the smoothing of the redundant variables is presented only. The input includes the task location

Algorithm 4 (Stage II only): Smoothing of redundant variables Input:

Task location matrices

, and the initial trajectory obtained from the Stage I (where the redundant variable profiles will be used). The output contains the smoothed trajectories of redundant variables and the robot joint coordinates

. A polynomial fitting is firstly applied to the initial trajectory of the redundant variable at the step IIa. Then, the smoothed

is used to update the robot trajectory taking into account the positioner direct kinematics and robot inverse kinematics (step II).

Performance Evaluation of the Enhanced Algorithms

In order to estimate computational efficiency of the developed algorithms and evaluate smoothness of the obtained motions, there were carried out a number of computational experiments that deal with the benchmark example adopted in this work (robotic system composed of 2-axis manipulator and 1axis positioner, 100 task locations). In these experiments, several algorithms were compared, including the conventional one (based on Dijkstra's shortest path), the DP-based algorithm with fine discretization and also enhanced algorithms with local fine optimization or with smoothing at the second stage. The computing time was estimated in the Matlab R2014b environment running in Win7SP1 (with Intel ® i5 CPU 2.67GHz 2.67GHz, 8G memory). It is clear that for C++ implementation the computing time is essentially smaller but it does not influence on the relation between the computing times for different algorithms.

tool frames are defined in slightly different way compared to the Section 2.2, which is imposed by the industrial partner. The robot and positioner kinematic parameters are presented in Table 17 where 1 R q ... 6 R q and P q denote corresponding joint coordinates. It should be mentioned that in contrast to other works where the kinematic models are based on the DH parameters; here another equivalent parameterization is used that is more convenient and computationally efficient. Definitions of the corresponding frames are presented in Section 2.2. Besides, it is worth mentioning that following the standard of the robot manufacturer, the manipulator configuration was defined using the indices S and T [START_REF] Kuka | [END_REF] where S allows to avoid ambiguities in the first, third and fifth joint angles in a usual way while T

Generation of Time-optimal Coordinated Motions of Robot and Positioner

The obtained task graph contains all possible postures of the robotic manipulator and positioner, which ensure following the laying path without collisions, without violation of the joint limits and keeping safe distance to the singularities. This graph includes 24089 vertices and is organized as a matrix with 361 rows and 159 columns where some of the cells are empty. Using the task graph and applying the developed algorithms presented in Chapter 3, it is possible to generate the desired timeoptimal motions of the robotic system for implementing the given laying path.

While generating the time-optimal motions, the algorithms require some additional constraints defining the actuators capacities to move the manipulator and positioner. These capacities are expressed in the form of the maximum allowable velocity and acceleration for each actuated joint that are presented in Table 20. It is worth mentioning that in practice the acceleration constraints are defined in the robot controller indirectly, via the special parameter (acceleration time) that it was assumed here to be equal to 0.25 s but it can be modified by the user depending on the payload. For this reason, during the motion planning the acceleration constraints were applied in "soft" manner, with different acceleration times. Another reason justifying this operation is related to approximate evaluation of the acceleration in the optimization algorithm that is based on the finite difference technique that uses three time instances only, in accordance with expression (103). As it has been observed, strict application of the acceleration constraints may produce some undesired peaks on the motion profiles that disappear if the acceleration constraints were slightly weaken (relaxed). However, detailed analysis of the final trajectories generated in such way shows that the obtained motions are smooth enough and satisfy all considered constraints. To find the desired time-optimal motions, the developed DP-based algorithm with the discretization 1 P q    was applied first and the trajectories were smoothed further using the cubic spline approximation (Algorithms 2 and 4, see Section 3.3). It took about one hour of computations for searching of the best path on the task graph and 0.1 sec for the smoothing. It is worth mentioning that further reduction of the discretization step is extremely time consuming, so the Algorithm 3 (with the local optimization) was not applied for this practical problem where the collision check requires more than 10 hours even for 1 P q    . Nevertheless, the results were quite acceptable for practice. In particular, for the generated trajectory the robotic system motion time is about 4.0 sec, which is much better compared to 14.0 sec that our industrial partner got using the software package Composicad.