Spatial Iterative Learning Control for Multi-material Three-Dimensional Structures

Afkhami, Zahra; Pannier, Christopher; Aarnoudse, Leontine; Hoelzle, David; Barton, Kira

doi:10.1115/1.4046576

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.

Fig. 1

Demonstrative multi-material 3D structure fabricated by AM. n1 and n2 are two arbitrary build materials.

View large Download slide

Demonstrative multi-material 3D structure fabricated by AM. n₁ and n₂ are two arbitrary build materials.

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.

The finite set of

Z_{n}

for an odd and positive integer n with

Z_{1} ≜ {0}

is defined as

Z_{n} ≜ {\frac{1 - n}{2}, \frac{3 - n}{2}, \dots, \frac{n - 3}{2}, \frac{n - 1}{2}}

A generic scalar function p(x, y) sampled at discrete values

(x, y) \in Z_{M} \times Z_{N}

can be combined in the following matrix form,

p (x, y) = [\begin{matrix} p (\frac{1 - M}{2} \frac{1 - N}{2}) & p (\frac{1 - M}{2} \frac{3 - N}{2}) & \dots & p (\frac{1 - M}{2} \frac{N - 1}{2}) \\ p (\frac{3 - M}{2} \frac{3 - N}{2}) & p (\frac{3 - M}{2} \frac{3 - N}{2}) & \dots & p (\frac{3 - M}{2} \frac{N - 1}{2}) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ p (\frac{M - 1}{2} \frac{1 - N}{2}) & p (\frac{M - 1}{2} \frac{3 - N}{2}) & \dots & p (\frac{M - 1}{2} \frac{N - 1}{2}) \end{matrix}]

A vectorization operator ν(.) can be applied to a matrix

P \in R^{m \times n}

to convert the matrix into a column vector form

P \in R^{m n \times 1}

given as

P ≜ ν (P) = v e c (P^{T})

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

The notation

⪯

and

⪰

are element-wise inequalities such that

A ⪯ B ⟹ [A]_{i j} \leq [B]_{i j} \forall i, j

The Frobenius norm of a matrix

A \in R^{M \times N}

is defined as

‖ A ‖ = \sqrt{\sum_{i = 1}^{M} \sum_{j = 1}^{N} | a_{i j} |^{2}}

3 General Iteration Varying Systems

Models describing the sequence of material addition and heightmap evolution (illustrated in Fig. 2) have been introduced in Refs. [5,13]. In these models, the printed topography at layer 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,

\begin{aligned} g_{l} (x, y) & = g_{l - 1} (x, y) \\ + \sum_{\binom{m \in Z_{M}}{n \in Z_{N}}} h_{l - 1}^{(x, y)} (x - m, y - n) * f_{l} (m, n) \end{aligned}

(1)

where f_l is the input signal,

h_{l - 1}^{x, y} (m, n)

is the response of the build material to an impulse applied at coordinate

(x, y) \in Z_{M} \times Z_{N}

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

h_{l}^{(x, y)} \in [{\bar{h}}_{l}, {\underline{h}}_{l}]

⁠. It should be mentioned that the model described in Eq. does not capture drop coalescence effects. The system defined in Eq. can be transferred into the lifted form Eq. 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. in the following equation,

g_{l} = g_{l - 1} + H (g_{l - 1}) f_{l}

(2)

with

f_{l} ≜ ν (f_{l} (x, y)) \in R^{M N \times 1}

⁠,

g_{l} ≜ ν (g_{l} (x, y)) \in R^{M N \times 1}

⁠, and

H (g_{l - 1})

\in R^{M N \times M N}

are the spatially varying plant matrix associated with the impulse response

h_{l - 1}^{(x, y)}

⁠, which is also a member of the following matrix set:

H^{I} = {H \in R^{M N \times M N} | \underline{H} ⪯ H ⪯ \bar{H}}

H and

\bar{H}

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. .

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)‖ < β, $\forall (x, y) \in R^{M \times N}$ ⁠.

Fig. 2

View large Download slide

Schematic of AM spatial dynamics described in Eq. (1)

Assumptions A₁ and A₂ 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₃ 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.

For a general additive manufacturing process, we are interested in controlling the heightmap increment with respect to the previous layer described as

