Ngoc-Yen-Phuong Cao
email: ngoc-yen-phuong.cao@ifp.fr

Benoit Celse
email: benoit.celse@ifpen.fr

Denis Guillaume
email: denis.guillaume@ifpen.fr

Isabelle Guibard
email: isabelle.guibard@ifpen.fr

Joris W Thybaut
email: joris.thybaut@ugent.be

Stabilization time modeling for hydroprocessing: identification of the dominant factors

Keywords: hydrocarbon, hydrotreating, stabilization, transient data

Hydroprocessing stabilization has been assessed using experimental data acquired at transient conditions. These data were obtained from a hydrotreating pilot plant in a wide range of operating conditions. It has been found that the stabilization evolution follows a first-order response and could be characterized by a stabilization time τ. A linear model was developed to relate τ to its most influential parameters. The model can then be combined with online transient data to predict the steady-state performance. By testing against new data with other feedstocks, the model has been found to provide a good prediction of the stabilization evolution and the ultimate steady-state hydrotreating performance. It is, hence, possible to "online" calculate the steady-state performance without the need to reach this steady state.

Introduction

Hydrotreating is a catalytic conversion process in petroleum refining, among others for removing impurities such as nitrogen and sulfur from hydrocarbon streams. It eliminates such heteroatoms in oil fractions from compounds by mixing them with hydrogen and treating this mixture in a fixed bed catalytic reactor at high temperature and pressure. The catalyst used in hydrotreating usually is molybdenum (group VI metal) based and promoted with cobalt or nickel (group VIII metal) and supported on alumina or silica-alumina. The most common catalysts used in hydrotreating are CoMo/Al 2 O 3 and NiMo/Al 2 O 3 . CoMo favors the sulfur removal and olefin saturation while NiMo is used for removing nitrogen and saturating aromatic compounds (Treese et al., 2015). Hydrotreating/hydrocracking catalysts require several days before achieving stable operation while the residence time in the reactor is of the order of maximum a few hours. The observed transient phenomena can, hence, at most partly be attributed to hydrodynamics phenomena in the equipment used and mainly to chemical phenomena occurring at the catalyst level. These are distinct from and, hence, are not to be confused with the long-term catalyst deactivation. There is very few literature on this subject. [START_REF] Sau | Effects of organic nitrogen compounds on hydrotreating and hydrocracking reactions[END_REF] investigated the effects of organic nitrogen on hydroprocessing reactions. The authors did some hydrocracking experiments using a zeolite-based catalyst and determined the evolution of the conversion with time. When the reactor temperature was increased from 390 to 405 °C while processing a hydrotreated vacuum gas oil feed containing 320 ppmw organic nitrogen, a slow stabilization of the conversion was observed. The phenomenon was explained by the slow rate of nitrogen desorption from the catalyst with increasing temperature. Thanks to this nitrogen desorption, the concentration of active acid sites increases and, hence does the conversion. As the nitrogen desorption is slow, the observed conversion takes long time to reach the steady state as well.

Similarly, these authors also investigated the effect of increasing nitrogen content in the feed from 5 to 125 ppmw. The observed time to reach the stabilization of conversion is also around eight days and is attributed to the adsorption of nitrogen on the active sites of the catalyst.

According to [START_REF] Elizalde | Dynamic modeling and simulation of a bench-scale reactor for the hydrocracking of heavy oil by using the continuous kinetic lumping approach[END_REF], the Liquid Hourly Space Velocities (LHSV) has an impact on the dynamic behavior of hydrocracking process. LHSV is the ratio of liquid volumetric flowrate and the catalyst volume. It is the inverse of the reactor space time. A lower LHSV is equivalent with a higher space time, so the time to reach the steady state is longer because the product concentration needs more time to reach the steady state. However, the authors reported that the temperature in the range of 380-400 °C does not seem to significantly affect the dynamic behavior. It means that no matter what the exact reactor temperature is, the time to reach the steady state remains similar.

A kinetic model is a significant asset, not to say essential, for the adequate design and simulation of processes, especially of large-scale operations with narrow profit margins. Since crude oil contains a lot of compounds with difficult structures, kinetic modeling of hydroprocessing is a challenging and time-consuming task [START_REF] Jarullah | Kinetic model development and simulation of simultaneous hydrodenitrogenation and hydrodemetallization of crude oil in trickle bed reactor[END_REF]. There are various approaches for hydroprocessing modeling, such as the detailed kinetic [START_REF] Schweitzer | A single events kinetic model for the hydrocracking of paraffins in a three-phase reactor[END_REF][START_REF] Oyekunle | Kinetic modeling of hydrodenitrogenation of pyridine[END_REF][START_REF] Charon-Revellin | Kinetic modeling of vacuum gas oil hydrotreatment using a molecular reconstruction approach[END_REF][START_REF] Nguyen | Kinetic modeling of quinoline hydrodenitrogenation over a NiMo(P)/Al2O3 catalyst in a batch reactor[END_REF][START_REF] Raghuveer | Pyridine hydrodenitrogenation kinetics over a sulphided NiMo/γ-Al2O3 catalyst[END_REF][START_REF] Nguyen | Indole hydrodenitrogenation over alumina and silica-alumina-supported sulfide catalysts-Comparison with quinoline[END_REF][START_REF] Doukeh | Kinetics of thiophene hydrodesulfurization over a supported Mo-Co-Ni catalyst[END_REF], lumping modeling [START_REF] Bonnardot | Kinetic modelling of hydro-treatment reactions by study of different chemical groups[END_REF][START_REF] García | In-depth modeling of gas oil hydrotreating: From feedstock reconstruction to reactor stability analysis[END_REF][START_REF] Lababidi | Modeling the hydrocracking kinetics of atmospheric residue in hydrotreating processes by the continuous lumping approach[END_REF][START_REF] Tang | Lumping kinetics of hydrodesulfurization and hydrodenitrogenation of the middle distillate from chinese shale oil[END_REF][START_REF] Becker | Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach[END_REF][START_REF] Esmaeel | 5-Lumps kinetic modeling, simulation and optimization for hydrotreating of atmospheric crude oil residue[END_REF] and black-box approach [START_REF] Elkamel | Modeling the hydrocracking process using artificial neural networks[END_REF][START_REF] Bahmani | Product Yields Prediction of Tehran Refinery Hydrocracking Unit Using Artificial Neural Networks[END_REF][START_REF] Sadighi | Comparison of kinetic-based and artificial neural network modeling methods for a pilot scale vacuum gas oil hydrocracking reactor[END_REF][START_REF] Sadighi | Evaluating the ability of R for modeling a commercial scale VGO hydrocracking plant using Artificial Neural Network (ANN)[END_REF]. The most common and widely used is the lumping approach which consists of regrouping chemical compounds with similar properties [START_REF] Oliveira | A Review of Kinetic Modeling Methodologies for Complex Processes[END_REF]. In these approaches, the model parameters are generally estimated by fitting the model to steady-state experimental data. Stabilization phenomena lead to long experimentation times in order to obtain sufficient steady-state data for kinetic modeling.

