Upgraded Kalman Filtering of Cutting Forces in Milling

1. Introduction

• sensitivity issues such as those affecting dynamometers based on strain gauges [5], capacitance sensors [6], optical sensors [7,8,9] and other sensors requiring large structural deformations at the inspected point in order to achieve a satisfactory signal to noise ratio; this is typically obtained by weakening the machining system, thus reducing the resonance frequency;
• load effects such as those introduced by rotating dynamometers [10,11] that may drastically reduce the first resonance frequency of the machine tool spindle (below 500 Hz) with respect to the unaltered configuration where no sensor (thus neither additional modal mass nor additional flexibility) is integrated into the machining system; and
• dynamic inertial disturbances deriving from the undesired but unavoidable oscillations of the modal mass placed in front of the piezoelectric load cells in case of plateform dynamometers; by so doing, even low-frequency vibration modes of the machine tool table where the dynamometer is clamped may cause undesired force fluctuations that cannot be avoided, as explained in [12].

2. Experimental Investigation on Dynamometer Dynamics

3. The Upgraded Augmented Kalman Filter (UAKF)

• Experimental modal analysis by means of pulse testing was carried out. The input hammer signal and $Q^{\prime }$ outputs were obtained from each measurement.
• $Q^{\prime }\times P$ nonparametric, empirical TFRFs $W_{qp}\left(\mathrm{j}{\omega }_k\right)$, coherence functions ${\gamma }_{qp}^2\left(\mathrm{j}{\omega }_k\right)$ ($k=1,2,\dots ,N/2$), and sampled Empirical Impulse Response Estimates (EIREs) $h_{qp}\left(t_n\right)$ with $n=1,\dots ,N$ were extracted from the experimental modal analysis data.
• All the EIREs were assembled into the matrix

  $$\mathbf{H}:=\left[\begin{array}{ccc}{\mathbf{h}}_{11}& \cdots & {\mathbf{h}}_{1P}\\ ⋮& \ddots & ⋮\\ {\mathbf{h}}_{Q^{\prime }1}& \cdots & {\mathbf{h}}_{Q^{\prime }P}\end{array}\right]$$

  (9)

  that has $(Q^{\prime }\times N)$ rows and P columns, as each ${\mathbf{h}}_{qp}$ is a $N\times 1$ column vector. The Singular Value Decomposition was applied to $\mathbf{H}$ in order to extract the most significant $P^{\prime }=6\ll P$ principal components carrying more than 99% of the information contained into $\mathbf{H}$. This allowed a considerable matrix compression that boosted the subsequent calculation.
• The discrete Green (Hankel) matrix ${\mathbf{H}}_{\mathrm{d}}$ and the velocity related retarded Green matrix ${\mathbf{H}}_{\mathrm{d}}^{\prime }$ were assembled from the compressed EIREs matrix according to the work in [29] based on the impulse dynamic subspace (IDS) approach. The Singular Value Decomposition was applied this time to ${\mathbf{H}}_{\mathrm{d}}$ and ${\mathbf{H}}_{\mathrm{d}}^{\prime }$. By a proper multiplication of the obtained principal components the natural pulsations ${\omega }_{\mathrm{n},m}$ and the damping ratios ${\xi }_m$ of the most important $m=1,\dots ,M$ vibration modes were obtained. In the current case, it was necessary to include $M\approx 50÷100$ vibration modes in order to achieve a satisfactory mathematical interpolation (though probably some modes were mostly mathematical artifacts rather than physical modes). The natural frequencies and damping factors were used to construct the state space matrix $\mathbf{A}$ as follows,

  $$\mathbf{A}=\left[\begin{array}{cccccc}& \mathbf{O}& & & \mathbf{I}& \\ -{\omega }_{\mathrm{n},1}^2& & & -2{\xi }_1{\omega }_{\mathrm{n},1}& & \\ & \ddots & & & \ddots & \\ & & -{\omega }_{\mathrm{n},M}^2& & & -2{\xi }_M{\omega }_{\mathrm{n},M}\end{array}\right]$$

  (10)

  where $\mathbf{O}$ and $\mathbf{I}$ are null and identity matrices of size $M\times M$.
