On the global nonlinear stochastic dynamical behavior of a class of exothermic CSTRs

Journal of Process Control 21 (2011) 1250–1264

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

On the global nonlinear stochastic dynamical behavior of a class of exothermic CSTRs Stefania Tronci, Massimiliano Grosso, Jesus Alvarez 1 , Roberto Baratti ∗ Dipartimento di Ingegneria Chimica e Materiali, Università di Cagliari, I-09123 Cagliari, Italy

a r t i c l e

i n f o

Article history: Received 16 December 2010 Received in revised form 26 July 2011 Accepted 26 July 2011

Keywords: Chemical reactor dynamics Stochastic estimation Chemical reactor stochastic dynamics Stochastic modeling Fokker–Planck equation

a b s t r a c t The problem of assessing the short and long time effects of stochastic fluctuations on the global-nonlinear dynamics of a class of closed-loop continuous exothermic reactors with temperature control and mono or bistable isothermal dynamics is addressed. The consideration of the problem within a Fokker–Planck (FP) stochastic framework yields: (i) the characterization of the global-nonlinear stochastic dynamics, and (ii) the connection between the deterministic and stochastic modeling approaches. The evolution of the state probability density function (PDF) is explained as the result of a complex interplay between deterministic dynamical features, initial PDF shape, and noise intensity. The correspondence between stationary PDF mono (or bi) modality and deterministic mono (or bi) stability is established, and the stochastic settling time is put in perspective with the deterministic, noise-diffusion, and escape times. The conditions for the occurrence of a retarded response, with respect to deterministic and noise-diffusion times, are identified. The proposed approach: (i) is illustrated with representative case class example, and (ii) constitutes an inductive step towards the development of a general-purpose stochastic modeling approach in chemical process systems engineering. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Industrial chemical process systems operate subjected to persistent fluctuations, which are due to variations around mean values of actuators, measurements, and high-frequency parasitic dynamics. Since the fluctuation frequencies are faster than the deterministic dynamics, they can be modeled as an exogenous stochastic noise input for the deterministic system. Depending on the system and its operating condition, the stochastic and deterministic modeling approaches can yield similar or considerably different results. For instance, in chemical reactive systems the presence of noise fluctuations can induce or eliminate a steady-state or a limit cycle [1,2] or may explain atypical long term behavior observed in experimental catalytic systems [3,4]. The global dynamical behavior of a stochastic nonlinear system is described by the evolution of the (possibly multimodal) state probability density function (PDF). The related stochastic modeling problem has been addressed with the simulation-oriented Monte Carlo (MC) and the analysis-oriented Fokker–Planck (FP) approaches.

∗ Corresponding author. Tel.: +39 0706755056; fax: +39 0706755067. E-mail address: [email protected] (R. Baratti). 1 On leave from Departamento de Ingenieria de Procesos e Hidraulica, Universidad Autonoma Metropolitana - Iztapalapa, 09340 Mexico, DF, Mexico. 0959-1524/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jprocont.2011.07.014

In the MC approach [5] a suitable sample of initial conditions and exogenous inputs (including time-invariant parameters) deviations are set, a bundle of discrete-time state motions is generated by running the nonlinear deterministic system over the data sample, and PDF identification-fitting techniques are applied to draw the evolution of the state PDF. The MC method is conceptually simple, can handle high-dimensional (up to ten-state) systems with a diversity of specialized numerical algorithms [5], but becomes cumbersome or intractable when the PDF is multimodal [6], and breaks down when the deterministic system evolves in close-tobifurcation condition [1,7]. The state of the art on the application of the MC method for chemical processes in general [5] and exothermic reactors [8] can be seen elsewhere, and here it suffices to mention that the MC method is one of the most widely applied stochastic modeling techniques in chemical process. In the FP approach the PDF dynamics are described by a multidimensional (with one dimension per state) linear partial differential equation (PDE), by a fundamental probability conservation principle with diffusive-convective transport mechanism [9,10]. The FP approach can a priori relate, in a qualitative way, important (mono/multimodality, settling time, and metastability) characteristics of the state (possibly multimodal) PDF dynamics with (mono/multistability, SS attractivity, dissipativity) features of the nonlinear deterministic dynamics. Due to the difficulty of solving the underlying FP PDE over an n-dimensional space (with one dimension per state), the quantitative evolution of the state PDF can be drawn for systems with up to four states [11–17]. The capabili-

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

ties and limitations of the FP and MC methods suggest that, as it is done in constructive control [18,19], the stochastic modeling task should be performed by adequately combining notions and tools of the two methods, in proportions which depend on the specific problem and modeling objective. From the FP approach theory and its applications to a diversity of problems in physics, electronics, chemistry, biology, and medicine [11–17], it is known that [9]: (i) the stationary PDF is unique and monomodal (or multimodal) if the deterministic system is monostable (or multistable), and (ii) metastability (extended permanence about a weak attractor) and escape time (duration of weak-tostrong attractor transition) are stochastic dynamic phenomena which do not exist in deterministic systems. The application of the FP approach to a stochastic system with on-line measurements yields the evolution of the conditioned state PDF and the corresponding FP equation with measurement injection is denominated nonlinear stochastic Kushner Filter (KF) [10]. When the PDF is approximated by a Gaussian (monomodal and symmetric) distribution, the KF becomes the nonlinear extended Kalman Filter (EKF), which is by far the most widely employed estimation technique in chemical process systems engineering. This kind of approximation has been employed to address the stochastic estimation problems for bioreactors described by gray models [20,21]. With the exception of a preliminary extension of the present work to two-state exothermic reactor with experimental data [22], the mechanistic model-based KF has not been applied to chemical process systems. The recent development of efficient on-line computational packages in computational fluid dynamics [23] enables the implementation of up to four-state nonlinear stochastic KF, and its application to an important class of estimation, control, fault detection and safe process design problems in chemical process, in the understanding that: (i) according to constructive control [24,25] and estimation [26] studies, only few-state subsystems of the process need to be estimated and/or controlled, and (ii) this observation applies to the development of the OF robust control version of the geometric control-based fault detection approach [8]. These considerations motivate the development of the FP-based nonlinear stochastic modeling approach in chemical process systems engineering. From previous studies in chemical process systems engineering it is known that in general when feedback is used in lumped [27] and distributed [28] deterministic systems the closed-loop dynamics is simpler (less nonlinear) than the open-loop dynamics, and the same is true for stochastic lumped [8,29,30] and distributed system. This simplified behavior feature offers the possibility of applying the stochastic FP-MC combined method to jointly perform the safe process and control designs over short and long time functioning horizons in the presence of fluctuations [29,30]. In chemical processes, the FP approach has been only applied to chemical reactors with one and two states. Excepting two studies on transient behavior [6,31], the studies have been circumscribed to the stationary behavior. Three studies [6,31,32] applied the nonlinear FP approach, while in the other cases the simulationbased local FP or MC methods have been performed. Regarding the stochastic stationary behavior of continuous exothermic reactors, the effect of the noise on the state variance in open [2,7] and closed [29] -loop regime has been assessed with numerical simulations, reporting that: (i) the difference between results with linear and nonlinear models grows with closeness to bifurcation condition [2,29], and (ii) the closed-loop concentration dynamics with PI control is almost decoupled from the fast nearly linear temperature dynamics [29]. Regarding reactor studies with nonlocal FP approach: (i) the stationary PDF of a one-state adiabatic combustion system has been obtained in analytic form, with a rigorous connection between mono (or bi) PDF modality and mono (or bi) deterministic stability [32], the numerically simulated stochastic

1251

steady-state transition (related to metastability) response of a bistable bioprocess has been detected by fitting the PDF evolution with a Two-Hump bimodal distribution [6], and (iii) the two (or six) dimensional FP equation resulting in a nonlinear sensitivity analysis of a two-state batch bioreactor (or two states with two uncertain initial values and two uncertain parameters) have been numerically solved using finite element method (two states) and a Gaussian-based particle method (two and six states) [31], demonstrating the feasibility of numerically solving the dynamic FP equations for reasonably complex systems. Summarizing, in the chemical process systems engineering field: (i) the FP stochastic modeling approach has been applied only to chemical reactors, (ii) with the exception of the adiabatic combustion reactor case [32], the stationary-stochastic and staticdeterministic behaviors have not been rigorously connected, (iii) in the context of a two-state bioreactor with nonmonotonic kinetics, the metastability and escape dynamic phenomena (which have important implications in safe process, control and fault detection designs) has been detected with simulations and bimodal timevarying PDF fitting techniques [6], but not rigorously connected with deterministic behavior features, (iv) excepting a preliminary extension of the present work to a two-state exothermic reactor with nonmonotonic kinetics and experimental data [22], the FP-based nonlinear stochastic KF has not been applied. In the systems theory field there are tractable numerical tools to address the FP equation-based stochastic modeling and estimation problems with up to four states, and the conceptual-oriented FP can assist the implementation of the MC method and the interpretation of its numerical results in higher dimensional systems. Thus, the combination of the FP and MC approaches offer the possibility of improving or developing methodologies for chemical process systems, including: (i) stochastic modeling of nonlinear systems, (ii) nonlinear stochastic Kushner filtering [10], (iii) nonlinear geometric control-based fault detection [8], and (iv) long term failure probability assessment in safe process design [33]. These considerations motivate the scopes of the present work: (i) the characterization of the global-nonlinear stochastic dynamical behavior of a class of temperature controlled exothermic reactors of practical interest [29], with emphasis on the formal connection between the stochastic and deterministic approaches, and (ii) the setting of an inductive step starting from a simple model towards the development of a general-purpose stochastic modeling approach for chemical processes, including a fruitful combination of the conceptual-oriented FP and numerical-oriented MC methods. First, the reactor concentration nonlinear deterministic dynamics is globally characterized, including a new approach in terms of potential-based dissipation recalled from mechanical and electrical system dynamics and control. Then, the stochastic nonlinear dynamics are described with the corresponding FP PDE equation. The stationary PDF is obtained in analytic form, relating PDF mono (or bi) modality with deterministic mono (or bi) stability. The transient PDF behavior is characterized in terms of dissipation mechanism, metastability and escape time phenomena, and the results are related to deterministic properties. The proposed approach is illustrated with representative reactor class example with deterministic monostability or bistability, depending on the kinetics parameters. 2. Stochastic modeling problem Consider a nonlinear chemical process (list of symbols in Table 2) ˜ x˙ = f (x, u, p) + f˜ [x, , u(t) + u(t), uu (t)],

x(0) = xo , t ∈ [0, ∞):=᏿ (1a)

1252

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

Table 1 Acronyms. OF SF SS QSS FP MC PDE PDF PP

