# Batch Model Predictive Control for Selective Laser Melting

Riccardo Zuliani, Efe C. Balta, Alisa Rupenyan, John Lygeros  R. Zuliani, E.
C. Balta, A. Rupenyan and J. Lygeros are with the Automatic Control
Laboratory, ETH Zurich, 8092 Zurich, Switzerland. A. Rupenyan is also with
Inspire AG, 8092 Zurich, Switzerland
$\\{$ebalta,ralisa,lygeros$\\}$<EMAIL_ADDRESS>Research supported by NCCR Automation, a National Centre of Competence in
Research, funded by the Swiss National Science Foundation (grant number
180545).

###### Abstract

Selective laser melting is a promising additive manufacturing technology
enabling the fabrication of highly customizable products. A major challenge in
selective laser melting is ensuring the quality of produced parts, which is
influenced greatly by the thermal history of printed layers. We propose a
Batch-Model Predictive Control technique based on the combination of model
predictive control and iterative learning control. This approach succeeds in
rejecting both repetitive and non-repetitive disturbances and thus achieves
improved tracking performance and process quality. In a simulation study, the
selective laser melting dynamics is approximated with a reduced-order control-
oriented linear model to ensure reasonable computational complexity. The
proposed approach provides convergence to the desired temperature field
profile despite model uncertainty and disturbances.

##  I Introduction

Selective laser melting (SLM) is a promising additive manufacturing technique
used to construct parts with highly complex, customizable geometries [1]. In
SLM, a part is constructed by sequentially depositing a layer of metallic
powder on the bed of a building chamber, which is then melted using a high
intensity laser to form a solid part. This process is repeated at each layer
to form a 3D part [2]. Ensuring and maintaining the quality and the
repeatability of produced parts through closed-loop control is a major
challenge in additive manufacturing, in particular in SLM [3]. The
microstructure and the mechanical properties of a part produced by SLM depend
on its thermal history. However, these properties are not directly measurable
in-situ. The mechanical properties are often correlated with the temperature
distribution on the surface, or with geometrical properties of the melt pool
[4]. Furthermore, the layer-to-layer nature of the process poses additional
challenges, as defects in the microstructure remain sealed and in most cases
cannot be resolved after the corresponding layer is completed [5]. Thus, to
ensure good quality parts, a suitable temperature profile must be tracked at
each layer, which can be controlled through a suitable laser power profile [6,
7].

Iterative learning control (ILC) has been proposed as a method to control the
temperature profile over successive layers by layer-to-layer adjustments of
the laser power [8, 9]. ILC aims at improving the performance of a system that
operates repetitively by learning from previous iterations, and minimizing an
error from a defined control objective. One of the advantages of ILC is that
it does not require extensive knowledge of the system dynamics [10]. ILC is
thus suitable for controlling the additive manufacturing processes, while only
utilizing reduced-order models of their dynamics [11, 12, 8]. ILC relies on an
open loop application of an input sequence during an iteration, which is
updated only in between successive iterations. Consequently, it is able to
successfully compensate for static model uncertainties and for disturbance
signals that repeat identically in each iteration (or layer in the SLM). The
major drawback of existing ILC techniques applied to SLM is their inability to
reject non-repetitive disturbances. To solve this problem and account for non-
repetitive disturbances, _model predictive control_ (MPC) can be used in
combination with iterative learning control to adjust the control input during
an iteration by utilizing run-time measurement feedback [13, 14, 15].

Batch model predictive control (B-MPC) is a receding horizon control approach
proposed initially for chemical process control [16, 17]. In B-MPC input
updates are carried out at each time step by minimizing a quadratic next-
iteration cost, which applies a penalty to both the error (or, more precisely,
the prediction of what the error is going to be) and the input update signal.
A model-based prediction of the error within an iteration due to the online
control updates is used to formulate a finite horizon MPC problem at each time
instant for optimal control updates, where input constraints may also be
enforced. The online nature of B-MPC provides robustness to non-repetitive
disturbances via the feedback updates during each iteration, whereas the
estimation scheme keeps track of the learned repetitive disturbances by warm
starting the error prediction at the beginning of each new iteration. The
estimator therefore represents the iterative part of the B-MPC method.

In this paper, we adapt and apply the B-MPC method of [16] to control in the
temperature distribution in an SLM process, in the presence of both repetitive
and non-repetitive disturbances as our main contribution. The proposed method
is validated in simulation and its performance is compared against (i) a
proportional controller that uses the same estimation scheme and (ii) a
conventional model predictive controller that does not have the estimation
updates between iterations to adapt for repetitive disturbances.

The paper is organized as follows: Section II introduces the control problem,
the control-oriented model is presented in Section III, Section IV presents
the control approach, and Section V demonstrates the tracking performance of
the proposed B-MPC approach.

##  II Problem description

A conceptual illustration of a typical SLM process is shown in Fig. 1. A
metallic powder is melted by a laser at each layer and then left to solidify
to form a solid 3D object.

Figure 1: Concept of SLM process

