A Novel Method for Control Performance Assessment with Fractional Order Signal Processing and Its Application to Semiconductor Manufacturing

The significant task for control performance assessment (CPA) is to review and evaluate the performance of the control system. The control system in the semiconductor industry exhibits a complex dynamic behavior, which is hard to analyze. This paper investigates the interesting crossover properties of Hurst exponent estimations and proposes a novel method for feature extraction of the nonlinear multi-input multi-output (MIMO) systems. At first, coupled data from real industry are analyzed by multifractal detrended fluctuation analysis (MFDFA) and the resultant multifractal spectrum is obtained. Secondly, the crossover points with spline fit in the scale-law curve are located and then employed to segment the entire scale-law curve into several different scaling regions, in which a single Hurst exponent can be estimated. Thirdly, to further ascertain the origin of the multifractality of control signals, the generalized Hurst exponents of the original series are compared with shuffled data. At last, non-Gaussian statistical properties, multifractal properties and Hurst exponents of the process control variables are derived and compared with different sets of tuning parameters. The results have shown that CPA of the MIMO system can be better employed with the help of fractional order signal processing (FOSP).


Introduction
Control loops are the most significant elements in automation systems.Product quality, operation safety, material and energy consumption, and thus the financial performance, are directly or indirectly linked to the performance of control systems.Control Performance Assessment (CPA) provides a benchmark for detecting and diagnosing root-causes of poor performance, reviewing and evaluating current performance, as well as improving control performance or avoiding performance degradation [1].
However, the mainstream assumption of Gaussian character of process signals contributes to the basis for CPA analysis algorithms and methods.On the contrary, a typical industrial process, which includes thousands of control loops, has witnessed many non-Gaussian behaviors, such as long range dependency (LRD), self-similarity, power-law of autocorrelation, infinite variance, spiky signals.Non-Gaussian signals and noises tend to produce large-amplitude fluctuations from the average value more frequently than Gaussian ones do [2].Process control systems are very complex, usually encompassing different hierarchical levels, therefore it is inconceivable for plant personals to maintain them on regular basis.For example, control loops are often multivariable and exhibit nonlinear dynamics stemming either from the plant, the transducers, the actuators, or even in some cases the controllers themselves in industrial applications.
To address the above issues, the fractional order signal processing (FOSP) techniques have been proposed to better characterize the control process in the recent years with the notion of fractional calculus [2].FOSP techniques include fractional order linear systems, autoregressive fractional integrated moving average (ARFIMA), Hurst parameter estimation, fractional order Fourier transformation (FrFT), fractals, Multifractal detrended fluctuation analysis (MFDFA), etc [3].In addition, fractional order thinking is inevitable to gain more insights to characterize complex objects [4].
CPA is an important asset-management technology to maintain highly efficient operation performance of automation systems in production plants [5].There are many classic performance assessment approaches, such as mean squared error (MSE), integral absolute error (IAE), statistical indexes, fractal indexes, etc.In addition, they can be various depending on systems under the different conditions, for example, with control action constraints, with deterministic disturbances, with setpoint changes.How to evaluate the quality of the control system performance is an essential issue for process engineers.In recent years, multifractal analysis has become the focus among the control engineering after some researchers gave applied it with CPA.The on-line control loop performance monitoring with non-Gaussian statistical and fractal measure has been proposed in [6].Doma ński presents results of the research on alternative CPA measures applied to control quality assessment for SISO loop with generalized predictive control (GPC) controller in [7,8].
However, as mentioned in [8], the above CPA analyses and observations are accomplished with a simple linear SISO case, which cannot capture the monofractal and multifractal properties, since the process complexity, cross-dependencies with varying delays, LRD and human factors are not reflected in simulations.Many intriguing questions are left behind in the above articles, such as multiple input multiple output (MIMO) and systems with significant delays.Moreover, the origins and meaning of multifractal properties with crossover phenomena are still needed to be analyzed.Therefore, it is worth investigating more complex scenarios, such as nonlinear, MIMO and systems with significant delay cases to verify method applicability and effectiveness.The aim of this paper is to cast more light on the new directions for process engineers and provide some novel techniques to assess the process performance.
More specifically, the purposes of this paper are to: 1.
Use MFDFA to analyze semiconductor data, derive the multifractal spectrum and select the characteristic parameters sensitive to changes of the control system.

