Long Short-Term Memory Neural Network for Temperature Prediction in
  Laser Powder Bed Additive Manufacturing
In context of laser powder bed fusion (L-PBF), it is known that the
properties of the final fabricated product highly depend on the temperature
distribution and its gradient over the manufacturing plate. In this paper, we
propose a novel means to predict the temperature gradient distributions during
the printing process by making use of neural networks. This is realized by
employing heat maps produced by an optimized printing protocol simulation and
used for training a specifically tailored recurrent neural network in terms of
a long short-term memory architecture. The aim of this is to avoid extreme and
inhomogeneous temperature distribution that may occur across the plate in the
course of the printing process.
  In order to train the neural network, we adopt a well-engineered simulation
and unsupervised learning framework. To maintain a minimized average thermal
gradient across the plate, a cost function is introduced as the core criteria,
which is inspired and optimized by considering the well-known traveling
salesman problem (TSP). As time evolves the unsupervised printing process
governed by TSP produces a history of temperature heat maps that maintain
minimized average thermal gradient.
  All in one, we propose an intelligent printing tool that provides control
over the substantial printing process components for L-PBF, i.e.\ optimal
nozzle trajectory deployment as well as online temperature prediction for
controlling printing quality.
Ashkan Mansouri Yarahmadi, Michael Breuß, Carsten Hartmann
Long Short-Term Memory Neural Network for Temperature Prediction in Laser Powder Bed Additive Manufacturing

###### Abstract

In context of laser powder bed fusion (L-PBF), it is known that the properties of the final fabricated product highly depend on the temperature distribution and its gradient over the manufacturing plate. In this paper, we propose a novel means to predict the temperature gradient distributions during the printing process by making use of neural networks. This is realized by employing heat maps produced by an optimized printing protocol simulation and used for training a specifically tailored recurrent neural network in terms of a long short-term memory architecture. The aim of this is to avoid extreme and inhomogeneous temperature distribution that may occur across the plate in the course of the printing process.

In order to train the neural network, we adopt a well-engineered simulation and unsupervised learning framework. To maintain a minimized average thermal gradient across the plate, a cost function is introduced as the core criteria, which is inspired and optimized by considering the well-known traveling salesman problem (TSP). As time evolves the unsupervised printing process governed by TSP produces a history of temperature heat maps that maintain minimized average thermal gradient.

All in one, we propose an intelligent printing tool that provides control over the substantial printing process components for L-PBF, i.e. optimal nozzle trajectory deployment as well as online temperature prediction for controlling printing quality.

Keywords:Additive manufacturing, laser beam trajectory optimization, powder bed fusion printing, heat simulation, linear-quadratic control

## 1 Introduction

In contrast to traditional machining, additive manufacturing (AM) builds objects layer by layer through a joining process of materials making the fabrication of individualized components possible across different engineering fields. The laser powder bed fusion (L-PBF) technique as an AM process, that we focus
on in this study, uses a deposited powder bed which is selectivity fused by a computer-controlled laser beam [17]. The extreme heating by the laser on the one hand, and on the other hand the influence of the degree of homogeneity of the heat distribution on the printing quality in L-PBF, make it highly challenging to conduct the printing process in an intelligent way that may guarantee high quality printing results. As explained in more detail when discussing related work, there has thus been a continuous effort to _(i)_ propose beneficial printing paths that help to avoid unbalanced heating and _(ii)_ to forecast the heat distribution in order assess the potential printing quality and terminate printing in case of foreseeable flaws.

In this paper, we propose to couple a laser beam trajectory devised on the basis of a heuristic control during the fabrication phase of L-PBF with prediction based on neural networks. The developed novel framework addresses both the abovementioned main issues in L-PBF and represents an intelligent printing tool that provides control over the printing process. To this end, we aim at conducting controlled laser beam simulation that approximately achieves _temperature constancy_ on a simulated melted power bed. In addition, we opt to perform temperature rate of change prediction as an important factor for microscopic structure of the final fabricated product.

The main novelty of the current paper is to adopt _long-short-term memory_ (LSTM) [8] prediction framework, which is introduced in Section 4 to predict the temperature distribution and its gradient during printing. This consequently can be used to avoid any overheating by taking necessary actions in advance, namely stopping the printing process to avoid the printer damage due to overheated deformed parts of the printing product. Based on this, we conjecture that our developed pipeline may provide a highly valuable step for practical printing that provides quality control of the printed product, while being efficient with regard to energy consumption and use of material. Finally, in Section 5, we present an effective numerical test concerning the predicted temperature gradients.

