Abstract

Principal component analysis (PCA) is a popular method for process monitoring. However, most processes are time-varying, thus older samples are not representative of the current process status. This led to the introduction of adaptive-PCA based monitoring, such as moving window PCA (MWPCA). In this study, near-infrared spectroscopy (NIRS) responses to digester failures were evaluated to develop a spectral data processing tool. Tests were performed with a spectroscopic probe (350–2,500 nm), using a 35 L mesophilic continuously stirred tank reactor. Co-digestion experiments were performed with pig slurry mixed with several co-substrates. Different stresses were induced by abruptly increasing the organic load rate, changing the feedstock or stopping the stirring. Physicochemical parameters as well as NIRS spectra were acquired for lipid, organic and protein overloads experiments. MWPCA was then applied to the collected spectra for a multivariate statistical process control. MWPCA outputs, Hotelling T2 and residuals Q statistics showed that most of the induced dysfunctions can be detected with variations in these statistics according to a defined criterion based on spectroscopic principles and the process. MWPCA appears to be a multivariate statistical method that could help in decision support in industrial biogas plants.

INTRODUCTION

Anaerobic digestion (AD) involves a complex series of degradation steps catalyzed by at least four trophic groups of microorganisms. The instability of the process is often cited as one of the main hurdles for AD technology implementation (Dupla et al. 2004). Such instability is mainly due to the interdependence of the different microbial groups involved in the degradation (Gujer & Zehnder 1983). The AD process is sensitive to environmental fluctuations (pH, alkalinity) and inhibitors (volatile fatty acids (VFAs), ammonia, etc.), which have a direct effect on the biological processes taking place inside the reactor (Chen et al. 2008). This makes it even more difficult to keep the biological process under optimal and stable operating conditions, especially in a context of mixing different types of organic waste depending on their availability, composition and price. Thus, real-time monitoring of the process is essential to avoid these disturbances and to improve biogas production.

So far, various methods have been used to monitor the process but none seem to be ideal (Madsen et al. 2011). Usually, a set of parameters, expected to be characteristic of the process status, are monitored and interpreted individually (Boe et al. 2010). A survey of 400 full-scale AD plants indicated that 95% of their inline instrumentation was limited to pH, temperature, water flow, biogas flow, level and pressure (Spanjers & van Lier 2006). Those parameters reflect the microbial behavior of the digester to a certain extent. However, ideal monitoring methods should be robust and give early indications of imbalances in the microbial status of the process. To this end, numerous studies have focused on rapid, reliable and inline monitoring tools as well as the development of efficient feedback alert and control strategies for AD (Madsen et al. 2011).

Among these new control strategies are spectroscopic techniques. Spectroscopic techniques, especially near-infrared spectroscopy (NIRS), have been strongly developed for the inline monitoring of biological media and thus of AD. NIRS has the advantage of being non-destructive, rapid and inexpensive. Some studies have been conducted to monitor the AD process using two main multivariate approaches combined with spectral measurements. The first is the use of partial least squares regression (PLS) to produce calibrations, with some success, for performance parameters such as total and volatile solids (TS/VS) (Holm-Nielsen et al. 2007), VFAs (Jacobi et al. 2009), alkalinity (Ward et al. 2011), biochemical methane potential (BMP) (Doublet et al. 2013). One of the weaknesses of the PLS approach is that it only provides the information requested. The second approach is the use of principle component analysis (PCA), which groups samples along axes maximizing their variability, thus highlighting existing clusters in data. It is therefore used as a tool to understand the spectra's variability over time and subsequently the digestate's variability. The PCA approach has been used with NIRS to highlight instabilities in digester operation with PCA score plots (Dias et al. 2008).

A study also investigated the use of spectra derived from Fourier transform-NIRS (FT-NIRS) and predictions of performance parameters within a PCA model as the basis of a process monitoring scheme for an anaerobic digester through the use of Hotelling's (T2) control charts (Reed et al. 2013). All these studies used a fixed PCA model. However, slow and normal changes often occur in real processes (Gallagher et al. 1997), especially with co-digestion, where the substrate quality and quantity could change over time. This leads to false alarms for a fixed-model monitoring approach. For the AD process, it also excludes the time-limited impact of input substrates on the digestate. An efficient multivariate statistical process control (MSPC) method for AD monitoring should therefore require the adaptation of the PCA model to accommodate the time-varying aspect of the process.

Multivariate statistical approaches based on PCA have been widely applied for fault diagnosis of chemical processes to improve their quality and productivity (MacGregor & Kourti 1995). A fixed (static) PCA model is first established based on collected process data under normal operating conditions (NOC) and then used to check new measurement data. For fault detection, differences between the new measurement data and their projections in the built model, the residuals, are then subjected to a statistical test to determine if they are significant. Usually, two statistics, squared prediction error (Q) and Hotelling distance (T2), are used to detect changes in the process (Jackson 1991). They represent the variability in the residuals subspace (Q) and principal components subspace (T2) (Qin 2003). An abnormal situation will cause at least one of the two indices to exceed control limits. As mentioned above, the main hypothesis of such a method is the existence of stable interrelationships between the variables. However, industrial processes are time-varying with characteristics such as (i) changes in the mean, (ii) changes in the variance and (iii) changes in the relationships structure among variables, including changes in the number of significant principal components (PCs) (Jeng 2010). When a time-invariant PCA model is employed, with the aforementioned normal changes, it may be unrepresentative of the status of the entire process, because the training set remains static. As a consequence, new observations may be wrongly classified as deviations from the reference conditions, leading to the reliability of the monitoring system being significantly compromised.

Several PCA-adaptive methods have been developed, such as dynamic PCA (DPCA) for auto-correlated data (Lu et al. 2005), recursive PCA (RPCA) (Li et al. 2000) and moving window PCA (MWPCA) (Jeng 2010) to improve the monitoring of dynamic processes. These methods include updating the PCA model with different datasets according to defined selection parameters (Schmitt et al. 2016). Recursive techniques update the model for an ever-increasing dataset by including new observations without discarding old ones. However, at some point older data may become unrepresentative of the varying process. The MWPCA method can address some of the weaknesses of RPCA by building a fitted model using its moving window. As the window slides along the data, a new process model is generated by including the newest samples, which are more representative of the current process operation, and discarding the oldest.