2.
Extract the characteristic parameters from the multifractal spectrum of reference data to form reference feature sets; 3.
Modify the standard single Hurst exponent estimation by the multiple Hurst exponent fitting method with crossover points; 4.
Select multifractal properties and modified Hurst exponents to distinguish different types of control actions (tunings).
Section 1 starts with a brief introduction of the current FOSP techniques in CPA.Section 2 provides the definition of the fractional Gaussian noise, Hurst exponent and α-stable distribution for readers with zero knowledge.Section 3 proposes the MFDFA algorithm and Hurst spline fit with crossovers.Data analysis with MFDFA for each loop of the real MIMO system is carried in Section 4. The analysis of the results are carried in Section 5.The conclusion is given in Section 6.

Fractal and Fractional Gaussian Noise
Fractal analysis is especially useful when the data show self-similarity, power-law, scale-invariant and nonlinear properties [9][10][11].From the fractal theory introduced by Mandelbrot in [9], the traditional concept of three-dimensional space can be extended to the fractal (fractional) dimension (FD).
Based on the definitions of fractional-order differential operators, many complex dynamic systems with complex memory behaviors can be more properly described by the fractional calculus.Some researchers found interesting the analytical results of the linear fractional-order differential equations are represented by the Mittag-Leffler function, which exhibits a power-law asymptotic behavior [12].Therefore, the fractional calculus is being widely used to analyze the random signals with power-law size distributions or power-law decay of correlations [13,14].
The fractional order integral of the function f (x) with α ∈ R+ is defined where Γ(•) is the Gamma function, a D −α t is the fractional integral of order α in [a, t] [14].Based on fractional order integral, the relationship of white Gaussian noise wGn ω(t), fractional Gaussian noise (fGn), α-stable noise ω α (t), fractional stable noise, multifractional Gaussian noise and multifractional stable noise can be described in Figure 1.Besides, the fGn can be expressed as the αth order integration of wGn [2,15]: Therefore, from the perspective of fractional signals and fractional-order systems the fGn can be simulated by the dth integrator with wGn as the input in the Figure 1.From the above figure, the major difference of monofractal and multifracal property is the time-depend variable H.The linear multifractional stable motion is obtained by generalizing the constant Hurst parameter H to a time-dependent local Hölder exponent H(t).

Hurst Parameter
LRD is characterized by the Hurst exponent H.It means that there is a strong coupling between values at different times.The coupling effect characterized by the autocorrelations functions (ACF) and power spectrum density (PSD) obey the power law, hyperbolically decaying, while the conventional assumptions are based on the exponentially decaying effects.To quantify the level of the coupling, the Hurst exponent and many valuable Hurst exponent estimators have been provided to more accurately characterize the LRD time series.
The rescaled range (R/S) method is one of the time-domain analysis of Hurst parameter defined as follows [16]: where E(•) denotes the expected value of the observations, R(n) is the range of the first n values, S(n) is their standard deviation, and C is a constant.Alternatively, Whittle's maximum likelihood estimator (MLE) and wavelet analysis are using periodogram based analysis in the frequency domain [17].LRD processes appear in many contexts, which is characterized by the Hurst parameter H (0 < H < 1).The phenomenon of LRD can be observed in hydrology, finance, economics and so on.Unlike with a stationary process, the autocorrelations between observations of a long memory series are slowly decaying to zero.However, according to the survey paper [18], Hurst estimators can be significantly affected by trend and seasonality.
LRD or long memory process can be defined by the autocovariance γ(k) in the time domain, or defined by the power spectrum P( f ) in the frequency domain.The power law is observed in the long memory ARFIMA(p, d, q) process {X t } with d ∈ (−1/2, 1/2) \ {0}, since the asymptotic behaviour of the autocovariance function γ(•) is given by [19] In particular, if 0 < d < 0.5, then {X t } is a long-memory process, or long range positive dependence, since . The process is said to exhibit intermediate memory (anti-persistence) (Also known as a mean-reverting process), or long-range negative dependence, for −0.5 < d < 0. In addition, the process is non-stationary for |d| ≥ 0.5, as it possesses infinite variance and first order difference is needed to get the stationary series.Therefore, the fractional order difference filter is the first step of processing long memory ARFIMA models [21].
Besides, the fractional order difference (integral) coefficient d has a closed relation with Hurst parameter and characteristic exponent α: where H is the Hurst parameter and α is from the α-stable distribution, which will be introduced in the next subsection.The Hurst parameter of fGn is related to d by the Equation (5).

