We present control strategies of a diffusion process for chemical vapor deposition for metallic bipolar plates. In the models, we discuss the application of different models to simulate the plasma-transport of chemical reactants in the gas-chamber. The contribution are an optimal control problem based on a PID control to obtain a homogeneous layering. We have taken into account one- and two-dimensional problems that are given with constraints and control functions. A finite-element formulation with adaptive feedback control for time-step selection has been developed for the diffusion process. The optimization is presented with efficient algorithms. Numerical experiments are discussed with respect to the diffusion processes of the macroscopic model.

1. Introduction

We motivate our studying on simulating a low-temperature low-pressure plasma that can be found in chemical vapor deposition (CVD) processes. In the last years the research and optimization in producing high-temperature films by depositing low-pressure processes have increased by using simulation tools, see [1, 2]. Theoretical models exist for deposition processes and can be modeled by coupled transport and flow equations, see [3, 4]. Further interest on standard applications to deposit titanium-nitrogen (TiN) and titanium-carbon (TiC) on metallic layers are immense, see [5]. Recently more and more focus on deposition with new material classes known as MAX-phases are becoming important, see [6, 7]. The MAX-phase are nanolayered terniar metal-carbides or -nitrids, where 𝑀 is a transition metal, 𝐴 is an A-group element (e.g., Al, Ga, In, Si, etc.) and 𝑋 is 𝐶 (carbon) or 𝑁 (nitride), see [8]. Such materials combine ceramic and metallic behavior and can be implanted in the metallic bipolar plates to obtain a new material with noncorrosive and good metallic conductivity behavior.

We discuss a model for low-temperature and low-pressure plasma that can be used to implant or deposit thin layers of important materials, see [9]. This model is used for the implantation process. First, we present the process in the plasma-reactor that transport the contaminants to the wafer surface, see [10]. We deal with a continuous flow model, while we assume a vacuum- and a diffusion-dominated processes. Second, the process at the wafer-surface is modeled by the heavy particles problem with underlying drift. This model deals more with the atomic behavior and we do not allow 𝑝=0, see [10].

To solve such optimization problems, we present a PID-controller (proportional, integral, differential) to control our deposition process, see [11]. We improved heuristic methods of deriving the PID parameters, while we discuss the posteriori error estimates respecting the time-step size control. Our contribution is a modified automatically step-size control and a best approximation is obtained with the time-dependent control method based on the Chien-Hrones-Reswick algorithm, see [12].

Numerical methods are described in the context of time- and spatial-discretization methods for the mesoscopic-scale model. We discussed different experiments and their convergence rates.

For the simulations we apply analytical and also numerical methods to obtain results to control the grow of thin layers.

The paper is outlined as follows. In Section 2 we present our mathematical model and a possible reduced model for the further approximations. In Section 3 we discuss the theoretical background for the simulation of CVD processes. The optimal control and their control paths based on the PID-control approach are discussed in Section 4. The software and program-tools are discussed in Section 5. The numerical experiments are given in Section 6. In the contents that are given in Section 7, we summarize our results.

2. Mathematical Model

In the following, the models are discussed in two directions of far-field and near-field problems:

(1)reaction-diffusion equations, see [13] (far-field problem);(2)Boltzmann-Lattice equations, see [9] (near-field problem). The modeling is considered by the Knudsen Number (Kn), which is the ratio of the mean free path 𝜆 over the typical domain size 𝐿. For small Knudsen numbers Kn0.011.0, we deal with a Navier-Stokes equation or with the convection-diffusion equation, see [5, 14], whereas for large Knudsen numbers Kn1.0, we deal with a Boltzmann equation, see [4].

2.1. Modeling with Partial Differential Equations

Dynamic processes with modifications in time and space will be reshaped by partial differential equations. There is (i) the PDE-formula itself which describes the physical laws of nature that influence the process and (ii) initial and boundary conditions in which specific characteristics of the process, like boundary behavior, can be coded.

There are two types of boundary conditions, namely, Dirichlet and Neumann boundary. With the Dirichlet type the exact value of the boundary is known, however, with Neumann boundaries the spatial derivation of the boundary values in normal direction is known, see an example of the boundary conditions in Figure 1.

2.2. Model for Optimal Control of the Layer

We will concentrate us on a continuum model of mass transportation and assume that the energy and momentum is conserved, see [13]. Therefore, the continuum flow of the mass can be described as diffusion reaction equation given as𝜕𝑡𝑐𝐷𝑐𝑅𝑔=0,inΩ×[0,𝑇],(2.1)𝑐(𝑥,0)=𝑐0(𝑥),onΩ,(2.2)𝜕𝑐(𝑥,𝑡)𝜕𝑛=𝑐1(𝑥,𝑡),on𝜕Ω×[0,𝑇],(2.3)where 𝑐 is the molar concentration, 𝐷 is the diffusion parameter, and 𝑅𝑔 is the reaction and source term.