Output feedback State feedback Steady state Quasi steady state Fokker–Planck Monte Carlo Partial differential equation Probability density function Phase portrait

˙ = f [x, , u(t), uu (t)],

(0) = o

(1b)

which evolves over the time interval ᏿, where x (or ) is the modeled (unmodeled) state, u(t) is the measured (exogenous and/or ˜ uu control) input with persistent actuator/measurement error u, is the unmeasured time-varying exogenous input with persistent fluctuations, p is the model parameter vector, and subsystem (1b) is the unmodeled parasitic (high frequency) dynamics associated with quasi steady-state (QSS) (list of acronyms in Table 1) assumptions (e.g., perfect mixing in a CSTR, static tray hydraulics and heat balance in a distillation column, fast species in a reaction network, etc.). Following the standard stochastic modeling approach in chemical process systems, let us consider that the initial state xo is a random variable with probability density function o (x), and that, from the perspective of the system dominant dynamics, the modeling error (f˜ ) manifests itself as an additive high-frequency disturbance which behaves as a normal (Gaussian) white noise w(t) with (possibly state dependent) intensity (variance) matrix Q, according to the following expressions xo ∼R[o (x)],

x ∈  = {x ∈ n |0 ≤ xi < ∞}

(2a)

(2b)

The enforcement of these assumptions upon system (1) yields the stochastic nonlinear dynamical model (x, t) : x ∈ ,

x˙ = f (x, u, p) + w(t), w(t)∼N[0, Q (x)],

x(to ) = xo ∼R[o (x)],

t ∈᏿

(3)

where (x, t) is the probability density function (PDF) at time t of the state stochastic process x(t). The case of given parameter uncertainty can be handled [10] by drawing a state-parameter dependent noise intensity q(x, p, t) with a first-order Taylor expansion of f about p [34]. The non-white (say colored) and/or non-Gaussian cases can Table 2 Nomenclature: nonlinear stochastic n-state system (Section 2). dij D f p qij Q t u uu w x   

᏿

(i, j) entry of the diffusion matrix D n × n diffusion matrix Model state function Model parameter vector (i, j) entry of the noise covariance matrix n × n noise covariance matrix Time Measured inputs Unmeasured inputs White noise Modeled state vector Unmodeled state vector Probability density function (PDF) State interval for n-dimensional system Time interval [0, ∞)

∂t (x, t) =

n n  

∂xi ∂xj [(x, t)dij (x)] −

i=1 j=1

n 

∂xi [(x, t)fi (x, u, p)],

i=1

0 < xi < ∞, 0 < t < ∞

(4a)

with initial (4b), conservation (4c), and boundary (4d) and (4e) conditions t=0:



(x, 0) = o (x) ∀x ∈ ;

(, t) d = 1 ∀t ∈ ᏿

(4b) (4c)



x=0: x = x∞ :

w(t) = f˜ [x(t), (t), u(t), uu (t)]∼N{0, Q [x(t)]}, t ∈ ᏿, Q (x) = [qij (x)], Q (x) = Q  (x)

be handled with a suitable dynamical extension driven by white noise or by adding white noises [35,36]. In the Monte Carlo (MC) method the PDF evolution is obtained as follows [5]: (i) first, a suitable sample of S(xo , w) of initial state-input pairs [xo , w(t)] is generated on the basis of the “data” (x, 0) = (xo ) and Q(x), (ii) then, a sample of state motions S(x) is generated by integrating the dynamical system (3) over the data sample S(xo , w), and (iii) finally, PDF identification-fitting techniques are applied to draw an analytic approximation of the PDF. The MC method functions well for high dimensional systems with deterministic monostability, but breaks down, even in one and two-state systems, with multistability or monostability in close-to-bifurcation condition [1,7,29]. This inherent limitation of the MC can in principle be remedied by resorting to the notions and tools from the FP method, by regarding the interplay between (multiplicity, stability and bifurcation) deterministic and (multimodality, response time and metastability) stochastic characteristics. In the Fokker–Planck approach, the evolution of the state PDF is given by the unique solution for the partial differential equation (PDE) [9]:

[d(0)][∂x (0, t)] = [f (0, p)](0, t) ∀t ∈ ᏿; ∂x (x∞ , t) = 0 ∀t ∈ ᏿

where dij (x) = qij (x)/2, and

 x2 x1

(4d) (4e)

(, t) d is the probability of having

the state x(t) in the n-dimensional hypercube delimited by x1 and x2 . The FP equation states the conservation of probability over the semi-infinite spatial domain : the rate of change of  (1st term in the left side of Eq. (4a)) is determined by diffusive plus advective flux, the diffusive flux (1st term in the right side of Eq. (4a)) is set by the noise intensity Q, and the advective flow (2nd term in the right side of Eq. (4a)) is set by the deterministic velocity drift f. According to the FP theory [9,10], when the deterministic version (x˙ = f ) of the stochastic system (3) has a unique SS with strong attractivity and robust stability, and the noise intensity matrix Q is independent of the state x, the multivariate PDF of the stochastic system is (nearly Gaussian) monomodal, and the PDF evolution can be obtained either: (i) by numerically solving (basically with computational fluid dynamics numerical techniques [23]) the FP PDE, when the system has up to four states, or (ii) by applying the MC method when the system has more states, with the assistance of the FP method in the design of the data sampling and PDF fitting stages in the light of the characteristics of the deterministic system nonlinear dynamics. Otherwise, when the deterministic behavior exhibits strongly nonlinear behavior such as multiplicity and structural instability, the stochastic system exhibits PDF multimodality and the possibility of the atypical retarded metastability associated with SS transition, which does not exist in deterministic modeling. In this case, as mentioned before, the MC breaks down even for low dimensional systems, and the way to address the stochastic modeling problem is by suitable combining ingredients from the FP and MC methods, according to the following rationale: (i) first, standard methods are applied to characterize the deterministic behavior,

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

(ii) then, FP theory is applied to qualitatively assess the stochastic behavior, and (iii) finally, quantitative PDF evolution results are obtained by either solving the FP PDF equation (when the system has few states), or by setting a MC method assisted by qualitative results drawn from the FP method (when the systems has a large number of states). These considerations motivate the study of the reactor class in the present work. The purpose is twofold: (i) the rigorous characterization of the stochastic global-nonlinear behavior of a rather simple reactor class which is of practical interest, including a formal deterministic-stochastic behavior connection as well as the identification of conditions for typical or atypical stochastic behavior, and (ii) the setting of a first step towards the development of the above discussed stochastic modeling approach for nonlinear chemical processes, through the combination of concepts and tools from nonlinear process dynamics, the FP approach, and the MC method. It must be pointed out that the reactor case example is a class of reactors over its kinetics parameter space, which includes: (i) as particular case the closed-loop exothermic reactor with PI temperature control and stochastic fluctuations of a previous study on stationary behavior [29], (ii) deterministic mono or bistability with structural stability or instability. Thus, according to the FP theory, the set of reactor stochastic responses should exhibit the basic possibilities that can be encountered in multi-state stochastic systems in general: stationary mono or multimodal PDF, and PDF transient with or without retardation by metastability.

Table 3 Nomenclature: closed-loop reactor (Section 2.1). Reactor area (m2 ) Input associated to observer Reactant concentration (mole m−3 ) Pure reactant concentration (mole m−3 ) Feed reactant concentration Neglected species concentration (mole m−3 ) Specific heat (kJ kg−1 K−1 ) Noise diffusion coefficient Damkholer number Kinetic constant (s−1 ) Adsorption constant (m3 mole−1 ) Dimensionless modeling error function Modeling error function (mole s−1 ) Heat transfer term Heat transfer term in dimensionless time units Agitator speed Parameters of the reaction rate R Volumetric flowrate (m3 s−1 ) Reaction rate (mole m−3 s−1 ) Dimensionless reaction rate Dimensionless time with respect to t Actual time (s) Reactor residence time (s) Reactor temperature (K) Ambient temperature (K) Jacket temperature (K) Feed temperature Temperature measurement Input vector Global heat transfer coefficient (kJ K−1 m−2 s−1 ) White noise ∼N[0, qT ] Reactor volume (m3 ) White noise ∼ N[0, qc ] Dimensionless concentration Dilution rate (s−1 ) Proportional gain of the PI control Proportional gain of the PI control Heat of reaction (kJ mole−1 ) Adiabatic temperature rise (K) Dimensionless concentration PDF Reaction rate in dimensionless time units Density of reacting mixture (kg m−3 ) Dimensionless adsorption constant Integral time of the PI control Dimensionless integral time of the PI control Unmodeled state Observer gain

A b C C0 Ce CN cp d Da kr ka f˜ F˜ H h N pR Qe R r t ta t T Ta Tc Te y u U

v

2.1. Reactor with PI temperature control In this subsection the continuous exothermic reactor employed in a previous simulation-based stochastic study on closed-loop stationary stochastic behavior with PI temperature control is recalled with two important differences: (i) the first order kinetics (with isothermal reactor monostability) is replaced by the more general inhibited Langmuir–Hinshelwood kinetics (with isothermal reactor mono or bistability), (ii) the system noise is related with parasitic dynamics, and (iii) the PI control is derived as a behavior approximation of an exact model-based robust nonlinear stabilizing state-feedback (SF) control. Consider a (possibly open-loop unstable) exothermic continuous reactor with volume V, where a reactant is fed with volumetric flow Qe at concentration Ce and temperature Te , and converted, at reactant concentration C and temperature T, into product via a reaction rate R according to the closed-loop reactor dynamics with temperature PI control (5d) (the nomenclature is listed in Table 3) o

C = −R(C, T, pR ) + (Ce − C) + F˜c (C, T, , u),

C(0) = Co

(5a)

o

T = R(C, T, pR ) + (Te − T ) − H(T − Tc ) + F˜T (C, T, , u), T (0) = To , y = T

o

 = F˜ (, C, T, u),

(5b)

(0) = o



Tc = T¯ c − a

(5c)





ta

(yT − T¯ ) + a −1

[yT () − T¯ ] d 0

(·) =

d(·) , dta

V w x  a H  

r 

a

 ω

ta is the actual time,  is the state of the unmodeled parasitic (high frequency) dynamics, including reactant concentration fluc˜ due to imperfect mixing and neglected species (C˜ N ),  tuations (C) is the heat of reaction (− H) divided by the volumetric specific heat capacity ( r cp ), H is the mixture-jacket heat transfer coefficient divided by the heat capacity (V r cp ), and a (or a ) is the proportional gain (or integral time) of the PI controller. The functions F˜c and F˜T denote the effect on the concentration-temperature dynamics of modeling errors as well as endogenous () and exogenous (u) fluctuations. The exogenous input (u) includes kinetics parameter errors (˜pR ), feed (T˜e ), jacket (T˜c ), and ambient (T˜a ) temperature, as ˜ fluctuations. well as volume (V˜ ), and agitator speed (N) In terms of dimensionless time t (with respect to the residence time t ) and concentration x (with respect to pure reactant concentration Co ) as well as reactor (T˜ ) and jacket (T˜c ) temperature deviations (from their steady-state values T¯ and T¯ c ), ·

