Study on the Nonlinear Dynamics of the Continuous Stirred Tank Reactors

: Chemical processes often exhibit nonlinear dynamics and tend to generate complex state trajectories, which present challenging operational problems due to complexities such as output multiplicity, oscillation, and even chaos. For this reason, a complete knowledge of the static and dynamic nature of these behaviors is required to understand, to operate, to control, and to optimize continuous stirred tank reactors (CSTRs). Through nonlinear analysis, the possibility of output multiplicity, self-sustained oscillation, and torus dynamics are studied in this paper. Speciﬁcally, output multiplicity is investigated in a case-by-case basis, and related operation and control strategies are discussed. Bifurcation analysis to identify di ﬀ erent dynamic behaviors of a CSTR is also implemented, where operational parameters are identiﬁed to obtain self-oscillatory dynamics and possible unsteady-state operation strategy through designing the CSTR as self-sustained periodic. Finally, a discussion on codimension-1 bifurcations of limits cycles is also provided for the exploration of periodic forcing on self-oscillators. Through this synergistic study on the CSTRs, possible output multiplicity, oscillatory, and chaotic dynamics facilitates the implementation of novel operation / control strategies for the process industry.


Introduction
It is well recognized that the natural world is complex, presenting an environment with chemical and ecological networks that interact on a global scale. Elements of many such systems always exhibit nonlinear dynamics, and continuous stirred tank reactors (CSTRs) are typical process units that exhibit nonlinear dynamics, presenting challenging operational problems due to complex behavior such as output multiplicities, oscillations, and even chaos [1]. For this reason, a complete knowledge of the static and dynamic nature of these behaviors is required to understand, to operate, to control, and to optimizea CSTR [2].
The presence of multiple steady states necessitates handy design of an efficient control system to regulate the operation states, and information about the transitions from state to state is critical [3]. Challenges would also arise when it is desirable to operate CSTRs under open-loop unstable conditions, where the reaction rate may yield good productivity while the reactor temperature is proper to prevent side reactions or catalyst degradation. Hence, it is important to consider the open/closed loop stability of the CSTRs, with analysis of output multiplicity providing practical guidelines for process design and operation.
Moreover, many studies show that CSTRs may exhibit a rich behavior in dynamic phenomena [4,5], with oscillation being the factor that has been subjected to intense research activity by both

Model Formulation and Analysis
To investigate the dynamics of the lumped system, we introduce a simple first-order, irreversible reaction A→B occurring at a jacket CSTR; the governing equations for mass and energy balance are provided as follows: The jacket outlet temperature Tj represents the characteristic temperature of the jacket for heat convection, and is assumed uniform in the jacket. This holds true when the jacket volume V is huge and well mixed and the thermal inertia of the metal walls is considered negligible.
For analytical purposes, further simplification steps could be applied to make the CSTR model two-dimensional, which is satisfied only when the dynamics related to the jacket temperature are much faster than that related to the reactor temperature. Then, the jacket temperature time constant is negligible and Tj can be derived by setting the right-hand side of the third formula at Equation (1) to zero: T j = q j ρ j c pj T j0 + U a T U a + q j ρ j c pj (2) Substituting Equation (2) into Equation (1), and the dimensionless form of the model is given by Vρc p (1+Kq j ) + kZ 1 where Z 1 = C A /C A0 ; Z 2 = Tρc p /(−∆H)C A Z f 0 = T f 0 ρc p /(−∆H)C A ; Z j0 = T j0 ρc p /(−∆H)C A K = ρ j c pj /U a ; k = k 0 exp(−E/RT) Setting the right-hand side of Equation (3) to zero obtains the equilibrium solutions: q(1+Kq j )ρc p /(qZ f 0 +kV+kVZ f 0 )+q j KU a Z j0 (q+kV) (q(1+Kq j )ρc p +q j KU a )(q+kV) (4) Note that the solutions in Equation (4) may not be stable and stability analysis could be performed using the Jacobian matrix of Equation (3), which is provided as follows: Since all eigenvalues of Equation (5) need to be negative to meet the stability criterion, a 11 a 22 − a 12 a 21 > 0 a 11 + a 22 < 0 (6) Other than the conditions prescribed in Equation (6), the CSTR is unstable. Singularity would cause complex dynamics to emerge, e.g., the necessary condition for self-oscillation, which is provided as follows: a 11 a 22 − a 12 a 21 > 0 a 11 + a 22 = 0 However, a general analytical investigation of instabilities of the CSTR would be infeasible because complexity of the terms in Equation (5). The presence of multiple solutions of a CSTR model indicates that each solution may have potentially different stability properties. One can expect that a multitude of coexisting steady states with a diversity of stability properties could lead to a complex global behavior. Moreover, when the rigorous three-dimensional model is considered, the analysis could be complicated further. Therefore, numerical bifurcation analysis is adopted to study the dynamics of CSTRs in a case-by-case manner, and on the basis of the analysis, related design, operation, and control issues are detailed.

