the airport before my first flight to France, our wedding, our trips,. . . but also the difficult moments that we have overcome together. I appreciate all of them as a part of our adventure.

AYERBE and Professor Didier DUMUR at the L2S laboratory, CentraleSupélec; Associate Professor Sylvain LAVERNHE and Professor Christophe TOURNIER at the LURPA laboratory, ENS Paris Saclay. Over my three-year PhD thesis, experiencing a lot of meetings and discussions, you all have given me precious knowledges and excellent supervisions, guiding me to accomplish the thesis. Thank you for your patience, thoroughness, belief, and encouragement. I appreciate deeply all of your supports, not only in the scientific point of view but also for my personal life.

Thank you so much! It is a great honor for me that Professor Kaan ERKORKMAZ at University of Waterloo, Canada and Associate Professor Richard BÉARÉE at Arts & Métiers ParisTech have accepted to be two reporters of my PhD thesis. Thank you for giving me many good comments, excellent proposes and recommendations to improve my thesis and to extend perspectives for future works. I would like also to thank Professor Emmanuel DUC at SIGMA Clermont, and Associate Professor Mohamed ABBAS-TURKI at ENS Paris-Saclay, for accepting to be the jury members of my thesis defense. Thank you so much for your very good questions and comments about my PhD thesis.

It is a great chance for me to thank the L2S and LURPA laboratories, where I have worked in a rich scientific environment, and have enjoyed friendly vivid activities. I would like also to acknowledge the Digiteo foundation for accepting and funding my 3-year PhD thesis.

I would like to say thank you for all of my friends and colleagues: Quang Binh, Lam Thanh, Minh Nguyen, Minh Tri, Son Nano, Dai Hai, Quang Bang, Duy Minh, Abid, Mohammed, Michele, Lorène, Kévin, Kamel, Laureen,... Thanks for your helps, for our good memories, picnics, trips together.

It is not only today that I would like to express my profound gratitude to my parents, who's giving me the best of their lives. Thank you so much! Mommy, Daddy, I love you! Finally, I dedicate this success for you, Phuong. You are my friend, my wife, my sources of motivation and love. I could not succeed this PhD thesis without you. At this moment, I remember not only our happy moments, our first favorite food "vit om sau", our phone call at En raison de la forte concurrence dans les industries manufacturières, il y a de plus en plus de demandes pour la mise en oeuvre des stratégies d'usinage à grande vitesse. Pour augmenter les performances de l'usinage à grande vitesse, ce travail de doctorat se concentre principalement xviii Variable gain contouring control for multi-axis machine tools sur l'amélioration de l'étape de commande d'axe, en particulier la commande en position des axes par la commande numérique. La motivation générale est que si l'action de commande fonctionne bien, la pièce usinée sera très proche de celle désirée. Par conséquent, les étapes de polissage manuelles après usinage peuvent être réduites. La qualité de coupe, la productivité et les avantages économiques sont donc évidents.

Le cadre de ces travaux de recherche est mis en évidence dans la Fig. 2. Le processus est constitué de multiples entrées et sorties.

Pièce:

-propriétés du matériau, -taille, géométrie, etc

Outil:

-propriétés du matériau, -géométrie, etc Processus d'usinage général en MO

Paramètres de coupe:

-vitesse de coupe, -profondeur de coupe

Machine-outil:

-spécification broche, -comportement dynamique -énergie, etc

CNC: -interpolation temporelle

Parties du processus d'usinage en MO considérées dans ce travail de recherche

Commande d'axe (CNC)

Perturbations:

-frottement -force de coupe, etc

Incertitudes:

-variation des paramètres, etc Taux d'enlèvement de matière Température de l'outil et de la pièce

Défauts de la pièce usinée:

-déformation, -rugosité, -erreur de contour, etc

Réponses de coupe:

-force de coupe, -énergie de coupe, etc Dans ce processus, les paramètres de coupe (vitesse de coupe, profondeur de coupe, etc), configuration et spécification de la machine-outil (spécification broche, comportement dynamique, etc), de la CN (interpolation temporelle), de la pièce (propriétés du matériau, taille, géométrie, etc) et de l'outil (propriétés du matériau, géométrie, etc) sont des données d'entrée. Tandis que les effets de sortie sont le taux d'enlèvement de matière, la température de l'outil et de la pièce, les défauts de la pièce usinée, les forces de coupe, etc. La partie centrale du processus d'usinage sur machine outil est la commande d'axe réalisée en boucle fermée, étant sujette à des incertitudes (par exemple la variation des paramètres d'entraînement, etc) et à des perturbations (par exemple frottements, force de coupe, etc). Dans ce travail de recherche, cinq facteurs d'entrée sont principalement pris en compte, notamment la vitesse de coupe, le comportement dynamique de la structure de la machine outil, la vitesse d'avance, la géométrie de la pièce/surface désirée et l'outil. Seul le frottement est considéré comme perturbation de la commande d'axe. En sortie, seul le profil d'erreur de contour est pris en compte.

Pour améliorer les performances de la CN industrielle, les techniques d'interpolation temporelle, le processeur, la motorisation et l'actionneur ont été considérablement améliorés. Toutefois, sa structure et sa stratégie de commande n'ont eu que peu d'améliorations. En fait, chaque entraînement d'axe est contrôlé à travers une structure classique en cascade, comprenant des boucles de position, de vitesse et de courant, respectivement des boucles externes aux boucles internes. Tandis que les deux boucles internes sont contrôlées par des régulateurs proportionnel -intégral (PI), la boucle de position externe est contrôlée par un régulateur Proportionnel (P) combiné avec une action d'anticipation en vitesse (Feed Forward -FFW) pour la boucle de vitesse [Altintas 2000a].

Au sein du processus de commande, il est évident que les petites discontinuités dans la trajectoire peuvent exciter les modes naturels de la structure mécanique, provoquant des vibrations de la machine. Par conséquent, la trajectoire de référence doit d'abord être lissée. Pour ce faire, la planification de la vitesse est réalisée en tenant compte des limites cinématiques de la machine, en particulier les limites de jerk [Erkorkmaz 2001a[START_REF] Barre | [END_REF], Beudaert 2012]. Par la suite, on s'attend à ce que le contrôleur produise la trajectoire résultante en respectant toutes les contraintes cinématiques de la trajectoire de référence.

Étant donné que l'usinage des surfaces de forme complexe est principalement concerné dans cette étude, il y a de nombreuses raisons qui expliquent qu'avec la commande classique mentionnée ci-dessus, il soit difficile d'obtenir une grande précision de contour. D'une part, dans l'usinage à grande vitesse, la vitesse de coupe élevée conduit à une vitesse d'avance élevée.

L'interpolation d'une telle vitesse d'avance augmente également la vitesse, l'accélération et le jerk de chaque axe lors de l'exécution de la trajectoire. Par conséquent, en suivant le profil libre, les axes de la machine-outil peuvent avoir des caractéristiques cinématiques élevées et des fréquences élevées. D'autre part, la structure de commande d'axe classique a une bande passante limitée. En outre, les sollicitations dynamiques de la machine-outil changent; elle est soumise à des perturbations, telles que des variations d'entraînement dues à la variation de la masse de la pièce et des conditions de lubrification, au frottement et aux forces de coupe [Altintas 2000a, Altintas 2000b, Erkorkmaz 2001b, Tang 2013, Khoshdarregi 2014]. Ces effets augmentent l'erreur de suivi et génèrent une erreur de contour, définie de manière classique comme la distance orthogonale entre le point de contact réel de l'outil et le profil souhaité.

Pour relever ces défis, les chercheurs ont proposé de nombreuses stratégies de commande avancées et intelligentes comme alternatives à la technique classique, telles que Cross Couple Control [Koren 1980[START_REF] Koren | [END_REF], Zhao 2013], Sliding Mode Control [Altintas 2000b], Adaptive Robust Control [Yao 1997, Davis 2015]. Ces approches sont plus ou moins validées par des résultats de simulation et/ou des essais expérimentaux utilisant une CN ouverte [Altintas 1994, Erkorkmaz 2001b, [START_REF] Beudaert | [END_REF]. Cependant, l'implémentation des contrôleurs avancés dans la version commerciale des CN n'est pas encore réalisée. Ceci est dû au fait que de façon générale, une CN ne permet jusqu'à présent qu'un accès limité au coeur des contrôleurs; elles restent donc "fermées". Une autre raison est que presque toutes les stratégies de commande proposées nécessitent une modification de la structure de commande classique.

En résumé, cette thèse de doctorat "Commande à gains variables de l'erreur de contour pour l'usinage multiaxes" vise à améliorer la précision de contour en usinage multi-axes grande vitesse pour des surfaces de forme complexe, en exploitant toutes les possibilités de la structure de commande classique de l'entraînement d'axes, et en proposant des solutions que les constructeurs pourraient intégrer dans leur CN commerciale. L'amélioration du suivi de la trajectoire signifie que l'erreur de contour doit être contrôlée aussi faible que possible. Dans cette étude, à l'exception du frottement, les autres sources de perturbations et d'incertitudes sont négligées.

De plus, compte tenu des vibrations de la machine, l'étape d'interpolation temporelle produit les mouvements d'axes, en respectant les contraintes cinématiques des axes, notamment les limites de jerk. Ces travaux proposent une méthode permettant de générer des gains de commande variables optimaux pour la réduction de l'erreur de contour dans le processus d'usinage.

L'optimisation proposée est ici soumise à diverses contraintes que sont les limitations cinématiques des axes (vitesse, accélération et jerk), la stabilité de l'entraînement des axes et les limites de courant du moteur.

Cette thèse de doctorat est réalisée à l'Université de Paris-Saclay, dans le cadre de la coopération Le second support de cette thèse est un algorithme d'interpolation temporelle VPOp (Velocity Profile Optimization), développé par X. Beudaert lors de sa thèse de doctorat conduite au LURPA de 2010 à 2013 [Beudaert 2013], permettant de générer les consignes de position de référence, contenant un profil de vitesse d'avance optimisé, et respectant toutes les contraintes cinématiques des axes mentionnées auparavant.

Stratégie principale

La stratégie est illustrée dans la Fig. 3.

Dans une optique d'optimisation du processus, on peut noter que le réglage de l'entraînement des axes pourrait être différent en fonction de la surface usinée. L'idée est donc d'optimiser le contrôleur d'axe en fonction de la géométrie de la trajectoire d'usinage, avec pour objectif la réduction des erreurs de contour.

En particulier, au lieu d'utiliser des gains de commande fixes dans le contrôleur de position, les gains variables seront calculés tout en préservant la structure de commande en cascade classique avec anticipation.

L'idée de départ de cette thèse de doctorat est de rechercher une séquence optimale de gains de Une autre spécificité est l'utilisation de "l'horizon glissant", approche présente dans la commande prédictive (MPC), pour résoudre le problème d'optimisation pour l'ensemble de la trajectoire.

Cette approche proposée, appelée "Offline Gain Adjustment (OGA)", permet de réduire l'erreur de contour pour l'usinage multi-axes grande vitesse, que ce soit en usinage 3 axes ou 5 axes.

Enfin, l'objectif est de développer les bases d'une méthode qui pourrait être mise en oeuvre dans une CN commerciale.

Architecture du manuscrit

L'architecture du manuscrit est organisée comme suit : Après l'introduction, le premier chapitre, "Géométrie et cinématique de trajectoire dans l'usinage à grande vitesse multi-axes", formule l'erreur de contour dans les cas d'usinage en bout et en flanc d'outil. Il met en évidence les contraintes cinématiques de la trajectoire de référence. Le deuxième chapitre, "Commande de contour adaptative pour l'usinage multi-axes", consiste à développer l'approche proposée nommée "Off-line Gain Adjustment"(OGA).

Cette approche permet de pré-compenser l'erreur de contour définie au chapitre précédent, qui se produit dans le processus d'usinage. Tout d'abord, la structure classique de commande des axes est présentée. La caractéristique non linéaire du modèle d'axe est également détaillée. Deuxièmement, l'état de l'art relatif aux techniques avancées de commande d'axes est passé en revue.

Ensuite, le coeur de ce chapitre consiste à formuler et à résoudre le problème d'optimisation sous contraintes lié à OGA.

Le but du troisième chapitre, "Simulations et discussions", est de mettre en évidence l'amélioration de la précision des contours grâce à l'approche OGA. Les effets du profil de gain variable optimal résultant, la vérification des contraintes et les effets de la configuration des paramètres d'OGA sont présentés et analysés à travers les résultats de simulation pour valider la pertinence de l'approche OGA.

Le quatrième chapitre, "Extensions OGA", consiste à explorer les applications d'OGA pour des objectifs diversifiés. Dans la première section, un gain global pour la commande en position est présenté et optimisé en utilisant l'approche OGA. Ensuite, la mise en oeuvre d'OGA dans le cas de l'utilisation d'une structure de commande prédictive généralisée est effectuée. Les résultats de simulation correspondants prouvent les possibilités offertes par OGA dans la pré-compensation de l'erreur de contour au sein de différents types de contrôleurs de position. Enfin, les "Conclusions et perspectives" visent à évaluer non seulement les contributions scientifiques mais aussi les perspectives de ces travaux.

xxiv Variable gain contouring control for multi-axis machine tools

Contributions de la thèse

Les principales contributions de cette thèse sont présentées ci-dessous.

La première contribution porte sur la formulation de l'erreur de contour en fraisage 5 axes en

Perspectives

Une perspective est de prendre en compte d'autres effets non linéaires, ainsi que des perturbations, comme par exemple les flexibilités de la structure machine, ou encore les efforts de coupe.

Introduction

Context

The advanced machining techniques are always the backbone of the manufacturing industries.

One of such techniques, High Speed Machining (HSM), is the main subject of this PhD thesis.

It replaces the Conventional Speed Machining (CSM) in sculpturing intricate parts in aerospace, automotive, die or mold industries, etc. This is because high speed machining increases material removal rates and enhances surface finish. As a result, an economical benefit can evidently be obtained.

Generally, high speed machining involves many processes, as shown in Fig. 5

Motivation and framework

Due to the high competition in the manufacturing industries, there are increasing demands for implementation of high speed machining strategies. To boost the performance of high speed machining, this PhD thesis mainly focuses on the improvement of the axis control stage, particularly the position controller of axis drive in the CNC process. The general motivation is that if the control action performs really well, the part finishing will be very close to the desired one.

Consequently, the hand finishing and polishing steps after machining can be eliminated or at least will require less effort. Hence, the cutting quality, productivity, and economical benefit can obviously be achieved.

The framework of the current research is highlighted within a multi-input and multi-output process of machining in Fig. 6.

Tool and workpiece temperature

Material removal rate

Cutting responses: -cutting forces -cutting energy

Machine tool:

-spindle specification -dynamic behavior -power -machine rigidity

Workpiece/desired profile:

-material properties -size, geometry

CNC: -feedrate interpolation

Tool/workpiece interface:

-deformation -roughness -contour error

Tool:

-material properties -geometry

Cutting parameters:

-cutting speed -depth of cut

Axis control (CNC)

Uncertainties: -parameter variations

Disturbances:

-friction -cutting force general machining procedure in MT specified machining procedure in MT within the current research In this process, the cutting parameters (e.g. cutting speeds, depth of cut, etc), configuration and specification on machine tool (e.g. spindle specification, dynamic behavior, etc), CNC (feedrate interpolation, axis control), workpiece (e.g. material properties, size, geometry, etc) and tool (e.g. material properties, geometry, etc) are designated as the inputs. Meanwhile, the output effects are such as material removal rate, tool and workpiece temperatures, deformation of machined part, cutting forces, etc. The core part of the machining process on machine tool is the axis control realized in a feedback manner, being subject to uncertainties (e.g. variation of drive parameters), and disturbances (e.g. friction, cutting force). In the limit of the current research, five input factors are mainly considered, including cutting speeds, dynamic behavior of machine tool (i.e. the axes), feedrate, geometry of desired part/surface and tool. Only friction is considered as the disturbance on axis control. At output, only contour error profile is taken into consideration.

To enhance the performance of the industrial CNC, interpolation techniques, computer processor, motor and actuator have been much improved. However, its control structure and strategy have undergone less improvement. In fact, each axis drive is controlled by CNC through a classical cascaded structure, including position, velocity and current loops ranging from outer to inner loops respectively. While the two inner loops are controlled by Proportional -Integral (PI) controllers, the external position loop is controlled by a Proportional (P) controller, combined with the Feed Forward (FFW) action for the velocity loop [Altintas 2000a].

During control process, it is known that small discontinuities in the trajectory can excite the natural modes of the mechanical structure, causing machine vibration. Therefore, the reference trajectory firstly needs to be smoothed. To do this, the trajectory planing is performed by considering the kinematic limitations of the machine, especially the limited jerk [Erkorkmaz 2001a[START_REF] Barre | [END_REF], Beudaert 2012]. Afterwards, the controller is expected to produce the resulting trajectory respecting all of the kinematic constraints of the reference one.

Considering that machining free-form surface is primarily concerned in this study, there are many reasons saying that with the above classical axis control, it is difficult to obtain high contouring accuracy. On the one hand, in high speed machining, the high cutting speed leads to high feedrate. By interpolating such a feedrate, the required axis velocity, acceleration and jerk are also increased. Thus, by following the free-form profile the machine tool axes may have high kinematic characteristics and high frequency. On the other hand, the classical axis control structure has a limited bandwidth. Furthermore, machine tool contains uncertainty characteristics and disturbances, such as variation in drive parameters due to the varying workpiece mass and lubrication condition, coulomb friction and cutting force [Altintas 2000a, Altintas 2000b, Erkorkmaz 2001b, Tang 2013, Khoshdarregi 2014]. These detrimental facts increase tracking error, and generate large contour error, which is classically defined as the orthogonal distance between the actual tool contact point and the desired profile.

To cope with the above challenges, researchers proposed many advanced and intelligent con-trol strategies as an alternative for the classical one, such as Cross Couple Control (CCC) technique [Koren 1980[START_REF] Koren | [END_REF], Zhao 2013], Sliding Mode Control [Altintas 2000b], Adaptive Robust Control (ARC) [Yao 1997, Davis 2015], Fuzzy Logic Control [Ngo 2013], Neural Network [Lin 2006], Model Predictive Control (MPC) [START_REF] Susanu | [END_REF], Tang 2012], etc. These approaches are more or less validated by simulation results and/or experimental tests using an open CNC [Altintas 1994, Erkorkmaz 2001b[START_REF] Beudaert | [END_REF]]. However, the implementation of the advanced controllers into the commercial version of CNC is not achieved yet. This is due to the fact that the CNC up to now permits only a limited access to the core of the controllers due to their close features. Another reason is that almost all the proposed control strategies require the modification of the classical control structure. That may be inconvenient and very expensive for the CNC manufacturer. Therefore, the previously mentioned classical control structure is still in use in the commercial CNC.

In summary, this PhD thesis "Variable gain contouring control for multi-axis machine tools" aims at improving the contouring accuracy in multi-axis high speed machining of free-form surfaces, by exploiting all possibilities of the classical control structure of axis drive, proposing solutions that the manufacturers can implement into their commercial CNC. The enhancement of contour following means that the contour error should be controlled to be as small as possible. In this research, except for the friction, other sources of disturbance and uncertainties are neglected. Moreover, considering machine vibration, the feedrate interpolation step produces the axis motions, respecting axis kinematic constraints, especially the limited jerk. This thesis refers to a method to generate optimal variable control gains for pre-compensating the contour error in machining process, using optimization techniques. The proposed optimization is here subject to the constraints, including the axis kinematic (velocity, acceleration and jerk) limitations, the stability of the axis drive and the motor current limits. Other factors causing the machine vibration are not taken into consideration. The second support of this thesis is a velocity optimization profile (VPOp) algorithm, developed by X. Beudaert during his PhD thesis conducted at LURPA from 2010 to 2013 [Beudaert 2013], allowing the generation of reference position setpoints, containing an optimized velocity profile, and respecting all the axis kinematic constraints.

This

Strategy

The main strategy of this PhD thesis is illustrated in Fig. 7. It is observed by CNC users that axis drive tuning can be different in function of the machined surface. Thus, the idea is to optimize the axis controller depending on the geometry of the machining trajectory, with the objective of contour error reduction.

Particularly, instead of using fixed proportional and feed forward control gains in the position controller, variable values of such gains will be computed while preserving the classical cascaded control structure with feed forward.

The original idea of this PhD thesis is to seek an optimal set of variable control gains that could excite appropriately the axis drive dynamic to reduce the contour error during machining process.

Control gains are generated off-line by solving an optimization problem and using the nonlinear axis model of the machine. Moreover, the optimization problem is subject to constraints, including axis kinematic limitations, stability criterion of servo drives and motor current limits.

Moreover, the contour error is considered in both 5-axis point milling and 5-axis flank milling.

The off-line execution is based on the nonlinear axis model, in which the disturbance coming from friction model is considered.

It is well known that when the friction is concerned, the feed forward friction compensation is widely used in CNC system to overcome the detrimental effects of friction at motion reversals, especially for contouring accuracy. However, the current nonlinear axis model, served as the fundamental base of this thesis, has not been equipped yet such a compensation technique.

Therefore, the control law for contour error reduction is proposed in this thesis without regard to this situation as well. It is assumed that the feed forward friction compensation is well designed in CNC and the proposed control law serves as an advanced supplementary approach to improve more the contouring accuracy in CNC machining.

Another specificity is the use of "receding horizon" approach of Model Predictive Control (MPC) to solve the optimization problem for the whole trajectory. This proposed approach, so-called Off-line Gain Adjustment (OGA), can be used to reduce the contour error for multi-axis high speed machining, either 3-axis or 5-axis machining.

Ultimately, the objective is to build the foundations of a method that could be considered for an implementation in a commercial CNC.

Thesis organization

The manuscript architecture is summarized as follows:

After the introduction, the first chapter, "Tool path geometry and kinematics in multi-axis high speed machining", is to formulate the contour error in multi-axis high speed machining and to highlight the kinematic constraints of the reference trajectory.

Firstly, the digital chain of machining is briefly introduced. The specifics in the generation of the programmed tool path is then discussed. After that, the contour error formulations for either point or flank milling are defined.

It is followed by the discussions of the feedrate interpolation techniques leading to presentation of the chosen technique to generate the tool trajectory and the joint motions.

The second chapter, "Adaptive contouring control in multi-axis high speed machining", consists in developing the proposed approach called "Off-line Gain Adjustment" (OGA).

This approach allows to pre-compensate the contour error defined in chapter two, that occurs in the machining process. Firstly, the introduction of the classical axis control structure is presented. The nonlinear feature of the axis model is also provided. Secondly, the state of the art related to the advanced techniques of contouring control is reviewed.

Then, the core section of this chapter is to formulate and solve the optimization problem under constraints related to OGA.

The intent of the third chapter, "Simulations and discussions", is to highlight the improvement of contouring accuracy thanks to the OGA approach proposed in chapter three. Impacts of the resulting optimal variable gain profile, constraint verification and effects of OGA's parameters configuration are shown and discussed through simulation results to validate the performance of the OGA approach.

The fourth chapter, "OGA extensions", is to explore the applications of OGA for diversified purposes. In the first section, the global gain of the position controller is introduced and then tuned by using the OGA approach. Next, the application of OGA in the case of using a Generalized Predictive Control structure is discussed. The corresponding simulation results prove the valuable possibilities of OGA in pre-compensating the contour error within different kinds of position controllers in multi-axis high speed machining.

Finally, "Conclusions and perspectives" aims at evaluating not only the scientific contributions but also the perspectives of this PhD thesis.

Chapter 1

Tool path geometry and kinematics in multi-axis high speed machining

The main goal of this chapter is to highlight tool path geometry and kinematic execution in multi-axis high speed machining, in order to support the formulation of axis control problem in the next chapter.

Firstly, the digital chain of machining is introduced. Multi-axis tool path computations, referring to tool geometry, tool positioning, geometrical errors, tool path description and milling context, are represented afterwards.

In order to propose a method to pre-compensate contour error in off-line process, this error must be estimated based on a machine simulator. Thus, the next section is dedicated to contour error estimation approach.

Lastly, relating to kinematic execution of trajectory, feedrate interpolation is carried out to evaluate its effects on the axis control and machining process.

Introduction of machining process

The digital chain of machining is shown in #2 Trajectory generation: From the nominal surface, CAM computes a set of reference contact points between the tool and the workpiece. Based on these points, a programmed tool path is determined. Most of time, such a tool path is represented by a set of discrete tool location points and tool orientation vectors. These data constitute the CL data, also called CL file.

#3 Post processor [Jung 2002, She 2007]: This is dedicated to transform the programmed tool path into the adapted language for numerical control (NC) in CNC, so-called "NC program".

#4 CNC [Schmitz 2001, Omirou 2005[START_REF] Mori | Remote Monitoring and Maintenance System for CNC Machine Tools[END_REF][START_REF] Yuen | [END_REF], Beudaert 2013, Besset 2017]:

This stage consists of two tasks: feedrate interpolation and axis control. The former allows to generate a trajectory having a feedrate profile, in such a way that the corresponding axis motions, which are obtained by solving the Inverse Kinematic Transformation (IKT), respect the machine constraints, e.g. the axes velocity, acceleration, and jerk limitations, to limit the machine vibration, etc. The later controls the axes motions through closed control loops, axes drive and motor. The combined motions of all axes allow to obtain the tool displacements.

#5 Machining [START_REF] Neugebauer | [END_REF]]: This task is performed in machining center. It is in fact the relative movements between the tool and the workpiece, controlled by the axis control in CNC, to sculpture the machined part.

Multi-axis tool path computation

The aim of this section is to analyze the way that tool path is programmed in the context of multi-axis machining. The main geometrical characteristics of the above tools are given as follows:

Tool geometry

• C E : Center point of the tool,

• C L : Cutter location described in the NC program, which is generally the extremity point of the tool,

• u: Unit vector of the tool axis orientation,

• r: Radius of the hemispheric tool and small radius of the toric tool,

• R: Radius of the cylindrical tool and large radius of the toric tool.

Tool positioning

The programmed tool path is characterized by a couple of values (C L , u), that represents the tool positioning.

To calculate the above components of the programmed tool path, CAM firstly computes the tool contact points, C C , onto the nominal surface. Using different tools for one specific milling technique, for example point milling illustrated in Fig. 1.4, leads to the fact that the tool location points are formulated differently from (1.1) to (1.4) [Duc 1998]. • for hemispheric tool:

a) c) b) C L C C n u C E r t v f R R v f f n n u u r C L C C C C C L C E
C E = C C + r.n C L = C C + r. (n -u) (1.1)
• for toric tool:

C E = C C + r.n + R.v C L = C C + r. (n -u) + R.v (1.2) with v = u ∧ n u ∧ n ∧ u (1.3)
• for cylindrical tool:

C L = C C + R.v (1.4) with v defined in (1.3).
in which f, n and t = f ∧ n, are the unit feed, surface normal and transverse vectors, respectively.

The tool axis direction u is defined by tilt angle θ t and yaw angle θ n in CAM, as seen in Fig.

1.5. In fact, u is obtained by two consecutive rotations: firstly a rotation of n by an angle of θ t around t and then the obtained vector is rotated by an angle of θ n around n [Lauwers 2003].

Tool path description

In multi-axis machining, there are two main formats for the description of tool path: either by linear segments (linear interpolation -G1) or by curves (polynomial interpolation -Bspline) [START_REF] Langeron | [END_REF]]. The former, which is more commonly used in the industrial community than the latter, is presented below. In principle, when the tool moves in the piecewise linear way, the tool contact point C C follows also the piecewise linear way, which approximates the nominal surface within an allowable chordal error. However, this is only true for 3-axis machining. As discussed in [Bohez 2002], in 5-axis machining, even if the linear interpolation is used in CAM, the real C C path is not formed by linear segments, but a curve, as seen in Fig. 1.11. This fact is due to the effect of tool orientation.

Thus, in 5-axis machining, the calculation of the chordal error is more complex than in 3-axis machining, because it should take the deviation generated by the tool orientation into account.

Milling contexts

Another factor that affects the configuration of the programmed tool path is the milling context, referring to point milling, flank milling, and combined milling, as illustrated in Fig. 1.12.

In point milling, the tool is in contact with the workpiece through a contact point C C . In flank milling, the tool intersects the workpiece through a contact segment on the tool trunk. The last milling context combines the above two kinds of intersection between the tool and the workpiece. Fig. 1.12 also highlights differences between the 3-axis and 5-axis configurations. The former refers to a machine having three linear axes, in which the tool orientation is fixed in vertical or horizontal direction depending on machine structure. The latter differs from its counterpart due to the addition of two rotary axes. Thus, the tool orientation can be changed relatively to the workpiece during machining process.

f C C C L u f u C L C L C C f n n u C L f C C f u u C L C C f n a) b) c) d) e) f)
Figure 1.12: Multi-axis machining contexts

Conclusion

According to the above discussions related to the multi-axis tool path computation, there are important remarks as follows:

• Different tool geometries lead to different formulas of the tool positioning.

• It is important to note that the tool does not remove all the material along its programmed tool path, because of the inherent existence of the chordal error and scallop height.

• The 5-axis machining center makes the milling more flexible than the 3-axis one, due to a variable tool orientation. However, it also places the challenges in computing and respecting the maximum chordal error and scallop height.

• The different milling contexts have their own characteristics in sculpturing the workpiece, especially the different contacts between the tool and the workpiece, that should be considered at the axis control stage.

3 Contour error estimation in multi-axis machining

Following, contour error formulation will be developed for 5-axis machining. This can obviously be applied for 3-axis machining as well. For simplicity, the superscripts m and w are denoted for the articular (machine axis) and workpiece spaces, respectively. The trajectory defined by the NC program is considered as the reference, while the trajectory generated by the effective tool movement during machining is considered as the actual one. The superscripts r and a are therefore used for the reference and actual parameters.

3.1 Contour error estimation in 5-axis point milling As it can be seen in Fig. 1.13.a, CAM generates a programmed tool path, characterized by a set of (C r L , u r ), in which the reference C L path can be represented by a series of reference tool location points C r L and tool axis orientation u r . If the tool follows exactly such a programmed tool path, it will sculpture the workpiece to obtain a desired machining profile, that approximates the nominal surface with the allowable geometrical errors, as discussed in § 2.3. This desired machining profile can be seen as the reference C C path, represented by a set of reference tool contact points C r C .

The above programmed tool path is interpreted into NC code, that is then fed into the CNC.

After that, tool axes are controlled to follow the axis position setpoints, that are generated from the programmed tool path through feedrate interpolation and the IKT, without knowledge of the reference tool contact points, but only the couple (C r L , u r ).

In fact, due to the limited bandwidth of the axis drives and the nonlinear characteristics, such as friction or cutting force, the actual tool positioning (C a L , u a ) lags behind the reference one (C r L , u r ), as illustrated in Fig. 1.13.b. Such lag distance is defined as the tool tracking error, e t , in (1.5).

e t = C r L -C a L (1.5)
As it can be seen, there is a deviation between the actual tool path and the programmed one.

Sets of the actual tool location points C a L and of the actual tool contact points C a C represent the actual C L and C C paths, respectively.

In 3-axis machining, many works in the literature [Koren 1983, Erkorkmaz 2001b[START_REF] Yao | [END_REF] defined contour error in point milling as the orthogonal distance (distance ⊥ ) between the actual tool location C a L and the reference C L path, denoted as ε L and expressed in (1.6).

ε L = distance ⊥ (C a L , Reference C L path) (1.6)
The purpose of such a contour error definition is to approximate the effective cutting error, caused by the actual tool contact point C a C . This is because the information of the tool contact points is not available in the axis control and machining stages. This error can be represented by the orthogonal distance between the actual tool contact point C a C and the reference C C path, denoted as ε C , and expressed in (1.7).

ε C = distance ⊥ (C a C , Reference C C path) (1.7)
In 5-axis machining, to improve the approximation accuracy of the effective cutting error ε C , [Altintas 2010, Yang 2015a] considered a complementary error related to the deviation of tool orientation. As it can be seen in Fig. 1.14, even though the tool location point can be tracked with high accuracy, the tool orientation deviation can still cause the cutting error ε C . The error related to the deviation of the tool orientation is denoted as ε u and defined in (1.8).

ε u = angle (u r , u a ) (1.8) C r L C r C u r C a C C a L u a ε u f r f a ε C Over cut

Literature review

The strategies for contour error estimation have some specificities as follows.

Firstly, they depend on the definition of the contour error. Whether contour error is defined by ε L or (ε L & ε u ) or ε C . As discussed in § 3.1.1, in 3-axis machining, contour error is usually defined by ε L [Ghaffari 2016, Yang 2015b, Altintas 2012[START_REF] Yao | [END_REF], Chiu 2001, Erkorkmaz 2006, Koren 1980], while (ε L & ε u ) are usually employed to define contour error in 5-axis machining [START_REF] Li | [END_REF], Yang 2015a, Khoshdarregi 2014, Zhang 2013, Altintas 2010]. To the best of the authors' knowledge, there is no report in which contour error is directly defined by ε C . This can be understood because the tool path information (C L , u) can be easily accessed, while the estimation of tool contact point C C requires further information from the CAM system.

Another factor affecting the estimation approaches of contour error is the purpose of use in axis control strategy: whether it is used for contour error compensation technique during machining process or for contour error pre-compensation technique in the off-line execution before the machining starts. Each control strategy has its own characteristics and requirements, thus the contour error estimation approach has to be consistent with the control strategy.

For contour error on-line compensation technique, contour error has to be estimated as fast as possible to respond to the real-time effect in on-line execution. Therefore, the estimation formula and assumption of contour error must be quite simple and/or the processor has to be powerful enough to overcome the computing burden. Within this context, some related works are proposed in the literature for both 3-axis and 5-axis machining, as follows.

In 3-axis machining, [Koren 1980] proved that contour error is proportional to the tool tracking error in machining linear paths. [Chiu 2001] proposed a task coordinate frame attached locally to the desired contour. By this way, some projections of the actual position tracking error in this moving task coordinate frame can be used to approximate contour error for feedback controller design. More accurate global task coordinate frame has been developed by [START_REF] Yao | [END_REF]]. It calculates contour error exactly to its first order approximation. [Ghaffari 2016] used a Newtonbased update law to estimate contour error.

In 5-axis machining, the influence of the actual tool orientation on contour error must be taken into account. [Altintas 2010] proposed two analytical models, for both tool tip contour error and tool orientation contour error, based on differential path geometry and the kinematics of the machine. [Yang 2015a] proposed that the contouring error components contributed by all axes drives are estimated through interpolated position commands and the generalized Jacobian function. [START_REF] Li | [END_REF]] used the second-order Taylor series expansion and an analytical root-seeking formula to obtain an analytical method rather than an iterative one to estimate contour error and to overcome the computing burden in the real-time process.

For contour error pre-compensation, the computation burden is not a constraint anymore. Thus, the algorithm for contour error estimation has enough time to perform. However, the precision of the estimated contour error much depends on the accuracy of the dynamic model. Within this context, the authors [Yang 2015b, Khoshdarregi 2014, Zhang 2013, Altintas 2012, Prévost 2011a, Erkorkmaz 2006] estimated contour error using the availability of the programmed tool path and the simulated one, in which the later is calculated by solving the Forward Kinematic Transformation (FKT) with the simulated axis position responses obtained from the dynamic model of the machine axes. The general idea of these works is quite simple. First, the closest segment on the programmed tool path to the simulated tool positioning is detected. Then, the contour error is formulated by ε L for 3-axis machining or (ε L & ε u ) for 5-axis machining.

In this work, the intended control strategy belongs to the class of the contour error precompensation technique in off-line execution. Therefore, the intended estimation approach for contour error should take the advantages of the existing approaches within this class. In particular, the calculation technique in [Erkorkmaz 2006] will be considered to formulate the intended contour error estimation approach.

According to the above survey, it can be noted that none of the existing works refer to the estimation of the effective tool contact point, in order to define the contour error by ε C . This fact gives a motivation for this work, that is to propose some assumptions to estimate contour error from the tool contact point C C in multi-axis machining within the context of off-line execution.

In summary, the above related studies and the intention of this work are classified by Table 1.1.

Table 1.1: Some related works for contour error estimation and the intention of this work 3-axis machining 5-axis machining The above issue gives a motivation for the following proposed contour error approximation approach. As a consequence, it is proposed to improve the contour error by controlling directly the C C point and the ε C error.

ε L ε C ε L & ε u ε C On-line [Ghaffari 2016] [ Yao
Contour error formulation. The mathematic formula of the tool contact point derives from the formulation of the tool positioning in § 2.2.

Take the case of toric tool for an example, according to (1.2) and (1.3), the tool contact point

C C is defined by (1.9), C C = C L -r. (n -u) -R. u ∧ n u ∧ n ∧ u (1.9)
in which the tool path information (C L , u) can be reformulated by solving the FKT based on the position setpoints and the position responses, as found in [Lavernhe 2006]. Clearly, the tool contact point C C will be determined if the unit surface normal vector n is known. Determination of n is the main difficulty to estimate the tool contact point.

As it can be seen in Fig. 1.5, as u is known, n can be reformulated if other parameters f, θ n and θ t are known or estimated. To obtain the latter values, it is proposed to make some assumptions as follows:

• For both the programmed and actual tool paths, the unit feed vector f at C C , denoted as f C C , is assumed to be equal to the unit tangential vector as a function of tool displacement at C L , denoted as f C L , if the two following conditions are satisfied:

-The tool diameter is much smaller than the trajectory curvature,

-The tool inclination variation and its angular velocity, ∆θ t and θt respectively, are very small, e.g. several micro-degrees.

This assumption is illustrated in Fig. 1.16.

C C C L f C L f C C _ u f C L f C C Figure 1.16: Estimation of the feed vector When C L is known, f C C and f C L are defined by the first derivative of C L coordination as a function of the tool displacement s, expressed in (1.10), f C C ≅ f C L = d ([C Lx , C Ly , C Lz ]) ds (1.10)
where the tool displacement at the k instant, denoted as s k , defined in (1.11).

s k = k i=1 C r,i L -C r,i-1 L (1.11)
In the following, only the symbol f is used for the feed vector, that is defined in (1.10).

• For θ n and θ t , their reference values are known in CAM configuration. Assume that the responses of θ n and θ t can be estimated by specified values Once all of the reference and actual values of u, f, θ n and θ t are determined, vector n can be calculated as follows.

According to Fig. 1.5, given the unit vectors

u (u x , u y , u z ), f (f x , f y , f z ) and n (n x , n y , n z ). A unit vector m (m x , m y , m z
) is obtained by a rotation of f with an angle of θ n around n, having the same direction with the projection of u onto the plane (t, f). Its coordinates are defined by (1.12).

[m x , m y ,

m z ] T = R n (θ n ) . [f x , f y , f z ] T (1.12) where R n (θ n ), a rotation matrix around n [Kovács 2012], is given in (1.13), R n (θ n ) =     c + n 2 x (1 -c) n x n y (1 -c) -n z s n x n z (1 -c) + n y s n y n x (1 -c) + n z s c + n 2 y (1 -c) n y n z (1 -c) -n x s n z n x (1 -c) -n y s n z n y (1 -c) + n x s c + n 2 z (1 -c)     (1.13) with c = cos(θ n ) and s = sin(θ n ).
Due to the fact that n is perpendicular to f, and the definition of θ t and θ n , the relations in

(1.14) are obtained.          n ⊥ f ∠ (n, u) = θ t ∠ (m, u) = 90 -θ t (1.14) From (1.14), leading to (1.15).          n.f = 0 n.u = n . u .cos (θ t ) m.u = m . u .cos (90 -θ t ) (1.15)
As u, f, n and m are unit vectors, from (1.15), (1.16) is deduced.

         n x .f x + n y .f y + n z .f z = 0 n x .u x + n y .u y + n z .u z = cos (θ t ) m x .u x + m y .u y + m z .u z = sin (θ t ) (1.16)
By solving (1.16), the coordinate components of n can be obtained.

To verify the above assumption, a test case is proposed with the yaw angle θ n sets to zero and the tilt angle θ t sets to 5 degrees. The previous equation (1.16) can be simplified by (1.17).

n = f ∧ u f ∧ u ∧ f (1.17)
Once n is estimated by (1.17), C C is estimated by (1.9).

Calculation technique. To fulfill the proposed contour error (ε C ) estimation, the orthogonal distance from the actual tool contact point C a C to the reference C C path must be found, as illustrated in Fig. 1.13.b.

The calculation technique of the contour error in off-line execution is based on [Erkorkmaz 2006],

whose main ideas are illustrated in Fig. 1.17 and summarized below.

Case 1

Case 2

Case 3

C C a C C a C C a C C i-1 C C i C C i+1 C C i-1 C C i C C i+1 C C i-1 C C i C C i+1 ε C ε C ε C Reference C C path I II III Normal surface of path vector C C i-1 C C i Normal surface of path vector C C i C C i+1 Bi-normal surface of angle C C i-1 C C i C C i+1 Figure 1
.17: Three cases for contour error estimation

• The actual tool contact point C a C can belong to one of the three regions, I, II and III, which are bounded by the normal surfaces of the path segment vectors and the bi-normal surface of the angle ∠C and given in (1.18). and given in (1.19). and given in (1.20).

i-1 C C i C C i+1 C . • Case 1: Contour error, ε C , is estimated using the previous closest C C segment, C i-1 C C i C ,
ε C = distance ⊥ C a C , C i-1 C C i C = ----→ C i-1 C C i C × ----→ C i-1 C C a C ----→ C i-1 C C i C (1.18) • Case 2: Contour error, ε C , is estimated using the next closest C C segment, C i C C i+1 C ,
ε C = distance ⊥ C a C , C i C C i+1 C = -----→ C i C C i+1 C × ----→ C i C C a C -----→ C i C C i+1 C (1.19) • Case 3: Contour error, ε C , is estimated using the closest reference tool contact point C i C ,
ε C = ----→ C a C C i C (1.20)
From (1.18) to (1.20), the absolute value of contour error, ε C , is obtained. To evaluate whether contour error represents an under cut or an over cut, the sign of contour error is a good indicator [Prévost 2011a]. It is proposed that the positive and negative signs of contour error represent for an under cut and an over cut, illustrated in Fig. 1.18.a and b, respectively. The sign of contour error is mathematically defined in (1.21).

C a C C a L u a f a ε C <0 Over cut C a C C a L u a f a ε C >0 Under cut Reference C C path Desired part C i C f i n i a) b) Actual C C path C i C f i n i
sign (ε C ) = sign n i . ----→ C i C C a C (1.21)
Validation. It is clear that the accuracy of the estimated contour error, from (1.18) to (1.20), belongs to the accuracy of the estimated tool contact point in (1.9). The latter belongs to the accuracy of n in (1.17).

To validate n, it is proposed to perform the two following steps:

• Firstly, from the estimated f and estimated n, calculate the estimated tilt angle θ t ,

• Then, compare the estimated θ t with the programmed one.

Note that in the case of study, the programmed tool path generated by CAM has the specificity of tilt and yaw angles as θ t = 5 deg and θ n = 0 deg, respectively.

Contour error estimation in 5-axis flank milling

The general 5-axis flank milling is illustrated in Fig. 1.21. As it can be seen, there exist the over cut and/or under cut due to the deviations of the higher and lower parts of the tool trunk over the workpiece.

C L C H = Hu (1.22)
where H is the distance between C L and C H .

Problem formulation

In 3-axis flank milling, due to the fixed tool orientation, only a limited number of surfaces can be machined, e.g. in flank milling with a 3-axis machine having a vertical tool axis, the walls with fixed orientations should be parallel to the machine tool axis. The definition of contour error in flank milling consists in approximating the above effective cutting errors between the nominal surface and the machined surface.

The literature review related to this topic is briefly presented below.

Literature review

Contour error in flank milling is usually determined by calculating the surface envelope of the tool movement, as illustrated in Fig. 1.7.

To determine the surface envelope, each author proposed their own approach, that generally belongs to one of the two following groups. The first one is the analytical approach, that refers to calculating the implicit equation of the envelope, such as using "Jacobian Rank

Deficiency" [START_REF][END_REF], Yang 2005], or "Sweep Differential Equation" [Blackmore 1994],

etc. The second one is the kinematic approach, that consists in estimating the surface envelope, such as using "Tangency fonction" [Chiou 2004], "Grazing Curve" [Bedi 2003], etc. While the analytical approach is generally complex, time-consuming and difficult to be implemented, the kinematic approach is quite simple to obtain quickly the approximated surface envelope.

That is the reason why the latter will be preferred in this thesis.

Because proposing a new estimation approach for the contour error of flank milling is not the objective of this thesis, the calculation technique of this error in the next paragraph is based on a kinematic approach proposed in the literature, that is [Pechard 2011].

Calculation technique

Similar to the calculation technique of contour error in point milling in § 3.1.3, it is possible to

find out the reference C L path segment C i L C i+1 L , that is closest to the actual tool positioning (C a L , C a H ), illustrated in Fig. 1.24.
From C i L , one can define a specified local plane P i , defined by (1.23),

P t = C i L , u i , f i (1.23)
where the unit feed vector f i is similarly defined in (1.10).

The estimated contour errors generated by the lower and higher tool location points, C L and C H , denoted by ε L and ε H , are provided by (1.24) and (1.25), respectively,

ε L = distance ⊥ C a L , P i = ---→ C a L N (1.24) ε H = distance ⊥ C a H , P i = ---→ C a H M (1.25) P i f i u i t i Over cut Zoom in Under cut M N u a f a t a C H a C L a C L i+1 C L i C H i C H i+1 u i u i+1 u a C H a C L a f a t a C L i ε H ε L Figure 1.24: Contour error definition in flank milling
where N and M are the projections of C a L and C a H on the plane P i , respectively. According to (1.24) and (1.25), contour error of 5-axis flank milling, denoted as ε F , is finally estimated by (1.26).

ε F = max (ε L , ε H ) (1.26)

Conclusion

To conclude, the estimated contour errors for both point milling and flank milling in multi-axis machining center have been formulated in this section.

For point milling, the primary idea of the proposed contour error estimation approach is based on the estimation of the tool contact point.

Meanwhile, in flank milling the estimation of the contour error is still based on the evaluation of the tool location points.

Feedrate interpolation

Another execution before controlling the machine axes is feedrate interpolation. This stage is important, because it generates the position setpoints, that are the input of the axis control stage.

It can be noted that the position response of each axis, which is the output of axis control stage, is expected to follow the position setpoint as close as possible. Thus, the behaviors of the position setpoint, i.e. its kinematic characteristics, will affect the behaviors of the position response as well. Consequently, the latter affects the tool and machining behaviors. For example, if the reference axis jerk exceeds the jerk limit of machine axis, its response may have a high possibility to exceed this limitation also. This will induce the tool vibration during machining, causing marks on the machined part.

Henceforth, in developing an advanced axis control law for the contour error reduction, it is important to take the constraints of the position setpoints into account, in the way that the resulting position responses should also respect such constraints.