In Section 3 of this paper a simulation framework is brought out by recalling the heat transfer model together with a cost function that consists of two terms aiming to maintain almost a constant temperature with a low spatial gradient across the power bed area. For simplicity, we confine ourselves to a 2-dimensional domain, which is still a realistic description of printing over the manufacturing plate. In Subsection 3.2, the idea of the travelling salesman problem (TSP) as a heuristics for the laser beam steering is explained; being one of the most fundamental and well-studied NP-hard problems in the field of combinatorial optimization (e.g. [4, 7]), we will use a stochastic optimization strategy (simulated annealing) to establish an optimal laser trajectory.

## 2 Related work in Laser Powder Bed Additive Manufacturing

In general, a variety of different laser beam parameters such as laser power, scan speed, building direction and laser thickness influence the final properties
of the fabricated product. Due to intensive power of laser during additive manufacturing, the printed product can have defects, such as deviations from the target geometry or cracks caused by large temperature gradients. For example, inhomogeneous heating may lead to unmelted powder particles that can locally induce pores and microscopic cracks [6]. At the same time, the cooling process determines the microstructure of the printed workpiece and thus its material properties, such as strength or toughness, which depend on the proportion of carbon embedded in the crystal structure of the material [1].

In a broader view, machine learning approaches may be deployed to provide monitoring capabilities over varying factors of L-PBF, namely the used metal powder and its properties both at the initial time of spread and during the printing process as well as the laser beam parameters, aiming to investigate and avoid any defect generation during the fabrication. See [2] for an survey.

Concerning the powder properties, different capturing technologies along with machine learning tools are used to automate the task of defect detection and avoidance during the printing process. In [19, 20], the k-means clustering [14] and convolution neural network (CNN) [12] respectively, were used to detect and classify defects at the time of initial powder spread and their probable consequences during the entire printing phase and based on captured grey images. In [11], high resolution temporal and evolving patterns are captured using a commercial EOS M270 system to find layer-wise heat inhomogeneities induced by the laser. In [10], an inline coherent imaging (ICI) system was used to monitor the defects and unstable process regimes concerning the morphology changes and also the stability of the melt pools. Here, the back scattered intensities from the melt pool samples are measured as a function of their heights called A-lines. Later, a Gaussian fitting of individual A-lines is performed to determine centroid height and amplitude of melt pools as a function of time corresponding to a range of different stainless steel powders with different properties.

About the laser beam and its parameter optimization task one can avoid conducting expensive real experiments, in terms of material and power usage, by simulating the printing process by means of finite element method [5] (FEM), Lattice Boltzmann method (LBM) or finite volume method (FVM) See [3, 18] for extensive surveys. Later the gathered simulated data may be used in a data-driven machine learning approach within a L-PBF framework. In this context, a prediction task of thermal history was performed in [13] by adopting a recurrent neural network (RNN) structure with a Gated Recurrent Unit (GRU) in a L-PBF process. A range of different geometries are simulated by FEM while accounting for different laser movement strategies, laser power and scan speed. A three-dimensional FEM is adopted in [24] to simulate the laser beam trajectory and investigate its effects on the residual stresses of the parts. The simulation results show modifications of the residual stress distributions and their magnitudes, that was validated through experimental tests, as a result of varying laser beam trajectory type. A parametric study [25] used the same FEM simulation setup as [24] with three varying factors namely the laser beam speed, the layer thickness and the laser deposition path width. While each factor value varies in its range
from low, medium to high the hidden relations among the factors and their affects on residual stresses and part distortions are revealed.

In context of FEM simulation with a steering source of heat to represent the laser movement, one can refer to the work developed in [21]. Here, the residual stresses during the printing is predicted though the laser nozzle steering rule is not revealed.

## 3 Heat transfer model and TSP formulation

As indicated, we first describe our heat simulation setting which is the framework for the TSP optimization protocol described in the second part of this section.

### Heat simulation framework

We set up a simulation environment, namely _(i)_ a moving source of heat (cf. (3)) to act as a laser beam on _(ii)_ an area \(\Omega\subset\mathbb{R}^{2}\) simulated as deposition of aluminium metal powder called a plate. We assume that the plate is mounted to a base plate with large thermal conductivity, which makes the choice of Dirichlet boundary conditions with constant boundary temperature appropriate; if the surrounding is an insulator, then a reflecting i.e. zero-flux or von Neuman boundary condition is more suitable. A sequence of laser beam movements, called a trajectory, is followed so that at each point the heat equation (1) is resolved based on FEM providing us a temperature map that varies on different plate locations as the time evolves.

