Optimal Control of Material Micro-Structures

In this paper, we consider the optimal control of material micro-structures. Such material micro-structures are modeled by the so-called phase field model. We study the underlying physical structure of the model and propose a data based approach for its optimal control, along with a comparison to the control using a state of the art Reinforcement Learning (RL) algorithm. Simulation results show the feasibility of optimally controlling such micro-structures to attain desired material properties and complex target micro-structures.


I. INTRODUCTION
In this paper, we consider the optimal control of complex and high DoF (Degree of Freedom) material micro-structure dynamics, by doing which we aim at developing a robust tool for exploring the fundamental limits of materials microstructure control during processing. The micro-structure dynamics process is governed by the so-called Phase Field (PF) model [1], which is a powerful tool for modeling micro-structure evolution of materials. It represents the spatial distribution of physical quantities of interest by an order parameter field governed by one or a set of partial differential equations (PDEs). The Phase Field method can naturally handle material systems with complicated interface, thanks to which it has been successfully applied to the simulation of a wide range of materials micro-structure evolution processes, such as e.g., alloy solidification [2] [3], alloy phase transformation [4] [5], grain growth [6] [7], and failure mechanism of materials [8] [9] etc. The phase field model is complex, nonlinear, and high DOF, and does not admit analytical solutions. Thus, data based approaches are natural to consider for the control of these systems. In recent work [10], we proposed a novel decoupled data based control (D2C) algorithm for learning to control an unknown nonlinear dynamical system. Our approach introduced a rigorous decoupling of the open-loop (planning) problem from the closed-loop (feedback control) problem. This decoupling allows us to come up with a highly efficient approach to solve the problem in a completely data based fashion. Our approach proceeds in two steps: (i) first, we optimize the nominal open-loop trajectory of the system using a 'blackbox simulation' model, (ii) then we identify the linear system governing perturbations from the nominal trajectory using random input-output perturbation data, and design an LQR controller for this linearized system. We have shown that the performance of D2C algorithm is approximately optimal, in the sense that the decoupled design is near optimal to second order in a suitably defined noise parameter [10]. In this work, we apply the D2C approach to the optimal control of the evolution of material micro-structures governed by the phase field model. For comparison, we consider RL techniques [11]. The past several years have seen significant progress in deep neural networks based reinforcement learning approaches for controlling unknown dynamical systems, with applications in many areas like playing games [12], locomotion [13] and robotic hand manipulation [14]. A number of new algorithms that show promising performance are proposed [15][16] [17] and various improvements and innovations have been continuously developed. However, despite excellent performance on a number of tasks, RL is highly data intensive. The training time for such algorithms is typically very large, and high variance and reproducibility issues mar the performance [18]. Thus, materials microstructure control is intractable for current RL techniques. Recent work in RL has also focused on PDE control using techniques like the Deep Deterministic Policy Gradient(DDPG) with modifications such as action descriptors to limit the large action spaces in infinite dimensional systems [19], as such algorithms tend to degrade in performance with an increase in dimensionality of the action-space. However, such an approach, in general, is not feasible for the Microstructure problem.
The control design of systems governed by PDEs are known to be computationally intractable [20]. Owing to the infinite dimensionality of the system, typically model reduction techniques are used to reduce the size of the state space. The most widely used method is the Proper Orthogonal Decomposition (POD) method that finds a reduced basis for the approximation of the nonlinear PDE using a simulation of the same [21]. However, in general, the basis eigenfunctions can be quite different around the optimal trajectory when compared to an initial trajectory, and an approximation using the latter's reduced basis can lead to unacceptable performance [22], [23]. In our approach, the nominal optimization problem is directly solved which does not need any reduced order modeling. Furthermore, the feedback design is accomplished by identifying a linear time-varying (LTV) model which automatically produces a reduced order model (ROM) of the system based purely on the input output data while facilitating the use of linear control design techniques such as LQG/ILQR on the problem as opposed to having to solve a nonlinear control problem in the methods stated previously. Furthermore, the entire control design for the problem is done purely in a data based fashion since the technique only queries a black box simulation model of the process. The contributions of the paper are as follows: we show how to model the dynamics of a multi-phase microstructure and unveil its underlying structure. We present the application of the D2C approach, and a state of the art RL technique (DDPG) to the control of such microstructure dynamics for the first time in the literature, to the best of our knowledge. Our results show that the local D2C approach outperforms the global RL approach such as DDPG when operating on higher dimensional phase field PDEs. We also show that the D2C approach is robust to noise in the practical regimes, and that global optimality can be recovered in higher noise levels through an open-loop replanning approach. Furthermore, we can exploit the peculiar physical properties of material micro-structure dynamics to optimize the control actuation architecture that leads to highly computationally efficient control that does not sacrifice much performance. We envisage that this work is a first step towards the construction of a systematic feedback control synthesis approach for the control of material micro-structures, potentially scalable to real applications.
The rest of the paper is organized as follows: In Section II, the micro-structure dynamics are expanded upon, describing the different PDEs tested. We propose the decoupled data based control (D2C) approach in Section III. In Section IV, the proposed approach is illustrated through customdefined test problems of varying dimensions, along with testing robustness to noise, followed by improvements to the algorithms which exploit the physics of the system.