α Stable Distribution
In the real system, non-stationary, non-Gaussian, spiky signals are usually regarded as outliers and thus discarded by engineers during the signal processing.In fact, however, each point of data from the real system could tell a "story" such as malfunctioning, inappropriate control philosophy, abnormal human intervention, etc. Non-Gaussian signals and noises tend to produce large-amplitude fluctuations from the average value more frequently than Gaussian ones do.Where did non-Gaussian signals come from?Many researchers have found that they are from fractional calculus or even more specifically, α-stable processes [22].The α stable distribution, which is the generalization of the Gaussian distribution, is used to model the non-Gaussian distribution.The α-stable distribution based techniques have been applied to describe many natural or man-made phenomena in various fields, such as physics, hydrology, biology, financial and network traffic.Stable distributions provide a useful theoretical tool for this type of signals and noises [23].The α-stable characteristic function (or distribution) is determined by four parameters: α, β, γ and δ.A univariate distribution function F(x) is stable if and only if its characteristic function has the form: where and α is called the characteristic exponent.A small value of α will imply considerable probability mass in the tails of the distribution.That is, the smaller α is, the heavier are the tails.α = 2 corresponds to the Gaussian distribution (for any β).γ is a scaling parameter called the dispersion.It is similar to the variance of the Gaussian distribution.β is a symmetry parameter.β = 0 indicates a distribution symmetric about δ.
There are two algorithms provide consistent estimators for all four parameters-Koutrouvelis's method, which is based on empirical characteristic function methods of Koutrouvelis [24] and McCulloch's method, which is a simpler, faster but less accurate method [25].In this paper, we apply Koutrouvelis method for estimation of these parameters in Section 4.

MFDFA Algorithm
MFDFA is an efficient FOSP technique to detect multifractality in a time series since it does not require any knowledge about the process time delay or other process parameters [26].MFDFA can estimate the multifractal spectrum of the generalized Hurst exponent from a time series and it does not require any knowledge about the process time delay or other process parameters [27].Currently, MFDFA has been successfully applied to analyze various data, such as hydrographic data [28], wind records [29], financial time series [30], traffic time series [31], mechanical vibration signals [32], etc.It has proven to be a powerful tool for uncovering the multifractality of non-stationary time series in the complex systems.In this study, MFDFA is applied to detect the presence of multifractal and monofratal properties in the MIMO systems with external disturbance and noise.Then, the origin of the multifractality of control signals, the generalized Hurst exponents of the original series are compared with shuffled data.

Basic MFDFA Algorithm
Detrended fluctuation analysis (DFA) was first proposed by Peng et al. on the DNA analysis in 1995 [33] and then MFDFA was proposed by Kantelhardt et al. in 2002 [34].For the calculation methodology, Ihlen developed MATLAB code in [35] and Doma ński modified the algorithm in [7].
The MFDFA method starts with a possibly non-stationary time series {e i } for i = 1, . . ., N, where N indicates its length.