Numerical Bifurcation Analysis Tools
For an effective illustration of the nonlinear dynamics in CSTRs, we sketch basics of bifurcation theory in this section. Consider a lumped system represented by ordinary differential equations (ODEs) depending on a parameter vector α, where f : R n × R→R n is smooth and nonlinear. The stability criterion is based on the equilibrium manifold f (x) = 0 with given parameter vector α. Numerical bifurcation analysis is implemented through (1) draw equilibriu mcurves with variation ofα; a numerical continuation algorithm, e.g., arch-length continuation, is available for solving such a problem; (2) detect bifurcation points along with these curves, where the qualitative or topological structure may change, leading the manifold to also change significantly. Limit point is correlative to the production of output multiplicity, which is identified as the Jacobian matrix being singular: where as Hopf bifurcation emerges when a pair of pure imaginary eigenvalues are detected, and Hopf bifurcation is related to the production of self-oscillatory behavior. A Hopf point can be viewed as amplitude zero limit cycle, and periodic forcing of a self-oscillator can bifurcate out complex dynamic trajectories such as torus. Hence, codimension-1 bifurcation of limit cycles, i.e., Hopf bifurcation of the limit cycles, is adopted to identify the generation of complex dynamics in CSTRs. Stemming from the fact that a rich dynamic behavior is observed in CSTRs, which disturbs the good development of the process, a complete knowledge of static and dynamic behavior is required to understand, to operate, to control, and to optimize the CSTRs. In the following sections, bifurcation analysis is implemented on the CSTRs, and for this purpose, we propose the analysis strategy shown in Figure 1. On the basis of the analytics, we discuss topics on periodic operation and optimization of the processes.

Numerical Bifurcation Analysis Tools
For an effective illustration of the nonlinear dynamics in CSTRs, we sketch basics of bifurcation theory in this section. Consider a lumped system represented by ordinary differential equations (ODEs) depending on a parameter vector α, where f: R n × R→R n is smooth and nonlinear. The stability criterion is based on the equilibrium manifold f(x) = 0 with given parameter vector α. Numerical bifurcation analysis is implemented through (1) draw equilibriu mcurves with variation ofα; a numerical continuation algorithm, e.g., arch-length continuation, is available for solving such a problem; (2) detect bifurcation points along with these curves, where the qualitative or topological structure may change, leading the manifold to also change significantly.
Limit point is correlative to the production of output multiplicity, which is identified as the Jacobian matrix being singular: where as Hopf bifurcation emerges when a pair of pure imaginary eigenvalues are detected, and Hopf bifurcation is related to the production of self-oscillatory behavior. A Hopf point can be viewed as amplitude zero limit cycle, and periodic forcing of a self-oscillator can bifurcate out complex dynamic trajectories such as torus. Hence, codimension-1 bifurcation of limit cycles, i.e., Hopf bifurcation of the limit cycles, is adopted to identify the generation of complex dynamics in CSTRs. Stemming from the fact that a rich dynamic behavior is observed in CSTRs, which disturbs the good development of the process, a complete knowledge of static and dynamic behavior is required to understand, to operate, to control, and to optimize the CSTRs. In the following sections, bifurcation analysis is implemented on the CSTRs, and for this purpose, we propose the analysis strategy shown in Figure 1. On the basis of the analytics, we discuss topics on periodic operation and optimization of the processes.

Control of Output Multiplicity
One interesting feature of nonlinear systems is their multiplicity behavior. For the case of input multiplicity [21], different feed flow rate q may lead to the same outcomes, which will affect CSTR control strategy because the "zero dynamics" are unstable. Hence, CSTRs are generally not controlled using q as the manipulate variable, and to attenuate the disturbances of q, feedward control is frequently applied. Output multiplicity in the CSTR is because the curve of generated heat and the line of removed heat versus temperature can intersect with multiple points.
A necessary condition to determine output multiplicity is provided [22], but direct employment of this condition to analyze a typical process is difficult because the maximum slope for the heat generation function is always difficult to solve explicitly. On the other hand, numerical continuation offers a method to analyze output multiplicity case-by-case. Taking the parameters in Table 1, and for the given input flow rates, multiple steady outputs may emerge. As shown in Figure 2, when q j = 10 m 3 /s, there are three equilibrium points-P1, P2, and P3-intersecting with the vertical line q = 16 m 3 /s. Hence, the necessary condition for output multiplicity is tangent of any curve to satisfy Table 1. Recipe from Douglas's model [22].

Values Parameter Values
Processes 20xx, x, x FOR PEER REVIEW 5 of 19