We modify our model equation (2.1) to a control problem with an additionally right-hand side source:𝜕𝑡𝑐𝐷𝑐=𝑐source,inΩ×[0,𝑇],𝑐(𝑥,0)=𝑐0(𝑥),onΩ,𝜕𝑐(𝑥,𝑡)𝜕𝑛=𝑐1(𝑥,𝑡),on𝜕Ω×[0,𝑇],(2.4)where 𝑐source(𝑥,𝑡)=𝑅𝑔 is the discontinuous or continuous source term of the concentration 𝑐 and we neglect a reaction term of this concentration.

We assume an optimal concentration at the layer𝑐opt(𝑥,𝑡),(2.5) where the layer is given as 𝑥Ωlayer and our constraints are given as𝑐source,min𝑐source𝑐source,max.(2.6)

Additionally, we have to solve the minimization problem:min(𝐽(𝑐,𝑐source1))=2𝑇Ωlayer||𝑐(𝑥,𝑡)𝑐opt||(𝑥,𝑡)2𝜆𝑑𝑥𝑑𝑡+2𝑇Ω||𝑐source||(𝑥,𝑡)2𝑑𝑥𝑑𝑡,(2.7)where 𝑇 is the time period of the process.

Remark 2.1. We choose the 𝐿2-error to control our minimization problem. In literature, see [15, 16], there exists further control-errors, which respect the time behavior.

In a first part, we only solve the transport equation with UG software-tool (unstructed grid software, see [17]) and try to find out the optimal control of the sources to obtain the best homogeneous layer.

In a second part, we consider the optimal control problem and solve also the backward problem.

3. Theoretical Background for Simulation of Diffusive CVD Processes

In what follows we discuss the approximation methods and errors for the simulation of the CVD processes.

3.1. Approximation and Discretization

For the numerical solutions we need to apply approximation methods, for example, finite-difference methods and iterative solver methods for the nonlinear differential equations, see [18, 19].

The finite-element discretization is based on Ω the variational boundary value problem reduces to find 𝑢𝑉 satisfying the initial condition 𝑢(0) such thatΩ𝜕𝑢𝑣𝜕𝑡+𝐷𝑢𝑣𝑑𝑥=0,𝑣𝑉.(3.1)We define the minimal length of triangle which we get from the spatial-discretization with Δ𝑥.

This leads to the following linear semidiscretized system of ordinary differential equations:𝑀𝑑𝑢𝑑𝑡+𝐴𝑢=0,(3.2)where 𝑀 is the mass and 𝐴 the M-matrix.

Here we have taken into account the Courant-Friedrichs-Levy- (CFL-) condition, which is given asCFL=2𝐷maxΔ𝑡min𝑒𝐸Δ𝑥2𝑒,(3.3)where 𝐷max is the maximal diffusion parameter, 𝐸 is the set of the edges of the discretization. We restrict the CFL-condition to 1, if we use an explicit time-discretization and can lower the condition, if we use an implicit discretization.

For the explicit time-discretization, we apply explicit Euler or Runge-Kutta methods.

We use the explicit lower-order Runge-Kutta methods:0121201(3.4)

Furthermore we use the following Heun method (third-order):0131323023014034(3.5)

The implicit time-discretization is done with implicit Euler or Runge-Kutta methods.

Here, we use the implicit trapezoidal rule:0112121212(3.6)

Furthermore we use the following Gauss Runge-Kutta method:123614143612+3614+36141212(3.7)

Remark 3.1. We apply implicit time-discretization methods for the pure diffusion part, where we apply explicit time-discretization methods for the pure convection part. Here we have to respect the CLF-condition, see [20].

3.2. Errors and Convergence Rate

For studying the errors and the convergence-rates in our test example, we have to define the following norm in two space-dimensions:

(i)discrete 𝐿max-norm:err𝐿max,Δ𝑥,Δ𝑡=𝑝max𝑖=1||𝑐num(𝑥𝑖,𝑇)𝑐ref(𝑥𝑖||,𝑇),(3.8)(ii) discrete 𝐿1-norm: err𝐿1,Δ𝑥,Δ𝑡=𝑝𝑖=1Δ𝑥2||𝑐num(𝑥𝑖,𝑇)𝑐ref(𝑥𝑖||,𝑇),(3.9)(iii) discrete 𝐿2-norm: err𝐿2,Δ𝑥,Δ𝑡=𝑝𝑖=1Δ𝑥2||𝑐num(𝑥𝑖,𝑇)𝑐ref(𝑥𝑖||,𝑇)2,(3.10) where Δ𝑥 is the spatial-step of the discretization, Δ𝑡 is the time-step of the discretization, and 𝑇 is the end-time of the computation. 𝑝 is the number of grid points in the discretization method. 𝑐num is the numerical solution and 𝑐ref is the reference solution, computed at fine spatial- and time-grids.

The numerical convergence rate are given as follows.

(i)For the spatial error, we define 𝜌𝐿2,Δ𝑥1,Δ𝑥2,Δ𝑡=logerr𝐿2,Δ𝑥1,Δ𝑡/err𝐿2,Δ𝑥2,Δ𝑡log(Δ𝑥1/Δ𝑥2),(3.11)where Δ𝑥1 is the coarse, Δ𝑥2 is the fine spatial grid-step, and Δ𝑡 is the time-grid step for both results.(ii) For the time error, we define 𝜌𝐿2,Δ𝑥,Δ𝑡1,Δ𝑡2=logerr𝐿2,Δ𝑥,Δ𝑡1/err𝐿2,Δ𝑥,Δ𝑡2log(Δ𝑡1/Δ𝑡2),(3.12) where Δ𝑡1 is the coarse, Δ𝑡2 is the fine-time-step, and Δ𝑥 is the spatial-grid step for both results.