1.
Define the "profile" E and transform original data into mean-reduced cumulative sums, where ē is the mean of series, such that the aggregated time series are with zero mean.

2.
Divide time series E j into N s = int(N/s) non-overlapping segments of equal length s, starting from the beginning.Since the length N of the series is often not a multiple of the considered time scale s, in order to not miss any piece of data, another set of segments starting from the end of data is made from the end coming to the beginning.As a result, 2N s segments are obtained covering the whole dataset.

3.
Calculate the local trend p k for each of the segments k = 1, . . ., 2N s by a least-square fit of the series.4.
Calculate the mean square error F 2 (k, s) for the estimate of each segment k of length s.

5.
Average over all segments to obtain the qth order variance (or fluctuation) function F q (s) for each size s: For q = 0 use Repeat steps (2)-(5) for different s evaluating new sets of variances F q (s).

7.
Plot F q (s) for each q in log-log scale and estimate the linear fit with least squares.If slope h(q) varies with q, multifractality is suspected.Single slope shows monofractal scaling.8.
Use Legendre transform to evaluate Hölder exponent α(q) and multifractal spectrum f (α): The slope H q of scaling function F q with different q order demonstrates the LRD behavior of signals.The multifractal spectrum indicates how much dominant are the various fractal exponents present in the series.The width of the singularity spectrum denotes the range of the generalized Hurst (Hölder) exponent, which is defined as It is often used to quantitatively measure the degree of multifractality of the series.Thus, the wider the spectrum the more multifractal series is.

Defining the Source of the Multifractality
The origin of multifractality of a time series can be distinguished as two different types, i.e., the multifractality due to: (i) the different long-range correlations of the number fluctuations; and (ii) the broadness of probability density function (PDF) of the distributions [34].
The easiest way to eliminate the correlations for (i) is shuffling the original series into random order, since the multifractality is due to the probability density, which is not affected by the shuffling procedure.For (ii), the surrogate process of data, defined as replacing the phase of discrete Fourier transform (DFT) coefficients of the original data with a set of pseudo independent and identically distributed quantities in (−π, π), can change the broad PDF of the original data into the Gaussian distribution and seldom destroys the intrinsic long-range correlations of the original data [34].

Plot Fitting Hurst Exponents with Crossovers
Generally, the conventional R/S plot fitting is performed by Equation ( 3).The standard estimation with least squares using QR factorization is used for single Hurst exponent estimation.In the current study, the curve is assumed to be piece-wise linear, which is achieved by solving the first-order spline least squares (LSQ) fitting problem with graph log(F q ) versus log(s) identifying crossovers.Each plot includes one, two, and three linear approximations.In the next section, the red line represents single line and estimation of one memory scale, which is the conventional single Hurst exponent H 1 .The green line represents two memory scales H 1 , H 2 with one crossover.The blue line represents three memory scales H 1 , H 2 , H 3 separated with two crossovers.Short-term is important from the perspective of the controller tuning, while the long-term scale indicates the control system structure.Therefore, the asymptotic and transient behavior of the system can be characterized with multiple Hurst exponents.

Case Studies
Semiconductor products are characterized by rapid changes in both improvement and deterioration of quality.Moreover, the embedded seasonal property of the high power plasma etching makes it difficult to analyze the real-time data since the process is nonlinear and non-stationary.
Thousands of runs are performed overnight for each recipe and they might be subject to unknown interventions.The precise control of coolant temperature is the key to the wafer uniformity during the etching process.Data from the process of dynamic coolant supply temperature are anonymously selected out of hundreds of control loops.The industrial part of the analysis is based on the real data stem from etching equipment in semiconductor manufacturing.
For thousands of data logs from different chambers, how to assess the control performance is a thorny problem for process engineers in semiconductor industry.In addition, the tuning of the MIMO system is a challenging task, because the optimization of PID parameters can be unachievable under coupled loops and nonlinearities.In Figure 2, r 1,2,3 , y 1,2,3 are set-points and outputs of the control system respectively.The three PID controllers C 1,2,3 are designed to control the temperature, flow and coolant level in the tank accordingly.Therefore, it is a MIMO system with internal coupling loops and external disturbance and noise with significant delays.In the article, we use Var1, Var2, Var3 to represent the error signals (r 1,2,3 − y 1,2,3 ) from each loop depicted in Figure 3.