To evaluate the constraints of the position setpoints, it is necessary to intrinsically understand the core idea of the feedrate interpolation.

Overview

The schematic diagram of feedrate interpolation is highlighted in Fig. 1.25, which is explained below.

Feedrate interpolation

NC program Kinematical limits, ... {V max ,A max ,J max } } (C L ,u) IKT (X m ,Y m ,Z m ,A,C) Axis motions /Geometry (X m ,Y m ,Z m ,A,C) Geometrical treatment

Position setpoints (C L ,u)

* Feedrate planning technique * Constraints: Firstly, the NC program contains the programmed tool path (C L , u) and the programmed fee-drate F prog . While the former is geometrically treated by the CNC, such as rounding off the discontinuities in the programmed tool path described by G1 code, the latter is always constant, determined by the cutting conditions between the tool and workpiece material. They are the two inputs for the IKT, to determine the joint movements.

Axis motions /Time q(s) q(t) s prog . (=F prog ) F prog s(t) .
Next, the programmed feedrate and the joints movements are evaluated as a function of geometry, namely denoted ṡprog and q (s), where s represents the tool displacement along the tool path.

The two data constitute the inputs of the feedrate interpolation. The outputs of this stage are the feedrate and the sampled axis position setpoints, evaluated as a function of time, and denoted ṡ (t) and q (t), respectively. The core part of the feedrate interpolation is primarily based on the feedrate planning technique with the respect of the concerned constraints, e.g. the axis kinematic limitations.

It can be said that CAM produces the programmed tool path typically without regards to the constraints of the machine, such as axis kinematic limits, machine vibration, etc. Therefore, the purpose of the feedrate interpolation is to generate a feedrate profile to follow the given trajectory, while complying with all of these constraints. Although the programmed feedrate is always constant, the one produced by the feedrate interpolation stage is often decreased to respect all of the concerned constraints.

Some studies related to this topic are briefly reviewed as follows.

Literature review

To perform the feedrate interpolation stage, there are many proposed techniques in the literature.

Each one solves its main problems and obtains its own objectives.

[ Bobrow 1985] solved the minimum time control problem based on the solution from the feedrate interpolation stage, in which the trajectory generation is considered as a dynamic system with two states, including path displacement and velocity; then the velocity along the tool path is limited by the actuator torque constraints. Moreover, the tool path acceleration is switched between its maximum and minimum limits at the identified path points to generate a bang-bang style trajectory.

Furthermore, there are problems associated to the feedrate planning along a curved path, be-ing subject to bounds on the acceleration to achieve minimum traversal time in machining [START_REF] Shin | [END_REF], Shiller 1990]. [START_REF] Timar | Algorithms for time-optimal control of CNC machines along curved tool paths[END_REF]] has addressed such problems for a 3-axis CNC machine, in which by using polynomial parametric curves, the optimal feedrate can be specified by a piecewise-rational function of the parameter. Consequently, a simple real-time interpolator algorithm, directly driving the machine based on the analytic path description, can be performed to obtain the optimized feedrate, without the need for piecewise-linear/circular G code approximations.

To remove feedrate fluctuations due to parametrization errors and to produce continuous profiles of position, velocity, and acceleration, [Erkorkmaz 2001a] proposed a quintic spline interpolation approach, in which fifth order polynomials are used to re-sample the reference trajectory generated with varying interpolation period at the servo loop closure period. This allows preservation of the original kinematic profiles. In addition, to counteract with the feedrate planning for long tool paths, a Linear Programming combined with a parallel windowing algorithm is proposed in [Erkorkmaz 2017], from which the interpolation can be performed over different portions along the tool paths in parallel, with minimal impact on the optimality of the obtained feedrate profile.

On top of this, it can be said that the axis kinematic constraints, including the velocity, acceleration and jerk limitations of the machine axes, are usually considered in the feedrate interpolation stage. This is due to the fact that such constraints allow preservation of the operational life of the drives and of course of the machine, avoidance of saturations, and the prevention of machine vibrations that can potentially lead to severe marks on the surface or instability of the axis control [START_REF] Sencer | [END_REF], Beudaert 2012]. The velocity and acceleration limitations of the drives in the feedrate interpolation are carried out in [Dong 2005, Renton 2000]. Indeed, jerk is also an important parameter, because it allows to limit the feedrate variations in high speed machining.

[Erkorkmaz 2001a, [START_REF] Barre | [END_REF]] have highlighted the effect of jerk limitation on the mechanical structure.

Chosen algorithm

Developing a new technique of feedrate interpolation is out of the scope of this work. Therefore, it is proposed to use one of the proposed feedrate interpolation technique in the literature to generate the position setpoint for the axis control. In fact, the chosen technique for the feedrate interpolation stage is the VPOp algorithm proposed in [Beudaert 2012], which has been developed at the laboratory LURPA. Main characteristics of VPOp are highlighted as follows:

• To reduce the frequency content of the trajectory and to avoid exciting the natural modes of the structure, the limitation of not only the tangential jerk but also of each axis jerk should be taken into consideration.

• The geometrical problem is separated from the temporal interpolation problem by a decoupled approach. The interpolation in sharp zones is based on a look ahead algorithm, in considering the geometrical analysis of the tool path.

• Its mathematical formulation allows processing both linear and rotary axes in the same manner and the algorithm is suitable for either serial or parallel kinematic machines.

The main equation representing the principle of VPOp algorithm, given in (1.27), expresses the upper limits of the feedrate ṡ subject to the minimum value between the programmed feedrate and the upper feedrate limits derived from the limitations of axis velocity, acceleration and jerk.

ṡ ≤ min F prog , V max |q s | , A max |q ss | , 3 J max |q sss | (1.27)
The equation (1.27) allows to predict the areas where the real feedrate will decrease, while respecting all of the kinematic constraints.

The resulting feedrate obtained by VPOp is illustrated in Fig. 1.26. As it can be seen, the feedrate respects all the constraints in (1.27).

Conclusions

The tool path geometry and the kinematic execution in multi-axis high speed machining have been analyzed in this chapter.

Concerning the geometry treatments, firstly the necessary knowledge in computing the multiaxis tool path have been given. Moreover, the issues of the real tool path derived from the changes of the tool orientation in 5-axis milling contexts have been discussed. The contour error estimation approaches in 5-axis point milling and 5-axis flank milling have been formulated. The former is proposed, based on the estimation of the tool contact point through some assumptions.

The latter is based on one of the proposed approach in the literature.

Relating to the kinematic execution, the important roles of feedrate interpolation for the axis control and for the machining are highlighted. Furthermore, the principle and purpose of this task are discussed. The chosen feedrate interpolation technique generates the position setpoint respecting the axis kinematic limitations of the machine. Therefore, the adaptive control law in the next chapter should be developed in the way that the position responses respect these constraints as well.

Finally, it can be concluded that this chapter has successfully provided the necessary elements (the position setpoint generation, the contour error formulation, the axis kinematic constraints) so that the intended adaptive control law for the contour error reduction can be developed.

Chapter 2

Adaptive contouring control in multi-axis high speed machining

The previous chapter has formulated the contour error estimation approach within the different contexts of multi-axis high speed machining. By using the nonlinear axis model, the contour error can be predicted by means of an off-line process.

This chapter aims at solving an optimization problem to produce the optimal variable control gains in pre-compensating the contour error, before it occurs in the machining process.

The introduction of the classical axis control structure is firstly presented. The nonlinear feature of the axis model is then highlighted. Finally, the core part of this chapter is to formulate and solve the optimization problem under constraints.

Introduction of the classical axis control and the nonlinear axis model

The axis control in CNC is built under the classical cascaded structure [Altintas 2000a]. It consists of the position, velocity and current loops, ranging from the outer to inner loops, respectively, illustrated in Fig. 2.1. The sampling time decreases from the position to current loops. The position loop is equipped with a Proportional (P) controller and a velocity Feed Forward (FFW) action. A Proportional-Integral (PI) controller with a torque FFW action is used in the velocity loop, while the current loop uses only a PI controller. The purpose of these controllers is to reduce the axis tracking error, with an expectation that the contour error is also reduced during the machining process.

In this thesis, due to the aim of contour error pre-compensation in an off-line phase, the nonlinear axis model of the above classical axis control is used to investigate the proposed control law. For the case of study, it is proposed to use the nonlinear axis model of the Mikron UCP 710 machine tool located at LURPA, that has been developed in [Susanu 2005], and enhanced and validated in [Prévost 2011a]. Obviously, the proposed control approach can be used for other axis models of other machining centers. It can be found in [Prévost 2011a], the chosen model has been proved to simulate the real behavior in 5-axis milling within Mikron UCP 710 machining center with high confidence. One validation example is for machining the impeller surface illustrated in Fig. 2.2, whose the simulated and measured axis tracking errors are superimposed almost at all portions of the machining profiles, as seen in Fig. In the selected axis model, only the nonlinear friction characteristic is taken into account. The friction law in each axis drive is studied and proposed in [Prévost 2011a] and illustrated in Fig. The fact is that the presence of the contour error definitely indicates the existence of the tool tracking error, the opposite is not always the case. This means that although the tool tracking error may exist, the contour error can be zero. It is well-known that obtaining small contour error is more important than small tool tracking error. This is because the contour error is usually referred to the conformity between the machined and desired parts. The control efforts for contour error reduction can be called contouring control [Erkorkmaz 2001b]. There are two main techniques in order to deal with such a control, including contour error compensation and pre-compensation techniques.

For the contour error compensation, advanced axis controllers are usually developed. There are two sub-philosophies in this way.

On the one hand, authors proposed advanced controllers for axis tracking. The motivation is that when the axis tracking error is reduced in all of the axes, it may indirectly lead to the contour error reduction.

Altintas et al. [Altintas 2000b] built a sliding mode controller to reduce the tracking error of each axis in the presence of external disturbances such as friction and cutting force, obtaining an indirect reduction of contour error. Erkorkmaz et al. [Erkorkmaz 2001b] used a pole placement controller with disturbance cancellation in the feedback loop, in order to deal with the detrimental effects of friction, cutting forces, and drive parameter variations. A combination of vibration avoidance, sliding mode control, torque ripple and friction compensation techniques to facilitate high bandwidth -high accuracy tracking in ball screw drives is proposed in [Kamalzadeh 2007].

The model predictive control is also proposed for the axis control in the context of high speed machining in [Rodriguez-Ayerbe 2014, Dumur 2008].

The benefit of this method is that the idea is simple, as it is based on the axis tracking error reduction to reduce the contour error. Moreover, the tracking error is easily obtained during the control process thanks to the measuring devices. However, it should be noted that the decrease of tracking error, which is still larger than zero, may be less effective in reducing the contour error or can even increase the contour error in some cases of multi-axis machining [START_REF] Koren | [END_REF], Tang 2013].

On the other hand, to directly reduce the contour error, especially in the case of multi-axis machining, authors proposed the advanced contouring controllers. In this way, all axis tracking error and position responses are collected and used to predict the contour error by one of the estimation methods. The control signal of each axis is then modified by one of the intelligent control laws to compensate the contour error [START_REF] Koren | [END_REF], Chiu 2001].

Following this philosophy, Koren [Koren 1980[START_REF] Koren | [END_REF] The benefit of the advanced contouring control approaches is that they can evaluate and control the real behaviors of the machine. The axis drive dynamics can be excited intelligently to compensate the contour error. In contrast, the drawback is the computing load in real time, especially in 5-axis machining. Moreover, its effectiveness much depends on the precision of the contour error estimation method. Even if the contour error can be estimated accurately, the advanced contouring controllers require modification of the classical control structure of the commercial CNC. That may be inconvenient and costly for machine tool manufacturers.

For the contour error pre-compensation techniques, the contour error is firstly predicted based on the dynamic axis model and one of the contour error estimation approaches. Then calibration efforts are performed to pre-compensate the predicted contour error.

To the best of the authors' knowledge, almost all of the works in the literature are dedicated to modifying the position setpoints in pre-compensating the contour error on the machining operation. Khoshdarregi et al. [Khoshdarregi 2014] proposed generating the shaping position commands in the FFW blocks that pre-compensate the resulting contour error and avoid structural vibrations in 5-axis CNC machine tools. Zhang et al. [Zhang 2013] modified the trajectory commands in pre-compensating the contour error and respecting axis kinematic constraints, which are axis velocity, acceleration and jerk limitations. Yang et al. [Yang 2015b] used the idea of MPC to adjust the position setpoints, preserving the axis velocity and acceleration constraints. Furthermore, the contouring accuracy in corner machining application is also improved by a trajectory planning strategy developed in [Erkorkmaz 2006].

The advantage of the contour error pre-compensation techniques is that the contour error can be calculated off-line through the nonlinear axis model and the computational load can be less critical. The modified position is easily implemented by users in CNC. However, this method requires strictly an accurate knowledge of the axis dynamic model obtained from the identification task.

According to the above discussions, it is proposed in this work to investigate a new strategy within the class of the contour error pre-compensation techniques. The proposed idea is that the control gains are adjusted in the off-line process through the nonlinear axis model, with an objective of the contour error reduction during the machining process. This proposed approach, so-called "Off-line Gain Adjustment" (OGA), is developed in detail below.

The works related to this topic in the literature and the intention of this thesis can be categorized by their different strategies through Table 2.1. Trajectory planning strategy [Erkorkmaz 2006], etc.

3 Off-line gain adjustment (OGA)

According to the pros and cons of the above two main contouring control techniques, it is proposed to take the advantages of both methods to compensate their drawbacks. The idea is to maintain the classical cascaded control structure and the position setpoints in the current commercial CNC, while using variable gains for the position controller, e.g. to adaptively excite the axis drive dynamic, to compensate the contour error during the machining process. The only modification on the commercial CNC is the extended functions, that are the gain update at each sampling time and the required memory to store the gain values for the given trajectory, as illustrated in Fig. 2.5.

Optimal variable gains

On-line process

Off-line process

Multi-axis High Speed Machining center Reference trajectory (setpoint)

Off-line Gain Adjustment

Gain modification influence in the contour error

In the proceeding, the cases of study are related to the 5-axis point milling on the 5-axis Mikron UCP 710 machine, having three linear axes {X, Y, Z} and two rotary axes {A, C}. For the simplicity of mathematic formulation, the following symbols are used:

• a: the axes {X, Y, Z, A, C},

• r and s: the reference and simulated values in view of the axis model simulation, respectively.

The core part of the proposed axis control is the OGA algorithm. Thus, the key question is how can the variable gains of position controllers generated by OGA reduce the contour error?

Answering this question, it is essential to analyze the procedure in Fig. 2.6 as follows. the different zones and between the different types of axis. For example: in the zone z 1 , Y axis is more excited than A axis, while in z 2 , the reverse is the case. These setpoints are then fed into the position closed loops of the axis drives to control the tool.

-3 0 ... V m Y V m A (rpm) (m/min) Δt t Δt t Δt Δt z 1 z 2 (d) (X w ,Y w ,Z w ) r z 1 z 2 F r (i,j,k) r (X w ,Y w ,Z w ) s F s (i,j,k) s (X m ,Y m ,Z m ,A,C) r Interpolation IKT (X m ,Y m ,Z m ,A,C) s
(S3) The axis tracking error depends partly on the reference dynamics (i.e spectral content).

For axes with references having different dynamics, this error would be different for each axis, having an impact in the contour error. To reduce the contour error, one possibility is to modify the axis drive dynamics in order to adapt it in function of the reference of each axis.

If the controllers use fixed control gains, the dynamic of axis drives is kept unchanged during the whole process of axis control. In the qualitative point of view, it can be said that the axis drive dynamic is too rigid to adapt appropriately the suitable positioning compliance. That causes the deviation between the resulting tool path and the reference one, that is responsible for the contour error presence.

Therefore, there is a need to excite more or less the axis drive dynamic during the axis control process to compensate the contour error. It can also be remarked that the axes are not needed to be permanently controlled with high gains, because for slow dynamics, a small gain is enough and the high frequencies of the machine are less excited. The OGA is derived based on such an idea.

Basically, when OGA adjusts the proportional gain of the position controller, K P a , it leads to adjusting of the dynamics of axis drives [Pritschow 1996], meaning that the bandwidth of the position loop B a is modified within the range in (2.1),

B min a ≤ B a ≤ B max a (2.1)
when OGA adjusts the feed forward gain of the position controller, K F a , it modifies the impact of feed forward action on the velocity loop.

According to the above effects of the two control gain adjustments, denoting K a as K P a and K F a , when OGA adjusts K a within a boundary as in (2.2), it leads to the fact that the axis tracking error, E a,m , is also modified within a boundary as in (2.3),

K min a ≤ K a ≤ K max a (2.2) E min a,m ≤ E a,m ≤ E max a,m (2.3)
where E X,m is for example the axis tracking error of X axis, given in (2.4).

E X,m = X r m -X s m (2.4) (S4
) By solving the FKT, a set of tool behavior responses in the workpiece space is determined.

Clearly, once the axis motions are adjusted more or less differently, the amplitude and direction of the resulting tangential feedrate vector in the workpiece space are also modified. This leads to the modification of the tool tracking error vector, e t , which is defined in (1.5). Its norm is also bounded by (2.5).

e min t ≤ e t ≤ e max t (2.5) Assume that between e min t and e max t , S represents a region where the tool can approach due to the effects of the above control gain adjustments. In this region, the desired profile (d) illustrates the tool contact positions corresponding to the zero contour error.

It can be said that the gain adjustments in OGA can finally modify the tangential feedrate vector so that the tool contact point can approach the desired profile (d), in order to remove or reduce the contour error.

OGA integrated control structure

The control structure with specificity of OGA is illustrated in Fig. 2.7.a and explained below. to modify the control gains in position loop, K P a and K F a . The torque feed forward action as seen in Fig. 2.1 is not considered in this work, in order to highlight the effects of the velocity feed forward control action and its modification based on the OGA.

Note that the unit of K P a is 1/s. However, in [Prévost 2011a] K P a is proposed in m/min/mm for the linear axes (X, Y and Z) and in rad/min/mrad for the rotary axes (A and C). Both of them are multiplied a factor of 1000/60 to obtain the original unit 1/s. In the proceeding, K P a in either m/min/mm or rad/min/mrad will be used. While, K F a is a constant without unit.

The inputs of OGA are the position setpoints and the position responses simulated by the nonlinear axis model. In OGA, firstly the FKT problem is solved to obtain the reference and simulated tool behaviors in the workpiece space. Then the tool contact points and the contour error are estimated, based on the proposed estimation method in Chapter 1, Section § 3.1. The constraints are checked. Solving the contour error optimization problem under constraints allows one to obtain the optimal variable gains at the output of OGA. The gains are afterwards used to update the control gains in the position loop of axis drive.

In the following, three main cases of study of OGA are considered.

• Case 1: OGA optimizes only K P a , illustrated in Fig. 2.7.b.

• Case 2: OGA optimizes only K F a , illustrated in Fig. 2.7.c.

• Case 3: OGA optimizes both K P a and K F a , illustrated in Fig. 2.7.d.

Optimization problem of OGA

To pre-compensate the contour error before it happens in the online machining process, the OGA generates the optimal variable gains in the simulated machining framework by using the nonlinear axis model. This has been done based on the idea of Model Predictive Control (MPC), that is to formulate the optimization problem within the prediction horizon.

Horizon definition

At instant k, the future trajectory ranging from k +1 to k +N k within a horizon h k is considered, as illustrated in Fig.

C L r,k k+N k O C L ,u r,1 r,1 Δt h,k Δs h,k F k+1 u r,k t k t k+1 t k+N k t 1 s 1 Figure 2.8: a) Horizon h k ; b) Temporal length of h k ; c) Geometrical length of h k
In time domain, as illustrated in Fig. 2.8.b, the temporal length of h k , ∆t h,k , is related to the sampling period of the position setpoint T e and N k through (2.6).

∆t h,k = N k -1 T e (2.6)
In the geometry domain, as illustrated in Fig. 2.8.c, the geometrical length of h k , ∆s h,k is defined by (2.7),

∆s h,k = s k+N k -s k+1 (2.7)
where s k is defined in (1.11). In addition, ∆s h,k can also be expressed through T e and the tangential feedrate F i , as in (2.8).

∆s h,k ∼ = T e   k+N k -1 i=k+1 F i   (2.8)
If ∆t h,k and ∆s h,k are chosen in advance, N k is obtained by solving (2.6) and (2.8), respectively.

There is a difference between the above two definitions of the horizon length. If all of the horizons over the trajectory are chosen as a function of time, having the same value of ∆t h,k , then each horizon has the same number of instants N k due to the constant value of T e . In contrast, if they are chosen as a function of tool displacement length, having the same value of ∆s h,k , then the number of instants N k in each horizon may be different, depending on the different values of feedrate F i . The two ways of the horizon length definition will be discussed in detail to solve the optimization of OGA in Chapter 3, Section § 2.1.

Optimization problem

Assume that at the k instant, the control gain values, the reference and simulated values of the axes motions and the tool behaviors at the present and past instants are known. The objective is to find out the optimal future values of the control gains from the instant k + 1 to the instant k + N k so that the future contour errors are minimized, over the horizon h k .

A matrix containing all position setpoints over the horizon h k , S r , is defined in (2.9):

S r =       X r,k+1 m Y r,k+1 m Z r,k+1 m A r,k+1 C r,k+1 X r,k+2 m Y r,k+2 m Z r,k+2 m A r,k+2 C r,k+2 • • • • • • • • • • • • • • • X r,k+N k m Y r,k+N k m Z r,k+N k m A r,k+N k C r,k+N k       (2.9)
Assume that within such a horizon, the values of K P a and K F a for all axes are given by the matrix K P in (2.10) and K F in (2.11), respectively.

K P =       K P,k+1 X K P,k+1 Y K P,k+1 Z K P,k+1 A K P,k+1 C K P,k+2 X K P,k+2 Y K P,k+2 Z K P,k+2 A K P,k+2 C • • • • • • • • • • • • • • • K P,k+N k X K P,k+N k Y K P,k+N k Z K P,k+N k A K P,k+N k C      
(2.10)

K F =       K F,k+1 X K F,k+1 Y K F,k+1 Z K F,k+1 A K F,k+1 C K F,k+2 X K F,k+2 Y K F,k+2 Z K F,k+2 A K F,k+2 C • • • • • • • • • • • • • • • K F,k+N k X K F,k+N k Y K F,k+N k Z K F,k+N k A K F,k+N k C       (2.11)
The instant variation amounts of the above gain values of the axis a are denoted as δK P a and δK F a , and defined in (2.12) and (2.13), respectively. Note that K P,k a and K F,k a are known at the k instant.

δK P a =       K P,k+1 a -K P,k a K P,k+2 a -K P,k+1 a • • • K P,k+N k a -K P,k+N k -1 a       (2.12) δK F a =       K F,k+1 a -K F,k a K F,k+2 a -K F,k+1 a • • • K F,k+N k a -K F,k+N k -1 a       (2.13)
The sums of squares of the instant gain variations of K P a and K F a are denoted as ∆K P a and ∆K F a , and defined in (2.14) and (2.15), respectively.

∆K P a = δK P a T δK P a (2.14) ∆K F a = δK F a T δK F a (2.15)
The vectors containing the above gain variations for all of the axes, denoted ∆K P and ∆K F , are defined in (2.16) and (2.17), respectively.

∆K P = ∆K P X , ∆K P Y , ∆K P Z , ∆K P A , ∆K P C (2.16) ∆K F = ∆K F X , ∆K F Y , ∆K F Z , ∆K F A , ∆K F C (2.17)
To simulate the machining from the k + 1 to k + N k instant over the horizon h k , each row of S r in (2.9), K P in (2.10) and K F in (2.11) are simultaneously sent to the machine simulator to generate the predicted position responses, denoted by Ŝs , and given by (2.18).