We often use Δ𝑥2=Δ𝑥1/2. In this case, we have 𝜌𝐿2,Δ𝑥1,Δ𝑥2,Δ𝑡=𝜌𝐿2,Δ𝑥1,Δ𝑡. Further we have to choose Δ𝑥1 with respect to the Δ𝑡𝐼=[0,Δ𝑡max], which have maximal 𝜌. Thus we define ArgMax(Δ𝑥):ArgMax(Δ𝑥)=argmaxΔ𝑡𝐼𝜌𝐿2,Δ𝑥,Δ𝑥/2,Δ𝑡.(3.13)

4. Optimal Control Methods

Here we discuss the control of a diffusion equation with a feedback based on a PID-controller.

4.1. Forward Controller (Simple P-Controller)

The first controller we discuss is the simple P-controller, see [11]. A first idea is to control linearly the error of the solved PDE.

In Figure 3, we present the P-controller.

Our control problem is given with the control of the error to the optimal concentration of the layer and correct the source-flux:𝜕𝑡𝑐𝐷𝑐=𝑐source,inΩ×[0,𝑇],𝑐(𝑥,0)=𝑐0(𝑥),onΩ,𝜕𝑐(𝑥,𝑡)𝜕𝑛=𝑐1(𝑥,𝑡),on𝜕Ω×[0,𝑇],(4.1)where 𝑐source(𝑥,𝑡) is a discontinuous source flow of the concentration 𝑐.

We assume an optimal concentration at the layer with the concentration 𝑐opt(𝑥,𝑡), where the layer is given as 𝑥Ωlayer and 0 elsewhere.

Our constraints are bounded as𝑐source,min𝑐source𝑐source,max.(4.2)

Remark 4.1. Taken into account the hysteresis of the deposition process, we apply a linear increase of 𝑐source,max in the optimal control with respect to time, see Figure 4.

4.2. PID-Controller

The PID-controller is used to control temperature, motion, and flow. The controller is available in analog and digital forms, see [16]. The aim of the controller is to get the output (velocity, temperature, position) in the area of the constraint output, in a short time, with minimal overshoot, and with small errors.

We have three elements in the PID-control, where 𝑃 is the proportional part, 𝐼 the integral part, and 𝐷 is the derivative part of the controller, see [16].

These terms describe three basic mathematical functions applied to the error signal, error=𝐶optimal𝐶computed.

The errors represented the difference between constraint (optimal set) and computed results in the simulation.

To accelerate a PID-controller means to adjust the three multipliers 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷 adding in various amounts of these functions to get the system to behave the way you want, see [11].

Figure 5 summarizes the PID terms and their effect on a control system.

Initialization of the PID-Controller
The algorithm of the initialization of the PID-control (i.e., search 𝐾𝑃,𝐾𝐼,𝐾𝐷) is given as in Algorithm 4.2 (see [15]).

Algorithm 4.2. (1) We initialize the P-controller: 𝐾𝐼=0.0,𝐾𝐷=0.0.
(2) The amplifying factor 𝐾𝑃 is increased till we reached the permanent oscillations as a stability boundary of the closed control system.
(3) We obtain for 𝐾𝑃 the critical value 𝐾𝑃,crit..
(4) The period-length of the permanent oscillation is given as 𝑇crit.
(5) We obtain the parameters from Table 1.

Further we compute the rest parameters as 𝐾𝐼=𝐾𝑃/𝑇𝑛,𝐾𝐷=𝐾𝑃𝑇𝑣, see [21].

4.3. Adaptive Time-Control

Often the heuristic assumptions of the PID-parameters are too coarse.

One can improve the method by applying an adaptive step-size control.

We discuss the step-size control with respect to our underlying error, that is, given by the computed and optimal output of our differential equation.

Based on the adaptive control, we can benefit to accelerate the control problem.

According to Hairer and Wanner [22], we apply the automatic control problem with a PID-controller.

The automatically step-size is given as (see [11])Δ𝑡𝑛+1=𝑒𝑛1𝑒𝑛𝐾𝑃tol𝑒𝑛𝐾𝐼𝑒2𝑛1𝑒𝑛𝑒𝑛2𝐾𝐷Δ𝑡𝑛,(4.3)where tol is the tolerance, 𝑒𝑛 is the error of the quantities of interest in time-step Δ𝑡.

We can control the step-size with respect to our heuristically computed 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷 parameters. Initialization of the adaptive control can be seen in Algorithm 4.3.

Algorithm 4.3. (1) Define Tolerance, Min and Max of the concentration.
(2) Apply the parameters: 𝐾𝑃,𝐾𝐷,𝐾𝐼 form a first run.
(3) Optimize the computations with a first feedback.

4.4. Identification of the Control Path

