TRAJECTORY GENERATION, CONTROL, AND GEOMETRIC ERRORCOMPENSATION FOR A 9-AXIS MICROMACHINING CENTERbyAlexander YuenB.A.Sc., The University of British Columbia, 2011M.A.Sc., The University of British Columbia, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)July 2018c© Alexander Yuen 2018The following individuals certify that they have read and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, a thesis/dissertation entitled:TRAJECTORY GENERATION, CONTROL, AND GEOMETRIC ERROR COMPENSATIONFOR A 9-AXIS MICROMACHINING CENTERsubmitted by Alexander Yuen in partial fulfillment of the requirements for the degree of Doctorof Philosophy in Mechanical EngineeringSupervisory Committee Members:Yusuf AltintasSupervisorTim SalcudeanSupervisory Committee MemberXiaodong LuSupervisory Committee MemberElizabeth CroftSupervisory Committee MemberiiAbstractThis thesis presents a trajectory generation algorithm, a control strategy, and a geometric errorcompensation methodology for a novel 9-axis micromachining center which combines a 3-axis micromill with a 6 degree of freedom magnetically levitated rotary table. The proposedtrajectory generation algorithm resolves redundant degrees of freedom by numerically solvingfor axes positions from desired tool positions and orientations. Differential axes positions arefound while ensuring the stroke limits of the drives are respected and singularities are avoided.The differential solution is numerically integrated to obtain the axes positions with respect todisplacement. The axes commands are then scheduled in time, while respecting the velocity,acceleration, and jerk limits of each of the drives, and traversing the toolpath as fast as possible.The experiments showed trajectories that resolved redundancies, avoided singularities, andrespected all physical limits of the drives.A control strategy which combines the capabilities of the micromill and the rotary table isintroduced. A sliding mode controller with a LuGre friction compensator is designed to controlthe position of the micromill, based on identified physical parameters. A lead-lag positioncontroller with an integrator and a notch filter is designed to control the rotary table. Sincethe translational axes of the micromill and rotary table are in parallel, the tracking error of themicromill is sent as a reference command to the rotary table, compensating the tracking errorsof the micromill with the higher bandwidth of the rotary table. In experiments, the dual stagecontrol law improved tracking error over the micromill alone.iiiThe geometric errors of the 3-axis micromill is compensated by using the precision motion ofthe 6 degree of freedom rotary table. The geometric errors of the 3-axis micromill are mea-sured with a laser interferometer, fit to quintic polynomials, and incorporated into the forwardkinematic model. The tooltip deviation is found by subtracting the ideal tooltip position fromthe tooltip position affected by geometric errors. Rotary table commands, from all 6 axes, thatcompensate for these deviations are found using a gradient descent algorithm. Experimentsshowed reductions in end effector deviations.ivLay SummaryTechnological advances in industries such as the electronics and biomedical field has driventhe demand for manufactured parts with high precision features, typically in the order of onethousandths of a millimeter. Subsequently there is a demand for manufacturing processesand machine tools capable of generating high precision features. A hybrid 9-axis machinetool, which combines a 3-axis micromilling machine and a 6 degree of freedom magneticallylevitated rotary table has been developed for this purpose.This thesis presents algorithms to plan the motion, control the precision, and control the ac-curacy for this 9-axis machine tool. The motion is planned to ensure all 9-axis are used asefficiently as possible given a predefined path. Furthermore the machine is controlled so theprecision and accuracy of the 3-axis micromill is improved with the rotary table. The developedalgorithms can be used with similar machines to manufacture high precision parts.vPrefaceChapter 3 and 4. A version of this material has been published in the CIRP Annals [YuenA., Altintas Y. Trajectory generation and control of a 9-axis CNC micromachining center.CIRP Annals. doi: 10.1016/j.cirp.2016.04.098] and IEEE/ASME Transactions on Mechatron-ics [Yuen A., Altintas Y. Constrained Trajectory Generation and Control for a 9-Axis Micro-machining Center With Four Redundant Axes. IEEE/ASME Transactions on Mechatronics.doi: 10.1109/TMECH.2017.2771260]. I was the lead investigator, responsible for the ma-jor development of the algorithms, experiment design, data collection and analysis, as wellas manuscript composition. Altintas Y. was the supervisory author on this project and wasinvolved throughout the project in concept formation and manuscript compositionChapter 5. A version of this material has been published in the ASME Journal of Manufac-turing Science and Engineering [Yuen A., Altintas Y. Geometric error compensation with a6 degree of freedom rotary magnetic actuator. ASME Journal of Manufacturing Science andEngineering. doi:10.1115/1.4040938]. I was the lead investigator, responsible for the ma-jor development of the algorithms, experiment design, data collection and analysis, as wellas manuscript composition. Altintas Y. was the supervisory author on this project and wasinvolved throughout the project in concept formation and manuscript compositionChapter 6. The results found in this chapter were previously included in the above three men-tioned works. Modifications include a more detailed analysis of the results. All figures andtables found in this thesis are used with permission from applicable sources.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Trajectory Generation for Redundant Manipulators and Machine tools . . . . . 52.3 Control laws for Feed drives and Dual Stage Actuators . . . . . . . . . . . . . 82.4 Geometric Error Modeling, Identification, and Compensation . . . . . . . . . 113 Trajectory Generation for a 9-axis Micromachining Center . . . . . . . . . . . . 133.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.2 Forward Kinematics of the 9-axis Micromachining Center . . . . . . . . . . . 14vii3.3 Redundancy Resolution of the 9-axis Micromachining Center . . . . . . . . . 193.4 Avoiding Stroke Limits of the Drives . . . . . . . . . . . . . . . . . . . . . . 213.5 Singularity Avoidance with augmentation of Jacobian Matrix . . . . . . . . . 253.6 Numerical Integration of Differential Solution to get Position Commands . . . 273.7 Closed loop correction of deviations from numerical integration . . . . . . . . 283.8 Full Redundancy Resolution Algorithm . . . . . . . . . . . . . . . . . . . . . 293.9 Feedrate Optimization with Redundancy Resolution . . . . . . . . . . . . . . 324 Control Design for Micromill Feed Drives, High Precision Rotary Table, and DualStage Feed Drives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.2 Modeling and Control Design of the Micromill . . . . . . . . . . . . . . . . . 394.2.1 Identification of Rigid Body Dynamics . . . . . . . . . . . . . . . . . 394.2.2 Identification of Friction Characteristics . . . . . . . . . . . . . . . . 434.2.3 Sliding Mode Position Controller . . . . . . . . . . . . . . . . . . . . 514.2.4 LuGre Feedforward Friction Compensator . . . . . . . . . . . . . . . 544.3 Modeling and Control Design of the Rotary Table . . . . . . . . . . . . . . . . 594.3.1 Lead-lag Position Controller . . . . . . . . . . . . . . . . . . . . . . . 594.3.2 Notch Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.4 Dual-stage Feed Drive Control . . . . . . . . . . . . . . . . . . . . . . . . . . 635 Geometric Error Modeling for a 3-axis Micromill and Compensation with a 6Degree of Freedom Rotary Table . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.2 Kinematic Model with Errors . . . . . . . . . . . . . . . . . . . . . . . . . . 685.3 Error Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71viii5.4 Geometric error compensation . . . . . . . . . . . . . . . . . . . . . . . . . . 745.4.1 Gradient Descent Optimization Algorithm Background . . . . . . . . 745.4.2 Geometric Error Compensation using Gradient Descent OptimizationAlgorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756 Simulation and Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . 806.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.2 Trajectory Generation Experimental Results . . . . . . . . . . . . . . . . . . . 806.3 Dual Stage Feed Drive Tracking Control Results . . . . . . . . . . . . . . . . 926.4 Geometric Error Compensation Results . . . . . . . . . . . . . . . . . . . . . 957 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.2 Future Research Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109AppendicesA Foward Kinematic Equation with Geometric Errors . . . . . . . . . . . . . . . . 117ixList of Tables4.1 Identified rigid body parameters of the micromill . . . . . . . . . . . . . . . . 434.2 Identified friction parameters of the micromill . . . . . . . . . . . . . . . . . . 504.3 Micromill parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.4 Rotary table parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 626.1 Circular contouring results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 926.2 Square contouring results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 926.3 Geometric errors with and without compensation . . . . . . . . . . . . . . . . 956.4 Mean and maximum tooltip errors with and without compensation for a multi-axis trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101xList of Figures1.1 UBC MAL’s hybrid 9-axis machine . . . . . . . . . . . . . . . . . . . . . . . 23.1 9-axis micromachining center and corresponding axes of the machine tool . . . 153.2 Generated trajectories without constraint cost function . . . . . . . . . . . . . 233.3 Generated trajectories with constraint cost function . . . . . . . . . . . . . . . 233.4 The norms of dqdswith and without constraint cost function . . . . . . . . . . . 243.5 Desired toolpath and generated tool path with w0 = µ0 = 1 . . . . . . . . . . . 263.6 Desired toolpath and generated tool path with w0 = µ0 = 10−6 . . . . . . . . . 263.7 Numerical errors of the redundancy resolution algorithm with closed loop cor-rection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.8 Numerical errors of the redundancy resolution algorithm without closed loopcorrection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.9 Real-time trajectory generation command from the result of feed optimizationalgorithm and 9-axis kinematic module . . . . . . . . . . . . . . . . . . . . . 364.1 Frequency rich signal for rigid body identification . . . . . . . . . . . . . . . . 424.2 Disturbance estimation for xc-axis at 10 [mm/s] . . . . . . . . . . . . . . . . . 484.3 Average disturbance estimation for xc-axis at all test speeds . . . . . . . . . . . 484.4 Disturbance estimation for xc-axis with corrected viscous friction and the fittedfriction curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.5 Graphical representation of LuGre friction model . . . . . . . . . . . . . . . . 55xi4.6 Block diagram of the sliding mode controller and LuGre friction compensatorfor the micromill . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.7 FRF of the af -axis including the nominal FRF and experimental FRFs withand without the notch filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.8 Block diagram of rotary table controller . . . . . . . . . . . . . . . . . . . . . 624.9 Dual stage feed drive controller . . . . . . . . . . . . . . . . . . . . . . . . . . 634.10 Simulated frequency response functions of micromill, rotary table, and mi-cromill with tracking error compensated with rotary table. . . . . . . . . . . . 655.1 Geometric errors of a general axis q . . . . . . . . . . . . . . . . . . . . . . . 685.2 zc-axis positioning errors and the resultant quintic polynomial fit . . . . . . . . 735.3 Tooltip errors with and without compensation for a circle on the x-y plane ofradius 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 796.1 Spiral toolpath . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 816.2 Reference commands for spiral trajectory . . . . . . . . . . . . . . . . . . . . 826.3 Velocity of reference commands for spiral trajectory . . . . . . . . . . . . . . . 826.4 Acceleration of reference commands for spiral trajectory . . . . . . . . . . . . 836.5 Jerk of reference commands for spiral trajectory . . . . . . . . . . . . . . . . . 846.6 Sinusoidal freeform surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . 856.7 Feedrate profile for freeform surface with limit set to 10 [mm/s] . . . . . . . . 866.8 Velocity, acceleration, and jerk profile of the translational axes for the freeformsurface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 876.9 Velocity, acceleration, and jerk profile of the rotational axes for the freeformsurface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88xii6.10 Velocity, acceleration, and jerk profile of the translational axes for the freeformsurface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 896.11 Velocity, acceleration, and jerk profile of the rotational axes for the freeformsurface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 906.12 Machined freeform surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . 916.13 a) Circular contour with sliding mode controller and proposed dual stage feeddrive control scheme. Radius: 30 mm, Tangential Feed Speed: 600 mm/min.Angular frequency of the traverse : ω(Hz) = f/R = 10[Hz]. b) Trackingerror of x and y axis during circular contour . . . . . . . . . . . . . . . . . . . 936.14 a) Square contour with sliding mode controller and proposed dual stage feeddrive control scheme. Tangential Feed Speed: 600 mm/min. b) Tracking errorof x and y axis during square contour . . . . . . . . . . . . . . . . . . . . . . 936.15 Geometric Errors for xc-axis positions -40mm to 40mm . . . . . . . . . . . . . 966.16 Geometric Errors for yc-axis positions -40mm to 40mm . . . . . . . . . . . . . 976.17 Geometric Errors for zc-axis positions -4mm to 76mm . . . . . . . . . . . . . 986.18 Compensating position commands when moving only zc axis . . . . . . . . . . 1006.19 Tooltip errors with and without compensation . . . . . . . . . . . . . . . . . . 101xiiiNomenclatureCc Micromill controllerCf Rotary table controllerEq Error matrix for the q axisFc Averaged Coulomb friction for LuGre friction modelFs Averaged static friction for LuGre friction modelF+/−coul Coulomb friction with forward/reverse velocityF+/−stat Static friction with forward/reverse velocityGc Plant dynamics of micromill feed driveGf Plant dynamics of rotary tableGc,cl Closed loop transfer function of micromill feed driveGf,cl Closed loop transfer function of rotary tableH Cost function constraining the stroke limits of all 9 drives in redundancy resolutionJ Jacobian matrix of the forward kinematic modelJ† Moore-Penrose inverse of the Jacobian matrix JJ c† Moore-Penrose inverse of the Jacobian matrix J cJ c Jacobian matrix of the tooltip positions as a function of rotary table positionsxivJx,y,z,a,b,c Inertia of x, y, z, a, b, and c axis of rotary tableKa Amplifier gainKf Gain of lead-lag controllerKi Gain of integratorKt Motor gainNa Notch filter to remove resonanceNb Notch filter to remove anti-resonanceNi,p B-spline basis function for feed profile with degree pOi Orientation of the tool about the X-axis of the workpiece coordinate frameOj Orientation of the tool about the Y-axis of the workpiece coordinate frameOk Orientation of the tool about the Z-axis of the workpiece coordinate framePi Control points for feed profilePx X position of tooltip in the work piece coordinate framePy Y position of tooltip in the work piece coordinate framePz Z position of tooltip in the work piece coordinate framePx,a X position of tooltip in the work piece coordinate frame considering geometric errorsPx,i Ideal X position of tooltip in the work piece coordinate framePy,a Y position of tooltip in the work piece coordinate frame considering geometric errorsxvPy,i Ideal Y position of tooltip in the work piece coordinate framePz,a Z position of tooltip in the work piece coordinate frame considering geometric errorsPz,i Ideal Z position of tooltip in the work piece coordinate frameS B-spline basis function withSc,cl Sensitivity transfer function of micromill feed driveT Constant to control frequency where phase is added with lead-lag controllerTs Sample time of servo systemVs Input signal∆P Error between the ideal tooltip position and tooltip position considering geometric er-rorsΩ+/−1 Input signalΩ+/−2 Input signalα Constant to control how much phase to add with lead-lag controllerα Time index...q Axis jerks...qmax Axis jerk constraintsq¨ Axis accelerationsq¨max Axis acceleration constraintsxviδx(q) Error in the x-direction for the q axisδy(q) Error in the y-direction for the q axisδz(q) Error in the z-direction for the q axisq˙ Axis velocitiesq˙max Axis velocity constraintss˙(s) Feed profile defined as a function of displacementǫ Measured general geometric errorǫˆ (q) Estimated general geometric error functionλ Bandwidth of sliding mode controllerq Vector of positions of all axesqf Vector of positions of rotary table axesqf Vector of rotary table positionsµ Singularity avoidance constantµ0 Singularity avoidance constant limitωa Resonant frequencyωb Anti-resonant frequencyωc Bandwidth of micromill feed driveωf Bandwidth of rotary tablexviiωm Design cross-over frequency of negative loop transmission of rotary tableρ Disturbance adaption gainσ Sliding surface for sliding mode controllerσ0 Stiffness of microstructures for LuGre friction modelσ1 Damping of microstructures for LuGre friction modelcq cos(q)sq sin(q)εx(q) Error about the x-direction for the q axisεy(q) Error about the y-direction for the q axisεz(q) Error about the z-direction for the q axisϑ Spline parameter of feed profile B-spline functionζa,1/2 Damping ratios to control size and width of Notch filter Naζb,1/2 Damping ratios to control size and width of Notch filter Nbaf Position of the a-axis of the rotary tableb Feed drive viscous frictionbf Position of the b-axis of the rotary tablecf Position of the c-axis of the rotary tablek Time indexxviiik1,2,3,4 Runge Kutta numerical integration constantsm Feed drive massn Displacement indexq General axis positions Tangential displacement along toolpathu Input signalw Manipulability of 9-axis micromachining centerw0 Threshold of manipulability to begin singularity avoidancexc Position of the x-axis of the micromillxf Position of the x-axis of the rotary tableyc Position of the y-axis of the micromillyf Position of the y-axis of the rotary tablezc Position of the z-axis of the micromillzf Position of the z-axis of the rotary tablexixAcknowledgmentsI am most grateful to my research supervisor, Dr. Yusuf Altintas for his support and guidanceover the course of the last 7 years. The insight and opportunities into the world of manufac-turing engineering provided by Dr. Altintas were second to none and through my experienceworking with Dr. Altintas, I have matured as an engineer and an academic researcher in a wayI never thought possible.I would also like to thank all members of the Manufacturing Automation Laboratory family.Whether it was a current member helping overcome research challenges or an alumni providingguidance, members of the MAL family never hesitated to offer a helping hand, and for that, Iam forever grateful.Finally, I would like to thank my mother, father, family, and friends. Without their moralsupport and sacrifices none of this would have been possible. I dedicate this thesis to mygrandfather, Dominic Leung, who lives in my heart.This research was sponsored by the Natural Science and Engineering Research Council ofCanadian Network for Research and Innovation in Machining Technology (NSERC CAN-RIMT).xxChapter 1IntroductionIn recent years there has been a decrease in the average size and increase in the complexityof manufactured parts due to technological progress in industries such as the biomedical andelectronics industry. At present, it is common to see manufactured parts with features in themicrometer scale. Though multiple manufacturing processes with micron level precision doexist, few are capable of producing parts with freeform features with precision. Wire electricaldischarge machining (EDM) and photolithography are capable of very high precisions but canonly be used to generate 2-D and 2.5-D parts, respectively. Though micro injection moldingand micro forming are capable of producing high precision free form parts, a higher precisionmanufacturing process must be used to produce the tooling. Given the limitations of the abovementioned manufacturing processes, multi-axis micromachining is one of the few viable meansfor producing freeform parts with micron level precision. This work contributes a developmentof novel multi-axis micromachine tool, which combines a traditional 3-axis micromill with a 6degrees of freedom (DOF) magnetically actuated rotary table shown in Fig. 1.1.The 6DOF rotary table was originally developed by the Precision Mechatronics Laboratory atthe University of British Columbia and has been presented in [1] and [2]. The rotary tableconsists of a Halbach magnetic array as the mover and a printed circuit board (PCB) as a stator.Actuation is achieved by passing current through the copper traces in the PCB, which createsan opposing magnetic field to the Halbach magnetic arrays on the mover. This configurationallows for a tetherless mover and a low form factor stator. The rotary table has a full 6 degrees1Figure 1.1: UBC MAL’s hybrid 9-axis machineof freedom with approximately 100 micron stroke limit in the x, y, and z directions, 0.5 degreestilt limit about the x and y axis, and unlimited rotation when rotating about the z-axis.By combining a lighter, and subsequently faster, actuator with the existing 3-axis micromill,it is possible to increase the precision and accuracy of the overall system. However, with thisnew configuration, the addition of multiple redundant axes introduce new research challenges.Typically, a 5-axis machine is capable of moving a cutting tool to a predefined cutter locationdefined by 3 positions and 3 orientations. Since a 5-axis machine tool has the necessary de-grees of freedom, analytic inverse kinematic solutions exist, allowing for unique mappings ofmachine tool axes positions to cutter locations. Conversely, the hybrid machine tool presentedin this work has 4 more axes than necessary. As a result, a range of machine tool axes posi-tions can correspond to a predefined cutter location. With respect to trajectory generation, themain challenge is selecting axes configurations that resolve the redundancies. Furthermore, the2trajectory must also avoid singularities and be scheduled in a way to avoid the physical limitsof the drives including the stroke, velocity, acceleration, and jerk limits. This thesis presentsa novel trajectory generation algorithm to avoid singularities while respecting the limits of allaxes.With the trajectory defined, it is necessary to ensure the machine tool follows these trajectoriesas close as possible. This is typically done by applying high bandwidth position controllers toeach axis of the machine tool. Though sophisticated modeling and control techniques can beused to increase the bandwidth to be as high as possible, the achievable bandwidth is limited bythe stability of the closed loop controller of the axis. Likewise, in the presented configuration,the 3-axis micromill is controlled with a high performance non-linear controller, but is limiteddue to stability issues. Due to its comparatively lighter mass and lack of mechanical contact,the rotary table is able to achieve a bandwidth that is an order of magnitude higher than the3-axis micromill. Since the translational axes of the rotary table and micromill run in parallel,it is possible to use the higher bandwidth of the rotary table to compensate the tracking errorscontributed by the low bandwidth of the micromill. The thesis presents a dual stage feed drivecontrol law that combines the strengths of the micromill and rotary table.In addition to increasing the tracking precision of the machine tool, it is possible to improve thevolumetric accuracy of the 3-axis micromill using the 6DOF rotary table. Though a feed driveshould be constrained to only move in the direction of actuation, due to errors in assembly it ispossible for the feed drive to deviate slightly in the orthogonal directions, in the translationaland rotational sense. Furthermore, these errors are typically not detected since the feed driveshave feedback on the direction of actuation. Though the deviations are small for an individualaxis, when multiple axes in series have assembly errors, it can correspond to an end effectordeviation that violates the tolerance of the parts to be machined. In this thesis, the rotary table3is used to compensate for the affect of these position dependent geometric errors. The errorsare measured with a laser interforemeter and mapped to end effector deviations. The trajectoryof the 6DOF rotary table is then modified with a novel algorithm to compensate for thesegeometric errors.Henceforth, the thesis is structured as follows: Chapter 2 discusses previous works reported inliterature specifically with regards to trajectory generation for redundant actuators, dual stagefeed drive control, and geometric error compensation. In Chapter 3 a novel trajectory genera-tion algorithm is presented. The trajectory generation method resolves the redundant degreesof freedom while ensuring that the generated motion commands are jerk continuous, time op-timal, and do not violate the stroke, velocity, acceleration, and jerk limits of the drive. Themodeling, identification and development of the control laws used to control the micromill, ro-tary table, and the combined efforts of the micromill and rotary table are presented in Chapter4. The higher bandwidth of the rotary table is used to compensate for the tracking errors fromthe lower bandwidth of the micromill. Chapter 5 provides an application of the rotary table,where the geometric errors of the micromill are first modeled through experimental measure-ments and compensated by modifying the trajectory of the rotary table. Chapter 6 presentssimulation and experimental results and Chapter 7 concludes the thesis and suggests futurereserach.4Chapter 2Literature Review2.1 OverviewThe configuration in the work presented is a 9-axis micromachining center. Due to its uniqueconfiguration, new trajectory generation algorithms, control laws, and accuracy enhancementtechniques must be developed. With regards to trajectory generation, since there are 4 moredegrees of freedom than necessary the main challenge is in resolving these redundant degreesof freedom. In the presented configuration, the translational axes of the micromill are parallelwith the translational axes of the rotary table. As a result, control laws must be developed thatexploit the strengths of the different actuators, which include the long stroke of the micromillaxis and the high bandwidth of the rotary table axis. Finally, the additional degrees of freedomcan also be used to compensate for geometric errors coming from the 3-axis micromill. Asa result, this literature survey evaluates existing work in trajectory generation of redundantmechanisms, control law design for machine tools and dual-stage servomechanisms, and themodeling, identification, and compensation of geometric errors.2.2 Trajectory Generation for Redundant Manipulators and Machine toolsTypically, a cutting tool only needs 5 degrees of freedom to define its position and orientationin space. Furthermore, it is the end effector position and orientation that is typically specifiedsince it is the element of the overall system that is interacting with its environment. As a result,one of the main objectives of trajectory generation, in the context of redundant manipulators, isresolving the redundant joints or axes so the end effector position and orientation corresponds5to the desired position and orientation. In machine tools, there are typically 5 or less axes, andas a result an analytic solution for this problem is typically available in the form of inversekinematics [3–5]. However, robotic devices typically have more joints than necessary as seenin various types of robotics arms. As a result, the same analytic inverse kinematic solutions cannot be used for these under constrained configurations. In literature, the earliest work for thetrajectory generation of redundant manipulators come from the field of robotics. This problemwas initially solved by Whitney [6], where given a redundant manipulator and specified endeffector velocities, the velocities of the joints are solved using the Moore-Penrose inverse ofthe robot’s Jacobian matrix. Since the Moore-Penrose inverse only considers the least normsolution of the redundancy resolution problem, an infinite range of solution still exists. Leigiosextended this work in [7], where the nullspace of the Jacobian matrix is used to minimize somecost function. By considering the nullspace of the Jacobian matrix, the desired end effectorposition is not affected, but the joint configuration is selected so the predefined cost function isminimized at each step.In Leigios’ original work, the cost function is specified so the solution selects a joint config-uration that always respects the stroke limits of the joints. Numerous authors have extendedLeigios’ original idea and used the cost function to achieve secondary goals such as the opti-mization of joint torques [8], minimization of energy consumption [9], or avoidance of singu-larities [10]. In these solutions, the differential solution minimizes the specified cost functionto a local optimal minimum at each solution step. Nakamura and Hanafusa [11] extended thesetypes of solution to consider global optimality of the cost function. In [11], Nakamura andHanafusa used the example of minimizing manipulability over the whole toolpath. It should benoted that due to the recursive nature of these globally optimal solution, the computation timewas orders of magnitude higher than with the locally optimal solution. It should be noted thatthe redundancy is resolved at the differential level, and as a result, there is a need to numerically6integrate the solution in order to obtain actual motion commands. In [12], the error differencesbetween the various numerical integration techniques were investigated. It was found that the4th Order Runge Kutta Method produced the best results. In order to improve the numericalintegration result, a closed loop corrective element is added in [13]. This way, the numericaldrift inherent in the use of numerical integration techniques is corrected at each step.Though the above methods are able to generate trajectories for redundant actuators, the resul-tant trajectories are not time optimal trajectories that consider the joint velocity, acceleration,and jerk limits. It should be noted that literature exists from the robotics field for optimaltrajectory generation [14] [15], but these works do not address kinematically redundant con-figurations. Since robotic applications do not require the same productivity demands as man-ufacturing applications, the need for time optimality is not as important and it is sufficient fora robot to operate under its maximum capacity. In contrast, trajectory generation algorithmsfor machine tools must be as time optimal as possible, at least acceleration continuous, andrespect the physical limits of all the drives. As a result, trajectory generation techniques inmachine tool literature are focused in fulfilling these criteria in different ways. In [16], Dong etal. proposed a solution to find a time optimal jerk limited trajectory with a bi-directional scanto optimize the trajectory only for jerk and a secondary acceleration-continuation algorithmto optimize the trajectory for acceleration. In [17] a solution is proposed where a tangentialdisplacement is selected at each servo time step in a way that saturates one constraint at alltimes. The solutions proposed in [16] and [17] are only acceleration continuous. In [18] and[5], similar constraints on velocity, acceleration and jerk are used. However, the tangential toolpath is jerk continuous. Due to this added complexity, non-linear optimization techniques wereused, where the feedrate is optimized while respecting the constraints by modifying segmentdurations in [18] and the control points of a spline that defines the feedrate profile in [5].7Based on this literature survey it can be seen that there exists a gap in literature. At present,techniques have not been developed where time optimal, smooth, and constrained trajectoriescan be generated for redundant configurations.2.3 Control laws for Feed drives and Dual Stage ActuatorsDue to the need for high precision positioning, a rich body of literature exist in the field ofsingle axis feed drive control for machine tools beyond the standard industrial Proportional-Integral-Derivative (PID) controllers. Erkorkmaz et al. presented a sophisticated feed drivecontroller in [19] which used a combination of zero phase error tracking control [20], poleplacement, Kalman filter, and feed forward friction compensation. Altintas et al. [21] devel-oped a high bandwidth, robust sliding mode controller with disturbance compensation showingsimilar results as those presented in [19] but with a simpler control structure. Okuwdire ex-tended the work in [22] and included the flexible modes of the ball screw in the design of thesliding mode controller. A switching gain scheduling controller is proposed in [23], whichaccommodates for position dependent dynamics and variations in mass. Hosseinabadi and Al-tintas [24] used an active damping network to damp the structural modes of the machine toolto increase the bandwidth of the sliding mode controller originally proposed in [19]. In [25],Kamalzadeh and Erkorkmaz compensated for the axial vibrations of the ballscrew drive byincluding the first axial mode into the sliding mode control law and showed superior perfor-mance over the use of a notch filter. Kamalzadeh and Erkorkmaz combined multiple controltechniques in [26] and [27], where excitation of torsional modes are avoided with notch fil-ters, control of rigid body dynamics is achieved with an adaptive sliding mode controller, andnon-linear friction and torque ripple are modeled and compensated in a feedforward fashion.It should be noted that the above works represent only a small subset of literature available inthe field of feed drive controls.8In contrast, due to its relative uniqueness, particularly in the field of manufacturing, there isless work for dual stage feed drives. Preliminary work was done in [28] where the trackingerror of the coarse stage is compensated with the fine stage. This work was extended in [29],where the estimated contouring error of the coarse stage is sent to the fine stage instead ofthe tracking error. In these works, a simple PID or PPI controller are used for the feedbackcontrol. Alfizy et al. [30] proposed a dual stage feed drive with a magnetic coarse actuationstage and a piezoelectric driven flexture fine actuation stage. In this work, the coarse stage iscontrolled with a simple PID controller and the fine stage compensates for the tracking error ofthe coarse stage. In [31], Choi et al. proposed a dual servo stage for UV lithography. Unlike theconfigurations presented in [28–30], where the position of the fine actuator is measured relativeto the coarse actuator, the position of the fine actuator is monitored with laser interferometerrelative to the machine tool’s base coordinate frame. Based on the ability to measure thefine actuator with a higher precision sensor relative to the base coordinate frame, the controlrelationship between the coarse and fine actuator is switched, and the coarse actuator followsthe tracking error of the fine actuator. Each axis is then controlled with a PID controller. Itshould be noted that the use of laser interferometer to measure positions in a CNC millingmachine tool would be impractical due to much larger actuating range resulting in difficultiesin mounting mirrors. In contrast to the multi-axis works presented above, in 2004, Kim etal. [32] developed a dual stage actuator and robust controller for camshaft turning. Given therepetitive nature of the input signal, the controller is designed to converge to an optimal designafter a certain amount of time.Beyond the field of multi-axis manufacturing machine tools, dual stage configurations are morecommonly found in 1 DOF read heads of hard disk drives. In these setups, a voice coil motordrives the coarse movement and a piezoelectric actuator drives the fine movement. As a result,the literature for the control of hard disk drive read heads is prevalent with examples of dual9stage control. It should be noted that machine tools and read heads of hard disk drives havedifferent control objectives and configurations. Unlike machine tools, where the axes typicallyfollow a smooth varying trajectory, read heads are typically given step commands to differenttracks located on the disk. As a result the objective is to react as fast as possible to a step com-mand, as opposed to following a smooth trajectory as close as possible. Furthermore, unlikemachine tools, where each individual axis has their own respective feedback sensor, the feed-back sensor detects the position of the read head, the positional sum of the voice coil motor andpiezoelectric actuator, allowing the closed loop control algorithm to control the actual positionof the end effector. As a result, the control laws developed address the control objectives anddynamic analysis in a different way. In [33] the authors implemented a linear controller forthe dual-stage system with the use of a zero phase error tracking controller and a feed forwardcompensator. Kobayashi et al. proposed a phase-stabilized servo controller in [34] where thestructural modes of the piezoelectric actuator are used to compensate for errors from windagedisturbance caused by the suspension vibration. Like the setup presented in this thesis, the finepiezoelectric actuator typically has stroke limits an order of magnitude lower than the coarsevoice coil motor actuator. In [35], Herrmann et al. addressed the problem of fine actuator sat-uration by implementing an anti-windup scheme. Other authors have implemented state spacebased controllers for dual-stage feed drives. In [36] and [37], basic implementations of H∞ andsliding mode controllers, respectively, were applied to the dual-stage configuration. She et al.[38], applied the equivalent input disturbance approach to a dual stage feed drive configuration.In [39] the coarse actuator was driven with an Adaptive Proximate Time Servomechanism andthe fine actuator was controlled by a Composite Nonlinear Feedback controller to reduce thesettling time. In this work the stiffness of the controller of the fine actuator was varied depend-ing on the read head’s proximity to the target. When the read head approaches the target thestiffness of the controller of the fine actuator increases, therefore increasing tracking precision.10Based on this literature survey it can be seen that many control laws exist for single stage feeddrives but few techniques exist for configurations with multi-axis dual stage configurations.The dual stage configurations and subsequently control laws that are found in manufacturingliterature are relatively limited, with no application of more sophisticated control algorithmson both the coarse and fine stages of actuation. Furthermore, the unique configurations ofthe works in multi-axis dual stage configurations do not necessarily apply to the configurationpresented in this work.2.4 Geometric Error Modeling, Identification, and CompensationThe modeling and identification of geometric errors began in the literature for the control ofcoordinate measuring machines (CMM). Since there are no process forces, software compen-sation of geometric errors in CMMs was a viable option. In one of the earliest examples ofsoftware compensation, Zhang et al. [40] modified feedback signals with a look up table inorder to compensate for geometric errors. Similar techniques have been extrapolated to 3-axisCNC machine tools [41] where homogenous transformation matrices are used to map the axesgeometric errors to the tooltip position errors. In [42] a 3D probe-ball and spherical test aredeveloped to measure the link errors in the rotary table of a 5-axis CNC machine tool. The tooldeviations caused by the axes geometric errors are then compensated in real-time by multiply-ing the inverse Jacobian matrix of the forward kinematic model with the tool deviation to obtaina compensating trajectory. Huang et al. also takes a similar approach in [43]. Alternatively, in[44], a method is proposed in which the trajectories for the A and C axis are modified to com-pensate for orientation geometric errors, then the X, Y, and Z axis are modified to compensatefor position geometric error, therefore avoiding the recursive nature of simultaneously solvingfor compensation commands on all axes. However this assumes that moving the translationalaxes will not incur any orientation errors, which may not be necessarily true. In [45] the geo-11metric errors are measured with ball-bar tests. The tooltip errors are mapped into the errors inthe axes positions then added back to the axes positions to compensate for the errors. In [46],the geometric errors are first measured with a ball-bar. The positional deviations are correctedusing the geometric error tables built into the CNC machine. However, positional deviationscaused by rotational geometric errors are corrected by modifying the NC code since these er-rors are not available for compensation in a standard geometric error table in CNC machines.Xiang et al. [47] modeled the geometric errors of a 5-axis CNC machine tool using screw the-ory and then compensated them using inverse kinematics. At present, there is ongoing researchon the geometric error compensation in 5-axis machine tools.As it can be seen, in the above mentioned strategies, the geometric errors are compensated bymodifying the toolpath with respect to the modeled geometric errors. However, comparativelyless work has been done on the compensation of geometric errors with a secondary actuator.This may be a desirable configuration as there may be scenarios in which the trajectory ofthe major actuator cannot be modified. In literature there are a few examples that attempt totackle this problem. A secondary magnetic bearing module has been used to compensate forthe straightness errors of a single axis [48]. In another work, a two dimensional PZT actuatoris used to compensate for the straightness error and positioning errors by sending the PZTactuator the inverted geometric error measured by a laser interferometer [49]. As it can beseen, the existing geometric error compensation techniques with the use of an external deviceare limited to a single axis and limited in general. Furthermore the standard geometric errorcompensation techniques mainly deal with configurations with analytic inverse kinematics,which allow for an analytic inverse differential solution. Since the rotary table has 6 degrees offreedom to compensate for the 3 translational errors, analytic inverse kinematic equations donot exist, and as a result, new techniques must be developed.12Chapter 3Trajectory Generation for a 9-axis Micromachining Center3.1 OverviewIn robotics literature, the trajectory generation techniques are capable of resolving redundantdegrees of freedom, but at present no works are able to do this in an optimal way, wherethe travel time is minimized while simultaneously respecting the velocity, acceleration, andjerk limits of the drives. Conversely, machine tool literature consists of trajectory generationtechniques where the travel time is optimal while respecting the physical constraints of alldrives involved, but have not been developed to resolve redundant degrees of freedom. Inthis chapter, a trajectory generation technique is presented where the strength of trajectorygeneration in robotics and machine tool literature are combined.Given a computer aided design (CAD) geometry of the part to be manufactured, the computeraided manufacturing (CAM) system will generate a corresponding toolpath for the tool to fol-low, which will typically consist of desired tooltip positions and tool orientations at varyingdisplacement intervals. In the case of the configuration presented, the trajectory generationalgorithm must first decompose the desired tooltip positions and orientation into 9-axis com-mands which correspond to the desired tooltip position and orientation. Next it must schedulethese commands in such a way that the velocity, acceleration, and jerk limits of all the axes arenot violated while traversing the toolpath as fast as possible.In order to do this, a forward kinematic model is first developed. Based on this kinematicmodel and the desired tooltip position and tool orientations, the redundancies are resolved13numerically and corresponding axes positions are found at fixed displacement intervals alongthe toolpath. These axes positions are then scheduled in time with a feedrate optimizer, to beas fast as possible without violating the limits of the axes. Finally the feed profile is resampled,resulting in real-time position commands to be sent to the drives of the machine tool.3.2 Forward Kinematics of the 9-axis Micromachining CenterGiven the position commands of each of the drives it is possible to find the tooltip position andtool orientation using the forward kinematic model of the machine tool. In order to generatethe forward kinematic model of a machine tool, first, two kinematic chains are formed withhomogenous transformation matrices (HTM) based on the sequence of the drives [50]. In theconfiguration presented, the workpiece kinematic chain starts from an arbitrary base coordinateframe and ends at the workpiece, BTw, and the tool kinematic chain starts from the samearbitrary base frame and ends at the tool holder, BTt, as shown in Fig. 3.1.With the kinematic chains defined, the inverse of the workpiece kinematic chain is multipliedwith the tool kinematic chain to obtain the final transformation matrix, wTt, which transformspositions defined in the tool coordinate into the workpiece coordinate frame. The first kine-matic chain consists of the x-y table of the micromill and all the axes of the rotary table. As aresult, based on this sequence, the transformation matrix to transform from the arbitrary base14Lsp,yLsp,zLtzc,0wtBTtLrzf,0 LsLy,zLx,zLwBTwzyBFigure 3.1: 9-axis micromachining center and corresponding axes of the machine tool15coordinate system to the workpiece coordinate frame is as follows:BTw,i =1 0 0 00 1 0 yc0 0 1 Ly,z0 0 0 1︸ ︷︷ ︸BTyc,i1 0 0 xc0 1 0 00 0 1 Lx,z0 0 0 1︸ ︷︷ ︸ycTxc,i1 0 0 00 1 0 00 0 1 Ls0 0 0 1︸ ︷︷ ︸xcTs,i×1 0 0 xf0 1 0 00 0 1 00 0 0 1︸ ︷︷ ︸sTxf ,i1 0 0 00 1 0 yf0 0 1 00 0 0 1︸ ︷︷ ︸xfTyf ,i1 0 0 00 1 0 00 0 1 zf + zf00 0 0 1︸ ︷︷ ︸yfTzf ,i×1 0 0 00 caf −saf 00 saf caf 00 0 0 1︸ ︷︷ ︸zfTaf ,icbf 0 sbf 00 1 0 0−sbf 0 cbf 00 0 0 1︸ ︷︷ ︸afTbf ,i×ccf −scf 0 0scf ccf 0 00 0 1 00 0 0 1︸ ︷︷ ︸bfTcf ,i1 0 0 00 1 0 00 0 1 Lr0 0 0 1︸ ︷︷ ︸cfTr,i1 0 0 00 1 0 00 0 1 Lw0 0 0 1︸ ︷︷ ︸rTw,i(3.1)where cq = cos(q) and sq = sin(q), q is a general axis position, and aTb,i denotes a transfor-mation matrix from coordinate frame a to coordinate frame b and the subscript i indicates anideal homogenous transformation. It should be noted that in Eq.(3.1), the coordinate frames s,r, and w, denote the stator, rotor, and workpiece coordinate frames, respectively. Furthermore,xc and yc are the position commands of the micromill, xf , yf , zf , af , bf , and cf are the positioncommands of the rotary table, Lx,z is the height of the xc-axis, Ly,z is the height of the yc-axis,Ls is the height of the stator, Lr is the height of the rotor, Lw is the height of the workpiece,and zf,0 is the initial floating distance of the rotary table. All the offsets are static and shown inFig. 3.1.16The second kinematic chain consists of the z-axis the micromill and the spindle offsets. As aresult, the tool kinematic chain is defined as:BTt,i =1 0 0 00 1 0 00 0 1 zc + L0,z0 0 0 1︸ ︷︷ ︸BTzc,i1 0 0 00 1 0 Lsp,y0 0 1 −Lsp,z0 0 0 1︸ ︷︷ ︸zcTt,i(3.2)where zc is the z-axis position command of the micromill, Lsp,y and Lsp,z are the linear offsetof the spindle, and L0,z is the initial offset of the z-axis of the micromill from the arbitrary basecoordinate frame. The initial offset, L0,z, can be defined as:L0,z = Lx,z + Ly,z + Ls + zf,0 + Lr + Lw + zc,0 + Lsp,z + Lt (3.3)where zc,0 is the initial position of the z-axis of the micromill and Lt is the length of the cuttingtool. With both transformation matrices originating at the arbitrary base coordinate frame, thetransformation matrix from the tool coordinate frame to the workpiece coordinate frame isfound by multiplying the inverse of Eq. (3.1) with Eq. (3.2) as follows:wTt,i =(BTw,i)−1 BTt,i (3.4)If the position of the tooltip is defined in the tool coordinate frame, then it is possible to findits position relative to the workpiece with the transformation defined in Eq. (3.4). Relative tothe tool coordinate frame, the tooltip position is the exposed length of the cutting tool and theorientation is always lined up with the z-axis of the tool coordinate frame. As a result, the toolcan be defined in the tool coordinate frame as ttp = [0 0 Lt]T and tto = [0 0 1]T for the tooltipposition and tool orientation, respectively. Given the transformation matrix in Eq. (3.4), thetooltip position and tool orientation can be transformed into the workpiece coordinate frame as17follows: Px OiPy OjPz Ok1 0 = T[ttptto1 0](3.5)Expanding the matrix equation defined in Eq. (3.5), the forward kinematic model gives thetooltip position and tool orientation defined as a function of the positions of the 9 drives asfollows:Px (q) = (−Lsp,y − yc − yf )(caf scf + ccf saf sbf )+(Lr + Lw + zc,0 + zc − zf )(saf scf − caf ccf sbf ) (3.6)+cbf ccf (−xc − xf )Py (q) = (−Lsp,y − yc − yf )(caf ccf − saf sbf scf )+(Lr + Lw + zc,0 + zc − zf )(ccf saf + caf sbf scf ) (3.7)+(xc + xf ))cbf scfPz (q) = (−xc − xf )sbf − Lwz − Lr+(Lr + Lw + zc,0 + zc − zf )caf cbf (3.8)+(Lsp,y + yc + yf ))cbf safOi(q) = saf scf − caf ccf sbf (3.9)Oj(q) = ccf saf + caf sbf scf (3.10)Ok(q) = caf cbf (3.11)where q = [xc, yc, zc, xf , yf , zf , af , bf , cf ]T , Px, Py, and Pz are the x, y, and z position of thetooltip in the workpiece coordinate frame, and Oi, Oj , and Ok are the orientation of the toolabout the x, y, and z axis of the workpiece coordinate frame. In the next section, the Jacobianof the forward kinematics, as defined in Eq. (3.6) to (3.11), along with the desired toolpath is18used to resolve the redundant degrees of freedom.3.3 Redundancy Resolution of the 9-axis Micromachining CenterTypically, machining toolpaths from the CAM system are specified as the tooltip position andtool orientation with respect to the workpiece, also known as cutter location (CL) data. The CLdata from the CAM system is first fit to a jerk continuous B-spline as a function of tangentialdisplacement , s, to preserve third order continuity [51] in order to avoid exciting the structuralmodes of the machine tool:R(s) = [Px(s), Py(s), Pz(s), Oi(s), Oj(s), Ok(s)]T (3.12)In traditional 3-axis and 5-axis machine tools, analytic inverse kinematic equations can be usedto find the corresponding axes positions given R(s). However, since the 9-axis micromachin-ing center has more axes than necessary, the system is under constrained and analytic inversekinematic solutions do not exist. In this work, numerical techniques from robotics literature [6][7] are adapted to resolve the under constrained system. In order to generate trajectories thatresolve the redundancy of the micromachining center, axes configurations, which correspondto the desired tooltip position and tool orientation, at fixed displacement intervals along thetoolpath are first found. With the forward kinematic equations, it is possible to find the Jaco-bian matrix, which is the differential of the tooltip position and tool orientation with respect toaxes positions:J =dRdq=dPxdxc· · · dPxdcf.........dOkdxc· · · dOkdcf6×9(3.13)19Given differential axes positions, it is possible to find the differential tooltip positions and toolorientations, as follows:dR = Jdq (3.14)However, since the objective is to find q from R(s), the inverse of J is required. Since J is nota square matrix it is not possible to simply multiple dR with the inverse of J . Alternatively, asoriginally proposed in [6], the Moore-Penrose inverse, defined as:J† = JT(JJT)−1 (3.15)can be used to find the differential axes positions. By multiplying the Moore-Penrose inversewith a differential tooltip position and tool orientation, it is possible to find the correspondingdifferential axes position. Since the spline describing the tooltip position and tool orientation,R(s), is a function of tangential displacement, s, it can be differentiated with respect to sand multiplied with J† to get the axes positions differentiated with respect to displacement asfollows:dqds= J†dRds(3.16)It should be noted that an infinite range of differential solutions exist and Eq. (3.16) representsthe least squares solution of the under constrained problem. This corresponds to the solutionwith the lowest average axes differential with respect to displacement at each step or moreformally:minimize∣∣∣∣∣∣∣∣dqds∣∣∣∣∣∣∣∣subject to dR = Jdq(3.17)However, since our trajectory has requirements beyond the least norm solution, Eq. (3.16) isaugmented to fulfill other requirements such as avoiding stroke limits and singularities. Fur-20thermore, the differential solution must be numerically integrated at each step in order to obtainposition commands at fixed displacement intervals.3.4 Avoiding Stroke Limits of the DrivesAs mentioned, Eq. (3.16) represents a single solution in an infinite range of solutions to Eq.(3.14). As a result Eq. (3.16) can be augmented to fulfill other criteria by using the nullspaceof the Jacobian matrix J in the following way:dqds= J†dRds+ β(I − JJ†)∇H (3.18)where(I − JJ†) is the nullspace of the Jacobian, β is a constant gain, and H is a cost functionto be minimized. Any projection onto the nullspace does not effect the position of the endeffector, so joint configurations can be selected that minimize the cost function defined in Hwithout affecting the end effector’s position. Due to the highly varying stroke limits of each ofthe two systems, it is important to ensure that the stroke limits of each of the drives are respectedwhen the trajectories are generated. In order to constraint the generated axes commands, thecost function is defined in a way that position commands close to the stroke limits of the drivesare penalized:H(q) =x2c(x2c,max − x2c)2 + · · ·+ c2f(c2f,max − c2f)2 (3.19)where qmax is the stroke limit of the respective axis. As it can be seen, the cost functionapproaches infinity as the axes positions approach the stroke limits. Equation (3.18) projectsthe gradient of H , which is defined as:∇H =[2x2c,maxxc(x2c,max − x2c)2 , . . . , 2c2f,maxcf(c2f,max − c2f)2](3.20)21on to the nullspace Jacobian, thereby selecting an axes configuration that gives the desiredtooltip position and tool orientation, while respecting the stroke limits of the drives. By substi-tuting Eq. (3.20) into Eq. (3.18), the solution works similar to a gradient descent optimization,where the solution selects the configuration that minimizes the cost function defined by H . Itshould be noted the solution divides its efforts between attempting to minimize the norm of theaxes differentials for all the axes and minimizing the cost function H .To show this, a trajectory for a circle of radius 4 [mm] on the x-y plane is generated using onlythe x-y axis of the micromill and rotary table. Figures 3.2 and 3.3 shows trajectories generatedby the algorithm without and with the stroke constraints, respectively. Figure 3.4 compares thenorm of the reference commands for the constrained and unconstrained configurations. As itcan be seen, without the constraint, as shown in Fig. 3.2, the fine and coarse actuator have thesame motion profile, resulting in the lowest average axes movements as shown in the norm ofthe motion commands in Fig. 3.4. However this would not be possible as the fine actuatorsonly have a stroke limit of 100 [µm]. When the constraints are included, the fine actuator isconstrained to its stroke limits while the coarse actuator takes up the remainder of the trajectory,as shown in Fig. 3.3, resulting in higher average axes movements as shown in Fig. 3.4.The strength of each portion of the solution is dictated by the gain β. A larger β results in asolution that prefers to keep the position commands away from the stroke limits but may notnecessarily have the lowest average axes movement. Conversely, a small β reduces the averageaxes movement but will have configurations that are closer to the stroke limits. It should benoted that if β is too small, Eq. (3.18) may result in solutions that violate the stroke limits.220 5 10 15 20 25-505Position Commands [mm]0 5 10 15 20 25Displacement [mm]-505Position Commands [mm]xc xf10xfstroke limits 10yfstroke limits ycyfFigure 3.2: Generated trajectories without constraint cost function0 5 10 15 20 25-505Position Commands [mm]0 5 10 15 20 25Displacement [mm]-505Position Commands [mm]10xfstroke limitsxc10yfyc10yfstroke limits10xfFigure 3.3: Generated trajectories with constraint cost function230 5 10 15 20 25Displacement [mm]012345678|q||q| without constraints|q| with constraintsFigure 3.4: The norms of dqdswith and without constraint cost function243.5 Singularity Avoidance with augmentation of Jacobian MatrixDue to the non-Cartesian movements of the 9-axis micromachining center, it is possible toend up in a singular configuration where the machine tool losses actuation capabilities in adirection. Mathematically, it results in the Jacobian matrix becoming degenerate, and as aresult, the Moore-Penrose inverse cannot be found and Eq. (3.18) will not have a solution. Inorder to avoid this scenario, the Jacobian matrix can be augmented by adding a diagonal matrixof small constants [10] so that a solution can always be found as follows:J∗ = JT(JJT + µI)−1 (3.21)The trade off for guaranteeing a solution is the introduction of an error into the solution, sincethe Moore-Penrose inverse has been modified. However, an augmented solution is only nec-essary in the vicinity of a singularity. As a result, the Moore-Penrose inverse is conditionallymodified by varying the magnitude of the small constant µ in the following way:µ =µ0(1− ww0)2, if w < w00, otherwise(3.22)where w =√det(JJT ) is manipulability, which corresponds to how close the machine isto a singular configuration. w0 is the threshold for manipulability in which the modificationconstant µ starts taking on a non-zero value. Furthermore, µ0 corresponds to the maximumamount of deviation the solution is allowed to have from the unmodified solution. If w0 andµ0, are set too low the solution may still be singular, resulting in a degenerate solution. Incontrast, if w0 and µ0 are too large, then the deviation from the desired path may be too large,violating user defined tolerances. The different effects of different sizes of w0 and µ0 can beseen in Fig. 3.5 and Fig. 3.6, where w0 = µ0 = 1 in Fig. 3.5 and w0 = µ0 = 10−6 in Fig. 3.6.As a result w0 and µ0 are tuning parameters which depend on the accuracy constraints.25-0--40 0--40-600Y -60Z-80 X -80 -100Reference ToolpathActual ToolpathFigure 3.5: Desired toolpath and generated tool path with w0 = µ0 = 1-2000-2020-40 0-20-40-600-60Z-axis [mm]-80-80 -100Y-axis [mm]X-axis [mm]200Reference ToolpathActual ToolpathFigure 3.6: Desired toolpath and generated tool path with w0 = µ0 = 10−6263.6 Numerical Integration of Differential Solution to get Position CommandsThough the solution is capable of avoiding the stroke limit of drives and singularities, thesolution only outputs a differential, dqds. In order to obtain axes positions, the solution must beintegrated over displacement, s, as follows:q(s) =∫ s0dqdsds (3.23)However, the analytic integration dqdsis not straightforward and the integration is performednumerically instead. The length of the whole toolpath, S, is first found with Simpson’s method[52], then divided into N intervals giving us the fixed displacement intervals ∆s = S/N .Through numerical integration, the axes positions q[n] are found at the displacements s[n] =∆s × n where n = 1, 2, 3, . . . , N − 1, N . Several works [12] have analyzed various ways toperform the numerical integration and as expected, the 4th Order Runge Kutta method providedthe best results. Since the initial joint configuration is known, and it is known that there is nodisplacement at the beginning of the tool path, the initial values can be set as q[0] = q0 ands[0] = 0. With these initial values, the numerical integration can be performed in the followingway:q[n+ 1] = q[n] + (k1 + 2k2 + 2k3 + k4)∆s6k1 =dq (s[n],q[n])dsk2 =dq(s[n] +∆s2,q[n] +∆s2k1)dsk3 =dq(s[n] +∆s2,q[n] +∆s2k2)dsk4 =dq (s[n] + ∆s,q[n] + ∆sk3)ds(3.24)which gives N axes positions, q[n], that correspond to the N desired tool position and orien-tation R(s[n]) at displacements s[n]. Unlike the majority of works in robotics literature that27employ a similar technique, the commands are generated with respect to displacement, allow-ing for further optimization over time. The 4th Order Runge Kutta integration is repeated untilthe end of the toolpath is reached, resulting in known axes positions at fixed displacementintervals.3.7 Closed loop correction of deviations from numerical integrationDue to the numerical nature of the integration algorithm, a small amount of numerical driftoccurs at each step, causing the corresponding tooltip position and tool orientation from theaxes positions q[n] to deviate from the desired tooltip position and tool orientation specified byR(s). As a result, a corrective action is introduced into the differential solution in Eq. (3.18)[13]. Given the solved axes position q[n] at a given displacement s[n], it is possible to find theresultant tooltip position and tool orientation with the forward kinematics as described in Eq.(3.6) to (3.11). For brevity, the forward kinematic equations are grouped into a single vectorfunction:f(q[n]) =Px (q[n])Py (q[n])Pz (q[n])Oi (q[n])Oj (q[n])Ok (q[n])(3.25)As result, the error between the desired tooltip position and orientation and the result of thedifferential solution can be found as follows:e[n] = R(s[n])− f(q[n]) (3.26)Given the error e[n], it is possible to include closed loop corrective action into Eq. (3.18) soaxes configurations can be selected to correct the numerical deviation caused by the numerical28integration algorithm as follows:dqds= J∗(Kee[n] +dRds)+ β(I− J∗J)∇H (3.27)where Ke is a scalar gain. As a result, similar to a controller, the differential solution seeks tominimize the error at each subsequent step. It should be noted that, the size of the scalar gainis limited by the stability of the solution. Should too high of a gain be set, the solution willbecome unstable. To demonstrate its effectiveness, Fig. 3.8 and Fig. 3.7 show the errors ofthe Runge Kutta numerical integration algorithm, with and without the closed loop correctionimplemented, respectively.From Fig.3.8 it should be noted that the 4th Order Runge Kutta integration method is able tokeep the numerical errors in the range of 10−11 [mm] and 10−10 [rads] which is well below thetolerance in most machining applications. However, it is important to note that since there isno corrective action for numerical drift, as the toolpath gets longer, numerical errors may beincurred at each step causing violation at some point on the toolpath. From Fig. 3.8, the closedloop corrective action corrects for the numerical drift resulting in a numerical error of 10−13[mm] and 10−12 [rads], at least an order of magnitude lower than without the corrective action.Furthermore, error is not incurred at each step, ensuring that the algorithm will not violate thetolerance for longer toolpaths.3.8 Full Redundancy Resolution AlgorithmThe full redundancy resolution algorithm incorporates numerous elements, including cost func-tions, singularity avoidance, numerical integration, and closed loop corrective action. Forbrevity:dqds(s[n],q[n]) = g (s[n],q[n]) (3.28)290 20 40 60 80 100 120 140-101NPfffiflffi 10!"#0 2$ 40 60 80 100 12$ 140D%&'()*+,+./ 0,,34202NOflfififlffi5ff 10!"#6y6x6z7k7i7jFigure 3.7: Numerical errors of the redundancy resolution algorithm with closed loop correc-tion0 20 40 60 80 100 120 140-10189:;<=>?@ABC=E=BFG<**?@`<=;FE?E=BFG<**** 0F−coul if v < 0(4.8)which must be modeled as an equivalent disturbance and included into the difference equa-tion as d(k). The equivalent disturbances, corresponding to the Coulomb friction, is found asfollows:d+f =F+coulKaKt(4.9)d−f =F−coulKaKt(4.10)Based on the friction model in Eq.(4.8) and the equivalent disturbances in Eq.(4.9) and Eq.(4.10),the disturbance from the Coulomb friction is defined as follows:df (v(k)) = PV (v(k)) · d+f +NV (v(k)) · d−f (4.11)40where PV (v(k)) and NV (v(k)) are the positive and negative velocity functions defined asfollows:PV (v(k)) =12· σ (v(k)) · (1 + σ (v(k))) (4.12)NV (v(k)) = −12· σ (v(k)) · (1− σ (v(k))) (4.13)and where σ is a deadband function to ensure encoder noise does not affect the identification:σ (v, Vd) =0 if |v| ≤ Vd1 if v > Vd−1 if v < −Vd(4.14)where Vd is the limit in which the measured velocity is considered to be non-zero. This ensuresthat noise from the velocity measurement is not mistaken for actual movement, which maypotentially affect the identification results. With the equivalent disturbance defined in Eq.(4.11), the difference equation can be rewritten as:v(k) = pvdv(k − 1) +Kvdu(k − 1)− [Kvdd+f PV (v(k − 1)) +Kvdd−f NV (v(k − 1))] (4.15)A frequency rich signal, u, as shown in Fig. 4.1 is sent to the machine and the velocity ismeasured. This type of signal is selected in order to excite the system at as many frequenciesas possible, resulting in better fitting results. Since the velocity, v, is measured by the encoder,it also defines the speed dependent sign functions, PV (v) and NV (v). The only unknowns inEq. (4.15) are the coefficients pvd, Kvd, d+f , and d−f . The problem can be redefined as a matrix410 lmn lmo lmp lmq 1 rmn rmostuv wx{|n-101nControl Signal [V]Figure 4.1: Frequency rich signal for rigid body identificationproblem as follows:v(2)v(3)...v(N)︸ ︷︷ ︸Y=v(1) u(1) PV (1) NV (1)v(2) u(2) PV (2) NV (2)............v(N − 1) u(N − 1) PV (N − 1) NV (N − 1)︸ ︷︷ ︸ΦpvdKvdKvdd+fKvdd−f︸ ︷︷ ︸θ(4.16)whereN is the total number of measurements in the identification experiments. The objective isto find coefficients in θ that model the relationship as accurately as possible. Since Y is known,and Φ is known, the coefficients in θ are found in order to minimize the error between thepredicted and measured value, which can be phrased in the following minimization problem:min12(Y − Φθ)T (Y − Φθ) (4.17)42Table 4.1: Identified rigid body parameters of the micromillParameter x y zmˆ [kg] 21.19 59.17 17.61bˆ [kg/s] 325.91 346.01 179.10Fˆ+coul [N] 6.29 19.98 17.47Fˆ−coul [N] -6.094 -18.93 -18.28which can be minimized by setting θ as follows:θˆ = (ΦTΦ)−1ΦTY (4.18)From θ the estimated values for mass (mˆ), viscous friction (bˆ), and coulomb friction (Fˆ+/−coul )can be found as follows:mˆ =(pˆvd − 1)KtKaTsKˆvd ln(pˆvd)(4.19)bˆ =(1− pˆvd)KtKaKˆvd(4.20)Fˆ+coul =ˆKvdd+fKˆvd(4.21)Fˆ−coul =ˆKvdd−fKˆvd(4.22)From these tests, the corresponding mass, viscous friction, and coulomb friction values foreach axis are found and outlined in Tab. 4.1. Based upon these results the friction model isimproved as outlined in the following section.4.2.2 Identification of Friction CharacteristicsIn addition to the Coulomb friction modeled and identified in Section 4.2.1, feed drives haveadditional non-linear friction phenomenon, which must be identified and compensated in order43to improve tracking performance. To identify this non-linear friction characteristic, a Kalmanfilter is designed to observe the disturbance when the axis is jogged back and forth at vari-ous speeds. In the absence of process forces, the disturbances observed by the Kalman filtercorresponds to non-linear friction phenomena. A simple pole placement controller is first de-signed, based on the rigid body parameters found in the previous section, allowing the axis tobe jogged at different speeds. It should be noted that the pole placement controller is used foridentification purposes, and a more sophisticated sliding mode controller is developed in thefollowing section.In order to design the Kalman filter, the feed drive must be modeled in discrete state spaceform, with disturbance as one of the states while considering the measurement and input noise.With the axis velocity defined in Eq. (4.3), the axis position can be found as follows:x(s) =1sKvs− pv [u(s)− d(s)] (4.23)which gives the following state space model:[x˙(t)v˙(t)]= Ac[x(t)v(t)]+Bc[u(t)d(t)](4.24)where Ac and Bc are the transition and input matrices defined as:Ac =[0 10 pv](4.25)Bc =[0 0Kv −Kv](4.26)From this state space model, the discrete time state is found:[x(k + 1)v(k + 1)]= Ad[x(k)v(k)]+Bd[u(k)d(k)](4.27)44where Ad and Bd are the discrete transition and input matrices defined as:Ad = eAcTs (4.28)Bd =∫ Ts0eAcλdλ ·Bc (4.29)In order to design a Kalman filter that can estimate disturbances in the presence of noise, thediscrete state space model in Eq.(4.27) is augmented in the following way: x(k + 1)v(k + 1)d(k + 1) = A x(k)v(k)d(k)+B [ u(k) ]+W[u˜(k)wd(k)](4.30)[xm(k)vm(k)]= C x(k)v(k)d(k)+V[x˜(k)v˜(k)](4.31)where u˜(k) is the input noise, x˜(k) is the feedback position noise, v˜(k) is the velocity feedbacknoise, and wd(k) is the disturbance noise. The matrices in the augmented state space model aredefined as:A =[Ad −Bd0 0 1]B =[Bd0]C =[1 0 00 1 0]W = Bd 000 1 V =[1 00 1](4.32)The augmented state space model treats the disturbance as a state, allowing the Kalman filterto observe it. Next, the noise of the feedback and input must be considered in order to selectoptimal gains for the Kalman filter. The position of the feed drives are measured with a linearencoder with a resolution of δx = 40 [nm], which introduces an error in measurement dueto quantization. As a result, the error between the actual and measured position is uniformly45distributed with a zero mean and values between +δx/2 and −δx/2. The variance can becomputed as:Rx˜ = E[(x˜− E[x˜]) = E[x˜2] =∫ ∞−∞p (x˜) · x˜2dx˜ =∫ δx/2−δx/2x˜2dx˜ =(δx)212(4.33)Since the noise is the digitally differentiated encoder signal, then the resolution is δv = 40/Ts[nm/s] resulting in a variance of:Rv˜ =(δv)212=(δx/Ts)212(4.34)Similar to the feedback, quantization of the input signal can introduce additional noise intothe system. The current amplifier is controlled with a PWM control signal with a frequency of50kHz. The real-time controller of dSpace is capable of varying the the PWM with a resolutionof 50ns. As a result, within a duty cycle, it is possible to have 400 unique values. Sincethe amplifier has a saturation limit of ±5[V ] then the resolution of the input signal is δu =10/400[V] = 0.025[V], which, similar to the feedback, gives a variance of:Ru˜ =(δu)212(4.35)Unlike the covariances Rx˜, Rv˜, and Ru˜, Rwd is a tuning parameter that adjusts the performanceof the Kalman filter with respect to predicting disturbance. A higher covariance results in aKalman filter that converges faster but outputs noisier results and vice versa for lower covari-ance. It was found the following values resulted in the best balance between performance andnoise:Rwd = 1× 10−6 (4.36)46Given these parameters, the Kalman filter takes the form: xˆ(k)vˆ(k)dˆ(k) = (I−KobsC)A xˆ(k − 1)vˆ(k − 1)dˆ(k − 1) = (I−KobsC)B [ u(k − 1) ]+Kobs[xm(k)vm(k)](4.37)where xˆ(k), vˆ(k), and dˆ(k) are the estimated position, velocity, and disturbance, respectively,at time step k, Kobs is the Kalman filter gain matrix and xm(k) and vm(k) are the measured axisposition and velocity, respectively. Based on these system parameters and noise variances, theKalman observer gain matrix, Kobs, is found [55] [19], giving us the following observer gainsfor the x, y, and z axis as follows:Kobs,x = 0.2227 2.8518E− 6285.1827 0.0061−73.5388 −0.0020 (4.38)Kobs,y = 0.1880 1.9796E− 6197.9680 0.0034−76.1335 −0.0017 (4.39)Kobs,z = 0.1993 2.2430E− 6224.3086 0.0042−75.3177 −0.0018 (4.40)With these Kalman filters, the axes are jogged back and forth at ±1 [mm/s], ±3 [mm/s], ±5[mm/s],±10 [mm/s],±20 [mm/s],±30 [mm/s],±40 [mm/s],±50 [mm/s],±75 [mm/s],±100[mm/s],±125 [mm/s], and±150 [mm/s] and the disturbance is measured over the whole strokeof the feed drive. An example of the disturbance prediction at 10 [mm/s] can be seen in Fig. 4.2.The disturbance over the whole stroke of the feed drive is averaged, resulting in a disturbanceestimate at each corresponding velocity as shown in Fig. 4.3.47-40 -30 -20 -10 0 10 20 30 40}~
Ł-40-2002040Disturbance Estimate [V]Positive VelocityNegative VelocityFigure 4.2: Disturbance estimation for xc-axis at 10 [mm/s] -100 0 100 -50050Disturbance Estimate [V] ¡¢£¡¤ ¥ ¦¢£¦§¨©ª¡ «¡¤ ¥ ¦¢£¦§Figure 4.3: Average disturbance estimation for xc-axis at all test speeds48From the velocity based disturbance data a Stribeck friction model of the form:Ff (v) =F+state−v/Ω+1 + F+coul(1− e−v/Ω+2)+ bv, for v > 0F−state−v/Ω−1 + F−coul(1− e−v/Ω−2)+ bv, for v < 0(4.41)is fit. In Eq.(4.41) F+stat and F−stat are the static frictions in the positive and negative directionsrespectively, F+coul and F−coul are the Coulomb frictions in the positive and negative directionsrespectively, and Ω+1 , Ω−1 , Ω+2 , and Ω−2 , are exponential coefficients which determine the ratein which the static friction converges to the Coulomb friction.Since the Kalman observer is based on the parameters identified in the rigid body model thedisturbance should consist of only the non-linear friction phenomenon. However, it is possiblethat the viscous friction, b, from the openloop rigid body identification may be inaccurate.Should this be the case, the disturbance will exhibit linear variation at the higher speeds, wherethe disturbance should only be defined by a constant Coulomb friction. An example of thiscan be seen in Fig. 4.3, where the observed friction has a negative slope at higher velocities,indicating a viscous friction that is higher than one obtained from the identification. If thisis the case, the originally identified viscous friction coefficient is adjusted by first identifyingthe slopes in the disturbance estimate in the positive,△b+, and negative direction,△b−, usinglinear regression. An average slope is found:△b¯ = △b+ +△b−2(4.42)and the original viscous friction is updated as follows:bˆ′ = bˆ+△b¯ (4.43)From Fig. 4.3, it can be seen that the slopes at higher velocities are removed once the viscousfriction has been corrected. Once corrected, the viscous friction in the Kalman filter is updated49Table 4.2: Identified friction parameters of the micromillParameter x y zb [kg/s] 25.05 20.55 102.9F+stat [N] 36.38 43.2902 30.6853F−stat [N] -29.03 -29.2014 -33.9802F+coul [N] 39.74 47.9513 23.6214F−coul [N] -26.74 -35.5574 -23.5320Ω+1 [mm/s] 1.31 8.81 6.2Ω−1 [mm/s] -0.34 -0.74 -6.8Ω+2 [mm/s] 1.79 10 4.2Ω−2 [mm/s] -0.45 -0.74 -4.1and the velocity tests are ran so the non-linear friction parameters can be identified more accu-rately. The disturbance with the corrected viscous friction can be seen in Fig. 4.4. In order tofit the parameters in the Stribeck model, F+stat and F−stat are first selected as the disturbance forceat the lowest velocities in the positive and negative direction, respectively. Since the non-linearfriction affects associated with static friction go to zero at high velocities, the Coulomb fric-tion, F+coul and F−coul, is selected as the average disturbance force of the four highest velocities.Finally, the exponential coefficients are found by minimizing the cost function:V =12N2∑k=N1[δF (k)]2 (4.44)where δF (k) is the difference between the measured disturbance and predicted disturbance:δF (k) = δF [v(k)] = F ′f [v(k)]− Fˆ+/−stat e−v/Ω+/−1 + Fˆ+/−coul(1− e−v/Ω+/−2)(4.45)Due to the non-linear and discrete nature of Eq. (4.45) the cost function is minimized in a bruteforce way. All permutations between 0 and 10 for are tested at increments of 0.1 [mm/s]. Basedupon this friction identification methodology, the parameters shown in Tab. 4.2 are found.50-100 ¬® 0 ® 100¯°±²³´µ¶ ·¸¸¹º»-50050Disturbance Estimate [V]¼½¾¿ÀÁÂÃÁÄ ¼ÂÅÆÇÅÈÉ ÊÅÇËÌÈÂÂÁÆÇÁÄ ÍÅÀÆÈÎÀ ¼ÂÅÆÇÅÈÉFigure 4.4: Disturbance estimation for xc-axis with corrected viscous friction and the fittedfriction curveWith the above coefficients, the friction behaviour, as described by the model, can be seen inFig. 4.4. Using the identified mass, corrected viscous friction, and non-linear friction param-eters, a sliding mode controller and a feedfoward friction compensator, based on the LuGremodel, is developed.4.2.3 Sliding Mode Position ControllerUnlike conventional ballscrews, where the pitch of the ball screw reduces the reflected distur-bance on the motor, a linear feed drive experiences a disturbance force completely. As a result,a controller must be selected which has good disturbance rejection properties. Furthermore,specific to the micromill in this work, the friction properties have significant variation over thelength of the travel. As a result, the sliding mode controller [56] was selected as the positioncontroller for the feed drives of the micromill. The design of the controller is based on work51presented in [21] but has been extended to linear feed drives. From Eq. (4.1) the correspondingtime domain state space model can be found:[x˙x¨]=[0 10 −be/me]·[xx˙]+[01/me]· u−[01/me]· d (4.46)where me = m/(KaKt) and be = b/(KaKt) and the indication for function of time (t) isdropped for brevity. Since the goal of the feed drive is to follow a trajectory, a sliding surfaceis designed to bring the tracking error to zero:σ = (x˙r − x˙) + λ (xr − x) = [λ 1]︸︷︷︸S·([xrx˙r]−[xx˙])= 0 (4.47)where λ is analogous to the bandwidth of the sliding mode controller, x and x˙ are the positionand velocity of the axis, respectively, and xr and x˙r are the reference position and referencevelocity of the axis respectively. The objective in sliding mode control is to bring the systemonto the sliding surface σ = 0, which then brings the tracking error to zero. In order to designa control law that drives the system to the sliding surface and is stable, a candidate Lyapunovfunction is first selected as a function of the sliding surface and a disturbance estimator asfollows:V (t) =12meσ2 +(d− dˆ)22ρ(4.48)where d is the actual disturbance, dˆ is the disturbance estimator, and ρ is the adaption gain ofthe disturbance estimator. The disturbance estimator acts similar to an integrator and is definedas a function of the sliding surface:dˆ[k + 1] = dˆ[k] + ρκσTs (4.49)52Ts is the sample time, and κ is a flag to ensure the disturbance estimator does not go beyondthe predetermined disturbance bound as follows:κ =0, if dˆ ≤ d− and σ ≤ 00, if dˆ ≥ d+ and σ ≥ 01, otherwise (4.50)As it can be seen the Lyapunov function is directly proportional to the size of sliding surfaceand the error in disturbance estimation. If the control law is selected so that the time derivativeof the Lyapunov function is less than zero, the sliding surface and error in disturbance estima-tion will be guaranteed to converge to zero. In order to do this, we take the time derivative ofEq. (4.48) and set the derivative to be less than 0:dV (t)dt= meσσ˙ +(d− dˆ)ρ˙ˆd < 0 (4.51)The derivative of the sliding surface and disturbance estimate is found as follows:σ˙ = [λ 1]︸︷︷︸S·([x˙rx¨r]−[x˙x¨])(4.52)˙ˆd = ρκσ (4.53)Substituting Eq. (4.46) into Eq. (4.52) into , the derivative of the sliding mode can be redefinedas a function of the rigid body model:σ˙ = [λ 1]︸︷︷︸S·([x˙rx¨r]−[0 10 −be/me]·[xx˙]+[01/me]· u−[01/me]· d)(4.54)Then by substituting the sliding surface (Eq. (4.47)), the disturbance estimate (Eq. 4.49), andtheir respective derivatives (Eq.(4.54) and Eq. (4.53)) into Eq. (4.51), the derivative of theLyapunov equation can be redefined as a function of known system parameters:dV (t)dt= meσ(λ (x˙r − x˙) + x¨r)− σ (u− d− bex˙)− σκ(d− dˆ) (4.55)53If Eq. (4.55) is always less than 0 then the sliding mode controller will be asymptoticallystable. In Eq. (4.55), σd − σκ(d − dˆ) can be set to σdˆ + σ(d − dˆ)(1 − κ). Given the boundsset in Eq. (4.50), then the condition σ(d − dˆ)(1 − κ) ≤ 0 is always guaranteed. To ensureasymptotic stability for the remainder of the Lyapunov equation, the following criteria is set:meσ(λ (x˙r − x˙) + x¨r)− σ(u− bex˙− dˆ)= −Ksσ2 (4.56)Since Eq. (4.56) is the criteria for asymptotic stability of the sliding mode controller, thecontrol signal u can be found that corresponds to this criteria. In order to do this Eq. (4.56) isrearranged to find u as follows:usmc = me [λ (x˙r − x˙) + x¨r] + bex˙+ dˆ+Ksσ (4.57)which ensures that the system will converge to the sliding surface, and subsequently bringtracking error to 0.4.2.4 LuGre Feedforward Friction CompensatorThough the sliding mode controller is capable of rejecting disturbance, it is only capable ofdoing so after some error has been incurred. Due to the non-linear friction phenomenon, atpoints where the velocity switches signs a disturbance occurs which causes a tracking error.Since this non-linear friction phenomenon can be modeled with relative accuracy, it is possibleto compensate for its effects before the occur. In this work, the LuGre friction model [57]is used to model the non-linear friction effects. The LuGre friction model is an extension ofthe traditional Stribeck friction model, as defined in Eq. (4.41), which is able to capture thehysteretic effects of friction. Unlike the Coulomb or Stribeck friction model, the LuGre frictionmodel does not treat the non-linear friction phenomenon as discontinuous force changes atpoints of velocity switching. Instead, it captures the friction dynamics in the stiction region,54x(t) ϵ(t)Figure 4.5: Graphical representation of LuGre friction modelwhich come from the deflection, ǫ(t), of micro structures between surfaces, before the static,Coulomb, and viscous friction affect become dominant. The microstructures are analogous tobristles and ǫ(t) can be seen as the average deflections of these bristles as shown in Fig.4.5.Based upon this model, the deflections of ǫ(t) are governed with the following differentialequation:dǫ(t)dt= x˙(t)− σ0g(x˙(t))ǫ(t) |x˙(t)| (4.58)where σ0 is analogous to the stiffness constant for the microstructures, which is found experi-mentally [58], and the speed dependent function, g(x˙(t)) is expressed as:g(x˙(t)) = Fc + Fsex˙(t)Ωs (4.59)where FC is coulomb friction, FS is the velocity (x˙(t)) dependent static friction, and Ωs is thevelocity in which the effect of static friction decays to 0. Equation (4.59) can be seen as a55simplified, averaged of Eq. (4.41), and as a result, the Stribeck friction effects are incorporatedinto the LuGre model through Eq. (4.59). The constants are found using the results of thefriction identification experiments in Section 4.2.2 as follows:Fs =12(F+stat + |F−stat|) (4.60)Fc =12(F+coul + |F−coul|) (4.61)and Ωs is found by minimizing Eq. (4.45) using Eq. (4.59). With the microstructure deflectionsmodeled in Eq. (4.58), the instantaneous friction force, Ff (t), can be modeled as follows:Ff (t) = σ0ǫ(t) + σ1dǫ(t)dt+ bx˙(t) (4.62)where σ1 is analogous to the damping constant [58] of the microstructures. The stiffness anddamping of the microstructures, σ0 and σ1 respectively, are first approximated with the ap-proach outlined in [58]. For σ0, a very small step xss, assumed to be in the stiction region, canbe given to the feed drive and the constant can be approximated as:σ0 = Fcsgn(x˙(t))/xss (4.63)In the stiction region, the dynamics of the system can be modeled as:mx¨+ (b+ σ1 + σ2) x˙+ σ0x = u(t) (4.64)and σ1 can be calculated as:σ1 = 2√mσ0 − b− σ2 (4.65)56With the constants defined in Eq. (4.63) and Eq. (4.65), the values were tuned in trackingexperiments to achieve the best performance. Since the sliding mode controller is able to sup-press errors contributed by viscous friction, the feed forward friction compensator only has tocompensate for the forces contributed by static friction, Coulomb friction, and microstructuredeflection ǫ(t). As a result, the feed forward compensator has the following structure:uff =1KaKt(σ0ǫ(t) + σ1dǫ(t)dt)(4.66)which represents the LuGre friction modeled with the viscous friction removed and normalizedby the axis amplifier and motor gain. The deflection ǫ(t) is found by numerically solving Eq.4.58 using the Forward Euler Approximation:x˙[0] = 0 ǫ[0] = 0dǫ[0]dt= 0dǫ[t]dt= x˙[k]− σ0g(x˙[k])ǫ[k − 1] |x˙[k]|ǫ[k] = ǫ[k − 1] + dǫ[k]dtTs(4.67)The full control law is the sum of sliding mode (Eq. 4.57) and feed forward compensator (Eq.4.66) commands:u(k) = usmc(k) + uff (k)= me [λ (x˙r(k)− x˙(k)) + x¨r(k)] + bex˙(k) + dˆ(k)+Ksσ(k) +1KaKt[σ0ǫ(k) + σ1dǫ(k)](4.68)The block diagram of the full micromill control law is shown in Fig. 4.6. With the tunedLuGre friction parameters, the sliding mode controller is tuned to achieve as high a bandwidthas possible before instabilities occurred and the control parameters are shown in Tab. 4.3.57Figure 4.6: Block diagram of the sliding mode controller and LuGre friction compensator forthe micromillTable 4.3: Micromill parametersAxis xc yc zcKa [A/V] 1 1 1Kt [N/A] 45.4 72.5 26.17Ks [Vs/mm] 0.05 0.25 0.08λ [1/s] 200 230 180ρ[V/mm] 50 50 50σ0 [N/mm] 900 900 1500σ1 [Ns/mm] 1.5 2.5 7.5Fc [N] 26.67 31.26 24.67Fs [N] 29.15 34.11 32.35Ωs [mm/s] 24.3 26.4 44.8584.3 Modeling and Control Design of the Rotary Table4.3.1 Lead-lag Position ControllerUnlike the feed drive, the rotary table is non-contact. As a result, the traditional rigid bodymodel for linear feed drives can be simplified by eliminating the viscous friction. This elimi-nates the need for sophisticated identification techniques, as it is possible to obtain the mass ofthe rotary table by simply weighing it. It was found that the mass of the rotary table was 2.3kgand as a result the mass in the translational axes can be defined as:Jx = Jy = Jz = 2.3[kg] (4.69)The rotational inertia about the x, y, and z axis were found based on the solid model of therotary table [59]:Ja = Jb = 20878[kg ·mm2] (4.70)Jc = 40508[kg ·mm2] (4.71)Since the rotary table is a free floating mass, the nominal transfer function for each axis of therotary table can be modeled as follows:Gf (s) =1Jqs2← q ∈ [x, y, z, a, b, c] (4.72)Due to its relatively simple dynamics, lead-lag controllers were selected to control the positionof the rotary table which take the form:Cf (s) = Kf1 + αTs1 + Ts(4.73)59where α is selected to add a desired phase, T is selected to add the phase at a specific location,and Kf is selected so Cf (s)Gf (s) crosses the 0dB line at a desired frequency. In order to dothis, the parameters are calculated as follows:α = 1+sinφm1−sinφm T =1ωm√αKf =1| 1+αTωmj1+Tωmj||Gf (ωm)| (4.74)where ωm is the desired crossover frequency of the openloop transfer function, which corre-sponds approximately to the bandwidth of the closed loop transfer function, and φm is thephase to be added. An integrator of the form:CI(s) = 1 +Kis(4.75)with Ki = ωm/10 is added to the control loop to reduce the steady state error. However, itwas found that unlike the nominal models in Eq. (4.72), the rotary table was not completelyrigid, hence its flexibility had to be considered while designing the controller. The resonancepeaks caused by the flexibilities of the rotary table created unstable conditions when the designbandwidth was set to 250Hz. As a result the gain is attenuated at the frequency of the resonancepeaks with notch filters so that the control action does not excite the low damped dynamics andits dynamic characteristics match the nominal model, as shown in Fig. 4.7.4.3.2 Notch FilterWhile designing the lead-lag controllers for the various axes of the rotary table, it was foundthat the nominal transfer function defined in Eq.(4.72) did not fully capture all the dynamicsof the rotary table. At high enough frequencies, the rotary table was not completely rigid andexhibited flexible characteristics. An example of this is shown in Fig. 4.7 where the rigidnominal model is compared with the actual frequency sweep of the af -axis.6010310ÏÐ10ÏÑ10Ï4ÒÓÔÕÖ×ØÙÚÛÙÜÝÞßàáâàãäå æßçèéêëìíîïãçð ÞñÞ íò óôçõïêÞïðöàßàè ÞñÞ íò óôçõïê÷ãòïðöàßàè ÞñÞ íò óôçõïêFigure 4.7: FRF of the af -axis including the nominal FRF and experimental FRFs with andwithout the notch filterAs it can be seen, though the experimental frequency response function (FRF) follows the nom-inal model at lower and higher frequencies, in between 1000 [rad/s] and 2000 [rad/s] resonantmodes from the flexibilities of the rotary table can be seen. With the increased magnitude atapproximately 1700 [rad/s], when the cross over frequency of the open loop transfer functionwas selected to be too high, the peak is pushed close to the 0dB reducing the gain marginand introducing oscillatory behaviour and potential instabilities. As a result, notch filters wereimplemented to remove the frequencies in the control signal that correspond to these resonantfrequencies. The notch filters take the following form:Na(s) =s2 + 2ωaζa,1 + ω2as2 + 2ωaζa,2 + ω2aNb(s) =s2 + 2ωbζb,1 + ω2bs2 + 2ωbζb,2 + ω2bN(s) = Na(s)Nb(s)(4.76)where Na and Nb are the notch filters to remove the peak and valley, respectively, caused by theresonant mode. In the notch filters Na and Nb, ωa, ωb are the resonance frequencies, selectedbased on peaks of the FRF, and (ζ1, ζ2) are the damping ratios, which are tuning parameters61Figure 4.8: Block diagram of rotary table controllerTable 4.4: Rotary table parametersAxis xf yf zf af bf cfJq 2.3[kg] 20878[kg/mm2] 40508 [kg/mm2]ωm [rad/s] 1570.8Kf 5.95E5 6.60E5 8.14E5 3.32E3 2.24E3 6.53E3α 1.31E2 1.31E2 1.31E2 1.31E2 1.31E2 3.75E2T 5.57E-5 5.57E-5 5.57E-5 5.57E-5 5.57E-5 2.73E-5Ki 157.1ωa [rad/s] - - 1992 1705 1429 2582ζa,1 - - 0.02 0.03 0.03 0.08ζa,2 - - 0.3 0.3 0.2 0.02ωb [rad/s] - - 1740 1533 1300 2421ζb,1 - - 0.3 0.3 0.2 0.05ζb,2 - - 0.02 0.04 0.03 0.262Figure 4.9: Dual stage feed drive controllerto achieve a desirable open loop transfer function. The full block diagram for the closed loopcontroller of the rotary table with the lead-lag compensator and notch filter is given in Fig. 4.8and the parameters of the lead lag controllers and notch filters are shown in Tab. 4.4.4.4 Dual-stage Feed Drive ControlBased on the above control designs, the achievable bandwidth of the micromill is 30Hz whilethe rotary table is 250Hz. Since the translational axes of the micromill and rotary table areparallel, the rotary table can be used to compensate for tracking errors coming from the lowerbandwidth of the closed loop controller of the micromill. As a result, in addition to its owntrajectory command, an additional input to the rotary table is the tracking error of the micromillas shown in Fig. 4.9.In this configuration, the summation of the two parallel drives has the bandwidth of the rotarytable, resulting in the long stroke of the micromill with the higher precision of the rotary table.The summation of the two drives gives the following transfer functions:x = xc + xf (4.77)63Based on the block diagram shown in Figure 4.9, Eq. (4.77) can be expressed as a function ofthe reference input as follows:x = xrCcGc(1 + CcGc)︸ ︷︷ ︸Gc,cl+xr1(1 + CcGc)︸ ︷︷ ︸Sc,clCfGf(1 + CfGf )︸ ︷︷ ︸Gf,cl(4.78)where xr is the reference command to the micromill feeddrive. As a result, as long as Gc,cl,Sc,cl, and Gf,cl are stable then x will be stable. Since the micromill, Gc,cl is stable, its trackingerrors, which are used as input to the rotary table, Gf,cl, will be bounded and stable. Further-more, the sensitivity transfer function of the coarse closed loop transfer function Sc,cl will alsobe stable since Gc,cl is designed to be stable. The two actuators are independently controlledand stable, and as a result the dual stage controller is also stable. The transfer function of thecoupled system can also be expressed as follows:xxr=CcGc + CfGf + CcGcCfGf(1 + CfGf )(1 + CcGc)(4.79)At frequencies below the closed loop bandwidth of the micromill, ωc, and rotary table, ωf ,|GcCc| ≫ 1 and |GfCf | ≫ 1. As a result the dual stage system has a closed loop gain of unity.However at frequencies greater than ωc but below ωf , |GcCc| ≤ 1 and |GfCf | ≫ 1 and as aresult the closed loop dynamics of the summed response is as follows:x ≈ GfCf1 +GfCfxr (4.80)As a result, at frequencies above the design bandwidth of the closed loop controller for themicromill, the closed loop dynamics of the summed response adopts the closed loop dynamicsof the rotary table. As a result, the summed system is able to have the stroke of the micromillbut the precision of the rotary table. This can be seen in Fig. 4.10.64Frequency [rad/s]Phase [degs]Magnitude [dB] Rotary Table FRFMicromill FRFMicromill + Rotary Table FRF500-50-100-1500-45-90-135-180101 102 103 104 105 106Figure 4.10: Simulated frequency response functions of micromill, rotary table, and micromillwith tracking error compensated with rotary table.65When the rotary table is used to compensate for the tracking errors of the micromill, the transferfunction of the axes adopts the bandwidth of the higher bandwidth controller. As a result, thedual stage feed drive system achieves the stroke of the micromill but with the high precision -high bandwidth control of the rotary table.66Chapter 5Geometric Error Modeling for a 3-axis Micromill and Compensation with a 6 Degree ofFreedom Rotary Table5.1 OverviewOften in the discussion of the effectiveness of a manufacturing process, the major focus is onthe precision of machine tool performing the manufacturing process. However, what is oftenless discussed is the accuracy of the manufacturing process. In most literature, it is assumedthat accuracy follows precision. However this is typically not the case and additional methodsmust be implemented in order to increase the accuracy of the machine tool to a sufficient level,even at sufficient levels of precision. The source of the problem lies in the fact that mostfeedback servomechanisms are only capable of detecting errors in the direction of actuation.As a result, errors that are in directions orthogonal to the direction of travel can go undetected.These errors can originate from assembly errors, deflection caused by loading, or deviationscaused by thermal expansion. Though these errors only create minimal problems for single axisof actuation, when multiple servomechanisms are combined and non-Cartesian movements areincluded, the geometric errors could result in large deviations from the desired end effectorposition and orientation when the whole kinematic model is considered.In this chapter, the fine movements of the rotary table are used to compensate for tooltip errorscaused by the assembly errors. Due to the invariant nature of errors caused by loading orthermal expansion, the rotary table would not be suitable for this task as these errors wouldchange with operating parameters or environmental conditions. The ideal kinematic model67from Chapter 3, is modified in Section 5.2 to include the effect of the geometric errors of eachof the axis. Since the geometric errors originate from assembly errors of the axes, it changesconsistently over the full actuating range of the axis. In Section 5.3, the geometric error ismeasured, then fit to a polynomial function of position. Finally, with the geometric errorsknown, a compensation algorithm is proposed in Section 5.4, where the rotary table positioncommands are modified to compensate for tooltip position errors.5.2 Kinematic Model with Errorsxyzx(q)y(q)εz(q)εy(q)εx(q)z(q)Figure 5.1: Geometric errors of a general axis qA single axis has a total of 6 geometric errors, which include one positioning error, two straight-ness errors, a roll error, a pitch error, and a yaw error. These errors can be seen on a general68axis, q, in Fig. 5.1.The positioning error is typically caused by inaccuracies in the metrology system, while theother geometric errors are caused by assembly errors. Though the errors are small in mag-nitude, with multiple axes stacked and rotational errors being amplified with linear offsets, itis possible to have large resultant errors at the end effector. In order to demonstrate this, thekinematic model defined in Eq. (3.6) to (3.8), is modified to include the geometric errors [41].The general error matrix is defined as follows:Eq =1 −εz(q) εy(q) δx(q)εz(q) 1 −εx(q) δy(q)−εy(q) εx(q) 1 δz(q)0 0 0 1 (5.1)where εx(q), εy(q), and εz(q) are the rotations about the x, y, and z axis of the current coordi-nate frame, respectively, and δx(q), δy(q), and δz(q) are displacements in the x, y, and z axis ofthe current coordinate frame, respectively, as shown in Fig. 5.1. It should be noted that Eq.(5.1)uses the small angle approximation for the cosine and sine of the rotational errors, ε, since therotational errors are typically sufficiently small. By multiplying Eq.(5.1) with the ideal HTM, itis possible to project the effect of errors in one axis onto the next axis, and subsequently, to thetooltip. For example, transforming from base coordinate frame to the y-axis of the micromillhas an ideal transformation as follows:BTyc,i =1 0 0 00 1 0 yc0 0 1 Ly,z0 0 0 1 (5.2)69However, when the error matrices are multiplied, the actual transformation becomes as follows:BTyc,a =BTyc,iEyc=1 −εz(yc) εy(yc) δx(yc) + ycεz(yc) 1 −εx(yc) δy(yc)−εy(yc) εx(yc) 1 δz(yc) + Ly,z0 0 0 1(5.3)As a result, in order to model the effect of the geometric errors of each axis on the positionof the tooltip, the ideal transformation matrices are multiplied with their respective error trans-formation matrices, making changes to the two kinematic chains in Eq.(3.1) and Eq.(3.2) asfollows:BTw,a =BTyc,iEycycTxc,iExcxcTs,isTxf ,i×sTxf ,iExf xfTyf ,iEyf yfTzf ,iEzf×zfTaf ,iEaf afTbf ,iEbf bfTcf ,iEcf×cfTr,irTw,i(5.4)BTt,a =BTzc,iEzczcTt,i (5.5)which will give the actual HTM from the tool coordinate frame to the workpiece coordinateframe:Ta =(BTw,a)−1 BTt,a (5.6)With the actual kinematic transformation defined in Eq.(5.6), the tooltip position with the effectof error is transformed to a position that is with respect to the workpiece with a method similarto Eq.(3.5): Px,aPy,aPz,a1 = Ta[T tp1](5.7)70which gives the actual tooltip position, Pa (q) = [Px,a (q) , Py,a (q) , Pz,a (q)]T . The result ofthe matrix multiplication in Eq. (5.7), with second order and higher terms error grouped asO(ε2, δ2, εδ), can be found in Appendix A. The forward kinematic model now also includesthe effect of the geometric errors on the tooltip as a function of axes positions. In this work,the forward kinematic equations consists of 9-axis commands and 18 position dependent errorterms. As a result, the equations can become very long, and subsequently, unmanageable andcomputationally expensive. In order to overcome this problem, the second order or higherterms involving geometric errors are approximated as follows:O(ε2, δ2, εδ) ≈ 0 for ε, δ → 0 (5.8)since the geometric errors are very small relative to the position commands and offsets. Giventhis approximation, the forward kinematic equations can be simplified and approximated byincluding only the zero and first order terms of geometric errors in the forward kinematicequations. With the kinematic equations modified to include the geometric errors, the positiondependent error equations must be identified through experimentation.5.3 Error ModelingFrom Eq. (A.1), Eq. (A.2), and Eq. (A.3), it can be seen that the geometric errors play a strongrole in determining the resultant position of the tooltip position. As a result these errors must bemeasured and fit to a function of position in order to be incorporated into the kinematic model.In order to do this, the errors of the machine tool are first measured with a laser interferom-eter using various mirror configurations to get geometric errors in different directions. Sincegenerated trajectories are on average jerk continuous [51], a quintic polynomial is selected asthe function to be fit to the error data in order to maintain the same jerk continuity. For eachgeometric error, there should be Nǫ ×Mǫ corresponding error measurements, ǫ, where Nǫ is71the total number of measurement points on the axis and Mǫ is the number of repeated mea-surements at each position. In this work, the errors are measured 6 times, every 1[mm] and themeasurement ǫ can take the form of a displacement error, δ, or a rotational error ε. The objec-tive is to find a quintic polynomial that best fits this data. First the Mǫ repeated measurementsat Nǫ locations are averaged:ǫ∗nǫ =ǫ1,nǫ + ǫ2,nǫ + · · ·+ ǫMǫ,nǫMǫ(5.9)where nǫ = 1, 2, . . . , Nǫ. With the averaged errors, the objective is to fit a quintic polynomialto the data points with as little error as possible. The quintic polynomial fitting problem can bephrased as follows:ǫ∗1ǫ∗2...ǫ∗Nǫ−1ǫ∗Nǫ︸ ︷︷ ︸Yǫ=q51 q41 q31 q21 q1 1q52 q42 q32 q22 q2 1.........q5Nǫ−1 q4Nǫ−1 q3Nǫ−1 q2Nǫ−1 qNǫ−1 1q5Nǫ q4Nǫq3Nǫ q2NǫqNǫ 1︸ ︷︷ ︸ΦǫϕAϕBϕCϕDϕEϕF︸ ︷︷ ︸θǫ(5.10)where Yǫ is a vector of averaged measured geometric errors, Φǫ is a matrix of the axis positionswhere the laser measurements took place, and θǫ is a vector of coefficients for the fitting poly-nomial . The objective is to find the coefficients in θǫ that minimizes the mean squared errorbetween the predicted geometric error and actual measured geometric error or more formally:min12(Yǫ − Φǫθǫ)T (Yǫ − Φǫθǫ) (5.11)This is done in a least squares fashion in the following way:θǫ = (ΦTǫ Φǫ)−1ΦTǫ Yǫ (5.12)720 10 20 30 40 50 60 70Z-Axis Position [mm]-20246810Positioning Errors [um]Average Laser Interferometer DataRaw Laser Interferometer DataQuintic Polynomial FitFigure 5.2: zc-axis positioning errors and the resultant quintic polynomial fitand as a result the equation of error becomesǫˆ (q) = ϕAq5 + ϕBq4 + ϕCq3 + ϕDq2 + ϕEq1 + ϕF (5.13)where ǫˆ indicates the approximated geometric error. This process is repeated for all 6 errorson each of the major axes giving a total of up to 18 quintic polynomial functions. An exampleof this process can be seen in Fig. 5.2 where the positioning errors of the zc-axis are fit to aquintic polynomial.Once all the geometric errors are measured and fitted, the errors can be included into Eq. (A.1),Eq. (A.2), and Eq. (A.3). As a result, the actual tooltip position is known and its errors can becompensated.735.4 Geometric error compensationOnce the geometric errors of each of the three axes are known, the actual tooltip position isknown. Subsequently it is possible to find the error of the tooltip position from the desiredtooltip position by subtracting the actual tooltip position from the desired tooltip position asfollows:∆P (q) = |Pi −Pa(q)| (5.14)The ideal tooltip position Pi can be provided from a computer aided manufacturing (CAM)toolpath. Alternatively, if only the axes positions are known Pi can also be the result of theideal forward kinematic model defined in Eq. (3.6) to Eq. (3.8), prior to the reference commandmodifications performed for geometric error compensation. The objective is to find rotary tablereference commands, qf = [xf , yf , zf , af , bf , cf ]T , that will reduce the geometric error to anacceptable minimum. In order to do this, two challenges exist. First, similar to the trajectorygeneration of 9-axis machine tool, there are more degrees of freedom available than necessary.Since only the tooltip position is defined with 3 degrees of freedom and the rotary table has 6degrees of freedom, there are 3 redundant degrees of freedom. As a result, redundancies mustbe resolved with the proposed geometric error compensation technique. Second, due to thenon-Cartesian kinematics of the overall 9-axis configuration, the compensation of geometricerrors is non-trivial and numerical techniques are necessary. To overcome these challenges, aniterative gradient descent algorithm which exploits the Jacobian matrix of the rotary table isproposed.5.4.1 Gradient Descent Optimization Algorithm BackgroundThe gradient descent optimization algorithm is a first order iterative optimization algorithmthat is used for finding the minimum of a function. The algorithm iterates by selecting function74variables that are proportional to the negative of the gradient of the function at the current point.This process is repeated until the gradient is sufficiently small, which corresponds to a localminimum of the function. If the variables to be modified are defined as η and the differentiablemultivariable cost function is defined as Ψ(η) then the cost decreases fastest if η goes in thedirection of the negative gradient of the cost function as follows:ηj+1 = ηj − γ∇Ψ(ηj) (5.15)where j is the current number of iteration for the gradient descent algorithm. For sufficientlysmall γ the cost will decrease at each iteration. As a result, starting with a guess of η0 the costwill decrease as follows:Ψ(η0) ≥ Ψ(η1) ≥ Ψ(η2) ≥ . . . (5.16)5.4.2 Geometric Error Compensation using Gradient Descent Optimization AlgorithmWith respect to geometric error compensation, the gradient descent algorithm is used to modifythe rotary table commands, qf , so that the tooltip position errors ∆P (qf ) are decreased to anacceptable minimum. As a result the variables to be modified are qf and the cost function isdefined to penalize tooltip errors as follows:Ψ(qf ) =12∆P (qf )T ∆P (qf ) (5.17)When applied to the problem of geometric error compensation, the gradient descent algorithmcan be phrased as follows:qf,j+1 = qf,j −∇Ψ(qf,j) (5.18)75where γ is set to 1, j is the number of iterations of the gradient descent algorithm and thegradient of the cost function Ψ is:∇Ψ(qf ) = dqfdPa∆P (qf ) (5.19)Since, ∆P (qf ) is known from Eq. (5.14), the differential term dqfdPa must be defined in order toobtain the gradient. Since there are more degrees of freedom than necessary, an infinite rangeof joint configurations are possible, or more formally, dqfdPais non unique.Similar to the redundancy resolution technique presented in Chapter 3, the Moore-Penroseinverse is used to resolve the redundancies. Since the geometric error compensation only usesthe rotary table commands, the Jacobian used for compensation will only consider the effectof the rotary table commands on the tooltip position and tool orientation. Furthermore, theJacobian will also incorporate the geometric errors in the forward kinematics equation. As aresult, the compensating Jacobian is defined as follows:J c =dPadqf=dPx,adxf. . . dPx,adcfdPy,adxf. . . dPy,adcfdPz,adxf. . . dPz,adcf3×6(5.20)where the superscript c differentiates the Jacobian defined in Eq. (5.20) from the Jacobiandefined in Eq. (3.13). With the compensating Jacobian, J c, a Moore-Penrose inverse can befound as a possible solution to dqfdPaJ c† =J cTJ cJ cT(5.21)As a result, the gradient of the objective function can be redefined as follows:∇Ψ(qf,j) = J c†∆P (qf,j) (5.22)76However, in addition to minimizing the error between the desired tooltip position and the actualtooltip position, the algorithm must ensure that the generated compensation commands do notviolate the stroke limits of the rotary table. Similar to the trajectory generation algorithmpresented in Chapter 3, an additional constraint can be projected onto the nullspace of J c†. Thecost function which constraints the magnitude of the stroke limits is considered as follows:Hc (qf ) =x2f(x2f,max − x2f)2 + · · ·+ c2f(c2f,max − c2f)2 (5.23)By projecting the gradient of this cost function:∇Hc (qf ) =[2x2f,maxxf(x2f,max − x2f)2 , . . . , 2c2f,maxcf(c2f,max − c2f)2](5.24)onto the nullspace of J c†, it is ensured that the solutions for∇Ψ(qf ) will not violate the strokelimits of the rotary table. It should be noted that the constraint equation, Eq. (5.23), and itsgradient, Eq. (5.24), only consider the limits of the rotary table. Since only the referencecommands of the rotary table are being modified, it is suffice to only consider the rotary tableposition commands. The gradient of the cost function is extended to consider this additionalconstraint in the following way:∇Ψ(qf ) = J c†∆P (qf ) +(I − J c†J c) βc∇Hc (qf ) (5.25)where βc determines the strength of the stroke limit constraint defined in Hc. The gradientdescent formulation can be redefined as follows:qf,j+1 = qf,j − J c†∆P (qf,j)−(I − J c†J c) βc∇Hc (qf,j) (5.26)Since the the cost, Ψ(qf,j), decreases the fastest when the rotary table commands, qf,j+1, arerecalculated in the direction of −∇Ψ(qf,j) then every iteration will minimize the cost. As a77result, Eq.(5.26) is iterated until the following condition is met:∆P(q∗f)< ν (5.27)Algorithm 2: Full compensation algorithm1 function Compensation (q);Input : Q = [q[0],q[1], . . . ,q[K − 1],q[K]]Output: Q∗ = [q∗f [0],q∗f [1], . . . ,q∗f [K − 1],q∗f [K]]2 k = 0;3 j = 0;4 for k ← 0 to K − 1 do5 ∆P = |Pi(q[k]0)−Pa(q[k]0)| ;6 while ∆P (q[k]j) > e do7 q[k]j+1 = q[k]j − J c† (q[k]j)∆P (q[k]j)−(I − J c†J c) βc∇Hc (qj) ;8 ∆P = |Pi(q[k]0)−Pa(q[k]j+1)| ;9 j = j + 110 end11 q∗[k] = q[k]j12 j = 013 k = k + 114 endwhere q∗f is the rotary table commands that reduce the geometric errors to an acceptable limitdefined by ν. Compensation commands are calculated for every point, q[k], which is knownfrom the trajectory algorithm. The full algorithm can be summarized in Algorithm 2.780 1 2 3 4 ø ù-101úûüýþßTþ10-40 1 2 3 4 ø ù-4-2024ùYûüýþßTþ10- 0 1 2 3 4ø ùD 024ùZûüýþßTþ10- CfffiflffiUfffiflffiCfffiflffiUfffiflffiUfffiflffiCfffiflffiFigure 5.3: Tooltip errors with and without compensation for a circle on the x-y plane of radius1Simulation results using this method can be seen in Fig. 5.3, where the algorithm eliminatesthe tooltip errors caused by the geometric errors of the machine tool when drawing a circle ofradius 1 [mm] on the x-y plane.79Chapter 6Simulation and Experimental Results6.1 OverviewExperiments were performed to validate the trajectory generation algorithm, dual stage feeddrive tracking control law, and geometric error compensation method presented in this thesis.All experiments were performed on the machine presented in the Introduction.6.2 Trajectory Generation Experimental ResultsIn order to validate the trajectory generation algorithm, two experiments were performed toensure that the generated trajectories were able to resolve redundancies, avoid singularities,respect the prescribed limits, and follow the desired end effector trajectory. First, the trajectorygeneration algorithm is used to generate axes commands for a spiral toolpath as shown in Fig.6.1The position commands for the spiral toolpath and the time derivatives toolpath are shown inFig. 6.2, Fig. 6.3, Fig. 6.4, and Fig. 6.5. As it can be seen, the position commands respect thestroke limits of all the drives, demonstrating the ability of the redundancy resolution componentof the trajectory generation algorithm to select axes configurations that respect the stroke limitsof the machine tool. In addition to this it can be seen that the generated axes commandsrespect the prescribed velocity, acceleration, and jerk limits of all the axes showing successfulintegration of the redundancy resolution algorithm with the feed optimization algorithm. Itshould also be noted that at certain portions of the trajectory, the velocity of the yc is saturated80020406 8 !40!2002040510152025X [mm]" #$$%Z [mm]0 5 100 15 200 25 !20020406 P&'()(&*+,,.Displacement [mm]PxPyPz0 5 100 15 200 25 !505x 10/3O0(3*)4)(&*+047.9:;<=>?@$@AB#$$%EiEjEk!1Figure 6.1: Spiral toolpath810 20 40 FG HG 100 120 140-1000100IJKLMLJNQRRSxcycVc0 20 40FG HG100 120 140WGX[0GX[IJKLMLJNQRRSxfyfVf0 20 40 FG HG 100 120 140WGXG[0GXG[IJKLMLJNQ\]^KSafbf0 20 40FG HG100 120 140Displacement [mm]-10010Position [rads]cfFigure 6.2: Reference commands for spiral trajectory0 50 100 150 200 250 300−1000100Velocity [mm/s] xcyczcxfyfzf0 50 100 150 200 250 300−40−2002040Displacement [mm]Velocity [rad/s] afbfcfFigure 6.3: Velocity of reference commands for spiral trajectory820 50 100 150 200 250 300−500005000Acceleration [mm/s2] xcyczcxfyfzf0 50 100 150 200 250 300−20020Displacement [mm]Acceleration [rad/s2 ] afbfcfFigure 6.4: Acceleration of reference commands for spiral trajectoryat its limits, and for portions of the trajectory where the velocity is not saturated, the jerk of theother axes are close to being saturated. This demonstrates the ability of the feed optimizationalgorithm to use all the kinematic limits of the feed drive to traverse the toolpath.In addition to respecting the limits of the axes, the generated axes commands should result ina toolpath that follows the desired toolpath. In the case of the trajectory generation algorithmpresented, the numerical nature of the trajectory generation algorithm and singularity avoid-ance components can cause the toolpath from the axes commands to deviate from the desiredtoolpath. As a result, the desired tooltip position and tool orientation are compared with thetooltip position and tool orientation from applying forward kinematics on the generated axesposition commands. The tool deviations for the spiral toolpath were shown earlier in Fig. 3.7showing negligible deviation, and as a result, successful implementation of the a higher order830 50 100 150 200 250 300−202x 104Jerk [mm/s3] xcyczcxfyfzf0 50 100 150 200 250 300−2000200Displacement [mm]Jerk [rad/s3 ] afbfcfFigure 6.5: Jerk of reference commands for spiral trajectory84Figure 6.6: Sinusoidal freeform surfacenumerical integration algorithm and closed loop corrective action.Next the trajectory generation algorithm is used to generate reference commands to machinethe freeform surface shown in Fig. 6.6.The freeform surface is a sinusoidal surface with a peak to peak amplitude of 500 microns anda working surface area of 10 [mm] x 10 [mm]. This surface was cut with a 1/64” (397[µm])Mitsubishi Carbide 2-fluted ball endmill. The finishing pass was done with a feedrate of 10[mm/s] with a spindle speed of 170000 [rev/min]. To avoid the ploughing affect from the zerocutting velocity from the tip of the ball end mill, the tool is commanded to have a constant tiltof 0.15 degrees using the magnetically levitate table’s six degrees of freedom. Since this was afinishing operation, the limitation of the feedrate will also come from the process itself, whichwas set at 10 [mm/s], in addition to the axes limits. The freeform surface was cut with a zig-850 2 4 _ ` 10 12024_`1012adeg hilDesire ateFigure 6.7: Feedrate profile for freeform surface with limit set to 10 [mm/s]860 2 4 m n 10 12opqq0pqqVelocity [mm/s]xcycrcxsysrstuvuw{0 2 |m8 10 12-101}~~
Ł2]100 2 | m 8 10 12 [s]-202Ł3105Figure 6.8: Velocity, acceleration, and jerk profile of the translational axes for the freeformsurface870 2 4 10 12-10010Velocity [rad/s]abc0 2 8 10 12-2000200 ¡ ¢£¤¥¦§¨¢£©ª«2]0 2 8 10 12¬®¯ °±²-200002000³ ¢´¨¢£©ª«3µFigure 6.9: Velocity, acceleration, and jerk profile of the rotational axes for the freeform surface880 2 4 ¶ · 10¸¹0¹Velocity [mm/s]ºcyc»cº¼y¼»¼0 2 ½ ¶ 8 10¾¿ÀÁ ÂÃÄ-100001000ÅÆÆÇÈÇÉÊËÌÍÎÏÐÐÑÒ2]Figure 6.10: Velocity, acceleration, and jerk profile of the translational axes for the freeformsurfacezag toolpath with a depth of cut of 100 [µm] and a step over of 10 [µm]. The planned feedrateis shown in Fig. 6.7 and as it can be seen, the feedrate is capped at 10 [mm/s]. However, atportions where the tool path must slow down and speed up, the limitation comes from the axes.This can be seen in the velocity, acceleration, and jerk of the reference commands in Fig. 6.8and Fig. 6.9 for the translational and rotational axes respectivelyFurthermore, it should be noted that when zoomed in, the velocity and acceleration profiles aresmooth and continuous as shown in Fig. 6.10 and Fig. 6.11 for the translational and rotationalaxes respectivelyThe manufactured surface can be seen in Fig. 6.12. As it can be seen, the produced partmatches the desired surface, showing the kinematic model is correct. Due to the finite step890 2 4 Ó Ô 10ÕÖ×Ø0Ö×ØVelocity [rad/s]aÙbÙcÙ0 2 Ú Ó 8 10ÛÜÝÞ ßàá020ÚÖâããäåäæçèéêëìæçíîï2]Figure 6.11: Velocity, acceleration, and jerk profile of the rotational axes for the freeformsurface90Figure 6.12: Machined freeform surface91Table 6.1: Circular contouring resultsError SMC DSFDMean | x-axis Tracking | [mm] 2.14E-4 1.98E-4Mean | y-axis Tracking | [mm] 1.06E-4 4.94E-5Mean Contour [mm] 1.58E-4 1.21E-4Max | x-axis Tracking | [mm] 2.67E-3 1.28E-3Max | y-axis Tracking | [mm] 1.72E-3 3.77E-4Max Contour [mm] 2.67E-3 8.58E-4Table 6.2: Square contouring resultsError SMC DSFDMean | x-axis Tracking | [mm] 2.58E-4 1.69E-4Mean | y-axis Tracking | [mm] 8.81E-5 7.18E-5Mean Contour [mm] 7.92E-5 9.93E-5Max | x-axis Tracking | [mm] 5.43E-3 2.26E-3Max | y-axis Tracking | [mm] 1.54E-3 5.56E-4Max Contour [mm] 7.59E-4 8.62E-4over length it can be seen that there are surface marks left by the scallop heights in betweenpasses of the zig-zag toolpath.From the above experiments, it can be said that the presented trajectory generation algorithmis able to fulfill all of its goals which include redundancy resolution, singularity avoidance, andoptimization6.3 Dual Stage Feed Drive Tracking Control ResultsTo verify that the dual stage feed drive controller improves tracking error performance, the dualstage feed drive controller is used to follow contours. In this set of experiments, the parametersused for the micromill and rotary table are as shown in Tab. 4.3 and Tab. 4.4.92a) b)−40 −20 0 20 40−30−20−10010203040~180nm~25nmY [mm]X [mm]ReferenceWith Compensation from Rotary TableWithout Compensation from Rotary Table0 0.5 1 1.5 2 2.5 3 3.5−101Without compensation from rotary tableWith compensation from rotary table0 0.5 1 1.5 2 2.5 3 3.5−50510Time [s]X Tracking Error [μm]Y Tracking Error [μm]215Rotary table supresses errors from directional changeFigure 6.13: a) Circular contour with sliding mode controller and proposed dual stage feeddrive control scheme. Radius: 30 mm, Tangential Feed Speed: 600 mm/min. Angular fre-quency of the traverse : ω(Hz) = f/R = 10[Hz]. b) Tracking error of x and y axis duringcircular contour0 1 2 3 4 ðñ4ñ2024 Without compensation from rotary tableWith compensation from rotary table0 1 2 3 4 5ñ101 Rotary table supresses errorsfrom directional changeX Tracking Error [òm]Y Tracking Error [òm]a) b)ó40 ó20 0 20 40ó30ó20ó10010203040Y [mm]X [mm]ReferenceWith Compensation from Rotary TableWithout Compensation from Rotary Table~400nmFigure 6.14: a) Square contour with sliding mode controller and proposed dual stage feed drivecontrol scheme. Tangential Feed Speed: 600 mm/min. b) Tracking error of x and y axis duringsquare contour93With these parameters, the micromill has a design bandwidth of 30Hz and the rotary table hasa design bandwidth of 250Hz. By having a bandwidth that is approximately 10 times higherthan micromill, it is expected that the rotary table will compensate the tracking errors causedby the low bandwidth of the micromill. Two experiments are performed where a circle, asshown in Figure 6.13, and a square, as shown in Figure 6.14, are tracked with the x-y axesof micromachining center. In these experiments xc and yc are commanded with the trajectorycommands while xf and yf compensate for the tracking errors of xc and yc respectively. Theresults can be seen in Fig. 6.13 and Fig. 6.14 and are summarized in Tab. 6.1 and Tab. 6.2,respectivelyFrom Tab. 6.1 and Tab. 6.2, it can be seen that the sliding mode controller is able to keep errorsunder 3[µm]. The rotary table assists in reducing all tracking errors as shown in Tab. 6.1 andTab. 6.2. Most notably, the higher bandwidth eliminates friction induced error peaks in thecircular contouring where the velocity reversals occur leading to improved contouring perfor-mance as shown in Figure 6.13 a) and Figure 6.13 b). Furthermore the dual stage configurationsees improvement in reducing fluctuations around zero tracking error. The fluctuations occurdue to quantization noise in the control action. Since these fluctuations occur at frequencieshigher than the bandwidth of the rotary table, they are unable to be completely compensatedfor.It should be noted that in the case of the square contouring results, as shown in Tab. 6.2, thereis a slight increase in contour error. Unlike a circular contour, for a square contour, one of theaxis is stationary at all times. While stationary, the axis error will converge to as close to zeroas possible since there is an integrator and very low level of disturbance. Since the rotary tablehas its own oscillatory movement caused by its own feedback noise, this will be added to thestationary coarse actuator resulting in precision that is worse than just the coarse actuator alone.94Table 6.3: Geometric errors with and without compensationErrorMean Mean Max MaxUncompensated Compensated Uncompensated Compensateddx (xc) [µm] 0.549 0.168 1.07 0.443dy (xc) [µm] 0.2 0.0571 0.629 0.265dz (xc) [µm] 0.396 0.162 1.1 0.37dx (yc) [µm] 1.16 0.222 3.06 0.626dy (yc) [µm] 1.65 0.44 4.18 0.883dz (yc) [µm] 0.19 0.0791 0.693 0.358dx (zc) [µm] 3.91 0.844 7.48 2.4dy (zc) [µm] 0.251 0.197 0.9 0.636dz (zc) [µm] 0.337 0.115 0.976 0.498As a result, for a trajectories where the coarse actuator already achieves very high precision onits own, the rotary table may not increase contouring performance. In contrast, for toolpathswhere multiple axes are moving there is a performance increase with respect to contouringerror as it can be seen from the circular contouring results.Based on these experiments, it can be concluded that the strength of the dual stage feed driveconfiguration lies in disturbance rejection for freeform toolpaths. It should be noted that theaddition of the rotary table is complementary to any control strategies used on the micromill.Though the sliding mode controller alone has good position tracking properties, the addition ofthe rotary table increases the overall precision on the respective axis particularly with regardsto disturbances. Should a different control strategy, besides the sliding mode controller be usedon the micromill, the addition of the rotary table would still be complementary to the overallprecision on the respective axis.6.4 Geometric Error Compensation Results95-40 -30 -20 -10 0 10 20 30 40-2-101d x (x c) [um]-40 -30 -20 -10 0 10 20 30 40-0.500.51d y (x c) [um]-40 -30 -20 -10 0 10 20 30 40xc [mm]-2-101d z (x c) [um]Interferometer Average Uncompensated Data Quintic Polynomial Fit to Uncompensated Average DataInterferometer Average Compensated Data Quintic Polynomial Fit to Compensated Average DataFigure 6.15: Geometric Errors for xc-axis positions -40mm to 40mm96-40 -30 -20 -10 0 10 20 30 40-505d y (y c) [um]-40 -30 -20 -10 0 10 20 30 40-505d x (y c) [um]-40 -30 -20 -10 0 10 20 30 40yc [mm]-0.500.51d z (y c) [um]Quintic Polynomial Fit to Compensated Average DataInterferometer Average Compensated DataInterferometer Average Uncompensated Data Quintic Polynomial Fit to Uncompensated Average DataFigure 6.16: Geometric Errors for yc-axis positions -40mm to 40mm970 10 20 30 40 50 60 70-50510d z (z c) [um]0 10 20 30 40 50 60 70-101d x (z c) [um]0 10 20 30 40 50 60 70zc [mm]-1012d y (z c) [um]Quintic Polynomial Fit to Uncompensated Average DataInterferometer Average Uncompensated DataQuintic Polynomial Fit to Compensated Average DataInterferometer Average Compensated DataFigure 6.17: Geometric Errors for zc-axis positions -4mm to 76mm98In addition experiments for the trajectory generation algorithm and dual stage feed drive con-trol, experiments were performed to verify the geometric error compensation algorithm. Inthese experiments, geometric errors of the 3-axis machine tool are compensated using the highprecision stroke of the 6DOF rotary table. First, geometric error measurements were performedwith a laser interferometer. It should be noted that due to the inability of the laser interferom-eter to measure roll errors along an axis, this error was set to 0 in the kinematic model. Themeasured errors are fit to quintic polynomials then included into the kinematic model, which isused with the gradient descent algorithm to calculate compensating commands for the 6DOFrotary table. To demonstrate the capabilities of the geometric error compensation algorithm,single axis displacement errors are measured with a laser interferometer then compensated withthe rotary table with reference commands generated by the algorithm. The objective of theseexperiments is to demonstrate that positional errors are compensated as the combined effort ofmultiple axes of the rotary table. The displacement errors, meant to mimic tooltip deviations,are measured with the laser interferometer with and without the compensating action of therotary table. The results are summarized in Figs. 6.15 to 6.17 and Table 6.3Based on the results shown in Table 6.3 there is a 64% and 60% improvement in the meanand max geometric error respectively, demonstrating the ability of the algorithm to reducegeometric errors. Furthermore, Figs. 6.15 to 6.17 show consistent improvement of geometricerror across the entire actuating range of each of the axis.To look at the results of the algorithm in greater detail, Fig. 6.18 show the rotary table com-mands used to compensate for the errors when actuating on the zc-axis. As it can be seen, thealgorithm favors the tilt commands of the rotary table to compensate for the straightness errorsin the x and y direction when actuating the zc-axis. The reason for this is that it is possible tomake large displacement corrections with relatively small tilts of the rotary table. Finally, it99-20 -10 0 10 20 30 40 ôõ024öCompensating ÷øùúûúøüCommands [um]-20 -10 0 10 20 30 40ôõýþßZ -100-ôõ0Compensating Rotational Commands [urad]xf, fzfbfcfafFigure 6.18: Compensating position commands when moving only zc axisshould be noted that the error is compensated as the combined effort of multiple axes, showingeffective use of the compensation Jacobian.Based upon the single axis experimental results and 2-axis simulation results shown in Fig.5.3 it can be seen that the tooltip error caused by the geometric errors is relatively small, in thesub-micron range. Due to the inconsistent nature of the cutting process itself, errors originatingfrom the cutting process would be the dominant source of error in the final machined piece.Measurements of a slot cut with a coordinate measuring machine found variations in the errorof nearly 30 microns. Furthermore, the configuration of the gantry type milling machine didnot allow for the use of a commercial ball bar, which would have allowed measurement ofcombined errors, independent of process forces.As an alternative, the full multi-axis compensation capabilities are demonstrated virtually. Inorder to do this, all 21 geometric errors of a three dimensional toolpath are modeled and the1000 1 2 3 4 5 60000500|Px|UncompensatedCompensated0 1 2 3 4 5 600.0050.01|Py| [um]UncompensatedCompensated0 1 2 3 4 5 6Time [s]00.050.1|Pz| [um]UncompensatedCompensatedFigure 6.19: Tooltip errors with and without compensationTable 6.4: Mean and maximum tooltip errors with and without compensation for a multi-axistrajectoryErrorMean Mean Max MaxUncompensated Compensated Uncompensated Compensated|∆Px| [µm] 0.97 0.19 6.5 3.7|∆Py| [µm] 0.79 0.31 6.5 7.1|∆Pz| [µm] 5.2 3.8 65.4 81.8101compensation commands are generated. The compensation commands are then used by themachine as reference commands. The encoder readings are fed into the kinematic model witherrors then compared with the ideal kinematic model. In these experiments, a spiral toolpathas shown in Fig. 6.1 is used as the original toolpath for the 3-axis micromill.Based on the results shown in Tab. 6.4 and Fig. 6.19 there is a 56% and 3% improvementin the mean and max geometric error respectively, demonstrating the ability of the algorithmto reduce geometric errors in multi-axis trajectories. Due to the presence of a tracking errorpeak at approximately 3.5 seconds this causes the maximum tooltip error be large. Since thecompensation algorithm only addresses geometric errors, tracking errors are not accounted forand can appear in the error results. This problem is further demonstrated in the zc-axis, wherethere is noisier error behavior than the other axes. This originates from the poor quantizationof the control signal of the zc-axis. Due to the limited number of digital-to-analog converters(DAC) on the dSPACE DS1103, the zc-axis is controlled with a PWM signal that is convertedinto voltage signal by a Axiomatic Universal Signal Converter. With the lower resolution of thePWM signal, approximately 9-bits, and further noise and distortion introduced by the universalsignal converter, the fluctuations of the tracking error will be higher than the other axes on themachine tool. However, on average, as shown in Tab. 6.4 and Fig. 6.19 the tooltip errors whichoriginate from the geometric errors are compensated for quite well.102Chapter 7Conclusions7.1 ConclusionsIn this thesis, a trajectory generation algorithm, controls strategy, and geometric error com-pensation technique have been developed for a novel 9 degree of freedom micromachiningcenter. The hybrid micro-machine tool combines a conventional 3-axis gantry type micromilland a 6DOF high-bandwidth, short stroke magnetic rotary table. Due to its unique 9-axis con-figuration, new CNC strategies have been developed in the thesis. A 9-axis novel trajectorygeneration algorithm, which can handle the four redundant drives while respecting the drivelimits, have been developed. In order to increase the precision on the translational axes, a con-trol strategy is proposed which combines the high bandwidth and precision of the rotary tableand the long stroke of the micromill. Finally, the rotary table is used to compensate for tooltiperrors caused by the geometric errors of the machine tool. The contributions are summarizedas follows:The proposed trajectory generation algorithm was developed in order to overcome the chal-lenges associated with generating trajectories for a configuration that has more degrees of free-dom than necessary. The position and orientation of a typical cutting tool can be defined by6 degrees of freedom and achieved by 5 degrees of freedom on a conventional 5-axis CNCmachine. Since the developed micro-machine tool has 9 degrees of freedom, 4 more degreesof freedom than necessary, traditional inverse kinematics and feed planning algorithms werenot applicable. In this thesis, a methodology was developed to overcome these challenges. A103forward kinematic model of the machine tool is first developed, which maps the 9-axis posi-tions to the the tooltip position and tool orientation with respect to the workpiece. A numericaltechnique is developed, which resolves the redundancies at the differential level with respectto displacement, using the Moore-Penrose inverse of the Jacobian of the forward kinematicmodel. The proposed differential solution ensures that singularities are avoided and the gen-erated trajectory does not violate the stroke limits of the axes. A corrective 4th Order RungeKutta numerical integration algorithm is used to extract the position commands from the dif-ferential solution, giving axes position commands that correspond to desired tool positions andorientations at fixed displacement intervals along the toolpath. The position commands arethen scheduled with respect to time, using a non-linear optimization algorithm to ensure thatthe toolpath is traversed as fast as possible without violating the velocity, acceleration, and jerkconstraints.In addition to a trajectory generation algorithm, a control strategy was developed which com-bines the long stroke of the 3-axis micromill with the high bandwidth tracking capabilitiesof the 6DOF rotary table. Prior to the control design, the rigid body dynamics of the mi-cromill were identified using a linear regression technique. The model is further refined, par-ticularly the non-linear friction characteristics, using disturbance observations with a Kalmanfilter. Based upon the rigid body model, a sliding mode controller with a bandwidth of approx-imately 30Hz was designed for position tracking. In order to improve tracking performance,a feedforward friction compensator, based on the LuGre friction model, is implemented. Dueto its comparatively simpler dynamics, sophisticated identification was not necessary for therotary table. Instead, the position controller was designed around a nominal model of the rotarytable, which consisted only of a free floating mass. To control position, a lead-lag controllerwas selected, cascaded with an integrator to remove steady state error. However, it was foundthat unmodeled flexibilities were limiting the potential bandwidth of the rotary table. Fre-104quency sweeps of the plant found resonance peaks at higher frequencies and notch filters oneach axis were implemented, resulting in an achievable bandwidth of approximately 250Hz. Inorder to combine the long stroke capabilities of the micromill with the high bandwidth track-ing capabilities of the rotary table, the tracking error of the micromill is sent as the referencecommand of the rotary table. Analysis of the transfer function showed that this configurationallowed the axis to adopt the bandwidth of the rotary table. Furthermore, results showed thatthe rotary table successfully compensated for tracking errors caused by the lower bandwidth ofthe micromill’s three Cartesian drives.Finally, the rotary table was used to compensate for geometric errors of the machine tool. Inorder to do this, the effect of the geometric errors of the 3-axis machine tool on the tooltipposition had to be modeled. The ideal transformation matrix of each moving axis is augmentedto account for the effect of six geometric errors, including positioning error, two straightnesserrors, roll error, pitch error, and yaw error. The errors are measured with a laser interferometerand fit to a quintic polynomial function of position to preserve jerk continuity, which are thenincorporated into the error augmented transformation matrices. The original forward kinematicmodel is reconstructed with the error augmented transformation matrices, resulting in a tooltipposition that accounts for the geometric errors of the 3-axis of the machine tool. By subtractingthe result of this forward kinematic model with the original ideal forward kinematic model, itis possible to model the tooltip deviations. Next, rotary table commands are generated tocompensate for these tooltip deviations. Due to the non-Cartesian kinematics of rotationalerrors and the presence of more degrees of freedom than necessary, a non-linear technique wasdeveloped in order to generate compensating commands for the rotary table. A gradient descentoptimization scheme was developed where the goal was the minimization of tooltip deviations.Since there were more degrees of freedom than necessary, the Moore-Penrose inverse wasused again. However, since only the rotary table commands were modified, the Jacobian in105this case only accounted for differential changes in tooltip position with respect to changes inrotary table positions. Furthermore, the solution is augmented to ensure the stroke limits of therotary table were not violated. Single axis laser interferometer experiments showed significantreductions in geometric errors as the combined effort of all 6 axes of the rotary table. Multi-axissimulations results also showed that tooltip deviations were minimized for multi-axis free-formtoolpaths.In summary, the thesis presents a novel 9-axis CNC micro machine tool with a new trajectorygeneration algorithm, dual axis control algorithm and geometric error compensation strategy.The proposed models can be applied to other multi-axes machine tools with redundant axes.7.2 Future Research DirectionsIn regards to the trajectory algorithm presented, further work can be done to select configu-rations that take advantage of the redundant degrees of freedom. At present, the trajectorygeneration algorithm only uses the nullspace of the Jacobian to select joint configurations thatavoid stroke limits. In reality, it may be possible to achieve multiple goals simultaneously,particularly with so many redundant degrees of freedom. The algorithm presented could beextended to fulfill multiple goals. To name a few, minimization of energy consumption or jointtorques can be considered. Specific to the machining process, if it is possible to model therelative stiffness between the tool and the workpiece as a function of the axes positions, thengiven an analytic gradient, it would be possible to select axes configurations in which the thestiffness is maximized. With the presence of multiple goals, it should be noted that additionalstrategies need to be developed in prioritizing the goals. With respect to the main redundancyresolution aspect of the trajectory generation algorithm, the Moore-Penrose inverse is used,which consequently minimizes the axes differential with respect to displacement locally. It hasbeen shown in literature that higher order versions of the presented solutions minimize the axes106differential at the respective orders. By performing the redundancy resolution at higher order itmay be possible to achieve more desirable acceleration and jerk responses when scheduling thecommands with respect to time. Finally, globally optimal variations of the presented solutionhave also been shown in literature and may be an interesting area of research and applicationfor the configuration presented in this work. However it should be noted that due to the glob-ally optimal criteria, even short toolpaths require unrealistic computation time, and optimizingfor a real free-form CNC toolpath may be outside the realm of practicality.Due to its relatively unique configuration, many research directions in the field of controlscould be pursued. In a gantry machine tool, one of the lowest frequency structural modesis from the structure of the z-axis itself. The rotary table could be used to provide activedamping for this structural mode, increasing the dynamic stiffness between the work pieceand the tool. Unlike traditional active damping devices, the rotary table could actively dampstructural modes in multiple directions simultaneously. It should be noted that the dual stagefeed drive control algorithm presented considers the micromill and rotary table as two separaterigid bodies. This is particularly true for the rotary table acting on the heavier feed drive.In reality, there may be flexibilities that when considered, could increase the precision of theoverall system. As a result, state-space control laws which consider the dynamic couplingbetween the coarse and fine actuator, with the purpose of increasing performance or robustnesswould be interesting to consider. With respect to geometric error compensation, there would bevalue in verifying the compensation algorithm for multi-axis trajectories. Though the algorithmhas shown effectiveness for single axis experimental results and multi-axis simulation results,verification via a ball-bar would be a possible research direction. Furthermore, the rotarytable can also be used for compensation of errors beyond geometric error. If it is possibleto model errors caused by thermal expansion or process forces, the tooltip deviations can bemodeled. These deviations can be incorporated into the compensation algorithm presented,107resulting in rotary table commands that compensate for the above mentioned sources of error.Alternatively, a simpler experimental approach could be taken. A part could be machined andthe tooltip deviations could be measured with a CMM. These deviations could be sent to thesame compensation algorithm which would generate compensating commands based on CMMmeasurements. This would require no modeling but a large amount of experimental data.Finally, in addition to further developing the presented algorithms, further research can be doneby applying the presented algorithms to different configurations. Since robotic arms typicallyhave more degrees of freedom than necessary, application of the trajectory generation algorithmpresented in this work to a different configuration could be considered novel. Furthermore,since the presented algorithm has an optimization aspect, it would make the most sense to applythis algorithm to robotic milling arms. Likewise, the dual stage feed drive control algorithmscan be applied to more coarse/fine actuator configurations found in manufacturing literature.In contrast to using the rotary table as an actuation device, it may be interesting to use therotary table as a sensing device. If the rotary table is commanded to be held stationary, then theoutput current to hold the rotary table would be proportional to the cutting force. If the sameconfiguration is used to machine a part, the feedrate of the overall system could be controlledto ensure a constant cutting force is held at all times, using the rotary table as force feedback.108Bibliography[1] X. Lu, M. Dyck, and Y. Altintas, “Magnetically levitated six degree of freedom rotarytable,” CIRP Annals, vol. 64, pp. 353–356, Jan. 2015.[2] M. Dyck, X. Lu, and Y. Altintas, “Magnetically Levitated Rotary Table With Six De-grees of Freedom,” IEEE/ASME Transactions on Mechatronics, vol. 22, pp. 530–540,Feb. 2017.[3] R.-S. Lee and C.-H. She, “Developing a postprocessor for three types of five-axis ma-chine tools,” The International Journal of Advanced Manufacturing Technology, vol. 13,pp. 658–665, Sept. 1997.[4] Psang Dain Lin and Ing Jyh Tsai, “The machining and on-line measurement of spatialcams on four-axis machine tools,” International Journal of Machine Tools and Manufac-ture, vol. 36, pp. 89–101, Jan. 1996.[5] B. Sencer, Y. Altintas, and E. Croft, “Feed optimization for five-axis CNC machine toolswith drive constraints,” International Journal of Machine Tools and Manufacture, vol. 48,pp. 733–745, June 2008.[6] D. E. Whitney, “Resolved Motion Rate Control of Manipulators and Human Prostheses,”IEEE Transactions on Man-Machine Systems, vol. 10, pp. 47–53, June 1969.[7] A. Liegeois, “Automatic Supervisory Control of the Configuration and Behavior ofMultibody Mechanisms,” IEEE Transactions on Systems, Man, and Cybernetics, vol. 7,pp. 868–871, Dec. 1977.109[8] J. Hollerbach and K. Suh, “Redundancy resolution of manipulators through torque opti-mization,” IEEE Journal on Robotics and Automation, vol. 3, pp. 308–316, Aug. 1987.[9] Y. Halevi, E. Carpanzano, G. Montalbano, and Y. Koren, “Minimum energy control ofredundant actuation machine tools,” CIRP Annals - Manufacturing Technology, vol. 60,no. 1, pp. 433–436, 2011.[10] Y. Nakamura and H. Hanafusa, “Inverse Kinematic Solutions With Singularity Robust-ness for Robot Manipulator Control,” Journal of Dynamic Systems, Measurement, andControl, vol. 108, pp. 163–171, Sept. 1986.[11] Y. Nakamura and H. Hanafusa, “Optimal Redundancy Control of Robot Manipulators,”The International Journal of Robotics Research, vol. 6, pp. 32–42, Mar. 1987.[12] E. Sanyildiz and H. Temelta, “A comparison study of the numerical integration meth-ods in the trajectory tracking application of redundant robot manipulators,” in 2011 7thInternational Conference on Electrical and Electronics Engineering (ELECO), pp. II–420–II–424, Dec. 2011.[13] L. Sciavicco and B. Siciliano, “A solution algorithm to the inverse kinematic problem forredundant manipulators,” IEEE Journal on Robotics and Automation, vol. 4, pp. 403–410,Aug. 1988.[14] K. J. Kyriakopoulos and G. N. Saridis, “Minimum jerk path generation,” in 1988 IEEEInternational Conference on Robotics and Automation Proceedings, pp. 364–369 vol.1,Apr. 1988.[15] A. Gasparetto and V. Zanotto, “Optimal trajectory planning for industrial robots,” Ad-vances in Engineering Software, vol. 41, pp. 548–556, Apr. 2010.110[16] J. Dong, P. M. Ferreira, and J. A. Stori, “Feed-rate optimization with jerk constraintsfor generating minimum-time trajectories,” International Journal of Machine Tools andManufacture, vol. 47, pp. 1941–1955, Oct. 2007.[17] X. Beudaert, S. Lavernhe, and C. Tournier, “Feedrate interpolation with axis jerk con-straints on 5-axis NURBS and G1 tool path,” International Journal of Machine Tools andManufacture, vol. 57, pp. 73–82, June 2012.[18] Y. Altintas and K. Erkorkmaz, “Feedrate optimization for spline interpolation in highspeed machine tools,” CIRP Annals, vol. 52, pp. 297–302, Jan. 2003.[19] K. Erkorkmaz and Y. Altintas, “High speed CNC system design. Part III: high speedtracking and contouring control of feed drives,” International Journal of Machine Toolsand Manufacture, vol. 41, pp. 1637–1658, Sept. 2001.[20] M. Tomizuka, “Zero phase error tracking algorithm for digital control,” Journal of Dy-namic Systems, Measurement, and Control, vol. 109, pp. 65–68, Mar. 1987.[21] Y. Altintas, K. Erkorkmaz, and W. H. Zhu, “Sliding mode controller design for high speedfeed drives,” CIRP Annals, vol. 49, pp. 265–270, Jan. 2000.[22] C. Okwudire and Y. Altintas, “Minimum tracking error control of flexible ball screwdrives using a discrete-time sliding mode controller,” Journal of Dynamic Systems, Mea-surement, and Control, vol. 131, pp. 051006–051006–12, Aug. 2009.[23] M. Hanifzadegan and R. Nagamune, “Switching gain-scheduled control design for flexi-ble ball-screw drives,” Journal of Dynamic Systems, Measurement, and Control, vol. 136,pp. 014503–014503–6, Sept. 2013.111[24] A. H. H. Hosseinabadi and Y. Altintas, “Modeling and active damping of structural vibra-tions in machine tools,” CIRP Journal of Manufacturing Science and Technology, vol. 7,pp. 246–257, Jan. 2014.[25] A. Kamalzadeh and K. Erkorkmaz, “Compensation of Axial Vibrations in Ball ScrewDrives,” CIRP Annals, vol. 56, pp. 373–378, Jan. 2007.[26] K. Erkorkmaz and A. Kamalzadeh, “High Bandwidth Control of Ball Screw Drives,”CIRP Annals, vol. 55, pp. 393–398, Jan. 2006.[27] A. Kamalzadeh and K. Erkorkmaz, “Accurate tracking controller design for high-speeddrives,” International Journal of Machine Tools and Manufacture, vol. 47, pp. 1393–1400, July 2007.[28] S. Staroselsky and K. A. Stelson, “Two-Stage Actuation for Improved Accuracy of Con-touring,” in 1988 American Control Conference, pp. 127–132, June 1988.[29] O. Masory and J. Wang, “Two-Stage Actuation Control System for CNC ContouringSystems,” in 1992 American Control Conference, pp. 375–378, June 1992.[30] A. T. Elfizy, G. M. Bone, and M. A. Elbestawi, “Design and control of a dual-stage feeddrive,” International Journal of Machine Tools and Manufacture, vol. 45, pp. 153–165,Feb. 2005.[31] Y. M. Choi and D. G. Gweon, “A High-Precision Dual-Servo Stage Using Halbach LinearActive Magnetic Bearings,” IEEE/ASME Transactions on Mechatronics, vol. 16, pp. 925–931, Oct. 2011.[32] B.-S. Kim, J. Li, and T.-C. Tsao, “Two-parameter robust repetitive control with applica-tion to a novel dual-stage actuator for noncircular Machining,” IEEE/ASME Transactionson Mechatronics, vol. 9, pp. 644–652, Dec. 2004.112[33] M. Kobayashi and R. Horowitz, “Track seek control for hard disk dual-stage servo sys-tems,” IEEE Transactions on Magnetics, vol. 37, pp. 949–954, Mar. 2001.[34] M. Kobayashi, S. Nakagawa, and S. Nakamura, “A phase-stabilized servo controllerfor dual-stage actuators in hard disk drives,” IEEE Transactions on Magnetics, vol. 39,pp. 844–850, Mar. 2003.[35] G. Herrmann, M. C. Turner, I. Postlethwaite, and G. Guo, “Practical implementation of anovel anti-windup scheme in a HDD-dual-stage servo-system,” IEEE/ASME Transactionson Mechatronics, vol. 9, pp. 580–592, Sept. 2004.[36] W. H. Yao and M. Tomizuka, “Robust controller design for a dual-stage positioning sys-tem,” in Proceedings of the IECON ’93 International Conference on Industrial Electron-ics, Control, and Instrumentation, pp. 62–66 vol.1, Nov. 1993.[37] T.-L. Tai and J.-S. Chen, “Discrete-time sliding-mode controller for dual-stage systemsAhierarchical approach,” Mechatronics, vol. 15, pp. 949–967, Oct. 2005.[38] J. H. She, X. Xin, and Y. Pan, “Equivalent-Input-Disturbance Approach-Analysisand Application to Disturbance Rejection in Dual-Stage Feed Drive Control System,”IEEE/ASME Transactions on Mechatronics, vol. 16, pp. 330–340, Apr. 2011.[39] J. Zheng and M. Fu, “Nonlinear Feedback Control of a Dual-Stage Actuator System forReduced Settling Time,” IEEE Transactions on Control Systems Technology, vol. 16,pp. 717–725, July 2008.[40] G. Zhang, R. Veale, T. Charlton, B. Borchardt, and R. Hocken, “Error Compensation ofCoordinate Measuring Machines,” CIRP Annals, vol. 34, pp. 445–448, Jan. 1985.113[41] A. C. Okafor and Y. M. Ertekin, “Derivation of machine tool error models and error com-pensation procedure for three axes vertical machining center using rigid body kinemat-ics,” International Journal of Machine Tools and Manufacture, vol. 40, pp. 1199–1213,June 2000.[42] W. T. Lei and Y. Y. Hsu, “Accuracy enhancement of five-axis CNC machines throughreal-time error compensation,” International Journal of Machine Tools and Manufacture,vol. 43, pp. 871–877, July 2003.[43] N. Huang, Y. Jin, Q. Bi, and Y. Wang, “Integrated post-processor for 5-axis machinetools with geometric errors compensation,” International Journal of Machine Tools andManufacture, vol. 94, pp. 65–73, July 2015.[44] Y. Y. Hsu and S. S. Wang, “A new compensation method for geometry errors of five-axis machine tools,” International Journal of Machine Tools and Manufacture, vol. 47,pp. 352–360, Feb. 2007.[45] S. Zhu, G. Ding, S. Qin, J. Lei, L. Zhuang, and K. Yan, “Integrated geometric errormodeling, identification and compensation of CNC machine tools,” International Journalof Machine Tools and Manufacture, vol. 52, pp. 24–29, Jan. 2012.[46] M. Tsutsumi, S. Tone, N. Kato, and R. Sato, “Enhancement of geometric accuracy offive-axis machining centers based on identification and compensation of geometric de-viations,” International Journal of Machine Tools and Manufacture, vol. 68, pp. 11–20,May 2013.[47] S. Xiang and Y. Altintas, “Modeling and compensation of volumetric errors for five-axis machine tools,” International Journal of Machine Tools and Manufacture, vol. 101,pp. 65–78, Feb. 2016.114[48] O.-S. Kim, S.-H. Lee, and D.-C. Han, “Positioning performance and straightness errorcompensation of the magnetic levitation stage supported by the linear magnetic bearing,”IEEE Transactions on Industrial Electronics, vol. 50, pp. 374–378, Apr. 2003.[49] Y. Deng, X. Jin, and Z. Zhang, “A macromicro compensation method for straightnessmotion error and positioning error of an improved linear stage,” The International Journalof Advanced Manufacturing Technology, vol. 80, pp. 1799–1806, Oct. 2015.[50] R. P. Paul, Robot manipulators: mathematics, programming, and control : the computercontrol of robot manipulators. MIT Press, 1981. Google-Books-ID: m99TAAAAMAAJ.[51] A. Yuen, K. Zhang, and Y. Altintas, “Smooth trajectory generation for five-axis machinetools,” International Journal of Machine Tools and Manufacture, vol. 71, pp. 11–19, Aug.2013.[52] M. Heng and K. Erkorkmaz, “Design of a NURBS interpolator with minimal feed fluctu-ation and continuous feed modulation capability,” International Journal of Machine Toolsand Manufacture, vol. 50, pp. 281–293, Mar. 2010.[53] L. Piegl and W. Tiller, The NURBS Book. Monographs in Visual Communication, BerlinHeidelberg: Springer-Verlag, 1995.[54] K. Erkorkmaz and Y. Altintas, “High speed CNC system design. Part II: modeling andidentification of feed drives,” International Journal of Machine Tools and Manufacture,vol. 41, pp. 1487–1509, Aug. 2001.[55] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” Journalof Basic Engineering, vol. 82, pp. 35–45, Mar. 1960.[56] V. I. Utkin, “Sliding mode control design principles and applications to electric drives,”IEEE Transactions on Industrial Electronics, vol. 40, pp. 23–36, Feb. 1993.115[57] C. C. De Wit, H. Olsson, K. J. Astrom, and P. Lischinsky, “A new model for control ofsystems with friction,” IEEE Transactions on Automatic Control, vol. 40, no. 3, pp. 419–425, 1995.[58] P. I. Ro, W. Shim, and S. Jeong, “Robust friction compensation for submicrometer posi-tioning and tracking for a ball-screw-driven slide system,” Precision Engineering, vol. 24,no. 2, pp. 160–173, 2000.[59] M. Dyck, “Magnetically levitated six degree of freedom micro-machining rotary table,”Master’s thesis, University of British Columbia, 2014.116Appendix AFoward Kinematic Equation with Geometric ErrorsThe actual tooltip position, Pa (q) = [Px,a (q) , Py,a (q) , Pz,a (q)]T are shown detail in thefollowing equations:Px,a (q) = −εx(xc)((L+ Lx,z)(caf scf + ccf saf sbf )+(Lsp,y − yc)(saf scf − caf ccf sbf ))+εx(zc)((Lsp,z + Lt)(caf scf + ccf saf sbf )+Lsp,y(saf scf − caf ccf sbf ))+εy(xc)((L+ Lx,z)cbf ccf + xc(caf ccf sbf − saf scf ))+εy(yc)((L)cbf ccf )+εz(xc)((Lsp,y − yc)cbf ccf + xc(caf scf + ccf saf sbf ))−εy(zc)((Lsp,z + Lt)cbf ccf )−εx(yc)(L(caf scf + ccf saf sbf )+(Lsp,y − yc)(saf scf − caf ccf sbf ))+εz(yc)((Lsp,y − yc)cbf ccf )+(−δy(xc)− δy(yc) + δy(zc))(caf scf + ccf saf sbf )+(−δz(xc)− δz(yc) + δz(zc))(saf scf − caf ccf sbf )+(Lsp,y − yc − yf )caf scf+(−δx(xc)− δx(yc) + δx(zc)−xc − xf − Lsp,yεz(zc))cbf ccf+(−L− Ls − Lx,z − zf − zf,0)saf scf+(L+ Ls + Lx,z + zf + zf,0)caf ccf sbf+(Lsp,y − yc − yf )ccf saf sbf−γx,zc((−L+ Ly,z)(caf scf + ccf saf sbf )+Lsp,y(caf ccf sbf − saf scf ))+γy,zc((−L+ Ly,z)cbf ccf )+γz,xc((Lsp,y − yc)cbf ccf )+O(ε2, δ2, εδ)(A.1)117Py,a (q) = εx(xc)((L+ Lx,z)(saf sbf scf − caf ccf )+(−Lsp,y + yc)(ccf saf + caf sbf scf ))+εx(zc)((Lsp,z + Lt)(caf ccf − saf sbf scf )+Lsp,y(ccf saf + caf sbf scf ))+εx(yc)(L(saf sbf scf − caf ccf )+(−Lsp,y + yc)(ccf saf + caf sbf scf ))−εz(xc)((Lsp,y − yc)cbf scf + xc(saf sbf scf − caf ccf ))+εy(zc)((Lsp,z + Lt)cbf scf )−εy(yc)(Lcbf scf )−εy(xc)((L+ Lx,z)cbf scf + xc(ccf saf + caf sbf scf ))−εz(yc)((Lsp,y − yc)cbf scf )+(−δy(xc)− δy(yc) + δy(zc))(caf ccf − saf sbf scf )+(−δz(xc)− δz(yc) + δz(zc))(ccf saf + caf sbf scf )+(Lsp,y − yc − yf )caf ccf+(−L− Ls − Lx,z − zf − zf,0)ccf saf+(δx(xc) + δx(yc)− δx(zc)+xc + xf + Lsp,yεz(zc))cbf scf+(−L− Ls − Lx,z − zf − zf,0)caf sbf scf+(−Lsp,y + yc + yf )saf sbf scf−γz,xc((Lsp,y − yc)cbf scf )−γy,zc((−L+ Ly,z)cbf scf )+γx,zc((L− Ly,z)(caf ccf − saf sbf scf )+Lsp,y(ccf saf + caf sbf scf ))+O(ε2, δ2, εδ)(A.2)118Pz,a (q) = εx(yc)(Lcbf saf + (−Lsp,y + yc)caf cbf )−Lw − εy(zc)(Lsp,zsbf + Ltsbf )−Lr + εz(yc)(Lsp,ysbf − ycsbf )+εx(xc)((L+ Lx,z)cbf saf + (−Lsp,y + yc)caf cbf )+(−δx(xc)− δx(yc) + δx(zc)− Lsp,yεz(zc))sbf+(−xc − xf )sbf+εy(xc)((L+ Lx,z)sbf − xccaf cbf )−εz(xc)((yc − Lsp,y)sbf + xccbf saf )+εy(yc)((L)sbf )−εx(zc)((Lsp,z + Lt)cbf saf − Lsp,ycaf cbf )+(−L− Ls − Lx,z − δz(xc)− δz(yc)+δz(zc)− zf − zf,0)caf cbf+(−Lsp,y + δy(xc) + δy(yc)−δy(zc) + yc + yf )cbf saf+γz,xc((Lsp,y − yc)sbf )+γy,zc((−L+ Ly,z)sbf )+γx,zc(Lsp,ycaf cbf + (−L+ Ly,z)cbf saf )+O(ε2, δ2, εδ)(A.3)where L = Lsp,z − L0,z + Lt + Lx,z − zc and O(ε2, δ2, εδ) are the second order and higherorder error terms.119**