Ŝs =       Xs,k+1 m Ŷ s,k+1 m Ẑs,k+1 m Âs,k+1 Ĉs,k+1 Xs,k+2 m Ŷ s,k+2 m Ẑs,k+2 m Âs,k+2
Ĉs,k+2 , ûs,k+1 ), as illustrated in Fig. 2.9.a.

• • • • • • • • • • • • • • • Xs,k+N k m Ŷ s,k+N k m Ẑs,k+N k m Âs,k+N k Ĉs,k+N k       (2.
From (1.1) to (1.17) and ( 2.19), a vector containing the predicted tool contact points, Ĉs C , is obtained in (2.20).

Ĉs C = Ĉs,k+1 C , Ĉs,k+2 C , ..., Ĉs,k+N k C (2.20)
As illustrated in Fig. 2.9.b, the contour error is calculated by the orthogonal distance from each element of Ĉs C to the reference C C path. To do this, as discussed in Chapter 1, Section § 3.1.3, a vector containing a set of the reference tool contact points, denoted as C r C , is employed. , due to the inherent tool tracking error. Therefore, the size of C r C must be greater than that of Ĉs C to ensure that the contour error can always be calculated. Consequently, C r C is proposed by (2.21),

X (mm) Y (mm) Z (mm) h k C C r,k+1 C C s,k+1 ^εC ^k+1 O C C r,1 C C s,k+N k ^CC r,k+N k X (mm) Y (mm) Z (mm) Reference C L path h k a) C L , u r,k+1 r,k+1 O ĈL , u s,k+1 s,k+1 C L ,u r,1 r,1 b) ĈL , u s,k+N k s,k+N k ĈL , u r,k+N k r,k+N k C C r,k-N b Simulated C L path Reference C C path Simulated C C path
C r C = C r,k-N b C , C r,k-N b +1 C , ..., C r,k+N k C (2.21) in which 0 ≤ N b ≤ k -1 is chosen in advance so that s r,k-N b < s s,k
, where s k defined in (1.11).

If