Control of Output Multiplicity
One interesting feature of nonlinear systems is their multiplicity behavior. For the case of input multiplicity [21], different feed flow rate q may lead to the same outcomes, which will affect CSTR control strategy because the "zero dynamics" are unstable. Hence, CSTRs are generally not controlled using q as the manipulate variable, and to attenuate the disturbances of q, feedward control is frequently applied. Output multiplicity in the CSTR is because the curve of generated heat and the line of removed heat versus temperature can intersect with multiple points.
A necessary condition to determine output multiplicity is provided [22], but direct employment of this condition to analyze a typical process is difficult because the maximum slope for the heat generation function is always difficult to solve explicitly. On the other hand, numerical continuation offers a method to analyze output multiplicity case-by-case. Taking the parameters in Table 1, and for the given input flow rates, multiple steady outputs may emerge. As shown in Figure 2, when qj = 10 m 3 /s, there are three equilibrium points-P1, P2, and P3-intersecting with the vertical line q = 16 m 3 /s. Hence, the necessary condition for output multiplicity is tangent of any curve to satisfy  In the above case, as qj decreases, a bifurcation point (qj ≈ 9.8 m 3 /s) appears and a collision of two steady statesoccurs, witha saddle-node bifurcation occurring when the critical equilibrium of the system has onezero eigenvalue. Detection of the limit points analytically could be difficult, but through tools of numerical bifurcation analysis, theyare easily detected, e.g., qj ≈ 10.0 m 3 /s in Figure 2 is identified with three output states. Table 1. Recipe from Douglas's model [22].

Parameter
Values In the above case, as q j decreases, a bifurcation point (q j ≈ 9.8 m 3 /s) appears and a collision of two steady statesoccurs, witha saddle-node bifurcation occurring when the critical equilibrium of the system has onezero eigenvalue. Detection of the limit points analytically could be difficult, but through tools of numerical bifurcation analysis, theyare easily detected, e.g., q j ≈ 10.0 m 3 /s in Figure 2 is identified with three output states.
When multiple output states are detected in the CSTR, controlling the working condition as desired is critical. A typical control strategy is to remove the multiplicity behavior between reactor temperature and cooling jacket flow rate by increasing the reactor feed temperature. As is shown in Figure 2, when the cooling water q j is below 9.8 m 3 /s, no limit point is detected, which indicates that no output multiplicity is generated. However, the dynamics may offer more information than the traditional steady state control. From the phase plot (q j = 10 m 3 /sand q = 16 m 3 /s), the convergent ridge (red arrow) of P2 separates the whole (Z1, Z2) domain into opposite trajectory directions, where the left-bottom zone is attracted to P1 point, but not all the initial points from the top-right zone converges to the P3 point because the red arrow almost crosses the P3 point. The stale zone that is attracted to P3 point is limited, and one can identify P3 stable criterion by the following equation: From Figure 3, the condition satisfying Equation (11) obtains a very small stable window, and there is a high possibility to converge at the P1 point. When q = 16 m 3 /s and holds constant, q j shifts across 9.8 m 3 /s (the limit point) can realize operation change between P1 and P3, approximately. Different q j would lead to different steady state, and the industry may expect different operational states at specific circumstances, i.e., failure to control the reactor properly would restrain it at P1, which has a lower conversion rate thanP3.Therefore, the following split-ranging control strategy shown in Figure 4 is adopted-at the start-up stage, steam valve B is opened to heat the reactor to around the P3 point. When the reaction heat is released, valve B is gradually closed and valve A is opened to remove the reaction heat, while at the shutdown stage, more heat is removed to force the working point to shift to P1. Therefore, on the basis of the phase dynamics, the CSTR operation/control could work at multiple states for different production load specification. When multiple output states are detected in the CSTR, controlling the working condition as desired is critical. A typical control strategy is to remove the multiplicity behavior between reactor temperature and cooling jacket flow rate by increasing the reactor feed temperature. As is shown in Figure 2, when the cooling water qj is below 9.8 m 3 /s, no limit point is detected, which indicates that no output multiplicity is generated. However, the dynamics may offer more information than the traditional steady state control. From the phase plot (qj = 10 m 3 /sand q = 16 m 3 /s), the convergent ridge (red arrow) of P2 separates the whole (Z1, Z2) domain into opposite trajectory directions, where the left-bottom zone is attracted to P1 point, but not all the initial points from the top-right zone converges to the P3 point because the red arrow almost crosses the P3 point. The stale zone that is attracted to P3 point is limited, and one can identify P3 stable criterion by the following equation: Figure 3, the condition satisfying Equation (11) obtains a very small stable window, and there is a high possibility to converge at the P1 point. When q = 16 m 3 /s and holds constant, qj shifts across 9.8 m 3 /s (the limit point) can realize operation change between P1 and P3, approximately. Different qj would lead to different steady state, and the industry may expect different operational states at specific circumstances, i.e., failure to control the reactor properly would restrain it at P1, which has a lower conversion rate thanP3.Therefore, the following split-ranging control strategy shown in Figure 4 is adopted-at the start-up stage, steam valve B is opened to heat the reactor to around the P3 point. When the reaction heat is released, valve B is gradually closed and valve A is opened to remove the reaction heat, while at the shutdown stage, more heat is removed to force the working point to shift to P1. Therefore, on the basis of the phase dynamics, the CSTR operation/control could work at multiple states for different production load specification. When multiple output states are detected in the CSTR, controlling the working condition as desired is critical. A typical control strategy is to remove the multiplicity behavior between reactor temperature and cooling jacket flow rate by increasing the reactor feed temperature. As is shown in Figure 2, when the cooling water qj is below 9.8 m 3 /s, no limit point is detected, which indicates that no output multiplicity is generated. However, the dynamics may offer more information than the traditional steady state control. From the phase plot (qj = 10 m 3 /sand q = 16 m 3 /s), the convergent ridge (red arrow) of P2 separates the whole (Z1, Z2) domain into opposite trajectory directions, where the left-bottom zone is attracted to P1 point, but not all the initial points from the top-right zone converges to the P3 point because the red arrow almost crosses the P3 point. The stale zone that is attracted to P3 point is limited, and one can identify P3 stable criterion by the following equation: Figure 3, the condition satisfying Equation (11) obtains a very small stable window, and there is a high possibility to converge at the P1 point. When q = 16 m 3 /s and holds constant, qj shifts across 9.8 m 3 /s (the limit point) can realize operation change between P1 and P3, approximately. Different qj would lead to different steady state, and the industry may expect different operational states at specific circumstances, i.e., failure to control the reactor properly would restrain it at P1, which has a lower conversion rate thanP3.Therefore, the following split-ranging control strategy shown in Figure 4 is adopted-at the start-up stage, steam valve B is opened to heat the reactor to around the P3 point. When the reaction heat is released, valve B is gradually closed and valve A is opened to remove the reaction heat, while at the shutdown stage, more heat is removed to force the working point to shift to P1. Therefore, on the basis of the phase dynamics, the CSTR operation/control could work at multiple states for different production load specification.

