Novel Robust Least-Squares Estimator for Linear ... - ACS Publications

Subscriber access provided by Hong Kong University of Science and Technology Library

General Research

A Novel Robust Least Squares Estimator for Linear Dynamic Data Reconciliation Rui Zhang, Yi Zhang, Yiping Feng, and Gang Rong Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b04853 • Publication Date (Web): 31 Dec 2018 Downloaded from http://pubs.acs.org on January 2, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

118x70mm (144 x 144 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. A simple process network 20x5mm (300 x 300 DPI)


Page 2 of 57

Page 3 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Flowchart of Robust Least Squares Dynamic Data Reconciliation 181x265mm (300 x 300 DPI)



Figure 3. Flowchart of the new robust least squares data reconciliation 262x550mm (300 x 300 DPI)


Page 4 of 57

Page 5 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Reconciliation results comparison of 313x209mm (72 x 72 DPI)


(no gross error)


Figure 5. Reconciliation results comparison of ( transient gross error) 313x204mm (72 x 72 DPI)


Page 6 of 57

Page 7 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure.6. Reconciliation results comparison of (continuous gross errors) 313x201mm (72 x 72 DPI)



Figure.7. Reconciliation results comparison of (continuous gross errors) 313x201mm (72 x 72 DPI)


Page 8 of 57

Page 9 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure.8. Reconciled bias of

by five methods (no gross errors)

312x149mm (72 x 72 DPI)



Figure.9. Reconciled bias of

by five methods (discontinuous gross errors)

379x206mm (72 x 72 DPI)


Page 10 of 57

Page 11 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure.10a. Reconciled bias of

(transient gross errors)

312x149mm (72 x 72 DPI)



Figure.10b. Reconciled bias of

(gross errors last for four time cells)

312x149mm (72 x 72 DPI)


Page 12 of 57

Page 13 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure.10c. Reconciled bias of

(gross errors last for six time cells)

312x149mm (72 x 72 DPI)



Figure 11. Network representation of steam system for a methanol synthesis unit 60x47mm (300 x 300 DPI)


Page 14 of 57

Page 15 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

A Novel Robust Least Squares Estimator for

2

Linear Dynamic Data Reconciliation

3

Rui ZHANG, Yi ZHANG, Yiping FENG, Gang RONG*

4

State Key Laboratory of Industrial Control Technology, College of Control Science

5

and Engineering, Zhejiang University, Hangzhou, 310027, China

6

7

ABSTRACT

8

Dynamic data reconciliation is used to reduce the uncertainties in process

9

measurement. Conventional data reconciliation theories and methods are based on least

10

squares estimation, whose conditions are hard to meet in real-world applications. Since

11

least squares estimators can be made robust by equivalent weight, many works

12

concentrated on robust estimators and their performance in data reconciliation.

13

However, robust least squares algorithm is still disable of detecting gross errors in low

14

redundancy variables. In this paper, our work takes a step further on improving the

15

robustness of Huber estimator by taking structural redundancy into consideration.

16

Modifications on influence function and weight function are made by involving local

17

redundancy and rejection interval. The detectivity of gross error in low redundant

18

variables is increased. Combined with above improvement and then taking the

1



1

advantage of online filtering, a novel robust dynamic data reconciliation algorithm is

2

proposed. Case studies demonstrate the superiority of this methodology.

3

KEYWORDS: Huber estimation; Robust least squares; Local redundancy; Rejection

4

interval;

5 6

1. INTRODUCTION

7

Linear dynamic data reconciliation (LDDR) is a way of estimating variables on

8

the basis of measurements and dynamic equilibrium equation. Dynamic material

9

balance model can be represented by continuous-state space equations or a sampled

10

input-output representation after discretization. Many other studies of LDDR have been

11

carried out, but most of them did not properly deal with the complexity of dynamic

12

systems and the demand for timeliness. Almasy(1990)1 transformed the dynamic

13

system balance problem into Kalman filtering problem by discretizing input and output

14

variables. Darouach and Zasadzinski(1991)2 transformed the LDDR into an algebraic

15

problem by using the forward difference approximation method, and obtained an

16

iterative algorithm to avoid the truncation error and matrix singularity problem. Integral

17

approach proposed by Bagajewicz and Jiang(1997)3 provides another choice for solving

18

the problem. But it is difficult to separate a signal from noise. Bagajewicz(2000)4 also

19

applies this method to steady state data reconciliation(SSDR) and DDR, proving that

20

for linear systems without tank flow, the performance of both approaches is similar in

21

the absence of biases and leaks.

22

Almasy(1990)1 studied DDR based on Kalman filter and proved that Kalman filter

23

equations are convenient for solving LDDR. Although Kalman filter has great 2


Page 16 of 57

Page 17 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

characteristics in theory, it has to suppose that the state variables are self-correlated. To

2

solve the general LDDR problem, the least squares principals of linear SSDR for

3

dynamic systems are extended. The new solution is based on the concept of generalized

4

linear dynamic model, which was first introduced to DR by Darouach and

5

Zasadzinski(1991)2. Generalized linear dynamic model refers to the model formulation

6

Dxk +1 = Bxk , which contains more state variables than constrains. Because D is

7

singular, the formulation cannot be written in a standard state space form. Obviously

8

standard Kalman Filter is not suitable for this situation, but we can utilize the recursive

9

form of Kalman filter, which is more convenient for real-time processing. Bai(2006)

10

5

11

many studies have been working on a new recursive filtering algorithm. However, the

12

model is quite huge and the derivation procedure is complex and tedious6. Rollins and

13

Devanathan(1993)7 adopted the same model formulation and proposed a constrained

14

least squares estimator,

15

instants. Xu(2010)8 proposed a new framework for generalized dynamic system, whose

16

core part is a simplified least squares formulation to represent the filtering problem.

pointed that DDR is actually a kind of filter technique. Based on the new model form,

making use of

measurements at continuous adjacent time

17

The least squares estimator(LSE) is the most widely-used parameter estimation

18

method. But the lack of robustness makes it not perfect in practice. Johnston and

19

Kramer(1995)9 reported the feasibility and better performance of robust estimators as

20

the objective function in data reconciliation especially when measurements are

21

contaminated by gross errors. In the literature, DDR problems were addressed using

22

Fair function(FF)10, three-part redescending estimator(TPRE)11 , Correntropy(CO)

23

function12, and so on. Subsequently, many studies have conducted on robust estimators

24

and their performance in data reconciliation. Including Arora and Biegler(2001)13, 3



1

many studies have demonstrated the potential of robust statistics developed by

2

Huber(1981)14, which attempts accurate estimation of statistical parameters. However,

3

according to Özyurt(2004)15, the influence function for weighted least square(WLS) is

4

proportional to measurement error. The occurrence

5

incorrect estimations, it is obviously not a robust estimator. But the research on Quasi-

6

weighted least squares estimator(QWLS)16 improved the robustness of WLS, which

7

retained the advantage of WLS, furthermore, enhanced the practicality in the presence

8

of gross errors. A unified view on robust data reconciliation is provided using

9

generalized objective functions, within a probabilistic framework17. M-estimators are a

10

maximum likelihood type estimator, which is one of the three main classes of robust

11

estimators. Based on the principle of equivalent weights, M-estimator can be

12

transformed into the form of classical least squares. In other words, classical LSE can

13

be made more robust by using equivalent weights18. This method not only improves the

14

robustness of classical estimators, but also preserves the advantages of simplicity in

15

mathematical forms and computation procedure. Particularly, many well-developed

16

mathematical models and computing methods based on least squares are still applicable.

17

The famous Huber estimator is the maximum likelihood estimation of Huber function19.

18

The Huber estimator is a kind of minimax estimation, which laid the foundation for

19

location robust estimation. Huber20 had made his first step towards attempting accurate

20

estimation of statistical parameters in the presence of gross errors. Since that, robust

21

estimators have been widely used in mathematical statistics, communication

22

engineering, exploration mapping and so on. From this point of view, data processing

23

derived from Huber estimation in this paper is essentially a new robust least

24

squares(RLS) methodology, which is named general robust least squares(GRLS). 4


of gross errors will lead to

Page 18 of 57

Page 19 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

After Huber, many researchers had paid much attention to improve other M-

2

estimators. Ni(2009)21 proposed two computational approaches for regression with

3

stochastically bounded noise (RSBN). Hampel(2011)22 focused on a smoothing

4

principle for the Huber and other location M-estimators. Directed towards the

5

correction of endogeneity problems in linear models, Kim(2007)23 proposed a new

6

robust estimator in the context of two-stage estimation methods. A new target function

7

was taken into consideration based on robust M-estimators by Jin(2012)24.

8

Zhang(2010)16 has proposed a quasi-weighted least squares(QWLS) estimator to

9

promote the performance of LSE. Claudia(2017)25 distinguished systemic errors into

10

outliers, biases, and drifts, and also studied the performance of different M-estimators.

11

However, few works have been carried out on the improvement of Huber estimator

12

strategies for linear processes. The effectiveness of conventional algorithm depends on

13

the reliability of approximated estimation solutions and reasonability of weight function

14

during iterations. Conventional weight function applies variances with prior probability

15

distributing measurement, disregarding the influence of robustness in structure space.

16

Zhou had pointed out that the Huber estimator has no limitations on location of spatial

17

distribution of observations18. In other words, gross errors may hardly be shown by

18

residuals, but the introduction of local redundancy would raise the detectivity of gross

19

errors especially for low redundancy variables.

20

The novel general robust least squares estimator(GRLSE) is proposed for linear

21

DDR, which has good performance in low redundancy system and detecting gross

22

errors. This algorithm focuses on making Huber estimator more robust based on the in-

23

depth discussion on advantages and disadvantages of classical LSE. In this paper, we

24

improved the most important function of Huber estimation, the influence function, and 5



Page 20 of 57

1

weight function by considering construction redundancy. The derivation process and

2

algorithm procedure are combined with online filtering for performing online

3

optimization easily. Due to the lack of accurate instruments and regular overhaul, low-

4

redundant variables commonly exist, which has been a big challenge for DDR. The

5

algorithm has good performance, especially for low-redundant variables. To avoid

6

miscalculation introduced from the newly used local redundancy, the procedure has

7

also been improved by combining constraints extracted from real process. The process

8

has been tested in several cases and dozens of situations, not only listed in the paper.

9

The industrial case presented here is persuasive.

10

2. PROBLEM STATEMENT

11

In this section, we focus on improving Huber estimator. The main part of Huber

12

distribution obeys normal distribution N (0, σ 2 ) , and the disturbance part obeys Laplace

13

distribution18.

14

function. Thus, the loss function of Huber estimator is

15

16

17

18

Huber estimator is the maximum likelihood estimation of Huber

 vi2 | vi |≤ cσ i  2  2σ i ρ (vi ) =  2 c | vi | − c | v |> cσ i i 2  σ i

(1)

The influence function of Huber estimator is −c vi / σ i < −c  ϕ (v= vi − c ≤ vi / σ i ≤ c i) c v / σ > c i i 

The weight function of Huber estimator is

6


(2)

Page 21 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


| vi / σ i |≤ c 1 w(vi ) =  cσ i sign(vi ) / vi | vi / σ i |> c

1

(3)

2

vi stands for residuals of the ith variable;

3

xi stands for measurements of the ith variable;

4

xˆi stands for reconciled values of the ith variable and vi =xi − xˆi ;

5

σ i stands for standard deviation of the ith variable;

6

c is a constant for Huber estimator, which is determined by trade-off between

7

robustness and estimation efficiency. Estimation efficiency is the ratio of Cramer

8

boundary and asymptotic variance. For Huber estimator, the efficiency is

9

= eff

1 [2Φ (c) − 1]2 = Var 2Φ (c) − 1 − 2 c φ (c) + 2 c 2 (1 − Φ (c))

( 4)

10

Φ (c) and φ (c) stand for distribution function and probability density value of

11

standard normal distribution at point c . When c = 1.4 , eff = 95.5% . In actual

12

calculation, it should be determined by the rate of gross error δ :

1 φ (c ) = 2Φ (c) − 1 + 2 1− δ c

13

14 15

( 5)

δ is always between 1% and 10%, so c is generally between 1 and 2. Mostly we choose c as 1.5 around.

16

The effect of Huber estimator is almost the same as LSE for tiny residuals. Thus, the

17

estimator has high efficiency derived from Gaussian distribution. When the residuals

18

are large, the performance of Huber estimator is approximately the same as least

19

absolute deviations estimator. If gross errors appear, bounded weight will be assigned

20

to the measurement biases to narrow the diffusion of gross error. However, for low7



1

redundant observation, Huber estimator could not reflect its gross error by residual,

2

which may lead to miscalculation. To preserve the advantage of Huber estimation, we

3

propose a novel RLS estimator based on the improvement of the Huber one.

4

Firstly, the influence function in Eq.(2) will be modified. Not only the priori

5

measurement variance was replaced by measurement residual variance, but also the

6

added local redundancy would be used to improve robustness in design space. The

7

modified function will increase the detection rate of gross error, especially for low

8

redundancy variables. Secondly, the weight function Eq.(3) of Huber estimator will be

9

reformulated to get rid of the contaminated measurements. Lastly, to weaken the

10

disadvantage of introducing local redundancy, a concrete iterative flow of the novel

11

robust LSE will be described in detail in section 4. The proposed improvements will be

12

elaborated and proved to be effective in the rest of this paper.

13

3. GENERAL FORMULATION

14

Robust least squares estimator (RLSE) combines the robust estimation principle with

15

the least squares method. In the robust least squares estimation, a reasonable unbiased

16

estimation of ρ function will be constructed to reflect the adopted scheme for

17

abnormal measurements. This paper uses the Huber function, which is not very

18

sensitive to measurement biases. Some improvements based on Huber function lead to

19

better robustness. The general formulation of modeling and reconciliation is described

20

in this section.

21

3.1 Local Redundancy

22

The intuitive concepts of observability and redundancy are explained in the literature.

23

The concepts were firstly introduced by Kretsovalis and Mah (1988)26. Carpentier et

24

al.(1991)27 split the variables into three classes according to redundancy, which is used 8


Page 22 of 57

Page 23 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

to formulate a well-conditioned least square problem and to compute the values of

2

redundant and calculable variables. Modron(1992)28 offered a more complex

3

classification of quantities. Maronna and Arcas(2009)29 showed that the linear

4

reconciliation problem can be represented by a standard multiple linear regression

5

model and the regression approach suggests a natural measure of the redundancy of an

6

observation. And Carpentier et al.(1991)27 first defined a measure of the local

7

redundancy, who called it error detectability of the i-th measurement. Here we call it

8

local redundancy in the following content. According to the footsteps, the concepts and

9

methods used in this article are presented below.

10

The general formulation of data reconciliation is always described as Eq.(6). Xˆ and

11

Z stand for reconciled variable vector and measured variable vector respectively.

12

Q = diag (σ 1 ,..., σ n ) denotes the covariance matrix. min ( Xˆ − Z )T Q −1 ( Xˆ − Z ) s.t. AXˆ = 0, A ∈ R m×n , X ∈ R n

13

14

The measurement is generally described as follows.

15

1,..., n zi =xi + ε i , i =

16 17

(6)

(7)

where random noise ε i ∼ N (0, σ i ) obeys normal distribution. The analytical solutions are listed below.

18

Xˆ= ( I − QAT ( AQAT ) −1 A) Z

(8)

19

Σ xˆ = cov( xˆ ) = ( I − QAT ( AQAT ) −1 A)Q

(9)

20 21

Given the following definition of sensitivity and redundancy respectively, H is redundancy. 9



Page 24 of 57

1

M = I − QAT ( AQAT ) −1 A

(10)

2

H = QAT ( AQAT ) −1 A

(11)

3

We can get Eq.(12) and Eq.(13).

4

Xˆ= MZ= ( I − H ) Z

(12)

5

Z − Xˆ = HZ

(13)

6

Here, M reflects how sensitive the results are to the measurements. When Q is

7

diagonal, the larger H ii is, the smaller M ii is. From Eq.(10), we know that the

8

smaller the covariance matrix is, the more precise the estimator is.

9

When Q is diagonal, M and H have the following characteristics:

10

M2 M = ,H2 H . M and H are idempotent matrices,=

11

0 ≤ H ii ≤ 1, 0 ≤ M ii ≤ 1 (i = 1, 2,...m) , where H ii is the diagonal element of H , M ii

12 13

14

is the diagonal element of M . = RT strace = (H )

trace( M )=

n

= H ii m ∑

n

∑M = i =1

(14)

i =1

ii

n−m

(15)

15

Here, RT is defined as system redundancy by Bagajewicz(1999)30, which is also

16

called overall redundancy by Wang(2001)31. H ii is local redundancy which was

17

firstly presented as comprehensive redundancy by Wang(2001)31. However, they have

18

been ignoring the influence of the precision of sensors on the system redundancy.

10


Page 25 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


S6 S1

U1 U1

S2

U2 U2

S3

U3 U3

S4

U4 U4

S5

S8 U5 U5 S9

S7

1 2

Figure 1. A simple process network

3

In the process of Fig.1, the precision of flowmeter of Si (i = 1...9) is defined as λi .

4

Suppose that λ3 is much higher than others in the measurement system, its

5

measurement could also be calculated by (λ2 , λ7 ) or (λ4 , λ6 ) . However, it is known

6

that the precision of the calculated measurement must be lower than any component, so

7

λ3 would be much lower than its real value, which cannot play the role of

8

reconciliation. Thus, the real precision of sensors should be concerned into redundancy

9

calculation. From Eq.(11), H ii involves both measurement precision Q and

10

network structure A , then the system redundancy and measurement precision both have

11

influence on comprehensive redundancy. H ii is influential in the adjacent variables,

12

which is also named as local redundancy.

13

It is noteworthy that the closer H ii is to 0, the smaller the local redundancy of the

14

variable is, and the greater its impact on the reconciliation accuracy is. While H ii

15

equals to 0, X i is not a redundant variable and cannot be reconciled.

16

3.2 Linear Dynamic Process Modeling

17

18

The general formulation of linear dynamic system is as follows  dW = AF   dt   CF = 0

11


(16)


1 2 3

W and F denote the measurements of inventory and flow; A and C denote

correlation coefficient matrix of capacity nodes and non-capacity nodes. The measurement vector is as follows:

 wi  zi =    fi 

4

5

Page 26 of 57

(17)

In the form of Kalman filter, Z denotes measurement vector for all variables

 z1  z  W  2 = Z = H  −v   F     zn 

6

(18)

7

Bagajewicz(1997)3 pointed out specific variable classification methods for linear

8

DDR. According to the projection matrix transformation proposed by Crowe(1983)32,

9

unmeasured variables in Eq.(16) can be eliminated. Then the general formulation of

10

the system is transformed into

11

 dWR = AR FR  BR dt    CR FR = 0

(19)

12

W  = Z R I  R  − vR  FR 

(20)

13

where all the matrices labeled with index R, represent the obtained measured variables

14

after matrix transformation.

15

After a series of transformations, the above formulations are transformed into the

16

following difference equations and the detailed process is described in the Supporting

17

Information.

18

EX (i + 1) = BX (i ) 12


(21)

Page 27 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Z R (i + = 1) DX (i + 1) − vR (v + 1)

1

2

(22)

Where the coefficient matrix E, B, D are defined as follows.

3

E= [ BR  − Ts ( AR 2 − AR1CR−11C2 )]

(23)

4

B = [ BR  0]

(24)

= D [ I1  I 22 − I 21CR−11CR 2 ]

5

(25)

C R = [CR1  CR 2 ] , remark that CR1 is square matrix and | CR1 |≠ 0 ;

6

where

7

Correspondingly AR = [ AR1  AR 2 ] ; Ts is sampling period.

8

3.3 Robust Data Reconciliation

9

Xu and Rong(2010)8 proposed a simplified least squares formulation to represent the

10

filtering problem in generalized linear dynamic systems. In our work, the least square

11

estimator was replaced by the improved Huber estimator and an iterative solution

12

procedure is designed. To initialize the problem, it is considered that only random errors

13

are present.

14

The generalized linear dynamic process model can be transformed into Eq.(21) and

15

Eq.(22) by deleting unmeasured variables, eliminating algebraic equations, and

16

discretization. Following the former section, the dynamic models will be combined with

17

robust least squares and filtering ideas.

18 19 20

Suppose the reconciliation results of time i were known as Xˆ (i | i ) , and the covariance matrix is Σ(i | i ) , introducing auxiliary variables Y , Y and ζ ,

= Y BX = (i ) EX (i + 1) 13


(26)


Page 28 of 57

1

Y = BXˆ (i | i )

(27)

2

ζ= Y − Y

(28)

3

The covariance matrix of ζ is BΣ(i | i ) B T .

4

Then, Eq.(21) and Eq.(22) can be transformed into Eq.(29).

 Y   E   ζ  =     X (i + 1) −   vR (i + 1)   Z R (i + 1)   D 

5

6 7

In this form, dynamic data reconciliation model could be reconstructed based on robust least squares theory18. min J =

8

(29)

n

∑Q j =1

−1 jj

1 2

ρ (vR , j (i + 1)) + ζ T [ BΣ(i | i ) BT ]−1ζ

 BXˆ (i | i )   E   ζ  s.t.  =    X (i + 1) −   vR (i + 1)   Z R (i + 1)   D 

(30)

9

where Z R (i + 1) is measurement variable, vR , j (i + 1) is the measurement bias of the

10

j th variable on ( i + 1 )th moment, Q jj is the j th diagonal element of covariance

11

12

13 14

matrix, ρ (vR , j (i + 1)) is the loss function. The adopted Huber function is as follows [vR , j (i + 1)]2 / 2, | vR , j (i + 1) |≤ c ⋅ Q jj  ρ (vR , j (i + 1)) =  1 2 c⋅ | vR , j (i + 1) | − c , | vR , j (i + 1) |> c ⋅ Q jj  2

The solution of problem (30) is obtained as follows. The derivative of Eq.(30) with respect to X (i + 1) is T

15

(31)

 E   BΣ(i | i ) BT ]−1    0 D 