In the forward problem, we computed a PDE to simulate the CVD process in the underlying domain. For our control problem, we have the behavior of two points in the underlying domain, 𝐶source our source and 𝐶𝑥 the source restricted to the response of the controlled system. To analyze the differences between given source the 𝐶source and the computed source in the control process 𝐶𝑥, we study the step response from our system:𝑥𝑎𝐶(𝑡)=𝑥(𝑡)𝐶Source(𝑡)=𝑥𝑒0.(4.4)

Algorithm 4.4. (1) We determine the model of the control path, for example, PT1 (Proportional time 1), PT2 (Proportional time 2), see [16]. For that purpose we investigate the step response 𝑥𝑎, his first and second derivate [16, pages 117 and 331].
(2) We determine the parameters of our model: Kp, T1, T2.
(3) Our goal is to control the system with a controller. Also we have to determine the control-parameter [16, page 405].

Identification of the Control Path
With the step response 𝑥𝑎, we identify the control path as a PT2 element.
Because lim𝑡𝑥𝑎(𝑡)=constant (also 0 nor ), the basic behavior of our element is 𝑃 (and neither 𝐷 nor 𝐼). Now we have to determine the number of time delays in the control path.

Identification of the Parameters
For such control path which has PT2 behavior, there exists a model as ordinary differential equation. Here a generalization can be obtained with four parameters and be divided by 𝑎0 a normed form with three parameters 𝜔0= angular frequency, 𝐴 = attenuation, and 𝐾𝑆 = proportional factor:𝑎2̈𝑥𝑎(𝑡)+𝑎1̇𝑥𝑎(𝑡)+𝑎0𝑥𝑎(𝑡)=𝑏0𝑥𝑒0,1𝜔20̈𝑥𝑎(𝑡)+2𝐴𝜔0̇𝑥𝑎(𝑡)+𝑥𝑎(𝑡)=𝐾𝑆𝑥𝑒0.(4.5)The solution of the differential equation results in the characteristic equation:𝛼2𝜔20𝛼+2𝐴𝜔0+1=0with𝛼1,2=𝐴𝜔0±𝜔0𝐴21.(4.6)
So that we have the analytical solution for the normed step response:𝑥𝑎(𝑡)𝑥𝑒0=𝐾𝑆𝑇11𝑇1𝑇2𝑒𝑡/𝑇1+𝑇2𝑇1𝑇2𝑒𝑡/𝑇2,with𝛼11=𝑇1,𝛼21=𝑇2.(4.7)We go now the inverse way to identify the system. We start with the step response. Using the inflected tangent at (𝑡𝑤,𝑥𝑤,̇𝑥𝑎(𝑡𝑤)) and the limit 𝑥𝑎(𝑡), we obtain our parameters.
Firstly, we have 𝐾𝑆=𝑥𝑎(𝑡)/𝑥𝑒0, 𝑡𝑔=𝐾𝑆(𝑥𝑒0/̇𝑥𝑎(𝑡𝑤)), and 𝑡𝑢=𝑡𝑤𝑡𝑔(𝑥𝑎(𝑡𝑤)/𝑥𝑎(𝑡)), see Figure 6. In the second step, we obtain 𝑇1 and 𝑇2 by the following formulas:𝑡𝑢𝑡𝑔=𝛼𝛼/(1𝛼)(𝛼ln(𝛼)+𝛼21)𝑡𝛼11,𝑔𝑇1=𝛼𝛼/(1𝛼)𝑇,𝛼=2𝑇1.(4.8)

Control Parameters
The parameter occurrence also in the transfer function 𝐺(𝑠), which is the Laplace transformation of the differential equation. Our first choice to determine the control parameter is here the PT2 element. Alternatively, we use an approximation to a PT1 element, but here we have additional a scan-time parameter:𝐾𝐺(𝑠)=𝑆(1+𝑇1𝑠)(1+𝑇2𝑠)PT2element,𝐺(𝑠)𝐾𝑆𝑒𝑇𝑢𝑠1+𝑠𝑇𝑔PT1element.(4.9)Independent from the behavior which is being wished, we obtain values for the parameters 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷, respectively. 𝐾𝑅=𝐾𝑃, 𝑇𝑛 and 𝑇𝑣 are of the PID-controller. In Table 2 are given the values in dependence to the percentage overshoot (PO). There are the values of PO=0, respectively, PO=20. The rest can be determined by linear interpolation. An alternative schema is given by Takahashi. Here we have in addition the parameter scan-time 𝑇. This is a general form of the Nichols-Ziegler method (see Table 1). Further are 𝐾𝑃=𝐾𝑅, 𝐾𝐼=𝐾𝑃/𝑇𝑛, and 𝐾𝐷=𝐾𝑃𝑇𝑣.

Table 3: MATLAB-toolbox pdetool.

5. Software and Program-Tools

MATLAB Functions
We use the MATLAB-toolbox pdetool for the time- and spatial-discretizations, where we have a finite-element method with 𝑃1-elements and an implicit Euler method for the time-discretization.

