Multiphysical modeling and optimal control of material properties for photopolymerization processes

Photopolymerization-based Additive Manufacturing (AM), a technique in which a product is built in a layerwise fashion by local curing of a liquid monomer, is increasingly being adopted by the high-tech sector. Nevertheless, industry still faces several challenges to improve the repeatability of product quality, as recognized by several authorities on AM standardization. It is commonly recognized that there is a need for an in-depth understanding, in-situ monitoring and real-time control of the curing process to work towards end-products of higher quality. This motivates the investigation on closed-loop control of the curing process and the build-up of material properties. This pioneering research contributes to the development of a control-oriented model in the form of a state-space description that describes the multiphysical photopolymerization process and connects curing kinetics, heat flow, strain and stress evolution. This work focuses on one spatial dimension and is extendable to higher dimensions. Moreover, an extension to existing control systems theory is proposed to anticipatively control the process through the quadratic tracking framework. The control strategy is based on sequential linearization of the nonlinear model obtained from multiphysical modelling. This theoreticalnumerical approach demonstrates the potential of model-based control of the material property build-up during vat photopolymerization processes such as stereolithography and serves as a proof of principle.


Introduction
Many methods and techniques for rapid prototyping have flourished over the past decade and have resulted in a substantial evolution of industrialized additive manufacturing (AM), also known as 3D printing. One of the more prominent types of AM is stereolithography, also known as vat photopolymerization [1,2]. The technique of vat photopolymerization comprises the processes where liquid photopolymers, usually contained in a vat, are selectively cured through targeted light. The field of stereolithography has progressed tremendously over the last years, but there are still issues that have to be resolved which currently prevent the vast adoption of photopolymerization for AM purposes. Among these issues are the insufficient repeatability and lack of consistency of the produced parts [3,4]. Other challenges lie in minimizing residual stresses, preventing defect formation and increasing the feasible product size [3,4]. These challenges originate from uncertainty in AM equipment, materials, and processes [4]. Examples of uncertainty include UV light source intensity variations, build environment disturbances, and resin composition perturbations as a result of aging and poor mixing [5,6]. Literature provides little information on the impact of these uncertainties on final product quality. However, cases have been reported in which a disturbance of 20 % UV power reduction led to severe quality issues [7]. Moreover, current working curve approaches involve several ex-situ processing steps, which tends to a trial and error method [8]. There is a commonly accepted need for closed-loop process control [3,4] in order to improve the robustness to uncertainty through in-situ sensor-based approaches. Hence, this work develops a model-based closedloop process controller for the curing process and the buildup of material properties.

Background
In the field of vat photopolymerization different methods exist to illuminate the resin. The most common types are either based on vector scanning or on digital mask projection. An illustration of both types is depicted in Fig. 1. The vector scanning method is based on a moving   actuator, illuminating a small portion of the resin. The mask projection systems project a two-dimensional array of beams.
The modeling of stereolithography has been described extensively in the literature, starting with Jacobs' cure depth model [1]. The physics have been described more detailed depending on the illumination source which influences the polymerized profile [9][10][11][12] and the effect of stacking layers has been investigated [13][14][15][16]. During the solidification of the resin, mechanical properties start to build up. From the gel point onward, the evolution of the elastic, viscous and even plastic properties can be described as a function of the level of conversion [1,[17][18][19][20][21][22][23][24]. Stress builds up as a result of the interplay between chemical shrinkage due to the solidification and thermal expansion due to the exothermic photopolymerization process. The multiphysical modeling of the vat photopolymerization process has been treated by means of a finite element framework [25,26]. The downside of such a finite element framework is that it is rather computationally expensive, and is therefore in its current state not suitable for process control purposes.
Even though control systems are inherently present in AM machines, the actual material and geometry transformations are controlled open-loop [27]. Few studies show feedback control on processes similar to vat photopolymerization. Among these studies, Zhao and Rosen [28,29] demonstrated closed-loop control of the cure height in a single-layer bottom-up vat photopolymerization process. Türeyen et al. [30] applied iterative learning control to improve the obtained geometry by measuring the difference between the actual geometry and reference geometry. Garanger controlled the stiffness of a 3D printed leaf spring by using different infill densities as the control variables [31,32]. Finally, Yebi [33], minimized cure level deviations throughout the layer's depth for thick resin-infused glass fiber composites.
Recently, closed-loop control of photopolymer conversion was demonstrated through proof-of-principle experiments at the sub-voxel scale with a slurry for ceramic vat photopolymerization [6]. This work can be considered as a first fundamental step towards real-time in-situ control of the polymerization reaction. However, several steps need to be taken before closed-loop control can be implemented in industrial, large-scale AM machines, especially when aiming to improve eventual product quality. Among these steps are the translation to actual material properties and moving towards product-scale.

