# Adaptive Step-Size Control in Simulation of Diffusive CVD Processes

Authors: Jürgen Geiser; Christian Fleck
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

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.


---

## Body

## 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 M is a transition
metal, A is an A-group
element (e.g., Al, Ga, In, Si, etc.) and X is C (carbon) or N (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 p=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 Section2 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 L. For small Knudsen numbers Kn≈0.01−1.0, we deal with a Navier-Stokes equation or with the
convection-diffusion equation, see [5, 14], whereas for large Knudsen numbers Kn≥1.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 Figure1.Dirichlet and
Neumann boundary conditions.
(a)(b)Spatial-discretization.
(a)(b)


### 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(2.1)∂tc−∇D∇c−Rg=0,inΩ×[0,T],(2.2)c(x,0)=c0(x),onΩ,(2.3)∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where c is the molar
concentration, D is the
diffusion parameter, and Rg is the reaction
and source term.We modify our model  equation (2.1) to a control problem
with an additionally right-hand side source:(2.4)∂tc−∇D∇c=csource,inΩ×[0,T],c(x,0)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where csource(x,t)=Rg is the
discontinuous or continuous source term of the concentration c and we neglect
a reaction term of this concentration.We assume an optimal concentration at the
layer(2.5)copt(x,t), where the layer is given as x∈Ωlayer and our
constraints are given as(2.6)csource,min≤csource≤csource,max.Additionally, we have to solve the minimization
problem:(2.7)min(J(c,csource)):=12∫T∫Ωlayer|c(x,t)−copt(x,t)|2dxdt+λ2∫T∫Ω|csource(x,t)|2dxdt,where T is the time
period of the process.Remark 2.1.
We choose theL2-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.


## 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 Figure1.Dirichlet and
Neumann boundary conditions.
(a)(b)Spatial-discretization.
(a)(b)


## 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(2.1)∂tc−∇D∇c−Rg=0,inΩ×[0,T],(2.2)c(x,0)=c0(x),onΩ,(2.3)∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where c is the molar
concentration, D is the
diffusion parameter, and Rg is the reaction
and source term.We modify our model  equation (2.1) to a control problem
with an additionally right-hand side source:(2.4)∂tc−∇D∇c=csource,inΩ×[0,T],c(x,0)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where csource(x,t)=Rg is the
discontinuous or continuous source term of the concentration c and we neglect
a reaction term of this concentration.We assume an optimal concentration at the
layer(2.5)copt(x,t), where the layer is given as x∈Ωlayer and our
constraints are given as(2.6)csource,min≤csource≤csource,max.Additionally, we have to solve the minimization
problem:(2.7)min(J(c,csource)):=12∫T∫Ωlayer|c(x,t)−copt(x,t)|2dxdt+λ2∫T∫Ω|csource(x,t)|2dxdt,where T is the time
period of the process.Remark 2.1.
We choose theL2-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Ωh the variational
boundary value problem reduces to find uh∈Vh satisfying the
initial condition uh(0) such
that(3.1)∫Ωh(∂uh∂tvh+D∇uh⋅∇vh)dx=0,∀vh∈Vh.We define the minimal length of
triangle which we get from the spatial-discretization with Δx.This leads to the following linear semidiscretized
system of ordinary differential equations:(3.2)Mdu*dt+Au*=0,where M is the mass and A the M-matrix.Here we have taken into account the
Courant-Friedrichs-Levy- (CFL-) condition, which is given as(3.3)CFL=2DmaxΔtmine∈EΔxe2,where Dmax is the maximal
diffusion parameter, E 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:(3.4)0121201Furthermore we use the following Heun method
(third-order):(3.5)0131323023014034The implicit time-discretization is done with implicit
Euler or Runge-Kutta methods.Here, we use the implicit trapezoidal
rule:(3.6)0112121212Furthermore we use the following Gauss Runge-Kutta
method:(3.7)12−361414−3612+3614+36141212Remark 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)
discreteLmax-norm:(3.8)errLmax,Δx,Δt=maxi=1p|cnum(xi,T)−cref(xi,T)|,(ii)
discreteL1-norm: (3.9)errL1,Δx,Δt=∑i=1pΔx2|cnum(xi,T)−cref(xi,T)|,(iii)
discreteL2-norm: (3.10)errL2,Δx,Δt=∑i=1pΔx2|cnum(xi,T)−cref(xi,T)|2, where Δx is the
spatial-step of the discretization, Δt is the
time-step of the discretization, and T is the end-time
of the computation. p  is the number
of grid points in the discretization method. cnum   is the
numerical solution and cref  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(3.11)ρL2,Δx1,Δx2,Δt=log(errL2,Δx1,Δt/errL2,Δx2,Δt)log(Δx1/Δx2),where Δx1 is the coarse, Δx2 is the fine
spatial grid-step, and Δt is the
time-grid step for both results.(ii)
For the time error, we define(3.12)ρL2,Δx,Δt1,Δt2=log(errL2,Δx,Δt1/errL2,Δx,Δt2)log(Δt1/Δt2), where Δt1 is the coarse, Δt2 is the
fine-time-step, and Δx  is the
spatial-grid step for both results.We often useΔx2=Δx1/2. In this case, we have ρL2,Δx1,Δx2,Δt=ρL2,Δx1,Δt. Further we have to choose Δx1 with respect to
the Δt∈I=[0,Δtmax], which have maximal ρ. Thus we define ArgMax(Δx):(3.13)ArgMax(Δx):=argmaxΔt∈IρL2,Δx,Δx/2,Δt.