where o

(5d)

 =

(·) =



(− H) , ( r cp )

˜ , u = (˜pR , Q˜ e , C˜ e , T˜e , T˜a , V˜ , N)

=

Qe , V

H=

˜ C˜ N , s )  = (C,

UA , (V r cp )

1253

d(·) , dt

t=

ta , t

t =

V , Qe

= a ,

(x, T˜ , pR ) = t R(C o x, T¯ + T˜ , pR ), x=

C , Co

T˜ = T − T¯ ,

 = T˜c ,

= T˜

=

a t

1254

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

˜ T˜˙ = −kT˜ − b,

the closed-loop reactor (5) is written as follows: x˙ = f (x, p) + f˜ (x, T˜ , , u), ˙ = f˜ (, x, T˜ , u),

x(0) = xo ;

(0) = o

T˜˙ = ˛(x, T˜ , pR ) + h + f˜T (x, T˜ , , u),







t

+ −1

 = −

T˜ (0) = T˜o ;

() d

T˜ (0) = T˜o

(11b)

(6a)

˜ , u, u), ˙ b˜˙ = −ωb˜ + f˜b (x, T˜ , b,

(6b)

˙ = f˜ (, x, T˜ , u),

(6c)

where

(6d)

˜ b(0) = bo ;

(0) = o

(11c) (11d)

˙ ˜ , u, u) ˜ ˙ = −ˇ(x, T˜ , T˜e , pR ) = −(1 + h − ∂T )(k ˜ T˜ + b) f˜b (x, T˜ , b, − (∂x )(f + f˜ ) − f˜˙ T − T˜˙ e

0

where f (x, p) = xe − x − r(x, pr ),

From previous constructive control studies for chemical processes [38], we know that, if the parasitic dynamics are monostable and the isothermal concentration dynamics:

r(x, pr ) = (x, 0, pR ),

(x, ˜ T˜ , pR ) = (x, T˜ , pR ) − (x, 0, pR )]

x˙ = f (x, p), ˛(x, T˜ , T˜e , pR ) =  (x, ˜ T˜ , pR ) − (1 + h)T˜ + T˜e ,

f˜ (x, T˜ , , u) = f˜c (x, T˜ , , u) − (x, ˜ T˜ , pR )

(f˜c , f˜T , f˜ )(x, T˜ , , u) = t (F˜c , F˜T , F˜ )(C o x, T¯ + T˜ , , u) According to the constructive chemical process control approach [24,37], the PI controller (6d): (i) can be realized in the observer-control form: ˙ = −ω − ω(h + ω ), =

−[(k + ω) h

+ ]

,

bˆ =  + ω ;

k + ω = h,

(7d)

kω =



(7e)

with estimation error dynamics: ˙ b˜ = −ωb˜ − b,

(8a)

˜ = b˜ o , b(0)

(8b)

b = ˛(x, T˜ , T˜e , pR ) + f˜T (x, T˜ , , u):=ˇ(x, T˜ , T˜e , pR )

(8c)

where bˆ is the (reduced observer-based) estimate of unknown input b of the closed-loop temperature dynamics (5d) written in the linear differential-nonlinear algebraic form: T˜˙ = h + b,

= T˜ + ε ,

T˜ (0) = T˜o ;

b = ˇ(x, T˜ , T˜e , pR )

and (ii) recovers the behavior (up to observer convergence) of the passivity-based robust nonlinear SF controller (9a), built the exact reactor model, which results from the enforcement of the prescribed closed-loop temperature dynamics (9b) upon the openloop temperature dynamics (5c), T˜c = T˜ −

[kT˜ + ˛(x, T˜ , pR ) + f˜T (x, T˜ , , u)] :=(x, T˜ , , u) h

T˜˙ = −kT˜

(9a) (9b)

in the understanding that the associated zero (i.e., isothermal) concentration-parasitic dynamics: x˙ = f (x, p),

x(0) = xo ;

˙ = f˜ (, x, 0, u),

(10a)

(0) = o

(10b)

represent the limiting behavior attainable with any robust nonlinear OF temperature controller. ˜ Accordingly, the application of the state coordinate (-to-b) change: b˜ = bˆ − b = ( + ωT˜ ) − ˇ(x, T˜ , T˜e , pR ) to system (5) yields the closed-loop reactor dynamics: x˙ = f (x, p) + f˜ (x, T˜ , , u),

x(0) = xo ;

(11a)

x(0) = xo ,

f (¯x, p) = 0

(12)

are mono stable with SS x¯ , with a suitable choice of the control gain pair ( , ) (according to a robust stability criterion) the closed-loop system (11) is monostable about the prescribed (possibly openloop unstable) SS associated with x¯ . If the isothermal dynamics are multistable and x¯ is one of the stable SSs, the closed-loop system is multistable and the nominal concentration x¯ is associated with one of the stable attractors, the controller functions only in the domain of attraction of the SS associated with x¯ , and this must be carefully accounted for in the control design. The preclusion of limit cycles in the closed-loop dynamics is ensured by [37,39]: (i) the absence of limit cycles in the zero-dynamics (10) of the nonlinear passive controller (9), and (ii) the behavior recovery property of the linear PI controller (6d) with respect to the nonlinear passive controller (9). In all cases, the concentration reaches a compact set about x¯ in about one settling residence time  ≈ 4t , and the temperature reaches (typically five to twenty times) faster its setpoint value (T˜ ≈ 0), without nothing important happening thereafter. Moreover, in the case of nonmonotonic kinetics (with respect to concentration), the deterministic isothermal dynamics (12) are the indistinguishable (unobservable) dynamics of the closed-loop reactor with temperature measurement [40]. The preceding observations evidence the central role that the isothermal concentration dynamics (12) play in the characteristics of the closed-loop reactor dynamics with conventional and advanced OF control. It must be pointed out that, as expected [25,27,28], the closed-loop deterministic dynamics (11) are a simpler (more linear) than their open-loop counterpart (Eqs. (6a)–(6c) with  = 0), and as we shall see: (i) in the case of monostable isothermal concentration deterministic dynamics, this will manifest itself as a simple stochastic dynamics with monomodal PDF, as reported in a previous four-state reactor study [8], and (ii) in the case of bistable isothermal concentration deterministic dynamics, the stochastic dynamics will be considerably more complex. 2.2. Stochastic modeling problem According to the preceding developments, the isothermal concentration-parasitic dynamics (10) is the limiting behavior attainable with robust OF temperature control, the isothermal concentration dynamics (12) determines the stability property of the closed-loop deterministic dynamics, and the closed-loop reactor motions reach a compact set of a SS in about one deterministic residence settling time (  ≈ 4t ), without nothing important happening thereafter. On the other hand, in stochastic modeling it is known that: (i) the stochastic dynamics depends strongly on the mono or multistability property of the deterministic dynamics, and (ii) the presence of stochastic fluctuations induces responses longer than the deterministic ones, and may induce a considerably retarded SS transition. These considerations suggest that fundamental aspects of the closed-loop reactor stochastic behavior can be

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

assessed by studying the behavior of the nearly isothermal dynamics (8a) in the presence of stochastic fluctuations (f˜ = w). For this aim, recall the closed-loop reactor dynamics, assume that the initial concentration xo is a random variable with (normalized) probability density function o (x), the isothermal value of f˜ is a white noise wc with intensity qc , and the disturbance f˜b of the fast observer error dynamics (11b) is a white noise v with intensity qv this is,

1255

rate R (in Eq. (5a)) with representation in terms of dimensionless concentration (x) and time (t) R(C, T, pR ) =

kr (T )C [1 + ka (T )C]

(17a)

2

(x, T˜ , pR ) = t R(C o x, T¯ + T˜ , pR )

(17b)

xo ∼R[o (x)],

(13a)

f˜ (x, 0, , u):=wc (t)∼N[0, qc ],

(13b)

where kr (or ka ) is the proportional (or equilibrium adsorption) constant kr (or ka ), and Co is the pure reactant concentration. Assuming pure reactant (xe = 1) is fed, the corresponding reactor isothermal dynamics are given by Eq. (10) with the rate of change function

f˜b (x, T˜ , , u):=v(t)∼N[0, qv ]

(13c)

f (x, p) = xe − x − r(x, pr )

(18a)

xe = 1

(18b)

and the resulting fast temperature-parasitic dynamics (11b) are in quasi stationary stochastic regime with respect to the slow concentration dynamics (11a). As shown in Appendix A, the enforcement of these assumptions upon system (11) yields the one-state stochastic nonlinear dynamical model: x˙ = f (x, p) + w(t), 0 ≤ x < ∞,

x(to ) = xo ∼R[o (x)],

w(t)∼N(0, q),

0=t<∞

(15a)

qT (k, ) =

(15b)

qv

T (x) > 0

(15c)

That the effect of the temperature stochastic fluctuations on the concentration dynamics of the nearly isothermal closed-loop reactor can be effectively lumped in a white noise term is a fact that has been demonstrated in a previous study with local nonlinear FP approach and numerical simulations [29]. The analysis of the effect of the PI control gain pair ( , ) on the concentration stochastic dynamics goes beyond the scope of the present study, and here it suffices to mention that the previous assessment [29] can be refined by applying the global-nonlinear approach of this study. Without restricting the approach, in this study the noise intensity q is assumed to be constant (independent of the state x). According to the stochastic model in the light of Eqs. (3) and (4), the evolution (x, t) of the concentration PDF is described by the FP equation ∂t (x, t) = ∂x [d∂x (x, t) − f (x, p)(x, t)], q 2

(16a)

(x, 0) = o (x)∀x ∈ [0, ∞)

(16b)

0 < x < ∞, 0 < t < ∞, d =

t=0:



∞

(, t)d = 1 ∀t ∈ [0, ∞)

(16c)

0

x=0: x=∞:

d(∂x )(0, t) = [f (0, p)](0, t) ∀t ∈ [0, ∞) (∂x )(∞, t) = 0 ∀t ∈ [0, ∞)

(16d) (16e)

For the purpose at hand let us consider a kinetics class which, depending on its parameters, yields mono or bistable isothermal reactor dynamics, in the understanding that for general n-state stochastic systems multistability is a necessary condition for the occurrence of the retarded stochastic SS transition phenomenon. For this purpose, let us consider a Langmuir–Hinshelwood reaction