The MWPCA approach has been tested on simulated data of continuously stirred-tank reactor processes (He & Yang 2008; Jeng 2010). However, the use of MWPCA on infrared data is still new, with few applications to processes such as crystallization in a pharmaceutical context (Taris et al. 2017). The window size determination can be empirical (He & Yang 2008) or standard (Schmidt et al. 2016) but an appropriate definition of alarm cases should be performed and discussed in relation to the studied system (Taris et al. 2017). The AD process is a perfect example of a varying process with complexity, which is its biological dimension. Indeed, due to the different biochemical steps encountered during the process, the definition of a reference state for these MSPC methods might be difficult. Moreover, the use of NIRS in AD is still challenging, as spectra are affected by the presence of water, particle size, and light scattering (Giordanengo et al. 2008; Xie et al. 2009).

However, the application of MSPC based on NIRS spectra for AD monitoring could reduce time-consuming analyses and enhance early detection of instabilities in the process. This study investigated the MSPC technique adapted for AD process monitoring based on NIRS data. Static PCA and MWPCA methods were tested with different experiments in normal or dysfunctional states. T2 and Q statistics control charts were used to detect instabilities in the process. Moreover, score plots and loadings helped identify which spectral variables could change and might be affected. The results section shows the comparison between both methods applied to infrared measurements for AD monitoring.

MATERIALS

Experimental setup

A pilot-scale continuously stirred tank reactor (CSTR, 35 L) was used for these experiments. The digester was operated at 38 ± 1 °C. The reactor was mounted on a scale in order to continuously monitor its weight, thereby controlling and recording feeding doses. The feedstock was prepared every week and stored at 4 °C in a hopper above the reactor. This allowed gravity feeding through an automated valve at the bottom of the hopper. The entire system was managed via an Agilent 34970a acquisition control unit with a specific program developed under Labview2014. The reactor was equipped with a gas flow meter connected to the acquisition program, allowing measurement and recording of kinetics and cumulative production of biogas. The exhaust biogas pipe was equipped with a sampling point to determine biogas quality. Biogas analysis was performed by gas chromatography equipped with a flame ionization detector (GC-FID) (Agilent Technologies, 6890N, USA) and a column Porapak Q according to the method described in Lucas et al. (2007).

Various experiments were conducted using this pilot-scale CSTR. Pig slurry was co-digested with different co-substrates at organic loading rate (OLR) varying between 1.5 and 5 kgCOD·m−3·d−1 and a hydraulic retention time (HRT) varying between 20 and 28 days (Table 1). For all experiments, the reactor was fed once a day. The co-substrates were horse feed residues, fruit waste, food fats and protein waste, and they were added to the pig slurry in different proportions (Table 1). The digester operated normally for 120 days before being subjected to different stresses. Feeding was usually stopped for several days when the methane concentration in the biogas dropped, to allow the digester to recover from each stress. It was sometimes necessary to empty the reactor to overcome restarting issues. Prior to these experiments, the digester was inoculated with digestate from a well-established laboratory-scale digester treating pig slurry and horse feed residues (HRT = 27 days, OLR = 1.01 kgVS·m−3·d−1) at a pH of 7.43, TS of 19.2 g/L and VS of 11.1 g/L. This same inoculum was used to restart the digester after each overload at the end of experiments Nº3, Nº4, Nº5, and Nº6.

Table 1

Operating conditions of the different experiments carried out