Δ g_{l} ≜ g_{l} - g_{l - 1} = H_{l - 1} f_{l}, H_{l - 1} \in H^{I}

(3)

where 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₁ and n₂, with repeated topology (see Fig. 3). It is assumed that the devices in Fig. 3 are printed on a thin pad of n₂ material to maintain consistency in the two spatial systems n₁ on n₂ and n₂ on n₁. 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₂ 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 $H_{l - 1, j}^{n_{1}, n_{2}}$ and $H_{l - 1, j}^{n_{2}, n_{1}}$ that describe spreading of the n₁ and n₂ materials on the corresponding previous layer topographies for the bi-material structure in Fig. 3.

Fig. 3

View large Download slide

Diagonal SILC: learning occurs in the horizontal direction from a previous device and vertical direction to incorporate topography. l denotes layer index and j denotes iteration (device) index.

Similar to Eq. , the spatial dynamics for the bi-material structure described in Fig. 3 can be expressed as

\begin{aligned} [\begin{matrix} Δ g_{l, j}^{n_{1}} \\ Δ g_{l, j}^{n_{2}} \end{matrix}] = [\begin{matrix} H_{l - 1, j}^{n_{1}, n_{2}} & 0 \\ 0 & H_{l - 1, j}^{n_{2}, n_{1}} \end{matrix}] [\begin{matrix} f_{l, j}^{n_{1}} \\ f_{l, j}^{n_{2}} \end{matrix}] & H_{l, j}^{n_{1}, n_{2}} \in H^{I_{n_{1}, n_{2}}} \\ H_{l, j}^{n_{2}, n_{1}} \in H^{I_{n_{2}, n_{1}}} \end{aligned}

(4)

with

Δ g_{l, j}^{i}

and

f_{l, j}^{i}

denoting heightmap increment and the control input of the build material i ∈ [n₁, n₂], 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,

H^{I_{n_{1}, n_{2}}}

and

H^{I_{n_{2}, n_{1}}}

are the interval sets corresponding to

H_{l, j}^{n_{1}, n_{2}}

and

H_{l, j}^{n_{2}, n_{1}}

⁠.

To aid in analysis and learning filter derivation, we define a nominal system plant matrix for each build material. These matrices are constructed from the impulse response of the desired topography at a given layer l,

H_{0}^{n_{1}, n_{2}} = H^{n_{1}, n_{2}} (g_{d_{l}}^{n_{2}}) H_{0}^{n_{2}, n_{1}} = H^{n_{2}, n_{1}} (g_{d_{l}}^{n_{1}})

(5)

where

g_{d_{l}}^{i}

is the desired topography of the build material i ∈ [n₁, n₂] at layer l. If

g_{d_{l}}^{i}

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₁ − A₃.

A4:
Plant spatial dynamics (e.g., $H_{l - 1, j}^{n_{1}, n_{2}}$ ⁠) is a function of the printing material (e.g., n₁), the build material of the previous layer (e.g., n₂), and previous layer topography (⁠ $g_{l - 1, j}^{n_{2}}$ ⁠).

Assumption A₄ 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.

To consider an arbitrary number of build materials, Eq. can be extended to a general multi-material model,

\begin{aligned} Δ g_{l, j} ≜ H_{l - 1, j} f_{l, j} \\ Δ g_{l, j} ≜ [\begin{matrix} Δ g_{1, j}^{n_{1}} \\ Δ g_{2, j}^{n_{2}} \\ ⋮ \\ Δ g_{l, j}^{n_{l}} \end{matrix}] f_{l, j} ≜ [\begin{matrix} f_{1, j}^{n_{1}} \\ f_{2, j}^{n_{2}} \\ ⋮ \\ f_{l, j}^{n_{l}} \end{matrix}] \end{aligned}

(6)

where n_l is the build material of layer l and

H_{l - 1, j}

is a block-diagonal plant matrix describing the spreading of the n_l material on the previous layer topography (⁠

g_{l, j}^{n_{l - 1}}

⁠) through the diagonal elements

H_{l - 1, j}^{n_{l}, n_{l - 1}}

⁠.

4 Diagonal SILC Design for Multi-Material Structures

