Accepted Manuscript

Title: Modifier Adaptation methodology based on transient and static measurements for RTO to cope with structural uncertainty

Authors: T. Rodríguez-Blanco, D. Sarabia, J.L. Pitarch, C. de Prada

PII: DOI:

Reference:

S0098-1354(17)30272-7

http://dx.doi.Org/doi:10.1016/j.compchemeng.2017.07.001 CACE 5857

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

26-9-2016 29-6-2017 1-7-2017

Please cite this article as: Rodríguez-Blanco, T., Sarabia, D., Pitarch, JL., & de Prada, C., Modifier Adaptation methodology based on transient and static measurements for RTO to cope with structural uncertainty.Computers and Chemical Engineering http://dx.doi.org/10.1016/jxompchemeng.2017.07.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Modifier Adaptation methodology based on transient and static measurements for RTO to cope with structural uncertainty

T. Rodríguez-Blanco3*, D. Sarabiab, J. L. Pitarcha, C. de Pradaa

a Dpt of Systems Engineering and Automatic Control, University of Valladolid, Doctor Mergelina s/n, 47011, Valladolid, Spain (tania.rodriguez@autom.uva.es, jose.pitarch@autom.uva.es, prada@autom.uva.es)

b Dpt of Electromechanical Engineering, University of Burgos, Avda. Cantabria s/n, 09006, Burgos, Spain (dsarabia@ubu.es)

Highlights

• Modifier Adaptation methodology uses plant measurements to bring the process to the real optimum despite the presence of uncertainty.

• A Modifier Adaptation approach is proposed to speed up the convergence to the optimum using transient measurements to estimate the process gradients.

• The described method is valid for parametric and structural uncertainty.

• The approach is illustrated through a simulated depropanizer distillation column.

* Corresponding author

Abstract

Optimal process operation is carried out by a Real-Time Optimization (RTO) layer which is not always able to achieve its targets due to the presence of plant-model mismatch. To overcome this issue, the economic optimization problem solved in the RTO is changed following the Modifier Adaptation methodology (MA), which uses plant measurements to find a point that satisfies the necessary optimality conditions (NCO) of an uncertain process. MA proceeds by iteratively adjusting the optimization problem with first and zeroth order corrections, calculated from steady state information at each RTO execution. This implies a long convergence time. This paper presents a new method based on a recursive identification algorithm to estimate process gradients from transient measurements to speed up the convergence of MA. The proposed approach is implemented in a simulated depropanizer column that incorporates a simplified model in the RTO, reducing by 8 the convergence time compared with traditional MA.

Keywords: Real-Time Optimization, Modifier Adaptation, uncertainty, distillation columns, MPC.

1. Introduction

The management of large scale systems, such as many in the petrochemical industry, consists of making decisions in order to satisfy process specifications and constraints on many variables. Ideally, these decisions should be optimal with respect to several aspects; such as efficiency, economy or environment, and are normally based on the use of large models and optimization methods.

RTO consists of an optimization layer that operates above the control layer and makes decisions on a time scale of hours by explicitly considering economic objectives. The optimum operating point obtained by the RTO layer is passed to lower-level controllers that include basic control and model predictive control.

However, optimal operation is not guaranteed since the process models used by RTO are always inaccurate, so the optimum computed from the model may not be the same as the optimum of the process. In addition, the layers of control structure are not always considered explicitly in the RTO model, creating inconsistencies between them that affect the final result.

Usually, the RTO layer uses a steady-state model of the process to make decisions, where the RTO problem is formulated as follows: min ^(u,P)

" (1) st gi(u,P) < 0 i = 1,...,ng

L ^ ^ U

u < u < u

where ^ is the cost function to be minimized, gi represents the constraints, P the model parameters, and u the decision variables which present lower and upper limits uL and uU.

In order to be able to compute the plant optimum, the NLP problem associated with the RTO must satisfy a couple of conditions:

(a) A necessary condition is that the Karush-Kuhn-Tucker (KKT) conditions of the model must coincide with those of the plant.

(b) A sufficient condition is that the model matches all the gradients of the plant at all points of interest. This results from the equality between the necessary optimality conditions of the process and the model and it is this which is the basis of the MA methodology, as will be explained later (Naysmith and Douglas, 1995).

Several proposals have been developed to cope with the uncertainty already mentioned and to drive the process to its real optimum point. The classical approach consists in an iterative two-stage algorithm involving a parameter estimation step (to update uncertain model parameters) followed by an economic optimization that is solved to obtain new decision variables. However, this formulation works well only if there is little structural plant-model mismatch and the changing operating conditions provide sufficient excitation to estimate the uncertain parameters (Yip and Marlin, 2004). Such conditions are rarely met in practice. For this reason, the classical approach of RTO will not satisfy the NCO and the real optimum will not be reached using this method.

From the two-step scheme, a new approach was developed by Roberts (Roberts, 1979) who incorporated information regarding plant gradients, adding a modifier to the economic cost function that emerges from the equality of the NCO for the real process and the static model used in the RTO layer. It implies that the modifier is computed from the difference between the gradient of the real cost and the model one. This method was called "integrated system optimization and parameter estimation" (ISOPE). In 2002, Tatjewski proved that the convergence to the optimum point does not depend on parameter estimation, but on the equality between the outputs of the process and the model at each RTO iteration (Tatjewski, 2002). For this reason, he introduced a new modifier that takes into account the difference between these outputs. New modifiers were also defined by Gao and Engell for process dependent constraints (Gao and Engell, 2005). The resultant modified RTO problem (2) in which and gM,i are the modified cost function and constraints, u*k-1 is the input applied in the previous steady state, that is, the optimal solution of the previous RTO, and the subscript "P" indicates that the gradient of the cost function, the gradient of the constraint and the constraint itself are evaluated from the process measurements:

min = ^(u, P) + ^kk (u -uk_i)

st gMi = gi (u,P) + Yk,i (u -uk-1) + sk,i < 0 uL < u < uU

i = 1,...nt„

where:

^T = Vu^P |u* . -vu4,*

Y ki =vu gp,i |uLi-Vu gi lu^ i = K../ig (3)

sk,i = gp,i (uk-1) - gi (uk-1)

Figure 1 presents the general formulation of the MA methodology.

These methods require the computation of model and experimental gradients. Whereas there are tools to compute the model gradients accurately, the plant gradients are more difficult to estimate. In particular, Brdys and Tatjewski, 1994 estimates experimental gradients from past operating points generated by the previous RTO iterations by using the definition of directional derivative and proposed a dual control approach that implicitly considers the requirement of the

gradient estimation in the set point optimization. To ensure that gradients are obtained accurately, a new constraint (K-1(Sk) > a) is added to the optimization problem, where k represents the condition number of Sk that is the matrix formed by the vectors of differences of the decision variables with respect to nu samples before, as (5) shows, and the parameter a indicates the minimum degree of excitation. This constraint represents the dual characteristic of the method: while the rest of the optimization aims to converge to the optimum of the modified model (primal objective), the dual constraint ensures that, in the next RTO iteration, the system will have enough excitation to estimate the process gradient adequately (dual objective). Marchetti et al., 2010 extended this approach by explicitly upper bounding the norm of the gradient estimation error formulating the Dual Modifier Adaptation methodology (DMA). Gao and Engell, 2005 used the dual constraint to decide whether to perturb the process additionally to estimate gradients accurately.

Assuming that there is a collection of nu+1 points Uk, Uk-i,..., Uk-nuapplied in the past RTO iterations, the vectors of differences with respect to previous points sk,; are defined by (4):

st, = U

Vi = 1...nu (4)

Supposing that the vectors sk,; are linearly independent, it is possible to formulate a nonsingular square matrix Sk e M(n x nu) whose condition number will give us the degree of excitation of the process:

Sk =kz .... sknu I1 (5)

The lower bound of the dual constraint (a) must be carefully chosen, since the convergence to the real optimum is very sensitive with this parameter and also because a big value implies an increase in the excitation of the system that can lead to an operating point that does not satisfy the NCO of the process.

This dual constraint will be also added to the novel MA approach presented in this paper, which is described in the section 2.2. This method is based on the direct estimation of the process gradients over the transient, for this reason, the dual constraint is used to ensure enough process excitation to accurately estimate the required gradients.

The calculation of experimental gradients can be avoided by using a different formulation called Nested-Modifier Adaptation (NMA) (Navia et al., 2015). This method uses a nested optimization architecture with a gradient-free optimization algorithm, for example, the Nelder-Mead algorithm, to directly update the modifiers, iterating with them over the modified optimization until the optimum of the process is found. In this way, the process gradient estimation is replaced by another method that takes into account the minimization of the cost function measured directly from the process.