Research objective
The main research objective of the present work is to demonstrate that the mechanical properties, which build up during photopolymerization, can be controlled in closed loop. The control of these mechanical properties can be considered at various spatial scales. To illustrate this, a roadmap to control the mechanical properties in AM is depicted in Fig. 2.
At the smallest scale, the chemical reaction is considered in an infinitesimally small 0D volume. This can be considered a single-input single-output (SISO) [34] problem where the input is the light intensity and the output the degree of conversion. At the largest scale, the entire three-dimensional AM machine is considered. This is a multi-input and spatially distributed or infinite-dimensional problem (MI∞O) [34]. The inputs are the light intensity and its spatial projection onto the resin. The degree of conversion has to be controlled at infinitely many points in the 3D build space. This research will simplify the spatial nature to a single dimension with the incident light intensity as a boundary input. The one-dimensional domain is considered as a first step to move towards full 3D MI∞O control. This work has a rather broad scope and develops a generic framework that can be used for different control objectives.

Outline
This paper is outlined as follows. The multiphysical photopolymerization modeling framework is developed in a modular approach in Section 2. The interconnection of the submodels is captured in a nonlinear control-oriented representation and a case study to cure a single layer is introduced. In Section 3, the nonlinear model is linearized and an extension to control theory on the tracking problem is formulated as an optimization problem. A state estimator is introduced to estimate variables that are not directly measured. To show the flexibility of the control strategy, the control of material properties in a single layer is demonstrated by means of simulations in Section 4. Finally, conclusions are drawn in Section 5.

Modeling of the photopolymerization process
Since photopolymerization is a multiphysical process spanning multiple research domains, a modular approach is taken, which allows for expansions at a later stage. Following this approach, the photopolymerization process is divided into four submodels as depicted in Fig. 3. The total model is defined as the interconnection of the respective submodels, denoted as Σ : In the remainder, a one-dimensional domain Ω is considered as illustrated in Fig. 4. The coordinate z ∈ [0, L] describes the depth coordinate of a material point. The incident light intensity is denoted I in , the top boundary is denoted Ω 0 and the bottom boundary is denoted Ω L .

Irradiation
The light intensity, which initiates the polymerization reaction, decreases as a function of the path length into the medium. In traditional photopolymerization, a homogeneous polymer resin is used, but (ceramic) particles can be included as well. In the absence of light scattering due to included particles, the light attenuation is absorption dominated and the Beer-Lambert law relates the optical attenuation to the attenuating species, their absorptivity and the path length in the medium [35,36].
For the one-dimensional problem with time-varying incident light intensity, the light intensity at coordinate z can be written in terms of the penetration depth D p [37,38] as where the time t ∈ R + , I in (t) ∈ R + and the penetration depth is the depth where the irradiance has reduced to 37% of its initial value [1,2]. Considering the domain Ω, the exponential irradiance attenuation is illustrated schematically in Fig. 4. Note that the incident light intensity, and thus also the actual light intensity, can not be negative. The cumulative energy provided to the resin is described by E = I in (t)dt. The input of this submodel, denoted by Σ I , is the incident light intensity I in (t) and the output the irradiance profile over depth and time I(z, t).

Reaction kinetics
The main principle of photopolymerization is the light initiated curing of a liquid monomer into solid polymer. The degree of conversion is often used to quantify the amount of monomer that is converted into polymer, defined as where [M ] is the concentration of monomer and [M ] 0 is the initial concentration of monomer [39]. The quantity p is also referred to as the degree of conversion.
To restrict the number of equations involved and retain a control-oriented framework, a phenomenological model is deemed suitable [40,41]. The typical form of such phenomenological model is where r(T ) is a temperature dependent rate constant and f (p) is a function of the degree of conversion. The constant exponent b ∈ [0.5, 1] [42]. In its simplest form, an m th order reaction kinetic model as a function of the z-coordinate and time is used to define the submodel where p max ≤ 1 denotes the asymptotic bound on degree of conversion and it is assumed that the rate constant r is temperature independent. The ultimate achievable degree of conversion, p max , is lower than one due to the reduced mobility of the reactive species upon curing [43]. With this relatively simple curing kinetics model, the dark curing reaction is ignored. The input of this submodel is the light intensity I(z, t) and the output is the degree of conversion p(z, t). The degree of conversion and the initial condition p 0 (z) are bounded by zero and p max .