In this work, we assume that the laser path is fixed and depends only on the
geometry of the part in construction. Additionally, following [8], we consider
a simplified thermal dynamics model of the process in which a layer directly
connected to the substrate at each iteration. With this approximation, the in-
layer process dynamics are identical for each layer, the dimensionality of the
model are fixed for each iteration $k$, and we have identical initial
conditions for each iteration representing a new layer. In other words, the
temperature distribution of a previous layer does not effect the subsequent
layer in this simplified model.

Our goal is to track a desired temperature field profile $T_{d}(\tau)$, where
$\tau\in[0,\tau_{f}]$ with $\tau_{f}$ as the layer print time, by adjusting
the laser power intensity $u(\tau)\in\mathbb{R}$. The desired temperature
field $T_{d}(\tau)$ is a time varying vector of temperatures for the whole
layer. For each layer $k$, we consider the temperature tracking control
problem

$\begin{split}\underset{u(\tau):[0,\tau_{f}]\to\mathcal{U}}{\mathrm{min.}}\quad&\int_{0}^{\tau_{f}}\|\mathbb{E}[T(\tau)-T_{d}(\tau)]\|_{Q(\tau)}^{2}d\tau,\\\
\text{s.t.}\quad&\dot{T}(\tau)=f(T(\tau),u(\tau),d(\tau)),\\\
\quad&T(0)=T_{0},\end{split}$ (1)

where $\mathcal{U}$ is a convex constraint set for the input,
$\mathbb{E}[\cdot]$ denotes the expectation of its argument, $T(\tau)$ is the
temperature state of the layer at time $\tau$, $\textstyle{d(\tau)}$ is a zero
mean i.i.d. multivariate Gaussian noise with covariance $\Sigma_{D}$, and
$f(\cdot,\cdot,\cdot)$ is a Lipschitz continuous function representing the
thermal dynamics of the layer. As we use a reduced order model of the process,
we include model uncertainty, in addition to process and measurement noise.
Additionally, we assume that the desired temperature field $T_{d}(\tau)$ is
identical for all layers (i.e., a fixed trajectory as a function of $\tau$ for
all layers and do not depend on the layer index $k$), and the initial
condition for each layer is identical, i.e., $\textstyle{T(0)=T_{0}}$, for all
$k$. Identical initial conditions is a reasonable assumption whenever there is
enough time between layers for sufficient cooling, or for cases when multiple
identical layers are printed sequentially over a build plate.

##  III Control-oriented Temperature Modeling for SLM

The temperature dynamics of an SLM process are nonlinear and complex,
involving phase transitions and temperature dependent material properties [1].
Finite element models, despite their fidelity, are not suited for control
purposes due to their computational complexity [18, 19]. In [8], it is shown
how the dynamics of the process below melting point can be well approximated
by a linear time-varying system. This linear model is better suited for
control purposes because of its simplicity and lower computational complexity.
According to [8], the powder layer can be represented as a grid of
$\textstyle{N=n_{x}\times n_{y}}$ nodes corresponding to a volume
$\textstyle{\delta_{x}\times\delta_{y}\times\delta_{z},}$ plus an additional
node representing the substrate, held at a constant temperature $T_{sub}$. The
heat transfer between nodes is expressed via a directed incidence matrix
$D\in\mathbb{R}^{N\times q}$, where $q$ is the number of links connecting the
nodes. Here, we assume the deterministic dynamics of $f$ in (1) is linear.
Using this formulation, the heat transfer equation, providing the thermal
evolution of the temperature of each node $T(\tau)\in\mathbb{R}^{N}$ in a
single layer, can be expressed as

$\displaystyle\begin{split}C_{c}\dot{T}(\tau)=&-DKD^{\top}T(\tau)+K_{sub}(T_{sub}-T(\tau))+\\\
&+B_{c}(\tau)u(\tau)~{}.\end{split}$ (2)

Here, $K=kI_{q}$ and $K_{sub}=k_{\text{sub}}I_{N}$ express the conductivity
between nodes and between nodes and substrate respectively, $C_{c}=cI_{N}$
describes the heat capacity of each node, and $I_{n}$ denotes the $n$
dimensional identity matrix. $B_{c}(\tau)\in\mathbb{R}^{N}$ describes, at each
time $\tau$, the fraction of the laser power $u(\tau)$ absorbed by each node
in the grid. Each node receives a different fraction of the laser power and,
overall, all the fractions must sum up to the laser transmission efficiency
$\textstyle{\eta\in(0,1]}$. Without loss of generality, we can set
$\sum_{i=1}^{n}B_{c,i}(\tau)=1$ for every $\tau$, i.e., take
$\textstyle{\eta=1}$ for simplicity. Equation (2) is a reasonably good
approximation of the SLM process provided that the melting point is not
reached and convection can be neglected.

Due to the modeling assumptions, the steady state temperature distribution is
the $T_{sub}$. Thus, we define all temperatures relative to $T_{sub}$ to
simplify the notation by absorbing the affine terms. By setting
$x(\tau)=T(\tau)-T_{sub}$ we can reformulate (2) into a continuous-time state
space form:

$\displaystyle\dot{x}(\tau)$ $\displaystyle=-A_{c}x(\tau)+B_{c}(\tau)u(\tau)$
$\displaystyle y(\tau)$ $\displaystyle=B_{c}(\tau)^{\top}x(\tau)\,.$

The output of the system is chosen as a weighted average of $x$ of the nodes
encompassed by the circumference of the laser beam. After discretizing the
system using exact discretization with sampling time $t_{s}$ we obtain the
following discrete-time state space system:

$\displaystyle\begin{split}x_{k}(t\\!+\\!1)&=Ax_{k}(t)+B(t)u_{k}(t)\\\
y_{k}(t)&=C(t)x_{k}(t)\,,\end{split}$ (4)

where $k$ refers to the layer, and $t$ is the time index of each iteration.
Assuming all laser scans have an identically finite duration $\tau_{f}$, we
have $t\in[0,\nu]$ with $\nu=\tau_{f}/t_{s}$ (assumed for simplicity to be a
positive integer). Thus, we utilize $y_{k}(t)$ for the minimization in (1)
written in $T(\tau)$, so that $Q(t)=C(t)^{\top}C(t)$. The error is defined as
$e_{k}(t)=y_{d}(t)-y_{k}(t)$, hence the objective term in (1), where
$y_{d}(t)=C(t)(T_{d}(t)-T_{sub})$ is the desired output. With input
disturbances and measurement noise present, system (4) becomes:

$\displaystyle\begin{split}x_{k}(t\\!+\\!1)&=Ax_{k}(t)+B(t)u_{k}(t)+B(t)v_{k}(t)\\\
y_{k}(t)&=C(t)x_{k}(t)+w_{k}(t)\,,\end{split}$ (5)

where we assume that $v_{k}(t)$ and $w_{k}(t)$ are Gaussian random variables
with zero mean and constant finite variance
$\textstyle{v_{k}(t)\sim\mathcal{N}\left(0,V\right)}$,
$\textstyle{w_{k}(t)\sim\mathcal{N}\left(0,W\right)}$, and are uncorrelated to
each other (both across iterations and across time). We construct a _lifted
representation_ of the system output as
$\mathbf{y}_{k}=[y_{k}(1)~{}y_{k}(2)~{}\dots~{}y_{k}(\nu)]^{\top}$. The same
can be done for the input, the noises, and the error. Assuming, without loss
of generality, that $x_{k}(0)=0$ for all $k$ (since any nonzero initial
condition can be seen as a disturbance), we have:

$\displaystyle\mathbf{y}_{k}$
$\displaystyle=G\mathbf{u}_{k}+G\mathbf{v}_{k}+\mathbf{w}_{k}$ (6a)
$\displaystyle\mathbf{e}_{k}$
$\displaystyle=\mathbf{y}_{d}-\mathbf{y}_{k}=\mathbf{e}_{k-1}-G\Delta\mathbf{u}_{k}+\bar{\mathbf{w}}_{k}\,,$
(6b)

where $\Delta\mathbf{u}_{k}=\mathbf{u}_{k}-\mathbf{u}_{k-1}$ is the input
update and

$G=\scalebox{0.7}{\mbox{$\displaystyle\begin{bmatrix}C(1)B(0)&\dots&0\\\
C(2)AB(0)&\dots&0\\\ \vdots&\ddots&\vdots\\\
C(T)A^{\nu-1}B(0)&\dots&C(\nu)B(\nu\\!-\\!1)\end{bmatrix}\,$}}.$

The term $\bar{\mathbf{w}}_{k}=-G\mathbf{v}_{k}-\mathbf{w}_{k}$ represents
layer-wise uncorrelated disturbances with
$\bar{\mathbf{w}}_{k}\sim\mathcal{N}(0,G^{\top}VG+W)$, i.e. disturbances that
have no correlation to disturbances in previous or next layers. We split the
error dynamics (6b) into its noise-free and noisy parts resulting in

$\displaystyle\begin{split}\bar{\mathbf{e}}_{k}&=\bar{\mathbf{e}}_{k-1}-G\Delta\mathbf{u}_{k}\\\
\mathbf{e}_{k}&=\bar{\mathbf{e}}_{k}+\bar{\mathbf{w}}_{k}\,.\end{split}$ (7)

In the next section, we propose a controller that iteratively solves (1) using
the formulation given in (7) by utilizing the Batch-MPC method proposed in
[16].

##  IV Proposed Batch MPC for SLM Processes

Our control approach leverages the repetitive nature of the process to
iteratively learn the optimal laser power profile. To account for model
uncertainty, a noise term that is correlated across the layers is introduced
in (7), next a Kalman filter is derived to construct an error estimate, and
finally, an optimization problem is formulated for the input update.

###  IV-A Dealing with Model uncertainty

The error dynamics (7) is derived from a simplified model (2), leading to a
mismatch between the complex physical process and its model used for control.
This mismatch leads to a systematic error that, unlike the stochastic noise
already present in (6) is correlated across iterations. To account for this
error, a noise term, $\bar{\mathbf{v}}_{k}$, that is correlated across layers
is introduced:

$\displaystyle\begin{split}\bar{\mathbf{e}}_{k}&=\bar{\mathbf{e}}_{k-1}-G\Delta\mathbf{u}_{k}+\bar{\mathbf{v}}_{k-1}\\\
\mathbf{e}_{k}&=\bar{\mathbf{e}}_{k}+\bar{\mathbf{w}}_{k}\,.\end{split}$ (8)

The statistical behaviour of $\bar{\mathbf{v}}_{k}$ depends on the accuracy of
the model $G$ as well as the magnitude of $\Delta\mathbf{u}_{k}$, hence, is
difficult to quantify. As suggested in [16], $\bar{\mathbf{v}}_{k}$ can be
modeled as an integrated white noise:

$\displaystyle\bar{\mathbf{v}}_{k}\sim\mathcal{N}(0,\bar{V})~{},~{}~{}~{}\bar{V}=\scalebox{0.7}{\mbox{$\displaystyle\begin{bmatrix}1&1&1&\cdots&1\\\
1&2&2&\cdots&2\\\ 1&2&3&\cdots&3\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
1&2&3&\cdots&\nu\end{bmatrix}$}}\sigma_{\bar{V}}^{2}\,.$ (9)

If $\bar{\mathbf{v}}_{k}$ is given as in (9), model (8) does not describe the
error dynamics exactly anymore, however, the statistical behaviour of the
noise term $\bar{\mathbf{v}}_{k}$ can be tuned by varying a single parameter
$\sigma^{2}_{\bar{V}}$. Guidelines on the choice of $\bar{\mathbf{v}}_{k}$,
which will impact the filter used to perform state estimation, are given in
[16].

With the model uncertainty, the statistical behaviour of
$\bar{\mathbf{w}}_{k}$ also becomes uncertain, hence, an additional tuning
parameter $\sigma_{\bar{W}}$ is introduced:

$\bar{\mathbf{w}}_{k}\sim\mathcal{N}(0,\bar{W})\>,\>\>\>\bar{W}=G^{\top}VG+W+I\sigma_{\bar{W}}^{2}\,.$
(10)

Note, that (10) does not use the same approach as in [16], instead we
leveraged the linear model knowledge ($G$) for approximating the covariance of
$\bar{\mathbf{w}}_{k}$.

To perform a state estimation at each time-step $t$, we partition $G$ as:

$G\triangleq[G(0),\>G(1),\>\cdots,\>G(\nu\\!-\\!1)],\>\>G(t)\in\mathbb{R}^{\nu}\,,$

let $\mathbf{e}_{k}(t)\in\mathbb{R}^{\nu}$ be the error sequence for the
$k$-th layer assuming that no input update is performed from time $t$ to time
$\nu$ (that is, until the end of the horizon), i.e., (6b) with

$\displaystyle\Delta{u}_{k}(t)=\cdots=\Delta{u}_{k}(\nu\\!-\\!1)=0\,.$

Then, an equivalent representation of system (8) in the in-layer domain, i.e.,
during the iteration, is given by the augmented system in (11):

$\begin{split}\begin{bmatrix}\bar{\mathbf{e}}_{k}(t\\!+\\!1)\\\
\mathbf{e}_{k}(t\\!+\\!1)\end{bmatrix}&=\begin{bmatrix}I&0\\\
0&I\end{bmatrix}\begin{bmatrix}\bar{\mathbf{e}}_{k}(t)\\\
\mathbf{e}_{k}(t)\end{bmatrix}-\begin{bmatrix}G(t)\\\ G(t)\end{bmatrix}\Delta
u_{k}(t)\,,\\\
e_{k}(t)&=\begin{bmatrix}0&H(t)\end{bmatrix}\begin{bmatrix}\bar{\mathbf{e}}_{k}(t)\\\
\mathbf{e}_{k}(t)\end{bmatrix}\,,\end{split}$ (11)

where we have

$H(t)=[~{}\underbrace{~{}0~{}~{}\cdots~{}~{}0~{}}_{t-1}~{}~{}1~{}~{}\underbrace{~{}0~{}~{}\cdots~{}~{}0~{}}_{\nu-t}~{}]\,.$

Notice that the output of (11) is the measured output of the SLM process
$e_{k}(t)$. The initial conditions of the state of the system (11), which also
convey the layer to layer dynamics of the system, are given in (12).

$\displaystyle\begin{bmatrix}\bar{\mathbf{e}}_{k}(0)\\\
\mathbf{e}_{k}(0)\end{bmatrix}\\!=\\!\begin{bmatrix}I&0\\\
I&0\end{bmatrix}\begin{bmatrix}\bar{\mathbf{e}}_{k-1}(\nu)\\\
\mathbf{e}_{k-1}(\nu)\end{bmatrix}+\begin{bmatrix}I\\\
I\end{bmatrix}\mathbf{\bar{v}}_{k-1}+\begin{bmatrix}0\\\
I\end{bmatrix}\mathbf{\bar{w}}_{k}\,.$ (12)

