3.1.1. Load-Carrying Capacity and Failure Modes

The failure modes of the MPSs are similar. The development of cracks was analyzed based on the MP25H16 specimen as an example. Firstly, there were no changes in the surface of these specimens at the initial loading stage. As the loading continued, the fragments at the shear plane began to exfoliate. After that, the initial cracks appeared on the shear plane and several small cracks appeared instantaneously. Herein, it should be noted that the crack width of 0.05 mm is adopted as the criterion of visible initial cracking [28]. With the increase in loading, fine cracks further spread, connected and penetrated to form a crack zone along the shear failure surface, and the fibers between the crack zones were gradually pulled out or pulled off. Finally, with the further increase in the load, along with a huge sound, specimens were sheared and damaged. The condition of the crack development is shown in Figure 8. In order to verify whether the UHPC specimens still possess the bearing capacity after the main crack occurs, the test machine was restarted to continue loading. It turns out the load could still reach 1/2~2/3 of the ultimate load. On the basis of the testing results, it can be seen that there are two main crack modes in the failure modes for these specimens, namely the single main crack (Figure 9a) and multiple diagonal cracks (Figure 9b), respectively.

Since the failure characteristics of the specimens are similar, we took the MP25H13 specimen as an example to describe the shear load–slip relationship of the MPSs. The shear load–slip curve of MP25H13 specimen is shown in Figure 10. The results show that the direct shear cracking failure characteristics of UHPC can be divided into three domains. The first domain is a linear development stage in which the relationship between load and slip is basically linear (ideal status). The second domain is a crack growth stage. With the initial stiffness degenerating, the slope of the load–slip curve begins to decline. After the micro-crack occurs, the stiffness gradually decreases, and then the load enters a linear deviation domain. As a whole, because the bridging action of the steel fibers suppresses the crack expansion, the load increases in the alternation as the slip increases. The third domain is the shear failure stage. At this stage, large amounts of steel fibers are broken or pulled out from the matrix, and the specimens are rapidly destroyed. *Materials* **2017**, *10*, x FOR PEER REVIEW 9 of 17

**Figure 8.** Failure modes and the crack propagation process. **Figure 8.** Failure modes and the crack propagation process.

(**a**) (**b**)

**Figure 9.** Two typical crack forms: (**a**) single main crack and (**b**) multiple diagonal cracks.

from the matrix, and the specimens are rapidly destroyed.

Single main crack

Since the failure characteristics of the specimens are similar, we took the MP25H13 specimen as an example to describe the shear load–slip relationship of the MPSs. The shear load–slip curve of MP25H13 specimen is shown in Figure 10. The results show that the direct shear cracking failure characteristics of UHPC can be divided into three domains. The first domain is a linear development stage in which the relationship between load and slip is basically linear (ideal status). The second domain is a crack growth stage. With the initial stiffness degenerating, the slope of the load–slip curve begins to decline. After the micro-crack occurs, the stiffness gradually decreases, and then the load enters a linear deviation domain. As a whole, because the bridging action of the steel fibers suppresses the crack expansion, the load increases in the alternation as the slip increases. The third domain is the shear failure stage. At this stage, large amounts of steel fibers are broken or pulled out

**Figure 8.** Failure modes and the crack propagation process.

(**a**) (**b**) **Figure 9.** Two typical crack forms: (**a**) single main crack and (**b**) multiple diagonal cracks. **Figure 9.** Two typical crack forms: (**a**) single main crack and (**b**) multiple diagonal cracks. *Materials* **2017**, *10*, x FOR PEER REVIEW 10 of 17

**Figure 10.** UHPC shear load–slip curve. **Figure 10.** UHPC shear load–slip curve.

## 3.1.2. Shear Strength Results

MP25H16

MP30H16

MP25S13

MP25H13

286.4 515.6

291.0 469.8

284.2 454.3

3.1.2. Shear Strength Results For the MPSs, the baroclinic bar and the steel fibers passing through the shear plane act together to form a truss to resist the shear force along the shear plane. The effect of steel fibers on the shear For the MPSs, the baroclinic bar and the steel fibers passing through the shear plane act together to form a truss to resist the shear force along the shear plane. The effect of steel fibers on the shear properties of UHPC was analyzed, and the test results of fifteen MPSs are listed in Table 5.

properties of UHPC was analyzed, and the test results of fifteen MPSs are listed in Table 5.

254.1 420.8 21.04 1.656 1.604

296.4 494.8 24.74 1.669 1.669

261.3 481.5 24.08 1.843 1.630

21.53 0.059

24.72 0.035

23.13 0.028

22.98 0.035

1.800

1.614

1.599

430.6

494.5

462.7

459.8

273.8 465.5 23.28 1.700

25.78

307.3 473.0 23.65 1.539

291.5 444.6 22.23 1.525

22.72

306.2 443.5 22.18 1.447

According to Table 5, the maximum shear strength of UHPC can reach 24.72 MPa (MP30H16), which is 1.37 times the shear strength of MP20H16. It can be seen that the direct shear strength of UHPC increases with the increase of the fiber volume fractions. In addition, the shear strength of MP25H16 is 19.6% higher than that of MP20H16, while the shear strength of MP30H16 is only 14.8% higher than that of MP25H16. This shows that the shear strength of UHPC increases more obvious at the initial phase of the fiber volume fraction increment. According to the ratio of ultimate load to initial crack load of the specimens and shear slip relationship along the shear plane under different fiber volume fraction (Figure 11), when the fiber volume fraction is fewer, cracking failure characteristics of UHPC are close to brittle failure. Based on the above analysis, it can be considered

**Table 5.** *Cont.*

According to Table 5, the maximum shear strength of UHPC can reach 24.72 MPa (MP30H16), which is 1.37 times the shear strength of MP20H16. It can be seen that the direct shear strength of UHPC increases with the increase of the fiber volume fractions. In addition, the shear strength of MP25H16 is 19.6% higher than that of MP20H16, while the shear strength of MP30H16 is only 14.8% higher than that of MP25H16. This shows that the shear strength of UHPC increases more obvious at the initial phase of the fiber volume fraction increment. According to the ratio of ultimate load to initial crack load of the specimens and shear slip relationship along the shear plane under different fiber volume fraction (Figure 11), when the fiber volume fraction is fewer, cracking failure characteristics of UHPC are close to brittle failure. Based on the above analysis, it can be considered that the ductile characteristic and direct shear strength of the monolithic placement specimens with appropriate steel fibers can be improved in direct shear load. *Materials* **2017**, *10*, x FOR PEER REVIEW 11 of 17 that the ductile characteristic and direct shear strength of the monolithic placement specimens with appropriate steel fibers can be improved in direct shear load.

**Figure 11.** Shear load–slip curves of UHPC specimens with a different fiber volume fraction: (**a**) MP20H16, (**b**) MP25H16 and (**c**) MP30H16. **Figure 11.** Shear load–slip curves of UHPC specimens with a different fiber volume fraction: (**a**) MP20H16, (**b**) MP25H16 and (**c**) MP30H16.

As shown in Figure 11, the initial stiffness of the specimens (Stage 1) is reduced with the increase in steel fibers. This can be explained as follows. On the one hand, the uniformity of fiber mixing will

is predominantly dependent on the overall modulus of the matrix. Only after micro-cracks occur, the shear force will be balanced by the steel fibers and the matrix. Thus, with the increase in steel fiber content, the overall modulus of the matrix itself is reduced, leading to a reduced shear stiffness at Stage 1. With the further increase in load, micro-cracks occur, and the influence of the internal defects is gradually eliminated. The subsequent stiffness (Stage 2) mainly depends on the comprehensive modulus of the steel fiber and matrix. It should be noted that these phenomena are only reflected

Table 5 also shows the effect of different types of fibers on the UHPC shear strength. Under the same volume fraction of fibers (2.5%), the shear strengths of MP25S13, MP25H13 and MP25H16 are almost at the same level. It shows that the shear strength is not significantly influenced by the shape of the steel fibers. In the case of a certain fiber length, whether or not the hooked-type steel fiber is used has little effect on the shear strength of UHPC. However, when using the same shape fibers, increasing the fiber length will slightly reduce the shear strength of UHPC. This is because the UHPC specimens with short fibers possess more fibers at the same volume fraction. From the previous analysis, increasing the fiber numbers is an important way to enhance the compressive strength of UHPC. In addition, according to the ratio of ultimate load to the initial crack load of the specimens

when the fiber content is more than 2.5%.

As shown in Figure 11, the initial stiffness of the specimens (Stage 1) is reduced with the increase in steel fibers. This can be explained as follows. On the one hand, the uniformity of fiber mixing will be affected with the increase of steel fiber content, and the probability of internal defects (voids) of the UHPC will also be increased. On the other hand, the shear stiffness at the early stage of loading is predominantly dependent on the overall modulus of the matrix. Only after micro-cracks occur, the shear force will be balanced by the steel fibers and the matrix. Thus, with the increase in steel fiber content, the overall modulus of the matrix itself is reduced, leading to a reduced shear stiffness at Stage 1. With the further increase in load, micro-cracks occur, and the influence of the internal defects is gradually eliminated. The subsequent stiffness (Stage 2) mainly depends on the comprehensive modulus of the steel fiber and matrix. It should be noted that these phenomena are only reflected when the fiber content is more than 2.5%.

Table 5 also shows the effect of different types of fibers on the UHPC shear strength. Under the same volume fraction of fibers (2.5%), the shear strengths of MP25S13, MP25H13 and MP25H16 are almost at the same level. It shows that the shear strength is not significantly influenced by the shape of the steel fibers. In the case of a certain fiber length, whether or not the hooked-type steel fiber is used has little effect on the shear strength of UHPC. However, when using the same shape fibers, increasing the fiber length will slightly reduce the shear strength of UHPC. This is because the UHPC specimens with short fibers possess more fibers at the same volume fraction. From the previous analysis, increasing the fiber numbers is an important way to enhance the compressive strength of UHPC. In addition, according to the ratio of ultimate load to the initial crack load of the specimens and shear load–slip curves under different fiber types (Figure 12), the use of hooked-type fibers has little effect on the ductile failure characteristics of UHPC specimens. *Materials* **2017**, *10*, x FOR PEER REVIEW 12 of 17 and shear load–slip curves under different fiber types (Figure 12), the use of hooked-type fibers has little effect on the ductile failure characteristics of UHPC specimens.

**Figure 12.** Shear load–slip curves of UHPC specimens with different fiber types: (**a**) MP25S13 and (**b**) MP25H13. **Figure 12.** Shear load–slip curves of UHPC specimens with different fiber types: (**a**) MP25S13 and (**b**) MP25H13.

