# An Energy Efficiency Evaluation Method Based on Energy Baseline for Chemical Industry

Authors: Dong-mei Yao; Xin Zhang; Ke-feng Wang; Tao Zou; Dong Wang; Xin-hua Qian
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

According to the requirements and structure of ISO 50001 energy management system, this study proposes an energy efficiency evaluation method based on energy baseline for chemical industry. Using this method, the energy plan implementation effect in the processes of chemical production can be evaluated quantitatively, and evidences for system fault diagnosis can be provided. This method establishes the energy baseline models which can meet the demand of the different kinds of production processes and gives the general solving method of each kind of model according to the production data. Then the energy plan implementation effect can be evaluated and also whether the system is running normally can be determined through the baseline model. Finally, this method is used on cracked gas compressor unit of ethylene plant in some petrochemical enterprise; it can be proven that this method is correct and practical.


---

## Body

## 1. Introduction

According to the “Statistical Review of World Energy” which is revised annually by BP company, the world’s primary energy consumption structure for last ten years is shown in Table1. Currently, the total annual energy consumption of the world is equivalent to about 13 billion tons of oil, and average growth rate of energy consumption is 2% per year; the fossil energy such as oil, natural gas, and coal occupies more than 86% of total consumption which shows that the world energy structure is based on fossil fuels. At the same time, the problems of energy depletion and environmental pollution have become increasingly serious. It will take a long time to discover and develop new energies and renewable energies, so the fossil fuels will be the energy base of human survival and development in a long period and the problems of fossil fuels can only be alleviated by increasing the energy efficiency now.Table 1
The world’s primary energy consumption structure during last decade.
Year Total primary energy/million tons of oil equivalent Share of each energy in primary energy structure/% Oil Natural gas Coal Nuclear energy Hydropower Renewable energy 2005 10557.1 36.1 23.5 27.8 6.0 6.3 / 2006 10878.5 35.8 23.7 28.4 5.8 6.3 / 2007 11099.3 35.6 23.8 28.6 5.6 6.4 / 2008 11294.9 34.8 24.1 29.2 5.5 6.4 / 2009 11164.3 34.8 23.8 29.4 5.4 6.6 / 2010 12002.4 33.6 23.8 29.6 5.2 6.5 1.3 2011 12225.0 33.4 23.8 29.7 4.9 6.5 1.7 2012 12476.6 33.1 23.9 29.9 4.5 6.7 1.9 2013 12730.4 32.9 23.7 30.1 4.4 6.7 2.2 2014 12928.4 32.6 23.7 30.0 4.4 6.8 2.5Therefore, energy efficiency in chemical industry organizations has been controlled and systematized by using a powerful instrument which is called energy management system (EnMS). It can assist the organizations to develop and implement an energy policy and establish quality objectives, targets, and action plans which take into account legal requirements and information related to significant energy use, improve energy efficiency, reduce costs, and improve environmental protection and benefits at last. And the EnMS conforming to international standard is based on ISO 50001 standard. The ISO 50001 standard is drawn up and published by the ISO international organization for standardization ISO/PC 242, energy management committee. The ISO 50001 standard can guide the organization to establish EnMS. Some researches have given the methods about how to establish the EnMS framework based on ISO 50001 in kinds of industries specifically [1–5]. For example, an EnMS framework proposed by Vikhorev et al. could implement the standards for energy data exchange, on-line energy data analysis, performance measurement, and display of energy usage [1]. Giacone and Mancò regard the whole energy system of a site as a single matrix equation model, use the model to quantify the specific energy consumption, and provide support for the energy benchmarking, budgeting, and targeting; finally the method is used in the glass and cast iron melting processes successfully [2]. Chiu et al. study the case of an electronic parts and components company which has established an integration-energy-practice model for introducing an ISO 50001 energy management system; the results of this study prove that the model could efficiently meet demands for energy performance indicators [4].However, ISO 50001 is only a framework; the researches mentioned above are also the extension of this framework in detail, it is not methodology-based, it cannot implement quantitative evaluation of energy efficiency by ISO 50001 only, and this shortage restricts the effective application of EnMS in the complex chemical industry. So it needs to implement the quantitative evaluation of energy efficiency, the EnMS will be effective, and then the EnMS can be used to produce a comprehensive understanding of the overall energy situation, main problems, and the energy consumption and energy-saving potential and eventually improve energy efficiency and economic benefit. Through combining the ISO 50001 standard framework and chemical business characteristics [6], this research puts forward a quantitative evaluation method of energy efficiency based on energy baseline for chemical industry, where the energy baseline is a curve that can react to the energy performance of system or equipment; it is quantitative reference providing a basis for comparison of energy performance. Therefore, the key of this method is to establish the energy baseline mathematical model that can accurately reflect the actual energy performance in a specified period of time. Then it can implement the energy efficiency evaluation by comparing the differences between the actual energy performance and the energy performance calculated by energy baseline model.


## 2. Methods

### 2.1. ISO 50001 Standard Framework

The purpose of ISO 50001 standard is to provide a standard framework which assists organizations to establish the energy management system and necessary energy management process to improve energy performance. According to ISO 50001 requirements, the energy management system structure is shown in Figure1.Figure 1
Energy management system architecture.The ISO 50001 standard is based on Plan-Do-Check-Act (PDCA) continual improvement principle, the details of the respective steps are shown in Figure1, and each step is divided into technical level and management level as shown in Figure 2.Figure 2
Technical level and management level of ISO 50001 standard.This study focuses on the technical level of ISO 50001 standard. Based on the technical level framework, this research will center on energy baseline and discuss how to establish energy baseline model that is suitable for different kinds of chemical production processes. And then energy efficiency evaluation will be implemented by using the energy baseline model and the technical level of the EnMS is ensured eventually to run as PDCA continual improvement principle.


### 2.2. Establishment of Energy Baseline Model

#### 2.2.1. Energy Performance Indicators and Impact Factors

According to the results of the energy review, it is ensured that the energy performance can be measured and impacts the energy efficiency, energy use, and energy consumption. Then, according to the measured energy input parameters, the energy performance will further be expressed by energy performance indicators (EnPIs) that can illustrate whether the energy performance is good or not. Then, the impact factors affecting EnPIs are analyzed and selected.


#### 2.2.2. The Form of Energy Baseline Model

The EnPIs and impact factors can be expressed as the following function equation:(1)Y=fX,θ+ε,where Y=(Y1,Y2,…,YM)T is an M-dimensional column vector; it represents M EnPIs; X=(X1,X2,…,XN)T is an N-dimensional column vector; it represents N impact factors of the EnPIs; f(·)=[f1(·),f2(·),…,fM(·)]T is a mapping relation from X to Y; θ is a coefficient vector in the mapping relations; ε=(ε1,ε2,…,εM); εi (i=1,2,…,M) is a random variable that obeys the normal distribution whose mean is zero and variance is σ2.Depending on the different energy objects, the energy baseline models are divided into three types: (A) empirical equation or mechanistic model is given directly in the form of (1); (B) mechanistic model is difficult to be obtained and the energy object is linear; then the energy baseline model is defined as a linear model; (C) mechanistic model is difficult to be obtained and energy object shows a strong nonlinear character, the energy baseline model is defined as a sum of nonlinear functions of each impact factor, and each nonlinear function of each impact factor is approximated by B-spline difference function. In the following, the energy baseline model forms are given as mentioned by (B) and (C).(1) Linear Model. Each EnPI is expressed as linear form of impact factors, as shown in the following:(2)Y1=θ01+θ11X1+θ21X2+⋯+θN1XN+ε1,Y2=θ02+θ12X1+θ22X2+⋯+θN2XN+ε2,⋮YM=θ0M+θ1MX1+θ2MX2+⋯+θNMXN+εM.(2) B-Spline Approximation Model. For many chemical industry processes, it is difficult to give an accurate empirical or mechanistic functional relationship, especially when there is a strong nonlinear relationship between the EnPIs and the impact factors. Therefore, general function forms are proposed to approximate various nonlinear relationships.The EnPIs are the sum of nonlinear relationships of each impact factor, as shown in the following:(3)Y1=f11X1+f21X2+⋯+fN1XN+ε1,Y2=f12X1+f22X2+⋯+fN2XN+ε2,⋮YM=f1MX1+f2MX2+⋯+fNMXN+εN.In (3), each fij(Xj) is approximated by a T-order B-spline function. And T is adjusted to adapt to different nonlinear degree models: the higher the nonlinear degree of the models is, the bigger T is. It generally can meet most of the needs when T takes 3 [7–9].Firstly, eachXj is divided into Lj segments evenly, the length of each segment is hj, the segmented point values are denoted from small to large as ξj0,ξj2,…,ξjLj, each segmented point corresponds to a function value denoted as βj0,βj2,…,βjLj, and these function values are the same as the parameters that need to be solved in the B-spline function. So hj, Lj, and the segmented point values are consistent with the relationship as shown in the following:(4)hj=max⁡Xj-min⁡XjLj,ξjl+1=hj+ξjl.The approximate function of B-spline forfij(Xj) in the lth segment range is as shown in the following equation:(5)fij,lXj=∑t=0Tβjl+t-3·Gt,TXj-ξjlhj,where i=1,2,…,N, j=1,2,…,M, l=3,4,…,(Lj-1), ξjl≤Xj≤ξj(l+1), and Gt,T((Xj-ξjl)/hj) is basis function of T-order B-spline function; make x=(Xj-ξjl)/hj, x∈[0,1]; then (6) can be deducted:(6)Gt,Tx=1T!∑s=0T-t-1sCT+1sx+T-t-sT.The recursive equation is shown as follows:(7)Gt,0x=10≤x≤10otherwise,Gt,Tx=1Tt+1-xGt,T-1x+1TT-t+xGt-1,T-1x.


