Dissipative Boundary Control for an Adiabatic Plug Flow Reactor With Mass Recycle

In this contribution the stability properties and regulation of a class of convective systems described by first order hyperbolic partial differential equations with boundary recycle is studied. The system’s setting is consistent with the first and second laws of thermodynamics, allowing to use the entropy functional as a storage function and the internal entropy production as the dissipation similarly to Hamiltonian systems, usually not well defined for systems with mass flows. It is found that the difference of the entropy evaluated at the boundaries is directly proportional to the supply rate, fulfilling the dissipation inequality. Furthermore, the dynamics of the entropy balance allow to define a saturated Proportional-Integral controller with a cascade structure: The inner loop tracks an entropy reference, while the outer loop regulates a process variable. The regulation is achieved with a lumped actuator, using continuous measurements at the boundaries. The controller is applied to an adiabatic plug flow reactor with a recycle of the output stream, a configuration known to be potentially unstable with dissociation reactions. Finally, the controller is tested to track a set point facing several disturbances using the recycle rate as the control variable.


I. INTRODUCTION
Reference tracking of distributed parameter systems based on nonlinear PDEs such as plug flow reactors (PFR) has been an active field since the advent of nonlinear control techniques [2]. Several techniques have been used to track references in distributed parameter systems: sum-of-squares and fuzzy control [3], [4], PI controllers for a class of semilinear hyperbolic systems [5] or linearized regulators [6]. Nonetheless, the regulation and stability analysis of infinitedimensional systems using passivity or thermodynamic properties is an open field and an opportunity to revisit problems with more insight of the phenomena [7], [8].
A class of dynamical systems using physical properties gained importance when [9] proposed the concepts of dissipation, storage function (or Lyapunov function) and supply rate. The criteria of stability and regulation using passivity or dissipativity properties is of wide use and has lead to several well known techniques in electrical and/or mechanical The associate editor coordinating the review of this manuscript and approving it for publication was Yiqi Liu . systems, such as the Interconnection Damping Assignment Passivity Based Control (IDA-PBC) for lumped parameter systems [10] where the energy domain can be well defined. However, energy flows and dissipation as defined by IDA-PBC are not equivalent when compared with systems with mass flows. If dissipation conditions are to be defined for systems with mass flows, the entropy balance should be be considered [11]- [13]. The internal entropy production as a stability criteria for physics-based dynamical systems was proposed in [14] for lumped parameter systems. The internal entropy production was tested as a Lyapunov function candidate in isolated systems in order to identify potential instabilities.
The extension of the dissipation concept to infinite dimensional systems is well established for the linear case with the port-Hamiltonian canonical structure, which has been adapted to represent several physical systems [15]. The canonical form of the latter has similar definitions respect to the lumped parameter counterpart, such as supply rate and a Hamiltonian function [13], [16], [17], where there have been contributions towards modeling using VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ thermodynamic properties [7], [8], [18]. Also, it is of paramount importance the numerical method of choice, since it might influence in the energy preservation, altering the stability properties [19]- [22]. For distributed parameter nonlinear systems under closed loops, is hard or even impossible to prove the uniqueness of the solutions. Nonetheless, there are novel applications towards formulations that satisfy the first law of thermodynamics, such as the development of boundary actuators with a nonlinear nature in [23]. Also, the dissipation properties postulated by the second law of thermodynamics have been used together with optimization techniques, yielding to novel applications in process engineering [24]- [27]. Systems with mass feedback are a classical process engineering case study [28] are still attractive due to their complexity, since these are potentially unstable in certain scenarios, even for simple processes such as heat exchange [29], [30]. The above motivates to revisit the closed loop problem of a constant velocity adiabatic PFR with mass feedback, now using thermodynamical concepts that lead to a dissipative formulation. By manipulating the recycle rate, and indirectly increasing the residence time of a fraction of the stream, it is possible to manipulate the supply rate and regulate a desired process variable without linearizing the model, using a boundary PI controller with a cascade structure, where the inner loop tracks an entropy rate reference while the outer loop tracks a process variable, providing robustness against strong disturbances [31]. Thus, in this paper a boundary controller is defined by using the passivity properties of the entropy function, that by construction allows to regulate a state variable of a class of first order hyperbolic systems. The organization of this document is as follows: In Section II the class of first order hyperbolic systems is defined. In Section III, the dissipation based on the second law is described, defining the controller. In Section IV an adiabatic PFR with a particular stoichiometry is analysed and the controller is tested. Finally, in Section V the performance of the controller, its formulation and future work is discussed.

