Simulation experiments for maximising the availability of a commercial octene production facility

Overall availability of a chemical process is of critical importance in industry. In this paper we evaluate the process design factors that influence the availability of a new chemical production facility by performing computer experiments on a stochastic simulation model. Experimental designs commonly used in the Design and Analysis of Computer Experiments (DACE) and Classical Design of Experiments (DOE) are evaluated and compared for application by means of simulation experiments. Furthermore, response surface and kriging models are evaluated for the approximation of the input-output relationships. The most accurate experimental design by approximation model combination is used to explore the design space, both in terms of the overall availability and the percentage time offline. We illustrate how the design and analysis of simulation experiments (DASE) are used for minimizing the risks in the design of a new 1-octene production facility in terms of maximising the overall availability and minimizing the percentage time offline simultaneously.


Introduction
Higher alpha olefins, such as 1-octene, are used commercially as co-monomers for the production of, amongst other things, linear low-density polyethylene.The problem with conventional linear alpha olefin technologies is that, in addition to 1-octene, a mathematical distribution of less useful olefins is also produced.Consequently the production of 1-octene via the selective tetramerisation of ethylene is of great commercial value.Perhaps the best known homogeneous catalyst system for the tetramerisation of ethylene is that developed by Phillips Petroleum Company which consists of a chromium source, 2,5dimethylpyrrole and an alkylaluminium as activator.A number of other chromium-based homogeneous catalyst systems have also shown promise in this regard [3].
Recently, Sasol Technology engaged in the design of a new commercial ethylene tetramerisation process for the production of octene.In the tetramerisation process various secondary products are also produced.The most problematic of these is a polyethylene polymer which fouls the reactor system and necessitates periodic shutdowns for cleaning, and consequently results in production losses due to plant unavailability.
Cleaning occurs by means of a hot solvent wash process during which the polymer is melted or dissolved in a different solvent to the reaction solvent at 180 • C. If fouling is too severe to be removed by this process, the unit must be opened and mechanically cleaned by means of hydroblasting, which results in extended downtime.The formation of polymer can be reduced by activating the catalyst at high Aluminium (Al) concentration [24].Once the catalyst is fully activated, the reaction can continue at higher efficiency and at a lower Al concentration.This equates to a reduction in cocatalyst consumption, which is crucial in achieving an economic total catalyst package cost.
The base case process design is to activate the catalyst in a comparatively small upfront activation reactor and feed the product from this reactor containing fully activated catalyst to a train of two larger main reactors in series.After exiting the main reactor train at 60 • C, the reaction product is heated in the post reactor heater before undergoing two flashes in the devolatilisation (Devo) train to recover the unreacted ethylene for recycling and to separate the desired product from the residual polymer.Both the heat exchanger and the Devo units can foul with moderate to excessive polymer.
The devolatilised, polymer-free product is then sent to the product work-up train to recover the desired 1-octene.The work-up train is not expected to foul.The base case flow sheet for both the reaction and hot wash processes is shown in Figure 1 [27].It includes: 1. Three parallel, equally sized batch activation reactors, two of which are in operation while the third is being defouled.The product from each activation reactor in operation is sent to a main reactor train.Each activation reactor can feed either of the two main reactor trains (one at a time).A batch charge from the activation reactor to the main reactor occurs at constant intervals.2. Two equally sized main reactor trains, each consisting of two reactors in series.If one of the main reactors in series needs to be shut down, the other reactor in the train will be shut down and cleaned as well.3. A spare post-reactor heater plus devolatilisation unit.Each post reactor and devolatilisation train services the output of both main reactor trains.If a Devo unit fouls, the associated post reactor heater must also be shut down for cleaning and vice versa.4.There is one hot wash system which can wash one main reactor train, one activation reactor and one post reactor heater plus Devo unit simultaneously.It can wash any of these sections separately or all three sections together.Thus it can simultaneously wash sections of the process flow diagram that are in series, but not sections that are in parallel.
The commercial plant must be designed to compensate for the unavoidable fouling, shutdown and cleaning processes.The impact of the fouling and cleaning of key process units on the availability of the commercial plant must therefore be quantified in order to better assess the risks involved and to identify areas that may require design optimisation so as to improve overall plant availability.
An industrial plant is typically part of a larger industrial complex.The plant receives feed streams from other units in the complex, and some of the products are feeds for other units.The consequence of this interdependence of the different units has the implication that the plant is under an "obligation" to process a predetermined volume of feed, and to produce the corresponding products.Not achieving this throughput has an impact on the overall stability of the industrial complex.
The throughput of the industrial plant is directly related to the availability and size of the plant.It is therefore possible to compensate for availability by increasing the size of the plant.The cost of the plant is, however, directly related to the size of the components.
From a cost perspective it is therefore beneficial to design the plant to be as small as possible.Consequently it is of critical importance to have a reliable estimate of the availability of the plant.This enables the process engineers to make informed decisions in designing the plant.
Real-world systems, such as the one discussed in this paper, are so complex that analytical models of these systems are virtually impossible to solve mathematically.Numerical computer-based simulation can be used to imitate the behaviour of the system over time.
The system depicted by the flowsheet in Figure 1 was modelled in the stochastic simulation software Arena (of Rockwell Software [10]).The simulation model was validated against the industrial system to prove that the model is a true representation of the physical system.For the remainder of this article the output from the simulation model will therefore be assumed to be representative of and interchangeable with the expected output from the physical system.
An estimate of overall availability was obtained by calculating the average availability from the simulation output over a four year period.The percentage time offline was calculated similarly.These responses were used as approximations of the true availability of the commercial plant.The simulation model was used to evaluate the effects of different design variables and operating philosophies on the availability.These experiments would have been impractical or impossible to perform on a commercial plant.
Even though the simulation model makes it possible to evaluate the effects of the design variables and operating philosophies, performing these experiments without formal statistical design and analysis procedures can waste considerable human and computer time since stochastic simulation models are very computing-intensive models.To ensure the proper utilisation of the time and money invested in creating the simulation model, the methodology of the Design and Analysis of Simulation Experiments (DASE) [13,14,27] was adopted to extract the maximum amount of information from the model with the minimum number of runs.
Using the DASE methodology an approximation model of the simulation model can be constructed that will solve virtually instantaneously.This approximation model can be used in lieu of the original model to explore the entire design space efficiently, as illustrated in Figure 2. Given a set of s input variables x 1 , . . ., x s the simulation model yields n output variables y 1 , . . ., y n .An approximation model of the input-output relationships is then constructed.The goodness-of-fit of the approximation model is evaluated with an additional set of validation points.
The Design and Analysis of Computer Experiments (DACE) has received a lot of attention in recent years [28,30].The majority of the work has, however, been performed on deterministic computer models in contrast to simulation models.A simulation model is more similar to a physical experiment in the sense that the outcome of each experiment is a random value.Therefore, replications are required in order to estimate experimental error.The simulation model is, however, still a computer model, and the benefits of the work performed on design and analysis of computer experiments may be realised in simulation models as well.
Therefore, in this study experimental designs commonly used in the DACE and Classical Design of Experiments (DOE) are compared for application with simulation experiments.Different metamodels are also evaluated for the approximation of the input-output relationships.
The optimal design by metamodel combination was used to explore the design space, both in terms of the overall availability and the percentage time offline, and valuable insight was gained about the industrial system.The methodology of multiple response optimisation was then used to define an operating envelope subject to constraints for both the overall availability and percentage time offline.
The paper is organised as follows.First we present the process variables under investigation.In §3 the experimental methodology is discussed.§4 and §5 consist of the design of experiments, as well as the approximation models evaluated are introduced in context of simulation models respectively.The evaluation criteria are then discussed in §6.The results from the evaluation of different designs and approximation models are further discussed in §7.1.In §7.2 the application of the DASE methodology to the case study is demonstrated and discussed.Finally some conclusions follow in §8.