Heat propagation
Polymerization is an exothermic process. This heat source and the heat transfer via conduction and convection is often modeled with the conservation law for energy [37] described by the partial differential heat equation where ρ is the density of the medium, c p its heat transfer coefficient, T the absolute temperature and k the thermal conductivity. The parameter ∆H p is referred to as the heat of polymerization. A convection boundary condition is assumed at the top of the specimen, Ω 0 , which can be described as where T ∞ is the ambient temperature. A Neumann boundary condition, resembling perfect insulation, is assumed at the bottom of the specimen Ω L and is described by A transformation to the deviation from the ambient temperature viaT (z, t) = T (z, t)−T ∞ completes the modeling of Σ T , which is defined as in Ω (8) where the input of this submodule is the degree of conversion p(z, t) and the curing rate ∂p(z,t) ∂t . The outputT (z, t) is the distribution of temperature over the medium, normalized with respect to the ambient temperature.

Stress and strain evolution
To model the evolution of mechanical properties during solidification, the specimen is envisioned to be quasi-threedimensional as depicted in Fig. 5. A slender specimen is considered where dx L, dy L and its constraints depend on the load case.
The total strain ε tot is decomposed into where ε el is the infinitesimal linear elastic strain, ε ch is the strain resulting from chemical shrinkage and ε th is the thermal strain. It is assumed that the material shrinks linearly and isotropically due to curing, described by where ε ch max is a negative scalar indicating the maximum shrinkage strain. Similarly, it is assumed that the material expands linearly and isotropically due to temperature changes. Therefore, the thermal strain is described by where α is the thermal expansion coefficient and T 0 is the initial temperature. If it is assumed that the specimen initially has temperature T (z, 0) = T 0 (z) = T ∞ , Eq. (11) can be written as ε th = αT .
Differentiating the chemical strain and the thermal strain with respect to time and writing it as an initial value problem givesε ch = ε ch maxṗ , where initially ε ch (z, 0) = ε ch 0 (z) and ε th (z, 0) = ε th 0 (z).
Linear elastic behavior is described by the generalized Hooke's law [44], described as where the stress σ := [σ xx , σ yy , σ zz , τ xy , τ xz , τ yz ] and the elastic strain ε el := ε el xx , ε el yy , ε el zz , γ el xy , γ el xz , γ el yz . The stiffness matrix is defined as where Lamé's constants λ := Eν (1+ν)(1−2ν) and µ := E 2(1+ν) . The solidification of the resin is captured in a conversiondependent Young's modulus E(p) [19,25]. It is assumed that the Poission's ratio, denoted by ν, is constant throughout polymerization. A relatively simple form to express such conversion-based elasticity is the piecewise linear function where E 0 1 and indicates that the mechanical properties are insignificant in the monomeric phase. The Young's modulus for the fully polymerized material is denoted by E pol . To omit the discontinuity of this function, it is approximated by an n th order polynomial such that E(p) = η(p)E pol , with the constraints η(0) = E 0 and dη(0) dp = 0 as The constants c 2 to c n are obtained via a least-squares curve fit. Depending on the application, different boundary conditions constrain the specimen. For now it is assumed that either the strain or the stress is fixed for each of the six components. The stress vector is rearranged such that the unknown and known subvectors are separated as σ =: [σ uk , σ k ] and similarly for the strains as ε =: [ε k , ε uk ] . The stiffness matrix is separated as C = η(p)C. The stress-strain relation can be written as Rearranging the unknowns to the left-hand side and writing as initial value problem gives in whichC = C 11 −C 12C −1 22C 21 and the initial conditions are denoted by ε uk (z, 0) = ε uk 0 (z) and σ uk (z, 0) = σ uk 0 (z). This completes the modeling of Σ M . For the onedimensional problem, it is defined as where f (·) denotes the corresponding function described in Eq. (20). The inputs of this module are the degree of conversion p(z, t) and temperature profileT (z, t) as well as their time derivatives. The boundary and initial conditions completely describe the mechanical model.