Letting \(u\) be the temperature across an open subset \(\Omega\subset\mathbb{R}^{2}\) as time \(t\) evolves in \(\in[0,T]\), the heat equation that governs the time evolution of \(u\) reads

\[\frac{\partial}{\partial t}u(x,y,t) =\alpha\nabla^{2}u(x,y,t)+\beta I(x,y)\,, (x,y,t)\in\Omega^{\circ}\times(0,T) \tag{1a}\] \[u(x,y,t) =\theta_{0} (x,y,t)\in\partial\Omega\times[0,T]\] (1b) \[u(x,y,0) =u_{0}(x,y) (x,y)\in\Omega \tag{1c}\]

where we denote by

\[\nabla^{2}\phi=\frac{\partial^{2}\phi}{\partial x^{2}}+\frac{\partial^{2}\phi }{\partial y^{2}} \tag{2}\]

the Laplacian of some function \(\phi\in C^{2}\), by \(\Omega^{\circ}\) we denote the interior of the domain \(\Omega\), and by \(\partial\Omega\) its piecewise smooth boundary; here \(u_{0}\) is some initial heat distribution, \(\theta_{0}\) is the constant ambient space temperature (\(20^{\circ}C\)), and we use the shorthands

\[\alpha\coloneqq\frac{\kappa}{c\rho}\quad\text{ and }\quad\beta\coloneqq\frac{1}{c\rho}\]

with \(\kappa\), \(c\) and \(\rho\) all in \(\mathbb{R}^{+}\), denoting thermal conductivity, specific heat capacity and mass density. Our power density distribution of choice to simulate a laser
beam is a Gaussian function :

\[I\left(x,y\right)=I_{0}\cdot\exp\!\left\{\left[-2\left[\left(\frac{x-x_{c}}{ \omega}\right)^{2}+\left(\frac{y-y_{c}}{\omega}\right)^{2}\right]\right]\right\} \tag{3}\]

with an intensity constant

\[I_{0}=\frac{2P}{\pi\omega^{2}} \tag{4}\]

by knowing \(\omega\) and \(P\) to be the radius of the Gaussian beam waist and the laser power, respectively. In our study we let \(x\in[-1,+1]\), \(y\in[-1,+1]\) with \((x,y)\in\Omega^{\circ}\), \(t\in\mathbb{R}^{+}\) and also \(u\left(x,y,0\right)=0\). The aluminium thermal properties are used to simulate the metal powder spread across the manufacturing plate.

We solved (1) using [15] by setting \(P=4200\left(\mathrm{W}\right)\) and \(\omega=35\) pixels while letting \((x_{c},y_{c})\) to take all possible trajectory points such that the domain \(\Omega^{\circ}\) always be affected by five consecutive heat source moves. In this way, we simulated the heat source movement across the board.

#### 3.1.1 Control objective.

By adoption of the TSP based protocol, we aim to minimize the value of a desired objective function :

\[J(m)=\frac{1}{2}\int_{0}^{T}\left(\int_{\Omega}|\nabla u_{m}(z,t)|^{2}+\left(u _{m}(z,t)-u_{g}\right)^{2}\mathrm{d}z\right)\mathrm{d}t \tag{5}\]

with

\[u_{m}=u_{m}(z,t)\,,\quad z=(x,y)\in\Omega^{\circ},\,t\in[0,T] \tag{6}\]

being the solution of the heat equation (1) on the interval \([0,T]\) under the control

\[m\colon[0,T]\to\Omega^{\circ}\,,\quad t\mapsto(x_{c}(t),y_{c}(t)) \tag{7}\]

that describes the trajectory of the center of the laser beam. Moreover, we have introduced \(u_{g}\) as the desired target temperature to be maintained over the domain \(\Omega^{\circ}\) as time \(t\) evolves.

The motivation behind (5) is to maintain a smooth temperature gradient over the entire plate for all \(t\in[0,T]\) as is achieved by minimizing the \(L^{2}\)-norm of the gradient, \(\nabla u\), while at the same time keeping the (average) temperature near a desired temperature \(u_{g}\) for any time \(t\).