Notice that the impact of the noise terms $\bar{\mathbf{v}}_{k}$ and
$\bar{\mathbf{v}}_{k}$ is accounted for in the transition between successive
iterations in a lifted form. Therefore, the iteration-wise error states
$\bar{\mathbf{e}}_{k}$ and $\mathbf{e}_{k}$ capture the effect of
disturbances.

###  IV-B State estimation for Batch Model Predictive Control

The optimal estimation of the state of the system (11) can be achieved with a
Kalman filter. The standard estimation procedure involves repetitive
application, at each time step, of the prior update (13) and the measurement
update (14).

(13) (14)

Where $K(t)$ is the steady-state Kalman gain of the form given in [20, chap.
16]. The initialization at the beginning of each iteration is:

$\displaystyle\begin{bmatrix}\bar{\mathbf{e}}_{k}(0|0)\\\
\mathbf{e}_{k}(0|0)\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\bar{\mathbf{e}}_{k-1}(\nu|\nu)\\\
\bar{\mathbf{e}}_{k-1}(\nu|\nu)\end{bmatrix}\,,$ (15)

with $\mathbf{\bar{e}}_{-1}=\mathbf{y}_{d}$. The repetitive nature of the
process is exploited in (15). The estimator implemented here represents the
iterative part of the B-MPC method. At the beginning of each new iteration,
the error sequence $\mathbf{e}_{k-1}(\nu|\nu)$ containing non-repetitive
disturbances is discarded and the new error sequence for the next iteration is
initialized considering only the noise with iteration-wise correlation, i.e.
the repetitive disturbances. This way, the controller adapts to non-repetitive
disturbances in each layer through state estimation while keeping track of the
previously learned repetitive disturbances.

The tuning of the Kalman filter should consider the fact that by giving more
weight to $\bar{\mathbf{w}}$ in relation to $\bar{\mathbf{v}}$, e.g.
decreasing $\sigma^{2}_{\bar{V}}$, the change from one iteration to the next
one becomes smoother, thus making the algorithm less aggressive and more
robust to model uncertainty, at the cost of a lower learning rate. Increasing
the value of $\sigma_{\bar{W}}$ produces a less oscillatory in-layer
behaviour.

Notice that the term $\bar{\mathbf{v}}_{k}$ is necessary to force the filter
to trade off between prediction and measurement. If
$\textstyle{\bar{\mathbf{v}}_{k}=0}$ then the filter would simply rely on the
prediction without considering the measured error and this would prevent the
B-MPC algorithm from learning from experience.

###  IV-C Batch Model Predictive Control (B-MPC) Formulation

Here, we present the derivation of a _receding horizon MPC controller_ , which
computes the optimal input by minimizing a quadratic cost over a finite
horizon. The strength of the B-MPC control law comes from its ability to
integrate both _online disturbance rejection_ (using $\mathbf{e}_{k}(t|t)$)
and _iterative learning_ to reject the constant components of the disturbances
as well as the model uncertainty (using a non-repetitive noise-free version of
the error $\bar{\mathbf{e}}_{k}(t|t)$ in state estimation). Thus, the state
estimation scheme presented in Section IV-B is an essential part of the B-MPC
and it is a distinguishing feature compared to conventional MPC.

The control algorithm uses the receding horizon principle, i.e. at time $t$
the MPC minimizes the predicted cost in the next $m$ time steps by choosing an
optimal input update. The time horizon $m$ of the MPC must be selected to
ensure that the computational complexity is low enough for real-time
implementation of the algorithm. After computing the input, which is a
sequence of $m$ elements, only the first element is applied. At the next time
step $t\\!+\\!1$ the procedure is repeated.

The B-MPC algorithm approximates an optimal solution to (1) by solving an
optimization problem with quadratic objective and input constraints over a
decision horizon.

$\displaystyle\begin{split}\underset{\Delta\mathbf{u}_{k}^{m}(t)}{\mathrm{min.}}&\mathbf{e}_{k}^{\top}(t\\!+\\!m|t)Q\mathbf{e}_{k}(t\\!+\\!m|t)+\Delta\mathbf{u}_{k}^{m,\top}R\Delta\mathbf{u}_{k}^{m}(t)\\\
\mathrm{s.t.}&\quad\mathbf{u}_{min}\leq\mathbf{u}_{k}^{m}(t)\leq\mathbf{u}_{max}\\\
&\quad\delta\mathbf{u}_{min}\leq\mathbf{u}_{k}^{m}(t)-\mathbf{u}_{k}^{m}(t-1)\leq\delta\mathbf{u}_{max}\,,\\\
&\quad\mathbf{e}_{k}(t\\!+\\!m|t)=\mathbf{e}_{k}(t|t)-G^{m}(t)\Delta\mathbf{u}_{k}^{m}(t)\,,\\\
&\quad\mathbf{u}_{k}^{m}(t)=\mathbf{u}_{k-1}^{m}(t)+\Delta\mathbf{u}_{k}^{m}(t),\\\
&\quad\mathbf{u}_{k}^{m}(t-1)=[u_{k}(t-1)~{}\mathbf{u}_{k}^{m-1,\top}(t)]^{\top},\end{split}$
(16)