Non-Gaussian Statistical Analysis
First, the histogram plot Var1, Var2, Var3 with a probabilistic distribution function (PDF) fitting is frequently carried out to detect if the data are impacted by non-Gaussian noises or subject to human interventions.It is often used to measure the heavy-tailedness underlying process during assessment and frequently the a suitable moment to call for further insight.Histogram fitting is significant since some parts of the process can be more impacted by the disturbances or can be more exposed to the process nonstationarity, can be cross-correlated or might be subject to human interventions.In real situations it is a challenge to have similar properties for all the loops.
The etching process requires uniform temperature behaviors from batch to batch, therefore Var1 is the major concern in the current MIMO system.As introduced in Section 2.3, the α-stable distribution is performed to detect whether the distributions are Gaussian, Lévy or else.In Figure 4, time series Var1 of the MIMO system has significant fat tail properties.Therefore, it can be inferred that there are non-Gaussian noises coming to the system.
Var2 is the variable of errors from flow loop.The higher flow rate can exchange heat efficiently whereas more turbulence could be introduced at the same time.Therefore, the flow loop is interacted with the temperature loop.From the Figures 4 and 5, we could see the clear correlations between Var1 and Var2.The histogram of the Var3 explicitly shows very poor or even not working control.The histogram is highly scattered, almost uniform.The loop is uncontrolled with possible oscillations and operating on the edge of the stability region or manually operated.

Hurst Exponents Fitting with Crossovers
Second, in Figures 6-8 the red line represents single memory scale estimation, i.e., one Hurst exponent fitting without crossover H 1 .Green line represents two memory scales H 1,2 with one crossover.Blue line represents three memory scales H 1,2,3 with two crossovers.All the crossover points are marked in circles.The phenomenon of the multiple memory scales in Var1 and Var2 are observed in Figures 6 and 7. Therefore, the variation of the memory scales with crossovers cannot be captured by the conventional Hurst R/S fitting method with the constant value.After processing data with MFDFA, the multiple Hurst exponents, along with crossover point position might carry on information of multifractality of the control performance.
Var3 is the variable of errors from coolant level loop, which in reality to keep the coolant away from the alarm levels.It should be noted that the control priority is comparatively lower in the design of the MIMO system.Compared with Var1, Var2, however, the monofractal behavior of Var3 can be validated through Figures 8-10.There might be done further conclusions on the control quality of the loops.Var1 and Var2 show multi-persistence behaviour with different Hurst exponents in different memory scales.In contrary Var3 is clearly mono-persistent.It means that this loop is rather not coupled with the other ones.Its Hurst exponent is approximately equal to 1.Its tuning is highly sluggish, almost not existent.It confirms initial observation done using the histogram.The other two loops are the most probably close coupled.Observation of the R/S plots and the shortest Hurst exponent also reveals dynamically sluggish control, however in the longer memory scale control improves and is back to the neutral one.

Multifractal Analysis
In order to clearly derive the origin of the multifractal property, the shuffled datasets of Var1, Var2, Var3 can remove correlations from data and any remaining scaling is caused by probability density function broadness.It is shown that the multifractal behaviors can be captured by the multifractal spectrum analysis with the shuffled datasets in Figure 11.
From multifractal spectrum in Figures 10-12, the multifractal spectrum width is significantly different.Notice the tiny arcs in Figure 10 demonstrate the constant H q for monofractal property in Var3, while the wider arcs of Var1, Var2 shows the multifractality in MIMO system.Moreover, the multifractality of the MIMO process have been detected and analyzed for the each loop with shuffled datasets Var1, Var2, Var3.