II. MICROSTRUCTURE DYNAMICS: CLASSES OF PDES IMPLEMENTED
This section provides a brief overview of the non-linear dynamics model of a general multi-phase micro-structure. We assume an infinitely large, 2-dimensional structure satisfying periodical boundary conditions.

A. The material system
The evolution of material micro-structures can be represented by two types of partial differential equations, i.e., the Allen-Cahn equation [24] representing the evolution of a non-conserved quantity, and the Cahn-Hilliard equation [25] representing the evolution of a conserved quantity. The Allen-Cahn equation has a general form of while the Cahn-Hilliard equation has the form where φ = φ(x, t) is called the 'order parameter', which is a spatially varying quantity. In Controls parlance, φ is the state of the system, and is infinite dimensional, i.e., a spatio-temporally varying function. It reflects the component proportion of each phase of material system. For a two-phase system studied in this work, φ = -1 represents one pure phase and φ = 1 represent the other, while φ ∈ (-1, 1) stands for a combination state of both phases on the boundary interface between two pure phases; M is a parameter related to the mobility of material, which is assumed constant in this study; F is an energy function with a non-linear dependence on φ; γ is a gradient coefficient controlling the diffusion level or thickness of the boundary interface. In essence, the Phase Field Method is one of gradient flow methods, which means the evolution process follows the path of steepest descent in an energy landscape starting from an initial state until arriving at a a local minimum. So, the behavior of micro-structures in phase field modeling highly depends on the selection of energy density function F. For instance, the double-well potential function owing to two minima at φ = ± 1 and separated by a maximum at φ = 0 (as plotted in Fig. 2) will cause the field φ to separate into regions where φ ≈ ± 1, divided by the boundary interface, while the single-well potential function owing to a single minimum at φ = 0 (see Fig. 2) predicts a gradual smoothing of any initial non-uniform φ field, yielding a uniform field with φ = 0 in the long-time limit.
Accordingly, the evolution of material micro-structure can be governed by selecting proper form of energy density function F. Controlling the evolution process of micro-structures represented by Allen-Cahn equation and Cahn-Hilliard diffusion equation are completely distinct. The former one is governed through generating or deleting order parameter straightforwardly, while the latter is done through guiding the transport and redistribution of the conserved order parameter across the whole domain. In this study, we adopt the following general form of energy density function F: Here, the parameters T and h are not governed by the A-C and C-H equations. Instead, they may be externally set to any value and may be spatially as well as temporally varying.