However, in the transient regime towards this steady state, the analyses of liquid effluent are already carried out at regular time intervals to verify whether the steady state has effectively been reached and to ensure that the reaction is under control. These transient data are currently not used for kinetic modeling because the stabilization behavior is not well understood. It leads to a true challenge, meaning that experimental data acquisition to 'calibrate' the model is a time and money consuming task.

The aim of this work is, first, to acquire a better understanding of the stabilization behavior during these transient conditions and secondly, to use these data in order to predict, from the first transient points, the steady-state performance. If this value is far from a target, the operators can change the operating condition without waiting for stabilization and without the use of a complex model.

Materials and methods

Pilot plant

The experimental data are acquired using a pilot plant located at IFP Energies Nouvelles, Solaize, France. It consists of four parallel fixed beds, operated in down-flow mode, which are used for hydrotreatment of feeds varying from gas oil to vacuum gas oil. The total catalyst volume in each reactor amounts to 50 cm 3 . The pilot plant operates at isothermal conditions with the temperature being controlled along the reactor. Once the unit is in continuous operation, the density, nitrogen content and refractive index of the liquid effluent are analyzed every 24 hours.

The steady state is ensured by having observed the stabilization of these characteristics.

Operating conditions are adjusted after having reached the steady state corresponding with the previous operating conditions.

Data representation

The acquired data cover 11 Vacuum Gas Oil (VGO) feeds over two similar catalysts. Liquid Hourly Space Velocities (LHSV) vary between 0.5 and 4 h -1 , the temperature from 350 to 410 °C and the total pressure between 50 and 140 bar. Performance data are measured in terms of different properties of the liquid effluent such as density, refractive index and nitrogen contents over time on stream, totaling 920 measurements. A series of consecutive points corresponding to one experimental run is denoted as an 'episode'. The latter characterizes the intended steady state kinetic analyses rather than the systematic, c.q., repeated, analysis of transient behavior. Figure 1 shows 7 episodes corresponding with 42 data points acquired in about 45 days. The first from these 7 episodes corresponds to the start of the test at specific operating conditions, while later episodes are initiated by a change of operating conditions such as LHSV or temperature. Other variables that may also change between consecutive episodes are pressure and feedstock composition. It is noted that the deactivation phenomena in the pilot plant are negligible for the test with duration less than 45 days.

Stabilization modeling

A computational study on the pilot plant has been carried out using a tracer technique to estimate the stabilization time purely due to hydrodynamic phenomena. The objective is to identify whether the experimentally observed stabilization time is determined by hydrodynamic or chemical phenomena. The corresponding dynamic simulation of the pilot plant (without chemical reaction) was carried out using Matlab & Simulink R2014b in order to follow the flow of an unreactive tracer component through the pilot plant.

Exploratory data analysis [START_REF] Tukey | Exploratory Data Analysis[END_REF] (see also Figure 1) shows that the stabilization phenomena follow a first-order transfer function as shown in Equation (1).

𝑦 𝑦 𝑦 𝑦 1 𝑒𝑥𝑝 𝑇𝑂𝑆 𝑇𝑂𝑆  (1)
Where y is the liquid effluent property at a specific time on stream (such as nitrogen content, density, refractive index); y init is the liquid effluent property corresponding to the first experimental point of episode (ppm); y final is the last experimental point of episode (ppm); TOS is the time on stream (h); TOS init is the time on stream corresponding to y init (h) and τ is the characteristic time of the episode (h). The characteristic time τ of each episode presented in the equation is estimated via nonlinear least-squares technique using Gauss-Newton algorithm [START_REF] Bates | Nonlinear Regression Analysis and Its Applications[END_REF].

The estimated characteristic time τ can logically be supposed to depend on some variables such as feed properties and operating conditions. A linear function accounting for possible interaction among the variables is defined to describe the characteristic time τ (Equation ( 2)).

𝜏 𝑎 𝑥 𝑎 𝑥 𝑎 𝑥 𝑥 ⋯ 𝑏

(2)

Where x 1 , x 2 are the input variables, which can be LHSV, temperature, pressure or feed properties such as organic nitrogen content, organic sulfur content, resin content in feed, etc.;

x 1 x 2 is the interaction term between x 1 and x 2 ; a 1 , a 2 , a 12 ,… and b are the coefficients.

The most influential input variables of the model are determined via variable selection technique called 'leaps' [START_REF] Furnival | Regressions by Leaps and Bounds[END_REF]. The main idea of this technique is based on an efficient branch-and-bound algorithm which was first proposed by [START_REF] Land | An Automatic Method of Solving Discrete Programming Problems[END_REF]. It searches for all possible subset solutions and stores them in the branches of a tree with the full set at the root. The algorithm runs an exhaustive search and shows the best variable subsets. It reveals the linear dependence of stabilization on the input variables. The quality of each subset is evaluated by a metric called adjusted R 2 . In terms of statistics, higher adjusted R 2 means a better quality. However, the number of variables in the model should not to be excessive. The selected model is a compromise between adjusted R 2 and the number of variables in the subset.