#### 2.2.3. Model Parameter Solution

The form of the energy baseline model is obtained in Section4; it needs to solve the model parameters according to the energy data. K group of energy data is given: (Y1k,Y2k,…,YMk,X1k,X2k,…,XNk), k=1,2,…,K. It can form M·K algebraic equations about the model parameters of θ by putting the energy data into the energy baseline model as shown in the following:(8)Y11=f1X11,X21,…,XN1,θ+ε11,Y21=f2X11,X21,…,XN1,θ+ε21,⋮YM1=fMX11,X21,…,XN1,θ+εM1,Y12=f1X12,X22,…,XN2,θ+ε12,⋮YMK=fMX1K,X2K,…,XNK,θ+εMK.The model parameter value is obtained by minimizing the objective function in the following equation:(9)Jθ=12∑k=1K∑m=1MYmk-fmX1k,X2k,…,XNk,θ-εmk2.Equation (9) can be simplified as M·K sum of error squares of θ as shown in the following equation:(10)Jθ=12∑j=1K·Mej2θ.Generally, (10) is a differentiable function of θ; then the necessary condition for J(θ) to obtain the minimum value is shown in the following equation:(11)∂Jθr∂θr=∑j=1K·Mejθejθ∂θr=0,where θ=(θ1,θ2,…,θR)T, r=1,2,…,R, and these two cases will be discussed according to the different forms of ej(θ) of θ.(1) θ as a Linear Relationship. In this case, ej(θ) is a linear expression of θ as shown in the following equation:(12)ejθ=Aj·θ+aj0+εj,where ej(θ)/∂θr=ajr; Aj=(aj1,aj2,…,ajR).Then (11) is a linear equation group as shown in the following equation:(13)ATAR×R·θ=ATb+ε,where (14)A=A1A2⋮AKMK·M×R,b=a10a20⋮aK·M0K·M×1,ε=ε1ε2⋮εK·MK·M×1.Generally,A is nonsingular; then θ has a unique solution. Because ε is an unknown random perturbation and there is (15), that is,(15)Eθ=EATA-1ATb-ε=EATA-1ATb-ATA-1ATEε=ATA-1ATb.Soθ^ which is the average value of θ is the estimate value of θ as shown in the following equation:(16)θ^=ATA-1ATb.The linear energy baseline model in (2) and the B-spline approximation model in (3) and (5) are both linear models of the parameters, so the models are solved by this method.(2) θ as a Nonlinear Relationship. When fm(X1k,X2k,…,XNk,θ) is a nonlinear expression of θ, (11) is nonlinear equation group, the solution would not be unique, there would be many pole points in (10), and the solution of (11) would not ensure making (10) minimized. Therefore, (10) is minimized directly. And the method combining PSO (Particle Swarm Optimization) and Gauss-Newton algorithm is proposed to solve this problem.PSO algorithm is an evolutionary computation method based on the movement of swarms looking for the most fertile feeding location. It uses a population-based, self-adaptive search optimization technique. Therefore, the global best solution can be found by simply adjusting the trajectory of each individual toward its own best location and toward the best particle of the entire swarm at each time step. So the PSO method is widely used due to its characteristics of global convergence and simplicity [10–12]. However, the convergence speed of PSO is slowed down soon when it is closed to the global optimal point. So the search is popular about the solution of this problem. A hybrid optimization algorithm (HOA) is put forward to deal with delamination identification in laminated composite structures. The HOA is based on combined PSO and SM (simplex method) [10]. Ratnaweera et al. discussed two methods to improve the performance of the PSO based on time-varying acceleration coefficients (TVAC): first by adding a small perturbation to a randomly selected modulus of the velocity vector of a random particle by predefined probability and second by discussing a kind of self-organizing hierarchical particle swarm optimizer with TVAC [11].Gauss-Newton algorithm is a well-known classical iterative algorithm to solve nonlinear least squares problems [13]. It uses Taylor series expansion to approximate the nonlinear regression model. Through iteration to correct the regression coefficient of the nonlinear regression model continuously, the regression coefficient will approximate the optimum regression coefficient of the nonlinear regression model continuously. Finally the sum of squared residuals to the original model is minimized. Meanwhile Gauss-Newton algorithm is local convergence, and the accuracy and speed of the algorithm depend on selection of the initial value to a large extent.So, according to the characteristics of the PSO and Gauss-Newton algorithm, the combination of the two methods is proposed to solve above nonlinear optimization problem [14]. Firstly, PSO is used to provide a rough solution. Then the rough solution is used as the iterative initial value of Gauss-Newton algorithm; the exact minimum value of (10) will be obtained by finite iterations. The steps using the PSO algorithm to acquire the rough solution are as follows.(1) Parameter Setting. The dimension of PSO is R, which is the number of the model parameters. The number of particles is D; the upper and lower bounds of the particle’s positions are (θrmin,θrmax); the upper and lower limits of particle velocities are (vrmin,vrmax), r=1,2,…,R. The position of d particle is expressed as θd=(θ1d,θ2d,…,θRd), the particle velocity is expressed as vd=(v1d,v2d,…,vRd), the history optimal position of the particle d is expressed as pb, and the optimal position of all the group is expressed as gb, d=1,2,…,D. The maximum weight is wmax, the minimum weight is wmin, and the maximum number of iterations is imax.(2) Initialization. Make the number of iterations i=0; the parameters are initialized according to the following:(17)θrd0=θrmin+θrmax-θrmin·rand1,vrd0=vrmin+vrmax-vrmin·rand1,pb=gb=θ10,where rand(1) is a function for generating a random value that meets standard normal distribution: r=1,2,…,R; d=1,2,…,D.(3) Particle Evaluation. Using (9) as the fitness function, the fitness of each particle can be calculated: Jd(i)=Jθd(i). If Jd(i)<J(pb), the history optimal position of particle pb is updated as pb=θd(i); if Jd(i)<J(gb), then the optimal position of all the particles is updated as gb=θd(i).(4) Particles Update. The position and speed of the particle are updated by the following:(18)vrdi+1=w·vvrdi+c1·rand1·pb-θrdi+c2·rand1·gb-θrdi,θrdi+1=θrdi+vrdi+1,w=wmax-iimaxwmax-wmin.(5) Stop Condition. When the fitness function reaches the setting accuracy or i reaches the maximum number of steps imax, the iteration stops, and the result of θ is gb.The rough solutiongb is within the smaller neighborhood of the exact optimal solution. Using gb as the iterative initial value of Gauss-Newton algorithm, the solution steps are as follows.(1) Make the iteration numberi=0, select the initial value point θ(0)=gb, and calculate the objective function value as shown in the following equation:(19)Jθ0=12∑j=1K·Mej2θ0.(2) Calculategjr(i)=∂ej[θ(i)]/∂θr, j=1,2,…,(K·M), and r=1,2,…,R; then gjr(i) can form the matrix G(i) as shown in the following:(20)Gi=g11ig12i⋯g1Rig21ig22i⋯g2Ri⋮⋮⋱⋮gK·M1igK·M2i⋯gK·MRiK·M×R.The gradient vector ofJ(θ) is b(i)=[b1(i),b2(i),…,bR(i)]T, where(21)bri=∑j=1K·Mejθi∂ejθi∂θr.(3) The linear equationsGT(i)·G(i)p=-b(i) of vector p are solved, where if the rank of GT(i)·G(i) is R, the equation has the unique solution p(i)=-GT(i)·G(i)-1·b(i); if the rank of GT(i)·G(i) is less than R, p(i)=-b(i).(4) Minimize the linear function groupFθ(i)+λ·p(i) of the single variable λ. The direct linear search method is used to obtain λ(k), and λ(k) meet Fθ(i)+λ(i)·p(i)=minλ⁡Fθ(i)+λ·p(i); then the new iteration point and the target function value are as shown in the following:(22)θi+1=θi+λi·pi,Jθi+1=12∑j=1K·Mej2θi+1.(5) If the termination condition is satisfied, the iteration will stop, andθ^=θ(i+1) is the final solution; otherwise make i=i+1, and return b step to continue the iteration.


#### 2.2.4. Energy Baseline Model Checking

