M Chèbre

Y Creff

N Petit

Feedback control and optimization for the production of commercial fuels by blending

This paper presents a control algorithm for blending systems. Such systems are used in refining to produce mixtures having specified properties from several components. The underlying control problem is multi-variable, with constraints on the inputs and outputs, and involves large uncertainties. To address this complexity, a constrained optimization problem is formulated, while the uncertainties are treated in closedloop by an estimator of the components properties. Besides a theoretical study of the main algorithm and a study of convergence, the paper presents numerous technical details that are needed to solve blending problems as they appear in refining operations. Among them are a general lack of measurements, variable delays, measurement synchronization, infeasibilities management, hydraulic constraints, and pre-blends. Industrial case-studies are provided and stress the relevance of the approach.

Introduction

The main contribution of this paper is a new method to control blending systems. These systems are used to produce a mixture having some desired properties. This control problem is very common in refineries. The method we propose can be used in various situations where non-reactive components are blended and linearly impact on the properties of the blend, provided nonlinear properties transformations are used (e.g. log scaling).

The proposed method has the advantage of dealing with uncertainties of the components properties. This robustness is of particular interest when time-varying upstream plants are considered. These upstream units produce components with time-varying, unmeasured, unknown or poorly known properties.

From an application point of view, the blending objectives are to produce a mixture with some prescribed properties. For gasoline, diesel or fuel production, these can be the octane number, sulfur concentration, Reid Vapor Pressure among others. Minimization of production costs is also an important topic.

Over the past 40 years, blending control systems have attracted much attention (see [START_REF] Perkins | Method for predictive intrinsic properties of a mixture[END_REF][START_REF] Bay | System for feed blending control[END_REF][START_REF] Feld | Fluid blending system. Patent[END_REF][START_REF] Le Febre | Process for controlling blending[END_REF][START_REF] Singh | Model-based real-time optimization of automotive gasoline blending operations[END_REF]). There has been significant research effort to propose closed-loop strategies using signals from on-line analysers located downstream the blender. It is worth mentioning that, usually, only downstream measurements are considered. The main reason for this choice is to minimize the number, and thus the cost, of required analysers. This point will be illustrated in Section 1. Basic strategies use single-variable controllers (mostly integral effect) in a single-input single-output modeling approach. A priori estimates of the components properties are used to assign the feedback loops. Following a common practice [START_REF] Walton | The Boss of blending[END_REF][START_REF][END_REF], we design a multi-variable strategy to improve performance. As will appear, we propose a genuine approach in which measurements are used to update knowledge on the components properties.

Another point of interest, in real applications, is the constraints. Constraints on the blend properties arise from production requirements. Others include flows limitations and pumping constraints. Simple rules can be used to take these into account one by one. More generally, in our approach, it is also possible to take them into account: we solve a constrained optimization problem.

The starting point of our study originates in the late 70's when, although commercial blending control packages were already available, the problem of robustness against uncertainties on the properties of components had not received any satisfactory answer. Various companies and teams of researchers have been involved since then. TOTAL has inherited a long experience in blending control, which can be traced back to the early 80's and the first online applications of the Anamel software. In the late 90's, it was decided to address the robustness problem: the result of this work is the development of the latest version of Anamel [START_REF] Bernier | Method and device for controlling production of a mixture of components, in particular a mixture with premix dead volumes[END_REF], which is now being used on 17 blenders located in 6 refineries within the TOTAL group and is presented here.

The main algorithm of this latest version of Anamel (in service since 2001) consists of two distinct, though connected, layers: an optimization problem and a feedback loop with an observer. The optimization problem permits to account for various constraints and production cost minimization. The observer is used to partially estimate the components properties in a spirit of adaptive control methods. Both layers are required to provide convergence and guarantee a successful blend. Other approaches based on estimates of the components properties perform a preliminary or in-line full identification of these properties (with laboratory sampling). In particular, one can refer to [START_REF] Perkins | Method for predictive intrinsic properties of a mixture[END_REF] for an estimation procedure performed before the recipe can be updated. On the other hand, one can refer to [START_REF] Feld | Fluid blending system. Patent[END_REF] for a strategy using direct in-line measurements of the components properties. Interestingly, to provide convergence of the blend properties to a prescribed target, Anamel's observer needs not to converge to the actual unknown values of the components properties. This might be the case though, but may not be so common in practice. Most of the time, the blend target is reached before accurate estimates of the components properties are obtained. This is not a concern, because blend properties are definitely the primary target, and should not be surprising to anyone familiar with adaptive control [START_REF] Ioannou | Robust adaptive control[END_REF][START_REF] Åström | A survey of adaptive control applications[END_REF]. Interestingly, this behavior can be analyzed using LaSalle's invariance principle (see e.g. [START_REF] Khalil | Nonlinear Systems[END_REF]) for the underlying dynamical system.

The paper is organized as follows. In Section 1, we present the blending process, the actuators and the available sensors. In Section 2, we expose the control problem. We introduce the notations to be used throughout the paper. In Section 3, we detail the solution and prove its theoretical convergence. Numerous important practical issues are discussed. We expose solutions for them. Finally, industrial results are presented and discussed in Section 4, while we draw some conclusions in Section 5.

Process description

In this section, we present the blending process which is used to obtain finished or semi-finished products from refined components or upstream units flows. All the key elements in the flow-sheet are detailed. We describe the available sensors and actuators. As will appear, the main operational problem is that the components properties are poorly known and usually not measured on line. In fact, only the blend properties are analyzed on line.

Obtaining semi-finished or finished products from refined components

In refineries, semi-finished or finished products (in particular commercial fuels) are not directly extracted as parts of crude oils, but are produced by blending several components. These components are transformed petroleum cuts, resulting from preliminary separations in atmospheric and vacuum columns. Mixing the various components with the right proportions provides the final blend, with properties required by the m specifications of interest. A typical blending process is pictured in Figure 1. Several components (n) are pumped to a blender (also referred to as static mixer) from intermediate storage tanks or pipes (run down). Downstream the blender, the product is either stored in a final product tank, or routed to another refining unit, or shipped off the refinery.

Several elements appear in the flow-sheet. These are set during the line-up procedure. Pipes connect tanks outlets to the blender, with the possibility of "pre-blends". This possibility, pictured in Figure 2, is used to reduce the total length of pipes. Additionally, a flushing tank and several product tanks can be found at the blender outlet. They enable alternative production modes such as sequenced batch (also called "rundown blending" which is frequent when very large productions are considered).

Actuators and sensors

The blend properties can be controlled with the n blender motorized inlet valves. Given a blender outlet total volume flow rate F , the valves openings define a control vector consisting of n volume ratios u = (u 1 , . . . , u n ) T , referred to as the recipe. Low-level flow controllers (named ratio control system) guarantee that this vector tracks any reference signal.

Let us now detail the sensors system. For cost and reliability reasons (sensors drift over time), the components properties are generally not measured on line. Yet, there exist large uncertainties on their values. Among the sources of uncertainty are the possible drifts in the operation of upstream units that may be inaccurately controlled during some periods of time, and slicing phenomena in storage tanks. From time to time, laboratory data provide samples measurements. These measurements have a very low time resolution, and can also be inaccurate and delayed. The main culprits are the semi-manual definition of the time stamp, the location of the sampling point, and faulty data input in the laboratory information system. The most interesting information is actually found downstream the blender. A limited number of analyzers are installed at the blender outlet. With these, the m properties of the blend are measured on-line. Usually, m is larger than n. Yet, depending on the type of sampling loop and the involved measurement technology, these measurements can suffer from large and variable delays.