We proceed, by dividing the entire plate \(\Omega^{\circ}\) into \(4\times 4\) sub-domains (see Fig. 1) and investigate our objective function (5) within each sub-domain as explained in Section 5.

### TSP-based Formulation

A common assumption among numerous variants of the TSP [4], as an NP-hard problem [7], is that a set of cities have to be visited on a shortest possible _tour_. Let \(\mathcal{C}_{n\times n}\) be a symmetric matrix specifying the distances corresponding to the paths connecting a vertex set \(\mathcal{V}=\{1,2,3,\cdots,n\}\) of cities to each other with
\(n\in\mathbb{N}\) to be the number of cities. A tour over the complete undirected graph \(\mathcal{G}\left(\mathcal{V},\mathcal{C}\right)\) is defined as a cycle passing through each vertex exactly once. The traveling salesman problem seeks a tour of minimum distance.

To adopt the TSP into our context, we formulate its input \(\mathcal{V}\) as the set of all \(16\times 16\) stopping points of the heat source over the board \(\Omega^{\circ}\), and the set \(\mathcal{C}\) as a penalty matrix with each element \(\mathcal{C}_{ij}\geq 0\) being the impact (i.e. cost) of moving the heat source from a \(i\in\mathcal{V}\) to \(j\in\mathcal{V}\). For every vertex \(i\in\mathcal{V}\), the possible movements to all \(j\neq i\) with the associated cost (5) is computed and assigned to \(\mathcal{C}_{ij}\) (see below for details). With this formulation, we want to remind the reader that \(\mathcal{C}\) elements are nonnegative and follow the triangle inequality :

\[\mathcal{C}_{ij}\leq\mathcal{C}_{ik}+\mathcal{C}_{kj} \tag{8}\]

with \(i,j,k\in\mathcal{V}\).

Note that, \(\mathcal{C}\) matrix is obtained based on a prior set of temperature maps produced using FEM without enforcing any particular protocol on them.

With this general formulation at hand, let us have a closer look at the discretized form of (5) that was used in current study to compute the elements of the penalty matrix:

\[\mathcal{C}_{ij}=\left|\sum_{l=1}^{4\times 4}\left(\left\|\boldsymbol{\Psi} \left(i,l\right)\right\|^{2}+\left(\boldsymbol{\Lambda}\left(i,l\right)-u_{g} \right)^{2}\right)-\sum_{l=1}^{4\times 4}\left(\left\|\boldsymbol{\Psi}\left(j,l \right)\right\|^{2}+\left(\boldsymbol{\Lambda}\left(j,l\right)-u_{g}\right)^{2 }\right)\right| \tag{9}\]

with \(l\) to be the sub-domain index. In addition, \(\boldsymbol{\Psi}\left(\cdot,l\right)=\sum_{z\in\Omega_{l}}\sum_{t\in t_{l}} \nabla u_{m}\left(z,t\right)\) represents the temperature gradient aggregation within each sub-domain, and

Figure 1: We divide the entire domain \(\Omega^{\circ}\) containing the diffused temperature values into \(4\times 4\) sub-domains separated by white lines. Within each sub-domain, (5) is computed to reveal how the temperature gradient \(\left|\nabla u_{m}(\cdot,t)\right|\) evolves as a function of time \(t\) evolves and how the average temperature \(\bar{u}\) is maintained near to a target value of \(u_{g}\). Note, the laser beam positions are in this image irrelevant.
\(\boldsymbol{\Lambda}\left(\cdot,l\right)=\frac{1}{|\Omega_{l}|}\sum_{z\in\Omega_{l}} \sum_{t\in t_{l}}u_{m}\left(z,t\right)\) is the average temperature value of each sub-domain, with \(t_{l}\) to be the time period on which the nozzle operates on \(\Omega_{l}\). Here, by \(|\Omega_{l}|\), we mean the number of discrete points in \(\Omega_{l}\subset\Omega^{\circ}\). In other words, (9) is the TSP cost of moving the nozzle from the \(i^{\text{th}}\) to the \(j^{\text{th}}\) stopping point that depends on (a) the mean square deviation of the temperature field from constancy and (b) on the mean square deviation from the global target temperature \(u_{g}\). In our simulation, the nozzle moves in the direction of the shortest (Euclidean) path connecting two successive stopping points. Thereby we assume the nozzle always adjusts its velocity so that the path between any arbitrarily chosen stopping points \(i\) and \(j\) always takes the same amount of time. The motivation behind this is to avoid heating up the entire domain \(\Omega^{\circ}\) as a result of keeping the nozzle velocity constant.