Through the above processes, the final energy baseline model is followed:(23)Y=fX,θ^.The adjusted decision coefficient is used as the evaluation criteria to check the obtained model.Firstly, the estimated values of EnPI are calculated on each measured point:(24)Y^mk=fX1k,X2k,…,XNk,θ^,k=1,2,…,K,m=1,2,…,M.Make(25)y=Y11,Y21,…,YM1,…,Y1K,Y2K,…,YMK,y^=Y^11,Y^21,…,Y^M1,…,Y^1K,Y^2K,…,Y^MK;the following statistics can be calculated.Total Square Sum is given as follows:(26)TSS=y-y-Ty-y-.Explain Square Sum is given as follows:(27)ESS=y^-y-Ty^-y-.Then the adjusted decision coefficient is as shown in the following equation:(28)R2¯=1-K×M-1K×M-R-11-ESSTSS2,where R2¯ is adjusted decision coefficient. The value range of R2¯ is [0,1]; the more close to 1, the better to express the relationship between the EnPI and the impact factors for the energy baseline model. The threshold can be set in advance according to the need; if the adjusted decision coefficient is greater than the threshold, it is considered that the model fitting is good; otherwise it is considered that the model fitting is poor; the model structure or impact factors need to be reselected.


### 2.3. Energy Efficiency Evaluation Based on Energy Baseline

According to the above sessions, the energy baseline can implement two functions in EnMS:(1)
Energy baseline can be used as energy performance evaluation benchmark after the implementation of energy plan.(2)
Energy baseline can be used as energy performance monitoring and system fault diagnosis on the base of its real-time estimation of energy performance indicators.


#### 2.3.1. Statistics Related to EnPIs

Firstly, three benchmark statistics are provided to judge whether energy performance achieves the desired standard:(1)
Expected reduced amplitude of EnPIs:eD=(e1D,e2D,…,eMD)(2)
The mean sample deviation of energy baseline model:e-B=(e-1B,e-2B,…,e-MB)(29)e-mB=1K∑k=1KYmk-Y^mk.(3)
The standard deviation of sample deviation of energy baseline model:(30)σB=σ1B,σ2B,…,σMB,σmB=1K-1∑k=1KYmk-Y^mk-e-mB2,m=1,2,…,M.Then a period of time to be evaluated is defined, and then the measured values of EnPIs and impact factors are(y1p,y2p,…,yMp,x1p,x2p,…,xNp), p=1,2,…,P. The measured values of EnPIs can be expressed as the following form: (31)ym=ym1,ym2,…,ymPT,m=1,2,…,M.Then the impact factors values are substituted into (21) of energy baseline model in order to obtain the estimated values of EnPIs: (32)y^m=y^m1,y^m2,…,y^mPT.Defineem=(em1,em2,…,emP)T=ym-y^m; then the mean value and standard deviation of em are (33)e-m=1P∑p=1Pemp,σm=1P-1∑p=1Pemp2.


#### 2.3.2. Energy Performance Evaluation

In order to evaluate the energy performance during a period of time after the implementation of the energy plan, the energy baseline model is established firstly according to the data before the implementation of the energy plan; then evaluated values of EnPIs through the model in this period are obtained. Through the comparison between the estimated values and the measured values, the implementation effect of the energy plan can be estimated quantitatively. Specifically, for eachYm (m=1,2,…,M), e-m is calculated through (33) and compared with emD. The evaluation criteria are as follows:(1)
If(-e-m)≥emD, the energy management system achieves or exceeds expected effect; it meets the requirements of indicators.(2)
If(-e-m)<emD, the actual reduced amplitude of EnPI is less than expected value, and it does not meet the requirements; the enterprise needs to review its energy again and make a new plan.


#### 2.3.3. Energy Performance Monitoring and System Fault Diagnosis

Energy efficiency monitoring and system fault diagnosis method based on energy baseline are used in the checking period of energy management system. The energy baseline model reflects the current state of the system operation; whether the system is normal operation can be judged by it. Specifically, for eachYm (m=1,2,…,M), e-m is calculated through (33), and the judging criteria are as follows:(1)
IfemB-3σmD≤e-m≤emB+3σmD, the energy management system runs normally.(2)
Ife-m<emB-3σmD, the energy baseline model is no longer suitable for estimating the current state of EnPI; the model needs to be rebuilt according to the current data.(3)
Ife-m>emB+3σmD, the energy consumption is too high or the system is failure; enterprise needs to make a new plan or system needs maintenance.


### 2.4. PDCA Principle Implementation of Energy Baseline

The EnMS is based on Plan-Do-Check-Act (PDCA) continual improvement principle, and the energy baseline will be adjusted in different phases to implement different functions:(1)
In Plan phase, energy baseline is established firstly based on energy history information and is as the benchmark of the next phase.(2)
In Do phase, the energy baseline is used to evaluate the difference of the energy performance before and after the implementation of the energy plan.(3)
In Check phase, the whole system runs smoothly, the energy baseline needs to be readjusted based on the energy data after the Do phase, and, at this time, the energy baseline can be used to monitor energy efficiency or diagnose system fault.(4)
In Action phase, the energy efficiency needs to be further improved to take actions; then it goes into the Plan phase; the present energy baseline is as the benchmark for next phase, and the EnMS forms a closed loop.


## 2.1. ISO 50001 Standard Framework

The purpose of ISO 50001 standard is to provide a standard framework which assists organizations to establish the energy management system and necessary energy management process to improve energy performance. According to ISO 50001 requirements, the energy management system structure is shown in Figure1.Figure 1
Energy management system architecture.The ISO 50001 standard is based on Plan-Do-Check-Act (PDCA) continual improvement principle, the details of the respective steps are shown in Figure1, and each step is divided into technical level and management level as shown in Figure 2.Figure 2
Technical level and management level of ISO 50001 standard.This study focuses on the technical level of ISO 50001 standard. Based on the technical level framework, this research will center on energy baseline and discuss how to establish energy baseline model that is suitable for different kinds of chemical production processes. And then energy efficiency evaluation will be implemented by using the energy baseline model and the technical level of the EnMS is ensured eventually to run as PDCA continual improvement principle.


## 2.2. Establishment of Energy Baseline Model

### 2.2.1. Energy Performance Indicators and Impact Factors

According to the results of the energy review, it is ensured that the energy performance can be measured and impacts the energy efficiency, energy use, and energy consumption. Then, according to the measured energy input parameters, the energy performance will further be expressed by energy performance indicators (EnPIs) that can illustrate whether the energy performance is good or not. Then, the impact factors affecting EnPIs are analyzed and selected.


### 2.2.2. The Form of Energy Baseline Model

The EnPIs and impact factors can be expressed as the following function equation:(1)Y=fX,θ+ε,where Y=(Y1,Y2,…,YM)T is an M-dimensional column vector; it represents M EnPIs; X=(X1,X2,…,XN)T is an N-dimensional column vector; it represents N impact factors of the EnPIs; f(·)=[f1(·),f2(·),…,fM(·)]T is a mapping relation from X to Y; θ is a coefficient vector in the mapping relations; ε=(ε1,ε2,…,εM); εi (i=1,2,…,M) is a random variable that obeys the normal distribution whose mean is zero and variance is σ2.Depending on the different energy objects, the energy baseline models are divided into three types: (A) empirical equation or mechanistic model is given directly in the form of (1); (B) mechanistic model is difficult to be obtained and the energy object is linear; then the energy baseline model is defined as a linear model; (C) mechanistic model is difficult to be obtained and energy object shows a strong nonlinear character, the energy baseline model is defined as a sum of nonlinear functions of each impact factor, and each nonlinear function of each impact factor is approximated by B-spline difference function. In the following, the energy baseline model forms are given as mentioned by (B) and (C).(1) Linear Model. Each EnPI is expressed as linear form of impact factors, as shown in the following:(2)Y1=θ01+θ11X1+θ21X2+⋯+θN1XN+ε1,Y2=θ02+θ12X1+θ22X2+⋯+θN2XN+ε2,⋮YM=θ0M+θ1MX1+θ2MX2+⋯+θNMXN+εM.(2) B-Spline Approximation Model. For many chemical industry processes, it is difficult to give an accurate empirical or mechanistic functional relationship, especially when there is a strong nonlinear relationship between the EnPIs and the impact factors. Therefore, general function forms are proposed to approximate various nonlinear relationships.The EnPIs are the sum of nonlinear relationships of each impact factor, as shown in the following:(3)Y1=f11X1+f21X2+⋯+fN1XN+ε1,Y2=f12X1+f22X2+⋯+fN2XN+ε2,⋮YM=f1MX1+f2MX2+⋯+fNMXN+εN.In (3), each fij(Xj) is approximated by a T-order B-spline function. And T is adjusted to adapt to different nonlinear degree models: the higher the nonlinear degree of the models is, the bigger T is. It generally can meet most of the needs when T takes 3 [7–9].Firstly, eachXj is divided into Lj segments evenly, the length of each segment is hj, the segmented point values are denoted from small to large as ξj0,ξj2,…,ξjLj, each segmented point corresponds to a function value denoted as βj0,βj2,…,βjLj, and these function values are the same as the parameters that need to be solved in the B-spline function. So hj, Lj, and the segmented point values are consistent with the relationship as shown in the following:(4)hj=max⁡Xj-min⁡XjLj,ξjl+1=hj+ξjl.The approximate function of B-spline forfij(Xj) in the lth segment range is as shown in the following equation:(5)fij,lXj=∑t=0Tβjl+t-3·Gt,TXj-ξjlhj,where i=1,2,…,N, j=1,2,…,M, l=3,4,…,(Lj-1), ξjl≤Xj≤ξj(l+1), and Gt,T((Xj-ξjl)/hj) is basis function of T-order B-spline function; make x=(Xj-ξjl)/hj, x∈[0,1]; then (6) can be deducted:(6)Gt,Tx=1T!∑s=0T-t-1sCT+1sx+T-t-sT.The recursive equation is shown as follows:(7)Gt,0x=10≤x≤10otherwise,Gt,Tx=1Tt+1-xGt,T-1x+1TT-t+xGt-1,T-1x.


