An optimization approach for a milling dynamics simulation based on Quantum Computing

Since the machining of complex aerospace components, like integral compressor-rotors (blade integrated disks), is very cost-intensive, optimizing the process by means of process simulations is an active field of research. With the rise of Quantum Computing, a new instrument with high optimization potential is moving into focus. In this paper, a possible application of Quantum Computing for the machining simulation of multi-axis milling of thin-walled aerospace components is discussed. For this reason, a simulation framework used for the milling simulation is analyzed and each component is evaluated separately in relation to Quantum Computing. Parts of the Harrow, Hassidim, and Lloyd algorithm are proposed to enhance the Finite-Element simulation-based component, like the modal analysis for dynamics simulation. This algorithm can solve linear system problems with exponential speed-up over the classical method. The paper presents a roadmap on how the classical steps of a modal analysis for dynamics simulation could be replaced by quantum algorithms based on quantum phase estimation. The implementation of the first working steps is presented to validate this approach. The linear system problem, arising from the dynamics simulation, is analyzed in detail and a minimal value problem of linear coupled oscillators is derived.


Introduction
In manufacturing applications, the potential of models and simulations for optimization based on digital twins (DT) is not fully exploited due to their computational requirements and the expert knowledge needed to operate them.Therefore, relevant physical effects in industrial practice are either neglected or only approximated by rough estimates.As a result, the quality of the DT and the insights and decisions derived from this suffer significantly, which leads in many cases to significant economic disadvantages.[1] In this context, the Fraunhofer IPT has been developing a domain-specific DT framework for machining, called dPart ® , for several years.It comprises data acquisition, processing, and analysis methods embedded in an Industrial-Internet-of-Things (IIoT) infrastructure (Fraunhofer Edge Cloud) for the targeted analysis and optimization of machining processes.[2] The following investigated simulation is embedded in the mentioned DT framework.With the advent of Quantum Computing (QC), a new tool with high optimization potential comes into focus.Mathematically it was shown that simulations with given accuracy allow polynomial and in certain cases even exponential acceleration (based on the Harrow, Hassidim, and Lloyd (HHL) algorithm) [3].A potential application for QC in manufacturing simulation is presented in the machining simulation of multi-axis milling of thin-walled aerospace components like integral compressor-rotors (blade integrated disks (blisks)).Blisks represent one of the most challenging components in turbomachinery manufacturing [4].Due to the high aspect ratio of the blade geometry resulting in low workpiece rigidity, various effects such as deflection and chatter can occur during the machining process [5].If vibrations occur during machining, it is often difficult for the operator to locate the cause.To avoid vibration issues, very conservative process parameters, such as reduced feed, cutting speed, cutting depth, and cutting width, are often used.The selection of an advantageous spindle speed for the present machining case is established based on stability diagrams.For the prediction of the process stability and thus for the generation of the stability diagrams, information about cutting force, dynamics of the total system (superposition of tool and workpiece dynamics), process parameters, and tool geometry is required.[6] In this paper, the workflow of the coupled technology models used in process simulation is presented in section 2.Then, two possible applications of QC out of the simulation workflow in subsection 2.1 and 2.2 are extracted.In section 3, the approach for enhancing the dynamics simulation of multi-axis milling of thin-walled aerospace components with quantum algorithms is described on the basis of linear coupled oscillators in detail.Moreover, in section 4 and 5 the quantum emulation as proof of concept is summarized so far and reveals the roadmap with the next steps to apply in the process design and optimization for an industrial thin-walled aerospace component.