The adjusted R 2 , see Equation ( 3), is a measure of model fitting quality taking into account the number of variables. It is usually employed to compare multiple linear models.

𝑅 1 1 𝑅 𝑛 1 𝑛 𝑝 1 (3)
Where R 2 is the 'coefficient of determination' or multiple R 2 ; n is the number of observations in data sample; p is the number of predictors, c.q., variables and/or parameters, in the model. Multiple R 2 is calculated via Equation (4).

𝑅 1 ∑ 𝑦 𝑦 ∑ 𝑦 𝑦 (4)
Where 𝑦 are observed values; 𝑦 are predicted values; 𝑦 are mean of sample.

Stabilization prediction

One of the properties of Equation ( 1) is that 95% of the value at steady state is reached after 3τ and 99% of the steady-state value after 5τ. It is a characteristic of first-order function which can be found in the literature [START_REF] Smith | A First Course in Differential Equations, Modeling, and Simulation[END_REF]. The time required to reach steady state can be estimated thanks to the τ predicted by the model.

The model for τ can then be used for stabilization performance prediction. For each episode, τ is calculated by the model. Knowing τ and the first two experimental points of the episode suffices to predict the episode evolution and, ultimately, the behavior at steady state using Equation ( 5), exemplified for the nitrogen content. If the predicted steady-state value is far from the target, the operating condition can be adjusted without waiting for stabilization. The prediction can be updated daily based on the new measurements.

𝑁 𝑁 𝑁

1 𝑒𝑥𝑝 𝑇𝑂𝑆 𝑇𝑂𝑆 𝜏

𝑁

(5)

Where N final is the liquid product nitrogen at steady state (ppm); N 2 is the liquid product nitrogen corresponding to the second point of episode (ppm); N 1 is the liquid product nitrogen corresponding to the first point of episode (ppm); TOS 2 is the time on stream corresponding to the second point of episode (h); TOS 1 is the time on stream corresponding to the first point of episode (h) and τ is the characteristic time predicted by model (h).

Results and discussion

This section exhibits the results of each step in the methodology shown in Figure 2. First part represents the hydrodynamic study of the pilot plant where the experimental data were acquired.

The analysis of stabilization behavior and its dependence on operating conditions, feedstock are displayed. Stabilization model is then constructed and the prediction against new data is performed.

Figure 2. Methodology adopted throughout this work.

Pilot plant hydrodynamics

A computational analysis of the pilot plant hydrodynamics has been carried out using Matlab & Simulink via the simulation of the behavior of an unreactive tracer component running through the pilot plant. The outlet concentration as a function of time provides insight into the stabilization behavior purely due to hydrodynamics. The pilot plant was simulated as several blocks representing each part of the pilot plant such as a heater, a reactor, a separator to remove hydrogen-rich gas from liquid stream and a stripping column to remove gas dissolved in the liquid stream such as H 2 S, NH 3 . Each element was modeled by considering it as a plug flow reactor (PFR) or as a continuous stirred-tank reactor (CSTR). Pipeline, heater and reactor were modeled as PFR while separator and stripping column were simulated as CSTR. Mass balance was then written for each block. In the simulation, the input stream is a tracer solution with a concentration set randomly at 1 mol/m 3 . The volumetric flowrate of the input stream was simulated at 50 cm 3 /h, which is similar to the one in the experimental data that would be used to compare. Temperature and pressure in the reactor were set at 370 °C and 140 bar respectively.

There was no simulation of the chemical reaction expected.

Exploratory data analysis

Stabilization was found to exhibit a first-order behavior, which is similar to the observations by [START_REF] Sau | Effects of organic nitrogen compounds on hydrotreating and hydrocracking reactions[END_REF]. The first-order characteristic time τ can be calculated from the evolution of effluent characteristics such as the nitrogen content, density and/or refractive index. Figure 4 compares the transient behavior of these properties, normalized between [0, 1], for the same experimental run. The stabilization behavior of these properties is similar and the values of τ are essentially the same. Following the evolution of stabilization by density and refractive index are more suitable than by the nitrogen content for the following reasons: (1) for episodes with low nitrogen level (< 10 ppm), it becomes difficult to observe the behavior in the nitrogen curve, while it remains feasible for density and refractive index; (2) density and refractive index measurements are typically more precise than for nitrogen and (3) density and refractive index measurements are easier (and less expensive) to carry out than nitrogen measurements. bar (τ = 60 h) is more rapid than the curve of 140 bar (τ = 106 h). For the later episodes, the difference in stabilization behavior between 115 and 140 bar seems to be less pronounced.

The impact of LHSV was analyzed using data from two tests with the same feed and identical pressure and temperature (P = 140 bar, T = 370 °C) with LHSVs amounting to 1 and 4 h -1 . As can be seen in Figure 5(c), a higher LHSV leads to a much quicker stabilization. The value of τ in the case of LHSV of 4 h -1 was 9 h which is lower than 65 h in the case of 1 h -1 .

The feedstock is also found to have an impact on the stabilization. Figure 5(d) compares the evolution of liquid effluent density as a function of time on stream between two feedstocks. The evolution with F3884 having a higher resin and nitrogen content stabilizes quicker than the one with F2588.

Other factors, such as temperature, catalyst were found to exert a less pronounced impact on the stabilization behavior.

Model for τ prediction

τ is modeled using linear regression accounting for interactions. As indicated before, two linear models for τ were built; one for first episodes and another for other episodes.

First episodes model

There are in total 21 values of characteristic time τ corresponding to first episodes. The possible input variables are operating conditions (LHSV, temperature, total pressure) and feed properties (organic nitrogen content N feed , organic sulfur content S feed , resin of feed res feed , density d feed and the weighted average temperature of simulated distillation TMP feed ). All input variables and the interaction terms which are products of operating conditions and feed properties are used in variable selection procedure. Figure 6 shows the evolution of the adjusted R 2 of the 5 best subsets of each size as a function of the size of subsets. It seems that 4 variables represents a good compromise between model accuracy and overfitting (bias/variance compromise), i.e., although the adjusted R 2 still slightly increases when including more than 4 variables, the increase seems so marginal that it is preferred to limit the number of variables in the model to 4.