### 2.2.3. Model Parameter Solution

The form of the energy baseline model is obtained in Section4; it needs to solve the model parameters according to the energy data. K group of energy data is given: (Y1k,Y2k,…,YMk,X1k,X2k,…,XNk), k=1,2,…,K. It can form M·K algebraic equations about the model parameters of θ by putting the energy data into the energy baseline model as shown in the following:(8)Y11=f1X11,X21,…,XN1,θ+ε11,Y21=f2X11,X21,…,XN1,θ+ε21,⋮YM1=fMX11,X21,…,XN1,θ+εM1,Y12=f1X12,X22,…,XN2,θ+ε12,⋮YMK=fMX1K,X2K,…,XNK,θ+εMK.The model parameter value is obtained by minimizing the objective function in the following equation:(9)Jθ=12∑k=1K∑m=1MYmk-fmX1k,X2k,…,XNk,θ-εmk2.Equation (9) can be simplified as M·K sum of error squares of θ as shown in the following equation:(10)Jθ=12∑j=1K·Mej2θ.Generally, (10) is a differentiable function of θ; then the necessary condition for J(θ) to obtain the minimum value is shown in the following equation:(11)∂Jθr∂θr=∑j=1K·Mejθejθ∂θr=0,where θ=(θ1,θ2,…,θR)T, r=1,2,…,R, and these two cases will be discussed according to the different forms of ej(θ) of θ.(1) θ as a Linear Relationship. In this case, ej(θ) is a linear expression of θ as shown in the following equation:(12)ejθ=Aj·θ+aj0+εj,where ej(θ)/∂θr=ajr; Aj=(aj1,aj2,…,ajR).Then (11) is a linear equation group as shown in the following equation:(13)ATAR×R·θ=ATb+ε,where (14)A=A1A2⋮AKMK·M×R,b=a10a20⋮aK·M0K·M×1,ε=ε1ε2⋮εK·MK·M×1.Generally,A is nonsingular; then θ has a unique solution. Because ε is an unknown random perturbation and there is (15), that is,(15)Eθ=EATA-1ATb-ε=EATA-1ATb-ATA-1ATEε=ATA-1ATb.Soθ^ which is the average value of θ is the estimate value of θ as shown in the following equation:(16)θ^=ATA-1ATb.The linear energy baseline model in (2) and the B-spline approximation model in (3) and (5) are both linear models of the parameters, so the models are solved by this method.(2) θ as a Nonlinear Relationship. When fm(X1k,X2k,…,XNk,θ) is a nonlinear expression of θ, (11) is nonlinear equation group, the solution would not be unique, there would be many pole points in (10), and the solution of (11) would not ensure making (10) minimized. Therefore, (10) is minimized directly. And the method combining PSO (Particle Swarm Optimization) and Gauss-Newton algorithm is proposed to solve this problem.PSO algorithm is an evolutionary computation method based on the movement of swarms looking for the most fertile feeding location. It uses a population-based, self-adaptive search optimization technique. Therefore, the global best solution can be found by simply adjusting the trajectory of each individual toward its own best location and toward the best particle of the entire swarm at each time step. So the PSO method is widely used due to its characteristics of global convergence and simplicity [10–12]. However, the convergence speed of PSO is slowed down soon when it is closed to the global optimal point. So the search is popular about the solution of this problem. A hybrid optimization algorithm (HOA) is put forward to deal with delamination identification in laminated composite structures. The HOA is based on combined PSO and SM (simplex method) [10]. Ratnaweera et al. discussed two methods to improve the performance of the PSO based on time-varying acceleration coefficients (TVAC): first by adding a small perturbation to a randomly selected modulus of the velocity vector of a random particle by predefined probability and second by discussing a kind of self-organizing hierarchical particle swarm optimizer with TVAC [11].Gauss-Newton algorithm is a well-known classical iterative algorithm to solve nonlinear least squares problems [13]. It uses Taylor series expansion to approximate the nonlinear regression model. Through iteration to correct the regression coefficient of the nonlinear regression model continuously, the regression coefficient will approximate the optimum regression coefficient of the nonlinear regression model continuously. Finally the sum of squared residuals to the original model is minimized. Meanwhile Gauss-Newton algorithm is local convergence, and the accuracy and speed of the algorithm depend on selection of the initial value to a large extent.So, according to the characteristics of the PSO and Gauss-Newton algorithm, the combination of the two methods is proposed to solve above nonlinear optimization problem [14]. Firstly, PSO is used to provide a rough solution. Then the rough solution is used as the iterative initial value of Gauss-Newton algorithm; the exact minimum value of (10) will be obtained by finite iterations. The steps using the PSO algorithm to acquire the rough solution are as follows.(1) Parameter Setting. The dimension of PSO is R, which is the number of the model parameters. The number of particles is D; the upper and lower bounds of the particle’s positions are (θrmin,θrmax); the upper and lower limits of particle velocities are (vrmin,vrmax), r=1,2,…,R. The position of d particle is expressed as θd=(θ1d,θ2d,…,θRd), the particle velocity is expressed as vd=(v1d,v2d,…,vRd), the history optimal position of the particle d is expressed as pb, and the optimal position of all the group is expressed as gb, d=1,2,…,D. The maximum weight is wmax, the minimum weight is wmin, and the maximum number of iterations is imax.(2) Initialization. Make the number of iterations i=0; the parameters are initialized according to the following:(17)θrd0=θrmin+θrmax-θrmin·rand1,vrd0=vrmin+vrmax-vrmin·rand1,pb=gb=θ10,where rand(1) is a function for generating a random value that meets standard normal distribution: r=1,2,…,R; d=1,2,…,D.(3) Particle Evaluation. Using (9) as the fitness function, the fitness of each particle can be calculated: Jd(i)=Jθd(i). If Jd(i)<J(pb), the history optimal position of particle pb is updated as pb=θd(i); if Jd(i)<J(gb), then the optimal position of all the particles is updated as gb=θd(i).(4) Particles Update. The position and speed of the particle are updated by the following:(18)vrdi+1=w·vvrdi+c1·rand1·pb-θrdi+c2·rand1·gb-θrdi,θrdi+1=θrdi+vrdi+1,w=wmax-iimaxwmax-wmin.(5) Stop Condition. When the fitness function reaches the setting accuracy or i reaches the maximum number of steps imax, the iteration stops, and the result of θ is gb.The rough solutiongb is within the smaller neighborhood of the exact optimal solution. Using gb as the iterative initial value of Gauss-Newton algorithm, the solution steps are as follows.(1) Make the iteration numberi=0, select the initial value point θ(0)=gb, and calculate the objective function value as shown in the following equation:(19)Jθ0=12∑j=1K·Mej2θ0.(2) Calculategjr(i)=∂ej[θ(i)]/∂θr, j=1,2,…,(K·M), and r=1,2,…,R; then gjr(i) can form the matrix G(i) as shown in the following:(20)Gi=g11ig12i⋯g1Rig21ig22i⋯g2Ri⋮⋮⋱⋮gK·M1igK·M2i⋯gK·MRiK·M×R.The gradient vector ofJ(θ) is b(i)=[b1(i),b2(i),…,bR(i)]T, where(21)bri=∑j=1K·Mejθi∂ejθi∂θr.(3) The linear equationsGT(i)·G(i)p=-b(i) of vector p are solved, where if the rank of GT(i)·G(i) is R, the equation has the unique solution p(i)=-GT(i)·G(i)-1·b(i); if the rank of GT(i)·G(i) is less than R, p(i)=-b(i).(4) Minimize the linear function groupFθ(i)+λ·p(i) of the single variable λ. The direct linear search method is used to obtain λ(k), and λ(k) meet Fθ(i)+λ(i)·p(i)=minλ⁡Fθ(i)+λ·p(i); then the new iteration point and the target function value are as shown in the following:(22)θi+1=θi+λi·pi,Jθi+1=12∑j=1K·Mej2θi+1.(5) If the termination condition is satisfied, the iteration will stop, andθ^=θ(i+1) is the final solution; otherwise make i=i+1, and return b step to continue the iteration.


### 2.2.4. Energy Baseline Model Checking