Periodic Exciting of the Inputs
Conventionally, a process is designed and operated at a desired steady state, and oscillations are taken as disturbances, which are not desirable; thus, attempts have been made to eliminate sources of variation by modifying process design specifications to ensure stability, or by adding a control system or a surge tank [23][24][25][26][27][28]. Recently, topics about designing for periodic operation began to resonate within the process systems engineering community, and nonlinearity of process systems may allow for τ-time averaged performance enhancement with periodic forcing of external parameters. The operation of CSTRs with time-varying feed conditions is an appealing prospect for the industry due to its relative ease in implementation. One can expect an optimal periodic operation problem that aims to maximize a scalar objective function as follows: where g(x, u) is the performance index, such as the profit. Let there be a stable state x 0 for a stationary input u 0 , and the corresponding scalar objective function is J 0 , which is equal to g(x 0 , u 0 ). Set the input in the form of where δu represents any continuous zero-mean τ-period perturbation vector, and the time-averaged objective function is given as J. If J > J 0 , then, the forced periodic operation performance surpasses that of the stationary operation. Essentially, the input perturbation δu, the dynamic system Equation (8), and the performance index g(x,u) all together determine performance of the periodic operation. The case introduced is the cascade-parallel reaction occurring in the adiabatic CSTR; industrial interest is to maximize the intermediate P, as well as to minimize the by-production R. The reaction kinetics is presented as follows: Assume the reaction volume is constant dV/dt = 0, the lumped model is provided as follows: where the flow rate is q, the concentration is C, and the subscript F represents the input condition. The objective is to maximize P production, but the separation cost of by-product R and reactants cost A and B are also taken into account, which are scaled by σ, σ 1 , and σ 2 , respectively.
The steady state performance is correlative to the extent of the reaction, η i = r i V (i = 1,2). The corresponding objective J 0 can be given as follows, where θ B is the mole ratio of B to A: When the design specifications σ, σ 1 , σ 2 , q A C AF , and q B C BF are held constant, the steady state operational problem becomes maximizing J 0 (criterion: ∂J 0 /∂x = 0, ∂ 2 J 0 /∂x 2 < 0), which is provided as follows: [29], then, the optimal outputs are C A0 = 0.3270, C B0 = 1.238, C P0 = 0.2083.
From an unsteady state operation viewpoint, one may inspect the overall performance if the input concentration C A and C B are excited by, say, cosine disturbances, where i 2 = −1, ω is the angular frequency, θ is the phase shift between two cosine disturbances, and Am is the amplitude. Therefore, the optimal periodic operation problem is to find proper ω and Am such that the τ-average objective function Equation (15) is maximized. Through Laplace-Borel [13] transform analysis, relation of the output concentrations with varying parameters are presented as follows: When the separation cost σ is neglected, Equation (15) retreats to the simple τ-time averaged production concentration δP. It is unfortunate that in the current case, periodic operation causes δP to decrease. As is shown in Figure 5, for different θ, J decreases with ω decrease. Hence, surge tanks are inserted before the reactor to eliminate possible oscillations of C A and C B .
From an unsteady state operation viewpoint, one may inspect the overall performance if the input concentration CA and CB are excited by, say, cosine disturbances, where i 2 = −1, ω is the angular frequency, θ is the phase shift between two cosine disturbances,and Am is the amplitude. Therefore, the optimal periodic operation problem is to find proper ω and Am such that the τ-average objective function Equation (15) When the separation cost σ is neglected, Equation (15) retreats to the simple τ-time averaged production concentration δP. It is unfortunate that in the current case, periodic operation causes δP to decrease. As is shown in Figure 5, for different θ, J decreases with ω decrease. Hence, surge tanks are inserted before the reactor to eliminate possible oscillations of CA and CB.