Spatial discretization
To simplify the model and facilitate the control of a spatially distributed system described by a partial differential equation, the model is discretized to obtain a set of ordinary differential equations. The z-coordinate z ∈ [0, L] is uniformly discretized as Here, N is the number of nodes in z−direction, including the boundary nodes; hence, the spacing between the nodes is equal to ∆z = L N −1 as shown in Fig. 6. For a stepwise derivation of the discretized model, the reader is referred to [45,46]. This discretized model is useful since it can be written as a control-oriented state-space representation of the formẋ(t) = f (x(t), u(t)), where x(t) denotes the state of the system and u(t) the input. Due to discretization, the state of the system x(t) is now of finite dimension.

Nonlinear state-space description
In control theory, systems are often described in a statespace representation which is a mathematical model of the system where the variables are related by first-order differential equations [47]. Casting the multiphysical photopolymerization into this framework and dealing with the distributed parameter nature of this system enables to apply well-established control techniques. To cast the system in this control-oriented form, the state vector is defined as The number of unknown strains in a node is denoted by n ε,uk and the number of unknown stresses by n σ,uk where the sum of both amounts to six. Then the unknown strains are defined as ε uk (t) : This means that the state x(t) ∈ X := P × R 9N . The initial state at t = 0 is defined as It is important to note that the light intensity cannot be negative which means that the input u(t) := I in (t) b ∈ I := R + . Now, the dynamics may be written in a state-space representation as where the nonlinear function f : X × I → R 10N . Considering that p(t) ∈ P, the nonlinear function g : and the constant j = ∆Hp ρcp . The matrix G relates the temperature flow between the nodes and is equal to where ζ = k ρcp(∆z) 2 and ξ = 2h ρcp(∆z) . The functions f ε uk (·) and f σ uk (·) describe the stress-strain differential equations for all nodes consecutively as and For an elaborate derivation, the reader is referred to [45].

Single, semi-infinite layer
Suppose that a semi-infinite layer of resin is equally illuminated as depicted in Fig. 7. Due to periodicity, the boundaries of the specimen are fixed in the (x, y)−plane, i.e., ε tot xx (z, t) = ε tot yy (z, t) = 0. Illumination of this specimen will not induce any shear. Hence γ tot xy (z, t) = γ tot xz (z, t) = γ tot yz (z, t) = 0. The material is free to deform in zdirection. If it is assumed that the material relaxes instantly, σ zz (z, t) = 0. Defined by the boundary conditions, the known strains comprise the elastic strain in x−direction, the elastic strain in y−direction and the shear strains. The stress is only known in the z−direction. The absence of shear enables to reduce the stress-strain relation from six-dimensional to three-dimensional. Furthermore, the problem is symmetric in the (x, y)-plane. This reduces the unknown in-plane stresses for each node, denoted σ xx,yy , to solving a scalar equation. The only unknown strain, ε el zz from which ε tot zz is derived, is also described by a scalar equation. The reduced state-space description for this particular case is derived and described in [45].
The behavior of the system, driven by a time-varying light intensity, can easily be analyzed by numerically solving the differential equation. The model can be used to solve a one-dimensional problem as illustrated for the case of a single layer, but there is sufficient potential to expand to the two-or three-dimensional case. Also, the model can be reduced to a zero-dimensional configuration by considering a single node.