### *3.2. Test Results and Analysis of the WJTSs 3.2. Test Results and Analysis of the WJTSs*

steel fibers.

3.2.2. Shear Strength Results

### 3.2.1. Load-Carrying Capacity and Failure Modes

3.2.1. Load-Carrying Capacity and Failure Modes The shear failure characteristic of the WJTSs is a typical brittle failure, which occurs abruptly and without obvious symptoms. Similarly, the shear failure of the specimens is accompanied by a loud noise, especially for the WJ25H16, while the sound of the specimens without steel fibers will not be harsh. From the failure of the bond interface of the specimens, the failure surface of the specimens without fibers is relatively smooth (Figure 13a). But parts of the poured UHPC are embedded together successively, and the peeling phenomenon between the UHPC blocks interfaces is not obvious. The bonding effect of steel fiber reinforced specimens (especially for the WJ25H16) is The shear failure characteristic of the WJTSs is a typical brittle failure, which occurs abruptly and without obvious symptoms. Similarly, the shear failure of the specimens is accompanied by a loud noise, especially for the WJ25H16, while the sound of the specimens without steel fibers will not be harsh. From the failure of the bond interface of the specimens, the failure surface of the specimens without fibers is relatively smooth (Figure 13a). But parts of the poured UHPC are embedded together successively, and the peeling phenomenon between the UHPC blocks interfaces is not obvious. The bonding effect of steel fiber reinforced specimens (especially for the WJ25H16) is excellent due to the bridging effect of steel fibers (Figure 13b), which is similar to the MPSs. When the

steel fibers can play an important role in enhancing the shear plane strength for the WJTSs.

(**a**) (**b**)

**Figure 13.** Direct shear failure of cast-in-place specimens: (**a**) without steel fibers and (**b**) doped with

excellent due to the bridging effect of steel fibers (Figure 13b), which is similar to the MPSs. When the specimens are damaged, obvious fiber pullout marks and UHPC fragments can be seen. Thus,

MP25H13.

*3.2. Test Results and Analysis of the WJTSs*

3.2.1. Load-Carrying Capacity and Failure Modes

specimens are damaged, obvious fiber pullout marks and UHPC fragments can be seen. Thus, steel fibers can play an important role in enhancing the shear plane strength for the WJTSs. excellent due to the bridging effect of steel fibers (Figure 13b), which is similar to the MPSs. When the specimens are damaged, obvious fiber pullout marks and UHPC fragments can be seen. Thus, steel fibers can play an important role in enhancing the shear plane strength for the WJTSs.

obvious. The bonding effect of steel fiber reinforced specimens (especially for the WJ25H16) is

*Materials* **2017**, *10*, x FOR PEER REVIEW 12 of 17

and shear load–slip curves under different fiber types (Figure 12), the use of hooked-type fibers has

(**a**) (**b**)

**Figure 12.** Shear load–slip curves of UHPC specimens with different fiber types: (**a**) MP25S13 and (**b**)

The shear failure characteristic of the WJTSs is a typical brittle failure, which occurs abruptly and without obvious symptoms. Similarly, the shear failure of the specimens is accompanied by a loud noise, especially for the WJ25H16, while the sound of the specimens without steel fibers will not be harsh. From the failure of the bond interface of the specimens, the failure surface of the specimens without fibers is relatively smooth (Figure 13a). But parts of the poured UHPC are embedded

little effect on the ductile failure characteristics of UHPC specimens.

**Figure 13.** Direct shear failure of cast-in-place specimens: (**a**) without steel fibers and (**b**) doped with steel fibers. **Figure 13.** Direct shear failure of cast-in-place specimens: (**a**) without steel fibers and (**b**) doped with steel fibers.


**4. Association of LUTD with UTI and VUR**

are often believed to originate from behavioural issues. These behavioural traits may evolve from adverse events that occur around or after the time of toilet training and/or personal stresses [13]. Severe emotional stressors, such as sexual abuse mainly in girls, has been asso‐ ciated with DV, and should be considered in a child, especially a girl, who presents with

Detrusor overactivity as a component of DV may represent a persistence of the normal in‐ fant voiding pattern after toilet training. It may be that a mild delay in the maturation of the central nervous system disrupts the ability of these children to learn true voluntary control over the micturition reflex [5, 14]. An alternative theory to the development of the urge syn‐ drome hypothesizes that detrusor overactivity is caused by the transient obstruction that oc‐ curs with DV [15]. In support of the belief that the roots of DV may be grounded in behavioural issues or central nervous system developmental delays, an association has been demonstrated between DV and attention deficit hyperactivity disorder (ADHD). Higher rates of enuresis, urinary incontinence, constipation, and other voiding symptoms have been

As a result of studies of DV diagnosed in infancy, a congenital or genetic component to the disorder has also been postulated. Small series of infants with signs and symptoms consistent with NNNB syndromes have been reported [18, 19]. Furthermore, DV has been linked to the Ochoa syndrome, a genetic disorder with an autosomal recessive in‐ heritance pattern. The gene locus at chromosome 10q23-q24 is identified as the defective gene in Ochoa syndrome. It is postulated to be the possible gene locus of the NNNB de‐ scribed by Hinman and Allan [20]. This information casts doubt on the commonly held belief that disturbance of behaviour is the sole cause of NNNB because these findings

At least some patients with DV represent occult neurogenic problems that will manifest provided these patients are followed longitudinally [21]. In all individuals with unex‐ plained severe DV, search for an unidentified neurological lesion must be made. A sub‐ tle neurological insult could present as DV. Such lesions may or may not be detectable with current imaging technologies. Routine magnetic resonance imaging in children with lower urinary tract problems without overt neurological signs and symptoms has a low yield of 7.5% but this may be improved by targeting children with abnormal cutaneous findings [22]. Tethered cord syndrome may be identified in some patients with subtle neurological signs and symptoms. Classical tethered cord has been diagnosed on the ba‐ sis of pathologically elongated conus or conus that lies below the L2 level. However, an‐ ecdotal successful outcomes following surgical division of the filum in children with apparently normally located cords suggest that the cord may sometimes be abnormally stretched without being at an abnormal location [23]. Such subtle lesions might also ex‐ plain the occasional presentation in infancy when the problem has its onset before toilet training has commenced [19]. Surgical division of the filum terminale in such patients is

Another possible indirect evidence for an unidentified neurological lesion in these pa‐ tients is the association of a peculiar facial expression in some children with DV [25].

controversial but may yield improvements in bladder dysfunction [24].

new onset DV and no other identifiable etiologic factor [4, 5].

described in children with ADHD [16, 17].

92 Recent Advances in the Field of Urinary Tract Infections

are present at or near the time of birth.

There is not sufficient reported data on the relationship between the type of LUTD with UTI and VUR. The earliest studies of LUTD mostly dealt with either OAB or DV. This is a some‐ what artificial distinction, as these conditions are often combined and sometimes difficult to separate. For that reasons this association will be presented in LUTD children as a whole.

The association between LUTD and UTI has been established, though the causal relationship is not clear. Recurrent UTI has been shown in any studies to be higher in VUR patients with LUTD than in those without such dysfunction [3, 27, 28]. It has been demonstrated that ade‐ quate management of LUTD not only decreases the rate of UTI but also increases resolution of the VUR [27, 29]. Traditionally, recurrent UTI and pyelonephritis have been recognized as potential causes of permanent renal damage [30]. Current opinion is that VUR alone is not sufficient to cause UTI or renal damage. Holland et al. reported that girls with primary VUR followed up for 10 years, with no recurrent UTI, did not develop renal scars [31]. Linshaw showed that VUR does not threaten the kidney as long as UTI is promptly treated [32]. The results of these studies suggest that the association between VUR and UTI is necessary for renal damage to occur, mainly in situations of low detrusor pressure. However, VUR may predispose invasion of the renal parenchyma by bacteria. It has been reported that LUTD is an important risk factor for VUR and renal damage [33]. In addition, current studies have showed that increased intravesical pressure associated with LUTD is a primary factor for in‐ ducing reflux and renal damage [12, 34]. In patients who had UTI the presence of reflux in‐ creased the rate of renal damage [12].

VUR in LUTD is theorized to be not the result of a short mucosal tunnel, but to be the consequence of the high filling and voiding pressures. In patients with LUTD, uninhibit‐ ed detrusor contractions and voluntary constriction of the sphincter, causing a functional obstruction, increase the intravesical pressure. Increased intravesical pressure can pro‐ mote VUR through a possible marginal competence in the valve mechanism [34, 35]. Re‐ cent studies have shown that the prevalence of VUR among children with idiopathic LUTD is between 14% and 46% [11, 12]. Some reports have emphasized bilateral reflux associated with LUTD [12, 36].

pertrophy can lead to an increase in risk factors that predispose the urinary tract to more UTIs and more severe VUR including anatomic bladder abnormalities, constipation, in‐ creased residual urine, higher voiding pressures, and increased urethral bacteria coloniza‐ tion secondary to turbulent flow [41, 42]. Turbulent flow is of particular importance as eddy currents formed by nonlaminar flow leads to reflux of periurethral bacteria to the proximal

Developments Regarding Dysfunctional Voiding and Urinary Tract Infections in Children

http://dx.doi.org/10.5772/53569

95

The pelvic floor musculature is closely related to bowel and bladder. Isolated contraction of the pelvic floor musculature was found to lead to spontaneous contraction of the bladder. The cycle of pelvic floor dysfunction contributing to recurrent UTIs and VUR can cause wor‐ sening bladder and bowel symptoms, and increased pelvic floor dysfunction with functional obstruction. Eventually this can lead to a hypertrophied, small-capacity bladder with highpressure voiding that will lead to renal damage. The integral relationship between the pelvic floor activity with UTIs and VUR can be used as a model for developing treatment strategies in affected children to address the cognitive and behavioural aspects of DV and prevent irre‐

Clinical symptoms may vary from mild incontinence to severe disorders with endpoints of irreversible bladder dysfunction with VUR, UTI and resulting nephropathy [8]. Children with DV voiding often present with urinary incontinence both during the day as well as at night. They may have urinary frequency, urgency, urge incontinence or nocturnal enuresis. Such storage symptoms may result from associated detrusor overactivity, urinary infection or reduced bladder capacity consequent to large residual urine and may be aggravated by constipation or behavioural disorders. A distinctive facial expression may be noted in some

Diagnosis relies heavily on a good history and physical examination, but also includes radi‐ ologic and urodynamic evaluation. The history should be directed towards the identification of children with neurologic or anatomic causes of their symptoms, and then distinguish be‐ tween which form or pattern of voiding dysfunction is present. Components of this history include maternal medical issues, perinatal history, developmental milestones, scholastic per‐ formance, behavioural history, specifics around toilet training, patterns of voiding and bow‐ el movements, history of UTI, and family history of voiding dysfunction. The use of a 3-day voiding diary is often helpful to identify the frequency of voiding, voided volumes, and tim‐