Process variables in the model
The system depicted by the flowsheet in Figure 1 was modelled in the simulation software Arena [10].To carry out a simulation using random inputs, such as time to failure and time to repair, the relevant probability distributions must be specified.Then, given that the input random variables to a simulation model follow particular distributions, the simulation proceeds through time by generating random values from these distributions [15] and uses these to evaluate various performance criteria.
The factors considered in the experimental study, together with their ranges, are depicted in Table 1.
The different probability distributions used to sample the time to failure and cleaning time are as follows.The Weibull distribution was used for the time to failure of the activation reactor, first main reactor in series and second main reactor in series.Exponential distributions were used for the time to failure for the post reactor heater plus Devo unit.Finally, for the cleaning by means of hot wash, a normal distribution was used.The selection of these distributions was guided by discussions with various process engineers.
In the case of the activation reactor or a main reactor fouling, it is assumed that, on average, one in thirty hot washes will be unsuccessful and that hydro-blasting will be required.
A normal probability distribution around the average is assumed.A constant turnaround time for hydro-blasting is assumed.This includes the time taken for attempting and abandoning the hot wash process.
In the case of the post reactor heater or Devo unit fouling, the same probability of a hot wash being unsuccessful and hydro-blasting being required is assumed.The reactors' failure distributions are specified as Weibull distributions.The Weibull distribution is defined by two parameters: α and β.The parameter α is a shape parameter.If α = 1, then the Weibull distribution coincides with the exponential distribution.The parameter β is a scale parameter.Possible applications of the Weibull distribution are the time to complete some task or the time to failure of a piece of equipment; it is used here as a crude model in the absence of data [15].A change in α generally alters a distribution's properties (e.g.skewness) more fundamentally than a change in location or scale parameter β [15].There are many other distributions available which could be used for the different units ( [15], Chapter 4).However, these distributions were applied with the current problem and yielded acceptable results.

