Automated Design and Implementation of Kalman Observer for Spindle Torque Estimation in CNC Machining by Zhao Wei Lu BASc, The University of British Columbia, 2018 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2020 © Zhao Wei Lu, 2020 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis entitled: Automated Design and Implementation of Kalman Observer for Spindle Torque Estimation in CNC Machining submitted by Zhao Wei Lu in partial fulfillment of the requirements for the degree of Master of Applied Science in Mechanical Engineering Examining Committee: Dr. Yusuf Altintas, Mechanical Engineering Supervisor Dr. Ryozo Nagamune, Mechanical Engineering Supervisory Committee Member Dr. Minkyun Noh, Mechanical Engineering Supervisory Committee Member iii Abstract Milling is a subtractive manufacturing technique, where the material is continuously removed from a workpiece until the desired part shape is obtained. Heavily used in the automotive and aerospace industries, Computer Numerical Control (CNC) milling machines occupy a large chunk in the manufacturing process. The current research focuses towards machining and machine tool monitoring systems that are more self-sufficient and self-adjusting to adapt the processes to machine tools. Cutting torque delivered by the machine tool spindle is one of the key sensory signals for machining process monitoring. This thesis presents a method that automatically reconstructs cutting torque from motor current commands generated by the servo controller of the machine tool. To estimate the cutting torque from commanded spindle current, the dynamics between torque to current relationship must be modeled and compensated. The thesis first presents an automated identification of spindle dynamics using data-driven system identification methods. The frequency response function (FRF) of the spindle dynamics is measured manually using CNC internal diagnostic tools. The identified FRF is then automatically converted to a state-space model using the Eigensystem Realization Algorithm (ERA). To reduce overfitting, an optimal threshold is applied to the ERA method to limit the identified system order. And to ensure the stability of the identified system, unstable eigenvalues of the system are removed using Schur decomposition. The identified system is then augmented such that the unknown torque input is modeled as a state changed by a random process noise. This augmented system is used to create a Kalman Observer, which compensates the spindle dynamics and estimates the torque from the spindle nominal iv current. The Kalman Observer is tuned automatically by estimating the noise covariance values using machining simulations. The method was eventually validated on a Quaser UX 600 industrial CNC system. v Lay Summary Milling is a process for which material is successively removed from a workpiece until the desired part geometry is achieved. This process can create high precision parts with surface finish qualities desirable to be used in aerospace, automotive, and die/cast industries. CNC (Computer Numerical Control) Milling employs control techniques to automate the execution of machine motions to cut the part, and the current research is focused on efficient self-reliant milling processes. The desire is to make CNC machines be able to self-diagnose machining problems and adjust the operating conditions accordingly. This is typically accomplished by measuring the cutting process states, which often require costly sensors or complicated dynamics models that require expert input. In this thesis, a sensor-less and automated solution to cutting process measurement is presented. vi Preface This thesis is the original, unpublished, independent work by the author, Zhaowei Lu. This work was supervised by Dr. Yusuf Altintas and based on the previous work done by former Ph.D. student Deniz Aslan and Yusuf Altintas. vii Table of Contents Abstract ...................................................................................................................................... iii Lay Summary .............................................................................................................................. v Preface........................................................................................................................................ vi Table of Contents ...................................................................................................................... vii List of Tables ............................................................................................................................. ix List of Figures ............................................................................................................................. x List of Symbols ........................................................................................................................ xiii List of Abbreviations ............................................................................................................... xvi Acknowledgments................................................................................................................... xvii Chapter 1: Introduction ............................................................................................................... 1 Chapter 2: Literature Review ...................................................................................................... 5 2.1 Overview ..................................................................................................................... 5 2.2 Direct Measurement Approaches ................................................................................ 5 2.3 Indirect Measurement Approaches ............................................................................. 6 2.4 Discussion ................................................................................................................... 9 Chapter 3: Automated Spindle Drive System Identification .................................................... 12 3.1 Spindle Disturbance Impulse Response Measurement ............................................. 12 3.2 Fitting to State-Space System ................................................................................... 15 3.3 Validation .................................................................................................................. 20 3.4 Discussion ................................................................................................................. 35 viii Chapter 4: Automated Kalman Filter Design ........................................................................... 36 4.1 State Space Transformation ...................................................................................... 36 4.2 Kalman Observer Design .......................................................................................... 39 4.3 Kalman Observer Automatic Tuning ........................................................................ 41 4.4 Results ....................................................................................................................... 47 4.5 Discussion ................................................................................................................. 62 Chapter 5: Conclusion............................................................................................................... 64 References ................................................................................................................................. 67 Appendix ................................................................................................................................... 71 ix List of Tables Table 1-1 Direct force/torque measuring methods ......................................................................... 9 Table 1-2 Indirect sensor based force/torque measuring methods .................................................. 9 Table 1-3 Indirect sensorless force/torque measuring methods .................................................... 10 Table 1-1 Kalman Observer Parameters ....................................................................................... 55 Table A.1-1 First Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 1-14 .. 72 Table A.1-2 First Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 14-28 73 Table A.1-3 Second Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 1-11....................................................................................................................................................... 74 Table A-1-4 Second Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 11-22....................................................................................................................................................... 75 x List of Figures Figure 1.1 Milling cutting schematic .............................................................................................. 2 Figure 1.2 Spindle disturbance torque response schematic ............................................................ 3 Figure 3.1: Simplified CNC Spindle Velocity Control Loop Block Diagram .............................. 14 Figure 3.2: Spindle disturbance torque to nominal current FRF at the first spindle speed range . 22 Figure 3.3 Spindle disturbance torque to nominal current FRF at the second spindle speed range....................................................................................................................................................... 23 Figure 3.4: Spindle impact test setup ............................................................................................ 24 Figure 3.5: Comparison of torque-current spindle FRF obtained from sine sweep and tap test .. 24 Figure 3.6: Steps for generating noisy simulated FRFs ................................................................ 26 Figure 3.7: Simulated Noisy FRF of 𝝉𝒅𝒋𝝎𝑰𝒏𝒐𝒎𝒋𝝎 for 10 State System with signal to noise ratio SNR = 10 dB ........................................................................................................................ 27 Figure 3.8: Fit results of 𝝉𝒅𝒋𝝎𝑰𝒏𝒐𝒎𝒋𝝎 for 10 State System with SNR = 10 dB with Normalized Root Mean Square NRMSE = 95.3 % ...................................................................... 28 Figure 3.9 Fit results of 𝝉𝒅𝒋𝝎𝑰𝒏𝒐𝒎𝒋𝝎 for 20 State System with SNR = 10 dB with NMRSE = 93.1 % ........................................................................................................................................... 28 Figure 3.10 Fit results of 𝝉𝒅𝒋𝝎𝑰𝒏𝒐𝒎𝒋𝝎 for 30 State System with SNR = 10 dB with NRMSE = 94.6 % ........................................................................................................................................... 29 Figure 3.11 Simulated Noisy FRF of 𝝉𝒅𝒋𝝎𝑰𝒏𝒐𝒎𝒋𝝎 for 10 State System with SNR = 0 dB ... 30 Figure 3.12 Fit results of 𝝉𝒅𝒋𝝎𝑰𝒏𝒐𝒎𝒋𝝎 for 10 State System with SNR = 0 dB with NRMSE = 95.8% ............................................................................................................................................ 30 Figure 3.13 Comparison between MATLAB ‘modalfit’ and ERA .............................................. 31 xi Figure 3.14 Fit result from measured first speed range torque-current spindle FRF with NRMSE = 83.2 % ........................................................................................................................................ 32 Figure 3.15: Fit result from measured second-speed range torque-current spindle FRF with NRMSE = 80.5 % ......................................................................................................................... 33 Figure 3.16 ERA Fit phase comparison for first spindle speed range .......................................... 34 Figure 3.17 ERA Fit phase comparison for second spindle speed range...................................... 34 Figure 4.1 State space system schematic ...................................................................................... 37 Figure 4.2 Expanded state-space model schematic....................................................................... 38 Figure 4.3 Kalman Observer Schematic ....................................................................................... 41 Figure 4.4 Uniform Distribution ................................................................................................... 42 Figure 4.5 Endmill Cutting Chip Thickness Diagram .................................................................. 44 Figure 4.6 Endmill Helix Angle Diagram..................................................................................... 45 Figure 4.7 Rotary dynamometer measurement dynamics FRF .................................................... 48 Figure 4.8 Rotary dynamometer FRF measurement input spectrum ............................................ 49 Figure 4.9 Simulated Cutting Torque ........................................................................................... 50 Figure 4.10 Simulated Cutting Torque Change per Sampling Interval ........................................ 51 Figure 4.11 Kalman Observer Estimated Torque (3000 RPM spindle speed) ............................. 53 Figure 4.12 Kalman Observer Estimated Torque (5000 RPM spindle speed) ............................. 56 Figure 4.13 Kalman Observer Estimated Torque (6000 RPM spindle speed) ............................. 57 Figure 4.14 Kalman Observer Estimated Torque (7000 RPM spindle speed) ............................. 58 Figure 4.15 Manual Kalman Observer Estimated Torque (3000 RPM Spindle Speed) ............... 60 Figure 4.16 Manual Kalman Observer Estimated Torque (5000 RPM Spindle Speed) ............... 60 Figure 4.17 Manual Kalman Observer Estimated Torque (6000 RPM Spindle Speed) ............... 61 xii Figure 4.18Manual Kalman Observer Estimated Torque (7000 RPM Spindle Speed) ................ 61 Figure A0.1 Real and Imaginary value comparision between fitted and measured FRF of First Spindle Speed Range .................................................................................................................... 77 Figure A0.2 Real and Imaginary value comparision between fitted and measured FRF of Second Spindle Speed Range .................................................................................................................... 78 xiii List of Symbols 𝑎, 𝑏 Upper and lower limits of the measurement resolution 𝑨,𝑩, 𝑪,𝑫 Current-torque dynamics state-space matrices ?̃?, 𝑩,̃ ?̃? Minimum realization of the identified state-space system ?̃?𝑑𝑒𝑐 , ?̃?𝑑𝑒𝑐 , ?̃?𝑑𝑒𝑐 Stable-unstable decomposed realization of the identified state-space system 𝑨𝒆𝒙𝒑, 𝑩𝒆𝒙𝒑, 𝑪𝒆𝒙𝒑 Expanded realization of the identified state-space system 𝒆 State estimation error 𝐺𝑜 Spindle velocity open loop dynamics 𝐺𝑙𝑜𝑜𝑝 Spindle velocity control loop dynamics 𝐺𝑠 Spindle mechanical structure dynamics 𝐺Ω Spindle velocity controller dynamics 𝐺𝑑𝜏 Spindle disturbance torque to current dynamics ℎ Chip thickness 𝐇 Impulse response Hankel matrix 𝐇′ Time shifted impulse response Hankel matrix 𝐼𝑎𝑐𝑡 Actual motor current 𝐼𝑛𝑜𝑚 Nominal spindle current 𝑗 Complex number operator 𝐽𝑒 Spindle inertia 𝐾𝑡 Motor torque constant 𝐾𝑡𝑐 Cutting force coefficient xiv 𝐾𝑡𝑒 Cutting edge coefficient 𝑲𝒔𝒔 Steady state Kalman Observer gain 𝑁 Number of teeth in cutting tool 𝑷 State estimation error covariance matrix 𝑷𝒌+𝟏′ A priori estimate of the state estimation error covariance 𝑸 Process noise covariance matrix 𝑹 Measure noise covariance 𝑼,𝑽 Eigen-time-delay co-ordinates in impulse response 𝑣 Measurement noise 𝑤 Process noise 𝒙 Current-torque dynamics state-space vector ?̂?𝒆𝒙𝒑𝒌+𝟏′ A priori estimate of the state vector 𝑦 Measured nominal current 𝜏𝑑 Disturbance torque Ω𝑎𝑐𝑡 Actual spindle velocity Ω𝑛𝑜𝑚 Nominal spindle velocity 𝓞 Observability matrix 𝓒 Controllability matrix 𝚺 Singular values matrix 𝓞 ̃, ?̃? Minimum realization of the observability and controllability matrices 𝜎2 Variance 𝜙𝑐 Tool instantaneous engagement angle xv 𝛽 Tool helix angle 𝜙𝑐𝛽 Helix lag angle 𝜙𝑝 Tool pitch angle xvi List of Abbreviations CNC Computer numerical control DARE Discrete algebraic Riccati equation FRF Frequency response function FFT Fast Fourier transform NRMSE Normalized root mean squared error SNR Signal to Noise Ratio SVD Singular value decomposition xvii Acknowledgments I would like to first and foremost express my gratitude to my supervisor Prof. Yusuf Altintas for his continued support and feedback throughout my research degree. I sincerely appreciate the opportunity that was given to me to study and be a part of the Manufacturing Automation Laboratory (MAL) and will treasure my time spent here. I would also like to thank Dr. Deniz Aslan and Dr. Hoai Nam Huynh for their discussions and feedback regarding my thesis. Their help has been invaluable in helping me navigate through the challenging parts of my degree. My gratitude also goes towards all my friends and peers at MAL, who made my journey through my master’s degree a thoroughly enjoyable one filled with laughter and joy. Special thanks to the NSERC CANRIMT for the financial support and the Industrial Technology Research Institute of Taiwan for providing the Quaser UX 600 CNC machining center. Finally, I would like to thank my family and loved ones for all the unconditional love and support that was given to me during my thesis. To my mother and father, thank you for believing in what I could accomplish, and encouraging me to push myself. To May, thank you for your unconditional support and optimism. I hope to one day be able to repay all the love and support that was given to me.1 Chapter 1: Introduction Milling is a subtractive manufacturing process in which material is successively removed from a workpiece using a rotating cutter. This machining process allows the creation of parts used in aerospace, automotive, and die/mold industries. Recently, there has been a push towards creating smarter and self-sufficient Computer Numerical Control (CNC) machine tools. Typically, the process and machine tool monitoring systems are operated based on the measure of several process states such as cutting force, vibrations, sound, and temperature. These process states are measured with machine-mounted sensors or internal machine control signals and used in applications such as machine health monitoring, chatter detection, tool breakage detection, and collision detection. In literature, the measured cutting forces and torque have been used extensively in monitoring machining operations since they can be directly correlated to the state of the milling process (Figure 1.1) directly [1]. 2 Figure 1.1 Milling cutting schematic The cutting torque generated during milling is directly transmitted through the spindle and tool structure, and finally to the spindle motor, which experiences it as a disturbance torque as shown in Figure 1.2. This disturbance torque causes variations in the spindle speed from the nominal speed, which is measured by the spindle encoder. Subsequently, the spindle’s servo controller suppresses the effects of the disturbance torque by generating extra current to the motor. By measuring the commanded spindle current from the spindle velocity controller, the cutting torque generated during milling can also be indirectly measured. Workpiece Chip Tool 3 Figure 1.2 Spindle disturbance torque response schematic However, to estimate the cutting torque from the spindle current, the dynamic distortion caused by the spindle mechanical assembly and motor dynamics must be compensated. This has been done in literature using a variety of techniques such as disturbance observers. However, these methods are often limited to research applications due to the modeling and control knowledge required for use and implementation. The knowledge requirements limit the application of these methods in industry and increases operating costs. To address this gap, this thesis presents a method of sensor-less cutting torque estimation requiring minimal expert input. The proposed method will provide a solution that is both low cost and easy to retrofit into existing machining centers. 4 The thesis presented is organized as follows. The literature related to the cutting force measurement and estimation techniques is reviewed in Chapter 2. A data-driven system identification method for spindles is applied in Chapter 3. Frequency Response Function (FRF) between the spindle’s torque and current command is measured using internal CNC commands, and the measured data is fitted to a state-space model using the data-driven method. The spindle’s model is used to build a Kalman Observer in Chapter 4 which is tuned automatically using machining torque simulations. The observer has been validated using an industrial CNC machine and tested for a variety of spindle speeds. The thesis is concluded with a summary of the research done and future directions in chapter 5. 5 Chapter 2: Literature Review 2.1 Overview The goal of this thesis is to create an automated and intelligent machining monitoring system that can predict cutting forces during machining operations. Many such systems exist in literature, with both direct and indirect approaches to measuring and monitoring the cutting forces. The direct approach relies on external measurement sensors installed in the cutting region. These methods provide accurate measurements of the cutting forces but require sensors and instrumentation that limit their applications to laboratory settings. On the other hand, the indirect methods employ controls and modeling techniques to estimate cutting forces using internal CNC measurements or sensors installed further from the cutting region. However, these methods are often limited as they require extensive control knowledge to apply and use. In this thesis, the cutting torque at the spindle will be predicted from servo control signals readily available within CNC systems. This chapter examines the methods in literature used to measure cutting torque and establish the study gap for an automated and sensor-less solution to cutting torque measurements. 2.2 Direct Measurement Approaches Many direct force measurement devices have been reported in the literature, using a variety of sensing techniques. Klocke et al. [2] and Totis et al. [3] used dynamometers mounted directly on the tool holder and table of the CNC respectively to measure and map cutting forces to the cutting tool path. These methods can provide accurate and reliable measurements of cutting forces but are very limited to laboratory environments. During cutting, the tool holder and table are subjected to flying chips, heat, coolant, and vibrations that are often hazardous to the sensors required by these 6 methods. Furthermore, the requirement of additional instrumentation limits the work volume of the CNC and adds additional flexibility to the structure of the machine tool. Some designs have attempted to address these issues by integrating force sensing elements in the tool holder. Inasaki et al. [4] developed a magnetostrictive torque sensor tool holder that can measure cutting torque during rotation. Mohring et al. [5] integrated piezoelectric sensors into the tool holder, and Haron et al. [6] created a tool holder with a strain-gage based rotating dynamometers. However, these methods still require instrumentation to measure and transmit the force data close to the cutting region, therefore decreasing the overall ruggedness of the design. 2.3 Indirect Measurement Approaches 2.3.1 Sensor-Based Indirect Measurement In this approach, the effects of cutting forces are measured indirectly and further away from the cutting region, such that the sensors and instrumentation are less exposed to the cutting environment and do not obstruct the machining work volume. Many different sensors placed on different locations of the machine tool spindle have been used in literature to indirectly measure force. In their literature review about the machine tool sensors, Tlusty and Andrews reported the use of strain gauges on the spindle bearings to measure the cutting forces indirectly by Promess [7]. Park and Altintas [8] integrated a ring-shaped force sensor between the spindle housing and the spindle flanges, and used piezoelectric force sensors to measure the forces experienced by the spindle. They applied a disturbance observer model to further increase the bandwidth of the measurement system, by accommodating for the effects of the spindle structure dynamics. Monno et al. [9] placed inductive displacement sensors to measure cutting forces and vibrations. Kalman filtering was used to fuse multiple sensors and estimate the cutting forces at the cutting tooltip. 7 However, these studies do not account for the speed [10] and temperature dependencies of the spindle structure dynamics. Kakino et al. [11] address this by using Eddy-current displacement sensors combined with thermocouples to accommodate the thermal displacement and change in the stiffness of the spindle. Similarly, Altintas et al. [12] measured the speed and load dependency of the spindle structural dynamics and used spindle mounted accelerometers to predict tool tip cutting forces and vibrations. These methods are capable of providing high-bandwidth force measurements of the cutting process, but are costly to implement and increases the maintenance cost of the spindle assembly. Furthermore, due to the indirect measurement nature, the selected sensors must be calibrated for each spindle/tool setup. These calibrations require modeling of the dynamic behavior from input force/torque to the measured values and also necessitate personnel knowledgeable in these topics. This further adds to the cost and complexity of the machine tool and is undesirable in industrial settings. 2.3.2 Sensor-less Indirect Measurement On the other hand, present-day CNC systems often allow access to internal signals in the controller to measure signals such as motor current commands [13]. Various studies in the literature have shown that the internal signals in a CNC machine can be used to estimate the cutting forces/torque. Altintas [14] demonstrated that by measuring the feed drive motor currents, the cutting forces can be predicted. This method required modeling of the feed-drive inertia and friction behavior and is limited by the bandwidth of the servo drive. Shinno et al. [15] designed a disturbance observer for a linear-driven aerostatic table system to estimate the cutting forces from the drive current. However, this was done via modeling only the rigid body dynamics of the system and cannot be applied in industrial machining centers as most are driven by ball-screw feed drives. Kakinuma 8 and Yamada [16] addressed this by using a Multi-Encoder-Based Disturbance Observer (MEDOB) to estimate cutting forces from ball-screw driven feed drives. A multi-encoder-based disturbance observer was used to model the feed-drive as a multi-inertia model and accommodate for the higher-order dynamics of the system. Similarly, Aslan and Altintas [17] used feed-drive currents to estimate cutting forces using a disturbance Kalman observer. The dynamics of the feed-drives were identified manually using hammer impact tests and the Kalman observer was used to compensate for the distortion it caused. Hence, they were able to increase the bandwidth of the force measurements past the bandwidth of the servo drive. To use the Kalman observer, the tuning parameters were manually adjusted such that the observer would provide good performance. Since these methods do not require an external sensory device, they can be directly retrofitted on existing machining centers. The main limitation of these methods is that they require accurate mathematical models of the force/torque to measurement dynamics. This can be first developed in laboratory settings using a variety of identification techniques, but given an industrial setting with many machining centers to identify it can become impractical. Furthermore, the disturbance observers used by these methods require manual tuning by someone knowledgeable in controls and dynamics, further increasing operating costs. There are some methods in literature that uses instead data-driven methods to estimate the cutting forces. Li et al. [18] used a neuro-fuzzy network to estimate the cutting force from feed-drive currents in a turning operation. Similarly, Kim et al [19] examined both disturbance observers and Artificial Neural Networks (ANN) models to estimate milling forces from feed-drive servo motor current. It was shown that the ANN system can provide higher accuracy force predictions when 9 compared against the disturbance observers. These data-driven methods allow for the estimation of cutting torque without tuning and has shown to be accurate. However, the main drawback of these systems is that they require large amounts of training data where the forces are known. This drastically increases the time required to implement these solutions and limits them to applications where the cutting force can be measured beforehand. 2.4 Discussion The methods used in literature for monitoring and predicting cutting forces/torque are summarized in the tables below: Table 1-1 Direct force/torque measuring methods Authors Method Benefits Drawbacks Klocke, Totis Directly mounted dynamometers Accurate, high bandwidth force measurements Delicated sensors, limited to lab applications Inasaki, Mohring, Haron Tool holder integrated sensors Sensors protected from cutting environment Bulky instrumentation near cutting area Table 1-2 Indirect sensor based force/torque measuring methods Authors Method Benefits Drawbacks Tlusty and Andrews Spindle mounted strain gages Sensors and instrumentation far away from cutting area Bandwidth limited by measurement dynamics Park and Altitnas, Monno, Kakino, Aslan and Altintas Spindle mounted sensors with compensated dynamics Further increases bandwidth of measurement solution Modelling of measurement dynamics needed, expensive to implement 10 Table 1-3 Indirect sensorless force/torque measuring methods Authors Method Benefits Drawbacks Altintas Estimation of force from feed drive current No sensors required Low bandwidth Shinmo Force from aerostatic feed drive current compensated using rigid body dynamics Higher bandwidth Limited to aerostatic tables, does not consider higher-order dynamics Kakinuma and Yamada, Alsan and Altintas Force from feed drive current using disturbance observers High bandwidth, considers higher-order dynamics of ball screw feed drives Modelling of feed drive dynamics needed, Observer tuning required Li, Kim Data-driven estimation of cutting force from feed drive current No manual modeling or tuning required Requires a large volume of training data with force measurements From this literature survey, there currently exist many solutions to the monitoring and measurement of cutting forces during machining. However, the direct measurement systems require expensive and unpractical instrumentation, whilst the indirect method requires modeling of the machine tool structure or the measurement system itself. Furthermore, indirect methods often employ methods that require tuning by experts in the control field. This increases the operating cost and complexity of the monitoring systems. Some data-driven methods attempt to address this issue, but still requires training data with pre-measured cutting forces. This makes these methods time consuming as the training data will have to be acquired for each machine. Therefore, it can be noted that there currently exists a gap in literature for an automated and industrial-friendly solution to machine force monitoring. In this thesis, this gap will be addressed 11 by introducing data-driven system identification techniques to the indirect force sensing techniques. to Furthermore, the tuning required by disturbance observer techniques will be automated by applying machining simulation estimate optimal observer parameters. 12 Chapter 3: Automated Spindle Drive System Identification This chapter presents the first step to automating the Kalman filter design, which is automating system identification. The model fitting concerns the identification of the spindle disturbance dynamics between spindle motor current and cutting torque. For that purpose, the spindle disturbance dynamics are first evaluated by measuring frequency response functions (FRFs) of the spindle velocity control loop. These measured FRFs are then combined to form the disturbance FRF from spindle current to disturbance torque. Using the combined FRF, the Eigen System Realization Algorithm (ERA) [20] is used to derive the state-space system model required to later predict the cutting torque. Improvements of the fitting ERA algorithm are proposed to minimize user input by supplementing an optimal threshold and ensuring system stability. Lastly, the fitting algorithm is validated in simulation with arbitrary generated systems to assess the quality of the derived model. 3.1 Spindle Disturbance Impulse Response Measurement During the cutting process, the spindle velocity control loop experiences perturbations in the form of disturbance torque due to a combination of friction and cutting loads: 𝜏𝑑 = 𝜏𝑑𝑐 + 𝜏𝑑𝑓 , (3.1) where 𝜏𝑑 is the disturbance torque, 𝜏𝑑𝑐 is the cutting torque and 𝜏𝑑𝑓 is the frictional torque. This is counteracted by the torque generated from the spindle current, which balances both the disturbance torque and the inertial loads of the spindle structure 𝐽𝑒𝑑Ω𝑎𝑐𝑡𝑑𝑡+ 𝐵𝑒Ω𝑎𝑐𝑡 = 𝐾𝑡𝐼𝑎𝑐𝑡 − 𝜏𝑑, (3.2) 13 where 𝐽𝑒 is the inertia of the spindle and Ω𝑎𝑐𝑡 is the actual spindle velocity. To predict the cutting torque from spindle current, the dynamic relationship between spindle current and disturbance torque must be first identified. After which the friction torque must be separated from the cutting torque. Since the spindle is often operated at a constant velocity, the resulting friction can be removed by simply measuring the mean spindle current when the tool is not cutting, and subtracting it from the cutting current. The torque to current dynamics is identified using the sine sweep technique. A sinusoidal signal with increasing frequency is given as input and the corresponding response is measured. The torque-current FRF is then calculated by taking the Fourier transform of both input and output signals and dividing the cross-power of the output and input by the auto-power of the input as given in Equation (3.4). 𝐺(𝑗𝜔) =𝐺𝑖𝑛𝑝𝑢𝑡(𝑗𝜔)𝐺𝑜𝑢𝑡𝑝𝑢𝑡∗ (𝑗𝜔)𝐺𝑖𝑛𝑝𝑢𝑡(𝑗𝜔)𝐺𝑖𝑛𝑝𝑢𝑡∗ (𝑗𝜔) (3.3) Where 𝐺𝑖𝑛𝑝𝑢𝑡(𝑗𝜔), 𝐺𝑜𝑢𝑡𝑝𝑢𝑡(𝑗𝜔) are the Fourier transform of the input and output signals respectively. Since it is difficult to directly perform a sweep of the torque-current dynamics in the frequency domain, the latter are reconstructed indirectly by combining frequency sweep data of various parts of the spindle velocity control loop. For a standard industrial CNC spindle, the velocity control loop block diagram can generally be simplified as shown in Figure 3.1. 14 Figure 3.1: Simplified CNC Spindle Velocity Control Loop Block Diagram Where Ω𝑛𝑜𝑚 is the commanded spindle velocity, which is fed to the velocity controller 𝐺Ω. The controller commands a nominal spindle current 𝐼𝑛𝑜𝑚, which is fulfilled by the current controller 𝐺𝐼. This produces the actual current 𝐼𝑎𝑐𝑡 in the motor, which generates a torque proportional to 𝐼𝑎𝑐𝑡 by torque constant 𝐾𝑡. For induction motors, the torque generated is usually also load and speed dependent. However, previous studies [21] have shown that this motor can be approximated as a DC motor for its nominal operating range. To construct the disturbance dynamics FRF, the following frequency sweeps must be identified: 1) Velocity open loop 𝐺𝑜 =𝛺𝑎𝑐𝑡𝛺𝑛𝑜𝑚= 𝐺𝛺(𝑗𝜔)𝐺𝑚(𝑗𝜔)𝐾𝑡𝐺𝑠𝑝(𝑗𝜔), (3.4) 2) Velocity control loop 𝐺𝑙𝑜𝑜𝑝 =𝛺𝑎𝑐𝑡𝐼𝑛𝑜𝑚= 𝐺𝑚(𝑗𝜔)𝐾𝑡𝐺𝑠𝑝(𝑗𝜔), (3.5) 3) Mechanical dynamics 15 𝐺𝑠 =𝛺𝑎𝑐𝑡𝐼𝑎𝑐𝑡𝑢𝑎𝑙= 𝐾𝑡𝐺𝑠𝑝(𝑗𝜔), (3.6) Given that the current control loop bandwidth is typically much higher than the mechanical dynamics, the motor voltage feedback 𝐾𝑏 is ignored for simplicity. Using measurable frequency responses 𝐺𝑜, 𝐺𝑙𝑜𝑜𝑝, and 𝐺𝑠, the dynamics of the velocity controller and spindle can be extracted as follows: • Velocity controller 𝐺Ω 𝐺Ω(𝑠) =𝐺𝑜(𝑠)𝐺𝑙𝑜𝑜𝑝(𝑠), (3.7) • Mechanical spindle dynamics 𝐺𝑠𝑝 𝐺𝑠𝑝(𝑠) =𝐺𝑠(𝑠)𝐾𝑡. (3.8) Having identified frequency responses 𝐺𝛺 and 𝐺𝑠𝑝, the torque-current dynamics is finally identified as follows: 𝐺𝑑𝜏(𝑠 = 𝑗𝜔) =𝐼𝑛𝑜𝑚(𝑗𝜔)𝜏𝑑(𝑗𝜔)=𝐺𝑠𝑝(𝑗𝜔)𝐺𝛺(𝑗𝜔)1 + 𝐺𝛺(𝑗𝜔)𝐺𝑚(𝑗𝜔)𝐾𝑡𝐺𝑠𝑝(𝑗𝜔)=𝐺𝑠𝑝(𝑗𝜔)𝐺𝛺(𝑗𝜔)1 + 𝐺𝑜(𝑗𝜔), (3.9) which is later used to develop a mathematical dynamics model. 3.2 Fitting to State-Space System Having obtained the frequency response of the spindle disturbance torque dynamics, the impulse response of the system is obtained via inverse Fourier Transform. Then the Eigensystem Realization Algorithm is applied to automatically fit the measurements into a discrete state-space system of the following form: 16 𝒙𝑘+1 = 𝑨𝒙𝑘 + 𝑩𝒖𝑘, 𝒚𝑘 = 𝑪𝒙𝑘 +𝑫𝒖𝑘 + 𝑣𝑘 . (3.10) Where 𝒖𝑘 represents the disturbance torque 𝜏𝑑 input to the system, and 𝒚𝑘 expresses the measured output of the spindle (Nominal spindle current 𝐼𝑛𝑜𝑚) which is corrupted by some white gaussian additive noise 𝑣. The dynamics of the system is described by the state matrix 𝑨 , the input matrix 𝑩 , the output matrix 𝑪, and the feedthrough matrix 𝑫 which are the unknowns to be identified. 𝒙𝑘 describes the state vector, which represents the process states of the system. Since there can be multiple realizations of a state-space system and that the system will be identified from data, the state vector has no direct physical meaning. Instead, it is more convenient to think of the system as an equivalent multiple mass spring damper system, and the states are the velocities and positions of such a system. Given Equation (3.10), the impulse input 𝒖𝑘 = {1 𝑘 = 00 𝑘 ≠ 0} , (3.11) would result in an output response of the following form 𝒚𝑘 = {𝑫 + 𝑣𝑘 𝑘 = 0𝑪𝑨𝑘−1𝑩+ 𝑣𝑘 𝑘 ≠ 0}. (3.12) By sampling time-shifted segments of the impulse response and populating the rows of a matrix with them, the Hankel matrix 𝐇 is formed as such: 𝐇 = [ 𝒚1 𝒚2 𝒚3 ⋯ 𝒚𝛽𝒚2 𝒚3 𝒚4 ⋯ 𝒚𝛽+1𝒚3 𝒚4 𝒚5 ⋯ 𝒚𝛽+2⋮ ⋮ ⋮ ⋱ ⋮𝒚𝛼 𝒚𝛼+1 𝒚𝛼+2 ⋯ 𝒚𝛼+𝛽] . (3.13) Substituting in Equation (3.12), the Hankel matrix 𝐇 can be shown to be the product of the observability 𝓞 and controllability matrix 𝓒 of the system, corrupted by some noise 𝑣, 17 𝐇 = [ 𝑪𝑩 + 𝑣1 𝑪𝑨𝑩 + 𝑣2 𝑪𝑨𝟐𝑩 + 𝑣3 ⋯ 𝑪𝑨𝛽−1𝑩+ 𝑣𝛽𝑪𝑨𝑩 + 𝑣2 𝑪𝑨𝟐𝑩 + 𝑣3 𝑪𝑨𝟑𝑩 + 𝑣4 ⋯ 𝑪𝑨𝛽𝑩 + 𝑣𝛽+1𝑪𝑨𝟐𝑩 + 𝑣3 𝑪𝑨𝟑𝑩 + 𝑣4 𝑪𝑨𝟒𝑩 + 𝑣5 ⋯ 𝑪𝑨𝛽+1𝑩+ 𝑣𝛽+2⋮ ⋮ ⋮ ⋱ ⋮𝑪𝑨𝛼−1𝑩+ 𝑣𝛼 𝑪𝑨𝛼𝑩 + 𝑣𝛼+1 𝑪𝑨𝛼+1𝑩+ 𝑣𝛼+2 ⋯ 𝑪𝑨𝛼+𝛽−2𝑩+ 𝑣𝛼+𝛽] = 𝓞𝓒 + 𝜎𝑽, (3.14) where 𝑽 is a matrix with independent and identically distributed zero-mean and unit variance random variables and 𝜎 is the variance of the measurement noise. To extract the noise-free minimal realization of the observability and controllability matrices, singular value decomposition (SVD) is applied. SVD decomposes the Hankel matrix into the following form: 𝐇 = 𝑼𝚺𝑽∗, (3.15) where 𝑼 and 𝑽 represent the eigen-time-delay coordinates of each pattern in the impulse response, and 𝚺 is the diagonal singular values matrix representing the energy of the corresponding pattern and notation * indicates matrix conjugate transpose. By applying a threshold to the energy values, the matrix is further broken down into truncated and non-truncated values as such: 𝑼𝚺𝑽∗ = [?̃? 𝑼𝒕] [?̃? 𝟎𝟎 𝚺𝒕] [?̃? 𝑽𝒕], (3.16) where ?̃?, ?̃? represent the dominant, high energy patterns to keep, and 𝑼𝒕, 𝑽𝒕 represent the patterns to be truncated. The asymptotic mean squared error (AMSE) optimal value of the threshold applied can be estimated as a function of the measurement noise [22]: 𝑟 =4√3√𝑛𝜎, (3.17) 18 where 𝑟 is the calculated optimal singular value hard threshold and 𝑛 represents the size of the Hankel matrix. Since in most cases the noise variance 𝜎 is unknown, Gavish and Donoho [23] showed that the ratio can be then estimated from the median of SVD energy values: 𝜏 =4√3𝜇𝛽Σ𝑚𝑒𝑑𝑖𝑎𝑛 (3.18) Where 𝜇𝛽 is a constant calculated by solving the following equation numerically; ∫√(4 − 𝑡)𝑡2𝜋𝑡𝑑𝑡 =12𝝁𝜷𝟎 (3.19) By calculating the optimal singular value hard threshold, the minimum realization of Observability ?̃? and Controllability ?̃? matrices are reconstructed using the truncated energy and eigen-values: ?̃? = ?̃?√?̃?, ?̃? = √?̃??̃?∗. (3.20) Now, system dynamics matrices 𝑨, 𝑩, 𝑪,𝑫 are then estimated from the observability and controllability matrices. First, to estimate the state dynamics matrix 𝑨, a second Hankel matrix is formed with the impulse response shifted one time-step forward. This Hankel matrix can be seen to be the product of the observability, state, and controllability matrices: 𝐇′ = [ 𝒚2 𝒚3 𝒚4 ⋯ 𝒚𝛽+1𝒚3 𝒚4 𝒚5 ⋯ 𝒚𝛽+2𝒚4 𝒚5 𝒚6 ⋯ 𝒚𝛽+3⋮ ⋮ ⋮ ⋱ ⋮𝒚𝛼+1 𝒚𝛼+2 𝒚𝛼+3 ⋯ 𝒚𝛼+𝛽+1] , 𝐇′ = [ 𝑪𝑨𝑩 𝑪𝑨𝟐𝑩 𝑪𝑨𝟑𝑩 ⋯ 𝑪𝑨𝛽𝑩𝑪𝑨𝟐𝑩 𝑪𝑨𝟑𝑩 𝑪𝑨𝟒𝑩 ⋯ 𝑪𝑨𝛽+1𝑩𝑪𝑨𝟑𝑩 𝑪𝑨𝟒𝑩 𝑪𝑨𝟓𝑩 ⋯ 𝑪𝑨𝛽+2𝑩⋮ ⋮ ⋮ ⋱ ⋮𝑪𝑨𝛼𝑩 𝑪𝑨𝛼+1𝑩 𝑪𝑨𝛼+2𝑩 ⋯ 𝑪𝑨𝛼+𝛽−1𝑩] = 𝓞𝐀𝓒. (3.21) 19 The state dynamics matrix 𝐀 can then be estimated by solving using the previously determined observability and controllability matrices. the estimate for the state matrix reads: ?̃? = ?̃?−𝟏𝐇′?̃?−𝟏 = ?̃?−12?̃?𝐇′?̃??̃?12. (3.22) The output matrix and the control matrix can be also extracted from the observability and controllability matrices as such: ?̃? = ?̃?12?̃?∗ [𝑰𝒑 𝟎𝟎 𝟎] , ?̃? = [𝑰𝒒 𝟎𝟎 𝟎] ?̃??̃?12. (3.23) Since the torque-current dynamics of the velocity control system in a CNC machine is always guaranteed to be stable, the identified state-space system must be stable as well. Therefore, to ensure the stability of the identified system, stable-unstable decomposition is applied. First, the Schur Decomposition [24] of the system is applied such that a unitary matrix 𝑻 is used to transform the state-space system as such; ?̃?𝒔𝒆𝒑 = 𝑻−𝟏?̃?𝑻 = [?̃?𝟏𝟏⏞𝑠𝑡𝑎𝑏𝑙𝑒?̃?𝟏𝟐𝟎 ?̃?𝟐𝟐⏟𝑢𝑛𝑠𝑡𝑎𝑏𝑙𝑒] ?̃?𝒔𝒆𝒑 = 𝑻−𝟏?̃? ?̃?𝑠𝑒𝑝 = 𝑻?̃? (3.24) Where ?̃?𝟏𝟏contains the stable eigenvalues of the system, ?̃?𝟐𝟐 contains the unstable eigenvalues, and ?̃?𝟏𝟐 is the coupling between the stable and unstable systems. To fully decouple the stable and unstable subsystems, a second transformation matrix is found by solving the following Lyapunov equation [25]: 20 ?̃?𝟏𝟏𝑺 − 𝑺?̃?𝟐𝟐 + ?̃?𝟏𝟐 = 𝟎 (3.25) The second transformation matrix is then constructed with the solution as follows: 𝑾 = [𝑰𝟏𝟏 𝑺𝟎 𝑰𝟐𝟐], (3.26) where 𝑰𝟏𝟏, 𝑰𝟐𝟐 are identity matrices with the size of ?̃?𝟏𝟏, ?̃?𝟐𝟐 respectively. By applying the transformation matrix, the decoupled system is found: ?̃?𝑑𝑒𝑐 = 𝑾−𝟏?̃?𝒔𝒆𝒑𝑾 =[ ?̃?𝒔𝒆𝒑𝟏𝟏⏞ 𝑠𝑡𝑎𝑏𝑙𝑒𝟎𝟎 ?̃?𝒔𝒆𝒑𝟐𝟐⏟ 𝑢𝑛𝑠𝑡𝑎𝑏𝑙𝑒] ?̃?𝑑𝑒𝑐 = 𝑾−𝟏?̃?𝒔𝒆𝒑 ?̃?𝑑𝑒𝑐 = 𝑾?̃?𝒔𝒆𝒑 (3.27) Then, the stable decoupled system is taken as the fitted system, and the unstable system is removed, thus guaranteeing the stability of the fitted system. The method discussed above provides a fully automated method of identifying the spindle torque-current dynamics requiring no expert input for use. This method is then tested in simulation in the next section to ensure its validity with a variety of system and noise conditions. 3.3 Validation The setup used for validation is a Quaser UX-600 machining center equipped with a Weiss 176039 – V7 spindle. It is controlled with a iTNC 530 Heidenhain CNC system. Using an external computer and the built-in diagnostic software TNCOpt, the required FRF measurements in Chapter 3.1 were obtained and the identified disturbance torque FRF is calculated using Equation (3.9) as 21 shown below in Figure 3.2. 22 Figure 3.2: Spindle disturbance torque to nominal current FRF at the first spindle speed range For this spindle motor, the induction motor consists of two operating ranges where the first one operates at 0-3300 RPM, and the second at 3300 RPM and above. Therefore, a second frequency sweep identification is done at 5000RPM, and the following FRF was measured. 23 Figure 3.3 Spindle disturbance torque to nominal current FRF at the second spindle speed range To validate the measured disturbance FRF, an impact test is conducted as shown in Figure 3.4. The spindle velocity is set to 0 rev/min, and an impact hammer is used to give a torque impulse to the spindle structure. The hammer force and the resulting spindle current are recorded, and the hammer force is converted to impact torque by multiplying the radius of the cutting tool hit. Then, the FRF is calculated using Equation (3.3). 24 Figure 3.4: Spindle impact test setup Using the obtained measurements, the tap test FRF and the sine swept FRF are compared against each other in Figure 3.5. Figure 3.5: Comparison of torque-current spindle FRF obtained from sine sweep and tap test 25 As shown. the tap test FRF and measured FRF matches closely with each other at lower frequencies. At the higher frequency region, there are larger discrepancies due to the bandwidth limitation of the hammer but the same trends are observed. Since the hammer cannot generate a perfect impulse impact on the tool, it could only excite the spindle structure up to 1 kHz. Thus the tap test FRF is not as reliable at the higher frequencies. On the other hand, the built-in sine – swept based FRF measurement function of CNC covers a wider frequency range. Since it is impossible to use an impact hammer on the spindle while it is in rotation, the second spindle speed FRF could not be tested. 3.3.1 Eigensystem Realization Algorithm Validation Simulation Setup To emulate the data collected in a real-world measurement, a typical stable system in state-space form is generated with an arbitrary number of states. This is done using the MATLAB function ‘rss’ [26] (Random State-Space). The number of states of the system is randomly selected between a range that is realistic for industrial spindles (between 10 – 30 states). This system is used to generate a noisy FRF measurement similar to what would be measured from an industrial CNC spindle. The steps to generating the noisy FRF is listed as follows: 1. A stable system with a specified number of states is randomly generated 2. The impulse response of the generated system is obtained using the MATLAB ‘impulse’ command 3. The impulse response is corrupted with additive white gaussian noise with a specified signal to noise ratio 4. The noise corrupted impulse response data is used in Equation (3.3) to calculate an FRF 26 Figure 3.6 summarizes of the process to generate the simulated noisy frequency response data. This FRF data is then fed to the system identification algorithm, which fits an estimated system to the data. The estimated system and the original system are then compared in terms of the frequency response function to confirm the validity of the algorithm discussed. The Normalized Root Mean Squared Error (NRMSE) of the frequency response is used to evaluate the fit and is defined as: 𝑓𝑖𝑡 = 100 (1 −‖𝐺𝑎𝑐𝑡𝑢𝑎𝑙(𝑗𝜔) − 𝐺𝑓𝑖𝑡(𝑗𝜔)‖‖𝐺𝑎𝑐𝑡𝑢𝑎𝑙(𝑗𝜔) − 𝑚𝑒𝑎𝑛(𝐺𝑎𝑐𝑡𝑢𝑎𝑙(𝑗𝜔)‖), where 𝐺𝑎𝑐𝑡𝑢𝑎𝑙(𝑗𝜔) represents the frequency response of the simulated system, and 𝐺𝑓𝑖𝑡(𝑗𝜔) represents the frequency response of the fitted system. Figure 3.6: Steps for generating noisy simulated FRFs 27 3.3.2 ERA Validation Simulation Results For a randomly generated 10-state system, the simulated noisy FRF is shown in Figure 3.7. Figure 3.7: Simulated Noisy FRF of 𝝉𝒅(𝒋𝝎)𝑰𝒏𝒐𝒎(𝒋𝝎) for 10 State System with signal to noise ratio SNR = 10 dB The simulated response was generated with a signal to noise ratio (SNR) of 10 dB. The fit achieved with the proposed algorithm in Chapter 3.3.2 is shown in Figure 3.8. The estimated system achieved an NRMSE fit of 95.3% when compared to the original system, where the main deviation from the real system occurs mostly in the high-frequency region. This is caused by the lower signal to noise ratio in the high-frequency regions of the FRF. 28 Figure 3.8: Fit results of 𝝉𝒅(𝒋𝝎)𝑰𝒏𝒐𝒎(𝒋𝝎) for 10 State System with SNR = 10 dB with Normalized Root Mean Square NRMSE = 95.3 % Varying the order of the system, it can be seen that the algorithm achieves a good estimation of the system regardless of system size. Where fitting the frequency response of a system with 20 states and 30 states resulted in a fit of 93.1% and 94.6%, respectively. Figure 3.9 Fit results of 𝝉𝒅(𝒋𝝎)𝑰𝒏𝒐𝒎(𝒋𝝎) for 20 State System with SNR = 10 dB with NMRSE = 93.1 % 29 Figure 3.10 Fit results of 𝝉𝒅(𝒋𝝎)𝑰𝒏𝒐𝒎(𝒋𝝎) for 30 State System with SNR = 10 dB with NRMSE = 94.6 % Furthermore, the SNR of the generated noisy frequency response was varied to ensure the algorithm handles a variety of measurement conditions. Figure 3.11 and Figure 3.12 show the simulated frequency response with 0 dB SNR and the resulting fit, respectively. 30 Figure 3.11 Simulated Noisy FRF of 𝝉𝒅(𝒋𝝎)𝑰𝒏𝒐𝒎(𝒋𝝎) for 10 State System with SNR = 0 dB Figure 3.12 Fit results of 𝝉𝒅(𝒋𝝎)𝑰𝒏𝒐𝒎(𝒋𝝎) for 10 State System with SNR = 0 dB with NRMSE = 95.8% 31 Lastly, the algorithm was compared against MATLAB’s ‘modalfit’ function which estimates the natural frequencies of a system from measured frequency-response functions. Figure 3.13 Comparison between MATLAB ‘modalfit’ and ERA As shown in Figure 3.13, in addition to requiring less manual tuning, the ERA method produces a closer fit than the MATLAB ‘modalfit’ [27] function. Furthermore, the MATLAB algorithm is much more sensitive to higher noise and results in higher levels of overfitting when compared to the method proposed. 3.3.3 Spindle FRF Fit Results Finally, the proposed algorithm was applied to the disturbance FRF measured on the Quaser UX600 machine tool center for validation purposes. A state-space model with a model order of 28 and 22 were fitted to the combined spindle current-torque frequency response of the first and 32 second speed range respectively. An NRMSE fit of 83.2% was achieved for the first speed range, and 80.5% was achieved for the second speed range as shown in Figure 3.14 and Figure 3.15 respectively. The normalized root-mean-square error was calculated using a manually noise reduced FRF to reduce over-inflating the error values due to the large noise content. The fitted state-space model parameters can be found in Appendix A.1. Figure 3.14 Fit result from measured first speed range torque-current spindle FRF with NRMSE = 83.2 % 33 Figure 3.15: Fit result from measured second-speed range torque-current spindle FRF with NRMSE = 80.5 % As shown, the magnitude fit result shows a very good fit at the lower frequency region of the measured FRF and can be seen as a good approximation. At the higher frequency region, some overfitting can be observed. This is due to the extremely low signal to noise ratio levels in this region. As well, the fitted system phase is compared against the measured system phase for the first and second speed range of the spindle shown in Figure 3.16 and Figure 3.17 respectively. 34 Figure 3.16 ERA Fit phase comparison for first spindle speed range Figure 3.17 ERA Fit phase comparison for second spindle speed range A similar result as the magnitude comparison can be seen where the phase matches closely at the lower frequencies, and at higher frequencies the high amount of noise in the measured FRF causes some overfitting in the identified system. The comparison for the real and imaginary parts 35 of the fitted and measured FRFs are shown in Appendix A.2. Overall, it can be seen that the identification method achieves a reasonably good fit for an automated system identification solution. 3.4 Discussion The algorithm proposed in this chapter demonstrates a method for identifying the spindle disturbance torque to current transfer function automatically with no expert input. The algorithm requires no tuning and achieves high accuracy fits for a variety of systems and measurement conditions. The steps to the identification method discussed are listed as follows: 1. Frequency sweeps are run using Heidenhain diagnostic software TNCOpt for the required measurements discussed in Chapter 3.1 2. Spindle torque-current frequency response data is obtained by combining measured data using Equation (3.9) 3. State-space model representation of the spindle torque-current dynamics are fitted using the Eigensystem realization algorithm 4. System order is determined automatically using Equation (3.17) 5. System stability is ensured using Schur decomposition In these steps, the only actions requiring a manual operator is the first step where the frequency sweeps for the spindle are conducted. When compared against pre-existing solutions in MATLAB, the proposed algorithm achieves both higher accuracy fits and requires less manual tuning of the fitting algorithm. Fitting is however sensitive to the signal-to-noise ratio which lowers its accuracy in the high noise FRF measurements. The method discussed provides a much more industrial friendly approach to identifying any spindle dynamics quickly and easily. 36 Chapter 4: Automated Kalman Filter Design In this chapter, the previously identified dynamics model is used to create a Kalman observer to estimate the cutting torque from spindle nominal current measurements. First, the model is augmented such that the unknown input torque is modeled as a state affected by process noise. Then a Kalman observer is designed using the expanded dynamics model to estimate the cutting torque. To tune the Kalman observer for optimal performance, the measurement and process noise covariance values are needed. The measurement noise covariance is obtained from the measurement resolution, whilst the process noise covariance is estimated via simulation. Using the obtained covariance values, the Kalman observer is tuned and applied to an industrial CNC machine for validation. 4.1 State Space Transformation From Chapter 3, the state-space model identified from the current-torque spindle dynamics are given as; 𝒙𝑘+1 = 𝑨𝒙𝑘 + 𝑩𝑢𝑘 , 𝑦𝑘 = 𝐼𝑛𝑜𝑚 = 𝑪𝒙𝑘 + 𝑫𝒖𝑘 + 𝑣𝑘 𝑢𝑘 = 𝜏𝑑(𝑘) (4.1) where 𝒙𝒌 is the discrete-time state vector of the spindle drive disturbance dynamics, 𝑦𝑘 is the measured nominal current, 𝑣𝑘 is the measurement noise and 𝑢𝑘 represents the unknown input cutting torque. The state-space model is represented in the following schematics in Figure 4.1. 37 Figure 4.1 State space system schematic To estimate the input using an observer, the state-space model is transformed via a linear transformation matrix such that the input torque 𝑢𝑘is modeled as another state as given; 𝒙𝒆𝒙𝒑𝒌 = {[𝒙𝒌]𝑢𝑘} (4.2) Subsequently, the state-space model becomes; 𝒙𝒆𝒙𝒑𝒌 = [𝑨 𝑩0 1]⏟ 𝑨𝒆𝒙𝒑{[𝒙𝒌]𝑢𝑘} + [𝟎1]𝑤𝑘 𝑦𝑘 = 𝐼𝑛𝑜𝑚 = [𝑪 0]⏟ 𝑪𝒆𝒙𝒑{[𝒙𝒌]𝑢𝑘} + 𝑣𝑘 (4.3) 1+ + = +1 38 In this model, the unknown input is modeled as a state that is changed by some unknown process noise 𝑤𝑘. The new model schematic can be seen in Figure 4.2. Figure 4.2 Expanded state-space model schematic In the expanded state-space model, the unknown input torque is modeled as a state of the system which is changed by some process noise 𝑤𝑘: 𝑢𝑘+1 = 𝑢𝑘 + 𝑤𝑘 (4.4) This model is slightly modified from the model presented by Altintas and Aslan [28] in that the process noise is modeled as the change to the unknown input as opposed to the input itself. For cutting processes, the cutting torque are typically harmonic signals. Taking the time derivative of it removes the DC signal component and shifts the mean towards zero. This better satisfies the assumption that the process noise has zero mean for designing a Kalman Observer. 1+ + + + + 39 4.2 Kalman Observer Design Using the model obtained in Section 4.1, the Kalman observer is then designed to estimate the cutting torque from the measured spindle current. This is done by estimating the expanded state vector (𝒙𝒆𝒙𝒑𝒌), which contains the unknown cutting torque. The Kalman observer is designed to minimize the error between the estimated state vector (?̂?𝒆𝒙𝒑𝒌) and the actual state vector [29]; 𝒆 = 𝒙𝒆𝒙𝒑𝒌 − ?̂?𝒆𝒙𝒑𝒌 (4.5) Where 𝒆 is the state estimation error vector. Since the expanded state vector includes the unknown cutting torque as shown in Equation (4.3), this Kalman observer will also find the minimum error estimate of the cutting torque. It is formulated in two steps where the first step finds the a priori estimate of the state vector and the state estimation error covariance; ?̂?𝒆𝒙𝒑𝒌+𝟏′ = 𝑨𝒆𝒙𝒑?̂?𝒆𝒙𝒑𝒌 𝑷𝒌+𝟏′ = 𝑨𝒆𝒙𝒑𝑷𝒌𝑨𝒆𝒙𝒑 + 𝑸 (4.6) where ?̂?𝒆𝒙𝒑𝒌+𝟏′ is the a priori estimate of the state vector at time 𝑘 + 1 and 𝑸 is the covariance of the process noise 𝑤𝑘. This step uses the identified system dynamics found in Chapter 3: to predict what the system state vector will be in the next time step. As well, it estimates the covariance of the estimation error in the next time step, which signifies the variability of the estimation error. The covariance for the estimation error is calculated as: 𝑷 = 𝑐𝑜𝑣(𝒆) =∑(𝒆𝒊 − ?̅?)2𝑁 − 1𝑁𝑖=1 (4.7) Having obtained an future estimate using the system dynamics, the a priori estimate is then corrected using measurement data via recursive least-squares: 40 ?̂?𝒆𝒙𝒑𝒌+𝟏 = 𝒙𝒆𝒙𝒑𝒌+𝟏′ + 𝑷𝒌+𝟏′ 𝑪𝒆𝒙𝒑𝑇 𝑹−𝟏 (𝑦𝑘 − 𝑪𝒆𝒙𝒑?̂?𝒆𝒙𝒑𝒌+𝟏′ ) 𝑷𝒌+𝟏 = 𝑷𝒌+𝟏′ − 𝑷𝒌+𝟏′ 𝑪𝒆𝒙𝒑𝑇 (𝑪𝒆𝒙𝒑𝑷𝒌+𝟏′ 𝑪𝒆𝒙𝒑𝑇 + 𝑹)−𝟏𝑪𝒆𝒙𝒑𝑷𝒌+𝟏′ (4.8) where 𝑹 is the covariance of the measurement noise 𝑣𝑘, and ?̂?𝒆𝒙𝒑𝒌+𝟏 and 𝑷𝒌+𝟏 are the posterior estimates of the state matrix and error covariance matrix respectively. The measurement noise covariance, like the estimation error covariance, describes the variability of the measurement noise. By calculating Equation (4.6) and (4.8) recursively using the measured current (𝑦𝑘), the state vector can be estimated. From the state vector, the input torque can be obtained by selecting the last state of the vector as follows; ?̂?𝑑𝑘 = [00⋮1] ?̂?𝒆𝒙𝒑𝒌 (4.9) Typically, the error covariance matrix converges to a constant quickly [29] Therefore, for ease of implementation, the error covariance matrix can be assumed to be time-invariant and its steady-state value can be calculated by solving the Discrete Algebraic Riccati Equation as such; 𝑨𝒆𝒙𝒑𝑷𝒔𝑨𝒆𝒙𝒑𝑻 − 𝑷𝒔 +𝑸 − 𝑨𝒆𝒙𝒑𝑷𝒔𝑪𝒆𝒙𝒑𝑻 (𝑹 + 𝑪𝒆𝒙𝒑𝑷𝒔𝑪𝒆𝒙𝒑𝑻 )−𝟏𝑪𝒆𝒙𝒑𝑷𝑨𝒆𝒙𝒑𝑻 = 𝟎 (4.10) This yields the steady-state value of the error covariance, and the Kalman Observer can be defined by the steady-state Kalman Observer gain: 𝑲𝒔𝒔 = 𝑷𝒔𝑪𝒆𝒙𝒑𝑻 𝑹−𝟏 (4.11) Then, by combining Equation (4.6), (4.8), and (4.11) the steady-state Kalman Observer simplifies to the following; ?̂?𝒆𝒙𝒑𝒌+𝟏 = 𝑨𝒆𝒙𝒑?̂?𝒆𝒙𝒑𝒌 +𝑲𝒔𝒔(𝑦𝑘 − 𝑪𝒆𝒙𝒑𝑨𝒆𝒙𝒑?̂?𝒆𝒙𝒑𝒌) (4.12) The schematic for the steady-state Kalman Observer can be shown below: 41 Figure 4.3 Kalman Observer Schematic Using Equation (4.9) and (4.12), the cutting torque can subsequently be estimated from the nominal current. 4.3 Kalman Observer Automatic Tuning To achieve optimal performance of the Kalman observer, the values of the measurement noise covariance 𝑹 and the process noise covariance 𝑸 must be determined. These covariances specify the relative relationship between the measurement noise variability and the variability of the unknown input torque. It is important to estimate these values correctly, as it determines the balance the trade off between compensating the spindle torque-current dynamics distortion versus noise rejection. 1+ + = + 1 + +1 + 42 4.3.1 Measurement Noise Covariance Estimation Firstly, the measurement noise covariance can be assessed by determining the additive noise present in the nominal current measurements. In this case, since the measured values from the CNC are the nominal commanded current values, the measurement noise present is due to the measurement resolution. The probability distribution of this noise can be modeled as a uniform distribution, as shown in Figure 4.4. The variance for this distribution can be calculated as; 𝜎2 =112(𝑏 − 𝑎)2 (4.13) where 𝑏 and 𝑎 are the upper and lower limits of the measurement resolution, respectively. By measuring/obtaining the measurement resolution of the spindle current commands, the measurement noise covariance can subsequently be determined using Equation (4.13). Figure 4.4 Uniform Distribution 43 4.3.2 Process Noise Covariance Estimation via Milling Torque Simulation On the other hand, the process noise specifies the rate of change in the input torque as shown in Equation (4.9). That is to say, the process noise 𝑤𝑘 is the incremental change in input torque 𝑢𝑘 between each time interval; 𝑤𝑘 = 𝑢𝑘+1 − 𝑢𝑘 (4.14) Therefore, by estimating cutting torque using milling torque simulations, the characteristics of the process noise 𝑤𝑘 can also be estimated. The steps to simulating milling torque are outlined below. The cutting torque from metal cutting results from the force needed to shear the workpiece material using the cutter. The torque that arises can be modeled as follows [30]: 𝑇𝑐 =𝐷2(𝐾𝑡𝑐𝑎ℎ + 𝐾𝑡𝑒𝑎) (4.15) Where 𝐷 is the cutting tool diameter, 𝑎 is the depth of cut, ℎ is the thickness of the chip cut and 𝐾𝑡𝑐 and 𝐾𝑡𝑒 are the cutting force and edge force coefficients respectively. 𝐾𝑡𝑐 represents the cutting forces resulting from material shearing, and 𝐾𝑡𝑒 represents cutting forces from friction at the edge of the cutting tool. These values are determined experimentally via cutting force measurements, and the depth of cut depends on the operator-determined cutting conditions. However, the chip cut thickness (ℎ) is a time-varying parameter that depends on the instantaneous position and angle of the endmill. 44 Figure 4.5 Endmill Cutting Chip Thickness Diagram As shown in Figure 4.5, the thickness of material cut by the endmill depends on the instantaneous engagement angle of the cutter 𝜙𝑐. For a rotating endmill, the chip thickness value varies periodically and can be approximated as: ℎ(𝜙𝑐) = 𝑐 sin(𝜙𝑐) (4.16) where 𝑐 is the feed-rate of the tool into the cutting material, measured in distance per tooth rotation. Moreover, an endmill typically has multiple cutting surfaces arranged in a helical pattern as shown in Figure 4.6 45 Figure 4.6 Endmill Helix Angle Diagram For any moment in time, the engagement angle 𝜙𝑐 is then also dependent on the cutting tool helix angle and varies as a function of length along the tool. The engagement angle change due to helix angle along the length of the tool can be found as; 𝜙𝑐𝛽 =2𝑧 tan(𝛽)𝐷 (4.17) where 𝑧 is the axial cutting depth, and 𝛽 is the helix angle of the cutting tool. Lastly, since there could be multiple endmill flutes engaged with the workpiece at any given time, the engagement angle for each tooth needs to be calculated. The pitch angle between each endmill tooth is found as a function of the number of flutes 𝑁 as follows: 𝜙𝑝 =2𝜋𝑁 (4.18) Combining Equations (4.16), (4.17), and (4.18), the engagement angle can be found for any given point and rotation angle of the cutting tool. To simulate the cutting torque, the rotation angle is increased at the desired spindle speed for small incremental time intervals. At each rotation angle, the cutter engagement angle is calculated for each point on the tool, and the cutting torque is 46 calculated accordingly. The cutting torque is then integrated along the length of the tool for each cutting flute. The steps of the simulation can be summarized as follows: 1. Find the tool rotation angle, and split the cutting tool into small segments along the length of the tool; 2. Calculate the engagement angle for each segment of the tool using Equation (4.17); 3. Calculate chip thickness and cutting torque using Equation (4.16) and (4.15) respectively for each tool segment; 4. Sum the cutting torque along the length of the tool; 5. Repeat for each cutting blade in the endmill; 6. Increment rotation angle and repeat. 7. Using this procedure, the cutting torque can be simulated provided that the cutting conditions such as depth of cut are known. Using the simulated cutting torque, the process noise is estimated by sampling the simulated cutting torque at the nominal current sampling frequency and applying Equation (4.14). From the simulated cutting torque values and Equation (4.3), the process noise covariance can be found as follows; 𝑄 =∑(𝑤𝑖 − ?̅?)2𝑁 − 1𝑁𝑖=1 (4.19) where 𝑤𝑖 is the estimated process noise, and ?̅? is the average estimated process noise value. Now that the process noise covariance 𝑄 and the measurement noise covariance 𝑅 are found, the 47 Discrete Algebraic Riccati Equation in Equation (4.12) can be solved and the steady-state Kalman Observer can be obtained. 4.4 Results The Kalman observer designed was validated using a Quaser UX 600 CNC Machine controlled via Heidenhain iTNC 530. Nominal spindle current values were measured using Heidenhain TNCScope at a 10 kHz sampling frequency. A Kistler 9123C1011 rotary dynamometer and tool holder was used to measure the cutting torque at 25600 kHz to compare and validate the cutting torque estimated from spindle current. To ensure that the measurement dynamics of the Kistler rotary dynamometer does not distort the measured torque, an impact test identification was carried out by applying an instrumented hammer impact to blade of the cutting tool held by the rotary dynamometer. The hammer impact force is then converted to an impact torque using the diameter of the cutting tool. The dynamometer measured torque and the hammer impacted torque were both measured and the following measurement dynamics FRF were found as shown in Figure 4.7. 48 Figure 4.7 Rotary dynamometer measurement dynamics FRF From the FRF measured, the bandwidth measured using the hammer impact test for the rotary dynamometer is seen to be around 800 Hz. As well, some resonance peaks can be seen at around 1100 Hz, however these peaks are most likely due to noise in the measurement since the impact hammer can only excite frequencies up to around 800 Hz as shown in Figure 4.8. 49 Figure 4.8 Rotary dynamometer FRF measurement input spectrum The experiments were carried out using a 12mm 2-Fluted high-speed steel tool, which was used to cut aluminum at a 2 mm depth of cut with a 0.2 mm/tooth feed-rate. Multiple slotting (full radial engagement) cuts were done using spindle speeds of 3000, 5000, 6000, and 7000 RPM. For these spindle speeds, the tooth passing frequencies are 100, 166, 200 and 233 Hz respectively, which are well within the measurement bandwidth of the rotary dynamometer. For the Quaser UX 600 spindle controller, the nominal current measurement resolution was 0.01A. Using Equation (4.13), the measurement covariance was calculated to be 0.0029. The cutting torque was simulated using the algorithm discussed in Chapter 3.3.1, using cutting coefficients obtained from the MAL in-house software CutPro. The cutting coefficients used for 𝐾𝑡𝑐 was 797.6 𝑁/𝑚𝑚2 and 25.7 𝑁/𝑚𝑚 for 𝐾𝑡𝑒. The simulated cutting torque for a 3000 RPM cut is found using the steps in Chapter 4.3.2 is shown below in Figure 4.9. 50 Figure 4.9 Simulated Cutting Torque Sampling the simulated cutting torque at the 10 kHz, the time differential of the cutting torque is found as shown in Figure 4.10. 51 Figure 4.10 Simulated Cutting Torque Change per Sampling Interval By taking the variance of the simulated torque differential, the process noise covariance is estimated to be 0.0035. Having determined the measurement and process noise covariance values, the Kalman Observer is designed by augmenting the model obtained in Chapter 3 using Equation (4.3). The steady-state Kalman Observer gain is found as the following by solving Equation (4.10) and (4.11) using the found covariance values: 𝑲𝒔𝒔 = [−0.86,0.44,0.21,0.16, −0.15,0.051,−0.054,0.067,0.028,0.036,0.01, −0.011, −0.013, −0.015,0.015,0.011,0.017,0.012,0.0036,−0.014,0.026,0.039,0.038,−0.01,0.0005,0.051] Using the steady-state Kalman Gain found above, the Kalman Observer is applied using Equation (4.12) to predict cutting torque from spindle current. The cutting torque estimated from nominal 52 spindle current at 3000 RPM is compared against the dynamometer reference torque as shown in Figure 4.11. 53 Figure 4.11 Kalman Observer Estimated Torque (3000 RPM spindle speed) 54 To evaluate the effectiveness of the estimated data, the Normalized Root Mean Squared Error is calculated using the range of the dynamometer measured torque as the normal reference: 𝑵𝑹𝑴𝑺𝑬 =√∑(𝜏𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑(𝑖) − 𝜏𝑑𝑦𝑛𝑜(𝑖))2𝑛𝑛𝑖=1max(𝜏𝑑𝑦𝑛𝑜) − min(𝜏𝑑𝑦𝑛𝑜)∗ 100% (4.20) Where 𝜏𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑 is the Kalman Observer estimated torque, and 𝜏𝑑𝑦𝑛𝑜 is the reference torque measured using the rotary dynamometer. Using Equation (4.20), the NRMSE was found to be 10.4%. Further experiments were conducted at the same cutting conditions as previously stated at 5000,6000, and 7000 RPM spindle speeds. For the cutting force measurements at other spindle speeds, the measurement noise covariance values stay constant since the measurement resolution stays constant. On the other hand, the process noise covariance changes based on the cutting conditions. For 5000,6000, and 7000 RPM spindle speeds, the cutting torque is simulated for each case using the methods discussed in Chapter 4.3.2. Using the simulated cutting torque, the process noise covariance, and subsequent Kalman Observer gain values are as found for each case as follows: 55 Table 1-1 Kalman Observer Parameters Process Noise Covariance Kalman Observer Gain 0.01 [1.17,0.21,0.66,0.3,0.31,0.028,0.006,−0.029,0.043, −0.019,0.026,−0.017 0.026 − ,0.0023, −0.0126,0, −0.016, −0.042,0.074,−0.12,−0.04,0.009, −0.024,−0.027,0.0002,1.61] 0.0142 [1.28,0.19,0.74,0.35,0.38,0.034,0.006,−0.036,0.051,−0.022,0.031,−0.02, 0.03,−0.0032,−0.015,0.0016,−0.019,−0.049,0.089,−0.14, −0.048,0.01, −0.029,−0.0319,0.0002,1.89] 0.0192 [1.39,0.18,0.81,0.39,0.45,0.04,0.006,−0.042,0.059,−0.026,0.036,−0.023, 0.035,−0.0042,−0.017,0.0036,−0.022,−0.056,0.013,−0.16,−0.056,0.013, −0.034,−0.037,0.0001,2.17 Using the identified Kalman Observers, the cutting torque for the experiment done at spindle speeds of 5000, 6000, and 7000 RPM are evaluated using Equation (4.12) and compared as shown in Figure 4.12, Figure 4.13, and Figure 4.14, respectively. 56 Figure 4.12 Kalman Observer Estimated Torque (5000 RPM spindle speed) 57 Figure 4.13 Kalman Observer Estimated Torque (6000 RPM spindle speed) 58 Figure 4.14 Kalman Observer Estimated Torque (7000 RPM spindle speed) 59 The errors for 5000, 6000 and 7000 RPM spindle speeds are 11.7%, 14.5% and 18.2% respectively. As can be seen, the designed Kalman observer performs well for a variety of input frequencies. The predicted torque matches relatively closely with the reference torque measured from the rotary dynamometer. At higher spindle speeds, the errors in the predicted torque are increased due to the effects of other unmodelled factors such as spindle runout become more apparent. The spindle runout causes uneven loading on the spindle bearings, which results in an additional harmonic friction force. As the spindle speed increases, the effects of spindle runout increase due to the higher centripetal forces whilst the cutting torque remains relatively the same. This results in distortion of the measured disturbance torque and increases the error observed. As well, occasional DC shifts in the predicted torque are observed likely due to temperature or power fluctuations in the spindle motor. This effect is reduced when the non-cutting spindle current mean is subtracted from the cutting spindle current values. Since the work presented is an automated procedure based off the work done by Altintas and Aslan [28], the results obtained are then also compared against observer estimated torque using the procedure outlined in their work. Since the spindle assembly used in Altintas and Aslan’s work are the same the spindle used in the previously presented results, the model for the spindle motor were copied over directly. Following the procedure outlined in their work, the following results were obtained for the same data measured at 3000, 5000, 6000 and 7000 RPM shown in Figure 4.15, Figure 4.16, Figure 4.17, and Figure 4.18 respectively. 60 Figure 4.15 Manual Kalman Observer Estimated Torque (3000 RPM Spindle Speed) Figure 4.16 Manual Kalman Observer Estimated Torque (5000 RPM Spindle Speed) 61 Figure 4.17 Manual Kalman Observer Estimated Torque (6000 RPM Spindle Speed) Figure 4.18Manual Kalman Observer Estimated Torque (7000 RPM Spindle Speed) 62 Where the errors measured between the dynamometer reference the estimated torque are 11.3%, 12.8%, 14.4% and 16.8% for the estimated torque at 3000, 5000, 6000 and 7000 RPM, respectively. As shown, the estimation error between the automated Kalman Observer and the manual Kalman Observer are similar at lower spindle frequencies. However, at higher spindle frequencies the manual Kalman Observer outperforms the automatically tuned observer. This is likely due to the high amounts of noise in the higher frequency regions of the measured FRF degrading the system fit accuracy at these higher frequencies. Whereas the manual observer is modelled based on physics of the spindle assembly, and therefore minimizing the potential for fitting error. At the same time, the modelling of the spindle assembly based on physics principles is a time consuming and complicated task and would therefore be difficult to apply to fleets of industrial CNC machine tools. From these comparisons, it is shown that the automatically designed Kalman Observer presented provides a streamlined observer design procedure at the expense of high frequency estimation accuracy. 4.5 Discussion This chapter presented a method of designing and implementing a Kalman Observer automatically requiring minimal user input of information. The only information required is that which are readily available on industrial CNC shop floors. A brief schematic of the Kalman Observer is shown in Figure 4.11. 63 Figure 4.11 Kalman Observer Schematic The Kalman Observer is tuned automatically via pre-simulation of the cutting process using known cutting parameters. Since cutting parameters are easily obtainable information in industrial CNC shops, the algorithm presented can be shown as an industrial friendly application for predicting the cutting torque. Furthermore, the sensorless approach allows for easy retrofitting of this monitoring solution onto pre-existing machining setups. Finally, the method presented was validated on an industrial CNC and shown that it could be practically applied to predict cutting torques during the machining process. Using the predicted torque obtained, the work presented in this thesis can then be applied to machine tool monitoring applications. 64 Chapter 5: Conclusion This thesis presented an automated system of predicting cutting torques from the CNC spindle commanded current. Two major steps are presented in this thesis: identifying and modeling the spindle dynamics, as well as automatically designing a Kalman observer. The torque to current spindle dynamics is identified using frequency sweeps of the velocity control loop measured via built-in function within CNC systems. The dynamics are indirectly identified by combining multiple frequency sweeps of different points in the velocity control loop. The identified dynamics are then fitted to a state-space model using the Eigensystem realization algorithm. This algorithm is automatically tuned using an asymptotic mean squared error (ASME) optimal thresholding technique to reduce overfitting, while at the same time also eliminating the need to manually adjust the fitting algorithm. Furthermore, Schur decomposition is used to remove unstable eigenvalues in the fitted model, since the spindle dynamics can always be assumed to be stable. The proposed algorithm was validated first in simulation, by fitting the simulated frequency response of random dynamics models and comparing the fitted model against the original simulated model. The algorithm showed above 90% normalized root mean squared error fit in a variety of system orders and for a variety of noise conditions. Finally, the proposed algorithm is applied to the spindle dynamics measured from a UX600 machining center and achieved a fit of above 80% Normalized Root Mean Square Error (NRMSE) for both spindle speed configurations. The identified state-space model is then applied to design a Kalman observer, such that the torque-current dynamics can be compensated. Since the goal of the observer is to estimate the unknown input, the state-space model is transformed such that the input is a state in the state vector. The input is subsequently modeled as a state that is changed by a random process noise. This 65 augmented state-space model is then used to design a steady-state Kalman observer by solving the discrete algebraic Riccati equation. To tune the observer, the covariance of the process noise was calculated by pre-simulating the cutting torque using known cutting conditions. This method is then validated on the UX 600 machining center and shown to estimate the cutting torque from spindle current accurately on a variety of cutting speeds. From the validation test estimation errors, it was shown that method presented can be reasonably applied for use in machine tool monitoring during high force and low spindle speed roughing operations. During these operations, the method presented could be applied to monitor and adaptively limit the cutting torque experienced by the cutting tool and spindle, thus minimizing the chances for machine tool spindle or cutting tool damage. Furthermore, if failures does occur, the resulting changes in cutting torque could be detected using the same algorithm, so as to warn the operator of the failure and prevent further damage. The observer presented combined with the potential monitoring tasks provide a framework for developing more intelligent machine tools. 5.1.1 Future Works This thesis serves as a basis for providing an automated and sensor-less cutting force monitoring solution. However, there are still some limitations and drawbacks in the system that could be addressed by future research: • The identified disturbance dynamics relies on the fact that the torsional stiffness of the cutting tool and spindle. If there are torsional dynamics below the frequencies that can be measured using CNC frequency sweeps (< 2500Hz), then these torsional dynamics should be identified and compensated. 66 • Spindle runout causes varying loads in the spindle ball bearings, and as a result, it creates additional disturbance torque at the spindle rotation speed harmonics. This effect decreases the accuracy of the disturbance torque estimations, and further modeling and identification can be applied to mitigate its effects. • The induction motor of the spindle is modeled using the spindle current frequency sweep and assumes a constant torque constant. This is shown to be a reasonable assumption at nominal spindle operating ranges, however, for high cutting load or high speed cutting the motor exhibits more non-linear speed and load dependent effects. Further modeling and application of non-linear observers could improve the accuracy of the system for high cutting load or high-speed cutting conditions. • The frequency response of the disturbance transfer function is measured during non-cutting and is assumed to be constant during cutting. However, temperature and load dependent behavior in the motor may cause the disturbance dynamics to change. Future research could be done to apply joint Kalman Observer techniques to simultaneously predict cutting torque and update model dynamics. 67 References [1] R. Teti, K. Jemielniak, G. O’Donnell, and D. Dornfeld, “Advanced monitoring of machining operations,” CIRP Ann. - Manuf. Technol., vol. 59, no. 2, pp. 717–739, 2010, doi: 10.1016/j.cirp.2010.05.010. [2] F. Klocke, S. Kratz, and D. Veselovac, “Position-oriented process monitoring in freeform milling,” CIRP J. Manuf. Sci. Technol., vol. 1, no. 2, pp. 103–107, 2008, doi: https://doi.org/10.1016/j.cirpj.2008.09.003. [3] G. Totis, O. Adams, M. Sortino, D. Veselovac, and F. Klocke, “Development of an innovative plate dynamometer for advanced milling and drilling applications,” Measurement, vol. 49, pp. 164–181, 2014, doi: https://doi.org/10.1016/j.measurement.2013.11.049. [4] H. Ohzeki, A. Mashine, H. Aoyama, and I. Inasaki, “Development of a Magnetostrictive Torque Sensor for Milling Process Monitoring,” J. Manuf. Sci. Eng., vol. 121, no. 4, pp. 615–622, Nov. 1999, doi: 10.1115/1.2833078. [5] H. C. Möhring, K. M. Litwinski, and O. Gümmer, “Process monitoring with sensory machine tool components,” CIRP Ann. - Manuf. Technol., vol. 59, no. 1, pp. 383–386, Jan. 2010, doi: 10.1016/j.cirp.2010.03.087. [6] M. Rizal, J. A. Ghani, M. Z. Nuawi, and C. H. C. Haron, “Development and testing of an integrated rotating dynamometer on tool holder for milling process,” Mech. Syst. Signal Process., vol. 52–53, no. 1, pp. 559–576, Feb. 2015, doi: 10.1016/j.ymssp.2014.07.017. [7] J. Tlusty and G. C. Andrews, “A Critical Review of Sensors for Unmanned Machining,” CIRP Ann. - Manuf. Technol., vol. 32, no. 2, pp. 563–572, Jan. 1983, doi: 10.1016/S0007-8506(07)60184-X. 68 [8] Y. Altintas and S. S. Park, “Dynamic compensation of spindle-integrated force sensors,” CIRP Ann. - Manuf. Technol., vol. 53, no. 1, pp. 305–308, Jan. 2004, doi: 10.1016/S0007-8506(07)60703-3. [9] P. Albertelli, M. Goletti, M. Torta, M. Salehi, and M. Monno, “Model-based broadband estimation of cutting forces and tool vibration in milling through in-process indirect multiple-sensors measurements,” Int. J. Adv. Manuf. Technol., vol. 82, no. 5–8, pp. 779–796, Feb. 2016, doi: 10.1007/s00170-015-7402-x. [10] H. Cao, B. Li, and Z. He, “Chatter stability of milling with speed-varying dynamics of spindles,” Int. J. Mach. Tools Manuf., vol. 52, no. 1, pp. 50–58, Jan. 2012, doi: 10.1016/j.ijmachtools.2011.09.004. [11] A. A. D. SARHAN, A. MATSUBARA, M. SUGIHARA, H. SARAIE, S. IBARAKI, and Y. KAKINO, “Monitoring Method of Cutting Force by Using Additional Spindle Sensors,” JSME Int. J. Ser. C, vol. 49, no. 2, pp. 307–315, Dec. 2006, doi: 10.1299/jsmec.49.307. [12] M. Postel, D. Aslan, K. Wegener, and Y. Altintas, “Monitoring of vibrations and cutting forces with spindle mounted vibration sensors,” CIRP Ann., vol. 68, no. 1, pp. 413–416, Jan. 2019, doi: 10.1016/j.cirp.2019.03.019. [13] J. F. G. Oliveira, F. Ferraz Júnior, R. T. Coelho, and E. J. Silva, “Architecture for machining process and production monitoring based in open computer numerical control,” Proc. Inst. Mech. Eng. Part B J. Eng. Manuf., vol. 222, no. 12, pp. 1605–1612, Dec. 2008, doi: 10.1243/09544054JEM1156. [14] Y. Altintas, “Prediction of cutting forces and tool breakage in milling from feed drive current measurements,” J. Manuf. Sci. Eng. Trans. ASME, vol. 114, no. 4, pp. 386–392, 69 Nov. 1992, doi: 10.1115/1.2900688. [15] H. Shinno, H. Hashizume, and H. Yoshioka, “Sensor-less monitoring of cutting force during ultraprecision machining,” CIRP Ann. - Manuf. Technol., vol. 52, no. 1, pp. 303–306, Jan. 2003, doi: 10.1016/S0007-8506(07)60589-7. [16] Y. Yamada and Y. Kakinuma, “Sensorless cutting force estimation for full-closed controlled ball-screw-driven stage,” Int. J. Adv. Manuf. Technol., vol. 87, no. 9–12, pp. 3337–3348, Dec. 2016, doi: 10.1007/s00170-016-8710-5. [17] D. Aslan and Y. Altintas, “Prediction of Cutting Forces in Five-Axis Milling Using Feed Drive Current Measurements,” IEEE/ASME Trans. Mechatronics, vol. 23, no. 2, pp. 833–844, Apr. 2018, doi: 10.1109/TMECH.2018.2804859. [18] X. Li, A. Djordjevich, and P. K. Venuvinod, “Current-sensor-based feed cutting force intelligent estimation and tool wear condition monitoring,” IEEE Trans. Ind. Electron., vol. 47, no. 3, pp. 697–702, 2000, doi: 10.1109/41.847910. [19] T. Y. Kim, J. Woo, D. Shin, and J. Kim, “Indirect cutting force measurement in multi-axis simultaneous NC milling processes,” Int. J. Mach. Tools Manuf., vol. 39, no. 11, pp. 1717–1731, Nov. 1999, doi: 10.1016/S0890-6955(99)00027-9. [20] J. N. Juang and R. S. Pappa, “An eigensystem realization algorithm for modal parameter identification and model reduction,” J. Guid. Control. Dyn., vol. 8, no. 5, pp. 620–627, May 1985, doi: 10.2514/3.20031. [21] D. Aslan and Y. Altintas, “On-line chatter detection in milling using drive motor current commands extracted from CNC,” Int. J. Mach. Tools Manuf., vol. 132, no. April, pp. 64–80, 2018, doi: 10.1016/j.ijmachtools.2018.04.007. [22] M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is 4/√3,” 70 IEEE Trans. Inf. Theory, vol. 60, no. 8, pp. 5040–5053, 2014, doi: 10.1109/TIT.2014.2323359. [23] M. Gavish and D. L. Donoho, “for Singular Values is 4 / 3,” pp. 1–14. [24] “Schur decomposition - MATLAB schur.” https://www.mathworks.com/help/matlab/ref/schur.html (accessed Oct. 26, 2020). [25] S. K. Nagar and S. K. Singh, “An algorithmic approach for system decomposition and balanced realized model reduction,” J. Franklin Inst., vol. 341, no. 7, pp. 615–630, Nov. 2004, doi: 10.1016/j.jfranklin.2004.07.005. [26] “Generate random continuous test model - MATLAB rss.” https://www.mathworks.com/help/control/ref/rss.html (accessed Oct. 26, 2020). [27] “Modal parameters from frequency-response functions - MATLAB modalfit.” https://www.mathworks.com/help/signal/ref/modalfit.html?s_tid=srchtitle (accessed Oct. 26, 2020). [28] D. Aslan and Y. Altintas, “On-line chatter detection in milling using drive motor current commands extracted from CNC,” Int. J. Mach. Tools Manuf., vol. 132, no. February, pp. 64–80, 2018, doi: 10.1016/j.ijmachtools.2018.04.007. [29] R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman Filtering with Matlab Exercises. Wiley, 2012. [30] Y. Altintas, Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design, 2nd ed. USA: Cambridge University Press, 2012. 71 Appendix A.1 Identified State-Space Matrix for Spindle Disturbance Dynamics Parameters72 Table A.1-1 First Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 1-14 0.973 -0.009 -0.007 0.022 0.003 -0.001 0.011 -0.009 -0.007 0.007 -0.002 -0.002 0.000 -0.004 0.000 0.808 -0.434 -0.058 -0.209 -0.078 0.012 0.041 0.028 -0.004 0.077 -0.043 0.014 0.017 0.000 0.160 0.808 -0.078 -0.297 -0.112 0.036 0.047 0.043 0.001 0.114 -0.063 0.020 0.022 0.000 0.000 0.000 0.588 0.605 -0.126 0.061 0.064 0.101 -0.009 0.148 -0.071 0.024 0.037 0.000 0.000 0.000 -0.765 0.588 -0.050 0.011 0.014 0.001 0.007 0.047 -0.030 0.009 0.006 0.000 0.000 0.000 0.000 0.000 0.550 -0.771 -0.003 0.077 0.012 0.080 -0.032 0.009 0.015 0.000 0.000 0.000 0.000 0.000 0.796 0.550 0.046 -0.007 -0.016 0.027 -0.021 0.008 0.011 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.487 0.822 0.013 -0.076 0.029 -0.010 -0.023 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.845 0.487 -0.018 -0.016 0.001 -0.001 -0.003 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.013 -0.950 -0.046 0.014 0.027 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.994 -0.013 -0.002 0.001 0.003 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.485 -0.875 -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.862 -0.485 0.013 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.544 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.833 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 73 Table A.1-2 First Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 14-28 -0.004 -0.011 0.013 0.005 -0.012 -0.005 0.017 0.001 0.002 0.007 0.003 0.010 0.006 0.004 -0.068 0.030 0.077 -0.029 -0.044 0.023 -0.013 0.060 -0.085 -0.006 0.002 -0.070 0.068 0.042 -0.102 0.040 0.119 -0.038 -0.069 0.028 -0.010 0.082 -0.120 -0.002 0.003 -0.093 0.100 0.060 -0.108 0.077 0.116 -0.049 -0.062 0.036 -0.036 0.085 -0.129 -0.006 -0.017 -0.123 0.104 0.056 -0.052 0.008 0.064 -0.016 -0.038 0.012 0.004 0.041 -0.059 0.002 0.010 -0.038 0.049 0.032 -0.056 0.035 0.068 -0.011 -0.039 0.005 -0.002 0.025 -0.050 0.018 -0.015 -0.041 0.050 0.020 -0.026 0.014 0.025 -0.022 -0.014 0.018 -0.015 0.037 -0.044 -0.021 0.007 -0.039 0.029 0.022 0.038 -0.052 -0.037 0.021 0.019 -0.014 0.026 -0.028 0.044 0.004 0.025 0.057 -0.037 -0.014 0.008 -0.011 -0.012 -0.007 0.005 0.011 -0.007 0.017 -0.011 -0.017 0.009 -0.011 0.003 0.012 -0.072 0.055 0.080 -0.029 -0.043 0.016 -0.019 0.047 -0.077 0.000 -0.014 -0.072 0.065 0.029 0.005 0.006 -0.009 -0.005 0.005 0.006 -0.010 0.006 -0.003 -0.010 -0.005 -0.011 0.000 0.003 0.009 -0.004 -0.010 0.002 0.005 -0.002 0.001 -0.006 0.009 -0.001 0.001 0.007 -0.007 -0.005 -0.039 0.022 0.043 -0.017 -0.024 0.010 -0.009 0.029 -0.044 -0.003 -0.002 -0.038 0.036 0.017 0.800 0.040 0.074 -0.026 -0.040 0.014 -0.012 0.044 -0.070 0.000 -0.005 -0.060 0.058 0.028 -0.544 -0.003 -0.027 0.005 0.016 -0.003 -0.003 -0.014 0.022 -0.004 -0.003 0.013 -0.020 -0.011 0.000 -0.765 -0.566 -0.033 -0.049 0.017 -0.019 0.053 -0.087 -0.001 -0.010 -0.077 0.072 0.032 0.000 0.613 -0.765 -0.009 -0.034 0.004 0.010 0.027 -0.043 0.008 0.007 -0.023 0.039 0.023 0.000 0.000 0.000 -0.829 -0.527 -0.013 0.010 -0.038 0.058 0.003 0.002 0.048 -0.047 -0.024 0.000 0.000 0.000 0.545 -0.829 -0.001 -0.007 -0.015 0.025 -0.007 -0.001 0.013 -0.025 -0.013 0.000 0.000 0.000 0.000 0.000 -0.889 0.437 0.022 -0.030 -0.009 -0.005 -0.034 0.020 0.012 0.000 0.000 0.000 0.000 0.000 -0.448 -0.889 -0.003 0.005 0.005 0.006 0.013 -0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.893 0.320 -0.009 -0.002 -0.092 0.087 0.050 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.395 -0.893 0.002 0.003 0.063 -0.060 -0.033 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.971 -0.216 0.012 -0.001 -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.215 -0.971 0.009 -0.005 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.950 -0.135 -0.042 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.052 -0.950 0.052 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.984 74 Table A.1-3 Second Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 1-11 0.967 -0.011 0.007 0.002 0.029 0.005 0.000 0.003 -0.006 0.001 -0.002 0.000 0.810 -0.183 -0.278 -0.084 0.006 -0.082 -0.020 -0.082 0.073 0.012 0.000 0.440 0.810 0.196 0.070 -0.001 0.058 0.016 0.053 -0.049 -0.009 0.000 0.000 0.000 0.605 -0.762 0.005 -0.036 -0.008 -0.037 0.031 0.004 0.000 0.000 0.000 0.612 0.605 0.000 -0.091 -0.026 -0.080 0.079 0.016 0.000 0.000 0.000 0.000 0.000 0.007 0.984 -0.011 -0.034 0.033 0.008 0.000 0.000 0.000 0.000 0.000 -0.999 0.007 -0.001 -0.004 0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.488 -0.856 -0.033 -0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.872 -0.488 -0.011 -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.524 0.849 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.832 -0.524 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 75 Table A-1-4 Second Spindle Speed Disturbance Dynamics State-Space A Matrix Columns 11-22 0.005 -0.011 0.028 0.000 -0.001 0.002 0.000 -0.001 0.012 -0.005 -0.002 -0.027 -0.054 0.076 0.024 0.107 -0.081 -0.225 0.015 -0.030 -0.043 -0.034 0.022 0.030 -0.032 -0.014 -0.072 0.060 0.153 -0.008 0.029 0.028 0.024 -0.009 -0.026 0.042 0.006 0.044 -0.037 -0.092 0.001 -0.006 -0.020 -0.018 -0.039 -0.043 0.040 0.030 0.116 -0.089 -0.248 0.022 -0.056 -0.041 -0.033 -0.015 -0.019 0.022 0.015 0.044 -0.030 -0.098 0.009 -0.020 -0.016 -0.009 0.001 -0.003 0.008 -0.001 0.001 -0.002 -0.002 -0.004 0.004 -0.001 -0.001 0.015 0.024 -0.033 -0.023 -0.049 0.025 0.099 -0.017 0.017 0.017 0.008 0.002 0.009 -0.019 -0.009 -0.016 0.003 0.029 -0.007 0.003 0.005 0.000 -0.004 -0.008 0.011 0.003 0.014 -0.010 -0.032 0.004 -0.004 -0.006 -0.004 -0.018 -0.024 0.026 0.030 0.047 -0.017 -0.108 0.029 -0.025 -0.015 -0.001 -0.791 -0.592 -0.027 -0.026 -0.040 0.012 0.089 -0.022 0.020 0.013 -0.001 0.602 -0.791 -0.029 -0.008 -0.015 0.002 0.030 -0.009 -0.002 0.006 0.001 0.000 0.000 -0.845 0.489 0.019 -0.012 -0.039 0.000 0.018 -0.011 -0.011 0.000 0.000 -0.460 -0.845 0.055 0.008 -0.135 0.056 -0.025 -0.016 0.016 0.000 0.000 0.000 0.000 -0.901 -0.371 -0.144 0.047 -0.033 -0.018 0.005 0.000 0.000 0.000 0.000 0.336 -0.901 0.126 0.004 0.021 0.025 0.027 0.000 0.000 0.000 0.000 0.000 0.000 -0.761 -0.143 0.126 0.079 0.008 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.925 -0.296 -0.015 0.020 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.235 -0.925 0.027 -0.013 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.988 -0.068 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.074 -0.988 76 First Spindle Speed Disturbance Dynamics State-Space B Matrix: [0.169,−0.040,−0.080, −0.103,−0.018, 0.005,−0.050, 0.028,−0.050,−0.028, −0.024 0.003 − 0.020,−0.023, 0.005,−0.025, 0.002, 0.016,−0.008,−0.048, 0.017,−0.073, 0.046, 0.017, 0.000, 0.080 − 0.061,−0.050]𝑇 First Spindle Speed Disturbance Dynamics State-Space C Matrix: [0.194, 0.150, 0.090, −0.024, 0.083, 0.030,−0.001, −0.017, 0.019, 0.001,−0.014, 0.016, 0.001, 0.014, 0.025, 0.041, −0.062,−0.018, 0.042,−0.012,−0.030,−0.031, 0.033, 0.008,−0.020, 0.033,−0.050,−0.053] Second Spindle Speed Disturbance Dynamics State-Space B Matrix: [ −0.186, 0.087,−0.047, 0.007, 0.135, 0.029,−0.009,−0.021, 0.000, 0.011, 0.018,−0.007, 0.007,−0.008 − 0.027, 0.022,−0.073,−0.121, 0.002,−0.050,−0.028,−0.028]𝑇 Second Spindle Speed Disturbance Dynamics State-Space C Matrix: [−0.200,−0.097, 0.156, −0.102, 0.019, 0.002,−0.020,−0.005,−0.019, 0.012,−0.001 0.009 − 0.027, 0.075,−0.021, 0.021,−0.045,−0.059,−0.012, 0.004,−0.020,−0.022] 77 A.2 Fitted and Measured System Real and Imaginary Comparison Figure A0.1 Real and Imaginary value comparision between fitted and measured FRF of First Spindle Speed Range 78 Figure A0.2 Real and Imaginary value comparision between fitted and measured FRF of Second Spindle Speed Range
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Automated design and implementation of Kalman observer...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Automated design and implementation of Kalman observer for spindle torque estimation in CNC machining Lu, Zhao Wei 2020
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Automated design and implementation of Kalman observer for spindle torque estimation in CNC machining |
Creator |
Lu, Zhao Wei |
Publisher | University of British Columbia |
Date Issued | 2020 |
Description | Milling is a subtractive manufacturing technique, where the material is continuously removed from a workpiece until the desired part shape is obtained. Heavily used in the automotive and aerospace industries, Computer Numerical Control (CNC) milling machines occupy a large chunk in the manufacturing process. The current research focuses towards machining and machine tool monitoring systems that are more self-sufficient and self-adjusting to adapt the processes to machine tools. Cutting torque delivered by the machine tool spindle is one of the key sensory signals for machining process monitoring. This thesis presents a method that automatically reconstructs cutting torque from motor current commands generated by the servo controller of the machine tool. To estimate the cutting torque from commanded spindle current, the dynamics between torque to current relationship must be modeled and compensated. The thesis first presents an automated identification of spindle dynamics using data-driven system identification methods. The frequency response function (FRF) of the spindle dynamics is measured manually using CNC internal diagnostic tools. The identified FRF is then automatically converted to a state-space model using the Eigensystem Realization Algorithm (ERA). To reduce overfitting, an optimal threshold is applied to the ERA method to limit the identified system order. And to ensure the stability of the identified system, unstable eigenvalues of the system are removed using Schur decomposition. The identified system is then augmented such that the unknown torque input is modeled as a state changed by a random process noise. This augmented system is used to create a Kalman Observer, which compensates the spindle dynamics and estimates the torque from the spindle nominal current. The Kalman Observer is tuned automatically by estimating the noise covariance values using machining simulations. The method was eventually validated on a Quaser UX 600 industrial CNC system. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2020-12-08 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NoDerivatives 4.0 International |
DOI | 10.14288/1.0395180 |
URI | http://hdl.handle.net/2429/76698 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2021-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2021_may_lu_zhaowei.pdf [ 2.6MB ]
- Metadata
- JSON: 24-1.0395180.json
- JSON-LD: 24-1.0395180-ld.json
- RDF/XML (Pretty): 24-1.0395180-rdf.xml
- RDF/JSON: 24-1.0395180-rdf.json
- Turtle: 24-1.0395180-turtle.txt
- N-Triples: 24-1.0395180-rdf-ntriples.txt
- Original Record: 24-1.0395180-source.json
- Full Text
- 24-1.0395180-fulltext.txt
- Citation
- 24-1.0395180.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0395180/manifest