Application of Active Set Method for Soft Sensor Model Identification of Crude Oil Distillation Process

Abstract

The problem of identification of the soft sensors is considered. Using the active set method the problem of least squares with inequality constraints on the variables has been solved. As a result of using soft sensors, obtained taking into account constraints on the model coefficients, their efficiency as compared with soft sensors obtained without constraints is shown. The approach to solve the soft sensor model obtaining problem, taking into account inequality constraints on the model coefficients, is described [1]. The obtained results are tested on industrial oil fractionation process for atmospheric distillation column. The soft sensor model coefficients for prediction of the key product quality are identified. The technological plant with several measured inputs u1,u2 ,,uN and one output y( ) is considered. The measured technological parameters (pressure, temperature, flow) are utilized as inputs. In practice, the amount of available process measured parameters to predict the quality of a product is much more than the required number of parameters. A priori knowledge of technological process allows to select necessary parameters. The problem of the soft sensor (SS) evaluation which is best predicting quality of products of crude distillation technological process is considered. The model of the soft sensor is obtained in the form of linear regression model for solution of the problem [2]: y( )  b0  b1u1    b2u2  ...  bNuN   ,

Keywords

• parametric identification
• optimization
• soft sensor
• constraints
• distillation column

-

Introduction and statement of the problem ( 1 ) Copyright © by the paper's authors. Copying permitted for private and academic purposes.

In: A. Kononov et al. (eds.): DOOR 2016, Vladivostok, Russia, published at http://ceur-ws.org where b j – j-th model coefficient, j  0,1,..., N , b0 – constant term, N – the number of input variables, 

- irregular timepoints of output measurement:  1, 2 , 3 ,... ,  i   i1  0   , i  2 ;  1   0   ;  0 - constant term;  - random component is limited by specific range.

The determination coefficient (the part of explained deviations variance of the dependent variable from its mean value): and root mean squared error (RMSE):

R 2  1  i( yi  yi )2 i( yi  y )

a 2 RMSE  iM1 ( yi  yi ) 2 / M 1/ 2 , ( 2 ) ( 3 ) are used as criteria of identification on a given time interval, where yi - the measured value of the output variable, yi - the value is obtained based on the SS, y a - the mean value of the measured output variable, M - the number of output measurement. The model is more consistent if the closer to unity the value of the coefficient of determination R 2 , or the closer to zero the value of the RMSE. 2

The proposed algorithm problem solution

Let u  [1, u1( ), u2 ( ),..., uN ( )]T be a combined vector of measured input variables, b = [b0 , b1,..., bN ]T - vector of coefficients of the same dimension, the components of which reflect the contributions of the respective input variables. Then the equation ( 1 ) takes the form:

We form the vector Y of dimension q from the output value y:

y  uT  b .

Y  [ y( 1 ), y( 2 ),..., y( q )]T   U       1 1 and we write the matrix equation:

We introduce error function:

Y  Ub . where Y is the actual measurement of output, and minimize the objective function: We obtain estimates of the parameters b by least squares method:

The multicollinearity case, which occurs when there is almost a linear relationship between inputs, is considered. In this case the matrix C  UT U is close to singular, so it is the smallest eigenvalue min  0 and the condition number is infinitely increased and causing the instability of the solution (5). If the min  0 then it corresponds the strict multicollinearity [3]. In order to obtain a stable solution of the equation (5) it is necessary to reduce the condition number of the matrix C, for example, by adding thereto a diagonal matrix B  kI (k> 0). Then the solution is found in a class of ridge parameter estimates: b  (UT U  kI)1UT Y . (6)

The quality, obtained by (5-6) models, depends on the number of available output measurements. The length of training sample is often insufficient to obtain reliable results. Also, the available data contain significant measurement error of inputs and outputs, unmeasured influences. Taking into account constraints on the model coefficients b j allows to avoid these problems. When taking into account constraints on coefficients at input, the problem of least squares with simple constraints on the variables is solved:

The solution of the problem (7) is obtained by the active set numerical method [4]. The given constraints are reduced to the form:

E  Y  Y  Y  Ub , Ψ  E2  (Y  Ub)2 .

b  (UT U)1 UT Y .

min  Y  Ub2 bmin  b  bmax .

Ab  bˆ , ( 4 ) (5) (7) where  1  0    A   0  1   0     0

Constraint aTi b  bˆi call active in acceptable point b if aTi b  bˆi , and inactive if aTi b  bˆi , a Ti - i-th row of A .