The PID-controller is also programmed in MATLAB. Our combined code is given in the DEPOSIT-PID toolbox and described in what follows. The DEPOSIT-PID toolbox is manipulated by a graphical user interface, by which simulation- and control-models are chosen and corresponding parameters can be manually adjusted.
The objective is to simulate the diffusion and deposition of the vapor in the apparatus and to obtain the optimal vapor concentration at the measuring point.
We can divide the process into three phases, namely, forward step, control step, backward step, which have a cyclic repetition. In the forward step, the diffusion takes place, which can be simulated by a time-step of the heat equation.
After that, in the control step, the actual concentration at the under boundary (0,1) can be measured (shown by computed in the graph) and compared with the optimal value (shown by optimal).
The control step is followed by the backward step: from the error in control step and the control model, the optimal alteration of the source can be computed. The vapor flows through the source point (0,0) in the apparatus and this value is shown by the SourceOutput in Figure 12. This will be simulated through the addition of SourceOutput to actual concentration at the source point (see Figure 8). Further, as a simplification, we set the concentration at the under boundary to zero, because here the gas transforms into solid matter.

In Section 4 we introduced some fundamentals in order to control the apparatus. The models and corresponding parameters can be altered by the Gui.

The Layout of the DEPOSIT-PID Gui
This Gui contains the following:
(1)short-time plot (2D) of computed, optimal, and SourceOutput; (2) long-time plot (2D);(3)listbox with names parameters:(a)textbox with actual value of the parameter chosen in [12], (b)textbox to change the value of the parameter chosen in [12]; (4)listbox with names parameters:(a)textbox with actual value of the parameter chosen in [23], (b)textbox to change the value of the parameter chosen in [23]; (5)3D plot of distribution;(6)3D grid plot of distribution;(7)listbox with names of parameters: (a)checkbox with actual value of the parameter chosen in [24];(8) push button: save;(9) push button: reset;(10)radio button: start.

Another software-tool, KONTOOL, is programmed to compute the numerical convergence rates of the applications. The software-tool has implemented the errors and convergence rates defined in Section 3.2. An error-analysis based on successive refinement of space and time is done and the resulting errors and convergence rates are computed. Optimal convergence rates with respect to balance the time- and spatial-grids are calculated, see Algorithm 5.2.

Remark 5.1. The software-tool can be modified and applied to arbitrary spatial- and time-discretization methods. The interface of KONTOOL needs at least the parameters of the spatial- and time-grid and the starting parameters of the underlying methods.

Algorithm 5.2 (The algorithm of the computation of the numerical convergence tableau). (1) We compute reference solutions: (a) numerically: fine time and spatial steps or (b) analytically (if there exists an analytical solution).
(2) We apply one spatial discretization of step Δ𝑥 and apply all time discretization with steps Δ𝑡, where the coarsest Δ𝑡 is given by the CFL condition or till first non-numerical results as oscillations. We compute the error 𝑐num𝑐ref in the 𝐿2-norm.
(3) We continue the next fine spatial steps, for example, Δ𝑥/2.
(4) We compute the convergence tableau with time and space.

In the next section, we discuss the numerical experiments.

6. Experiment for the Plasma Reactor

In this section, we present our numerical experiments for the CVD processes in a plasma reactor.

6.1. Simulation of a Diffusion Equation with Analytical Solution (Neumann Boundary Conditions)

Here we simulate a diffusion equation with Neumann boundary conditions and right-hand side 0. Our control problem has only the forward problem to solve and we consider the accuracy of our simulations.

We have the following equation:𝜕𝑡𝑐𝛽2(𝜕𝑥𝑥+𝜕𝑦𝑦𝜕)𝑐=0,inΩ×[0,𝑇],𝑐(𝑥,𝑦,0)=cos(2𝑥)+cos(2𝑦),onΩ,𝑛𝑐(𝑥,𝑦,𝑡)=0,on𝜕Ω×[0,𝑇],(6.1)where 𝑐 is the molar concentration, Ω=[1,1]×[1,1], and 𝑡(0,𝑇). 𝐷=𝛽2 is the diffusion parameter of the diffusion equation.

We have the following analytical solution:𝑐ana(𝑥,𝑦,𝑡)=sin(2)+𝑛=1𝐴𝑛exp(𝛽2𝑛2𝜋2𝑡)(cos(𝑛𝜋𝑥)+cos(𝑛𝜋𝑦)),(6.2)where 𝐴𝑛=(1)𝑛(4sin(2)/(𝑛2𝜋24)). We apply the diffusion coefficient 𝐷=0.01, respectively, 𝛽=0.1. We obtain the following result after time 𝑇=12.0, see Figure 11.

We see that for large Δ𝑥 the numerical solution converge faster to the stable constant endsolution (offset) than solution for smaller Δ𝑥 and the analytical solution. The error between the constant analytical endsolution (offset, sin(2)=0.9093) is also greater than for smaller Δ𝑥. The 𝐿2-error is given in Table 4.

Remark 6.1. We test for the pure diffusion equation our underlying discretization methods and apply finite elements for the spatial-discretization and implicit Runge-Kutta methods for the time-discretization. In the results, we obtain decreasing errors for the different time- and spatial-steps.

6.2. Simulation of an Optimal Control of a Diffusion Equation with Heuristic Choise of the Control Parameters

Here we simulate a first example of a diffusion equation and control the concentrations in the deposition process.