Control of the photopolymerization process
Capturing the vat photopolymerization process behavior in a control-relevant model in state-space form opens up new applications in the well-established field of control techniques. In the current literature there is little documentation and little consensus about efficient strategies to achieve beneficial material properties. For this reason, the aim is to opt for a flexible control strategy in the sense that it should accommodate for different control objectives and it should be tunable based on engineering intuition. Finally, in view of the short reaction time scale in industrial applications, the real-time computational cost should be low.
The quadratic tracking framework fits these aforementioned requirements [48]. First of all, the framework is able to cope with nonlinear dynamics by means of a sequential linearization strategy. The framework can be used for a variety of objective functions and is intuitively tunable. Also, the control gains can be implemented as look-up tables, resulting in low real-time computational demands. Solving the quadratic tracking problem for nonlinear systems is, however, not trivial. Ç imen and Banks approached this problem by describing the system as an input-affine pseudolinear system of the formẋ(t) = A(x(t))x(t) + B(x(t))u(t) [49,50]. The nonlinear optimal tracking problem has been solved for Linear-Time-Varying (LTV) approximations of the true nonlinear system. Alternatively, this work proposes a recursive linearization strategy. Although convergence properties are yet to be proven, the proposed strategy is likely to converge to the optimal control input for the true nonlinear system, similar to the work of Ç imen and Banks.
Typically not every variable is measured in the vat photopolymerization process, but unmeasured variables can be estimated from measured ones even in the presence of measurement noise. Ideally, one would be able to measure the material properties at every (sub-)voxel in the product. However, the state-of-the-art measuring systems are not capable yet to do so. It is possible to measure variables at the boundary of the layer, such as the degree of conversion or temperature, and use these measurements to infer material properties throughout the product. Previous work provides a detailed discussion on measurement, actuation and control in vat photopolymerization [27]. A framework to estimate unmeasured variables from noise corrupted measured data incurred from a nonlinear process model is the extended Kalman filter [51].
The to-be-controlled nonlinear system is described by Eq. (23). In contrast to traditional actuation methods where the control variable is the exposure time, the present approach proposes to exploit the hardware's actuation capabilities by modulating the UV light intensity. As exemplified by soft-start cure strategies [52], exploiting such actuation capabilities is expected to improve the material properties with respect to simple UV on/off control strategies which merely control the exposure time [5,33]. Literature is inconclusive regarding the difference in performance between tracking complete reference trajectories throughout the curing reaction and merely steering towards a final cure level target. However, cure strategies such as soft-start suggest that different irradiance trajectories I in (t) but with the same exposure E can lead to different material properties. Hence, reference tracking provides a mechanism to control the time evolution of the curing process, rather than aiming to reach a final degree of cure alone. Note that in either case, the current machine software typically does not allow such modulation as it is closed and proprietary.
In the remainder of this section, the system is first linearized. Then the optimal control problem is formulated for the approximated system with time-varying weights and finally a state estimator is introduced in the form of the extended Kalman filter.

Linearized system
The function f (x, u), can be approximated with a Taylor expansion around a time-varying nominal trajectory (x * (t), u * (t)) on the interval t ∈ [0, T ]. The nonlinear system is approximated by a first-order Taylor expansion aṡ where and B(t) := ∂f ∂u (x * (t),u * (t)) .
Equation (28) is rearranged intȯ . (31) where d(t) is an additive term solely dependent on the nominal trajectory. Similarly, the time-varying matrices A(t) and B(t) are solely dependent on x * (t) and u * (t) and are therefore known.

Quadratic tracking problem with time-varying weights for nonlinear systems
The aim is to track a time-varying reference x ref (t) on the interval t ∈ [0, T ] with the system described by a nonlinear state-space representation of the forṁ with initial condition x 0 := x(0). This nonlinear system is approximated with Eq. (31) along a nominal trajectory (x * (t), u * (t)) on the interval t ∈ [0, T ]. Let us define the cost criterion J as and where the time-varying scalars θ Q (t) and θ R (t) define how the weight matricesQ andR vary over time for t ∈ [0, T ]. Given an initial condition x 0 ∈ X 0 := Rn, wheren is the state dimension, and an a-priori known reference signal where the input signal u ∈ U := {u : [0, T ] → Rm} and J : X 0 × X ref × U → R is defined in Eq. (33). The parameterm denotes the input dimension. For the specific state-space description considered in Section 2.6,n = 10N andm = 1.

Theorem I:
Given the initial condition x 0 ∈ X 0 and the reference trajectory x ref ∈ X ref . Then the cost Eq. (33) of the optimization problem Eq. (36) satisfies where the norm ||φ(t)|| 2 R(t) := φ(t) R(t)φ(t) and the scalar function (38) in which the symmetric matrix function P (t), the vector function r(t) and the scalar function c(t) are the unique solutions of the differential equations with boundary conditionṡ Moreover, the optimal control that solves the optimization problem Eq. (36) is given by with 0 ≤ t ≤ T and the optimal value J(x 0 , x ref , u opt ) = V (x 0 , 0).
It is interesting to point out that the optimal control u opt anticipates on the reference trajectory Theorem II: Suppose that the nominal trajectory (x * , u * ) converges to a solution (x, u) ofẋ = f (x, u) in a L 2 sense, i.e., In particular, the optimal cost of the optimization problem subjected to the approximated dynamics with x(0) = x 0 and u ∈ U converges to the optimal cost of the optimization problem subjected to the true nonlinear model with x(0) = x 0 and u ∈ U.
The proofs of Theorem I and Theorem II are provided in Appendix B and Appendix C, respectively. In case a perfect linearization is achieved, i.e., the nominal trajectory coincides with the solution, the optimal control law for the actual nonlinear system is obtained. The control scheme is visualized in Fig. 8. Note that the computation of the anticipative signal u AN at time t requires knowledge of the reference signal x ref (τ ) for all τ ∈ [0, T ]. Key to the control strategy is that the nominal trajectory is close to the true state trajectory. If these differ too much, the system approximation will be poor and performance will be suboptimal. Estimates of the optimal control law can be iteratively improved via a sequential linearization sequence. The solution of the previous iteration is recursively used as nominal trajectory for the subsequent iteration until the deviation between nominal trajectory and true state trajectory converge to zero [45,53].