In practice no polynomial-time algorithm is known for solving the TSP [7], so we adopt a _simulated annealing algorithm_[23] that was first proposed in statistical physics as a means of determining the properties of metallic alloys at a given temperatures [16]. In the TSP context, we adopt [23] to look for a good (but in general, sub-optimal) tour corresponding to the movement of the heat source leading to the minimization of (9).

In Section 5, we reveal our prediction results obtained by adopting a TSP based heuristic along with the LSTM network. Before moving to the next section, let us observe a subset of temperature maps obtained based on TSP shown as Fig. 2.

Figure 2: A subset of heat maps produced by FEM as the solution to the heat equation (1). One clearly observes the effect of the previous laser positions on current status of the map, in terms of diffused temperature. The TSP as a heuristics steers the heat source across the plate aiming to keep temperature constancy. Note that all temperatures are in Celsius.
## 4 The LSTM Approach

Let us start discussion of our deep learning framework structure by investigating its LSTM [8] cell building blocks shown as Fig. 2(a) used to comprise a stack of three LSTM layers (see Fig. 2(b)) followed by a fully connected layer.

Here, we use temperature gradient values of \(\mu=14\) previous (i.e. from previous time) heat maps to predict the gradients values of the current heat map. By letting \(\zeta\) to be the current heat map, its history feature values formally lie in a range of \([\zeta-\mu,\zeta-1]\) heat maps with \(\zeta>\mu\). By considering each heat map to have 16 sub-domains and the same number of gradient features \(\Psi\left(\cdot,l\right)\), each corresponding to one sub-domain, we obtain in total \(\nu=\mu\times 16\) number of gradient feature history values that we vectorise to frame the vector \(\mathcal{X}\in\mathbb{R}^{\nu}\). Our aim is to use sub-sequences from \(\mathcal{X}\) to train the stacked of LSTMs and forecast a sequence of 16 number of gradient feature values corresponding to the sub-domains of a heat map of interest \(\zeta\).

Let us briefly discuss the weight and bias matrix dimensions of each LSTM cell. Here, we use \(q\in\mathbb{N}\) as the number of hidden units of each LSTM cell and \(n\in\mathbb{N}\) to represent the number of features that we obtain from FEM based heat maps and fed to the LSTM cell. More specifically, we have only one feature \(\Psi\left(\cdot,l\right)\) per sub-domain, i.e. \(n=1\). In practice, during the training process and at a particular time \(t^{\prime}\), a batch of input feature values \(\mathcal{X}\supset\mathcal{X}^{\left\langle t^{\prime}\right\rangle}\in \mathbb{R}^{b\times n}\) with \(b\in\mathbb{N}\) to be the batch size, are fed to each LSTM cells of the lowest stack level in Fig. 2(b). Here LSTM learns to map each feature value in \(\mathcal{X}^{\left\langle t^{\prime}\right\rangle}\) to its next adjacent value in \(\mathcal{X}\) as its label. The mapping labels are applied during the training and to the only neuron \(\mathcal{R}_{\eta}\in\mathbb{R}\) of the last fully connected layer with \(\eta=1\).

In addition to \(\mathcal{X}^{\left\langle t^{\prime}\right\rangle}\), each LSTM cell accepts two others inputs, namely \(h^{\left\langle t^{\prime}-1\right\rangle}\in\mathbb{R}^{b\times q}\) and \(c^{\left\langle t^{\prime}-1\right\rangle}\in\mathbb{R}^{b\times q}\), the so called the _hidden state_ and _cell state_ both of which are already computed at the time \(t^{\prime}-1\). Here, the cell state \(c^{\left\langle t^{\prime}-1\right\rangle}\) carries information from the intervals prior to \(t^{\prime}\).

A few remarks on how the cell state \(c^{\left\langle t^{\prime}\right\rangle}\) (10) at time \(t^{\prime}\) is computed by the formula (10) below are in order: A more precise look at (10) reveals that it partially depends on \(\Gamma_{f}\odot c^{\left\langle t^{\prime}-1\right\rangle}\) with \(c^{\left\langle t^{\prime}-1\right\rangle}\) being the cell state at the previous time \(t^{\prime}-1\). The term \(\Gamma_{f}\) satisfies