Clinical examination must include an assessment of higher mental functions and their ageappropriateness and basic neurological evaluation. On physical examination, attention should be paid to the back and lower spine for cutaneous manifestations of an occult spinal dysraphism and/or sacral agenesis. Neurologic examination should include assessment of lower extremity function, rectal tone, perineal/anal sensation, and intactness of bulbocaver‐

urethra and bladder (milk back phenomenon), causing recurrent infections.

versible damage to the upper urinary tract [40].

**5. Diagnosis and evaluation**

of these patients [8].

ing of incontinent episodes [9].

Van Gool et al. addressed the relationship between DV and reflux for the first time with a retrospective questionnaire in 1992 [37]. The prevalence of uncoordination between the de‐ trusor and urethral sphincter approached 18% and included voiding pattern abnormalities such as urge syndrome, staccato voiding, fractionated and incomplete voiding, and voiding postponement. They also found that, in those children who had spontaneous resolution of their VUR, there was a lower prevalence of DV.

It is important to carefully assess all children with reflux for subtle signs of DV. Chil‐ dren with DV are more likely to have recurrent UTIs, have mild bilateral reflux with less spontaneous resolution, and are less likely to have success with surgical management. The treatment of DV in such children can improve the chances of spontaneous resolution of the reflux and may also reduce recurrent urinary infection. Koff et al. [3] reported on their series of children with VUR, who either resolved spontaneously or were surgically treated, and found that DES was observed in 43% of children with primary VUR and in 77% of a subset of these children who had breakthrough UTI. The presence of DES pat‐ terns was associated with a longer time for spontaneous resolution of low-grade reflux and with unsuccessful surgical outcomes. Of children in the surgically treated group, on‐ ly those with DES developed recurrent and/or contralateral reflux. Children with untreat‐ ed DV undergoing ureteral re-implantation may be at a higher risk for developing recurrent reflux or a new bladder diverticulum [38].

Contrary to these results, Chen and colleagues reported that VUR and UTI are not inde‐ pendent of DV [39]. Their findings are new observations contradicting the previous belief that both UTI and VUR are independently associated with DES. Chen et al. performed a multivariate analysis on 2759 paediatric urology patients, further examining the relationship between DES, VUR, and UTI. Their data demonstrated a higher rate of DES in girls than in boys: 43.7% compared with 23.8%. This group also found that there was no difference in the presence of DES in patients with unilateral and bilateral VUR. Surprisingly, they observed no association of VUR or UTI individually with DES but rather DES was only noted when both of these issues were present. Although this large-scale multivariate analysis is a more statistically powerful study than previous smaller retrospective analyses, the data could be potentially skewed in that all of the patients were recruited from a paediatric urology popu‐ lation, which is not representative of the general paediatric population. Patients with known reflux could potentially be protected from DES by prophylactic antibiotics or continued clin‐ ical follow-up in a specialty setting. On the other hand, it is suggested that there may be two types of reflux: one that is primary or congenital in nature and another that is secondary and in part due to DES and UTI. This study is in line with the neuroplastic theory that postulates that hypertrophy of the bladder and bowel musculature is caused by trophic factors re‐ leased during pelvic floor hyperactivity secondary to central nervous disturbances [40]. Hy‐ pertrophy can lead to an increase in risk factors that predispose the urinary tract to more UTIs and more severe VUR including anatomic bladder abnormalities, constipation, in‐ creased residual urine, higher voiding pressures, and increased urethral bacteria coloniza‐ tion secondary to turbulent flow [41, 42]. Turbulent flow is of particular importance as eddy currents formed by nonlaminar flow leads to reflux of periurethral bacteria to the proximal urethra and bladder (milk back phenomenon), causing recurrent infections.

The pelvic floor musculature is closely related to bowel and bladder. Isolated contraction of the pelvic floor musculature was found to lead to spontaneous contraction of the bladder. The cycle of pelvic floor dysfunction contributing to recurrent UTIs and VUR can cause wor‐ sening bladder and bowel symptoms, and increased pelvic floor dysfunction with functional obstruction. Eventually this can lead to a hypertrophied, small-capacity bladder with highpressure voiding that will lead to renal damage. The integral relationship between the pelvic floor activity with UTIs and VUR can be used as a model for developing treatment strategies in affected children to address the cognitive and behavioural aspects of DV and prevent irre‐ versible damage to the upper urinary tract [40].


3.1.1. Time-Dependent Density Functional Theory—Fundamentals

The standard approach of DFT is based on the existence of a one-to-one mapping between the exact ground-state density and the external potential applied to a quantum system [48]. Furthermore, one can extend this mapping between two systems, the many-body interacting "real" system and a many-body non-interacting "Kohn-Sham" system [39] where we introduce an effective external potential (dubbed KS potential). This second mapping allows for calculating quantities, i.e., the ground-state density and energy, belonging to the real system by using the KS system, bringing about a large numerical simplification. DFT is therefore the actual method of choice for calculating electronic structure of materials as well as the energies of complex atoms and molecules [42]. The KS mapping was later extended to the dynamics of the single-particle density by Runge and Gross [55]. The existence of this second mapping lays the foundations for performing, e.g., numerical spectroscopy with so-called ab-initio methods, i.e., without—in principle—any fitting parameter [41]. TDDFT can theoretically be used for studying electrical transport (see for example Ref. [56–60]), but the information it can provide is only partial [61,62].

It appears natural therefore, from the Runge and Gross' theorem (RG), to establish a mapping between the single-particle current density and the external (vector) potential applied to the real system [63]. In principle, this Time-Dependent Current-Density Functional Theory (TDCDFT) should be the workhorse for charge transport calculations, but the lack of suitable approximations for the KS's potential renders the theory of little use at the moment. One the other hand, addressing thermal transport phenomena within DFT would require instead a non-equilibrium theory based on either a local temperature or a local energy density. Recently, a functional theory based on the excess energy density as the basic variable has been presented [64,65] which is suited to study thermoelectrical phenomena in the static and time-dependent case. In the following, we will look only at standard TDDFT.

The RG theorem proves that for a general time-dependent Hamiltonian *H* that describes the dynamics of a many-body system,

$$H(t) = T + V\_{\text{int}} + V\_{\text{ext}}(t) \tag{35}$$

there exists a one-to-one mapping between the single-particle density *n*(**r**, *t*) and the external potential *V*ext(**r**, *t*), given the initial conditions. Here, *T* is the kinetic energy, *V*int the particle-particle interaction, and *V*ext the time-dependent external potential. The proof of the theorem assumes that both density and potential can be expanded in series around the initial time: a detailed proof and the physical assumption on which this is based can be found in [40,41]. It has since been shown that there exists a system of non-interacting particles whose single-particle density is identical at each time to *n*(**r**, *t*) of the interacting system [66]. The existence of this mapping then allows investigation of the dynamics of a system by looking at its non-interacting Doppelgänger. Clearly, this brought about the same reduction of computational requirements as the standard DFT. However, as was the case with DFT, we are exchanging some of the physical information about the system to reduce the dimensionality of the problem. It should then come as no-surprise that TDDFT cannot for example reproduce the exact single-particle current density **j**(**r**, *t*) since some of its contributions are not derivable from the knowledge of the density alone [67]. Indeed, if starting from the continuity equation,

$$
\partial\_l n(\mathbf{r}, t) = -\nabla \cdot \mathbf{j}(\mathbf{r}, t), \tag{36}
$$

we need to conclude that given the density we can determine only the longitudinal part of **j**(**r**, *t*) and any term **j***T* such that **j**(**r**, *t*) = **j***L*(**r**, *t*) + **j***T*(**r**, *t*) and ∇ · **j***T* ≡ 0 cannot be obtained from Equation (36). However, some information about the total current can still be obtained from the continuity equation [68].

As with the static DFT, we study the KS Hamiltonian *H*KS,

$$H\_{\rm KS} = T\_{\rm S} + V\_{\rm KS} \tag{37}$$

where *TS* = ∑*i p*2 *i* /2*m* is the kinetic energy of the non-interacting particles of mass *m* and momenta {**p***i*}, and *VKS* = *V*ext + *V*Hxc. The potential *V*Hxc is the sum of the Hartree and exchange–correlation (xc) potentials that replace the particle–particle interaction. By the RG theorem *V*Hxc is a functional of the density only, *V*Hxc = *V*Hxc[*n*](**r**, *t*), and the dynamics of the single-particle density of the real system evolving with Hamiltonian *H* is identical to the dynamics of the single-particle density evolving with *H*KS. Notice that in principle, given the particle–particle interaction of the original many-body problem, *V*Hxc is a universal functional and does not depend on the external potential *V*ext. This implies for example that if we are interested in the dynamics of an electronic system, *V*Hxc is the same either if we are studying an atom, a molecule, a slab, or a bulk, i.e., this potential is transferrable to whatever system we want to study. This high transferability makes finding reliable approximation for *V*Hxc difficult. Indeed, the universal *V*Hxc for the electrons contains the physical information needed to go from the standard weakly interacting Landau's Fermi liquid, to the strongly correlated Wigner's crystal or superconducting phases. Notice that for the electrons, a large part of *V*Hxc is given by the Hartree's (mean-field) interaction, and only a relatively small contribution of the overall energy determines all these interesting phases of matter. Therefore, standard approaches to build some approximation to *V*Hxc are based on interpolating the numerical solution of the many-body problem with some known high- and low-density limits [69,70]. Many of these approximations are static. To apply them to the dynamics, a common method is the so-called adiabatic local density approximation (ALDA). In this approximation, we take a static *V*Hxc[*n*] and replace the static density with the instantaneous density *n*(**r**, *t*). Indeed, this neglects the history of the system—while we generally expect that the xc potential to be history dependent. For example, one can show that in this approximation there is not relaxation induced by particle–particle interaction, in contrast with observation [71,72].


**EMERGING POLLUTANTS IN THE ENVIRONMENT - CURRENT AND FURTHER IMPLICATIONS**

Edited by **Marcelo L. Larramendy** and **Sonia Soloneski**

#### **Emerging Pollutants in the Environment - Current and Further Implications**

http://dx.doi.org/10.5772/59332 Edited by Marcelo L. Larramendy and Sonia Soloneski

#### **Contributors**