Generation of Self-Sustained Oscillations
It is well known that a CSTR without control can generate self-sustained oscillations, and outputs from such a process vary periodically even if there is no disturbance in input streams. It is interesting to analyze the conditions under which oscillations arise in a CSTR, and Hopf bifurcations can give birth to limit cycles, when the equilibrium changes stability via a pair of purely imaginary eigen values. Suppose f(x, u) is a general vector field on R n and satisfies f (0, u 0 ) = 0 for some constant input u 0 , and J(u 0 ) denotes the associated Jacobian matrix evaluated at x = 0. For a neighborhood of u 0 , the eigenvalues of J(u) split invariantly into three disjoint eigen spaces W s , W c , and W u containing stable, central, and unstable eigenvalues, respectively. Hopf bifurcation is generated by a pair of imaginary eigenvalues (c = ±iω(u)), and by central manifold theorem, f(x, u) can be reduced to a two-dimensional equivalence along with the zero manifold. Therefore, the study ofHopf bifurcations can be reduced to a two-dimensional form, .
where α(u 0 ) = 0 and the residual is the 2-norm. Equation (10) can be transformed to the normal form Hopf bifurcation as follows: .
where β(u = 0) = 0. One can verify that the quadratic terms will not affect the structural properties of the dynamic system and can be truncated off. By introducing the complex variables z = y 1 + iy 2 and z* = y 1 − iy 2 , Equation (11) is re-formulated as the following complex form: Using z = ρe iϕ and substituting to Equation (12) gives the polar form: The first formula of Equation (13) has the equilibrium point ρ = 0 for all values of β. The equilibrium is linearly stable if β < 0; it remains stable at β = 0 but nonlinearly, and the rate of solution convergence to zero is no longer exponential; for β > 0, the equilibrium becomes linearly unstable. Moreover, there is an additional stable equilibrium point ρ 0 (β) = sqrt(β) for β > 0. The Hopf bifurcation can be supercritical or subcritical due to the stability criterion of the limit cycles, which is evaluated by the Lyapunov function, as is detailed in Appendix A.
Previous studies have shown that a self-sustained oscillation is generated from the Oxo reaction system [30], and this process can be simplified as a first-order irreversible reaction system with the modeling parameters displayed in Table 2. A holistic picture of the possible bifurcation mechanisms that the CSTR model may exhibit is best achieved by showing both the loci of the saddle nod points and the Hopf points. The parameters q and q j are easy to manipulate in practice-a two-parameter bifurcation diagram is provided with varying q and q j . A two-parameter bifurcation diagram with varying q and q j is given in Figure 6, where two loci of Hopf curves are detected and shown as the dash line. The Bogdanov-Taken (BT) point bifurcates the solid curve with the Hopf branch, which is labeled as dashed curve 2, and the solid curve is the saddle node branch. Together with another Hopf curve 1, these three curves separate the whole region into six sub-sections, each exhibiting different global behavior. Within this diagram, the operational point P mentioned above falls into section B, as represented by the large sized dot, meaning that a unique stable limit cycle is generated. A two-parameter bifurcation diagram with varying q and qj is given in Figure 6, where two loci of Hopf curves are detected and shown as the dash line. The Bogdanov-Taken (BT) point bifurcates the solid curve with the Hopf branch, which is labeled as dashed curve 2, and the solid curve is the saddle node branch. Together with another Hopf curve 1, these three curves separate the whole region into six sub-sections, each exhibiting different global behavior. Within this diagram, the operational point P mentioned above falls into section B, as represented by the large sized dot, meaning that a unique stable limit cycle is generated. Figure 6. Two-parameter bifurcation diagram of the Oxo reaction model. BT stands for a Bogdanov-Taken point, and CP stands for a cusp point. Any points in section A generate a unique stable point thatcoverts to the low conversion region; points in section E generate three equilibriums, the high and low conversion ones are stable while the middle one is a saddle; points in section B generate one stable limit cycle; points in section C generate a unique stable point thatcoverts to the low conversion region; points in section D generate three equilibriums, the high conversion region is a stable limit cycle, then a saddle, and a stable point at the low conversion region; points in section F generate a saddle that is incorporated by a limit cycle.
From the analysis diagram provided in Figure 6, one can find that designing at q = 0.004923 m 3 s −1 and qj = 0.01 m 3 s −1 generates oscillatory dynamics. If the oscillation is not desirable, one can redesign the operation point to, say, q = 0.01 m 3 s −1 and qj = 0.01 m 3 s −1 , which falls into section E. If self-oscillation is on demand, i.e., the steady-state operating point is set as q = 0.004923 m 3 s −1 and qj = 0.01 m 3 s −1 , a stable limit cycle with frequency 0.00167 s −1 is generated, as is shown in Figure 7 (black line). Since points in section E generate three equilibriums, and the high and low conversion ones are stable while the middle one is a saddle, the phase portrait obtained is similar to the subplot in Figure 2. One can implement the strategy provided in Figure 4 to control the CSTR at a high conversion equilibrium point. Because the mean value of the process output of a process operating at oscillatory state will not be the same as that of a steady state process output at the mean values of operating variables, one can also utilize the oscillatory dynamics for process intensification purpose, which is incorporated in the following section.  Figure 6. Two-parameter bifurcation diagram of the Oxo reaction model. BT stands for a Bogdanov-Taken point, and CP stands for a cusp point. Any points in section A generate a unique stable point thatcoverts to the low conversion region; points in section E generate three equilibriums, the high and low conversion ones are stable while the middle one is a saddle; points in section B generate one stable limit cycle; points in section C generate a unique stable point thatcoverts to the low conversion region; points in section D generate three equilibriums, the high conversion region is a stable limit cycle, then a saddle, and a stable point at the low conversion region; points in section F generate a saddle that is incorporated by a limit cycle.
From the analysis diagram provided in Figure 6, one can find that designing at q = 0.004923 m 3 s −1 and q j = 0.01 m 3 s −1 generates oscillatory dynamics. If the oscillation is not desirable, one can redesign the operation point to, say, q = 0.01 m 3 s −1 and q j = 0.01 m 3 s −1 , which falls into section E. If self-oscillation is on demand, i.e., the steady-state operating point is set as q = 0.004923 m 3 s −1 and q j = 0.01 m 3 s −1 , a stable limit cycle with frequency 0.00167 s −1 is generated, as is shown in Figure 7 (black line). Since points in section E generate three equilibriums, and the high and low conversion ones are stable while the middle one is a saddle, the phase portrait obtained is similar to the subplot in Figure 2. One can implement the strategy provided in Figure 4 to control the CSTR at a high conversion equilibrium point. Because the mean value of the process output of a process operating at oscillatory state will not be the same as that of a steady state process output at the mean values of operating variables, one can also utilize the oscillatory dynamics for process intensification purpose, which is incorporated in the following section. A two-parameter bifurcation diagram with varying q and qj is given in Figure 6, where two loci of Hopf curves are detected and shown as the dash line. The Bogdanov-Taken (BT) point bifurcates the solid curve with the Hopf branch, which is labeled as dashed curve 2, and the solid curve is the saddle node branch. Together with another Hopf curve 1, these three curves separate the whole region into six sub-sections, each exhibiting different global behavior. Within this diagram, the operational point P mentioned above falls into section B, as represented by the large sized dot, meaning that a unique stable limit cycle is generated. Figure 6. Two-parameter bifurcation diagram of the Oxo reaction model. BT stands for a Bogdanov-Taken point, and CP stands for a cusp point. Any points in section A generate a unique stable point thatcoverts to the low conversion region; points in section E generate three equilibriums, the high and low conversion ones are stable while the middle one is a saddle; points in section B generate one stable limit cycle; points in section C generate a unique stable point thatcoverts to the low conversion region; points in section D generate three equilibriums, the high conversion region is a stable limit cycle, then a saddle, and a stable point at the low conversion region; points in section F generate a saddle that is incorporated by a limit cycle.
From the analysis diagram provided in Figure 6, one can find that designing at q = 0.004923 m 3 s −1 and qj = 0.01 m 3 s −1 generates oscillatory dynamics. If the oscillation is not desirable, one can redesign the operation point to, say, q = 0.01 m 3 s −1 and qj = 0.01 m 3 s −1 , which falls into section E. If self-oscillation is on demand, i.e., the steady-state operating point is set as q = 0.004923 m 3 s −1 and qj = 0.01 m 3 s −1 , a stable limit cycle with frequency 0.00167 s −1 is generated, as is shown in Figure 7 (black line). Since points in section E generate three equilibriums, and the high and low conversion ones are stable while the middle one is a saddle, the phase portrait obtained is similar to the subplot in Figure 2. One can implement the strategy provided in Figure 4 to control the CSTR at a high conversion equilibrium point. Because the mean value of the process output of a process operating at oscillatory state will not be the same as that of a steady state process output at the mean values of operating variables, one can also utilize the oscillatory dynamics for process intensification purpose, which is incorporated in the following section.