Results and Discussion
This paper presented a CPA method with non-Gaussian statistical factors.First, we extracted four parameters α, β, γ, δ with α stable distribution and calculated the multifractal spectrum width ∆(α).
It is well known that the PID tuning of the MIMO system is a time consuming and arduous work.Two sets of real data Var1, Var2, Var3 and Var1', Var2', Var3' with different PID parameters are compared in the end of article, in which the former settings are conservative and latter ones are more aggressive.In Table 1, the proposed FOSP method has been conducted according to the previous sections.The aggressive tuning with small α depicts the broadness and fat tails of the distribution of the underlying process.Moreover, the absolute values of β, γ, δ indicate the worse behavior with the latter tuning.In addition, the multi-scale Hurst spline fit with crossovers are located, compared with single Hurst estimation.Using multiple Hurst exponent fitting at different scales representing the dynamic behavior of the system with noise in Table 2. Review of the persistence properties of the data reveals visible effects of tuning and we observe changes in H 1 .For the first loop it moved from more sluggish value of 0.682 towards neutral tuning, while the second loop Var2 is slightly more sluggish.Loop Var3 is unchanged.However, the effect of tuning (whatever the set is) is evident in comparison with original MIMO system presented in R/S plots (Figures 6-8).Var1, Var2 require disturbance decoupling compared with aggressive tuning while Var3 requires possible redesign of the control philosophy and fine tuning.
As it is well known there is no single universal measure that handles and detects all possible features that may be encountered in the industrial control.Each method well suites to the specific loop properties.Classical methods, such as integral, statisitical or minimum variance methods cover wide are of the applications.However, they may have difficulties in case of nonlinear and complex processes.FOSP approach enables to address controller performance, which depends on the controller tuning and control philosophy suitability.It allows to incorporate into the research the stable properties of the control error histogram.As distribution fitting may be inefficient, the multifractal spectrum can disclose further information about control properties.Complex cascaded properties of the overall system, human impact or superposition of various impacts and oscillations (non-concave spectra) may be additionally assessed.

Conclusions
The multifractal properties of process control variables exhibit fractional order system dynamics.In this paper, we have described the origin of multifractal properties and how the FOSP techniques relate to complex phenomena in the control system.More significantly, CPA of the MIMO control system is first proposed with the help of FOSP through the real semiconductor industry data.Indications for the control improvement are defined and the conclusions upon the control performance of the MIMO system are subsequently derived.The FOSP analysis makes it possible to perform analysis in different scales, allowing to address controller dynamics (short scale) and system cross-dependencies (mostly reflected in longer scales).Similar phenomena have been already analyzed in time-delayed systems and networks, allowing to assess both system parametrization and its structure.Therefore, it is believed that FOSP will find a unique place that conventional signal processing methods cannot even reach the answer to.

Figure 10 .
Figure 10.Multifractal spectrum for shuffled Var3.The parallel behavior with different scales (left plot) and the tiny arch (right plot) indicate the monofractal behaviors of Var3 in the MIMO system.

Figure 11 .
Figure 11.Multifractal spectrum for shuffled Var1.The non-parallel behavior with different scales (left plot) and the arch with large width (right plot) indicate the multifractal behaviors of Var1 in the MIMO system.

Figure 12 .
Figure 12.Multifractal spectrum for shuffled Var2.The non-parallel behavior with different scales (left plot) and the arch with large width (right plot) indicate the multifractal behaviors of Var2 in the MIMO system.

Table 1 .
Indexes for variables with different PID tunings.

Table 2 .
Hurst exponents at different memory scales.