Experimental methodology
The methodology that was used to compare the different experimental design by approximation model combinations is presented in Figure 3.The methodology is presented as a tree diagram with branches and nodes.Each experimental design and approximation model is presented as a node in the tree.Phase 2 was performed for each screening branch, but due to space constraints it is only shown for the D-Optimal branch.The approxima-tion models were fitted to each of the second phase designs.Therefore, the response surface and kriging models were fitted to 25 design combinations.The objective was to determine the most accurate model in the final stage of the tree diagram for predicting the availability of the plant.

Design of experiments
Computer-based simulation and analysis are used extensively in engineering to predict the performance of a system or product.These computer models can be either deterministic or random in nature.The design and analysis of (deterministic) computer experiments is currently undergoing extensive research [28,29,30].Stochastic simulation models are, however, random in nature.The DASE is discussed in [13].In this paper four designs from the DOE (fractional factorial designs, Plackett-Burman designs, central composite designs and D-optimal designs), and two space filling designs from the DACE literature (uniform designs and Latin hypercube designs) were compared for application in the DASE.
As the number of factors k increases, the number of runs in a 2 k factorial design required for a complete replicate of the design rapidly outgrows the resources of most experimenters.If the experimenter can reasonably assume that certain high-order interactions are negligible, information from the main effects and low-order interactions may be obtained by running only a fraction of the complete factorial experiment.These are denoted as 2 k−p designs, where p is the fraction (2 p fraction).These fractional factorial designs are among the most widely used types of designs for product and process design [22].
Plackett-Burman designs, attributed to [25], are two-level fractional factorial designs for studying k = N − 1 variables in N runs, where N is a multiple of 4. If N is a power of 2, these designs are identical to two-level fractional factorial designs.However, for N = 12, 20, 24, 28 and 36, Plackett-Burman designs are sometimes of interest.More specifically, these designs cannot be represented as cubes; they are sometimes called non-geometric designs ( [16], [22, p. 319]).
The D-optimal criterion minimises the generalised variance of the parameter estimates in a linear model [1].Specifically, it minimises the determinant of (X T X) −1 , where X is the expanded design matrix having one column for each coefficient to be estimated in the model.For any candidate design X * , the D-efficiency of X * is defined as |X * T X * |/|X T X| relative the optimal design X.
The central composite design (CCD) consists of three points: a 2 k−p two-level factorial or fractional factorial design, a set of 2k star points at a distance α = √ k from the centre of the design, and n 0 centre runs.Lucas [19] stated that the use of the central composite design is the routine 'production run' of the applied statistician.He also showed in [20] that composite designs perform well in terms of the D-efficiency measure.Due to its practical usage, and the flexibility it offers in design, the CCD will be considered for evaluating design efficiencies.DACE was developed from the methodology of the design and analysis of physical experiments.However, the theories of DOE [22] and of response surfaces [23] are based on the fact that an observation in a physical experiment is affected by variability due to the effect of a number of independent factors and random variability.In contrast, computer experiments are deterministic and yield the same result from repeated runs.Therefore, space filling experimental designs are employed for computer experiments.Latin hypercubes and uniform designs are examples of different types of space filling designs [30].
Since uniform designs were introduced in the early 1980s by Fang [5], it has become very popular.The design points of a uniform design are uniformly scattered on the experimental domain.It is a type of fractional factorial design with an added uniformity property.According to Fang [6], the uniform design is superior to other designs because many other design criteria are simultaneously optimised together with minimisation of the discrepancy criterion.A considerable advantage of uniform designs is that less information is required of the underlying model.A large number of different uniform designs are available at the website [7], and they are specified by the notation U n (q k ), where U denotes uniform design, n the number of runs, k the number of factors and q the number of levels.See [6,8] for the theory and application of uniform designs.