II. THERMODYNAMICS OF CONVECTIVE SYSTEMS
Consider a first order hyperbolic system where diffusion/dispersion phenomena can be neglected and a macroscopic thermodynamic process is taking place. These class of systems are also known as convective systems. The dynamics is represented by where η (z, t) , ζ (z, t) ∈ H n [(0, L) , R n ] are the extensive and intensive (or conjugated) properties state (co-states) vectors, respectively., being H n [(0, L) , R n ] the infinite dimensional Hilbert space on n−dimensional functions defined on the interval [0, L], the spatial coordinate is z ∈ [0, L] ⊂ R + and t ∈ [0, ∞) respectively, v = F/A ∈ R + is the constant linear velocity, F is a constant volumetric flow rate and A is a constant cross-sectional area of the system, the matrix M ∈ R m×n may contain stoichiometric coefficients and/or heat transfer coefficients for transport phenomena and f : H n [(0, L) , R n ] → R n is the vector field of the thermodynamic phenomena. The system has the following initial conditions: where α (z) ∈ H n [(0, L) , R n ]. The nonlinear boundary condition with mass feedback is where η in (t) ∈ R n are the inlet conditions and u (t) ∈ (0, 1) is the recycle rate, which can be used as the control variable, assumed as smooth functions of time.

A. ENTROPY BALANCE
The entropy density function or local entropy, s (t, z), is a real valued function and a state variable, therefore it can be computed from a map of the extensive variables, s = φ (η), where φ (η) might depend on the set or a subset of the state variables that define the state of the system [32]. Assuming local equilibrium [33], the Gibbs relation that describes the differential behavior of the entropy density at a given point z is given by where µ is the chemical potential vector, defined as µ = µ o + ln (a) with µ o as the reference chemical potential vector and ln (a) is the logarithm of the activities vector, defined as a i = x i γ i , where x i = c i /c T is the molar fraction, c T = 1 T c is the total concentration, and γ i is an activity coefficient. Notice that the Gibbs relation considers an isochoric system where P/Td (v) = 0, wherev is a volume fraction, thus the entropy does not depend on a complete set of extensive variables. The extensive properties are volumetric, then the functional J i (t) = A z 0 η i (λ, t) dλ is the total (extensive) state variable along the spatial domain and λ is a dummy variable, where the upper limit is usually the boundary z = L.
The local entropy gradient of a convective system, or intensive properties vector, is defined as being the entropy at least two times differentiable, reaching a maximum value at the thermodynamic equilibrium in isolated systems. The concavity must hold Since in (4) can be written locally as ds = ζ T dη (Gibbs relation), the entropy balance is given by ζ T ∂η/∂t, and using (1), it leads to 30940 VOLUME 10, 2022 where the first product on the rhs is the entropy transport, ζ T M T is the vector of local driving forces and is the local internal entropy production, which is equal to zero if the driving forces are zero (thermodynamic equilibrium). The total entropy of the system is computed by integrating the entropy density through the spatial domain, and by multiplying by the constant cross-sectional area A: which leads to the entropy balance and is the total internal entropy production, which is not necessarily a bilinear expression as defined by the Onsager relations [33]. Notice that the total internal entropy production in (10) is positive definite and zero only at the thermodynamic equilibrium.
The initial condition of (10): where φ (·) is a map of the extensive properties, 1 and the entropy density at the boundaries are 2