Through the above processes, the final energy baseline model is followed:(23)Y=fX,θ^.The adjusted decision coefficient is used as the evaluation criteria to check the obtained model.Firstly, the estimated values of EnPI are calculated on each measured point:(24)Y^mk=fX1k,X2k,…,XNk,θ^,k=1,2,…,K,m=1,2,…,M.Make(25)y=Y11,Y21,…,YM1,…,Y1K,Y2K,…,YMK,y^=Y^11,Y^21,…,Y^M1,…,Y^1K,Y^2K,…,Y^MK;the following statistics can be calculated.Total Square Sum is given as follows:(26)TSS=y-y-Ty-y-.Explain Square Sum is given as follows:(27)ESS=y^-y-Ty^-y-.Then the adjusted decision coefficient is as shown in the following equation:(28)R2¯=1-K×M-1K×M-R-11-ESSTSS2,where R2¯ is adjusted decision coefficient. The value range of R2¯ is [0,1]; the more close to 1, the better to express the relationship between the EnPI and the impact factors for the energy baseline model. The threshold can be set in advance according to the need; if the adjusted decision coefficient is greater than the threshold, it is considered that the model fitting is good; otherwise it is considered that the model fitting is poor; the model structure or impact factors need to be reselected.


## 2.2.1. Energy Performance Indicators and Impact Factors

According to the results of the energy review, it is ensured that the energy performance can be measured and impacts the energy efficiency, energy use, and energy consumption. Then, according to the measured energy input parameters, the energy performance will further be expressed by energy performance indicators (EnPIs) that can illustrate whether the energy performance is good or not. Then, the impact factors affecting EnPIs are analyzed and selected.


## 2.2.2. The Form of Energy Baseline Model

The EnPIs and impact factors can be expressed as the following function equation:(1)Y=fX,θ+ε,where Y=(Y1,Y2,…,YM)T is an M-dimensional column vector; it represents M EnPIs; X=(X1,X2,…,XN)T is an N-dimensional column vector; it represents N impact factors of the EnPIs; f(·)=[f1(·),f2(·),…,fM(·)]T is a mapping relation from X to Y; θ is a coefficient vector in the mapping relations; ε=(ε1,ε2,…,εM); εi (i=1,2,…,M) is a random variable that obeys the normal distribution whose mean is zero and variance is σ2.Depending on the different energy objects, the energy baseline models are divided into three types: (A) empirical equation or mechanistic model is given directly in the form of (1); (B) mechanistic model is difficult to be obtained and the energy object is linear; then the energy baseline model is defined as a linear model; (C) mechanistic model is difficult to be obtained and energy object shows a strong nonlinear character, the energy baseline model is defined as a sum of nonlinear functions of each impact factor, and each nonlinear function of each impact factor is approximated by B-spline difference function. In the following, the energy baseline model forms are given as mentioned by (B) and (C).(1) Linear Model. Each EnPI is expressed as linear form of impact factors, as shown in the following:(2)Y1=θ01+θ11X1+θ21X2+⋯+θN1XN+ε1,Y2=θ02+θ12X1+θ22X2+⋯+θN2XN+ε2,⋮YM=θ0M+θ1MX1+θ2MX2+⋯+θNMXN+εM.(2) B-Spline Approximation Model. For many chemical industry processes, it is difficult to give an accurate empirical or mechanistic functional relationship, especially when there is a strong nonlinear relationship between the EnPIs and the impact factors. Therefore, general function forms are proposed to approximate various nonlinear relationships.The EnPIs are the sum of nonlinear relationships of each impact factor, as shown in the following:(3)Y1=f11X1+f21X2+⋯+fN1XN+ε1,Y2=f12X1+f22X2+⋯+fN2XN+ε2,⋮YM=f1MX1+f2MX2+⋯+fNMXN+εN.In (3), each fij(Xj) is approximated by a T-order B-spline function. And T is adjusted to adapt to different nonlinear degree models: the higher the nonlinear degree of the models is, the bigger T is. It generally can meet most of the needs when T takes 3 [7–9].Firstly, eachXj is divided into Lj segments evenly, the length of each segment is hj, the segmented point values are denoted from small to large as ξj0,ξj2,…,ξjLj, each segmented point corresponds to a function value denoted as βj0,βj2,…,βjLj, and these function values are the same as the parameters that need to be solved in the B-spline function. So hj, Lj, and the segmented point values are consistent with the relationship as shown in the following:(4)hj=max⁡Xj-min⁡XjLj,ξjl+1=hj+ξjl.The approximate function of B-spline forfij(Xj) in the lth segment range is as shown in the following equation:(5)fij,lXj=∑t=0Tβjl+t-3·Gt,TXj-ξjlhj,where i=1,2,…,N, j=1,2,…,M, l=3,4,…,(Lj-1), ξjl≤Xj≤ξj(l+1), and Gt,T((Xj-ξjl)/hj) is basis function of T-order B-spline function; make x=(Xj-ξjl)/hj, x∈[0,1]; then (6) can be deducted:(6)Gt,Tx=1T!∑s=0T-t-1sCT+1sx+T-t-sT.The recursive equation is shown as follows:(7)Gt,0x=10≤x≤10otherwise,Gt,Tx=1Tt+1-xGt,T-1x+1TT-t+xGt-1,T-1x.


## 2.2.3. Model Parameter Solution

The form of the energy baseline model is obtained in Section4; it needs to solve the model parameters according to the energy data. K group of energy data is given: (Y1k,Y2k,…,YMk,X1k,X2k,…,XNk), k=1,2,…,K. It can form M·K algebraic equations about the model parameters of θ by putting the energy data into the energy baseline model as shown in the following:(8)Y11=f1X11,X21,…,XN1,θ+ε11,Y21=f2X11,X21,…,XN1,θ+ε21,⋮YM1=fMX11,X21,…,XN1,θ+εM1,Y12=f1X12,X22,…,XN2,θ+ε12,⋮YMK=fMX1K,X2K,…,XNK,θ+εMK.The model parameter value is obtained by minimizing the objective function in the following equation:(9)Jθ=12∑k=1K∑m=1MYmk-fmX1k,X2k,…,XNk,θ-εmk2.Equation (9) can be simplified as M·K sum of error squares of θ as shown in the following equation:(10)Jθ=12∑j=1K·Mej2θ.Generally, (10) is a differentiable function of θ; then the necessary condition for J(θ) to obtain the minimum value is shown in the following equation:(11)∂Jθr∂θr=∑j=1K·Mejθejθ∂θr=0,where θ=(θ1,θ2,…,θR)T, r=1,2,…,R, and these two cases will be discussed according to the different forms of ej(θ) of θ.(1) θ as a Linear Relationship. In this case, ej(θ) is a linear expression of θ as shown in the following equation:(12)ejθ=Aj·θ+aj0+εj,where ej(θ)/∂θr=ajr; Aj=(aj1,aj2,…,ajR).Then (11) is a linear equation group as shown in the following equation:(13)ATAR×R·θ=ATb+ε,where (14)A=A1A2⋮AKMK·M×R,b=a10a20⋮aK·M0K·M×1,ε=ε1ε2⋮εK·MK·M×1.Generally,A is nonsingular; then θ has a unique solution. Because ε is an unknown random perturbation and there is (15), that is,(15)Eθ=EATA-1ATb-ε=EATA-1ATb-ATA-1ATEε=ATA-1ATb.Soθ^ which is the average value of θ is the estimate value of θ as shown in the following equation:(16)θ^=ATA-1ATb.The linear energy baseline model in (2) and the B-spline approximation model in (3) and (5) are both linear models of the parameters, so the models are solved by this method.(2) θ as a Nonlinear Relationship. When fm(X1k,X2k,…,XNk,θ) is a nonlinear expression of θ, (11) is nonlinear equation group, the solution would not be unique, there would be many pole points in (10), and the solution of (11) would not ensure making (10) minimized. Therefore, (10) is minimized directly. And the method combining PSO (Particle Swarm Optimization) and Gauss-Newton algorithm is proposed to solve this problem.PSO algorithm is an evolutionary computation method based on the movement of swarms looking for the most fertile feeding location. It uses a population-based, self-adaptive search optimization technique. Therefore, the global best solution can be found by simply adjusting the trajectory of each individual toward its own best location and toward the best particle of the entire swarm at each time step. So the PSO method is widely used due to its characteristics of global convergence and simplicity [10–12]. However, the convergence speed of PSO is slowed down soon when it is closed to the global optimal point. So the search is popular about the solution of this problem. A hybrid optimization algorithm (HOA) is put forward to deal with delamination identification in laminated composite structures. The HOA is based on combined PSO and SM (simplex method) [10]. Ratnaweera et al. discussed two methods to improve the performance of the PSO based on time-varying acceleration coefficients (TVAC): first by adding a small perturbation to a randomly selected modulus of the velocity vector of a random particle by predefined probability and second by discussing a kind of self-organizing hierarchical particle swarm optimizer with TVAC [11].Gauss-Newton algorithm is a well-known classical iterative algorithm to solve nonlinear least squares problems [13]. It uses Taylor series expansion to approximate the nonlinear regression model. Through iteration to correct the regression coefficient of the nonlinear regression model continuously, the regression coefficient will approximate the optimum regression coefficient of the nonlinear regression model continuously. Finally the sum of squared residuals to the original model is minimized. Meanwhile Gauss-Newton algorithm is local convergence, and the accuracy and speed of the algorithm depend on selection of the initial value to a large extent.So, according to the characteristics of the PSO and Gauss-Newton algorithm, the combination of the two methods is proposed to solve above nonlinear optimization problem [14]. Firstly, PSO is used to provide a rough solution. Then the rough solution is used as the iterative initial value of Gauss-Newton algorithm; the exact minimum value of (10) will be obtained by finite iterations. The steps using the PSO algorithm to acquire the rough solution are as follows.(1) Parameter Setting. The dimension of PSO is R, which is the number of the model parameters. The number of particles is D; the upper and lower bounds of the particle’s positions are (θrmin,θrmax); the upper and lower limits of particle velocities are (vrmin,vrmax), r=1,2,…,R. The position of d particle is expressed as θd=(θ1d,θ2d,…,θRd), the particle velocity is expressed as vd=(v1d,v2d,…,vRd), the history optimal position of the particle d is expressed as pb, and the optimal position of all the group is expressed as gb, d=1,2,…,D. The maximum weight is wmax, the minimum weight is wmin, and the maximum number of iterations is imax.(2) Initialization. Make the number of iterations i=0; the parameters are initialized according to the following:(17)θrd0=θrmin+θrmax-θrmin·rand1,vrd0=vrmin+vrmax-vrmin·rand1,pb=gb=θ10,where rand(1) is a function for generating a random value that meets standard normal distribution: r=1,2,…,R; d=1,2,…,D.(3) Particle Evaluation. Using (9) as the fitness function, the fitness of each particle can be calculated: Jd(i)=Jθd(i). If Jd(i)<J(pb), the history optimal position of particle pb is updated as pb=θd(i); if Jd(i)<J(gb), then the optimal position of all the particles is updated as gb=θd(i).(4) Particles Update. The position and speed of the particle are updated by the following:(18)vrdi+1=w·vvrdi+c1·rand1·pb-θrdi+c2·rand1·gb-θrdi,θrdi+1=θrdi+vrdi+1,w=wmax-iimaxwmax-wmin.(5) Stop Condition. When the fitness function reaches the setting accuracy or i reaches the maximum number of steps imax, the iteration stops, and the result of θ is gb.The rough solutiongb is within the smaller neighborhood of the exact optimal solution. Using gb as the iterative initial value of Gauss-Newton algorithm, the solution steps are as follows.(1) Make the iteration numberi=0, select the initial value point θ(0)=gb, and calculate the objective function value as shown in the following equation:(19)Jθ0=12∑j=1K·Mej2θ0.(2) Calculategjr(i)=∂ej[θ(i)]/∂θr, j=1,2,…,(K·M), and r=1,2,…,R; then gjr(i) can form the matrix G(i) as shown in the following:(20)Gi=g11ig12i⋯g1Rig21ig22i⋯g2Ri⋮⋮⋱⋮gK·M1igK·M2i⋯gK·MRiK·M×R.The gradient vector ofJ(θ) is b(i)=[b1(i),b2(i),…,bR(i)]T, where(21)bri=∑j=1K·Mejθi∂ejθi∂θr.(3) The linear equationsGT(i)·G(i)p=-b(i) of vector p are solved, where if the rank of GT(i)·G(i) is R, the equation has the unique solution p(i)=-GT(i)·G(i)-1·b(i); if the rank of GT(i)·G(i) is less than R, p(i)=-b(i).(4) Minimize the linear function groupFθ(i)+λ·p(i) of the single variable λ. The direct linear search method is used to obtain λ(k), and λ(k) meet Fθ(i)+λ(i)·p(i)=minλ⁡Fθ(i)+λ·p(i); then the new iteration point and the target function value are as shown in the following:(22)θi+1=θi+λi·pi,Jθi+1=12∑j=1K·Mej2θi+1.(5) If the termination condition is satisfied, the iteration will stop, andθ^=θ(i+1) is the final solution; otherwise make i=i+1, and return b step to continue the iteration.