Consider the multi-layer structure shown in Fig. 3 with two different spatial plant dynamics. The components of the SILC framework can be found below. For this example, we assume that all devices in the same layer (and hence same material) have the same desired heightmap increment,

Δ g_{d_{l}}

⁠.

\begin{aligned} Δ e_{l, j} ≜ [\begin{matrix} Δ e_{l, j}^{n_{1}} \\ Δ e_{l, j}^{n_{2}} \end{matrix}], Δ g_{d_{l}} ≜ [\begin{matrix} Δ g_{d_{l}}^{n_{1}} \\ Δ g_{d_{l}}^{n_{2}} \end{matrix}], H_{0} = [\begin{matrix} H_{0}^{n_{1}, n_{2}} & 0 \\ 0 & H_{0}^{n_{2}, n_{1}} \end{matrix}] \\ Δ g_{l, j} ≜ [\begin{matrix} Δ g_{l, j}^{n_{1}} \\ Δ g_{l, j}^{n_{2}} \end{matrix}], H_{l - 1, j} = [\begin{matrix} H_{l - 1, j}^{n_{1}, n_{2}} & 0 \\ 0 & H_{l - 1, j}^{n_{2}, n_{1}} \end{matrix}] f_{l, j} ≜ [\begin{matrix} f_{l, j}^{n_{1}} \\ f_{l, j}^{n_{2}} \end{matrix}] \end{aligned}

with

Δ e_{l, j} \in R^{2 M N \times 1}

⁠,

Δ g_{d_{l}} \in R^{2 M N \times 1}

⁠,

Δ g_{l, j} \in R^{2 M N \times 1}

⁠,

f_{l, j} \in R^{2 M N \times 1}

⁠, and

H_{l - 1, j} \in R^{2 M N \times 2 M N}

⁠.

Leveraging the compact form presented above, the error in heightmap increment can be defined as

Δ e_{l, j} ≜ Δ g_{d_{l}} - Δ g_{l, j} = Δ g_{d_{l}} - H_{l - 1, j} f_{l, j}

(7)

where Δg_l,j is replaced using the structure in Eq. .

Combining the feedforward signal and corresponding error illustrated in Fig. 4 yields a SILC update law of the form,

f_{l, j + 1} = L_{f_{l, j}} f_{l, j} + L_{e_{l, j}} Δ e_{l, j}

(8)

with

L_{f_{l, j}}

and

L_{e_{l, j}}

block diagonal matrices of the input and error filters defined as

L_{f_{l, j}} = [\begin{matrix} L_{f_{l, j}}^{n_{1}} & 0 \\ 0 & L_{f_{l, j}}^{n_{2}} \end{matrix}], L_{e_{l, j}} = [\begin{matrix} L_{e_{l, j}}^{n_{1}} & 0 \\ 0 & L_{e_{l, j}}^{n_{2}} \end{matrix}]

L_{f_{l, j}}^{i}

and

L_{e_{l, j}}^{i}

are the input and error filters corresponding build materials i ∈ [n₁, n₂]. Note that the update law in Eq. 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.

Fig. 4

View large Download slide

Diagonal SILC, each layer has its own SILC loop and each device is one iteration. l is the layer index and j is the iteration (device) index. w is the iteration shift operator, f_j+1 = wf_j, n₁ and n₂ are the build materials corresponding to layer l and l + 1. The update filters are iteration varying.

Substituting Eq. into Eq. yields the SILC algorithm that will be used throughout the paper.

f_{l, j + 1} = (L_{f_{l, j}} - L_{e_{l, j}} H_{l - 1, j}) f_{l, j} + L_{e_{l, j}} Δ g_{d_{l}}

(9)

The term diagonal learning stems from the inclusion of information from both previous iterations and previous layers within the update law. This has been illustrated in Figs. 3 and 4.

4.1 Design of Learning Filters.

There are many different methods that can be employed from temporal ILC to design the input and error filters in Eq. such as proportional type ILC, model inversion, and Q-filter design [15]. Here, we implement a norm optimal-SILC (NO-SILC) algorithm [5] to design the

L_{f_{l, j}}

and

L_{e_{l, j}}