One of the main disadvantages of MA is the necessity to wait for the steady state of the process before updating the modifiers. In many real applications, such as the operation of distillation columns or similar process units, the transients can last for several hours, so the convergence of MA can be very slow and the real optimum may only be achieved after several days of operation. During this period of time, the operating conditions or the plant-model mismatch may change and the method will not converge to an operating point that satisfies the NCO of the process. This issue makes the application of this methodology impractical in these cases. Besides, and not less important, this kind of method does not guarantee a feasible path during the convergence to the real optimum, only when this point is achieved are all the process constraints satisfied. So, the longer the time to achieve the optimal operation the bigger the probability of violating operational constraints that may lead to the loss of benefits or to the violation of safety and environmental constraints that may result in a dangerous operation.

In order to speed up the convergence of MA methodology for slow dynamic processes and reduce the risk of process constraint violations, several researchers have suggested the use of transient measurements to estimate the variables required by the steady-state optimization. This idea was pursued by Zhang and Roberts in 1990 (Zhang and Roberts, 1990), who combined the ISOPE scheme with a linear dynamic model identification to compute process gradients for the steady state optimization of nonlinear constrained processes with slow dynamics, but this work did not address the problem of shortening the gradient estimation time. In contrast, in 2014, François and Bonvin (François and Bonvin, 2014) proposed an approach that uses transient measurements to compute process gradients by the neighbouring extremal method which relies on the accuracy of the linearization resulting from a variational analysis of the nominal model. However, none of these techniques work well in the presence of strong structural plant-model mismatch.

This paper tries to extend the idea of using transient measurements to speed up the convergence to the optimum of the MA methodology, since waiting for the steady state at each RTO iteration is no longer necessary, estimating the process cost and constraint gradients directly by means of a recursive identification method. This identification method does not require any assumption about the type of uncertainty of the model, parametric or structural, or the knowledge of the parameter responsible for the plant- model mismatch, for this reason, it is able to deal with both parametric and structural uncertainties. The key here is to estimate the plant gradients at the steady state from the measurements during the transient of the process. The estimated gradients are then employed to update the value of the modifiers and to solve the modified economic optimization problem every sample time during the transient to drive the operation of the process to a point that satisfies the NCO of the plant faster than traditional MA.

The performance of the proposed method is illustrated through a case study corresponding to an industrial depropanizer distillation column of a petrol refinery. Besides the rigorous model of the column, a simplified one has been used in the implementation of the RTO for MA. One of the main advantages of the use of simplified models in the RTO layer is that the computation time in the optimization is significantly reduced, especially in large scale systems. In addition, a reduced model is easier to understand and easier to update and maintain if necessary. However, using reduced models implies greater plant-model mismatch appearing as structural uncertainty, which makes it even more difficult to achieve a point that satisfies the NCO of the process. The case study also considers the problem of RTO implementation when the decision variables in the RTO do not correspond to the ones of the process due to an intermediate MPC layer not being taken into account in the RTO model, which contributes to the structural uncertainty of the problem.

The paper is structured as follows. Section 2 reviews the MA approaches based on transient measurements. The next section describes the case study used to illustrate the advantages and performance of these methods, the depropanizer column. In section 4, a reduced model of this distillation column has been developed for use in the RTO layer. The RTO economic optimization problem is presented in section 5. Finally, the obtained results are shown in section 6, followed by brief conclusions in section 7.

2. Modifier Adaptation methodology based on transient measurements

MA is normally applied based on static information, as happens in the gradient based method called DMA (Brdys and Tatjewski, 1994) or the gradient free method, NMA (Navia et al., 2015). However, the implementation of these methods is sometimes impractical, especially in processes with a long settling time, as the process needs to reach the steady-state at each RTO execution to estimate the process gradients. Over these long periods of time, the operating conditions or differences between process and model may change and the method may not converge to the real optimum, that is, an operating point that satisfies the NCO of the process.

To overcome this problem, MA approaches that use transient measurements have been developed to estimate process gradients, speeding up the convergence to the plant optimum since the number of required steady states is reduced because the modifiers are updated each sample time during the transient. One technique based on this idea that is valid for parametric uncertainty (François and Bonvin, 2014) is presented first in this section, followed by the proposed method that extends the applicability to plant-model structural mismatch.

Both approaches use transient measurements to estimate the steady-state value of plant gradients and therefore the value of modifiers, thus allowing convergence faster than traditional MA based on static information. For this purpose, RTO is applied at regular time intervals over the transient. At every interval, measurements are used to estimate the modifiers of the optimization problem, the solution of which provides the new set of constant inputs u to be applied to the plant until the next optimization instant, proceeding in such a way that, after the transient, the real plant optimum is achieved. If enough measurements are taken during the transient, and therefore, the estimated gradients have been enough corrected during this time, the true value will be achieved at the steady-state. This implies that an operating point that satisfies the NCO of the process will be reached in only one steady-state. However, this steady-state will be extended over time because the process inputs are manipulated very frequently, every sample time during the transient.

2.1. Modifier Adaptation approach based on NEC to estimate process gradients during the transient

This method was developed by Francois and co-workers in 2014 (François and Bonvin, 2014). The philosophy behind this framework is inspired from Neighbouring-Extremal techniques (NEC), which use transient information for steady-state optimization. The main difference is that the control update is not obtained by computing a control law, but rather by solving a modified optimization problem.

NEC attempts to maintain process optimality in the presence of disturbances by implementing a control law that forces the satisfaction of the process NCO. It relies on linear approximations around the nominal optimum, namely the first-order variations of the NCO, to estimate the variation of the parametric uncertainty ôfi from the changes in the output measurements ôy and the inputs ôu. For this reason, one condition that must be satisfied is ny > np, that is, there are at least as many output measurements y as there are uncertain parameters p.

Consider the following static optimization problem:

min ^(u,P)

5./ gi(u,ß) < 0 i = l,...,ng

where $ is a smooth function that represents the cost depending on the unknown parameters ß and the decision variables u and gi the inequality constraints. The model includes the output equation (7) where h represents the set of equality constraints:

y = h(u, ß) (7)

Here, the variables u and y represent the inputs and outputs at steady state.

From the variational analysis of the first-order variations of the NCO, considering the parametric variations Sß around the nominal values of the parameters, ßnom, a first-order approximation to the cost gradient can be expressed as (8):

Vu^ = ASu + BSp (8)

with A being the Hessian matrix (nu x nu) and B the (nu xnp)matrix given by (9):

A = V2uu^ B = V2up^ (9)

Output variables can be linearized with respect to u and ft:

Sy = QSi + PSp (10)

with the (ny x nu) matrix Q and the (ny x np) matrix P computed by (11): Q = Vuy P = Vpy (11)

Let us assume ny >n^, that is, there are at least as many output measurements as there are uncertain parameters. Using (12), the parametric variations S ft can be inferred from Sy and Su as follows:

Sp = D(Sy - QSi) (12)

where D is a (np x ny ) pseudoinverse of P.

Equation (8) provides a first-order approximation to the cost gradient, which can be estimated from Sy and Su, using (12) to eliminate Sft:

Vu^ = GySy + GuSi (13)

Gy = BD = V2up^(Vpy)+ (14)

Gu = A - BDQ=V2uu^-V2up^(Vpy)+Vuy (15)

RTO modifiers X can be expressed by the following equations, where yp and y are the vector of output measurements from the process (measured during the transient) and the model, up and u the vector of inputs applied to the process and the model, and finally, ynom the vector of model outputs for the nominal solution unom (the solution of RTO without modifiers).

^ = v>P -vu^ = Gy^yp + Guaip -(Gy^y + Guai) = Gy • (¿yp -¿y) + Gu • (<Sip -5i) (16)

where:

¿yp = yp (t) - ynom

¿y = y - ynom (17)

¿Up = Up - Unom

¿1 = u - u

The same estimation procedure is followed to compute the constraint modifiers ji, formulating (13) as the estimation of constraint gradients and computing Gy and Gu as a function of the constraint value instead of the cost value.

Vu gt = Gy^ + Gu£i i = 1,...,ng (18)

Gy = BD = V2uP gt (Vpy)+ (19)

Gu = A - BDQ= V2uugi - V2upgi(Vpy)+Vuy (20)

Yi = VugPi -Vugt = Gy^yp + Gu^ip -(Gy^y + Gu^i) = Gy • (¿yp -¿y) + Gu • (¿up -¿u) (21)

The main advantage of this method is that the process cost gradient is only estimated by using the model offline, through the expression for Gy and Gu, whereas ¿yp is measured from the process. Nevertheless, this method also presents strong limitations, such as:

■ The method works well only if there is parametric uncertainty or little structural plant model mismatch.

■ The user has to know what the uncertain parameters are, although their exact values are unknown.

■ There has to be at least as many available output measurements as uncertain parameters

( ny ^ np ).