Da(1 + ) x

(18c)

(1 + x)2

p = (Da, ) ∈ P where

q ≈ qc + T (x)qT (k, ) 2hk2

r(x, p) =

(14)

for the nearly isothermal closed-loop reactor, where (T is defined in Appendix A)



2

Da = t



kr (T¯ ) (1 + )2

(18d)

 ,

 = C o ka (T¯ )

(18e)

p ∈ P = {p ∈ 2 |0 ≤ Da ≤ 2, 0 ≤  ≤ 30}

(18)

Da is the Damkohler dimensionless number (i.e., the ratio of residence to reaction time) which acts as proportionality constant for the stabilizing first-order rate mechanism,  is the dimensionless adsorption constant associated to the destabilizing quadratic reaction inhibition mechanism, p is the Damkohler number-inhibition constant parameter pair, and P is the associated parameter set. With suitable modifications, the reaction rate function R includes Monod (isotonic) and Haldane (nonmonotonic) kinetics in bioreactor engineering [41]: the Haldane kinetics has a denominator in standard three-term quadratic monic form with two parameters, the LH kinetics corresponds to a particular choice of the two parameters of the Haldane kinetics, and the Monod kinetics has a linear inhibition term. Our reactor stochastic modeling problem consists in assessing the short and long time behavior of the concentration (possibly multimodal) PDF evolution (x, t), over the reactor parameter set P and in the light of the noise intensity and the (multiplicity and stability) characteristics of the deterministic behavior. We are interested in formally: (i) relating the PDF settling time with the deterministic and noise diffusion times, (ii) establishing the correspondence between deterministic and stochastic behavior, and (iii) identifying the conditions for the occurrence of PDF retarded response due to SS-to-SS transition phenomenon. 3. Deterministic behavior In this section, the deterministic behavior of the reactor (8) is characterized within a global dynamics framework (list of symbols in Table 4), using bifurcation [42,43], geometric dynamics [27,44,45], and dissipation [46] methods. Recall the reactor stochastic model (10), and write its deterministic dynamics x˙ = f (x, p),

x(to ) = xo ,

x ∈ X = [0, 1],

p∈P

(19)

Due to the presence of the quadratic inhibition in the denominator of the reaction rate (Eq. (17a)), this equation can have one or several SSs, depending on the parameter p, with bifurcation phenomena for some parameter values.

1256

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

B B+ (or B− ) e g ns P Pm (or Pb ) Pbe x¯ x¯ − (or x¯ + ) x¯ o x¯ s X X− (or X+ )

Bifurcation curve Low (or high) concentration branch of the bifurcation curve Lyapunov function Force related to the potential  Number of steady states Parameter space Parameter region of monostability (or bistability) Parameter curve with equipotential bistability Attractor for the monostability case Low (or high) concentration attractor for the bistability case Intermediate concentration repulsor for the bistability case Concentration value at bifurcation point Domain of attraction for the monostability case Domain of attraction of the attractor x¯ − (or x¯ + ) for the bistability case Potential energy Convergence rate

 

3.1. Multiplicity and stability

Fig. 1. Limit-point bifurcation diagram over the -Da parameter space (solid line) and equipotentiality locus for the two stable solutions (dotted line). Points A, B, C, and D correspond to steady state behaviors illustrated in Fig. 2a, 2b, 2c and 2d, respectively.

In the boundary of X = [0, 1] we have that x˙ = f (0, p) = xe = 1 > 0

(20a)

x˙ = f (1, p) = −r(1, p) < 0

(20b)

meaning that a motion x(t) started in the boundary of X moves towards the interior of X. This in turn implies that: (i) X is an invariant set in the sense that any trajectory born in X stays in X [27,45], and (ii) the ns SSs of the reactor are in the interior XI of X, or equivalently, the real roots (¯x1 , . . . , x¯ s ) of the algebraic equation f(x, p) = 0 are in XI , this is, f (¯xi , p) = 0,

i = 1, . . . , ns (p)

⇒

Moreover, the set of bifurcation pairs (x, p) in the threedimensional state-parameter space x − p satisfy the algebraic equation pair [43] f (x, p) = 0

(22a)

∂x f (x, p) = 0

(22b)

and the solution for (x, Da) of this equation pair yields the twobranch curve x± = ˛± ()

(23a)

Da± = ˇ± ()

(23b)

8 <  ≤ 30

(23c)

in the state-parameter x − p, where ±() −  , ˛ () = 4 =

2

,

() =



 2 − 8

The projection of this curve in the parameter plane P (18d) yields the v-shaped bifurcation curve (plotted in Fig. 1) B = B− ∪ B+ ⊂ P

(24a)

P = Pm ∪ B ∪ Pb

(24b)

B = {p ∈ P|ˇ () = Da,  ≥ 8} B−

B+ )

X − = [0, x¯ o ),

X + = (¯xo , 1]

and the PP on top of Fig. 2b and d. (iii) B (saddle-node bifurcation) where system (19) undergoes mono-bistability transition, with one attractor (¯xa with domain of attraction Xa ) and one saddle (¯xs with domain of attraction Xs ), concentration domain partition X = Xa ∪ Xs and the PP on top of Fig. 2c for p ∈ B. The PP for p ∈ B+ is the mirror image of the one in Fig. 2c. In branch B− (or B+ ) of the bifurcation curve, the saddle x¯ s results from the coalescence () of the repulsor x¯ o with the high (or low) concentration attractor x¯ + (or x¯ − ), and consequently for p ∈ B− :

xa = x¯ − ∈ X a = [0, x−s ),

X s = [¯xo+ , 1],

for p ∈ B+ :

x¯ s = x¯ o  x¯ +

xa = x¯ + ∈ X a = (¯xo− , 1], x¯ s = x¯ o  x¯ +

ˇ ()

(1 + )

±

X = X − ∪ x¯ o ∪ X + ,

X s = [0, x¯ o− ],

±

{−1 ± () + 5/2 −  2 [−1 ± ()]/8}

±

(i) Pm (monostability) where system (19) has one attractor x¯ in X, and the phase portrait (PP) on top of Fig. 2a. (ii) Pb (bistability) where system (19) has two attractors (¯x− with domain of attraction X− , and x¯ + with domain of attraction X+ ), one repulsor (¯xo ) in between, state interval partition

x¯ i = f −1 (p) ∈ XI = (0, 1) ∀p ∈ P (21)

±

σ

Table 4 Nomenclature: reactor deterministic behavior (Section 3).

(24c)