## 2.2.4. Energy Baseline Model Checking

Through the above processes, the final energy baseline model is followed:(23)Y=fX,θ^.The adjusted decision coefficient is used as the evaluation criteria to check the obtained model.Firstly, the estimated values of EnPI are calculated on each measured point:(24)Y^mk=fX1k,X2k,…,XNk,θ^,k=1,2,…,K,m=1,2,…,M.Make(25)y=Y11,Y21,…,YM1,…,Y1K,Y2K,…,YMK,y^=Y^11,Y^21,…,Y^M1,…,Y^1K,Y^2K,…,Y^MK;the following statistics can be calculated.Total Square Sum is given as follows:(26)TSS=y-y-Ty-y-.Explain Square Sum is given as follows:(27)ESS=y^-y-Ty^-y-.Then the adjusted decision coefficient is as shown in the following equation:(28)R2¯=1-K×M-1K×M-R-11-ESSTSS2,where R2¯ is adjusted decision coefficient. The value range of R2¯ is [0,1]; the more close to 1, the better to express the relationship between the EnPI and the impact factors for the energy baseline model. The threshold can be set in advance according to the need; if the adjusted decision coefficient is greater than the threshold, it is considered that the model fitting is good; otherwise it is considered that the model fitting is poor; the model structure or impact factors need to be reselected.


## 2.3. Energy Efficiency Evaluation Based on Energy Baseline

According to the above sessions, the energy baseline can implement two functions in EnMS:(1)
Energy baseline can be used as energy performance evaluation benchmark after the implementation of energy plan.(2)
Energy baseline can be used as energy performance monitoring and system fault diagnosis on the base of its real-time estimation of energy performance indicators.


### 2.3.1. Statistics Related to EnPIs

Firstly, three benchmark statistics are provided to judge whether energy performance achieves the desired standard:(1)
Expected reduced amplitude of EnPIs:eD=(e1D,e2D,…,eMD)(2)
The mean sample deviation of energy baseline model:e-B=(e-1B,e-2B,…,e-MB)(29)e-mB=1K∑k=1KYmk-Y^mk.(3)
The standard deviation of sample deviation of energy baseline model:(30)σB=σ1B,σ2B,…,σMB,σmB=1K-1∑k=1KYmk-Y^mk-e-mB2,m=1,2,…,M.Then a period of time to be evaluated is defined, and then the measured values of EnPIs and impact factors are(y1p,y2p,…,yMp,x1p,x2p,…,xNp), p=1,2,…,P. The measured values of EnPIs can be expressed as the following form: (31)ym=ym1,ym2,…,ymPT,m=1,2,…,M.Then the impact factors values are substituted into (21) of energy baseline model in order to obtain the estimated values of EnPIs: (32)y^m=y^m1,y^m2,…,y^mPT.Defineem=(em1,em2,…,emP)T=ym-y^m; then the mean value and standard deviation of em are (33)e-m=1P∑p=1Pemp,σm=1P-1∑p=1Pemp2.


### 2.3.2. Energy Performance Evaluation

In order to evaluate the energy performance during a period of time after the implementation of the energy plan, the energy baseline model is established firstly according to the data before the implementation of the energy plan; then evaluated values of EnPIs through the model in this period are obtained. Through the comparison between the estimated values and the measured values, the implementation effect of the energy plan can be estimated quantitatively. Specifically, for eachYm (m=1,2,…,M), e-m is calculated through (33) and compared with emD. The evaluation criteria are as follows:(1)
If(-e-m)≥emD, the energy management system achieves or exceeds expected effect; it meets the requirements of indicators.(2)
If(-e-m)<emD, the actual reduced amplitude of EnPI is less than expected value, and it does not meet the requirements; the enterprise needs to review its energy again and make a new plan.


### 2.3.3. Energy Performance Monitoring and System Fault Diagnosis

Energy efficiency monitoring and system fault diagnosis method based on energy baseline are used in the checking period of energy management system. The energy baseline model reflects the current state of the system operation; whether the system is normal operation can be judged by it. Specifically, for eachYm (m=1,2,…,M), e-m is calculated through (33), and the judging criteria are as follows:(1)
IfemB-3σmD≤e-m≤emB+3σmD, the energy management system runs normally.(2)
Ife-m<emB-3σmD, the energy baseline model is no longer suitable for estimating the current state of EnPI; the model needs to be rebuilt according to the current data.(3)
Ife-m>emB+3σmD, the energy consumption is too high or the system is failure; enterprise needs to make a new plan or system needs maintenance.


## 2.3.1. Statistics Related to EnPIs

Firstly, three benchmark statistics are provided to judge whether energy performance achieves the desired standard:(1)
Expected reduced amplitude of EnPIs:eD=(e1D,e2D,…,eMD)(2)
The mean sample deviation of energy baseline model:e-B=(e-1B,e-2B,…,e-MB)(29)e-mB=1K∑k=1KYmk-Y^mk.(3)
The standard deviation of sample deviation of energy baseline model:(30)σB=σ1B,σ2B,…,σMB,σmB=1K-1∑k=1KYmk-Y^mk-e-mB2,m=1,2,…,M.Then a period of time to be evaluated is defined, and then the measured values of EnPIs and impact factors are(y1p,y2p,…,yMp,x1p,x2p,…,xNp), p=1,2,…,P. The measured values of EnPIs can be expressed as the following form: (31)ym=ym1,ym2,…,ymPT,m=1,2,…,M.Then the impact factors values are substituted into (21) of energy baseline model in order to obtain the estimated values of EnPIs: (32)y^m=y^m1,y^m2,…,y^mPT.Defineem=(em1,em2,…,emP)T=ym-y^m; then the mean value and standard deviation of em are (33)e-m=1P∑p=1Pemp,σm=1P-1∑p=1Pemp2.