B. ENTROPY PRODUCTION AND ITS DYNAMICS
In distributed parameter systems, the internal entropy production is a functional of the local internal entropy production, whose dynamic balance is and, by expanding the gradient 1 If a subset of the total extensive properties is used, then s = ζ T η + q (η), where q (·) depends on a state equation [32], therefore the product ζ T η does not necessarily yield an entropy. 2 The subscripts 0 and L will indicate values at the boundaries, where DFOI stands for ''Dissipative First Order Interaction'', negative according to the entropy Hessian, and DSOI for ''Dissipative Second Order Interaction'' as defined in [14], which contains the gradient of the vector field, and also the driving forces, whose magnitude may change its sign for nonlinear phenomena far from the equilibrium, (e.g. chemical reactions).

III. DISSIPATIVITY PROPERTIES AND BOUNDARY CONTROL
It is postulated that the entropy function is concave for isolated systems [34]. By considering the local equilibrium postulate of classical non-equilibrium thermodynamics [33], by construction the entropy has the same properties of a storage function for a thermodynamically consistent system [11], [35]. The following proposition will give the stability properties of the system.

Proposition 1: Consider a first order hyperbolic thermodynamic system as defined in (1) with constant flow/velocity. If the vector field considers transport expressions such that only at the thermodynamic equilibrium their flux is equal to zero, then the entropy functional is a storage function and the internal entropy production functional defines a dissipation.
Proof: We define a differential with the opposite sign of the entropy with the auxiliary variable wheres andS = ALs are constant and arbitrarily large values of entropy, hence w (z, t) ≥ 0 ∀t ≥ 0, z ∈ [0, L] and W (t) ≥ 0 ∀t ≥ 0, respectively. Then, the entropy negative balance in (10) is and since (t) ≥ 0, in (20) is dissipative and satisfies the dissipation inequality, with (t) as the dissipation function, with integral form The rhs of (15), describes the supply rate, which is bounded and directly proportional to the system velocity which can be seen as a bounded supply rate, since for a closed system s L (t) = s 0 (t) and dW (t) dt = − (t). The variable VOLUME 10, 2022 W (t) allows then to formulate the dissipativity in its original definition [36]. A significant difference with respect to the Hamiltonian systems is that the supply rate ω is not the product of a conjugated force and an effort, rather the entropy function is a map of the extensive variables. From (17) it can be seen that the minimum supply rate ω min = 0, which means full recycle u (t) = 1, while the maximum supply rate is The entropy functional itself cannot provide stability information, nor reveal potentially unstable trajectories. Instead, the internal entropy production functional has proved to be the good criteria for stability analysis [14].