Vladimira Hanzelova, Rafael Grossi Grossi Botelho, Sérgio Monteiro, Valdemar Tornisielo, Ismail Saadoun, Jacek Namiesnik, Błażej Kudłak, Natalia Szczepańska, Katarzyna Owczarek, Zofia Mazerska, Sajal Ray, Soumalya Mukherjee, Niladri Sekhar Bhunia, Anindya Sundar Bhunia, Mitali Ray, Mohamed Faycal Atitar, Ralf Dillert, Detlef Bahnemann, Hamza Belhadj, Juan Jose Alava, Frank Gobas, K. Paola Calle, Monica Riofrio, Mercy J. Borbor-Cordova, Nastenka Calle, Michael G. Ikonomou, Mandy R. R. McDougall

#### **© The Editor(s) and the Author(s) 2015**

The moral rights of the and the author(s) have been asserted.

All rights to the book as a whole are reserved by INTECH. The book as a whole (compilation) cannot be reproduced, distributed or used for commercial or non-commercial purposes without INTECH's written permission. Enquiries concerning the use of the book should be directed to INTECH rights and permissions department (permissions@intechopen.com).

Violations are liable to prosecution under the governing Copyright Law.

Individual chapters of this publication are distributed under the terms of the Creative Commons Attribution 3.0 Unported License which permits commercial use, distribution and reproduction of the individual chapters, provided the original author(s) and source publication are appropriately acknowledged. If so indicated, certain images may not be included under the Creative Commons license. In such cases users will need to obtain permission from the license holder to reproduce the material. More details and guidelines concerning content reuse and adaptation can be foundat http://www.intechopen.com/copyright-policy.html.

#### **Notice**

Statements and opinions expressed in the chapters are these of the individual contributors and not necessarily those of the editors or publisher. No responsibility is accepted for the accuracy of information contained in the published chapters. The publisher assumes no responsibility for any damage or injury to persons or property arising out of the use of any materials, instructions, methods or ideas contained in the book.

First published in Croatia, 2015 by INTECH d.o.o. eBook (PDF) Published by IN TECH d.o.o. Place and year of publication of eBook (PDF): Rijeka, 2019. IntechOpen is the global imprint of IN TECH d.o.o. Printed in Croatia

Legal deposit, Croatia: National and University Library in Zagreb

Additional hard and PDF copies can be obtained from<EMAIL_ADDRESS>
Emerging Pollutants in the Environment - Current and Further Implications Edited by Marcelo L. Larramendy and Sonia Soloneski p. cm. ISBN 978-953-51-2160-2 eBook (PDF) ISBN 978-953-51-5404-4


*Article* **Predicting Primary Energy Consumption Using Hybrid ARIMA and GA-SVR Based on EEMD Decomposition**

#### **Yu-Sheng Kao 1, Kazumitsu Nawata 1 and Chi-Yo Huang 2,\***

Received: 6 September 2020; Accepted: 27 September 2020; Published: 7 October 2020

**Abstract:** Forecasting energy consumption is not easy because of the nonlinear nature of the time series for energy consumptions, which cannot be accurately predicted by traditional forecasting methods. Therefore, a novel hybrid forecasting framework based on the ensemble empirical mode decomposition (EEMD) approach and a combination of individual forecasting models is proposed. The hybrid models include the autoregressive integrated moving average (ARIMA), the support vector regression (SVR), and the genetic algorithm (GA). The integrated framework, the so-called EEMD-ARIMA-GA-SVR, will be used to predict the primary energy consumption of an economy. An empirical study case based on the Taiwanese consumption of energy will be used to verify the feasibility of the proposed forecast framework. According to the empirical study results, the proposed hybrid framework is feasible. Compared with prediction results derived from other forecasting mechanisms, the proposed framework demonstrates better precisions, but such a hybrid system can also be seen as a basis for energy management and policy definition.

**Keywords:** ensemble empirical mode decomposition (EEMD); autoregressive integrated moving average (ARIMA); support vector regression (SVR); genetic algorithm (GA); energy consumption; forecasting

#### **1. Introduction**

Research on energy supply and demand has become critical since the 1973 oil crisis. In the past decades, the average annual worldwide energy consumption grew due to the rapid economic growth of major economies. Based on the forecast by BP [1], the worldwide energy consumption will increase by 34% between 2014 and 2035. Over 90% of the world's energy consumption comes from coal, oil, natural gas, and nuclear sources [1]. Furthermore, energy consumption always plays a dominant role in countries' long-term sustainability. For most industries, e.g., heavy industries, much more energy will be required for sustainable growth. Therefore, the understanding and prediction of energy consumption in general, and of a specific economy in particular, are critical from an economic perspective.

Recently, scholars have started to forecast the demand and supply of energy by integrating various models into a hybrid one [2–4]. In general, such hybrid forecasting methods can be divided into two: causal models and time series models. On the one hand, causal models are mainly constructed based on one or more independent variables. Then the dependent variables can be predicted. Strict assumptions and theoretical bases are required for constructing such causal models. On the other hand, time series models are based on historical data. Whether linear or nonlinear, such models are used for estimating future values. These approaches are always regarded as the most feasible ways to predict energy consumption. The fundamental purpose of time series is to derive trends or patterns that can be modeled by econometric methods such as the autoregressive integrated moving average (ARIMA). However, the models for nonlinear time series predictions are always difficult to realize due to the uncertainty and volatility of the time series. Since the linear models cannot be used to predict complex time series, nonlinear approaches are more suitable for such purpose. Thus, this study aims to predict energy consumption by using integrated methods that incorporate linear and nonlinear methods. Due to issues that arise in time series forecasting, accurate predictions are essential. Conventional linear approaches are effective in the event of forecasting issues. However, more studies are finding that, compared to nonlinear methods, such as support vector machines (SVMs) and tree-based algorithms, linear methods do not perform well in the event of various time series problems, especially complex time sequences. This is because linear methods cannot be used to detect the complex implicit patterns in time series. This study adopts a hybrid model by incorporating linear and nonlinear methods to predict energy consumption and overcome this problem.

However, in accordance with previous studies, no single prediction model is applicable to all scenarios. Therefore, many researchers have introduced hybrid models for predicting time series; such models incorporate both linear and nonlinear models or combine two linear models [5]. Earlier works have also revealed that such hybridization of prediction frameworks not only shows the complementary nature of the frameworks with respect to predictions but also enhances the accuracy of predictions. Thus, models in hybrid forms have become a common practice in forecasting. However, noise and unknown factors exist in time series. These factors influence the volatility of time series and cannot be easily solved by hybridizing linear and nonlinear patterns only. Such hybridization probably produces an overfitting problem; thus, the optimal parameters required to model a prediction framework cannot be derived. Fortunately, such difficulties can be partially solved by leveraging the ensemble empirical mode decomposition (EEMD) proposed by Wu and Huang [6] because the EEMD can solve the noise problems and enhance the prediction performance. Noises, such as trend, seasonality, and unknown factors, which often exist in time series, influence forecasting performance. To make the prediction more precise, noise problems should be carefully dealt with. There are several ways to approach noise problems and enhance forecasting performance. One way is to tune the hyperparameter for algorithms (such as the support vector regression, SVR). As a result, the prediction's performance will improve. Another way is to first deal with the time sequence using a decomposition method, such as the EEMD method. If a decomposition approach is used, the time series can be split into several stable sequences for prediction. The more stable the sequences, the better the model's prediction performance. The feasibility of EEMD has been verified by various works (e.g., [7,8]) in solving nonstationary time series and complex signals.

Given the abovementioned advantages regarding the advantages of the hybrid prediction system as well as the EEMD for time series, this research proposes a novel hybrid framework integrating the EEMD, ARIMA, genetic algorithm (GA), and the SVR to predict primary energy consumption. The concept of the proposed model comes from several sources. First, in conventional time series methods, ARIMA is popular and powerful. It has been extensively used to deal with various forecasting issues. Therefore, the ARIMA method is suitable for this study. Second, based on past research, univariate or single methods used to deal with forecasting problems cannot yield a high forecasting performance when compared to hybrid methods. Therefore, this study simultaneously uses another nonlinear method, SVR, to enhance the prediction performance. In past decades, SVR has been useful in a wide variety of prediction domains. Therefore, SVR is suitable for this study. Although SVR has obtained several important prediction records, it has a significant problem: hyperparameter tuning. Hyperparameter tuning will affect forecasting performance. It is thus necessary to find a way to select the ideal hyperparameter for SVR. Due to computing speed, the authors cannot spend much time searching for the ideal hyperparameters. Consequently, the greedy algorithm (grid search) will be abandoned. Instead, the authors intend to adopt a heuristic algorithm, such as the GA approach, to find the best hyperparameters in SVR. Finally, decomposition methods' effectiveness in enhancing models' forecasting performance has been verified by past literature. The EEMD, Wavelet, and other

methods have been hugely successful in the signal processing field. Due to the innovativeness and power of EEMD, recent studies have widely adopted the method in signal processing research.

Based on above-mentioned reasons, the authors attempt to combine the EEMD, ARIMA, GA, and SVR into a prediction model for energy consumption. First, the time series of energy consumption is divided into several intrinsic mode functions (IMFs) and a residual term. Here, IMFs stand for different frequency bands of time series, which range from high to low, while each IMF represents a series of oscillatory functions [9,10]. That way, each time series component can be identified and modeled accordingly. The characteristics of the time series can also be captured in detail. The ARIMA is then introduced to predict the future values of all extracted IMFs and the residue independently. Since the accuracy of the nonlinear time series derived using the ARIMA may be unacceptable, the SVR is utilized based on the nonlinear pattern to further improve prediction performance. In addition, the accuracy of the SVR-based prediction models completely depends on the control parameters, and the parameters should be optimized. Therefore, the GA is leveraged to derive the optimal parameters. In general, the prediction model fuses EEMD, ARIMA, SVR, and GA into a hybrid prediction framework. The predicted IMFs and the residue will be split into a final ensemble result. The proposed hybrid framework will be verified by a prediction of Taiwanese primary energy consumption for the next four years. Meanwhile, the accuracy of the prediction results will be compared with the ones derived by other forecast methods, which include ARIMA, ARIMA-SVR, ARIMA-GA-SVR, EEMD-ARIMA, EEMD-ARIMA-SVR, and EEMD-ARIMA-GA-SVR. Based on the empirical study results, the effectiveness of the prediction framework can be verified.

The remaining part of this work is structured as follows. Section 2 reviews the related literature regarding the consumption and forecasting methods of energy. Section 3 introduces the methods used in this paper, which include ARIMA, SVR, and GA. Section 4 describes the background of the empirical study case, the dataset, and the empirical study process. Section 5 concludes the whole work, with major findings and opportunities for future studies.

#### **2. Literature Review**

For decision-makers to effectively understand the trend of energy fluctuation, which may generate far-reaching implications, the precise forecasting of primary energy consumption is indispensable. On the one hand, with an accurate forecasting of energy consumption, the government can establish energy plans for fluctuations in oil supply. On the other hand, energy predictions can be useful for the investment of firms. Albeit important, predicting energy consumption is not always simple. Therefore, a robust forecasting model will be necessary.