ExperimentsFeedstock compositionDigester operationsOLR (kgCODm−3·d−1)HRT (days)
Pig slurry + horse feed residues (20 g·L−1Normal 1.4 20.6 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Normal 2.6 20.8 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1) + food fats (20 g·L−1Overload occurred 4.9 21 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Normal 1.4 28.2 
Pig slurry + horse feed residues (40 g·L−1) + fruit waste (540 g·L−1Overload occurred 3.01 23.7 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Stirring stopped 1.7 27.2 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Normal 1.6 26.4 
Pig slurry + fruit waste (135 g·L−1) + protein waste (20 g·L−1Overload occurred 3.8 23.7 
ExperimentsFeedstock compositionDigester operationsOLR (kgCODm−3·d−1)HRT (days)
Pig slurry + horse feed residues (20 g·L−1Normal 1.4 20.6 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Normal 2.6 20.8 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1) + food fats (20 g·L−1Overload occurred 4.9 21 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Normal 1.4 28.2 
Pig slurry + horse feed residues (40 g·L−1) + fruit waste (540 g·L−1Overload occurred 3.01 23.7 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Stirring stopped 1.7 27.2 
Pig slurry + horse feed residues (20 g·L−1) + fruit waste (270 g·L−1Normal 1.6 26.4 
Pig slurry + fruit waste (135 g·L−1) + protein waste (20 g·L−1Overload occurred 3.8 23.7 

In addition to the experiments carried out under normal conditions, four stress cases will be thoroughly discussed in this study. They are overload induced through feedstock mainly composed of carbohydrates, lipids or protein and a stirring defect case.

Sampling and reference analysis

Digestate samples (∼1.7 L) were taken daily immediately before feeding the reactor. Physicochemical analyses of digestate and feedstock (∼1 L) were performed once a week. TS/VS, total Kjeldhal nitrogen (TKN), total ammonium nitrogen (NH4+), chemical oxygen demand (COD) and pH were analyzed using standard methods (APHA 2012). VFAs (VFAHPLC) were analyzed by high performance liquid chromatography (HPLC, Varian©, U3000). VFAtit/alkalinity ratio (VFA/alk) was determined by titrimetry (Titralab AT1000, Hach Lange). Long-chain fatty acids (LCFA) were extracted with hexane/isopropanol (3/2) using an accelerated solvent extractor (Dionex, ASE 350). After evaporation of the solvent, LCFA content VFAs was determined by gas chromatography/mass selective detection (GC/MS, Agilent Technologies, 7890B/5977A) equipped with a DB-FFAP column (30 m × 0.32 mm i.d., 0.25 mm film thickness) as described in Zhang et al. (2015). Biogas production and quality (CH4/CO2) were monitored daily for all experiments. During overload experiments, parameters such as pH, VFAHPLC, VFA/alk and LCFA were also monitored daily on digestate samples.

Spectral acquisition

All digestate samples were maintained at the reactor bath temperature and scanned using a near infrared instrument in reflectance mode while fresh (i.e. without any pretreatment, such as drying). Spectra were taken at a frequency of 1–7 spectra per week. Acquisition frequency was increased (daily acquisition) for overload experiments to pick up all information on the process state. For spectra acquisition, an off-line spectroscopic probe was used on all collected digestates. The system (Figure 1) consisted of two fibers, one fiber for lighting and a second fiber for signal collection. The core diameter of the two fibers was 1,000 μm with a numerical aperture of 0.39 (BFY1000, Thorlabs). The principle of this probe was based on diffuse optical spectroscopy. The device was linked to a tungsten-halogen light source (Ocean Optics HL-200-FHSA) and a LabSpec1 spectrometer (ASD, Boulder).

Figure 1

Schematic of the NIRS two-fiber probe.

Figure 1

Schematic of the NIRS two-fiber probe.

The spectral range extended from 350 nm to 2,500 nm with a step of 1 nm, with resolutions of 3 nm (350 nm–1,000 nm range) and 10 nm (1,000 nm–2,500 nm range). The intensity of the reflected light (I(λ)) was measured for each sample. Dark current (In (λ)) was recorded for all measured spectra and subtracted. To overcome the spectral responses to the various optical components of the instrumentation (light source, fibers and spectrometer), a white reflectance reference (SRS99, Spectralon®) was always measured (I0 (λ)). For each sample, the spectrum was measured three times. From these measurements, a reflectance (R (λ)) was calculated for each sample as follows: 
formula
(1)

For subsequent calculations, the spectra were normalized by applying the logarithm (absorbance) and smoothed using the Savitzky–Golay algorithm (Savitzky & Golay 1964) in order to reduce the scattering effect and delete the baseline (Zeaiter et al. 2005). Finally, the spectra were truncated to the 450–2,200 nm range in order to exclude noise present in the spectrum.

METHODS: MULTIVARIATE DATA ANALYSIS

Static PCA method

The primary purpose of PCA is to explore or describe a large multivariate dataset. As a multivariate method, PCA allows one to reduce the size and extract pertinent information from experimental datasets while preserving the maximum variance of the original space. The size reduction is an indicator of the data structure and of the underlying chemical, physical or biological phenomena. Considering a generic spectra data matrix X (n × p), the PCA model is written according to Equation (2) below: 
formula
(2)
where n is the number of observations and p the spectral variables (e.g. wavelengths). T is the score's matrix with a (a = 1, 2 …) being the number of retained PCs which are orthogonal and define the new subspace of the spectral variables. PT is the transposition of the loadings matrix which explains the relationship between spectral variables and PCs; E is the residual matrix. Variables are in general mean centered and scaled to unity variance before PCA implementation. Having established a PCA model, the process monitoring is generally reduced to multivariate control charts as T2 and Q. Therefore, the monitoring and fault detection system is implemented by means of two major steps:
  • Calibration: a model is built on training data acquired in the NOC region. From this model, the scale parameter vectors, the mean x and standard deviation s, the covariance and the loading matrix P, are computed.

  • Prediction: new observations are projected onto the PCA model to evaluate whether they behave accordingly to the training set. The new observation xi is scaled using the scale parameter vectors x and s. T2 and Q statistics of the new observation are also evaluated using the PCA model (Equations (3) and (4)). 
    formula
    (3)
     
    formula
    (4)
    where Λ is the covariance matrix of the a scores column and is the i-th observation predicted through the PCA model. To detect abnormal process behavior, T2 and Q statistic upper control limits are determined. If one of these exceeds the upper limits, this measurement is considered as an alarm. For the T2 statistic, the upper limit is expressed by the means of Equation (5) (MacGregor & Kourti 1995), where is the limiting value based on the Fisher–Snedecor distribution with (A, N − A) degrees of freedom and α (90% ≤ α ≤ 95%), the level of significance. 
    formula
    (5)
With the Q statistic describing how well the PCA model predicts the xi vector, abnormal events can also be detected by its upper limit, expressed in Equation (6) (Jackson & Mudholkar 1979) below: 
formula
(6)
with 
formula
where cα is the value of the normal distribution with α, the level of significance. Some studies have suggested that some consecutive number of alarms should be detected before establishing that an uncommon event has occurred.

MWPCA method

In MWPCA, the basic algorithm stipulates that a data window of fixed length is moved in real time to update the PCA model once a new sample is available. With each new observation, this window excludes the oldest observation and includes the newest one. One challenge in implementing MWPCA is to select the window length. The window size should be large enough to cover sufficient sample data for modeling and monitoring. A strict lower bound for the window size is at one observation more than the number of retained components used in the model. There are different criteria for the selection of the window size. Schmitt et al. (2016) proposed a selection method based on the minimization of the sum of square predicted errors (SSPE) of a dataset (either calibration or prediction), chosen from the NOC region. The goal is to explore MWPCA model performance for different window sizes on the validation set, rather than fault detection. Then, given a validation set Xv, the sum of squared prediction errors of the adaptive PCA model for each observation in Xv is: 
formula
(7)
where c is the calibration set length from Xc; v the validation set length and is the one-step-ahead prediction of xi given the PCA model from time i1. After choosing the proper window size, a training set is selected to build a PCA model; new data are projected on the PCA subspace. Thus, for a window size k, the data matrix at time t is Xt = (xt−k+1, xt−k+2, …, xt) and at a time t+1 it is Xt+1 = (xt−k+2, xt−k+3, …, xt+1). Scale parameter vectors can then be calculated using observations in the new window. T2 and Q statistics are evaluated in the new observation for fault detection. If the T2 or Q statistic exceeds the limits because the observation is a fault or an outlier, the model is not updated; otherwise it is included in the new training set and the oldest observation is removed (Jeng 2010).

MWPCA method for AD process monitoring

In this work, NIR spectra collected during anaerobic co-digestion experiments represented the varying data. A MWPCA-based approach was proposed to detect changes that could lead the process to dysfunction and then give an early alarm. Since feedstock, inhibitors and environmental factors influence the AD process, the detected changes or faults will be related to the imbalances of the process. An deviation of spectra from those collected before could therefore be considered as an out-of-control state. However, AD is a biological process. Therefore, an out-of-control observation will not be excluded from the next window, contrary to the usual algorithm. This is due to the fact that this ‘new state’ of the process will affect the next observations and could even result in a new normal steady state of the digester operation, especially for the AD process, which has an evolving steady state.

The method implemented here does not require any prior knowledge of the evolving peaks or spectral regions. The adapted MWPCA algorithm is shown in detail in Figure 2. A ‘k’ optimal window size was first defined for MWPCA. Xk spectra were then considered for the first training set and a PCA model was built with determination of the loading matrix P, the number of PCs, T2 and Q limits (Equations (5) and (6)). The next collected spectrum was projected onto the PCA subspace, and T2 and Q were estimated according to Equations (3) and (4). The criterion adopted for fault detection according to the AD process was as follows: the new observation had to exceed the control limits simultaneously for both statistics before declaring the beginning of a dysfunctional state. Then, the window was still moved forward for new detection. Different cases of alarm were analyzed to assess the validity of this criterion based on the contribution plot allowing the identification of the spectral regions that significantly changed. The analysis was carried out with in-house procedures written in MatlabR2013b equipped with the Statistics toolbox.

Figure 2

Moving window PCA algorithm implemented to detect AD dysfunctions.

Figure 2

Moving window PCA algorithm implemented to detect AD dysfunctions.

RESULTS AND DISCUSSION

Process operation

Table S1 (Supplementary Materials) provides a summary of the results for concentrations of parameters measured on digestate samples during the experiments. Parameters such as TS, VS, COD and TKN show variations mostly between the beginning and the end of an experiment. Such differences are mainly substrate dependent. The reactor showed good degradation performance over the first 100 days of the study without VFA or LCFA accumulation, corresponding to experiments N°1 and N°2 (Figure 3). After those two experiments, various failures occurred in the digester after 2–3 weeks at a normal load rate as illustrated in Figure 3, which shows NH4+, VFAHPLC and/or LCFA accumulation and biogas production decrease.

Figure 3

OLR, biogas production, pH, VFAHPLC, LCFA and NH4+ variations for experiments N°1–6. Pig slurry (PS); horse feed residues (HF); fruit waste (FW); food fats (F) and protein waste (PW).

Figure 3

OLR, biogas production, pH, VFAHPLC, LCFA and NH4+ variations for experiments N°1–6. Pig slurry (PS); horse feed residues (HF); fruit waste (FW); food fats (F) and protein waste (PW).

For experiment N°3, the digester encountered lipid overload leading to LCFA accumulation. LCFA has been reported to be inhibitory to methanogen activity (Lalman & Bagley 2002). Due to the amphiphilic properties, LCFA is easily adsorbed onto a microbial surface, thereby impeding the passage of essential nutrients through the cell membrane (Pereira et al. 2005). By adding fats, the biological media was not only exposed to such inhibition, the OLR also increased (from 3 to 5). However, the digester worked fine for 20 days at this OLR before biogas production dropped (day 123); which showed that the main cause of this dysfunction was LCFA accumulation. As suggested by some studies, results obtained in this study confirmed that LCFA started accumulating (day 113) before VFA (day 120). It was also shown that inhibitory values for LCFA are in the range of a few milligrams (biogas production dropped around 0.7 g/L). This makes LCFA a good indicator for lipid overload detection. However, such parameters are not easily available in industrial units.

For experiment N°4, by increasing the OLR of the feedstock from 1.5 to 3, the methanogenesis rate became unable to successfully deal with the uncoupling of the acidogens/acetogens and methanogens, therefore leading to acidosis. VFAHPLC accumulation allows quick identification of the failure (starting from day 303) even if the concentration values were still under 1 g/L. There was a peak of biogas production on day 304, before it started dropping due to organic overload and the associated acidosis. However, field parameters such as VFA/alk and pH reacted weakly and later during the experiment, with an increase of VFA/alk (to 0.5) and a decrease of pH (to 7.2). The stirring defect experiment (experiment N°5) had similar dynamics and results, given that stirring improves contact between biomass and the feedstock. This absence of contact disturbed the process, which was highlighted by VFAHPLC accumulation and later by the drop in biogas production.

In experiment N°6, protein overload provided interesting results because of its capacity to maintain an acceptable biogas production during the experiment. VFA accumulated up to 13.5 g/L, while the VFA/alk ratio tended to remain at 1 due to an established balance between the acidity of the medium (VFA) and the buffer capacity (alkalinity). This buffer capacity is due to ammonium ions (NH4+) produced by protein degradation increasing their concentration from 2 gN/Kg in previous experiments to 4 gN/Kg. In this experiment, pH decreased slowly and stayed around 7.5. However, the digester was no longer working in optimal conditions due to the lower biogas production rate and VFA accumulation.

Therefore, as shown by all these results, potential indicators often used in industrial units to detect failure in the AD processes (pH, VFA/alk) might not be sufficient to prevent dysfunctions in the process. As for biogas production, it can be seen in Figure 3 that acidosis was preceded by an increase of production in overload experiments. This peak of production was reached when the process was already accumulating inhibitory compounds (VFA, LCFA, NH4+). This would certainly create confusion in biogas plants when the production drops quickly. Consequently, additional sensors are required for a better management of the process, particularly in industrial units. Those sensors could be fitted into a decision support system to quickly prevent process dysfunctions.

Static PCA applied to NIRS-derived data for AD process monitoring

Static PCA was first implemented to assess its ability to detect faults. As stated above, two steps are necessary for the monitoring of a process by a static PCA: calibration and prediction. First a static PCA model has to be built using a training dataset in the NOC regions. To meet these requirements, experiment N°2 was selected for the training step as it worked well and gave satisfying results with biogas production maintained at a good flow rate of 42.5 ± 4.5 NL/day (Figure 3) and CH4 at 60.2 ± 1.3% of biogas. The collected spectra during this experiment were baseline corrected and preprocessed through mean-centering, and a PCA model was built. The first component explained 88.44% of the variance, the second 6.31%, the third 2.91% and the fourth 1.65%, while each of the subsequent components added less than 0.5%. Therefore, four PCs were believed to describe most of the variance present in the evolving system. As prediction sets experiments N°3 and N°4 were used to simulate two cases: one where the prediction set had a slightly different feedstock (experiment N°3 with food fats) and another where it was the same feedstock as in the experiment used in the training step (experiment N°4). Those two experiments also corresponded to lipid and carbohydrate overloads respectively. After applying scale parameters and projecting these new spectra on the built model, T2 and Q statistics were estimated for both training and prediction sets (Figure 4). As expected, the statistics values calculated for the training set fell under the control limits.

Figure 4

Static PCA monitoring results: T2 and Q with control limits for experiments N°3 (a) and (b) and N°4 (d) and (e). Comparison between the process and control charts for each experiment (c) and (f).

Figure 4

Static PCA monitoring results: T2 and Q with control limits for experiments N°3 (a) and (b) and N°4 (d) and (e). Comparison between the process and control charts for each experiment (c) and (f).

For the first prediction set (experiment N°3), T2 and Q (Figure 4(a) and 4(b)) were shown to exceed the control limits for some consecutive observations, mostly the last ones. To determine the accuracy of these alarms, the dynamics of these statistics were compared to the biogas production and LCFA accumulation. As shown in Figure 4(c), there was an alarm on day 114 corresponding to the start of LCFA accumulation in the digester. This first alarm occurred before VFAHPLC started accumulating, with high biogas production and VFA/alk at 0.16. This event could be considered to be an early sign of imbalance. There was also a good accuracy with the last observations control charts when biogas production started dropping. T2 and Q were above control limits from day 125, with: 1.3 g/L of VFAHPLC; 0.25 of VFA/alk and 0.92 g/L of LCFA in the digester (Figure 4(c)). The system could then be considered to be in an out-of-control state. An early warning can be established after a sensitivity study of the first alarms detected.

For the second prediction set (experiment N°4), Q (Figure 4(e)) was above the control limit for all spectra. For T2 (Figure 4(d)), one out of two observations was found above the control limit and there was not consistent trend. This suggested that almost all acquired spectra in this experiment were in an out-of-control state. However, comparing these statistics dynamics to the digester history (Figure 4(f)), it can be seen that the digester worked fine for the first 20 days of the experiment before being subject to acidosis. There was no VFA accumulation and nothing suggested a dysfunction in the process during those 20 days. This clearly shows that the new steady state defined in the experiment N°4 was not considered by the static PCA model. The training set was no longer representative of the process.

The first established model in the training set was still representative of the spectra in experiment N°3 because of the continuity between the calibration set and the prediction set. We could assume that the evolution took place gradually in the digester and did not affect the robustness of the model. However, when the same model was used to monitor a new similar experiment with no continuity with the training set, it was no longer representative. These results demonstrate that static PCA led to questionable outcomes when used to monitor dynamic systems. Indeed, in the case of the AD process, a new steady state will be considered to be inconsistent with the training set behavior. This justifies the use of PCA-adaptive methods for AD monitoring as they can update the training set and resolve the ‘new steady state’ issue.

MWPCA APPLIED TO NIRS-DERIVED DATA FOR AD PROCESS MONITORING

Optimal window size determination

MWPCA was then implemented according to the procedure shown in Figure 2; spectra preprocessing through mean, centering and unity variance standardization were performed each time the window was moved. Also, the algorithm was restarted from the beginning each time the digester had to be totally or partially emptied and replaced with a new inoculum to overcome restarting issues.

To choose the window size k, an SSPE test (Schmitt et al. 2016) was performed with the dataset from experiment N°2 (12 spectra) in the NOC region. SSPE was then calculated for a varying window length k; 6 ≤ k ≤ 10. The validation set was the last two spectra. The PCA model was built with four PCs. As shown in Figure 5, SSPE decreased as the window grew; so there was no real minimum. We then added the representativeness of the modeled space T2* (sum of the corresponding T2) in order to choose the window that had a small SSPE and a high representativeness. It was noted that at k = 7, the model considered the largest amount of information (Figure 5), unlike for other window sizes. Along with the decreasing SSPE, which is the main criterion, T2* allowed discrimination between all tested window sizes. Then the moving window length was defined for k= 7.

Figure 5

SSPE and T2* calculated while varying the window size. The optimum window value was observed for k = 7.

Figure 5

SSPE and T2* calculated while varying the window size. The optimum window value was observed for k = 7.

Cases of fault detection according to the AD process

As mentioned above, a criterion of T2 and Q simultaneously above control limits was selected to validate a dysfunctional state detection. To better understand this criterion, this section explains and details the possible cases of variation of T2 and Q statistics with examples based on the process. MWPCA has a short-term memory, thus it is important to understand NIRS spectra deviation based on spectroscopic principles (light scattering, optical path …) specific to this type of medium and to the process. This first step will enable the elimination of false alarms.

First, all digestate spectra (Figure S1, Supplementary Materials) have water spectral signature with absorption peaks at 970 nm, 1,200 nm, 1,440 nm and 1,940 nm (Roggo et al. 2007). The difference of several spectra will be a function of the absorbance intensity as a result of a combination of digestate chemical composition and the light's optical path, related to the light scattering, in the digestate. These characteristics are mostly related to physical particles in the digestate. This suggests that the physical aspect and light scattering of the digestate have a great influence on the infrared signal but water absorption bands could bring some non-negligible chemical information. Another characteristic of these spectra is the occasional presence of a feature at 2,200 nm, which is a dip between two peaks of water on an absorbance spectrum. This dip is most likely due to the appearance of another water adsorption peak at 2,250 nm which signals a change in the light's optical path.

Loadings will therefore represent global intensity variation: displacement or apparition of water band absorption or light scattering by the presence of a slope between 400 and 1,000 nm. The associated PCs are orthogonal and can all highlight the same information. Changes in the loading associated with a PC might occur between two PCs but it will not affect the T2 or Q statistic since they include all loadings in their determination. Table 2 highlights those possible changes in the spectra through their loadings and scores and how they could lead to T2 and Q variations, or not. It is important to note that the model, once updated, can adapt itself to those variations. This adaptation makes it possible to take into account new steady states. A model that will not readjust itself will suggest continuous and rapid changes or events occurring in the digester and therefore will highlight imbalances of the process.

Table 2

Monitoring cases for T2 and Q statistics, score plots and loadings variations of the built PCA model

CasesT2QScore plotsLoadingsObservations
↓ ↓   The eighth spectrum projected is well represented by the PCA subspace regardless of the axis and practically superimpose with the mean of the seven previous spectra. Loadings 1 and 2 are characteristic of the water absorption dip at 1,440 nm and the global spectra intensity respectively. There is no indication of change, so the process is operating normally. 
↑ ↓   The eighth spectrum projected is away from the other spectra in the PCA subspace. The shapes of loadings 1 and 3, similar to the spectra mean, show that their variations are due to the spectra intensity. Loading 2 is characteristic of the media scattering light (slope over smaller wavelengths). This indicates the presence of more particles in the sample. The feedstock probably increased, inducing a ‘new state’ in the digester, but has not cause dysfunction. 
↓ ↑   This case is similar to the first case regarding the eighth spectrum projection and the characteristics of loadings 1 and 2. There might be some light scattering (loading 4) in the media, but there is no indication of change in the process. A high residual Q mainly implies that the current model was not able to address all the sample information. The next observation might be a fault if the next window does not adjust its model. 
↑ ↑   This case is similar to the second case in terms of the eighth spectrum projection and the forms and characteristics of loadings 1 and 2. In addition, loadings 3 and 4 respectively show a flattening (1,500–2,000 nm) and a displacement of the 1,440 nm water peak. This not only indicates the presence of more particles in the sample but also a change in the light optical path (digestate whitening). The feedstock increase has caused dysfunction. 
CasesT2QScore plotsLoadingsObservations
↓ ↓   The eighth spectrum projected is well represented by the PCA subspace regardless of the axis and practically superimpose with the mean of the seven previous spectra. Loadings 1 and 2 are characteristic of the water absorption dip at 1,440 nm and the global spectra intensity respectively. There is no indication of change, so the process is operating normally. 
↑ ↓   The eighth spectrum projected is away from the other spectra in the PCA subspace. The shapes of loadings 1 and 3, similar to the spectra mean, show that their variations are due to the spectra intensity. Loading 2 is characteristic of the media scattering light (slope over smaller wavelengths). This indicates the presence of more particles in the sample. The feedstock probably increased, inducing a ‘new state’ in the digester, but has not cause dysfunction. 
↓ ↑   This case is similar to the first case regarding the eighth spectrum projection and the characteristics of loadings 1 and 2. There might be some light scattering (loading 4) in the media, but there is no indication of change in the process. A high residual Q mainly implies that the current model was not able to address all the sample information. The next observation might be a fault if the next window does not adjust its model. 
↑ ↑   This case is similar to the second case in terms of the eighth spectrum projection and the forms and characteristics of loadings 1 and 2. In addition, loadings 3 and 4 respectively show a flattening (1,500–2,000 nm) and a displacement of the 1,440 nm water peak. This not only indicates the presence of more particles in the sample but also a change in the light optical path (digestate whitening). The feedstock increase has caused dysfunction. 

↑ above control limits ↓ under control limits.

For each case shown in Table 2, a PCA model is built on seven spectra (as the determined window size is k = 7) in the NOC region with four PCs. An eighth spectrum is projected on the PCA subspace to simulate variations of T2 and Q statistics. Score plots indicate how the new spectra behave in relation to the others and loadings come to relate the process to the variations highlighted by the model. The position of this new spectrum is discussed based on the principle of fault detection related to the center of the model and to other samples as described, for example, by Yamanaka et al. (2012) (see Figure S3, Supplementary Materials). Standard normal variate (SNV) (Barnes et al. 1989) was performed on the mean, loadings and the eighth spectrum for a better visualization. All the studied cases are based on spectra from samples used in this study and processed by the MWPCA routine.

Lipid overload

Experiments N°1, N°2, and N°3 were carried out one after another, without interruption, after a classical load increase. The MWPCA routine started at the beginning of experiment N°1. The first PCA model was built using four latent variables explaining a cumulative variance of 99.82%. The window was then slid along experiment N°2. As previously illustrated (Figure 4), a new steady state (day 55) was established with the addition of fruit waste to the feedstock. This new state was properly taken into account by the two statistics, which remained under control limits (Figure 6). On day 71, Q was above the control limit, which shows the model was not representative at this point but has readjusted itself in the next window. The results for the monitoring of experiment N°3 showed that after encountering a ‘new state’ alarm on day 106 (only T2 above limit), there was an early warning (compared to existing indicators on biogas plants) starting on day 114 (T2 and Q above limits) that corresponded to the beginning of LCFA accumulation, as previously shown by static PCA. The last observations were also above control limits for both statistics, which related to the acidosis. As explained in Table 2 above, continuous and rapid changes were occurring in the digester that prevented the model from readjusting, thus highlighting a dysfunctional state. At some point the spectra returned under control limits but the digester had not recovered yet (VFAHPLC 7.8 g/L, LCFA 1.7 g/L and VFA/alk 1.8). This was due to the fact that the current model was made using dysfunctional (high variance) spectra. Then, a potential recovery could also be detected by T2 and Q being above control limits given that this MWPCA algorithm applied to the AD process is somehow a state change detector.

Figure 6

Comparison of biogas production, LCFA and T2 and Q statistics for experiments N°1, N°2 and N°3.

Figure 6

Comparison of biogas production, LCFA and T2 and Q statistics for experiments N°1, N°2 and N°3.

Organic overload

Experiment N°4 took place after a total emptying of the digester. The algorithm was then restarted and a first training set was defined with the same window length. The first PCA model was developed also using four PCs for 99.60% of cumulative variance explained. The monitoring results are depicted in Figure 7(a). As a reminder, a regular load increase was performed until a full normal load was reached on day 295 and the OLR increased abruptly on day 301. There was an early warning on day 304 when VFA/alk was 0.2 and VFAHPLC 0.7 g/L and biogas production had just slightly decreased. This slight decrease of production would not be remarkable in biogas plants as the production is buffered by gas storage. To recover more rapidly from the acidosis, the feedstock was restricted to pig slurry only (low OLR) on day 306. This resulted in a drop of the biogas production on day 307 and another warning was found, which was in good agreement with the process.

Figure 7

Comparison of biogas production, VFAHPLC and T2&Q statistics for experiments N°4 (a) and N°5 (b).

Figure 7

Comparison of biogas production, VFAHPLC and T2&Q statistics for experiments N°4 (a) and N°5 (b).

Stirring defect

For experiment N°5, 30% of the digester was changed before restarting feeding. The collected spectra were monitored through MWPCA and the first training set explained 99.80% of cumulative variance. After a regular load increase, the feeding was increased to a normal charge on day 323 causing an adjustment of the model (Figure 7(b)). Stirring was turned off starting from day 328 and biogas production decreased slightly. From this moment, a settling gradually occurred in the digester, which became substantially free of suspended matter on day 330. The settling caused a change of the light's optical path and the spectra's global form changed mostly between 1,400 and 2,300 nm (Figure S2, Supplementary Materials). Obviously, both T2 and Q were above control limits at this point (VFA/alk 0.21 and VFAHPLC 0.9 g/L) and returned under control on day 332 after the stirring was switched on. This confirmed that the digestate's physical aspect has a high impact on the monitoring but still remains correlated with chemical events occurring in the digester.

Protein overload

This last experiment started after 50% of the digester was emptied. For the first 10 days the digester was fed with a mix of pig slurry, horse feed and fruit waste. Protein waste was abruptly added from day 356. To apply the MWPCA algorithm, a first training set was defined and the PCA model built explained 99.78% of cumulative variance. Looking at the control charts depicted in Figure 8, there was an alarm on day 356 just before the beginning of the protein addition. This suggests that there was already a potential dysfunction in the process but it was not correlated or confirmed by any physicochemical parameters. This alarm is therefore solely attributable to change in the spectral characteristics with interferences that might be related to changes in the particle size in the digestate or to the fouling of the probe. VFAHPLC started accumulating the day after protein waste was added to the feedstock (day 357). Biogas production dropped slightly, after a peak on day 359 (48.6 NL/day), but maintained a flow rate of 32 NL/day until the end of the experiment. The control charts reacted at the same moment and several alarms were identified until the end of the experiment. This confirms that even though the process apparently had a good biogas production rate, several parameters like VFAHPLC (0.2–13 g/L), VFA/alk (0.2–1) and NH4+ (2–4 gN/kg) showed that the digester was no longer working in the initial optimal conditions. The control charts highlighted continuous change, which indicated that the apparent balance was not real and confirmed that the digester was in a dysfunctional state.

Figure 8

Comparison of biogas production, VFAHPLC and T2 and Q statistics for experiment N°6.

Figure 8

Comparison of biogas production, VFAHPLC and T2 and Q statistics for experiment N°6.

Given the various cases studied, it was shown that false alarms (only Q or T2 above control limits) were not always included on control charts by new models in the MWPCA routine, which led to some interpretation difficulties. Besides infrared interferences (probe fouling, light scattering, etc), the changes of the feedstock was often at the root of this problem. However, in biogas plants where feedstock varies more often, several other types of false alarm could arise. Thus, one of the weaknesses of MWPCA is its short-term memory, the length of which is the window size of the model. Combining RPCA and MWPCA is very effective for adaptive process monitoring and is expected to have broad applicability in industrial processes (Jeng 2010). This combination could increase the memory of the control charts, while simultaneously updating the calibration datasets. A mix of both methods can be more useful and easier to implement in an expert system for AD process monitoring.

It is also important to note that these pilot-scale experiments were conducted with highly biodegradable substrates leading to a low average TS value (18.8 gTS/Kg) and therefore to a low absorbance of the acquired spectra. Working on industrial plants could increase the TS and subsequently the signal intensity perceived by the NIRS probe. Spectra acquisition was made off-line for analysis convenience in this study, but a directly connected probe would increase sampling frequency, which might increase the MWPCA model's robustness.

CONCLUSION

MWPCA was performed on NIRS spectra data for the detection of overloads or mechanical disturbances in an anaerobic digester. Using spectroscopic principles and aspects, we showed that the spectra change according to the chemical (intensity of the absorption) and physical (scattering light, optical path) properties of the digestate. An analysis of T2 and Q variations allowed the definition of a criterion of fault detection when both T2 and Q were above control limits. This criterion, defined according to the digester dynamic, showed great prediction accuracy in the process monitoring. By combining the digester operation and spectroscopic principles, the adopted criterion seems to answer the questions raised by giving early warnings compared to other existing field indicators. The MWPCA model can be introduced to provide decision support that will allow an expert or expert system to monitor the AD process.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/wst.2020.117.

REFERENCES

REFERENCES
APHA
2012
Standard Methods for the Examination of Water and Wastewater
, 22nd edn.
American Public Health Association, American Water Works Association, Water Environment Federation
,
Washington, DC, USA
.
Barnes
R.
Dhanoa
M.
Lister
J.
1989
Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra
.
Appl. Spectrosc.
43
,
772
777
.
https://doi.org/10.1366/0003702894202201
.
Boe
K.
Batstone
D. J.
Steyer
J. P.
Angelidaki
I.
2010
State indicators for monitoring the anaerobic digestion process
.
Water Res.
44
(
20
),
5973
5980
.
https://doi.org/10.1016/j.watres.2010.07.043
.
Chen
Y.
Cheng
J. J.
Creamer
K. S.
2008
Inhibition of anaerobic digestion process: a review
.
Bioresour. Technol.
99
(
10
),
4044
4064
.
https://doi.org/10.1016/j.biortech.2007.01.057
.
Dias
A. M. A.
Moita
I.
Páscoa
M.
Alves
M. M.
Lopes
J. A.
Ferreira
E. C.
2008
Activated sludge process monitoring through in situ near-infrared spectral analysis
.
Water Sci. Technol.
57
(
10
),
1643
1650
.
https://doi.org/10.2166/wst.2008.147
.
Doublet
J.
Boulanger
A.
Ponthieux
A.
Laroche
C.
Poitrenaud
M.
Cacho Rivero
J. A.
2013
Predicting the biochemical methane potential of wide range of organic substrates by near infrared spectroscopy
.
Bioresour. Technol.
128
,
252
258
.
https://doi.org/10.1016/j.biortech.2012.10.044
.
Dupla
M.
Conte
T.
Bouvier
J. C.
Bernet
N.
Steyer
J. P.
2004
Dynamic evaluation of a fixed bedanaerobic digestion process in response to organic overloads and toxicant shock loads
.
Water Sci. Technol.
49
(
1
),
61
68
.
https://doi.org/10.2166/wst.2004.0019
.
Gallagher
V. B.
Wise
R. M.
Butler
S. W.
White
D. D.
Barna
G. G.
1997
Development and benchmarking of multivariate statistical process control tools for a semi- conductor etch process; improving robustness through model updating
.
Proc. ADCHEM
97
,
78
83
.
https://doi.org/10.1016/S1474-6670(17)43143-0
.
Giordanengo
T.
Charpentier
J. P.
Roger
J. M.
Roussel
S.
Brancheriau
L.
Chaix
G.
Baillères
H.
2008
Correction of moisture effects on near infrared calibration for the analysis of phenol content in eucalyptus wood extracts
.
Ann. For. Sci.
65
(
8
),
803
.
https://doi.org/10.1051/forest:2008065
.
Gujer
W.
Zehnder
J. B.
1983
Conversion processes in anaerobic digestion
.
Water Sci. Technol
15
(
8–9
),
127
167
.
https://doi.org/10.2166/wst.1983.0164
.
He
X. B.
Yang
Y. P.
2008
Variable MWPCA for adaptive process monitoring
.
Ind. Eng. Chem. Res.
47
(
2
),
419
427
.
https://doi.org/10.1021/ie070712z
.
Holm-Nielsen
J. B.
Andree
H.
Lindorfer
H.
Esbensen
K. H.
2007
Transflexive embedded near infrared monitoring for key process intermediates in anaerobic digestion/biogas production
.
J. Near Infrared Spectrosc.
15
(
2
),
123
135
.
https://doi.org/10.1255/jnirs.719
.
Jackson
J. E.
1991
A Users Guide to Principal Components
.
John Wiley & Sons Inc
.
https://doi.org/10.1002/0471725331.scard
.
Jackson
J.
Mudholkar
G.
1979
Control procedures for residuals associated principal component analysis
.
Technometrics
21
(
3
),
341
349
.
https://doi.org/10.1080/00401706.1979.10489779
.
Jacobi
H. F.
Moschner
C. R.
Hartung
E.
2009
Use of near infrared spectroscopy in monitoring of volatile fatty acids in anaerobic digestion
.
Water Sci. Technol.
60
(
2
),
339
346
.
https://doi.org/10.2166/wst.2009.345
.
Jeng
J. C.
2010
Adaptive process monitoring using efficient recursive PCA and moving window PCA algorithms
.
J. Taiwan Inst. Chem. Eng.
41
(
4
),
475
481
.
https://doi.org/10.1016/j.jtice.2010.03.015
.
Lalman
J.
Bagley
D. M.
2002
Effects of C18 long chain fatty acids on glucose, butyrate and hydrogen degradation
.
Water Res.
36
(
13
),
3307
3313
.
https://doi.org/10.1016/S0043-1354(02)00014-3
.
Li
W.
Yue
H.
Valle-Cervantes
S.
Qin
S.
2000
Recursive PCA for adaptive process monitoring
.
J. Process Control
10
(
5
),
471
486
.
https://doi.org/10.1016/S0959-1524(00)00022-6
.
Lu
N.
Yao
Y.
Gao
F.
Wang
F.
2005
Two-dimensional dynamic PCA for batch process monitoring
.
AIChE J.
51
(
12
),
3300
3304
.
https://doi.org/10.1002/aic.10568
.
Lucas
T.
Le Ray
D.
Peu
P.
Wagner
M.
Picard
S.
2007
A new method for continuous assessment of CO2 released from dough baked in ventilated ovens
.
J. Food Eng.
81
(
1
),
1
11
.
https://doi.org/10.1016/j.jfoodeng.2006.09.024
.
MacGregor
J. F.
Kourti
T.
1995
Statistical process control of multivariate process control
.
Control Eng. Pract.
3
(
3
),
403
414
.
https://doi.org/10.1016/0967-0661(95)00014-L
.
Madsen
M.
Holm-Nielsen
J. B.
Esbensen
K. H.
2011
Monitoring of anaerobic digestion processes: a review perspective
.
Renew. Sust. Energ. Rev.
15
(
6
),
3141
3155
.
https://doi.org/10.1016/j.rser.2011.04.026
.
Qin
S. J.
2003
Statistical process monitoring: basics and beyond
.
J. Chemom.
17
(
8–9
),
480
502
.
https://doi.org/10.1002/cem.800
.
Reed
J. P.
Devlin
D.
Esteves
S. R. R.
Dinsdale
R.
Guwy
A. J.
2013
Integration of NIRS and PCA techniques for the process monitoring of sewage sludge anaerobic digester
.
Bioresour. Technol.
133
,
398
404
.
https://doi.org/10.1016/j.biortech.2013.01.083
.
Roggo
Y.
Chalus
P.
Maurer
L.
Lema-Martinez
C.
Edmond
A.
Jent
N.
2007
A review of near infrared spectroscopy and chemometrics in pharmaceutical technologies
.
J. Pharm. Biomed. Anal.
44
(
3
),
683
700
.
https://doi.org/10.1016/j.jpba.2007.03.023
.
Savitzky
A.
Golay
M. J.
1964
Smoothing and differentiation of data by simplified least squares procedures
.
Anal. Chem.
36
(
8
),
1627
1639
.
https://doi.org/10.1021/ac60214a047
.
Schmitt
E.
Rato
T.
De Ketelaere
B.
Reis
M.
Hubert
M.
2016
Parameter selection guidelines for adaptive PCA-based control charts
.
J. Chemom.
30
,
163
176
.
https://doi.org/10.1002/cem.2783
.
Spanjers
H.
van Lier
J. B.
2006
Instrumentation in anaerobic treatment – research and practice
.
Water Sci. Technol.
53
(
4–5
),
63
76
.
https://doi.org/10.2166/wst.2006.111
.
Taris
A.
Hansen
T. B.
Rong
B.-G.
Grosso
M.
Qu
H.
2017
Detection of nucleation during cooling crystallization through moving window PCA applied to in situ infrared data
.
Org. Process Res. Dev.
21
(
7
),
966
975
.
https://doi.org/10.1021/acs.oprd.7b00076
.
Ward
A. J.
Hobbs
P. J.
Holliman
P. J.
Jones
D. L.
2011
Evaluation of near infrared spectroscopy and software sensor methods for determination of total alkalinity in anaerobic digesters
.
Bioresour. Technol.
102
(
5
),
4083
4090
.
https://doi.org/10.1016/j.biortech.2010.12.046
.
Xie
L.
Ye
X.
Liu
D.
Ying
Y.
2009
Quantification of glucose, fructose and sucrose in bayberry juice by NIR and PLS
.
Food Chem.
114
(
3
),
1135
1140
.
https://doi.org/10.1016/j.foodchem.2008.10.076
.
Yamanaka
O.
Yoshizawa
N.
Hiraoka
Y.
Sano
K.
Hamada
T.
2012
A monitoring technique using multivariate statistical process control method for performance improvement with application to wastewater treatment plant operation
.
IFAC Proc.
45
(
16
),
619
624
.
https://doi.org/10.3182/20120711-3-BE-2027.00223
.
Zeaiter
M.
Roger
J. M.
Bellon-Maurel
V.
2005
Robustness of models developed by multivariate calibration. Part II: the influence of pre-processing methods
.
TrAC, Trends Anal. Chem.
24
,
437
445
.
https://doi.org/10.1016/j.trac.2004.11.023
.
Zhang
H.
Wang
Z.
Liu
O.
2015
Development and validation of a GC-FID method for quantitative analysis of oleic acid and related fatty acids
.
J. Pharm. Anal.
5
(
4
),
223
230
.
https://doi.org/10.1016/j.jpha.2015.01.005
.

Supplementary data