\[c^{\left\langle t^{\prime}\right\rangle}=\Gamma_{u}\odot\tilde{c}^{\left\langle t ^{\prime}\right\rangle}+\Gamma_{f}\odot c^{\left\langle t^{\prime}-1\right\rangle} \tag{10}\]

with \(\odot\) representing the element-wise vector multiplication. As (10) further shows, \(c^{\left\langle t^{\prime}\right\rangle}\) also depends on \(\tilde{c}^{\left\langle t^{\prime}\right\rangle}\) that itself is computed based on the feature vector \(\mathcal{X}^{\left\langle t^{\prime}\right\rangle}\) and the previous hidden state \(h^{\left\langle t^{\prime}-1\right\rangle}\) as:

\[\tilde{c}^{\left\langle t^{\prime}\right\rangle}=\tanh\left(\left[\left(h^{ \left\langle t^{\prime}-1\right\rangle}\right)_{b\times q}\right]\left( \mathcal{X}^{\left\langle t^{\prime}\right\rangle}\right)_{b\times n}\right] \times\mathcal{W}_{c}+\left(b_{c}\right)_{b\times q}\right) \tag{11}\]

with \(\times\) visualizing in this work standard matrix multiplication, and \(b_{c}\) and \(\mathcal{W}_{c}\) to be the corresponding bias and weight matrices, respectively.
Equation (10) contains two further terms, \(\Gamma_{u}\) and \(\Gamma_{f}\), called the _update gate_ and _forget gate_ defined as

\[\Gamma_{u}=\sigma\Bigg{(}\bigg{[}\Big{(}h^{\langle t^{\prime}-1\rangle}\Big{)}_{ b\times q}\bigg{]}\left(\mathcal{X}^{\langle t^{\prime}\rangle}\right)_{b \times n}\bigg{]}\times\mathcal{W}_{u}+\left(b_{u}\right)_{b\times q}\Bigg{)} \tag{12}\]

and

\[\Gamma_{f}=\sigma\Bigg{(}\bigg{[}\Big{(}h^{\langle t^{\prime}-1\rangle}\Big{)} _{b\times q}\bigg{]}\left(\mathcal{X}^{\langle t^{\prime}\rangle}\right)_{b \times n}\bigg{]}\times\mathcal{W}_{f}+\left(b_{f}\right)_{b\times q}\Bigg{)} \tag{13}\]

that are again based on the feature vector \(\mathcal{X}^{\langle t^{\prime}\rangle}\) and the previous hidden state \(h^{\langle t^{\prime}-1\rangle}\) with \(b_{u}\), \(b_{f}\), \(\mathcal{W}_{u}\) and \(\mathcal{W}_{f}\) to be the corresponding biases and weight matrices. Let us shortly conclude here, that the feature vector \(\mathcal{X}^{\langle t^{\prime}\rangle}\) and the previous hidden state \(h^{\langle t^{\prime}-1\rangle}\) are the essential ingredients used to compute \(\tilde{c}^{\langle t^{\prime}\rangle}\), \(\Gamma_{u}\) and \(\Gamma_{f}\), that are all used to update the current cell state \(c^{\langle t^{\prime}\rangle}\) (10).

The motivation to use the _sigmoid function_\(\sigma\) in structure of the gates shown in (12) and (13) is its activated range in \([0,1]\), leading them in extreme cases to

Figure 3: (a) A graphical representation of the LSTM cell accepting the hidden state \(h^{\langle t^{\prime}-1\rangle}\) and the cell state \(c^{\langle t^{\prime}-1\rangle}\) from the previous LSTM and the feature vector \(\mathcal{X}^{\langle t^{\prime}\rangle}\) at the current time. (b) A schematic representation of the adopted stack of LSTMs comprised of three recurrent layers processing the data. The upper LSTM layer is followed by a fully connected layer. The network performs a regression task being trained based on half-mean-square-error loss function (16). Note, the fully connected layer is established between the output of LSTM stack \(h^{\langle t^{\prime}\rangle}_{3}\) and the only neuron of the last layer \(\mathcal{R}^{\langle t^{\prime}\rangle}_{\eta}\) with \(\eta=1\).
be fully on or off letting all or nothing to pass through them. In non-extreme cases they partially contribute the previous cell state \(c^{\langle t^{\prime}-1\rangle}\) and the on the fly computed value \(\tilde{c}^{\langle t^{\prime}\rangle}\) to the current cell state \(c^{\langle t^{\prime}\rangle}\) as shown in (10).