B. Discretized model for PDE control
The phase-field is discretized into a 2D grid of dimension N × N . Let φ i,j t denote the phase field of the cell defined by grid location (i, j) for all i, j = {1, ..., N }, at time t.
From Eq. (1) and (3), we have the A-C Equation of the form Now, we apply the above discretization to Eq. (4), and use a central-difference scheme to get the high-dimensional 2D phase-field model. Re-arranging the terms gives the timevariation of the local phase-field: C. Minimization of energy function F for controlling microstructure dynamics The control variables T and h govern the evolution process of micro-structure by influencing the slope and minimums of energy density function F as shown in Fig. 3. Suppose that h = 0, a positive T yields a double-well F while negative T gives a single-well F , as seen in Fig. 3(a). As explained earlier, under h = 0, T determines the Φ phase will go to a uniform state or phase-separating state in the long-time limit. When h = 0, the symmetric energy profile shown in Fig. 3(a) becomes skewed as shown in Fig. 3(b). Under this condition, the parameter h determines which phase will be favored during evolution. If h > 0, the minimum locates on the right hand side, and therefore the phase represented by Φ = 1 is favored while, if h < 0, then the Φ field favors the phase Φ = −1. In most applications, T and h represent specific input conditions, i.e., representing temperature and an external force field (e.g. electromagnetic or mechanical), respectively. Since we aim at investigating the control of material microstructure evolution during processing, we suppose that T and h can be adjusted to any value instantaneously, and that the feedback to the controller is the entire phase field Φ (please see the Simulations section for details).

III. DECOUPLED DATA BASED CONTROL (D2C) ALGORITHM
In this section, we propose a data-based approach to solve the material phase-separation control problem. In particular, we apply the so-called decoupled data-based control (D2C) algorithm that we recently proposed [10], and extend it for high dimensional applications. The D2C algorithm is a highly data-efficient Reinforcement Learning (RL) method that has shown to be much superior to state of the art RL algorithms such as the Deep Deterministic Policy Gradient (DDPG) in terms of data efficiency and training stability while retaining similar or better performance. In the following, we give a very brief introduction to the D2C algorithm (please see [10], [26], [27] for the relevant details).
Let us reconsider the dynamics of the system given in Eq. (1) and (2) and rewrite it in a discrete time, noise perturbed state space form as follows: where Φ t is the state, U t is the control, w t is a white noise perturbation to the system, and is a small parameter that modulates the noise in the system. Suppose now that we have a finite horizon objective function: where c(Φ, u) denotes a running incremental cost, Ψ(Φ) denotes a terminal cost function and E[·] is an expectation operator taken over the sample paths of the system. The objective of the control design is then to design a feedback policy U t (Φ t ) that minimizes the cost function above and is given by The D2C algorithm then proposes a 3 step procedure to approximate the solution to the above problem.
First, a noiseless open-loop optimization problem is solved to find an optimal control sequence,Ū * t , that generates the nominal state trajectoryΦ * t : subject to the zero noise nominal dynamics: Second, the perturbed time varying linear system around the nominal given by δΦ t+1 = A t δΦ t + B t δU t is identified in terms of the time varying system matrices A t , B t .
Third, an LQR controller for the above time varying linear system is designed whose time varying gain is given by K t . Finally, the control applied to the system is given by U t (Φ t ) =Ū * t + K t δΦ t . The following subsections provide details for each of the above-mentioned steps of the D2C algorithm.

