A New Homotopy-Based Strategy for the Robust Determination of All the Feasible Solutions for CSTR Systems

A new homotopy-based strategy is presented that can be used in the robust determination of multiple steady-state solutions for continuous stirred tank reactor (CSTR) systems. The strategy relies on the features of homotopy parameter and variables bounding, and requires only that one feasible solution of the system is either known beforehand or can be solved with an existing solving algorithm. The strategy systematically results in all the multiple solutions, or alternatively confirms that the problem does not have multiple solutions, within the predefined problem domain. The strategy was successfully demonstrated with CSTR cases gathered from the literature. Finding all the feasible solutions was verified in simple CSTR systems by applying tools available in the literature. Variables bounding constrained the homotopy path to travel only within the pre-defined variable domain. The strategy is applicable for determining multiple steady states for a variety of chemical engineering systems.


Introduction
The steady-state and dynamic behaviour analyses of continuous stirred tank reactors (CSTRs) have received considerable interest in the field of chemical engineering.Even though the search for complete analysis methods of CSTR behaviour has been fairly exhaustive, especially during the 1980s, the quotation from Macbeth Act V in the title of the study by Farr and Aris [1]: "Yet who would have thought the old man to have had so much blood in him?"describes well how intriguing and challenging CSTRs have proven to be.The field continues to receive some attention through the introduction of more complex CSTR systems (see e.g.Mjalli et al. [2]; Švandová et al. [3]; Kiss et al. [4]; Brooks [5]; Yermakova and Anikeev [6]; Waschler et al. [7]).However, the emphasis in the current studies in the literature has shifted towards the investigation of CSTR systems from the perspective of numerical solving methods [8][9][10][11][12].The reason for this is that CSTR models exhibit characteristics that pose challenges to practically all the available general-purpose numerical solving methods.Hence, they are naturally interesting case study examples for different novel solving methods.This is also the focus of this study.
CSTRs are typically modelled by relying on the material balances of the components, the enthalpy balance of the reactor and the enthalpy balance of the cooling, or heating, medium.In the case of a two-or multi-phase reactor system, the CSTR model also includes phase equilibrium relations for the component system under investigation.
When comparing the CSTR model to models of other process units in chemical engineering, such as distillation or absorption column models, it can be seen that the total number of equations in a CSTR model is relatively small.This, however, does not signify that the operation and behaviour of a CSTR is straightforward and simple.In fact, a CSTR is the simplest type of chemical reactor that exhibits both multiple steady states and oscillations [13].The complexity of CSTR behaviour is due to the combined effect of chemical reactions and interaction between the chemical reactions and phase equilibrium.Therefore, CSTR models are generally highly non-linear.The nonlinearity is reflected in the form of multiple steady states and a 1 Chemical Process Engineering, Faculty of Technology, University of Oulu, PO Box 4300, FI-90014 University of Oulu, Finland multiform dynamic behaviour including instability, sustained oscillations, and phenomena such as strange attractors and chaotic behaviour (see e.g.Jayakumar et al. [14]; Kiss et al. [4]; Pérez-Polo and Pérez-Molina [15]; Razón and Schmitz [16]; Russo and Bequette [17]).
The non-linearity and subsequent occurrence of multiple steady states of a CSTR model may arise even when the modelled reactor system seems to be extremely limited in terms of sources of non-linear behaviour, i.e. even in the absence of thermal effects (see Horn and Jackson [18]).Uppal et al. [19][20] analysed the steady-state operation for the system of an irreversible first-order, exothermal, Arrhenius-type reaction in a non-adiabatic CSTR and observed that as many as three steady states may exist.Similarly, Uppal et al. [19][20] showed that six different multiplicity patterns and as well as a unique steady state exist for the system when residence time is used as the bifurcation parameter, i.e. the system has in total seven bifurcation diagrams.This analysis was later confirmed by Golubitsky and Keyfitz [21] using singularity theory.The number of possible multiplicities and especially the number of possible bifurcation diagrams increases rapidly with an increase in the complexity of the system.For example, in the case of consecutive first-order Arrhenius-type reactions in a non-adiabatic CSTR, seven steady states may exist and 23 different bifurcation diagrams can be formulated [16].
Even though a model describing a CSTR may have multiple solutions, only some of the solutions correspond to conditions where the state of CSTR is stable and thus feasible for practical operation [23].The stability of each solution with respect to practical operation can be ascertained, for example by dynamic simulation.The dynamic simulation of the system is definitely necessary for system control, management and safety analysis purposes.However, without extensive knowledge about the steady-state behaviour of the system, the dynamic simulation might yield insufficient operation analysis results.Hence, it is of great value if we can define at least the maximum number of steady-state solutions of a CSTR model and furthermore determine parameter regions corresponding to a specific number of solutions.As a consequence, the prediction of multiplicity patterns for different CSTR systems has received considerable attention, especially in the studies of Balakotaiah and Luss [24][25][26][27][28][29].
Even though the CSTR systems were investigated in depth during the 1980s, attention has been paid mainly to the prediction of multiplicity patterns and stability analysis of CSTR dynamics, and the formulation of robust solving strategies to obtain the corresponding steady-state profiles has not been as extensive.In addition, it is worth noting that the majority of the studies concentrate on the investigation of simplified CSTR systems, e.g. the effect of the cooling jacket temperature dynamics is not typically incorporated in the model [17].Therefore, there is a clear need for a robust method for the determination of all the feasible steady states of more complex CSTR systems, and evidently for more general chemical engineering systems with given specifications.In addition, at the moment, there is no robust routine in commercial steady-state simulation packages to perform this analysis.
In this paper, a solving strategy is proposed that aims at tackling the weakness of the commercial steady-state simulation packages mentioned above.The strategy exploits the fact that generally one feasible solution for a process model can be found with some existing solving algorithm at a reasonable calculation cost.The solution found is then used as a starting point for the Newton homotopy-based solving method, which in turn is used to find all the remaining feasible solutions of the model.
In this study, it will be shown that the target of finding all the feasible solutions for the investigated CSTR systems can be reached with the Newton homotopy-based solving method if it is equipped with the features of homotopy parameter and variables bounding presented in Malinen and Tanskanen [30][31].The solving strategy of the present work extends the applicability of the method presented in Malinen et al. [32] in the determination of multiple steady states of chemical engineering systems.In addition, the method proposed by Malinen et al. [32] has a shortcoming when applied to systems where the trivial solution is located in the vicinity of the variable domain.This shortcoming is tackled in this study.On the other hand, the focus in this work is to demonstrate the characteristics and robustness of the proposed solving strategy in various CSTR cases gathered from the literature.The cases investigated are ordered in terms of increasing complexity.In the simpler examples, the number of solutions can be verified through the analysis methods available in the literature for CSTR systems.For the two latest CSTR systems this kind of conclusive method is not available and hence applying a robust solving strategy, like the one presented in this work, to find the solutions is the only feasible way to analyse the multiplicity behaviour of the system thoroughly.

