Phase Stability Analysis with Equations of State—A ... - ACS Publications

Article pubs.acs.org/IECR

Phase Stability Analysis with Equations of StateA Fresh Look from a Different Perspective Boyan B. Ivanov, Anatolii A. Galushko, and Roumiana P. Stateva* Institute of Chemical Engineering, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria S Supporting Information *

ABSTRACT: The paper introduces a novel approach to the solution of the phase stability problem under constant temperature and pressure with equations of state (EoSs) as the thermodynamic model. An effective strategy that is designed to consecutively transform the tangent-plane distance function (hyper)surface and guarantees a nonrepetitive determination of all its stationary points is interwoven within the method’s paradigm. The performance profile, efficiency, and effectiveness of the method are illustrated on fifteen systems, which comprise as many as 12 components and are at a wide range of temperature and pressure conditions. The new approach is a first step toward the development of a general solver that can be used with a complete guarantee not only in phase equilibria modeling and calculations but in problems requiring the determination of all minima (maxima) of multimodal functions as well.

■

the recent excellent and extensive review of Zhang et al.23 who examine the state-of-the art of global optimization methods. The review offers an analytical, systematic, and in-depth analysis of the strong points and shortcomings, and of a number of other significant aspects of global optimization methods’ application to the solution of a range of important problems in phase equilibria modeling and calculations such as phase stability analysis, Gibbs free energy minimization, parameter estimation, etc. Additional concise enumeration of the merits and limitations of deterministic and stochastic global optimization methods can be found in Bhargava et al.,4 who conclude that finding an optimization approach capable of solving the different problems in phase equilibria modeling is still a very hot research topic in the current literature. The present contribution is a continuation of our research devoted to the phase stability analysis problem3,12,32−34 but outlines a different approach to its solution. It is a fresh first step toward the development of a general solver that can be used with confidence and with a complete guarantee not only in phase equilibria modeling and calculations but in problems requiring the determination (location) of all minima (maxima) of multimodal functions as well. We address here the phase stability problem under constant temperature and pressure (i.e., PT-stability), and we apply equations of state (EoSs) as the thermodynamic model. In the particular examples that we have studied, cubic EoSs are used, but the solver is designed to work with any EoS regardless of its complexity. The application of EoSs to the phase stability problem results in a more difficult formulation than when activity-coefficient equations are used. This is because the objective function in this case is multivariable, nonconvex, and highly nonlinear. Furthermore, the TPDF often has unfavorable

INTRODUCTION The phase equilibrium problem has been around for many years, but it is still an object of continual interest. Among the numerous reasons for this interest is the fact that phase equilibria calculations are executed nowadays not only as standalone calculations but also as a part of a more general simulation, synthesis, or design case, performed by a process simulator, where they are run thousands of times. In these cases, the accent is on the reliability and robustness of phase equilibria calculations since their failure at any stage will disrupt the solution process completely. At the heart of the phase equilibria problem is the phase stability analysis, which, after the publication of the seminal work of Michelsen1 who introduced the tangent plane distance function (TPDF), following Baker et al.2 and Gibbs before them, has dramatically changed the course of research toward the direction of advocating new reliable and efficient numerical technique to solve the problem.3 The robustness of the phase stability test is the key to successful phase behavior and phase composition prediction. Over the years, the phase stability analysis problem has been tackled using different solution methods. These methods can be divided into two big groups: local and global. Local methods include reformulations of Newton’s method. Global methods are either stochastic (based on the application of exploration and exploitation strategies, which have been developed from principles of natural phenomena and artificial intelligence) or deterministic optimization methods and guarantee, at least theoretically, that the optimum is obtained irrespective of the starting point.3,4 The number of papers and reviews devoted to the application of different optimization methods to the solution of the phase stability problem, in the past decade only, is considerable.3−31 The references cited here are given rather as an illustration of the amount of work performed and of the wide variety of optimization methods used by the different research groups than with the ambition to be exhaustive. Interested readers are prompted to acquaint with © XXXX American Chemical Society

Received: April 4, 2013 Revised: June 26, 2013 Accepted: July 17, 2013

A

dx.doi.org/10.1021/ie401072x | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

k* = ln φi(y*) + ln yi* − hi