A. Open Loop Trajectory Optimization
We utilize an ILQR based method to solve the openloop optimization problem. ILQR typically requires the availability of analytical system Jacobian, and thus, cannot be directly applied when such analytical gradient information is unavailable (much like Nonlinear Programming software whose efficiency depends on the availability of analytical gradients and Hessians). In order to make it an (analytical) model-free algorithm, it is sufficient to obtain estimates of the system Jacobians from simulations, and a sample-efficient randomized way of doing so is described in the following subsection. Since ILQR is a well-established framework, we skip the details and instead present pseudocode in algorithm 1.
1) Estimation of Jacobians: Linear Least Squares by Central Difference (LLS-CD): Using Taylor's expansions of the dynamics 'h', where h(Φ, U ) = f (Φ) + g(Φ)U is the non-linear model of Section 2), about the nominal trajectory (Φ t ,Ū t ) on both the positive and the negative sides, we obtain the following central difference equation: . Multiplying by δΦ t T δU t T on both sides to the above equation and apply standard Least Square method: where 'n s ' be the number of samples for each of the random variables, δΦ t and δU t . Denote the random samples as δX t = δΦ 1 t δΦ 2 t . . . δΦ ns t , δU t = δU 1 t δU 2 t . . . δU ns t and δY t = δX t δU t . We are free to choose the distribution of δΦ t and δU t . We assume both are i.i.d. Gaussian distributed random variables with zero mean and a standard deviation of σ. This ensures that δY t δY T t is invertible. Let us consider the terms in the matrix δY t δY T . From the definition of sample variance, for a large enough n s , we can write the above matrix as

B. Linear Time Varying System Identification
Closed-loop control design in step 2 of D2C requires the knowledge of the linearized system parameters A t and B t for 0 ≤ t ≤ T −1. Here we use the standard least square method to estimate these parameters from input-output experiment data.
First start from the perturbed linear system about the nominal trajectory and estimate the system parameters A t and B t from: δΦ t+1 =Â t δΦ t +B t δU t , where δΦ (n) t is the state perturbation vector and δU (n) t is the control perturbation vector we feed to the system at step t, n th simulation. All the perturbations are zero-mean, i.i.d, Gaussian noise with covariance matrix σI. σ is a o(U ) small value selected by the user. δΦ (n) t+1 denotes the deviation of the output state vector from the nominal state after propagating for one step.
Run N simulations for each step and collect the inputoutput data: Y = [Â t |B t ]X and write out the components: Finally, using the standard least square method, the linearized system parameters are estimated as

C. Closed Loop Control Design
Given the estimated perturbed linear system, we design a finite horizon, discrete time LQR [28] along the trajectory for each time step to minimize the cost function, and use the standard Riccati equations to solve for the feedback gains.

D. D2C Algorithm: Summary
The Decoupled Data Based Control (D2C) Algorithm is summarized in Algorithm 2.

IV. SIMULATION RESULTS
In this section, we compare the training and performance of the data based D2C approach, with the RL-based Deep Deterministic Policy Gradient(DDPG), on a material with control inputs as the temperature and external field inputs on subsets of grid points.

A. Structure and Task
We simulated the phase-separation dynamics in Python, through calling an explicit, second-order solver subroutine in FORTRAN. The system and its tasks are defined as follows: 2) Linearize and identify the LTV system parameters (Â t ,B t ) via least square (Section III-B).
3) Solve the Riccati equations for each step along the nominal trajectory for feedback gain {K t } T −1 t=0 . 4) Apply the closed-loop control policy, while t < T do end while a) Material Micro-structure: The material model simulated consists of a 2-dimensional grid of dimensions 10x10, 20X20 and 50X50 i.e., 100, 400 and 2500 fully observable states respectively. Furthermore, the model is solved using an explicit, second order, central-difference based scheme. The control inputs (T, h) are actuated at each grid point separately; thus, the number of control variables is twice the total number of observable states(200, 800, 5000 resp.). The control task is to influence the material dynamics to converge at a 1) Banded phase-distribution ( Fig.4(a)), with 100 state variables and 200 control channels 2) Checkerboard phase-distribution ( Fig.4(b)), with 400 state variables and 800 control channels 3) A custom phase-distribution (Fig.4(c)), with 2500 state variables and 5000 control channels starting from a uniform initial condition φ 0 = 0.
B. Algorithms Tested a) Deep Deterministic Policy Gradient: Deep Deterministic Policy Gradient (DDPG) is a policy-gradient based off-policy reinforcement learning algorithm that operates in continuous state and action spaces. It relies on two function approximation networks one each for the actor and the critic. The critic network estimates the Q(s, a) value given the state and the action taken, while the actor network engenders a policy given the current state. Neural networks are employed to represent the networks.