To give a bigger picture, let us visualize the role of the \(\Gamma_{u}\) and \(\Gamma_{f}\) gates concerning the cell state \(c^{\langle t^{\prime}\rangle}\). In Fig. 3a, a direct line connecting \(c^{\langle t^{\prime}-1\rangle}\) to \(c^{\langle t^{\prime}\rangle}\) carries the old data directly from time \(t^{\prime}-1\to t^{\prime}\). Here, one clearly observes the \(\Gamma_{u}\) and \(\Gamma_{f}\) gates both are connected by \(+\) and \(\times\) operators to the passed line. They linearly contribute, as shown in (12) and (13), the current feature value \(\mathcal{X}^{\langle t^{\prime}\rangle}\) and the adjacent hidden state \(h^{\langle t^{\prime}-1\rangle}\) to update the current cell state \(c^{\langle t^{\prime}\rangle}\). Meanwhile, \(\Gamma_{u}\) shares its contribution through \(\times\) operator with \(\tilde{c}^{\langle t^{\prime}\rangle}\) to the passing line.

Finally, to make the current LSTM activated we need the cell state value at the time \(t^{\prime}\), namely \(c^{\langle t^{\prime}\rangle}\), that we obtain from (10) and also the so called _output gate_ obtained from

\[\Gamma_{o}=\sigma\Bigg{(}\bigg{[}\Big{(}h^{\langle t^{\prime}-1\rangle}\Big{)} _{b\times q}\bigg{|}\Big{(}\mathcal{X}^{\langle t^{\prime}\rangle}\Big{)}_{b \times n}\bigg{]}\times\mathcal{W}_{o}+\left(b_{o}\right)_{b\times q}\Bigg{)} \tag{14}\]

with \(\Gamma_{o}\in[0,1]\), and \(b_{o}\) and \(\mathcal{W}_{o}\) to be the corresponding bias and weight matrices. The final activated value of the LSTM cell is computed by

\[h^{\langle t^{\prime}\rangle}=\Gamma_{o}\odot\tanh\Big{(}c^{\langle t^{\prime }\rangle}\Big{)}. \tag{15}\]

Here, the obtained activated value \(h^{\langle t^{\prime}\rangle}\) from (15) will be used as the input hidden state to the next LSTM cell at the time \(t^{\prime}+1\).

Let us also mention that all the biases \(b_{c},b_{u},b_{f},b_{o}\in\mathbb{R}^{b\times q}\) and the weight matrices are further defined as

\[\mathcal{W}_{c}\coloneqq\Big{[}\left(\mathcal{W}_{ch}\right)_{q \times q}\!\!\big{|}\!\left(\mathcal{W}_{cx}\right)_{n\times q}\Big{]}^{\top}, \mathcal{W}_{u}\coloneqq\Big{[}\left(\mathcal{W}_{uh}\right)_{q \times q}\!\!\big{|}\!\left(\mathcal{W}_{ux}\right)_{n\times q}\Big{]}^{\top}\] \[\mathcal{W}_{f}\coloneqq\Big{[}\left(\mathcal{W}_{fh}\right)_{q \times q}\!\!\big{|}\!\left(\mathcal{W}_{fx}\right)_{n\times q}\Big{]}^{\top}, \mathcal{W}_{o}\coloneqq\Big{[}\left(\mathcal{W}_{oh}\right)_{q \times q}\!\!\big{|}\!\left(\mathcal{W}_{ox}\right)_{n\times q}\Big{]}^{\top}\]

leading both the \(\tilde{c}^{\langle t^{\prime}\rangle},c^{\langle t^{\prime}\rangle}\in \mathbb{R}^{b\times q}\).

Finally, we have a fully connected layer that maps the output \(h^{<:>}_{b\times q}\) of the stacked LSTM to the only neuron of the output layer \(\mathcal{R}_{\eta}\). This is achieved during the training process while the weight matrix \(\hat{\mathcal{W}}\in\mathbb{R}^{q\times\eta}\) and bias vector \(\hat{b}\in\mathbb{R}^{b\times\eta}\) corresponding to the fully connected layer are updated based on a computed loss \(\mathcal{L}\) using half-mean-square-error (16) between the network predictions and target temperature gradient values obtained from the heat maps produced by FEM.

\[\mathcal{L}=\frac{1}{2\eta b}\sum_{i_{1}=1}^{b}\sum_{i_{2}=1}^{\eta}\left(p_{i _{1}i_{2}}-y_{i_{1}i_{2}}\right)^{2} \tag{16}\]