Milling dynamics simulation workflow
The FEM simulation workflow (Figure 1), originally published by Rudel et.al., is briefly summarized below for a better understanding.It shows the detailed sequence of inputs, methods, and variables involved in calculating the virtual representation of the dynamic process behavior.The user programs the toolpath for the milling operation in the first step.Necessary inputs for the dexel-based engagement simulations are the NC path, the stock geometry of the workpiece, and the tool geometry such as diameter D, number of teeth N th , and helix angle λ.The NC path contains all relevant cutting parameters like axial stepover a p , radial stepover a e , feed per tooth f z , and spindle speed n.It is also necessary to select discrete Cutter Locations (CLs) for which the Finite-Element (FE) modal analysis should be performed in the following steps.Since a single toolpath consists of several thousand points, the simulation of all CLs of the toolpath is not realistic in terms of computational time and intensity.[6] In order to calculate the engagement angles between tool and workpiece on the one hand, and to extract CLdependent In-Process-Workpieces (IPWs) on the other hand, a dexel-based engagement simulation was then executed.Moreover, the user needs to define material properties, e.g.density ρ; Young's Modulus E; Poisson's ratio ν, and boundary conditions, e.g.constraints or friction; in order to run FE modal analysis and obtain CL-dependent Frequency Response Functions (FRFs).In addition to the number of selected CLs, for the FE modal analysis of the selected CL resolution, the computational effort also scales exponentially with the number of FE mesh nodes, the selected node resolution, and model size [7,8].Additionally, it is also decisive how many of the relevant eigenmodes are considered during the FE modal analysis.In the presented case, the simulation was limited to the first six eigenmodes [9].Furthermore, specific cutting force coefficients need to be defined to execute a dual-mechanistic cutting force simulation and to obtain the process-related cutting force F c .Fig. 1.FEM simulation workflow based on Rudel et.al. [6] Knowing the position-dependent cutting force and FRFs, stability diagrams can be generated to select suitable spindle speeds to avoid vibrations during the milling process.[6,9,10] The FEM simulations workflow (Figure 1) includes analytical and numerical models.Determination of the relative time required for each step, for a finishing operation at a simplified blade geometry (Figure 2 (a)) showed that the analytical models were negligible in terms of their computing time (see Table 1, computer specifications: Intel ® Xeon ® Gold-5118 2.30 GHz CPU, 24 cores, 380 GB RAM).Note that the computational time of simulation tasks related to the IPW, e.g. the FE workflow, highly scales with the defined number of cutter locations.As an example, selecting 200 CLs for one operation results in a time effort of 4400 seconds.Both analytical methods, force and stability simulation in the milling dynamics simulation workflow, are less computationally extensive compared to the dexel-based engagement simulation and the FE workflow.This disadvantage of numerical simulation has been elaborated on in previous scientific works [11].For the use case of finishing in blisk milling, the volume of the removed chip is very small.High resolution is required to calculate the engagement precisely leading to a very small dexel size, which negatively affects the computational time.However, with the help of suitable quantum algorithms (e.g. the HHL-Algorithm), it is assumed that the computation time will be exponentially accelerated compared to classical algorithms on classical computers [12].Analogous behavior applies to the FE modal analysis.Due to the complexity of numerical methods, only discrete cutter locations are currently modeled and considered in the simulation of manufacturing processes.Both methods show the potential to bring advantages in the field of milling dynamics simulation by using quantum algorithms.The two methods are further discussed in the following subsections 2.1 and 2.2.

Dexelsimulation
Within the dPart ® framework, the engagement simulation is a dexel-based simulation.The main components of a dexel model are a structured two-dimensional grid on a plane that serves as a nail coordinate system (1), nails, perpendicular to the coordinate system's plane located at the intersection points of the grid (2), and a subset defined by the intersection points between the nail and the volume that is represented by the nails (3).A three-dimensional volume is discretized using three orthogonal depth fields aligned with the axes of a cartesian coordinate system, an example is shown in Figure 2 (b).The engagement simulation provides the tool's contact area, contact length, and width.It also allows tracking removed volume and tool cutting direction to enable an evaluation of the workpiece cut by cut.The algorithm calculates all segments of any dexel (for all three grids) in the nail model that need to be changed.The intersections are performed individually for each of the three dexel grids to take advantage of memory closeness.The calculations performed by the intersectors are based on floatingpoint calculations.As stated above, the engagement simulation was investigated for optimization potential.In the investigated use case, the order of magnitude between tool size and cutting precision is very different (a big tool and compared to a small nail model), leading to a vast amount of intersection calculation per cut.A speed-up was reached using parallelization concluding the bottleneck remaining in access operations and therefore out of scope for quantum-based optimization.

Classical modal analysis
The FE workflow is performed for each IPW in subsequent calculations.Starting after the dexel-based engagement simulation in CAM, the main steps are IPW conversion from dexel to solid, meshing process, and the execution of modal analysis.Considering the condition for quantum acceleration, the workflow was analyzed step by step both in terms of computational time and with regards to extract numerical problems that can be transformed in quantum algorithms.The conversion steps from the dexel representation to the solid geometry and the meshing process consumed the most time.However, the conversion and meshing processes are highly affected by internal software algorithms and were therefore not focused subsequently.In contrast, the step 'modal analysis', i.e. solving the eigenvalue problem, is a well-known mathematical numerical problem and therefore found to be suitable for computation with quantum algorithms.Within the classical step 'modal analysis', the equation of motion for a free and undamped system is consulted in the physical domain on a local simulation computer.Here, the mass matrix and the stiffness matrix are extracted for each respective IPW.With the equation of a harmonic oscillation and the second derivative of the harmonic oscillation, the equation of motion can be transformed into an eigenvalue problem.By knowing eigenvalues and eigenvectors (and damping ratio), the dynamics for each simulated IPW can be described in the state space representation, usually visualized as FRF.The results can be further utilized for the calculation of stability diagrams for optimization of the spindle speed to avoid process vibrations.In the future approach, the quantum-based modal analysis will be connected to a QC via a new microservice in the batch layer, within the Core Domain of the framework dPart ® [6].It is planned, that the calculated solutions can be sent back to the local simulation computer via the mentioned microservice.There, the results can be further utilized as described above for optimization of the spindle speed to avoid process vibrations.