## 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Ωh the variational
boundary value problem reduces to find uh∈Vh satisfying the
initial condition uh(0) such
that(3.1)∫Ωh(∂uh∂tvh+D∇uh⋅∇vh)dx=0,∀vh∈Vh.We define the minimal length of
triangle which we get from the spatial-discretization with Δx.This leads to the following linear semidiscretized
system of ordinary differential equations:(3.2)Mdu*dt+Au*=0,where M is the mass and A the M-matrix.Here we have taken into account the
Courant-Friedrichs-Levy- (CFL-) condition, which is given as(3.3)CFL=2DmaxΔtmine∈EΔxe2,where Dmax is the maximal
diffusion parameter, E 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:(3.4)0121201Furthermore we use the following Heun method
(third-order):(3.5)0131323023014034The implicit time-discretization is done with implicit
Euler or Runge-Kutta methods.Here, we use the implicit trapezoidal
rule:(3.6)0112121212Furthermore we use the following Gauss Runge-Kutta
method:(3.7)12−361414−3612+3614+36141212Remark 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)
discreteLmax-norm:(3.8)errLmax,Δx,Δt=maxi=1p|cnum(xi,T)−cref(xi,T)|,(ii)
discreteL1-norm: (3.9)errL1,Δx,Δt=∑i=1pΔx2|cnum(xi,T)−cref(xi,T)|,(iii)
discreteL2-norm: (3.10)errL2,Δx,Δt=∑i=1pΔx2|cnum(xi,T)−cref(xi,T)|2, where Δx is the
spatial-step of the discretization, Δt is the
time-step of the discretization, and T is the end-time
of the computation. p  is the number
of grid points in the discretization method. cnum   is the
numerical solution and cref  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(3.11)ρL2,Δx1,Δx2,Δt=log(errL2,Δx1,Δt/errL2,Δx2,Δt)log(Δx1/Δx2),where Δx1 is the coarse, Δx2 is the fine
spatial grid-step, and Δt is the
time-grid step for both results.(ii)
For the time error, we define(3.12)ρL2,Δx,Δt1,Δt2=log(errL2,Δx,Δt1/errL2,Δx,Δt2)log(Δt1/Δt2), where Δt1 is the coarse, Δt2 is the
fine-time-step, and Δx  is the
spatial-grid step for both results.We often useΔx2=Δx1/2. In this case, we have ρL2,Δx1,Δx2,Δt=ρL2,Δx1,Δt. Further we have to choose Δx1 with respect to
the Δt∈I=[0,Δtmax], which have maximal ρ. Thus we define ArgMax(Δx):(3.13)ArgMax(Δx):=argmaxΔt∈IρL2,Δx,Δx/2,Δt.


## 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 Figure3, we present the P-controller.Figure 3
P-controller
for the solutionC.Our control problem is given with the control of the
error to the optimal concentration of the layer and correct the
source-flux:(4.1)∂tc−∇D∇c=csource,inΩ×[0,T],c(x,0)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where csource(x,t) is a
discontinuous source flow of the concentration c.We assume an optimal concentration at the layer with
the concentrationcopt(x,t), where the layer is given as x∈Ωlayer and 0 elsewhere.Our constraints are bounded as(4.2)csource,min≤csource≤csource,max.Remark 4.1.
Taken into account the hysteresis of
the deposition process, we apply a linear increase ofcsource,max in the optimal
control with respect to time, see Figure 4.Figure 4
Linear
constraintcopt for the
deposition process. x-axis: Time, y-axis: c (concentration).


### 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, whereP is the
proportional part, I the integral
part, and D is the
derivative part of the controller, see [16].These terms describe three basic mathematical
functions applied to the error signal,error=Coptimal−Ccomputed.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 multipliersKP, KI, and KD adding in
various amounts of these functions to get the system to behave the way you
want, see [11].Figure5 summarizes  the PID terms and their effect on a
control system.Figure 5
PID-control. Effect on control system: the main
influence in a control loop,KP reduces a large
part of the overall error. KI reduces the
final error in a system. Summing even a small error over time produces a drive
signal large enough to move the system toward a smaller error. KD counteracts the KP and KI terms when the
output changes quickly. This helps reduce overshoot and ringing.Initialization of the PID-Controller
The algorithm of the initialization of the PID-control (i.e.,
searchKP,KI,KD) is given as in
Algorithm 4.2  (see [15]).Algorithm 4.2.
(1)
We initialize the P-controller:KI=0.0,KD=0.0.
(2) The amplifying factorKP is increased
till we reached the permanent oscillations as a stability boundary of the
closed control system.
(3) We obtain forKP the critical
value KP,crit..
(4) The period-length of the permanent oscillation is given
asTcrit.
(5) We obtain the parameters from Table1.Table 1
Heuristic derivation of the control parameters
(Nichols-Ziegler).
ControllerKPTnTvP0.5KP,critPI0.45KP,critPID0.6KP,critFurther we compute the rest parameters asKI=KP/Tn,KD=KP*Tv, 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])(4.3)Δtn+1=(en−1en)KP(tolen)KI(en−12enen−2)KDΔtn,where tol is the
tolerance, en is the error of
the quantities of interest in time-step Δt.We can control the step-size with respect to our
heuristically computedKP, KI, and KD 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:KP,KD,KI 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,Csource our source and Cx the source
restricted to the response of the controlled system. To analyze the differences
between given source the Csource and the
computed source in the control process Cx, we study the step response from our
system:(4.4)xa(t):={Cx(t)∣CSource(t)=xe0}.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 xa, 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 responsexa, we identify the control path as a PT2 element.
Becauselimt→∞xa(t)=constant (also ≠0 nor ∞), the basic
behavior of our element is P (and neither D nor I). 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 bya0 a normed form
with three parameters ω0= angular
frequency, A = attenuation,
and KS = proportional
factor:(4.5)a2⋅x¨a(t)+a1⋅x˙a(t)+a0⋅xa(t)=b0⋅xe0,1ω02⋅x¨a(t)+2Aω0⋅x˙a(t)+xa(t)=KS⋅xe0.The solution of the differential
equation results in the characteristic
equation:(4.6)α2ω02+2A⋅αω0+1=0withα1,2=−A⋅ω0±ω0⋅A2−1.
So that we have the analytical solution for the normed
step response:(4.7)xa(t)xe0=KS⋅[1−T1T1−T2⋅e−t/T1+T2T1−T2⋅e−t/T2],withα1=−1T1,α2=−1T2.We go now the inverse way to
identify the system. We start with the step response. Using the inflected
tangent at (tw,xw,x˙a(tw)) and the limit xa(t→∞), we obtain our parameters.
Firstly, we haveKS=xa(t→∞)/xe0, tg=KS⋅(xe0/x˙a(tw)), and tu=tw−tg⋅(xa(tw)/xa(t→∞)), see Figure 6. In the second step, we obtain T1 and T2 by the
following formulas:(4.8)tutg=αα/(1−α)⋅(α⋅ln(α)+α2−1)α−1−1,tgT1=αα/(1−α),α=T2T1.Figure 6
Step response
and inflection tangent.Figure 7
Empirical determination of the model of the control
path. The three main cases: (1)xa(t=0)≠0:P,PDT1,PPT1,PIDT1, (2) xa(t=0)=0,x˙a(t=0)≠0:I,PT1,PDT2,  (3) xa(t=0),x˙a(t=0)=0,x¨a(t=0)≠0:IT1,PT2,….,  (1) xa(0)=0, (2) x˙a(0)=0, (3)  x¨a(0)≠0.Control Parameters
The parameter occurrence also in the transfer functionG(s), 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:(4.9)G(s)=KS(1+T1s)(1+T2s)PT2element,G(s)≈KS⋅e−Tus1+s⋅TgPT1element.Independent from the behavior which is
being wished, we obtain values for the parameters KP, KI, and KD, respectively. KR=KP, Tn and Tv 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 T. This is a general form of the Nichols-Ziegler method
(see Table 1). Further are KP=KR, KI=KP/Tn, and KD=KP*Tv.Table 2
Derivation of
the control parameters (Chien-Hrones-Reswick and Takahashi).
PID(PO)=0(PO)=20TakahashiKR0.6⋅tgtu⋅KS0.95⋅tgtu⋅KS1,2⋅tgKs⋅(tu+T)Tn1.0⋅tg1.35⋅tg2⋅(tu+T/2)2tu+TTv0.5⋅tu0.47⋅tu0.5(tu+T)Table 3
MATLAB-toolbox pdetool.
MATLAB functionDescriptionEnhancementParabolicSolve parabolic PDE (heat equation)In the future we plan to use
academical code to compute the PDE-solutions (e.g.,θ-method). Furthermore we will use
Comsol/FemlabPdetoolMATLAB toolbox to create the geometry
of the FEM-structure and the boundary conditionsWe have also implemented an alternative
mesh, with orthogonal triangles. MATLAB always uses equilateral trianglesRefineThis function refines the geometry of
the mesh. All triangles, were replaced with four new trianglesFor convergence rate calculations we
must guarantee that the geometry near the points, which are changed in the
backward step (source), are similarGuideMATLAB toolbox to create graphical user
interfaces (gui)


