- Email: [email protected]
Machine Learning-Based Model Predictive Control of Distributed Chemical Processes⁎
3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by 3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by Partial 3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by Partial Available online at www.sciencedirect.com Differential Equation, and XI Workshop Control of Distributed 3rd IFAC/IEEE CSS Workshop on Control of Systems Governed by Partial Differential Equation, and XI Workshop Control of Distributed Parameter Systems Partial Differential Equation, and XI Workshop Control of Distributed Parameter Systems Oaxaca, Mexico, May 20-24, Parameter Systems Differential Equation, and XI 2019 Workshop Control of Distributed Oaxaca, Mexico, May 20-24, 2019 Parameter Systems Oaxaca, Mexico, May 20-24, 2019 Oaxaca, Mexico, May 20-24, 2019 IFAC PapersOnLine 52-2 (2019) 120–127
ScienceDirect
Machine Learning-Based Model Predictive Machine Learning-Based Model Predictive Machine Learning-Based Model Predictive Control of Distributed Chemical Processes MachineofLearning-Based Model Processes Predictive Control Distributed Chemical Control of Distributed Chemical Processes Control Chemical Processes ∗ Distributed ∗ ∗ ∗ Zhe Wuof ∗ Anh Tran ∗ Yi Ming Ren ∗ Cory S. Barnes ∗
Zhe Wu ∗∗ Anh Tran Ren ∗∗ Cory S. Barnes ∗ ∗ ∗ ∗∗ ∗ Yi Ming D. ∗ Siyao Chen Christofides Zhe Wu Anh Tran Yi Ming D. Ren Cory S. Barnes ∗ Panagiotis ∗∗ Siyao Chen Panagiotis Christofides ∗ ∗ ∗ ∗ ∗∗ ∗ ∗∗ Zhe Wu Anh TranPanagiotis Yi Ming D. Ren Cory S. Barnes ∗ Siyao Chen Christofides ∗ Siyao Chen ∗ Panagiotis D. Christofides ∗∗ ∗ Department of Chemical and Biomolecular Engineering, University of of Chemical and Biomolecular Engineering, University of ∗ ∗ Department California, Los Angeles, CA, 90095-1592, USA (email: Department of Chemical and Biomolecular Engineering, University of California, Los CA, USA ∗ Department of Chemical and Biomolecular Engineering, University of California, Los Angeles, Angeles, CA, 90095-1592, 90095-1592, USA (email: (email: [email protected]) [email protected]) ∗∗ California, Los Angeles, CA, 90095-1592, USA (email: Department of Chemical and Biomolecular Engineering and [email protected]) ∗∗ and Biomolecular Engineering and ∗∗ Department of Chemical ∗∗ [email protected]) Department of Electrical and Computer Engineering, University Department of Chemical and Biomolecular Engineering and of Department of Electrical and Computer Engineering, University ∗∗ Department of Chemical and Biomolecular Engineering and of California, Los Angeles, CA 90095-1592, USA (e-mail: Department of Electrical and Computer Engineering, University of California, Los Angeles, CA 90095-1592, USA (e-mail: Department of Electrical and Computer Engineering, University of California, Los Angeles, CA 90095-1592, USA (e-mail: [email protected]) [email protected]) California, Los Angeles, CA 90095-1592, USA (e-mail: [email protected]) [email protected]) Abstract: This work proposes a general framework for linking of a state-of-the-art compuAbstract: This work proposes a general framework for linking of a state-of-the-art computational dynamics (CFD) ANSYS Fluent, and other computing platforms using Abstract: This work proposes a general framework linking a state-of-the-art computational fluid fluid dynamics (CFD) solver, solver, ANSYS Fluent,for and otherof computing platforms using Abstract: This work proposes a general forand linking ofcomputing a state-of-the-art computhe lock fluid synchronization mechanism inANSYS anframework effortFluent, to extend the utilities of CFD solvers from tational dynamics (CFD) solver, other platforms using the lock synchronization mechanism in an effort to extend the utilities of CFD solvers from tational fluid dynamics (CFD) solver, Fluent, and the other computing platforms using the lockmodeling synchronization mechanism inANSYS an and effort to extend utilities of CFD solvers from strictly and design to also control optimization applications. Specifically, phthalic strictly and to also optimization applications. Specifically, phthalic the lockmodeling synchronization mechanism in this an and effort to extend the utilities of CFD solvers from strictly modeling and design design also control control and optimization applications. Specifically, phthalic anhydride (PA) synthesis is to chosen for investigation because of its industrial significance anhydride (PA) synthesis is chosen for this investigation because of its industrial significance strictly modeling and design to also control and optimization applications. Specifically, phthalic and its extreme high exothermicity. Initially, a high-fidelity two-dimensional axisymmetric anhydride (PA) synthesis is chosen for this investigation because of its industrial significance and its extreme high exothermicity. Initially, a high-fidelity two-dimensional axisymmetric anhydride (PA)CFD synthesis isfor chosen forInitially, this investigation of its industrial significance heterogeneous model an industrial-scale is because developed in ANSYS Fluent. Next, and its extreme high exothermicity. a FBR high-fidelity two-dimensional axisymmetric heterogeneous CFD model for an industrial-scale FBR is developed in ANSYS Fluent. Next, and its extreme exothermicity. Initially, a FBR high-fidelity two-dimensional axisymmetric heterogeneous model for an industrial-scale is developed in ANSYS Fluent. Next, the CFD modelCFD ishigh used to explore a wide operating regime of the FBR to create a database, the CFD model is to a operating regime of the to create aa database, heterogeneous CFD model fornetwork an industrial-scale FBR is developed in ANSYS Fluent. Next,a the model is used used to explore explore a wide wide regime oftechniques the FBR FBR to fromCFD which recurrent neural and operating ensemble learning arecreate used todatabase, derive from which recurrent neural network and ensemble learning techniques are used derive a the CFD model is usedneural to explore a wide operating regime oftechniques the FBR to create ato database, homogeneous ensemble regression model using a state-of-the-art application program interface. from which recurrent network and ensemble learning are used to derive a homogeneous ensembleneural regression model using a state-of-the-art application program interface. from recurrent network and ensemble learning areprogram used the tointerface. derive Then,which a model predictive control (MPC) formulation that istechniques designed to drive processa homogeneous ensemble regression model using a state-of-the-art application Then, a predictive control (MPC) formulation that designed to drive the process homogeneous ensemble regression model using state-of-the-art application Then, a model model predictive control (MPC) formulation that is isdeactivation designed toisprogram drive theinterface. process performance to the desired set-point and to aavoid catalyst developed using performance to the desired set-point and to avoid catalyst deactivation is developed using Then, a model predictive control (MPC) that isdeactivation designed toisdrive the process the ensemble regression model. Subsequently, CFD model, the ensemble regression model performance to the desired set-point and formulation to the avoid catalyst developed using the ensemble to regression model. Subsequently, CFD model,deactivation the ensemble regression model performance the combined desired set-point and to the avoid catalyst is developed using and the MPC are to create aa closed-loop system by linking ANSYS Fluent to the ensemble regression model. Subsequently, the CFD model, the ensemble regression model and the MPC are combined to create closed-loop system by linking ANSYS Fluent to the ensemble regression model. the lock CFDsynchronization model, the ensemble regression model SciPy viaMPC a message-passing interface (MPI) with mechanism. Finally, the and the are combined toSubsequently, create a closed-loop system by linking ANSYS Fluent to SciPy via a message-passing interface (MPI) with lock synchronization mechanism. Finally, the and the are combined to create a closed-loop system by demonstrate linking ANSYS Fluent the to SciPy viaMPC adata message-passing (MPI) with lock are synchronization mechanism. simulation generated byinterface the closed-loop system used to the Finally, robustness simulation generated byinterface the closed-loop system used to demonstrate the Finally, robustness SciPy via adata message-passing (MPI) with lock are synchronization mechanism. the and effectiveness of the proposed simulation data generated by theapproach. closed-loop system are used to demonstrate the robustness and effectiveness of the proposed simulation data generated by theapproach. closed-loop system are used to demonstrate the robustness and effectiveness of the proposed approach. © 2019, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. and effectiveness of the proposed approach. Keywords: Phthalic anhydride synthesis, Machine learning-based modeling, Model predictive Keywords: Phthalic anhydride synthesis, Machine learning-based modeling, Model predictive control, CFD modeling, Recurrent neuralMachine network learning-based modeling, Model predictive Keywords: Phthalic anhydride synthesis, control, CFD modeling, Recurrent neuralMachine network Keywords: Phthalic anhydride synthesis, control, CFD modeling, Recurrent neural network learning-based modeling, Model predictive control, modeling, Recurrent neural networkfrom a computational standpoint is well documented in 1. CFD INTRODUCTION 1. INTRODUCTION from a computational standpoint is well documented in literature. Specifically, itstandpoint begins with 1. INTRODUCTION from a computational is the welldevelopment documentedofina literature. Specifically, itstandpoint begins with ofina INTRODUCTION a computational is the welldevelopment documented detailed first-principles model for the FBR followed the literature. Specifically, it begins with the development Modeling, control 1.and optimization of fixed-bed reactors from first-principles model forwith the the FBR followed by by of thea Modeling, control and optimization of fixed-bed reactors detailed literature. Specifically, it begins development of development of a corresponding computationally efficient first-principles model for the FBR followed by thea (FBRs) hascontrol been anand active research area for both academia Modeling, optimization of fixed-bed reactors detailed development of aa corresponding computationally efficient (FBRs) hascontrol been anand active research area for both academia detailed first-principles model for the FBR followed by the development of corresponding computationally efficient Modeling, optimization of fixed-bed reactors model, which is used in the development and industries asan FBRs are building blocks of the petro- reduced-order (FBRs) has been active research area for both academia reduced-order model, which is used in the development and industries as FBRs are building blocks of the petrodevelopment of a corresponding computationally efficient (FBRs) has been an active research area for both academia and industries as FBRs are building blocks of the petroof a model-based control algorithm. Subsequently, the reduced-order model, which is used in the development chemical and refining industries. It is found that more of a model-based control algorithm. Subsequently, the chemical and refining industries. It blocks is found that more of model, which is used inSubsequently, the development and industries as gas-solid FBRsindustries. arecatalytic building ofin theexisting petrofirst-principles model, the reduced-order model and the a model-based control algorithm. than 90-95% chemical and of refining It reactors is found that more reduced-order first-principles model, the model and the than 90-95% of gas-solid catalytic reactors in existing a model-based control algorithm. Subsequently, first-principles model, the reduced-order reduced-order model and the chemical and refining industries. It reactors is found that more control algorithm are integrated to generate a closed-loop production lines are FBRs. The prevalent application of of than 90-95% of gas-solid catalytic in existing control algorithm are integrated to generate a closed-loop production lines are FBRs. The prevalent application of first-principles model, the reduced-order model and the control based algorithm are integrated to generate a closed-loop than of are gas-solid catalytic reactors in existing production lines FBRs.industries The prevalent application of system on the of is FBRs 90-95% in chemical process are because of their system on which which the performance performance of the theaalgorithm algorithm is FBRs in process industries are of their control based algorithm are integrated to generate closed-loop production are FBRs. The prevalent application of system FBRs in chemical chemical process industries are because because of other their evaluated. For that reason, the first-principles model that based on which the performance of the algorithm is simplicity inlines design and operation compared with evaluated. For that reason, the first-principles model that simplicity in design and operation compared with other system based onthat which thedynamics performance of the algorithm is FBRs in chemical process industries are because of other their accurately describes the of the physical FBR evaluated. For reason, the first-principles model that configurations, e.g., fluidized-bed reactors. However, desimplicity in design and operation compared with describes the of configurations, e.g., fluidized-bed reactors. However, de- accurately reason, the first-principles model FBR that accurately describes the dynamics dynamics of the the physical physical FBR simplicity in design and operation compared with other is essentialFor to that the development of control algorithms as spite the simple appearance, FBRs reactors. with highly exothermic configurations, e.g., fluidized-bed However, de- evaluated. to of algorithms as spite the appearance, FBRs with highly exothermic accurately the dynamics ofidea theinphysical FBR is essential essential to the the development development of control control algorithms as configurations, e.g., However, de- is spite the simple simple appearance, FBRs with exothermic well as the describes transition from an novel academia to reactions give rise tofluidized-bed some of the reactors. mosthighly challenging conwell as the transition from an novel idea in academia to reactions give rise to some of the most challenging conessential to the development of industry. control algorithms as spite the simple appearance, FBRs highly exothermic reactions give in rise to someengineering of the with mostowing challenging con- is aa cost-effective implementation in The need for well as the transition from an novel idea in academia to trol problems chemical to extreme cost-effective implementation in industry. The need for trol problems in chemical engineering owing to extreme well as the transition from that an novel ideaused in The academia to reactions givenonlinear risechemical to some of the dynamics, mostowing challenging concomputer-aided toolbox can be to expedite a cost-effective implementation in industry. need for nonlinearity, distributed moving hot trol problems in engineering to extreme computer-aided toolbox that can be to expedite nonlinearity, nonlinear distributed dynamics, moving hot a cost-effective implementation in industry. need has for computer-aided toolbox that first-principles can be used usedThe to expedite trol problems in ofchemical engineering owing moving to extreme development of high-fidelity models spots, high risk thermal runaway, constraints on manonlinearity, nonlinear distributed dynamics, hot athe the development oftoolbox high-fidelity first-principles models has spots, high risk of thermal runaway, constraints on maa computer-aided that can be used to expedite nonlinearity, nonlinear distributed dynamics, moving hot the spots, highinputs risk of thermal runaway, constraints on mabeen recognized by many and has led to the developdevelopment of high-fidelity first-principles models has nipulated and state variables, and limited on-line been recognized by many and has led to the developnipulated inputs and state variables, and limited on-line the development of high-fidelity first-principles models has spots, highinputs riskwith of thermal runaway, and constraints on mament of powerful computing platforms, i.e., computational recognized by many and has led to the developmeasurements high uncertainty long delays Wu been nipulated and state variables, and limited on-line of powerful computing platforms, i.e.,tocomputational measurements with high uncertainty and long delays Wu ment been recognized by many and hasas led the developnipulated inputs and state variables, and long limited on-line fluid dynamic (CFD) solvers, such ANSYS Fluent and ment of powerful computing platforms, i.e., computational and Huang (2003). A standard framework to address the measurements with high uncertainty and delays Wu dynamic (CFD) solvers, such as Fluent and Huang (2003). A standard framework to address the fluid of powerful computing platforms, i.e., computational fluid dynamic (CFD) solvers, such as ANSYS ANSYS Fluent and and measurements with high uncertainty and long delays Wu OpenFoam, which allow process engineering researchers to control aspect of FBRs to improve the process and Huang (2003). A standard framework toperformance address the ment OpenFoam, which allow process engineering researchers to control aspect of FBRs to improve the process performance fluid dynamic (CFD) solvers, such as ANSYS Fluent and Huang (2003). A standard framework to address the control aspect of FBRs to improve the process performance overlook the nitty-gritty details of numerical methods for OpenFoam, which allow process engineering researchers to to suppress the magnitude of the hot-spot temperature overlook the nitty-gritty details of numerical methodsand for and to suppress the magnitude of the hot-spot temperature OpenFoam, which allow process engineering researchers to control aspect of FBRs to improve the process performance the solution of the first-principles model and to devote overlook the nitty-gritty details of numerical methods for and to suppress the magnitude of the hot-spot temperature Financial support from the National Science Foundation and the the solution of the first-principles model and to devote Financial overlook the nitty-gritty details of numerical methods for and to suppress thefrom magnitude of the hot-spot temperature their complete attention to crafting more optimal reactor support the National Science Foundation and the the solution of the first-principles model and to devote Department of Energy is gratefully acknowledged. their complete to more Financial support from the National Science Foundation and the the the first-principles and to reactor devote Department of Energy is gratefully acknowledged. theirsolution completeofattention attention to crafting crafting model more optimal optimal reactor Financial support from the National Science Foundation and the Department of Energy is gratefully acknowledged. their complete attention to crafting more optimal reactor Department of Energy is gratefully acknowledged.
2405-8963 © © 2019 2019, IFAC IFAC (International Federation of Automatic Control) Copyright 120 Hosting by Elsevier Ltd. All rights reserved. Copyright 2019 IFAC 120 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 2019 IFAC 120 10.1016/j.ifacol.2019.08.021 Copyright © 2019 IFAC 120
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
designs and developing algorithms to identify more costeffective operating conditions. However, as ANSYS Fluent and OpenFoam are developed to be stand-alone computing platforms which means that the integration of CFD solvers with other computing platforms is often difficult and yet to be standardized as the source codes are often encrypted, and consequently, CFD solvers are often used solely as the design toolbox. Recently, the community has made several efforts to expand the utilities of CFD solvers by studying the communication between CFD solvers and other computing platforms. For example, in Vaquerizo and Cocero (2018), the CFD solver was connected with Aspen Plus to obtain physical and chemical properties of compounds that do not exist in the ANSYS database. Additionally, in our previous works (Wu et al. (2017)), an in-house optimization subroutine that solves a specific class of quadratic programming problems was developed via user-defined functions (UDF) such that ANSYS Fluent was able to compile the program of model predictive controller (MPC) and integrated it within the CFD model to form a closed-loop model. However, recognizing that MPC optimization problems may be non-convex and NP-hard for complex chemical reaction systems, a better solution is to integrate existing robust optimization solvers for solving large-scale nonlinear optimization problems with CFD solvers. For that reason, this work focuses on developing a general framework for linking ANSYS Fluent CFD solver to various computing platforms using the lock synchronization mechanism that enables the development of a realistic closed-loop system in which a high-fidelity CFD model is used to represent the physical system.
121
by a centralized controller as shown in Fig. 1. Specifically, the FBR is fed with a mixture of air and o-xylene that undergoes a highly exothermic conversion to produce PA leading to the formation of hot spots, which are regulated by the centralized controller to prevent catalyst deactivation and thermal runaway. Previous studies have used computational fluid dynamics (CFD) tools to create high-fidelity model for the FBR of this type to study the sensitivity of the hot-spot temperature and product yield to operating parameters such as the concentration of o-xylene in the feed, the feed flow rate and the jacket temperature profile. We have found that it is reasonable to assume that the FBR is symmetrical in the azimuthal direction, which allows for the generation of a computationally affordable two-dimensional (2D) axisymmetric in space heterogeneous CFD model with the same output of meaningful data as a computationally expensive 3D model. In remainder of this section, the development of a two-dimensional (2D) axisymmetric in space heterogeneous CFD model of the FBR via ANSYS Fluent will be outlined.
Fig. 1. Two-dimensional axisymmetric reactor geometry with outer cooling jacket (blue) divided into four noninteracting discrete zones along the tube.
2. CHEMICAL REACTOR MODEL Phthalic anhydride (PA) synthesis is chosen for this investigation because of its industrial significance since PA is one of the most important intermediates for a variety of polymer products such as polyesters, resins and plasticizer and is a key pharmaceutical ingredient of cellulose acetate phthalate, in addition to its extreme high exothermicity. Specifically, PA is commonly produced by partial oxidation of o-xylene in great excess air over V2 O5 /T iO2 catalysts housed in a multi-fixed-bed reactor, inside which the fixedbed reactors (FBRs) are submerged in cooling molten salt. The reactor design and the operating parameters of the production of PA are optimized to suppress the magnitude of the hot-spot temperatures to prevent explosion, catalyst deactivation and thermal runaway. Despite the off-line efforts, the formation of hot spots in the first section of the FBRs is a known problem due to insufficient cooling, and therefore, the need of high-fidelity models that accurately simulate transient and steady-state response of the reactor is evident. This investigation outlines the modeling (Sections 3 and 4.2) and control (Section 4.3 and 4.4) of a FBR inside the industrial-scale multi-fixed-bed reactor. In the existing production line of PA, the FBR is typically 4 m in length and 2.5 cm in diameter and is packed with V2 O5 /T iO2 catalyst particles. Additionally, the FBR is equipped with an outer cooling jacket, which is divided into four noninteracting discrete equidistant zones along the FBR, each of which jacket temperature is regulated 121
3. COMPUTATIONAL FLUID DYNAMICS MODELING 3.1 Mesh Generation The mesh quality plays an important role in CFD modeling, where an accurate CFD model with high mesh quality can significantly save computational resources and obtain a converged solution with robustness. Due to the axisymmetric geometry property of the fixed-bed catalytic reactor, a two-dimensional (2D) axisymmetric structured mesh is constructed in ANSYS ICEM software as shown in Fig. 1. Specifically, the mesh discretizes the reactor volume into 85211 control volumes, within each of which the first principles model (as shown in Eq. 12) of the FBR is numerically solved via second-order upwind finite volume methods. Additionally, high mesh density is applied near the boundary of the reactor to capture large spatial gradients in that region. With the high-quality structured mesh, evaluated based on the minimum orthogonal factor of 1 and maximum ortho skew of 0, the CFD simulations are not expected to encounter convergence difficulty, and numerical solutions generated from the CFD model are expected to carry small numerical error (Tran et al. (2018)). Additionally, a mesh independent study demonstrates that further increasing the mesh size does not improve the simulation results of CFD models but leads to higher computation time to convergence.
2019 IFAC CPDE-CDPS 122 Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
3.2 2D Axisymmetric Heterogeneous CFD Model
−En ) where n = 1, 2, · · · , 5 (5a) Rg T s r1 = K1 αPAs ; r2 = K2 αPBs ; r3 = K3 αPAs (5b) r4 = K4 αPAs ; r5 = K5 αPCs (5c) KC POs 2 α= KC POs 2 + (K1 + 6.5K3 + 3K4 )PAs + K2 PBs + K5 PCs (5d) where A, B, C and D represent o-XL, o-tolualdehyde (o-TA), phthalide (PH), and PA, respectively, En (kJ kmol−1 ), K0n (kmol m2 kg−1 s−1 N−1 ), Kn (kmol m2 kg−1 s−1 N−1 ) and rn (kmol kg−1 s−1 ) are the activation energy, the pre-exponentionl factor, the rate constant and the intrinsic reaction rate of the n reaction, respectively, α is the fraction of unoccupied oxidized sites, Pis (kPa) is the partial pressure of species i on the catalyst (solid) surface, Ts (K) is the catalyst temperature and KC (kmol m2 kg−1 s−1 N−1 ) is the catalyst reoxidation rate constant. In this kinetic model, r6 is set to zero because within the operating regime optimized for the yield of PA, the oxidation rate of PA to COx has been found to be insignificant (Anastasov (2003)). Furthermore, it is important to note that the partial oxidation of o-xylene does not take place spontaneously as the reactants enter the FBR. In fact, the reactants must be transported from the bulk flow to the catalyst surface, and subsequently, diffuse into catalyst pores and bind to catalyst active sites at which the partial oxidation of o-xylene occurs. Therefore, the effectiveness factor (η) of 0.34 is introduced to account for the effects of transport resistances between the free flowing fluid and the catalyst surface. In addition, it is equally important to recognize that the production of PA only occurs in the catalyst particles, and the reaction rates (Eq. 5) are expressed in terms of the rate of change per unit mass of catalyst; therefore, prior to be integrated into the material and energy conservation equations of the FBR, the reaction rates are converted into the rate of change per unit volume of the FBR, i.e., kmol m−3 s−1 , as follows, Kn = K0n exp(
In this subsection, we outline the modeling strategy to customize ANSYS Fluent to create the 2D axisymmetric heterogeneous CFD model for the FBR, and we begin with the bed porosity profile. It has been found in previous experiments that a bed porosity profile exhibits an oscillatory pattern with a decreasing amplitude as the distance to the reactor wall increases due to different packing structure induced by the wall effect; and this pattern is especially more pronounced in FBRs with small tube-to-pellet ratio N < 10. Therefore, in an effort to obtain an accurate simulation of the flow field profile in the FBR, the radial variation of the catalytic bed porosity must be accounted for, and in this CFD model, the bed porosity profile in the radial direction encoded in a user-defined function (UDF) is modeled as follows (De Klerk (2003)), γ(r) = 2.14r2 − 2.53r + 1 when r ≤ 0.637 (1) = γb + 0.29 exp(−0.6r)[cos(2.3π(r − 0.16))] (2) + 0.15 exp(−0.9r) when r > 0.637 where γ(r) and γb are the bed porosity at a particular position inside the bed and the bed porosity in the absence of wall effect, respectively and r is the dimensionless distance from the wall and is computed as follows, R−y (3) r= dp where R, y and dp represent reactor radius, radial position relative to the reactor axis, and particle diameter, respectively. Specifically, at regions close to the reactor wall (r ≤ 0.637), in which the catalyst packing pattern is more structured, the bed porosity profile is modeled by Eq. 1. While at regions further from the wall (r > 0.637), in which the catalyst packing pattern is more random, the bed porosity profile with the oscillatory pattern is modeled by Eq. 2. By encoding an accurate description of the bed porosity profile in the CFD model, we have made the first step toward building a model that can simulate a representative flow field profile in the FBR. In addition, when the reactants flow into the FBR with a nonuniformly packed catalytic bed and a low average porosity of 0.4, the flow field profile is expected to deviate from the plug flow profile and to have a complex pattern due to coupling effects with other transport phenomena, the flow superficial velocity is expected to decrease, and the pressure drop across the catalytic bed is expected to be significant. To address this problem, ANSYS Fluent porous zone function is enabled, which adds additional source terms into the Naiver-Stokes equation to simulate the effects of the catalytic bed on the flow field profile, and its parameters, i.e., the permeability (αb ) and the inertial loss coefficient (C2 ), encoded in a UDF are computed as follows, 3
dp (γ(r)) 3.5 1 − γ(r) ; C2 = , (4) 150 (1 − γ(r))2 dp (γ(r))3 and the Naiver-Stokes equation can be written as shown below in Eq. 12b. αb =
In this study, the intrinsic partial oxidation rate of oxylene to PA is modeled with a global heterogeneous kinetic model (Anastasov (2003)), in which the reaction rate expressions are written as follows, 122
rno = ηρs (1 − γ(r)) rn (Pif ).
(6)
where ρs (kg m−3 ) is the density of the catalyst particles. Subsequently, the kinetic model is used to develop the source term (SY,i ) for the material conservation equation for the species i as follows, SY,i =
5 Mi δni rno γ(r) n=1
(7)
where Mi (kg kmol−1 ) is the molecular weight of the species i and δni is the stoichiometric coefficient of species i in the reaction n. It is important to note that when a source term is introduced a CFD model in which ANSYS Fluent porous function is utilized, the source term is automatically multiplied by the bed porosity and subsequently integrated in the corresponding equation as shown in Eqs. 12c, 12d and 12e; therefore, a correction factor, 1/γ(r), must be used in the formulation of all source terms developed in this work. To develop the heterogeneous CFD model for the FBR, we need to explicitly develop the first-principle model for the solid phase. In this work, we assume that the partial oxidation of o-xylene is a transport limited process, and the catalyst particles are isothermal. As a result, the first-
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
principle model for the solid phase is reduced to the energy conservation equation as shown, ∂ ((1 − γ(r))ρs CP,s Ts ) = ∂t 5 rio ∆Hi + (1 − γ(r))hf s (Tf − Ts )
(8)
i=1
where CP,s =836.0 (J kg−1 K−1 ) and =540.0 (m2 m−3 ) is the specific constant pressure heat capacity and the external surface per unit volume of catalyst particles, Tf (K) is the bulk temperature, hf s (W m−2 K−1 ) is the overall heat transfer coefficient between the bulk flow and the catalyst surface and ∆Hi (J kmol−1 ) is the enthalpy of reaction of the ith reaction. Specifically, the right hand side of Eq. 8 accounts for the rate of change in the solid temperature due to the partial oxidation of o-xylene and the heat transfer between the free flowing fluid and the catalyst particles, respectively. Next, to integrate the energy conservation equation for the solid phase into the CFD model, the solid temperature is defined as a user-defined scalar (UDS). Upon this, ANSYS Fluent automatically generates a generic convective-diffusive transport equation for this UDS in a porous media of the following form, ∂γρφ + ∇ · (γρvφ) = ∇ · (γΓ∇φ) + γSφ ∂t
(9)
where φ is a UDS, Γ is a diffusion coefficient, and Sφ is the source term, so that our task is to make use of UDFs to transform Eq. 9 into Eq. 8. Specifically, the convective and diffusive terms are zero out, two additional source terms accounting for the two mechanisms that change the solid temperature as discussed are formulated as follows, 5 ro ∆Hi STs ,1 = i=1 i (10) γ(r) (1 − γ(r))hf s (Tf − Ts ) STs ,2 = , (11) γ(r) and the transient term in Eq. 9 is replaced by that of Eq. 8. As a result of the above analysis, the 2D heterogeneous CFD model used to describe the behavior of a catalytic process taking place in the FBR, i.e., the continuity, momentum, energy and species material balances, are given by the following equations: ∂ (γρf ) + ∇ · (γρf v) = 0 ∂t ∂ (γρf v) + ∇ · (γρf vv) = ∂t − γ∇P + ∇ · (γτ ) − γ
(12a) (12b) 2
γµ γ C2 ρf |v|v v+ αb 2
∂ (γρf CP,f Tf ) + ∇ · (γv(ρf CP,f Tf + P )) = (12c) ∂t ∇ · γkf ∇T − hi J i + τ · v − γSTs ,2 i
∂ ((1 − γ)ρs CP,s Ts ) = γSTs ,1 + γSTs ,2 ∂t ∂ (γρf Yi ) + ∇ · (γρf vYi ) = −∇ · (γJ i ) + γSY,i ∂t
(12d) (12e)
123
boundary conditions: z = 0; Yi = Yi,0 , Tf = Tf,0 ∂Tf ∂γρf Yi = 0; =0 r = 0; ∂r ∂r ∂γρf Yi ∂Tf = 0; − λr = hf j (Tf − Tj ) r = R; ∂r ∂r
123
(13a) (13b) (13c)
where v (m s−1 ), P (kPa), ρf (kg m−3 ), CP,f (J kg−1 K−1 ), kf (W m−1 K−1 ) and µ (m2 s−1 are the velocity, static pressure, density, specific constant-pressure heat capacity, thermal conductivity and molecular viscosity of the bulk flow, respectively, ∇ ( i hi J i ) represents the effect of enthalpy transport due to species diffusion, Yi and J i (kg m−2 s−1 ) are the mass fraction and mass diffusion flux, respectively, Tj is the jacket temperature profile, hf j the lumped heat transfer coefficient between the bulk flow and the cooling jacket and τ is the stress tensor. 4. FEEDBACK CONTROL DESIGN In this section, the feedback control system is developed for the phthalic anhydride synthesis reactor to drive the to its setoutlet concentration of phthalic anhydride xoutlet D point xset D and maintain the fluid temperature xT inside the reactor below the maximum allowable temperature xmax . Specifically, we derive a recurrent neural network T model capturing the relationship between fluid temperature xT , species concentration xc in the reactor and the cooling jacket temperature Tji , i = 1, 2, 3, 4. The RNN models are obtained based on training, validation and testing dataset generated from a large number of open-loop CFD simulations with various step changes of Tji , i = 1, 2, 3, 4. Subsequently, MPC is developed using the RNN model fR (·) that will be discussed in the following section as the prediction model with constraints on both jacket temperature and fluid temperature to calculate optimal control actions (i.e., cooling jacket temperature Tji , i = 1, 2, 3, 4) to regulate xoutlet . Finally, the structure D of the entire closed-loop system that integrates MPC solver and CFD dynamic simulation is presented. 4.1 Open-loop Dynamics To build a recurrent neural network model to represent the CFD model of a fixed-bed catalytic reactor for a phthalic anhydride synthesis, open-loop simulations are first conducted to generate the dataset that captures the system dynamics within a wide range of operating conditions. Specifically, by manipulating cooling jacket temperatures Tji , i = 1, 2, 3, 4 of four zones, the dataset is created for fluid temperature xjT and species concentration xjc in the reactor, where j = 1, ..., 8 represents 8 nonequidistant surfaces set in ANSYS Fluent to monitor both xc and xT along the length of the reactor (i.e., 5 of them are set within the first meter of reactor since the maximum fluid temperature occurs in this region). The open-loop CFD simulations are repeated under various step-changes on jacket temperatures to drive xoutlet to D different steady-states of xoutlet under a disturbance-free D environment. The dataset of fluid temperature and species concentration profiles in the reactor are obtained with the initial steady-state as follows: the outlet PA concentration
2019 IFAC CPDE-CDPS 124 Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
= 7.550 × 10−3 and jacket temperatures Tj1 = xoutlet D 640 K, Tj2 = 640 K, Tj3 = 640 K and Tj4 = 640 K. The CFD model is simulated with a sampling period ∆CF D = 0.1 s which represents the integration step size for simulations of the CFD model. Subsequently, timeseries data from open-loop simulations are separated into a large number of time-series data samples with shorter period via a sliding window such that the dimension of time-series data is reduced and more data is generated for the development of an RNN model. Additionally, data preprocessing is employed to remove data samples with repeated or similar dynamic behavior to avoid over-fitting during the training process of RNN. 4.2 Ensemble Regression Modeling In this study, a recurrent neural network (RNN) learning algorithm (Lipton et al. (2015)) and a 10-fold cross validation are used to construct a homogeneous ensemble regression that uses the feed operating parameters, the jacket temperature profile and the real-time state measurements at eight discrete locations to predict the nonlinear distributed dynamics of the FBR. The motivation for using the ensemble regression is four-fold. First, the RNN learning algorithm, which is also known as the backpropagation through time (BPTT), can easily derive RNNs that memorize the training data well, however, without the generalization capability, and therefore, by exposing the learning algorithm to different subsets of the training data, generated by 10-fold cross validation, effectively prevents it from overfitting the RNNs. Second, the learning algorithm is known to be non-convex and is an infamous example of NP-hard problems, and hence, by using different starting initial weight matrices, the learning algorithm might be able to avoid getting trapped in a local minimum and arrive at the set of weights that allows the RNN to accurately approximate the latent function that transforms input sequences to corresponding output sequences. Third, an ensemble regression is known to often outperform its individual constituents and simultaneously accounts for uncertainty in model selection. Fourth, the proposed approach creates a unique opportunity for us to demonstrate our general framework for linking ANSYS Fluent to several powerful APIs that are typically incompatible. For those reasons, the ensemble regression structured to consist of eight sub-ensembles, and each of which is a 10-fold cross-validated committee of RNNs as depicted in Fig. 2 and 3 and is derived to predict the nonlinear distributed dynamics of the FBR. Specifically, Fig. 2 demonstrates the RNN structure where the first column is the input layer, the second and the third columns are hidden layers and the last column is the output layer. In the input layer, Xkij , k = 1, ..., 6 represents the realtime state measurements at eight discrete locations i (i = 1, ..., 8) for j th (j = 1, ..., 10) RNN model, Xkij , k = 7, ..., 12 represents the estimated states from previous location at current time, and Xkij , k = 13, ..., 16 represents jacket temperatures for four zones. In the hidden layer, the nodes calculate the weighted sum of input nodes with nonlinear activation functions and also introduced feedback loops to maintain information in ‘memory’ over time. At the end, the nodes in the output layer ykij , k = 1, ..., 6 report 124
the estimation of future states, and will be used as the predicted states in the MPC in the next section.
Fig. 2. The structure of the recurrent neural network. In this study, each sub-ensemble of the ensemble regression is designed to predict the trajectories of the state variables at a fixed location over a time interval equal to the sampling rate of the centralized controller, which can be an integer multiple of the CFD integration time step (in this work it is set as 50 CFD integration time steps). Specifically, as shown in Fig. 3, the sub-ensemble uses the on-line measurements of the state variables at the designated location along the FBR, the estimated trajectories of the state variables at an upstream location over the elapsed time and the jacket temperature profile to calculate its predictions, which are the average trajectories of those generated by its constituent RNNs. This ensemble integration strategy appears to be rudimentary; nonetheless, the performance of the sub-ensemble measured in terms of the absolute mean error is statistically better than that of its constituent RNNs. The constituent RNNs of individual sub-ensembles are known as the 10-fold cross-validated committee; this is because that 10-fold cross validation is used to generate 10 different subsets of the training data, and each of which is subsequently used to train a RNN. Although the state-ofthe-art API, Keras, can be readily used to create the 10fold cross-validated committee, training RNNs remains a nontrivial task because the learning algorithm is known to suffer from the vanishing and exploding gradient problems. For that reason, the RNNs is designed to have two hidden recurrent layers consisting of 64 recurrent units that use the rectified linear unit function as the activation function (i.e., ReLu function f (x) = max{0, x}) and are initialized as proposed in (Le et al. (2015)). Specifically, Le et al. (2015) demonstrates that this initialization in which the recurrent and bias weight matrices are set to the identity matrix and zero, respectively, allows the trained RNNs to have the consistently comparable performance with the standard long short term memory (LSTM) networks in various applications. On a last note of this Section, we would like to detail the procedure encoded in the ensemble regression which is used to predict the nonlinear distributed dynamics of the FBR over a finite time interval in the future. Specifically, the first sub-ensemble normalizes the input signals, i.e., the feed operating parameters, the real-time state measurements at the first location and the jacket temperature profile, using the statistics of the first training set and, subsequently, transmits the normalized input signals to its cross-validated committee of RNNs. It is critical to note that all input and output signals of individual RNNs are
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
normalized based solely on the statistics of its training set to prevent data leakage from its test set. Next, the constituent RNNs simultaneously compute the trajectories of the normalized state variables over a finite time interval in the future, which are pooled to compute the respective average trajectories. Then, the first sub-ensemble generates its output signals by rescaling the average trajectories of the normalized state variables using the corresponding statistics. The output signals of the first sub-ensemble are used to infer the dynamics of the FBR at the corresponding location and as input signals to the second sub-ensemble in addition to the real-time state measurements at the second location and the jacket temperature profile. All downstream sub-ensembles execute the same procedure, and the outputs of the sub-ensembles are used to estimate the nonlinear distributed dynamics of the FBR. 4.3 Model Predictive Controller In this work, the primary objectives of the feedback control scheme is to regulate xoutlet while mitigating issues D related to catalyst deactivation and thermal runaway. Based on the RNN model derived in Section 4.2, a model predictive controller is developed to drive xoutlet to its D desired set-point xset and ensure the fluid temperature D xT does not exceed the maximum allowable hot-spot max temperature xT . We use x = [xD xT xO ]T ∈ R6 to represent the states of the weighted-averaged crosssectional mole fraction of PA xD (termed averaged PA concentration in the following text), the weighted-averaged cross-sectional fluid temperature xT (termed averaged fluid temperature in the following text) and all the other species concentrations in the reactor. Additionally, fluid temperature and jacket temperatures are bounded by := {0 < xT ≤ xmax } ∈ R and U bound := {umin ≤ xbound i T T max 4 ui ≤ ui } ∈ R (i = 1, 2, 3, 4), respectively to avoid hotspot formation. The optimization problem of MPC is of the following form: tk+N
J = min
u∈S(∆)
s.t.
tk
2 (t) − xset |˜ xoutlet D | D
(14a)
x ˜(tk+i+1 ) = fR (˜ x(tk+i ), u(tk+i )), i = 0, 1, ..., N − 1 (14b) (14c) x ˜(tk ) = x(tk ) u(tk+i ) ∈ U bound , i = 0, 1, ..., N − 1
(14d)
x ˜T (tk+i ) ∈ xbound , i = 0, 1, ..., N − 1 (14e) T where x ˜D (t), x ˜T (t) are the predicted state trajectories for the averaged PA concentration at the outlet and the maximum averaged fluid temperature, respectively. fR (·) represents the RNN model, S(∆) is the set of piecewise constant functions with period ∆, and N is the number of sampling periods in the prediction horizon. In the optimization problem of Eq. 14, the objective func2 tion of Eq. 14a is the sum of |˜ xoutlet (t) − xset D D | over the prediction horizon. The constraint of Eq. 14b is the RNN model fR (·) used to predict the evolution of states including xD (t) and xT (t). Eq. 14c takes the measurements of states at t = tk from the CFD simulation as the initial conditions for the RNN model of Eq. 14b. Eq. 14d defines the input constraints over the entire prediction horizon. Eq. 14e defines state constraints over the prediction hori125
125
zon. Additionally, we assume that the states of the closedloop system are measured at each sampling time. After the optimal solution u∗ (t) of the optimization problem of Eq. 14 is obtained, only the first control action of u∗ (t) is sent to the control actuators to be applied over the next sampling period. Then, at the next instance of time tk+1 := tk + ∆, the optimization problem is solved again, and the horizon is rolled forward one sampling time. Under the MPC of Eq. 14, the closed-loop state xoutlet can D be ultimately driven to its desired set-point xset D in the absence of model mismatch or disturbances (e.g., reactor feed disturbances) if there is sufficient energy to satisfy the constraints of Eq. 14d and Eq. 14e. To account for model mismatch and eliminate the offset between the final steady-state value of xoutlet and xset D D , an additional deviation term e(tk−1 ) = x ˜(tk−1 )−x(tk−1 ) which represents the deviation between the actual state and the predicted state in the last step t = tk−1 is applied in the predictive models of Eq. 14b such that the model mismatch is accounted for in the predictive models (i.e., x ˜(tk+1 ) = fR (˜ x(tk ), u(tk )) + e(tk−1 )). Alternatively, offset can also be eliminated by integrating an integral feedback control scheme with MPC such that the overall control action can be obtained as: uM P C+I (tk ) = uM P C (tk ) + uI (tk ) (Wu et al. (2017)). In our work, the MPC is developed with RNN models via the scientific computing package SciPy in Python. The integration of the optimization problem of MPC and RNN models is demonstrated by the algorithm shown in Fig. 4. Specifically, at t = tk the RNN model and the state measurements are first imported to MPC. Subsequently, the MPC control actions are calculated by solving the optimization problem of Eq. 14. It is noted that in Eq. 14b, the RNN model fR (·) is utilized to predict future states (i.e., state profiles over the prediction horizon t = [tk , tk+N )) by taking state measurement x(tk ) and control actions u(t), t = [tk , tk+N ) as the inputs. In Fig. 4, it is further demonstrated that following the RNN structure shown in (i) Fig. 3, each sub-RNN model (i.e., the RNN model fR (·), i = 1, ..., 8 for each surface in the reactor) not only utilizes state measurements and control actions, but also takes the predicted outputs x(tk+1 ) from the previous RNN model (i−1) fR (·) (i.e., at (i − 1)th surface) as the inputs to predict the system dynamics at ith surface. 4.4 Integrating MPC and Dynamic CFD Simulation In this section, an integrated framework of MPC and dynamic CFD simulation is developed via message passage interface (MPI) to implement MPC within CFD dynamic simulations. Since UDF supports C language only, the co-compilation of UDF and MPC solver depends on the specific programming language used. For example, the optimization problem of MPC is solved in Python environment in our work, and therefore, the MPC solver and UDF are compiled in Python and C environments, respectively, and a communication bridge between the MPC and ANSYS Fluent UDF needs to be developed through a lock synchronization function embedded in the UDF. As shown in Fig. 5, the lock synchronization mechanism is utilized to exchange data between CFD simulation and MPC solver. Specifically, at the end of each sampling period of CFD dynamic simulation, the measurements of
2019 IFAC CPDE-CDPS 126 Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
Fig. 3. The schematic of ensemble learning with RNN models.
Fig. 4. The algorithm of integrating MPC with RNN models where N is the length of prediction horizon and i is the index of RNN models.
process states are first obtained via UDF (e.g., functions C T, F YI are utilized to derive fluid temperature and PA concentration in ANSYS Fluent via user-defined functions (UDF)) and then sent to the MPC solver. Meanwhile, the CFD simulation is forced to wait for the control actions from MPC. Then, the MPC solver is invoked to solve the optimization problem to derive control actions u(t), for the next sampling period t ∈ [tk , tk+1 ) based on received process state measurements at t = tk . Finally, the CFD simulation receives u(tk ) and continues dynamic simulation until the next sampling time step. The above procedure is repeated with new measurements at the end of each sampling period. Therefore, under the integration of MPC and dynamic CFD simulation, jacket temperatures (i.e., manipulated inputs calculated by MPC) can be adjusted at each sampling step such that xoutlet can D be ultimately driven to its set-point in the absence of model mismatch or disturbances. Additionally, it should be noted that actual process states are measured in CFD simulations while normalized process states are utilized in RNN models and MPC solvers. Therefore, data processing function is also embedded in MPC solver to convert it to normal value before sending it back to CFD simulator. 126
Fig. 5. Basic structure of the proposed integrated MPC and CFD dynamic simulation via a lock synchronization mechanism. 5. CLOSED-LOOP SIMULATION RESULTS In this section, the closed-loop performance of a fixedbed catalytic reactor for the phthalic anhydride synthesis under MPC is investigated. All simulation settings of the closed-loop simulations are the same as those used in the previously studied open-loop system with ∆M P C = 5.0 s and ∆CF D = 0.1 s, respectively. We demonstrate the case of set-point tracking under a disturbance-free environment. The initial steady-state is at xoutlet = 7.550 × 10−3 , D and the jacket temperatures are Tj1 = 640 K, Tj2 = 640 K, Tj3 = 640 K and Tj4 = 640 K, respectively. The new set-point for averaged PA concentration at the outlet is −3 xset . It is shown in Fig. 6 that averaged D = 7.573 × 10 PA concentration at the outlet xoutlet is successfully driven D −3 to xset under MPC, and it takes around D = 7.573 × 10 170 s for xoutlet to reach its final value. The dynamic perD formance of the closed-loop system under MPC improves significantly compared to open-loop control, while it takes −3 around 300 s for the system to settle to xset . D = 7.573×10 The corresponding manipulated input profiles are given in Fig. 7, where it is demonstrated that the jacket temperatures of the four zones stay at high temperature (around
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
quickly, and 685 K) for the first 20 s to increase xoutlet D then the jacket temperatures decrease to their steady-state values such that xoutlet can approach its desired set-point D set xset D and ultimately settle down to xD . Additionally, it is demonstrated in Fig. 8 that the fluid temperature profile is maintained below the maximum allowable temperature 773 K under the constraints of Eq. 14e, which avoids the formation of hot-spots that could deactivate the catalyst. 7.58
×10
-3
Closed-loop system under MPC Desired xset D
7.575
xoutlet D
7.57 7.565 7.56 7.555 7.55
0
50
100
150
200
250
300
350
Time (s)
Fig. 6. The outlet PA concentration xoutlet profile of the D closed-loop system under MPC with xset D = 7.573 × 10−3 .
highly exothermic process was utilized to demonstrate the effectiveness of the proposed approach, under which the product yield can be maximized while suppressing the hot-spot temperature to avoid catalyst deactivation and thermal runaway. Initially, a high-fidelity two-dimensional axisymmetric heterogeneous CFD model for an industrialscale FBR was developed in ANSYS Fluent. Then, the open-loop simulation of CFD model was conducted to create a database for a wide operating regime of the FBR, from which recurrent neural network and ensemble learning were used to derive a homogeneous ensemble regression model using the application program interface (API) Keras. Subsequently, based on the RNN models derived from ensemble regression, an MPC control scheme was developed to drive the process outputs to the desired set-points and suppress the magnitude of the hot-spot temperature to avoid catalyst deactivation. The CFD model, the ensemble regression and the MPC were integrated to create a closed-loop system by linking ANSYS Fluent to MPC optimizer in Python via a message-passing interface (MPI) with lock synchronization mechanism. Finally, the closed-loop simulation results demonstrated the robustness and effectiveness of the proposed approach. 7. ACKNOWLEDGMENTS
Jacket temperature: Zone 1
680
127
Jacket temperature (K)
660
Financial support from the National Science Foundation and the Department of Energy is gratefully acknowledged.
640 680
Jacket temperature: Zone 2
660 640 680
REFERENCES
Jacket temperature: Zone 3
660 640
680 Jacket temperature: Zone 4
660 640
0
50
100
150
200
250
300
350
Time (s)
Average cross-sectional temperature of fluid (K)
Fig. 7. The jacket temperature profiles of the closed-loop −3 system under MPC with xset . D = 7.573 × 10 720 700 680 660 640 300 200
Time (s)
100 0
0
1
2
3
4
Distance along the PA reactor (m)
Fig. 8. The fluid temperature profile of the closed-loop −3 system under MPC with xset . D = 7.573 × 10 6. CONCLUSION In this work, a general framework for linking of a CFD solver, ANSYS Fluent and other computing platforms using the lock synchronization mechanism was proposed to extend the utilities of CFD solvers from strictly modeling and design to also control and optimization applications. Specifically, a fixed-bed reactor (FBR) with a 127
Anastasov, A.I. (2003). An investigation of the kinetic parameters of the o-xylene oxidation process carried out in a fixed bed of high-productive vanadia–titania catalyst. Chemical Engineering Science, 58, 89–98. De Klerk, A. (2003). Voidage variation in packed beds at small column to particle diameter ratio. AIChE journal, 49, 2022–2029. Le, Q.V., Jaitly, N., and Hinton, G.E. (2015). A simple way to initialize recurrent networks of rectified linear units. arXiv preprint arXiv:1504.00941. Lipton, Z.C., Berkowitz, J., and Elkan, C. (2015). A critical review of recurrent neural networks for sequence learning. arXiv preprint arXiv:1506.00019. Tran, A., Pont, M., Crose, M., and Christofides, P.D. (2018). Real-time furnace balancing of steam methane reforming furnaces. Chemical Engineering Research and Design, 134, 238–256. Vaquerizo, L. and Cocero, M. (2018). CFD–Aspen Plus interconnection method. Improving thermodynamic modeling in computational fluid dynamic simulations. Computers & Chemical Engineering, 113, 152–161. Wu, W. and Huang, M.Y. (2003). Nonlinear inferential control for an exothermic packed-bed reactor: Geometric approaches. Chemical Engineering Science, 58, 2023– 2034. Wu, Z., Aguirre, A., Tran, A., Durand, H., Ni, D., and Christofides, P.D. (2017). Model predictive control of a steam methane reforming reactor described by a computational fluid dynamics model. Industrial & Engineering Chemistry Research, 56, 6002–6011.
ScienceDirect
Machine Learning-Based Model Predictive Machine Learning-Based Model Predictive Machine Learning-Based Model Predictive Control of Distributed Chemical Processes MachineofLearning-Based Model Processes Predictive Control Distributed Chemical Control of Distributed Chemical Processes Control Chemical Processes ∗ Distributed ∗ ∗ ∗ Zhe Wuof ∗ Anh Tran ∗ Yi Ming Ren ∗ Cory S. Barnes ∗
Zhe Wu ∗∗ Anh Tran Ren ∗∗ Cory S. Barnes ∗ ∗ ∗ ∗∗ ∗ Yi Ming D. ∗ Siyao Chen Christofides Zhe Wu Anh Tran Yi Ming D. Ren Cory S. Barnes ∗ Panagiotis ∗∗ Siyao Chen Panagiotis Christofides ∗ ∗ ∗ ∗ ∗∗ ∗ ∗∗ Zhe Wu Anh TranPanagiotis Yi Ming D. Ren Cory S. Barnes ∗ Siyao Chen Christofides ∗ Siyao Chen ∗ Panagiotis D. Christofides ∗∗ ∗ Department of Chemical and Biomolecular Engineering, University of of Chemical and Biomolecular Engineering, University of ∗ ∗ Department California, Los Angeles, CA, 90095-1592, USA (email: Department of Chemical and Biomolecular Engineering, University of California, Los CA, USA ∗ Department of Chemical and Biomolecular Engineering, University of California, Los Angeles, Angeles, CA, 90095-1592, 90095-1592, USA (email: (email: [email protected]) [email protected]) ∗∗ California, Los Angeles, CA, 90095-1592, USA (email: Department of Chemical and Biomolecular Engineering and [email protected]) ∗∗ and Biomolecular Engineering and ∗∗ Department of Chemical ∗∗ [email protected]) Department of Electrical and Computer Engineering, University Department of Chemical and Biomolecular Engineering and of Department of Electrical and Computer Engineering, University ∗∗ Department of Chemical and Biomolecular Engineering and of California, Los Angeles, CA 90095-1592, USA (e-mail: Department of Electrical and Computer Engineering, University of California, Los Angeles, CA 90095-1592, USA (e-mail: Department of Electrical and Computer Engineering, University of California, Los Angeles, CA 90095-1592, USA (e-mail: [email protected]) [email protected]) California, Los Angeles, CA 90095-1592, USA (e-mail: [email protected]) [email protected]) Abstract: This work proposes a general framework for linking of a state-of-the-art compuAbstract: This work proposes a general framework for linking of a state-of-the-art computational dynamics (CFD) ANSYS Fluent, and other computing platforms using Abstract: This work proposes a general framework linking a state-of-the-art computational fluid fluid dynamics (CFD) solver, solver, ANSYS Fluent,for and otherof computing platforms using Abstract: This work proposes a general forand linking ofcomputing a state-of-the-art computhe lock fluid synchronization mechanism inANSYS anframework effortFluent, to extend the utilities of CFD solvers from tational dynamics (CFD) solver, other platforms using the lock synchronization mechanism in an effort to extend the utilities of CFD solvers from tational fluid dynamics (CFD) solver, Fluent, and the other computing platforms using the lockmodeling synchronization mechanism inANSYS an and effort to extend utilities of CFD solvers from strictly and design to also control optimization applications. Specifically, phthalic strictly and to also optimization applications. Specifically, phthalic the lockmodeling synchronization mechanism in this an and effort to extend the utilities of CFD solvers from strictly modeling and design design also control control and optimization applications. Specifically, phthalic anhydride (PA) synthesis is to chosen for investigation because of its industrial significance anhydride (PA) synthesis is chosen for this investigation because of its industrial significance strictly modeling and design to also control and optimization applications. Specifically, phthalic and its extreme high exothermicity. Initially, a high-fidelity two-dimensional axisymmetric anhydride (PA) synthesis is chosen for this investigation because of its industrial significance and its extreme high exothermicity. Initially, a high-fidelity two-dimensional axisymmetric anhydride (PA)CFD synthesis isfor chosen forInitially, this investigation of its industrial significance heterogeneous model an industrial-scale is because developed in ANSYS Fluent. Next, and its extreme high exothermicity. a FBR high-fidelity two-dimensional axisymmetric heterogeneous CFD model for an industrial-scale FBR is developed in ANSYS Fluent. Next, and its extreme exothermicity. Initially, a FBR high-fidelity two-dimensional axisymmetric heterogeneous model for an industrial-scale is developed in ANSYS Fluent. Next, the CFD modelCFD ishigh used to explore a wide operating regime of the FBR to create a database, the CFD model is to a operating regime of the to create aa database, heterogeneous CFD model fornetwork an industrial-scale FBR is developed in ANSYS Fluent. Next,a the model is used used to explore explore a wide wide regime oftechniques the FBR FBR to fromCFD which recurrent neural and operating ensemble learning arecreate used todatabase, derive from which recurrent neural network and ensemble learning techniques are used derive a the CFD model is usedneural to explore a wide operating regime oftechniques the FBR to create ato database, homogeneous ensemble regression model using a state-of-the-art application program interface. from which recurrent network and ensemble learning are used to derive a homogeneous ensembleneural regression model using a state-of-the-art application program interface. from recurrent network and ensemble learning areprogram used the tointerface. derive Then,which a model predictive control (MPC) formulation that istechniques designed to drive processa homogeneous ensemble regression model using a state-of-the-art application Then, a predictive control (MPC) formulation that designed to drive the process homogeneous ensemble regression model using state-of-the-art application Then, a model model predictive control (MPC) formulation that is isdeactivation designed toisprogram drive theinterface. process performance to the desired set-point and to aavoid catalyst developed using performance to the desired set-point and to avoid catalyst deactivation is developed using Then, a model predictive control (MPC) that isdeactivation designed toisdrive the process the ensemble regression model. Subsequently, CFD model, the ensemble regression model performance to the desired set-point and formulation to the avoid catalyst developed using the ensemble to regression model. Subsequently, CFD model,deactivation the ensemble regression model performance the combined desired set-point and to the avoid catalyst is developed using and the MPC are to create aa closed-loop system by linking ANSYS Fluent to the ensemble regression model. Subsequently, the CFD model, the ensemble regression model and the MPC are combined to create closed-loop system by linking ANSYS Fluent to the ensemble regression model. the lock CFDsynchronization model, the ensemble regression model SciPy viaMPC a message-passing interface (MPI) with mechanism. Finally, the and the are combined toSubsequently, create a closed-loop system by linking ANSYS Fluent to SciPy via a message-passing interface (MPI) with lock synchronization mechanism. Finally, the and the are combined to create a closed-loop system by demonstrate linking ANSYS Fluent the to SciPy viaMPC adata message-passing (MPI) with lock are synchronization mechanism. simulation generated byinterface the closed-loop system used to the Finally, robustness simulation generated byinterface the closed-loop system used to demonstrate the Finally, robustness SciPy via adata message-passing (MPI) with lock are synchronization mechanism. the and effectiveness of the proposed simulation data generated by theapproach. closed-loop system are used to demonstrate the robustness and effectiveness of the proposed simulation data generated by theapproach. closed-loop system are used to demonstrate the robustness and effectiveness of the proposed approach. © 2019, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. and effectiveness of the proposed approach. Keywords: Phthalic anhydride synthesis, Machine learning-based modeling, Model predictive Keywords: Phthalic anhydride synthesis, Machine learning-based modeling, Model predictive control, CFD modeling, Recurrent neuralMachine network learning-based modeling, Model predictive Keywords: Phthalic anhydride synthesis, control, CFD modeling, Recurrent neuralMachine network Keywords: Phthalic anhydride synthesis, control, CFD modeling, Recurrent neural network learning-based modeling, Model predictive control, modeling, Recurrent neural networkfrom a computational standpoint is well documented in 1. CFD INTRODUCTION 1. INTRODUCTION from a computational standpoint is well documented in literature. Specifically, itstandpoint begins with 1. INTRODUCTION from a computational is the welldevelopment documentedofina literature. Specifically, itstandpoint begins with ofina INTRODUCTION a computational is the welldevelopment documented detailed first-principles model for the FBR followed the literature. Specifically, it begins with the development Modeling, control 1.and optimization of fixed-bed reactors from first-principles model forwith the the FBR followed by by of thea Modeling, control and optimization of fixed-bed reactors detailed literature. Specifically, it begins development of development of a corresponding computationally efficient first-principles model for the FBR followed by thea (FBRs) hascontrol been anand active research area for both academia Modeling, optimization of fixed-bed reactors detailed development of aa corresponding computationally efficient (FBRs) hascontrol been anand active research area for both academia detailed first-principles model for the FBR followed by the development of corresponding computationally efficient Modeling, optimization of fixed-bed reactors model, which is used in the development and industries asan FBRs are building blocks of the petro- reduced-order (FBRs) has been active research area for both academia reduced-order model, which is used in the development and industries as FBRs are building blocks of the petrodevelopment of a corresponding computationally efficient (FBRs) has been an active research area for both academia and industries as FBRs are building blocks of the petroof a model-based control algorithm. Subsequently, the reduced-order model, which is used in the development chemical and refining industries. It is found that more of a model-based control algorithm. Subsequently, the chemical and refining industries. It blocks is found that more of model, which is used inSubsequently, the development and industries as gas-solid FBRsindustries. arecatalytic building ofin theexisting petrofirst-principles model, the reduced-order model and the a model-based control algorithm. than 90-95% chemical and of refining It reactors is found that more reduced-order first-principles model, the model and the than 90-95% of gas-solid catalytic reactors in existing a model-based control algorithm. Subsequently, first-principles model, the reduced-order reduced-order model and the chemical and refining industries. It reactors is found that more control algorithm are integrated to generate a closed-loop production lines are FBRs. The prevalent application of of than 90-95% of gas-solid catalytic in existing control algorithm are integrated to generate a closed-loop production lines are FBRs. The prevalent application of first-principles model, the reduced-order model and the control based algorithm are integrated to generate a closed-loop than of are gas-solid catalytic reactors in existing production lines FBRs.industries The prevalent application of system on the of is FBRs 90-95% in chemical process are because of their system on which which the performance performance of the theaalgorithm algorithm is FBRs in process industries are of their control based algorithm are integrated to generate closed-loop production are FBRs. The prevalent application of system FBRs in chemical chemical process industries are because because of other their evaluated. For that reason, the first-principles model that based on which the performance of the algorithm is simplicity inlines design and operation compared with evaluated. For that reason, the first-principles model that simplicity in design and operation compared with other system based onthat which thedynamics performance of the algorithm is FBRs in chemical process industries are because of other their accurately describes the of the physical FBR evaluated. For reason, the first-principles model that configurations, e.g., fluidized-bed reactors. However, desimplicity in design and operation compared with describes the of configurations, e.g., fluidized-bed reactors. However, de- accurately reason, the first-principles model FBR that accurately describes the dynamics dynamics of the the physical physical FBR simplicity in design and operation compared with other is essentialFor to that the development of control algorithms as spite the simple appearance, FBRs reactors. with highly exothermic configurations, e.g., fluidized-bed However, de- evaluated. to of algorithms as spite the appearance, FBRs with highly exothermic accurately the dynamics ofidea theinphysical FBR is essential essential to the the development development of control control algorithms as configurations, e.g., However, de- is spite the simple simple appearance, FBRs with exothermic well as the describes transition from an novel academia to reactions give rise tofluidized-bed some of the reactors. mosthighly challenging conwell as the transition from an novel idea in academia to reactions give rise to some of the most challenging conessential to the development of industry. control algorithms as spite the simple appearance, FBRs highly exothermic reactions give in rise to someengineering of the with mostowing challenging con- is aa cost-effective implementation in The need for well as the transition from an novel idea in academia to trol problems chemical to extreme cost-effective implementation in industry. The need for trol problems in chemical engineering owing to extreme well as the transition from that an novel ideaused in The academia to reactions givenonlinear risechemical to some of the dynamics, mostowing challenging concomputer-aided toolbox can be to expedite a cost-effective implementation in industry. need for nonlinearity, distributed moving hot trol problems in engineering to extreme computer-aided toolbox that can be to expedite nonlinearity, nonlinear distributed dynamics, moving hot a cost-effective implementation in industry. need has for computer-aided toolbox that first-principles can be used usedThe to expedite trol problems in ofchemical engineering owing moving to extreme development of high-fidelity models spots, high risk thermal runaway, constraints on manonlinearity, nonlinear distributed dynamics, hot athe the development oftoolbox high-fidelity first-principles models has spots, high risk of thermal runaway, constraints on maa computer-aided that can be used to expedite nonlinearity, nonlinear distributed dynamics, moving hot the spots, highinputs risk of thermal runaway, constraints on mabeen recognized by many and has led to the developdevelopment of high-fidelity first-principles models has nipulated and state variables, and limited on-line been recognized by many and has led to the developnipulated inputs and state variables, and limited on-line the development of high-fidelity first-principles models has spots, highinputs riskwith of thermal runaway, and constraints on mament of powerful computing platforms, i.e., computational recognized by many and has led to the developmeasurements high uncertainty long delays Wu been nipulated and state variables, and limited on-line of powerful computing platforms, i.e.,tocomputational measurements with high uncertainty and long delays Wu ment been recognized by many and hasas led the developnipulated inputs and state variables, and long limited on-line fluid dynamic (CFD) solvers, such ANSYS Fluent and ment of powerful computing platforms, i.e., computational and Huang (2003). A standard framework to address the measurements with high uncertainty and delays Wu dynamic (CFD) solvers, such as Fluent and Huang (2003). A standard framework to address the fluid of powerful computing platforms, i.e., computational fluid dynamic (CFD) solvers, such as ANSYS ANSYS Fluent and and measurements with high uncertainty and long delays Wu OpenFoam, which allow process engineering researchers to control aspect of FBRs to improve the process and Huang (2003). A standard framework toperformance address the ment OpenFoam, which allow process engineering researchers to control aspect of FBRs to improve the process performance fluid dynamic (CFD) solvers, such as ANSYS Fluent and Huang (2003). A standard framework to address the control aspect of FBRs to improve the process performance overlook the nitty-gritty details of numerical methods for OpenFoam, which allow process engineering researchers to to suppress the magnitude of the hot-spot temperature overlook the nitty-gritty details of numerical methodsand for and to suppress the magnitude of the hot-spot temperature OpenFoam, which allow process engineering researchers to control aspect of FBRs to improve the process performance the solution of the first-principles model and to devote overlook the nitty-gritty details of numerical methods for and to suppress the magnitude of the hot-spot temperature Financial support from the National Science Foundation and the the solution of the first-principles model and to devote Financial overlook the nitty-gritty details of numerical methods for and to suppress thefrom magnitude of the hot-spot temperature their complete attention to crafting more optimal reactor support the National Science Foundation and the the solution of the first-principles model and to devote Department of Energy is gratefully acknowledged. their complete to more Financial support from the National Science Foundation and the the the first-principles and to reactor devote Department of Energy is gratefully acknowledged. theirsolution completeofattention attention to crafting crafting model more optimal optimal reactor Financial support from the National Science Foundation and the Department of Energy is gratefully acknowledged. their complete attention to crafting more optimal reactor Department of Energy is gratefully acknowledged.
2405-8963 © © 2019 2019, IFAC IFAC (International Federation of Automatic Control) Copyright 120 Hosting by Elsevier Ltd. All rights reserved. Copyright 2019 IFAC 120 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 2019 IFAC 120 10.1016/j.ifacol.2019.08.021 Copyright © 2019 IFAC 120
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
designs and developing algorithms to identify more costeffective operating conditions. However, as ANSYS Fluent and OpenFoam are developed to be stand-alone computing platforms which means that the integration of CFD solvers with other computing platforms is often difficult and yet to be standardized as the source codes are often encrypted, and consequently, CFD solvers are often used solely as the design toolbox. Recently, the community has made several efforts to expand the utilities of CFD solvers by studying the communication between CFD solvers and other computing platforms. For example, in Vaquerizo and Cocero (2018), the CFD solver was connected with Aspen Plus to obtain physical and chemical properties of compounds that do not exist in the ANSYS database. Additionally, in our previous works (Wu et al. (2017)), an in-house optimization subroutine that solves a specific class of quadratic programming problems was developed via user-defined functions (UDF) such that ANSYS Fluent was able to compile the program of model predictive controller (MPC) and integrated it within the CFD model to form a closed-loop model. However, recognizing that MPC optimization problems may be non-convex and NP-hard for complex chemical reaction systems, a better solution is to integrate existing robust optimization solvers for solving large-scale nonlinear optimization problems with CFD solvers. For that reason, this work focuses on developing a general framework for linking ANSYS Fluent CFD solver to various computing platforms using the lock synchronization mechanism that enables the development of a realistic closed-loop system in which a high-fidelity CFD model is used to represent the physical system.
121
by a centralized controller as shown in Fig. 1. Specifically, the FBR is fed with a mixture of air and o-xylene that undergoes a highly exothermic conversion to produce PA leading to the formation of hot spots, which are regulated by the centralized controller to prevent catalyst deactivation and thermal runaway. Previous studies have used computational fluid dynamics (CFD) tools to create high-fidelity model for the FBR of this type to study the sensitivity of the hot-spot temperature and product yield to operating parameters such as the concentration of o-xylene in the feed, the feed flow rate and the jacket temperature profile. We have found that it is reasonable to assume that the FBR is symmetrical in the azimuthal direction, which allows for the generation of a computationally affordable two-dimensional (2D) axisymmetric in space heterogeneous CFD model with the same output of meaningful data as a computationally expensive 3D model. In remainder of this section, the development of a two-dimensional (2D) axisymmetric in space heterogeneous CFD model of the FBR via ANSYS Fluent will be outlined.
Fig. 1. Two-dimensional axisymmetric reactor geometry with outer cooling jacket (blue) divided into four noninteracting discrete zones along the tube.
2. CHEMICAL REACTOR MODEL Phthalic anhydride (PA) synthesis is chosen for this investigation because of its industrial significance since PA is one of the most important intermediates for a variety of polymer products such as polyesters, resins and plasticizer and is a key pharmaceutical ingredient of cellulose acetate phthalate, in addition to its extreme high exothermicity. Specifically, PA is commonly produced by partial oxidation of o-xylene in great excess air over V2 O5 /T iO2 catalysts housed in a multi-fixed-bed reactor, inside which the fixedbed reactors (FBRs) are submerged in cooling molten salt. The reactor design and the operating parameters of the production of PA are optimized to suppress the magnitude of the hot-spot temperatures to prevent explosion, catalyst deactivation and thermal runaway. Despite the off-line efforts, the formation of hot spots in the first section of the FBRs is a known problem due to insufficient cooling, and therefore, the need of high-fidelity models that accurately simulate transient and steady-state response of the reactor is evident. This investigation outlines the modeling (Sections 3 and 4.2) and control (Section 4.3 and 4.4) of a FBR inside the industrial-scale multi-fixed-bed reactor. In the existing production line of PA, the FBR is typically 4 m in length and 2.5 cm in diameter and is packed with V2 O5 /T iO2 catalyst particles. Additionally, the FBR is equipped with an outer cooling jacket, which is divided into four noninteracting discrete equidistant zones along the FBR, each of which jacket temperature is regulated 121
3. COMPUTATIONAL FLUID DYNAMICS MODELING 3.1 Mesh Generation The mesh quality plays an important role in CFD modeling, where an accurate CFD model with high mesh quality can significantly save computational resources and obtain a converged solution with robustness. Due to the axisymmetric geometry property of the fixed-bed catalytic reactor, a two-dimensional (2D) axisymmetric structured mesh is constructed in ANSYS ICEM software as shown in Fig. 1. Specifically, the mesh discretizes the reactor volume into 85211 control volumes, within each of which the first principles model (as shown in Eq. 12) of the FBR is numerically solved via second-order upwind finite volume methods. Additionally, high mesh density is applied near the boundary of the reactor to capture large spatial gradients in that region. With the high-quality structured mesh, evaluated based on the minimum orthogonal factor of 1 and maximum ortho skew of 0, the CFD simulations are not expected to encounter convergence difficulty, and numerical solutions generated from the CFD model are expected to carry small numerical error (Tran et al. (2018)). Additionally, a mesh independent study demonstrates that further increasing the mesh size does not improve the simulation results of CFD models but leads to higher computation time to convergence.
2019 IFAC CPDE-CDPS 122 Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
3.2 2D Axisymmetric Heterogeneous CFD Model
−En ) where n = 1, 2, · · · , 5 (5a) Rg T s r1 = K1 αPAs ; r2 = K2 αPBs ; r3 = K3 αPAs (5b) r4 = K4 αPAs ; r5 = K5 αPCs (5c) KC POs 2 α= KC POs 2 + (K1 + 6.5K3 + 3K4 )PAs + K2 PBs + K5 PCs (5d) where A, B, C and D represent o-XL, o-tolualdehyde (o-TA), phthalide (PH), and PA, respectively, En (kJ kmol−1 ), K0n (kmol m2 kg−1 s−1 N−1 ), Kn (kmol m2 kg−1 s−1 N−1 ) and rn (kmol kg−1 s−1 ) are the activation energy, the pre-exponentionl factor, the rate constant and the intrinsic reaction rate of the n reaction, respectively, α is the fraction of unoccupied oxidized sites, Pis (kPa) is the partial pressure of species i on the catalyst (solid) surface, Ts (K) is the catalyst temperature and KC (kmol m2 kg−1 s−1 N−1 ) is the catalyst reoxidation rate constant. In this kinetic model, r6 is set to zero because within the operating regime optimized for the yield of PA, the oxidation rate of PA to COx has been found to be insignificant (Anastasov (2003)). Furthermore, it is important to note that the partial oxidation of o-xylene does not take place spontaneously as the reactants enter the FBR. In fact, the reactants must be transported from the bulk flow to the catalyst surface, and subsequently, diffuse into catalyst pores and bind to catalyst active sites at which the partial oxidation of o-xylene occurs. Therefore, the effectiveness factor (η) of 0.34 is introduced to account for the effects of transport resistances between the free flowing fluid and the catalyst surface. In addition, it is equally important to recognize that the production of PA only occurs in the catalyst particles, and the reaction rates (Eq. 5) are expressed in terms of the rate of change per unit mass of catalyst; therefore, prior to be integrated into the material and energy conservation equations of the FBR, the reaction rates are converted into the rate of change per unit volume of the FBR, i.e., kmol m−3 s−1 , as follows, Kn = K0n exp(
In this subsection, we outline the modeling strategy to customize ANSYS Fluent to create the 2D axisymmetric heterogeneous CFD model for the FBR, and we begin with the bed porosity profile. It has been found in previous experiments that a bed porosity profile exhibits an oscillatory pattern with a decreasing amplitude as the distance to the reactor wall increases due to different packing structure induced by the wall effect; and this pattern is especially more pronounced in FBRs with small tube-to-pellet ratio N < 10. Therefore, in an effort to obtain an accurate simulation of the flow field profile in the FBR, the radial variation of the catalytic bed porosity must be accounted for, and in this CFD model, the bed porosity profile in the radial direction encoded in a user-defined function (UDF) is modeled as follows (De Klerk (2003)), γ(r) = 2.14r2 − 2.53r + 1 when r ≤ 0.637 (1) = γb + 0.29 exp(−0.6r)[cos(2.3π(r − 0.16))] (2) + 0.15 exp(−0.9r) when r > 0.637 where γ(r) and γb are the bed porosity at a particular position inside the bed and the bed porosity in the absence of wall effect, respectively and r is the dimensionless distance from the wall and is computed as follows, R−y (3) r= dp where R, y and dp represent reactor radius, radial position relative to the reactor axis, and particle diameter, respectively. Specifically, at regions close to the reactor wall (r ≤ 0.637), in which the catalyst packing pattern is more structured, the bed porosity profile is modeled by Eq. 1. While at regions further from the wall (r > 0.637), in which the catalyst packing pattern is more random, the bed porosity profile with the oscillatory pattern is modeled by Eq. 2. By encoding an accurate description of the bed porosity profile in the CFD model, we have made the first step toward building a model that can simulate a representative flow field profile in the FBR. In addition, when the reactants flow into the FBR with a nonuniformly packed catalytic bed and a low average porosity of 0.4, the flow field profile is expected to deviate from the plug flow profile and to have a complex pattern due to coupling effects with other transport phenomena, the flow superficial velocity is expected to decrease, and the pressure drop across the catalytic bed is expected to be significant. To address this problem, ANSYS Fluent porous zone function is enabled, which adds additional source terms into the Naiver-Stokes equation to simulate the effects of the catalytic bed on the flow field profile, and its parameters, i.e., the permeability (αb ) and the inertial loss coefficient (C2 ), encoded in a UDF are computed as follows, 3
dp (γ(r)) 3.5 1 − γ(r) ; C2 = , (4) 150 (1 − γ(r))2 dp (γ(r))3 and the Naiver-Stokes equation can be written as shown below in Eq. 12b. αb =
In this study, the intrinsic partial oxidation rate of oxylene to PA is modeled with a global heterogeneous kinetic model (Anastasov (2003)), in which the reaction rate expressions are written as follows, 122
rno = ηρs (1 − γ(r)) rn (Pif ).
(6)
where ρs (kg m−3 ) is the density of the catalyst particles. Subsequently, the kinetic model is used to develop the source term (SY,i ) for the material conservation equation for the species i as follows, SY,i =
5 Mi δni rno γ(r) n=1
(7)
where Mi (kg kmol−1 ) is the molecular weight of the species i and δni is the stoichiometric coefficient of species i in the reaction n. It is important to note that when a source term is introduced a CFD model in which ANSYS Fluent porous function is utilized, the source term is automatically multiplied by the bed porosity and subsequently integrated in the corresponding equation as shown in Eqs. 12c, 12d and 12e; therefore, a correction factor, 1/γ(r), must be used in the formulation of all source terms developed in this work. To develop the heterogeneous CFD model for the FBR, we need to explicitly develop the first-principle model for the solid phase. In this work, we assume that the partial oxidation of o-xylene is a transport limited process, and the catalyst particles are isothermal. As a result, the first-
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
principle model for the solid phase is reduced to the energy conservation equation as shown, ∂ ((1 − γ(r))ρs CP,s Ts ) = ∂t 5 rio ∆Hi + (1 − γ(r))hf s (Tf − Ts )
(8)
i=1
where CP,s =836.0 (J kg−1 K−1 ) and =540.0 (m2 m−3 ) is the specific constant pressure heat capacity and the external surface per unit volume of catalyst particles, Tf (K) is the bulk temperature, hf s (W m−2 K−1 ) is the overall heat transfer coefficient between the bulk flow and the catalyst surface and ∆Hi (J kmol−1 ) is the enthalpy of reaction of the ith reaction. Specifically, the right hand side of Eq. 8 accounts for the rate of change in the solid temperature due to the partial oxidation of o-xylene and the heat transfer between the free flowing fluid and the catalyst particles, respectively. Next, to integrate the energy conservation equation for the solid phase into the CFD model, the solid temperature is defined as a user-defined scalar (UDS). Upon this, ANSYS Fluent automatically generates a generic convective-diffusive transport equation for this UDS in a porous media of the following form, ∂γρφ + ∇ · (γρvφ) = ∇ · (γΓ∇φ) + γSφ ∂t
(9)
where φ is a UDS, Γ is a diffusion coefficient, and Sφ is the source term, so that our task is to make use of UDFs to transform Eq. 9 into Eq. 8. Specifically, the convective and diffusive terms are zero out, two additional source terms accounting for the two mechanisms that change the solid temperature as discussed are formulated as follows, 5 ro ∆Hi STs ,1 = i=1 i (10) γ(r) (1 − γ(r))hf s (Tf − Ts ) STs ,2 = , (11) γ(r) and the transient term in Eq. 9 is replaced by that of Eq. 8. As a result of the above analysis, the 2D heterogeneous CFD model used to describe the behavior of a catalytic process taking place in the FBR, i.e., the continuity, momentum, energy and species material balances, are given by the following equations: ∂ (γρf ) + ∇ · (γρf v) = 0 ∂t ∂ (γρf v) + ∇ · (γρf vv) = ∂t − γ∇P + ∇ · (γτ ) − γ
(12a) (12b) 2
γµ γ C2 ρf |v|v v+ αb 2
∂ (γρf CP,f Tf ) + ∇ · (γv(ρf CP,f Tf + P )) = (12c) ∂t ∇ · γkf ∇T − hi J i + τ · v − γSTs ,2 i
∂ ((1 − γ)ρs CP,s Ts ) = γSTs ,1 + γSTs ,2 ∂t ∂ (γρf Yi ) + ∇ · (γρf vYi ) = −∇ · (γJ i ) + γSY,i ∂t
(12d) (12e)
123
boundary conditions: z = 0; Yi = Yi,0 , Tf = Tf,0 ∂Tf ∂γρf Yi = 0; =0 r = 0; ∂r ∂r ∂γρf Yi ∂Tf = 0; − λr = hf j (Tf − Tj ) r = R; ∂r ∂r
123
(13a) (13b) (13c)
where v (m s−1 ), P (kPa), ρf (kg m−3 ), CP,f (J kg−1 K−1 ), kf (W m−1 K−1 ) and µ (m2 s−1 are the velocity, static pressure, density, specific constant-pressure heat capacity, thermal conductivity and molecular viscosity of the bulk flow, respectively, ∇ ( i hi J i ) represents the effect of enthalpy transport due to species diffusion, Yi and J i (kg m−2 s−1 ) are the mass fraction and mass diffusion flux, respectively, Tj is the jacket temperature profile, hf j the lumped heat transfer coefficient between the bulk flow and the cooling jacket and τ is the stress tensor. 4. FEEDBACK CONTROL DESIGN In this section, the feedback control system is developed for the phthalic anhydride synthesis reactor to drive the to its setoutlet concentration of phthalic anhydride xoutlet D point xset D and maintain the fluid temperature xT inside the reactor below the maximum allowable temperature xmax . Specifically, we derive a recurrent neural network T model capturing the relationship between fluid temperature xT , species concentration xc in the reactor and the cooling jacket temperature Tji , i = 1, 2, 3, 4. The RNN models are obtained based on training, validation and testing dataset generated from a large number of open-loop CFD simulations with various step changes of Tji , i = 1, 2, 3, 4. Subsequently, MPC is developed using the RNN model fR (·) that will be discussed in the following section as the prediction model with constraints on both jacket temperature and fluid temperature to calculate optimal control actions (i.e., cooling jacket temperature Tji , i = 1, 2, 3, 4) to regulate xoutlet . Finally, the structure D of the entire closed-loop system that integrates MPC solver and CFD dynamic simulation is presented. 4.1 Open-loop Dynamics To build a recurrent neural network model to represent the CFD model of a fixed-bed catalytic reactor for a phthalic anhydride synthesis, open-loop simulations are first conducted to generate the dataset that captures the system dynamics within a wide range of operating conditions. Specifically, by manipulating cooling jacket temperatures Tji , i = 1, 2, 3, 4 of four zones, the dataset is created for fluid temperature xjT and species concentration xjc in the reactor, where j = 1, ..., 8 represents 8 nonequidistant surfaces set in ANSYS Fluent to monitor both xc and xT along the length of the reactor (i.e., 5 of them are set within the first meter of reactor since the maximum fluid temperature occurs in this region). The open-loop CFD simulations are repeated under various step-changes on jacket temperatures to drive xoutlet to D different steady-states of xoutlet under a disturbance-free D environment. The dataset of fluid temperature and species concentration profiles in the reactor are obtained with the initial steady-state as follows: the outlet PA concentration
2019 IFAC CPDE-CDPS 124 Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
= 7.550 × 10−3 and jacket temperatures Tj1 = xoutlet D 640 K, Tj2 = 640 K, Tj3 = 640 K and Tj4 = 640 K. The CFD model is simulated with a sampling period ∆CF D = 0.1 s which represents the integration step size for simulations of the CFD model. Subsequently, timeseries data from open-loop simulations are separated into a large number of time-series data samples with shorter period via a sliding window such that the dimension of time-series data is reduced and more data is generated for the development of an RNN model. Additionally, data preprocessing is employed to remove data samples with repeated or similar dynamic behavior to avoid over-fitting during the training process of RNN. 4.2 Ensemble Regression Modeling In this study, a recurrent neural network (RNN) learning algorithm (Lipton et al. (2015)) and a 10-fold cross validation are used to construct a homogeneous ensemble regression that uses the feed operating parameters, the jacket temperature profile and the real-time state measurements at eight discrete locations to predict the nonlinear distributed dynamics of the FBR. The motivation for using the ensemble regression is four-fold. First, the RNN learning algorithm, which is also known as the backpropagation through time (BPTT), can easily derive RNNs that memorize the training data well, however, without the generalization capability, and therefore, by exposing the learning algorithm to different subsets of the training data, generated by 10-fold cross validation, effectively prevents it from overfitting the RNNs. Second, the learning algorithm is known to be non-convex and is an infamous example of NP-hard problems, and hence, by using different starting initial weight matrices, the learning algorithm might be able to avoid getting trapped in a local minimum and arrive at the set of weights that allows the RNN to accurately approximate the latent function that transforms input sequences to corresponding output sequences. Third, an ensemble regression is known to often outperform its individual constituents and simultaneously accounts for uncertainty in model selection. Fourth, the proposed approach creates a unique opportunity for us to demonstrate our general framework for linking ANSYS Fluent to several powerful APIs that are typically incompatible. For those reasons, the ensemble regression structured to consist of eight sub-ensembles, and each of which is a 10-fold cross-validated committee of RNNs as depicted in Fig. 2 and 3 and is derived to predict the nonlinear distributed dynamics of the FBR. Specifically, Fig. 2 demonstrates the RNN structure where the first column is the input layer, the second and the third columns are hidden layers and the last column is the output layer. In the input layer, Xkij , k = 1, ..., 6 represents the realtime state measurements at eight discrete locations i (i = 1, ..., 8) for j th (j = 1, ..., 10) RNN model, Xkij , k = 7, ..., 12 represents the estimated states from previous location at current time, and Xkij , k = 13, ..., 16 represents jacket temperatures for four zones. In the hidden layer, the nodes calculate the weighted sum of input nodes with nonlinear activation functions and also introduced feedback loops to maintain information in ‘memory’ over time. At the end, the nodes in the output layer ykij , k = 1, ..., 6 report 124
the estimation of future states, and will be used as the predicted states in the MPC in the next section.
Fig. 2. The structure of the recurrent neural network. In this study, each sub-ensemble of the ensemble regression is designed to predict the trajectories of the state variables at a fixed location over a time interval equal to the sampling rate of the centralized controller, which can be an integer multiple of the CFD integration time step (in this work it is set as 50 CFD integration time steps). Specifically, as shown in Fig. 3, the sub-ensemble uses the on-line measurements of the state variables at the designated location along the FBR, the estimated trajectories of the state variables at an upstream location over the elapsed time and the jacket temperature profile to calculate its predictions, which are the average trajectories of those generated by its constituent RNNs. This ensemble integration strategy appears to be rudimentary; nonetheless, the performance of the sub-ensemble measured in terms of the absolute mean error is statistically better than that of its constituent RNNs. The constituent RNNs of individual sub-ensembles are known as the 10-fold cross-validated committee; this is because that 10-fold cross validation is used to generate 10 different subsets of the training data, and each of which is subsequently used to train a RNN. Although the state-ofthe-art API, Keras, can be readily used to create the 10fold cross-validated committee, training RNNs remains a nontrivial task because the learning algorithm is known to suffer from the vanishing and exploding gradient problems. For that reason, the RNNs is designed to have two hidden recurrent layers consisting of 64 recurrent units that use the rectified linear unit function as the activation function (i.e., ReLu function f (x) = max{0, x}) and are initialized as proposed in (Le et al. (2015)). Specifically, Le et al. (2015) demonstrates that this initialization in which the recurrent and bias weight matrices are set to the identity matrix and zero, respectively, allows the trained RNNs to have the consistently comparable performance with the standard long short term memory (LSTM) networks in various applications. On a last note of this Section, we would like to detail the procedure encoded in the ensemble regression which is used to predict the nonlinear distributed dynamics of the FBR over a finite time interval in the future. Specifically, the first sub-ensemble normalizes the input signals, i.e., the feed operating parameters, the real-time state measurements at the first location and the jacket temperature profile, using the statistics of the first training set and, subsequently, transmits the normalized input signals to its cross-validated committee of RNNs. It is critical to note that all input and output signals of individual RNNs are
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
normalized based solely on the statistics of its training set to prevent data leakage from its test set. Next, the constituent RNNs simultaneously compute the trajectories of the normalized state variables over a finite time interval in the future, which are pooled to compute the respective average trajectories. Then, the first sub-ensemble generates its output signals by rescaling the average trajectories of the normalized state variables using the corresponding statistics. The output signals of the first sub-ensemble are used to infer the dynamics of the FBR at the corresponding location and as input signals to the second sub-ensemble in addition to the real-time state measurements at the second location and the jacket temperature profile. All downstream sub-ensembles execute the same procedure, and the outputs of the sub-ensembles are used to estimate the nonlinear distributed dynamics of the FBR. 4.3 Model Predictive Controller In this work, the primary objectives of the feedback control scheme is to regulate xoutlet while mitigating issues D related to catalyst deactivation and thermal runaway. Based on the RNN model derived in Section 4.2, a model predictive controller is developed to drive xoutlet to its D desired set-point xset and ensure the fluid temperature D xT does not exceed the maximum allowable hot-spot max temperature xT . We use x = [xD xT xO ]T ∈ R6 to represent the states of the weighted-averaged crosssectional mole fraction of PA xD (termed averaged PA concentration in the following text), the weighted-averaged cross-sectional fluid temperature xT (termed averaged fluid temperature in the following text) and all the other species concentrations in the reactor. Additionally, fluid temperature and jacket temperatures are bounded by := {0 < xT ≤ xmax } ∈ R and U bound := {umin ≤ xbound i T T max 4 ui ≤ ui } ∈ R (i = 1, 2, 3, 4), respectively to avoid hotspot formation. The optimization problem of MPC is of the following form: tk+N
J = min
u∈S(∆)
s.t.
tk
2 (t) − xset |˜ xoutlet D | D
(14a)
x ˜(tk+i+1 ) = fR (˜ x(tk+i ), u(tk+i )), i = 0, 1, ..., N − 1 (14b) (14c) x ˜(tk ) = x(tk ) u(tk+i ) ∈ U bound , i = 0, 1, ..., N − 1
(14d)
x ˜T (tk+i ) ∈ xbound , i = 0, 1, ..., N − 1 (14e) T where x ˜D (t), x ˜T (t) are the predicted state trajectories for the averaged PA concentration at the outlet and the maximum averaged fluid temperature, respectively. fR (·) represents the RNN model, S(∆) is the set of piecewise constant functions with period ∆, and N is the number of sampling periods in the prediction horizon. In the optimization problem of Eq. 14, the objective func2 tion of Eq. 14a is the sum of |˜ xoutlet (t) − xset D D | over the prediction horizon. The constraint of Eq. 14b is the RNN model fR (·) used to predict the evolution of states including xD (t) and xT (t). Eq. 14c takes the measurements of states at t = tk from the CFD simulation as the initial conditions for the RNN model of Eq. 14b. Eq. 14d defines the input constraints over the entire prediction horizon. Eq. 14e defines state constraints over the prediction hori125
125
zon. Additionally, we assume that the states of the closedloop system are measured at each sampling time. After the optimal solution u∗ (t) of the optimization problem of Eq. 14 is obtained, only the first control action of u∗ (t) is sent to the control actuators to be applied over the next sampling period. Then, at the next instance of time tk+1 := tk + ∆, the optimization problem is solved again, and the horizon is rolled forward one sampling time. Under the MPC of Eq. 14, the closed-loop state xoutlet can D be ultimately driven to its desired set-point xset D in the absence of model mismatch or disturbances (e.g., reactor feed disturbances) if there is sufficient energy to satisfy the constraints of Eq. 14d and Eq. 14e. To account for model mismatch and eliminate the offset between the final steady-state value of xoutlet and xset D D , an additional deviation term e(tk−1 ) = x ˜(tk−1 )−x(tk−1 ) which represents the deviation between the actual state and the predicted state in the last step t = tk−1 is applied in the predictive models of Eq. 14b such that the model mismatch is accounted for in the predictive models (i.e., x ˜(tk+1 ) = fR (˜ x(tk ), u(tk )) + e(tk−1 )). Alternatively, offset can also be eliminated by integrating an integral feedback control scheme with MPC such that the overall control action can be obtained as: uM P C+I (tk ) = uM P C (tk ) + uI (tk ) (Wu et al. (2017)). In our work, the MPC is developed with RNN models via the scientific computing package SciPy in Python. The integration of the optimization problem of MPC and RNN models is demonstrated by the algorithm shown in Fig. 4. Specifically, at t = tk the RNN model and the state measurements are first imported to MPC. Subsequently, the MPC control actions are calculated by solving the optimization problem of Eq. 14. It is noted that in Eq. 14b, the RNN model fR (·) is utilized to predict future states (i.e., state profiles over the prediction horizon t = [tk , tk+N )) by taking state measurement x(tk ) and control actions u(t), t = [tk , tk+N ) as the inputs. In Fig. 4, it is further demonstrated that following the RNN structure shown in (i) Fig. 3, each sub-RNN model (i.e., the RNN model fR (·), i = 1, ..., 8 for each surface in the reactor) not only utilizes state measurements and control actions, but also takes the predicted outputs x(tk+1 ) from the previous RNN model (i−1) fR (·) (i.e., at (i − 1)th surface) as the inputs to predict the system dynamics at ith surface. 4.4 Integrating MPC and Dynamic CFD Simulation In this section, an integrated framework of MPC and dynamic CFD simulation is developed via message passage interface (MPI) to implement MPC within CFD dynamic simulations. Since UDF supports C language only, the co-compilation of UDF and MPC solver depends on the specific programming language used. For example, the optimization problem of MPC is solved in Python environment in our work, and therefore, the MPC solver and UDF are compiled in Python and C environments, respectively, and a communication bridge between the MPC and ANSYS Fluent UDF needs to be developed through a lock synchronization function embedded in the UDF. As shown in Fig. 5, the lock synchronization mechanism is utilized to exchange data between CFD simulation and MPC solver. Specifically, at the end of each sampling period of CFD dynamic simulation, the measurements of
2019 IFAC CPDE-CDPS 126 Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
Fig. 3. The schematic of ensemble learning with RNN models.
Fig. 4. The algorithm of integrating MPC with RNN models where N is the length of prediction horizon and i is the index of RNN models.
process states are first obtained via UDF (e.g., functions C T, F YI are utilized to derive fluid temperature and PA concentration in ANSYS Fluent via user-defined functions (UDF)) and then sent to the MPC solver. Meanwhile, the CFD simulation is forced to wait for the control actions from MPC. Then, the MPC solver is invoked to solve the optimization problem to derive control actions u(t), for the next sampling period t ∈ [tk , tk+1 ) based on received process state measurements at t = tk . Finally, the CFD simulation receives u(tk ) and continues dynamic simulation until the next sampling time step. The above procedure is repeated with new measurements at the end of each sampling period. Therefore, under the integration of MPC and dynamic CFD simulation, jacket temperatures (i.e., manipulated inputs calculated by MPC) can be adjusted at each sampling step such that xoutlet can D be ultimately driven to its set-point in the absence of model mismatch or disturbances. Additionally, it should be noted that actual process states are measured in CFD simulations while normalized process states are utilized in RNN models and MPC solvers. Therefore, data processing function is also embedded in MPC solver to convert it to normal value before sending it back to CFD simulator. 126
Fig. 5. Basic structure of the proposed integrated MPC and CFD dynamic simulation via a lock synchronization mechanism. 5. CLOSED-LOOP SIMULATION RESULTS In this section, the closed-loop performance of a fixedbed catalytic reactor for the phthalic anhydride synthesis under MPC is investigated. All simulation settings of the closed-loop simulations are the same as those used in the previously studied open-loop system with ∆M P C = 5.0 s and ∆CF D = 0.1 s, respectively. We demonstrate the case of set-point tracking under a disturbance-free environment. The initial steady-state is at xoutlet = 7.550 × 10−3 , D and the jacket temperatures are Tj1 = 640 K, Tj2 = 640 K, Tj3 = 640 K and Tj4 = 640 K, respectively. The new set-point for averaged PA concentration at the outlet is −3 xset . It is shown in Fig. 6 that averaged D = 7.573 × 10 PA concentration at the outlet xoutlet is successfully driven D −3 to xset under MPC, and it takes around D = 7.573 × 10 170 s for xoutlet to reach its final value. The dynamic perD formance of the closed-loop system under MPC improves significantly compared to open-loop control, while it takes −3 around 300 s for the system to settle to xset . D = 7.573×10 The corresponding manipulated input profiles are given in Fig. 7, where it is demonstrated that the jacket temperatures of the four zones stay at high temperature (around
2019 IFAC CPDE-CDPS Oaxaca, Mexico, May 20-24, 2019
Zhe Wu et al. / IFAC PapersOnLine 52-2 (2019) 120–127
quickly, and 685 K) for the first 20 s to increase xoutlet D then the jacket temperatures decrease to their steady-state values such that xoutlet can approach its desired set-point D set xset D and ultimately settle down to xD . Additionally, it is demonstrated in Fig. 8 that the fluid temperature profile is maintained below the maximum allowable temperature 773 K under the constraints of Eq. 14e, which avoids the formation of hot-spots that could deactivate the catalyst. 7.58
×10
-3
Closed-loop system under MPC Desired xset D
7.575
xoutlet D
7.57 7.565 7.56 7.555 7.55
0
50
100
150
200
250
300
350
Time (s)
Fig. 6. The outlet PA concentration xoutlet profile of the D closed-loop system under MPC with xset D = 7.573 × 10−3 .
highly exothermic process was utilized to demonstrate the effectiveness of the proposed approach, under which the product yield can be maximized while suppressing the hot-spot temperature to avoid catalyst deactivation and thermal runaway. Initially, a high-fidelity two-dimensional axisymmetric heterogeneous CFD model for an industrialscale FBR was developed in ANSYS Fluent. Then, the open-loop simulation of CFD model was conducted to create a database for a wide operating regime of the FBR, from which recurrent neural network and ensemble learning were used to derive a homogeneous ensemble regression model using the application program interface (API) Keras. Subsequently, based on the RNN models derived from ensemble regression, an MPC control scheme was developed to drive the process outputs to the desired set-points and suppress the magnitude of the hot-spot temperature to avoid catalyst deactivation. The CFD model, the ensemble regression and the MPC were integrated to create a closed-loop system by linking ANSYS Fluent to MPC optimizer in Python via a message-passing interface (MPI) with lock synchronization mechanism. Finally, the closed-loop simulation results demonstrated the robustness and effectiveness of the proposed approach. 7. ACKNOWLEDGMENTS
Jacket temperature: Zone 1
680
127
Jacket temperature (K)
660
Financial support from the National Science Foundation and the Department of Energy is gratefully acknowledged.
640 680
Jacket temperature: Zone 2
660 640 680
REFERENCES
Jacket temperature: Zone 3
660 640
680 Jacket temperature: Zone 4
660 640
0
50
100
150
200
250
300
350
Time (s)
Average cross-sectional temperature of fluid (K)
Fig. 7. The jacket temperature profiles of the closed-loop −3 system under MPC with xset . D = 7.573 × 10 720 700 680 660 640 300 200
Time (s)
100 0
0
1
2
3
4
Distance along the PA reactor (m)
Fig. 8. The fluid temperature profile of the closed-loop −3 system under MPC with xset . D = 7.573 × 10 6. CONCLUSION In this work, a general framework for linking of a CFD solver, ANSYS Fluent and other computing platforms using the lock synchronization mechanism was proposed to extend the utilities of CFD solvers from strictly modeling and design to also control and optimization applications. Specifically, a fixed-bed reactor (FBR) with a 127
Anastasov, A.I. (2003). An investigation of the kinetic parameters of the o-xylene oxidation process carried out in a fixed bed of high-productive vanadia–titania catalyst. Chemical Engineering Science, 58, 89–98. De Klerk, A. (2003). Voidage variation in packed beds at small column to particle diameter ratio. AIChE journal, 49, 2022–2029. Le, Q.V., Jaitly, N., and Hinton, G.E. (2015). A simple way to initialize recurrent networks of rectified linear units. arXiv preprint arXiv:1504.00941. Lipton, Z.C., Berkowitz, J., and Elkan, C. (2015). A critical review of recurrent neural networks for sequence learning. arXiv preprint arXiv:1506.00019. Tran, A., Pont, M., Crose, M., and Christofides, P.D. (2018). Real-time furnace balancing of steam methane reforming furnaces. Chemical Engineering Research and Design, 134, 238–256. Vaquerizo, L. and Cocero, M. (2018). CFD–Aspen Plus interconnection method. Improving thermodynamic modeling in computational fluid dynamic simulations. Computers & Chemical Engineering, 113, 152–161. Wu, W. and Huang, M.Y. (2003). Nonlinear inferential control for an exothermic packed-bed reactor: Geometric approaches. Chemical Engineering Science, 58, 2023– 2034. Wu, Z., Aguirre, A., Tran, A., Durand, H., Ni, D., and Christofides, P.D. (2017). Model predictive control of a steam methane reforming reactor described by a computational fluid dynamics model. Industrial & Engineering Chemistry Research, 56, 6002–6011.