Here, \(p\) and \(y\) values represent the predicted and the target gradient temperature values, respectively.
## 5 Results

To begin with, we consider a set of computed root-mean-square-error measures (RMSE) between the predicted and the target gradient values corresponding to the nozzle moves as shown in Fig. 4. More precisely, each curve value represents a computed RMSE between all \(4\times 4\) sub-domains gradient feature values of predicted heat map \(\zeta\) and their ground truth counterpart. Since we use a history of \(\mu=14\) previous gradient heat maps, the first prediction can be performed for the \(15^{\text{th}}\) nozzle move. Among all the measured RMSE values, we highlight four of them as can be seen in Fig. 4, that correspond to the \(25^{\text{th}}\), \(50^{\text{th}}\), \(75^{\text{th}}\) and \(100^{\text{th}}\) percentiles.

As one observes in Fig. 4, a relatively low RMSE measure is obtained almost across all nozzle moves on horizontal axis, though there exist some outliers. We further visualize the corresponding prediction results of the percentiles as Fig. 5. Specifically, let us take the \(25^{\text{th}}\) RMSE percentile computed between the black and its overlapping part of the pink curve shown in Fig. 4(a). The black curve in Fig. 4(a) is comprised of \(4\times 4\) forecasted vectorized gradient feature values of heat map sub-domains produced by \(54^{\text{th}}\) nozzle move with a RMSE equal to \(0.009\), compared to its overlapping pink curve. In this case, we let the \(i\) as the nozzle move number to range in \([\zeta-\mu,\zeta-1]\) to produce a history of feature gradient values corresponding to the heat map with \(\zeta=54\). This consequently means, the non-overlapping part of the pink curve in Fig. 4(a) represent the vectorized

Figure 4: Each curve value represents a computed RMSE between all \(4\times 4\) sub-domains gradient feature values of predicted heat map \(\zeta\) and their ground truth counterpart. The RMSE computation can be started from \(15^{\text{th}}\) nozzle move onward, since we use a history of \(\mu=14\) previous gradient maps. Those RMSE measures, highlighted as \(\times\) in ascending order correspond to the \(25^{\text{th}}\), \(50^{\text{th}}\), \(75^{\text{th}}\) and \(100^{\text{th}}\) percentiles, respectively.
history feature values of \(40^{\text{th}}\) to \(53^{\text{th}}\) heat maps that comprise \(\mu=14\) number of preceding heat maps of \(\zeta=54\), each of them with \(4\times 4\) sub-domains. The black curves in Figs. 5b, 5c and 5d are also comprised of the predicted gradient feature values of the heat maps \(\zeta\) equal to \(129\), \(220\) and \(244\), respectively, that are forecasted based on \(\mu\) number of their previous heat maps.

A closer look at the four prediction samples shown in Figs. 5 reveals that the even the \(100^{\text{th}}\) percentile, that marks some kind of outlier, is accurately predicted in that the shape of the black curve tracks the pink curve (ground truth). Concerning other three RMSE percentile values, the synchronicity among the black and its pink curves is preserved equally well though in some part we do not have a full overlap.

Finally, the parameters used during the training phase are revealed to be the Adam optimizer [9] applied on batch data of size \(6\). The epoch value is chosen to be \(350\) that results to a meaningful reduction of RMSE and loss measures within each batch. The initial learning rate was also chosen to be \(0.008\) with a drop factor of \(0.99\) concerning each \(12\) epochs. To avoid the overfitting phenomenon, the flow of data within the network structure is randomly adjusted by setting the LSTMs outputs with a probability of \(0.25\) to zero [22].

## 6 Conclusion

We developed a novel and practical pipeline and mathematically justified its comprising components. Our proposed model consists of two major components, namely the simulation part of a laser power bed fusion setup based on FEM and an intelligent agent based on LSTM network that actively judges the simulation results based on a proposed cost function. The FEM simulation can be robustly applied before conducting expensive real-world printing scenarios so that the intelligent component of the pipeline can decide on early stopping of the printing process. The LSTM based network predicts the forthcoming temperature rate of change across the simulated power bed based on previously seen temperature history leading us to have a means of control to achieve a final optimal printing process as visualized by our obtained results.

## Acknowledgements

The current work was supported by the European Regional Development Fund, EFRE 85037495.