## 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 Figure3, we present the P-controller.Figure 3
P-controller
for the solutionC.Our control problem is given with the control of the
error to the optimal concentration of the layer and correct the
source-flux:(4.1)∂tc−∇D∇c=csource,inΩ×[0,T],c(x,0)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where csource(x,t) is a
discontinuous source flow of the concentration c.We assume an optimal concentration at the layer with
the concentrationcopt(x,t), where the layer is given as x∈Ωlayer and 0 elsewhere.Our constraints are bounded as(4.2)csource,min≤csource≤csource,max.Remark 4.1.
Taken into account the hysteresis of
the deposition process, we apply a linear increase ofcsource,max in the optimal
control with respect to time, see Figure 4.Figure 4
Linear
constraintcopt for the
deposition process. x-axis: Time, y-axis: c (concentration).


## 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, whereP is the
proportional part, I the integral
part, and D is the
derivative part of the controller, see [16].These terms describe three basic mathematical
functions applied to the error signal,error=Coptimal−Ccomputed.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 multipliersKP, KI, and KD adding in
various amounts of these functions to get the system to behave the way you
want, see [11].Figure5 summarizes  the PID terms and their effect on a
control system.Figure 5
PID-control. Effect on control system: the main
influence in a control loop,KP reduces a large
part of the overall error. KI reduces the
final error in a system. Summing even a small error over time produces a drive
signal large enough to move the system toward a smaller error. KD counteracts the KP and KI terms when the
output changes quickly. This helps reduce overshoot and ringing.Initialization of the PID-Controller
The algorithm of the initialization of the PID-control (i.e.,
searchKP,KI,KD) is given as in
Algorithm 4.2  (see [15]).Algorithm 4.2.
(1)
We initialize the P-controller:KI=0.0,KD=0.0.
(2) The amplifying factorKP is increased
till we reached the permanent oscillations as a stability boundary of the
closed control system.
(3) We obtain forKP the critical
value KP,crit..
(4) The period-length of the permanent oscillation is given
asTcrit.
(5) We obtain the parameters from Table1.Table 1
Heuristic derivation of the control parameters
(Nichols-Ziegler).
ControllerKPTnTvP0.5KP,critPI0.45KP,critPID0.6KP,critFurther we compute the rest parameters asKI=KP/Tn,KD=KP*Tv, 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])(4.3)Δtn+1=(en−1en)KP(tolen)KI(en−12enen−2)KDΔtn,where tol is the
tolerance, en is the error of
the quantities of interest in time-step Δt.We can control the step-size with respect to our
heuristically computedKP, KI, and KD 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:KP,KD,KI 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,Csource our source and Cx the source
restricted to the response of the controlled system. To analyze the differences
between given source the Csource and the
computed source in the control process Cx, we study the step response from our
system:(4.4)xa(t):={Cx(t)∣CSource(t)=xe0}.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 xa, 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 responsexa, we identify the control path as a PT2 element.
Becauselimt→∞xa(t)=constant (also ≠0 nor ∞), the basic
behavior of our element is P (and neither D nor I). 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 bya0 a normed form
with three parameters ω0= angular
frequency, A = attenuation,
and KS = proportional
factor:(4.5)a2⋅x¨a(t)+a1⋅x˙a(t)+a0⋅xa(t)=b0⋅xe0,1ω02⋅x¨a(t)+2Aω0⋅x˙a(t)+xa(t)=KS⋅xe0.The solution of the differential
equation results in the characteristic
equation:(4.6)α2ω02+2A⋅αω0+1=0withα1,2=−A⋅ω0±ω0⋅A2−1.
So that we have the analytical solution for the normed
step response:(4.7)xa(t)xe0=KS⋅[1−T1T1−T2⋅e−t/T1+T2T1−T2⋅e−t/T2],withα1=−1T1,α2=−1T2.We go now the inverse way to
identify the system. We start with the step response. Using the inflected
tangent at (tw,xw,x˙a(tw)) and the limit xa(t→∞), we obtain our parameters.
Firstly, we haveKS=xa(t→∞)/xe0, tg=KS⋅(xe0/x˙a(tw)), and tu=tw−tg⋅(xa(tw)/xa(t→∞)), see Figure 6. In the second step, we obtain T1 and T2 by the
following formulas:(4.8)tutg=αα/(1−α)⋅(α⋅ln(α)+α2−1)α−1−1,tgT1=αα/(1−α),α=T2T1.Figure 6
Step response
and inflection tangent.Figure 7
Empirical determination of the model of the control
path. The three main cases: (1)xa(t=0)≠0:P,PDT1,PPT1,PIDT1, (2) xa(t=0)=0,x˙a(t=0)≠0:I,PT1,PDT2,  (3) xa(t=0),x˙a(t=0)=0,x¨a(t=0)≠0:IT1,PT2,….,  (1) xa(0)=0, (2) x˙a(0)=0, (3)  x¨a(0)≠0.Control Parameters
The parameter occurrence also in the transfer functionG(s), 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:(4.9)G(s)=KS(1+T1s)(1+T2s)PT2element,G(s)≈KS⋅e−Tus1+s⋅TgPT1element.Independent from the behavior which is
being wished, we obtain values for the parameters KP, KI, and KD, respectively. KR=KP, Tn and Tv 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 T. This is a general form of the Nichols-Ziegler method
(see Table 1). Further are KP=KR, KI=KP/Tn, and KD=KP*Tv.Table 2
Derivation of
the control parameters (Chien-Hrones-Reswick and Takahashi).
PID(PO)=0(PO)=20TakahashiKR0.6⋅tgtu⋅KS0.95⋅tgtu⋅KS1,2⋅tgKs⋅(tu+T)Tn1.0⋅tg1.35⋅tg2⋅(tu+T/2)2tu+TTv0.5⋅tu0.47⋅tu0.5(tu+T)Table 3
MATLAB-toolbox pdetool.
MATLAB functionDescriptionEnhancementParabolicSolve parabolic PDE (heat equation)In the future we plan to use
academical code to compute the PDE-solutions (e.g.,θ-method). Furthermore we will use
Comsol/FemlabPdetoolMATLAB toolbox to create the geometry
of the FEM-structure and the boundary conditionsWe have also implemented an alternative
mesh, with orthogonal triangles. MATLAB always uses equilateral trianglesRefineThis function refines the geometry of
the mesh. All triangles, were replaced with four new trianglesFor convergence rate calculations we
must guarantee that the geometry near the points, which are changed in the
backward step (source), are similarGuideMATLAB toolbox to create graphical user
interfaces (gui)