Table 1 presents an overview of the selected variables in the best subsets with the number of variables in subset from 1 to 7. The 5 best subsets with 4 variables are detailed as well. Among these 5 best subsets, the subset (LHSV, P, res feed , LHSV*res feed ) giving the highest adjusted R 2 is selected. This model contains only 1 interaction term LHSV*res feed , which also appears in the 4 left. It is a robust and less complicated model.

The model presents an adjusted R 2 of 0.79 and a multiple R 2 of 0.83. The corresponding parity diagram is given in Figure 7. The results are quite good for all the feedstocks. Figure 6. Variable selection using leaps for linear model with interaction (first episodes). The coefficients of the model as well as the statistical values are presented in Table 2. The global p value amounts to 5.3×10 -6 and the F value for the global significance of the regression exceeds the 95% quantile by far, which both indicate that the model is statistically significant.

All variables are individually significant as well as their p value is sufficiently low. As can be seen, a negative correlation was determined between LHSV and τ. The modeling confirmed that temperature is not a dominant factor. These two results are coherent with those acquired by [START_REF] Elizalde | Dynamic modeling and simulation of a bench-scale reactor for the hydrocracking of heavy oil by using the continuous kinetic lumping approach[END_REF]. A positive proportionality between pressure and τ was observed. However, the impact of pressure is less pronounced than that of other variables. The interaction term shows that the impact of LHSV on τ depends on the value of feed resin, i.e., the polar components with high molecular weight. An overview of the diagnostics for the multiple linear regression with interaction yielding the combination (LHSV, P, res feed , LHSV* res feed ) as best performing one, is displayed in Figure 8.

The analysis is as below:

 Residuals vs. fitted shows that there is linear trend.

 The normal Q-Q indicates that the residuals are almost normally distributed.



The assumption of equal variance is also valid by checking the scale-location.

 Residuals vs. Leverage show that there is no influential observation. Even the observation #13 is on the Cook's distance line (red dashed line on the Residual vs Leverage plot), it is not an influential observation. Because when this point is excluded from the regression, multiple R 2 changes from 0.8301 to 0.8335, which means a small impact.

Figure 8. Diagnostic plots for linear model with interaction (first episodes model).

Later episodes model

There are 33 observations of τ in the whole data set for later episodes. The same procedure for variable selection was applied. The input variables are the feed properties, operating conditions of the current episode as well as of the previous episode. Variables relating to the previous episode are, e.g., temperature and LHSV of this previous episode: T pre , LHSV pre ; and ratio of temperature and LHSV of the episode and the previous one. Only a selection of interaction terms was allowed based on the inspiration acquired by modeling the first episodes, to keep the number of considered variables and, hence, adjustable parameters, within reasonable constraints. Table 3 shows the 17 input variables taken into account in variable selection procedure. Figure 9 show the evolution of the adjusted R 2 of the best subsets as a function of the number of variables. On similar grounds as for the first episodes modeling (Section 3.3.1), for the later episodes regression 7 variables is determined to be a good compromise between accuracy and over fitting for these later episodes. Diagnostics of model for later episodes is discussed as below (see also Figure 11).

 Residuals vs. fitted shows that there is a linear trend.

 The normal Q-Q indicates a normal distribution.

 Scale-location shows that the assumption of equal variance is valid.

 Residuals vs. Leverage show no influential observation.

Figure 11. Diagnostic plots for linear model with interaction (later episodes model).

Prediction

The model was used to predict τ and the steady-state performance was calculated via Equation

(5) assuming that the first two points of the episode have already been measured. The prediction was carried out on a new experimentation with another feedstock operating on another pilot plant but the same catalyst.

First episodes

First episodes model is tested using new data from other feedstocks. Figure 12

Later episodes

Later episodes model is tested using new data from other feedstocks as well. Figure 13 compares the prediction (solid line) and the experimental data. Again, the first two points of episode which are used to estimate the steady-state value are shown as circles. The remaining data in the episode are represented by triangles. The model predicts well the evolution of episodes 2, 3, 5, 7 and 8. Since nitrogen values are low in episode 7 and 8, the prediction is considered not far from the experimental value, even if the evolution is not perfectly in line with the experimentally observed one. Moreover, during the first point of an episode, hydrodynamic effects still may have a certain impact, indicating that we should be careful with the reliability of this first point.

Regarding episode 4 and 6, the stabilization behavior did not follow the first-order response since the reactor temperature was slightly adjusted to reach an intended nitrogen level. The reactor temperature profile is shown in Figure 14. A target of 30 ppm was fixed for the value of liquid effluent nitrogen in episode 4. The reactor temperature was firstly set at 378.7 °C. The prediction using the first two points was 25 ppm (see solid line in episode 4, Figure 13). The reactor temperature was, hence, gradually decreased to 377.1 °C. The temperature profile was considered constant from the fourth point. The steady-state value prediction using the fourth and fifth point was 32.8 ppm, which is not far from the target (see dashed line in Figure 13). The third point was not used for prediction since the temperature was adjusted during that time.

Similar to episode 4, the liquid effluent nitrogen target of episode 6 was 50 ppm. The temperature was firstly set at 374.7 °C for episode 6 and the steady-state forecast using the first two points was 80 ppm. The temperature was gradually modified to 376.7 °C. The steady-state value prediction using the fourth and fifth point was 52 ppm, which was close to the target (see also dashed line in episode 6, Figure 13). As the liquid product nitrogen values seem to be sensitive to the reactor temperature value, it is needed to carefully control the temperature profile. As can be seen, the steady-state performance prediction depends not only on the predicted value of τ but also on the first two "online" measurements. This 'online' prediction is, Figure 15 shows the steady-state prediction based on the first two points and different values of predicted τ. In fact, a variation of ±15% and ±30% was added on τ predicted by model in order to study the impact of predicted τ on the steady-state prediction. As can be seen, the nitrogen value comes to its steady state at about 14 -20 ppm for ±30% variation of predicted τ. The steady-state value of nitrogen is 15.5 -18.5 ppm for ±15% variation. It is shown that a small variation of τ does not impact significantly on the prediction. Moreover, the model can be employed at the experimentation scheduling step. The predicted τ (in hours) can be used to estimate the order of magnitude of the time required to reach steady-state (in days). As explained in Section 2.4, the steady state can be considered to reach after about 3τ -5τ.