Forced Periodic Inputs on the Self-Sustained Oscillations
From process intensification viewpoint, regulating the set point as periodic may achieve better performance than controlling at steady state, but monitoring delay or process model mismatch may cause interference between manipulate and control variables. Hence, it is necessary to study the dynamics through external forcing of the chemical oscillator.
In current study, q j is excited rather than q on two reasons: (1) for a typical operational problem, varying the reactant inlet flow rate in large scale is not easy to realize in practice; (2) q is linearly correlated with state variable Z 1 in the dynamic system, while nonlinearity is the reason for τ-time averaged performance improvement. On the other hand, varying q j is advisable because coolant is easy to manipulate, and because the flow rate of q j usually determines heat removal from the reactor. For example, a deviation of q j by 10% will have an impact up to 5.5% for the heat removal term in Equation (1), hence havinga direct influence on the process itself.
Taking q j as forced oscillation variable, with the form of shown as follows: where q j0 is the nominal coolant flow rate, A is the forcing amplitude, f is the forcing frequency. A sinusoidal function by itself is an orbit if it is cast to polar coordinate plane. Setting u = sin(2πft + ϕ), where the derivative is du/dt = cos(2πft + ϕ) 2πf, then defining v = cos(2πft + ϕ), for analytical convenience, Equation (22) is written in the formal form of the two dimensional Hopf bifurcation; the sinusoidal forcing of the chemical oscillator is presented as follows: For the system above, the first two formulas in combined forms is a limit cycle while the remaining two forms are the other; how the two orbits correlate with each other to intensify overall behavior of the process system is the focus of this discussion.
As the intrinsic frequency affects the dynamics of the studied process, it can be calculated according to the forcing rotation velocity by casting the normal form of the unforced process to the polar coordinates. The method chosen in this work is to integrate Equation (3) at the nominal state by a fourth-order Runge-Kutta algorithm, and a frequency of 0.00167 s −1 is extracted from its trajectory. Apparently, when the difference between the forcing frequency f and the intrinsic frequency f 0 is large, the influence of an oscillation with a lower frequency would fade away; when f ≈ f 0 , resonance takes place. Neither scenario is expected from the process intensification point of view. Therefore, the forcing frequencies in the range 0.1 < f/f 0 < 1 plus 1 < f/f 0 < 10 are of interest [31], and f/f 0 = 3 would be a preferable choice.
The impact of process inlet flow rate forcing on state variables is expected to be insignificant when the forcing amplitude is very small. Yet reasonably large forcing amplitudes provide an opportunity to modify the dynamics of the process. As shown in Figure 7, with A = 0.006 m 3 /s, f = 0.006 s −1 , the output frequency shifts from the intrinsic frequency to the forcing frequency of 0.006 s −1 . Complex dynamics would emerge because of codimension-1 bifurcations, i.e., bifurcation of the limit cycles, as is shown in Figure 8a periodic doubling (PD) and Figure 8b Neimark-Sacker (NS) bifurcation.
analytical solution of PD and NS bifurcations in Appendix B, and from Equation (A2), it is obvious that μ > 0 causes Equation (A1) to diverge, while for μ< 0 it shrinks to a point. Periodic or pseudoperiodic solutions are obtained when Re(μ) = 0, which is the singular point that may causecodimension-1 bifurcation of limit cycles. NS bifurcation emerges when μ1,2 = 0 ± iφ, where a two-dimensional invariant torus appears, while the fixed point changes stability by a Hopf bifurcation, as is shown in Figure 8b.  Taking amplitude as the bifurcation parameter for Equation (15), the continuation diagram of the limit cycle is provided in Figure 9. A supercritical NS point is generated when the amplitude is about 0.00266 s −1 , with the normal form coefficient of this NS point being negative. One can find the analytical solution of PD and NS bifurcations in Appendix B, and from Equation (A2), it is obvious that µ > 0 causes Equation (A1) to diverge, while for µ< 0 it shrinks to a point. Periodic or pseudo-periodic solutions are obtained when Re(µ) = 0, which is the singular point that may causecodimension-1 bifurcation of limit cycles. NS bifurcation emerges when µ 1,2 = 0 ± iϕ, where a two-dimensional invariant torus appears, while the fixed point changes stability by a Hopf bifurcation, as is shown in Figure 8b. Taking amplitude as the bifurcation parameter for Equation (15), the continuation diagram of the limit cycle is provided in Figure 9. A supercritical NS point is generated when the amplitude is about 0.00266 s −1 , with the normal form coefficient of this NS point being negative. One can find the analytical solution of PD and NS bifurcations in Appendix B, and from Equation (A2), it is obvious that μ > 0 causes Equation (A1) to diverge, while for μ< 0 it shrinks to a point. Periodic or pseudoperiodic solutions are obtained when Re(μ) = 0, which is the singular point that may causecodimension-1 bifurcation of limit cycles. NS bifurcation emerges when μ1,2 = 0 ± iφ, where a two-dimensional invariant torus appears, while the fixed point changes stability by a Hopf bifurcation, as is shown in Figure 8b.  Therefore, when the orbits are on the left-hand side of NS (A > A NS ), the output is a stable limit cycle with the same frequency as the forcing frequency. When the orbits are on the right-hand side of NS (A < A NS ), an invariant two-dimensional torus is generated in the corresponding phase space, as shown in the subplot of Figure 10. The emergence of torus is viewed as coexistence of two cycling with different frequency, but torus is to be avoided from apractical perspective because the dynamics of the torus are complicated, the large area of working states makes it difficult to control, and lower average conversion ratio is obtained in this process. Therefore, forcing amplitude should be larger than the NS point, and application constrains such as operation cost should also be considered to decide proper forcing amplitude. When the periodic forcing is A = 0.006, f = 0.006, better outcome is generated in at least three aspects: (1) the process itself is dominated by forcing frequency, which means as long as the forcing input is controlled, the objective is well controlled, irrelevant of the internal oscillation; (2) the process is intensified by periodic forcing; (3) smaller output amplitude in comparison with the internal oscillation makes the system more robust to various disturbances. Therefore, when the orbits are on the left-hand side of NS (A > ANS), the output is a stable limit cycle with the same frequency as the forcing frequency. When the orbits are on the right-hand side of NS (A < ANS), an invariant two-dimensional torus is generated in the corresponding phase space, as shown in the subplot of Figure 10. The emergence of torus is viewed as coexistence of two cycling with different frequency, but torus is to be avoided from apractical perspective because the dynamics of the torus are complicated, the large area of working states makes it difficult to control, and lower average conversion ratio is obtained in this process. Therefore, forcing amplitude should be larger than the NS point, and application constrains such as operation cost should also be considered to decide proper forcing amplitude. When the periodic forcing is A = 0.006, f = 0.006, better outcome is generated in at least three aspects: (1) the process itself is dominated by forcing frequency, which means as long as the forcing input is controlled, the objective is well controlled, irrelevant of the internal oscillation; (2) the process is intensified by periodic forcing; (3) smaller output amplitude in comparison with the internal oscillation makes the system more robust to various disturbances. To numerically show how the periodic forcing can improve the performance of a process, the performance index calculated as follows: where τ is the forcing period; g(x,u) is the performance function, whose unforced equivalence is represented as ys = g(xs); Z is the state variable vector; and u is the forcing input.  Figure 4, indicating that, for various f values, the dynamics of the process system is about to change when A decreases to NS point. The periodic forcing parameters in Equation (26) are recommended to be A = 0.006 m 3 /s, f = 0.006 s −1 , and better outcome is obtained in at least three aspects: (1) the process itself is dominated by the forcing frequency, which means that the objective is well controlled as long as the forcing input is controlled, regardless of its internal oscillation; (2) the process is intensified by periodic forcing; (3) a smaller output amplitude in comparison with the internal oscillation makes the system more robust to various disturbances. To numerically show how the periodic forcing can improve the performance of a process, the performance index calculated as follows: where τ is the forcing period; g(x,u) is the performance function, whose unforced equivalence is represented as y s = g(x s ); Z is the state variable vector; and u is the forcing input. In this case, the state variable vector is [Z 1 , Z 2 ], with the performance function g(Z(u)) as the outlet concentration Z 1 . For the case where the forcing is given by A = 0.006 m 3 /s, f = 0.006 s −1 , the average outlet concentration is Z 1 forced = 0.067, compared with the unforced average value of Z 1 unforced = 0.078, thus showing an increase in conversion of 16.4%. As shown in Figure 9, the center of the orbit would move left as the forcing amplitude A increases, but A cannot be increased in an unbounded manner, i.e., the jacket volume, coolant control valve, or the reactor itself would restrict the upper limit of A. The same situation happens with the forcing frequency f ; when f is too high, even the feasibility of the process model, especially the assumption of remark 2, is questionable. For different f values, the location of the NS point does not migrate significantly, as shown in Figure 4, indicating that, for various f values, the dynamics of the process system is about to change when A decreases to NS point. The periodic forcing parameters in Equation (26) are recommended to be A = 0.006 m 3 /s, f = 0.006 s −1 , and better outcome is obtained in at least three aspects: (1) the process itself is dominated by the forcing frequency, which means that the objective is well controlled as long as the forcing input is controlled, regardless of its internal oscillation; (2) the process is intensified by periodic forcing; (3) a smaller output amplitude in comparison with the internal oscillation makes the system more robust to various disturbances.