## 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 withP1-elements and
an implicit Euler method for the time-discretization.DEPOSIT-PID Toolbox
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.Backward step.
(a)(b)Figure 9
Layout of the
DEPOSIT-PID Gui.Figure 10
Convergence diagram tool.In Section4 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.KONTOOL
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 Section3.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Δx and apply
all time discretization with
steps Δt, where the coarsest Δt
is given by the CFL condition
or till first non-numerical results
as oscillations. We compute
the error cnum−cref in the
L2-norm.
(3)  We continue the next fine
spatial steps, for example,Δx/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 side0. Our control problem has only the forward problem to
solve and we consider the accuracy of our simulations.We have the following
equation:(6.1)∂tc−β2(∂xx+∂yy)c=0,inΩ×[0,T],c(x,y,0)=cos(2x)+cos(2y),onΩ,∂nc(x,y,t)=0,on∂Ω×[0,T],where c is the molar
concentration, Ω=[−1,1]×[−1,1], and t∈(0,T). D=β2 is the
diffusion parameter of the diffusion equation.We have the following analytical
solution:(6.2)cana(x,y,t)=sin(2)+∑n=1∞Anexp(−β2n2π2t)(cos(nπx)+cos(nπy)),where An=(−1)n(−4sin(2)/(n2π2−4)). We apply the diffusion coefficient D=0.01, respectively, β=0.1. We obtain the following result after time T=12.0, see Figure 11.2D experiment of the diffusion equation at the end-timeT=12.0, β=0.1. x-axis: x-coordinate,
y-axis: y-coordinate and z-axis: c (concentration).
(a)(b)(c)(d)2D experiment of the diffusion equation and
control of a single point.
(a)(b)(c)(d)We see that for
largeΔx the numerical
solution converge faster to the stable constant endsolution (offset) than
solution for smaller Δx and the
analytical solution. The error between the constant analytical endsolution
(offset, sin(2)=0.9093) is also
greater than for smaller Δx. The L2-error is given
in Table 4.Table 4
Offset convergence (Δt=0.1, time-step=100, β=0.61644).
Δxoffset(analy.)max(num.)min(num.)L1-errorL2-error0.200.9092974268260.9021084076430.9021083812957.189e–35.17e–50.100.9092974268260.9075149940360.9075149870081.782e–33.18e–60.050.9092974268260.9088585115000.9088585096664.389e–41.93e–7Table 5
Numerical results for the P-controller for
different spatial-steps withD=0.1 and λ=1.0 as P-value
for the controller.
ΔxΔterr=|unum,Δx,Δt−unum,fineΔx,Δt|Convergence rate0.10.10.0770073.25310.10.050.0161534.00770.10.0250.00200853.94470.10.01250.000260871.82760.10.006250.0001469900.050.10.278732.85910.050.050.0768333.74810.050.0250.0114374.0520.050.01250.0013793.77760.050.006250.0002011100.0250.10.65642.45520.0250.050.239393.24490.0250.0250.0505054.0080.0250.01250.00627813.99990.0250.006250.0007848200.01250.11.042.14240.01250.050.47112.71730.01250.0250.143273.72510.01250.01250.0216684.05160.01250.006250.002613300.006250.11.20921.5860.006250.050.805542.41790.006250.0250.301483.21790.006250.01250.0648034.0050.006250.006250.00807260Table 6
The associatedΔt, which have maximal convergence rate  (ArgMax), for Δx=0.1, 0.05, 0.025, 0.0125.
ΔxArgMax−log(Max).10000.650000.62.05000.190002.4.025000.044.64.01250.010006.64Table 7
The associatedΔt, which have maximal convergence rate  (ArgMax), for Δx=0.1, 0.05, 0.025, 0.0125.
ΔxArgMax−log(Max).100000.07003.84.050000.01705.88.025000.00407.97.012500.00109.97Table 8
Convergence ofλ with
bisection.
λa(max)a(min)Time-step790.0000.0312676128758890.0312622187118722592.5000.0179287914475040.0180126082683413493.7500.0097555203576400.0096631156512394444.3750.0050116366334120.0050070357573585419.6870.0024642094388410.0025049709735506407.3430.0011316284141310.0011825337734157401.1710.0004669074552840.0004657945761918398.0850.0001092600682670.0001896372556049397.3140.0000803123993140.00004320189361311396.9280.0000416282413610.00001211849267412396.7350.0000029233842930.00000604731877613396.72918396.723−0.000017536337301−0.00000044452143517396.711−0.000016718130883−0.00000491279461716396.687−0.000002204413679−0.00002227953598115396.639−0.000015370659931−0.00001674376709314396.542−0.000042906107641−0.00001511764100210395.000−0.000163960451645−0.0001649989271891Remark 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:(6.3)∂tc−∇D∇c=f(t),inΩ×[0,T],c(x,0)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where c is the molar
concentration, D is the
diffusion parameter of the diffusion equation, and f(t) is the
right-hand side or diffusion source. We have the following constraint: coptimal(xpoint,ypoint)=0.5, where (xpoint,ypoint) is the control
point in our domain. The parameters are given as D=0.01 and f(t)=at+b, so we deal with a linear source, a=0.2 and b=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 Figure12.