attributes such as discontinuity and nondifferentiability (particularly when cubic CEoS are used for modeling thermodynamic properties). Additional complexities arise near the phase boundaries, in the vicinity of critical points or saturation conditions, and when the same model is used for determining the thermodynamic properties of the mixture.3,7,23 Consequently, TPDF has several local minima including trivial and nonphysical solutions. The paper is organized as follows: First, the modification, introduced by us,32 of the well-known TPDF1 and the importance of finding all its stationary points (s.p.) are briefly discussed. Next, the paradigm of the new approach, which requires the completion of four stages, is presented. We continue with a detailed description of each of these stages. Then, the robustness, efficiency, and effectiveness of the solver are demonstrated on the example of 12 systems, widely used in the literature as bench-mark tests. We conclude with a brief outline of the future developments of the present method with the view to extend its application to a wide class of problems. Tangent Plane Distance Function and Modified Tangent Plane Distance Function. The tangent plane criterion was formulated and solved theoretically (though not computationally) more than 100 years ago by Gibbs.3 A practical numerical realization of the Gibbs tangent plane stability test has been developed (derivations and figure illustrations) and given by Baker et al.,2 while a first implementation was presented by Michelsen.1 These papers have revolutionized the way phase equilibrium calculations are performed and have spawned a steady series of papers and algorithms refining the basic ideas.3 Minimization of the TPDF has been found to be more reliable and robust than minimizing the Gibbs function directly, since theoretically to guarantee a global minimum for G, using global optimization, the formulation must contain a large number of variables, making a global optimization algorithm computationally expensive, particularly for larger systems. For these reasons, the minimization of the TPDF presents many advantages. Modified TPDF. The new approach to the solution of the phase stability problem seeks to determine all minimum points of a modified TPDF, Φ(y), as advocated originally by Stateva and Tsvetkov:32

This quantity can be either positive or negative. When all calculated k* are positive, the system is stable. The original TPDF1 and the modified one32 are completely equivalent, as demonstrated in refs 32−34. In principle, it is very difficult to give a complete guarantee that all minima of the TPDF (all zeros of Φ(y) in our case) are located because their number N is not known a priori. From this point of view, the most important performance measure of a solver’s reliability is whether or not it is capable of finding all minima of a TPDF (or of any other multimodal function to that matter), which is directly connected with the reliability of the stopping criterion/criteria of the optimization methods applied. Different authors accept different criteria how to consider the search complete. The most obvious, but very demanding computationally because of the unnecessarily large increase of problem dimension, is based on thermodynamic considerations and requires setting as an upper bound the phase rule for nonreacting systems. Wasylkiewicz et al.35 proposed a topological criterion for stationary points relating the number of minima, maxima, and saddles in the case of liquid−liquid equilibrium. Unfortunately, this criterion is not valid for systems involving vapor and liquid phases.20,36 Sofyan et al.11 state that the number of stationary points is always an odd number. It was, however, later demonstrated that the number of stationary points can be even.20 A common argument is that it is not necessary to find all minima; actually, once a minimum of the TPDF to which a negative k* corresponds is located, the process can be abandoned, as the system is unequivocally identified as unstable. Still, stability analysis on its own cannot determine which is the stable phase configuration of a system (identified as unstable at the given temperature and pressure), corresponding with the minimum of the Gibbs energy. Also, use of compositions corresponding to the TPDF global minima to initialize flash calculations does not guarantee convergence to correct phase behavior and compositions. A second stability test should be performed to validate the correctness of flash calculations. Consequently, the focus of phase stability analysis has recently shifted to providing good starting points to improve the reliability and convergence of phase equilibrium calculations. It has been argued and demonstrated that localization of all zeroes of the TPDF is the key element of the operation of multiphase algorithms.3,32−34 This is because the zeros of the TPDF provide very good approximations to the compositions to which the consequent flash calculations that have to be preformed in order to determine the phase configuration corresponding with the minimum Gibbs energy, will converge. Another issue of the ongoing debate in the literature is devoted to the question: which is the most appropriate optimization method to be applied to solve the phase stability problem? Since there is no theoretically based guarantee of reliability (owing to the nonconvex and nonlinear nature of the thermodynamic models used to describe the Gibbs free energy in the TPDF and because of the difficulties that may arise in solving such problems by standard methods), it is widely accepted that global optimization methods are the better alternative for the solution of the problem. Yet, some of the most important aspects of a global solver, among which are (i) how it balances global vs local searches since emphasis on global search increases the number of function evaluations