• Each $W_{qp}\left(\mathrm{j}{\omega }_k\right)$ was then expressed as the sum

  $$W_{qp}\left(\mathrm{j}{\omega }_k\right)=\frac{F_{\mathrm{dyn},q}\left(\mathrm{j}{\omega }_k\right)}{F_p\left(\mathrm{j}{\omega }_k\right)}=\sum _{m=1}^M\frac{G_{m,qp}+G_{m,qp}^{\prime }\mathrm{j}{\omega }_k}{{\left(\frac{\mathrm{j}{\omega }_k}{{\omega }_{\mathrm{n},m}}\right)}^2+2{\xi }_m\frac{\mathrm{j}{\omega }_k}{{\omega }_{\mathrm{n},m}}+1},$$

  (11)

  where the static gains $G_{m,qp}$ were associated to elastic force contributions due to device deformation, whereas $G_{m,qp}^{\prime }$ were associated to viscous damping forces that are proportional to the modal velocities. These coefficients were determined by stepwise linear regression and weighted linear regression procedures that were carried out in the frequency domain. As a result, a total of $\left(4MP{Q}^{\prime }\right)$ modal parameters have been determined so far. This would have implied a very large state space vector of size $\left(4MP\right)\times 1$, that would have caused numerical issues if not properly handled. A considerable dimension reduction was then achieved by noting that the P column vectors composing the matrices $G_{m,qp}$ and $G_{m,qp}^{\prime }$ (for any fixed m) had to be all proportional to the same eigenmode shape inspected at the $Q^{\prime }$ available outputs. The constants of proportionality were expected to be in general different from each other since they corresponded to the participation factors associated to each of the P inputs. However, the two matrices $G_{m,qp}$ and $G_{m,qp}^{\prime }$ had to be rank-one matrices. In other words, their first principal component obtained from SVD should have provided a good approximation of their structure, as follows,

  $${\mathbf{G}}_m:=\left[G_{m,qp}\right]=\underset{\mathrm{from}\mathrm{SVD}}{\underbrace{{\mathbf{U}\mathbf{\Sigma }\mathsf{\Theta }}^{\mathrm{H}}}}\approx \underset{\mathrm{into}\mathrm{matrix}\mathbf{C}}{\underbrace{\left[\begin{array}{c}{u}_{m,1}^{\left(1\right)}\\ ⋮\\ u_{m,Q^{\prime }}^{\left(1\right)}\end{array}\right]{\sigma }_m^{\left(1\right)}}}\underset{\mathrm{into}\mathrm{matrix}\mathbf{B}}{\underbrace{\left[\begin{array}{ccc}{\vartheta }_{m,1}^{\left(1\right)}& \cdots & {\vartheta }_{m,P}^{\left(1\right)}\end{array}\right]}}$$

  (12)

  and similarly for $G_{m,qp}^{\prime }$

  $${\mathbf{G}}_m^{\prime }:=\left[G_{m,qp}^{\prime }\right]=\underset{\mathrm{into}\mathrm{matrix}\mathbf{B}}{\underbrace{{{\mathbf{V}\mathbf{S}\mathbf{\Lambda }}^{\mathrm{H}}}_{\mathrm{from}\mathrm{SVD}}\approx \underset{\mathrm{into}\mathrm{matrix}\mathbf{C}}{\underbrace{{\mathbf{v}}_m^{\left(1\right)}{s}_m^{\left(1\right)}}}\underbrace{{\lambda }_m^{\left(1\right)}}}}$$

  (13)
  Thus, only the modes satisfying this approximation were kept in the final model.
• Then, the matrix $\mathbf{B}$ could be assembled as follows,

  $$\mathbf{B}=\left[\begin{array}{ccc}& \mathbf{O}& \\ {\omega }_{n,1}^2{\vartheta }_{1,1}^{\left(1\right)}& \cdots & {\omega }_{n,1}^2{\vartheta }_{1,P}^{\left(1\right)}\\ ⋮& & ⋮\\ {\omega }_{n,M}^2{\vartheta }_{M,1}^{\left(1\right)}& \cdots & {\omega }_{n,M}^2{\vartheta }_{M,P}^{\left(1\right)}\\ & \mathbf{O}& \\ {\omega }_{n,1}^2{\lambda }_{1,1}^{\left(1\right)}& \cdots & {\omega }_{n,1}^2{\lambda }_{1,P}^{\left(1\right)}\\ ⋮& & ⋮\\ {\omega }_{n,M}^2{\lambda }_{M,1}^{\left(1\right)}& \cdots & {\omega }_{n,M}^2{\lambda }_{M,P}^{\left(1\right)}\end{array}\right]$$

  (14)

  whereas matrix $\mathbf{C}$ was written as follows,

  $$\mathbf{C}=\left[\begin{array}{cccccccc}{u}_{1,1}^{\left(1\right)}{\sigma }_1^{\left(1\right)}& \cdots & u_{M,1}^{\left(1\right)}{\sigma }_M^{\left(1\right)}& & & v_{1,1}^{\left(1\right)}{s}_1^{\left(1\right)}& & v_{M,1}^{\left(1\right)}{s}_M^{\left(1\right)}\\ ⋮& & ⋮& \mathbf{O}& \mathbf{O}& & & \\ u_{1,Q^{\prime }}^{\left(1\right)}{\sigma }_1^{\left(1\right)}& \cdots & u_{M,Q^{\prime }}^{\left(1\right)}{\sigma }_M^{\left(1\right)}& & & v_{1,Q^{\prime }}^{\left(1\right)}{s}_1^{\left(1\right)}& & v_{M,Q^{\prime }}^{\left(1\right)}{s}_M^{\left(1\right)}\end{array}\right]$$

  (15)
  Matrix $\mathbf{B}$ was of size $4M\times P$, while matrix $\mathbf{C}$ was of size $Q^{\prime }\times 4M$. Accordingly, the global state space matrix $\mathbf{A}$ had to be doubled, i.e.,

  $$\mathbf{A}:=\left[\begin{array}{cc}\mathbf{A}& \mathbf{O}\\ \mathbf{O}& \mathbf{A}\end{array}\right]$$

  (16)
  Thus, the associated state space vector was of size $4M\times 1$, which was still large but less problematic. In detail, the final state space vector size was in the range $(200÷400)\times 1$ instead of $(3200÷6400)\times 1$.