2.2. Modifier Adaptation approach based on the direct estimation of the process gradients over the transient.

An alternative approach based on transient measurements considers that process gradients can be estimated from input-output data during the transient using adaptive estimation techniques. The proposed technique uses transient measurements to estimate the plant gradients, and therefore the value of modifiers, thus allowing faster convergence to the steady-state plant optimum than traditional MA.

An approximation for the variation of the process cost function Afa is given by (22), where Auk-1 is the column vector of moves of the input variables (decision variables of the RTO problem) and k represents the RTO execution that in transient methods does not coincide with the steady-state number, as happens in static MA:

A<f>k = Au[-19k-! (22)

The estimator for the variation of the cost function is given by (23) where the hat "A" indicates that the variable is estimated from the process measurements:

Af = AuT_A-1 (23)

where 6 kis the column vector of estimated parameters which contains the required process gradients to compute the modifiers in the first element, i.e., the gradients of the process cost function with respect to the decision variables uk-1.

In this case, we have approximated the variation of the process cost function as a quadratic Taylor polynomial that depends on the input applied until the present sample time and the one applied in the previous sample instant. Other simpler functions can be supposed, for example, a first-order approximation. However, this approximation would be realistic only for linear systems which are rarely met in practice.

Af = k-1 = Vu^u^ + V^-f^u k-1AuT-1 (24)

= Auk1/2Auk-1AuT-1 J (25)

6T-1 =[Vuk-1f V2uk-1uk-1fJ (26)

The gradients contained in 0 k are estimated by employing the recursive extended least squares algorithm (RELS) with forgetting factor a. The concept of forgetting means that older data are gradually discarded in favour of more recent information (Vahidi et al., 2005).

It is based on the difference between the current input uk-1 and the previous one and the difference between the measured A^k and the predicted Ajk change in the cost function for the gradient estimation (Guay, 2014). Then, the parameter estimation update approach is given as follows:

So =11 (27)

Au k-1 = u k-1 - u k - 2

Ajk = <Pk-i0 t-i (29)

ek =Ajk -Aj (30)

(Sk)-1 = 1 (Sk-i )-1 - -r (Sk-1 )-1Pk-1 fi +1 (ph (Sk-1 )-Vk-1 1 vh (S3-1 )-1

. _. ? V—k-1/ ~rk-11 - ^

a a2 v a

0k = 0k-1 + - (Sk-1 )-1Pk-1 (1 + 1 PkT-1 (Sk-1 )-1Pk-1 | (ek) a V a

where E is the covariance matrix of the estimated error whose initial value is E0, ek is the output prediction error, and uk is the solution of the modified RTO problem for each sampling instant during the transient k. On the other hand, A^k = ^k - ^k-1 is the difference between the current process cost function and the cost function measured one sampling time before.

An important drawback of using passive identification with transient measurements is the degree of excitation that the system must have to estimate the gradients accurately. For this reason, a dual constraint is added to the RTO problem (Brdys and Tatjewski, 1994), (Marchetti et al., 2010), (K-1(Sk) > a) that is described by (4) and (5). This constraint represents the dual characteristic of the method: while the rest of the optimization tries to converge to the optimum of the modified model (primal objective), the dual constraint ensures that in the next RTO iteration the system will have enough excitation to estimate the process gradient again (dual objective). However, this technique also identifies second order terms included in 0 k from second order data included in fk whose excitation is not guaranteed by the dual constraint. For this reason the persistent excitation condition (33) based on the theory of system identification must be checked in order to guarantee that all the plant gradients (included the second-order terms) have been properly estimated (Goodwin and Sin, 1984), (Ljung, 1987). This condition, given by (33), indicates that the square matrix fk fkT must be positive definite to guarantee that there is enough process excitation. c (a positive constant) fixes a lower limit for the excitation.

1 k T - ZPkPk ^ CiI

k k=1 (33)

This condition can be satisfied if there is at least as many data sets measured from the process as estimated parameters, that is, k > size[ 0 k ] .

The persistent excitation condition can be checked after the end of the algorithm, once the steady-state of the process has been achieved (as we have made in this paper), or by adding a constraint to the RTO problem at the kh execution, being k = size[ 0k ], to guarantee that there is enough excitation. If this constraint is satisfied at the kh iteration, this condition will be satisfied

during the whole experiment since the more information we get, the better the estimation of the plant gradients.

This constraint would be formulated as shown in (33) where fk is equal to fk-i plus one column vector that includes the decision variables of the RTO, in this way, the optimizer will decide the next inputs to be applied uk forcing that the process is enough excited.

Taking into account that sometimes the inputs of the model are not the same as the inputs of the process, for example, when there is an intermediate control layer that is not considered in the RTO model, or the fact that applying RTO during the transient implies that the input variable may not be able to reach the fixed set points, the plant gradients will have to be estimated with respect to the set points of the process. Therefore, the variable uk that represents the solution of the RTO problem would indicate the set points computed for the input variables. This issue will be shown in the case study described in the section 3.

One advantage of this technique compared to that described in section 2. i is that the computation of second derivatives is not required to estimate the plant gradients. In addition, notice that the formulation does not require any assumption about the type of uncertainty of the model, parametric or structural, or the knowledge of the parameter responsible for the plant-model mismatch. Consequently, the proposed method can be applied to both parametric and structural uncertainty without modifying the RTO model or identifying what the uncertain parameters are.

A summary of the method is given next:

Step 0. Set k = 0 and initialize the vector of estimated parameters 0 0, the value of the estimation error eo and the covariance matrix So whose initial value is given by (27). Then, choose the tuning parameters a (the forgetting factor), the lower limit a for the dual constraint and the lower limit for the persistent excitation c. To start the algorithm it is necessary to apply two operating points, uk-i and Uk-2 to make the first estimation for the process gradients. Only two previous points are necessary to estimate the first gradient, regardless the number of decision variables.

Step 1. At the sample instant k, the process cost function and the constraints gp,i are measured or calculated directly through other measurements. From these measurements and the ones taken in the previous sample time k- i, the variation terms A0p,k and Agpi,k are computed, and the data vector is formed by:

i = i,..'ng (34)

Figure 4 shows how the time domain is divided into several sample periods, which are the time intervals between two consecutive RTO executions. The yellow points indicate the time instants when the measurement was taken and the blue points the instants when the new inputs obtained by solving the RTO are applied.

Step 2. The estimation error for the variation of the process cost function e$,k and for the constraints eg i,k is computed as follows:

, t 1 (35)

egi,k - AgPi,k -Pk-iVg,k-i

1 k - 1 k - ■ 0p, k-i

¿g Pt, k - gPi,k - gPi, k

Au k-i - Uk-i - " U k - 2

Step 3. Estimate the process cost and constraint gradients, Vu<P I and VugPi

luP,i-i ■ uP,i-i

contained in 0<k and 9g k using the RELS algorithm described by (27)-(32).

Step 4. Compute the modifiers Xk, jt,k and et,k of the RTO problem applying the following expressions:

Xk = Vu<P |u -Mu.

k uP,k-1 uk-1

Y k , i =Vu gp, i L -Vu gi |u* i =1,...ng (36)

uP, k-1 uk-1 *

si, k = gp, i (uP, k-1) - gi (uk-1)

Step 5. Solve the RTO modified optimization problem (37) to obtain uk , which is applied to the process until the next RTO execution at k+1.

min = <(u, 9) + XT(u - u^k-1)

st gM,i = gi(u,9) + YTk(u -uP,k-1) + £i,k <0 i = 1,...,ng (37) uL < u < uU K-'(Sk) > «

Step 6. Check if the convergence criterion is satisfied, for instance, u* - u< tol that is

there are no changes in the RTO decision variables higher than a fixed tolerance. Other stop criteria can be considered, for example, there are no changes in the process cost function

higher than a specified value k -$"PkA < tol . If any of these criteria is satisfied, stop the

algorithm, otherwise, k = k + 1, and return to step one when the next sampling time is reached.

2.3. BIBO stability of the closed- loop system

The gradients estimated during the transient using RELS will be equal to the steady-state gradients only when the system converges to a steady state. This fact together with the manipulation of the control set points during the transient may imply that the closed-loop system (controlled plant + RTO) becomes unstable. However, according to the small gain theorem (SGT) for nonlinear systems (Haddad and Chellaboina, 2008), (Khalil, 1996), the following derivations can be made to ensure closed-loop bounded-input, bounded-output (BIBO) stability. Figure 5 shows the basic scheme that concerns the interconnection of two systems, in this case, the controlled plant (input u, output <frP ) and the RTO (input <P , output u). The condition provided by the SGT to ensure BIBO stability of the closed-loop system G x H, is given by (38) for any norm of the involved signals:

GI HI < 1 (38)

In order to fulfil (38) for the induced Lra ^ Lra norm (denoted by the system norm below), the following assumptions are made.

Assumption 1:

According to (37), u is bounded, so Ju^ - a, where a > 0 . Then, it can be concluded that the modified optimization problem described by H is BIBO stable, that is, ||H||i < x.

Assumption 2:

The controlled plant, represented by G, is properly designed, stable and well posed, such that the system is BIBO stable with |<P ^ - %|uiL, for some 1 > % > 0. So, |G||; <i. Then, two possible scenarios appear:

1. The output of the modified RTO problem, i.e. u, is not on the constraints. Then, nothing can be concluded about |H[ , that is, its value may vary between (0,®), so the overall

closed-loop system behaviour may be unstable (<P increases as u does). If this worst case happens, then we will enter to the scenario 2.

2. u is on the constraints. Then, JuL - a independently of how large is <P . Therefore, \\H\[ < i applies for large <P , since - ^lA«,, with 0 < £ < i. In this case, considering Assumption 2, the closed loop system is BIBO stable by the Small Gain Theorem.

Remark 1:

The behaviour of the closed-loop system may lie on the stability limit in the worst case ( G|| • |H| -1), i.e., the system may be continuously jumping between scenario i and 2, but, without becoming unstable in the BIBO sense.

Remark 2:

Note that the assumption 2 about |G|| < i is strong, but it can be relaxed indeed to be |G|| <

still fulfils (38). Furthermore, the filter parameters and the sampling time of the MA algorithm can act as tuning parameters to smooth the system response, thus reducing |H|.

3. Case study

The proposed algorithm for gradient estimation has been applied in an industrial case study and the results obtained have been compared with those obtained with other MA approaches.

3.1 Process description

The case study considers the optimal operation of a depropanizer distillation column of a petrol refinery which isolates propane (C3H8) from a mixture containing butane (C4H10) and other components. In general terms, continuous distillation columns use the variation of temperature and pressure conditions along the height of the column to get more volatile component at the top of the column, propane in this case, and less volatile component at the bottom of the column, among them, butane. The control objective for the distillation column is to maintain the composition of propane in distillate and bottom streams within their desired specifications.

A set of basic controllers maintains liquid levels in the accumulator and the bottom of the column at their desired values, there are also steam and reflux flow controllers and another one

for the head pressure. In addition, an MPC controller, which manipulates the references of steam and reflux flow controllers, tries to maintain head and bottom temperatures in the column close to its set points, as an indirect way of controlling distillate and bottom composition. Later on, an RTO layer is added to improve the economic performance of the column while maintaining the composition of propane in distillate and bottom streams within the desired specification. This control structure is shown in Figure 6.

The decision variables of the RTO layer are the head and bottom temperatures, which are passed as set points to the MPC controller (a DMC one).

The simulated depropanizer column is based on an industrial example located in the Repsol YPF Tarragona Refinery (Spain), made up of a total condenser, a partial reboiler, 37 equilibrium stages, and operating at 1.57 106 Pa. The feed mixture enters at stage 19 at a flow rate of 468 kmol/h and 330.42 K. The composition of the feed in the example is 45.55 mol% propane, 44.67 mol% butane and 9.77 mol% ethane, while the tray efficiency is 60 %. A rigorous first principles model was developed to represent the process whose main equations are described below.

3.2 Material and energy balances

The column mass balances for each tray are given by equations (39)-(43), with reference to Figure 7, where dmol/dt is the molar accumulation rate (kmol/h), i is the tray number, l, v, f stand for liquid, vapour and feed molar flows respectively (kmol/h), lref is the molar reflux flow (kmol/h), and b and d are the molar flow of bottom and distillate streams, also in kmol/h. (39) presents the overall mass balance around the ith tray, whereas (40), (41), (42) and (43) show the mass balances in the feed tray, the top of the column (tray n), the bottom (tray 1) and the top accumulator. In (43), lacum is the flow that comes from the condenser and loverflow is the excess liquid that overflows from the accumulator (Coulson and Richardson, 2002) and (Rodríguez-Blanco et al., 2015):

dmol, dt dmol

- = l/+i + V-1 " h - V

i = 2,...n -1

n _ feed dt

f + ln feed+1 + Vn_feed-1 ln feed Vn_feed

dmoln dt dmol dt dmol

= lref + Vn-1 - ln - Vn

1 = I2 + V - b

~acum = l - l - d - l

acum ref overflow

From the molar flow rates, the mass flow rates (kg/h) of the different streams, reflux R, distillate D, bottom B and feed F, can be obtained by applying the expressions (44)-(47), where Pm (kg/kmol) indicates the molecular weight of each stream:

R = lrefPm,lref (44)

D = dPmM (45)

B = bPmb (46)

F = fPm,f (47)

Component mass balances are expressed by (48)-(50), where Xj,i and yj,i are the composition of component j (butane, propane and ethane) in the liquid and vapour streams through the ith tray (°/1mol):

Ji = l,+ixj,,+i + vi-iyj,i-i " l,xj,, -v,yj,,

dx j,n feed dt

2 xj, = 1

i = 2,...n -1

fzj + ln feed+1xj,n feed+1 + vn_feed-1yj,n_feed-1 ln _feedxj,n_feed Vn _feedyj ,n feed

2 yj,i =1

i = 1,...,n

Due to their speed of response, energy balances around the ith tray have been modelled as steady-state equations and are expressed by (51)-(53), where Hv is the vapour enthalpy (kJ/kg), hi is the liquid enthalpy (kJ/kg), and hj is the specific enthalpy of each component (kJ/kg) that depends on the tray temperature T (K):

Hv,,V, = hl,,+1ln+1 + Hv,,- 1V,-1 -hlßi

i = 2,... n -1

v,n_feedVn_feed = hff + hl,n_feed +1ln_feed +1 + Hv,n_feed- 1Vn_feed-1 hl,n_feedln

Hv,, =2 yjihj(Tl)

, =2 xj,ihj(Ti)

i = 1, ...n

i = 1, ...n

The key temperatures in the operation of distillation columns are the head temperature (Thead) and the bottom temperature (Tbottom). These temperatures correspond to the sensitive trays of the column where sensors are located. In this case, the sensitive trays are the number 34 for Thead and 4 for Tbottom.

3.3 Vapour-liquid equilibrium:

The concentration of vapour in contact with liquid at equilibrium yeqj, i in each tray is expressed by Raoult's law (55), where the equilibrium constant Kj,i is given by the ratio between the component j vapour pressure Psatj(Pa) and the total pressure in the ith tray Pi (Pa). The parameter Ef indicates the tray efficiency:

yeqj, = ef kj ,xj i

_ ' sat,j j p

i = 1,...n i = 1,...n

3.4 Energy balance of the reboiler:

A static energy balance is considered in the reboiler, where Qb (kJ/h) is the heat flow added to the reboiler, and Ah is the latent heat of vaporization (kJ/kg) that depends on the reboiler pressure Preboiier (bar), which has been considered as a constant value.

Hv,1V1 = Qb-hl,1b Qb = S ^h(Preboiler)

3.5 Energy balance of the condenser:

The steady-state energy balance corresponding to the condenser is expressed by the equation (59), where Qc is the heat removed in the condenser (kJ/h):

Qc - Vn (Hv ,n (Thead) " hl ,n (Thead))

The nonlinear dynamic model described above has been simulated using the software EcosimPro (EA Int, 2013), a modern object oriented simulation environment. The model is composed by 2056 equations (129 differential equations and 1927 algebraic ones) and 2152 variables (1976 explicit, 129 derivative, 40 algebraic and 7 boundary variables). This rigorous model was also used to identify the linear dynamic model of the DMC controller, which uses a step response model that mainly relates the head and bottom temperatures with the set points of the reflux and steam controllers (see Figure 6).

4. RTO model

While the previous dynamic model will be used to represent the real process, the RTO layer uses a steady-state model to make decisions. A rigorous steady- state model may be developed directly by removing derivatives from the dynamic column described above, but the resulting model has a large number of variables, equations and algebraic loops that make the initialization and convergence of the optimization difficult. In addition, this complexity makes the optimization problem solved in the RTO run slowly.

For this reason, a simplified model has been developed to be used inside the RTO layer. The reduced model takes into account the fact that distillation columns can be divided into three sections: namely, the rectifying or enriching section, the stripping section and the feeding tray, as Figure 8 shows. The feeding tray is the tray of the distillation column where the feed is introduced. The enriching section consists of the trays above the feed tray and is so named because, in this section, rising vapours are enriched in more volatile component. Since liquid is stripped off the more volatile component in the section below the feeding tray, it is called the stripping section.

In order to obtain a fairly small reduced model, the rigorous column model has been simplified by considering only three trays and a global efficiency. The equations that describe the steady-state simplified model used by the RTO layer are described below.

The overall mass balance of a distillation column is given by (60), where f d and b are the feed, distillate and bottom molar flow rate (kmol/h); whereas the mass balances on the enriching and stripping sections are given by (61) and (62) respectively, where v and Vs are the vapour molar flow rate in the enriching and stripping sections (kmol/h), lS is the liquid molar flow rate in the stripping section (kmol/h), and lref is the reflux molar flow rate (kmol/h).

f - d + b

V - lref + d

lS - vS + b

(60) (61) (62)

The component mass balances are expressed by (63)-(65), where xSj and ys . are the composition of component j (butane, propane and ethane) in the liquid and vapour streams through the stripping section, and Xjandy,- are the same compositions for the enriching section. These equations are applied for two components, in this case butane and ethane, which are those present in lower concentrations; while the composition of propane is obtained by (66). The feed enters the column with a composition of Xj,njeed and distillate and bottom streams leave the column with a composition equal to xDj and xBj respectively.

f xj,n _feed - d XD,j + b XB,j

lSxS - vSyS + bx

yj- x S ,x,

-S,xS-S,yS, -S,

Xj -S jyj -S j xH ,■ -1

J j xj,n_feed - S j xj - S j yj - S j xj - S j yj - S j xB,j

The energy balances are expressed by (67)-(69), where Qc is the heat removed in the condenser (kJ/h), Qb is the heat flow added to the reboiler (kJ/h), Ah is the latent heat of vaporization (kJ/kg) that depends on the reboiler pressure Preboiier(bar) , Hs and hs are the enthalpies for liquid and vapour streams in the stripping section (kJ/kmol), hB is the enthalpy of the bottom stream (kJ/kmol), hDis the enthalpy of the distillate stream (kJ/kmol), and HD is the enthalpy input vapour stream to the condenser (kJ/kmol). Every value of enthalpy is a function of head, Thead (°C), and bottom, Tbottom (°C), temperatures.

Qc - v(HD (Thead) " hD (Thead))

Qb - S 3l{Preboiler)

lSh-S (bottom) + Qb - v HS (Tbottom) + b hB (Tbottom)

(68) (69)

Along the column, there is a pressure drop that is a function of the vapour stream passing through the column and the number of trays ntray present in the real column. As Vand Vs have approximately the same value, both of them could be used to compute this pressure drop and, therefore, the value of the bottom pressure Pbottom- The head pressure Phead is mantained by a controller, while c is a constant value to compute the pressure drop through the column.

bottom head

+ (vS/c)2 • n,

The concentration of vapour and liquid streams leaving each zone are at equilibrium, which is expressed by Raoult's law (71), which depends on the equilibrium constants K. The parameter Ef indicates the tray efficiency, which varies from 0 to 1, which is the same as saying that it varies from 0 to 100%. At the reboiler, the efficiency is equal to 1, since it is considered an ideal tray:

yS K (PbottomtTbottom)xS

yj - Ef (K (PheaJhead >j " yj ) + yS

Due to the absence of energy balances for each tray of the column, head and bottom temperatures have been obtained from experimental data using the ALAMO (Automated Learning of Algebraic Models for Optimization) software for model building, (Cozad et al., 2014) and (Cozad et al., 2015). ALAMO generates algebraic models of simulations or experiments, in this case, models for head and bottom temperatures have been generated as functions of the inputs of the model, that is, the steam flow to the reboiler S (kg/h) and the reflux flow R (kg/h). These expressions are obtained from the simulation of the rigorous dynamic model of the distillation column, where f represents a second order polynomial function. For this task, 82 experiments on different operating points have been made to obtain the functions that better adjust all the data taken from the simulation.

Thead - f (S'R) Tbottom - f (S>R)

Obtaining R from the following equation, where Pm,R (kg/kmol) is the molecular weight of lre/

R = lrefPm,R (75)

It is assumed that all the light key component, ethane, will be removed completely from the column in the distillate stream.

xB (C2H6) = 0 (76)

The main differences between the rigorous model that can be obtained directly from the dynamic one and the simplified steady-state model are presented in Tables 3 and 4:

As shown in Table 3, the model size has been reduced by 28 times. The execution of one RTO problem is 32 times faster using the reduced model, which allows the use of RTO in real applications where decisions must be implemented within a time scale of a few seconds. In this example, the use of the rigorous model does not imply a problem with respect to the optimization time, but it is important to note that this issue could happen by using large scale models, where the optimization time could be prohibitive and the use of reduced models is the only way to deal with them.

In order to test the reduced model, several changes in the input variables, reflux and steam flow rates, have been applied to the reduced and full models in parallel, as shown in Figure 9, starting from a steady-state to check that the behaviour of the output variables is similar for both models. Some of the output variables have been represented in Figures 10, 11, 12, 13 and 14, where it can be observed that the achieved steady-states are not the same, due to the strong structural uncertainty.

5. RTO problem formulation

The RTO problem is based on the static reduced model described in the previous section that can be considered as a grey-box model, since it combines a partial theoretical structure with experimental data to complete the model. This model reduction implies a strong structural plant-model mismatch; in addition, the real process is controlled by an MPC layer which is not considered in the steady-state model, so the inputs of the real process, DMC set points for bottom and head temperatures, are not the same as the model inputs, that is, steam and reflux flows.

Another source of uncertainty is due to the identification of the linear dynamic models that are required by the DMC controller. This identification is carried out around a certain operating point, so these models could be incorrect if the process is operated away from this point, resulting in more uncertainty the further the decision variables are moved from the identification point. The presence of this strong uncertainty implies that the optimal solution of the RTO model will not be the same as the process optimum, making it necessary to apply MA methodology to achieve the optimal operating point.

The economic objective of RTO is to maximize the operation profit considering the benefits obtained from producing distillate D (kg/h) with a given purity specification, and minimizing the steam consumption S (kg/h) for a given feed (Porru et al., 2015). Prices Ps and Pd have been fixed for both streams, 20 €/ton for the steam consumed and 80 €/ton for propane in distillate, assuming that it is above the composition specification. Below this target, the price decreases as a function of the composition following a sigmoid behaviour. The objective function is represented by the value of $ (€/h) calculated by (77) and the constraint over the composition of propane in the distillate stream g is expressed by (78), whereas the RTO problem is formulated as:

max < = PD (xD (C3H8 )D) - PsS

Stationary reduced mode/ of the process g = Xd (C3H8 ) > 0.80

Lower and upper limits are fixed for the model inputs, that is, steam and reflux set points.

4000 < S < 6000 6000 < R < 11000

The value of the temperatures (°C) Thead and Tbottom shall fall within the following intervals:

40 < Thead < 60 80 < Tbottom < 110

Other constraints have been taken into account; however, for simplicity, they will not be modified following MA methodology. These constraints are the maximum loss of pressure AP to prevent the flooding of the column and the specification for propane composition in the bottom stream:

AP < 0.3

xB(C3H8 ) < 0.02

As mentioned before, in this example, the inputs of the model are not the same as the inputs of the process due to the presence of the MPC layer. The RTO model does not consider the DMC controller working in closed loop, so the set of model inputs (S, R) does not include the closed-loop plant degrees of freedom (Tbottom, Thead), as shown in Figure 15.

The modified problem is given by (82) and (83). The subscript "P" indicates that the variable is measured from the process and the subscript "k-1" is the measurement taken in the previous steady state:

max<M = PD (xD (C3H8))D - PsS + ^k (c(u) - c(uP,k-1))

gM = -xd(C3H8) + 0.80 + yT(c(u)-c^-O + ek <0 Stationary mode/ of the process

4000 < S < 6000 6000 < R < 11000

40 < Thead < 60

80 < Tbottom < 110

AP < 0.3

xB (C3H8) < 0.02

(82) (83)

The number of modifiers required to adapt the RTO problem nk is given by (84), where nu is the number of decision variables which are head and bottom temperatures (Thead, Tbottom), and ng is the number of constraints (distillate composition); so, in our problem, nk = 5.

nk = ng + nu (ng + 1)

Modifiers X, y and e are given by (85)-(87) and represent the difference between experimental and model gradients. Although several techniques to compute process gradients have been presented in the previous sections, another key point is the calculation of model gradients, but in this case they cannot be computed directly, since Tbottom and Thead are not considered as independent variables of the RTO model, so the modifiers are computed as follows (Costello et al., 2013):

! T 1 k = VA lc(up,k _1) / * c( uk _ l) = V^l

T Y k = Vcgp c(u ) _Vcg c(uP,k _1) / * c(u*k_ 1) =Vc gP

sk = gp (c(up,k -1)) _ g(uk - 1)

c(u ) "Vu^(VuC)+ .

c(up,k-1) uk _1

lc(uP,k _1)

_Vu g(Vuc)+

(86) (87)

where (+) indicates the Moore-Penrose pseudo-inverse.

As usual in MA methodology, the modifiers are applied filtered by using a first-order exponential filter where Ki, Kr and K are diagonal matrices with eigenvalues in the interval (0, 1], (88)-(90) to smooth changes in the decision variables and to ensure convergence. I k = (I _ KJX k_i + K ^ k (88)

Y k = (I _ K y )y k_i + K yy k (89)

8k = (I _ K£ )Ek _1 + K£8k (90)

As has been shown in the figures of the previous section, the reduced model is not perfect so the model optimum is not the same as the real optimum (see Table 5).

To solve the nominal RTO problem without modifiers implies a decrease of 8.77 % in the profit obtained. To solve this issue, MA techniques can be implemented like the ones explained below. The NLP optimization problems have been solved using a sequential approach with a sequential quadratic programming (SQP) algorithm implemented in SNOPT library (Gill et al., 2008) and executed in EcosimPro software (EA Int, 2013).

6. Results

The rigorous dynamic model (described in Section 3.1) has been simulated in EcosimPro software to be used as a substitute for the real process in order to study the performance of the different MA approaches. In particular, we aim to test the proposed method based on the direct estimation of the process gradients, comparing it with standard DMA and NMA formulations, as well as with the one based on NEC and transient measurements.

6.1 Static Modifier Adaptation methodology

The DMA (Brdys and Tatjewski, 1994) , (Marchetti et al., 2010) and NMA (Navia et al., 2015) methods have been implemented on the depropanizer column solving the RTO optimization problem every 6 hours, the time required for the process to achieve a new steady-state. In both approaches, DMA and NMA, the first two iterations were used to estimate the process gradients by finite differences; in the second approach this is necessary to provide a good starting point for the upper optimization layer. The third operating point corresponds to the nominal solution, that is, the result obtained by solving the RTO problem without modifiers.

The DMA methodology requires the addition of a new constraint on the optimization problem (S > a), which takes into account the grade of excitation of the process so as to ensure sufficient information in the measurements to guarantee gradient accuracy. This is equivalent to choosing the parameter a (lower limit of K-1(Sk)), whose selected value in this case is 0.10.

Figure 16 shows the evolution of the cost function during the process operation, whereas Figures 17 and 18 show the bottom and head temperatures. The evolution of the composition of the distillate is shown in Figure 19, where it can be observed that the achieved operating point corresponds to an active constraint represented by the dotted red line. Figures 20, 21 and 22 represent the steam, reflux and distillate flows, respectively.

Fixing a tolerance band of 0.5% with respect to the optimal value of the cost function, the graphs (in particular, Figure 16) show that the optimum of the process is achieved after approximately 60 hours, using both DMA and NMA approaches. For the DMA approach, the number of steady-states needed was 10, the first two to initialize the computation of the process gradients (n = 2) and 8 related to RTO solutions. NMA also requires 10 steady-states, first, two initial steady-states (nu = 2) to estimate plant gradients, then five in order to construct the first simplex corresponding to the number of first- order modifiers plus 1, and finally 3 RTO executions.

DMA only requires one tuning parameter a, the bigger this parameter the greater the excitation of the system, driving the process to the optimum following a less optimal path. However, the Nelder Mead algorithm used in NMA to estimate the modifiers has many tuning parameters which strongly affect the speed convergence and the path followed to achieve the real optimum. It is a difficult task to tune these parameters, since the effect of each one is not really known, so a bad choice could mean reaching the optimum slowly, implying many steady-states. In addition, the estimation of the initial modifiers also affects the behaviour of the method.

An overshoot is observed in Figures 20 and 21, this is due to the saturation of the valves that manipulate the reflux flow and the steam flow. This happens because the reference given by the RTO layer to the DMC cannot be reached by the process (see Figures 17 and 18). This solution corresponds to the nominal one, that is, the resultant from solving the RTO problem without modifiers. In spite of this unfeasible first solution, MA is able to bring the process to the real optimum, giving feasible solutions in the following iterations.

One important disadvantage of this kind of method is that they can work well with few decision variables but, when the number of inputs increases, the number of modifiers also grows, involving a bigger number of steady-states to converge to the optimal solution, which supposes a problem, especially for processes with long settling times like in the depropanizer column.

As shown in Figures 16 to 22, the convergence to the optimum using static MA techniques takes two and a half days, which is a problem if the methodology is implemented in real problems. To deal with this issue, transient measurements have been incorporated in the MA methodology to speed up the convergence to the real optimum, as described previously.

6.2 Modifier Adaptation methodology using transient measurements

In this section, the results obtained by applying the MA techniques based on transient measurements to estimate process gradients are presented. Measurements from the process are taken every hour during the transient and the RTO problem is updated at the same rate. This sample time corresponds to one sixth of the settling time of the process.

The first approach described in section 2, based on NEC to estimate process gradients, needs to know where the parametric uncertainty between process and model is located. The use of the

reduced steady-state model implies a strong structural plant model mismatch; for this reason, this mismatch has been adjusted to the rigorous one that represents the real process by adding some offsets that take into account the structural mismatch. In this case, the steady-states achieved by the dynamic process can be matched with the ones of the model by changing three uncertain parameters that have been considered as parametric uncertainty: two constants (chead and cbottom) added in the equations used to compute head and bottom temperatures, (73) and (74) with a nominal value of 0, being the easiest way to deal with it, but not the most efficient, and the global efficiency (Ef) of the column (72), whose nominal value is equal to 0.98, a value which has been adjusted from simulation results. However, these parameters are time variant, so the implementation of this technique is possible and works effectively within a certain operating range. This task, to translate the structural uncertainty into parametric, could be a difficult task in some cases, so it is a clear disadvantage of this technique:

Thead = f (S,R) + Chead (91)

Tbottom = f (S,R) + cbottom (92)

Therefore, in Algorithm 1, the vector of uncertain parameters is ft = [chead, Cbottom, Ef], the decision variables (the model inputs) are u = [R, S] and the vector of output variables y = [B, D, xs(C^Hs)]. It is possible to choose a different set of output variables, but the selected ones are the most sensitive to the uncertain parameters. The results using this method are shown in the figures as TMA.

In the same way, this section presents the results of the proposed method based on estimating process gradients using a recursive parameter identification algorithm, to be precise, Recursive Extended Least Squares (RELS), with a forgetting factor a equal to 0.99. The lower bound for the degree of excitation, C1, of the system is fixed to 0.05. The results obtained using this approach are shown in the figures as EstMA. In this example the size of the data vector and the estimated parameters, given by (93) and (94), is five, so, after the fifth iteration there is enough collected data to check if the persistent excitation condition is satisfied.

Vk = \^Tbottomk AThead,k ATbottomk AThead,k ATbottomkAThead,k J (93)

^Tbommdt$ ^TheadTbotbm,kThbotbmk,$ 2 ^ Thead,kTheadk,$ ^ Tbottom,kTheadk

The process starts from a certain operating point and, after one hour, the nominal solution is applied to the process, that is, the RTO solution without modifiers. Thus, both methodologies have the same transient information when they start working. In this example, modifiers have been applied without filtering.

Figure 23 presents the evolution of the process cost function, whereas Figures 24 and 25 show the changes in the reflux and steam flow rates obtained by applying the described methodologies. The evolution of the composition of the distillate is shown in Figure 26, observing that the achieved operating point corresponds to an active constraint indicated by the dotted red line. Figures 27, 28 and 29 represent the steam, reflux and distillate flows respectively.

Observing the previous graphs (in particular Figure 23), one can see that the convergence of both methods to a region near the optimum (a tolerance band of 0.50% with respect to the process cost function has been considered) is faster than in the previous section, around 8 hours in both cases, which is a clear advantage that the approaches based on transient information present. However, notice that, due to the presence of structural uncertainty, there is an offset

with respect to the real process optimum in the solution provided by the method based on NEC. By contrast, when estimating the process gradients directly with the proposed recursive identification algorithm, the real optimum is achieved after one steady-state. Nevertheless, it must be pointed out that the settling time of the process has grown, since the DMC set-points were being changed continuously (every hour) during the transient.

By applying the RELS algorithm to update the modifiers during the transient, the optimum operating point is achieved after 8 hours, which involves 8 RTO solutions (one per hour), reducing by approximately 8 times the time required to achieve the optimum as compared with traditional static MA techniques, such as DMA or NMA. At iteration five, when the identification method has enough information, the persistent excitation condition has been checked observing that this condition is satisfied and therefore, concluding that the process gradients have been properly estimated.

6.3 Evaluation of the dual constraint in MA based on the direct estimation of the process gradients over the transient.

A dual constraint has been added to the MA technique based on the estimation of the process gradients using RELS to ensure that in the next RTO iteration, the system will have enough excitation to estimate the process gradients again. Figure 30 shows the evolution of the dual constraint in the implementation described in the section before. In this case the lower limit for the dual constraint is 0.02.

As view of the results, it can be concluded that the system has enough excitation during the whole operation. In some iterations the dual constraint becomes active. This is because, the RTO solution does not provide enough excitation, so, the dual constraint acts to guarantee that the gradients will be properly estimated in the next RTO.

A second example is shown below to proof the effect of the dual constraint in a system that not is adequately excited. A very small constant filter has been considered, Kx = Kr = K = 0.15, to force that situation, that is, the changes in the input variables proposed by the RTO do not provide enough excitation to estimate the gradients accurately. In this situation, if excitation is forced by the dual constraint, an operating point that satisfies the NCO of the process can be achieved. In the following figures, RTO filtered, indicates the results obtained without adding the dual constraint and RTO dual, presents the results adding the dual constraint to ensure excitation. Figure 31 presents the evolution of the process cost function, whereas, Figure 32 and 33 show the evolution of head and bottom temperatures respectively. Figure 34 shows the value of the dual constraint whose lower limit has been fixed in 0.10.

The previous results show that the addition of the dual constraint, with an appropriate lower limit, provides enough excitation to the system. In this way, the transient MA method presented in this paper will be able to achieve an operating point that satisfies the NCO of the process, even if the RTO does not induce important changes in the decision variables by itself.

It is important to note that the convergence rate is strongly affected by the constant of the filters. In this last example, the optimum is achieved after 60 hours, this is due to the modifiers have been strongly filtered, slowing down the convergence of the method.

7. Conclusions

In this paper, a new method to speed up the convergence of RTO-MA to the real plant optimum has been proposed. It is based on transient information, obtaining process gradients directly

from truncated Taylor expansions of the process cost and gradients combined with adaptive filtering estimation techniques. The method has been tested in a realistic case study corresponding to a depropanizer distillation column of a petrol refinery and the results obtained show that it is possible to effectively perform the optimization of the distillation column in the presence of structural plant-model mismatch and the presence of an MPC layer not considered explicitly in the RTO model. This reduces by a factor of 8 the time required to achieve the process optimum as compared with traditional static MA techniques, such as DMA or NMA.

Nevertheless, there are some open problems that deserve further research. For instance, the size of the problem with respect to the number of decision variables and the number of constraints is still an open issue, because it increases the number of experimental gradients that must be estimated and the number of modifiers that are needed.

This paper also shows the clear advantages of using simplified steady-state models in the RTO layer; the solution time is significantly reduced during the optimization of the depropanizer column operation (approximately 28 times in this case), which allows the implementation of RTO in real processes. In addition, the reduced model is easier to understand and easier to update and maintain if necessary, without requiring an extensive knowledge of the process.

Therefore, the combination of modifier adaptation methodology, the least squares algorithm to use transient measurements to compute process gradients, and the use of simplified models in the RTO layer is really a powerful approach to achieve the optimal operating point of processes, reducing the convergence time.

8. Acknowledgments

The authors wish to express their gratitude to project DPI2015-70975-P of the Spanish MINECO for the financial support for this study under FPI grant BES-2013-062737.

9. References

Brdys, M. and Tatjewski, P., 1994. An algorithm for steady-state optimizing dual control of

uncertain plants. In 1st IFAC Workshop on New Trends in Design of Control Systems, pp. 249-254, Smolenice, Slovakia.

Costello, S., François, G., Bonvin, D., 2013. Real-time optimization when the plant and the

model have different inputs, in: Proc. IFAC Symp. DYCOPS 2013, Mumbai, pp. 39-44.

Coulson, J.M., Richardson, J.F., 2002. Coulson and Richardson's Chemical Engineering, Volume 2, 5th edition, Butterworth-Heinemann.

Cozad, A., Sahinidis, N.V., Miller, D.C., 2014. Learning surrogate models for simulation- based optimization. AIChE Journal, 60, 2211- 2227.

Cozad, A., Sahinidis, N.V., Miller, D.C., 2015. A combined first- principles and data-driven approach to model building. Computers and Chemical Engineering, 73, 116- 127.

EA Int, 2013. EcosimPro 5.2. User Manual (http://www.ecosimpro.com).

François, G., Bonvin, D., 2014. Use of Transient Measurements for the Optimization of Steady-State Performance via Modifier Adaptation. Industrial & Engineering Chemistry Research, 53 (13), pp 5148-5159.

Gao, W., Engell, S., 2005. Iterative set-point optimization of batch chromatography. Computers & Chemical Engineering, 29, 1401-1409.

Gill, P. E., Murray, W., Saunders, M. A., 2008. User's Guide for SNOPT Version 7: Software for Large-Scale Nonlinear Programming.

Goodwin, G. C., Sin, K. S., 1984. Adaptive filtering prediction and control. Prentice-Hall, Information and system sciences series.

Guay, M., 2014. A time-varying extremum-seeking control approach for discrete-time systems. Journal of Process Control 24, pp. 98-112.

Haddad, W. M. and Chellaboina, V. Nonlinear Dynamical Systems and Control: A Lyapunov-based approach. Princeton Universtity Press, ISBN-9780691133294, 2008.

Khalil, H.K. Nonlinear systems. Prentice-Hall, ISBN-0130673897, 1996.

Ljung, L., 1987. System Identification. Theory for the user. Prentice-Hall, Information and system sciences series.

Marchetti, A., Chachuat, B., Bonvin, D., 2010. A dual modifier-adaptation approach for realtime optimization. Journal of Process Control, 20, 1027-1037.

Navia, D., Briceño, L., Gutiérrez, G., de Prada, C., 2015. Modifier- adaptation methodology for real-time optimization reformulated as a nested optimization problem. Ind. Eng. Chem. Res 2015; 54:12054-71.

Naysmith, M.R., Douglas, P.L., 1995. Review of Real Time Optimization in the Chemical Process Industries. Dev. Chem. Eng. Mineral Process. 3: 67-87. doi: 10.1002/apj.5500030202.

Porru, M., Baratti, R., Alvarez, J., 2015 Energy saving through control in an industrial

multicomponent distillation column. 9th International IFAC-ADCHEM (7-10): 11391144.

Roberts, P.D., 1979. An Algorithm for Steady-State System Optimization and Parameter-Estimation. International Journal of Systems Science, 10, 719-734.

Rodríguez-Blanco, T., Sarabia, D., Navia, D., de Prada, C., 2015. Modifier- Adaptation

methodology for RTO applied to distillation columns. ADCHEM 2015, Whistler, British Columbia, Canada.

Tatjewski, P., 2002. Iterative Optimizing Set-Point Control- The Basic Principle Redesigned. 15th Triennial World Congress, IFAC. Barcelona, Spain.

Vahidi, A., Stefanopoulou, A., Peng, H., 2005. Recursive least squares with forgetting for online estimation of vehicle mass and road grade: Theory and experiments. Vehicle System Dynamics, 2005: 43, pp.31-55.

Yip, W.S., Marlin, T.E., 2004. The effect of model fidelity on real-time optimization performance. Computers & Chemical Engineering, 28, 267-280.

Zhang, Y., Roberts, P., 1990. On-line steady-state optimisation of nonlinear constrained

processes with slow dynamics. Transactions of the Institute of Measurement and Control 12 (5), 251 - 261.

$ P,k-l gp,i k-l

YES STOP

Figure 1. General formulation of Modifier Adaptation methodology $ P,k-1 gP i,k-1

Upper optimization layer (Nelder Mead Algortihm)

Filtering

^k, Yk,i ek,i

Nested Modified Optimization

k = k +l

YES STOP

Figure 2. Schematic NMA methodology

< Process Gradient Estimation (RELS)

Modifiers computation

\k, Vi,k £¡,k r

Modified Economic Optimization

MATHEMATICAL MODEL (Taylor Polynomial)

Figure 3. Schematic of MA based on the direct estimation of process gradients over the transient

Unom 1 Uk-2 Ufr-l u* -

Au/f-i

Qp, QP k-2 k- 1 <t> P,k ____<

°,nom ——- j r -—• 1 <k, ___(f)p,t k-2 ' -1 .

V,nam ,-2 C-l c k+1

Figure 4. Sample periods for the estimation ofplant gradients

u N CONTROLLED PLANT (G)

RTO (#) s

Figure 5. Closed-Loop Setup

Figure 6. Control structure of a depropanizer distillation column

Vi+1 ' yj.i+1

Vi > yji li+1 K xj,i+1 \ I

Vi-i y yj.i-i ir~ xj,i \ \ t

1,-1 xj,i-1 \ f

Tray i+\

Tray i

Tray i-1

Figure 7. Schematic of material balances for each tray

Condenser

y = XD

Vs, yS

(g) Qc

L, x LS, xS

Enriching section Stripping section

Reboiler

Figure 8. Schematic diagram of a simplified model of a continuous distillation column

~ 9000 m

— 8000 §

^ 6000 5000

0 5 10

Time (h)

Figure 9. Changes in the reflux and steam flows(R and S)

Steam Flow • Reflux Flow

0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5

10 Time (h)

• Rigorous Model

• Reduced model

Figure 10. Evolution of the composition ofpropane in the distillate (xd(C3Hs))

13000 12000 11000 10000 9000 8000 7000 6000 5000 4000

Figure 11. Evolution of the bottom flow rate (B)

■M ■M

10 Time (h)

h) 14000

Bfl 13000

S O 12000

e 11000

= 10000

Q 9000

Time (h)

Figure 12. Evolution of the distillate flow rate (D)

• Rigorous Model

• Reduced model

Rigorous Model Reduced model

■M ■M

105 104.5 104 103.5 103 102.5 102 101.5

Time (h)

Rigorous Model Reduced Model

Figure 13. Evolution of the bottom temperature (Tbottom)

u ( 55

3 ■M 50

Q. m 45

■a (U 40

Time (h)

Figure 14. Evolution of the head temperature (Thead)

• Rigorous Model

• Reduced Model

u = [S, R]7

* MODEL

u) = \Tbottom., TheadV

c(u*) = \Tbottom

Figure 15. Comparison between model and plant inputs and outputs

Time (h)

Figure 16. Evolution of the process cost function using static MA techniques

Tbottom NMA SP Tbottom NMA Model Optimum Real Optimum Tbottom DMA SP Tbottom DMA

Time (h)

Figure 17. Evolution of the bottom temperature (Tbottom) using static MA techniques

Figure 18. Evolution of the head temperature (Thead) using static MA techniques

0.835 0.83 0.825 _ 0.82 ~ 0.815 jç 0.81 * 0.805 0.8 0.795 0.79

0 20 40 60

Time (h)

Figure 19. Evolution of the composition ofpropane in the distillate (xD(C3H8))using static MA techniques

Active Constraint

E 5900

Steam Flow DMA (kg/h) ■ Steam Flow NMA (kg/h) Real Optimum

Time (h)

Figure 20. Evolution of the steam flow (S) using static MA techniques

!c 13500 m

* 12500 o 11500

J 10500 St 9500 8500 7500

Reflux Flow DMA (kg/h) Real Optimum Reflux Flow NMA (kg/h)

Time (h)

Figure 21. Evolution of the reflux flow (R) using static MA techniques

11900 11400 10900 10400 9900 9400 8900

Distillate Flow DMA (kg/h) Distillate Flow NMA (kg/h) Real Optimum

Time (h)

Figure 22. Evolution of the distillate flow (D) using static MA techniques

850 800 750

Si 700

650 600

Time (h)

TMA EstMA

Model Optimum Real Optimum

Figure 23. Evolution of the process cost function ($p) using MA based on transient measurements

Time (h)

Thead EstMA SP Thead EstMA Model Optimum Real Optimum Thead TMA SP Thead TMA

Figure 24. Evolution of the head temperature (Thead) using MA based on transient measurements

104 102 100

94 92 90

Tbottom EstMA SP Tbottom EstMA Model Optimum Real Optimum Tbottom TMA

Time (h)

Figure 25. Evolution of the bottom temperature (Tbottom) using MA based on transient measurements

TMA EstMA

Time (h)

Figure 26. Evolution of the composition ofpropane in the distillate (xofCHsJJusing MA based on transient measurements

BO 5700

s o 5500

m m 5300

■ Steam Flow TMA (kg/h) Steam Flow EstMA (kg/h) Real Optimum

Time (h)

Figure 27. Evolution of the steam flow using MA (S) based on transient measurements

BO 10900

o 9900

M- 8900

Reflux Flow TMA (kg/h)

Real Optimum

Reflux Flow EstMA (kg/h)

0 2 4 6 8 10 12 14

Time (h)

Figure 28. Evolution of the reflux flow (R) using MA based on transient measurements

Time (h)

Figure 29. Evolution of the distillate flow (D) using MA based on transient measurements

Figure 30. Evolution of the dual constraint

0 10 20 30 40 50 60

Time (h)

Figure 31. Evolution of the process cost function ($p)

52 50 48

H 46 e

44 42 40

SP Thead RTO dual Thead RTO dual SP Thead RTO filtered Thead RTO filtered Model Optimum Real Optimum

Time (h)

Figure 32. Evolution of the head temperature (Thead) 106

104 102 0 100 o 98

0 10 20 Time (h) 30

Figure 33. Evolution of the bottom temperature (Tbottom) 0.45

SP Tbottom RTO dual Tbottom RTO dual SP Tbottom RTO filtered Tbottom RTO filtered Model Optimum Real Optimum

Lower limit a RTO dual RTO filtered

Ih_TT_n-

Time (h)

Figure 34. Evolution of the dual constraint

Table 1. Nomenclature used in the dynamic model

Variable Meaning Units

l Liquid molar flow rate kmol/h

v Vapour molar flow rate kmol/h

f Feed molar flow rate kmol/h

dmol/dt Molar accumulation rate kmol/h

lref Reflux molar flow rate kmol/h

lacum Condenser outlet molar flow rate kmol/h

loverflow Liquid molar flow rate overflowing from the accumulator kmol/h

b Bottom molar flow rate kmol/h

d Distillate molar flow rate kmol/h

R Reflux mass flow rate kg/h

D Distillate mass flow rate kg/h

B Bottom mass flow rate kg/h

F Feed mass flow rate kg/h

S Steam flow to the reboiler kg/h

Pm Molecular weight kg/kmol

xj,i Molar fraction of component j in the liquid stream leaving the tray i °/1mol

yi.i Molar fraction of component j in the vapour stream leaving the tray i °/1mol

yeqj,i Equilibrium concentration °/1mol

Kj,i Equilibrium constant of component j in the tray i

Hv Vapour enthalpy kJ/kg

hi Liquid enthalpy kJ/kg

hj Specific enthalpy of component j kJ/kg

T. Tray temperature K

n feed Feed tray #

Thead Head temperature K

Tbottom Bottom temperature K

Psat,j Vapour pressure of component j Pa

P, Tray temperature Pa

Ef Tray efficiency °/1

Qb Heat flow added to the reboiler kJ/h

Qc Heat removed in the condenser kJ/h

Ah Latent heat of vaporization kJ/kg

Preboiler Reboiler pressure bar

Table 2. N omenclature used in the RTO model

Variable Meaning Units

v Vapour molar flow rate in the enriching section kmol/h

Vs Vapour molar flow rate in the stripping section kmol/h

lS Liquid molar flow rate in the stripping section kmol/h

V Vapour mass flow rate in the enriching section kg/h

Vs Vapour mass flow rate in the stripping section kg/h

L Liquid mass flow rate in the enriching section kg/h

LS Liquid mass flow rate in the stripping section kg/h

j Molar fraction of component j in the liquid stream through the stripping section °/1mol

x Molar fraction of component j in the liquid stream through the enriching section °/1mol

yf Molar fraction of component j in the vapour stream through the stripping section °/1mol

y Molar fraction of component j in the vapour stream through the enriching °/1mol

section

xj,n feed Molar fraction of component j in the feed stream °/1mol

XD,j Molar fraction of component j in the distillate stream °/1mol

Molar fraction of component j in the bottom stream °/1mol

H Enthalpy for liquid stream in the stripping section kJ/kg

hs Enthalpy for vapour stream in the stripping section kJ/kg

Hb Enthalpy of the bottom stream kJ/kg

hD Enthalpy of the distillate stream kJ/kg

Hd Enthalpy of the input vapour stream to the condenser kJ/kg

Pbottom Bottom pressure Pa

Phead Head pressure Pa

ntray Number of trays #

Table 3. Size comparison between steady-state models, rigorous ^ and reduced model.

Rigorous Reduced

Equations 1076 39

Explicit variables 921 29

Algebraic variables 148 3

Boundary variables 7 7

Table 4. Time comparison between steady- state models, rigorous and reduced model

Rigorous Reduced

Optimization time (s) 7.594 0.235

Number of iterations 56 50

Time per iteration (s) 0.135 0.0047

Table 5. Comparison between process and

model optimum

Real optimum Model optimum

S (kg/h) 5558.20 5000.00

R (kg/h) 9389.67 9265.70

D (kg/h) 11179.53 9965.52

Thead (°C) 50.98 42.58

Tbottom ( C) 102.02 96.33

Xd(C3H) (°/1) 0.800 0.819

$ (€/h) 783.19 714.50