The proposed methodology can be applied to any reaction where the chemical stabilization requires significantly more time than the hydrodynamic stabilization and analysis time to forecast the evolution towards the steady state. The latter can help the operator to follow the trend of the stabilization behavior. It enhances the decision making in process operation to reach an intended target.

Conclusions

Stabilization is crucial for hydroprocessing. In this work, transient hydrotreating data were used to probe the details of this stabilization, which was found to follow a first-order response.

Two models, i.e., one for the first episodes, another for the later episodes, for stabilization time 'τ' prediction were constructed. It was found that a higher LHSV leads to a quicker stabilization.

The extent of the impact of LHSV on τ depends on the feed resin content. A proportional relationship between pressure and stabilization was found. Temperature is not a dominant factor.

Stabilization of later episodes depends not only on the feedstock and operating conditions but also on the operating conditions of the previous episode.

The availability of two initial measurements during stabilization typically suffices to predict the steady-state hydrotreating performance. Good prediction results were obtained, particularly for the first episodes. The stabilization of later episodes is more complicated to predict. The Table 1. Selected variables corresponding to each subset obtained from variable selection in Figure 6 (only show the best selected subsets of 'from 1 to 7' variables and 5 best subsets of 4 variables)

Index

Selected variables Adjusted R 2

Figure 1 .

 1 Figure 1. Data illustration -Liquid effluent nitrogen as a function of time on stream (TOS).

Figure 3

 3 Figure 3 compares the transient behavior of tracer concentration in pure hydrodynamics mode with the effluent properties in reactive mode. The figure on the left represents the tracer outlet concentration. The figure on the right shows the transient behavior of liquid product nitrogen in reactive mode. As can be seen, the hydrodynamic stabilization of the pilot plants occurs within 25 hours while the observed stabilization of the pilot plant in reactive mode requires around 200 hours. The hydrodynamic response to a tracer step reached steady state significantly faster than the stabilization of the hydrodenitrogenation behavior. Hence, stabilization is mainly driven by chemical phenomena.

Figure 3 .

 3 Figure 3. Calculated effluent concentration profiles of tracer (a) and experimentally observed liquid effluent nitrogen content as a function of time on stream in the pilot plant (b) (LHSV = 1 h -1 , T = 370 °C, P = 140 bar).

Figure 4 .

 4 Figure 4. Comparison of transient behavior of nitrogen/density/refractive index (outputs are normalized between [0, 1]). Values of τ are 39.8, 36.3, 39.2 h calculated respectively from density, nitrogen content and refractive index evolution. There were 54 values of τ which have been classified into 21 first episodes and 33 later episodes, which are shown in Figure 5(a). It can be seen, indeed, that first episodes, on average, require more time to reach steady state as compared to later episodes. It was, hence, decided notto describe the τ of first and later episodes with a single model, but to develop dedicated models, i.e., one for the τ from first episodes and one for the τ of later episodes.

Figure 5 .

 5 Figure 5. (a) Classification of τ regarding the type of episode; (b, c, d) Impact of pressure (solid line: 115 bar, dashed: 140 bar, feed F2588); LHSV (feed F4384, P = 140 bar, T = 370 °C) and feedstock (LHSV = 1 h -1 , P = 140 bar, T = 370 °C) respectively on transient behavior of liquid effluent density.

Figure 5

 5 Figure 5(b), (c) and (d) show the transient behavior of liquid effluent density as a function of time on stream. Figure 5(b) displays the evolution of stabilization at two different pressures, i.e., 115 and 140 bar. It can be observed that for the first episodes, the evolution of the curve of 115

Figure 7 .

 7 Figure 7. Parity plot with 95% confidence interval (R 2 = 0.83) (first episodes model).

Figure 9 .

 9 Figure 9. Variable selection using leaps for later episodes model.

Figure 10 .

 10 Figure 10. Parity plot with 95% confidence interval (R 2 = 0.66) (later episodes model). The regression is globally statistically significant since the global p-value is sufficiently small and the F value exceeds the critical value. Table 4 displays the parameter values obtained for the model and their individual statistical significance. The p values for the individual variables indicate that S feed and LHSV pre are at the margin of statistical significance (P = 0.11 and P = 0.07, respectively) and that the other variables are significant. LHSV, sulfur and resin content in the feed represent the operating conditions and feed properties. Other variables represent the impact of previous conditions on the

  shows two examples of such predictions by the model. The first two points of the episode which are used to predict the steady-state value are shown as circles. The remaining data in the episode are represented by triangles. The predicted nitrogen values are very close to the measured ones. The model has been found to provide a good prediction of the stabilization evolution and the steadystate value. Using the model for τ prediction enables the calculation of the steady-state performance without the need to wait for reaching steady state. The advantage is that if the steady-state value could be far from target, operators can change the operating condition without waiting for stabilization.

Figure 12 .

 12 Figure 12. Two test predictions using the first two points of episode (points: experimental data, solid line: model prediction; Test (a): LHSV = 1.34 h -1 , P =132 bar, T = 370 °C; Test (b): LHSV = 1.71 h -1 , P =140 bar, T = 386 °C).

  hence, flexible and the operators should consider the certainty of experimental data and adapt to what really happen (modification of operating conditions, measurement error and equipment problem).

Figure 13 .

 13 Figure 13. Prediction with the first two points of other episodes (points: experimental data, solid line: model prediction using the first two points, dashed line: prediction while discarding the first three points).

Figure 14 .

 14 Figure 14. Reactor temperature profile corresponding to each episode

Figure 15 .

 15 Figure 15. Sensitivity of predicted τ by model on the steady-state prediction (points: experimental data; line: prediction based on the first two points and different values of τ)

  Figures & Tables