A. BOUNDARY CONTROL
In many applications, distributed actuators are not available or are discrete, due to the infinite-dimensional nature of the spatial domain [2], [3]. The following derivations assume that the boundaries are continuously measured. This implies that potential disturbances coming through the boundaries, d (t) := η in (t) are known. Due to the constant velocity, the mean residence time of the reactants can be modified with the recycle rate. Then, the extensive properties at the boundary hence (11) is We assume that when the fresh stream at z = 0, d (t), and the recycle from z = L are well mixed and together they create a new boundary condition. Nonetheless, since there is mixing at z = 0, this implies a s mix , and the entropy at the boundary z = 0 is not a linear combination of the properties at the boundaries. The entropy at z = 0 is a function of the input properties (disturbances), the recycle rate and the properties at z = L, namely s 0 (t) = φ d (t) , u (t) , η L (t) . By using the negative of the entropy, W (t), considering the input properties and the control variable (recycle rate) can be rewritten aṡ (20) and the controlled variable tends to the unity, the supply rate tends to zero, as shown in (15).
Note that (20) is not an explicit function of the control variable, since the supply rate is a map of u (t). Nonetheless, when (20) is derived with respect to time yields where˙ (t) is the total entropy production rate, with unknown sign. Besides, notice that the control variable u (t) and its time derivativeu (t) appear. In order to formulate the dynamics of the control variable, the terms that do not includė u (t) in (21) are grouped into a new auxiliary variable which might be unknown at some extent due to the dynamics of the internal entropy production or the rate of change of the disturbances.
In addition, the dissipation function defines a variable χ (t) :=Ẇ , such thatẄ in (21) together with the uncertain terms in (22) becomė By construction, a cascade structure regulates χ (t) with an inner loop around the dynamic reference χ ref (t), defining (23) as a reference and its erroṙ where e 1 (t) = χ (t) − χ ref (t) and κ 1 > 0 is a tunable constant. This defines the error dynamics asė 1 (t) = −κ 1 e 1 (t), and the integral form of (24) is 1 t) , and for a sufficiently large κ 1 guarantees that The dynamics of the control variable is a function of the error and the conditions at the boundaries such thaṫ where It is important to note that the intensive properties at the z = 0 boundary, ζ 0 (t), are not a linear combination of the disturbance and the recycle conditions at the exit of the PFR. On the other hand, the extensive properties η L (t) − d (t) are the difference between the properties at both boundaries before mixing. Even though this term is the inner product of the vector of intensive and extensive properties, its result is not an entropy function. The reference dynamicsχ ref (t), is considered to include the PI controller in the dynamics of the manipulated variable as followṡ where θ sat (t) is a saturation function and the error: 30942 VOLUME 10, 2022 depends on the measured variable y (t) and its time varying reference, y ref (t); K p is a proportional gain and τ I is the reset time of the integral action. The dynamics of the manipulated variable in (26) becomė Since the entropy function is concave, W (t) must be convex. As mentioned, the sign of˙ (t) cannot be determined for all driving forces, hence˙ (t) must be computed for each physical system [32]. However, when a process is near the thermodynamic equilibrium,˙ (t) is negative definite, being a Lyapunov functional [7], [8].
The regulation of an extensive variable at the boundary z = L is accomplished with (29), where the estimation of uncertain terms, such as˙ , is carried out by the integral action. Note that the controlled variable, u (t), is obtained by integrating (29), therefore the PI controller is in fact implemented on the dynamics of the controlled variable. On the other hand, the use of reference (28) allows by construction to cancel the dynamics of a inner loop in (24), meaning that in order to follow a process variable, an entropy is being tracked.

SATURATION CONSTRAINTS
PI control offers robustness as an advantage, however, in order to avoid performance issues, it is convenient to define a saturation constraint. For this particular reacting dissipative system with mass recycle, a saturation is proposed such that upper and lower limits are. Thus, the controller structure in (29) is modified by multiplying the controller inu by a saturation function θ sat , where The tracked reference, either concentration or temperature (proportional to internal energy) must be bounded between two physically meaningful values: the minimum conversion reached when u = 0 (open loop), and the equilibrium conditions when u = 1 and the PFR operates as a batch reactor (full recycle). In order to avoid recycle rates that fall into these last scenarios, the restriction imposed in (30) establishes an upper limit and a lower limit, namely u + and u − .

IV. CASE STUDY
Let consider a reversible reaction system taking place in an adiabatic PFR, where the boundary condition at z = 0 is a mixture of a fresh stream and a fraction of the flow at z = L. The system is depicted in Figure 1 and is assumed to have no transport delays, that is, the mass from z = L has the same properties at the mixture point at z = 0. The saturated regulator is applied using the parameters given in Table 1.

