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.
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 . 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 . 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) . Through AM, a 3D structure can be generated by the sequential addition of materials on the surface . 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 , biological sensors , and optical filters . 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 , 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 . The proposed iteration varying system in Ref.  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  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.
In this section, preliminary notations and definitions that will be used in the upcoming sections will be defined.
3 General Iteration Varying Systems
The spreading behavior on a flat surface is different from that on a non-flat surface .
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].
The plant matrix is considered bounded input, bounded output stable, meaning that there exist positive finite scalars α and β such that given a bounded input, ‖fl(x, y)‖ < α, the resulted output will always be bounded, ‖gl(x, y)‖ < β, .
Assumptions A1 and A2 denote that the spatial dynamics of a given plant (Hl−1) are a function of previous layer topography (gl−1) and the surrounding environment. Assumption A3 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 α . Given the interval plant (positive and bounded), Hl−1, and bounded input, the output will be always bounded.
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, n1 and n2, with repeated topology (see Fig. 3). It is assumed that the devices in Fig. 3 are printed on a thin pad of n2 material to maintain consistency in the two spatial systems n1 on n2 and n2 on n1. 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 A2 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 and that describe spreading of the n1 and n2 materials on the corresponding previous layer topographies for the bi-material structure in Fig. 3.
In multi-material and multi-layer structures, the following assumption holds for the plant matrix along with the assumptions described in A1 − A3.
Plant spatial dynamics (e.g., ) is a function of the printing material (e.g., n1), the build material of the previous layer (e.g., n2), and previous layer topography ().
Assumption A4 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.
4 Diagonal SILC Design for Multi-Material Structures
4.1 Design of Learning Filters.
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 . 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 n1 and 120 nm for n2. Each layer is printed in a single printing pass, so the build material alternates between n1 and n2 at each new layer.
For the first (bottom) layer in the simulation, composed of material n1, the underlaying surface (l = 0) is assumed to be a pre-layer of cured n2 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, 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 gl dependence on the convolution kernel is modeled using the multivariate regression model (method M2) of Ref. . 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 for the 3 × 3 pixel crop of the heightmap gl centered at the pixel coordinates (x, y). Note that the magnitude of gl does not affect the impulse response ; only the local variation in gl about (x, y) affects .
Since there are two material combinations (n1 printed on n2 and n2 printed on n1), two regression models must be specified. For n1 printed on n2, an equilibrium contact angle of 10 deg is used in the numerical simulations of drop spreading, whereas 5 deg is used for n2 printed on n1.
The nominal models and , corresponding to the printing of n1 and n2, are calculated using the regression models’ prediction of spreading on reference topographies gd, denoted and , 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. , we take a 2D array of cube roots of drop volumes to be the input denoted fl,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 n1 printed on n2, the deterministic term is , and the additive uncertainty is used, giving , and for printing n2 on n1: . The additive uncertainty is normally distributed as , where σi is the standard deviation of the elements of the nominal model , 1(3, 3) is a 3 × 3 matrix of 1’s for i ∈ [n1, n2].
From the impulse responses h, the plant matrices Hl,j in Eq. (4) are calculated from the lifted domain conversion defined in Sec. 3. The impulse response bounds are chosen such that and for i ∈ [n1, n2]. Similarly, the plant matrix bounds H and are calculated from hi and by the BCCB construction method described in Ref. .
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 n1 material, and three layers of the n2 material is considered as shown in Fig. 1. For analysis and ease in the implementation, we consider the learning filters in Eqs. (12a)–(12d) based on the nominal dynamic model ( and ), while the update law Eq. (9) will still use the iteration varying plant dynamics Hl−1,j.
It is important to note that the inputs of the first iteration at layers 1 and 2 are zero, and , 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 fl,1 = fl−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 n1 and 26 for the n2 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 n1 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.
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.
The authors would like to gratefully thank the support from National Science Foundation Grant Nos. GRFP-1256260 and CMMI-1727894.