(or is the low (or high) concentration branch. Thus, where the bifurcation curve B induces the partition of the parameter space P in three sets (see Fig. 2):

At the bifurcation curve B (Fig. 1), system (19) is not structurally stable, in the sense that a small parameter change produces a different steady state structure (number of states and stability properties). System is robustly stable when its parameter p is sufficiently away from B. Given the one-dimensionality of the state space X in the light of the topological restrictions imposed by the SS structure and stability properties [27,44,45], the deterministic system (19) over P does not have limit cycles. This rules out the presence of limit cycles in the actual closed-loop dynamics (11) with PI control, and consequently, in their nearly isothermal stochastic approximation (14). From chemical reaction arguments [27], we know that a motion x(t) started at xo in the domain of attraction Xa of the attractor x¯ a

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

(a)

1257

X

X

+

X

a

X

s

X

-

X

+

g(x)

X

(d)

(c)

(b) -

-0

f (x)

x

0.0

0.2

0.4

0.6

0.8

1.0

0.2

x

0.4

0.6

0.8

1.0

0.2

0.4

x

0.6

x

0.8

1.0

0.2

0.4

0.6

0.8

1.0

x

Fig. 2. Phase portraits, restoring force g(x) and potential function (x) for the following cases: (a) monostability ( = 0 and Da = 2.0, linear model), (b) bistability with different potential minima ( = 20 and Da = 0.2350), (c) structural instability with one attractor and a not hyperbolic point ( = 20 and Da = 0.28), and (d) bistability with equipotential minima ( = 20 and Da ≈ 0.22903). Attractors, repulsors and not hyperbolic points are represented with black circles, white circles and gray squares, respectively.

reaches within a 98% deviation (referred to xo − x¯ a ) in about four residence times, this is

r ≈ 4tr ,

tr =

(V/Qe ) =1 t

(25)

without nothing important happening thereafter. 3.2. Dissipation For the deterministic-stochastic connection purpose at hand, the dissipation characterization will not be performed with the standard quadratic Lyapunov function approach, but instead a potential-based approach [46] recalled from mechanical systems will be employed. For this aim, recall the Newtonian concentration dynamics (19) and write it in Hamiltonian form: x˙ +   (x, p) = 0,



x(0) = xo

(26a)

x

(x, p) =

g(, p) d;

(26b)

0

g(x, p) = −f (x, p)

(26c)

where  is the “potential” (energy) referred to x = 0, and g = −f is the related “force”. For an attractor xa with domain of attraction Xa , the dissipation (rate of energy change) is given by x ∈ Xa :

e˙ = −f 2 (x, p) ≤ 0,

x = x¯ a :

e = (x, p) − (¯xa , p):=(x, ˜ p) ≥ 0;

e˙ = 0, e = 0

meaning that e is a Lyapunov function. The corresponding convergence rate  (in deterministic residence time units) and 98.7% deterministic settling time r are given by the expressions (x, p) =



−e˙ f 2 (x, p) > 0 ∀x ∈ X a = e (x, ˜ p)

r

{[x(), p]} d = ln 0

e(x )

o

e(x)

=4

(27a)

(27b)

The convergence rate  at the state x is equal to the squared force f2 at x divided by the energy content ˜ (with respect to the attractor x¯ a ) at x. In the monostability case (p ∈ Pm ) with SS x¯ and domain of attraction X (see PP on top of Fig. 2a): (i) the force g vanishes at the SS x¯ , and is restoring over X [i.e., g(x)(x − x¯ ) > 0 for x = / x¯ ], and (ii) consequently, the potential energy is positive ( ≥ 0) over X, and has a unique minimum, at x¯ , over X (see Fig. 2a). In the bistability case (p ∈ Pb ) with stable SSs x¯ − and x¯ + (with domains of attraction X− and X+ ), respectively, and repulsor x¯ o (see PPs on top of Fig. 2b and d): (i) g is a restoring force over X− (or X+ ), and vanishes at the SS x¯ − (or x¯ + ), and (ii) consequently, the potential  is positive ( ≥ 0) over X− (or X+ ), and has a minimum, at x¯ − (or x¯ + ), over X− (or X+ ). Globally speaking (i.e., over the entire domain X), the potential (x) has a double well shape (see Fig. 2b and d) with two local minima (corresponding to the attractors x¯ − and x¯ + ), and a local maximum (corresponding to the repulsor x¯ o ) between the two minima. When the two attractors have different potentials, the underlying deterministic system is said to have alopotential bistability, and the SS with the smallest (or largest) potential is referred to as strong (or weak) attractor. Otherwise, when the attractors have the same potential, the system is said to have equipotential bistability. Accordingly, parameter set Pb with bistability is partitioned as follows (see in Fig. 1) Pb = Pba ∪ Pbe ,

Pbe = {p ∈ P|[ + (p), p] = [ − (p), p]} ⊂ Pb

(28)

Pba

where is the two-dimensional set with alopotential bistability, Pbe is the curve with equipotential bistability, and [ − (p), + (p)] is the solution pair for x of the equation: (¯x− , p) = (¯x+ , p) which states that the two attractors (¯x− and x¯ + ) have the same potential. As it can be seen in Fig. 1, the equipotential curve Pbe separates the v-shaped bistability set Pb in two regions. As we shall see later, this dissipation characterization, which is not standard in chemical process systems engineering, will become a central point in the connection between the deterministic and stochastic modeling approaches.

1258

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

Table 5 Nomenclature: reactor stochastic behavior (Section 4). E F N td te tf tr t  

d

e

r

 ω

4.2. Stationary behavior The enforcement of the condition ∂t  = 0 upon the FP PDE (16a) yields the linear second-order ordinary differential equation (ODE) [(·) : = d(·)/dx]

Lyapunov functional for the PDF Linear probability flux Normalization constant Noise diffusion characteristic time Mean escape characteristic time Convection characteristic times Deterministic settling time Stochastic characteristic time Stochastic potential Convergence rate Diffusion settling time Escape settling time Deterministic settling time Stochastic settling time Weight function

  ¯ ¯ − f (x, p)(x)] = 0, [d(x)

¯ ¯ (0) = f (0, p)(0) d  (∞, t) = 0

¯ = Ne−(x,p) (x)

N=

In this section, the reactor stochastic nonlinear dynamics are globally characterized by assessing, with functional analysis and Lyapunov tools [47,48], the solution of Fokker–Planck’s (FP) PDE (16a) in the light of the features of the global deterministic dynamics. The deterministic and stochastic modeling approaches are rigorously connected, and the theoretical developments are verified with numerical testing (list of symbols in Table 5).

4.1. Analysis based on transport phenomena According to Eq. (16a), the linear PDE with concentrationdependent coefficient f(x, p) undergoes the combined effect of: (i) spread by diffusive transport with diffusion d, and (ii) drift by advective transport with “velocity” f(x, p) (see flows in the PPs in Fig. 2). The PDF reaches asymptotically a unique stationary value ¯ (x), regardless of the shape of the initial PDF o (x). In the monostability case with attractor x¯ in X, the velocity f over X/¯x points towards the stable SS x¯ , and vanishes at x¯ , implying that the PDF ¯ should asymptotically reach a monomodal PDF (x) with maximum value at x¯ . In the bistability case, the velocity f over X/¯x− (or X + /¯x+ ) points towards the stable SS x¯ − (or x¯ + ), and vanishes at x¯ − (or x¯ + ), implying that the PDF should reach asymptotically a bimodal PDF ¯ with two maxima (at x¯ − and x¯ + ), (x). The diffusion (td ) and convection (tf ) characteristic times (expressed in deterministic time units) are given by tf ≈ tr ,

tr = 1;

d = 4td ,

(30c)

¯ of this equation yields the stationary and the unique solution for  PDF

4. Stochastic behavior

1  tr , d

(30a) (30b)



(x, p) = −

td =

0


 ∞ 0

1 d



(31a) x

[f (, p)] d = 0

1 e−(,p)

d

(x, p) d

(31b)

(31c)

where N is a normalization constant, and  is the stochastic potential, which is proportional to the deterministic one (26b), with proportionality constant 1/d. ¯ According to Eq. (31): (i) the stationary PDF (x) is mono (or bi) modal if and only if the deterministic system is mono (or bi) stable, ¯ has two maxima at the attractors (¯x− and (ii) in the bimodal case,  + and x¯ ), and a minimum at the repulsor (¯xo ) between the attractors. In the case of equipotential bistability, the two maxima have the same value, or equivalently, the two attractors have comparable likelihood. In the case of alopotential bistability, the two maxima have the different values, and the stronger attractor exhibits more likelihood than the weaker one. In Fig. 3 are presented, for large (dh = 10−2 ) and small (ds = 10−3 ) diffusion values (29) and inhibition constant  = 20: (i) the first moment  of the stationary PDF as a function of the Damkohler number (Da), (ii) the concentration-Da deterministic solution diagram, and (iii) the stationary PDF at three Da values. As it can be seen in Fig. 3, at small diffusion value ds : (i) in the alopotential bistability case (Da = 0.246), the PDF is practically monomodal about the strong attractor (see the left and right panels in Fig. 3), (ii) in the

r = 4tr ,

where r (or d ) is the deterministic (or diffusion) 98.7% settling time. Assuming the fluctuation f˜ in the concentration dynamics (11a) has typical small (|f˜ |1 ≈ 0.04), intermediate (|f˜ |m ≈ 0.1), and large (|f˜ |h ≈ 0.14) standard deviation values, the corresponding diffusion values and times (in deterministic time units) are given by (ds , dm , dl ) = (10−3 , 5 × 10−3 , 10−2 ) (tds , tdm , tdl ) = (1000, 500, 100)

(29)

While in the deterministic model (without noise), the reactor concentration settles in about 4 residence time (i.e., r = 4), in the stochastic model (with noise) the PDF evolution has a settling time

 ( r , d ) which depends (in a way to be defined) on the deterministic ( r ) and diffusion ( d ) settling times as well as on the mono or bistability property of the deterministic system.

Fig. 3. (a) First moment  of the stationary solution of the FP model as a function of the Damkohler number ( = 20) at d = 10−2 (solid line) and d = 10−3 (dashed line). The solution diagram for the deterministic model is also reported (dotted line); (b) details of the stationary solution in the case of bimodality.

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

equipotential bistability case (Da ≈ 0.22903), the PDF is bimodal with two nearly nonoverlapping monomodal components (see the central panel in Fig. 3), and (iii) the PDF moment exhibits an abrupt change with Da, because, as d decreases, the PDF bimodality feature becomes less evident. 4.3. Dissipation Following functional analysis tools [47,48], the “energy” E of the PDF (x, t) is given by its weighted squared deviation with respect ¯ to the stationary PDF (x) (derivation in Appendix B):



∞

˜ 2 (x, t) dx ≥ 0 ω(x)

E=

(32a)

0

¯ ˜ t) = (x, t) − (x) (x, ω(x) = Ne

(32b)

(x,p)

(32c)

with weight function ω(x) determined by the stochastic potential ¯ , or equivalently by the stationary PDF . The application of Lyapunov’s second method yields the dissipation mechanism (derivation in Appendix B).



∞

˜ ω(x)F 2 [(x, t)] dx ≤ 0

E˙ = −d 0

f (x, p)

˜ x (x, t) − ˜ t)] =  F[(x,

d

(33a)

˜ t) (x,

(33b)

where F is a linear probability flux and E (32a) is a Lyapunov functional for the FP equation (16a). The corresponding convergence rate ω (in deterministic residence time units) and 98.7% energy settling time  are given by the expressions

 ∞

˜ t)] dx d 0 ω(x)F 2 [(x, ˙ ∞ ω (t) = −E/E >0 =







ω () = ln 0

0

E(0) E( p )

˜ 2 (x, t) dx ω(x)

(34a)

 =4

(34b)

The time-varying convergence rate ω depends in a complex ¯ (26a), way on: (i) the time-invariant shape of the stationary PDF  and (ii) the time-varying shape of the probability flux F (33b). As the deterministic force f vanishes, we have that f →0

⇒

(x, p) → 0,

ω(x) → 1/2,

˜ → ˜x F()

and as expected, the dissipation-convergence mechanism (34) becomes the familiar one



E˙ = −d

∞



0

1/2

2

˜ ∂x (x, t)

d −E˙ = = E

 ∞ 0

∞ 0

dx ≤ 0

2

˜ ∂x (x, t)

˜ 2 (x, t) dx 

(35a)

dx >0

(35b)

of the heat equation over the semi-infinite domain [0, ∞) [49,50], ˜ and d standing for temwith constant convergence rate 1/2 , and  perature deviation and thermal diffusion, respectively. According to the preceding dissipation mechanism ((32) and (33)), the time-varying PDF (x, t) reaches asymptotically the mono ¯ (or bi) modal stationary PDF (x) determined (through the deterministic potential ) by the deterministic mono (or bi) stability property, with convergence rate (34a) and settling time  (34b) depending in a complex way (34) on the deterministic (tr ) and diffusion (td ) characteristics times as well as on the shapes of the initial ¯ PDFs. [o (x)] and stationary [(x)]

1259

4.4. Transient behavior without metastability In the cases of deterministic monostability from any initial PDF o (x), and of bistability with initial PDF in the domain of attraction of the strong attractor, the transient PDF (x, t) reaches the station¯ ary PDF (x) with a settling time  that is bounded from below (or above) by the deterministic residence (or stochastic diffusion) settling time r (or d ), this is,



¯ (x, t)−→(x),

4 = r <  < d =

8 q

(36)

In the monostability case, the stationary PDF is monomodal with maximum height, at the attractor x¯ , determined by the deepness of the deterministic potential. In the bistability case, the stationary PDF is bimodal with two maxima, at the attractors x¯ − and x¯ + , with heights determined by the deepness of the wells of the deterministic potential. In both cases, the stochastic settling time  decreases with the attractiveness (potential deepness) of the SS, and the noise intensity q. In the case of deterministic bistability, there is an important difference between the transient behaviors of the deterministic (19) and stochastic (14) models: while the deterministic motion x(t) reaches either the strong or weak attractor, depending on the initial condition xo (see phase portrait at the top of Fig. 2b and d), the stochastic PDF motion (x, t) always reaches a bimodal PDF with its largest maximum at the strong attractor, meaning that, regardless of the ¯ o (x), the strong attractor is the more likely shape of the initial PDF  stationary steady-state. 4.5. Transient behavior with metastability According to the FP theory [9], while in the case of determinis¯ o (x) confined to the strong attractor tic bistability with initial PDF  domain, the transient PDF (x, t) reaches asymptotically a bimodal ¯ stationary PDF (x) with a settling time  < d (36) smaller than ¯ o (x) over the diffusion settling time d , in the case with initial PDF  the weak attractor domain the PDF reaches the same stationary PDF with a settling time  > d larger than the diffusion settling time

d . This retardation feature of the stochastic response is associated to the escape, metastability, and SS-to-SS transition phenomena [11,51]. From the perspective of the dissipation characterization performed in Section 4.3, this stochastic retardation phenomenon is explained by the complex dependency of the convergence rate on the initial and stationary PDF shapes, and on their interplay. In the stochastic metastability phenomenon the PDF undergoes, driven by noise, a weak-to-strong attractor domain transition with a stochastic settling time  = 4t  d related to the characteristic time te of the escape phenomenon, with te being the average time necessary for the concentration to move from the weak (say x¯ − ) to strong (¯x+ ) attractor by overcoming the repulsor minus weak attractor potential barrier  (see Fig. 2b and Eq. (37b)). This escape phenomenon underlies tunneling effects in quantum physics as well as chemical reaction by collision plus energetic barrier overcoming in molecular dynamics [52]. The exact formula to calculate the escape time te and its mean value can be seen in [52], from which it is possible to obtain the Arrhenius-type estimate te ≈ [ e (¯xo , x¯ + , d, p)]e (¯x

o ,¯x+ ,p)

:=fe (¯x+ , x¯ o , p),

x ∈ X − , e = 4te (37a)

(¯xo , x¯ + , p) =

[(¯xo , p) − (¯x+ , p)] >0 d

e (¯xo , x¯ + , d, p) =

2 1/2 [| (¯xo , p)|| (¯x+ , p)|] d

(37b) (37c)

1260

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

instance, in industrial practice a strongly exothermic polymerization reactor operation is said to be safe if its probability of exploding (moving towards a strong ignition attractor) in one year is less than one in a million [53].

Escape time t E

1e+9 1e+8

d = 10-3

1e+7

d = 5 10 -3 d = 10-2

4.6. Verification with PDF evolution tests

1e+6 1e+5 1e+4 1e+3 1e+2 1e+1 -0.032

-0.028

-0.024

-0.020

-0.016

Da - Da + Fig. 4. Escape time as estimated by expression (37) vs the distance between the actual Da parameter and the one corresponding to the limit point bifurcation (Da − Da+ ), at different d values.

on the basis of a parabolic approximation of the stochastic potential  = /d over [¯x− , x¯ + ]. Here,  and e are activation energy and collision frequency-like terms, respectively, and e is the escape settling time. Accordingly, when there is deterministic bistability and the initial PDF o (x) is in the domain of attraction of the weak attractor and the noise intensity has a typical low-noise value (say with td = 103 ), ¯ the transient PDF (x, t) reaches the stationary bimodal PDF (x) with a settling time  which is larger than diffusion settling time, according to the expressions



¯ (x, t)−→(x),

4 = r <  ≈ e > d

(38)

As the noise intensity decreases (i.e., d grows) and the strongweak potential barrier ( ) grows, the difference ( e > d ) between the escape and diffusion times grows by orders of magnitude. In Fig. 4 is presented the dependency of the escape characteristic time te on the Damkohler number difference (Da − Da+ ) (with Da+ being the bifurcation value in the B+ curve of Fig. 1) according to formula (32), for (xe , ) = (1, 20) and the three diffusion values (ds , dm , and dl ) listed in Eq. (29), showing that: (i) as expected, the escape time te decreases with the approach of Da to its bifurcation value Da+ , (ii) at large noise intensity (dl ), te grows moderately with Da, with values from te ≈ 4 deterministic characteristic times to ≈102 diffusion characteristic times, and (iii) at small noise intensity (ds ), te grows significantly with Da, with values from 102 to 107 diffusion characteristic times td . In the case of alopotential deterministic bistability with initial monomodal distribution in the domain of attraction of the weak attractor and sufficiently low noise intensity the escape time becomes similar or larger than the diffusion time (i.e., te ≤ td ): (i) the initial PDF transient response appears to remain for a “long time” (with respect to r ) nearly monomodal over the weak attractor domain, and this behavior is called metastability [11,51], (ii) then, at some point, the PDF starts moving “slowly” towards its stationary nearly monomodal PDF over the strong attractor domain, and (iii) at the “long” settling time  ≈ e (37a) associated with the escape phenomenon, the PDF reaches, its stationary bimodal shape with largest likelihood about the strong attractor. In stochastic system dynamics, this weak-to-strong domain PDF shape evolution is known as (mean) SS-to-SS transition. The characterization of this weak-to-strong attractor long-term transient due to the accumulation of energy injected by the exogenous fluctuations (depending on equipment and control system) has applicability in safe reactor control and process designs. For

On the basis of the deterministic (Section 3) and stochastic (Sections 4.1–4.5) behavior characterizations performed over the set P (18d) of Damkohler-adsorption parameter pairs (Da, ), a set of three test runs were designed to illustrate the manifestation of the deterministic behavior characteristics on the main short and long term stochastic behavior features. For this aim, the adsorption constant ( = 20) and the noise intensity (d = 10−3 : typical noise level) were fixed, and three Damkohler numbers were set in order to have three case examples with prescribed deterministic stability and initial PDF shape (characteristics listed in Table 6): • Case 1 (Da = 0.160) with deterministic monostability, and without possibility of stochastic metastability for any initial PDF. • Case 2 (Da ≈ 0.22903) with deterministic equipotential bistability, possibility of metastability, and with initial PDF (about the repulsor) so that metastability is not at play. • Case 3 (Da = 0.260) with deterministic alopotential bistability, possibility of stochastic metastability, and with initial PDF (about the weak attractor) so that metastability is at play. In all cases the deterministic settling time is r = 4. In Table 6 are listed: (i) the deterministic SSs and their stability properties, (ii) the estimates, according to Eq. (37), of the escape settling times for the cases of equipotential ( e = 1014 ) and alopotential ( e = 3.2 × 104 ) bistability, and (iii) the mono/bimodal shape feature of the initial PDF o (x). The difference of 10 orders of magnitude is due to the fact that in the equipotential case the difference between the maximum and the minimum value of the potential is rather large (see potential vs x plot in Fig. 2b–d and the te vs Da plot in Fig. 4). In Case 2 with equipotential bistability and rather large escape time, the initial monomodal PDF o (x) will be set with a per-domain area ¯ similar to the one of the stationary PDF (x), in such a way that the metastability phenomenon is not activated. Differently, in Case 3 with alopotential bistability and smaller escape time, the initial PDF will be confined in the weak domain of attraction, in such a way that the metastability phenomenon is activated. According to the theoretical developments presented in previous subsections: (i) in Cases 1 and 2 without metastability, the stochastic settling time should be bounded from below (or above) by the deterministic (or stochastic diffusion) settling time r (or d ) (36), according to Eq. (39a), and (ii) in Case 3 with metastability, the stochastic dynamics should exhibit a retarded response set by the associated escape settling time e (37), according to Eq. (39b), this is, Cases 1 and 2 ( e ≈ 1014 ) : 3

Case 3 ( d = 4 × 10 ) :

4 = tr <  < d = 4 × 103 4

4 = r <  ≈ e ≈ 3.2 × 10 > d

(39a) (39b)

It must be pointed out that the preceding qualitative and qualitative behavior assessments for the three case examples have been drawn a priori (i.e., before simulation-based analysis of FP’s equation) from theoretical arguments, on the basis of: (i) the qualitative and quantitative characterization of the deterministic behavior, (ii) the assessment of the stochastic behavior in the light of the diffusion mechanism and the interplay between the diffusion and drift mechanisms associated with FP’s equation, and (iii) the escape time estimates according to Eq. (37). Thus, the simulation-based stage of our study has two purposes: (i) the illustration and verification of the preceding theoretical

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

1261

Table 6 Data and results for the PDF transient tests (Section 4.6). Inhibition constant  = 20 Noise intensity q = 2d, d = 10−3

Damkholer Point (, Da) in the bifurcation diagram (Fig. 1) Stability property Attractor x¯ − Repulsor x¯ o Attractor x¯ + Escape settling time e Initial PDF Stationary PDF Stochastic settling time 

Deterministic settling time r = 4 Diffusion settling time d = 4 × 103 Case 1

Case 2

Case 3

0.160 A Monostability

0.22903 D Equipotential bistability 0.018 0.206 0.676 4 × 1014 Monomodal about x¯ + Bimodal equimaxima at: x¯ − , x¯ + 20

0.260 B Alopotential bistability 0.014 (strong) 0.303 0.583 (weak) 3.2 × 104 Monomodal about x¯ o Bimodal large maximum at: x¯ .small maximum at: x¯ + 3.2 × 104

0.81 Bimodal Monomodal maximum at: x¯ 20

developments and results, and (ii) the refinement of the transient behavior characterization in terms of PDF shape, stochastic settling time, and onset of the metastability phenomenon. In Fig. 5 is presented the PDF response (x, t) for Case 1 (with deterministic monostability and initial bimodal PDF), showing that, as expected from the theoretical developments the PDF: (i) becomes monomodal with settling time close to the deterministic one ( r = 4); (ii) evolves with monomodal shape without ¯ metastability, reaching the stationary function (x) in a settling time  ≈ 20 (around five times the deterministic settling time

r = 4) between the deterministic and noise diffusion settling times (4 = r <  ≈ 20 < d = 4 × 103 ) as predicted by expressions (36) and (39a). In Fig. 6 is presented the PDF response (x, t) for Case 2 (with deterministic equipotential bistability and initial monomodal PDF about the deterministic repulsor x¯ o ) showing that, as in Case 1, the PDF: (i) becomes bimodal with settling time close to deterministic one ( r = 4), and (ii) reaches, without metastability, the stationary ¯ shape (x) with a settling time  ≈ 20 (around five times the deterministic settling time r = 4) between the deterministic and noise diffusion settling times (4 = r <  ≈ 20 < d = 4 × 103 ) as predicted by equations (36) and (39a). In Fig. 7 is presented the PDF response (x, t) for Case 3 (with deterministic alopotential bistability and initial monomodal PDF confined in the domain of attraction X+ of the high-concentration weak attractor x¯ + ), showing that, in agreement with the theoretical ¯ developments, the PDF: (i) reaches its stationary PDF function (x) with a settling time  ≈ e quite similar to the a priori estimate (37) of the escape settling time e ≈ 3.2 × 104 , in agreement with expres-

Fig. 5. Time evolution, in logarithmic scale, of (x,t) for  = 20, Da = 0.160 and d = 10−3 . The initial condition is a bimodal PDF. Continuous line: distribution at r ; dashed line: distribution at  .

Fig. 6. Time evolution, in logarithmic scale, of (x,t) for  = 20, Da ∼ 0.229 (equipotentiality) and d = 10−3 . The initial condition is a Gaussian PDF. Continuous line: distribution at r ; dashed line: distribution at  .

sions (36) and (39b), with largest stationary likelihood about the strong attractor x¯ − and appreciable span over its domain of attraction X− , and (ii) exhibits slow SS-to-SS transition likelihood with a metastability period from 1m ≈ 20 (i.e., the stochastic settling time

 = 20 when metastability is not at play) to 1m ≈ ≈ 3.2 × 104 (i.e., the escape settling time predicted by Eq. (37)). In Cases 1 (Fig. 5) and 2 (Fig. 6) without metastability, the initial bimodal (Case 1) or monomodal (Case 2) PDF becomes monomodal

Fig. 7. Time evolution, in logarithmic scale, of the (x,t) for  = 20, Da = 0.260 and d = 10−3 . The initial condition is a Gaussian PDF. Continuous line: distribution at  , dashed line: distribution at e .

1262

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

(Case 1) or bimodal (Case 2) in one deterministic settling time ( r = 4), and then the PDF undergoes a gradual change until the stationary monomodal (Case 1) or bimodal (Case 2) PDF is reached at time  = 20 (between the deterministic and diffusion settling times, and closer to the deterministic one). In Case 3 with metastability at play the response is longer and more complex: (i) first, the initial monomodal PDF over the weak attractor domain X+ (about the weak attractor x¯ + ) remains monomodal over X+ , with small change, along the settling time period [0,  = 20] of Cases 1 and 2 without metastability at play, (ii) then, over the diffusion-to-escape long settling time interval [  = 20,  ≈ e ≈ 3.2 × 104 ] the PDF becomes bimodal with growing (or decreasing) maximum about the strong (or weak) attractor x¯ − (or x¯ + ), and (iii) the PDF reaches its stationary PDF shape at  ≈ e ≈ 3.2 × 104 . For this (typical) noise intensity value (q = 2 × 10−3 in deterministic characteristic time units), the response of the PDF with metastability is one thousand times slower than the PDF response without metastability. 4.7. Discussion The stochastic closed-loop reactor dynamics has a PDF response determined by three characteristic times: the deterministic (residence) time tr , the stochastic diffusion time td fixed by the noise intensity q, and the stochastic escape time te determined by the strong-minus-weak potential difference. While the deterministic and stochastic-diffusion times are always present, the escape time is at play only when there is deterministic bistability with initial PDF about the weak attractor. Typically, the deterministic settling time is about four residence times. Depending on the noise intensity: (i) the stochastic diffusion time ranges from td = 102 (ample fluctuations |f˜ | ≈ 0.14) to 103 (narrow fluctuations |f˜ | ≈ 0.04), and (ii) the stochastic escape time, close to the bifurcation point (right side in Fig. 4), ranges from 20 (ample fluctuations) to 5 × 103 (small fluctuations) residence times, and from 50 (ample fluctuations) to 4 × 108 (small fluctuations), away from the bifurcation (left side in Fig. 4). For large fluctuations the three time responses may overlap, for small fluctuations they may be separated by orders of magnitude. These considerations must be accounted for, depending on the specific reactor and stochastic modeling objective: while for standard online observer and control emphasis is placed on the transient along the deterministic-to-diffusion (hours to days) horizon, for off or on-line quality, productivity, and safety assessments with respect to equipment, control and monitoring schemes the longterm diffusion-to-escape time (weeks to months or years) horizon must be regarded. The FP-based global-nonlinear stochastic characterization of the nearly isothermal closed-loop reactor are in agreement with the fact that the nonisothermal reactor with isothermal mono (or bi) stability is globally observable (or detectable) [40], in the sense that the measured temperature determines one (or two) concentration trajectories, and in the case of detectability the two concentration converge with dilution rate speed. Thus, in the presence of noise: (i) the nearly isothermal reactor exhibits a stationary mono (or bi) modal concentration PDF when there is deterministic mono (or bi) stability, and (ii) the stochastic stationary concentration PDF is attained with a settling time determined by the interplay of deterministic, diffusion and escape times. While deterministic detectability establishes the existence of two concentration values for the same temperature, in the stochastic approach the maxima of the bimodal PDF correspond to the most probable concentration values. These considerations suggest that the deterministic estimator and output feedback control design of nonisothermal reactor with nonmnonotonic kinetics [40] can be improved with the employment of the nonlinear stochastic filter [22] developed as an extension of the present study.

In spite of its simplicity, the reactor class case example of the present study captures the basic feature that underlie the stochastic behavior of a wider class of more complex systems, and consequently, the general-purpose conceptual framework associated with the present study can explain results and assist the test design for stochastic studies with MC method. For instance, the numerically drawn statements that in chemical reactive systems the presence of fluctuations can introduce considerable changes in dynamical behavior [1,2] and that differences between stochastic modeling results with linear and nonlinear approaches grows with closeness to bifurcation condition [2,29]. In the two-state bistable bioreactor study [6] the PDF bimodality identification and features can be performed in a more systematic way with the twostate extension of the present study [22]. Moreover, in a recent study on the nonlinear control-based fault detection problem for a four-state exothermic polymer reactor [8], from the application of model-based MC simulation with T2 statistics-based PDF fitting it was found that the open-loop stationary PDF was monomodal and asymmetric (nongaussian), and the closed-loop stationary PDF was more narrow and nearly Gaussian. According to the FP theory, these results can be explained as follows: (i) the open-loop PDF is monomodal and asymmetric because the reactor is openloop stable with an asymmetric deterministic potential, and (ii) the exact model-based nonlinear geometric SF controller induces a partially linear-decoupled closed-loop dynamics with a unique strong attractor and a nearly parabolic (symmetric) deterministic potential. This statement exemplifies the way in which the MC and FP methods complement each other on this particular polymer reactor case: following the FP approach, standard tools employed in chemical reaction engineering are applied to conclude the monostability and attractivity properties of the open and closed-loop deterministic dynamics, and these results establish a priori the feasibility of performing a monomodal PDF fitting through MC simulation.

5. Conclusions The global-nonlinear dynamical behavior of the closed-loop reactor with Langmuir–Hinshelwood (inhibited) kinetics has been characterized in terms of deterministic dynamical features and noise level intensity, including stationarity mono/bi modality, transient response (shape and speed) features in terms of deterministic (residence) and stochastic (diffusion and escape) characteristic times, initial PDF shape, and presence of long-term metastability phenomenon (not present in deterministic systems) due to deterministic bistability and initial PDF about the weak attractor. Two key connecting points between the deterministic and stochastic approaches are the deterministic global-nonlinear dynamical characterization in terms of a potential-based dissipation, and the correspondence of this potential with the one of the stationary PDF. The theoretical findings were illustrated and verified with numerical run tests. This theoretical-drawn global-nonlinear characterization: (i) rigorously connects the global-nonlinear deterministic and stochastic modeling approaches, (ii) explains or agrees with the results of previous open and closed-loop reactor studies with local-nonlinear MC and FP approaches, and (iii) enables the connection between the deterministic and stochastic nonlinear state estimation approaches. The selected reactor class case study, in spite of its simplicity, captures the stochastic behavior of a wider class of more complex systems, meaning that the conceptual framework associated with the present study: (i) can employed to assist the sample and PDF fitting designs of the MC method, and (ii) constitutes an inductive step towards the development of a general-purpose stochastic modeling approach in chemical process systems engineering, by combining the FP and MC approaches. The combined FP-MC method

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

approach offers the possibilities of exploiting the existing deterministic behavior characterization methods in the assessment of the short and long term stochastic behavior of open and closedloop dynamics in the presence of fluctuations, and of developing more systematic means to address the joint safe process-control design problem with emphasis on the interplay between control gains and safety measures. The present study has conceptual and technical points of departure for the consideration of several problems in chemical reactors in particular and processes in general, among them are: (i) nonlinear stochastic filtering for nonisothermal reactors [22] and bioreactors [6,20,21], (ii) connection through the notion of potential with power shaping nonlinear control [54], (iii) medium and long term stochastic behavior assessment for off or on-line quality, productivity, and safety assessments with respect to equipment, control, monitoring and maintenance schemes, and (iv) improvement of constructive estimation [26] and control [24,25] schemes where the multi-state MIMO problem is reduced to a set of SISO problems for one to three-state subsystems, with the possibility of behavior recovery of nonlinear robust OF control via linear PI control with antiwindup protection, and (v) the application of the constructive approach for the development of the robust OF control version of the nonlinear SF control-based fault detection approach [8].

Appendix B. PDF dissipation Here, the dissipation mechanism (33) of FP’s equation (16) is obtained by applying functional analysis [47,48] and Lyapunov [49] tools. Recall the stationary solution (31) of the FP equation (16), and set the coordinate change

Regard closed-loop reactor dynamics (11), take the time derivative of the temperature dynamics (11b) and substitute the estimation error dynamics (11c), solve for b˜ the temperature dynamics (11b) to get (A1b), and recall (A1c) the stochastic modeling assumption (13c) T˜¨ + kT˜˙ = −b˜˙ = ωb˜ − f˜b

(A1a)

b˜ = −kT˜ − T˜˙

(A1b)

b˙ = −f˜b = v(t)∼N[0, qv ]

(A1c)

The substitution of (A1b) and (A1c) into (A1a) yields the temperature-estimation error dynamics in Langevin’s equation form (A2a) with stationary covariance qT for the temperature fluctuations [36]: T˜¨ + (k + ω)T˜˙ + (kω)T˜ = v

(A2a)

qv 2(k + ω)(kω)

(A2b)

By virtue of Assumption (13b), the disturbance f˜ of the closedloop concentration dynamics (11a) can be approximated, through a first-order Taylor series expansion about T˜ = 0, by (A3a) and the covariance q of f˜ by (A3b) and (A3c): f˜ (x, T˜ , , u) ≈ wc + [∂T f˜ (x, 0, , u)]T˜

(A3a)

q ≈ qc + T (x)qT

(A3b) 2

¯ u)] ¯ >0 T (x) = [∂T f˜ (x, 0, , According to Eqs. (7d)–(7e) (k+

(A3c)

ω, kω) = (h , / ), and the substitution of this equation into (A3b) and (A3c) yields Eqs. (15a)–(15c)

¯ ˜ t) = (x, t) − (x) (x,

(B1)

to take FP equation (16a) into the selfadjoint form ¯ t = (wx )x :=Lw (), 



¯ w(x) = (x),

¯ 0 x0 = 0, d

∞

¯ ∞ x∞ = 0, 

¯ t) dx = 0 (x)(x, 0

where Lw is a Sturm–Liouville operator with weight w. Multiply the domain equation by , and substitute t = (2 /2)t to obtain



2 (x, t) ¯ (x) 2



= (x, t)[w(x)x ]x t

Integrate over [0, ∞)



∞

¯ (x) 0

2 (x, t) 2



 dx =

t

∞

(x, t)[w(x)x ]x dx:=R 0

The exchange of x-integration and t-differentiation operations in L yields (B2b), and the integration by parts of R followed by the application of the boundary conditions yields (B2c), and the substitution of (B2b) and (B2c) into (B2a) yields the dissipation (B2d), L = E˙

Appendix A. Stochastic model

˜ t) (x, , ¯ (x)

(x, t) =

L:=

The authors kindly acknowledge Regione Sardegna for the financial support (CRP2 370) and J. Alvarez also for the support through the program “Visiting Professor 2008”, for the realization of this work at the Dipartimento di Ingegneria Chimica e Materiali of the University of Cagliari.

qT =

for the noise intensity q of the nearly isothermal closed-loop concentration stochastic dynamics (14).



Acknowledgements

1263

E=

1 2

(B2a)

 

∞

(B2b)

w(x)2x (x, t) dx

(B2c)

w(x)2x (x, t) dx

(B2d)

∞

R=−



2 ¯ (x, t)] dx (x)[ 0

0

E˙ = −

∞

0

where E is the “energy” associated to the PDF relative deviation . Finally, the application of the coordinate change (B1) to Eq. (B2d) ˜ yields the dissipation (33) in PDF deviation coordinate (x, t). References [1] L.K. Doraiswamy, B.D. Kulkarni, Relevance of stochastic modeling in chemically reacting systems, Ind. Eng. Chem. Fundam. 25 (1986) 511–517. [2] N.J. Rao, D. Ramkrishna, J.D. Borwanker, Nonlinear stochastic simulations of stirred tank reactors, Chem. Eng. Sci. 29 (1974) 1193–1204. [3] W.R.C. Graham, D.T. Lynch, CO oxidation on Pt: model discrimination using experimental bifurcation behavior, AIChE J. 33 (1987) 792–800. [4] M.P. Harold, P. Luss, An experimental study of steady-state multiplicity features of two parallel catalytic reactions, Chem. Eng. Sci. 40 (1985) 39–52. [5] A. Doucet, N. de Freitas, N. Gordon, Sequential Monte Carlo Methods in Practice, Springer, New York, 2001. [6] T. Vesterinen, R. Ritala, Bioprocesses and other production processes with multi-stability for method testing and analysis, in: European Symposium on Computer Aided Process Engineering (ESCAPE 15) 2005, Barcelona, Spain, 2005. [7] T.M. Pell, R. Aris, Some problems in chemical reactor analysis with stochastic features, Ind. Eng. Chem. Fundam. 8 (1969) 339–345. ˜ ˜ J.F. Davis, P.D. Christofides, Enhancing datade la Pena, [8] B.J. Ohran, D. Munoz based fault isolation through nonlinear control, AIChE J. 54 (2008) 223–241. [9] H. Risken, The Fokker–Planck Equation: Methods of Solutions and Applications, Springer-Verlag, Berlin, 1996. [10] A.H. Jazwinski, Stochastic Processes and Filtering Theory, Academic Press Inc., New York, 1970. [11] P. Hanggi, P. Jung, Bistability in active circuits: application of a novel Fokker–Planck approach, IBM J. Res. Dev. 32 (1988) 119–126.

1264

S. Tronci et al. / Journal of Process Control 21 (2011) 1250–1264

[12] S. Chandrasekhar, Stochastic problems in physics and astronomy, Rev. Mod. Phys. 15 (1943) 1–89. [13] T.K. Soboleva, A.B. Pleasants, Population growth as a nonlinear stochastic process, Math. Comput. Modell. 38 (2003) 1437–1442. [14] D.C. Mei, C.W. Xie, L. Zhang, The stationary properties and the state transition of the tumor cell growth model, Eur. Phys. J. B 41 (2004) 107–112. [15] C.F. Lo, Stochastic Gompertz model of tumor cell growth, J. Theor. Biol. 248 (2007) 317–321. [16] J. Kuipers, P. Hoyng, J. Wicht, G.T. Barkema, Analysis of the variability of the axial dipole moment of a numerical geodynamo model, Phys. Earth Planet. In. 173 (2009) 228–232. [17] D.W. Huang, H.L. Wang, J.F. Feng, Z.W. Zhu, Modelling algal densities in harmful algal blooms (HAB) with a stochastic dynamics, Appl. Math. Modell. 32 (2008) 1318–1326. [18] R. Sepulchre, M. Jancovic, P.V. Kokotovic, Constructive Nonlinear Control, Springer-Verlag, New York, 1977. [19] M. Krstic, I. Kanellakopoulos, P.V. Kokotovic, Nonlinear and Adaptive Control Design, Wiley, New York, 1995. [20] N.R. Kristensen, H. Madsen, S.B. Jorgensen, A method for systematic improvement of stochastic black-box models, Comp. Chem. Eng. 28 (2004) 1431–1449. [21] N.R. Kristensen, H. Madsen, S.B. Jorgensen, Parameter estimation in stochastic grey-box models, Automatica 40 (2004) 225–237. [22] S. Tronci, A. Balzano, J. Alvarez, R. Baratti, Global-nonlinear stochastic estimation of exothermic reactors with temperature measurement, in: Proceedings of the 9th International Symposium on Dynamics and Control of Process Systems 2010 (DYCOPS 2010), July 5–7, Leuven, Belgium, 2010. [23] A. Balzano, S. Tronci, R. Baratti, Accurate and efficient solution of distributed dynamical system models, in: 20th European Symposium on Computer Aided Process Engineering – Escape 20, Ischia, NA, 6–8 June 2010, Comp.-Aided Chem. Eng. 28 (2010) 421–426. [24] E. Castellanos-Sahagun, J. Alvarez, J. Alvarez-Ramírez, Two-point compositiontemperature control of binary distillation columns, Ind. Eng. Chem. Res. 45 (2006) 9010–9023. [25] J. Alvarez, P. Gonzalez, Constructive control of continuous polymer reactors, J. Process Control 17 (2007) 463–476. [26] J. Alvarez, C. Fernandez, Geometric estimation of nonlinear process systems, J. Process Control 19 (2009) 247–260. [27] J. Alvarez, J. Alvarez, R. Suàrez, Nonlinear bounded control for a class of continuous agitated tank reactors, Chem. Eng. Sci. 46 (1991) 3235–3249. [28] I. Karafyllis, P.D. Christofides, P. Daoutidis, Dynamic of a reaction-diffusion system with Brusselator kinetics under feedback control, Phys. Rev. E 59 (1999) 372–380. [29] M. Ratto, A theoretical approach to the analysis of PI-controlled CSTRs with noise, Comput. Chem. Eng. 22 (1998) 1581–1593. [30] O. Paladino, M. Ratto, Robust stability and sensitivity of real controlled CSTRs, Chem. Eng. Sci. 55 (2000) 321–330. [31] I. Horenko, S. Lorenz, C. Schutte, W. Huisinga, Adaptive approach for nonlinear sensitivity analysis of reaction kinetics, J. Comput. Chem. 26 (2005) 941–948. [32] M. Oberlack, R. Arlitt, N. Peters, On stochastic Damkohler number variations in a homogeneous flow reactor, Combust. Theor. Model. 4 (2000) 495–509.

[33] J.E. Berryman, D.M. Himmelblau, Stochastic modelling of reactors for process analysis and design – I. Isothermal steady state well mixed and plug flow reactors, Chem. Eng. Sci. 28 (6) (1973) 1257–1272. [34] G. Leu, R. Baratti, An extended Kalman filtering approach with a criterion to set its tuning parameters, Comp. Chem. Eng. 23 (2000) 1839–1849. [35] H. Kwakernaak, R. Sivan, Linear Optimal Control Systems, Wiley, New York, 1972. [36] A. Papoulis, Probability, Random Variables and Stochastic Processes, McGrawHill, New York, 1965. [37] P. Gonzalez, J. Alvarez, Combined proportional/integral-inventory control of solution homopolymerization reactors, Ind. Eng. Chem. Res. 44 (2005) 7147–7163. [38] J. Alvarez, R. Suárez, C. Sánchez, Nonlinear decoupling control of free-radical polymerization continuous stirred tank reactor, Chem. Eng. Sci. 45 (1990) 3341–3354. [39] J. Alvarez-Ramírez, J. Alvarez, A. Morales, An adaptive cascade control for a class of chemical reactors, Int. J. Adapt. Control 16 (2002) 681–701. [40] A. Schaum, J.A. Moreno, J. Diaz-Salgado, J. Alvarez, Dissipativity-based observer and feedback control design for a class of chemical reactors, J. Process Control 18 (2008) 896–905. [41] J.E. Bailey, D.F. Ollis, Biochemical Engineering Fundamentals, McGraw Hill Higher Educational, Tokio, 1986. [42] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, New York, 1990. [43] L. Perko, Differential Equations and Dynamical Systems, Springer-Verlag, New York, 1991. [44] M.W. Hirsch, S. Smale, Differential Equations Dynamical Systems and Linear Algebra, Academic Press Inc., London, 1974. [45] S. Lefschetz, Differential Equations: Geometric Theory, Dover Publications, New York, 1977. [46] J.E. Slotine, W. Li, Applied Nonlinear Control, Prentice Hall, Upper Saddle River, NJ, 1991. [47] P.A. Markowich, C. Villani, On the trend to equilibrium for the Fokker–Planck equation: an interplay between physics and functional analysis, Mat. Contemp. 19 (1999) 1–29. [48] T.D. Frank, Lyapunov and free energy functionals of generalized Fokker–Planck equations, Phys. Lett. A 290 (2001) 93–100. [49] H.F. Weinberger, A First Course in Partial Differential Equations, John Wiley and Sons, New York, 1965. [50] I. Prigogini, From Being to Becoming, W.H. Freeman and Company, San Francisco, 1980. [51] R. Bonifacio, L. Lugiato, J.D. Farina, L.M. Narducci, Long time evolution for a onedimensional Fokker–Planck process: application to absorptive optical stability, IEEE J. Quantum Electron 17 (1981) 357–365. [52] C.W. Gardiner, Handbook of Stochastic Methods, Springer-Verlag, Germany, 1997. [53] D.A. Crowl, J.F. Louvar, Chemical Process Safety. Fundamentals with Applications, Prentice Hall PTR, Upper Saddle River, US, 2002. [54] A. Favache, D. Dochain, Thermodynamics and chemical systems stability: the CSTR case study revisited, J. Process Control 9 (2009) 1324–1332.