Control problem

The primary goal of any blending system is to produce a mixture having some specified properties. In other words, the blending system has to find a n-dimensional recipe u such that the m properties of the mixture satisfy some objective.

The instantaneous blend properties are considered as the output of the system. They are denoted by the m-dimensional vector y. Alternatively, if a product tank is considered, it may be desired to control the average (m-dimensional) properties z of this tank.

The components properties are grouped in a m × n matrix B.

The two ends of the flow-sheet are the components inlets, and the blender. Flows through the pipes of the flow-sheet are assumed to satisfy a one dimensional plug-flow model. When they are considered, "pre-blends" consist of several pipes merging at the same location (see again Figure 2).

In this paper, the blending models are assumed to be linear. This assumption is not restrictive, because, up to some change of variables, numerous properties actually satisfy this linearity assumption. Therefore, preliminary vector coordinate wise non linear mappings can be used to validate this assumption.

In particular, at steady-state, the following relation holds

y = Bu (1)
Balance equations can be considered for non steady cases.

Several constraints on the recipe u need to be considered. For mathematical consistency, the recipe vector u coordinates must all lie in [0, 1] and satisfy n i=1 u i = 1. From an operational and economical point of view, u should remain close to a recipe of interest u opt . Further, hydraulic constraints (physical limitations of the pumps and pipes) and components availability impose upper and lower bounds on the coordinates u i , i = 1, ..., n.

Various constraints on the blend properties y need to be considered as well. A reference y r and/or upper and lower bounds are associated to each coordinate of y. From a practical point of view, these bounds can be considered as hard bounds (related to commercial specifications) or soft bounds, which can be violated at the expense of profit losses (also referred to as "give-away").

The matrix of the components properties B is poorly known. Yet, B an initial estimate for it (most frequently given by laboratory samples) and some on-line blend properties measurements are available.

Before we present our control solution and discuss implementation issues, let us summarize our goal.

Problem 1 Given B an initial estimate of B, given real time measurements of the output blend properties y, find a closed-loop control scheme, acting on u, such that y converges to y r and remains between pre-specified bounds. At all times, u must satisfy the operational constraints, and preferably be close to a recipe of interest u opt .

Proposed solution

To solve the control problem discussed in Section 2, we propose a twofold approach. The constraints and the various control objectives are formulated in an optimization problem. Simultaneously, an observer reconstructs an estimate of the components properties. These two parts of the control law closely interact. Under some mild simplifying assumptions, theoretical convergence of this strategy is studied. Interestingly, the reconstructed components properties need not to converge to their actual values to guarantee a successful blend.

Optimization problem

It is considered here, that B an estimate of the components properties matrix B is given. Later, in § 3.2 this matrix will be updated.

As discussed in Section 2, not every blend property needs to match a specified reference. Some of them are simply asked to remain within some given bounds. Values of blend properties associated to specified references y r (dim y r = r ≤ m) can be estimated, using [START_REF] Anamel | Post benefit study. Real-time blending optimization[END_REF], through an (r × n) sub-matrix Br of B. Typically, r ranges from 2 to 5, while n ranges from 5 to 10. Similarly, the blend properties associated to hard and soft bounds can be computed using the sub-matrices Bh and Bs . Lower and upper bounds vector on the hard and soft constraints are noted y h,lb , y h,ub , y s,lb , and y s,ub , respectively. Vector lower and upper bounds on the control vector are noted u lb and u ub . Taking into account the consistency equation n i=1 u i = 1, the recipe of interest u opt , and, most importantly, the blending objectives, we formulate the control problem under the form of the following optimization problem

min u u -u opt 2 Q ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 ≤ u lb ≤ u ≤ u ub ≤ 1 n i=1 u i = 1 Br u = y r y h,lb ≤ Bh u ≤ y h,ub y s,lb ≤ Bs u ≤ y s,ub (2)
where a symmetric definite matrix Q is used to weight the Euclidian norm, i.e. u 2 Q = u T Qu. This matrix is chosen to promote or to penalize the use of some components (further details are given in §3.7.2). The optimization problem (2) is a quadratic programming problem. It can be handled with various software packages such as IMSL [START_REF]User's guide[END_REF]. Its solution gives an open-loop control u.

Feedback

As mentioned earlier, on-line blend properties measurements can be used to update the open-loop control law. These measurements are usually used to monitor the blend, and we wish to take them into account under the form of a feedback into the optimization problem (2).

In the presented context of linear multi-variable control, a classic way to proceed (e.g. in Model Predictive Control [START_REF] Froisy | Model predictive control-building a bridge between theory and practice[END_REF][START_REF] Garcia | Model predictive control: Theory and practice-a survey[END_REF] or Internal Model Control [START_REF] Garcia | Internal model control. 2. Design procedure for multivariable systems[END_REF]), is to introduce an integral term fed by the difference between the measurement and the setpoint, or a filter of the difference between the measurement and the prediction of a model. Such classic procedures fail here. A main reason is the equality constraint on the control vector, which reduces the number of free control parameters by one. The actual gains of the (n -1) remaining control variables are in fact unknown, because they consist of differences between entries of the B matrix which are also unknown. For sake of illustration, a very simple test case is given in Appendix A. Therefore, another approach must be considered.

First, we relate the measurements, which are assumed to be done continuously, to the current values of the control variable by y = Bu. Then, the estimation B of B, which is assumed to be constant, is updated as follows. Considering its j th row Bj , the continuous-time update law is

d BT j dt = -β j Hu Bj u -y j , ( 3
)
where H is the following diagonal scaling matrix (ū being a reference recipe, e.g. a constant value close to u opt )

H = 1 ū ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 ū1 0 . . . 0 1 ūn ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
and β j is a strictly positive parameter. This update law is analogous to those found in adaptive control [START_REF] Ioannou | Robust adaptive control[END_REF][START_REF] Åström | A survey of adaptive control applications[END_REF]. Considering the output relation (1), the essential idea is the comparison of the observed system response Bu with the model output Bu.

Finally, our feedback control law is designed as follows.

Solution to problem 1 (Simple version) Sequentially, solve the optimization problem [START_REF] Åström | A survey of adaptive control applications[END_REF] and update the estimate of the components properties B when new measurements are available. The discrete-time update is the sampled version of the continuous-time update law [START_REF] Åström | Adaptive control[END_REF].

We now prove convergence of this strategy under the assumption that the optimization problem (2) is always feasible, i.e. admits a (necessarily unique because of strict convexity) solution.

Convergence

Solution 1 combines an on-line parameter estimator (3) and a control law which is defined as the solution of the optimization problem [START_REF] Åström | A survey of adaptive control applications[END_REF]. From this description, it can be viewed as an (indirect) adaptive controller [START_REF] Ioannou | Robust adaptive control[END_REF]. As will now appear, B is continuously adjusted so that Bu(t) approaches Bu(t) as t → +∞. Yet, no particular effort is made to design the input u(t) so that B converges toward B, as would normally be desired in an on-line parameter estimation technique. This is not one of the objectives, as it could cause large variations of the input signal u(t) (e.g. to satisfy some persistency of excitation property [START_REF] Ioannou | Robust adaptive control[END_REF][START_REF] Åström | Adaptive control[END_REF][START_REF] Khalil | Nonlinear Systems[END_REF]).

Consider, for property j = 1, . . . , m, the scalar function (Lyapunov function candidate)

Ψ( Bj ) = 1 2 ( Bj -B j )H -1 ( Bj -B j ) T
This function is strictly positive away from B j , where it equals 0. Its timederivative along the trajectories of (3) is