Figure 1 .

 1 Figure 1. Data illustration -Liquid effluent nitrogen as a function of time on stream (TOS).

Figure 2 .

 2 Figure 2. Methodology adopted throughout this work.

Figure 3 .

 3 Figure 3. Calculated effluent concentration profiles of tracer (a) and experimentally observed liquid effluent nitrogen content as a function of time on stream in the pilot plant (b) (LHSV = 1 h -1 , T = 370 °C, P = 140 bar).

Figure 4 .

 4 Figure 4. Comparison of transient behavior of nitrogen/density/refractive index (outputs are normalized between [0, 1]). Values of τ are 39.8, 36.3, 39.2 h calculated respectively from density, nitrogen content and refractive index evolution.

Figure 5 .

 5 Figure 5. (a) Classification of τ regarding the type of episode; (b, c, d) Impact of pressure (solid line: 115 bar, dashed: 140 bar, feed F2588); LHSV (feed F4384, P = 140 bar, T = 370 °C) and feedstock (LHSV = 1 h -1 , P = 140 bar, T = 370 °C) respectively on transient behavior of liquid effluent density.

Figure 6 .

 6 Figure 6. Variable selection using leaps for linear model with interaction (first episodes).

Figure 7 .

 7 Figure 7. Parity plot with 95% confidence interval (R 2 = 0.83) (first episode model).

Figure 8 .

 8 Figure 8. Diagnostic plots for linear model with interaction (first episodes model).

Figure 9 .

 9 Figure 9. Variable selection using leaps for later episode model.

Figure 10 .

 10 Figure 10. Parity plot with 95% confidence interval (R 2 = 0.66) (later episodes model).

Figure 11 .

 11 Figure 11. Diagnostic plots for linear model with interaction (later episodes model).

Figure 12 .

 12 Figure 12. Two test predictions using the first two points of episode (points: experimental data, solid line: model prediction; Test (a): LHSV = 1.34 h -1 , P =132 bar, T = 370 °C; Test (b): LHSV = 1.71 h -1 , P =140 bar, T = 386 °C).

Figure 13 .

 13 Figure 13. Prediction with the first two points of other episodes (points: experimental data, solid line: model prediction using the first two points, dashed line: prediction while discarding the first three points).

Figure 14 .

 14 Figure 14. Reactor temperature profile corresponding to each episode 568

Table 1 .

 1 Selected variables corresponding to each subset obtained from variable selection in

	Figure 6 (only show the best selected subsets of 'from 1 to 7' variables and 5 best subsets of 4
	variables)
	Index	Selected variables	Adjusted R 2
	1	P	0.369
	2	P, LHSV*S feed	0.528
	3	T, P, LHSV*S feed	0.625
	4a	LHSV, P, res feed , LHSV*res feed	0.788
	4b	LHSV, res feed , LHSV*res feed , P*TMP feed	0.782
	4c	res feed , LHSV*res feed , LHSV*TMP feed , P*TMP feed	0.780
	4d	LHSV, P, LHSV*res feed , T*res feed	0.763
	4e	LHSV, LHSV*res feed , T*res feed , P*TMP feed	0.756
	5	LHSV, T, P, res feed , LHSV*res feed	0.819
	6	P, res feed , LHSV*res feed , LHSV*TMP feed , T*res feed , P*N feed	0.835
	7	LHSV, P, N feed , LHSV*Res feed , T*TMP feed , T*N feed , P*N feed	0.828

Table 2 .

 2 Coefficient and statistical values for first episodes model

		Interce pt	LHSV	P	res feed	LHSV*res feed
	Coefficient	63.339	-69.055	0.693	-7.081	5.481
	p-value	3.33×10 -3	3.15×10 -6	7.1×10 -6	5.31×10 -5	1.9×10 -4
			Global p-value	F	F 95%
			5.3×10 -6	19.55	3.01

Table 3 .

 3 Variables used in variable selection technique for later episodes regression

		Operating
	Feed		Condition switching	Interaction terms
		conditions
	N feed	LHSV	T T	LHSV*res feed
	S feed	T	LHSV LHSV	LHSV*N feed
	res feed	P	T T	LHSV*S feed
	d feed	L H S V pre	LHSV LHSV
	TMP feed	T pre

Table 4 .

 4 Coefficient and statistical values for later episodes model

		Intercept		LHSV	S feed	LHSV pre
	Coefficient	177.349		-10.943	2.589	8.851
	p-value	8.32×10 -5		0.036	0.114	0.071
		T pre	LHSV*res feed	LHSV LHSV	LHSV LHSV
	Coefficient	-0.379		-0.828	8.832	-11.713
	p-value	8.34×10 -4		3.85×10 -3	5.78×10 -3	0.020
		Global p-value	F	F 95%
		1.4×10 -4		6.83	2.40

Table 2 .

 2 N feed , LHSV*Res feed , T*TMP feed , T*N feed , P*N feed 0.828 Coefficient and statistical values for first episodes model

	1	P

Table 3 .

 3 Variables used in variable selection technique for later episodes regression

		Operating
	Feed		Condition switching	Interaction terms
		conditions
	N feed	LHSV	T T	LHSV*res feed
	S feed	T	LHSV LHSV	LHSV*N feed
	res feed	P	T T	LHSV*S feed
	d feed	L H S V pre	LHSV LHSV
	TMP feed	T pre

TOS: time on stream (h) ppm: parts-per-million VGO: vacuum gas oil


Jean-Christophe Burnel

Kilian Fatras

Nicolas Courty

Generating Natural Adversarial Hyperspectral examples with a modified Wasserstein GAN

Keywords: Adversarial Examples, Remote sensing, Deep Learning

, it reweights the true data empirical distribution to encourage the classifier to generate adversarial examples. We provide a proof of concept of our method by generating adversarial hyperspectral signatures on a remote sensing dataset.

Introduction