State estimation
If not all states x(t) in Eq. (42) are measured, then the unmeasured states may be estimated with a state estimator. It is assumed that measurements inferred from sensors are available and corrupted with additive Gaussian white noise N (µ, σ 2 ), with mean µ and co-variance σ 2 . For nonlinear systems, the states can be estimated with an extended Kalman filter [51]. Specifically, consider the nonlinear dynamic system described bẏ where y(t) is a measured output variable. An additive process noise is denoted with w(t) and an additive measurement noise with v(t). The continuous time extended Kalman filter has the estimator dynamicṡ with gain matrix in which t ∈ [0, T ]. The time-varying matrix Φ(t) is the symmetric positive definite solution to the algebraic Riccati differential equatioṅ where t ∈ [0, T ] and The control scheme, augmented with the state estimator, is visualized in Fig. 9.
anticipative Figure 9: Quadratic tracking control scheme of Fig. 8, extended with a state estimator.

Simulation of the Controlled System
To illustrate the capabilities of the closed-loop control system, the semi-infinite thin layer, introduced in Section 2.7, c.f. Fig. 7, is considered. The state of the system is defined as x(t) := [p(t),T (t), ε ch (t), ε th (t), ε tot zz (t), σ xx,yy (t)] ∈ P × R 5N . The control gains are determined via the recursive linearization strategy and are implemented to control the nonlinear system. The main aim is to show the versatility of the model as well as its potential to be used for various purposes. First, two case studies are shown with full state feedback. This is followed by a case study with a limited set of measured quantities that consist of the degree of conversion and the measured temperature both located at the top Ω 0 . The control goal in this last study is to track a desired degree of conversion at the bottom Ω L . The physical parameters for all simulations are given in Table 1. The simulation settings and initial conditions are given in Appendix A. The differential equations are solved using the explicit Runge-Kutta (4,5) method [56,57].

4.1.
Tracking of an arbitrary degree of conversion reference In the first case study, the aim is to track an arbitrary degree of conversion reference throughout the entire thickness of the specimen. The designed reference contains multiple discontinuities such that the anticipative nature of the controller can be observed clearly. The weights are chosen such that the error throughout the layer is equally penalized. Note that this reference is not typically used in photopolymerization.
To illustrate the possibility for time-varying weights, Q(t) as defined by Eq. (34) is chosen to be linearly increasing such that the error is penalized more over time. Therefore, θ Q (t) = 100+450t is chosen to linearly increase and the weight matrix is chosen to be diagonal of the form Q = diag (I N ×N , 0 N ×N , . . . , 0 N ×N ). The terminal cost in Eq. (33) is chosen as Q T =Q. The weight on the control input R(t), as defined by Eq. (35), is chosen constant over time where θ R (t) = 10 −2 andR = 1. The system is simulated in closed-loop with the optimal control law described by Eq. (42). The control law is obtained from the sequential linearization algorithm, with initial nominal trajectories x * 0 (t) = 0 N and u * 0 (t) = 0. After five iterations, the nominal trajectory has converged sufficiently to the true state trajectory. The degree of conversion at the top of the layer, z = 0, in the middle of the layer, z = 0.5L and at the bottom of the layer, z = L, are shown in Fig. 10. The Figure 10: Tracking of an arbitrary degree of conversion reference ( ) throughout the entire specimen. The degree of conversion at z = 0 m ( ), z = 2.5 · 10 −5 m ( ) and z = 5 · 10 −5 m ( ) are depicted, as well as the optimal control input ( ) and cumulative supplied energy ( ). The reference trajectory is followed on average; however, deviations from the reference are inevitable due to the spatial variation in light intensity.
reference for each of these states is equal and also depicted.
It can be observed that the effect of light attenuation with depth (the Beer-Lambert law) is unavoidable, which means that there always is a certain deviation from the reference signal at some depth z. However, this control strategy enables to minimize this state error optimally over the entire layer in the sense of minimizing the pre-defined cost. For this specific case, the penetration depth is chosen equal to the layer thickness, i.e., D p = L. In case L Dp 1 there is a lot of deviation throughout the layer since only the top nodes are sensitive to the input. In case L Dp 1, all nodes receive approximately the same light intensity and the trajectories over the layer thickness almost coincide.
The anticipative nature of the controller can be observed clearly at t ≈ 0.2 s and at t ≈ 0.9 s. The controller initiates the polymerization reaction before the ramp and step in the reference occurs in order to minimize the cost over the entire time window.
The temperature profile of the specimen is shown in Fig. 11. The temperature throughout the specimen is al- Figure 11: Temperature evolution at different z-coordinates, corresponding with tracking an arbitrary degree of conversion reference. The temperature at z = 0 m ( ), the temperature at z = 2.5 · 10 −5 m ( ) and the temperature at z = 5 · 10 −5 m ( ) are depicted. Due to the relatively high heat conduction, the temperature deviates little throughout the specimen. most equal due to the relatively high heat conduction. After t = 1.2 s the effect of convection is clearly visible and the specimen eventually cools down to the ambient temperature. This implies that the thermal strain will eventually converge to zero.
The chemical and thermal strain follow the same trend as the degree of conversion and temperature depicted in Figures 10 and 11. The total strain in z-direction is shown in Figure 12. In the time window t ∈ [0, 2] s, the thermal expansion dominates the chemical shrinkage resulting in a positive total strain in z−direction.
The principal stresses in the (x, y)−plane are depicted in Fig. 13. The stresses start to emerge as a result of the chemical and thermal strains that will result in a mechanical strain. Since on this time window the thermal strain dominates the chemical strain, the material has the tendency to expand, resulting in negative stresses.
In a similar fashion, any other state variables may be tracked by altering the time-varying weight matrix Q(t),   i.e., tracking temperature, strain or stress profiles. Quantities can be tracked at a specific position in the specimen, but also on average as illustrated by this example.