matrices. NO-SILC offers more controllability in terms of the convergence behavior and final error. The quadratic cost function that is minimized to solve for the NO-SILC filters is defined as

\begin{aligned} J & = Δ e_{l, j + 1}^{T} [\begin{matrix} Q_{n_{1}} & 0 \\ 0 & Q_{n_{2}} \end{matrix}] Δ e_{l, j + 1} + f_{l, j + 1}^{T} [\begin{matrix} S_{n_{1}} & 0 \\ 0 & S_{n_{2}} \end{matrix}] f_{l, j + 1} \\ + (f_{l, j + 1} - f_{l, j})^{T} [\begin{matrix} R_{n_{1}} & 0 \\ 0 & R_{n_{2}} \end{matrix}] (f_{l, j + 1} - f_{l, j}) \end{aligned}

(10)

where

Q_{n_{1}}

and

Q_{n_{2}} \in R^{M N \times M N}

⁠,

S_{n_{1}}

and

S_{n_{2}} \in R^{M N \times M N}

⁠,

R_{n_{1}}

and

R_{n_{2}} \in R^{M N \times M N}

are symmetric positive definite matrices that penalize the layer error, input signal, and change in the input signal from iteration to iteration, respectively. We will consider the same weighting matrices for both systems such that

Q_{n_{1}} = Q_{n_{2}} = Q

⁠,

R_{n_{1}} = R_{n_{2}} = R

⁠, and

S_{n_{1}} = S_{n_{2}} = S

⁠. The weighting matrices are defined here as identity matrices multiplied by positive scalars q, s, and r such that Q = qI, S = sI, and R = rI.

Given the assumption of iteration invariant desired reference trajectories, we can relate two successive errors (iteration to iteration) within the same layer using the following equation

Δ e_{l, j + 1} = Δ e_{l, j} + H_{l - 1, j} f_{l, j} - H_{l - 1, j + 1} f_{l, j + 1}

(11)

Substituting Eq. into Eq. and setting the partial derivatives of

J

with respect to

f_{l, j + 1}^{n_{1}}

and

f_{l, j + 1}^{n_{2}}

equal to zero, the following relationships can be obtained for the filters,

\begin{aligned} L_{f_{l, j}}^{n_{1}} & = ((S + R) + H_{l - 1, j + 1}^{T_{n_{1}, n_{2}}} Q H_{l - 1, j + 1}^{n_{1}, n_{2}})^{- 1} \\ \times (R + H_{l - 1, j + 1}^{T_{n_{1}, n_{2}}} Q H_{l - 1, j}^{n_{1}, n_{2}}) \end{aligned}

(12a)

L_{e_{l, j}}^{n_{1}} = ((S + R) + H_{l - 1, j + 1}^{T_{n_{1}, n_{2}}} Q H_{l - 1, j + 1}^{n_{1}, n_{2}})^{- 1} (H_{l - 1, j + 1}^{T_{n_{1}, n_{2}}} Q)

(12b)

L_{f_{l, j}}^{n_{2}} = ((S + R) + H_{l - 1, j + 1}^{T_{n_{2}, n_{1}}} Q H_{l - 1, j + 1}^{n_{2}, n_{1}})^{- 1} \times (R + H_{l - 1, j + 1}^{T_{n_{2}, n_{1}}} Q H_{l - 1, j}^{n_{2}, n_{1}})

(12c)

L_{e_{l, j}}^{n_{2}} = ((S + R) + H_{l - 1, j + 1}^{T_{n_{2}, n_{1}}} Q H_{l - 1, j + 1}^{n_{2}, n_{1}})^{- 1} (H_{l - 1, j + 1}^{T_{n_{2}, n_{1}}} Q)

(12d)

The learning filters in Eqs. to are iteration varying matrices that depend on previous layer topography as well as existing material combinations.

If the update law in Eq. converges to some finite input value and the iteration varying spatial dynamics converge to the nominal model such that

\begin{aligned} lim_{\binom{l \to \infty}{j \to \infty}} f_{l, j} = f_{\infty}, lim_{\binom{l \to \infty}{j \to \infty}} H_{l, j} = H_{0}, lim_{\binom{l \to \infty}{j \to \infty}} L_{f_{l, j}} = L_{f_{0}}, lim_{\binom{l \to \infty}{j \to \infty}} L_{e_{l, j}} = L_{e_{0}} \end{aligned}