Quantum based modal analysis
In the beginning, the calculation of a simplified blisk blade geometry was started.The selected geometry and simulation parameters are shown in Figure 2 (a) and Table 2. Despite the highly simplified geometry, the matrix size of 1005 by 1005 was not computable with existing quantum hardware.Therefore, the solution of the modal analysis and the eigenvalue problem using quantum algorithms was investigated within this paper based on a simplification approach of linear coupled oscillators.The corresponding matrix was of size 2 N × 2 N , where N is the number of oscillators.However, this approach can be scaled up again as soon as the appropriate quantum hardware is available.

Linear coupled oscillators
The system of linear coupled oscillators (Fig. 3) was studied with the aim of simulating its behavior under the influence of an external force F = ( f 0 , . . ., f N−1 ) T .This is described by its equation of motion (Eq.1), where X = (x 0 , . . ., x N−1 ) T contains the displacement x i of the N oscillators, M is the mass matrix, K the stiffness matrix, and C the damping.
This equation can be solved by switching into frequency space using X(t) = X(ω) exp(iωt)dω, thus resulting in where the frequency dependent response function χ(ω) between the amplitude of the applied force F(ω) and the displacement X(ω) is introduced with the eigenfrequency ω k and the eigenvector V k for the matrix M − 1 2 KM − 1 2 .The damping is represented by the two coefficients α and β, based on Rayleigh Damping and depending on experimental measurements [13].A full description of the response function depends on the accurate computation of the eigenfrequencies and eigenvectors.Those can be computed using classical algorithms within polynomial time (O(N 3 ) with the size of the Matrix N [14]).The big drawback of those algorithms is, that they have to deconstruct an equation of vectors into their components and solve those one at a time.A quantum system can be described directly with linear algebra, which makes it optimal for linear algebra computations Fig. 3. Linear coupled oscillators as a toy model.The system is fully described by the displacement x i of the i-th oscillator, its external forces f i , its mass m i and the stiffnesses {k i , k i+1 } of the surrounding connections.in which quantum states represent vectors and quantum operations or gates represent matrices.However, efficient methods are needed to prepare those quantum states and gates without losing any speed-up.One example of quantum algorithms that suits the problem described above is the quantum phase estimation (QPE), which can be used to estimate eigenvalues and eigenvectors and promise exponential speed-up in comparison to classical algorithms [15,16].This speed-up is only given for a certain family of problems, as it depends on assumptions.This exact algorithm is further described in subsection 3.2.

Quantum phase estimation
QPE is a central quantum algorithm that was introduced by Kitaev [15].Its main function is to find the eigenvalue of a unitary operator U, which finds use in many physics [17] and chemistry [18] problems as well in Shor's algorithm [19].In the standard setting for QPE, an eigenstate |Ψ⟩ of the unitary U is given and, in addition, the unitary U is given in the form of a quantum gate.More precisely, the goal is to estimate the phase φ ∈ [0, 1) of the eigenvalue e 2πiφ associated with |Ψ⟩.For a classical algorithm, this task scales exponentially with respect to the number of qubits necessary to represent U |Ψ⟩.
It is further possible to use QPE to determine the eigenvalues of a hermitian matrix H if there exists a quantum circuit for U = e 2πiH .This is not always possible without losing all quantum speed-up.However, QPE can be used without losing the quantum advantage as long as the matrix can be decomposed in a linear superposition of elementary quantum gates [20] or if the rule of how to compute the entries of the matrix is known 4. The QPE starts with the preparation of the eigenstate Ψ in the lower register.This is followed by a series of controlled U, which prepares the desired phase φ.The inverse QFT probes for any possible phase, such that the eigenvalues φ and true eigenstates |Ψ⟩ are stored with high probability in the two quantum registers.[21].The central element of the original version is the quantum Fourier transform (QFT) [22,23,24,25] as shown in Fig. 4. For a detailed description of the QPE algorithm, see the studies of Cleve et al. (1998), Nielsen et al. (2010), and Scarani (2012) [16,25,26].
In its original version, QPE needs a register of m ancilla qubits in which it stores the phase with high probability.This phase can only be exact up to a finite error ε < 2 −(m+1) due to its binary representation [25].The probability distribution in the ancilla register is peaked around the desired phase φ and its width increases with the error ε.The minimum probability to measure the correct phase never falls below 4/π 2 [16].The time efficiency of the algorithm depends mostly on the number of queries of U, which scales as O(ε −1 ) while the number of additional operators in the QFT scales as O( log(1/ε) 2 ).Hence, the algorithm scales inversely proportional with respect to the error ε.Notice that this is the scaling obtained by assuming that the unitary U k is implemented by applying the unitary U k times.The scaling further reduces if there exists a quantum circuit for the direct implementation of U k which scales similarly, in the number of gates, as U.This makes it an algorithm that scales in the number of U k calls logarithmically with the inverse of the error O(log(1/ε)) and hence, the QFT dominates the computation time [16].
A basic assumption in QPE is that one is capable of implementing the unitary operator U and preparing a state Ψ , which overlaps significantly with the eigenstate |Ψ⟩, within polynomial time on a quantum device.As already mentioned, this is not always the case for arbitrary U = e iH .However, if H can be decomposed into a superposition of elementary gates, one can use a truncated Taylor Series expansion to bypass the exponential function [20].Another case in which the implementation of the unitary U is efficient is when the matrix H has sparsity properties and its entries can be efficiently evaluated [27].The preparation of the eigenstate |Ψ⟩ is also a non-trivial task that needs to be accomplished efficiently in the number of quantum gates to preserve the speed-up of the algorithm.Methods to prepare an approximate eigenstate efficiently are known [28,29].Here, a good approximation of Ψ is sufficient due to the probabilistic nature of the algorithm, which suppresses small errors.The desired eigenvalues φ are hidden in the ancilla register as superposition (see Fig. 4).Extracting them is a non-trivial problem, as a measurement destroys the superposition and yields only one eigenvalue with a certain probability.QPE amplifies the probability of measuring the true eigenvalues, which can be picked out by repeating the whole process multiple times.The corresponding true eigenvectors |Ψ⟩ are in the bottom register (see Fig. 4) and can be further used in other algorithms such as the Hadamard test [30] or extracted by measuring it.