Tracking of an arbitrary stress reference
In the following case study, an arbitrary stress profile is tracked at the bottom of the layer. The weight on the error, Q(t) as described with Eq. (34), is chosen to be constant over time, hence θ Q (t) = 10 −14 . The weight matrix is defined asQ = diag(0 N ×N , . . . , 0 N ×N , 0 (N −1)×(N −1) , 1) and the terminal weight in Eq. (33) is chosen as Q T = 10 −3Q . Hence, only the stress deviation from the reference at the bottom of the specimen is penalized. The weight on the control input R(t) as defined in Eq. (35) is chosen as θ R (t) = 1 andR = 10 −3 . The difference in magnitude of the weights is due to the high magnitude of stresses. To prevent ill-conditioning, the weights are adjusted and stresses are computed in a scaled unit, i.e., MPa in this case. The simulated stress profiles and corresponding control input are depicted in Fig. 14, illustrating that this arbitrary feasible stress profile can be tracked.

Tracking of an arbitrary conversion reference with state estimation
The aim of this case study is to show that unmeasured variables can also be controlled on the condition that these are observable from the measurements. Only two state variables are measured: the degree of conversion and temperature at the top of the specimen. Based on these measurements and the state estimator model, the other state variables are estimated. This is a relevant case since it is typically easier to measure at a boundary of the specimen rather than measuring inside the specimen. The tracking error is penalized at the bottom boundary to illustrate the control of inferred and non-measured variables. The aim is to track a sigmoid-shaped degree of conversion profile at the bottom of the specimen. Hence, only the error in degree of conversion at the bottom of the specimen is penalized. It is assumed that the nonlinear system is perturbed with an additive, normally-distributed, process noise w(t) viȧ The process noise has the characteristics w(t) ∈ N (0, W ), where the covariance matrix W = (I N ×N · 10 −2 , I N ×N , I N ×N · 10 −4 , I N ×N · 10 −4 , I N ×N · 10 −4 , I N ×N · 10 16 ). The output equation is given by :=C x(t)+v(t), (51) where the output noise v(t) ∈ N (0, V ) and the matrix V = diag(10 −2 , 1). The variances for the noise signals w(t) and v(t) are based on the magnitude of the signalṡ x(t) and y(t) without noise and are chosen such that the noise signal contributes between 1% and 10% to the total signal. For example the degree of conversion measurement, i.e., the first output in y(t), is perturbed with an additive noise signal with characteristics N (0, 10 −2 ). Hence, the mean of the perturbation µ = 0 and the standard deviation σ = 0.1. Similarly, the temperature measurement is perturbed by a signal with mean µ = 0 K and standard deviation σ = 1 K. Based on experience with FTIR spectrometry [6], a 10 % measurement uncertainty in degree of conversion is considered large. Nevertheless, the proposed strategy is able to handle larger uncertainties at the cost of a deteriorated estimate and as a result a larger overall error. Often, in practice little is known about the noise signal and the engineer uses the matrices W and V to tune how much to rely on the model with respect to the measurements. The degree of conversion is shown in Fig. 15 as well as its estimate and the control input. The corresponding stress evolution and estimate are depicted in Fig. 16.  It can be concluded that the state-estimator estimates the actual state well and that the control of the degree of conversion at the bottom of the layer is achieved with limited measurement knowledge. In case the magnitude of the noise is increased, the quality of the estimates deteriorate. The quality of the estimate will deteriorate most in the stress variable, since it is based on a number of other estimated variables.