A. FIRST LAW BALANCES
Consider the following single reversible reaction 3 : where E 1,2 are the chemical species and T = −1 2 = B − A is the vector of stoichiometric coefficients with sign, being A T = 1 0 and B T = 0 2 the vectors of reactants and products, respectively. The component balance is: where r is a specific reaction rate and C T (z, t) = C 1 C 2 is the concentration vector. The total mass is conserved and if (32) is multiplied by the molecular mass M w , the total mass balance m T is given by the transport equation due to and since other types of energy are neglected. The internal energy is the total energy, thus where U (z, t) = C T u is the total internal energy density and is the internal energy. By assuming constant heat capacities, the temperature dynamics is nonetheless, we prefer balances of the extensive properties. The boundary conditions of the system are For simplicity, it is assumed that all chemical species behave ideally, i.e. incompressible liquids and diluted nonionic chemical species, such that no excess functions are  In a), the peaks in the entropy production due to the chemical reaction together with the mass recycle define a different steady state profile, for instance, in the mole fraction X 1 as seen in b) [18]. required and the activities are equal to the mole fraction.These assumptions lead to a mass action rate obtained from mesoscopic non-equilibrium thermodynamics [37] is where T (z, t) is the local absolute temperature, µ o (z, t) is the reference chemical potential, K eq = exp − T µ o /RT is the equilibrium constant at a given temperature, ε 0 is a mesoscopic kinetic parameter, E a0 ε 0 T µ o is an activation energy and D R is an effective diffusion coefficient. The expression r (z, t) = 0 holds only at T µ = A = 0, i.e. at thermodynamic equilibrium.

B. ENTROPY, ENTROPY PRODUCTION AND STABILITY
The driving force of a chemical reaction is the affinity, A/T = T µ/T and the flux is given by −r A , therefore the FIGURE 3. Closed loop of a system under relatively moderate ramp disturbances. In a), there is a ramp for the reactant concentration with positive slope in t d 1 ≤ t ≤ t d 2 , then the ramp has a negative slope in t d 2 ≤ t ≤ t d 3 . In b), the manipulated variable counteracts this disturbance, and in c), the magnitude of t decreases as the manipulated variable gets closer to full recycle.
negative of the entropy balance is The Hessian of the PFR under study is wherec v = c T v C is the local mean heat capacity, with The entropy production is σ = −rA/T > 0 and its dynamic equation is ∂σ/∂t = −v∂σ/∂z + DFOI + DSOI, being In a), distributed behavior of C 2 z, t for a system with relatively moderate ramp disturbances. The controller smoothly follows the reference C 2,ref at the boundary rejecting the disturbances. In b) the behavior of σ z, t shows mild variations. , the reaction of the manipulated variable is kept within the limits defined by u + and u − (red lines) in order to prevent batch operation, perhaps reaching the thermodynamic equilibrium, c) which can make t /F tend to zero, since the supply rate ω t (not shown, but qualitatively similar to t ) also tends to zero.
where r u = T u is the local internal energy change, E a0 is an activation energy, A is the reaction affinity and 1 is a vector of ones of proper dimensions. DFOI is a negative definite function, and DSOI is a function whose sign can be positive for large driving forces (far from thermodynamic equilibrium) and its magnitude may be greater than that of the DFOI thus being dominant, which is the case for chemical reactions. The mesoscopic conservative term E a0 RT is the only term in the DSOI that causes this sign shift in the entropy production dynamics [14], [37].
First order hyperbolic systems are expected to be stable. Nonetheless, the stoichiometry of the dissociation reaction along with the mass recycle can lead to high parametric sensibility with respect to the recycle rate. For instance, in Figure 2 it is shown how a tiny (with practical effects, infinitesimal) change in the open loop recycle rate generates a new steady state profile along the PFR. On this sense, a linearized version of a reacting system with mass recycle was studied in [38], in which an irreversible dissociation reaction with the same stoichiometry was used, leading to a FIGURE 6. In a), the distributed behavior of C 2 z, t for a periodic disturbance between times t d 1 < t < t d 2 . In a), the controlled variable reaches the set point after the disturbances pass and b), the local internal entropy production has a couple of peaks when the chemical reaction is reactivated.
spectral first order hyperbolic system with eigenvalues on the right half complex plane.