The sufficient minimum conditions for simple constraints are as follows: 1. bmin  b*  bmax , b mFRin  b FR  b mFRax

* 2. UTFR (Y  Ub* )  0 3. λ min  Umin (Y  Ub * ) , λimin  0 , i  1,, t min ,

T λmax  UTmax (Y  Ub* ) , λ miax  0 , i  1,, t max 4. UTFR U FR is positive definite, where b* - the minimum point of the solution of problem (7); subscript FR indicates that in the vector and matrices the elements and columns with index numbers corresponding to the index numbers of b elements, that have not met the boundary values (7), are used; subscript min, max, indicates that the in matrix only the columns with index numbers corresponding to the index numbers of b elements, taking the appropriate minimum or maximum boundary value, are used. tmin,, tmax.- number of active upper and lower limits respectively; λ min , λ max - vectors of Lagrange multipliers corresponding to the lower and upper active constraints.

To start the method of the active set it is necessary to determine the starting point using (6).

The minimum point b* for the search algorithm for iteration k is: 1. Performance verification of the stop conditions. (Reaching the performance errors of conditions (8), constraints on the number of iterations). 2. Selection of a logic branch. Does it make sense to remove any constraint of the set of active constraints list. The condition of performance of a condition 3 in (8) is checked. If the condition is not satisfied for some of the vector element, constraint is excluded from the list of active constraints. 3. The calculation of the search direction pk . Like ( 4 ) solves the problem min Y  Ubk  UFRpFR 2 . Calculate the non-zero  N 1 tk  - dimensional vector (8) pFR and the direction of search pk  AT FR p FR , where tk - the number of active constraints on k iteration. 4. Calculate the step length  k . From  b FR   Ψ p FR   bˆ FR diagonal matrix  b FR   p FR  Ψ is calculated, bˆ FR - consists of the elements bˆ which aren't active constraints, elements bˆ , opposite boundary values in (7) for constraints in the active set, are excluded from bˆ FR . The  k  min Ψii is an available minimum positive step from bk along pk . The index j of minimum positive diagonal element Ψ is remembered. If  k  1, then  k  1, otherwise  k   k . 5. Constraint is added in the list of active constraints. If  k   , then j constraint ˆ bFR becomes active, it is necessary to add to the list. 6. Recalculation approximation. bk1  bk  kpk is calculated, and return to step 1 of the algorithm is carried out.

The influence of the process dynamics on the quality of the products is taken into account by the dynamic SS. The predictive model is represented as a sum of convolutions of plant inputs and a finite impulse response (FIR) hi (discrete analogues of the first degree Volterra kernels): y( )  h0  kn101h1 (k 1)u1 (  k)  kn201h2 (k 1)u2 (  k)  ...  knN01 hN (k 1)u N (  k), (9) where h0 – constant term, 

- irregular timepoints of output measurement:  1 , 2 , 3 ,... ,  i   i1  0   , i  2 ;  1  0   ;  0 - constant component;  random component is limited by the specific range.

Let u  [1, u1 ( ),..., u1 (  n1  1),..., uN ( ),..., uN (  nN  1)]T - the combined vector of measured input variables of dynamic SS (DSS) with dimensionality q  1  kN1 nk where nk - is a number of values of k-th input, h  h0 , h1 ( 1 ), ..., h1 (n1 ), ..., hN ( 1 ), ..., hN (nN ) т - vector FIR of the same dimension, the components of which reflect the contributions of the respective input variables of DSS. Then the equation (9) takes the form:

y  uT  h .

We form the vector Y of dimension q from the output value y: and the matrix U, containing the measured inputs uj, corresponding to output value y from (9): 1 u1(1)  u1(1  n1 1)  uN(1)  uN(1  nN 1) 1 u1( 2)  u1( 2  n1 1)  uN( 2)  uN( 2  nN 1). U             1 u1( q)  u1( q  n1 1)  uN( q)  uN( q  nN 1) Then write the matrix equation: Y  Uh .

We introduce the error function:

E  Y Y  Y Uh , whereY is the actual measurement of output, and minimize the objective function: The constraints on transient response components are written as: Ψ  E2  (Y  Uh)2.

smin  s  smax , where s  s1( 1 ),..., s1(n1),..., sN( 1 ),..., sN (nN )т , smin s1min,...,smNin т , smax s1max ,...,smax т.

N

The transient response components s are related with the components of the impulse response h by the relations:

sjk ik1hj(i), j 1, 2, , N, k 1, , nj.

The constraints (11) are reduced to:

Ah~  sˆ , where h~   h1(n1)  , sˆ   smin  .

   smax   hN( 1 )  (10) (11) (12) (13) where h* - is the minimum point, the solution of problem (10) with constraints (13); subscript ACT indicates that in vector, matrix only the elements, rows with index numbers corresponding to the elements index numbers of active constraint in (13) are used; t - number of active constraints, λ - the vectors of Lagrange multipliers corresponding to the active constraints, Z - matrix the columns of which is basis of the feasible direction of search for equality constraints (13). The matrix Z is formed by the variable-reduction technique. [4].

In order to start the method of the active set it is necessary to determine the starting point (using a solution of the problem (10) without any constraints, with subsequent correction of coefficients hi that does not fall under the constraints (13)).

The search algorithm of minimum point h* for iteration k is: 1. Performance verification of the stop conditions (reaching the performance errors of conditions (14), constraints on the number of iterations). 2. Selection of a logic branch. Does it make sense to remove any constraint of the set of active constraints list. The condition of performance of a condition 3 in (14) is checked. If the condition is not satisfied for some of the vector element, constraint is excluded from the list of active constraints and it is need to recalculate Zk . 3. The calculation of the search direction pk . Like ( 4 ) solves the problem min Y  Uh k  UZ k p Z 2 . Calculate the non-zero 1  kN1 nk  tk  - dimensional vector p Z and the direction of search pk  ZkpZ , where tk - the number of active constraints on k iteration. 4. Calculate the step length  k . From Ah~  Ψp~k   sˆ diagonal matrix Ψ is calculated. Calculated  k  minΨii  - a minimum non-negative available step from hk along pk , where i is the index number of element, which is not active constraints in (13) and is not element of opposite boundary values (11) for constraints in the active set, p~ k consists from elements of pk without first element. The index j of minimum positive diagonal element of Ψ ii is remembered. If  k  1 , then  k  1, otherwise  k   k . 5. Constraint is added to the list of active constraints. If  k   , then j constraint sˆ FR becomes active and recalculate Zk . 6. Recalculation approximation. The hk1  hk  p k k is calculated and return to step 1. The influence of constraints of parameters of soft sensors model

A priori knowledge about the process and industrial step test (when the value of the one control variables is changed at fixed others) allow to define the value of constraints.

In order to investigate the influence of constraints on the quality of the obtained static soft sensor model we compare solutions of the equations (6) and constrained optimization problem (7). One and the same value of the ridge coefficient is used.

The Fig. 1 and Table 1 show the results of the performance of the static models on the test sample for the bubble-point temperatures of the target product (BP) when a model obtained on the training sample, consisting of a number of measurements specified in the Table 1. In the verification sample of models the number of measurements is equal to 420.

In order to investigate the influence of constraints on the quality of the obtained dynamic soft sensor model we compare solutions of the optimization problem (10) and optimization problem (10) with constrains (13). One and the same value of the ridge coefficient is used.

The Fig. 2 and Table 2 show the results of the performance of the dynamic models on the test sample for the dew-point temperatures of the target product (DP) when a model obtained on the training sample, consisting of a number of measurements specified in the Table 2. In the verification sample of models the number of measurements is equal to 650.

The estimation of improvements of the prediction quality by the criterion RMSE of static model obtained with the constraints on the parameters of SS is 100( 1,79 1,249 ) / 1,79  30% compared to the model without constraints. The estimation of improvements of the prediction quality by the criterion RMSE of dynamic model obtained with the constraints on the parameters of SS is 100 · ( 1,456– 1,074 )/ 1,456  25 % compared to the model without constraints.

Conclusion

The using the method of the active set, taking into account constraints on the model coefficients can improve quality of the evaluated SS models.

The test of the proposed approach to solving the problem of obtaining a soft sensor model for industrial crude oil distillation unit is showed that the decrease root mean square error on the test sample can be not less than 25%.

References

1. Bakhtagze. N.N. : Virtual Analyzers: Identification Approach . In: Automation and Remote Control . Vol. 65 , issue 11, pp. 1691 - 1709 . ( 2004 )

2. Draper N. , Smith H. : Applied regression analysis . M .: Finance and Statistics ( 1986 )

3. Bolshakov A.A. , Karimov R.N. : Methods of processing of the multidimensional given and time numbers . M .: Hotline Telecom ( 2007 )

4. Gill , P.E. , Murray , W. , and Wright , M.H. : Practical Optimization, London: Academic, 1981 . Translated under the title Prakticheskaya optimizatsiya , Moscow: Mir ( 1985 )