where, $\mathbf{u}_{k-1}^{m}(t)$ is the $m$ control inputs applied starting
from time $t$ in the previous layer, $u_{k}(t-1)$ is the most recent input
applied in the current layer, $\Delta\mathbf{u}_{k}^{m}(t)$ is the control
update in the next $m$ time steps ($t,t+1,\dots,t\\!+\\!m\\!-\\!1$):

$\Delta\mathbf{u}_{k}^{m}(t)=\begin{bmatrix}\Delta u_{k}(t)&\cdots&\Delta
u_{k}(t\\!+m\\!-1)\end{bmatrix}^{\top}\,,$

$\textstyle{\mathbf{e}_{k}(t\\!+\\!m|t)\in\mathbb{R}^{\nu}}$ is the prediction
of the output error at time $\textstyle{t+m}$ considering the future control
inputs from $t$ to $\textstyle{t\\!+\\!m\\!-\\!1}$, and
$\textstyle{G^{m}(t)=\begin{bmatrix}G(t)&\cdots&G(t\\!+\\!m\\!-\\!1)\end{bmatrix}}$.
$Q$ and $R$ are tunable positive semi-definite cost matrices. The problem (16)
is a quadratic program (QP) that can be solved efficiently. Let
$\textstyle{\bm{u}_{k}^{m,*}(t)}$ denote the optimizer of (16). Then the first
optimal input $\Delta u_{k}^{*}(t)$ of the optimal solution is used for the
input update, i.e., $\Delta u_{k}(t)=\Delta u_{k}^{*}(t)$, and the procedure
is repeated in a receding horizon. The control horizon $m$ is fixed as long as
$t+m\leq\nu$, otherwise we apply a shrinking horizon of $m=\nu-t$. The
resulting algorithm to implement the B-MPC is given in Algorithm 1.

Computation of the optimal input using (16) accounts for non-repetitive
disturbances during an iteration, while the state estimation through the
Kalman update keeps track of the learned repetitive disturbances throught the
error state. A block diagram of the system is presented in Figure 2.

Algorithm 1 Batch MPC

1:Init Compute $G,K(t)$ for $t\in[0,\nu]$

2: Initialize estimator state using $\mathbf{y}_{d}$

3:for $k=1:$ Maxit do

4: Initialize estimator state using (15)

5: for $t=1:\nu$ do

6: Solve QP given in (16) and compute $\Delta{u}_{k}(t)$

7: Apply the input: $u_{k}(t)=u_{k\\!-\\!1}(t)+\Delta u_{k}(t)$

8: Measure the error: $e_{k}(t)=y_{k}(t)-y_{d}(t)$

9: Kalman update using (13), (14)

10: end for

11:end for

B-MPC$+$$+$Plant$+$$v_{k}(t)$$\Delta
u_{k}(t)$$w_{k}(t)$$y_{k}(t)$Estimator$u_{k}(t)$$+$$e_{k}(t)$$\bar{\mathbf{e}}_{k}(t|t)$$y_{d}(t)$$\mathbf{u}_{k-1}=\left\\{u_{k-1}(0),\dots,u_{k-1}(N)\right\\}$
Data block $u_{k-1}(t)$

Figure 2: Block diagram of the B-MPC in closed-loop.

##  V Case Studies

We demonstrate the B-MPC approach on a simulated SLM process, with the
corresponding model and process parameters provided in Table I. A spiral-in
scan path (see Figure 3) has been chosen with a raster spacing of $60\mu m$.
The desired output $\textstyle{y_{d}(t),t\in[0,\nu]}$ has been computed on a
system with no model uncertainty by applying a constant power of $20$ W while
moving the laser along the diagonal of the grid with a constant velocity of
$0.5\,m/s$. The noise covariances $W,V$ have been chosen as $2\%$ of the
maximum magnitudes of the desired output and the input respectively.

TABLE I: Physical parameters used in the simulations Parameter | Symbol | Value
---|---|---
Dimensions of node | $\delta_{x},\delta_{y}$ | $2\cdot 10^{-5}\,[m]$
Powder layer thickness | $\delta_{z}$ | $5\cdot 10^{-5}\,[m]$
Conductivity with substrate | $k_{sub}$ | $20\delta_{z}\,[m]$
Heat capacity of a node | $c$ | $8.5\cdot 10^{-8}\,[J/K]$

The process model (5) is used for simulations. To emulate the effect of the
model mismatch, multiplicative uncertainty has been added to the diagonal
entries of the heat capacity matrix $C_{c}$, to the matrix describing the
conductivity with the substrate $K_{\text{sub}}$, and to each entry of matrix
$B(t)$, $t\in[0,\nu]$ in the simulated model. The random uncertainty is
sampled uniformly from the interval $[0,0.3]$ ($[-0.3,0]$ for the uncertainty
on $C_{c}$). For example, the heat capacity matrix is given by
$\textstyle{C_{c}=C_{c,\text{nominal}}(I+\text{diag}(\mathbf{r}))}$, where
$\mathbf{r}\in\mathbb{R}^{N}$ is a random vector whose entries are uniformly
sampled from the interval $[-0.3,0]$.

