## Abstract

Iterative learning control (ILC) is a powerful technique to regulate repetitive systems. Additive manufacturing falls into this category by nature of its repetitive action in building three-dimensional structures in a layer-by-layer manner. In literature, spatial ILC (SILC) has been used in conjunction with additive processes to regulate single-layer structures with only one class of material. However, SILC has the unexplored potential to regulate additive manufacturing structures with multiple build materials in a three-dimensional fashion. Estimating the appropriate feedforward signal in these structures can be challenging due to iteration varying initial conditions, system parameters, and surface interaction dynamics in different layers of multi-material structures. In this paper, SILC is used as a recursive control strategy to iteratively construct the feedforward signal to improve part quality of 3D structures that consist of at least two materials in a layer-by-layer manner. The system dynamics are approximated by discrete 2D spatial convolution using kernels that incorporate in-layer and layer-to-layer variations. We leverage the existing SILC models in literature and extend them to account for the iteration varying uncertainties in the plant model to capture a more reliable representation of the multi-material additive process. The feasibility of the proposed diagonal framework was demonstrated using simulation results of an electrohydrodynamic jet printing (e-jet) printing process.

## 1 Introduction

Iterative learning control (ILC) is a powerful technique that has been widely used in systems with repetitive characteristics, even those that lack real-time feedback signals, in order to achieve (near) perfect output tracking of a reference trajectory over a short number of iterations [1]. Temporal ILC uses past information in the time domain in order to build an appropriate feedforward control signal with the aim of ensuring convergence of the tracking error from iteration to iteration. Previous studies [2–4] have considered ILC architectures that address bounded iteration varying model parameters and provide convergence guarantees to a bounded neighborhood of a nominal system over finite iterations [3,4]. However, these methods have primarily been considered for temporal rather than spatial dynamics.

Recent work in the literature has focused on extending temporal ILC to the spatial domain such that the system parameters are all defined based on spatial coordinates [5]. In these systems, the ILC algorithm aims to decrease the spatial tracking error *e*(*x*, *y*) from iteration to iteration. Spatial ILC (SILC) has been demonstrated for topography control in additive manufacturing (AM) [5]. Through AM, a 3D structure can be generated by the sequential addition of materials on the surface [6]. Due to the wide range of materials that can be used in AM processes, AM-fabricated devices have been made in diverse applications including flexible electronics [7], biological sensors [8], and optical filters [9]. These devices can be fabricated layer-by-layer, depositing uniform layers of multiple materials as shown in Fig. 1.

Importantly, the performance of such devices depends on the uniformity and consistency of the layers; as such, the fabrication process must be able to provide strict adherence to the desired design requirements through robust control of the process. However, the lack of real-time monitoring devices that can capture in situ measurements has been a challenge for most AM systems, and in particular (*μ*-AM) systems [10], making online measurement of the output signal difficult in real-time. The output and subsequent error measurement are only available after the material is placed on the substrate using distributed sensing devices. Due to the challenges in real-time monitoring and control, SILC provides an appealing option for recursive control during the additive process.

In recent years, the SILC framework has been extended to consider iteration varying dynamics as well as model uncertainties [11]. The proposed iteration varying system in Ref. [11] enables robust monotonic convergence for a 2D, run-to-run e-jet printing process based on Lyapunov stability criteria. Despite the progress in the area of spatial dynamics, current SILC algorithms do not consider multiple sets of spatial dynamics due to multi-materials nor surface variations in multi-layer structures. Multiple layers lead to initial condition variability due to previous layer dynamics. Furthermore, in multi-material structures, the change in spatial dynamics dues to the different materials requires a MIMO SILC approach.

This work presents a new spatial learning control framework for multi-material 3D constructs that incorporate model uncertainties and spatially varying dynamics for multiple plant models. We propose a MIMO configuration that incorporates **vertical learning** through consideration of previous layer spatial dynamics [12] and **horizontal learning** from part to part to derive a **diagonal** SILC framework. Simulation results using a model of an electrohydrodynamic jet (e-jet) printing process is used to demonstrate the feasibility of the proposed diagonal SILC framework.