dΨ dt ( Bj ) = -β j Bj u -y j 2 ≤ 0
Therefore, Ψ( Bj ) is a Lyapunov function for system (3) (see [START_REF] Khalil | Nonlinear Systems[END_REF]). From LaSalle's invariance principle, for any initial condition, the solution of system (3), Bj (t), converges when t → +∞ towards the largest invariant set of (3), included in the subset { Bj s. t. dΨ/dt( Bj ) = 0}. Therefore, Bj (t) converges in a way such that B j u = Bj u. Yet, by definition of the optimal control problem [START_REF] Åström | A survey of adaptive control applications[END_REF], which is assumed to possess a solution (which is necessarily unique), Bj u satisfies the blend objectives. Therefore, so does y j (t) = B j u(t), in the limit as t → +∞. The same reasoning applies to all the blend properties. In summary, the blend is successful, even though Bj does not converge to B j . In details, we have

lim t→+∞ Br (t)u(t) = y r , y h,lb ≤ lim t→+∞ Bh (t)u(t) ≤ y h,ub , y s,lb ≤ lim t→+∞ Bs (t)u(t) ≤ y s,ub
while the equality lim t→+∞ B = B might not hold.

Handling infeasibility with cascaded optimization problems

The (preliminary) solution to problem 1 we propose converges provided that the optimization problem ( 2) is always feasible. Certainly, it is a reasonable assumption to require that the blend is feasible using the components at our disposal. Yet, one should notice that the optimization problem is not written in terms of the components properties matrix B, but in terms of its estimate B. Thus, the feasibility assumption bears on an estimate, which, as it changes over time, might reveal troublesome.

For real-time applications, it is required that the control algorithm always provides an answer, even a partly disappointing one when it seems impossible to satisfy every control objective. When the optimization problem (2) is not feasible, it remains possible to solve easier problems, with a smaller number of constraints. As already discussed, the output constraints are sorted with respect to their relative importance. This ranking is a common formulation of blending operators. When the convex optimization problem (2) is not feasible, the constraints can be relaxed according to the following procedure (see Appendix B for details)

• Step 1: Hard constraints management. Lower and upper bounds on u are considered, along with the consistency equation n i=1 u i = 1. The equality constraints for the references and the inequality constraints corresponding to the soft constraints are totally forgotten. Hard constraints are substituted with (positive) slack variables s appearing in the modified objective function. These slack variables are used as penalty functions to be minimized. The target recipe u opt does not appear in this new objective function. This problem is always feasible. Its solution (u * , s * ) gives achievable values for the hard constraints. The resolution is achieved sequentially. The typical problem solved at iter-ation k + 1 is given by min u,s

Bh k+1 u + s -y h,ub k+1 2 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 ≤ u lb ≤ u ≤ u ub ≤ 1 s ≥ 0 n i=1 u i = 1 Bh 1k u ≤ y hach,ub 1k
where Bh k+1 and y h,ub k+1 correspond to the hard constraints treated at iteration k + 1; Bh 1k correspond to the hard constraints relaxed at iterations 1 to k, the values of which being given by y hach,ub 1k . Only maximum bounds have been considered here to simplify the formulation.

• Step 2: References management. Lower and upper bounds on u are considered, along with the consistency equation n i=1 u i = 1. The previously updated hard constraints, now achievable, are also considered, but the soft constraints are totally forgotten. References simply appear in the objective function, while the target recipe u opt does not appear in the formulation. This problem is also always feasible. Its solution u * * gives achievable values for the references. The typical problem solved in this step is given by

min u Br u -y r 2 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 ≤ u lb ≤ u ≤ u ub ≤ 1 n i=1 u i = 1 Bh u ≤ y hach,ub
where the hard constraints updated in Step 1 are given by y hach,ub . Only maximum bounds have been considered here to simplify the formulation.

• Step 3: Soft constraints management. Lower and upper bounds on u are considered, along with the consistency equation n i=1 u i = 1. The previously updated hard constraints are still considered, as well as the updated references that are now achievable. Soft constraints are substituted with (positive) slack variables s appearing in the modified objective function. The target recipe u opt does not appear in the formulation. This problem is always feasible. Its solution u * * * is applied to the blending process. The resolution is achieved sequentially. The typical problem solved at iter-ation k + 1 is given by min u,s

Bs k+1 u + s -y s,ub k+1 2 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 ≤ u lb ≤ u ≤ u ub ≤ 1 s ≥ 0 n i=1 u i = 1 Bh u ≤ y hach,ub Br u = y rach Bs 1k u ≤ y sach,ub 1k
where Bs k+1 and y s,ub k+1 correspond to the soft constraints treated at iteration k + 1; Bs 1k correspond to the soft constraints relaxed at iterations 1 to k, the values of which being given by y sach,ub 1k ; y rach correspond to the achievable references determined in Step 2. Only maximum bounds have been considered here to simplify the formulation.

In practice, hard and soft constraints are subdivided into sets of decreasing priorities or ranks. Steps (1) and ( 3) are treated iteratively, each iteration involving a set of constraints with equal priorities. In summary, these steps always produce a recipe u which satisfies as many of the constraints appearing in the original optimization problem as possible. This is interesting because it is consistent with the requirements of real-time control.

Accounting for delays

As had appeared earlier (in §3.3), the feedback law uses an error between predicted measurements and actual ones (see Equation ( 3)). The delays, that are the main sources of mis-synchronization, must be carefully accounted for, for the error to be correctly computed. These delays either come from the measurement process, or from dead volumes (pre-blends) that induce transportation delays. Interestingly, the delays are either known or can be estimated with a good accuracy. In practice, it appears that we do not have to use specific control methods robust to delay uncertainty (one can refer to [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF] for an overview of such methods).

Measurement delays

The measurement delay may induce important mismatches and inconsistencies in the above error calculation. Fortunately, measurements come along with sampling dates (time stamps) and it is possible to compare them against the right (delayed) predictions. The formula (3) is modified accordingly.

Variable transportation delays

In practice, transport delays must be accounted for in the preceding methodology. Interestingly, these transport delays are not constant. They actually depend on the control variables, which are themselves updated by the control algorithm. Under a plug-flow assumption, this dependance can be fully determined.

We now expose means to calculate the transport delays in a very straightforward manner which is eventually implemented.

First, let us consider a simple example, pictured in Figure 3. Consider two connected pipe sections. Their volumes are noted V 1 and V 2 , respectively. We note b 1 , b 2 and b 3 the composition of the fluid flowing through this pipe. These variables have distinct values when the fluid contained in the two pipe sections is not homogeneous (e.g. due to upstream transients). We note F the time varying flow rate. Under the assumption of incompressibility, the volumic flow rate is spatially uniform. The compositions are constrained by the following

b 1 (t) b 2 (t) b 3 (t) V 2 V 1 F (t)
(t) = b 1 (t -δ 1 (t)), with V 1 = t t-δ 1 (t) F (τ )dτ b 3 (t) = b 2 (t -δ 2 (t)), with V 2 = t t-δ 2 (t) F (τ )dτ
The preceding integral relations (see also [START_REF] Petit | Motion planning for two classes of nonlinear systems with delays depending on the control[END_REF][START_REF] Petit | Systèmes à retards, platitude en génie des procédés et contrôle de certaines équations des ondes[END_REF]) implicitly define the varying delays δ 1 and δ 2 . Therefore, the upstream and downstream properties satisfy the following equality

b 3 (t) = b 1 (t -δ 2 (t) -δ 1 (t -δ 2 (t)))
These relations can be generalized to cases of interest in the context of blending processes such as the one depicted in Figure 4.