The input magnitude and rate constraints are defined as $u_{\text{max}}=20$ W,
$u_{\text{min}}=0$ W, and $\delta u_{\text{max}}=2$ W, $\delta
u_{\text{min}}=-2$ W, respectively, with the control horizon $m=20$. The
Kalman filter has been tuned using the heuristic guidelines given in [16]. The
resulting parameter values are
$\textstyle{\sigma_{\bar{V}}=0.8,\sigma_{\bar{W}}=70}$, and the corresponding
sample time is $\textstyle{T_{s}=10^{-5}\,s}$. The cost matrices $Q$ and $R$
are tuned heuristically.

We now compare the tracking performance of B-MPC to the performance achieved
with a proportional controller utilizing the estimation scheme explained in
Section IV-B. The input update commanded by the proportional controller at
every time-step is given by:

$\displaystyle\Delta u_{k+1}^{c}(t)=k_{p}H(t)\mathbf{e}(t|t)\,,$ (17)

The proportional gain $k_{p}$ has been determined heuristically, starting with
$\textstyle{k_{p}=0}$ and increasing until the system reaches instability.
Then, the gain has been set to the value that produced the smallest normed
error. $\Delta u_{k+1}^{c}(t)$ is then projected on the constraint set and
applied to the system.

Figure 3 shows the corresponding error norm $||\mathbf{e}_{k}||$ for the two
controllers as a function of the iteration index. The convergence rate of the
B-MPC is better than that of the proportional controller, and the steady state
error half of the one obtained with the proportional controller.

Figure 3: Comparison of the 2-norm errors of the B-MPC and the proportional
controller across different iterations.

In Figure 4 we compare the error norm across different iterations as obtained
by the B-MPC with the error norm obtained using a traditional MPC. The MPC
update rule is the same as in Algorithm 1, however, at the beginning of each
new iteration the input $\mathbf{u}_{k}$ is reset to the zero vector, and
$\mathbf{e}_{k}$ is reset to $\mathbf{y}_{\text{des}}$. The main difference
between B-MPC and MPC is that the latter does not benefit from the estimation
of the non-repetitive disturbances, since any information coming from previous
iterations is discarded once the iteration ends. Indeed, as can be seen in
Figure 4, the error norm obtained using the MPC does not improve over
successive iterations.

In Figure 4 we also compare the error norm across different iterations for
different choices of the parameter $\sigma_{\bar{V}}$. By decreasing
$\sigma_{\bar{V}}$ (e.g. to $0.1$), the KF relies more on the prediction,
which is subject to model uncertainty. As a result the convergence rate of the
B-MPC decreases and the error norm increases. If $\sigma_{\bar{V}}$ becomes
too large (e.g. equal to $25$), the KF relies too much on the noisy
measurement, as a result, the B-MPC produces more oscillatory input
trajectories, which lead to larger error.

Figure 4: Comparison of the 2-norm errors of the B-MPC for different values of
$\sigma_{\bar{V}}$ and with the MPC.

Figure 5 shows the error trajectories corresponding to the two controllers in
the final iteration, i.e., ${e}_{10}(t)$ with $t\in[0,\nu]$. The proportional
controller, which does not consider the constraints, produces larger errors
when the laser is in the vicinity of the corners of the spiral path. Notice,
however, that due to the tightness of the rate constraints, the B-MPC also
incurs errors along the corners of the spiral path.

Figure 5: Comparison of error trajectories in the final iteration.

Figure 6 shows the input applied in the final iteration, i.e., ${u}_{10}(t)$
with $t\in[0,\nu]$. Notice that the commanded input by the proportional
controller exceeds the input constraints and is saturated by the system. The
lower panel of Figure 6 shows the magnitude of the updates at each time step,
$\Delta{u}_{k}(t)$. The proportional controller, which does not consider the
constraints explicitly, does not always make full use of the entire range of
$[\delta\mathbf{u}_{min},\delta\mathbf{u}_{max}]$, due to the saturation.

Figure 6: Comparison of input trajectories in the final iteration.

##  VI Conclusion

In this paper we propose a B-MPC method for closed-loop temperature control of
an SLM process. The B-MPC controller enables both in-layer and between-layers
input updates, and achieves good performance even with non-repetitive
disturbance. After the end of an iteration, the information about the error
and the input of the system is the starting point for the next iteration. The
proposed controller is able to obtain high-performance tracking despite
repetitive and non-repetitive disturbances and model uncertainty. Furthermore,
the performance of the B-MPC outperforms the proportional controller in our
case study, both in terms of convergence rate and steady-state error.