where

L_{f_{0}}

and

L_{e_{0}}

are based on the nominal plant model, the converged input would tend toward,

f_{\infty} = (I - L_{f_{0}} + L_{e_{0}} H_{0})^{- 1} L_{e_{0}} Δ g_{d_{l}}

(13)

Under these conditions, the performance of the system can be obtained by substituting Eq. into Eq. as following

Δ e_{\infty} = (I - H_{0}^{l} [(H_{0}^{l^{T}} Q)^{- 1} S + H_{0}^{l}]^{- 1}) Δ g_{d_{l}}

(14)

The convergence behavior and performance of NO-SILC can be tuned by the weighting coefficients in the filters [5,11].

Extending the learning filters in Eqs. – to the general linear iterative model described in Eq. results in

\begin{aligned} L_{f_{l, j}}^{n_{l}} & = ((S + R) + H_{l - 1, j + 1}^{T_{n_{l}, n_{l - 1}}} Q H_{l - 1, j + 1}^{n_{l}, n_{l - 1}})^{- 1} \\ \times (R + H_{l - 1, j + 1}^{T_{n_{l}, n_{l - 1}}} Q H_{l - 1, j}^{n_{l}, n_{l - 1}}) \end{aligned}

(15a)

L_{e_{l, j}}^{n_{l}} = ((S + R) + H_{l - 1, j + 1}^{T_{n_{l}, n_{l - 1}}} Q H_{l - 1, j + 1}^{n_{l}, n_{l - 1}})^{- 1} (H_{l - 1, j + 1}^{T_{n_{l}, n_{l - 1}}} Q)

(15b)

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.

Fig. 5

View large Download slide

Schematic of multi-material e-jet system. High-voltage pulse actuation can be independently applied to either nozzle to eject a drop of material with drop volume related to pulsewidth.

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).

Fig. 6

View large Download slide

Heightmap evolution system. (a) Controlled device fabrication follows the print-cure-scan cycle, with a change in build material at each new layer. Subscript k represents a layer and material combination. (b) Evolution of the solid surface topography as a heightmap signal is modeled as the integration of a convolution of input f with kernel h. h is a function of local variation in g for each pixel in g.

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