C. CLOSED LOOP
The product concentration at z = L, C 2,L , was regulated using the recycle rate u (t) as the control variable using controller in (29), that for the particular case study has the form, being the error e (t) =

1) CASE I: DISTURBANCE REJECTION WITHOUT USING THE SATURATION
The magnitude of this rejection is such that it causes relatively moderate changes in the inlet concentrations, mainly in the reactant concentration. Notice how the behaviour of the control variable u (t) resembles to the total entropy S (t). This is visible in Figure 3, where the qualitative behavior of both variables is quite similar. The distributed behavior of C 2 (z, t) and σ (z, t) is shown in Figure 4, where there are no noticeable peaks.

2) CASE II: DISTURBANCE REJECTION USING THE SATURATION
This Case has a periodic disturbance with a magnitude such that the saturation is required at some points. Figure 5 shows that when the system is close to full recycle, the manipulated variable reaches the saturation conditions preventing batch operation, which is reflected in c) of the same Figure. Recall that the PI controller is divided by the function (t). Furthermore, Figure 6 shows the distributed behavior of the regulated variable and the internal entropy production with the expected behavior, that is, higher conversion, higher entropy production. When the recycle rate is close to full recycle, the supply rate ω is close to zero, thus the entropy production is close to zero. When u (t) decreases again, allowing fresh stream into the PFR, another peak in σ (z, , t) can be observed. Once the disturbance ends, the systems reaches the set point again without offset. In both I and II Cases, the qualitative behavior of the entropy functional, S (t), was similar to the disturbance and the manipulated variable behavior, while the entropy production functional, (t), has a similar qualitative behavior to (t). See Figure 7 (only shown for Case II).

V. CONCLUSION
A dissipative PI boundary controller with saturation constraints is proposed and tested for a class of convective reacting systems with mass recycle, known to be a control challenge due to potential instabilities. The entropy allows to define a dissipation functional and the entropy convective flows can be associated to a supply rate, such that a process variable can be tracked while dissipating an entropy potential until a new steady state is reached. This property can be easily exploited for systems that are shown to be dissipative for any driving force, but it can become difficult when stability is compromised with DSOI, which is the case for chemical reactions with large activation energies. The latter can be overcome by using the robust properties of a PI structure together with a saturation function, allowing the system to smoothly reach a desired set-point while avoiding trivial operating conditions (batch) caused by strong disturbances.
Future work aim to analyse systems with linear fluxes, thus being able to use more powerful estimation and control laws without losing the thermodynamic foundations, and also describing simultaneous phenomena (e.g., reacting systems with heat transfer) for more diverse configurations. ARMANDO CAMPOS-RODRÍGUEZ was born in Jalisco, Mexico. He received the B.S., M.S., and Ph.D. degrees in chemical engineering from the Universidad de Guadalajara, Mexico. He is currently with the Universidad Autónoma de Guadalajara as a Research Professor. His research interests include wastewater treatment, non-linear control, and Tequila process.
EFRÉN AGUILAR-GARNICA was born in Jalisco. He received the B.S. and Ph.D. degrees in chemical engineering from the Universidad de Guadalajara, Mexico. He is currently with the Universidad Autónoma de Guadalajara as the Head of the Research Office and the Chemical Engineering Laboratory. His research interests include catalysis and process control.
J. PAULO GARCÍA-SANDOVAL was born in Mexico City. He received the B.S. and M.S. degrees in chemical engineering from the Universidad de Guadalajara and the Ph.D. degree in electrical engineering from CINVESTAV Guadalajara, Mexico. He has been awarded with the Premio Estatal de Ciencia of the Jalisco State. He is currently the Head of the Process Control Laboratory, CUCEI Campus, Universidad de Guadalajara, where he is a Research Professor. His research interests include fractional calculus, thermodynamics, and nonlinear control.