Conclusions
This work presents the development of a modular control-oriented model describing the full multiphysical vat photopolymerization process. The main focus of this work is on a one-dimensional problem describing the spatio-temporal evolution of the degree of conversion, which is related to the build-up of mechanical properties. A nodal discretization is applied to deal with the infinite dimensionality of the problem and a quasi-three-dimensional domain is considered to describe the build-up of stresses and strains. The solidification process is described by a linear elastic material model with conversion-dependent Young's modulus. The entire control-oriented model is captured in a single nonlinear state-space description enabling to close the control loop at a new level in AM: the level of material properties.
To control the system, a tracking control strategy is proposed, based on the sequential linearization of the nonlinear model and the minimization of a quadratic cost criterion. An optimal control input is derived for the approximated dynamics and optimality of the strategy is proven for the true nonlinear system, provided that the nominal trajectory coincides with the true state trajectory. The control algorithm's implementation consists of a look-up table and therefore requires little real-time computational effort.
The feasibility to control the nonlinear dynamics that describe the evolution of material properties during photopolymerization is shown on the basis of multiple simulation-based case studies. It is shown that arbitrary state variables of the model can be tracked in an anticipative fashion. In addition, it is shown that nonmeasured states may be controlled by augmenting the control scheme with an extended Kalman filter. Evidently, the model is a simplified representation of reality. The modular modeling approach leaves room to expand the model and include more phenomena, such as dark cure or inhibition of, e.g., oxygen, in order to better represent reality. The arguably oversimplified material model can be improved by introducing more rate dependent effects, such as viscoelasticity. Even though the real-time control algorithm itself is computationally effective, computing the control gains is computationally expensive. This means that in order to make the transition to a 3D or a significantly more complex model, the scaling of the computational load needs to be addressed. In spite of these limitations the presented results are promising and can be considered as a step towards a new process control paradigm where material properties are controlled in closed-loop and in real-time at the full machine scale. This control paradigm is therefore expected to become increasingly important in improving AM part quality.

(B.1)
For the remainder of this proof, the function arguments are left out to improve readability. The optimization problem is subjected to the systemẋ = Ax + Bu + d which after transposing may be written asẋ = x A + u B + d . Substitution of the latter and the finite horizon Riccati equation, described by Eq. (39a), into Eq. (B.1) giveṡ V = x P BR −1 B P − Q x + u B P x + x P Bu + 2x A r + 2u B r + 2x ṙ +ċ + d P x + x P d + 2d r.

(B.2)
Via completion of the squares, the norm u + R −1 B P x 2 R = u B P x + x P Bu + u Ru + x P BR −1 B P x can be used to rewriteV intȯ V = −x Qx − u Ru + u + R −1 B P x 2 R + 2x A r + 2u B r + 2x ṙ +ċ + d P x + x P d + 2d r.  Note that the left-hand side is equal to J(x 0 , x ref , u), c.f. Eq. (33). Since x(0) P (0)x(0) + 2x(0) r(0) + c(0) is constant, the optimal input, u opt (t), minimizing the cost J is provided by Eq. (42).

Appendix C. Proof Theorem II
In case the nominal trajectory is equal to the actual state trajectory, the approximation of the nonlinear system is perfect. Revisiting the Taylor approximation of the nonlinear system given by Equation (28) shows that if the nominal trajectory, in the limit case, has converged to the actual state trajectory in a square integrable sense, the perturbations x * (t) − x(t) and u * (t) − u(t) vanish. This means that a perfect approximation of the nonlinear statespace form is achieved in the limit case.