Nc

Φ(y) =

∑ (ki+ 1(y) − ki(y))

(1)

i=1

where ki(y) = ln φi(y) + ln yi − hi

i = 1, 2, ..., Nc

hi = ln zi + ln φi(z)

(1a) (1b)

and kNc+1(y) is equated to k1(y). From eq 1, it follows directly that when y = y*, where k1(y*) = k 2(y*) = ... = k Nc(y*)

(3)

(2)

then Φ(y*) = 0, which is its minimum value. The zeros y* of the function Φ(y) are its minima, and from a computational point of view, it is more advantageous to search for minima with known (zero) minimum value. The zeros correspond with points on the Gibbs energy (hyper)surface, where the local tangent (hyper)plane at each y* is parallel to that at z. To each y*, a quantity k* corresponds, such that B

dx.doi.org/10.1021/ie401072x | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Stage 1. Determination of Starting Points. The solver uses two sets of starting points. The first set comprises starting points corresponding to “almost pure phases”. As noted by many authors, these starting points are not sufficient particularly when a local solver is used; see, for example, ref 37. The second set comprises a vector y0, which consists of Nc individuals (y0i ) according to

significantly while emphasis on local search decreases the number of function evaluations but decreases the reliability of the resulting solution, (ii) how effectively the repetitive determination of the same minimum is avoided, and most importantly, (iii) how the reliability of the stopping criterion/a is guaranteed, by no means automatically find answers. Thus, it is our belief that a conceptually novel approach, which does not apply global methods but addresses robustly and effectively the above issues, can prove to be a very successful alternative. The new method will serve as a basic element and will mark the avenue toward the design of a general solver framework that will fit in with the requirements of a different perspective to the solution of a wider range of problems. Outline of the New Approach. Paradigm of the Approach. The approach advocated by us requires the completion of four stages: (1) determination of starting points; (2) localization of a zero of the TPDF; (3a) determination of an auxiliary function corresponding with the zero located in stage 2; and (3b) successive generation of a hybrid surrogate function; (4) evaluation of the stopping criterion and a check whether it is satisfied. A schematic representation of the paradigm of the solver is shown in Figure 1. In what follows each of the above stages of the solver will be outlined briefly.

yir = random[0.0, 1.0],

i = 1, Nc

(4)

and normalized according to yi0 =

yir N

∑i =c 1 yir

y 0 = {yi0 },

i = 1, Nc

,

(5)

i = 1, Nc

(6)

where random[0.0,1.0] is a uniformly distributed random value between 0.0 and 1.0 generated automatically by the machine. Thus, the total number of starting points used by the solver is M = 2Nc. Stage 2. Localization of a Zero of the TPDF. In this stage of the solver, any of the local optimization methods known gradient and gradient free alikecan be used. As pointed out previously, the specific form of Φ(y) (its zeros are its minima), and the fact that it is easily differentiated analytically allows the application of a wide spectrum of nonlinear minimization techniques for locating its stationary points.3 In some of our previous studies, for example,32 we have implemented with success the BFGS with a line searcha Quasi-Newton method with a superlinear order of convergence. In the current version of the solver, we use the Pattern Search Method (PSM), a representative of one of the most popular classes of methods to minimize functions without the use of derivatives or of approximations to derivatives. PSM is a local method, based on generating search directions that positively span the search space. It is conceptually simple, can be designed to rigorously identify points satisfying stationarity for local minimization (from arbitrary starting points), and demonstrates good convergence properties. Moreover, PSM can be easily imbedded in a general solver paradigm without jeopardizing the convergence for local stationarity. A detailed description of the PSM can be found in refs 38−40, and hereunder, a very concise outline of the method will be presented. The following problem is under consideration: min(Fn(y))

s.t. Fn(y) ≤ εFun 0 ≤ yi ≤ 1,

i = 1, Nc

(7) (7a)

Nc

∑ yi = 1 i=1

(7b)

where εFun is a small positive real number. Taking into consideration that yNc = 1 − Σi =Nc 1− 1yi, the equality condition, eq 7b, can be replaced by the following inequality: Nc − 1