The SLM process has been modeled using a reduced-order linear time-varying
model, which captures the dynamics below melting point. A possible extension
could consider a nonlinear model which takes into account the melt pool
dynamics for the process dynamics. Formal convergence analysis of the proposed
method with model mismatch and extensions of the control framework with multi-
layer models are subject for future work.

## References

  * [1] S. Bremen, W. Meiners, and A. Diatlov, “Selective laser melting: A manufacturing technology for the future?” _Laser Technik Journal_ , vol. 9, no. 2, pp. 33–38, 2012.
  * [2] C. Y. Yap, C. K. Chua, Z. L. Dong, Z. H. Liu, D. Q. Zhang, L. E. Loh, and S. L. Sing, “Review of selective laser melting: Materials and applications,” _Applied physics reviews_ , vol. 2, no. 4, p. 041101, 2015.
  * [3] V. Renken, A. von Freyberg, K. Schünemann, F. Pastors, and A. Fischer, “In-process closed-loop control for stabilising the melt pool temperature in selective laser melting,” _Progress in Additive Manufacturing_ , vol. 4, no. 4, pp. 411–421, 2019.
  * [4] P. Lott, H. Schleifenbaum, W. Meiners, K. Wissenbach, C. Hinke, and J. Bultmann, “Design of an optical system for the in situ process monitoring of selective laser melting (slm),” _Physics Procedia_ , vol. 12, pp. 683–690, 2011.
  * [5] S. Di Cataldo, S. Vinco, G. Urgese, F. Calignano, E. Ficarra, A. Macii, and E. Macii, “Optimizing quality inspection and control in powder bed metal additive manufacturing: Challenges and research directions,” _Proceedings of the IEEE_ , vol. 109, no. 4, pp. 326–346, 2021.
  * [6] G. Tapia and A. Elwany, “A review on process monitoring and control in metal-based additive manufacturing,” _Journal of Manufacturing Science and Engineering_ , vol. 136, no. 6, 2014.
  * [7] D. Liao-McPherson, E. C. Balta, R. Wüest, A. Rupenyan, and J. Lygeros, “In-layer thermal control of a multi-layer selective laser melting process,” in _2022 European Control Conference (ECC)_ , 2022.
  * [8] M. J. Spector, Y. Guo, S. Roy, M. O. Bloomfield, A. Maniatty, and S. Mishra, “Passivity-based iterative learning control design for selective laser melting,” in _2018 Annual American Control Conference (ACC)_. IEEE, 2018, pp. 5618–5625.
  * [9] A. Shkoruta, W. Caynoski, S. Mishra, and S. Rock, “Iterative learning control for power profile shaping in selective laser melting,” in _2019 IEEE 15th International Conference on Automation Science and Engineering (CASE)_. IEEE, 2019, pp. 655–660.
  * [10] H.-S. Ahn, Y. Chen, and K. L. Moore, “Iterative learning control: Brief survey and categorization,” _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ , vol. 37, no. 6, pp. 1099–1121, 2007.
  * [11] I. Lim, D. J. Hoelzle, and K. L. Barton, “A multi-objective iterative learning control approach for additive manufacturing applications,” _Control Engineering Practice_ , vol. 64, pp. 74–87, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0967066117300692
  * [12] Z. Wang, C. P. Pannier, K. Barton, and D. J. Hoelzle, “Application of robust monotonically convergent spatial iterative learning control to microscale additive manufacturing,” _Mechatronics_ , vol. 56, pp. 157–165, 2018.
  * [13] S.-K. Oh and J. M. Lee, “Iterative learning model predictive control for constrained multivariable control of batch processes,” _Computers & Chemical Engineering_, vol. 93, pp. 284–292, 2016.
  * [14] J. H. Lee, K. S. Lee, and W. C. Kim, “Model-based iterative learning control with a quadratic criterion for time-varying linear systems,” _Automatica_ , vol. 36, no. 5, pp. 641–657, 2000.
  * [15] U. Rosolia and F. Borrelli, “Learning model predictive control for iterative tasks. a data-driven control framework,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 7, pp. 1883–1896, 2017.
  * [16] K. S. Lee, I.-S. Chin, H. J. Lee, and J. H. Lee, “Model predictive control technique combined with iterative learning for batch processes,” _AIChE Journal_ , vol. 45, no. 10, pp. 2175–2187, 1999.
  * [17] Z. K. Nagy and R. D. Braatz, “Robust nonlinear model predictive control of batch processes,” _AIChE Journal_ , vol. 49, no. 7, pp. 1776–1786, 2003\.
  * [18] P. Bondi, G. Casalino, and L. Gambardella, “On the iterative learning control theory for robotic manipulators,” _IEEE Journal on Robotics and Automation_ , vol. 4, no. 1, pp. 14–22, 1988.
  * [19] S.-R. Oh, Z. Bien, and I. H. Suh, “An iterative learning control method with application to robot manipulators,” _IEEE Journal on Robotics and Automation_ , vol. 4, no. 5, pp. 508–514, 1988.
  * [20] Z. Bien and J.-X. Xu, “Iterative learning control: analysis, design, integration and applications,” 1998.