## 2 Preliminaries

In this section, preliminary notations and definitions that will be used in the upcoming sections will be defined.

*n*with $Z1\u225c{0}$ is defined as

*p*(

*x*,

*y*) sampled at discrete values $(x,y)\u2208ZM\xd7ZN$ can be combined in the following matrix form,

*vec*(.) is the conventional column-wise vectorization operator.

## 3 General Iteration Varying Systems

*l*,

*g*

_{l}(

*x*,

*y*) is the sum of previous layer topography

*g*

_{l−1}(

*x*,

*y*) and newly added material Δ

*g*

_{l}. It is assumed that the added material will spread on the previous layer to a finite extent [14]. The simplified system model can be described using the following 2D convolution format,

*f*

_{l}is the input signal, $hl\u22121x,y(m,n)$ is the response of the build material to an impulse applied at coordinate $(x,y)\u2208ZM\xd7ZN$ and layer

*l*− 1. It is assumed that

*h*

_{l}is an uncertain, spatially varying interval parameter, bounded by invariant upper and lower bounds such that $hl(x,y)\u2208[h\xafl,h_l]$. It should be mentioned that the model described in Eq. (1) does not capture drop coalescence effects. The system defined in Eq. (1) can be transferred into the lifted form Eq. (2) through the use of a vectorization operator defined in Sec. 2. A full description of the lifted domain conversion can be found in Ref. [12]. For brevity, we present the lifted form of Eq. (1) in the following equation,

**H**and $H\xaf$ are invariant lower and upper bounds associated with

**H**(

**g**

_{l−1}), respectively. We will use

**H**

_{l−1}instead of

**H**(

**g**

_{l−1}) for brevity. The following assumptions are considered for the system spatial dynamics described in Eq. (2).

- A1:
The spreading behavior on a flat surface is different from that on a non-flat surface [12].

- A2:
The plant spatial dynamics are causal in the temporal and noncausal in the spatial domain, meaning that the applied input at a given position will affect the output in the advanced layers and surrounding coordinates [5,12].

- A3:
The plant matrix is considered bounded input, bounded output stable, meaning that there exist positive finite scalars

*α*and*β*such that given a bounded input, ‖**f**_{l}(*x*,*y*)‖ <*α*, the resulted output will always be bounded, ‖**g**_{l}(*x*,*y*)‖ <*β*, $\u2200(x,y)\u2208RM\xd7N$.

Assumptions **A**_{1} and **A**_{2} denote that the spatial dynamics of a given plant (**H**_{l−1}) are a function of previous layer topography (**g**_{l−1}) and the surrounding environment. Assumption **A**_{3} holds for the additive system described in Eq. (2), given that material addition to the substrate is bounded by a pre-defined volume of available material. Furthermore, because of actuator constraints, the input in Eq. (2) is limited by an upper bound *α* [11]. Given the interval plant (positive and bounded), **H**_{l−1}, and bounded input, the output will be always bounded.

**H**

^{I}is the interval set associated with

**H**

_{l−1}.

### 3.1 Multi-plant System Dynamics.

Equation (3) can be extended to combine multi-plant dynamics into a single MIMO architecture. As an example, we will consider the fabrication of a two-material construct, *n*_{1} and *n*_{2}, with repeated topology (see Fig. 3). It is assumed that the devices in Fig. 3 are printed on a thin pad of *n*_{2} material to maintain consistency in the two spatial systems *n*_{1} on *n*_{2} and *n*_{2} on *n*_{1}. Although we assume that the spatial dynamics for a given material are invariant to the layer index, we must still apply the noncausality assumption of **A**_{2} to guarantee that the spatial dynamics of a given layer are a function of previous layer topography and build material. As a result, we consider two plant matrix $Hl\u22121,jn1,n2$ and $Hl\u22121,jn2,n1$ that describe spreading of the *n*_{1} and *n*_{2} materials on the corresponding previous layer topographies for the bi-material structure in Fig. 3.