∑ Figure 1. Schematic representation of the solver’s paradigm.

i=1

C

yi ≤ 1

(7c) dx.doi.org/10.1021/ie401072x | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

In ym = {ym,i}, i = 1, 2, ..., Nc are the coordinates (in mole fractions) of a point in the thermodynamic phase space and An is a constant positive number. The distance R (eqs 8a and 8b) is calculated by Euclidean norm as follows:

Equations 7a and 7c constitute the composition space, which is a convex simplex for the Nc − 1 independent variables. Obviously, the nonconvexity of the problem arises from and is inherent to the thermodynamic model applied. By introducing the inequality constraint (eq 7c), a 2-fold goal is achieved: first, the dimension of the optimization problem is reduced from Nc to (Nc − 1); second, the constrained optimization problem given by eqs 7, 7a, and 7c is now a bit easier to solve than the one given by eqs 7, 7a, and 7b because, generally, inequality constraints are more straightforward to handle by any optimization method (including the PSM for bound constrained minimization we are using). The tolerance-related parameters that directly influence how and when convergence of the PSM is achieved and that also affect the solution are εCon, εFun, and εy. Tolerance on constraints εCon is the termination tolerance placed on constraint violations and represents the maximum value by which parameter estimates can violate a constraint and still allow successful convergence. Tolerance on function εFun is the termination tolerance placed on the log-likelihood objective function. Successful convergence occurs when the log-likelihood function value changes by less than εFun. Tolerance on variable εy is the termination tolerance placed on the estimated parameter values. Similar to εFun, successful convergence occurs when the parameter values change by less than εy. The values of the tolerance-related parameters of the PSM are embedded in the commercial software used by us, namely MATLAB version 7.2.0.232.41 For example, the value of the tolerance on constraints by default is 10−7, while the tolerance on function and the tolerance on variable have the same default value of 10−6.41 Stage 3a. Auxiliary Functions: Generation and Application. A crucial point for any solver is how to guarantee that a minimum already located will not be determined again; in other words how to prevent the solver from being trapped into the current promising area and, hence, from losing the ability to search for other minima in an effective manner. A common approach is to “destroy” the minimum already determined by placing a pole at the minimum point and generating new search directions away from it.20 Our approach is different and is based on the idea of a consecutive transformation of the TPDF (hyper)surface in such a way as to provide a complete assurance that there will be no repetitive determination of the same minima. This is achieved by generating, for each zero located with coordinates yn* = {yn,i}, i = 1, 2, ..., Nc, an auxiliary function (AF) which is nonnegative, smooth, and twice differentiable. An AF is a “substitute” of a zero and modifies the TPDF in such a way as to “eliminate” the s.p. located from further consideration. A function that satisfies the above requirements is of the form A = AFn(A n , y *n , ym)

i = (Nc − 1)

R=

if

if

m≠n (9)

In the general Nc dimensional case, an auxiliary function is represented by an Nc − 1-dimensional convex body of revolution, which is characterized by the two parameters An and Rext n . The parameter An determines the maximal value of an auxiliary function (in geometrical terms it represents the height of the body of revolution) and should satisfy the following condition: A n > εAF (10) where εAF is a small positive real number. Rnext defines an area, denoted AA, in the vicinity (proximity) of a current zero located that should completely encompass that zero, according to R next =

R+i

max (R i+ , R i−)

i = 1,(Nc − 1)

(11)

R−i

where and are the distances from the current s.p. located to a point on the TPDF, for which Φ(y +) = Φ(y −) = εAF

(12a)

and Φ(y) < εAF

(12b)

holds. Figure 2 depicts also the area AA for a two-dimensional case. Details of an auxiliary function determination are shown in Figure 3 for the methane + hydrogen sulphide system with an overall composition in mole fractions z = (0.9813; 0.0187) at T = 190 K and p = 40.53 bar. It should be noted that the auxiliary function determined is the one corresponding with the respective feed (trivial solution, which is a zero of the TPDF known a priori) of the system.

(8a)

R < R next

and

Figure 2. Illustration of the determination of the parameters R and Rnext, characterizing an auxiliary function for a two-dimensional case. The area AA is the area bounded within the red circle.