In the literature on energy prediction, several researchers have completed accurate predictions [2–4,11,12]. Some researchers adopted economic indicators by mixing various energy indicators for predicting energy consumption. Other researchers used only time series data for forecasting. While these forecast methods are different, the prediction results of the two categories of models can serve as solid foundations for further investigations on energy consumption.

Furthermore, to enhance the accuracy of predictions, some authors employed hybrid models. Of such models, linear and nonlinear ones were integrated. The fusion of linear and nonlinear models can overcome the shortage of adopting only one kind of method and provide more accurate results [3,5,13]. In addition to the hybrid methods involving both linear and nonlinear models, some studies have attempted to transform the data by integrating data preprocessing and post-processing procedures [14,15]. By doing so, the forecasting capabilities of the hybrid models with data preprocessing and postprocessing procedures can show superior performance in energy predictions.

Further, several studies have proposed machine learning methods for predicting energy consumption. Al-Garni, Zubair, and Nizami [16] used weather factors as explanatory variables of a regression model for predicting the consumption of electric energy in eastern Saudi Arabia. Azadehet al. [12] modeled Turkey's electricity consumption on the sector basis by utilizing the

artificial neural network (ANN). Wong, Wan, and Lam [17] developed an ANN-based model for analyzing the energy consumption of office buildings. Fumo and Biswas [11] employed simple and multiple regression models as well as the quadratic regression model to predict residential energy consumption. Ahmadet al. [13] attempted to review the applications of ANN and the SVM for energy consumption prediction and found that both seem to show better performance in energy forecasting [13]. Ardakani and Ardehali [3] applied regressive methods consisting of linear, quadratic, and ANN models by incorporating an optimization algorithm into the model to achieve better performance in predicting long-term energy consumption.

Although many scholars have empirically verified the effectiveness of machine learning methods for dealing with time series problems, no single prediction model seems applicable to all scenarios. That is, even if the machine learning model outperforms other traditional linear methods, using a single machine learning model to address all time series issues would be problematic and unrealistic.

Many researchers have thus employed hybrid time series forecasting models, which incorporate linear with nonlinear models or combine two kinds of linear models [5]. Previous studies have also revealed that such hybrid frameworks not only complement each other in prediction but also enhance prediction accuracy. Thus, models in hybrid forms have become a common practice in forecasting.

For example, Yuan and Liu [2] proposed a composite model that combined ARIMA and the grey forecasting model, GM (1,1), to predict the consumption of primary energy in China. Based on their findings, the results obtained when using the hybrid model were far superior to those obtained when only using the ARIMA or GM (1,1) models. Zhang [18] developed a hybrid prediction model consisting of both ARIMA and ANN. Zhu et al. [4] developed a hybrid prediction model of energy demands in China, which employed the moving average approach by integrating the modified particle swarm optimization (PSO) method for enhancing prediction performance. Wang, Wang, Zhao, and Dong [19] combined the PSO with the seasonal ARIMA method to forecast the electricity demand in mainland China and obtain a more accurate prediction. Azadeh, Ghaderi, Tarverdian, and Saberi [20] also adopted the GA and ANN models to predict energy consumption based on the price, value-add, number of customers, and energy consumption. Lee and Tong [21] proposed a model that combined ARIMA and genetic programming (GP) to improve prediction efficiency by adopting the ANN and ARIMA models.

Further, Yolcuet et al. [5] developed a linear and nonlinear ANN model with the modified PSO approach for time series forecasting. They achieved prediction results superior to those of conventional forecasting models. According to the analytical results, the hybrid model is more effective because it adopts a single prediction method and can thus improve the prediction accuracy.

Though hybrid models based on ARIMA and ANN have achieved great success in various fields, they have several limitations. First, this hybrid approach requires sufficient data to build a robust model. Second, the parameter control, uncertainties in weight derivations, and the possibility of overfitting must often be discussed when using ANN models. Because of these limitations, more researchers started to adopt SVR in forecasting since it can mitigate the disadvantages of ANN models. SVR is suitable for forecasts based on small datasets.

Pai and Lin's [22] work is a representative example of adopting SVR methods in hybrid models for forecasting. They integrated ARIMA and SVR models for stock price predictions. Patel, Chaudhary, and Garg [23] also adopted the ARIMA-SVR for predictions and derived optimal results based on historical data. Alwee, Hj Shamsuddin, and Sallehuddin [24] optimized an ARIMA-SVR-based model, using the PSO for crime rate predictions. Fan, Pan, Li, and Li [25] employed independent component analysis (ICA) to examine crude oil prices and then used an ARIMA-SVR-based model to predict them.

Based on the results of the literature review, hybrid models, including both the SVR and the ANN, have achieved higher prediction accuracies than traditional prediction techniques. However, the invisible and unknown factors which can influence the volatility of time series cannot be addressed easily by hybridizing linear and nonlinear patterns. The problem of overfitting can emerge; thus, the optimal parameters cannot be derived.

Fortunately, such difficulties can be partially resolved by leveraging the EEMD proposed by Wu and Huang [6]. The method has been feasible and effective in solving problems consisting of nonstationary time series and complex signals [7,8]. Wang et al. [7] integrated the EEMD method with the least square SVR (LSSVR) and successfully predicted the time series of nuclear energy consumption. Prediction performance has increased significantly and outperformed some well-recognized approaches based on level forecasting and directional prediction. Zhang, Zhang, and Zhang [26] predicted the prices of crude oil by hybridizing PSO-LSSVM and EEMD decomposition. The work demonstrated that the EEMD technique can decompose the nonstationary and time-varying components of times series of crude oil prices. The hybrid model can be beneficial to model the different components of crude oil prices and enhance prediction performance.

The previous studies in the literature review section aimed to develop a model that could effectively and accurately predict energy consumption and demand. In their methodologies, these works being reviewed attempted to use linear or nonlinear methods to predict energy consumption. Furthermore, they tried to use the parameter search algorithm in their model to enhance its prediction accuracy. Based on the review results, complex time series can be split by the EEMD into several relatively simple subsystems. The hidden information behind such complex time series can be explored more easily. Thus, in the following section, a hybrid analytical framework consisting of ARIMA, SVR, and GA will be proposed. The framework will be adopted to predict primary energy consumption.

#### **3. Research Methods**

This section first introduces the data processing method. Next, the individual models including ARIMA and SVR will be introduced. Afterward, the optimization approach based on GA will be introduced. Finally, the analytical process of the proposed hybrid model will be described.

#### *3.1. EEMD*

Empirical mode decomposition (EMD), an adaptive approach based on the Hilbert–Huang transformation (HHT), is often used to deal with time series data including ones with nonlinear and nonstationary forms [8]. Since such time series are complicated, various fluctuation modes may coexist. The EMD technique can be used to decompose the original time series into several simple IMFs, which correspond to different frequency bands of the time series and range from high to low; each IMF stands for a series of oscillatory functions [9,10]. Moreover, the IMFs must satisfy two conditions [6]: (1) in the whole data series, the number of extrema and zero crossings must either be equal or differ at most by 1; and (2) at any point, the mean value of the envelope (envelope, in mathematics, is a curve that is tangential to each one of a family of curves in a plane) defined by the local minima is 0.

Based on the above definitions, IMFs can be extracted from the time series *y*(*t*) according to the following shifting procedures [27]: (1) identify the local maxima and the minima; (2) connect all local extrema points to generate an upper envelope *emax*(*t*) and connect all minima points to generate a lower envelope *emin*(*t*) with the spline interpolation, respectively; (3) compute the mean of the envelope, *a*(*t*), from the upper and lower envelopes, where *a*(*t*)=(*emax*(*t*) + *emin*(*t*))/2; (4) extract the mean from the time series and define the difference between *y*(*t*) and *a*(*t*) as *c*(*t*), where *c*(*t*) = *y*(*t*) − *a*(*t*); (5) check the properties of *c*(*t*): (i) if *c*(*t*) satisfies the two conditions illustrated above, an IMF will be extracted and replace *y*(*t*) with the residual, *r*(*t*) = *y*(*t*) − *c*(*t*); (ii) if *c*(*t*) is not an IMF, then *y*(*t*) will be replaced by *c*(*t*); and (6) the residue *r*1(*t*) = *y*(*t*) − *c*1(*t*) is regarded as the new data subjected to the same shifting process, which was described above for the next IMF from *r*1(*t*). When the residue *r*(*t*) becomes a monotonic function or at most has one local extrema point from which no more IMF can be extracted [27], the shifting processes can be terminated.

Through the abovementioned shifting process, the original data series *y*(*t*) can be expressed as a sum of IMFs and a residue, *y*(*t*) = *m i*=1 *ci*(*t*) + *rm*(*t*), where *m* is the number of IMFs, *rm*(*t*) is the final residue, and *ci*(*t*) is the *i*th IMF. All the IMFs are nearly orthogonal to each other, and all have nearly zero means.

Although the EMD has been widely adopted in handling data series, the mode-mixing problem still exists. The problem can be defined as either a single IMF consisting of components of widely disparate scales or a component of a similar scale residing in different IMFs. To overcome this problem, Wu and Huang [6] proposed the ensemble EMD (EEMD), which adds white noise to the original data, and thus the data series of different scales can be automatically assigned to proper scales of reference built by the white noise [7]. The core concept of the EEMD method is to add the white noise into the data processing. White noise can be viewed as a sequence with zero mean value; this sequence does not fall under any distribution. Based on different algorithms, the white noise can be assigned to specific distributions for calculation. In EEMD, the purpose of this method is to make the original sequence the stable sequence. Hence, this method employs simulation, using the original sequence to generate various sequences of normal distributions—this is the white noise concept. With this method, the original sequence can be split into several different sequences. Meanwhile, the sum of the decomposed sequences equals the original sequence. These decomposed sequences are called IMFs. This way, the mode-mixing problem can be easily solved. The EEMD procedure is developed as follows [6]: (1) add a white noise series to the original data; (2) decompose the data with added white noise into IMFs; (3) repeat steps 1 and 2 iteratively, but with different white noise each time; and (4) obtain ensemble means of corresponding IMFs as the final results.

In addition, Wu and Huang [6] established a statistical rule to control the effect of added white noise: *en* = ε/ √ *n*, where *n* is the number of ensemble members, ε represents the amplitude of the added noise, and *en* is the final standard deviation of error, which is defined as the difference between the input signal and the corresponding IMFs. Based on previous studies, the number of ensemble members is often set to 100, and the standard deviation of white noise is set to 0.2.