Quantum emulation as proof of concept
As proof of concept, we emulated the QPE for a system of N = 4 linear coupled oscillators of equal mass m = 4 kg, which are connected by springs with stiffness k = 1 N m −1 , on a classical computer.The resulting FRF is shown in Fig. 5  Frequency ω (Hz) arbitrarily chosen pair of input and output nodes (2 and 0).In here, the damping was neglected.The blue graph represents the result of the numerically solved eigenvalues, eigenvectors of the matrices and vectors, and their FRF.Conversely, the orange one shows the previous emulation of the FRF based on the presented circuit in Qiskit on a classical computer.The accuracy of the quantum emulation depends directly on the number of qubits and with it on the computation-time of the algorithm.To classify the readiness of the individual simulation steps, it should be noted that the steps of state preparation, gate preparation, and data extraction are implemented, but work inefficiently in terms of computation time.An improved implementation is currently being worked on.

Summary & Outlook
The presented paper introduces an optimization approach for a dynamics milling simulation based on QC algorithm.The FEM simulation workflow of the dynamics simulation of multiaxis milling of thin-walled aerospace components, like integral compressor rotors, was presented, and individual steps of it were evaluated for quantum readiness.The FE modal analysis was selected as a suitable use case.A scalable quantum algorithm utilizing QPE was presented based on the minimum value problem of linear coupled oscillators.Nevertheless, the state-of-the-art quantum devices, so-called noisy intermediate scale quantum (NISQ) devices, are characterized by non-trivial qubit numbers (∼ 100), but also by relatively high error rates between 10 −2 and 10 −3 in the implementation of gates and measurements [31].This limits the depth of the quantum circuit that can be simulated on current quantum hardware without having a large error on the output to a few thousand gates [32].QPE generally requires high-depth circuits for exploiting its full potential, and as such it is not thought to be an algorithm that is suited for NISQ devices, but it requires a full fault-tolerant QC [33].However, further developed algorithms can be used to es-timate the phase even on NISQ technology using the evolution operator for various times [34,35].
Continuous further development of the presented quantum algorithm will take place.The main emphasis will be on troubleshooting and debugging specific parts of the algorithm.Additionally, efforts will be directed toward scaling up the minimum value problem of linear coupled oscillators to a threedimensional approach.The ultimate goal is to achieve a practical industrial use case, such as the milling dynamics simulation of integral compressor rotors like a blisk blade.Moreover, QC could support more simulations in the future in the manufacturing industry.For this purpose, a backend, to fill the interface between classical simulation and QC, will be developed in the dPart ® framework.

Fig. 5 .
Fig. 5. Numerically simulated and quantum emulated frequency response function for a system of 4 linear coupled oscillators with mass m = 4 kg and stiffness k = 1 N m −1 as proof of concept.

Table 1 .
Measured relative time effort for each simulation step

Table 2 .
Geometry and simulation parameters for an