One of the weakness of DNNs is the so called adversarial example sensibility. Adversarial examples are samples which are not classified in the correct class. These examples can occur even for very accurate classifiers and can be problematic for their real-world applications, especially when security is an issue. Therefore one may want to have access to examples which can fool a classifier in order to measure how robust is the DNN. Adversarial examples have received a lot of attention from the machine learning community. There are two main strategies to generate such examples. The first strategy is to use the gradients from classifiers. The generated adversarial examples are images drawn from the true distribution where we add a small perturbation following the gradients. The intuition is that the classifier will change its prediction because the gradients are the steepest classifier variation directions. This is the main idea of the Fast Gradient Sign Method (FGSM) algorithm [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF]. One problem from this approach is that one needs to know the full architecture to compute the classifier's gradients. The second paradigm does not require access to gradients. It looks for adversarial examples only by having access to the classifier's prediction. Our new method follows the second strategy as we detailed later.

To generate adversarial examples without requiring to classifier's gradients, [START_REF] Zhao | Generating natural adversarial examples[END_REF] uses a GAN and an invertor. The generator is a function which goes from the latent space to the true data space while the invertor is an inverse function. From the true data, they use the invertor to go to the latent space and look for adversarial examples. They Contactauthor:jean-christophe.burnel@irisa.fr add small perturbations to the latent vector z * which represents a data x * until they find an adversarial example. This method allows them to have very realistic adversarial examples without using the classifier's gradients. In [START_REF] Song | Constructing unrestricted adversarial examples with generative models[END_REF] authors use an AC-GAN to model the data distribution with the assumption that an ideal model could generate all the set of legitimate data. With such a model they search in the latent space for all the adversarial examples. Their search is done by minimizing the confidence score of a classifier while having the auxiliary classifier still predicting the correct class. So the adversarial examples look like true images and are adversarial for the attacked classifier.

We propose a different approach to generate natural adversarial examples from those methods. Like them, our method only needs to know the output of a given classifier, seen as it a "black box" classifier. From a re-weighted distribution of the true data based on a classifier's predictions, a generator is learnt to generate adversarial examples. The idea is to create a map between the latent space and the set of adversarial images which are present in the true data.

The paper is structured as follows: after a brief introduction on GAN, we present our main contributions in section 3. We will detail our modified loss function for the WGAN architecture which is specialised in adversarial examples generation. In section 4, we will test our method on hyperspectral data from the Data Fusion Contest 2018 [START_REF]IEEE GRSS Data Fusion Contest[END_REF], by generating data which are misclassified by a classifier.

Proposed methods

Gerative Adversarial Networks. The principle of a GAN is to learn a generator of realistic data. Intuitively, it tries to minimize the distance between the distributions of generated data and target data. This notion of distance between distribution is critical and has led to several variants of GAN. Formally, the generator takes z ∼ P Z , as input to generate data G(z) ∼ P G and then distribution their distance to real data distribution P r . P Z is usually a gaussian or a uniform distribution. GANs improve the generated data quality by minimizing the distance between the distributions P G and P r . However as we do not know the true distributions P G and P r , we only consider their empirical distributions from the available true data and the generated ones. Wasserstein GAN [START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF] is a GAN variant which uses the Wasserstein distance to minimize the distance between distributions. Thanks to the Kantorovich-Rubinstein duality [START_REF] Peyré | Computational optimal transport[END_REF] for calculating this Wasserstein distance, the objective function of WGAN is :

min θ max φ E x∼Pr [D φ (x)] -E z∼Pz [D φ (G θ (z))]
where D is within the set of 1-Lipschitz functions parameterized by θ. Analogous to GANs, we still call D the "discriminator" although it is actually a real-valued function. In order to respect the 1-Lipschitz constraint, [START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF] used a weight clipping trick for the discriminator's weights during the optimization procedure. However this practice is far from being optimal and might lead to unstability during optimization and poor minima. Another method [START_REF] Petzka | On the regularization of wasserstein GANs[END_REF] involves a gradient penalty and enforces the gradient to be less or equal to one. This is the variant of WGAN that we will consider in the rest of the paper:

min θ max φ E x∼Pr [D φ (x)] -E z∼Pz [D φ (G θ (z))] + E x∼P x [(max{0, ∇D φ (x) 2 -1}) 2 ]
where P x is the distribution of samples along the straight lines between a pair of points from P r and P G .

Generating adversarial data. For a fixed trained classifier on the true data, we aim to find a way to generate adversarial data with GANs. Instead of encouraging the generator to generate data close to the true data, we will encourage it to generate adversarial data which follows the true data distribution P r . Intuitively, the generator will be a map between the latent space and the adversarial data from the true distribution. We opted for a modification of the loss function of WGAN by considering a reweighting strategy of the empirical true distribution. Indeed, The generator generates data similar to the true data because of the true data empirical distribution which is 1 N n i=1 δ xi . Our main idea is to reweight a data weight according to the classifier's prediction. If a data is misclassified, the resulting weight will be bigger than a correctly classified data. The resulting distribution must have the following form: n i=1 p i δ xi , where:

n i=1 p i = 1.

Now we will describe and review the impact of different reweighting methods:

Uniform weight Here we are in the case of traditional WGAN, all data are given the same weight so there is no particular reason for the generator to produce adversarial examples, we can visualize it with figure 1a where we have a distribution where all points got the same weight 1/N.

Hard weighting An intuitive way to generate adversarial example is to only consider misclassified data. In figure 1b, we see the data given the following weighting : 1 -c(x) * x, where c(x) is the classifier output. The downside of this weighting strategy is that we are left with few data for very accurate classifier. Unfortunately, GANs are data hungry and then this method is not efficient to generate adversarial examples.

Soft weighting Another weighting strategy is to use the prediction vector from the classifier. This can be done through a softmax function, defined as:

softmax(x, w) exp(w * [x j -x max ]) K i=0 exp(w * [x i -x max ]) (1)
where x is our prediction vector and w is a temperature coefficient that controls the entropy of the rsulting distribution. In the figure 1c we can see an example of weighting using the softmax function with a temperature w = 8.

To demonstrate the relevance of the different reweighting methods, we use a toy dataset to illustrate our approach from Figures 1a to 1f, where we represent data from a certain class with points, the classifier's decision boundary with a line, so that points under the line are correctly classified and points over the line are misclassified.