We have the following equation:𝜕𝑡𝑐𝐷𝑐=𝑓(𝑡),inΩ×[0,𝑇],𝑐(𝑥,0)=𝑐0(𝑥),onΩ,𝜕𝑐(𝑥,𝑡)𝜕𝑛=𝑐1(𝑥,𝑡),on𝜕Ω×[0,𝑇],(6.3)where 𝑐 is the molar concentration, 𝐷 is the diffusion parameter of the diffusion equation, and 𝑓(𝑡) is the right-hand side or diffusion source. We have the following constraint: 𝑐optimal(𝑥point,𝑦point)=0.5, where (𝑥point,𝑦point) is the control point in our domain. The parameters are given as 𝐷=0.01 and 𝑓(𝑡)=𝑎𝑡+𝑏, so we deal with a linear source, 𝑎=0.2 and 𝑏=0.1 are constants. In the following tests, we propose the 3 possibilities to control the optimal temperature:

(i)P-control with constant optimal constraint,(ii)PID-control with constant optimal constraint,(iii)PID-control with Linear optimal constraint.

The results for the control methods are given in Figure 12.

6.3. Preliminary Remark for Simulations for the Convergence Order

To determine the function ArgMax, firstly, we have to determine a practicable interval 𝐼 forArgMax(Δ𝑥)=argmaxΔ𝑡𝐼𝜌𝐿2,Δ𝑥,Δ𝑥/2,Δ𝑡.(6.4)One possibility is the interval𝐼CFL=(0,Δ𝑡CLF],(6.5)with Δ𝑡CFL=Δ𝑠2/2𝐷max from the CFL-condition. In the experiments, we find another interval, where the convergence rate function is relative stable convex:𝐼stable=(0,Δ𝑡stable],(6.6)whereΔ𝑡stable=maxΔ𝑡{𝜌(0,Δ𝑡)isconvex}.(6.7)It is clear that we can in 𝐼stable take some restrictions to a subinterval 𝐼sub=[Δ𝑡𝑠min,Δ𝑡𝑠max]𝐼stable, when we know that some Δ𝑡𝐼sub with 𝜌(Δ𝑡)>𝜌(Δ𝑡𝑠min),𝜌(Δ𝑡𝑠max).

The function of the numerical convergence rate is discrete in the spatial-discretization variable Δ𝑥 since we get a finer discretization with a bisection of Δ𝑥. With a finer discretization, a triangle is replaced by four subtriangles.

In the temporal discretization variable Δ𝑡, we are in contrast not restricted to such conditions and Δ𝑥 could be chosen to any arbitrary value above 0. To get a first glance, we have selected the methods of bisection. Subsequently, we consider finer discretizations in intervals which are of special interest.

6.4. Simulations for the Convergence Order

We consider the convergence rate 𝜌 as a function in dependent of Δ𝑥 and Δ𝑡. This function also depends on some parameters, for instance, the diffusion coefficient 𝐷 and the control parameter 𝜆. We now present the results of two chosen experiments. For the first experiment, the parameters are 𝐷=0.1 and 𝜆=1, while we increase the propagation velocity in such a way that we increase 𝐷 to 1 for the second experiment.

We observe that, for instance, in the case of spatial-discretization Δ𝑥=0.05 and 𝑥=0.25, the maximal convergence rate lies between Δ𝑡=0.05 and Δ𝑡=0.0125. The maximum itself lies close to 0.025. To determine the precise value of Δ𝑡=ArgMax(Δ𝑥), we refine our method in these interesting intervals.

In Figure 13, we observe that there is a relatively stable convex area starting at 0. This area then switches over to an area with strong fluctuations. We now compare the stable area with the CFL-condition. The area where CFL<1 lies in the stable convex area of the function (Δ𝑡CFL<Δ𝑡stable). For Δ𝑥=0.05,𝐷=0.1, the CFL-area ends at Δ𝑡CFL=0.0125, the stable convex area range to about Δ𝑡stable=0.2. For our convergence diagram, we try to find for every Δ𝑥 the value Δ𝑡, where the convergence rate becomes maximal (ArgMax). It is clear that we only use Δ𝑡 at which the convergence rate function is stable: Δ𝑡<Δ𝑡stable.

Remark 6.2. The experiment shows the linear convergence rate of the P-controller with different 𝜆 values. So we obtain a stable method with respect to the P-controller. In the examples, we apply heuristic methods to derive the control parameters for the P- and PID-controller. We show that we have reached the linear order of the underlying finite element discretization method. We have higher control errors if we did not compute the correct control parameters and the numerical errors are smaller than our control error. To prohibit this problem, we have to compute in the next example the control parameters by a feedback equation, see [25].

6.5. Simulation of an Optimal Control of a Diffusion Equation with Adaptive Control

In the second example, we simulate the diffusion equation and control the temperature with and adaptive control based on a PID-controller, see [11].