#### *3.2. The ARIMA Model*

The ARIMA model for forecasting time series was proposed by Box, Jenkins, and Reinsel [28]. The model consists of the autoregressive (AR) and the moving average (MA) models. The AR and MA models were merged into the ARMA model, which has already become matured in predictions. The future value of a variable is a linear function of past observations and random errors. Thus, the ARMA can be defined as

$$y\_t = \phi\_1 y\_{t-1} + \phi\_2 y\_{t-2} + \dots + \phi\_p y\_{t-p} + \varepsilon\_t - \theta\_1 \varepsilon\_{t-1} - \theta\_2 \varepsilon\_{t-2} - \dots - \theta\_q \varepsilon\_{t-q} \tag{1}$$

where *yt* is the forecasting value; φ*i* is the coefficient of the *i*th observation; *yt*−*i* is the *i*th observation; θ*i* is the parameter associated with the *i*th white noise; ε*t* is the white noise, whose mean is zero; and ε*t*−*i* is the noise terms.

The ARMA model can satisfactorily fit the original data when the time series data is stationary. However, if the time series are nonstationary, the series will be transformed into a stationary time series using the *d*th difference process, where *d* is usually set as 0, 1, or 2. ARIMA is used to model the differenced series. The process is called ARIMA (*p*, *d*, *q*), which can be expressed as

$$w\_t = \phi\_1 w\_{t-1} + \phi\_2 w\_{t-2} + \dots + \phi\_p w\_{t-p} + \varepsilon\_t - \theta\_1 \varepsilon\_{t-1} - \theta\_2 \varepsilon\_{t-2} - \dots - \theta\_q \varepsilon\_{t-q} \tag{2}$$

where *wt* is denoted as ∇*d yt*. When *d* equals zero, the model is the same as Equation (1) of ARMA.

The ARIMA model was developed by using the Box–Jenkins method. The procedure of the ARIMA involves three steps: (1) Model identification: Since the stationary series is indispensable for the ARIMA model, the data needs to be transformed from a nonstationary one to a stable one. For the stability of the series, the difference method is essential for removing the trends of the series. This way, the *d* parameter is determined. Based on the autocorrelation function (ACF) and the partial autocorrelation function (PACF), a feasible model can be established. Through the parameter estimation

and diagnostic checking process, the proper model will be established from all the feasible models. (2) Parameter estimation: Once the feasible models have been identified, the parameters of the ARIMA model can be estimated. The suitable ARIMA model based on Akaike's information criterion (AIC) and Schwarz's Bayesian information criterion (BIC) can be further determined. (3) Diagnostic checking: After parameter estimation, the selected model should be tested for statistical significance. Meanwhile, hypothesis testing is conducted to examine whether the residual sequence of the model is a white noise. Based on the above procedures, the forecasting model will be determined. The derived model will be appropriate as the training model for predictions.

In this research, the ARIMA models will be built by feeding each decomposed data. The separated ARIMA models will then be integrated with an SVR model for further analysis. Fitting performance is expected to be enhanced further.

#### *3.3. SVR*

SVR was proposed by Vapnik [29], who thought that theoretically, a linear function *f* exists to define the nonlinear relationship between the input and output data in the high-dimensional-feature space. Such a method can be used to solve the function with respect to fitting problems. Based on the concepts of SVR, the basis function can be described as

$$f(\mathbf{x}) = \mathbf{w}^T \cdot \boldsymbol{\varphi}(\mathbf{x}) + b,\tag{3}$$

where *f*(*x*) denotes the forecasting values, *x* is the input vector, *w* is the weight vector, *b* is the bias, and ϕ(*x*) stands for a mapping function to transform the nonlinear inputs into a linear pattern in a high-dimensional-feature space.

Conventional regression methods take advantage of the square error minimization method for modeling the forecasting patterns. Such a process can be regarded as an empirical risk in accordance with the loss function [29]. Therefore, the ε − insensitive loss function (*T*ε) is adopted in the SVR and can be defined as

$$T\_{\varepsilon}(f(\mathbf{x}), y) = \left\{ \begin{array}{c} \left| f(\mathbf{x}) - y \right| - \varepsilon \\ 0 \end{array} \right\} \qquad \text{if} \left| f(\mathbf{x}) - y \right| \ge \varepsilon \\ \text{otherwise} \end{array} \tag{4}$$

where *y* is the target output and ε is expressed as the region of ε − insensitive. When the predicted value falls into the band area, the loss is equal to the difference between the predicted value and the margin [30]. *T*ε(*f*(*x*), *y*) is leveraged to derive an optimum hyperplane on the high-feature space to maximize the distance which can divide the input data into two subsets. The weight vector (w) and constant (*b*) in Equations (4) can be estimated by minimizing the following regularized risk function:

$$P(c) = c \frac{1}{n} \sum\_{i=1}^{n} T\_{\varepsilon} \left( f(x\_i), y\_i \right) + \frac{1}{2} \|\mathbf{w}\|^2,\tag{5}$$

where *T*ε(*f*(*xi*), *yi*)is the ε−insensitive loss function in Equation (5). Here,1/2 w 2 plays the regularizer role, which tackles the problem of trade-off between the complexity and approximation accuracy of the regression model.

Equation (5) above aims to ensure that the forecasting model has an improved generalized performance. In the regularization process, *c* is used to specify the trade-off between the empirical risk and the regularization terms. Both *c* and ε can be defined by hyper-parameter search algorithms and users. These parameters significantly determine the prediction performance of the SVR.

In addition, based on the concept of the tube regression, if the predicted value is within the ε − tube, the error will be zero. However, if the predicted value is located outside the ε − tube, the error will be produced. Such an error, the so-called ε − insensitive error, is calculated in terms of the distance between the predicted value and the boundary of the tube. Since some predicted values exist outside the tube, the slack variables (ξ + ξ∗ *i* ) are introduced and defined as tuning parameters. These variables stand for the distance from actual values to the corresponding boundary values of the tube. Given the synchronous structural risk, Equation (5) is transformed into the following constrained form by using the slack variables:

$$\begin{aligned} \text{Minimize } & \frac{1}{2} \|w\|^2 + c \sum\_{i=1}^n \left(\xi + \xi\_i^\*\right) \\ \text{Subject to} & \frac{1}{2} \end{aligned}$$

$$\left\{ \begin{array}{ll} y\_i - \left( \mathbf{w}^T \cdot \boldsymbol{\varphi}(\mathbf{x}\_i) \right) - b \le \varepsilon + \xi\_i \\ \left( \mathbf{w}^T \cdot \boldsymbol{\varphi}(\mathbf{x}\_i) \right) + b - y\_i \le \varepsilon + \xi\_i^\* \\ \xi\_{i\prime} \xi\_j^\* \ge 0, \qquad \text{for} \quad i = 1, \dots, n \end{array} \right\} \tag{6}$$

The constant *c* determines the trade-off between the flatness of *f* and the amount up to which deviations larger than ε are tolerated. To solve the above problem, the Lagrange multiplier and the Karush–Kuhn–Tucker (KKT) conditions will be leveraged. After the derivation, the general form of the SVR function can be expressed as

$$f(\mathbf{x}, \mathbf{w}) = f(\mathbf{x}, a\_i, a\_i^\*) = \sum\_{i=1}^n \left(a\_i - a\_i^\*\right) \times K(\mathbf{x}, \mathbf{x}\_i) + b \tag{7}$$

where α*i* and α∗ *i* are the Lagrange multipliers, and *K*(x, x*i*) is the inner product of two vectors in the feature space, ϕ(x*i*) and ϕ(x*j*). Here, *K*(x, x*i*) is called the kernel function. In general, the most popular kernel function is the Gaussian radial basis function (RBF), which can be defined as *K*(x*i*, x*j*) = exp(− x*i* − x*j* 2/2σ2) [29]. In this work, the RBF is employed for predictions. In addition, since σ is a free parameter, the RBF kernel can be described as *K*(x*i*, x*j*) = exp(−γ x*i* − x*j* 2), where γ is a parameter of the RBF kernel.

#### *3.4. Optimization by GA*

While defining a prediction model, the enhancement of prediction accuracy and the avoidance of overfitting are the most important tasks. By doing so, the training model can achieve far better performance in predictions when testing data are inputted. In the SVR-based models, the *c*, ε, and γ parameters play a dominant role in determining modeling performance. That is, if these parameters can be defined correctly and appropriately, the forecasting model will be efficient. To select the best parameters, the method for searching these parameters will be indispensable. In this work, GA will be adopted in selecting the optimal values for these three SVR parameters.

GA, a concept first proposed by Holland [31], is a stochastic search method based on the ideas of natural genetics and the principle of evolution [32,33]. GA works with a population of individual strings (chromosomes), each of which stands for a possible solution to a given problem. Each chromosome is assigned a fitness value based on the result of the fitness function [34]. GA allows more opportunities to fit chromosomes to reproduce the shared features originating from their parent generation. It has been regarded as a useful tool in many applications and has been extensively applied to derive global solutions to optimization problems. The algorithm is also applicable to large-scale and complicated nonlinear optimal problems [35]. The GA procedure is summarized below based on the work by [36]:

Step 1: Randomly generate an initial population of *n* agents; each agent is an *n* − bit genotype (chromosome).

Step 2: Evaluate the fitness of each agent.

Step 3: Repeat the following procedures until *n* offspring has been created.

Step 4: Replace the current population with the new one.

Step 5: The whole process is finished when the stopping condition is satisfied. Then the best solution is returned to the current population. Otherwise, the process will go back to Step 2 until the terminating condition can be satisfied.

#### *3.5. The Proposed Hybrid Model*

The abovementioned prediction approaches including ARIMA and SVR deliver good performance in general since these methods deal with regression problems effectively. For linear models, ARIMA performs quite well in forecasting time series. However, its prediction performance is limited since it cannot appropriately predict highly nonlinear and nonstationary time series. Therefore, the data stabilization technique based on EEMD decomposition is introduced to handle nonlinear data; the nonstationary process is introduced as well. Complexities such as randomness and intermittence often exist in the time series. Even if the series have been proceeded by the EEMD, unknown factors that influence the series remain. To enhance the forecasting accuracy of the EEMD-ARIMA, the SVR method is integrated. SVR based on the nonlinear method is good at coping with small data and unstable data series. Thus, the EEMD-ARIMA, integrated with the SVR, will be useful in predicting nonlinear time series generally and energy consumption specifically. Meanwhile, GA is adopted to derive the best parameters for SVR, which can improve the performance of the hybrid model.

In general, the proposed EEMD-ARIMA-GA-SVR prediction framework (Figure 1) is composed of the following steps:

**Figure 1.** General procedure of the proposed ensemble empirical mode decomposition- autoregressive integrated moving average-genetic algorithm-support vector regression (EEMD)-(ARIMA)-(GA)-(SVR) modeling framework for primary energy consumption forecasting.

Step 1: The original time series of primary energy consumptions, *yt*, *t* = 1, 2, ... , *n*, is decomposed into *m* IMF components *ci*(*t*), *i* = 1, 2, ... , *m*, and a residual component *rm*(*t*) using the EEMD method.

Step 2: ARIMA and SVR are introduced as stand-alone prediction methods to extract the IMFs and residual of the time series, respectively. Meanwhile, GA is introduced to optimize the parameters associated with SVR. Accordingly, the corresponding prediction results for all components can be obtained.

Step 3: The independent prediction results of all the IMFs and the residual are aggregated as an output, which can be regarded as the final prediction results of the original time series *y*(*t*).

Thus, the fitted values can be accordingly derived from this proposed hybrid prediction framework. Further, to demonstrate the effectiveness of the proposed hybrid EEMD-ARIMA-GA-SVR framework, the time series of Taiwanese primary energy consumption will be adopted to verify the feasibility of the proposed framework. Meanwhile, the prediction results based on ARIMA, ARIMA-SVR, ARIMA-GA-SVR, EEMD-ARIMA, and EEMD-ARIMA-SVR will be introduced for comparison.

#### *3.6. Performance Measures for Predictions*

Different measures for prediction errors will be adopted to evaluate the accuracy of the prediction models. In this research, the mean absolute error (MAE), mean absolute percentage error (MAPE), mean square error (MSE), and root mean square error (RMSE) will be adopted. The four metrics are used here to evaluate the forecasting performance. The MAPE and RMSE metrics can be useful to explain the performance of predictions. To yield more accurate evaluation results, we further provide results being derived by MAE and MSE as references. These performance measures are defined below:

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} \left| y\_t - \hat{y}\_t \right| \tag{8}$$

$$\text{MAPE} = \frac{1}{n} \sum\_{t=1}^{n} \left| \frac{y\_t - \hat{y}\_t}{y\_t} \right| \times 100 \tag{9}$$

$$\text{MSE} = \frac{1}{n} \sum\_{i=1}^{n} \left( y\_t - \hat{y}\_t \right)^2 \tag{10}$$

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left( y\_t - \hat{y}\_t \right)^2} \tag{11}$$

where *n* stands for the size of the test data and *yt* as well as *y*ˆ*t* denote the actual value and the predicted value. Based on these measures, the lower values of all performance measures represent superior forecasting. The MAE reveals how similar the predicted values are to the observed values while the MSE and RMSE measure the overall deviations between the fitted and predicted values. Generally, the MAPE value should be less than 10%. To prove the effectiveness of the proposed hybrid EEMD-ARIMA-GA-SVR model in forecasting, alternative methodologies consisting of ARIMA, ARIMA-SVR, ARIMA-GA-SVR, EEMD-ARIMA, and EEMD-ARIMA-SVR will be used as benchmark models to compare with proposed model.

To test the differences of means between fitted and actual values, the paired-sample Wilcoxon signed-rank test is introduced in this research, where μ1 stands for the mean of the actual data while μ2 represents the mean of the fitted data. The null hypothesis is defined as *H*0: μ1 − μ2 = 0 while the alternative hypothesis is defined as *H*1: μ1 − μ2 - 0. The null hypothesis cannot be rejected while the result of the subtraction between the mean value of the actual data and the mean value of the fitted data is 0. Such a case means that there is no mean difference between the fitted data and the actual data. In other words, the training model is suitable as a prediction model since it is consistent with the real situation.

#### **4. Forecasting Taiwan's Primary Energy Consumption**

In this section, an empirical case based on Taiwanese primary energy consumption is presented to verify the feasibility and effectiveness of the proposed hybrid framework. Comparisons with other benchmark models will be provided to demonstrate the forecasting capabilities of the proposed framework. The background of the Taiwanese primary energy consumption will be presented first. Then the raw data and analytic process will be introduced in Section 4.2. Afterward, the time series will be decomposed by the EEMD for the predictions in Section 4.3. Modeling via ARIMA will be introduced in Section 4.4. The predictions of energy consumptions using the EEMD-ARIMA-SVR model optimized by the GA method will be discussed in Section 4.5. Finally, the evaluations of hybrid models and forecasting results will be described in Section 4.6.

#### *4.1. Background*

Because of its shortage of natural resources, Taiwan relies heavily on energy imports. In managing energy imports, energy consumption predictions are indispensable. Such energy predictions can help the government sector define relevant energy policies for sustainable development.

#### *4.2. Raw Data and Modeling Method*

In this study, the time series of the Taiwanese primary energy consumption (Figure 2) was adopted to verify the effectiveness of the prediction models. Annual primary energy consumption from 1965 to 2014 was provided by BP [1]. The time series was separated into two subsets where 90% (46 samples) of the dataset were chosen as the training set while the remaining 10% (4 samples) were selected as the test set for verifying the prediction efficiency.

**Figure 2.** Taiwan's primary energy consumption (million barrels of oil equivalent). Source: BP Statistical Review of World Energy 2015.

For accurate predictions, the original dataset will be transformed by the EEMD method. After the decomposition of the time series, the results of decompositions can be further predicted by the ARIMA method. Then, the prediction results derived from the ARIMA method will be aggregated as the final prediction results for energy consumption.

Meanwhile, such procedures can be extended by adopting SVR. Moreover, to help build the training model and avoid overfitting risks, the *k*-fold cross-validation was introduced into model construction within the prediction process. The *k*-fold cross-validation is used in SVR prediction error examinations based on selected hyperparameters from GA. In this study, the authors adopt five-fold cross-validation. A five-fold validation entails that all samples are divided into five portions; four of the five portions are used for training, and one of the five portions is for testing. The process is repeated five times. At first, the GA algorithm generates a set of hyperparameters. Through the five-fold cross-validation, we can obtain a performance. Next, the GA algorithm will continue to generate different hyperparameters based on cross-validation, thus obtaining a performance. Finally, we can select the best hyperparameters in terms of the best performances. The training data was randomly divided into *k* subsamples. Among the *k* samples, the *k* − 1 subsamples were selected as the training data, and the remaining subsample was considered as validation data for testing the model. This research adopted 5-fold cross-validation into the model construction process.

#### *4.3. Data Preprocessing Using EEMD Decomposition*

Before conducting the prediction using the proposed hybrid framework, the time series of energy consumptions will be processed using the EEMD decomposition method, which separates the original time series into several IMFs and a residue. The results are depicted in Figure 3. The four independent IMF components correspond to different frequency bands of the time series, which range from high to low.

**Figure 3.** Decomposition results for primary energy consumption via EEMD.

In Figure 3, the IMFs stand for changing frequencies, amplitudes, and wavelengths. IMF1 represents the highest frequency, maximum amplitude, and the shortest wavelength. The frequencies and amplitudes associated with the rest of the IMFs are lower while the wavelengths are longer. The residue represents a mode slowly varying around the long-term average.

#### *4.4. Forecasting with the ARIMA Model*

To establish the prediction model of primary energy consumption in terms of the historical time series dataset, the ARIMA method is introduced. Since the ARIMA model must be built into the series stationarity, the *d* times difference needs to be obtained to have an ARIMA (*p*, *d*, *q*) model with *d* as the order of differencing used. To test the stationarity of *d* times differencing, the augmented Dickey–Fuller (ADF) test is utilized.

The construction of the ARIMA model depends on model identification. Here, differencing will be important for solving non-stationarity. Moreover, the order of AR (*p*) and MA (*q*) needs to be identified. Through ACF and PACF, the order of AR (*p*) and MA (*q*) can further be determined [28]. However, the ACF and PACF method may not be useful when performing hybrid ARMA processes. Commonly, AIC or BIC measures can be used to easily inspect the appropriateness of the ARMA model. In this research, the best forecasting model is determined based on AIC measures and statistically significant results. Once the forecasting model has been selected, the nonlinear method based on the optimization approach will be integrated into this model.

First, the stationary test for the decomposed data series and the first difference of the original time series derived from the ADF test is implemented and shown in Table 1. According to the results of the ADF test, all the data series belonged to the stationary. That is, all the transformed data can be used to model constructions.

**Table 1.** Stationary analysis of primary energy consumption in Taiwan.

To further determine the parameters of the ARIMA order, ACF and PACF will be utilized. Likewise, such a procedure can also be followed for the decomposed data set using EEMD. To simplify the analytical procedure of EEMD-ARIMA, the corresponding ACF and PACF diagrams are not presented here.

Through the AIC and the statistical significance test, the suitable ARIMA models can be derived for the first difference time series and decomposed time series via the EEMD. Then the optimal form is specified as ARIMA(1,1,1). Table 1 shows the test results in the first difference of the original sequence. The original sequence is not stationary. Based on the first difference in the original difference, the sequence shows a stable status. Therefore, the difference *d* will be set as 1. IMF1 ~ IMF4 and the residual are derived from the original sequence using the EEMD method. The sum of the IMFs and the residual is equal to the original sequence. The rationality has already been explained in the fourth paragraph of Section 3.1, where the original sequence was split into several different sequences. The sum of decomposed sequences equals to original sequence. Finally, after determining the proper parameters of the ARIMA models, whether the residual of the selected model possesses the autocorrelation problem should be confirmed. Therefore, the ACF and PACF tests were conducted to verify the selected model. Figure 4 demonstrates the estimated residuals using the ACF and PACF tests. According to the test results, no autocorrelation and partial autocorrelation exist within the residuals.

**Figure 4.** Residual error by autocorrelation function (ACF) and partial autocorrelation function (PACF) for ARIMA(1,1,1).

#### *4.5. Forecasting with the EEMD-ARIMA-SVR Model Optimized by the GA Method*

After building the EEMD-ARIMA model, the SVR method will be introduced to reduce the errors produced by ARIMA and enhance forecasting accuracy. According to earlier works, the prediction performance of the SVR method is outstanding; further, it can be fused with other nonlinear or linear methods successfully. Thus, after building the EEMD-ARIMA model, the SVR method will be introduced to reduce the errors produced by ARIMA and enhance forecasting accuracy.

Regarding the hybrid models in this research, the ARIMA initially served as a preprocessor to filter the linear pattern of the decomposed data series. Then, the error terms of the ARIMA model were fed into the SVR model to improve prediction accuracy. Generally, three parameters, *c*, ε, and γ, can influence the accuracy of the SVR model. Currently, no clear definition and standard procedure are available for determining the above three parameters [21]. However, the improper selection of the three parameters will cause either overfitting or underfitting. To prevent these, the GA method with cross-validation will be introduced to derive the best parameters for constructing the forecasting model. Meanwhile, some studies have pointed out that the utilization of the RBF can yield better prediction performance [21,22]. Thus, the RBF kernel with the 5-fold cross-validation based on the RMSE measure is adopted to help derive the best parameters of the SVR using GA. The above procedures will be applied to the ARIMA-SVR and EEMD-ARIMA-SVR models.