## 2.3.2. Energy Performance Evaluation

In order to evaluate the energy performance during a period of time after the implementation of the energy plan, the energy baseline model is established firstly according to the data before the implementation of the energy plan; then evaluated values of EnPIs through the model in this period are obtained. Through the comparison between the estimated values and the measured values, the implementation effect of the energy plan can be estimated quantitatively. Specifically, for eachYm (m=1,2,…,M), e-m is calculated through (33) and compared with emD. The evaluation criteria are as follows:(1)
If(-e-m)≥emD, the energy management system achieves or exceeds expected effect; it meets the requirements of indicators.(2)
If(-e-m)<emD, the actual reduced amplitude of EnPI is less than expected value, and it does not meet the requirements; the enterprise needs to review its energy again and make a new plan.


## 2.3.3. Energy Performance Monitoring and System Fault Diagnosis

Energy efficiency monitoring and system fault diagnosis method based on energy baseline are used in the checking period of energy management system. The energy baseline model reflects the current state of the system operation; whether the system is normal operation can be judged by it. Specifically, for eachYm (m=1,2,…,M), e-m is calculated through (33), and the judging criteria are as follows:(1)
IfemB-3σmD≤e-m≤emB+3σmD, the energy management system runs normally.(2)
Ife-m<emB-3σmD, the energy baseline model is no longer suitable for estimating the current state of EnPI; the model needs to be rebuilt according to the current data.(3)
Ife-m>emB+3σmD, the energy consumption is too high or the system is failure; enterprise needs to make a new plan or system needs maintenance.


## 2.4. PDCA Principle Implementation of Energy Baseline

The EnMS is based on Plan-Do-Check-Act (PDCA) continual improvement principle, and the energy baseline will be adjusted in different phases to implement different functions:(1)
In Plan phase, energy baseline is established firstly based on energy history information and is as the benchmark of the next phase.(2)
In Do phase, the energy baseline is used to evaluate the difference of the energy performance before and after the implementation of the energy plan.(3)
In Check phase, the whole system runs smoothly, the energy baseline needs to be readjusted based on the energy data after the Do phase, and, at this time, the energy baseline can be used to monitor energy efficiency or diagnose system fault.(4)
In Action phase, the energy efficiency needs to be further improved to take actions; then it goes into the Plan phase; the present energy baseline is as the benchmark for next phase, and the EnMS forms a closed loop.


## 3. Results and Discussion

In this section, the cracked gas compression units of the ethylene plant in some petrochemical companies are regarded as the energy management object in order to establish the energy efficiency monitoring and evaluation system according to the ISO 50001 standard and decrease the energy consumption. And the specific indicators are to reduce energy consumption of per unit output more than 5%.


### 3.1. Cracked Gas Compression Unit Process Principle

The main equipment of cracked gas compressor unit includes five stage compressors whose numbers are K-13001, K-13002, K-13003, K-13004, and K-13005; alkali washing tower whose number is C-1340; and high pressure depropanizer tower whose number is C-1365.In order to prevent compressor surge since the entry of the compressor flow is too small, 3 antisurge lines are designed in the cracked gas compressor:(1)
“Three return one”: return a portion of the gas discharged from the third stage compressor to the suction inlet of the first stage compressor to ensure that the first three compressors’ flow is greater than required minimum flow.(2)
“Four return four”: return a portion of the gas discharged from the fourth stage compressor to the suction inlet of the fourth stage compressor to ensure that the fourth compressor’s flow is greater than required minimum flow.(3)
“Five return five”: return a portion of the gas discharged from the fifth stage compressor to the high pressure propane tower; then the gas of the top of the high pressure propane tower returns into the fifth stage compressor.


### 3.2. Energy Baseline of Cracked Gas Compression Unit

To implement advanced process control (APC) system, it is planned to reduce the return flow of the 4 surge lines to improve the output of unit energy consumption and achieve the purpose of energy-saving [15–19].


#### 3.2.1. EnPI and Impact Factors

Energy input for cracked gas compressor unit is ultrahigh pressure steam; product output for cracked gas compressor unit is compressed cracked gas, so EnPI is defined as(34)EnPI=turbineinletflowK-13005outletflow.And the impact factors areX1 to X15 shown in Table 2.Table 2
Energy performance related variables and parts of measured data.
Name Number Unit Data (per minute) 10:00 10:01 10:02 10:03 10:04 10:05 ⋯ The compressor speed X 1 RPM 4078 4067 4025 4070 4059 4012 ⋯ K-13001 inlet pressure X 2 KPA 33.3 32.8 32.9 33.2 33.0 33.2 ⋯ K-13001 outlet temperature X 3 °C 35.0 34.9 35.4 35.4 35.3 35.2 ⋯ K-13004 outlet pressure X 4 °C 75.2 74.8 73.9 75.2 74.9 74.7 ⋯ K-13004 outlet pressure X 5 MPA 1.555 1.549 1.540 1.547 1.548 1.549 ⋯ K-13004 outlet flow X 6 T/H 241.9 241.9 242.7 245.2 243.2 241.4 ⋯ K-13005 outlet temperature X 7 °C 70.9 71.8 70.8 71.3 70.8 71.2 ⋯ K-13005 outlet pressure X 8 MPA 3.593 3.634 3.660 3.654 3.639 3.638 ⋯ Turbine inlet pressure X 9 MPA 11.48 11.59 11.58 11.54 11.65 11.56 ⋯ Turbine inlet temperature X 10 °C 518.2 515.9 515.5 514.8 513.5 515.0 ⋯ C-1365 tower top temperature X 11 °C −21.5 −21.9 −21.6 −21.5 −21.7 −21.8 ⋯ C-1365 tower top pressure X 12 MPA 1.461 1.440 1.451 1.458 1.463 1.442 ⋯ K-1300 three return one X 13 % 7.1 7.2 7.1 7.2 7.1 7.2 ⋯ K-1300 four return four X 14 % 6.5 6.5 6.4 6.4 6.5 6.5 ⋯ K-1300 five return five X 15 % 6.6 6.6 6.6 6.6 6.6 6.5 ⋯ K-13005 outlet flow X 16 T/H 223.7 223.8 223.1 222.8 226.5 226.5 ⋯ Turbine outlet flow X 17 T/H 323.1 322.8 320.0 324.0 323.5 322.3 ⋯ EnPI (X17/X16) Y — 1.444 1.442 1.434 1.454 1.428 1.423 ⋯


#### 3.2.2. Energy Baseline Model

Since the cracked gas compression unit is working near the set point, it presents linear relationship between EnPI and impact factors, so the energy baseline model is defined as a linear model as follows:(35)Y=θ0+∑i=115θiXi,where θ0 and θi are model parameters, Y is an EnPI, Xi is an impact factor, and i=1,2,…,15. Variables defined in detail and parts of the measured data are shown in Table 2.The model parameters are solved according to the linear model solution method; the results are as follows:(36)θ^=0.52440.00030.0310-0.02220.0045-0.23590.0016-0.0218-0.01910.1481-0.00010.0289-0.74380.1478-0.13230.0545T,R2¯=1-120-1120-16-11-0.2720.2832=0.904.R 2 ¯ is 0.904, which is close to 1 (generally, for the time series, it is considered to have a higher fitting accuracy that R2¯ achieves 0.9.), so the model has higher fitting accuracy.


### 3.3. Energy Performance Evaluation

Through the application of the advanced process control on the cracking gas compressor unit, the return flow ratio of each surge line is reduced. The compression unit is working stably on the point close to the surge line, and energy consumption of per unit of output is reduced. The effect analysis of energy-saving is as follows.The measured data is obtained in two hours after implementation of APC. And the comparison result between the actual measured value and estimated values of EnPI after implementing APC is shown in Figure3.Figure 3
EnPI contrast curve before and after implementation of APC.The actual measured values of EnPI are expressed asy=(y1,y2,…,y120)T, the estimated values of EnPI calculated by the energy baseline model are expressed as y^=(y^1,y^2,…,y^120)T, e=y-y^ is defined, and statistics are calculated as follows: (37)y^-=1120∑p=1120yp=1.3686,e-=1120∑p=1120ep=-0.1213.So the average of reduced amplitude of EnPI is(38)EnPI  reduced  amplitude  ratio=-e-y^-·100=8.86%.Therefore, according to the evaluation criteria of energy efficiency, the establishment of the energy management system and the used technical method complete the expected indicator to save energy by 5% for the cracked gas compression unit.


## 3.1. Cracked Gas Compression Unit Process Principle

The main equipment of cracked gas compressor unit includes five stage compressors whose numbers are K-13001, K-13002, K-13003, K-13004, and K-13005; alkali washing tower whose number is C-1340; and high pressure depropanizer tower whose number is C-1365.In order to prevent compressor surge since the entry of the compressor flow is too small, 3 antisurge lines are designed in the cracked gas compressor:(1)
“Three return one”: return a portion of the gas discharged from the third stage compressor to the suction inlet of the first stage compressor to ensure that the first three compressors’ flow is greater than required minimum flow.(2)
“Four return four”: return a portion of the gas discharged from the fourth stage compressor to the suction inlet of the fourth stage compressor to ensure that the fourth compressor’s flow is greater than required minimum flow.(3)
“Five return five”: return a portion of the gas discharged from the fifth stage compressor to the high pressure propane tower; then the gas of the top of the high pressure propane tower returns into the fifth stage compressor.