In the setup presented in Figure 4, five storage tanks are connected to a single blender through a network of pipes using three pre-blends. Three hold-ups

(b 1 , f 1 ) (b 2 , f 2 ) (b 3 , f 3 ) (b 4 , f 4 ) (b 5 , f 5 ) (b 6 , f 6 ) (b 7 , f 7 ) (b 8 , f 8 ) V 1 V 2 V 3 y
f 8 = f 1 + f 2 + f 3 .
Following the calculations presented in the introductory case, transport delays can be derived from the following equations

V 1 = t t-δ 1 (t) f 6 (τ )dτ, V 2 = t t-δ 2 (t) f 7 (τ )dτ, V 3 = t t-δ 3 (t) f 8 (τ )dτ
Then, intermediate blends properties are easily determined by

b 6 (t) = f 4 (t -δ 1 (t)) f 6 (t -δ 1 (t)) b 4 + f 5 (t -δ 1 (t)) f 6 (t -δ 1 (t)) b 5 b 7 (t) = f 1 (t -δ 2 (t)) f 7 (t -δ 2 (t)) b 1 + f 2 (t -δ 2 (t)) f 7 (t -δ 2 (t)) b 2 b 8 (t) = f 7 (t -δ 3 (t)) f 8 (t -δ 3 (t)) b 7 (t -δ 3 (t)) + f 3 (t -δ 3 (t)) f 8 (t -δ 3 (t)) b 3
Eventually, the end-product property is easily deduced, through a cascade of mixing laws, by combining the time-varying transport delays. After some easy algebra, the sought-after relation is

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (f 1 (t) + f 2 (t) + f 3 (t) + f 4 (t) + f 5 (t))y(t) = f 8 (t) f 7 (t -δ 3 (t)) f 8 (t -δ 3 (t)) f 1 (t -δ 2 (t -δ 3 (t))) f 7 (t -δ 2 (t -δ 3 (t))) b 1 + f 8 (t) f 7 (t -δ 3 (t)) f 8 (t -δ 3 (t)) f 2 (t -δ 2 (t -δ 3 (t))) f 7 (t -δ 2 (t -δ 3 (t))) b 2 + f 8 (t) f 3 (t -δ 3 (t)) f 8 (t -δ 3 (t)) b 3 + f 6 (t) f 4 (t -δ 1 (t)) f 6 (t -δ 1 (t)) b 4 + f 6 (t) f 5 (t -δ 1 (t)) f 6 (t -δ 1 (t)) b 5
(4) Interestingly, in the case of constant flow rates from the storage tanks, this last formula can be simplified down to the following usual relation

(f 1 + f 2 + f 3 + f 4 + f 5 )y(t) =f 1 b 1 + f 2 b 2 + f 3 b 3 + f 4 b 4 + f 5 b 5
More generally (see Appendix C), Equation ( 4) takes the form

y(t) = i=1,...,n U i (t)b i , ( 5
)
where U i (t) can be expressed in terms of (ratios of) values of present and past values of components of the vector recipe. Vectors b i are the columns of the properties matrix B. Denoting U(t) the n-dimensional column vector with entries U i (t), we simply have y(t) = BU(t). This expression is a substitute to Equation ( 1) for systems with pre-blends. It is used in estimator (3), which is transformed into:

d BT j dt = -β j HU Bj U -y j , ( 6
)
It can be easily shown that this modified estimator still converges to the same invariant set [START_REF] Barraud | Commande de procédés à paramètres variables[END_REF].

Actual control algorithm

The proposed solutions to the issues raised in § 3.4 and in § 3.5 are incorporated in the control algorithm. This yields the following solution.

Solution to problem 2 (Actual version)

Sequentially, solve the optimization problem (2), using whenever required the procedure described in § 3.4 to handle infeasibilities, and then update the estimate of the components properties B when new measurements are available. The discrete-time update is the sampled version of the continuous-time update law [START_REF] Bernier | Method and device for controlling production of a mixture of components, in particular a mixture with premix dead volumes[END_REF].

Further enhancements could be considered. It can be noted that, in the optimization problem (2), the expressions y = Bu are still used in all cases, with or without pre-blends. This can be viewed as solving a predictive control problem with a control horizon equal to 1 ("single move"), and a single prediction point at infinity, as for a given τ , U i (t + τ ) = u i (t), ∀i ∈ {1, . . . , n}.

It is possible to extend problem (2), by using ( 5) to express the output constraints at various times from current time t to t + τ . This can be done at the expense of introducing strong nonlinearities in the optimization problem. Such an approach has been tested through simulations [START_REF] Barraud | Commande de procédés à paramètres variables[END_REF], but has not been implemented, as it did not provide any significant improvement for problems of practical interest.

Further functionalities

A blending controller has to propose more functionalities than those described so far in this paper. For instance, the total flowrate F (t) is generally handled, for sake of output maximization. Additives have to be managed for diesel fuels. Two other important features are the quality control of storage tanks, and the recipe optimization (the way to choose u opt is equation ( 2)). These two features are briefly presented in the next sections. All the functionalities stated above are included in Anamel. The control strategy presented in the previous sections focuses on tracking a given (instantaneous) blend properties reference. This reference need not be constant, and can be updated over time. It is possible to cascade the main control algorithm with a trajectory tracking strategy. This can be done to manage blends in storage tanks. A simple (generic) scenario with one storage tank is reported in Figure 5. Note V the current volume of fluid contained in the storage tank, z its composition. Consider H a finite time horizon over which the total blender flow rate F and the blend properties y are assumed constant. The storage tank properties values satisfy the following balance equation

(V + F H)z(t + H) = V z(t) + F Hy.
It is desired that z(t) tracks a reference signal z r (t). Then, we set z(t + H) = z r (t) to define the reference y r (t) as

y r (t) = z r (t) + V F H (z r (t) -z(t))
In practice, z(t) is not measured, but frequently estimated through the balance equation given above.

In the case when the storage tank is not empty at the beginning of the blend, a heel management strategy is used. Based on the same balance equation, it consists of defining a reference trajectory originating at this known value.

Recipe optimization

When problem ( 2) is feasible, the objective function can be used to optimize the recipe, i.e. to promote or penalize the usage of some specific components, which also impacts on the cost of the blend. Using several techniques detailed in [START_REF] Bernier | Method and device for controlling production of a mixture of components, in particular a mixture with premix dead volumes[END_REF], it is possible to push the blend properties toward their specified hard constraints. We briefly sketch an overview of such techniques.

When no particular recipe optimization is desired, u opt is either set equal to the initial value provided by the scheduling department of the considered refinery, or equal to the filtered current recipe u f , to avoid oscillations. When recipe optimization is an objective, u opt equals u f + ζ, where the entries of ζ are positive (resp. negative) for components to promote (resp. to limit) in the blend. When the cost of the blend must be minimized, a price is associated to each component. Then, the prices of the components with respect to the cost of the current blend, are used to set ζ entries. Details can be found in [START_REF] Bernier | Method and device for controlling production of a mixture of components, in particular a mixture with premix dead volumes[END_REF].

Industrial results

The algorithm detailed in the previous section, has been implemented in C, and installed in numerous refineries of TOTAL under the name "Anamel". Thanks to its relatively low complexity, the whole control scheme can be run on standard control systems at a rate consistent with the process dynamics (e.g. every 1 to 5 minutes).

Statistical study