Characteristics of some present solving methods
When the CSTR is considered to operate at a steady state, the mathematical description is reduced to a set of equations: Basically, the solving of this non-linear equation set can be carried out with locally convergent, Newton-Raphson-based solving methods.Locally convergent solving methods usually work satisfactorily if a feasible initial guess is supplied.In addition, when an initial guess close to the solution is available, the solving methods converge superlinearly to a solution of the CSTR model.However, the locally convergent solving methods are able to converge only to a single solution from an initial (1) guess.Since a CSTR model may have several feasible solutions, achieving at least some of them robustly would require a large set of initial guesses over the entire feasible domain.Thus, it is evident that this kind of 'trial and error' method is incapable of resulting in all the feasible solutions with full certainty.
In contrast, homotopy continuation methods form an attractive set of solving methods, which offer a means of reaching multiple steady-state solutions for a CSTR model from one starting point.Homotopy continuation methods can be roughly divided into problem-dependent and problem-independent homotopies.The division can be made on the basis of the role of the continuation parameter.
In problem-dependent homotopies, the continuation parameter has a physical meaning.The homotopies can be used to constitute bifurcation diagrams where for example the conversion of a reactant in a CSTR is presented as a function of residence time.Problem-dependent homotopies have been used in several CSTR studies (see e.g.Vadapalli and Seader [12]; Wang et al. [33]; Yermakova and Anikeev [6]).Since problemdependent homotopies are inherently strongly problem-tailored, they are not directly applicable as such to the solving of various chemical engineering problems.
Problem-independent homotopies, on the other hand, can be applied directly to the solving of various chemical engineering problems, because the continuation parameter is an artificial parameter without physical meaning.Problem-independent homotopies are based on the homotopy function: The basic principle in problem-independent homotopies is to track a homotopy path defined by h(x, θ) from a starting point, (x 0 ,0), to a solution, (x * ,1).Depending on the selection of the auxiliary function, g(x), different homotopies are obtained.As an example, by defining g(x) = f(x) − f(x 0 ), the traditional Newton homotopy is obtained.The homotopy continuation methods have been discussed in detail, e.g. in Allgower and Georg [34], Gritton et al. [35], Kovach and Seider [36], Kuno and Seader [37], Lin et al. [38], Seader et al. [39], Seydel [40], Seydel and Hlavacek [41], Wayburn and Seader [42], and Watson [43].Hence, the theory behind homotopy continuation methods is discussed only when appropriate in the present work.
One advantage of the homotopy continuation methods is that several solutions for Eq.(1) can be reached on one homotopy branch.However, the methods cannot guarantee that all the feasible steady-state solutions of a chemical engineering problem are approached with certainty on one continuous homotopy branch.The main reason for this is that the homotopy path may travel outside the feasible variable domain.It is also due to the absence of real space connections between separate homotopy path branches, thus preventing the attainment of solutions on isolated branches.Recently, Rahimian et al. [8][9] combined the features of fixed-point (FP) and Newton (N) homotopies to form a new homotopy called as the FPN homotopy.Rahimian et al. [8][9] applied FPN homotopy in the solving of a multitude of different chemical engineering related problems to reach all the feasible solutions of the problems successfully.However, the formulation of FPN homotopy does not include variables bounding, and hence the path may travel outside the feasible variable domain.
The challenges mentioned above can also be encountered in the solving of a CSTR model with traditional homotopies.Herein, the adiabatic CSTR with a first-order exothermal reaction presented in Chapter 4.2 is used as an illustrative example.As shown in Fig. 1a, in some cases the starting point is such that the traditional Newton homotopy successfully results in a homotopy path which passes through all the solutions on the θ =1 plane.However, as Figs.1b-1d show, separate homotopy path branches are formed when the starting point is changed.In these cases, some of the solutions of the original equation set exist on a different homotopy path branch than the starting point.Thus, the wrong conclusion about the number and existence of multiple solutions might be drawn.
It is noteworthy that when a starting point is also a solution of the original equation set (see Fig. 1d), the traditional Newton homotopy results in as many separate homotopy path branches as there are solutions in the original equation set.As Fig. 1d illustrates, the branches formed exist parallel with each other and travel from negative infinity to positive infinity with respect to the homotopy parameter.
As Fig. 1e indicates, in the worst case, none of the solutions are on the same homotopy path branch as the starting point.Therefore, the solutions on the isola branch would be completely missed, and the wrong conclusion might be made, namely that the problem has no feasible solutions with the given specifications.Thus, as Fig. 1 illustrates, the traditional Newton homotopy is not able to find all the steady-state solutions of a CSTR system robustly.
Fig. 1 also illustrates another fundamental problem with the traditional problem-independent homotopy methods.Similarly to locally convergent solving methods, they do not prevent the solution path from travelling outside the feasible variable domain.As a consequence, negative mole fraction and temperature values for example might be encountered along the solution path (see Fig. 1c).Negative values, when substituted into logarithm or square root functions, result in complex numbers, which are not usually tolerated in the thermodynamic subroutines of process simulation packages.The solving may also be completely interrupted by a fatal error if the variable value of pure zero is substituted into a logarithm function.
The above-mentioned weaknesses, i.e. the impact of the starting point on the accessibility of the solutions and the challenge of keeping the homotopy path within the feasible problem domain, were also recently demonstrated by Giunta et al. [44] in their study of fixed-point homotopy to find multiple steady states for (2) the system of CO oxidation in a catalyst pellet.Evidently, there is a general need for a problem-independent homotopy method, which would allow the achievement of all the feasible solutions robustly and would also keep the variable values strictly within the feasible problem domain boundaries throughout the solving.
Efforts aiming at developing such a method include the bounded homotopies initially proposed by Paloschi [45][46][47], and later revised as modified bounded homotopies by Malinen and Tanskanen [29,48].Bounded homotopies prevent the homotopy path from exceeding the feasible problem domain and thus prevent failures originating from the substitution of unfeasible variable values in thermodynamic subroutines.
In order to improve the capabilities of homotopy methods to approach solutions existing on separate homotopy path branches, Malinen and Tanskanen [31] proposed the concept of homotopy parameter bounding.However, as was concluded in Malinen [49], when used separately, both the concept of homotopy parameter bounding and bounding with respect to problem variables are incapable of finding all the solutions if the solutions exist on isolated branches.
Fig. 1 The homotopy paths formed with the traditional Newton homotopy method for the set of Eqs. ( 11) and ( 12) describing the adiabatic CSTR with a first-order exothermal reaction when the starting point is a) ( , ) ., c T , .c T ) and e) ( , )   ., .c T 0 0 0 2 1 3
Recently, Malinen et al. [32] presented a Newton homotopybased method, which was equipped with the features of both homotopy parameter and variables bounding.It was noticed that the method was capable of robustly determining all the stationary points for the tangent plane distance function (TPDF).Since the method presented in Malinen et al. [32] utilises the properties of the TPDF problem to determine all the solutions of the problem, it is not as such directly usable in the solution of a general chemical engineering problem.In this study, the method presented in Malinen et al. [32] is modified to extend its applicability to CSTR systems.