or ⎛ πR ⎞ A = A n cos⎜ ext ⎟ ⎝ 2R n ⎠

(ym, i − yn, i )2 )

where n is the number of the current zero of the TPDF located and m is an arbitrary point in the thermodynamic phase space. For a two-dimensional case the determination of the distance R is illustrated in Figure 2.

(8)

R > R next

∑ i=1

where either

A=0

(

(8b) D

dx.doi.org/10.1021/ie401072x | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Figure 3. Illustration of the auxiliary function AFz determination for the feed with coordinates z = (0.9813; 0.0187) of the methane + hydrogen sulphide system at T = 190 K and p = 40.53 bar. Blue line, the TPDF shown for the interval 0.974 < y1 < 0.99; red dashed-line, the auxiliary function AFz; Pink dot-dashed line, εFun; green dot-dashed line, εAF. The high values of the two parameters εFun and εAF shown are chosen intentionally for illustrative purposes only.

Stage 3b. Surrogate FunctionGeneration. A pivotal point of the solver’s performance is the consecutive generation of a surrogate function (SF), which, when in its entirety, will serve as a substitute of the initial TPDF. The general form of the surrogate function is

Int =

i

⎧ εAF ⎪ Fint(y) = ⎨ ⎪ ⎩ Fn(y)

j=n

+

∑ j=2

y *j , y),

∀ n,

∏

dyi

i=1

(14)

Fn(y) > εAF ⎫ ⎪ ⎬ Fn(y) < εAF ⎪ ⎭

(15)

where