After an initial installation phase of this Anamel software in TOTAL refineries from 2001 to 2004, a statistical study of the performance of the blends was conducted at the TOTAL Leuna (Germany) refinery, which was then the latest of this installation program. The study was done so as to compare the results of the new on-line blending optimizer that was installed in 2004 in Germany: all the blends of the year 2003 was thus compared to the blends of year 2005 (with 405 gasoline blends, 425 diesel blends, and 173 fuel blends in 2005). Enhancements are detected on gasoline blends with a higher proportion (42% of improvement for Reid Vapor Pressure RVP and 30% for octane number MON) of the blends closer to the hard contraint (lower give-away for monitored qualities). On diesel blenders, a better convergence with less oscillations was obtained. In particular, the cloud point regulation has been improved significantly. Anamel users in this refinery are convinced that the new release of the software has contributed to improve the performance of the overall blend operations.

Besides these general comments, we now briefly present two representative test cases. The first one is a continuous blend upstream a process unit for diesel production, the second one is a batch production of gasoline. Today (2008), Anamel is installed on 17 distincts blending systems (5 for gasoline production, 7 for diesel production, and 5 for fuel production).

Industrial case 1: continuous blending upstream an HDS unit

In the industrial case presented here, the product is diesel fuel. The blender outlet is directly connected to an hydro-desulphurization unit (HDS). In this context, the instantaneous properties of the blend are the variables that must be controlled.

In details, 6 components are used out of which 3 (representing more than 20% of the total volume flow rate) are set to constant flow rates and cannot be manipulated. There remains a n = 3 dimensional recipe vector. One control objective is to maximize the incorporation of component 1.

As before, Anamel is run (according to Solution 2) every 5 minutes. The recipe is subjected to several contraints including hydraulic contraints, min. and max. bounds. Also, its rate of change is limited.

Three properties are of particular interest. These are the flash point, and the cloud point which must track their setpoints and the sulphur content which must remain below some max. value. Further, three other properties must be monitored. These are the density (D15), the cetane index, and the ASTM360.

The results are presented in Figure 6. The recipe histories, the flash point histories, the cloud point histories, and the sulphur rate are presented in (a), (b), (c) and (d), respectively. The D15, cetane index and ASTM360 are easily kept inside their bounds. They are not reported in the plots. The time horizon reported here is approximately 24 hours (one day).

The results stress that the flash point and the cloud point are kindly controlled. Over time, the controller tends to maximize the rate of component 1, while keeping the sulphur rate below its maximum. At the beginning of the trends, the flash point is below its minimum. This induces a strong control action which vanishes once the problem is solved.

The case reported here is representative of the behavior of the proposed control scheme used on continuous blends. In the presented case of blending upstream an HDS unit, several beneficial aspects have been noted. Gradually, it has been possible to increase the overall feed rate of the unit by 60% over 4 years. This is mainly due to a revamping and a debottlenecking of the unit. But according to end-users in the refinery, it can also be considered as an indirect effect of the controller which keeps the properties of interest within the prescribed bounds and thus minimizes the input disturbances of the HDS unit. It also allowed to maximize the throughput without causing any give-away, and to reduce the variability of the feed rate. Together with a better cloud point control, the blending optimization and the increase of the feed rate have yielded substantial profits which are approximately of 4.3 million Euros/year for a single HDS unit [START_REF] Anamel | Post benefit study. Real-time blending optimization[END_REF].

Industrial case 2: batch gasoline blend

The second case presented here is a gasoline blend. In the considered setup, the blender outlet is connected to a product tank, the properties of which are of interest. Before the blend is started, the product tank already contains about 20000 cubic meters of gasoline. A remaining amount of 15000 cubic meters has to be added during the blend. In this situation, the heel management technique discussed in §3.7.1 is of importance. A total of 7 components are used, out of which 2 (representing 25% of the total volume flow rate) are set to some given ratios (in percentage) that are to remain unchanged during this operation. The recipe dimension is n = 5. To maximize profits, the overall flow rate must be optimized, the target being 1000 cubic meters per hour. No particular recipe optimization is desired. As discussed in § 3.7.2, the target is the initial recipe. A total of m = 9 properties are considered. Two properties need to be controlled to their references. These are the motor octane number (MON) and the vapor pressure (VP). During transients, their values (in the storage tank) must remain between some min. and max. bounds. Besides, other properties such as the research octane number (RON), the gravity, the sulphur content, the benzene content, the olefins contents, the aromatics contents, and the 70% distillation point must be monitored. They correspond to hard and soft bounds.

Anamel is run (according to Solution 2) every 5 minutes. Further, as discussed in Section 2, the recipe is required to satisfy several contraints such as: hydraulic constraints (which can be formulated using incorporation ratios depending upon the current total volume flow rate), min. and max. values, and bounded rate of change (from one run to the next).

The results are presented in Figure 7. The recipe u histories, the MON histories, and the VP histories are presented in (a), (b), and (c), respectively. The total blend time is approximately 15 hours. Early, both MON and VP converge to their setpoints. To achieve these very good results, it is necessary to use varying instantaneous MON and VP setpoints. The recipe varies significantly during the batch, while remaining strictly within its bounds. Noises on the measurements do not impact on the convergence of the algorithm.

The control method is quite robust. This is particularly visible when, after 2.75 h, one of the analyser fails and provides unusable information. Simply, Anamel "freezes" the value of the recipe until the analyser recovers, which happens at t = 3.25 h.

During the first 6 hours (approximately), the system is rather difficult to control because the tank has been filled with off-specifications products and the MON and VP are rather far from their objectives. In particular, the VP in the product tank violates the allowed maximum at the start of the blend. This infeasibility is treated using the technique presented in § 3.4. After about 6 hours, the properties have almost reached their targets. Then, the recipe optimization methods starts, and smoothly brings the current recipe close to the initial recipe.

At the end of the batch, the end-user is left with a tank filled with the desired amount of product with desired specifications and a recipe close to the target one. During the batch all the constraints have been satisfied, and despite the sensor failure, the controller has provided a smooth and efficient transient.

Conclusion and perspectives

The presented algorithm addresses the problem of blending control in a way which is new in this field. The algorithm consists of two connected layers forming an adaptive controller: one is an estimator which partially estimates the plant parameters, and one is an optimization problem solver. Available measurements are used in a feedback loop passing through the estimator, which does not need to converge to the actual components properties to guaran-tee the success of the blend. This fact is of great importance in applications.

In the presented context of industrial blending, measurement synchronization has appeared as an important issue, and the algorithm incorporates a solution to compute input-dependent delays in an effective and accurate way. Handling the infeasibilities of the optimal control problem formulated to solve the multi-variable control problem has also revealed an important topic, especially for real-time implementation. The software package resulting from this work, "Anamel" is now used in numerous refineries. It is constantly improved further and is the subject of active research efforts [START_REF] Chèbre | On line blend monitoring tools in TOTAL refineries[END_REF][START_REF] Chèbre | Feed control for an hydrodesulphurization unit using Anamel blend optimizer[END_REF].

A A test case for simple feedback loops This section proposes an example to discard the classic methods discussed in §3.2. As indicated in this section, a classic control technique is to introduce an integral term fed by the difference between the measurement and the setpoint (Strategy 1), or a filter of the difference between the measurement and the prediction of a model (Strategy 2). Due to the uncertainties in the components properties, the gains of the systems are poorly know. As will appear, it is difficult to distinguish situations in which closed-loop stability can be guaranteed. To illustrate this point, let us consider the following simple example. Assume

C Varying delays

Consider a blend network with p pre-blends. Note