In GA, the number of iterations is set as 100, the population size is defined as 50, and the maximum number of iterations is defined as 50. The search boundaries for *c*, ε, and γ are within the intervals [ 10−4, 102 ], [ 10−4, 2 ], and [ 2−4, 22 ], respectively. The optimal values for *c*, ε, and γ can thus be 0.441, 3.813, and 0.219, respectively. Further, the *c*, ε, and γ parameters of the decomposed time series belonging to the four independent IMFs and the residual (Figure 3) are summarized in Table 2. Once the parameters have been derived using GA, the optimal hybrid prediction model can be established.

**Table 2.** The parameter selection of forecasting model optimized by GA.

#### *4.6. Evaluations of Hybrid Models and Forecasting Results*

After the construction of the hybrid models, the effectiveness of predictions will be compared further with those of different models including ARIMA, ARIMA-SVR, ARIMA-GA-SVR, EEMD-ARIMA, and EEMD-ARIMA-SVR. The superiority of the proposed EEMD-ARIMA-GA-SVR model of forecasting capability will be verified accordingly. The parameters derived using the GA with the fivefold cross-validation will be utilized to construct the SVR model. The parameters can yield better forecasting performances in related ARIMA-SVR models.

Based on the proposed hybrid framework and the five models for comparisons, the primary energy consumption in Taiwan for 2010–2014 is predicted and summarized in Table 3. According to the prediction results, EEMD-ARIMA-SVR and EEMD-ARIMA-GA-SVR outperformed the other four models. Meanwhile, the prediction results derived using ARIMA and ARIMA-SVR were unsatisfactory from the aspect of inconsistencies between predicted versus actual values. The four prediction performance measures, MAE, MAPE, MSE, and RMSE, derived from the six training models versus the actual values are illustrated in Figure 5 and summarized in Table 4. Based on the results of comparisons, the proposed model outperformed ARIMA. MAE and MSE decreased by 70.43% and 93.28% in the training stage, respectively; further, the two measures decreased by 71.89% and 88.51% in the testing stage, respectively. From the aspects of MAPE and RMSE, both measures improved by 64.68% and 74.07% in the training stage, respectively; they also improved by 71.85% and 66.10% in the testing stage, respectively.

**Table 3.** Data Test by four sample ranging from 2011 to 2014.

**Figure 5.** Fittings of the six models.

The proposed model outperformed ARIMA-SVR. MAE and MSE decreased by 67.85% and 91.59% in the training stage, respectively; further, the two measures decreased by 61.52% and 80.31% in the

testing stage, respectively. Meanwhile, MAPE and RMSE improved by 62.26% and 71.01% in the training stage and by 61.44% and 55.63% in the testing stage, respectively.

**Table 4.** Comparison of forecasting indices.

Remark: \*: numbers in percentage.

Compared with ARIMA-GA-SVR, the proposed model performed better as well. MAE and MSE decreased by 67.14% and 91.02%, in the training stage and by 55.08% and 74.79% in the testing stage, respectively. Further, MAPE and RMSE improved by 61.61% and 70.04% in the training stage and by 54.96% and 49.79% in the testing stage, respectively. Based on the above comparison results, the hybrid model integrated with the EEMD method showed better forecasting performance than other models without the data decomposition process.

The proposed model also outperformed EEMD-ARIMA. MAE and MSE decreased by 17.17% and 37.57% in the training stage and by 34.19% and 45.87% in the testing stage, respectively. Meanwhile, MAPE and RMSE were enhanced by 8.41% and 20.98% in the training stage and by 34.27% and 26.43% in the testing stage, respectively.

Compared with EEMD-ARIMA-SVR, the proposed framework performed better. From the aspect of MAE and MSE, the proposed framework outperformed EEMD-ARIMA-SVR by reducing both measures by 2.76% and 1.80% in the training stage and by 3.57% and 10.33% in the testing stage, respectively. At the same time, MAPE and RMSE were enhanced by 3.19% and 0.91% in the training stage and by 3.58% and 5.30% in the testing stage, respectively.

More specifically, from the aspect of training models, the hybrid models including EEMD-ARIMA-SVR and EEMD-ARIMA-GA-SVR performed better, with relatively smaller residual errors in comparison with results derived from any stand-alone or hybrid models. The predictions based on the models without the EEMD decomposition show the limited forecasting capability of the training models. Similarly, from the perspective of the testing model, the proposed model achieved better performance in predicting primary energy consumption. In this study, the comparison results show that the hybrid models with EEMD decomposition can significantly reduce overall forecasting errors. That is, the EEMD method is useful in manipulating the nonstationary time series; thus, better prediction results can be derived by integrating other forecasting methods. Further, GA can reduce forecasting errors by ARIMA-SVR. Based on the analytical results, the proposed EEMD-ARIMA-GA-SVR is a powerful tool and model for energy consumption prediction.

Finally, to identify the significant differences between the prediction results of any two models adopted in this work, the Wilcoxon signed-rank test is performed. This test was adopted extensively in examining the prediction results of two different models and justifying whether these results are significantly different based on small samples [37,38]. The null and alternative hypotheses are described as follows:

*H*0: μ1 − μ2 = 0 (null hypothesis) and

*H*1: μ1 − μ2 -0 (alternative hypothesis),

where μ1 represents the mean of the actual data and μ2 stands for the mean of the predicted data. The Wilcoxon test can easily determine whether the mean differences are significant or not. All the

*p*-values derived from the testing of the six pairs of methods were higher than 0.05. That is, no mean differences are observed between the test and predicted values derived from each method.

#### **5. Conclusions**

This study presented a hybrid prediction model consisting of ARIMA, SVR, EEMD, and GA. The empirical results have verified the feasibility of the proposed method. Such a hybrid model that combines linear and nonlinear patterns based on the ARIMA and the SVR models as well as adopts data preprocessing and parameter optimization for time series predictions can produce more precise prediction results. From the aspect of limitations and future research possibilities, different datasets on energy consumption can be used at the same time to evaluate whether the forecasting performance of the proposed model will be the best among all prediction models. The data post-processing procedure can be integrated; the differences of the prediction results derived from the proposed model and the framework consisting of the post-processing procedure can be compared. In the future, a dynamic procedure in terms of multiple-step-ahead forecasting can be adopted to replace the one-step-ahead forecasting techniques being used in this work. Much more meaningful and valuable information can be derived for decision-makers in energy predictions.

**Author Contributions:** Y.-S.K. designed, performed research, analyzed the data, and wrote the paper. K.N. advised on the research methods. C.-Y.H. advised on the research methods, re-wrote, and proof-read the whole article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This article was subsidized by the Taiwan Normal University (NTNU), Taiwan and the Ministry of Science and Technology, Taiwan under Grant numbers T10807000105 and MOST 106-2221-E-003-019-MY3.

**Conflicts of Interest:** The authors declare no conflict of interests.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).


**5. Genetic basis**

Short QT syndrome occurs mainly in an autosomal dominant pattern of inheritance with high phenotypic and genetic heterogeneity. To date, potential deleterious rare variants located in nine genes (*CACNA1C*, *CACNA2D1*, *CACNB2*, *KCNH2*,

*AE3: Anion exchanger; BrS: Brugada syndrome; CDSP: Systemic primary carnitine deficiency; GOF: gain-of-function; ICa: Voltage-gated calcium currents; IKr: Rapidly activating potassium currents; IKs: Slowly activating potassium currents; IK1: Inward rectifier potassium currents; LOF: loss-of-function; LQTS: Long QT syndrome; SQTS: Short QT syndrome.*

#### **Table 1.**

*Genes associated with Short QT Syndrome or Shorter than normal QT interval and its phenotypic overlap with the main arrhythmogenic syndromes.*

#### **Figure 3.**

*Diagram of the overlap between the genes associated with the short QT syndrome (SQTS) and the main channelopathies: Brs: Brugada syndrome; LQTS: long QT syndrome and CPVT: catecholaminergic polymorphic ventricular tachycardia.*

*KCNJ2*, *KCNQ1*, *SLC22A5*, *SLC4A3* and *SCN5A*) have been associated with SQTS (**Table 1**) [22]. However, only three genes (*KCNH2*, *KCNQ1* and *KCNJ2*) have been shown to cause SQTS definitively so far, and the *SLC4A3* gene presents moderate evidence [23]. The association of the other genes with SQTS remains controversial (**Figure 3**) [24].


**Acknowledgements**

The authors gratefully acknowledge the support provided by national research council (NRC), Sri Lanka from the research grant (no 11-067). The authors would also like to extend their sincere gratitude to Mr. Nilhan Niles for his cooperation.


**Validation of an Inertial Sensor Algorithm to Quantify Head and Trunk Movement in Healthy Young Adults and Individuals with Mild Traumatic Brain Injury**

**Lucy Parrington 1,2,\*, Deborah A. Jehu 1,2, Peter C. Fino 1,2,3, Sean Pearson 4, Mahmoud El-Gohary 4 and Laurie A. King 1,2,5**

Received: 20 October 2018; Accepted: 16 December 2018; Published: 19 December 2018

**Abstract:** Wearable inertial measurement units (IMUs) may provide useful, objective information to clinicians interested in quantifying head movements as patients' progress through vestibular rehabilitation. The purpose of this study was to validate an IMU-based algorithm against criterion data (motion capture) to estimate average head and trunk range of motion (ROM) and average peak velocity. Ten participants completed two trials of standing and walking tasks while moving the head with and without moving the trunk. Validity was assessed using a combination of Intra-class Correlation Coefficients (ICC), root mean square error (RMSE), and percent error. Bland-Altman plots were used to assess bias. Excellent agreement was found between the IMU and criterion data for head ROM and peak rotational velocity (average ICC > 0.9). The trunk showed good agreement for most conditions (average ICC > 0.8). Average RMSE for both ROM (head = 2.64◦; trunk = 2.48◦) and peak rotational velocity (head = 11.76 ◦/s; trunk = 7.37 ◦/s) was low. The average percent error was below 5% for head and trunk ROM and peak rotational velocity. No clear pattern of bias was found for any measure across conditions. Findings suggest IMUs may provide a promising solution for estimating head and trunk movement, and a practical solution for tracking progression throughout rehabilitation or home exercise monitoring.

**Keywords:** concussion; inertial motion units (IMUs); vestibular exercises; validation; motion capture