Screening experiments
The importance of factors depends on the experimental domain.The users should supply information on this domain, including realistic ranges of the individual factors and limits on the admissible factor combinations (e.g.some factor values must add up to 100%).Therefore, in practice, user involvement is crucial for the application of screening methods [13].
For this study the 14 variables and ranges are shown in Table 1.The following designs were evaluated for the screening phase.Designs with approximately 20 runs were chosen.
1. 2  (4 14 ) uniform design.4. Latin hypercube design with 20 runs.The design was created using the lhs package in R [26]. 5. D-optimal design with 15 runs.The design was created using Design Expert [32].
For each of the designs the simulation model was executed only once at each design point.Another approach would have been to replicate each design point, for instance, five times and to use the average response value of the designs.Only one replication was used because the goal of this study was specifically to assess the most efficient (smallest number of model runs) experimental design strategy for simulation experiments.During the screening phase there can literally be hundreds of factors, and adding replications on each design point can increase the computer time of the study prohibitively.
The outputs for all the designs are summarised in an asterisk were found to be statistically significant for each design from the analysis of variance [22, p. 60].
The following may be observed from Table 2.The order of highest impact of the first three values are mostly similar, but for the rest of the factors the designs differ.The number of significant variables also differ for all five of the designs in this screening study.For example, the D-optimal design yielded only 4 significant factors, compared to the Resolution III fractional factorial design which yielded 8 significant factors.Note the R 2 values, i.e. the percentage variance explained by the model, are very high for all the designs.

Interaction/second-order experiments
During this phase experimental designs were evaluated on the sets of variables which were identified as significant from the different screening designs (Phase One).Note that for each branch in Figure 3 different variables were found to be significant.Therefore, different variables and numbers of experiments were used for the second phase designs of the different branches.The designs evaluated during this phase are shown in Table 3. Twenty-five design combinations were investigated during this phase (See Figure 3).
Generally, during the second phase of a sequential design process, interaction designs are employed, and the data are checked for curvature [23].If curvature is present, a third phase follows where second-order models are fitted.In this study, the curvature was found to be insignificant for all the designs and models evaluated.A central composite design was, however, included in this phase, because it is a well-known and widely used design.
The results are discussed in §7.1.

Approximation models
The response surface and kriging models used in this study are outlined in this section.