• f i (t) the volume flow rate from the storage tank i at time t, with i ∈ {1, . . . , n}; • f n+i (t) the total volume flow rate through the dead volumes i at time t, with i ∈ {1, . . . , p}; • F (t) the total volume rate of the blender at time t, F (t) = i=1,n f i (t); • V i the dead volume associated to pre-blend i; • b i the properties of component i.

To each component i, one can associate a path Π i defined by the sequence of the p i dead volumes that are passed by from the storage tank to the blender. This sequence consists of p i distinct integers corresponding to the numbering of pre-blends, Π i = {π 1 i , π 2 i , . . . π p i i }, with π j i ∈ {1, . . . , p} for all j ∈ {1, . . . , p i }. Π i = ∅, i.e. p i = 0 means that component i directly feeds the blender.

For each pre-blend i, one can define the set Γ i of its q i input flow rates. This set consists of q i distinct integers corresponding to the numbering of the volume flow rates, Γ i = {γ 1 i , γ 2 i , . . . γ q i i }, with γ j i ∈ {1, . . . , n+p} for all j ∈ {1, . . . , q i }.

When Π i = ∅, the properties b i appear under the form f i (t)/F (t)b i in the blend equation. The case Π i = ∅ is different. Consider pre-blend π j i , its total flow rate is f π j i (t) = k∈Γ π j i f k (t). The transportation delay in this pre-blend is given by δ π j i (t) with

V π j i = t t-δ π j i (t) f π j i (τ )dτ
Let us define Δ j i (t) : t → tδ π j i (t), for all π j i in Π i . For a given i, the composition of these functions is

Δ k,j i (t) Δ k i (Δ j i (t)) : t → t -δ π j i (t) -δ π k i (t -δ π j i (t))
and Δ l,k,j i (t) Δ l i (Δ k,j i (t)). Note

U i (t) = f π p i i (t) F (t) f π p i -1 i (Δ p i i (t)) f π p i i (Δ p i i (t)) f π p i -2 i (Δ p i -1,p i i (t)) f π p i -1 i (Δ p i -1,p i i (t)) • • • f π 1 i (Δ

Fig. 1 .Fig. 2 .

 12 Fig.1. Typical blending process. Components from storage tanks are pumped to a blender to be stored in a final product tank.

Fig. 3 .

 3 Fig. 3. Transport phenomenon through two connected pipe sections.

delayed equalities b 2

 2

Fig. 4 .

 4 Fig. 4. Transport phenomenon with pre-blends. play a role in the corresponding transport phenomena. These are noted V 1 , V 2 , and V 3 . Besides the storage tanks compositions b i=1,...,5 , one has to consider intermediate blends properties. These are noted b 6 , b 7 , and b 8 . Flow rates from the storage tanks are noted f i=1,...,5 . Intermediate flow rates are f 6 = f 4 + f 5 , f 7 = f 1 + f 2 , and f 8 = f 1 + f 2 + f 3 .

Fig. 5 .

 5 Fig. 5. Reference trajectories management in the case of storage tank composition control.

1 ]

 1 that 3 components are mixed to produce a blend with a single property of interest. Note b 1 , b 2 , b 3 the components properties. It is desired to track a setpoint y r . In the optimization problem (2), no hard nor soft bounds are considered. Finally, u lb = 0 and u ub = 1. Let y r = 5.7 and u opt = [0.3, 0.4, 0.3]. By going through the iterations of the control system, various behaviors can be observed, depending on the value of the above mentioned parameters. With [b 1 , b 2 , b 3 ] = [4, 5, 6] and an initial estimate [ b1 , b2 , b3 ] = [4, 6, 7], both Strategy 1 and Strategy 2 applied repeatedly provide convergence. With the exact same setup except [b 1 , b 2 , b 3 ] = [6.5, 5, 6], both Strategy 1 and Strategy 2 fail. It is not a feasibility issue. To get some insight into this phenomenon, consider Strategy 1. As already discussed, the equality constraint reduces the number of free control parameter to 2. Consider that u 1 , and u 2 are used. Their actuals gains are [b 1b 3 , b 2b 3 ] while they are estimated to [ b1 -b3 , b2 -b3 ]. In the first case, these quantities evaluate to [b 1b 3 , b 2b 3 ] = [-2, -1], [ b1 -b3 , b2 -b3 ] = [-3, -1] while, in the second case [b 1b 3 , b 2b 3 ] = [0.5, -1], [ b1 -b3 , b2 -b3 ] = [-3, -It can be shown, through a simple analysis of the optimization problem, that the solution [u 1 , u 2 , u 3 ] is such that [u 1u 3 , u 2u 3 ] = α[ b1 -b3 , b2 -b3 ],where α ∈ R. The induced improvement in the tracking error is[u 1u 3 , u 2u 3 ].[b 1 -b 3 , b 2 -b 3 ] t instead of the expected [u 1 -u 3 , u 2 -u 3 ].[ b1 -b3 , b2 -b3 ] t . If [b 1b 3 , b 2b 3 ].[ b1 -b3 , b2-b3 ] t < 0, the closed loop strategy diverges instead of converging.

  2,...,p i

		i	(t))	f i (Δ 1,2,...,p i i	(t))
		f π 2 i (Δ 2,...,p i i	(t))	f π 1 i (Δ 1,2,...,p i i	(t))
	We have U i (t) = u i (t) = f i (t)/F (t) for Π i = ∅ and
	n
	y(t) =	U
	i=1

i (t)b i .

Acknowledgements

The authors wish to thank the following persons: Claude Jablon, Philippe Zaccagnino, Erik Michelsen, Patrick Caillon, Jean-Claude Bonneau, Jean-Louis Frin, Edith di-Crescenzo, Guillaume Pitollat, Bruno Romeyer, Jean-Claude De Wit, Fari Djenab, Jean-Marie Gorszczyk, David Vissière, Julien Barraud, Vincent Bachelot, Michel Bernier, Valérie Guyard, Tarik Idaoubella, Christophe Legros, Stéphane Bourgeon, Cyrille Mirgain, Roland Forterre, Pierre Hoornaert, Bertrand Ecale, Elena Braceras, Lionel Armand, and the TOTAL Grandpuits refinery west-side operators.

B Cascaded optimization

Consider the following simple optimization problem min u J(u)

where y 1 and y 2 are scalar constraints. Infeasibility means that the two inequality constraints can not be both enforced for u lb ≤ u ≤ u ub . It can be important however to provide a solution which satisfies them "at best". Assume that the first inequality constraint is more important than the second one. Introducing the slack variable s, we can consider the problem min

It is always feasible, with solution (u * , s * ). s * > 0 means that satisfying b T 1 u ≤ y 1 is possible with u lb ≤ u ≤ u ub . But s * = 0 means that b T 1 u * , equal to or greater than y 1 , is the nearest achievable value from y 1 for u lb ≤ u ≤ u ub . Accordingly, we set the achievable value y 1,ach as y 1,ach = y 1 if s * > 0 and

This problem is always feasible, with solution (u * * , s * * ). u * * can be seen as an acceptable relaxed solution of the initial problem. The cascaded optimization problems described in section 3.4 are solved in a similar way, from highest priority hard constraints to lowest priority soft constraints. Notice that the slack variable s is not required to deal with the equality constraints used for the references. . Initially, the storage tank is more than half-filled with an off-specifications product. This situation is handled by the heel-management strategy. The total blend time is 15 hours. Convergence of both the MON and the VP to their setpoints is achieved early. Despite large noises, the proposed algorithm is quite robust and computes control values (a) that remain within their bounds. At the end of the blend, the storage tank is filled with the desired amount of product with the desired properties, while the recipe is close to the target.


Michael Geiger

Lamissa Diabaté

Lionel Ménard