N b = k -1, the first element of vector C r C is always C r,1
C , the contour error can always be calculated, however it induces a computation burden. This is because the size of C r C will be increased when k increases. To overcome such an inconvenience, depending on the dynamic response of the axis drive, N b can be chosen as a specific value to reduce the size of C r C and so reduce the computation time. For example, during simulations of OGA with the axis drive in this thesis, N b is chosen as in (2.22).

   N b = k -1 if 1 ≤ k ≤ 100 N b = 100 if k > 100 (
J = ÊC ÊC T + λ P ∆K P T + λ F ∆K F T (2.24)
where λ P , λ F are the weighting vectors for the gain variations of K P a and K F a , given in (2.25) and (2.26) respectively.

λ P = λ P X , λ P Y , λ P Z , λ P A , λ P C (2.25) λ F = λ F X , λ F Y , λ F Z , λ F A , λ F C (2.26)
From (2.24), the optimization problem of OGA is formulated by (2.27),

KP , KF = argmin (J) (2.27)
where KP and KF are the optimal values of K P in (2.10) and K F in (2.11).

In solving (2.27), the meaning of the weighting factors is that smaller λ P and λ F , greater ∆K P a and ∆K F a can be accepted, and vice-versa. The optimization problem in (2.27) is subject to constraints, that are discussed below.

Constraints

In order to assure machine kinematic limitations, stability and feasibility, the optimization of OGA in (2.27) is subject to the three following constraints.

Stability criterion

The first constraint is the stability criterion of the servo drive. It is proposed that the variable values of proportional gain in the position controller, generated by OGA, need to respect the limitations of phase margin and gain margin of the position open loop of the axis drive, as in (2.28) [Younkin 2003]. (2.28)

Axis kinematic limitations

The second constraint refers to the axis kinematic limitations. They are the axis velocity, acceleration and jerk limitations, given in V max (m/min -rpm) 30 30 30 15 20

A max (m/s 2 -rad/s 2 ) 2.5 3 2.1 0.83 × 2π 0.83 × 2π J max (m/s 3 -rad/s 3 ) 5 5 50 5 × 2π 100 × 2π          -V max a ≤ V s,i a ≤ V max a -A max a ≤ Âs,i a ≤ A max a -J max a ≤ Ĵs,i a ≤ J max a with k + 1 ≤ i ≤ k + N k (2.29)
where V s,i a for example are the predicted velocity values of the axis a from k

+ 1 to k + N k instant. V s,k+1 X , Âs,k+1
X and Ĵs,k+1

X of X axis for example are calculated by (2.30).

                                 V s,k+1 X = Xs,k+1 m -X s,k m T e
Âs,k+1

X = V s,k+1 X -V s,k X T e = Xs,k+1 -2X s,k + X s,k-1 T 2 e Ĵs,k+1 X = Âs,k+1 X -A s,k X T e = Xs,k+1 -3X s,k + 3X s,k-1 -X s,k-2
T 3 e (2.30)

Motor current limitations

Relating to the motor current response, it is noted that the sampling time of the position loop T e is different from that of the current loop T ′ e . The k instant of the position response corresponds to the k ′ instant of the current response. The relation between T e , T ′ e , k and k ′ is evaluated through the machining time ∆t m , as represented in (2.31).

∆t m = (k -1) T e = k ′ -1 T ′ e
(2.31)

There are two important values related to the motor current. The first one is the instant value of motor current I k ′ a (A). The second one is its increment ∆I k ′ i (A), which is defined in (2.32).

∆I k ′ a = I k ′ a -I k ′ -1 a (2.32)
The two values are subject to their constraints in (2.33),

     I k ′ i a ≤ I k ′ i ,max a ∆I k ′ i a ≤ ∆I k ′ i ,max a (2.33) in which given k ′ in the prediction horizon h k , ∆t m (k + 1) ≤ ∆t m k ′ i ≤ ∆t m k + N k ; I k ′ i ,max a
is proportional to the admissible maximum of the motor torque C max m (N.m) by a factor of K t (N.m/A), as in (2.34).

I k ′ ,max a = C max m /K t (2.34)
Depending on the type of motor, e.g. the Siemens motor 1FT6084 -8WF71 for Mikron Machine [Seimens 2005], the motor torque limit, C max m , is a function of the motor angular velocity, V m (rpm), through the "torque-velocity" diagram in Fig. 2.10, that is represented in (2.35). Therefore, from (2.34) and (2.35), the motor current limit is a function of K t and V m , as in (2.36).

C max m = f (V m ) (2.
I k ′ ,max a = f (V m ) /K t (2.36)
This means that during the control operation, each motor drive has an angular velocity profile, leading to a corresponding profile of the motor current limit I k ′ ,max a , as illustrated in Fig. 2.11 (black dashed curve).

I(A) t(s)

I a k',max I a k',max - k'
ΔI a I(OGA) 0 In addition, respecting the limit profile of motor current is not sufficient if the response of the motor current is too noisy and contains a large variation, because this fact can potentially damage the motor and further the machine. Thus, it is proposed that the maximum current increment is equal to k I (%) of the nominal value of motor current, I n a , which corresponds to the value of the constant torque constraint, C c m , in Fig. 2.10.

Consequently, I n a and ∆I k ′ ,max a are given in (2.37) and (2.38), respectively.

I n a = C c m /K t (2.37) ∆I k ′ ,max a = k I I n a (2.38)
To summarize, at the k instant the OGA aims at solving the optimization problem given in (2.27), being subject to the constraints formulated by (2.28), (2.29), and (2.33), to generate a set of optimal variable control gains used in the position controller of one of the three cases of study, illustrated in Fig. 2.7.b, c and d, in pre-compensating the contour error predicted by (2.23). The developments of OGA are valid for both the temporal horizon length in (2.6) and the geometrical horizon length in (2.7).

Admissible gain range

For each variable gain value at each instant, OGA always checks its feasibility based on the three above constraints. If one certain gain value does not satisfy the constraints, it will not be chosen. Therefore, the control gain range, in which OGA can be performed, can be theoretically unlimited. However, it is clear that if OGA is performed within an unlimited gain ranges, the computation burden turns to be too serious to finish the OGA optimization.

To address this issue, OGA should be performed within one limited range for each control gain, so-called admissible gain range. It is proposed to define the admissible gain range in the way that the classical axis control can use any fixed gain values in this range without violating any constraints in § 3.4.

For the three cases of study illustrated in Fig. 2.7.b, c and d, the admissible gain ranges of both K P a and K F a are established as follows.

• For the admissible range of K P a : The upper limit of K P a , K P,max a , is chosen respecting the stability criterion of the position loop in (2.28), expressed in (2.39).

K P a ≤ K P,max a s.t stability criterion in (2.28) (2.39)

The lower limit of K P a , K P,min a , is globally proposed as in (2.40).

0 < K P,min a < K P,max a (2.40)

One specific value of K P,min a will be chosen by users.

• For the admissible range of K F a : The lower limit of K F a , K F,min a , should be zero, that is the case having only P controller activated, and is described by (2.41).

K F,min a = 0 (2.41)

In fact, K F a is outside the position closed loop, its maximum value does not affect the stability criterion of this loop. However, the variation of K F a affects the behavior of the reference velocity signal in the velocity loop. A too large value of K F a can violate the axis kinematic limitations and/or the current limitations.

The upper limit of K F a , K F,max a , depends on two following factors. The first one is the reference axis kinematic profiles, including the reference axis velocity, acceleration and jerk profiles. One example for the reference jerk profile is that, as it can be seen in Fig.

2.12, case 1 has the reference jerk profile much below its limitations than case 2. Thus, there are more spaces for the variations of the actual axis jerks in case 1 than in case 2. Consequently, the potential K F,max a in the former can be larger than in the later,

K F,max,1 a > K F,max,2 a .
-J max The second factor is the chosen value of K P a . For example, as it can be seen in Fig. 2.13, the reference velocity signal U r a,m is a sum of two control signals U P a and U F a , derived from K P a and K F a respectively. Obviously, when U r a,m is subject to a limit value, the chosen K P a value will affect the upper limit of K F a , and vice versa. In addition, the choice of K F,max a is subject to the two constraints, including the axis kinematic limitations in (2.29) and the current limitations in (2.33).

According to the above facts related to K F,max a , it is proposed two steps for obtaining its value as follows.

-Firstly, the admissible gain range of K P a is discretized into numerous discrete values. With each discrete value of K P a , increasing the value of K F a from zero till one limit value, over which the axis kinematic limitations and the current limitations are not satisfied anymore. This first step is represented by (2.42),

               K P a = K P,1 a , K F a ≤ K F,1 a K P a = K P,2 a , K F a ≤ K F,2 a . . . K P a = K P,N a , K F a ≤ K F,N a s.t.    axis kinematic limitations in (2.29) current limitations in (2.33) (2.42) with 1 ≤ i ≤ N , K P,min a ≤ K P,i a ≤ K P,max a .
-Secondly, obtain K F,max a by (2.43).

K F,max a = min K F,1 a , K F,2 a , . . . , K F,N a (2.43)
To summarize, the control gains K P in (2.10) and K F in (2.11) are bounded within the admissible gain ranges, represented in (2.44),

   K P,min a ≤ K P,i a ≤ K P,max a K F,min a ≤ K F,i a ≤ K F,max a with k + 1 ≤ i ≤ k + N k (2.44)
in which it can be highlighted that K P,min a is chosen by users, respecting (2.40); K P,max a is determined depending on the chosen stability criterion for the servo drive, as in (2.39); K F,min a is obviously zero, as in (2.41); and K F,max a is determined depending on the reference trajectory and the chosen value of K P a , being subject to the two constraints, as expressed in (2.43).

Solve the optimization problem

To solve the OGA optimization problem, two questions need to be answered. How to optimize the variable gains for the future trajectory over the horizon h k in (2.10) and ( 2.11)? Which technique allows to solve the optimization problem in (2.27)? Following, two solving techniques are proposed to solve the optimization problem of OGA for the whole trajectory.

At the first instant, k = 1, the initial values of control gains are freely chosen respecting the admissible gain ranges in § 3.5 to start the simulated machining, as in (2.45).

   K P,1 a = K P,initial a K F,1 a = K F,initial a (2.45)
This leads to the fact that the first simulated/reference values of the axes motions and tool behaviors are achieved. Then, the prediction horizon starts at k = 1 and the OGA is solved to obtain optimal variable gains for k > 1.

First solution technique: OGA(a)

The first solution technique of OGA is denoted as OGA(a). Its motivation is illustrated in Fig.

and explained as follows.

Over the horizon h k in Fig.

K k+1:k+l 1 a = - 1 2 H k a sin (t s ) + 1 2 H k a + K k a (2.46) (l 2 ) : K k+l 1 +1:k+N k a = K k+l 1 a (2.47)
where

• The horizon length characterized by N k can be a fixed value or adaptive values along the trajectory, depending on user's choice for either the temporal horizon length in (2.6) or the geometrical horizon length in (2.7).

• l 1 is the length of the sinusoidal part. This is an user-defined value, proposed in (2.48).

l 1 = k l N k with 0 < k l ≤ 1 (2.48)
If k l = 1, the constant part defined in (2.47) is removed, there is only the sinusoidal part in • t s is a time array of the sinusoidal part, defined by (2.49).

t s = [π/2 : 1/ (l 1 -1) : 3π/2] (2.49)
• H k a represents the gain amplitude modification of such a gain function. It can have a positive, zero or negative value. When the value of H k a is changed, different behaviors of the gain function are obtained for tuning, as illustrated in Fig. 2.14.b.

In respecting the admissible gain range in (2.44), the constraint of H k a is given in (2.50).

K min a -K k a ≤ H k a ≤ K max a -K k a (2.50)
There are many possibilities to propose such a gain function, e.g. using the parabola function instead of the sinusoidal one. However, only the above gain function is chosen, in order to create smooth gain transition between the two consecutive horizons, by choosing an appropriate value of l 1 in (2.48). Other possibilities will be checked in future works.

To solve the optimization problem of OGA, a gridding approach is used. In this approach, l 1 is chosen constant in advance and kept fixed for all of the horizons, the criteria in (2.24) is evaluated for H k a equal to each element of the vector U k a in (2.51),

U k a = -mδ a -(m -1) δ a . . . 0δ a . . . (n -1) δ a nδ a (1,m+n+1) (2.51)
where

• The element {0δ a = 0} of U k a means that if H k a = 0, all of the variable gains in the horizon h k , defined by (2.46) and (2.47), are equal to the current gain value K k a .

• δ a is an amplitude step of the gain function in (2.46). It is subject to (2.52).

0 < δ a ≤ K max a -K min a (2.52)
• Denote p 0 as the total number of potential gain values for K a . From (2.51), p 0 is defined in (2.53).

p 0 = m + n + 1 with 0 ≤ {n, m} ≤ p 0 -1 (2.53)
• The constraint in (2.50) leads to the condition in (2.54).

   0 ≤ nδ a ≤ K max a -K k a 0 ≤ mδ a ≤ K k a -K min a (2.54)
It is proposed that the maximum gain variation, (p 0 -1) δ a , is always less than or equal to the admissible gain range, represented in (2.55),

(p 0 -1) δ a = T 0 K max a -K min a with 0 < T 0 ≤ 1 (2.55)
in which T 0 , denoted for the variation factor, is an user-defined parameter. It allows to regulate the maximum gain variation. If T 0 = 1, the maximum gain variation is equal to the admissible gain range. If 0 < T 0 < 1, the maximum gain variation is smaller than the admissible gain range.

From (2.55), the gain amplitude step, δ a , is defined in (2.56).

δ a = T 0 (K max a -K min a ) (p 0 -1) (2.56)
If p 0 is a fixed value, the value of T 0 also allows to adjust the amount of gain step.

To determine n and m in (2.51), and respect the constraint in (2.53), n is firstly proposed by (2.57) and illustrated in Fig. 2.17. The meaning of (2.57) is that if

(q) : n = p 0 -1 K max a -K min a -K k a + K max a with K k a ∈ K min a : K max a ( 2
K k a = K min a , then n = p 0 -1. If K k a = K max a , then n = 0.
Therefore, when K k a ranges from K min a to K max a , n ranges from p 0 -1 to 0.

Afterwards, from (2.53), m is defined in (2.58).

m = p 0 -n -1 (2.
> K k a -K min a , then mδ a is set as K k a -K min a .
In the case of T 0 = 1, Fig. • In Fig. 2.18.b, if K k a = K max a -δ a , from (2.57) and (2.58), leading n = 1 and m = 3. Thus,

U k a in (2.51) is [-3δ a , -2δ a , -1δ a , 0δ a , 1δ a ]. • In Fig. 2.18.c, if K k a = K max a
, from (2.57) and (2.58), leading n = 0 and m = 4.

Thus, U k a in (2.51) is [-4δ a , -3δ a , -2δ a , -1δ a , 0δ a ].
Similarly, in the case of T 0 = 0.5, Fig. 2.19 shows three examples of the tuning cases of U k a in (2.51). It can be seen that in this case the tuning gain range, the maximum gain variation and the gain step are half of the ones in the previous case with T 0 = 1. From Fig. 2.19.a to Fig. 2.19.c, it can be noticed a receding effect of the tuning gain range, that does not exist when T 0 = 1.

To summarize, the OGA(a) has three main properties: (1) generate the variable gains using a predefined gain transition function as a template, (2) tune and optimize the gain amplitude of the gain function within the admissible gain range to generate the optimal variable gains and

(3) recede the horizon without overlap.

The gain modifications over the trajectory are defined by parameters N k , k l , T 0 and H k a . The parameters N k , k l and T 0 are fixed in advance and are equal for each axis. Only parameter H k a , one for each axis, is obtained by optimization.

Second solution technique: OGA(b)

The second solution technique of OGA is denoted as OGA(b). It takes the advantages of MPC into account, that is to find the optimal gains for the whole horizon h k at the k instant but only the first optimal gain of the horizon is kept for the instant k + 1, as illustrated in Fig. 2.20.a and b. Then the horizon is receded to the instant k + 1 and the optimization is repeated, as shown in Fig. 2.20.c. The method of receding horizon allows one to apply the OGA optimization for the whole trajectory. During the prediction horizon h k , the gain can take different values for each sampling time. In order to simplify the optimization problem and as receding horizon is used, it is proposed that the gain is constant during the prediction horizon.

Assume that δK a (2.59), representing for either δK P a in (2.12) or δK F a in (2.13), is the gain modifications within the horizon h k .

δK a =       K k+1 a -K k a K k+2 a -K k+1 a • • • K k+N k a -K k+N k -1 a       (2.59)
Solving the OGA optimization problem with this method consists in tuning δK a in (2.59) within the admissible gain ranges in (2.44). Applying a gridding approach to solve this optimization problem, the criteria (2.24) is evaluated for δK a in (2.59) equal to the columns of matrix W a in (2.60), to find out all of the possibilities of K P in (2.10) and K F in (2.11), (2.60) in which as same as OGA(a), it is proposed to introduce the total number of potential gain values p 0 in (2.53), and the variation factor T 0 in (2.55). Then, the gain step δ a , the values of n and m are defined in (2.56), (2.57) and (2.58), respectively.

W k a =       -mδ a -(m -1) δ a . . .
One example is that if δK a in (2.59) is equal to the last column of W k a in (2.60), the variable gains become (2.61):

   K k+1 a = K k a + nδ a K k+2:k+N k a = K k+1 a (2.61)
Then, by checking all of the possibilities of K P in (2.10) and K F in (2.11) by means of the machine simulator, the optimal variable gains satisfying (2.27) are found.

The illustration for the values in the first row of W a in (2.60) is similar to Fig. 2.18 and Fig.

4 Conclusion

To conclude, the main part of this chapter has been dedicated to the formulation of the optimization problem of OGA. The optimization problem of OGA is proposed in (2.27), being subject to the stability, axis kinematic, and current constraints, represented in (2.28), (2.29), and (2.33), respectively.

The optimization problem of OGA is basically built for the future trajectory over a horizon, in order to deal with the challenges in machining the free-form geometry. After solving this optimization problem, the horizon is receded so that the OGA optimization can cover the whole trajectory. Finally, an optimal variable gain profile is obtained, pre-compensating the contour error.

For solving such an optimization problem of OGA, there are many techniques that can be performed, such as genetic algorithm, heuristic method, etc. However, the two proposed solving methods for the OGA optimization problem have been preferred here, based on a gridding technique. The primary reason for this choice is its simplicity. That makes the OGA more practical and easily implementable in the machining community. However, such a technique also causes a problem related to computation burden. This drawback will be discussed in Chapter 3, Section § 2.2. The optimization problem under constraints of OGA and the two solving techniques of OGA have been proposed in the previous chapter. The aim of this chapter is to analyze the efficiency of the above developments, especially in the contouring accuracy improvement and the constraint verification.

To do this, the necessary starting points and parameters configurations are performed. Many cases of study are carried out. Simulation results are shown and discussed to prove the performance of the proposed approach.

Starting points

Before performing and evaluating the OGA, the following tasks are considered:

• Generate the reference trajectory and position setpoints.

• Find out the admissible gain range, based on the analyses in Chapter 2, Section § 3.5.

• Find out the best responses of the classical axis control, consider it as the reference case to compare with the performance of OGA.

• Analyze the sensitivity of K P and K F in the position loop.

Reference trajectory and position setpoints

Admissible gain range

An important starting point of OGA is the admissible gain range. Based on the analysis in Chapter 2, Section § 3.5, the upper limit of K P a in the position controller should firstly be determined by respecting the stability criterion. To do so, the linear axis model of the position open loop should firstly be identified from the nonlinear one, that comes from the complete axis model in Fig. 2.7.a and simply illustrated in Fig. 3 Firstly, the nonlinear velocity model is excited by a step signal. The velocity step response of the nonlinear model depends on the amplitude of the step, principally because of the friction in the machine. The identification has been done using a step amplitude equal to the half of the maximum velocity of the axes, that is because the nonlinear axis model is previously identified through this value as well.

Indeed, the numerical model needs to provide the same response as the nonlinear one. This precision depends on the chosen order of the transfer function. Thus, an appropriate choice for the order of the model should be achieved. As it can be seen, the third and fourth order linear models approximate reasonably well the response of the nonlinear velocity loop of Y axis. In the following, the third order model is finally chosen as it is simpler than the fourth one and has the same response. Its mathematical representation in discrete time is given by (3.1),

V s Y,m (q -1 ) U r Y,m (q -1 ) = 1.16q -1 -0.3656q -2 1 -0.4007q -1 + 0.1411q -2 + 0.05416q -3 (3.1)
with q -1 the backward-shift operator.

Secondly, the discrete time transfer function of the integrator action is given in (3.2).

Y s m (q -1 ) V s Y,m (q -1 ) = 0.006 1 -q -1 (3.2)
As a result, the linear model of the position open loop is determined by combining the linear velocity model in (3.1) with the integrator action (3.2), and given in (3.3).

Y s m (q -1 ) U r Y,m (q -1 ) = 0.00696q -1 -0.002193q -2 1 -1.401q -1 + 0.5418q -2 -0.08693q -3 -0.05416q -4 (3.3)
From the obtained linear model of the position open loop, the upper limit of K P Y can be obtained, by tuning this gain, illustrated in Fig. 3.4, respecting the stability criterion in (2.28).

U Y,m (rad/s) (rad) Y m r s K P Y 1000 60 (rad) Y m r Y m s U Y,m
r (q -1 ) (q -1 ) Consequently, the admissible ranges of K P a and K F a for OGA simulations are obtained and given in Table 3.1. Note that, the K P,min a can be chosen as a very small positive value, satisfying (2.40). However, this increases the computation burden of OGA. To have a reasonable computation time for the simulations of OGA in many cases of study, only K P,min a = 0.9 is checked and shown.

Reference response in the classical axis control

To evaluate the benefit of variable gains obtained by OGA in axis control, the contouring accuracy of the classical axis control using fixed control gains is firstly discussed. Table 3.2 shows three cases of study using fixed control gains, in which the K P values are the same as K P * for the three cases, the K F values are slightly different. (3) K F * 0.9 0.8 0.9 14.4

Note that the third case, K P * combined with K F * are the best fixed gains, which are obtained as follows. Testing all combinations of K P and K F in the admissible range from Table 3.1, with a grid step of ∆K P = 0.1 for K P and ∆K F = 0.01 for K F , as illustrated by Fig. 3.6. Such values of ∆K P and ∆K F are chosen considering the sensitivity of these two gains on the reference velocity signal of the velocity loop, that will be discussed in the next section § 1.4. The smallest mean contour error have been obtained with K P * and K F * . This test takes about three days to finish with a PC Intel(R) Core(TM) i7-4790 CPU 3.60 GHz, RAM 16 Go and the testing algorithm is performed on the nonlinear axis model implemented in Matlab R2014b. The contouring accuracy is quantitatively evaluated through the mean contour error, ε, defined in (3.4),

ε = 1 N N i=1 ε i (3.4)
with N is the number of tool positionings on the trajectory.

As it can be seen from Table 3.2, even small changes in the K F value can produce different mean contour error values. Based on the curvature profile, the trajectory can be divided into nine portions, from S1 to S9, having the equivalent curvature characteristic in each portion. S2 and S8 are concave portions with the increasing curvatures. S3 and S9 are concave portions with the decreasing curvatures.

S5 and S6 are convex portions having the increasing and decreasing curvatures respectively. S1, S4 and S7 have curvatures close to zero and can be seen as linear portions. In each portion, the feedrate can be increased and/or decreased, depending on the strategy in the feedrate interpolation task. As it can be seen, the feedrate in the three linear portions is increased quickly. • The contour error in the high curvature portions, S2, S3, S5, S6, S8 and S9, is much bigger than in the linear portions, S1, S4 and S7.

• In the same curved portions, the increasing curvatures (S2, S5 and S8) and the decreasing curvatures (S3, S6 and S9) produce different contour error profiles.

• In S2, S5, and S6, there are slight differences in feedrate tracking, however, their corre-sponding contour error profiles are greatly different. It means that the resulting contour error profile is very sensitive to the feedrate.

• One specified combination of the fixed gain values may produce a small contour error in some portions, but may cause a large contour error in other ones. For example, case (3) has smaller contour error than case ( 2) in S2 and S5; while case (2) has smaller contour error than case (3) in S3, S6 and S9.

From the above facts, there are three general remarks in using the fixed control gains, as follows:

• Firstly, large contour error usually occurs in the curved portions (S2, S5, S6, S9). Meanwhile, the linear portions usually have small contour error (S1, S4, S7).

• Secondly, in the same curved portions, when the curvature increases or decreases, the contour error profile may behave differently (S2 and S3, S5 and S6, S8 and S9).

• Lastly, the modification of control gains leads to modifying the resulting feedrate tracking.

A small adjustment of the resulting feedrate can produce a large change in the contour error.

According to the above discussions, it can be concluded that the curvature characteristics and the control gain values are two main factors that are responsible for the presence and variations of contour error.

In addition, it would be significant to compare the contouring accuracy between the case (3) using the best fixed gain values K P * a and K F * a , and the case using the nominal gain values, K P,n a and K F,n a , proposed by the machine tool manufacturer. In [Prévost 2011a], it can be found the nominal gain values in the position controller of the Mikron UCP 710 machine, as given in Table 3.3. (m/min/mm or rad/min/mrad) 1.5 3.5 1

K F,n a 1 1 1 Figs 3
.8 shows that the obtained contour error of the case (3) using the best fixed gain is clearly much smaller than the one obtained by the case using the nominal gain values. In contrast, for the tool tracking error, the reverse is the case. This can be explained that the nominal gain values are proposed by the machine tool manufacturer for mainly reducing the tool tracking error, expecting to indirectly reduce the contour error. However, the above comparison once again proves that the reduction of tool tracking error does not much reduce the contour error in machining free-form profile.

According to the above results and discussions, the response of the case (3), using the best fixed gain values K P * a and K F * a , is considered as the reference one, in order to compare with the performance of OGA.

Sensitivity of different control gains

The OGA is in fact applied to modify the values of K P a and K F a in the position controller, yielding the changes in the reference signal of the velocity loop and in the resulting feedrate profile after that, in order to compensate the contour error.

Thus, the two following questions come naturally: which control gain allows to achieve more effectivenesses with OGA? And the effectiveness occurs in which conditions? To answer these questions, it is important to analyze the sensitivity of these control gains on the reference signal of the velocity loop.

The reference velocity signals from P and FFW actions, denoted as U P and U F respectively, are compared when one of the two control gains is kept fixed as its best fixed value, K P * a or K F * a , and the other is tuned with different values. For such a analysis, Fig. 3.9 compares U P and U F for Y axis. On top, Fig. 3.9 indicates that when K P Y is kept fixed, the modification of K F Y leads to clear modifications in both U P and U F . It can be seen that greater K F Y , greater U F , smaller U P ; and vice-versa. The contribution of U P to the reference velocity signal is much less than U F . At bottom, Fig. 3.9 shows that when K F Y is kept fixed, the modification of K P Y does not change U F , but slightly modifies U P . This is obvious because U F only depends on the position setpoint and the value of K F Y . In these three cases of study, greater K P Y , greater U P and vice-versa.

According to the above discussions, it can be said that K F a modification is more sensitive than K P a in adjusting the reference velocity signal. It can be predicted that K F a only needs to be slightly adjusted to achieve the benefit of OGA; while K P a needs to be much modified to have similar effect as the K F a adjustment. These predictions will be checked in § 3.

Configurations

To perform the cases of study of OGA(a) and OGA(b), it is needed to configure their parameters, as follows:

• The sampling period of the position loop of Mikron UCP 710 machine is T e = 6 ms.

• The axis kinematic limitations are taken from Table 2.2.

• The initial gain values (2.45) can be chosen as whatever value in the admissible gain range.

In this work, it is proposed to start with the best fixed gains.

• The nominal motor current value I n a (2.37) is calculated based on the documentation of Siemens motor 1FT6084 -8WF71 [Seimens 2005]. ∆I max a (2.38) is chosen as 2% of I n a (k I = 2%), in order to limit noise and large variations of the motor current.

• Other user-defined parameters, including the horizon length, the number of potential gain values, the length of the sinusoidal part in the pre-defined gain function of OGA(a), and the weighting factors, are discussed in detail below.

Horizon length

Chapter 2, Section § 3.3.1 proposed two ways to define the horizon length, either as a function of time or as a function of the trajectory geometry.

For OGA(a), the horizon length is defined considering the following two factors:

• Firstly, because the horizon in OGA(a) completely recedes without overlap, the optimal variable gains generated within a horizon are totally applied for controlling the trajectory over that horizon.

• Secondly, according to the three remarks in § 1.3, there is the fact that the free-form trajectory can be divided into many different portions, and each portion has its own curvature characteristics, e.g. either convex or concave or linear features. Moreover, one specified combination of the fixed control gains is suitable for controlling only one specified portion. Consequently, the horizon lengths should correspond to the above portions. Like this, the variable gains are optimized to be suitable with the curvature characteristics in those portions.

Due to the above two facts, it is proposed to define the horizon length of OGA(a) as a function of the trajectory geometry, which is illustrated in Fig. 3.10 and realized by Algorithm 1.

h1 h2 h3 h4 h5 h6 h7 h8 h9 C 0 -C 0 i 1 i 2 i 3 i 4 i 5 i 6 i 7 i 8 i 9 i 10 Y Figure 3
.10: Illustration of the definition of geometry horizon for OGA(a)

Algorithm 1 Find the horizon length of OGA(a)

1: Calculate the curvature values of the reference trajectory, as illustrated in Fig. 3.10.

2:

Choose the curvature value C 0 so that [-C 0 : C 0 ] represents linear portions. In this work, C 0 = 0.02 (mm -1 ) is chosen.

3: Determine the total instants over the reference trajectory, denoted as n i .

4: for k = 1 : n i do 5:

Determine the first and last indexes on the trajectory, e.g. i 1 and i 10 respectively.

6:

Determine the indexes, e.g. {i 2 , i 4 , i 5 , i 7 , i 8 }, that stay on the boundaries of the above linear portions.

7:

Determine the critical indexes, e.g. {i 3 , i 6 , i 9 }, that are the negative and positive curvature peaks.

8: end for 9: From the above obtained indexes, determine the horizon lengths as a function of tool displacement length as in (2.7).

For OGA(b), at each instant the objective is to find the optimal gain for the first future instant, while considering the future curvature characteristics of the trajectory.

If the horizon length of OGA(b) is chosen as adaptive as OGA(a), it may cause the wrong prediction effect. For example, in Fig. 3.11, the length of linear portion is much longer than that of the curved one. At the k instant, the tool stays at the linear portion. If OGA(b) takes the length of linear portion to evaluate the future trajectory containing mainly the curved portion, then the chosen gain, K k+1 a , applied at the k + 1 instant is not the optimal one. This induces the wrong effect of OGA(b). This value means that if T e = 0.006 s and the horizon has a constant feedrate of 3 m/min, then the geometrical length of such a horizon is 7.2 mm. The tuning effect of different horizon lengths in OGA(b) will be discussed in § 3.7.

Number of potential gain values

For both OGA(a) and OGA(b), the number of potential gain values is given by p 0 (2.53), that is in relation with the value of gain step in δ a (2.56). For greater p 0 , smaller δ a , more possibilities are evaluated to obtain the optimal variable gains, but it is more time consuming. For example, if all of the five axes, {X, Y, Z, A, C} and both K P a and K F a are involved into the optimization of OGA and if p 0 = 20, the number of tuning cases in each horizon will be 20 10 = 1.024e 13 . This takes too huge computation time to finish the optimization of OGA. As a consequence, the OGA implementation would be intractable.

For the simulation of OGA(a) in § 3, p 0 will be chosen as 10. This value means that if all of the three axes {Y, Z, A} and both of the two control gains are involved into the optimization, there are 10 6 tuning cases in each horizon, that is large enough to verify the effect of OGA(a).

Because the horizon in OGA(a) is receded without overlap and the Algorithm 1 allows to generate 9 horizons over the reference trajectory, as seen in Fig. 3.10, thus there are 9.10 6 simulating iterations over the whole trajectory.

For OGA(b), the horizon is receded only one instant. Assume that the trajectory has N instants, there are N -1 horizons. In fact, the reference trajectory in Fig. 3.1.a has 852 instants, meaning that there are 851 iterations. Because the number of horizons in OGA(b) is much more than in OGA(a), p 0 in OGA(b) must be chosen less than in OGA(a), in order to have a reasonable computation time. Finally, p 0 = 5 will be chosen to simulate OGA(b) in § 3, that corresponds a total number of 851.5 6 ≈ 13.10 6 simulations.

Length of the sinusoidal part in the predefined gain function

For OGA(a) only, the length of the sinusoidal part, l 1 (2.46), allows to adjust the shape of gain profile, as seen in Fig. 2.16. Because the horizon length is adaptive, depending on the curvature characteristics of the free-form trajectory, as discussed in § 2.1, the length of the sinusoidal part should be adaptive also. However, how to choose its adaptive values is an open question. To firstly analyze the responses of OGA(a), the length of the sinusoidal part is chosen as 30% of the horizon length. According to (2.48), this means l 1 = 0.3N k (k l = 0.3). The tuning effect of different values of l 1 will be discussed in § 3.6.

Weighting factor

The weighting factors, λ P and λ F (2.24), allow to manage the gain variation amount. Smaller λ P and λ K , greater gain variations ∆K P a and ∆K F a can be accepted respectively, and vice-versa.

In the OGA, the main objective is to obtain the contour error as small as possible, as long as its constraints are respected. If a large gain variation can satisfy this motivation, then it can be chosen. In the sequel, λ P = λ F = 0 is firstly chosen to simulate both OGA(a) and OGA(b).

However, the effect of different weighting factor values on the gain variations and the resulting contour errors would also be worth discussion in § 3.8.

Summary of the configurations

To summarize, the configurations of both OGA(a) and OGA(b) are given in Table 3.4. Particularly, T 0 (2.55) will be a tuning parameter in the next section § 3. ( + ) For the desired trajectory in Fig. 3.1.a, there are 852 instants, corresponding to the total length of the tool displacement of s end = 101 mm. The nine adaptive horizons generated by the Algorithm 1 have the adaptive N k and ∆s h,k as follows: 77,56,111,116,63,84,114,87,144] instant. 8.3126, 10.7076, 12.9685, 17.1182, 13.5307, 11.2934, 9.5413, 8.8023, 8.7577] mm.

N k (h1:h9)=[
∆s h,k (h1:h9) =[

Discussions

This section consists of analyzing not only the effects of the parameter configurations on the performance of OGA, but also the responses of OGA, including the obtained variable gain profile, the axis kinematic responses, the resulting feedrate and contour error profiles.

The cases of study of OGA(a), OGA(b) and the reference case are denoted and explained in Table 3.5. When only one control gain, K P a or K F a , is modified, the remaining control gain, K F a or K P a , is kept constant as its best fixed value, K F * a or K P * a in Table 3.2, respectively.

Impact of the variation factor

In the following, the effects of the variation factor T 0 (2.55) on the performance of both OGA(a)

and OGA(b) are discussed. In fact, if a very high number of potential gains, p 0 (2.53), can be checked by OGA without suffering serious computation burden, T 0 has less meanings. However, as discussed in § 2.2, due to the inherent drawback of computation burden in the gridding solving techniques of OGA(a) and OGA(b), p 0 must be limited, e.g. as 10 and 5 for OGA(a)

and OGA(b), respectively, to assure a reasonable computation time for the OGA optimization.

Due to this fact, the value of T 0 is needed to regulate the gain step δ a and the maximum gain variations. As a result, it allows to compensate inconveniences due to the limited number of potential gains. It can be argued that adjusting δ a can also have the same effect as T 0 . However, T 0 is preferred because it has more physical meanings than δ a , when it directly shows the ratio between the maximum allowable gain variation and the admissible gain range.

The cases of study related to this parameter in OGA(a) and OGA(b) are illustrated in Fig. With the admissible gain range in Table 3.1 and the number of potential gains in Table 3.4, the corresponding gain steps of the above cases of study of T 0 are calculated in (2.56) and given in Table 3.6 and Table 3.7. It can be seen that ranging from Ra1 to Ra5 and from Rb1 to Rb5, the variation factor T 0 is decreasing, consequently the gain steps, δ P and δ F , are decreasing also. Table 3.6: Different gain steps in OGA(a)

δ P δ F Ra1 Ra2 Ra3 Ra4 Ra5 Ra1 Ra2 Ra3 Ra4 Ra5
Y 0.0889 0.0178 0.00889 0.0044 0.000889 0.12 0.0240 0.012 0.0060 0.0012 Z 0.1111 0.0222 0.01111 0.0056 0.001111 0.2167 0.0433 0.02167 0.0108 0.00433 A 0.0222 0.0044 0.00222 0.0011 0.000222 0.1689 0.0338 0.01689 0.0084 0.00338 Table 3.7: Different gain steps in OGA(b)

δ P δ F Rb1 Rb2 Rb3 Rb4 Rb5 Rb1 Rb2 Rb3 Rb4 Rb5
Y 0.04 0.02 0.004 0.0004 0.00004 0.054 0.027 0.0054 0.00054 0.000054 Z 0.05 0.025 0.005 0.0005 0.00005 0.098 0.0488 0.0098 0.00098 0.000098 A 0.01 0.005 0.001 0.0001 0.00001 0.076 0.038 0.0076 0.00076 0.000076

According to the configurations in § 2, the OGA(a) and OGA(b) are performed with the above different cases of study of T 0 . The obtained mean contour errors, ε (3.4), are given in Table 3.8. There are important remarks based on the results in Table 3.8, as follows:

• Almost all of the cases of study with OGA(a) and OGA(b) produce smaller contour error than the reference case r0, except for OGA(b1) and OGA(b2) in Rb1. Based on these quantitative evaluations, it is hard to explain such low performances. However, when the obtained gain profiles and jerk responses are analyzed in the section § 3.4, these worse results of OGA will be clarified.

• OGA(a3) and OGA(b3) are the best cases of OGA(a) and OGA(b), respectively. It can be explained that the OGA(a3) and OGA(b3) have more possibilities in tuning the gains, due to the combination of all potential gains of both K P and K F .

• In OGA(a), OGA(a1) has smaller contour error than OGA(a2) with Ra1; while in all of the remaining cases the reverse is the case. In OGA(b), OGA(b1) is better than OGA(b2) in the cases of Rb1, Rb2 and Rb3; while in Rb4 and Rb5, OGA(b2) is better than OGA(b1).

The above phenomenons can be explained by the sensitivity of the K P and K F values over the reference velocity signal, that has been discussed in § 1.4, i.e. Fig. 3.9. It can be said that in order to see the effects of the K P adjustment, the gain variations of K P should be large enough to be able to regulate the reference velocity signal and of course the feedrate, in order to further compensate the contour error. In the opposite, the small gain variations of K F are suitable to adjust the axis velocity and the feedrate in a delicate manner to compensate the contour error.

• If the gain variations are too small, as in Ra5 and Rb5, they generate too small effects on the axis velocity and the feedrate. As a consequence, the contour error is not much compensated in these cases.

• For OGA(a), the best case of OGA(a1) is Ra1, while the best cases of OGA(a2) and OGA(a3) are Ra3. For OGA(b), the best cases of OGA(b1), OGA(b2) and OGA(b3) are namely Rb2, Rb4 and Rb3. Based on these facts, it is hard to find out a clear method to choose T 0 , in order to achieve the best performances of OGA(a) and OGA(b).

According to the above results and for the given axis drives, it can only be pointed out two following facts:

In OGA(a), for adjusting K P , T 0 = 1 in Ra1 should be chosen. It means that the maximum gain variation of K P should be equal to the admissible gain range. For adjusting K F , its maximum gain variation should be 10% of its admissible gain range, as in Ra3.

In OGA(b), the maximum gain variations of K P and K F should be chosen as within 5% and 0.1% of their admissible gain range, respectively.

Because OGA(a) and OGA(b) have different ways and characteristics in generating the optimal gain profile, it is hard to compare them within the same condition. It is proposed to compare between their best cases in Table 3.8 (in bold). The improvement and the computation time of these best cases are given in Table 3.9. There are some important remarks as follows.

• In general, the best performance of OGA(b) is better than that of OGA(a). However, the computation time of the former is much longer than the latter.

• For OGA(a), the best cases of OGA(a1), OGA(a2) and OGA(a3) improve the contouring accuracy with a contour error reduction of around 67%, 90% and 95%, respectively, compared with the reference case r0.

With the same computation time, OGA(a2) is enhanced around 71% as compared with OGA(a1). Although having 50% of improvement compared with OGA(a2), OGA(a3) takes much more computation time than OGA(a2), 1.5 days compared with 12.5 minutes. • It can be seen that the best cases of OGA(b1), OGA(b2) and OGA(b3) can respectively reduce around 70%, 95% and 98% the contour error of the reference case r0. The OGA(b2)

takes the same computation time as the OGA(b1), but the former reduces around 83% the contour error of the latter. Meanwhile, the contouring accuracy of OGA(b3) is improved around 66% compared with OGA(b2), but the former takes 7 times longer to finish the computation than the latter.

• The best case of OGA(b), OGA(b3) in Rb3, has improved 67.14% the contouring accuracy of the best case of OGA(a), OGA(a3) in Ra3.

From the above results and discussions, it can be said that to have a reasonable computation time and a good contouring accuracy, it should be chosen OGA(a2) or OGA(b2). However, if in off-line execution the computation time is normally not limited, and the contouring accuracy has the highest priority, it is clear that OGA(a3) and OGA(b3) are preferred.

Impact of variable gains on contour error profile

This section is to illustrate the variable gains and contour error profiles, obtained from the best performances of the cases of study of T 0 in § 3.1, represented in Table 3.9. Theses responses are plotted in function of the trajectory in the workpiece space.

Fig. 3.13 gives a general view for the best variable gains and contour error profiles of OGA(a1) in Ra1, OGA(a2) and OGA(a3) in Ra3.

Zoom in

Fixed gains Variable gains Figure 3.13: Resulting variable gains and contour error profiles of OGA(a)

For OGA(a), the horizon boundaries are also plotted to easily check the variable gain formulation due to the predefined gain function within each horizon. It can be seen that all of the variable gains respect the admissible gain range in Table 3.1. The gain variations of OGA(a3) and OGA(b3) are slightly less than the remaining cases. It can be understood that, when only one control gain is adjusted, it must be much changed in order to compensate the contour error. While, if both K P and K F are adjusted, a small modification of each gain is enough to have the effect of contour error compensation.

These variable gains affect the contour error profiles, as seen in the three plots at the bottom of Fig. 3.13 and Fig. 3.14. From OGA(a1) to OGA(a3) and from OGA(b1) to OGA(b3) the contour error profile is improving in a way that it is approaching the zero line over the whole trajectory.

In OGA(a1) and OGA(b1), there are contour error peaks at the region around Y w = 60 mm.

While, in OGA(a2) and OGA(b2) there are smaller contour error peaks at the region around Y w = 20 mm. These contour error peaks are all smaller than the ones in the reference case. In OGA(a3) and OGA(b3), all of the contour error peaks appeared in the previous cases have been removed. The zoomed views indicate that the contour error of OGA(a3) and OGA(b3) only fluctuates several micrometers and several tenth of micrometers, respectively, around the zero line of contour error.

It can generally be said that the optimal variable gains of either OGA(a2) or OGA(a3) or OGA(b2) or OGA(b3) have mostly overcome the challenges due to the high curvature along the reference trajectory, in order to generate a small contour error over the whole free-from trajectory.

Impact of variable gains on axis kinematic and motor current responses

In this section, OGA(a3) in Ra3 and OGA(b3) in Rb3, presented in § 3.1, are employed to highlight the effect of the obtained variable gains on the axis kinematic and current responses. 3.15 and Fig. 3.16 show that all of the axis kinematic responses respect their limitations in Table 2.2. With a general view, when the control gains are adjusted, there are small variations of the axis position, velocity and acceleration, but a large variation of the axis jerk, especially in OGA(b3).

A comparison of the axis kinematic responses between OGA(a3) and OGA(b3), e.g. Y axis, is illustrated in Fig. 3.17. It can be noticed that the gain profile of OGA(b3) fluctuates more than the one of OGA(a3). As a consequence, the axis kinematic responses of OGA(b3) is much more excited than that of OGA(a3), especially the axis jerk. While the axis position, velocity and acceleration responses are quite the same as the reference ones, the response of axis jerk is more excited.

For Y axis For Z axis

For A axis

Relation between variable gains, jerk response and contour error

This section consists in answering the question: why OGA(a1) in Ra1 (T 0 = 1) is better than in Ra2 (T 0 = 0.2), while OGA(b1) in Rb1 (T 0 = 0.1) is worse than in Rb2 (T 0 = 0.05), in which OGA(b1) in Rb1 is even worse than r0? Fig. 3.19 shows the obtained variable gains, axis jerk responses (e.g. Y axis) and the resulting contour error profiles of OGA(a1) in Ra1 and Ra2. There are two important remarks as follows.

Firstly, large gain variation of K P a in OGA(a1) does not much excite the jerk response. As it can be seen, the gain variations in Ra1 are much larger than in Ra2. This is because T 0 in Ra1 is five times more than in Ra2. This leads to the fact that the jerk variations in Ra1 are more significant than in Ra2. However, the axis jerk variations in Ra1 are still much less than its limitation and they only occur at the sinusoidal gain part in each horizon.

Secondly, due to the sensitivity contribution of K P a in the position controller discussed in § 1.4, it can be stated that it would be good to have a large variation of K P a to adjust the reference velocity signal of the velocity loop and the feedrate after that, in order to compensate the contour error.

Clearly, the above two facts prove that for the contour error compensation in OGA(a1) the large gain variation of K P a in Ra1 is better than the small one in Ra2. Consequently, if a large gain variation for the minimum contour error does not satisfy the constraints, especially the axis jerk limitation, the OGA will remove this optimal choice. Therefore, the obtained gain profile of OGA(b1) in Rb1, containing many constant gain portions and several large gain jumps, is not the optimal one for the contouring following anymore and even worse than the classical case. Meanwhile, the small gain variations in Rb2 can both compensate the contour error and still respect the jerk limitation. These facts explain why the contour error compensation of OGA(b1) in Ra1 is worse than in Ra2 and even worse than r0.

Relation between variable gains, feedrate, tool tracking error and contour error

This section is dedicated to analyze the relation between the variable gains, the feedrate, the tool tracking error and the contour error in the workpiece space. While the contour error is much reduced by OGA, the resulting tool tracking error of OGA is sometimes increased as compared with r0. This confirms the fact that a zero contour error can be obtained even if the tool tracking error still exists.

The zoomed views in Fig. 3.21 and Fig. 3.22 prove that both OGA(a) and OGA(b) make the tool displacement very accurate over the peak of the convex curve zone, which is one of the critical regions in machining the free-form trajectory.

The comparison of the contouring accuracy between OGA(a3) and OGA(b3) is illustrated in detail through Fig. 3.23. Over the trajectory, the contour error of OGA(b3) mostly fluctuates within the limit of 1 µm, while the OGA(a3) has some contour error peaks from 2 to 4 µm. 3.6 Impact of the sinusoidal part of OGA(a)

The effects in tuning different lengths l 1 (2.48) of the sinusoidal part in the predefined gain function (2.46) of OGA(a), are now discussed. Note that l 1 is kept the same over all of the adaptive horizons of OGA(a).

Three cases of study related to the length l 1 of OGA(a2) in Ra3 from § 3.1 are given in Table 3.10. Remember that OGA(a2) in Ra3 of Table 3.8 has used l 1 = 0.3N k . As it can be seen, the above three cases of study have compensated the contour error much more than the reference case r0. L1 reduces about 60% the contour error of r0. L2 and L3 similarly reduce about 90% and 75 % the contour error of r0 and L1, respectively.

In addition, the effects of different lengths of the sinusoidal part on the obtained variable gain profile, the axis jerk responses and the contour error profile are illustrated in Fig. 3.25.

As shown in Fig. 2.16, longer l 1 , smoother gain profile, and vice-versa. Because L1 has the smallest l 1 in the three cases, its gain profile contains sharp changes. As it can be noticed, several sharp gain variations in L1, even with small amplitudes, can also lead to the corresponding large peaks in the axis jerk responses. This implies that other large and sharp changes of the optimal gains in L1 have been eliminated from the OGA optimization due to the constraints of axis jerk.

Therefore, the gain profile of L1 is less flexible and has less optimal variable gains, leading to the fact that its contour error is not much compensated. This means that the last two cases can generate more optimal gains than the first case. As a result, the contour error in L2 and L3 is reduced much more than in L1. Moreover, a smaller contour error peak of L2 around Y w = 20 mm has been removed in L3.

According to the above discussions, it can be concluded that longer length of the sinusoidal part, more flexible profile of the variable gains, less excitations of the axis jerk responses, and more contour error reductions. However, how to globally choose the optimal length of the sinusoidal part for the different adaptive horizons is still an open question. Up to now, it is still needed to do some tests and choose the best fixed length of the sinusoidal part for all of the adaptive horizons of OGA(a).

Impact of the horizon length of OGA(b)

This section refers to the tuning effect of different values of the fixed horizon length on the performance of OGA(b).

Four cases of study related to the horizon length of OGA(b1) in Rb2 from § 3.1 are provided in Table 3.11. Remember that the previous value of the horizon length of OGA(b1) in Rb2 in Table 3.8 is ∆t k = 0.144 s (or N k = 25 instants). As it can be seen, Hb1 has the largest contour error peaks at the critical positions on the trajectory, being the curvature peaks on the convex and concave portions. Even its large gain variations at these critical positions cannot remove the contour error peaks. This implies that the reactions of the variable gains in Hb1 are too slow to deal with the changes in the curvature at the critical positions, which usually cause large contour errors. This is due to the shortest horizon length of Hb1. Therefore, it can be said that if the horizon length is too short, the OGA(b) will lose the prediction effect at the critical positions on the trajectory.

In Hb4, due to having the longest horizon length, its variable gain profile behaves too early as compared with the other cases, e.g. as seen in the K P Y profile at Y w = [30 : 50] or the K P A profile at Y w = [50 : 70]. Besides, it can be noticed that at Y w = [50 : 60], Hb4 has the largest contour error in the four cases of study. This can be understood since the detrimental effect of a too long horizon length of OGA(b) discussed in § 2.1, illustrated in Fig. 3.11, that leads to the fact that the obtained variable gain is not optimal for the closest future tool position.

According to the above discussions, the horizon length should be chosen appropriately not only to have the prediction effect, but also to generate the optimal gains. In this way, the horizon length in Hb2 is finally chosen. However, how to choose explicitly the horizon length for OGA (b) is still an open issue. It is proposed that some tests must be performed to choose the best horizon length for OGA(b).

Impact of the weighting factor

The tuning effects of different values of the weighting factor on the OGA(a) and OGA(b) responses are discussed in this section. The cases of study related to this parameter for OGA(a2) in Ra3 and OGA(b2) in Rb3 from § 3.1 and their quantitative results are given in Table 3.12. It can be seen that when the weighting factor λ F is increasing, the mean contour error of OGA(a2) is increasing; while the mean contour error of OGA(b2) slightly decreases. This can be explained through the fact that OGA with large weighting factor permits only small or maybe too small gain variations. This means that for OGA(a) a small gain amplitude of the pre-defined gain function leads to small gain variations for each horizon, that cannot sufficiently adjust the axis velocity, feedrate, and then cannot much compensate the contour error. In contrast, each optimization cycle of OGA(b) only produces one single optimal gain. Therefore, small gain variation of each single gain due to large weighting factor does not much affect on its performance. Moreover, similar to the cases of study in § 3.1, for good contouring accuracy, the gain variations of OGA(a) are usually larger than that of OGA(b). Henceforth, the choice of weighting factor depends on the chosen solving technique of OGA, either OGA(a) or OGA(b).

In addition, the effects of the different weighting factors on the obtained gain profiles, the axis jerk responses, and the resulting contour error profiles are illustrated by Fig. 3.27 and Fig. 3.28.

For OGA(a2), when the weighting factor value is zero, as in Wa1, the variable gain profiles undergo larger adjustments, but the axis jerk variations are quite similar to the other cases.

Consequently, the contour errors are more reduced. For OGA(b2), when the weighting factor value is zero, as in Wb1, the variable gain profiles also undergo more adjustments, and the axis jerk responses are more varied than the other cases. Consequently, the optimal gains having large variations may be removed by the OGA constraints, especially the limited jerk constraint, leading to the fact that their contour errors are not much reduced.

Figure 3.28: Impact of the weighting factor on the OGA(b2)

The above discussions highlight the fact that the weighting factor is a good parameter to regulate the behaviors of not only the variable gain profiles, but also the axis kinematic responses, i.e.

the axis jerk response, and the resulting contour error profile.

Conclusion

To conclude, the configuration and the performance of OGA are discussed through the simulation results of the two solving techniques, OGA(a) and OGA(b). Each solving technique has its own advantages and disadvantages or challenges, as follows.

For OGA(a)

The advantages include: * The adaptive horizon lengths allow OGA(a) to consider the challenges of the different curvature characteristics over the free-form trajectory. * The adaptive horizon lengths are easily obtained by the Algorithm 1, that can be applied for whatever free-form trajectory. * The principle of receding horizon without overlap makes the computation time of OGA(a) shorter than OGA(b). * The behaviors of the variable gain profile and the axis kinematic responses can be managed by the formulation of the gain function, i.e. by choosing an appropriate length of the sinusoidal part. * The obtained variable gains are shaped into the gain function. They can be computed on-line using the gain function and its optimal parameters. Hence, the optimal variable gains obtained by OGA(a) will need less memory to be stored if OGA(a) is implemented into the real machine.

The challenges include: * There are many possibilities to define such a gain function of OGA(a). Finding an optimal gain function is an open issue. * Besides, how to choose the optimal length of the sinusoidal part within the different adaptive horizons is also an open question.

For OGA(b)

The advantages include: * The solving technique is basically based on tuning/searching. The idea is simple, so it can be easily understood and used by new users.

Their common challenges include: * There are many parameters for choice in the configurations of OGA(a) and OGA(b). How to choose optimal values for some of the user-defined parameters, e.g. the length of the sinusoidal part in OGA(a) or the horizon length in OGA(b), are still open questions. * The optimizations in the cases of OGA(a3) and OGA(b3) can significantly improve the contouring accuracy, but they have a high computation burden.

According to the above discussions, it can be concluded that if the computation burden can be regardless or overcome within the off-line execution phase, tuning both K P a and K F a is the best choice for the contouring accuracy. If not only a high contouring accuracy but also a short computation time are needed, then tuning only K F a is preferred.

In general, the proposed OGA has significantly improved the contouring accuracy of the given free-form profile, while preserving all of the concerned constraints. It could be a practical solution for the axis control in the future CNC.

The last two chapters have presented, analyzed and proved the core idea of this PhD thesis, that is to keep the classical axis control structure being used in the current CNC and to modify only the control gain inside the position loop by the OGA approach to reduce contour error. In this chapter, the applicable possibilities of OGA are extended.

In the first part, one global control gain is introduced into the classical axis control structure and adjusted by OGA. In the second part, the OGA is also examined in the case using advanced controllers for axis control.

Simulation results are shown and discussed to prove the possibilities of these OGA extensions.

Finally, an OGA-GUI interface is created to facilitate the use of OGA.

Introduction of the variable global gain

Generally, the OGA is proposed to adjust the control gains used in the position controller This chapter assumes that machine tool manufacturers can develop the axis position control with either classical or advanced controllers. Although the number of control gains or parameters in the position loop can be numerous depending on which kind of controller is used, there is only one reference signal of the velocity loop. Thus, the adjustments of several or all of the control gains result in the modification of the reference signal of the velocity loop, as illustrated in Fig. The more gains adjusted, the more degrees of freedom of the adjustment. However, this leads to a higher computation burden, especially if the optimization is performed for all of the axes. If only one or several control gains are optimized by OGA, it is important to choose the appropriate gains to adjust. The sensitivity and effectiveness of such gains need to be discussed firstly. This is due to the fact that their contribution to the global control action is not the same. For example, as shown in the previous chapter, the velocity reference produced through K F is much more important than through K P . Thus, adjusting the K F gain is more sensitive and more effective than adjusting K P . However, K P is essential to maintaining feedback loop stability.

According to the above analysis, to reduce the computation burden and to adequately maintain the efficiency of OGA, a new idea is proposed, in which the control gains inside the position controller, either a classical or advanced one, are kept unchanged, then one variable global gain is introduced just after the position controller, ahead of the velocity loop, as represented in Fig. 4.2. This gain is without unit, as its function is just to modify the amplitude of the reference velocity signal. The OGA is then performed to adjust this new gain. Before analyzing the efficiency of this idea related to the contouring accuracy, it is obvious that the computation time in this case of study is much reduced compared to the global case where all the control gains are optimized, and it is equivalent to the case where only one classical control gain, e.g. either K P or K F in the classical axis control, is adjusted.

Another benefit of the global gain adjustment is that it can easily be applied into any type of position controller without regard to the complexity of its control structures. Therefore, CNC technicians or specialists can easily implement this new gain for the application of OGA into their commercial CNC versions.

In the following, such a global gain is namely integrated with the classical axis control structures and with one of the advanced axis controllers.

Admissible range for the global gain

In fact, the principle of OGA formulated in the previous chapter is kept unchanged. The only difference is that the optimization task of OGA is now to optimize the global gain, instead of the P and/or FFW control gains. Therefore, the only requirement is to establish constraints for this new gain.

It can be noticed that the adjustment of K global a is equivalent to the adjustments of both K P a and K F a by the same factor. Therefore, the adjusted K global a values will affect the concerned constraints proposed in Chapter 2, Section § 3.4, including the stability of the axis drive, the axis kinematic and the motor current limitations. As a result, the bounds of K global a must be chosen respecting all of the constraints. They are proposed as follows.

Firstly, without the variable global gain in the position loop, K P a and K F a are chosen as fixed values. In this case of study, it is proposed to choose the two gains as their best fixed values in view of contouring accuracy, K P * a and K F * a , given in Table 3.2, and represented by (4.1).

   K P a = K P * a K F a = K F * a (4.1)
The preserving the axis kinematic limitation in Table 2.2 and the motor current limitations in (2.33). A set of equations in (4.3) must be obtained.

       K global,1 a = K P,max a K P * a K global,2 a = K F,max a K F * a (4.3)
Consequently, K global,max a is in relation with K P,max a and K F,max a though (4.4).

     K global,max a = min K global,1 a , K global,2 a = min K P,max a K P * a , K F,max a K F * a (4.4)
Then, the lower bound of K global a , which is obviously greater than zero and less than K global,max a , is chosen by users.

Finally, the admissible range of K global a is given in (4.5).

K global a ∈ 0 : K global,max a (4.5)

Simulations and discussions

The applications of the two solving techniques of OGA on K global a are denoted as OGA(a4) and OGA(b4).

Configurations of OGA

The configurations of OGA(a4) and OGA(b4) are the same as the ones given in Table 3. The variation factor T 0 in (2.55), which is examined as a tuning parameter, is similar to the cases of study illustrated in Fig. 3.12.

Responses of OGA

The corresponding responses of OGA(a4) and OGA(b4) for each case of study of T 0 are discussed as follows. It should be kept in mind that T 0 in (2.55) is now applied for K global a , while the two classical control gains, K P a and K F a , are kept fixed as K P * a and K F * a , respectively. The obtained quantitative results are compared with the reference case r0 in Table 4.2. It can be seen that the best performances of OGA(a4) and OGA(b4) correspond to the cases of Ra3 and Rb4 where they reduce the contour error of r0 about 77% and 79%, respectively.

Although the contouring accuracy of OGA(b4) in Rb4 is slightly improved as compared with OGA(a4) in Ra3, the computation time of the former is 20 times longer than the latter.

To highlight the differences in adjusting the global gain, K global , and other control gains, either K P or K F , a comparison between Table 3.8 and Table 4.2 is performed.

It can be seen that the best cases of OGA(a4) in Ra3 is better than the best case of OGA(a1) in Ra1 and worse than the best case of OGA(a2) in Ra3. Similarly, the best cases of OGA(b4) in Rb4 is better than the best case of OGA(b1) in Rb2 and worse than the best case of OGA(b2) in Rb4.

The contouring performances of OGA(a4) and OGA(b4) are respectively worse than OGA(a3)

and OGA(b3) in all of the cases of study. This is due to the fact that the former cases have less degrees of freedom in tuning the gains than the latter cases. Because of this, the former cases take much less computation time than the latter cases.

Based on the above comparisons, it can be said that the performances of OGA(a4) and OGA(b4) are more or less the same as OGA(a1) or OGA(a2) and OGA(b1) or OGA(b2), respectively.

In addition, although the adjustment of global gain reduces the contour error less than the adjustment for both classical control gains, the former achieves more reasonable computation time than the latter. Similar to the adjustment of K P a and/or K F a , the efficiency of the K global a adjustment much depends on the parameter configurations, such as the gain variation factor T 0 . The contour error profiles of OGA(a4) and OGA(b4) have a peak of about 30 µm at the critical regions just over Y w = 20µm and Y w = 60µm, respectively. However, in general their contour errors have been much reduced as compared with the one of the reference case.

The contour error histogram shows that the contour errors of OGA(a4) and OGA(b4) have a centered normal distribution, mainly ranging from -10 to 10 µm, while r0 induces the contour errors mainly between 5 to 20 µm.

For OGA(a4), over all of the five cases of study of T 0 , it has much improved the contouring accuracy with at least 40 % improvement as compared with the reference case r0. As it can be seen in Table 4.2, it should be good to have the variation factor T 0 ranging from 0.05(Ra4) to 0.2(Ra2) for OGA(a4).

For OGA(b4), from Rb2 to Rb4, the contour error has been reduced as compared with r0, however the mean of contour error in Rb1 is much larger than in r0. This reduced contouring performance of OGA(b4) in Rb1 is due to its largest variation factor, which is explained in detail as follows. It can be seen that the gain profile in Rb4 is much smoother than in Rb1. In fact, the gain profile of Rb1 consists of straight lines and large gain jumps, which are explained by the effect of the kinematic constraints, i.e. the axis jerk limitation, in solving the optimization of OGA.

As it can be seen, each gain variation of Rb1 induces a large peak of jerk. Therefore, it can be understandable that if one optimal gain value of Rb1 in term of the smallest contour error corresponds to a large gain variation, that violates the axis jerk limitation, this gain is not chosen by OGA. Consequently, the gain profile of Rb1 remains unchanged in many segments. At the end of the optimization procedure along the trajectory, the obtained variable global gain profile

is not an optimal one for the purpose of contour error compensation.

According to the above relation between the gain variation, kinematic constraints and the chosen gain value, the obtained gain profile in Rb1 causes a poor CE profile, as shown in Fig. 4.6. As it can be seen, at the beginning of the trajectory, OGA(b4) in Rb1 works well when it allows the smallest contour errors within the three cases of study. To obtain this, the gain profile of Rb1 started with large gain variations. For example, the initial variations of K global Y are as large as 1.2. Then, because the global gain has been increased too much, its next gain variations may violate the jerk limitations. Therefore, after a good start its gain profile kept unchanged until around Y w = 40 mm. This explains why the obtained gain profile in Rb1 is not optimal, and consequently the contour error in Rb1 is not reduced, but even increased along the trajectory.

Zoom in

Zoom in

Other analyses of OGA(a4) and OGA(b4) can be performed the same as the ones in previous chapter, such as the impacts of the sinusoidal gain part, of the horizon length or of the weighting factor on the OGA responses(Chapter 3, § 3.6, § 3.7, § 3.8). However, for simplicity, only the axis kinematic responses of both OGA(a4) and OGA(b4) are discussed as follows. As it can be seen, all of the axis kinematic responses respect their limitations in Table 2.2.

It can be seen that in general the amplitude variation of the global gain in OGA(b4) is smaller than in OGA(a4). However, the gain profile of the former varies nearly every instant along the trajectory, while that of the latter has many constant gain segments due to effect of the constant part in its predefined gain function.

Conclusion

Indeed, it is hard to conclude that the adjustment of the K global a is more interesting than the previous adjustment of the classical control gains, K P a and/or K F a . It can only be said that it has the advantages of the one single gain adjustment, especially the reasonable computation time. In addition, the above simulation results prove one important thing that the OGA can be applied correctly as a means of modifying only one single global gain, in the context of the classical axis control structure. Therefore, it is inspired to test this idea with the advanced axis control structure. It makes sense because like in this case the applicable possibilities of OGA can be further extended.

OGA integrated advanced axis controller

In this section, the application of OGA for the global gain is examined within the structure of an advanced position controller.

It should be noted that although the development of advanced controller for feed drives is out of the main scope of this PhD thesis, it is also valuable to highlight the generic possibilities of OGA, not only for the classical axis control but also for advanced structures. Moreover, this idea is developed with an assumption that the advanced position controllers can be implemented in the future in a commercial CNC version, which is not obvious.

A case of study for advanced feed drive control is the GPC (Generalized Predictive Control), which has been extensively studied in the literature [START_REF] Susanu | [END_REF], Dumur 2008, Rodriguez-Ayerbe 2014].

Note that in the previous chapters, the sampling period in the framework of machine tools for the position loop and the setpoint interpolation period were set to 6 ms, that is consistent with the Mikron UCP 710 machining center. Furthermore, the trajectory for the desired part in Fig. 3.1.a is interpolated to machine with two linear axes Y and Z, and one rotary axis A.

Following, the machining of the same part is also simulated. Nevertheless, to evaluate the response of OGA in both different types of position controller and different contexts of multiaxis high speed machining, other configurations for the reference trajectory are proposed, as follows:

• Firstly, the sampling period in the framework of machine tools for the position loop and the setpoint interpolation period were set to 1 ms, instead of 6 ms. The reason for this configuration is that the GPC provides better performance with 1 ms than 6 ms [Rodriguez-Ayerbe 2014]. It can be understood that each advanced position controller has its own specificities, and its own applicable domain. This makes sense, because it allows to see the adaptation of OGA for whatever condition of axis control.

• The above modification related to the sampling period leads to the re-calculation in the feedrate interpolation stage, using VPOp algorithm with a new interpolation period of 1 ms. Unfortunately, this task for the case of 5-axis machining, i.e. for Y, Z, and A axes in the previous case of study, takes too long to finish. However, it can be finished with an reasonable computation time for a 3-axis machining, i.e. for Y and Z axes for this case of study. Due to the time limit in this PhD framework, finally it is proposed that instead of using 5-axis machining center, machining the above desired part with 3-axis machining center is now simulated.

One example for the reference trajectories obtained by the different configurations in CAM and in the feedrate interpolation is illustrated in Fig. 4.9. As it can be seen, the trajectory interpolated with T e = 1 ms has more programmed tool path data than with T e = 6 ms. Furthermore, it is clear that the tool axis orientation of the 3-axis machining is always fixed vertically, while that of 5-axis machining is changed, respecting the configuration of tilt angle θ t = 5 degrees and yaw angle θ n = 0 degree, as illustrated in Fig. 1.19.

Tool axis directions

Generalized predictive control (GPC)

Following, the synthesis of GPC controller is briefly reminded.

GPC belongs to the class of model-based advanced controllers. Therefore, the design of the GPC controller for the axis control requires the knowledge of a numerical model of the machine axis.

In this work, this is obtained through two successive identification procedures as shown in Fig. q -1 B(q -1 ) A(q -1 ) θ U ref Figure 4.11

4.11.

: Identification procedures of the axis models

The first identification process is to obtain the nonlinear axis model. Based on this nonlinear model, a linear discrete time transfer function of the axis dynamics, in other words the numerical model, is derived through a standard least square identification method [Landau 1990]. Once the appropriate linear axis model is obtained, the design of the GPC can be performed, as stated below.

In the GPC theory [Clarke 1987], the plant is classically modeled by the input/output CARIMA form as in (4.6), A q -1 y (t) = B q -1 u (t -1) + C (q -1 ) ξ (t) ∆ (q -1 ) (4.6) with q -1 the backward-shift operator, u(t) the input and y(t) the output of the system. ξ (t) is a zero mean non-correlated white noise, C (q -1 ) models the noise influence [Clarke 1989]. The introduction of the difference operator ∆ (q -1 ) = 1 -q -1 in the disturbance model helps to find an integral action in the controller and so eliminate the static errors. A (q -1 ) and B (q -1 ) are the two polynomials obtained by the least square identification:

   A (q -1 ) = 1 + a 1 q -1 + . . . + a na q -na B (q -1 ) = b 0 + b 1 q -1 + . . . + b n b q -n b .
The block diagram of the above CARIMA model is represented in Fig. 4.12.

u(t) y(t) B(q -1 ) q -1 A(q -1 ) 1 + + Δ(q -1 ) C(q -1 ) ξ(t) The control signal of GPC is obtained by minimization of a quadratic cost function in (4.7),

J = N 2 j=N 1 [y ref (t + j) -ŷ (t + j)] 2 + λ Nu j=1 ∆u (t + j -1) 2 (4.7)
in which

• N 1 and N 2 define the lower and upper output prediction horizons, respectively

• N u is the control horizon,

• λ is the control weighting factor,

• y ref is the setpoint,

• ŷ (t + j) is the optimal j-step ahead predictor, given in (4.8), ŷ (t + j) = F j q -1 y (t) + H j q -1 ∆u (t -1) free response + G j q -1 ∆u (t + j -1)

forced response (4.8) where the F j , H j , G j polynomials are unique solutions of a set of Diophantine equations [Clarke 1987],

• u is the control signal.

The receding horizon principle assumes that only the first value of the optimal control sequence resulting from the minimization of (4.7) is applied to the system, so that at the next sampling period the same procedure is repeated. This control strategy leads to a two-degree of freedom RST controller, implemented through a difference equation (4.9) by means of 3 polynomials R, S and T , as illustrated in Fig. 4.13. One major characteristic to be noted is that T (q) has a non-causal structure and assures the anticipation effect in the closed loop.

∆(q -1 )S q -1 u (t) = -R q -1 y (t) + T (q) y ref (t) (4.9) R(q -1 )

u(t) y(t) y ref (t) +

-Δ(q -1 )S(q -1 ) 1 B(q -1 ) q -1 A(q -1 ) 1 T(q) + + Δ(q -1 ) C(q -1 ) ξ(t) Once the polynomials R (q -1 ), S (q -1 ) and T (q) of the GPC controller are found, they are used to control the nonlinear axis model, e.g. of Y axis, shown in Fig. 4.14. R Y (q -1 ) Δ(q -1 )S Y (q -1 ) 1 + - without the OGA, while verifying the axis kinematic limitations in Table 2.2 and the motor current limitations in (2.33).

Equivalent polynomial controller

GPC controller

U Y,m (rad/s) (rad) Y m r s K global Y Y m s U Y,m r (q -1 ) (q -1 ) R Y (q -1 ) Δ(q -1 )S Y (q -1 ) 1 (rad/s) U Y,m

Simulations and discussions

The main objective of this section is to primarily prove that the OGA can be applied into diversified control structures, either classical or advanced ones and can reduce more or less the contour error as compared with the control structures without OGA.

Due to the above intention, and for the simplicity, only OGA(a4) is performed. Following, four cases of study are considered, as in Table 4.3.

Linear model

The linear model for the position loop is needed not only to determine the stability criterion of the axes but also to design the GPC controller. To obtain this, the same identification procedure as Chapter 3, Section § 1.2 is performed, but not represented here for the simplicity. It is noted that because the only difference is the sampling time of 1 ms, instead of 6 ms, a third-order linear velocity model is also obtained. Following, only the obtained mathematical representation of the linear model of the position open loop, illustrated in Fig. 4.16, is given in (4.10).

Y s m (q -1 ) Û r Y,m (q -1 ) = 0.000224q -1 -0.0001933q -2 1 -2.729q -1 + 2.444q -2 -0.6697q -3 -0.04502q -4 (4.10)

Reference responses

To perform the cases of study in Table 4.3, the reference cases of P+FFW and GPC controllers must be determined firstly.

Reference response of the classical axis controller (P + FFW)

As discussed in Chapter 3, Section § 1.3, the reference case of the classical axis controller corresponds to the response based on the best fixed gains, K P * and K F * . These gains are obtained by testing many values of K P and K F in order to find out the best ones having the smallest mean contour error, respecting all of the constraints.

In this case of study, the above searching approach is employed again to find out K P * and K F * .

These best fixed gains and the corresponding mean of contour error are given in Table 4.4 with T e = 1 ms are different from the ones in 5-axis machining with T e = 6 ms. This is due to the fact the position setpoints and its reference kinematic characteristics of these two cases are produced differently by the feedrate interpolation algorithm, as discussed and illustrated in Fig. 4.9 and Fig. 4.10. Therefore, it can be understood that although the trajectory curvature is the same, the control gains need to be consistent with the different kinematic references, i.e.

different feedrate profiles, to obtain a high contouring accuracy.

Reference response of the advanced controller (GPC)

The synthesis of GPC polynomials is based on the obtained numerical model of the position loop in (4.10), along with the pre-defined parameters including N 1 , N 2 , N u , λ, in (4.7) and the polynomial C in (4.6). These parameters are chosen respecting the stability constraints in (2.28).

In particular, it is proposed that the reference response of GPC corresponds to the best value of λ in (4.7), denoted as λ * a , which is achieved by the same manner as the best fixed gains, K P * and K F * , of the P + FFW controller. In this case of study, λ Y and λ Z are tuned within a set of values [0.0001 : 0.0001 : 0.01], in order to find out the best values, λ * Y and λ * Z , corresponding to the smallest mean contour error. Table 4.5 represents the values of all of the predefined parameters in the reference case of the GPC controller, along with its mean of contour error. N u 1

C (q -1 ) (1 -q -1 ) (1 -0. The synthesis of GPC leads to the following R, S and T polynomials, e.g. of Y axis, in (4.11)

                      
R Y (q -1 ) = 90.8526 -160.3818q -1 + 67.1472q -2 + 4.2234q -3 ∆ (q -1 ) S Y (q -1 ) = 1.0000 -1.9246q -1 + 1.1778q -2 -0.2262q -3

T Y (q) = 0.2534q -1 + 0.1279 + 0.0346q 1 -0.0475q 2 -0.1162q 3 -0.1712q 4 -0.2124q 5 -0.2404q 6 -0.2558q 7 -0.2596q 8 -0.2531q 9 -0.2373q 10 -14.2170q 11 + 17.4360q 12 (4.11)

Frequency analysis Fig. 4.17 shows the Black-Nichols diagrams of the position open loop to verify the stability criterion of the P + FFW and GPC controllers, using the above configurations for the reference cases. For both Y and Z axes, the gain margin and phase margin of the P + FFW controller are larger than the ones of the GPC controller. However, both of them respect the stability constraint in (2.28). For the OGA configurations, the parameters in Table 3.4 are kept unchanged, while the variation factor T 0 is chosen as 1, corresponding to case Ra1 in Fig. 3.12. With the interpolation period of 1 ms, the trajectory of the desired part in Fig. 3.1.a includes 4965 instants. The total length of the tool displacement is s end = 101.5 mm. The nine adaptive horizons generated by the Algorithm 1 have the adaptive N k and ∆s h,k as follows:

N k (h1 : h9) = [531,540,747,726,621,515,354,637,294] instant.

∆s h,k (h1 : h9) = [8. 5383, 10.7139, 12.8409, 17.2418, 12.7771, 11.9358, 9.5926, 12.2614, 5.1457] mm.

Based on the analysis in § 3.2, the admissible range of the global gain is obtained and given in Table 4.6. Only K global,min = 0.9 is checked to have a reasonable computation time, that is about several tenths of minutes. [0.9 : 1.16] [0.9 : 1.25] GPC+OGA(a4) [0.9 : 1.62] [0.9 : 1.24] As it can be seen, the global gain is adjusted within its admissible gain range, and all of the axis kinematic responses respect the kinematic limitations in Table 2.2. Moreover, it can be seen that for GPC and GPC +OGA(a4), the variation of jerk response is the least in all of the cases.

The comparison of the mean contour error between the above cases of study is given in Table 4.7. It can be seen that the application of OGA in the two cases f1b and f2b has reduced the contour error in the two reference cases f1a and f2a, with around 17 % and 39%, respectively.

Another remark is that the development of the advanced axis controller, GPC, can reduce 50% the contour error of the classical axis controller, P+FFW. While, if GPC is combined with OGA(a4), then it can even reduce around 70% the contour error of the classical controller. Moreover, Fig. 4.20 illustrates the resulting profiles of the variable global gain, the resulting feedrate, tool tracking error and contour error profiles as a function of trajectory.

It can be seen that due to the optimal adjustment of the global gain value, although the feedrate and the tool tracking error is slightly changed, the contour error has been much compensated.

Furthermore, as it can be seen in the two histograms in Fig. 4.21, the OGA(a4) allows to obtain a centered normal distribution of contour error, ranging from -0.8 to 0.8 µm for P+FFW+OGA(a4)

and from -0.2 to 0.2 µm for GPC + OGA(a4). While, without OGA, the norm of contour error ranges from 0.4 to 1.2 µm and from 0.2 to 0.4 µm for P+FFW and GPC, respectively.

Zoom in Zoom in

Tool TE

Conclusion

According to the above simulation results and discussions, it can be said that the OGA can be combined with the advanced controller by means of the global gain, in order to much more compensate the contour error.

It is important to remind that the above performances of OGA in these cases of study may not be the best ones, because the parameter tuning is not discussed, it only aims at proving that the OGA can apply and have its effect in different control structures, either classical or advanced one. In this aspect, the above simulation results have satisfied the initial objective.

OGA-GUI

To facilitate the end user in configuring and performing OGA, an OGA-GUI interface is developed and illustrated in Fig. 4.22. The working principle of this interface is explained below.

(3) (4) (1) (2)

(3) (1) (2)

Conclusions

This chapter has introduced a new idea related to the global gain and the application of OGA into this new control gain. The advantages and disadvantages of this proposed idea are given as follows.

• Advantages

The application of OGA for the global gain can lead to a high contouring accuracy, and the same computation time as when the OGA is applied into one single control gain K P or K F .

The general principle of OGA is not changed in this case of study. Only the admissible range for this new control gain is needed to define, but the idea and the constraints for such an admissible range is not changed.

OGA can be applied not only with a classical controller but also with an advanced one, e.g. the GPC controller in this case of study.

• Disadvantages A drawback of this method is that it is needed to slightly modify the software part in CNC, in order to implement this new control gain.

Although the synthesis of the advanced controller is out of the scope of this PhD thesis, it is also interesting to see that the GPC, in the framework of the advanced model-based controller, can give a better contouring accuracy than the classical axis controller for the position loop in CNC.

The OGA-GUI interface is also built to facilitate the end users to tune the parameters, to switch between the techniques and the control gains, and to start the OGA algorithm. The visual comparisons and analysis are also integrated in this interface. Finally, the user can easily achieve the resulting variable gain profiles for use in the future CNC, to pre-compensate the contour error in the machining process.

Conclusion and perspectives

Conclusions

This PhD thesis focuses on improving the contouring accuracy in machining free-form part in the context of multi-axis high speed machining, for either 3-axis or 5-axis high speed machining center, based on the means of an off-line process. Its main strategy is to exploit all possibilities of the classical control structure of axis drive, proposing the elementary base of a method that could be considered for an implementation in a commercial CNC.

Inspiring from the above motivations, the main contributions of this thesis are summarized as follows.

The first contribution is related to the contour error formulation in 5-axis point milling, based on ε C , using the proposed estimation approach for the tool contact points C C . In the literature, almost all of works proposed the contour error estimation approach based on ε L and/or (ε L , ε u ), using the tool path information (C L , u). This is because this information is available and can be accessed during machining process, while the determination of the tool contact points C C during machining stage is still an open issue. In fact, the use of ε L and/or (ε L , ε u ) is to estimate the effective cutting error (under cut or over cut), expressed by ε C . Therefore, the formulation of contour error by ε C , using directly the estimated tool contact points C C , is proposed to predict the effective cutting error more accurate than the one based on the tool path information, ε L and/or (ε L , ε u ).

In addition to this, the kinematics taken into account during in the feedrate interpolation stage have also been analyzed. They are the essential elements supporting generally the axis control and particularly fulfilled the development of the proposed adaptive control law. Especially, the axis kinematic constraints, including the axis velocity, acceleration, and jerk limitations, in the feedrate interpolation are used as one of the constraints in solving the optimization problem of the proposed axis control approach.

The thesis has scoped with the challenges in machining with high speed machining center, being the high feedrate and the axis kinematic constraints of the machine. Besides, the nonlinear characteristic derived from the static friction is also considered through the friction model.

The motor current limits have also been taken into account to assure a good behavior of the machine drives. Moreover, the closed structure of commercial CNC and axis drives have been also discussed and taken into account, in order to proposed an easily implementable solution in actual CNC systems.

An adaptive axis control law has been proposed for the commercial CNC based on an Off-line Gain Adjustment (OGA) approach, to compensate the contour error during machining. The main idea is to determine the optimal variable gain values used in the position loop of the classical axis control to reduce the contour error during the machining phase. This optimal gain calculation is done off-line by solving an optimization problem using a machine virtual environment. The optimization is formulated subject to the axis kinematic constraints, the stability criterion of servo drive, and the limits of motor current. The optimization result is the gain values of the position loop and feed forward for the considered machining profile. In this way the gains of the position controller and the feed forward are changed on-line in real time to reduce the contour error in the machining.

One of the important contributions of the current research is the development of two solving techniques for the optimization of OGA. They are based on the receding horizon principle.

The gain tuning is performed within a trajectory horizon, that is receded until the end of the trajectory.

In the first solving technique, OGA(a), the gain is obtained based on the pre-defined geometry horizon, considering the curvature characteristic of the free-form profile and based on the predefined gain function. When the horizon is receded without overlap, the optimization is repeated for the next horizon.

In contrast, the second solving technique, OGA(b), predicts the optimal constant gains over one horizon length but only the first optimal gain is kept and the horizon is receded one step ahead.

This solving technique has a more important computation burden than the first one but allows the control gain to vary more flexibly. That leads to fine-tune the resulting feedrate and as a consequence, this second solving technique allows to reduce the contour error much more than the first one. In the concerned cases of study, the best case of OGA(b) has improved 67.14% the contouring accuracy, compared with the best one of OGA(a). However, the former takes more than 9 hours to finish the OGA optimization than the later.

The two solving techniques are different in their own principles therefore their configurations are also set differently. The analysis of the parameter configuration on the performance of OGA have been discussed. The choices of the maximum number of potential gains, the gain variation factor, the horizon lengths and the weighting factor are important for the high performance of OGA.

Another contribution of this work is the proposition of a global variable gain in the axis controller. The influence of the position control gain variation and the feed forward gain variation in the contour error reduction has induced this new idea related to a global control gain. This gain can be integrated into whatever kinds of controllers, either classical or advanced ones. Afterwards, this new gain is used in the application of the proposed OGA approach. This method guarantees not only the improvement of the contouring accuracy but also a reasonable computation time. That extents the tractable possibilities of the proposed OGA approach in an industrial environment.

Simulation results based on a validated machining simulator showed that the precision level of the contour following and the computation time are different depending on the different solving techniques and the number of axes and control gains involved in the optimization of OGA. The results also verified that the responses of OGA have respected all of the concerned constraints. The relations between the obtained variable gains, axis jerk response, resulting feedrate and contour error profiles have been highlighted. Due to the high confidence of the machine simulator, it can be said that the contour error has been reduced significantly thanks to the proposed OGA approach.

In addition, a machine interface has been built (OGA-GUI) to make the use of OGA easier to users and increase the practical applications of the OGA in the machining community.

It can be concluded that the proposed OGA approach satisfies the initial objective of the PhD thesis, that is to reduce the contour error while keeping unchanged the classical axis control structure inside the current commercial CNC of machine tools. The only proposed modification on the future CNC is the extended functions, in which the gains are updated at each sampling time and the required memory to store the gain values for a given trajectory.

Perspectives

In perspective, the OGA will consider more nonlinear characteristics and disturbances effects, e.g. machine flexibilities, into the optimization problem under constraints of OGA. That can extend the effectiveness of the OGA.

Figure 1 :

 1 Figure 1: Processus d'usinage à grande vitesse simplifié

Figure 2 :

 2 Figure 2: Usinage général en machine outil [Dagiloke 1995] et paramètres d'usinage considérés dans ce travail de recherche

  entre deux laboratoires de recherche français: le Laboratoire des Signaux et Systèmes (L2S) de CentraleSupélec et le Laboratoire Universitaire de Recherche en Production Automatisée (LURPA) de l'ENS Paris-Saclay. Ce travail est enregistré sous N o 2014-812D -Projet OMEGA et soutenu par la fondation DIGITEO. La méthode proposée est appliquée au centre d'usinage 5 axes Mikron UCP 710, composé de trois axes linéaires et deux axes rotatifs. Cependant, elle peut être appliquée quelle que soit la structure du centre d'usinage 5 axes considéré. Deux briques logicielles élémentaires précédemment développées au LURPA en partenariat avec le L2S constituent le début du processus: La première est le modèle non linéaire des axes du centre UCP 710 Mikron, qui a été initialement proposé par M. Susanu dans le cadre de sa thèse de doctorat au L2S de 2002 à 2005 [Susanu 2005], puis enrichi et validé par D. Prévost dans sa thèse de doctorat au LURPA de 2008 à 2011 [Prévost 2011a]. Le modèle a fait ses preuves pour simuler le comportement réel du centre d'usinage Mikron UCP 710 avec une grande fiabilité.

Figure 3 :

 3 Figure 3: Vue globale de la stratégie proposée

  Tout d'abord, la chaîne d'usinage numérique est brièvement présentée. Les particularités de la génération de la trajectoire sont ensuite relevées. Les techniques d'interpolation temporelle classiques sont expliquées et mènent à la présentation de la technique choisie pour générer la trajectoire de l'outil et les mouvements articulaires.

  bout d'outil et en flanc d'outil. Pour obtenir de telles formulations d'erreur de contour, la géométrie de la trajectoire FAO est détaillée et certaines hypothèses sont proposées. En outre, la cinématique prise en compte lors de l'étape d'interpolation temporelle est également exploitée. Une loi de commande d'axe adaptative est proposée fondée sur l'approche Offline Gain Adjustment (OGA) pour compenser l'erreur de contour pendant l'usinage. Le diagramme général d'OGA est illustré dans la Fig. 4. Ajustement des gains hors-ligne Consigne de position + Utiliser le modèle non-linéaire d'axe + Estimer l'erreur de contour + Résoudre le problème sous contraintes

Figure 4 :

 4 Figure 4: Idée générale d'OGA

Figure 5 :

 5 Figure 5: Simplified high speed machining procedure

Figure 6 :

 6 Figure 6: Machining process on machine tool [Dagiloke 1995] and specified concerns in the research framework

  PhD work is carried out at the University of Paris-Saclay, France, within the cooperation between two French research laboratories: Laboratoire des Signaux et Systèmes (L2S), Centrale-Supélec and Laboratoire Universitaire de Recherche en Production Automatisée (LURPA), ENS Paris-Saclay. This work is registered under N • 2014-812D -Projet OMEGA and supported by DIGITEO foundation. The proposed method is specifically built for the 5-axis machining center Mikron UCP 710, consisting of three linear and two rotary axes. However, it is expected to be applicable to the structures of 5-axis machining centers. There are two starting points favoring this work. The first one is the nonlinear model of the translation axes of the Mikron UCP 710, that has been proposed by M. Susanu within her PhD thesis at L2S from 2002 to 2005 [Susanu 2005], enhanced and validated by D. Prévost in his PhD thesis at LURPA from 2008 to 2011 [Prévost 2011a].The model was proven to simulate the real behavior of the Mikron UCP 710 machining center with high confidence.

Figure 7 :

 7 Figure 7: Global view of this PhD thesis

Figure 1

 1 Figure 1.1: Digital chain of machining

Figure 1 . 3 :

 13 Figure 1.2: APT parameterization

Figure 1 . 4 :

 14 Figure 1.4: Point milling with a) hemispheric tool; b) toric tool; c) cylindrical tool

  Figure 1.5: Tool positioning parameters

Figure 1

 1 Figure1.6: The tool path planning[START_REF] Duc | [END_REF]

Figure 1

 1 Figure1.7: Tool movement envelope and C L cross-section[Kim 1995]

Figure 1 . 9 :

 19 Figure1.9: Scallop height of a flat surface machined by a hemispheric tool[Lin 1996]

  Linear interpolation allows to describe tool path through a series of tool positionings. The tool location points are under the form of linear segments (C L point -ISO 3592), while tool axis orientations are interpolated in workpiece frame, as illustrated in Fig. 1.10. Each tool positioning constitutes a block of machining program, including three position components (x, y, z), and three orientation components (i, j, k).

Figure 1

 1 Figure1.10: Tool path description by linear interpolation[Lavernhe 2006]

Figure 1 .

 1 Figure 1.11: Convex/concave part -C C and C L points [Bohez 2002]

Figure 1

 1 Figure 1.13: a) Reference tool positioning; b) Actual tool positioning with contour error