*i*∈ [

*n*

_{1},

*n*

_{2}], in an layer

*l*and iteration (device)

*j*. It is assumed that the total number of layers and iterations are limited such that

*l*= 1, 2, …,

*L*and

*j*= 1, 2, …,

*J*. In addition, $HIn1,n2$ and $HIn2,n1$ are the interval sets corresponding to $Hl,jn1,n2$ and $Hl,jn2,n1$.

*l*,

*i*∈ [

*n*

_{1},

*n*

_{2}] at layer

*l*. If $gdli$ is a flat topography, then the nominal matrices are BCCB [5]. The BCCB property of a matrix makes fast Fourier transforms (FFT) possible which has been shown to be computationally less expensive in calculating matrix products and norms [5].

In multi-material and multi-layer structures, the following assumption holds for the plant matrix along with the assumptions described in **A**_{1} − **A**_{3}.

- A4:
Plant spatial dynamics (e.g., $Hl\u22121,jn1,n2$) is a function of the printing material (e.g.,

*n*_{1}), the build material of the previous layer (e.g.,*n*_{2}), and previous layer topography ($gl\u22121,jn2$).

Assumption **A**_{4} ensures that the surface energy of different materials is captured through the spatial dynamics that describe the spreading of a sessile drop as a function of the interactions induced by a particular material combination and orientation.

*n*

_{l}is the build material of layer

*l*and $Hl\u22121,j$ is a block-diagonal plant matrix describing the spreading of the

*n*

_{l}material on the previous layer topography ($gl,jnl\u22121$) through the diagonal elements $Hl\u22121,jnl,nl\u22121$.

## 4 Diagonal SILC Design for Multi-Material Structures

**g**

_{l,j}is replaced using the structure in Eq. (3).

*i*∈ [

*n*

_{1},

*n*

_{2}]. Note that the update law in Eq. (8) is based on the error in each layer, Δ

**e**

_{l,j}, not the total error. Layer error was selected as the desired goal to ensure layer repeatability and consistency, which is often a design objective in additively manufacturing constructs.

### 4.1 Design of Learning Filters.

*q*,

*s*, and

*r*such that

**Q**=

*q*

**I**,

**S**=

*s*

**I**, and

**R**=

*r*

**I**.

*a*)

*b*)

*c*)

*d*)

*a*) to (12

*d*) are iteration varying matrices that depend on previous layer topography as well as existing material combinations.

*a*)–(12

*d*) to the general linear iterative model described in Eq. (6) results in

*a*)

*b*)

## 5 Simulation Setup

In this section, we investigate the proposed diagonal SILC framework through a simulation study using a model of an electrohydrodynamic jet (e-jet) printing process. E-jet achieves material deposition using an electrostatic field, allowing for superior resolution and build material diversity. Drop-on-demand e-jet is achieved using synchronized substrate motion and high-voltage pulses applied to either nozzle of a bi-material e-jet printer, with a schematic shown in Fig. 5. Varying the rectangular wave pulsewidth allows for variation in printed drop size.

This simulation assumes a known relationship between pulsewidth and drop size so that the control input can be taken as drop size rather than pulsewidth. Controlled topography evolution requires topography updates, so printing is performed in a cycle of (1) printing an array of drops of varying sizes at discretized coordinates, (2) curing (solidifying) the drops, and then (3) scanning the solid surface to obtain a heightmap measurement of topography [12]. This print-cure-scan cycle is depicted in Fig. 6(a).

For this simulation, the device structure has the topology of Fig. 1 with layer heights of 100 nm for *n*_{1} and 120 nm for *n*_{2}. Each layer is printed in a single printing pass, so the build material alternates between *n*_{1} and *n*_{2} at each new layer.

For the first (bottom) layer in the simulation, composed of material *n*_{1}, the underlaying surface (*l* = 0) is assumed to be a pre-layer of cured *n*_{2} so that first-layer surface interactions with other substrate materials need not be considered.