Response surface models
Response surface methodology (RSM) is a collection of mathematical and statistical techniques useful for the modelling and analysis of problems in which a response of interest is influenced by several variables and the objective is to optimise the response [22].
Response surface modeling postulates a model of the form where y(x) is the unknown function of interest, and f (x) is a known polynomial function of x ∈ χ, where χ is the design space.Here is a random error, which is assumed to be normally distributed with zero mean and common variance σ 2 [22].
In most RSM problems, the form of the relationship between the response and the independent variables is unknown.Thus, the first step in RSM is to find a suitable approximation for the true functional relationship between y and the set of independent variables.Usually, a lower-order polynomial in some region of the independent variables is employed.
If the response (y) has a linear relationship with the independent variables (x 1 , . . ., x k ), then the approximating function is the first-order model If there is curvature in the system, then a polynomial of higher degree must be used, such as the second-order model Obviously, it is unlikely that a polynomial model will be a reasonable approximation of the true functional relationship over the entire space of the independent variables, but for a relatively small region they usually work quite well [22].
The parameters, β0 , βi , βii , and βij , i = j = 1, 2, . . ., k, of the polynomials in ( 2) and ( 3) are determined by means of least squares regression which minimises the sum of squares of the deviations of the predicted values, ŷ(x), from the actual values, y(x).The least squares estimates are obtained from where X is the model matrix of sample data points, X T is its transpose, and y is a column vector containing the values of the response at each sample point [22].
If the model is sufficiently accurate, it is applied to search the design space for improved or optimal solutions.Myers and Montgomery [23, Chapter 1, p. [10][11] give a detailed discussion on the phases of response surface modelling.

Kriging
Kriging is an interpolation technique that was originally developed in the field of geostatistics.Sacks et al. [28] initiated the application of kriging to DACE where it postulates a combination of a regression part and a stochastic part, Here x is a vector of design variables, y(x) is the unknown response function, f T (x)β the known (usually polynomial) regression function of x, β the p × 1 vector of parameters to be estimated, f (x) the vector containing the functions of x expanded for the polynomial model, and Z(x) the realisation of a stochastic process with zero mean, variance σ 2 , and non-zero covariance.While f T (x)β approximates the design space, Z(x) creates "localize" deviations or departures so that the kriging model interpolates the n t data points [9].
The covariance matrix of Z(x) is given by where R is the correlation matrix, and R(x i , x j ) the correlation function between any two of the sampled n t data points x i and x j .Here R is an n t × n t symmetric matrix with ones along the main diagonal.The correlation function R(x i , x j ) needs to be specified by the user.Several correlation functions may be used and are discussed in [30].Throughout this paper, the Gaussian correlation function is employed, where n d is the number of design variables, θ k is the unknown correlation parameters used to fit the model, and x ik and x jk are the k th components of sample data points x i and x j .
Predicted values of the response y(x) at untried values of x are given by ŷ where y(x) is the column vector of length n t containing the sample values of the response, and F (x) the n t × p matrix of design points.Here r T (x) is the correlation vector between an untried x and the sampled data points {x 1 , . . ., x nt }, and is defined as The vector β = (F (x is used in (8).The variance between the underlying model β and y is estimated as The maximum likelihood estimates (MLE's) of θ k in (7) used to fit the model are found by Any values for θ k will create an interpolative model but the "best" kriging model is found by solving the k-dimensional unconstrained non-linear optimisation problem in ( 12) [31].

Evaluation criteria
The comparisons between the different designs and approximation models were performed by sampling additional validation points to assess the accuracy over the region of interest.
For each set of validation points, the maximum absolute error (MAX) and the root mean square error (RMSE) were computed as MAX = max i∈{1,...,ne} and respectively, where n e is the number of additional validation points.While RMSE provides good estimates of the "global" error over the region of interest, MAX gives a good estimate of the "local" error by measuring the worst error within the region of interest.A good approximation will have low RMSE and MAX values [17].

Results and discussion
Two approximation models for the input-output relationship are compared in this section, namely response surface models (see §5.1) and kriging models (see §5.2).

Design and model comparison
For each of the designs in §4.2 a response surface model and kriging model were fitted to the data from the simulation experiments.The response surface model was fitted using Design Expert [32].The kriging model was fitted using the DACE toolbox [18] in Matlab [21].A twenty-run random Latin hypercube design was used as validation grid, created with the lhs [2] package in R [26].
For each approximation model in Phase 3 (Figure 3), the plant availabilities were predicted by means of the response surface and kriging models for the validation grid.The root mean square error ( 14) and maximum absolute error (13) values for the predictions were calculated.
The general trends for the root mean square error and maximum absolute error were found to be the same.Therefore, due to space constraints, only the results for the root mean square error are discussed here.
The calculated root mean square errors are shown in Figures 4 and 5 for the response surface model and kriging model results, respectively.The vertical axis represents the criterion values, and the horizontal axis represents the interaction designs.From Figures 4 and 5 it may be observed that the lowest criterion values for both the response surface models and the kriging models originate from the Plackett-Burman and maximin Latin hypercube screening designs.Using effective and efficient screening designs are of great practical importance.In an industrial experimental design study, the screening results are the building blocks for the follow-up experimental phases.If the incorrect factors are used in the follow-up phases, it may lead to misleading and suboptimal findings.
From Figures 4 and 5 it may also be observed that the lowest criterion values for the response surface models originate from the fractional factorial and D-optimal second phase designs for both the Plackett-Burman and maximin Latin hypercube screening branch.The central composite design also yields comparable results.For the kriging models it may be observed that the lowest criterion values are obtained from the uniform second phase designs for the Plackett-Burman screening branch, and all the second phase designs yield comparable results for the maximin Latin hypercube design.Therefore, it may be concluded that the space filling designs (the uniform design and the maximin Latin hypercube design) are more effective for fitting kriging models than for fitting response surface models.
To illustrate the consistency of the results we replicated the experiment for the Plackett-Burman screening design five times.More specifically, response surface models were constructed for each of the five replications for each of the Phase 2 designs.The validation grid was then predicted with each of the five response surface models.The results are summarised in Table 4.
The trend in the root mean square errors between the different designs is consistent with  that illustrated in Figure 4, although the confidence intervals do overlap.Furthermore, the narrow width of the confidence intervals illustrates that the comparisons for each of the interaction designs will yield the same conclusions, i.e. the fractional factorial design is the most accurate interaction design in the second phase.Therefore, replicating the experiment will result in the same conclusion regarding the designs and models.
Kriging models and linear regression models share a common mathematical framework consisting of regressors and errors, but the emphasis is quite different.Linear regression focuses entirely on the regressors and their coefficient estimates, and makes simplistic assumptions about the errors (independence).In contrast, kriging makes simplistic assumptions about the regressors (just a constant term) and focuses entirely on the correlation structure of the errors.Thus, regression and kriging are probably best thought of as diametric opposites.Regression is about estimating regression coefficients that (together with the assumed functional form) completely describe what the function is.Kriging is about estimating correlation parameters that describe how the function typically behaves.Kriging makes predictions by interpolating and extrapolating from the data in a way most consistent with the estimated typical behavior [9].
Overall the difference between using the response surface models ( §5.1) and the kriging models ( §5.2) are not as pronounced.It would, however, seem that the response surface models give slightly better results.This could possibly be due to the random nature of stochastic simulation models.Each design point of the simulation model is a sample from a population.Therefore, there exists uncertainty around the actual value of the response at the specific design point.As discussed earlier, kriging models are interpolation models and are built on the assumption that the values of the design points are deterministic, i.e. two replications of the design point should yield the same response value.
To conclude, the overall best branch of the tree in Figure 3 is the Plackett-Burman design for Phase 1, the fractional factorial design for Phase 2, and the response surface model for Phase 3, both in terms of the root mean square error ( 14) and the maximum absolute error (13) criteria.This result might be specific to the case study.However, it is provided as a general recommendation for the design and analysis of simulation experiments in industry.

Overall plant availability
Typical questions in a study such as the one presented here, centres around the risk of decisions.The engineers need to compare the risks associated with different scenarios.
Risk is extremely difficult to quantify if the design space is not well understood.The application of the response surface model in (15) for quantifying risk is discussed in this section.
One way to quantify the risk in designing for a desired availability is to calculate the area of the design space for which the availability is at least as large as the design availability over the total design space.For the current process it is not possible to calculate this area analytically, but it can be approximated by means of random sampling techniques.Specifically, uniform random data were generated for all the variables in the design space, and the availability was predicted by means of the approximation model at each point.More specifically, one million uniform random values were generated for each factor (Table 1).For the "Hotwash Parallel or Series" factor the values below 0.5 were combined to create "Hotwash Series", coded as 0, and values greater than 0.5 were combined to create "Hotwash Parallel" coded as 1.The response surface model in (15) was used to predict the plant availability for each of the factor combinations.
Let L denote the desired percentage availability.Then the probability that the availability for a given set of random input variables x, A(x), is greater than or equal to L% availability may be defined as where I is the indicator function and χ is the design space (Table 1).In other words, the number of values above a certain availability divided by the total number of simulations is an approximation of the percentage of design space that satisfies the criterion.The mean variable values which satisfies a given L% availability may be calculated as where The plot in Figure 6 was generated by identifying the number of data points that yield an availability greater than or equal to the specified availability (L%), specified on the horizontal axis, and the number of points that satisfy this criterion divided by the total during periods of plant unavailability.If the output from the operational plant is processed by another unit, the receiving unit will have to receive alternative feed, or source the feed from a buffer.If no such flexibility is available, these units will have to shut down during this period.Buffers are expensive equipment, and the cost is related to the size of the buffer.The percentage of time that the total plant is offline should therefore be minimised.The approximation model for the percentage time offline is given by where the variables have the same meanings and are encoded in the same way as before.
Let L denote the desired percentage of time offline.Then the probability that the time offline for a given set of random input variables x is less than or equal to L % is defined as The mean variable values which satisfies a given L % time offline value may be defined as where Figure 7 depicts the percentage of time offline below the specified time offline (L %) on the horizontal axis, and the percentage of design space that satisfies the criterion on the vertical axis.This gives an approximation of the percentage of the design space that is covered by the data for a percentage of time offline lower than or equal to the specified value.The means of the variables which satisfy the specified target are summarised in Table 6.For example, at an offline target of L = 1% the means of the factor values that yield percentage time offline values lower than or equal to 1% are depicted.The "Percentage of Total Design Space" column gives an approximation of the percentage of the design space that is covered at the specified percentage time offline target value.

Multiple response optimisation
For the process design it was required to find the region of design space where both the percentage time offline and overall availability are within an acceptable range.As an illustration, the following criteria are specified: • Percentage time offline between 0 and 1%, and • Total availability between 90 and 100%.One approach is to construct a graphical overlay plot.The optimised results depicted in Tables 5 and 6 were used to set the factor levels for the graphical overlay plot.Figure 8 depicts the overlay plot for the average time to hotwash and the minimum uptime for the main reactor satisfying the different criteria.Alternative methods such as the Derringer and Suich desirability function approach may be used for multiple response optimisation [4].
This plot may be used to aid design decisions regarding the average time to hotwash, given a minimum uptime for the main reactor.For instance, given a minimum uptime for the main reactor of 7 days, the target average time to hotwash must be approximately between 12 and 14 hours in order to satisfy the criteria.It is always beneficial to have a range of options available, as this allows for flexibility in deciding on the most practical and cost effective factor levels for the process design.

Conclusions
In this study experimental designs commonly used in the DACE and in DOE were compared for application by means of simulation experiments.Different metamodels were also evaluated for the approximation of the input-output relationships.
The optimal design by metamodel combination was used to explore the design space, both in terms of the overall availability and the percentage time offline, and valuable insight was gained about the industrial system.The models in (15) and (19) were then used to define an operating envelope subject to constraints on both the overall availability and percentage time offline.The Plackett-Burman and maximin Latin hypercube designs were found to be the best screening designs for the case study discussed; both these designs are easy and convenient designs to implement.Most modern statistical design packages should include these designs.Choosing between these two designs would mostly be a matter of personal preference and modelling considerations.A possible benefit of the Latin hypercube family of designs is that it can be created for any number of runs and factor levels.The Plackett-Burman design would, however, be beneficial if setting the factor levels in a simulation model is not trivial.For the Latin hypercube designs each design point is at a different factor level, whereas for the Plackett-Burman designs the factor levels are only varied at their upper and lower levels.In some simulations a large amount of modelling effort and coding is involved in setting the factors at each level of the design.For these kinds of models the Plackett-Burman design could be advantageous.
During the second phase of the design process the fractional factorial design, the D-optimal design and the central composite design yielded comparable results when fitting response surface models.If curvature is observed in the design space, the central composite design or second-order D-optimal design should be used.For the kriging model all the designs gave comparable results.An interesting observation is that the two space filling designs (the maximin Latin hypercube design and the uniform design) performed better with the kriging model compared to the response surface model.
Response surface models and kriging models differ in their ease of use.Kriging models are much more difficult to fit and to interpret.Interpreting the optimal θ values for the model is also very difficult.Furthermore, it would be problematic to release the model to a client, or use the model as 'n "black box" in a current simulation model.Response surface models are, however, easy to fit with most statistical packages.The coefficients are easy to interpret, and it is trivial to build the approximation model into a simulation model, or any other software to release to a client.
The response surface model constructed on the fractional factorial design was used to analyse the case study, and to make recommendations to the design engineers.The case study consists of 14 factors, and it would have been extremely difficult to develop a good understanding of such a complex system without the use of simulation and design of experiments.The resulting approximation model was used to explore the design space effectively and efficiently, and to gain valuable insight about the system under investigation.These results were used to decide on what factors to focus on in further experimental work in the pilot plant, and to make informed decisions about the sizing of the reactors.
The overall availability and percentage time offline criteria were evaluated simultaneously in order to identify a feasible operating region of the average time to hotwash and the minimum uptime for the main reactor which satisfies a dual criterion.This is a powerful technique which may be used to find an optimal compromise between multiple response criteria [4,11,12].Building a simulation model is a difficult and time consuming process.It is therefore only natural to expect the maximum amount of value to be gained from the model.Therefore, it is recommended that experimental design techniques must be used to enhance the quality and quantity of information and understanding that can be extracted from simulation models.This ensures that the returns on the time and money invested in building the model are maximised.As was demonstrated in the industrial case study in this paper, once the approximation models are obtained, they can be employed for design exploration and optimisation.An experimental design and analysis strategy should form a critical part of any simulation study, if possible.

Figure 1 :
Figure 1: Reactor flow sheet for an ethylene tetramerisation process.

Figure 2 :
Figure 2: Computer model and experiments for an industrial system.

Figure 3 :
Figure 3: Methodology used for comparing different design and approximation model combinations.

Figure 7 :
Figure 7: Percentage design space below percentage time offline targets.

Figure 8 :
Figure 8: The feasible area for percentage time offline between 0 and 1%, and availability between 90 and 100%.

Table 1 :
Input variables and ranges used in the experimental design process.
Hotwash Parallel or SeriesIndicates whether the hotwash can treat everything in parallel or only the equipment that are in series.This is a categorical variable, i.e. it can only have discrete values.Series is coded as 0, and parallel as 1. 0 1 [32]0Fractional factorial design.The effects of 14 factors are evaluated in only 16 runs.The design was created using Design Expert[32]. 2. Placket-Burman design with 20 runs.In this case study an N = 20 run Plackett- [13,14,15]ign in 19 factors was applied during the screening phase.This is the smallest Plackett-Burman design equal to or larger than 14 factors.After constructing the design in the 19 factors, columns 15 to 19 were removed (these columns are not required because there are only 14 factors present in this study).Plackett-Burman designs are often applied as screening designs in the design and analysis of simulation experiments[13,14,15].3. U 20

Table 2 .
The numbers in the table indicate the relative sizes of the absolute values of the effects of the factors[22, p. 68].This gives an indication of the relative importance of the factors.The values marked with

Table 2 :
Summarised results from the screening experiments.

Table 3 :
Table of designs used in the second phase experiments.Column names are phase 1 designs, and row names are phase 2 designs.In the row names the number in brackets denotes the number of significant variables, and in the table the number in brackets denotes the number of experiments performed.

Table 4 :
Summary of the root mean square error (RMSE) values of the response surface models for the Plackett-Burman screening design.

Table 5 :
Mean values of factors for different availability targets.

Table 6 :
Mean values of factors for different percentage time offline targets.