### 6.3. Preliminary Remark for Simulations for the Convergence Order

To determine the functionArgMax, firstly, we have to determine a practicable interval
I for(6.4)ArgMax(Δx):=argmaxΔt∈IρL2,Δx,Δx/2,Δt.One possibility is the
interval(6.5)ICFL:=(0,ΔtCLF],with ΔtCFL:=Δs2/2Dmax from the
CFL-condition. In the experiments, we find another interval, where the convergence
rate function is relative stable convex:(6.6)Istable:=(0,Δtstable],where(6.7)Δtstable:=maxΔt{ρ:(0,Δt)→ℝisconvex}.It is clear that we can in Istable take some
restrictions to a subinterval Isub:=[Δtsmin,Δtsmax]⊂Istable, when we know that some Δt∈Isub with ρ(Δt)>ρ(Δtsmin),ρ(Δtsmax).The function of the numerical convergence rate is
discrete in the spatial-discretization variableΔx since we get a
finer discretization with a bisection of Δx. With a finer discretization, a triangle is replaced by
four subtriangles.In the temporal discretization variableΔt, we are in contrast not restricted to such conditions
and Δx 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 Δx and Δt. This function also depends on some parameters, for
instance, the diffusion coefficient D and the control
parameter λ. We now present the results of two chosen
experiments. For the first experiment, the parameters are D=0.1 and λ=1, while we increase the propagation velocity in such a
way that we increase D to 1 for the second
experiment.We observe that, for instance, in the case of
spatial-discretizationΔx=0.05 and x=0.25, the maximal
convergence rate lies between Δt=0.05 and Δt=0.0125. The maximum itself lies close to 0.025. To determine the precise value of Δt=ArgMax(Δx), we refine our method in these interesting intervals.In
Figure13, 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 (ΔtCFL<Δtstable). For Δx=0.05,D=0.1, the CFL-area ends at ΔtCFL=0.0125, the stable convex area range to about Δtstable=0.2. For our convergence diagram, we try to find for
every Δx the value Δt, where the convergence rate becomes maximal (ArgMax).
It is clear that we only use Δt at which the
convergence rate function is stable: Δt<Δtstable.ρ for a
P-controller, D=0.1.
(a)(b)Figure 14
Convergence diagram (KONTOOL). We can see in the loglog-plotArgMax(Δx) linear
dependence.Figure 15
Convergence diagram (KONTOOL). The loglog-plotArgMax(Δx)
with linear dependence.ρ for a
P-controller, D=1.
(a)(b)(c)(d)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:(6.8)∂tc−∇D∇c=f(t),inΩ×[0,T],c(x,t)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where c is the molar
concentration, D is the
diffusion parameter of the diffusion equation, and f(t) is the right
hand side or source.We have the following constraint:(6.9)coptimal(xpoint,ypoint)=0.5, where (xpoint,ypoint)=(0,−1) is the control
point in our domain.The
automatically step-size is given as (see [11])(6.10)Δtn+1=(en−1en)KP(tolen)KI(en−12enen−2)KDΔtn,where tol=1 is the
tolerance, en is the error of
the quantities of interest in time-step Δt.The errors are given as(6.11)en=||un−un−1||||un||,where un is the result
at time-step tn.The parameters are given asD=0.1, Pkrit=15, Tkrit=5, Δt (PID-control), tol=1 (adaptive
PID-control), Δtmax=0.1 (adaptive
PID-control), Δtmin=0.01 (adaptive
PID-control) (see Figure 17).Figure 17
2D experiment with and without the adaptive time-step control.Furthermore, we change the parametertol=1 to tol=0.1 and Δtmin=0.01 to Δtmin=0.0001 (see Figure 18).Adaptive PID
with modified parameters.
(a)(b)Remark 6.3.
In  Figures17 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 KP, KI, and KD before the
step-size control for the whole time-interval 0,T. The optimal control parameters are given to the
whole time-interval 0,T. A modified algorithm to compute the control
parameters for the initialization time-interval 0,tbegin and the whole
time-interval 0,T 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 parametersKP,global, KI,global and KD,global for the
time-interval 0,T.
(2) We apply the automatically step-size control for
the global control parameters withtolglobal, Δtmax,global and Δtmin,global, which can by chosen large.
(3) We stop the computation till we reach the optimal
solution and mark remember the timetoszill.
(4) We compute the local control parametersKP,local, KI,local and KD,local for the
time-interval 0,toszill.
(5) We restart the computation with the local control
parameters and smaller step-size parameterstollocal, Δtmax,local and Δtmin,local till we reach toszill and continue
the computation with the global parameters.
(6) We stop the computation if we reacht=T. 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
parametersKP, KI, and KD 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 Algorithm6.4.To be more precise the computation of global control
parametersKP,global=KP(Δt), KI,global=KI(Δt), and KD,global=KD(Δt) can be
automatized in an interval Δt∈(0,T).We improve the automatically step-size, which is given
in Section4.3 to(6.12)(Δt)n+1=(en−1en)KP(Δt)(tolen)KI(Δt)(en−12enen−2)KD(Δt)(Δt)n,where tol is the
tolerance, en is the error of
the quantities of interest in time-step n.With approximations to each subintervals of(0,T), we derive the control parameters KP(Δt), KI(Δt), and KP(Δt).The improved globalized control parameters are
computed in the following algorithm.Algorithm 6.6.
Computing of the control parametersKP,global=KP(t), KI,global=KI(t) and KD,global=KD(t) for the
time-interval 0,T.
(1) Compute the critical control function with the
minima and maxima. Based on the optimal oscillated function, we can derive theKP(T/3) for the
interval (0,T/3).
(2) Redo the step (1) for the intervals(T/3,2T/3) and 2T/3,T and approximate
the function KP(t). Based on the KP(t) we can derive
the KI(t) and KD(t) function.
(3) The optimal control parameters are used in
Algorithm6.4.Remark 6.7.
Further ideas to determine the time-dependent parametersKP(Δt), KI(Δt), and KD(Δt) 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
P-value.