Fn(y) = Φ(y) + AFz (A z , R zext , z , y) AFj(Aj , R jext ,

∭y =0,1

(Nc − 1)

Fint(y)

∈ 2, N

if if

Thus, there is a mathematical guarantee that all stationary points (zeros) of the TPDF are located

(13)

where Φ(y) is the TPDF (eq 1), AFz(Az, Rext z , z, y) is the auxiliary function for the trivial solution, and AFj(Aj, Rext j , y* j , y) is the jth auxiliary function corresponding with the jth stationary point (zero) determined. The surrogate function (eq 13) is nonnegative for all y and differentiable without points of discontinuity of the second kind. The surrogate function is a hybrid function that is represented geometrically by a combination of several elements: parts of the TPDF (eq 1); the AFz, superimposed on the trivial solution (see Figure 3); and the auxiliary functions superimposed on the (N − 1) stationary points located. Thus, the total number of the auxiliary functions is N and is equal to the number of zeros of the TPDF. Needless to mention that since N is not known a priori, the number of the AFj elements of the SF is not known in advance as well. At this point a brief comment highlighting further the organizational details of the solver’s paradigm should be made: Since the trivial solution is the only zero of the TPDF known in advance, the actual optimization process begins af ter an auxiliary function for the trivial solution (the feed), namely AFz, is determined. The surrogate hybrid function (eq 13) is thus generated successively during the iteration process by adding elements to the third term on the RHS of eq 13. Figure 4 depicts the consecutive generation of the hybrid surrogate function for a mixture of methane + hydrogen sulphide, with an overall composition in mole fractions (0.9813; 0.0187) at T = 190 K and p = 40.53 bar. Stage 4. Stopping Criterion. The stopping criterion is evaluated from

If and only if If

Int = 1.0.

(14a) (14b)

Int < 1.0

that means that there exists at least one zero in the range {yMIN , i yMAX }, i = 1, Nc , which has not been yet located by the solver. i In the above, {yimin = 0.0}, ∀ i = 1, Nc

∀ i = 1, Nc ,

{yimax = 1.0}, (16)

Figure 4b illustrates the case when the stopping criterion of the solver is satisfied, while Figure 5 illustrates the case when the stopping criterion is still not satisfied for the system CO2 + methane, with an overall composition (0.6; 0.4), at T = 220 K and p = 60.8 bar. The stopping criterion advocated (eqs 14) is robust, completely reliable, and ensures sufficient exploration search power of the solver. The iteration process should be terminated if and only if Int = 1.0. (Figure 4b). At this point, some further elaborations are called for. The inherent connection between the surrogate function (eq 13) and the stopping criterion (eq 14) is realized via εAF, which, in principle, is the parameter responsible for the correct generation of the surrogate function. Thus, if the surrogate function is generated correctly on each step and in its entirety, and if the stopping criterion is satisfied then there is a complete mathematical guarantee that either all s.p. (zeros) of the TPDF are located by the solver or that none exist (the system under consideration is globally stable). E

dx.doi.org/10.1021/ie401072x | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

Figure 5. Hybrid surrogate function for the system carbon dioxide + methane with an overall composition (0.6; 0.4) at T = 220 K and p = 60.8 bar. Illustration of the case when the stopping criterion is not yet satisfied. Black line, the SF; green dot-dashed line, εAF. The high value of the parameter εAF shown is chosen intentionally for illustrative purposes only.

Still, the following should be made clear: there are numerous contributions published in the specialized literature that advocate robust and efficient algorithms and methods designed to reduce the cost of the numerical evaluation of multidimensional integrals, because that problem is of particular interest to physics, quantum chemistry, etc. Hence, if it is important that the computational time should be reduced then any of those algorithms can be adopted and interwoven in the software’s environment of the solver. That discussion, however, is very much outside the scope of the present contribution and hence will not be elaborated further.

■

Figure 4. TPDF and the auxiliary functions generated for the feed and the three s.p. located for a mixture of methane + hydrogen sulphide (a) and hybrid surrogate function in its entirety for a mixture of methane + hydrogen sulphide, with an overall composition in mole fractions (0.9813; 0.0187) at T = 190 K and p = 40.53 bar. Illustration of the case when the stopping criterion is satisfied. Color key: Black line, the SF; blue line, the TPDF; red dashed line, the auxiliary functions generated for the feed and for each of the three zeros of the TPDF located; green dot-dashed line, εAF. The high value of the parameter εAF shown is chosen intentionally for illustrative purposes only.

RESULTS AND DISCUSSION The performance of the new solver was tested using a great number of examples, quite different in naturefrom highly difficult theoretical problems to large practical problems in petroleum and chemical engineering. Some of the systems examined were chosen primarily for the high degree of complexity, which provided a rigorous test for the solver’s capabilities. From all the systems examined by us we have selected to illustrate the performance profile, efficiency, and effectiveness of the solver on just 12 examples, all of which are very well-known bench-mark tests from the literature. The systems are at a wide range of temperature and pressure conditions and comprise as many as 12 components. In order to ensure a complete and reliable basis for comparison of our results with those of the other authors, we show the thermophysical properties of the pure components used in Table S1 (Supporting Information). As mentioned previously, the PT-stability problem is addressed here applying cubic EoSs as the thermodynamic model, namely either the SRK or the PR EoS. Furthermore, the solution process is organized in such a way that when the EoS has three real roots only the one corresponding with the lowest Gibbs energy is accepted, the other two are rejected. Lastly, the EoSs binary interaction parameters (BIPs) employed (Table S2a,b, Supporting Information) were those of the original references.

If, however, any element of the SF is generated incorrectly but the stopping criterion is satisfied then it cannot be guaranteed that the solver has located all zeros of the TPDF. It must be underlined that the stopping criterion (the integral, eq 14) is evaluated numerically. For systems with low to medium dimension, this is not computationally demanding and the numerical integration can be realized with a high degree of accuracy, at a very modest computational cost. In order to guarantee high accuracy of the approximation the adaptive Simpson quadrature method of numerical integration, as advocated in MATLAB,41 is interwoven in the solver’s environment. For systems with high dimension, however, a large number of calculations will be required, which, not unexpectedly, will increase computational effort. F

dx.doi.org/10.1021/ie401072x | Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

Article

If a too high value of εAF is chosen then, for a zero located, the corresponding auxiliary function (area in the twodimensional case/body of revolution in the Nc dimensional case) could overlay and thus “eliminate” another zero that might be located in the proximity of the former. As a result, the particular “wrong” element of the SF generated at this step will distort the representation of the TPDF (hyper)surface in its entirety, regardless of the fact that all other elements of the SF were obtained correctly. Furthermore, in the worst case scenario, it might so happen that the “eliminated” zero was the only one to which a negative k* corresponds and consequently the system will be identified as a “stable” rather than as “unstable” one. This is illustrated in Figures 6 and 7. Figure 6 shows the TPDF (in the interval 0.61