As mentioned before, an important hyperparameter of our softmax function is the temprature coefficient. Its value allows controling the weight for data with low prediction score. In the next figure, we use a factor of 5 for figure 1d and we can see that while misclassified data are given bigger weights the major part of the points still have a decent amount of weight. In figure 1e the weights for data too far from the classification line start fading due to their weights being low and finally in figure 1f we are vastly using the points along the classification line.

Experiments and Results

Unlike traditional RGB images, hyperspectral imagery divides the color spectrum in multiple contiguous bands that can outreach the visible spectrum. In those images a pixel is a spectrum, and we can identify materials in a pixel by analysing it because each materials have its own spectral signature. To test our algorithm we first considered a pixel-based classification task on hyperspectral images. The classifier takes a spectrum as input and will output a probability vector of belonging to one particular material. The classifier we used is based on [START_REF] Hu | Deep convolutional neural networks for hyperspectral image classification[END_REF] and have an Overall Accuracy of 0.52 with κ = 0.43 on disjoint train/test, which seem coherent with the result observe in recent reviews [START_REF] Audebert | Deep learning for classification of hyperspectral data: A comparative review[END_REF].

Our objective here is to generate natural adversarial spectra. Dataset. We use the DFC2018 dataset [START_REF]IEEE GRSS Data Fusion Contest[END_REF] which was acquired over the University of Houston campus and its neighbourhoods. It has hyperspectral data covering a 380-1050 nm spectral range with 48 bands at a 1-m GSD. Here the classes do not only cover urban classes (buildings, cars, railways, etc.) but it also has some vegetation classes (healthy or stressed grass, artificial turf, etc.). A view of the data is given with figure 2, where we only take 3 of the 48 bands for visualization purpose, and the ground truth labels are shown with figure 3 Experimentation. For this experiment we want to generate spectra of healthy grass, we know that the classifier has an accuracy of 82% for that particular class on the whole dataset. For the WGAN architecture we choose 1-dimensional convolutions for both the generator and the critic as it keeps coherence between close spectral bands. For the generator we have the followings layers (as seen in figure 4):

(Linear, 384)(Reshape, (6,64))(Conv1D, (12,32), 3)(Conv1D, (24,16), 3)(Conv1D, (48,1), 5), and for the critic we have : We use all the available data. In this experiment we choose to use w = 20 Results. We consider the healthy grass for target class, and we want to generate spectra similar to it but misclassified by the classifier. For this we used our ground truth to replace the pixels of the target class by the closest generated spectra. To verify that our generated spectra belong to the target class, we choose to plot the statistics of the spectra, in plain the mean of the spectra and dotted the standard deviation. To interpret the results we choose to plot in the same figure the target class from the image, the adversarial examples generated and the two most predicted classes by our classifier. The laters are evergreen trees and stressed grass figure 5. In both cases we see that the adversarial spectra match the targeted class more than the predicted class. With figure 6 we first show classification of healthy grass from the original image, here the classifier have 88.77% accuracy, and the classification of the adversarial spectra.

To go further we show some results with two other classes with different accuracy. We first choose Cars class where the attacked classifier have an accuracy of 0.39 and Crosswalks where the classifier have an accuracy of only 0.05. The results are visible with figure 7, here the comparison is made with the most predicted class, and we can see that it works in both cases. Moreover we show in table 1 that the more accurate our classifier is the less there are adversarial examples per batch. This can be explained by the fact that we do not have access to an infinite number of data so a more accurate classifier leads to less adversarial data to train the model. Qualitative comparison with state of the art. We now compare our results with [9]1 . We adapted this method to the same considered architectures considered in this paper. We used untargeted attacks, which means that we want our spectra to be misclassified but we do not want it to be classified as a specific label. However, it is to note that their algorithm is capable of doing targeted attack. The results are visible with figure 8. We did not include class Car for it was not possible to obtain any adversarial example over 100 runs of 10 000 iterations. For the Healthy grass class, the produced adversarial spectrum exhibits a lot of noise. For the class Crosswalks, it looks like their algorithm was not able to produce spectra close the real distribution, while our method worked in both cases. This somehow illustrates the benefits of our generation scheme. Further studies will try to confirm this view on a more extensive validation, involving practitioners from the field.

Conclusion and discussion

In this work we propose a new method to generate natural adversarial examples for a pretrained classifier, using a modified loss function for the WGAN with a re-weighted distribution of our data based on the pretrained classifier's prediction. Our method departs from other generative networks in the sense that we do not generate adversarial samples from existing ones by addition of a small perturbation. We rather consider the case where we train a generator to produce adversarial samples. As an immediate consequence, our method can treat the initial pretrained classifier as a black-box, which is a more realistic scenario. We applied our method to the generation of hyperspectral data and show that it is possible to generate spectra that fool a classifier while keeping a similar content. One possible extension of our work would be to adapt our method to generate other modalities frequently encountered in remote sensing, such as point cloud data, but also to generate entire portions of images rather than solely focusing on spectra (spatial-spectral cubes). Here again, apart from the complexity of the needed neural architectures for such tasks, we do not see any major issue that would prevent our method from working. Complementarily, extensive tests with photo-interprets will be conducted to further assess the qualities of our algorithm.

Fig. 1 :

 1 Fig. 1: Visualization of differents weighting methods

Fig. 2 :

 2 Fig. 2: False RGB Fig. 3: Labels

Fig. 4 :

 4 Fig. 4: Architecture used for Generator and Critic

Fig. 5 :

 5 Fig. 5: [Best viewed in color] Comparison with Evergreen trees and Stressed grass

Fig. 8 :

 8 Fig. 8: [Best viewed in color] Results of [9] for classes Healthy grass and Crosswalks

Table 1 :

 1 Succes rate of our method (over 10 runs)

		Mean Std Clf. Acc
	Healthy Grass 0.66 0.06 0.82
	Car	0.81 0.05 0.39
	Crosswalks	0.99 0.009 0.05

following their online implementation : https://github.com/ermongroup/