Conclusions
The nonlinear dynamics discussed above shows that the design, operation, and control of a CSTR could be very complicated. The presence of the multiple steady states necessitates handy design of an efficient control system to regulate the operational conditions to desired states. The exploration on oscillatory CSTRs has evolved into two distinct directions: one is the elimination of the oscillations and the other is to take advantage of the process dynamics. Chaotic dynamics would also generate through coupling of two oscillatory CSTRs or CSTR forced. Therefore, a complete knowledge of static and dynamic behavior of these behaviors is required to understand, to operate, to control, and to optimize CSTRs.

Conflicts of Interest:
The authors declare no conflict of interest.
Likewise, x 1 x 2 3 and x 1 3 x 2 can vanish to make Equation (A5) simpler. This determines e and k, which are functions of a, b, c, d, and h, while h can be arbitrary. Since the rest of the deduction does not include e and k, their explicit expression will not be included here. From Equation (A3), dV/dt is sign definite when α(u) 0. However, when α(u) = 0, conditions for dV/dt being positive definite are determined by parameters of x 1 4 , x 2 4 , and x 1 2 x 2 2 , Observe that β 1 , β 2 , and β 3 are independent of g and j, while g and j are confidents need to be determined. Instead of determining these two confidents, we introduce σ (x 1 2 + x 2 2 ) 2 , which makes the sign of σ < 0 equal to Equation (A7) when Equation (A9) is satisfied, Solving Equation (A9) and the Lyapunov coefficient, σ, is obtained, Then, σ in Equation (A10) determines the stability of an internal oscillation.

Appendix B Bifurcaitons of the Limit Cycles
To study bifurcations of limit cycles, we introduce the Floquet theory. Considering Equation (3)  Floquet theorem provides solution structure for Equation (A11), which can be used to study bifurcations; when oscillatory q j at Equation (11) is included, one can verify that the solution v of Equation (A11) need not to be periodic, but should be the form as follows: where p(t) has period T 0 , and µ i is the ith-Floquet exponents, plus ρ i (T 0 ) = e µiT0 is the Floquet multiplier, satisfying