Lucien Wald
email: lucien.wald@ensmp.fr.

Controlling the quality of solar irradiation data by means of a web service

The control of the quality of the irradiation data is often a prerequisite to the further processing of the data. Though data are usually controlled by meteorological offices, the sources are so numerous that the user often faces time-series of measurements containing questionable values. As customers of radiation data, we established our own procedures to screen time-series of measurements. Since this problem of quality control is of concern to many researchers and engineers and since it is often a lengthy and tedious task, we decided to make this screening procedure available to everyone as a web service. This service is the purpose of this paper. The objective is not to perform a precise and fine control, an objective out of reach without details on the site and instruments, but to perform a likelihood control of the data and to check their plausibility. This is achieved by comparing observations with some expectations based upon the extraterrestrial irradiation and a simulation of the irradiation for clear skies. This service is available to everyone on the Web site www.helioclim.net. It offers a very convenient mean to check time-series of irradiation: data are input in a HTML page by a copy and paste procedure and the return is also a HTML page that can be analyzed in details for the data flagged as suspicious.

Controlling the quality of irradiation measurements

The measurements of solar global radiation on the ground level are made primarily with pyranometers. These hemispherically sensitive sensors measure the downward solar radiation flux density on a horizontal surface. The properties of pyranometers which are of concern when evaluating the accuracy and quality of radiation measurements are: sensitivity, stability, response time, cosine response, azimuth response, linearity, temperature response, and spectral response [START_REF] Esra | European Solar Radiation Atlas[END_REF]ISO-9060 1990). The deviations are primarily caused by the directional sensitivity of the individual receiver surfaces and by effects of the glass dome which covers the upwardfacing sensor. Some errors that appear are due to a dust layer or snow or ice on the dome. Other errors may occur from inaccurate horizontal adjustment of the receiver surface, moisture inside the dome, etc. Achievable relative uncertainties at the 95 per cent confidence level are about 3, 8 and 20 per cent for operational pyranometers of high, good and moderate quality in case of hourly sums. In case of daily sums, relative uncertainties of 2, 5 and 10 per cent can be achieved (WMO, 1981).

When taking observations and processing their results, errors may be made which should be detected and, as far as possible, corrected. A reliable estimate of errors is only possible by the final user, when information about the pyranometer is available. This is seldom the case. The data from radiation measurements should be checked in two stages: technical control, and quality control. When observations are acquired from national meteorological services, it is usually assumed that the measurements first went through a technical control by these services. The process of quality control is distinct from technical control. Not only is the correctness of arithmetic calculations checked but also agreement of every radiation parameter with the sun elevation and other parameters. The quality control of radiometric observations is guided by basic physical principles which determine the relations between individual radiation parameters; it may also take into account statistical knowledge on the spatial and temporal variability of these parameters [START_REF] Esra | European Solar Radiation Atlas[END_REF].

Quality control may be a lengthy task. Most often, customers do not have the necessary information to perform a precise and fine control. Accordingly, they are trusting the procedures used by the meteorological offices and are ready to use the data in confidence. However, the sources of information are very numerous and do not exhibit the same reliability. It often happens that the user faces time-series of measurements containing questionable values. The reasons are numerous, ranging from non-standard values for the absence of measurements (e.g., 0 instead of -999) to typing errors if the processing chain of measurements is not entirely digital.

As customers of radiation data, we felt the need to establish our own procedures to screen time-series of measurements. The objective is not to perform a precise and fine control, an objective out of reach without details on the site and instruments, but to flag the data that are obviously questionable. Examples are null values when -999 is expected, values much larger than the extraterrestrial irradiation or positive values when the sun is well below the horizon. The aim of these procedures are to perform a likelihood control of the data and to check the plausibility of the measurements. This is performed by comparing observations with those expected. Since this problem of quality control is of concern to many researchers and engineers, we decided to make this screening procedure available to everyone as a web service. This service is the purpose of this paper.

Quality control procedures

The quality control procedure adopted in the present paper meets the guidelines of the World Meteorological Organisation (WMO 1981) and implements some of the procedures used in the construction of the European Solar Radiation Atlas [START_REF] Esra | European Solar Radiation Atlas[END_REF]. The procedure applies to daily or hourly sums of global irradiation. It applies to individual values or time-series. Because of the limited amount of information, the procedure is restricted to the analysis of each measurement with respect to the extraterrestrial irradiation and the extreme that can be expected for this particular instant. The procedure is exactly the same for daily and hourly irradiation.

Quality control for maximum values of daily sums of global solar irradiation

First of all, each observed daily sum G d (resp. hourly sum G h ) should be less than the daily (resp. hourly) extraterrestrial irradiation received by a horizontal plane located at the same geographical co-ordinates. Otherwise it is flagged as questionable.