## 3.2. Energy Baseline of Cracked Gas Compression Unit

To implement advanced process control (APC) system, it is planned to reduce the return flow of the 4 surge lines to improve the output of unit energy consumption and achieve the purpose of energy-saving [15–19].


### 3.2.1. EnPI and Impact Factors

Energy input for cracked gas compressor unit is ultrahigh pressure steam; product output for cracked gas compressor unit is compressed cracked gas, so EnPI is defined as(34)EnPI=turbineinletflowK-13005outletflow.And the impact factors areX1 to X15 shown in Table 2.Table 2
Energy performance related variables and parts of measured data.
Name Number Unit Data (per minute) 10:00 10:01 10:02 10:03 10:04 10:05 ⋯ The compressor speed X 1 RPM 4078 4067 4025 4070 4059 4012 ⋯ K-13001 inlet pressure X 2 KPA 33.3 32.8 32.9 33.2 33.0 33.2 ⋯ K-13001 outlet temperature X 3 °C 35.0 34.9 35.4 35.4 35.3 35.2 ⋯ K-13004 outlet pressure X 4 °C 75.2 74.8 73.9 75.2 74.9 74.7 ⋯ K-13004 outlet pressure X 5 MPA 1.555 1.549 1.540 1.547 1.548 1.549 ⋯ K-13004 outlet flow X 6 T/H 241.9 241.9 242.7 245.2 243.2 241.4 ⋯ K-13005 outlet temperature X 7 °C 70.9 71.8 70.8 71.3 70.8 71.2 ⋯ K-13005 outlet pressure X 8 MPA 3.593 3.634 3.660 3.654 3.639 3.638 ⋯ Turbine inlet pressure X 9 MPA 11.48 11.59 11.58 11.54 11.65 11.56 ⋯ Turbine inlet temperature X 10 °C 518.2 515.9 515.5 514.8 513.5 515.0 ⋯ C-1365 tower top temperature X 11 °C −21.5 −21.9 −21.6 −21.5 −21.7 −21.8 ⋯ C-1365 tower top pressure X 12 MPA 1.461 1.440 1.451 1.458 1.463 1.442 ⋯ K-1300 three return one X 13 % 7.1 7.2 7.1 7.2 7.1 7.2 ⋯ K-1300 four return four X 14 % 6.5 6.5 6.4 6.4 6.5 6.5 ⋯ K-1300 five return five X 15 % 6.6 6.6 6.6 6.6 6.6 6.5 ⋯ K-13005 outlet flow X 16 T/H 223.7 223.8 223.1 222.8 226.5 226.5 ⋯ Turbine outlet flow X 17 T/H 323.1 322.8 320.0 324.0 323.5 322.3 ⋯ EnPI (X17/X16) Y — 1.444 1.442 1.434 1.454 1.428 1.423 ⋯


### 3.2.2. Energy Baseline Model

Since the cracked gas compression unit is working near the set point, it presents linear relationship between EnPI and impact factors, so the energy baseline model is defined as a linear model as follows:(35)Y=θ0+∑i=115θiXi,where θ0 and θi are model parameters, Y is an EnPI, Xi is an impact factor, and i=1,2,…,15. Variables defined in detail and parts of the measured data are shown in Table 2.The model parameters are solved according to the linear model solution method; the results are as follows:(36)θ^=0.52440.00030.0310-0.02220.0045-0.23590.0016-0.0218-0.01910.1481-0.00010.0289-0.74380.1478-0.13230.0545T,R2¯=1-120-1120-16-11-0.2720.2832=0.904.R 2 ¯ is 0.904, which is close to 1 (generally, for the time series, it is considered to have a higher fitting accuracy that R2¯ achieves 0.9.), so the model has higher fitting accuracy.


## 3.2.1. EnPI and Impact Factors

Energy input for cracked gas compressor unit is ultrahigh pressure steam; product output for cracked gas compressor unit is compressed cracked gas, so EnPI is defined as(34)EnPI=turbineinletflowK-13005outletflow.And the impact factors areX1 to X15 shown in Table 2.Table 2
Energy performance related variables and parts of measured data.
Name Number Unit Data (per minute) 10:00 10:01 10:02 10:03 10:04 10:05 ⋯ The compressor speed X 1 RPM 4078 4067 4025 4070 4059 4012 ⋯ K-13001 inlet pressure X 2 KPA 33.3 32.8 32.9 33.2 33.0 33.2 ⋯ K-13001 outlet temperature X 3 °C 35.0 34.9 35.4 35.4 35.3 35.2 ⋯ K-13004 outlet pressure X 4 °C 75.2 74.8 73.9 75.2 74.9 74.7 ⋯ K-13004 outlet pressure X 5 MPA 1.555 1.549 1.540 1.547 1.548 1.549 ⋯ K-13004 outlet flow X 6 T/H 241.9 241.9 242.7 245.2 243.2 241.4 ⋯ K-13005 outlet temperature X 7 °C 70.9 71.8 70.8 71.3 70.8 71.2 ⋯ K-13005 outlet pressure X 8 MPA 3.593 3.634 3.660 3.654 3.639 3.638 ⋯ Turbine inlet pressure X 9 MPA 11.48 11.59 11.58 11.54 11.65 11.56 ⋯ Turbine inlet temperature X 10 °C 518.2 515.9 515.5 514.8 513.5 515.0 ⋯ C-1365 tower top temperature X 11 °C −21.5 −21.9 −21.6 −21.5 −21.7 −21.8 ⋯ C-1365 tower top pressure X 12 MPA 1.461 1.440 1.451 1.458 1.463 1.442 ⋯ K-1300 three return one X 13 % 7.1 7.2 7.1 7.2 7.1 7.2 ⋯ K-1300 four return four X 14 % 6.5 6.5 6.4 6.4 6.5 6.5 ⋯ K-1300 five return five X 15 % 6.6 6.6 6.6 6.6 6.6 6.5 ⋯ K-13005 outlet flow X 16 T/H 223.7 223.8 223.1 222.8 226.5 226.5 ⋯ Turbine outlet flow X 17 T/H 323.1 322.8 320.0 324.0 323.5 322.3 ⋯ EnPI (X17/X16) Y — 1.444 1.442 1.434 1.454 1.428 1.423 ⋯


## 3.2.2. Energy Baseline Model

Since the cracked gas compression unit is working near the set point, it presents linear relationship between EnPI and impact factors, so the energy baseline model is defined as a linear model as follows:(35)Y=θ0+∑i=115θiXi,where θ0 and θi are model parameters, Y is an EnPI, Xi is an impact factor, and i=1,2,…,15. Variables defined in detail and parts of the measured data are shown in Table 2.The model parameters are solved according to the linear model solution method; the results are as follows:(36)θ^=0.52440.00030.0310-0.02220.0045-0.23590.0016-0.0218-0.01910.1481-0.00010.0289-0.74380.1478-0.13230.0545T,R2¯=1-120-1120-16-11-0.2720.2832=0.904.R 2 ¯ is 0.904, which is close to 1 (generally, for the time series, it is considered to have a higher fitting accuracy that R2¯ achieves 0.9.), so the model has higher fitting accuracy.


## 3.3. Energy Performance Evaluation

Through the application of the advanced process control on the cracking gas compressor unit, the return flow ratio of each surge line is reduced. The compression unit is working stably on the point close to the surge line, and energy consumption of per unit of output is reduced. The effect analysis of energy-saving is as follows.The measured data is obtained in two hours after implementation of APC. And the comparison result between the actual measured value and estimated values of EnPI after implementing APC is shown in Figure3.Figure 3
EnPI contrast curve before and after implementation of APC.The actual measured values of EnPI are expressed asy=(y1,y2,…,y120)T, the estimated values of EnPI calculated by the energy baseline model are expressed as y^=(y^1,y^2,…,y^120)T, e=y-y^ is defined, and statistics are calculated as follows: (37)y^-=1120∑p=1120yp=1.3686,e-=1120∑p=1120ep=-0.1213.So the average of reduced amplitude of EnPI is(38)EnPI  reduced  amplitude  ratio=-e-y^-·100=8.86%.Therefore, according to the evaluation criteria of energy efficiency, the establishment of the energy management system and the used technical method complete the expected indicator to save energy by 5% for the cracked gas compression unit.


## 4. Conclusions

This study analyzes the energy structure of the current world and the problems faced and expounds the necessity of energy management in chemical enterprises. Based on ISO 50001 standard, the structure of energy management system is established, and the energy performance evaluation method based on energy baseline is proposed. The method firstly establishes the baseline models for different kinds of systems and gives the general solution method for each kind of models and then illustrates how to use the energy baseline model for energy performance evaluation and system fault diagnosis. Finally, the application of this method in the cracked gas compression unit of certain petrochemical enterprise makes the object that energy performance is decreased by 5% come true, and the production costs are greatly reduced at the same time. It can be concluded that this method is correct and practical.



---

*Source: 7087393-2016-11-02.xml*