• In order to improve the global model adequacy, the coefficients of the $\mathbf{B}$ matrix were recalculated by means of another (weighted) linear regression procedure that was carried out in the frequency domain, by exploiting the relation

  $$\underset{\begin{array}{l}\mathrm{known}\mathrm{term}\\ \mathrm{from}\mathrm{experiment}\end{array}}{\underbrace{\mathbf{W}\left(\mathrm{j}{\omega }_k\right)}}\cong \underset{\mathrm{coefficient}\mathrm{matrix}}{\underbrace{\mathbf{C}{\left(\mathrm{j}{\omega }_k\mathbf{I}-\mathbf{A}\right)}^{-1}}}\underset{\mathrm{unknown}}{\underbrace{\mathbf{B}}}.$$

  (17)
• The obtained dynamic system may be used for accurately estimating the measured outputs from known inputs. Nevertheless, here we were interested in the inverse problem: estimating all the inputs from the outputs. This is in general not possible when two or more distinct inputs generate very similar outputs. In order to extract the maximum amount of information regarding the input structure, the Singular Value Decomposition was applied on matrix $\mathbf{B}$, thus yielding

  $$\mathbf{B}=\underset{\mathrm{from}\mathrm{SVD}}{\underbrace{{\mathbf{U}}_{\mathbf{B}}{\mathbf{\Sigma }}_{\mathbf{B}}{\mathsf{\Theta }}_{\mathbf{B}}^{\mathrm{H}}}}\Rightarrow \mathbf{B}\mathbf{F}\left(t\right)=\underset{{\mathbf{B}}_{\mathrm{pca}}}{\underbrace{{\mathbf{U}}_{\mathbf{B}}{\mathbf{\Sigma }}_{\mathbf{B}}}}\underset{{\mathbf{F}}_{\mathrm{pca}}\left(t\right)}{\underbrace{{\mathsf{\Theta }}_{\mathbf{B}}^{\mathrm{H}}\mathbf{F}\left(t\right)}}\cong {\tilde{\mathbf{B}}}_{\mathrm{pca}}{\tilde{\mathbf{F}}}_{\mathrm{pca}}\left(t\right),$$

  (18)

  where $\tilde{\square }$ denotes a reduced form including only the most important principal components. In the current approximation, $\tilde{P}=6\ll 16$ principal components were sufficient to represent ~99% of the information contained in $\mathbf{B}$. Thus, the final state space form was

  $$\left\{\begin{array}{l}\dot{\mathbf{x}}\left(t\right)=\mathbf{A}\mathbf{x}\left(t\right)+{\tilde{\mathbf{B}}}_{\mathrm{pca}}{\tilde{\mathbf{F}}}_{\mathrm{pca}}\left(t\right),\\ {\mathbf{F}}_{\mathrm{dyn}}\left(t\right)=\mathbf{C}\mathbf{x}\left(t\right),\end{array}\right.$$

  (19)

  having $\tilde{P}$ “orthogonal” inputs producing distinguishable outputs, whence they could be reconstructed by applying the Kalman filter. This was a key step in order to obtain a successful filter.
• The Upgraded Augmented Kalman Filter was obtained from the expanded state space form as

  $$\left\{\begin{array}{c}{\dot{\mathbf{x}}}_{\mathrm{e}}\left(t\right)=\underset{{\mathbf{A}}_{\mathrm{e}}}{\underbrace{\left[\begin{array}{cc}\mathbf{A}& {\tilde{\mathbf{B}}}_{\mathrm{pca}}\\ \mathbf{O}& \mathbf{O}\end{array}\right]}}\underset{{\mathbf{x}}_{\mathrm{e}}\left(t\right)}{\underbrace{\left[\begin{array}{c}\mathbf{x}\left(t\right)\\ {\tilde{\mathbf{F}}}_{\mathrm{pca}}\left(t\right)\end{array}\right]}}+\underset{{\mathbf{G}}_{\mathrm{e}}}{\underbrace{\left[\begin{array}{c}\mathbf{O}\\ \mathbf{I}\end{array}\right]}}{\mathbf{w}}_{proc}\left(t\right),\hfill \\ \mathbf{F}\left(t\right)=\underset{{\mathbf{C}}_{\mathrm{e}}}{\underbrace{\left[\begin{array}{cc}\mathbf{C}& \mathbf{O}\end{array}\right]}}{\mathbf{x}}_{\mathrm{e}}\left(t\right)+{\mathbf{w}}_{meas}\left(t\right),\hfill \end{array}\right.$$

  (20)

  where the input force ${\tilde{\mathbf{F}}}_{\mathrm{pca}}$ was incorporated into an expanded state space vector ${\mathbf{x}}_{\mathrm{e}}$, and it was artificially generated as the integral of a Gaussian process noise ${\mathbf{w}}_{\mathrm{proc}}$ of size $\tilde{P}\times 1$. A Gaussian measurement noise ${\mathbf{w}}_{\mathrm{meas}}$ of size $Q^{\prime }\times 1$ was also added to the output. $\mathbf{O}$ and $\mathbf{I}$ were null and identity matrices of adequate sizes.
• Process and measurement noises were described by the covariance matrices

  $$\left\{\begin{array}{lll}\mathrm{E}\left[{\mathbf{w}}_{\mathrm{proc}}{{\mathbf{w}}_{\mathrm{proc}}}^{\mathrm{T}}\right]& =& \left({10}^{10}÷{10}^{12}\right)\mathbf{I},\\ \mathrm{E}[{{\mathbf{w}}_{\mathrm{meas}}{\mathbf{w}}_{\mathrm{meas}}}^{\mathrm{T}}]& =& \left({10}^0÷{10}^2\right)\mathbf{I},\\ \mathrm{E}\left[{\mathbf{w}}_{\mathrm{proc}}{{\mathbf{w}}_{\mathrm{meas}}}^{\mathrm{T}}\right]& =& \mathbf{O},\end{array}\right.$$

  (21)

  where the covariances’ orders of magnitude were selected after a preliminary trial and error phase.
• A stationary Kalman filter was eventually determined as follows,

  $$\left\{\begin{array}{c}{\dot{\mathbf{x}}}_{\mathrm{e},\mathrm{obs}}\left(t\right)=\left({\mathbf{A}}_{\mathrm{e}}-{\mathbf{L}}_{\mathrm{e}}{\mathbf{C}}_{\mathrm{e}}\right){\mathbf{x}}_{\mathrm{e},\mathrm{obs}}\left(t\right)+{\mathbf{L}}_{\mathrm{e}}{\mathbf{F}}_{\mathrm{dyn}}\left(t\right),\hfill \\ {\tilde{\mathbf{F}}}_{\mathrm{pca}}\left(t\right)=\left[\begin{array}{cc}\mathbf{O}& \mathbf{I}\end{array}\right]{\mathbf{x}}_{\mathrm{e},\mathrm{obs}}\left(t\right),\hfill \end{array}\right.$$

  (22)

  where the observer (obs) gain matrix ${\mathbf{L}}_{\mathrm{e}}$ was computed by the classical lqe Matlab routine. Finally, the inverse PCA transformation matrix was applied in order to reconstruct the input force vector, as follows,

  $$\mathbf{F}\approx {\tilde{\mathsf{\Theta }}}_{\mathbf{B}}{\tilde{\mathbf{F}}}_{\mathrm{pca}}$$

  (23)

  and therefore the resultant force components $R_i$ could be derived by applying Equation (1). It is worth noting that the Upgraded Augmented Kalman Filter (UAKF) proposed here is practically a classic Augmented Kalman Filter based on a much more general model of dynamometer dynamics.

4. Filter Performance Evaluation

5. Conclusions