In addition, the daily irradiation should only exceed by a small amount the daily irradiation G cd (resp. hourly irradiation G ch ) likely to be observed under exceptionally clear skies with a high atmospheric transparency. The clear sky irradiation G cd (resp. G ch ) is assessed by the means of the ESRA model for clear sky. This model is based upon parameterization formula by Kasten (in ESRA 1984 and[START_REF] Kasten | The Linke turbidity factor based on improved values of the integral Rayleigh optical thickness[END_REF] and has been detailed in [START_REF] Rigollier | On the clear sky model of the 4 th European Solar Atlas with respect to Heliosat method[END_REF]. Corrections were brought to this model as proposed by John Page (personal communication). The formula for the Rayleigh optical thickness in Rigollier et al. behaves incorrectly with terrain altitude in the original model. This was not evidenced in this paper because the sites used for validation have altitudes less than 500 m. Europe as a whole offers a few areas of altitudes larger than 1000 m. Other continents such as Africa or Asia have large areas of larger altitude and such corrections are necessary. Using the notations of Rigollier et al., the beam transmittance under cloudless skies is given by their Equation 2:

exp [-0.8662 T L (AM2) m  R (m)]
(1)

where T L (AM2) is the Linke turbidity factor for an air mass equal to 2, m is the relative optical air mass and  R (m) is the integral Rayleigh optical thickness. The air mass m is given by (their Equation # 3):

m( S ) = (p/p 0 ) / [ sin  S + 0.50572 ( S + 6.07995

) -1.6364 ] (2)
where  s is the solar altitude angle, possibly corrected for refraction effects, and the station height correction is given by their Equation # 4: (p/p 0 ) = exp (-z / z h ), where z is the site elevation and z h is the scale height of the Rayleigh atmosphere near the Earth surface, equal to 8434.5 meters. The parameterization of the Rayleigh optical thickness is given by their Equation # 7:

if m  20, ( S  1.9°), 1/ R (m) = 6.6296 + 1.7513m -0.1202m 2 +0.0065m 3 -0.00013m 4 (3) if m > 20, ( S < 1.9°), 1/ R (m) = 10.4 + 0.718m
In the modified ESRA model, the equations become: A VIRER !!!!!! m( S ) = 1 / [ sin  S + 0.50572 ( S + 6.07995) -1.6364 ] (4)

Let corr_ R (p/p 0 ) be the correction of the integral Rayleigh optical thickness due to the elevation of the site. John corr_ R (0.75) = 1.248174 -0.011997 m( S ) + 0.00037 m 2 ( S )

(5) corr_ R (0.50) = 1.68219 -0.03059 m( S ) + 0.00089 m 2 ( S )

Given that corr_ R (1)=1 and assuming that corr_ R (p/p 0 ) = corr_ R (0.5) for (p/p 0 )>0.5, corr_ R (p/p 0 ) can be determined for any (p/p 0 ) by piecewise linear interpolation. The integral Rayleigh optical thickness is thus given by:

if m  20, ( S  1.9°) 1/ R (m) = corr_ R (p/p 0 ) / [6.625928 + 1.92969m -0.170073m 2 +0.011517m 3 -0.000285m 4 ] if m > 20, ( S < 1.9°), 1/ R (m) = 10.4 + 0.718 m (p/p 0 ) (6)
The observed value G d (resp. G h ) is compared to this clear sky value G cd (resp. G ch ). The second condition is that G d should be less than (1.1 G cd ). Otherwise the observation for this day is flagged as questionable.

The clear-sky model makes use of the Linke turbidity factor, a parameter quantifying the atmospheric visibility (aerosol plus water vapor) under clear skies. Values of this factor are available on climatological basis [START_REF] Angles | A Web server for accessing a database on solar radiation parameters[END_REF][START_REF] Angles | A climatological database of the Linke turbidity factor[END_REF]; the SoDa web site: http://soda.jrc.it). It varies from site to site and from month to month. Typical values can be found in WMO (1981). A value of 1 stands for very optically clean atmosphere and may be adopted to control the upper limit of the daily irradiation. It is recommended to perform this test twice with two different Linke turbidity factors (e.g., 1 and 3 or 2 and 4) to assess the sensitivity of the results to this factor.

Quality control for minimum values of daily sums of global solar irradiation

The third condition to be fulfilled is that the daily sum G d (resp. hourly sum G h ) is greater than a minimum value expected in continuous overcast conditions. This is equivalent to saying that the clearness index KT d (resp. KT h ) should be greater than a minimum value. This minimum value was set up according to the analysis of the collected data and to the minima of KT d (resp. KT h ) found in the last and previous editions of the European solar radiation atlas (ESRA 1984[START_REF] Esra | European Solar Radiation Atlas[END_REF]. Accordingly, the smallest acceptable value for KT d and KT h is 0.03, which represents a heavily overcast sky.

The case of very high latitude sites

Special attention was paid to the observed solar irradiation data for the sites at very high latitudes, when the elevation of the sun above horizon at noon is less than 2°. When the extraterrestrial irradiation on the horizontal plane G 0d is greater than 1 J cm -2 (2.78 Wh m -2 ), the upper acceptance limit is set to (2 G cd ) and the lower acceptance limit is set to (0.015 G 0d ). When the extraterrestrial irradiation is less than 1 J cm -2 , values of G d are still accepted if they are less than 10 J cm -2 (27.78 Wh m -2 ) and no minimum value is set. The reasons for doing this are that on the one hand measuring instruments have an optimum precision of 5 W m -2 and that on the other hand atmospheric phenomena like refraction are present [START_REF] Esra | European Solar Radiation Atlas[END_REF].

Web service

The quality control procedure is part of an on-going effort of the Group 'Télédétection & Modélisation' of the Ecole des Mines de Paris to provide free of charge valuable tools and information to the solar radiation community through the most know and used media, namely the world wide web.

The web server HelioClim (http://www.helioclim.net) hosts the quality control procedure and other solar radiation oriented sizing tools. It is built from open source software. Apart the fact that those software can be used free of charge, they are rated among the most used, reliable, fast and secure tools available on the market (Netcraft 2001). Installed on a Linux box an Apache web server is delivering result by the mean of CGI (Common Gateway Interface) based scripts (see references on tools Linux and Apache). The interesting feature of the CGI-based quality control procedure is the re-use of source code (mainly C code), procedures and libraries already written, used and validated by the scientific community in various projects like [START_REF] Esra | European Solar Radiation Atlas[END_REF], Helioserve [START_REF] Angles | A Web server for accessing a database on solar radiation parameters[END_REF] and SoDa [START_REF] Wald | SoDa: a project for the integration and exploitation of networked solar radiation databases[END_REF]. The HelioClim web site use also extensively MySQL a relational database to store, extract and manipulate solar radiation parameters and PHP, a server-side, cross-platform, HTML embedded scripting language (see references on tools MySQL and PHP).

The user interface has been kept simple for a maximum reliability over the network bandwidth. The user is requested to provide information to compute the quality control procedure: geographical co-ordinates, elevation and dates (if single value). Detailed help HTML pages are available to better understand and correctly fill the forms of each quality control procedure. Also in line are a document explaining the algorithm used in the calculation and references to articles, web site of interest and publications on solar radiation topic.

The quality control procedure has been divided in four HTML documents:

1. Daily irradiation for a single day

Example

This web service was used to perform a quality screening of time-series of daily irradiation for the period 1994-1998 and for sites in Africa. Data were provided by the networks belonging to the World Meteorological Organization. These data are prominent in constructing global maps of the Linke turbidity factor within the SoDa project [START_REF] Wald | SoDa: a project for the integration and exploitation of networked solar radiation databases[END_REF] as well as in assessing the accuracy of the method Heliosat-II, aiming at retrieving solar irradiations from an appropriate processing of satellite data (Rigollier, Wald 1999). Hence, it is important that these data are properly screened. The number 10 means that the daily irradiation is greater than the extraterrestrial value. An example is the day 24 November 1994. When clicking on a cell flagged as incorrect, another page appears, similar to Figure 1b.

This new page gives the observed irradiation for this instant (day or hour), the extraterrestrial irradiation, the irradiation for clear sky and the solar elevation at noon. It permits to see how large is the discrepancy between the observed value and the threshold under concern. Values may be rehabilitated or discarded.

In August 1998, two days (13 and 21) exhibit values that are well exceeding the expected value for very clear skies. Here, the Linke turbidity factor was set to 2, a value well below the typical values for this site. It follows that the irradiations for these two days are overestimated and should be rejected. They are flagged with the number 11.

Missing values are set to 0. Such days (e.g., 31 November 1994, 29, 30 and31 February 1998) are flagged with the number 12.

The other flagging values 21 to 24 are used for the tests for low solar elevations at noon.

Conclusion

This web service is available to everyone on a free basis. Data are input in the HTML page by a copy and paste procedure. The return is also a HTML page that can be saved and further processed if necessary. This service calls upon up-to-date algorithms that have proven quality. They may be changed as progresses are made. This service offers a very convenient mean to users to check time-series of irradiation without writing and testing themselves the procedures.

Acknowledgements

The content of the present paper was influenced by the chapter II of the ESRA handbook.

2. Daily irradiation for several days 3 .

 3 Hourly irradiation for a single hour 4. Hourly irradiation for several hours Simple forms are used inside the HTML document allowing a one-click per page result. Figure 1 exhibits a snapshot of the input and output forms for screening a single day. The service was tested by the means of several time-series of hourly and daily irradiations available in the CD-ROM of the European solar radiation atlas (ESRA 2000).

Figure 1 .

 1 Figure 1. Screening daily irradiation. a) the form for screening a single day. b) the form of the screening results.

Figure 2

 2 Figure 2 exhibits an extract of the HTML page containing the results of the test. The site under concern is Casablanca (Senegal). Here are only reported three months of 1994 and nine of 1998 for illustration. Each line represents a month. For each line, there are 31 cells, one per day. The character "V" stands for "verified" and means that the value is OK with respect to the procedure. The days not passing the test are flagged with a number taking values 10, 11, 12 or 21 to 24.

  Fruitful discussions were held with John Page (Sheffield, United Kingdom) and Gerhard Czeplak (Hamburg, Germany). The program Leonardo of the Commission of the European Communities offered a grant for the stay of Michael Geiger in France. This work was also partly supported by the Commission under Contract Number IST-1999-12245.