Figure 1 .

 1 Figure 1.14: Cutting error caused by tool orientation deviation

3 Figure 1

 31 Figure 1.15: Zero cutting error in the presence of tool path deviation

Figure 1 .

 1 Figure 1.18: Illustration of contour error sign

Fig. 1 .Figure 1

 11 Fig. 1.19 shows the programmed tool path, the estimated C C path and the estimated vectors.

Figure 1

 1 Figure 1.21: Over cut and under cut in flank milling

Fig. 1 .Figure 1 Figure 1

 111 Figure 1.22: Impact of tool location error on the effective cutting errors of 3-axis flank milling

Figure 1 .

 1 Figure 1.25: Schematic diagram of feedrate interpolation

Figure 2 .

 2 Figure 2.1: Axis control structure [Susanu 2005]

  2.3.

Figure 2 . 2 :

 22 Figure 2.2: Definition of impeller surface and the machining trajectory [Pechard 2011]

2. 4 .

 4 It is a combination of both viscous and Coulomb frictions. The detailed descriptions of the friction model as well as the motor drive model are given in Appendix A.

Figure 2 . 3 :Figure 2 . 4 :

 2324 Figure 2.3: Simulated and measured axis tracking errors of the impeller surface [Prévost 2011a]

+

  Use nonlinear axis model + Calculate/predict the contour error + Solve the contour error optimization problem under constraints Commercial CNC with extended functions