For the first (bottom) layer in the simulation, composed of material n₁, the underlaying surface (l = 0) is assumed to be a pre-layer of cured n₂ 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, $g_{d_{l}}$ 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 $h_{l - 1}^{(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 $h_{l - 1}^{(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 $h_{l - 1}^{(x, y)}$ ⁠; only the local variation in g_l about (x, y) affects $h_{l - 1}^{(x, y)}$ ⁠.

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

The nominal models $h_{0}^{n_{1}}$ and $h_{0}^{n_{2}}$ ⁠, corresponding to the printing of n₁ and n₂, are calculated using the regression models’ prediction of spreading on reference topographies g_d, denoted $h_{0}^{n_{1}} = h (g_{d}^{n_{2}})$ and $h_{0}^{n_{2}} = h (g_{d}^{n_{1}})$ ⁠, 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.

Fig. 7

View large Download slide

Nominal model. At left is a 3 × 3 pixel crop of a layer’s reference topography signal away from the edges so that it is uniform in height. At right is the model’s corresponding impulse response for printing the alternate material on the flat surface.

Fig. 8

View large Download slide

Effect of non-flat surface. At left is a 3 × 3 pixel crop of a g signal, indexed by material. The crop is taken near an edge to show the effect of significant surface variation. At right is the spatially varying impulse response local to a certain pixel (x, y) for printing the alternate material at the location corresponding to the crop at left.

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₁ printed on n₂, the deterministic term is $h^{n_{1}, n_{2}} (g_{l - 1, j}^{n_{2}})$ ⁠, and the additive uncertainty ${\hat{h}}^{n_{1}}$ is used, giving $h_{l, j}^{n_{1}, n_{2}} = h^{n_{1}, n_{2}} (g_{l - 1, j}^{n_{2}}) + {\hat{h}}^{n_{1}}$ ⁠, and for printing n₂ on n₁: $h_{l, j}^{n_{2}, n_{1}} = h^{n_{2}, n_{1}} (g_{l - 1, j}^{n_{1}}) + {\hat{h}}^{n_{2}}$ ⁠. The additive uncertainty is normally distributed as ${\hat{h}}^{i} = N (0, σ_{i}^{2}) 1 (3, 3)$ ⁠, where σ_i is the standard deviation of the elements of the nominal model $h_{0}^{i}$ ⁠, 1(3, 3) is a 3 × 3 matrix of 1’s for i ∈ [n₁, n₂].

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 ${\underline{h}}^{i} = h_{0}^{i} - σ^{i} 1 (3, 3)$ and ${\bar{h}}^{i} = h_{0}^{i} + σ^{i} 1 (3, 3)$ for i ∈ [n₁, n₂]. Similarly, the plant matrix bounds H and $\bar{H}$ are calculated from hⁱ and ${\bar{h}}^{i}$ 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₁ material, and three layers of the n₂ 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 (⁠ $L_{f_{0}}$ and $L_{e_{0}}$ ⁠), 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, $f_{1, 1}^{n_{1}} = 0$ and $f_{2, 1}^{n_{2}} = 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₁ and 26 for the n₂ 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.

Fig. 9

Average height increments for both materials approach a nominal value over the iterations (devices). n1 and n2 are the build materials in the odd and even layers, respectively (r = 0.2, s = 0.001). Iteration j = 0 refers to the first device.

View large Download slide

Average height increments for both materials approach a nominal value over the iterations (devices). n₁ and n₂ are the build materials in the odd and even layers, respectively (r = 0.2, s = 0.001). Iteration j = 0 refers to the first device.

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.

Fig. 10

View large Download slide

The normalized surface roughness converges over the iterations (devices) for all layers (r = 0.2, s = 0.001). Iteration j = 0 refers to the first device.

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₁ 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.

Fig. 11

SILC convergence for different values of the NO tuning parameters for the first layer with n1 material. Iteration j = 0 refers to the first device.

View large Download slide

SILC convergence for different values of the NO tuning parameters for the first layer with n₁ material. Iteration j = 0 refers to the first device.

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.

References

1.

Wang

,

Y.

,

Gao

,

F.

, and

Doyle

,

F. J.

,

2009

, “

Survey on Iterative Learning Control, Repetitive Control, and Run-to-Run Control

,”

J. Process Control

,

19

(

10

), pp.

1589

–

1600

. 10.1016/j.jprocont.2009.09.006

Google Scholar

Crossref

2.

Norrlöf

,

M.

, and

Gunnarsson

,

S.

,

2002

, “

Time and Frequency Domain Convergence Properties in Iterative Learning Control

,”

Int. J. Control

,

75

(

14

), pp.

1114

–

1126

. 10.1080/00207170210159122

Google Scholar

Crossref

3.

Meng

,

D.

, and

Moore

,

K. L.

,

2017

, “

Robust Iterative Learning Control for Nonrepetitive Uncertain Systems

,”

IEEE Trans. Automat. Contr.

,

62

(

2

), pp.

907

–

913

. 10.1109/TAC.2016.2560961

Google Scholar

Crossref

4.

Ahn

,

H.-S.

,

Moore

,

K. L.

, and

Chen

,

Y.

,

2007

,

Iterative Learning Control: Robustness and Monotonic Convergence for Interval Systems

, Vol.

1

,

Springer-Verlag, London

.

5.

Hoelzle

,

D. J.

, and

Barton

,

K. L.

,

2016

, “

On Spatial Iterative Learning Control Via 2-D Convolution: Stability Analysis and Computational Efficiency

,”

IEEE Trans. Control Syst. Technol.

,

24

(

4

), pp.

1504

–

1512

. 10.1109/TCST.2015.2501344

Google Scholar

Crossref

6.

Gibson

,

I.

,

Rosen

,

D.

, and

Stucker

,

B.

,

2015

,

Additive Manufacturing Technologies: 3D Printing, Rapid Prototyping, and Direct Digital Manufacturing

, 2nd ed.,

Springer-Verlag, New York

.

7.

Wang

,

K.

, and

Stark

,

J. P.

,

2010

, “

Direct Fabrication of Electrically Functional Microstructures by Fully Voltage-Controlled Electrohydrodynamic Jet Printing of Silver Nano-Ink

,”

Appl. Phys. A Mater. Sci. Process.

,

99

(

4

), pp.

763

–

766

. 10.1007/s00339-010-5701-5

Google Scholar

Crossref

8.

Park

,

J.-U.

,

Lee

,

J. H.

,

Paik

,

U.

,

Lu

,

Y.

, and

Rogers

,

J. A.

,

2008

, “

Nanoscale Patterns of Oligonucleotides Formed by Electrohydrodynamic Jet Printing With Applications in Biosensing and Nanomaterials Assembly

,”

Nano Lett.

,

8

(

12

), pp.

4210

–

4216

9.

Theiler

,

P. M.

,

Lütolf

,

F.

, and

Ferrini

,

R.

,

2018

, “

Non-Contact Printing of Optical Waveguides Using Capillary Bridges

,”

Opt. Express

,

26

(

9

), pp.

11934

–

11939

. 10.1364/OE.26.011934

Google Scholar

Crossref

PubMed

10.

Huang

,

Y.

,

Leu

,

M. C.

,

Mazumder

,

J.

, and

Donmez

,

A.

,

2014

, “

Additive Manufacturing: Current State, Future Potential, Gaps and Needs, and Recommendations

,”

ASME J. Manuf. Sci. Eng.

,

137

(

1

), p.

014001

.

Google Scholar

Crossref

11.

Altin

,

B.

,

Wang

,

Z.

,

Hoelzle

,

D. J.

, and

Barton

,

K.

,

2018

, “

Robust Monotonically Convergent Spatial Iterative Learning Control: Interval Systems Analysis Via Discrete Fourier Transform

,”

IEEE Trans. Control.

,

27

(

6

), pp.

2470

–

2483

. 10.1109/TCST.2018.2868039

Google Scholar

Crossref

12.

Pannier

,

C. P.

,

Ojeda

,

L.

,

Wang

,

Z.

,

Hoelzle

,

D.

, and

Barton

,

K.

,

2018

, “

An Electrohydrodynamic Jet Printer With Integrated Metrology

,”

Mechatronics

,

56

, pp.

268

–

276

. 10.1016/j.mechatronics.2018.01.001

Google Scholar

Crossref

13.

Pannier

,

C.

,

Wu

,

M.

,

Hoelzle

,

D.

, and

Barton

,

K.

,

2019

, “

LPV Models for Jet-Printed Heightmap Control

,”

2019 American Control Conference

,

Philadelphia, PA

,

July

, pp.

5402

–

5407

.

14.

Wang

,

Z.

,

Sammons

,

P. M.

,

Pannier

,

C. P.

,

Barton

,

K.

, and

Hoelzle

,

D. J.

,

2018

, “

System Identification of a Discrete Repetitive Process Model for Electrohydrodynamic Jet Printing

,”

2018 Annual American Control Conference (ACC)

,

Milwaukee, WI

,

June

.

15.

Bristow

,

D. A.

,

Tharayil

,

M.

, and

Alleyne

,

A. G.

,

2006

, “

A Survey of Iterative Learning Control

,”

IEEE Control Syst. Mag.

,

26

(

3

), pp.

96

–

114

. 10.1109/MCS.2006.1636313

Google Scholar

Crossref

Spatial Iterative Learning Control for Multi-material Three-Dimensional Structures

Abstract

1 Introduction

2 Preliminaries

3 General Iteration Varying Systems

3.1 Multi-plant System Dynamics.

4 Diagonal SILC Design for Multi-Material Structures

4.1 Design of Learning Filters.

5 Simulation Setup

6 Simulation Results

7 Conclusion

Acknowledgment

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Spatial Iterative Learning Control for Multi-material Three-Dimensional Structures

Abstract

1 Introduction

2 Preliminaries

3 General Iteration Varying Systems

3.1 Multi-plant System Dynamics.

4 Diagonal SILC Design for Multi-Material Structures

4.1 Design of Learning Filters.

5 Simulation Setup

6 Simulation Results

7 Conclusion

Acknowledgment

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only