We have the following equation:𝜕𝑡𝑐𝐷𝑐=𝑓(𝑡),inΩ×[0,𝑇],𝑐(𝑥,𝑡)=𝑐0(𝑥),onΩ,𝜕𝑐(𝑥,𝑡)𝜕𝑛=𝑐1(𝑥,𝑡),on𝜕Ω×[0,𝑇],(6.8)where 𝑐 is the molar concentration, 𝐷 is the diffusion parameter of the diffusion equation, and 𝑓(𝑡) is the right hand side or source.

We have the following constraint: 𝑐optimal(𝑥point,𝑦point)=0.5,(6.9) where (𝑥point,𝑦point)=(0,1) is the control point in our domain.

The automatically step-size is given as (see [11])Δ𝑡𝑛+1=𝑒𝑛1𝑒𝑛𝐾𝑃tol𝑒𝑛𝐾𝐼𝑒2𝑛1𝑒𝑛𝑒𝑛2𝐾𝐷Δ𝑡𝑛,(6.10)where tol=1 is the tolerance, 𝑒𝑛 is the error of the quantities of interest in time-step Δ𝑡.

The errors are given as𝑒𝑛=||𝑢𝑛𝑢𝑛1||||𝑢𝑛||,(6.11)where 𝑢𝑛 is the result at time-step 𝑡𝑛.

The parameters are given as 𝐷=0.1, 𝑃krit=15, 𝑇krit=5, Δ𝑡 (PID-control), tol=1 (adaptive PID-control), Δ𝑡max=0.1 (adaptive PID-control), Δ𝑡min=0.01 (adaptive PID-control) (see Figure 17).

Furthermore, we change the parameter tol=1 to tol=0.1 and Δ𝑡min=0.01 to Δ𝑡min=0.0001 (see Figure 18).

Remark 6.3. In Figures 17 and 18, we see an oscillating time interval at the beginning of the automatically step-size control. In the first experiments, we had only taken into account a previous heuristically computation of the control parameters 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷 before the step-size control for the whole time-interval 0,𝑇. The optimal control parameters are given to the whole time-interval 0,𝑇. A modified algorithm to compute the control parameters for the initialization time-interval 0,𝑡begin and the whole time-interval 0,𝑇 overcome the oscillation problems.

A modified automatically step-size control, which minimize the oscillations is given in the following algorithm.

Algorithm 6.4. (1) We compute the reference control parameters 𝐾𝑃,global, 𝐾𝐼,global and 𝐾𝐷,global for the time-interval 0,𝑇.
(2) We apply the automatically step-size control for the global control parameters with tolglobal, Δ𝑡max,global and Δ𝑡min,global, which can by chosen large.
(3) We stop the computation till we reach the optimal solution and mark remember the time 𝑡oszill.
(4) We compute the local control parameters 𝐾𝑃,local, 𝐾𝐼,local and 𝐾𝐷,local for the time-interval 0,𝑡oszill.
(5) We restart the computation with the local control parameters and smaller step-size parameters tollocal, Δ𝑡max,local and Δ𝑡min,local till we reach 𝑡oszill and continue the computation with the global parameters.
(6) We stop the computation if we reach 𝑡=𝑇. If we obtain also high oscillation with the local parameters, we refine the local interval and go to step (3).

Remark 6.5. The modified automatically step-size control had taken into account the local behavior of the control problem. We could adapt the control parameters 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷 with respect to the local time-intervals. This modified algorithm considers a local time behavior more accurate and reduces oscillations at the initialization process.

6.6. Simulation of an Optimal Control of a Diffusion Equation with Adaptive Control and Recovering of the Control Parameters

In what follows, we present the adaptive control based on Algorithm 6.4.

To be more precise the computation of global control parameters 𝐾𝑃,global=𝐾𝑃(Δ𝑡), 𝐾𝐼,global=𝐾𝐼(Δ𝑡), and 𝐾𝐷,global=𝐾𝐷(Δ𝑡) can be automatized in an interval Δ𝑡(0,𝑇).

We improve the automatically step-size, which is given in Section 4.3 to(Δ𝑡)𝑛+1=𝑒𝑛1𝑒𝑛𝐾𝑃(Δ𝑡)tol𝑒𝑛𝐾𝐼(Δ𝑡)𝑒2𝑛1𝑒𝑛𝑒𝑛2𝐾𝐷(Δ𝑡)(Δ𝑡)𝑛,(6.12)where tol is the tolerance, 𝑒𝑛 is the error of the quantities of interest in time-step 𝑛.

With approximations to each subintervals of (0,𝑇), we derive the control parameters 𝐾𝑃(Δ𝑡), 𝐾𝐼(Δ𝑡), and 𝐾𝑃(Δ𝑡).

The improved globalized control parameters are computed in the following algorithm.

Algorithm 6.6. Computing of the control parameters 𝐾𝑃,global=𝐾𝑃(𝑡), 𝐾𝐼,global=𝐾𝐼(𝑡) and 𝐾𝐷,global=𝐾𝐷(𝑡) for the time-interval 0,𝑇.
(1) Compute the critical control function with the minima and maxima. Based on the optimal oscillated function, we can derive the 𝐾𝑃(𝑇/3) for the interval (0,𝑇/3).
(2) Redo the step (1) for the intervals (𝑇/3,2𝑇/3) and 2𝑇/3,𝑇 and approximate the function 𝐾𝑃(𝑡). Based on the 𝐾𝑃(𝑡) we can derive the 𝐾𝐼(𝑡) and 𝐾𝐷(𝑡) function.
(3) The optimal control parameters are used in Algorithm 6.4.