Figure 2 .

 2 Figure 2.5: Proposed offline execution

Figure 2 . 6 :

 26 Figure 2.6: Motivation of gain adjustment

Figure 2 . 7 :

 27 Figure 2.7: Control structure with OGA, e.g. for Y axis

  2.8.a. The horizon h k is characterized by N k instants or N k tool positions, represented by (C L , u). Its length is physically defined either as a function of time or as a function of tool displacement length, as proposed below.

Figure 2 .

 2 Figure 2.9: a) Predicted tool behaviors; b) Predicted contact points and predicted contour error



  Gain margin ≥ 15 (dB) P hase margin ≥ 70 (Deg)

Figure 2 .

 2 Figure 2.10: "Torque-Velocity" diagram of Siemens motor 1FT6084-8WF71

Figure 2 .

 2 Figure 2.11: Motor current response of OGA

Figure 2 .

 2 Figure 2.12: Reference and actual axis jerks

Figure 2 .

 2 Figure 2.13: Impacts of K P a and K F a on the reference velocity signal

Figure 2 .

 2 Figure 2.14: OGA(a): a) Horizon length; b) Tuning cases; c) Receding horizon without overlap

Figure 2 .

 2 Figure 2.15: Gain prediction by predefined gain function

( 2 .Figure 2 .

 22 Figure 2.16: Effects of different values of l 1

Figure 2 .

 2 Figure 2.17: Value of n as a function of K k a

Figure 2 .

 2 Figure 2.18: Three examples for the tuning case of T 0 = 1

Figure 2 .

 2 Figure 2.19: Three examples for the tuning case of T 0 = 0.5

Figure 2 .

 2 Figure 2.20: OGA(b): a) Horizon length; b) Tuning the constant gains; c) Receding horizon with overlap

  0δ a . . . (n -1) δ a nδ a