Output heightmaps are measured relative to the top of the flat pre-layer, and the domain is 256 × 256 pixels. The pitch size is considered 1 *μ*m for both materials. The desired output, $gdl$ at layer *l* is uniform except for the four outer rings of pixels, which are reduced by half to better represent material dropoff at edges.

Heightmap evolution from layer to layer is simulated according to Eq. (1), as depicted in Fig. 6. The *g*_{l} dependence on the convolution kernel $hl\u22121(x,y)$ is modeled using the multivariate regression model (method M2) of Ref. [13]. In this model, numerical simulations of drop spreading to equilibrium on non-flat surfaces are pre-computed for a given equilibrium contact angle of the build material. Subsequently, an ordinary least squares multivariate linear regression is performed, taking the variation in the nine elements of each 3 × 3 pixel crop of the heightmap *g* as the predictor variables and nine elements of each measured 3 × 3 pixel impulse response *h* as the response variables. The fitted regression model is used to evaluate the spatially varying impulse response $hl\u22121(x,y)$ for the 3 × 3 pixel crop of the heightmap *g*_{l} centered at the pixel coordinates (*x*, *y*). Note that the magnitude of *g*_{l} does not affect the impulse response $hl\u22121(x,y)$; only the local variation in *g*_{l} about (*x*, *y*) affects $hl\u22121(x,y)$.

Since there are two material combinations (*n*_{1} printed on *n*_{2} and *n*_{2} printed on *n*_{1}), two regression models must be specified. For *n*_{1} printed on *n*_{2}, an equilibrium contact angle of 10 deg is used in the numerical simulations of drop spreading, whereas 5 deg is used for *n*_{2} printed on *n*_{1}.

The nominal models $h0n1$ and $h0n2$, corresponding to the printing of *n*_{1} and *n*_{2}, are calculated using the regression models’ prediction of spreading on reference topographies *g*_{d}, denoted $h0n1=h(gdn2)$ and $h0n2=h(gdn1)$, and shown in Fig. 7. In contrast, examples of impulse responses predicted on non-flat surfaces are shown in Fig. 8. The nominal plants in Eq. (5) can be calculated based on the BCCB matrix construction. Following the model of Ref. [13], we take a 2D array of cube roots of drop volumes to be the input denoted *f*_{l,j}.

Both the input *f* and output *g* are heightmaps, but only *g* represents a scanned topography. Choosing the same height units for both *f* and *g* makes *h*, and thus the system, dimensionless.

In addition to the deterministic spatial variation in the system described above, stochastic noise is added to the simulation.

For *n*_{1} printed on *n*_{2}, the deterministic term is $hn1,n2(gl\u22121,jn2)$, and the additive uncertainty $h^n1$ is used, giving $hl,jn1,n2=hn1,n2(gl\u22121,jn2)+h^n1$, and for printing *n*_{2} on *n*_{1}: $hl,jn2,n1=hn2,n1(gl\u22121,jn1)+h^n2$. The additive uncertainty is normally distributed as $h^i=N(0,\sigma i2)1(3,3)$, where *σ*_{i} is the standard deviation of the elements of the nominal model $h0i$, **1**(3, 3) is a 3 × 3 matrix of 1’s for *i* ∈ [*n*_{1}, *n*_{2}].

From the impulse responses *h*, the plant matrices **H**_{l,j} in Eq. (4) are calculated from the lifted domain conversion defined in Sec. 3. The impulse response bounds are chosen such that $h_i=h0i\u2212\sigma i1(3,3)$ and $h\xafi=h0i+\sigma i1(3,3)$ for *i* ∈ [*n*_{1}, *n*_{2}]. Similarly, the plant matrix bounds **H** and $H\xaf$ are calculated from h^{i} and $h\xafi$ by the BCCB construction method described in Ref. [5].

## 6 Simulation Results