The applied solving strategy
The core of the solving strategy presented here is a modified Newton homotopy method having features of both homotopy parameter and variables bounding.The method applied in this paper is a revised formulation of the method presented in Malinen et al. [32].

The modified Newton homotopy method with both homotopy parameter and variables bounding
The modified Newton homotopy formulation applied in this paper can be written as ) is a penalty function resulting in a scalar.The penalty function annihilates the magnitude of the function, f(x inf ), whenever the homotopy path runs outside the predefined problem domain.Correspondingly, the auxiliary functions, M(θ − θ b )e and f (x x x ) ' )( inf i nf ,inf 0 − b , compensate for the annihila- tion, thus resulting in a path that is bounded with respect to the homotopy parameter, θ, and mapped variables, x inf .As was noticed in Malinen et al. [32], since there are n + 1 variables, i.e. x inf and θ , but only n equations, the homotopy parameter and variables bounding features in Eq. ( 3) are concurrently unable to keep the path completely bounded.Basically, the path may exceed the problem domain with respect to one of the problem variables or the homotopy parameter.However, by using a small numerical value for parameter M, such as 0.001, the effect of the auxiliary function, f (x x x ) ' )( inf i nf ,inf 0 − b , will be stronger than the effect of the other auxiliary function, M(θ − θ b )e, on the course of the homotopy path.Thus, the boundaries of the feasible problem variable domain are neither reached nor crossed along the homotopy path.The significance of the value of parameter M is investigated later on in this study.
Due to the fact that the homotopy parameter has no physical meaning, the boundaries of the homotopy parameter can in essence be defined arbitrarily.However, in practice, the curvature of the homotopy path should be reasonable from the perspective of a path tracking method and overall solving efficiency.The selection of boundaries with respect to solution efficiency is not, however, the aim of the present work.In the present study, the boundaries of the homotopy parameter have been set as −1 and 1.
It has to be emphasized that mapped problem variables must be used when applying Eq. ( 3).On the basis of the maximum b i max and minimum b i min values that are set for every variable x i , mapping from a finite space into an infinite space is carried out as follows: The details of the symbols and functions in Eq. ( 3) can be found in Malinen et al. [32].It is worth emphasizing that Eq.
(3) and the formulation presented in Malinen at al. [32] differ in the way the Jacobian matrix term is determined in the auxiliary . In Malinen et al. [32], the Jacobian matrix term was determined at the trivial solution of the In general, however, chemical engineering problems have no trivial solutions.In addition, if the Jacobian matrix term is determined at a starting point close to the problem domain boundary, it may cause the numerical values of some non-zero elements of the Jacobian matrix term to become nearly zero.This in turn will pose numerical challenges for homotopy path tracking inside the bounding zone.
In addition, if some non-zero element of the Jacobian matrix is evaluated incorrectly as zero, the homotopy path will not be bounded with respect to the variables.As a result, the homotopy path will exit the feasible domain with respect to the variables rather than the homotopy parameter.
In order to tackle these challenges, in this paper, a generic way to determine the Jacobian matrix term in Eq. ( 3) is imple- is determined at the centre of the mapped variable space, i.e. x 0 With this selection, the problem of evaluating the Jacobian matrix term elements close to the problem domain boundaries is avoided, as it is common that the values of some non-zero elements approach machine precision, and pure zero near the domain boundaries.The values typically become sufficiently significant at the centre of the mapped variable space.Thus, successful path tracking can also be guaranteed, assuming non-singular f'(0) inside a narrow bounding zone with respect to the problem variables.
To distinguish clearly the starting point and the point where the Jacobian matrix f (x ' ) inf 0 is evaluated, it is stressed that herein x 0 inf x is not the actual starting point for homotopy path tracking with Eq. ( 3).Instead, one solution of the original problem (Eq.( 1)) in the mapped variable space, x inf* , acts as the starting point for path tracking.Thus, f(x inf* , 0) = 0, x inf* = x b,inf and θ = θ b at the starting point.Therefore, in order to apply the proposed strategy, one solution of the original problem must be either known beforehand or must be solved with some existing solving algorithm.When any feasible solution of the original problem is available, the solution can be used as a starting point for homotopy path tracking based on Eq. ( 3).It is worth noting that in the present solving method, Eq. ( 3), the solutions of the original problem, Eq. ( 1), are situated on the θ = 0 plane and not on the θ =1 plane as is the situation with the traditional problemindependent homotopy Eq. ( 2).

The solving strategy for multiplicity determination
Fig. 2 represents a two-phase solving strategy for the determination of multiple solutions of an equation set, f(x) = 0, within a feasible problem domain.In Phase I, one solution of the problem is solved, if it is not already available.If the solution has to be solved, the solving can be performed with any existing solving method.Usually, local solving methods, e.g.Newton-Raphson-based solving methods, result in a solution with reasonable calculation expenses.If a local solving method is not able to result robustly in a solution, some (global) solving method, e.g. a homotopy method, can be used instead.When investigating CSTR systems, the problem-dependent homotopies are applicable in Phase I.This is due to the fact that the starting point for the homotopy path is easily formulated by applying e.g. the residence time in CSTR as the homotopy parameter.As a consequence, the CSTR system has only one easily obtainable steadystate solution at the starting point when the residence time has a value of zero.At that point, the reactants have zero conversion and the outlet stream temperature equals the feed flow temperature.If no solution has been found, not even with global solving methods, it is highly likely that the problem has no feasible solution with the specifications.The existence of at least one feasible solution for the problem can be ascertained before solving in certain problems by applying physically realizable values for the system parameters, as in the case of a CSTR in which a single exothermic, first-order reaction occurs.Furthermore, in the system even the multiplicity pattern, i.e. the number of feasible solutions, can be predicted [24]., is determined.Finally, the solution determined in Phase I of the solving strategy is mapped into the infinite variable space, and then the homotopy path defined by Eq. ( 3) is tracked forward and backward from the θ = 0 plane, i.e. from the starting point (x inf* ,0).All the points where the homotopy path intercepts the θ = 0 plane are collected.However, only the points that do not lie inside the bounding zone represent the feasible solutions of the original problem f(x inf ) = 0.The path is tracked until the path intercepts the θ = −1 and θ = 1 planes.Approaching the θ = −1 and θ = 1 planes indicates that all possible solutions within the feasible problem domain have been approached.As will be illustrated with the test cases, the mapped problem variables have the same values with opposite signs at the interception points on the θ = −1 and θ = 1 planes.This behaviour was also noticed for the TPDF problems in Malinen et al. [32].
The multiplicity determination strategy proposed in Fig. 2 can be flexibly implemented as a complementary part of an existing NAE solving algorithm.The strategy requires only that one feasible solution of the problem is either known beforehand or can be solved with an existing solving algorithm.Then, the actual multiplicity determination strategy systematically results in all the solutions, or alternatively the information that the problem does not have multiple solutions, within the predefined problem domain.As illustrated in Chapter 4, the strategy is highly suitable for solving different CSTR system models, where problem variables must have physically meaningful values throughout the solving procedure.

Test problems
To illustrate the performance of the solving strategy introduced in Chapter 3, the solving of four CSTR problems gathered from the literature are considered here.In all the cases, the width of the bounding zone with respect to the unmapped problem variables is kept extremely narrow.The lower l i inf and upper u i inf inner boundaries have been selected as -10 and 10, respectively.As a result, the bounding zone width is 5 • 10 -11 of the width of the problem domain in the unmapped variable space.
The CSTR problems and the proposed solving strategy have been implemented in the MATLAB environment.The actual homotopy path tracking has been carried out with CL_ MATCONT, which is a continuation toolbox in MATLAB [50].

BioCSTR
Following the representation of Problem 8.5 in Finlayson [51], the growing of cells in a CSTR can be described as where Q is the volumetric flow rate, S is the substrate concentration in the reactor, S in is the substrate concentration in the feed flow, V R is the reactor volume, μ m is the maximum specific growth reaction rate, K S is the Monod constant and K I The number of real solutions of this cubic polynomial in Eq. ( 8) can be evaluated with the help of its discriminant ∆: a a a a a a a a a a a a , where If the discriminant • is negative, the polynomial has only one real-valued solution; the other solutions being complex, • equals zero, the polynomial has a multiple solution and all the solutions are real-valued, • is positive, all the solutions of the polynomial are distinct and exist in the real-value space.
Based on the evaluation of the discriminant with respect to changes in the Damköhler number while the other parameters remain unchanged, the BioCSTR system has seven distinct domains: The differences between the domains with respect to the real-valued solutions can be seen in Fig. 3.With the specifications Da = 1.19, ω = 0.00356 and ε = 2.53, the discriminant is positive, as can be observed in Eq. (10).Thus, the problem has three distinct real solutions within the parameter set.In addition, as can be seen in Fig. 3, all of the solutions have positive real values.The substrate concentration inside the reactor cannot exceed the inlet flow substrate concentration, i.e. σ ≤ 1, and naturally the concentration cannot have negative values, i.e. σ ≥ 0. Thus, all of them are feasible since the solutions are bounded in the range of 0 ≤ σ ≤ 1.Thus, for the investigation of the solving strategy characteristics, we define the problem domain boundary values such that b σ max = 1 and 0 min = σ b .In Phase I of the proposed solving strategy, the first solution of the BioCSTR problem is solved.The MATLAB fsolve routine was used as the solver.By using σ inf = 0 as the initial guess for the solving, a solution of the BioCSTR model was found at σ inf* = −0.0584,which corresponds in the unmapped variable space to the solution σ * = 0.4371.In Phase II of the proposed solving strategy, the Jacobian at the centre of the variable space is first calculated.In the mapped variable space . By using the solution (10) σ inf* = −0.0584as the starting point for homotopy path tracking, the homotopy path is tracked both in forward and backward directions and solutions are collected.The path is represented in Fig. 4. As can be seen, all the feasible solutions of Eq. ( 7) are achieved along a single homotopy path.The solutions on the θ = 0 plane with the mapped variables are σ inf* = [−0.0584− 0.5359 −1.3542], which correspond in the unmapped variable space to σ * = [0.43710.1456 0.0221].It can also be noticed that the path intercepts the planes θ = −1 and θ = 1 at the points where the mapped problem variable σ inf has the same numerical value, but with opposite signs, i.e. (θ, σ inf ) = (−1 10.6429) and (θ, σ inf ) = (−1 10.6429).The interception takes place inside the bounding zone where both the homotopy parameter and variables bounding features are in force.The interception location is dependent on the value of parameter M. The effect of the parameter on the course of the homotopy path and the interception location can be seen in Fig. 5.
As can be seen in Fig. 5a, all the homotopy paths, irrespective of the M parameter value, intercept the domain boundaries with respect to the homotopy parameter.However, the location of the crossing moves closer to the variable outer boundary as the M parameter value increases.Thus, from this perspective, the value of the M parameter should be small, as stated previously.Instead, based on Fig. 5b, we can state that the course of the homotopy path is shorter if we apply a higher value for the M parameter.On the other hand, the path may become more difficult to follow if we apply a high M parameter value.Nevertheless, regardless of the parameter value all the feasible solutions of the problem are found.The primary objective in this study is to find all the solutions of the CSTR systems robustly and not to seek for an optimal value for M. Thus, in all the cases in this study a small value of 0.001 for M has been applied.

Exothermal reaction in an adiabatic CSTR
Following the representation in Finlayson [51], steady-state material and energy balances for an adiabatic CSTR in which a first-order reaction (e.g.A → B) occurs can be represented as where c A is the normalized concentration of A in the reactor outlet c A with respect to the feed concentration of A c A,0 , T is the normalized temperature of the reactor outlet T with respect to the reactor inlet temperature T 0 , and the problem parameters γ and β can be represented as where E a is the activation energy of the reaction, R is the gas constant, ∆H r is the heat of reaction, k is the kinetic constant of the reaction, ρ is the reaction mixture density and C p is the heat capacity of the mixture.In the investigated case, k(T 0 ) = 1.Let us also define the Damköhler number Da in this example as 0 , and parameter B, whose value can be applied as an indicator of the appearance of multiplicities in the system as proven by Balakotaiah and Luss [24], as . ρ γβ Applying these parameters, Eqs. ( 11) and ( 12) can be combined to give γ from which c A can be solved and substituted either to Eq. ( 11) or ( 12) to give us also the value of T .On the other hand, according to the analysis of Balakotaiah and Luss [24], the exact uniqueness-multiplicity boundary can be presented as where c A * is the dimensionless concentration of A at the uniqueness-multiplicity boundary.As can be observed from Eq. ( 18), a multiplicity within the feasible domain exists for some Da if, and only if, Hence, in the investigated case the criterion can be presented in the form where γ > 0 by applying the model parameters.However, it is worth noting that the criterion defined in Eq. ( 21) is sufficient for a system to have multiple solutions, if Da is within the domain defined by the uniqueness-multiplicity boundary in Eqs.(18) and (19).Before implementing the solving method, we also need to define the variable boundaries.Since Eqs. ( 11) and ( 12

Single solution
The adiabatic CSTR model has only one feasible solution unless the condition of Eq. ( 21) is fulfilled.Let us apply the following set of parameters: Q/V R =8.7, β = 0.15 and γ = 30.Hence, the right-hand side of Eq. ( 21) is given the value of 0.15385, which is bigger than the value of β.Thus, the model has only one feasible solution.Let us investigate the properties of the proposed modified Newton homotopy-based strategy with this simple problem.Use of the parameters results in the Jacobian matrix .By using this solution as a starting point for homotopy path tracking, Eq. ( 3) results in the path shown in Fig. 6.
As shown in Fig. 6, the path intercepts the θ = 0 plane at the starting point only.This means that with this set of parameters, Eqs. ( 11) and ( 12) have only one solution, as stated above.The path intercepts the θ = −1 and θ = 1 planes at the points where the variables have the same values but with opposite signs.where the path intercepts the θ = −1 and θ = 1planes.

Multiple steady-state solutions
Now, let us investigate the same adiabatic CSTR model, but with different parameters, to evaluate the capabilities of the proposed solving strategy when the model has solutions.We will use the parameter set Q/V R = 25, γ = 30 and β = 0.25 .Hence, β has a value bigger than 0.15385.Thus, the multiplicity condition given by Eq. ( 21) is fulfilled.As can be seen in Fig. 7, with these parameter values, Da = V R /Q = 0.04, the CSTR system has three feasible solutions.
As illustrated in Fig. 7, the system exhibits an S-shape solution curve, which was to be expected based on the analysis of Balakotaiah and Luss [24].In addition, the multiplicity of the solutions occurs within the range of Da Î [0.02630, 0.06032].These uniqueness-multiplicity boundaries can also be calculated from Eqs. ( 18) and (19).It is worth noting that the system has two feasible solutions only where the solution curve crosses the uniqueness-multiplicity boundaries.The feasible solutions with the applied parameter set are represented in Table 1.
The Jacobian matrix evaluated at the centre of the variable space with mapped variables is f' f 0 ( , ) .By using this solution as a starting point for homotopy path tracking, usage of Eq. ( 3) results in the path represented in Fig. 8.
As can be seen in Fig. 8, all the solutions listed in Table 1 are achieved on a continuous homotopy path.It can also be noticed that the path intercepts both the θ = −1 and θ = 1 planes at the points where the variables have the same values but with opposite signs.

Exothermal reaction in a cooled CSTR
Following the representation in Shacham et al. [52], steadystate material and energy balances for a cooled, perfectly mixed CSTR (see Fig. 9) with an irreversible exothermal first-order reaction with respect to component A can be written as In Eqs. ( 22)-( 24), U is the overall heat transfer coefficient, A is the overall heat transfer area, F j is the volumetric flow rate cooling media, ρ j is the density of the cooling media, C j is the heat capacity of the cooling media, λ is the heat of reaction, and T j0 and T j are the inlet and outlet temperatures of the cooling media.The reaction rate coefficient, k, follows the representation where α is the reaction frequency factor.In the process model, Eqs. ( 22)-( 25), negligible heat losses, constant densities and perfect mixing both inside the reactor (reacting mixture) and cooling jacket (cooling media) have been assumed.
Table 1 The feasible solutions of the adiabatic CSTR model with an exothermal reaction.Q/V R = 25, β = 0.25 and γ = 30 T [-]   0   Cooled CSTR systems with simple exothermal kinetic reaction, i.e.Eqs. ( 22)-( 23), have been investigated in a multitude of studies (see e.g.Golubitsky and Keyfitz [21]; Balakotaiah and Luss [24,27]; Russo and Bequette [17].However, the effects of the cooling medium energy balance, i.e.Eq. ( 24), on the multiplicities have received considerably less attention.Russo and Bequette [17] observed that the addition of the cooling jacket dynamics along with the energy balance changes the multiplicity patterns of the system.
Utilising the parameter values presented in Table 2, the equation set, Eqs. ( 22)- (24), has three unknowns: T , c A and T j .Before the solution of the equation set, the feasible variable domain needs to be defined.First, the values of the unknowns must be positive.On the other hand, the unknowns do not have any generally defined maximum for their values.Thus, the upper domain boundary values must be artificially specified.In addition, in order to make the problem behave better numerically, the lower domain boundary for temperatures is specified to have a small positive value instead of 0. In this case, the problem domain boundary values b i max and b i min have been specified as presented in Table 3. Again, the MATLAB fsolve routine is used in Phase I of the solving strategy.The fsolve routine converges from the initial guess presented in Table 3 to the solution x inf* = [0.0574−0.0233 0.0565], which corresponds to the solution x * = [537.16410.4739 536.6157] in the unmapped space.By tracking the homotopy path determined by Eq. (3) forward and backward with respect to the homotopy parameter, the homotopy path illustrated in Fig. 10 is formed.
The continuous homotopy path that is formed passes through all the feasible solutions of the cooled CSTR model.The solutions are presented in Table 4.The solutions are the same as summarized in Shacham et al. [52].

Adiabatic CSTR with consecutive reactions
The following adiabatic CSTR has been presented in Seader et al. [39] and can also be found in the test problem library accessed through the website: http://www.polymath-software.com/library[53].The system consists of an adiabatic CSTR wherein two reactions take place.The first reaction is catalytic and irreversible: and the second reaction is reversible and non-catalytic: The rate expression [kmol m -3 s -1 ] of the first reaction is , and the rate expression [kmol m -3 s -1 ] of the second reversible reaction is The kinetic [s -1 ] and adsorption equilibrium [m 3 kmol -1 ] terms as a function of temperature are where T is in Kelvins and R 8.314 J mol -1 K -1 .The steadystate mass balance equations for the CSTR reactor in terms of the residence time, Θ, can be written as and the steady-state enthalpy balance, taking into account the sensible enthalpy and reaction enthalpies, as As will be noticed, the equation set to be solved constitutes a set of four equations (Eqs.( 34)-( 37)) and four variables, if the reactor inlet conditions (c A0 ,c B0 , c C0 and T 0 ) and residence time Θ have been specified.Thus, the unknown variables are the reactor effluent quantities (c A ,c B and c C ) and the reactor temperature T. The model of the adiabatic CSTR reactor, i.e.Eqs. ( 34)- (37), may have multiple solutions depending on the given specifications.Use of the specifications c A0 = 3 kmol m -3 , c B0 = c C0 = 0, T 0 = 298 K and Θ = 300 s results in as many as five feasible steady-state solutions.
All the problem variables have a physical meaning.In addition, they are feasible only when they are positive.The domain boundary values and initial guess utilised are illustrated in Table 5.   3), all the five solutions listed in Table 6 are approached on a single continuous path.In addition, as Fig. 11 illustrates, the path also intercepts the θ = 0 plane at two additional points.Because these points lie inside the bounding zone, they do not represent valid solutions of the adiabatic CSTR system model and are thus not taken into account in the set of feasible solutions.The additional points do not indicate that the problem would have feasible or unfeasible solutions outside the defined problem domain or inside the bounding zone.Rather, the additional points are a result of the mathematical formulation of Eq. ( 3).
The probability of having feasible solutions inside the problem bounding zone can be further decreased by using higher absolute numerical values for the lower, l i inf , and upper, inf i u , inner boundaries.In addition, the capability of the method to approach feasible solutions very close to the problem domain boundaries b i min and b i max is improved.However, because the homotopy path is tracked in the mapped variable space in the applied solving method, the length of the path would inevitably increase, thus increasing the solving time.Knowledge about the behaviour of the system can be exploited in order to select numerical values for the inner domain boundaries so that the method would be able to approach all the feasible solutions of the problem robustly but still restrict the solving time from becoming unnecessarily high.As can be seen in Table 6, solutions 1 and 2 are located close to the domain boundary with respect to c C .Let us examine the importance of the correct numerical evaluation of the term in Eq. ( 3), which is the main difference between the previous attempt of Malinen et al. [32] and the current work to formulate a robust solution strategy for chemical engineering systems.Let us assume that we evaluate the Jacobian close to the domain boundary, i.e. at solution 1, following the formulation of Malinen et al. [32] instead of the centre of the feasible variable domain, i. has a considerable effect on the course of the homotopy path.Depending on the evaluation location, the homotopy path crosses the domain boundaries with respect to the homotopy parameter or c C inf .Luckily, also when applying solution 1 as the starting point, all the feasible solutions are found.It is also important to acknowledge that when applying solution 1 as the starting point, the domain boundary crossing occurs in a similar fashion as in the proposed strategy, i.e. the path interception occurs when the bounding is in force for all variables.In addition, the path intercepts the Thus, even in this case, we can continue path tracking at (θ, x inf ) in the same direction from (−θ, −x inf ).
The main reason for the differences between the homotopy paths of different f (x ' ) inf 0 evaluation locations is the numerical procedure applied in the evaluation.In this study Jacobian was evaluated within the fsolve routine of Matlab.In the investigated case, the evaluated Jacobian pattern is different when evaluated at solution 1 rather than at the centre of the domain.The Jacobian patterns are for the proposed strategy and the evaluation at solution 1 as proposed in Malinen et al. [32], respectively.The Jacobian patterns differ with respect to one element.The element at (2,3) is evaluated zero at solution 1, whereas evaluation of the Jacobian with the proposed solving strategy results in the correct Jacobian pattern, where the element is non-zero.
As a consequence of the incorrect Jacobian pattern evaluation at solution 1, the term f ( does not compensate for the annihilation of f(x inf ), thus resulting in a path that is not bounded with respect to the mapped variables, x inf , as shown in Fig. 12.The eventual outcome of this is that the solution is interrupted as the path enters the non-physical variable domain.Thus, the shortcoming of the method proposed in Malinen et al. [32] can be tackled by appropriate selection of the Jacobian evaluation location.However, it is worth noting that the incorrect Jacobian pattern evaluation may occur anywhere within the variable domain, but applying evaluation in the centre of the domain increases the probability of yielding  where the path intercepts the domain boundaries.
the correct Jacobian pattern.Hence, it may be useful to verify that we have evaluated the Jacobian pattern correctly by evaluating it in multiple locations along the variable domain before starting the solving.Nevertheless, assuming that the Jacobian pattern is correct and the non-zero elements of the Jacobian have numerically significant values, robust solving of multiple solutions for CSTR systems is enabled.

Conclusions
In the present paper, a solving strategy has been presented, which can be used for the robust determination of all the feasible solutions of CSTR systems under given specifications, or alternatively to provide the information that the problem has no multiple solutions, presuming only that one feasible solution of the system is available.With this respect, the strategy complements the existing NAE solving algorithms by offering a robust problem-independent tool for the determination of multiple solutions.
The proposed solving strategy utilizes the previously presented Newton homotopy-based method [32].All the feasible solutions of the studied CSTR models were obtained along a single homotopy branch without exiting the feasible variable domain.Thus, the solving method applied in the proposed solving strategy is highly suitable for chemical engineering problem solving, where problem variables must have physically meaningful values.In addition, it was shown that the robustness of the method presented in Malinen et al. [32] can be improved by evaluating the Jacobian pattern in the centre of the domain rather at a solution near the domain boundary.On the other hand, the value of M has a considerable effect on the length of the homotopy path.Thus, the selection of the value of M should be investigated in future work.
As a whole, the applicability and performance of the proposed solving strategy should also be considered for other chemical engineering systems, such as various separation systems and reactive systems with several phases in thermodynamic equilibrium, in future research.

Nomenclature
max m in Variables mapping requires that the maximum, b i max , and minimum, b i min , values are given for every problem variable in the unmapped variable space.The values can be realized as natural or artificial domain boundary values.For example, in the case of mole fractions, the natural selection for the domain boundary values is b σ min = 0 and b i max = 1, whereas in the case of molar flows or concentrations the upper domain boundary value is strongly case-sensitive and must be set artificially or based on case-dependent information.

Fig. 2
Fig. 2 Solving strategy for the determination of all the feasible solutions of a CSTR system at the steady state.
enable the investigation of number of real solutions, into the form

Fig. 3
Fig. 3 BioCSTR system solutions determined with the proposed solving strategy as a function of the Damköhler number.The markers indicate solutions found for the BioCSTR model at different Damköhler numbers.All the feasible solutions were obtained regardless of the value of Da.Also, all the unfeasible real solutions were located in all cases by changing min b σ from zero to −0.5.

Fig. 4
Fig. 4 Homotopy path in the mapped variable space with the solving strategy of the present study.Starting point of the homotopy path (○), solutions (×), and points (•) where the path intercepts the θ = −1 and θ = 1 planes.

9 Fig. 5
Fig. 5 The course of the homotopy path in the BioCSTR model solving with different M parameter values a) throughout the path course and b) in the vicinity of the solutions.The starting point of the homotopy path is indicated by (○) and the solutions by (×).Da = 1.19, ω = 0.00356 and ε = 2.53.
) have been normalized by the inlet concentration and inlet temperature, the reactant concentration c A cannot obtain values above 1 and below 0. Therefore it is justified that b c A max = 1 and b c A min = 0 .Even though T cannot obtain values below 0, values above 1 are feasible when an exothermal reaction takes place in an adiabatic reactor.In this case, the maximum value T can be determined by evaluating the temperature, which the reactor reaches for complete conversion of A. However, here the domain boundary values for T as b T max .= 1 5 and b T min .= 0 5 are applied for illustration purposes of the solving progress,

Fig. 6
Fig. 6 Homotopy path in the mapped variable space with the method of the present study by using the parameters: Q/V R = 8.7, β = 0.15 and γ = 30.Starting point of the homotopy path (○) and points (•)

Fig. 7 Fig. 8
Fig.7 The adiabatic CSTR system solutions determined with the proposed solving strategy as a function of the Damköhler number defined in Eq. (15).The meshes indicate the location of the uniqueness-multiplicity boundaries.

Fig. 9 A
Fig. 9 A schematic representation of a cooled CSTR with an exothermal reaction.

Fig. 10
Fig. 10 Homotopy path with respect to the homotopy parameter and inf A c with the method of the present study.Starting point of the homotopy path (○), solutions (×), and points (•) where the path intercepts the θ = −1 and θ = 1 planes.

Fig. 11
Fig. 11 Homotopy path with respect to the homotopy parameter, c A inf and c

e. x 0 0 inf=
, as proposed in the solving strat- egy of this study.The change in the homotopy path course is shown in Fig.12.It is evident based on Fig.12that the selected location of evaluating f c C inf = −20 and c C inf = 20 2 planes at the points where the homotopy parameter, c A inf and c C inf have the same values but with opposite signs.

Fig. 12
Fig.12 Homotopy paths with respect to the homotopy parameter, c A inf

Table 4
The feasible solutions of the cooled CSTR model.

Table 3
Domain boundary and initial guess values.

Table 5
Domain boundary values and initial guesses.

Table 6
Feasible solutions for the adiabatic CSTR.