C. Training and Testing
D2C implementation is done in three stages corresponding to those mentioned in the previous section, and a black box phase field model is simulated in Python.
Training: The open-loop training plots in Fig. 5 show the cost curve during training. After the cost curves converge, we get the optimal control sequence that could drive the systems to accomplish their tasks. The training parameters and outcomes are summarized in (Tables I, II, III). We note that the training time is dominated by the Python simulation code, and can be made significantly faster with an optimized and parallellized code. Thus, we envisage this work as an initial proof of the concept for controlling material microstructures.
The open-loop optimal trajectory learned by the D2C algorithm after convergence (for goal state-III, Allen-Cahn PDE) is shown in Figure 6.   Testing Criterion: For the closed-loop design, we proceed with the system identification and feedback gain design step of the D2C algorithm mentioned in the previous section to get the closed-loop control policy. For testing, we compare the performance between the open-loop D2C control policy and the closed-loop D2C control policy under different noise levels. The open-loop control policy is to apply the optimal control sequence solved in the open-loop training step without any feedback. Hence, the perturbation drives the model off the nominal trajectory and increases the episodic cost as the noise level increases. Zero-mean Gaussian i.i.d. noise is added to every control channel at each step. The standard deviation of the noise is proportional to the maximum control signal in the designed optimal control sequence for D2C. The feedback policy obtained appears to be near-optimal and highly robust for low noise levels (Fig 7), which is as expected from the O( 4 )-optimality of the algorithm used [27]. But when operating in high noise regimes, the robustness of the feedback policy quickly degenerates. Thus, to reduce the variance of the policy in these domains, we advocate the use of a rapid open loop solver along with replanning at every step, i.e. a Model-Predictive Control approach, to solve the RL problem, since MPC should be able to recover the global optimality for these problems [29].
The recursive MPC algorithm implemented, which uses a fast and reliable local-planner, is summarised in Algorithm 3, and Figure 8 compares the robustness to noise between the stochastic policy obtained from the recursive MPC with the decoupled closed-loop feedback policy at varying noise regimes. It must be noted that the range of noise levels (at least up until about 100 % of the maximum nominal control signal) that we are considering here is far beyond what is typically encountered in practical scenarios. As for the criterion for performance, we use episodic cost at each noise level. 100 rollouts are simulated at every noise level tested. Tables I and  II shows that D2C learns the optimal policy significantly faster than DDPG. The primary reason for this disparity is the feedback parametrization of the two methods: the DDPG deep neural nets are complex parametrizations that are difficult to search over, when compared to the highly compact open-loop + linear feedback parametrization of D2C, as is apparent from Table III. We observe from experiments that the DDPG algorithm is not able to converge to the goal state for the material system of state-dimension over algorithm can quickly generate the closed-loop feedback solution even for the more complex systems.

V. CONCLUSIONS
In this article, we have provided an overview of the modeling of multi-phase micro-structures in accordance with the Allen-Cahn and Cahn-Hiliard equations. We have also compared learning/ data-based approaches to the control of such structures, outlining their relative merits and demerits. In future work, we shall further develop a hybrid modeldata based method which can exploit a priori knowledge of the dynamics to significantly reduce the computation time. We seek to implement a smarter actuator-selection scheme to reduce action space, which may improve viability of DDPG for higher order systems, while also significantly improving performance of D2C. We also plan to extend the methodology for higher-dimensional materials, leading up to an implementation on a realistic full-scale model.

VI. ACKNOWLEDGEMENTS
This work was supported by the NSF CDS&E program under grant number 1802867.