In this section, the simulation results of the system described in Sec. 3 with repeated topology such as Fig. 1 using diagonal NO-SILC are investigated. A multi-layer structure with four layers of the *n*_{1} material, and three layers of the *n*_{2} material is considered as shown in Fig. 1. For analysis and ease in the implementation, we consider the learning filters in Eqs. (12*a*)–(12*d*) based on the nominal dynamic model ($Lf0$ and $Le0$), while the update law Eq. (9) will still use the iteration varying plant dynamics **H**_{l−1,j}.

It is important to note that the inputs of the first iteration at layers 1 and 2 are zero, $f1,1n1=0$ and $f2,1n2=0$, implying that there is no prior knowledge on the appropriate inputs for these materials. This results in no material deposition during the first iteration of layers 1 and 2. However, the input of the first device (*j* = 1 at other layers, *l* > 2, comes from the last device in the previous layer with the same material, such that **f**_{l,1} = **f**_{l−2,J}, where *J* is the total number of iterations per layer as shown in Fig. 3. A normally distributed white noise signal has been added to the input signal, with the mean and variance of 0.00 and 0.01 *μ*m, to better represent the experimental environment. The signal-to-noise ratio is approximately 12 for the *n*_{1} and 26 for the *n*_{2} material for the mean value of the converged input. Furthermore, the input will be constrained to positive definite values to ensure an additive process. The set of simulations was run for 30 iterations (devices) using matlab for different values of *s* and *r*. In all simulations, the parameter *q* was set equal to one (*q* = 1).

The average height increment and the corresponding desired height increment for each material class are presented in Fig. 9. The outputs from the layers of a given material all converged to roughly the same offset from the desired height increment where better convergence is observed in higher layers. In additive manufacturing, the design goal often focuses on achieving consistent layers with a repeatable thickness distribution, which is a highly desirable characteristic in most sensory applications where uniformity and periodicity of layers are of great importance.

The normalized surface roughness, *Ra*, at layer *l* is calculated from the standard deviation of the total heightmap divided by the average of the desired total heightmap at layer *l*. The surface roughness, which is the representation of layer flatness, converges to its final values as shown in Fig. 10. The upper layers of a single material class exhibit lower surface roughness and are more uniform across the layer. This uniformity results in smaller deviations from the nominal plant convergence to a flat layer. Small deviations from the nominal model combined with limited noise in the system results in a SILC controller that converges to a finite state. For example, in Fig. 9, the average height increment has reached to its final value after three iterations.

To compare the effects of the tuning parameters, the performance of the SILC system in Frobenius norm is demonstrated in Fig. 11 for different values of the penalty terms *r* and *s* in layer one with *n*_{1} material as a representation of all other layers. The weighting coefficients enable NO-SILC to control the rate of convergence, the final converged error, and the converged output (Figs. 9 and 11). Increasing the value of the input penalty *s* would increase the final error. On the other hand, increasing the value of *r* would reduce the convergence speed. Due to the nature of the additive process, it is also very important to use appropriate weighting coefficients to avoid accumulation and propagation of errors in the successive layers.

## 7 Conclusion

In this paper, we present a novel spatial ILC framework for non-repetitive systems that have multiple spatial dynamics with application to microscale additive manufacturing of multi-material 3D structures. To address the combined challenges of multiple dynamics due to multiple build materials and varying initial conditions due to roughness of the previous layer surface, a new diagonal SILC algorithm is proposed. Under the assumption of bounded iteration varying uncertainty, diagonal SILC learns from layer to layer and from device to device. To demonstrate the performance of the proposed framework, a bi-material structure is considered in simulation where the corresponding MIMO configuration involves two subsystems with distinct SILC control loops that are internally connected due to the layer-wise nature of the AM process. Simulation results show that diagonal SILC can be successfully employed to regulate the input of the iterative system to improve the heightmap reference tracking and the corresponding surface roughness. Future work will focus on extending the proposed diagonal SILC strategy in this paper to estimate the design boundaries for the convergence of the tracking error for iteration varying systems with multi-plant dynamics. In addition, future work is focused on the experimental validation of the diagonal SILC framework.

## Acknowledgment

The authors would like to gratefully thank the support from National Science Foundation Grant Nos. GRFP-1256260 and CMMI-1727894.