Remark 6.7. Further ideas to determine the time-dependent parameters 𝐾𝑃(Δ𝑡), 𝐾𝐼(Δ𝑡), and 𝐾𝐷(Δ𝑡) can be done by using our identification method from Algorithm 4.4 in Section 4.4. There we have discussed two methods to determine the parameters from the control path, namely, Chien-Hrones-Reswick and Takahashi. Alternatively, we have the method from Algorithm 4.2 (Nichols-Ziegler) in Section 4.2. For the identification of the control path, we have only look to the step response. The alternative method demanded to find the critical 𝑃-value.

6.7. Determinate of Critical P with Exponential Regression

To initialize the PID-controller, we have to derive the 𝐾𝑃, 𝐾𝐷, and 𝐾𝐼 parameters.

An idea is derived with the step function response, see [15].

This idea has also be included into the time-step control of the PID-controller.

In Figure 19, we can see that the maximum and minimum of the error function, done with the step function response, have a exponential behavior.

So we use the exponential regression𝑦=𝑏exp(𝑎𝑡)+𝜖,(6.13)with 𝜖𝑁(0,𝜎2) normal distributed. Via log the exponential regression can be formed in a linear regression:𝑦=𝑏exp(𝑎𝑡),log(𝑦)=log(𝑏)+𝑎𝑡.(6.14)Before we have to convert the minima into the positive area.

Let 𝜆(𝑛) be the actual lambda value, 𝜆(𝑛)+ the next greater 𝜆, and 𝜆(𝑛) the next smaller one. Firstly, we use 𝜆(𝑛+1)=(𝜆(𝑛)+𝜆(𝑛)+)2𝑎(max,𝑛)<0,𝑎(min,𝑛)>0(𝜆tosmall),(𝜆(𝑛)+𝜆(𝑛))2𝑎(max,𝑛)>0,𝑎(min,𝑛)<0(𝜆togreat).(6.15) In the first step and in limit case, we use 𝜆(𝑛+1)=2𝜆(𝑛), respectively, 𝜆(𝑛+1)=𝜆(𝑛)/2. In the first experiment, we start with 𝜆=359, which we have found with manual tests. The automatic found that the value 396.729 the 𝑎-value of 𝜆=3.59 have 3 zeros, the new less than five.

Remark 6.8. The adaptive step-size control of the PID-controller had taken into account the optimization in each local time-interval. We consider a larger time-interval and derive the optimal control parameters 𝐾𝑃,global, 𝐾𝐼,global, and 𝐾𝐷,global. Based on this parameterization, we apply the adaptation in time to be optimal in such a large time-interval.

6.8. Results of Our Experiments to Find the Time-Dependent Parameters

To automatize the adaptive time-step algorithm, we have to compute the step response and derive the 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷 parameters of this process. Here, we present the experiments in Figure 21, which had taken into account the critical values 𝑃crit and 𝑇crit to derive the optimal 𝐾𝑃, 𝐾𝐷, and 𝐾𝐼 parameters based on Algorithm 6.6.

A further algorithm to derive the control parameters is done by Chien et al., see [12].

The results are given in Figure 23.

The 𝐾𝑃,𝐾𝐷, and 𝐾𝐼 parameters are derived with Algorithm 6.6.

6.9. Experiments with Time-Dependent Parameters

In these experiments, we concentrate on the Chien-Hrones-Reswick method, see [12].

To automatize the adaptive time-step algorithm, we have to compute the step response and derive the 𝐾𝑃, 𝐾𝐼, and 𝐾𝐷 parameters of this process. The minimization problem is presented in Figure 24, which had taken into account the critical parameters.

In Figure 25, we see the step function response and the best fit of this function with the adaptive time control method. Because of the sensitivity to the initial problem, we can overcome the oscillation with the Chien-Hrones-Reswick method.

In Figure 26, we see the comparison of our different adaptive methods. Here the best method is the benefit of the time adaptive control method based on Algorithm 6.6.

Remark 6.9. The benefit of the time-dependent control parameters 𝐾𝑃(Δ𝑡), 𝐾𝐷(Δ𝑡), and 𝐾𝐼(Δ𝑡) are important in the initialization of the control problem. Based on the fine-time scales globally chosen, parameters neglect the local time behavior. In the experiments, we found out the accelerated relaxation behavior and the fast control to the optimal value. Such adaptive algorithms to optimize the control parameters can improve the control methods and had taken into account the local-time scales.

7. Conclusions and Discussions

We present a continuous or kinetic model, due to the fare-field or near-field effect of our deposition process. We discuss the PID-controller to automatize our deposition process. Due to heuristic methods of deriving the PID parameters, we discuss aposteriori error estimates to automatize the time-stepping methods. A modified automatically step-size control is discussed and the best approximations are obtained with the time-dependent control method based on the Chien-Hrones-Reswick algorithm. For the mesoscopic-scale model, we discussed different experiments and their convergence rates. In future, we will analyze the validity of the models with physical experiments.