4 . 3

 43 Chosen algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Fig. 3 .

 3 Fig. 3.1.a shows the desired part having free-form profile for milling with 5-axis Mikron UCP 710 machine center.

Figure 3

 3 Figure 3.1: a) Desired part; b) Feedrate cartography from feedrate planning

. 2 .Figure 3 . 2 :

 232 Figure 3.2: Nonlinear model of the position open loop without control gain (T e = 6 ms)

Fig. 3 .

 3 Fig. 3.3 shows the comparison between the step velocity responses of the different order linear models and that of the nonlinear one, for example in the case of Y axis.

Figure 3 . 3 :

 33 Figure 3.3: Step response of different order linear models of velocity loop

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Position open loop with the control gain K P Y (m/min/mm) and its linear plan

Figure 3 . 6 :

 36 Figure 3.6: Discretized admissible gain range by the gain steps: ∆K P = 0.1 and ∆K F = 0.01

Fig. 3 .

 3 Fig. 3.7 shows the curvature, feedrate (F) and contour error (CE) profiles as a function of the trajectory geometry for the above three cases of study.

Figure 3 . 7 :

 37 Figure 3.7: Responses with the fixed control gain values

Figure 3 . 8 :

 38 Figure 3.8: Comparing contour error profile of C3 and the nominal case

Figure 3 . 9 :

 39 Figure 3.9: Comparing velocity reference of FFW and P actions, for Y axis

Figure 3 .

 3 Figure 3.11: Illustration of the effect of the bad horizon length in OGA(b)

3. 12 .

 12 In fact, T 0 is tuned over the range [0.0001 : 1], however only five cases of study of OGA(a), including Ra1, Ra2, Ra3, Ra4 and Ra5, and five cases of study of OGA(b), including Rb1, Rb2, Rb3, Rb4 and Rb5, are shown, because they are enough to highlight the effect of the variation factor T 0 on the performance of OGA(a) and OGA(b).

Figure 3 .

 3 Figure 3.12: Illustration of the cases of study for the variation factor of OGA(a) and OGA(b)

Fig. 3 .Figure 3 .

 33 Fig. 3.14 does the same but for the best cases of OGA(b), meaning OGA(b1) in Rb1, OGA(b2) in Rb4 and OGA(b3) in Rb3.

Fig.

  Fig. 3.15 and Fig. 3.16 show that all of the axis kinematic responses respect their limitations in

Figure 3 .Figure 3 .Figure 3 .

 333 Figure 3.15: Axis kinematic responses of OGA(a3)

Figure 3 .Fig. 3 .

 33 Figure 3.19: Analysis of OGA(a1) in Ra1 and Ra2

Figure 3 .

 3 Figure 3.20: Analysis of OGA(b1) in Rb1 and Rb2

Fig. 3 .Figure 3 .Figure 3 .

 333 Figure 3.21: Impact of variable gains of OGA(a3) on the resulting feedrate, tool tracking error and contour error

Figure 3 .Figure 3 .

 33 Figure 3.23: Comparison of the contouring accuracy between OGA(a3) and OGA(b3)

Fig. 3 .

 3 Fig.3.24.b illustrates that the contour error of the reference case r0 mainly ranges from 0.5 to 50 µm. While, OGA(a3) and OGA(b3) have a centered normal distribution of contour error, ranging from -3 to 3 µm and from -0.5 to 0.5 µm, respectively.

Figure 3 .

 3 Figure 3.25: Impact of the sinusoidal part of gain function on OGA(a2)

Fig. 3 .

 3 Fig.3.26 shows the effects of the obtained gain profiles on the axis jerk responses and the resulting contour error profiles. In general, from Hb1 to Hb4, when the horizon length is increasing, the gain variations are generally decreasing and therefore the variations of the axis jerk responses are decreasing also.

Figure 3 .

 3 Figure 3.26: Impact of the prediction horizon length on the OGA(b1)

Figure 3 .

 3 Figure 3.27: Impact of the weighting factor on the OGA(a2)

(

  P+FFW) of the classical axis control in the current commercial CNC. This method can be easily implemented by the machine tool manufacturers, because they do not need to change their control structure. The only requirement is to change the way in which the control gain values are stored in CNC, a set of variable gain values instead of a fixed gain value.

Figure 4 .

 4 Figure 4.1: Adjustment of the classical control gains in position controller

Figure 4 . 2 :

 42 Figure 4.2: Adjustment of the global gain

2Figure 4 . 3 :

 43 Figure 4.3: Axis control structure with specificities of global gain and OGA

  4, except for the admissible gain range of K global a . There are many possibilities to choose the K global,min a (4.5), but only K global,min a = 0.5 is checked considering the reasonable computation time in the simulations. The K global,max a is obtained according to the discussion in § 2.1. Finally, the admissible range of the global gain is obtained and provided in Table 4.1.Table 4.1: Admissible range of the global gain

Fig. 4 .Figure 4 . 4 :

 444 Fig. 4.4 compares the resulting feedrate, tool tracking error and contour error profiles between OGA(a4) in Ra3 and OGA(b4) in Rb4, in the relation with the trajectory and the obtained variable global gains.

Fig. 4 .

 4 Fig. 4.5 compares the obtained variable global gains and the jerk responses of OGA(b4) in Rb1 (T 0 = 0.1) and Rb4 (T 0 = 0.001).

Figure 4 . 5 :

 45 Figure 4.5: Resulting variable global gains and jerk responses of OGA(b4) in the cases of Rb1 and Rb4

Figure 4 . 6 :

 46 Figure 4.6: Contouring accuracy of OGA(b4) in the cases of Rb1 and Rb4

Fig. 4 .Figure 4 . 7 :Figure 4 . 8 :

 44748 Fig. 4.7 and Fig. 4.8 show the resulting variable global gain profiles and the axis kinematic responses as a function of trajectory in cases of Ra3 and Rb4, which are the two best cases of OGA(a4) and OGA(b4), respectively.

Figure 4 Figure 4 .

 44 Figure 4.9: 3-axis trajectories with different interpolation periods

Figure 4 .

 4 Figure 4.12: Block diagram of CARIMA model

Figure 4 .

 4 Figure 4.13: Two-degree of freedom RST controller

Figure 4 .Figure 4 .

 44 Figure 4.14: GPC controller for the position loop of the nonlinear Y-axis model

  open loop with GPC controller and linear model

Table 4 . 5 :

 45 Parameters and ε in the reference response of the GPC

  8q -1 ) (1 -0.9q -1 ) above polynomial C is chosen considering[Rodriguez-Ayerbe 2014], in which the obtained GPC controller for axis control has been experimentally validated. Such a polynomial C has been chosen due to the two following reasons. Firstly, it includes a root at q = 1, in order to remove the integral action of the GPC controller. In fact, the static friction in the axis produces oscillations in the output when an integral action is included in the predictive controller [Rodriguez-Ayerbe 2014]. Secondly, the other two roots of the C polynomial permits to obtain good robustness margins [Rodriguez-Ayerbe 2005].

Figure 4 .

 4 Figure 4.17: Frequency analysis of f1a and f2a: a) Y axis; b) Z axis

3. 3 . 4

 34 Responses of the OGA integrated GPCThe resulting variable global gain profiles and the axis kinematic responses of all of the cases of study in Table4.3 are firstly shown in Fig.4.18 and Fig. 4.19.

Figure 4 .Figure 4 .

 44 Figure 4.18: Kinematic responses of f1a and f1b

Figure 4 .Figure 4 .

 44 Figure 4.20: Comparisons of feedrate, tool tracking error, contour error between f1a, f1b, f2a and f2b

Figure 4 .

 4 Figure 4.22: OGA-GUI interface

Figure 4 .Figure 4 .

 44 Figure 4.23: OGA-GUI: Select machine structure

Figure 4 .

 4 Figure 4.25: OGA-GUI: Reference trajectory and reference axis kinematic

Figure 4 .

 4 Figure 4.26: OGA-GUI: Select the position controller and design it

Figure 4 .

 4 Figure 4.27: OGA-GUI: Select adjustable gain(s)

Figure 4 .

 4 Figure 4.28: OGA-GUI: Select the axes for optimization by OGA

Figure 4 .

 4 Figure 4.29: OGA-GUI: Select OGA(a) or OGA(b) and configure it

Figure 4 .

 4 Figure 4.30: OGA-GUI: Add the constraints

Figure 4 .

 4 Figure 4.31: OGA-GUI: Run OGA

Figure 4 .

 4 Figure 4.32: OGA-GUI: OGA results

Tool displacement s (mm) Upper limits by acceleration Upper limits by jerk Feedrate (m/min) Actual feedrate F prog

Figure 1.26: Resulting feedrate by VPOp [Beudaert 2013]

  was a pioneer of Cross Couple Control (CCC) for compensating the contour error online. Afterwards, a large number of extended and/or improved versions of CCC are proposed[Huo 2012, Zhang 2015]. Cheng et al.[Cheng 2009] proposed an integrated motion control scheme combined with a fuzzy logic-based feedrate regulator, which adjust the value of the desired feedrate, to reduce the contour error. In addition, Khal-

ick et al.

[START_REF] El-Khalick | Discrete-time model predictive contouring control for biaxial feed drive systems and experimental verification[END_REF]

added the real time estimated contour error into the optimization problem of Model Predictive Control (MPC), in yielding the contour error reduction. Adaptive robust control is also proposed for biaxial contouring control in

[START_REF] Yao | [END_REF], Davis 2015

].

Table 2

 2

		Controller	Trajectory
		Tracking following	Contouring following	(C L , u)
			Cross couple control
	On-line	Sliding mode control [Kamalzadeh 2007], [Altintas 2000b], Zero Phase Error Tracking Controller [Erkorkmaz 2001b], etc.	[Koren 1980, Koren 1991] [Huo 2012, Zhang 2015], Fuzzy logic control [Cheng 2009], Model predictive control [El-Khalick 2011], Adaptive robust control
			[Yao 2012],[Davis 2015], etc.
				MPC's idea[Yang 2015b],
				FFW shaping[Khoshdarregi 2014],
	Off-line	Model predictive control [Rodriguez-Ayerbe 2014] [Dumur 2008], etc.	Off-line gain adjustment (OGA)	Analytical method and kinematic constraints [Zhang 2013],

.1: Contouring control methods and the intention of this work

  18)From(2.18), solving the FKT problem allows one to obtain a vector representing the predicted tool behaviors, Ps , given in(2.19),

			Ps = P s,k+1 w	, P s,k+2 w	, ..., P s,k+N k w	(2.19)
	where P s,k+1 w	represents for (	Ĉs,k+1 L

  2.22) From (1.18) to (1.20),(2.20) and (2.21), the predicted contour error values over the horizon h k are calculated and represented by ÊC in(2.23).

	ÊC = εk+1

C , εk+2 C , ..., εk+N k C (2.23)

From (2.16), (2.17), and (2.23), the cost function, J, is proposed in (2.24),

Table 2 .

 2

	2, denoted by V max a	, A max a	, and J max a	respectively,

Table 2 .

 2 2: Axis kinematic constraints of Mikron UCP 710 machine

	X Y Z	A	C

  2.14.a, a gain modification profile is defined to generate variable values for the control gain K a of the axis a, which are the gain values over the columns of either

	K P in (2.10) or K F in (2.11), as illustrated in Fig. 2.14.b. This gain profile is a function of several
	parameters. The principle of OGA(a) is to tune the parameter values of such a gain profile. By
	this way, different parameter values generate different responses of the gain profile, yielding

different values of the control gains. Then, by evaluating the optimization problem in (2.27), the best parameter value of the gain function is chosen, corresponding to the optimal variable gains

Table 3 .

 3 1: Admissible gain range of K P a (m/min/mm or rad/min/mrad) and K F

	a

Table 3 .

 3

		2: Three cases of study for using fixed gains
	Cases of study	Gains K P * (m/min/mm or rad/min/mrad) 1.2 1.8 Y Z	A 1	ε (µm)
	(1)	K F 1	0.95 0.85 0.95	33.3
	(2)	K F 2	0.85 0.75 0.85	18.5

Table 3 .

 3

	3: Nominal gain values in position controller of Mikron UCP 710 machine
	Y	Z A
	K P,n a

Table 3 .

 3 4: Configurations of OGA(a) and OGA(b)

	Type	Parameters		Data source	Value	Unit
		T e				Mikron machine	0.006	s
		[I n Y , I n Z , I n A ] (2.37)	Motor documentation [26.2, 26.2, 42.7] A
	General	[∆I max Y	, ∆I max Z	, ∆I max A ](2.38) User-defined	2% [I n Y , I n Z , I n A ]	A
		λ P , λ F (2.24)		User-defined	0	-
		K P,1 a , K F,1 a (2.45)	User-defined	K P * a , K F * a	see Chap. 2, § 3.2
		K P,min,max a	, K F,min,max a	User-defined	see Table 2.2	see Table 2.2
		N k , ∆s h,k (2.7)		Algorithm 1	( + ) see below	instant, mm
	OGA(a)	l 1 (2.48)			User-defined	0.3N k	instant
		p 0 (2.53)			User-defined	10	-
	OGA(b)	N k , ∆t h,k (2.6) p 0 (2.53)	User-defined User-defined	25, 0.144 5	instant, s -

Table 3 .

 3 5: Cases of study for OGA with the classical control structure

	Control strategy Notation Description
	Classical control r0	Using the best fixed gains
		OGA(a1) Tuning K P a
	OGA(a)	OGA(a2) Tuning K F a
		OGA(a3) Tuning both K P a and K F a
		OGA(b1) Tuning K P a
	OGA(b)	OGA(b2) Tuning K F a
		OGA(b3) Tuning both K P a and K F a

Table 3 .

 3 8: The mean contour errors of OGA(a) and OGA(b)

	ε(µm)

Table 3 .

 3 9: The best performances of OGA(a) and OGA(b)

		ε(µm) Improvement Computation time
	r0	14.4	-	-
		OGA(a1) 4.9	65.97% r0	≃ 12.5 min
		OGA(a2) 1.4	90.28% r0	≃ 12.5 min
	OGA(a)		71.43%a1
		OGA(a3) 0.7	95.14% r0	≃ 2260 min (≃ 1.5 days)
			50%a2
		OGA(b1) 4.2	70.83% r0	≃ 6.5 h
		OGA(b2) 0.68	95.28% r0	≃ 6.5 h
	OGA(b)		83.81%b1
		OGA(b3) 0.23	98.40% r0	≃ 45.5 h (≃ 1.89 days)
			66.18%b2
			67.14%a3

Table 3 .

 3 10: Cases of study for the sinusoidal part

	Control strategy Notation l 1	ε(µm) Improvement
	Classical control r0	-	14.4	-
		L1	0.1N k 5.7	60.42%r0
	OGA(a2)	L2	0.3N k 1.4	90.28%r0 75.44%L1
		L3	0.8N k 1.4	90.28%r0

Table 3 .

 3 11: Cases of study for the horizon length of OGA(b)

	Control strategy Notation Horizon length	ε(µm) Improvement
	Classical control r0	-	14.4	-
		Hb1	0.024 s (5 instants)	4.7	67.36%r0
	OGA(b1)	Hb2 Hb3	0.144 s (25 instants) 4.2 0.294 s (50 instants) 4.8	70.83%r0 66.67%r0
		Hb4	0.594 s (100 instants) 6.2	56.94%r0

Table 3 .

 3 11 indicates that in general the above four cases of study produce smaller contour error than the reference case r0. The best case is Hb2, which has reduced around 70% the contour error of r0.

Table 3 .

 3 12: Cases of study for the weighting factor of OGA(a) and OGA(b)

	Control strategy Notation λ F	ε(µm) Improvement
	Classical control r0	-	14.4	-
		Wa1	0	1.4	90.28%r0
	OGA(a2)	Wa2	0.0001 3	79.17%r0
		Wa3	0.001 5.8	59.72%r0
		Wb1	0	6.3	56.25%r0
	OGA(b2)	Wb2	0.0001 4.7	67.36%r0
		Wb3	0.001 4.7	67.36%r0

  4.1.

		Variable gain 1
	θ ref (rad)	Variable gain 2	...	U ref (rad/s)	Velocity loop (Current loop: Motor + Friction model)	V (rad/s)	1 s	θ (rad)
		Variable gain n
		Position controller
						Position loop

Table 4 .

 4

		2: Performance of OGA(a4) and OGA(b4)
	r0	ε(µm)			14.4
			Ra1	Ra2	Ra3	Ra4	Ra5
		T 0	1	0.2	0.1	0.05	0.01
	OGA(a4)	ε(µm)	7.2	3.6	3.3	4.2	8.7
		Improvement (%)	50% r0	75% r0 77.08% r0 70.83% r0 39.58% r0
		Computation time (min)			≈ 12.5
			Rb1	Rb2	Rb3	Rb4	Rb5
		T 0	0.1	0.05	0.01	0.001	0.0001
	OGA(b4)	ε(µm) Improvement (%)	25.6 -77.78% r0 9.03% r0 60.42% r0 79.17% r0 18.75% r0 13.1 5.7 3 11.7
		Computation time (min)			≈ 400

Table 4 .

 4 3: Cases of study for different axis controllers (T e = 1 ms)

	Control strategy	Notation Tuning Gain
	P+FFW	f1a	None (Fixed gains)
	P+FFW+OGA(a4) f1b	K global a
	GPC	f2a	None (Fixed gains)
	GPC+OGA(a4)	f2b	K global a

Table 4 .

 4 4: K P * (m/min/mm), K F * and ε in 3-axis machining with T e = 1 ms

	Gain	Y	Z
	K P *	1.5 2.3
	K F *	0.72 0.57
	ε (µm)		0.46
	Comparing

Table 3

 3 .2 and Table4.4, it can be noticed that K P * and K F * in 3-axis machining

Table 4 .

 4 6: Admissible range of the global gain with T e = 1 ms

	Gain	Case	Y	Z
	K global	r0+OGA(a4)

Table 4 .

 4 7: Contouring performances comparison

	Control strategy	Notation Mean contour error Improvement	Computation time
			(µm)	%	(min)
	P+FFW	f1a	0.46	-	-
	P+FFW+OGA(a4) f1b	0.38	17.39%f1a	≈ 12.5
	GPC	f2a	0.23	50%f1a	-
	GPC+OGA(a4)	f2b	0.14	69.57%f1a, 39.13%f2a ≈ 12.5

Acknowledgments

First, I would like to express my special thanks to my supervisors: Professor Pedro RODRIGUEZ-

* The prediction effect of OGA(b) is more global than that of OGA(a). It is not necessary to analyze the curvature characteristics of the free-form trajectory before performing OGA (b). * According to the above results, it can be said that the obtained variable gains in OGA(b) are more flexible than in OGA(a). Moreover, the best cases in OGA (b) produce smaller contour error than in OGA(a).

The challenges include: * Because the horizon is receded by only one instant, the computation time in OGA(b) is much longer than in OGA(a). It can be argued that if the horizon is receded by more than one instant, then the computation time can be reduced. This is correct. However, initially the OGA(b) is inspired from the idea of MCP, in which one-step receding horizon is chosen. Moreover, if the receding horizon length can be varied, it will be another parameter for tune. Therefore, this will be a part of future works. * If the variable gains obtained by OGA(b) is implemented into the real machine, it takes more memory space than the gain implementation of OGA(a). * How to explicitly choose optimal horizon length of OGA(b) is still an open issue.

In addition, the above two solving techniques, OGA(a) and OGA(b), share the same characteristics as follows.

Their common advantages include: * Both techniques reduce much more contour error than the reference case using the best fixed gain in the classical axis control.

The application of OGA for the multi-axis flank milling would also be performed in future works.

In addition, other techniques of optimization, e.g. heuristic method or genetic algorithm, etc, should be examined. The objective is to reduce the computation burden in solving the OGA optimization problem and to produce more optimal control gains for the contour error reduction.

The gain modification pattern depends on all axis trajectories. This dependency has to be deeper analyzed to reduce the computation time. Ideally, it is expected to analyze the obtained gain variation profiles to eliminate or reduce the optimization stage. Moreover, using optimization results of similar trajectories is to produce gain modification profiles in function of trajectories.

Last but not least, the experimental tests of OGA are expected to be performed in the framework of an OPEN CNC.

Appendices

Nonlinear axis model 1 Motor model The back EMF voltage E(t) is proportional to the angular velocity Ω(t) (rad/s) of the rotor in the motor, expressed as (A.1),

where K e (V/rad/s) is the back EMF constant.

According to the Kirchhoff's voltage law, the electrical equation of the DC motor is described by (A.2),

where i(t) (A) is the armature current. This current passing through the rotor windings in the magnetic region of the stator generates the Laplace force. As a result, this force produces a torque C m (N.m), leading to the rotation of the rotor, given as (A.3),

where K t (N.m/A) is the torque constant. Such a torque drives the external torque through the motor shaft, expressed by the dynamic equation in (A.4),

where C res (t) (N.m) represents the resistant torque, which is due to the viscous and dry friction or the gravity of the Z axis; J (kg.m 2 ) the total equivalent rotor moment of inertia, that is obtained by (A.5),

where J mot , J red , J screw and J axis are the equivalent inertia of the motor, the gear reducer, the ball screw and the driven axis; N is the reduction factor.

Considering all of the equations from (A.1) to (A.5) in the Laplace domain s, the electric, eletro-mechanic and mechanic equations are obtained in (A.6), (A.7) and (A.8), respectively.

The equations from (A.6) to (A.8) can be illustrated by the block diagram in Fig.

The above model is used for a rotary motor. For the linear motor, the detailed description can be found in [Prévost 2011b]. The obtained friction model is characterized by (A.9),

where i f r (A) is the current caused by the friction, V (m/min) is the axis velocity, the parameters a, b, c, d and the value of the static friction current i 0 (A) are given in Table A.1.

Experimental points Tendency curve

Results (X axis) Double exponential law C m (s) = s.J.Ω(s)

The motor model with the friction model (A.9) is illustrated in Fig Title : Variable gain contouring control for multi-axis machine tools Keywords : adaptive control laws; contour error; high speed machining Abstract : The advanced machining techniques are always the backbone of the manufacturing industries. Among such techniques, high speed machining is the main subject of this PhD thesis. Indeed, the main objective of this work is to improve the contouring accuracy in multi-axis high speed machining of free-form surfaces, directly acting inside the axis control loops. To do that, a first step aims at elaborating a strategy to estimate as accurately as possible the contour error for different tool configurations. This contour error is then minimized by means of an off-line adaptation for a given profile of the proportional and feedforward gains of the axis position loop controllers. This gain adaptation is performed via an optimization algorithm that considers a nonlinear model of the machine behaviour, in particular including friction related to each axis. This optimization leading to the controllers gains takes into account several constraints, including the axis kinematic (velocity, acceleration and jerk) limitations, the stability of the controlled loops and the motor current limits. Finally, to help their integration within an industrial framework, the developed strategies can be directly implemented in commercial CNC, by exploiting all possibilities of the classical control structure of axis drive.