### 6.7. Determinate of Critical P with Exponential Regression

To initialize the PID-controller, we have to derive
theKP, KD, and KI 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 Figure19, we can see that the
maximum and minimum of
the error function, done with the step function response, have a exponential
behavior.Figure 19
Error of the
time-step control.Exponential
regression of min and max. Behavior ofλ.
(a)(b)So we use the exponential regression(6.13)y=b⋅exp(a⋅t)+ϵ,with ϵ~N(0,σ2) normal
distributed. Via log the exponential
regression can be formed in a linear regression:(6.14)y=b⋅exp(a⋅t),log(y)=log(b)+a⋅t.Before we have to convert the
minima into the positive area.Letλ(n) be the actual
lambda value, λ(n)+ the next
greater λ, and λ(n)− the next
smaller one. Firstly, we use  (6.15)λ(n+1)={(λ(n)+λ(n)+)2:a(max,n)<0,a(min,n)>0(λtosmall),(λ(n)+λ(n)−)2:a(max,n)>0,a(min,n)<0(λtogreat). In the first step and in limit case, we use λ(n+1)=2⋅λ(n), respectively, λ(n+1)=λ(n)/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 a-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 parametersKP,global, KI,global, and KD,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 theKP, KI, and KD parameters of
this process. Here, we present the experiments in Figure 21, which
had taken into account the critical values Pcrit and Tcrit to derive the
optimal KP, KD, and KI parameters
based on Algorithm 6.6.Δt dependence
parameters determined by Nichols-Ziegler.
(a)(b)(c)Δt  dependence
parameters determined by Chien-Hrones-Reswick with PO=0,…,20.
(a)(b)(c)A further algorithm to derive the control parameters
is done by Chien et al., see [12].The results are given in Figure23.Δt dependence
parameters determined by Takahashi.
(a)(b)(c)TheKP,KD, and KI 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 theKP, KI, and KD parameters of
this process. The minimization problem is presented in Figure 24, which
had taken into account the critical parameters.Parameters
chosen by Chien-Hrones-Reswick: Summands of the minimization problem.
(a)(b)(c)(d)In Figure25, 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.Figure 25
Parameters chosen by Chien-Hrones-Reswick: Actual
value.In Figure26, 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.Figure 26
Parameters
chosen by Chien-Hrones-Reswick: Source.Remark 6.9.
The benefit of the time-dependent
control parametersKP(Δt), KD(Δt), and KI(Δt) 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.


## 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 side0. Our control problem has only the forward problem to
solve and we consider the accuracy of our simulations.We have the following
equation:(6.1)∂tc−β2(∂xx+∂yy)c=0,inΩ×[0,T],c(x,y,0)=cos(2x)+cos(2y),onΩ,∂nc(x,y,t)=0,on∂Ω×[0,T],where c is the molar
concentration, Ω=[−1,1]×[−1,1], and t∈(0,T). D=β2 is the
diffusion parameter of the diffusion equation.We have the following analytical
solution:(6.2)cana(x,y,t)=sin(2)+∑n=1∞Anexp(−β2n2π2t)(cos(nπx)+cos(nπy)),where An=(−1)n(−4sin(2)/(n2π2−4)). We apply the diffusion coefficient D=0.01, respectively, β=0.1. We obtain the following result after time T=12.0, see Figure 11.2D experiment of the diffusion equation at the end-timeT=12.0, β=0.1. x-axis: x-coordinate,
y-axis: y-coordinate and z-axis: c (concentration).
(a)(b)(c)(d)2D experiment of the diffusion equation and
control of a single point.
(a)(b)(c)(d)We see that for
largeΔx the numerical
solution converge faster to the stable constant endsolution (offset) than
solution for smaller Δx and the
analytical solution. The error between the constant analytical endsolution
(offset, sin(2)=0.9093) is also
greater than for smaller Δx. The L2-error is given
in Table 4.Table 4
Offset convergence (Δt=0.1, time-step=100, β=0.61644).
Δxoffset(analy.)max(num.)min(num.)L1-errorL2-error0.200.9092974268260.9021084076430.9021083812957.189e–35.17e–50.100.9092974268260.9075149940360.9075149870081.782e–33.18e–60.050.9092974268260.9088585115000.9088585096664.389e–41.93e–7Table 5
Numerical results for the P-controller for
different spatial-steps withD=0.1 and λ=1.0 as P-value
for the controller.
ΔxΔterr=|unum,Δx,Δt−unum,fineΔx,Δt|Convergence rate0.10.10.0770073.25310.10.050.0161534.00770.10.0250.00200853.94470.10.01250.000260871.82760.10.006250.0001469900.050.10.278732.85910.050.050.0768333.74810.050.0250.0114374.0520.050.01250.0013793.77760.050.006250.0002011100.0250.10.65642.45520.0250.050.239393.24490.0250.0250.0505054.0080.0250.01250.00627813.99990.0250.006250.0007848200.01250.11.042.14240.01250.050.47112.71730.01250.0250.143273.72510.01250.01250.0216684.05160.01250.006250.002613300.006250.11.20921.5860.006250.050.805542.41790.006250.0250.301483.21790.006250.01250.0648034.0050.006250.006250.00807260Table 6
The associatedΔt, which have maximal convergence rate  (ArgMax), for Δx=0.1, 0.05, 0.025, 0.0125.
ΔxArgMax−log(Max).10000.650000.62.05000.190002.4.025000.044.64.01250.010006.64Table 7
The associatedΔt, which have maximal convergence rate  (ArgMax), for Δx=0.1, 0.05, 0.025, 0.0125.
ΔxArgMax−log(Max).100000.07003.84.050000.01705.88.025000.00407.97.012500.00109.97Table 8
Convergence ofλ with
bisection.
λa(max)a(min)Time-step790.0000.0312676128758890.0312622187118722592.5000.0179287914475040.0180126082683413493.7500.0097555203576400.0096631156512394444.3750.0050116366334120.0050070357573585419.6870.0024642094388410.0025049709735506407.3430.0011316284141310.0011825337734157401.1710.0004669074552840.0004657945761918398.0850.0001092600682670.0001896372556049397.3140.0000803123993140.00004320189361311396.9280.0000416282413610.00001211849267412396.7350.0000029233842930.00000604731877613396.72918396.723−0.000017536337301−0.00000044452143517396.711−0.000016718130883−0.00000491279461716396.687−0.000002204413679−0.00002227953598115396.639−0.000015370659931−0.00001674376709314396.542−0.000042906107641−0.00001511764100210395.000−0.000163960451645−0.0001649989271891Remark 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:(6.3)∂tc−∇D∇c=f(t),inΩ×[0,T],c(x,0)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where c is the molar
concentration, D is the
diffusion parameter of the diffusion equation, and f(t) is the
right-hand side or diffusion source. We have the following constraint: coptimal(xpoint,ypoint)=0.5, where (xpoint,ypoint) is the control
point in our domain. The parameters are given as D=0.01 and f(t)=at+b, so we deal with a linear source, a=0.2 and b=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 Figure12.


## 6.3. Preliminary Remark for Simulations for the Convergence Order

To determine the functionArgMax, firstly, we have to determine a practicable interval
I for(6.4)ArgMax(Δx):=argmaxΔt∈IρL2,Δx,Δx/2,Δt.One possibility is the
interval(6.5)ICFL:=(0,ΔtCLF],with ΔtCFL:=Δs2/2Dmax from the
CFL-condition. In the experiments, we find another interval, where the convergence
rate function is relative stable convex:(6.6)Istable:=(0,Δtstable],where(6.7)Δtstable:=maxΔt{ρ:(0,Δt)→ℝisconvex}.It is clear that we can in Istable take some
restrictions to a subinterval Isub:=[Δtsmin,Δtsmax]⊂Istable, when we know that some Δt∈Isub with ρ(Δt)>ρ(Δtsmin),ρ(Δtsmax).The function of the numerical convergence rate is
discrete in the spatial-discretization variableΔx since we get a
finer discretization with a bisection of Δx. With a finer discretization, a triangle is replaced by
four subtriangles.In the temporal discretization variableΔt, we are in contrast not restricted to such conditions
and Δx 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 Δx and Δt. This function also depends on some parameters, for
instance, the diffusion coefficient D and the control
parameter λ. We now present the results of two chosen
experiments. For the first experiment, the parameters are D=0.1 and λ=1, while we increase the propagation velocity in such a
way that we increase D to 1 for the second
experiment.We observe that, for instance, in the case of
spatial-discretizationΔx=0.05 and x=0.25, the maximal
convergence rate lies between Δt=0.05 and Δt=0.0125. The maximum itself lies close to 0.025. To determine the precise value of Δt=ArgMax(Δx), we refine our method in these interesting intervals.In
Figure13, 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 (ΔtCFL<Δtstable). For Δx=0.05,D=0.1, the CFL-area ends at ΔtCFL=0.0125, the stable convex area range to about Δtstable=0.2. For our convergence diagram, we try to find for
every Δx the value Δt, where the convergence rate becomes maximal (ArgMax).
It is clear that we only use Δt at which the
convergence rate function is stable: Δt<Δtstable.ρ for a
P-controller, D=0.1.
(a)(b)Figure 14
Convergence diagram (KONTOOL). We can see in the loglog-plotArgMax(Δx) linear
dependence.Figure 15
Convergence diagram (KONTOOL). The loglog-plotArgMax(Δx)
with linear dependence.ρ for a
P-controller, D=1.
(a)(b)(c)(d)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:(6.8)∂tc−∇D∇c=f(t),inΩ×[0,T],c(x,t)=c0(x),onΩ,∂c(x,t)∂n=c1(x,t),on∂Ω×[0,T],where c is the molar
concentration, D is the
diffusion parameter of the diffusion equation, and f(t) is the right
hand side or source.We have the following constraint:(6.9)coptimal(xpoint,ypoint)=0.5, where (xpoint,ypoint)=(0,−1) is the control
point in our domain.The
automatically step-size is given as (see [11])(6.10)Δtn+1=(en−1en)KP(tolen)KI(en−12enen−2)KDΔtn,where tol=1 is the
tolerance, en is the error of
the quantities of interest in time-step Δt.The errors are given as(6.11)en=||un−un−1||||un||,where un is the result
at time-step tn.The parameters are given asD=0.1, Pkrit=15, Tkrit=5, Δt (PID-control), tol=1 (adaptive
PID-control), Δtmax=0.1 (adaptive
PID-control), Δtmin=0.01 (adaptive
PID-control) (see Figure 17).Figure 17
2D experiment with and without the adaptive time-step control.Furthermore, we change the parametertol=1 to tol=0.1 and Δtmin=0.01 to Δtmin=0.0001 (see Figure 18).Adaptive PID
with modified parameters.
(a)(b)Remark 6.3.
In  Figures17 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 KP, KI, and KD before the
step-size control for the whole time-interval 0,T. The optimal control parameters are given to the
whole time-interval 0,T. A modified algorithm to compute the control
parameters for the initialization time-interval 0,tbegin and the whole
time-interval 0,T 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 parametersKP,global, KI,global and KD,global for the
time-interval 0,T.
(2) We apply the automatically step-size control for
the global control parameters withtolglobal, Δtmax,global and Δtmin,global, which can by chosen large.
(3) We stop the computation till we reach the optimal
solution and mark remember the timetoszill.
(4) We compute the local control parametersKP,local, KI,local and KD,local for the
time-interval 0,toszill.
(5) We restart the computation with the local control
parameters and smaller step-size parameterstollocal, Δtmax,local and Δtmin,local till we reach toszill and continue
the computation with the global parameters.
(6) We stop the computation if we reacht=T. 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
parametersKP, KI, and KD 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 Algorithm6.4.To be more precise the computation of global control
parametersKP,global=KP(Δt), KI,global=KI(Δt), and KD,global=KD(Δt) can be
automatized in an interval Δt∈(0,T).We improve the automatically step-size, which is given
in Section4.3 to(6.12)(Δt)n+1=(en−1en)KP(Δt)(tolen)KI(Δt)(en−12enen−2)KD(Δt)(Δt)n,where tol is the
tolerance, en is the error of
the quantities of interest in time-step n.With approximations to each subintervals of(0,T), we derive the control parameters KP(Δt), KI(Δt), and KP(Δt).The improved globalized control parameters are
computed in the following algorithm.Algorithm 6.6.
Computing of the control parametersKP,global=KP(t), KI,global=KI(t) and KD,global=KD(t) for the
time-interval 0,T.
(1) Compute the critical control function with the
minima and maxima. Based on the optimal oscillated function, we can derive theKP(T/3) for the
interval (0,T/3).
(2) Redo the step (1) for the intervals(T/3,2T/3) and 2T/3,T and approximate
the function KP(t). Based on the KP(t) we can derive
the KI(t) and KD(t) function.
(3) The optimal control parameters are used in
Algorithm6.4.Remark 6.7.
Further ideas to determine the time-dependent parametersKP(Δt), KI(Δt), and KD(Δt) 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
P-value.


## 6.7. Determinate of Critical P with Exponential Regression

To initialize the PID-controller, we have to derive
theKP, KD, and KI 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 Figure19, we can see that the
maximum and minimum of
the error function, done with the step function response, have a exponential
behavior.Figure 19
Error of the
time-step control.Exponential
regression of min and max. Behavior ofλ.
(a)(b)So we use the exponential regression(6.13)y=b⋅exp(a⋅t)+ϵ,with ϵ~N(0,σ2) normal
distributed. Via log the exponential
regression can be formed in a linear regression:(6.14)y=b⋅exp(a⋅t),log(y)=log(b)+a⋅t.Before we have to convert the
minima into the positive area.Letλ(n) be the actual
lambda value, λ(n)+ the next
greater λ, and λ(n)− the next
smaller one. Firstly, we use  (6.15)λ(n+1)={(λ(n)+λ(n)+)2:a(max,n)<0,a(min,n)>0(λtosmall),(λ(n)+λ(n)−)2:a(max,n)>0,a(min,n)<0(λtogreat). In the first step and in limit case, we use λ(n+1)=2⋅λ(n), respectively, λ(n+1)=λ(n)/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 a-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 parametersKP,global, KI,global, and KD,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 theKP, KI, and KD parameters of
this process. Here, we present the experiments in Figure 21, which
had taken into account the critical values Pcrit and Tcrit to derive the
optimal KP, KD, and KI parameters
based on Algorithm 6.6.Δt dependence
parameters determined by Nichols-Ziegler.
(a)(b)(c)Δt  dependence
parameters determined by Chien-Hrones-Reswick with PO=0,…,20.
(a)(b)(c)A further algorithm to derive the control parameters
is done by Chien et al., see [12].The results are given in Figure23.Δt dependence
parameters determined by Takahashi.
(a)(b)(c)TheKP,KD, and KI 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 theKP, KI, and KD parameters of
this process. The minimization problem is presented in Figure 24, which
had taken into account the critical parameters.Parameters
chosen by Chien-Hrones-Reswick: Summands of the minimization problem.
(a)(b)(c)(d)In Figure25, 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.Figure 25
Parameters chosen by Chien-Hrones-Reswick: Actual
value.In Figure26, 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.Figure 26
Parameters
chosen by Chien-Hrones-Reswick: Source.Remark 6.9.
The benefit of the time-dependent
control parametersKP(Δt), KD(Δt), and KI(Δt) 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.



---

*Source: 728105-2009-05-10.xml*