0 ζ  =0   P  v R (i + 1) 

14


(32)

Page 29 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

 p1 0 P =  0   0 0

where

 1  Q  jj pj =   c ⋅ sign[vR , j (i + 1)]  Q jj ⋅ vR , j (i + 1) 

2

3

0 0  pn 

(33)

| vR , j (i + 1) |≤ c Q jj

(34) | vR , j (i + 1) |> c Q jj

The solution of Eq.(32) is

4

| i ) [ E T ( BΣ(i − 1| i − 1) BT ) −1 E + DT Q −1 D]−1 Xˆ (i = × [ E T ( BΣ(i − 1| i − 1) BT ) −1 BXˆ (i − 1| i − 1) + DT Q −1Z R (i )]

5

To be similar with filtering formulations, the covariance matrix of Xˆ (i + 1| i + 1)

6

should be iteratively calculated. Eq.(30) is essentially a robust least squares problem.

7

According to the covariance matrix formulation of RLSE, the covariance matrix of

8

Xˆ (i + 1| i + 1) is as follows.

(35)

Σ(i + 1| i + 1) = [ BΣ(i | i ) BT ]−1 0  [ BΣ(i | i ) BT ]−1 0   E  −1 M [E  D ]  T   M 0 0 Q −1   Q −1   D   −1

9

10

11

12

where

T

T

(36)

T 0  T = 1   0 T2 

(37)

 E {Φ '(ζ )}2  E {Φ '(ζ 1 )Φ '(ζ m )} 1      T1 =     2 E {Φ ''(ζ m )}   E {Φ '(ζ m )Φ '(ζ 1 )} 

(38)

 E { ρ '(v )}2  R ,1     T2 =   2  E { ρ '(vR ,n )}  15


(39)


Page 30 of 57

G  E M = [ E T  DT ]  1 G2   D  

(40)

3

−1  [ BΣ(i | i ) BT ]1,1 E {Φ ''(ζ 1 )}  [ BΣ(i | i ) BT ]1,−1m E {Φ ''(ζ m )}    G1 =      [ BΣ(i | i ) BT ]m−1,1 E {Φ ''(ζ 1 )}  [ BΣ(i | i ) BT ]m−1,m E {Φ ''(ζ m )}  

(41)

4

−1 Q1,1  E { ρ ''(vR ,1 )}    G2 =   −1  Qn ,n E { ρ ''(vR ,n )}  

(42)

1

2

5

and

Here Φ (⋅) uses ζ 2 / 2 , so

6

T1= BΣ(i | i ) BT

(43)

7

 E{ρ '(vR ,1 )}2     T2 =   2  { '( )} E v ρ R ,n  

(44)

8

9

S1 , S 2 are defined as follows: −1  [ BΣ(i | i ) BT ]1,1 E {Φ ''(ζ 1 )}  [ BΣ(i | i ) BT ]1,−1m E {Φ ''(ζ m )}       S1 =   [ BΣ(i | i ) BT ]−m1,1 E {Φ ''(ζ 1 )}  [ BΣ(i | i ) BT ]−m1,m E {Φ ''(ζ m )}  

(45)

= [ BΣ(i | i ) BT ]−1

10

11

−1 Q1,1  E { ρ ''(vR ,1 )}    S2 =   −1  Qn ,n E { ρ ''(vR ,n )}  

So M is stated as 16


(46)

Page 31 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


M =Σ E T [ B (i | i ) BT ]−1 E + DT S 2 D

1

2

(47)

The covariance matrix becomes Σ(i + 1| i + 1)

3

= [ E T ( BΣ(i | i ) BT ) −1 E + DT S 2 D]−1[ E T ( BΣ(i | i ) BT ) −1 E + DT Q −1T2Q −1 D] T −1

(48)

−1

× [ E ( BΣ(i | i ) B ) E + D S 2 D] T

T

4

As deduction above, the robust least squares data reconciliation formulation is Eq.(35)

5

and Eq.(48). The calculation of relevant parameters P, T2 , S 2 are listed in Eq.(33),

6

(34), (44), (46).

7 8

3.4 Procedure of Robust Data Reconciliation using Huber estimator

9

The DDR method eliminates the unmeasured variables by projection matrix and

10

transforms it into a discrete model. We define T as the reconciliation period and

11

sampling frequency is presented by f . The algorithm process is given as following

12

and the flowchart of the algorithm is presented in Fig.2.

13 14

STEP 1: Initialize the sampling parameters period T and frequency f ,

= i 0,= k 1，Xˆ (0) = X ， = c 1.5 ;

15

STEP 2: Transform and discretize the linear dynamic process model as

16

presented in section 3.2, calculate coefficient matrix E , B, D following

17

Eq.(23), Eq.(24) and Eq.(25);

18

STEP 3: Calculate the reconciled value Xˆ (i | i ) , according to the Eq.(35);

19

v(i ) DXˆ (i ) − Z R (i ) ; STEP 4: Calculate the residual= 17



| v(i ) |

Page 32 of 57

1

STEP 5: Calculate the test statistics ε i =

2

c , assign the weight according to Eq.(3);

3

STEP 6: Calculate P (i ) , according to the obtained residual;

4

STEP 7: Recalculate Xˆ (i | i ) ;

5

( k +1) STEP 8: If the Euclidean distance θ= xˆ

6

constant e , go on to the next step, otherwise k= k + 1 and go back to step 4;

7

STEP 9: Calculate S 2 (i ), T2 (i ) , and the covariance matrix of Xˆ (i ) ;

8

STEP 10: i = i + 1 , calculate the dynamic data on the next moment until

9

exceeding given period.

σi

, if it is greater than the constant

xˆ ( k +1) − xˆ k

18


is smaller than a tiny

Page 33 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Initialize the sampling parameters T and f i=0,k=1 , Initialize Xˆ (0)=X (0)，c = 1.5 Transform and discretize model Calculate coeffient matrix E , B, D

Calculate Xˆ (i | i )

Calculate residual = v(i ) DXˆ (i ) − Z R (i )

Calculate P (i )

Recalculate Xˆ (i | i )

N

θ xˆ( k +1)= xˆ ( k +1) − xˆ k ≤ e

k= k + 1

N

Y Calculate S 2 (i ), T2 (i ) and the covariance matrix of Xˆ (i ) i = i +1

i > Tf

Y End

1 2 3

Figure 2. Flowchart of Robust Least Squares Dynamic Data Reconciliation 4. SOLUTION STRATEGY

4

In this section, an improved algorithm is proposed to deal with disadvantages of

5

RLSE by taking both observation and structural redundancy into account. The modified

6

influence function and weight function can help us recognize lower redundant variables 19



Page 34 of 57

1

with gross errors. Meanwhile, the disadvantage from introducing local redundancy

2

would be weakened by a redesigned iterative flow.

3

4.1 Improved Influence Function of Huber Estimator

4

Leverage is a measure of the distance between the independent variable values of an

5

observation and the others. High-leverage points are observations made at extreme or

6

outliers of the independent variables. These points are lack of neighboring observations,

7

which means that the appropriate regression model will approach to that special point.

8

If the errors are homoscedastic, an observation's leverage score determines the degree

9

of noise in the model's misprediction of that observation. Conventional Huber estimator

10

applied variances from prior probability distribution of measurement disregarding the

11

importance of space robustness. It has no limitations on location of spatial distribution

12

of observations, which leads to poor robustness in structure space. In this section, we

13

introduced local redundancy to take spatial characteristics into consideration.

14 15

The definition of H has been given in Eq.(11). The relationship of system redundancy RT

16 = (H ) RT strace =

and local redundancy

H ii

has been defined in Eq.(14)

n

= H ∑ i =1

ii

m . Because

H ii

value just takes part of system

17

redundancy, it couldn’t be criterion of the redundancy of system. Here average local

18

redundancy is taken as a threshold. If the local redundancy of variable is lower than

19

threshold, it is supposed to be weakly-redundant.

20

21

= H m ,ii

1 m trace( H ii ) < n n

For linear system, the least square analytical solution is listed in Eq.(8) and 20


(49)

Page 35 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


v = Z − Xˆ =QAT ( AQAT ) −1 AZ =HZ

1

(50)

2

v stands for residuals, Z stands for measurements, Xˆ stands for reconciled

3

values, The diagonal element of matrix Q is σ i 2 , and σ i stands for standard

4

deviation.

5

From Eq.(50), ii

6

vi =∆xi =H ii xi + ∑ H ii x j

(51)

j ≠i

7

H ii is the local redundancy of ith variable. If H ii is close to zero, the ith variable

8

is a weak redundancy variable. The gross error of the ith measurement couldn’t be

9

reflected from ∆xi , so the weak redundancy variables would hide the gross errors.

10 11

The local robustness properties of general penalized M-estimators are studied via the influence function. The influence function of Huber estimation is Eq.(2).

12

For variables in the [–c, c] range, function values would be proportional to residual

13

errors. Variables beyond this range are regarded as suspect data. If there are low

14

redundancy variables, the absolute value of residuals is much smaller than actual errors.

15

The gross errors would hide themselves in the normal value section. Thus, the influence

16

and weight functions should be reformulated and improved.

17 18

In this paper, local redundancy is introduced into the calculation of residuals. T = Q∆x QA = ( AQAT ) −1 AQ HQ

21


(52)


1

Thus, we let σ r ,i = H ii σ i replace σ i . Since 0 ≤ H ii ≤ 1 , v j / σ r , j ≥ vi / σ i ,

Page 36 of 57

the

2

lower the redundancy is, the smaller H ii is, the larger v j / σ r , j is. Then the influence

3

function is formulated as

−c H ii vi / H ii σ i < −c  ϕ (vi ) vi = − c ≤ vi / H ii σ i ≤ c  vi / H ii σ i > c c H ii

4

5

(53)

By analyzing Eq.(53), the influence of the improvement is to shorten the linear part

6

of former influence function.

7

4.2 Improved Weight Function of Huber estimator

8

The weight function is the key to robust least squares. The improvement of the weight

9

function should consider both robustness and efficiency. Observations would be

10

divided into normal, available and gross error. Based on the three types of observation,

11

the interval should be divided into reservation, descending weight, and rejection

12

intervals respectively.

13

The efficiency and reliability of robust least squares data reconciliation is mostly

14

based on the descending weight interval of influence function w(vi ) . Furthermore, for

15

reliable w(vi ) , the residuals should reflect measurement errors as much as possible,

16

so local redundancy H ii should be large. The reliability of initial iterative solution will

17

be quite important. When k = 0 , x (0) = x , the initial value is just measurement value.

18

The residual is zero, so the first circulation is just conventional least squares estimation.

19

The existence of gross error would lead to bias. Meanwhile, conventional Huber 22


Page 37 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

function has no rejection interval, which leads to worse robustness. Considering above

2

all, the weight function is revised as below.

3

4 5

 vi max(| Yi = |) 0 σ i H ii   vi  = w(vi ) cσ i H ii sign(vi )= / vi Yi | |> c σ i H ii   vi 1 = Yi | |≤ c σ i H ii 

(54)

In this weight function, the largest error statistics Yi will be weighted zero in rejection interval, that means that would be marked and eliminated.

6

In addition, for actual industrial process, measurement values usually have

7

boundaries. These boundaries are always used as constraints in optimizations. Thus,

8

according to the actual production, boundaries are also introduced to decide whether

9

the eliminated gross error in the last circulation is reasonable. If not, the second largest

10

value would replace the maximum to be weighted zero.

11

4.3 Procedure of Novel Robust Least Squares Data Reconciliation

12

The procedure of novel robust least squares data reconciliation is as following:

13

STEP 1: Initialize the reconciliation period T and sampling frequency f , and

14

= i 0,= k 1 ; When k = 1 , the reconciled value xˆ (1)= ( I − QAT ( AQAT ) −1 A) x ,

15

Q = diag (σ i2 ) ; residual v (1) =xˆ (1) − x =− Hx

16

STEP 2: Transform and then discretize the linear dynamic process model,

17

calculate the coefficient matrix E , B, D, H ;

18

STEP 3: Calculate the reconciled value Xˆ (i | i ) , according to Eq.(35); 23



Page 38 of 57

1

v(i ) DXˆ (i ) − Z R (i ) ; STEP 4: Calculate the residual=

2

STEP 5: Calculate P (i ) , according to the obtained residual;

3

STEP 6: Recalculate Xˆ (i | i ) ;

4

STEP 7: Calculate Yi =

5

decrease its weight according to Eq.(54);

6

STEP 8: Check if the value satisfies xL ,i ≤ xi(k) ≤ xU ,i . If not, remove the latest

7

added one in gross error set S , and add the second suspicious one into S ;

8

otherwise k= k + 1 and go back to step 4;

9

( k +1) STEP 9: If the Euclidean distance η= xˆ

| v(i ) | H ii( k ) σ i

, if it is greater than a constant c , then

xˆ ( k +1) − xˆ k

is smaller than a tiny

10

constant e , go on to the next step, otherwise k= k + 1 and go back to step 4;

11

STEP 10: Calculate S 2 (i ), T2 (i ) and the covariance matrix of Xˆ (i ) ;

12

STEP 11: i = i + 1 , calculate the dynamic data on the next moment until

13

exceeding given period.

24


Page 39 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Initialize the sampling parameters T , f and i=0,k=1, initialize the value of x, Q, v Transform and discretize model and calculate matrix E , B, D, H

Calculate Xˆ (i | i )

Calculate residual = v(i ) DXˆ (i ) − Z R (i )

Calculate z (i ) and compare it with constant c to calculate w(i )

Calculate P (i ) k= k + 1

Recalculate Xˆ (i | i ) Y N xL ,i ≤ xi(k) ≤ xU ,i

N

N

Update gross error set S

η x(ˆ k +1)=

xˆ ( k +1) − xˆ k ≤ e

Y Calculate S 2 (i ), T2 (i ) and the covariance matrix of Xˆ (i ) i = i +1

i > Tf

Y End

1 2

Figure 3. Flowchart of the new robust least squares data reconciliation 25



1

Page 40 of 57

5. CASE STUDIES

2

The three methods, least squares estimator, Huber estimator and GRLS estimator are

3

applied to the same process. The following proposed performance indices are used to

4

compare the effectiveness and robustness of the above three algorithms. The smaller

5

the performance index η

6

is, the better is the performance of the algorithm.

The sum of square error (SSE) is defined as below. = SSE

7

n

∑ ( x − xˆ )

2

( 55)

i =1

8

The overall power (OP) is an index for correct detection rate of biased variable.

OP =

9

10 11

12

13

14

number of bias identified correctly number of bias simulated

( 56)

The average type I error (AVTI) is an index for wrong identification of a biased variable.

AVTI =

number of bias identified wrongly number of bias simulated

( 57)

The average type II error (AVTII) is a measure for deviation failed to be identified.

AVTII =

number of bias failed to be identified number of bias simulated

( 58)

15 16

Accumulative reconciled bias: n

17

η1 = ∑ ( Xˆ i − X i,true )T Q −1 ( Xˆ i − X i,true ) i =1

26


(59)

Page 41 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Accumulative measured bias: n

η2 = ∑ [Z R (i) − X i,true ]T Q −1[Z R (i) − X i,true ]

2

(60)

i =1

3 Performance index: n

4

η1 η = = η2

∑ ( Xˆ i =1

i =1

5

− X i ,true )T Q −1 ( Xˆ i − X i ,true )

(61)

n

∑[Z

i

R

−1

(i ) − X i ,true ] Q [ Z R (i ) − X i ,true ] T

5.1 Numerical Case

6

The selected process is described in Fig.1, which was first introduced by Bagajewicz

7

and Jiang3. It shows a simple process which contains 9 flows and 5 nodes, among which

8

U 1 ,U 2 ,U 3 ,U 4 are capacity nodes, U 5 is non-capacity node. W and F represent

9

volume and flow respectively. Sampling interval is 1 second. Measurements are

10

generated from real data adding normal noise. The process of rejection of unmeasured

11

variables and algebraic equations by projection matrix transformation is simplified here.

12

Suppose all the variables are measured. In the process of simulation, c = 1.5 for Huber

13

function and tiny constant e = 1e−6 . Suppose the initial state are

14

W (0) = [100 150 100 200]T ,

15

F (0) = [5 15 20 10 5 10 5 3 2]T .

16

The standard deviation of each variable is 10% of the initial state.

17

The local redundancy of this system is 27



1

H ii = [0.9664 0.9723 0.9854 0.9852 0.0024 0.0310  0.0268 0.0123 0.0025 0.0121 0.0035 0 0]T

.

2

As stated in 3.1, while H ii =0 , i is non-redundant variable and couldn’t be

3

reconciled. So S 8 and S 9 are non-redundant. S1 and S 5 have the greatest

4

influence on the reconciliation of the system. The effect of H ii in the robust

5

reconciliation would be demonstrated afterwards.

6

5.1.1

Scenario 1

7

In this scenario, we suppose that there is no gross error. Fig.4 shows the reconciled

8

data of W2 by both LS and GRLS. The two curves almost coincide. Tab.1 also shows

9

that the effects of both methods are almost the same when there is no gross error.

10 11

12

Figure 4. Reconciliation results comparison of W2 (no gross error) Table 1. Performance comparison in scenario 1

28


Page 42 of 57

Page 43 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


η1

η2

η

LS

150.92

245.81

0.614

GRLS

148.74

245.81

0.605

1 2

5.1.2

Scenario 2

3

The process network is the same as Fig.1. When the measurement of W2 has gross

4

errors from T5 to T7, and the magnitude of bias is 20, the reconciliation results are

5

shown as Fig.5. The performance of conventional RLS are disappointing. Obviously,

6

the GRLS shows much more robustness against gross error. And we can see within

7

three time-adjacent outliers, the GRLS algorithm has good performance as expected.

8 9

Figure 5. Reconciliation results comparison of W2 ( transient gross error)

29



Page 44 of 57

1

The performance index of different method are listed below, the GRLS shows better

2

performance.

3

Table.2. Performance comparison in scenario 2

η1

η2

η

RLS

178.22

293.02

0.6082

GRLS

164.26

293.02

0.5606

4 5

5.1.3

Scenario 3

6

When errors of W1 begin from T45 and last at least 20 time cells, the magnitude of

7

the bias is still 20, the results are shown in Fig.6 and Fig.7. Both W1 and W2 are

8

contaminated, the validity of the traditional robust least squares algorithm has been

9

affected. Our improved algorithm shows better robustness which could also be revealed

10

from Table 3. But from the trend of line chart of GRLS, almost after the gross error last

11

over 10 hours, the gap between GRLS and true data begins to be wider, but the gap

12

between GRLS and RLS begins to be smaller. And for W1 the line of GRLS coincides

13

with RLS near 15 hours later; for W2 the line of GRLS coincides with RLS near 17

14

hours later. So the algorithm maybe break down when gross error lasts over 15 time

15

periods.

30


Page 45 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 2

Figure.6. Reconciliation results comparison of W1 (continuous gross errors)

3 4

Figure.7. Reconciliation results comparison of W2 (continuous gross errors)

31



1

Page 46 of 57

Table.3. Performance comparison in scenario 3

η1

η2

η

RLS

457.25

783.79

0.5834

GRLS

260.74

783.79

0.3327

2 3

5.2 Comparative Case

4

A performance comparative analysis has been conducted. It does not only include

5

redescending M-estimators like Correntropy and adaptive Hampel, but also some

6

filtering methods. The process is the same with above case in Fig.1.

7

5.2.1

Scenario 1

8

In this scenario, we suppose that there is no gross error. Fig.8 shows the reconciled

9

bias of W2 by the listed five methods, including Correntropy, Kalman filter, GRLS,

10

adaptive Hampel and Median filter. The results show that when there is no gross error,

11

the four methods have almost the same performance, except for Kalman filter. The

12

performance is judged by accumulative reconciled bias and list in Tab.4. The simulation

13

time is 200 time cells.

14

Table.4. Accumulative reconciled bias of five methods

η1

Correntropy

Kalman

GRLS

570.2244

604.1043

568.7501

Adaptive Hampel 553.5601

32


Median 561.6709

Page 47 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 2

3

Figure.8. Reconciled bias of W2 by five methods (no gross errors) 5.2.2

Scenario 2

4

In this scenario, we suppose that there are four discontinuous gross errors. The time

5

and margin of gross errors are (40, 45), (90, 25), (120, -35) and (180, -13). Fig.9 shows

6

the reconciled bias of W2 by the same five methods. The results show that when gross

7

error occurs, the five methods show different performance. The accumulative

8

reconciled biases are listed in Tab.5. Kalman filter has poor robustness with single gross

9

error. Other method had almost the same performance, except for adaptive Hampel. It

10

demonstrates that GRLS has good performance against gross error, and at least has the

11

same robustness with Correntropy and median filter.

33



Page 48 of 57

1 2

3

Figure.9. Reconciled bias of W2 by five methods (discontinuous gross errors) Table.5. Accumulative reconciled bias of five methods

η1

Correntropy

Kalman

GRLS

568.5634

642.2739

566.1528


Median 566.9222

4 5

5.2.3

Scenario 3

6

In this scenario, we suppose that there are continuous gross errors. Firstly, deviation

7

with amplitude around 25, and lasting for three time cells are added to measurements

8

of W2 . Fig.10a shows the reconciled bias of W2 by five methods. The accumulative

9

reconciled biases are listed in Tab.6. When the disturbance lasts for just three time cells,

10

Kalman filter has performed worse, but others perform well.

34


Page 49 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 Figure.10a. Reconciled bias of W2 (transient gross errors)

2 3

Table.6. Accumulative reconciled bias of five methods

η1

Correntropy

Kalman

GRLS

570.3373

672.9898

564.3034


Median 562.7660

4 5 6

When the disturbance lasts longer, adaptive Hampel and Median filter meet their break-down point as shown in Fig.10b.

7 8

Figure.10b. Reconciled bias of W2 (gross errors last for four time cells) 35



1

Page 50 of 57

Table.7. Accumulative reconciled bias of five methods Correntropy

Kalman

GRLS

573.4068

699.7748

565.4687

η1


Median 643.9940

2

When the distance lasts much more longer, the robustness of adaptive Hampel is

3

unstable. Hampel couldn’t detect all the gross errors sometimes. But the tendency of

4

GRLS is still stable as shown in Fig.10b and 10c. And the performance of GRLS is a

5

little better than Correntropy.

6 Figure.10c. Reconciled bias of W2 (gross errors last for six time cells)

7 8

Table.8. Accumulative reconciled bias of five methods Correntropy

Kalman

GRLS

577.1355

776.1355

572.8728

η1


Median 716.4984

9 10 11

Above all, most of the time, GRLS has better robustness than listed M-estimators and filters. 36


Page 51 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

5.3 Industrial Case

2

The performance of the proposed GRLS algorithm has been tested using a classical

3

industrial model first proposed by Serth(1986)33, which was a methanol synthesis unit

4

of a large chemical plant. The size and structure of this case are similar with many

5

industrial systems. A network representation of the system consisting of 12 nodes and

6

28 streams is shown in Fig. 11.

v4

v3

e10

e4

e1

e5 e9

v2

e6

v9

e12 e11

e17

e16

e22

v1

e3

v5

e15

v8

e2

e13 e18

v10

e27

e19

e7 v7 e8

e24

e23 e14

v6

e25

e21

e26

e20 v11

e28

7 8

Figure 11. Network representation of steam system for a methanol synthesis unit

9

Assuming that e7 , e8 , e11 , e26 are contaminated by 50% gross errors. True values,

10

measurements, variances, redundancy, and local redundancy are listed respectively.

11

Reconciliation results of LS, RLS, RLS-LR(robust least squares considered local

12

redundancy) and GRLS are compared.

13

The detailed results are shown in the supporting information. It can be concluded that

14

the LSE algorithm has the worst robustness. Our novel robust least squares iterative

15

algorithm performs much better than other two. Traditional MT test based on LSE

16

recognized 24 gross errors, but only four of them are right. Thus, the MT test based on

17

LSE cannot locate the gross error correctly. After comparing the last three robust 37



Page 52 of 57

1

iterative algorithms, it is found that the RLSE detected that e8 , e26 contain gross errors,

2

but e7 , e11 are missing. We found that the local redundancy of e7 , e11 are relatively

3

lower, and the test statistics cannot reflect the gross error. The robust least squares

4

algorithm

5

e1 , e7 , e8 , e11 , e22 , e26 , e27 . Compared with the actual gross errors, there was no missing

6

report, but the false report increased first kinds of errors. Because 0 ≤ H ii < 1 , the

7

introduction of H ii would decrease the threshold of error detection and increase the

8

detective rate of gross error, but at the same time, it improves the error statistics of other

9

variables to some extent. The last column shows that our improved algorithm has the

10

considering

space

robustness

has

detected

gross

errors

in

best performance. It located the four gross errors accurately.

11

After analyzing convergence, RLS iterated 15 times, RLS-LR iterated 36 times and

12

GRLS iterated 9 times, our algorithm costs the least computation resources and has

13

good convergence.

14

To prove the statistical significance of the algorithm, the industrial case has been

15

simulated for 100 times. And the average AVTI and AVTII are listed below. From the

16

statistical results, since the conventional robust least square method has no

17

consideration of space construction and redundancy, it has high possibility of

18

misjudgment. Due to the introduction of local redundancy, RLS-LR has higher AVTI.

19

After improving the algorithm flow, GRLS has added bounding constraints for

20

variables into the iteration, and the failure rate of gross error detection is relatively much

21

lower.

22

Table.9. AVTI and AVTII of four methods

38


Page 53 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


AVTI

AVTII

LS

83.33%

0.25%

RLS

0.13%

49.90%

RLS-LR

12.50%

0.08%

GRLS

0.083%

0.042%

1 2

6. CONCLUSION

3

In this work, we propose a novel robust least squares estimator for linear DDR. The

4

algorithm takes the structural robustness into consideration via using local redundancy

5

in influence function and rejection interval in weight function. We redesign an iterative

6

flow to weaken the disadvantage from introducing local redundancy. Combined with

7

the improvement and taking the advantage of online filtering, a novel robust dynamic

8

data reconciliation algorithm is proposed. The novel estimator has almost the same

9

effectiveness with RLSE without gross errors, but much better performance when the

10

measurements are contaminated by gross errors. Compared to the conventional LSE or

11

RLSE, the GRLSE improves the accuracy of DDR, and its application in industrial

12

systems could improve the quality and accuracy of large-scale measurement. More

13

challenges about implementing the method in complex industrial processes will be

14

addressed in our future works.

15 16 39



Page 54 of 57

1

7. SUPPORTING INFORMATION

2

8. The Supporting Information is available free of charge on ACS Publications website

3

at http://pubs.acs.org/.

4

1. Mathematical formulation: Some formula deduction extended from Section 3.2

5

2. Data sheet of industrial case study, including

6

Table.1. Comparison of different methods for methanol synthesis unit

7

Table.2. Comparison of different methods using SSE, OP, AVTI

8

3. Nomenclature for important symbols

9 10

9. TABLE OF CONTENTS

Composition

Improvement

M-Estimators

Local redundancy

Influence function Hampel Estimator Robust Least Square Estimator

Elimination interval

Weight function

Online filtering

Correntropy Estimator

. . .

Iterative process

Improved iterative process

11 12

40


General Robust Least Square Estimator

Page 55 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

AUTHOR INFORMATION

2

Corresponding Author

3

*Tel: +86 571 87953145. Email: [email protected].

4 5

ACKNOWLEDGEMENTS

6 7 8 9 10 11

The authors gratefully acknowledge financial support by Smart Manufacturing Comprehensive Standardization and New Mode Application Project of MIIT (No.2017ZJ-003) and Science Fund for Creative Research Groups of NSFC (Grant No.61621002). REFERENCES

12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

1. Almasy, G. A., Principles of dynamic balancing. AIChE Journal 1990, 36 (9), 1321-1330. 2. Darouach, M.; Zasadzinski, M., Data reconciliation in generalized linear dynamic systems. AICHE Journal 1991, 37 (2), 193-201. 3. Bagajewicz, M. J.; Jiang, Q., Integral approach to plant linear dynamic reconciliation. AIChE Journal 1997, 43 (10), 2546-2558. 4. Bagajewicz, M. J.; Jiang, Q., Comparison of steady state and integral dynamic data reconciliation. Computers & Chemical Engineering 2000, 24 (11), 2367-2383. 5. Bai, S.; Thibault, J.; McLean, D. D., Dynamic data reconciliation: Alternative to Kalman filter. Journal of Process Control 2006, 16 (5), 485-498. 6. Veverka, V. V.; Madron, F., Material and energy balancing in the process industries: From microscopic balances to large plants. Elsevier: 1997; Vol. 7. 7. Rollins, D. K.; Devanathan, S., Unbiased estimation in dynamic data reconciliation. AIChE journal 1993, 39 (8), 1330-1334. 8. Xu, H.; Rong, G., A new framework for data reconciliation and measurement bias identification in generalized linear dynamic systems. AIChE Journal 2010, 56 (7), 1787-1800. 9. Johnston, L. P.; Kramer, M. A., Maximum likelihood data rectification: Steady‐ state systems. AIChE Journal 1995, 41 (11), 2415-2426. 10. (a) Fair, R. C., On the robust estimation of econometric models. In Annals of Economic and Social Measurement, Volume 3, number 4, NBER: 1974; pp 667-677; (b) Albuquerque, J. S.; Biegler, L. T., Data reconciliation and gross‐error detection for dynamic systems. AIChE Journal 1996, 42 (10), 2841-2856. 11. (a) Hampel, F. R., The influence curve and its role in robust estimation. Journal of the american statistical association 1974, 69 (346), 383-393; (b) Nicholson, B.; López-Negrete, R.; Biegler, L. T., On-line state estimation of nonlinear dynamic systems with gross errors. Computers & Chemical Engineering 2014, 70, 149-159. 41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

12. (a) Chen, J.; Peng, Y.; Munoz, J. C., Correntropy estimator for data reconciliation. Chemical Engineering Science 2013, 104, 1019-1027; (b) Zhang, Z.; Chen, J., Correntropy based data reconciliation and gross error detection and identification for nonlinear dynamic processes. Computers & Chemical Engineering 2015, 75, 120-134. 13. Arora, N.; Biegler, L. T., Redescending estimators for data reconciliation and parameter estimation. Computers & Chemical Engineering 2001, 25 (11), 1585-1599. 14. Huber, P. J., Robust statistics. 1981. 15. Özyurt, D. B.; Pike, R. W., Theory and practice of simultaneous data reconciliation and gross error detection for chemical processes. Computers & Chemical Engineering 2004, 28 (3), 381-402. 16. Zhang, Z.; Shao, Z.; Chen, X.; Wang, K.; Qian, J., Quasi-weighted least squares estimator for data reconciliation. Computers & Chemical Engineering 2010, 34 (2), 154-162. 17. Wang, D.; Romagnoli, J., A framework for robust data reconciliation based on a generalized objective function. Industrial & Engineering Chemistry Research 2003, 42 (13), 3075-3084. 18. Zhou, J.-w.; Huang, Y.-c.; Yang, Y.-x., Theory of Robustified Least Squares Method [M]. Wuhan: The Press of Huazhong University of Science and Technology 1997, 71 (95), 104-118. 19. Huber, P. J., Robust estimation of a location parameter. The Annals of Mathematical Statistics 1964, 35 (1), 73-101. 20. Huber, P. J., Robust statistics. In International Encyclopedia of Statistical Science, Springer: 2011; pp 1248-1251. 21. Ni, X. S.; Huo, X., Another look at Huber's estimator: A new minimax estimator in regression with stochastically bounded noise. Journal of Statistical Planning and Inference 2009, 139 (2), 503-515. 22. Hampel, F.; Hennig, C.; Ronchetti, E., A smoothing principle for the Huber and other location M-estimators. Computational Statistics & Data Analysis 2011, 55 (1), 324-337. 23. Kim, T.-H.; Muller, C., Two-stage Huber estimation. Journal of Statistical Planning and Inference 2007, 137 (2), 405-418. 24. Jin, S.; Li, X.; Huang, Z.; Liu, M., A New Target function for robust data reconciliation. Industrial & Engineering Chemistry Research 2012, 51 (30), 1022010224. 25. Llanos, C. E.; Sánchez, M. C.; Maronna, R. A., Classification of Systematic Measurement Errors within the Framework of Robust Data Reconciliation. Industrial & Engineering Chemistry Research 2017, 56 (34), 9617-9628. 26. Kretsovalis, A.; Mah, R., Observability and redundancy classification in generalized process networks—I. Theorems. Computers & Chemical Engineering 1988, 12 (7), 671-687. 27. Carpentier, P.; Cohen, G., State estimation and leak detection in water distribution networks. Civil Engineering Systems 1991, 8 (4), 247-257. 28. Madron, F.; Veverka, V., Optimal selection of measuring points in complex plants by linear models. AIChE journal 1992, 38 (2), 227-236. 29. Maronna, R.; Arcas, J., Data reconciliation and gross error diagnosis based on regression. Computers & Chemical Engineering 2009, 33 (1), 65-71. 42


Page 56 of 57

Page 57 of 57 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 2 3 4 5 6 7 8

30. Bagajewicz, M. J.; Sanchez, M. C., Design and upgrade of nonredundant and redundant linear sensor networks. AIChE Journal 1999, 45 (9), 1927-1938. 31. Xiuping, W.; Gang, R.; Shuqing, W., Key-variable Analysis in Data Reconciliation. Chinese Journal of Scientific Instrument 2001, 22 (3), 235-239. 32. Crowe, C.; Campos, Y.; Hrymak, A., Reconciliation of process flow rates by matrix projection. Part I: linear case. AIChE Journal 1983, 29 (6), 881-888. 33. Serth, R.; Heenan, W., Gross error detection and data reconciliation in steam‐ metering systems. AIChE Journal 1986, 32 (5), 733-742.

9 10

43


Novel Robust Least-Squares Estimator for Linear ... - ACS Publications

Recommend Documents