Comprehensive Framework for Controlling Nonlinear Multi-Species Water Quality Dynamics

1 Comprehensi v e Frame work for Controlling Nonlinear Multi-species W ater Quality Dynamics Salma M. Elsherif † , ¶ , Ahmad F . T aha † , ∗∗ , Ahmed A. Abokifa ⋄ , and Lina Sela ♠ Abstract T racing disinfectant (e.g., chlorine) and contaminants ev olution in water networks requires the solution of 1-D adv ection- reaction (AR) partial differential equations (PDEs). With the absence of analytical solutions in many scenarios, numerical solutions require high-resolution time- and space-discretizations resulting in large model dimensions. This adds complexity to the water quality control problem. In addition, considering multi-species water quality dynamics rather than the single-species dynamics produces a more accurate description of the reaction dynamics under abnormal hazardous conditions (e.g., contamination events). Y et, these dynamics introduce nonlinear reaction formulation to the model. T o that end, solving nonlinear 1-D AR PDEs in real time is critical in achieving monitoring and control goals for various scaled networks with a high computational burden. In this work, we propose a nov el comprehensive framew ork to ov ercome the large-dimensionality issue by introducing different approaches for applying model order reduction (MOR) algorithms to the nonlinear system follo wed by applying a real-time water quality regulation algorithm that is based on an advanced model to maintain desirable disinfectant le vels in water networks under multi-species dynamics. The performance of this framew ork is validated using rigorous numerical case studies under a wide range of scenarios demonstrating the challenges associated with regulating water quality under such conditions. Index T erms Multi-species water quality dynamics, water quality regulation and control, model predictive control, McCormick relaxation, linear/nonlinear model order reduction. I . I N T RO D U C T I O N A N D L I T E R A T U R E R E V I E W W ater quality dynamics are widely modeled by the one-dimension advection-reaction (1-D AR) partial dif ferential equations (PDEs). These AR PDEs allow the tracing of the disinfectant and other chemical substances’ e volution throughout the components of water distribution networks (WDNs). In most cases, analytical solutions are non-existent to solve these PDEs network-wide. Nonetheless, PDEs can be solved using numerical techniques, yet, they require high-resolution time- and space- discretization. This results in high-dimension models that add computational burden to the problem of regulating water quality in drinking networks. That leads to physical-dri ven models that are intractable when considering constrained control and water quality (WQ) re gulation algorithms. Moreov er , in water quality simulations the most widely used decay and reaction model is the single-species model. In this model, disinfectant (i.e., chlorine) is assumed to decay at a constant rate that only accounts for puriﬁed water contamination lev els. Y et, contamination sources v ary from microbial, non-microbial components in the bulk ﬂow , attached to the pipe walls, or contamination ev ents that get intruded into the system [ 1 ]. This driv es the need for a more accurate representation of these scenarios which can be achiev ed by the multi-species reaction dynamics. The multi-species dynamics enable the model to simulate chlorine ev olution with the existence of another reactive component in the system. This representation duplicates the number of variables to be traced network-component-wide while unfortunately adding comple xity to the model by introducing nonlinear reaction dynamics. T o this end, model order reduction (MOR) is an essential step to mov e forward in achie ving a compact formulation of the multi-species water quality dynamics to be integrated into a model-based control framework. MOR techniques transform the full-or der model (FOM) to a r educed-or der model (R OM) in a way that preserves the structure, properties, and the closed-form representation of the FOM while achie ving the pre-speciﬁed lev el of accuracy and reducing computational time. Eventually , the goal is to control chlorine injections dosed by rechlorination stations to maintain residual lev els that meet water quality standards. That can be achieved by applying an effecti ve control algorithm on the deriv ed R OM. Our group has been interested in v arious dimensions of this research area. A summary of our work and the prior literature is gi ven ne xt. † Department of Civil and En vironmental Engineering, V anderbilt University , Nashville, TN, USA. Emails: salma.m.elsherif@vanderbilt.edu, ah- mad.taha@vanderbilt.edu ¶ Secondary: Department of Irrigation and Hydraulics Engineering, Faculty of Engineering, Cairo University . ⋄ Department of Ci vil, Materials, and En vironmental Engineering, The Univ ersity of Illinois Chicago. Email: abokifa@uic.edu ♠ Department of Ci vil, Architecture, and En vironmental Engineering, The Univ ersity of T exas at Austin. Email: linasela@utexas.edu ∗∗ Corresponding author . This work is partially supported by National Science Foundation under grants 1728629, 2015603, 2015671, 2151392, and 2015658. 2 Fig. 1: Conceptual framework of the paper . A. Liter ature Re view Hereinafter , we survey the literature on the topics of MOR for dynamic systems, in general, and water systems, in particular , and water quality regulation and control while highlighting the gaps and dra wbacks moti vating this paper’ s contributions. MOR for Dynamic Systems: Se veral studies hav e proposed and implemented dif ferent MOR algorithms in v arious disciplines (e.g., electromagnetics, electro-mechanics, structural and ﬂuid dynamics) where the lar ge-dimensionality issue is faced [ 2 ]–[ 5 ]. Most of these studies ha ve applied either singular v alue decomposition (SVD) approach [ 6 ], [ 7 ] or Krylo v subspace methods [ 8 ], [ 9 ]. Whilst, combining SVD and Krylov methods is in vestigated and implemented in [ 10 ]. In infrastructure networks preserving the system’ s properties including stability , controllability , and observ ability is a major concern while applying MOR with the aim of applying a post-reduction control. Ne vertheless, Krylo v methods do not preserve such properties which limit their suitability to be used in our study [ 3 ]. Sev eral SVD-based model reduction methods have been proposed for linear systems and more realizations and extensions hav e been inv estigated and integrated to tackle the reduction of nonlinear systems; a revie w of linear and nonlinear models order reduction can be found in [ 11 ], [ 12 ]. The Balanced T runcation (BT) method [ 2 ] is b uilt based on both the controllability and observability Gramians for stable, linear systems. The study by [ 13 ] extends BT to be applied to nonlinear systems, while the authors in [ 14 ], [ 15 ] build extensions for unstable systems. Howe ver , BT becomes computationally intractable for large- scale systems. Nev ertheless, the famous widely used Proper Orthogonal Decomposition (POD) in ﬂuid dynamics community [ 16 ] is considered tractable at the expense of accuracy compared to BT . Y et, in some cases where relativ ely lower accuracy is acceptable, POD may result in an unstable system even near stable equilibrium points depending on the actual formulation of the full-order model. Therefore, the methods that balance between the BT and POD methods hav e been propounded to integrate the adv ance of both methods into one. For example, the authors in [ 17 ] hav e proposed a balanced method but it has failed to successfully reduce models when the number of outputs of the system is large. Con versely , the balanced POD (BPOD) [ 6 ] is tractable with an overall computation time similar to POD, but it computes adjoint snapshots to combine and balance controllability and observability similarly to BT , which is not raised in POD. Furthermore, POD can be extended to reduce the order of nonlinear systems by approximately projecting the nonlinearity term in the system to a subspace of the dynamics [ 3 ], [ 18 ]. Therefore, the nonlinear term is ev aluated separately and approximated at on l y a small set of interpolation points (hyper-parameter) using a combination of projection and interpolation methods such as the discrete empirical interpolation method (DEIM) [ 18 ], the Gappy POD method [ 19 ], [ 20 ], and the Gauss-Newton with approximated tensors (GN A T) method [ 21 ], refer to the revie w paper [ 22 ] for details. MOR for W ater Systems: W ater systems model order reduction has been broadly inv estigated for network hydraulics over the past decades with a limited number of studies looking into MOR for w ater quality dynamics. These studies adopt different approaches to reduce the hydraulic model dimension by applying methods varying between performing nodal Gaussian elimination [ 23 ], Gaussian elimination on the linearized form of the model and recov ering the nonlinearity of the system as a post-reduction step [ 24 ], genetic algorithm [ 25 ], and system aggregation [ 26 ]. Perelman and Ostfeld [ 27 ] consider a coupled model that combines both hydraulics and water quality dynamics of the network and apply systems aggregation. Lately , two studies have applied different approaches to cover the MOR for water quality dynamics. Authors in [ 28 ] have proposed reducing the order of the water quality model by formulating a bi-linear spatially-discretized but a temporally- continuous representation of the dynamics. This formulation augments the input vectors in a way that preserves the system’ s stability . The induced error between the actual and reduced-order models is minimized by the reduction of the H 2 -norm. In this study , water dynamics transport and reaction are simulated using the Advection-Dif fusion-Reaction Partial Differential Equations that include the dif fusion term in comparison to our work that neglects its ef fect. Nonetheless, studies [ 29 ], [ 30 ] state that diffusion is dominated in network branches with signiﬁcantly low velocities. T o that end, it is an acceptable assumption to neglect the diffusion effect in networks with limited dead-end branches, higher velocities, and changing demands. On the other hand, augmenting and transferring the model into nonlinear formulation results in a more comple x one when considering the multi-species nonlinear water quality dynamics and does not preserve the stability of the system. Secondly , study [ 31 ] applies dif ferent SVD-based projection algorithms to reach a reduced-order w ater quality model including BT , POD, and BPOD in addition to preserving the stability of the BPOD method. Results have prov en that the BPOD method is more usable while being computationally tractable and robust for zero and non-zero initial conditions. 3 Howe ver , their model only includes single-species linear reaction dynamics where chlorine is assumed to be decaying at a constant rate resulting in a linear state-space formulation. Therefore, our work allows ﬁlling the gap in applying MOR for multi-species nonlinear dynamics. Moreover , in their model the explicit central Lax-W endroff discretization scheme is used. Howe ver , the Upwind schemes giv e a more accurate physical description of the advection-reaction problem. In our study , we apply the explicit and implicit Upwind discretization schemes while highlighting the differences and the level of dif ﬁculty . Notice that, on the contrary to studies [ 32 ]–[ 34 ] where MOR is performed for compositional simulation, the authors in [ 31 ] state that it is considered as a pre-step to apply an efﬁcient control algorithm which also applies to our work in this paper . W ater Quality Control: The topic of controlling chlorine has been cov ered in se veral studies with various algorithms, objectiv es, and constraints [ 35 ]–[ 37 ]. Objectives vary between minimizing the cost of injecting chlorine into the system, maintaining minimal deviations from chlorine setpoint concentrations, minimizing the formation of the excess DBPs, and minimizing computational time [ 38 ]. The problem formulation is either a single-objective optimization problem or a multi- objectiv e one with more of the aforementioned objectiv es. Howe ver , such studies do not build a closed-form representation of all inputs, states, and outputs that updates ev ery speciﬁed time-step ov er the simulation period and allows network-wide control. Whereas, studies [ 31 ], [ 39 ] apply Model Predictiv e Control (MPC) on the full-order and reduced-order single-species models in both studies with no clear explanation/extension for scenarios where multi-species dynamics take place. Our prior W ork: W e have been focusing on tackling and covering the water quality modeling and control in WDNs. First, the problem of modeling and controlling single-species water quality dynamics is thoroughly in vestigated in [ 39 ], followed by reducing this model’ s order and verifying the validity of controlling the reduced order model in [ 31 ]. Moreover , as a ﬁrst state-of-the-art attempt, study [ 40 ] has identiﬁed single-species water quality models using only input-output experimental data and, accordingly , data-dri ven system identiﬁcation algorithms. Lastly , a survey study on how to accurately simulate multi- species water quality dynamics has been conducted in [ 41 ]. This study has built a closed-form, network- and control-theoretic representation of all system inputs, variables, and output measurements under such dynamics that giv e a more realistic WQ formulation. The performance of this formulation has been validated using the widely-used simulation tool, EP ANET and its multi-species water quality simulation extension, EP ANET -MSX [ 42 ], [ 43 ]. Howe ver , controlling chlorine under multi-species dynamics, based on a control-theoretic explicit model, is to the authors’ knowledge has not been inv estigated—a gap that is ﬁlled in this paper . B. P aper Contributions This paper’ s major objective is to in vestigate the implementation and complexity of regulating and controlling chlorine lev els under multi-species water quality . The detailed paper contributions are: • Construct and propose a comprehensive framework to ov ercome the large-dimensionality issue associated with discretizing the 1-D AR PDEs and the complexity associated with the nonlinearity of the multi-species water quality dynamics. Different paths can be taken, starting by linearizing the system and applying MOR for linear systems (MOR-LS). Another path is to consider the nonlinear MOR (MOR-NLS) algorithm on the original FOM. • Utilizing the reduced-order models in an MPC algorithm. W e apply it to the formulated R OMs and compare them to each other and to the original FOMs. Also, we compare it with basic scenarios with single-species dynamics demonstrating the challenges associated with controlling chlorine lev els under multi-species water quality dynamics. • Position the framework in a generalized scalable form in the sense that simpliﬁcations are included to consider single-species water dynamics and differentiations are suggested to consider chlorine linear/nonlinear decay and reaction models that have been developed in the literature to simulate various ev ents/scenarios. • V alidate the performance of the framework using thorough numerical case studies to test accuracy , computational burden, and robustness to the system hydraulics changes. Our proposed frame work is illustrated in Fig. 1 . As shown, different approaches can be followed to formulate a reduced- order model to be controlled for the multi-species water quality model. Each step to be taken and each path to be chosen are explained in the following sections of the paper . The paper’ s sections are organized as follows, Section ” State-space Multi- species W ater Quality Model ” provides the formulation of the state-space representation of the multi-species water quality model (MS-WQM). This formulation is based on the transport and reaction model in pipes, mass balance for the other network components, and the multi-species dynamics expression. Section ” Model Order Reduction and T ransformation of MS-WQM ” provides full descriptions of the methods used in our framew ork to reach a compact reduced-order model. Section ” Real-T ime Regulation of MS-WQM via Model Predictiv e Control and McCormick Relaxations ” introduces the control problem and its implementation on the linear and nonlinear ROM. Section ” Case Studies ” showcases the framew ork performance on different networks under a wide range of scenarios. Section ” Conclusion, Paper’ s Limitations, and Recommendations for Future W ork ” comes last. I I . S TA T E - S PAC E M U LT I - S P E C I E S W A T E R Q UA L I T Y M O D E L W e model WDN by a directed graph G = ( N , L ) . The set N deﬁnes the nodes and is partitioned as N = J ∪ T ∪ R where sets J , T , and R are collections of junctions, tanks, and reservoirs. Let L ⊆ N × N be the set of links, and deﬁne 4 the partition L = P ∪ M ∪ V , where sets P , M , and V represent the collection of pipes, pumps, and valv es. T otal number of states is n x = n L + n N , where n L and n N are numbers of links and nodes. The number of reservoirs, junctions, tanks, pumps, valv es, and pipes are n R , n J , n TK , n M , n V , and n P . Each pipe i with length L i is spatially discretized and split into s L i segments. Hence, number of links is expressed as n L = n M + n V + P n P i =1 s L i while n N = n R + n J + n TK is the number of nodes. In this paper, the state-space representation is formulated for multi-species dynamics with two chemicals: chlorine and a ﬁctitious reactant. The system representation of the two-species which is able to capture chemicals evolution, booster stations injections, and sensors measurements, is expressed by an NDE as follows   E 11 ( t ) 0 0 E 22 ( t )   | {z } E ( t )   x 1 ( t + ∆ t ) x 2 ( t + ∆ t )   | {z } x ( t +∆ t ) =   A 11 ( t ) 0 0 A 22 ( t )   | {z } A ( t )   x 1 ( t ) x 2 ( t )   | {z } x ( t ) +   B 11 ( t ) 0 0 B 22 ( t )   | {z } B ( t )   u 1 ( t ) u 2 ( t )   | {z } u ( t ) + f ( x 1 , x 2 , t ) , (1a)   y 1 ( t ) y 2 ( t )   | {z } y ( t ) =   C 11 ( t ) 0 0 C 22 ( t )   | {z } C ( t )   x 1 ( t ) x 2 ( t )   | {z } x ( t ) +   D 11 ( t ) 0 0 D 22 ( t )   | {z } D ( t )   u 1 ( t ) u 2 ( t )   | {z } u ( t ) (1b) where v ariable t represents speciﬁc time in a simulation period [0 , T s ] ; ∆ t is the time-step or sampling time; vectors x 1 ( t ) and x 2 ( t ) ∈ R n x depict the concentrations of chlorine and the other ﬁctitious reactant (two species model) in the entire network; vector u 1 ( t ) ∈ R n u 1 represents the dosages of injected chlorine; vector u 2 ( t ) ∈ R n u 2 accounts for planned or unplanned injection of the ﬁctitious component; vector f ( x 1 , x 2 , t ) encapsulates the nonlinear part of the equations representing the mutual nonlinear reaction between the two chemicals; vector y 1 ( t ) ∈ R n y 1 denotes the sensor measurements of chlorine concentrations at speciﬁc locations in the network while y 2 ( t ) ∈ R n y 2 captures the ﬁctitious reactant measurements by sensors in the network if they exist. The state-space matrices { E , A , B , C , D } • are all time-varying matrices that depend on the network topology and parameters, hydraulic parameters, decay rate coefﬁcients for the disinfectant, and booster stations and sensors locations. It is customary to assume that these matrices evolv e at a slower pace than the states x ( t ) and control inputs u ( t ) . On another note, matrices E 11 , E 22 are changing ev ery hydraulic time-step allowing them to be represented at time t not t + ∆ t of the water quality simulation horizon. The concentration ev olution throughout network components is covered by the conservation of mass law , transport, decay , and reaction models of the substances. A full description of how the models are deriv ed for each type of the components is provided in [ 41 ]. Howev er, for the reader to be able to follow up with the dev elopments of this paper, some material from [ 41 ] need to be reproduced and altered. W e list a brief overvie w of the governing equations formulating our model and its state-space representation in the following sections. A. T ransport and Reaction in Pipes Conservation of mass during transport and reaction in pipes is simulated by the one-dimension advection-reaction (1-D AR) partial differential equation, which for Pipe i is expressed as ∂ t c P i = − v i ( t ) ∂ t c P i + R P MS ( c P i ( x, t )) , (2) where c P i ( x, t ) is concentration in pipe at location x and time t ; v i ( t ) is the mean ﬂo w v elocity; and R P MS ( c P i ( x, t )) is the multi-species reaction rate in pipes expression (more explanation is given in Section ” Multi-species Reaction and Decay Model ”). Eq. ( 2 ) is discretized over a ﬁxed spatio-temporal grid, that for a Pipe i with length L i is split into a number of segments s i = j L i v i ( t )∆ t k of length ∆ x i = L i s i . In the considered 1-D AR model, the main two processes are the adv ection where the concentration at a certain location and time is affected by upstream concentrations, and reaction where chemicals decay and/or mutually react. That being said, Upwind discretization schemes are more descriptive to the actual physical process considered among other schemes [ 44 ]. Applying the Eulerian Finite-Difference based Implicit Upwind scheme on the multi- species water quality dynamics representation adapted in this paper has shown reliable results that trace chemicals contractions within different networks with various scales, according to [ 41 ]. In this paper we consider both Explicit and Implicit Upwind schemes to in vestigate their performance from a control-theoretic perspectiv e (See Fig. 2 ). 1) Explicit Upwind Scheme: F or segment s of Pipe i except for the ﬁrst segment, the concentration is calculated as c P i ( s, t + ∆ t ) = (1 − λ i ( t )) c P i ( s, t ) + λ i ( t ) c P i ( s − 1 , t ) + R P MS ( c P i ( s, t ))∆ t, (3) 5 Fig. 2: Implicit and Explicit Upwind discretization schemes for Pipe i connecting Junctions 1 and 2. Each scheme calculates concentration c P i ( s, t + ∆ t ) at segment s (colored in maroon) depending on concentrations at the segments/nodes included in its frame. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.) where λ i ( t ) = v i ( t )∆ t ∆ x i is the Courant number and according to Courant-Friedrichs-Lewy condition (CFL), Courant number (CN) is maintained to be in the range of 0 < λ i ( t ) ≤ 1 so that the scheme is stable. Moreover , the concentrations in the ﬁrst segment is expressed as in ( 4 ) assuming that the connecting upstream node is Junction j . c P i (1 , t + ∆ t ) = (1 − λ i ( t )) c P i (1 , t ) + λ i ( t ) c J j ( t ) + R P MS ( c P i ( s, t ))∆ t, (4) 2) Implicit Upwind Scheme: The difference is that the concentration at the upstream segment/node is taken at the current time-step we are calculating at. That is, Equ. ( 5a ) calculates the concentration for segment s of Pipe i . As well, the concentration of the ﬁrst segment with Junction j as the upstream node is expressed in Equ. ( 5b ). (1 + λ i ( t )) c P i ( s, t + ∆ t ) − λ i ( t ) c P i ( s − 1 , t + ∆ t ) = c P i ( s, t ) + R P MS ( c P i ( s, t ))∆ t, (5a) (1 + λ i ( t )) c P i (1 , t + ∆ t ) − λ i ( t ) c J 1 ( t + ∆ t ) = c P i (1 , t ) + R P MS ( c P i (1 , t ))∆ t. (5b) B. Mass Balance at Network Components For components other than pipes, conservation of mass is applied to formulate expressions for concentrations calculation. 1) Mass Balance at Reservoir: Reservoirs are assumed to have constant concentrations. For each Reservoir i concentration is expressed as c R i ( t + ∆ t ) = c R i ( t ) . 2) Mass Balance at Pumps and V alves: The model deals with pumps and valves as transmission links with concentration equals the concentration of the node upstream them. That being said, for Pump i or V alve j installed after Reservoir k (as an example), concentrations are expressed as c M i ( t + ∆ t ) = c R k ( t + ∆ t ) , and c V j ( t + ∆ t ) = c R k ( t + ∆ t ) . 3) Mass Balance at J unctions: Chemicals are assumed to ha ve complete and instantaneous mixing in junctions with no storage time. Thus, chemical concentration at each Junction i is expressed as c J i ( t ) = P j ∈ L in q j in ( t ) c j in ( t ) + q B J i ( t ) c B J i ( t ) q D J i ( t ) + P k ∈ L out q k out ( t ) , (6) where j and k are the counters for total L in links ﬂowing into the junction and L out links extracting ﬂow from the junction; q j in ( t ) and q k out ( t ) are the inﬂows and outﬂows from these links connected to the junction; c j in ( t ) is the concentration in the inﬂow solute; q B J i ( t ) is the ﬂow injected to the junction with concentration c B J i ( t ) by booster station if located; and q D J i ( t ) is demand. 4) Mass Balance at T anks: Mass conservation in tanks assumes complete instantaneous mixing of all inﬂows, outﬂo ws, and stored water following the continuously stirred tank reactor (CSTR) model. V TK i ( t + ∆ t ) c TK i ( t + ∆ t ) = V TK i ( t ) c TK i ( t ) + X j ∈ L in q j in ( t ) c j in ( t )∆ t + V B TK i ( t + ∆ t ) c B TK i ( t + ∆ t ) − X k ∈ L out q k out ( t ) c TK i ( t )∆ t + R TK MS ( c TK i ( t )) V TK i ( t )∆ t, (7) 6 where V B TK i ( t + ∆ t ) is the v olume injected to the tank with concentration c B TK i ( t + ∆ t ) by booster station if located. R TK MS ( c P i ( x, t )) is the multi-species reaction rate in tanks expression (refer to Section ” Multi-species Reaction and Decay Model ”). C. Multi-species Reaction and Decay Model Dividing the model into decay and mutual reaction dynamics allows it to consider a substance with relati vely different reaction rates than the decay rate and for the model to be less sensitiv e to the other reactants’ concentrations. Decay model is a ﬁrst-order model that depends on only chlorine concentration and constant decay rate. Hence, the chlorine decay reaction rates for Pipe i and T ank j are k P i = k b + 2 k w k f r P i ( k w + k f ) , k TK j = k b , where k b is the bulk reaction rate constant; k w is the wall reaction rate constant; k f is the mass transfer coefﬁcient between the bulk ﬂow and the pipe wall; r P i is the pipe radius. The mutual reaction model is expressed by a second-order nonlinear ODEs which are discretized using Forward Euler method c ( t + ∆ t ) − c ( t ) = − k r ∆ t ( c ( t )˜ c ( t )) , ˜ c ( t + ∆ t ) − ˜ c ( t ) = − k r ∆ t ( c ( t )˜ c ( t )) , where c ( t ) , ˜ c ( t ) are the concentrations for chlorine and ﬁctitious reactant; and k r is the mutual reaction rate between them. Eventually , reaction expressions for pipes and tanks are R P M ( c P i ( s, t )) = − k r c P i ( s, t ) ˜ c P i ( s, t ) , R TK M ( c TK j ( t )) = − k r ˜ c TK j ( t ) c TK j ( t ) , (8a) R P M (˜ c P i ( s, t )) = − k r c P i ( s, t ) ˜ c P i ( s, t ) , R TK M (˜ c TK j ( t )) = − k r ˜ c TK j ( t ) c TK j ( t ) . (8b) A full description of the state-space matrices construction for the Upwind discretization schemes and an example on a simple three-node network (consists of reserv oir, a pump, a junction, a pipe, and a tank—Fig. 6) is included in [ 41 ] for reader’ s reference on how to formulate the representation for different network component. It is worth mentioning that this study validates the utilization of these EFD discretization schemes and the model performance as mentioned in comparison to EP ANET and its extension, EP ANET -MSX (WQ multi-species simulation tool). The comparison is considered reliable as the governing laws and equations are the same for all network components in both models. It should be noted that EP ANET+EP ANET -MSX employs the Lagrangian time-dri ven method, dividing each pipe into changing-sized segments. Whilst, the adopted EFD schemes in our study work within a ﬁxed grid, facilitating the construction of a state-space representation with ﬁnite dimensions. The drawback associated with these discretization schemes is the large dimensionality of the model. Ho wev er, main objectiv es of this study are to address this challenge by employing model order reduction techniques and to integrate the reduced-order multi-species model effecti vely into time-efﬁcient real-time feedback control algorithm, which are outlined and presented in detail in the next sections. On the other hand, coupling the EP ANET+EP ANET -MSX model with a real-time control algorithm is complex and presents challenges due to the need to handle changes in segment count and size per pipe at each simulation time-step, as well as being familiar with and able to leverage and use their toolkits in the coding language used (i.e., MA TLAB and Python). In the next section, we inv estigate different MOR algorithms for ( 1 ). I I I . M O D E L O R D E R R E D U C T I O N A N D T R A N S F O R M A T I O N O F M S - W Q M The state-space representations formulated in the previous section are in forms of nonlinear difference equations (NDEs) ( 1 ) with large numbers of variables resulted from high resolution spatio- temporal-discretization. T o reach the end-goal of this paper , which is controlling chlorine lev els for ( 1 ), we propose dif ferent methodologies to reduce the model order and showcase their limitations, accuracy , computational time, and robustness/sensiti vity to initial conditions and ﬁctitious reactant type. That being said, we list full descriptions of the methods cov ered in our framework. W e start with linearizing ( 1 ), then explain model order reduction and transformation for linearized and original nonlinear systems. A. Model Linearization The mutual reaction is expressed as a nonlinear term that can be linearized using T aylor series approximations [ 45 ]. By linearizing around operating points c o , ˜ c o , the nonlinear term R M ( c ( t ) , ˜ c ( t )) for both chemicals is expressed as: R M ( c ( t ) , ˜ c ( t )) = − k r ( c o ˜ c o + c o (˜ c ( t ) − ˜ c o ) + ˜ c o ( c ( t ) − c o )) , = − k r ( c o ˜ c o + c o ˜ c ( t ) − c o ˜ c o + ˜ c o c ( t ) − ˜ c o c o ) , = − k r ( c o ˜ c ( t ) + ˜ c o c ( t ) − ˜ c o c o ) . (9) For each of the chemicals, the mutual reaction after linearization is broken down to a term that depends on its concentration, a term that depends on the other chemical’ s concentration, and a constant. The general state-space representation ( 1 ) has a 7 (a) (b) Fig. 3: (a) Linear and (b) nonlinear MOR methods conﬁguration. block-diagonal matrix of A matrices with no dependency between the chemical except in the f function. That is, by applying linearization to the model the state-space representation is updated to linear difference equations (LDEs):   E 11 ( t ) 0 0 E 22 ( t )   | {z } E ( t )   x 1 ( t + ∆ t ) x 2 ( t + ∆ t )   | {z } x ( t +∆ t ) =   ˘ A 11 ( t ) ˘ A 12 ( t ) ˘ A 21 ( t ) ˘ A 22 ( t )   | {z } ˘ A ( t )   x 1 ( t ) x 2 ( t )   | {z } x ( t ) +   B 11 ( t ) 0 0 B 22 ( t )   | {z } B ( t )   u 1 ( t ) u 2 ( t )   | {z } u ( t ) + Φ , (10) where ˘ A 11 ( t ) and ˘ A 22 ( t ) are the modiﬁed diagonal matrices; ˘ A 12 ( t ) and ˘ A 21 ( t ) are the matrices gathering the dependency between the two species concentrations; and Φ is the vector containing the constants. Note that the changes in ˘ A 11 ( t ) and ˘ A 22 ( t ) from the original matrices are only in the sub-matrices/elements representing pipes and tanks only (i.e., A P P and A TK TK ). B. Model Order Reduction and T ransformation Algorithms In our study , we in vestigate two SVD-based projection methods; POD and BPOD. The reason behind not applying BT method is that it has been proven to be computationally impractical for the linear water quality model [ 31 ]. Both POD and BPOD are applied on the linearized MS-WQM, while an e xtension to the POD method is applied to reduce the nonlinear model where the nonlinear term is directly ev aluated. Before explaining the detailed approach of the aforementioned methods, we start by explaining the general approach of SVD-based methods where a snapshot of the original space is taken. For the general nonlinear state-space representation ( 1 ) which can concisely formulated as follows E ( t ) x ( t + ∆ t ) = A ( t ) x ( t ) + B ( t ) u ( t ) + f ( x ( t )) , y ( t ) = C ( t ) x ( t ) + D ( t ) u ( t ) , (11) the ﬁrst step is to map the representation states x ∈ R n x to another space state w ∈ R n x . This mapping aims to re-order the states according to their inﬂuence in the preserved property . Driven by the goal of applying control algorithm on our model, we care to capture the most controllable and observ able snapshots of the original space. Transformation is performed through constructing a non-singular matrix V ∈ R n x × n x , so that x = V w . That is, Eq. ( 11 ) is expressed in terms of w as follows E w ( t ) w ( t + ∆ t ) = A w ( t ) w ( t ) + B w ( t ) u ( t ) + V − 1 f ( V w ( t )) , y w ( t ) = C w ( t ) w ( t ) + D ( t ) u ( t ) , (12) where E w = V − 1 E V , A w = V − 1 AV , B w = V − 1 B and C w = C V . Next, the reduced-order model is captured from the transformed mapping with number of states n r ≪ n x donated by x r ∈ R n r . A snapshot is taken of x equal to V r x r where V r is the matrix comprised of the ﬁrst n r columns of V . Similarly , we deﬁne L r as the ﬁrst n r rows of V − 1 . Finally , the reduced-order model is expressed as E r ( t ) x r ( t + ∆ t ) = A r ( t ) x r ( t ) + B r ( t ) u ( t ) + f r ( x r ( t )) , y r ( t ) = C r ( t ) x r ( t ) + D ( t ) u ( t ) , (13) where E r = L r E V r , A r = L r AV r , B r = L r B and C r = C V r . 8 The choice of n r can be done arbitrarily as a ﬁxed number or to conserve a speciﬁed level of energy between ROM and FOM. The energy of a system is determined by the summation of its eignv alues, hence n r can be chosen to keep a certain energy percentage of FOM in R OM [ 46 ]. Ho wev er, we in vestigate choosing different numbers of n r for each case study where the energy persevered is increased with larger n r . Additionally , majority of MOR methods deal with original systems with zero initial conditions which does not align with the nature of water quality dynamics. Previously , the authors in [ 31 ] have dealt with that by recognizing the non-zero initials network-wide as inputs for the system and setting ˆ x ( t ) = x ( t ) − x (0) in the original model. W e follow same approach with further analysis for the nonlinear term of the mutual dynamics. The mutual reaction dynamics as stated in Section ” Multi- species Reaction and Decay Model ” take place in pipes and tanks. That is, vector f contains zeros except for states of pipes’ segments and tanks. W e deﬁne x MS 1 ( t ) := { c TK ( t ) , c P ( t ) } and x MS 2 ( t ) := { ˜ c TK ( t ) , ˜ c P ( t ) } . Accordingly , f ( x MS 1 ( t ) , x MS 2 ( t )) = α · x MS 1 ( t ) · x MS 2 ( t ) , (14) where α := { α TK , α P } ; α TK j = − k r ∆ t V TK j ( t ) V TK j ( t +∆ t ) ∀ j = 1 , . . . , n TK ; and α P l = − k r ∆ t ∀ l = 1 , . . . , P n P i =1 s L i . Henceforward, by setting ˆ x MS 1 ( t ) = x MS 1 ( t ) − x MS 1 (0) and ˆ x MS 2 ( t ) = x MS 2 ( t ) − x MS 2 (0) and substituting in ( 14 ) we get f ( ˆ x MS 1 ( t ) , ˆ x MS 2 ( t )) = α · ˆ x MS 1 ( t ) · ˆ x MS 2 ( t ) = α · ( x MS 1 ( t ) − x MS 1 (0)) · ( x MS 2 ( t ) − x MS 2 (0)) = α · ( x MS 1 ( t ) · x MS 2 ( t ) − x MS 2 (0) · x MS 1 ( t ) − x MS 1 (0) · x MS 2 ( t ) + x MS 1 (0) · x MS 2 (0) ) , (15) which proves that considering f ( ˆ x MS 1 ( t ) , ˆ x MS 2 ( t )) can be utilized by updating A ( t ) in the original model to eliminate the negati ve terms (in blue) for pipes and tanks, while the positiv e constant term (in green) encapsulates the nonlinear term at the initial concentrations which is already considered (refer to the online version of the paper for the actual colors). Subsequently , the full-order model is formulated as E ( t ) ˆ x ( t + ∆ t ) = ˆ A ( t ) ˆ x ( t ) + ˆ B ( t ) ˆ u ( t ) + f ( ˆ x ( t )) , y ( t ) = C ( t ) ˆ x ( t ) + ˆ D ( t ) ˆ u ( t ) , (16) where ˆ A ( t ) =   A 11 ( t ) ˆ A 12 ( t ) ˆ A 21 ( t ) A 22 ( t )   , , ˆ B ( t ) = [ B ( t ) A ( t ) x (0)] , ˆ D ( t ) = [ D ( t ) C ( t ) x (0)] , and ˆ u ( t ) = [ u ⊤ ( t ) 1 ⊤ ] ⊤ . On the other hand, for the linearized full-order model in ( 10 ) same approach as in [ 31 ] is followed and the ﬁnal model formulated as E ( t ) ˆ x ( t + ∆ t ) = ˆ A ( t ) ˆ x ( t ) + ˆ B ( t ) ˆ u ( t ) + Φ , y ( t ) = C ( t ) ˆ x ( t ) + ˆ D ( t ) ˆ u ( t ) . (17) where ˆ A ( t ) = ˘ A ( t ) , ˆ B ( t ) = [ B ( t ) ˘ A ( t ) x (0)] , ˆ D ( t ) = [ D ( t ) C ( t ) x (0)] and ˆ u ( t ) = [ u ⊤ ( t ) 1 ⊤ ] ⊤ . Lastly , we judge the performance of the MOR methods by calculate the root-mean-square error (RMSE) metric, RMSE = v u u t 1 N p N p X j =1 || y ( j ) − y ( j ) || 2 2 . (18) The error is calculated for a speciﬁc simulation period of N p time-steps through which we apply same system inputs u ( j ) to the two models. In the following sections, we give full description of the utilized methods. W e start with applying POD and BPOD for the linearized formulation of the system, followed by integrating and handling the nonlinearity in the original representation of the system (Eq. ( 1 ) for case of zero initial conditions and Eq. ( 16 ) for case of non-zero initial conditions). The basic and the balanced POD methods are considered data-driv en SVD methods. The main idea is to build empirical Gramians based on snapshots of the original system. These empirical Gramians avoid solving complicated, intractable in many case, L yapunov equations. POD method relies on constructing controllability Gramian while BPOD constructs ﬁnite horizon controllability and observability Gramians. Notably , POD method fav ors highly controllable states ov er highly observable but less controllable ones which BPOD averts by reﬂecting observability in the captured snapshot. It is important to highlight that in our system the concepts of contr ollability and observability for the two chemicals are different in what the y do reﬂect. While the input vector u 1 ( t ) depicts chlorine injections into the system by source or rechlorination stations, vector u 2 ( t ) enables simulating the intrusion of the contaminant to the system [ 47 ]. Henceforward, controllability for the second chemical is indicating which network components get exposed/affected by the contamination ev ent. On the other hand, typically water quality sensors are located to measure chlorine levels and from here comes the abstract concept of the system being observ able for water quality measurements. This is a main reason for chlorine monitoring to be a solid proxy of the water quality state in a speciﬁc network. Howe ver , no sensors are placed for contaminants detection speciﬁcally with their wide range. That is, their observability is reﬂected in chlorine lev els and not quantiﬁable in the matrix 9 Procedur e 1: POD for general MS-WQM 1 Construct snapshot X m as in ( 19 ) 2 if n x ≪ m then 3 Calculate W C m = X m X ⊤ m 4 Obtain transformation matrix V by applying eigen value decomposition W C m V = V Λ 5 else 6 Calculate f W C m = X ⊤ m X m 7 Obtain matrices forms of eigen vector and eigen value of f W C m ; Q and Λ 8 Calculate transformation matrix as V = X m Q Λ − 1 2 9 end if 10 Specify n r 11 Deﬁne V r as the ﬁrst n r columns of V 12 Deﬁne L r as the ﬁrst n r rows of V − 1 13 Calculate E r , A r , B r , and C r 14 if FOM is nonlinear then 15 Follo w Procedure 2 16 end if C 22 of ( 1b ) (i.e., a zero matrix). That puts a limitation on applying BPOD method as it will overlooks this contaminant because it is not observable. In Section ” Balanced Proper Orthogonal Decomposition (BPOD) ”, we propose a special approach to solve this issue. In addition, with no output measurement for that chemical, the RMSE metric in Eq. 18 only measures the error for chlorine. In fact, the main purpose of this work is to control and monitor chlorine under contamination events which makes it valid to focus on the output of measuring its concentrations that are accurately representing the real-time state. Nev ertheless, to ev aluate the performance of the applied MOR methods we assume the existence of ”imaginary” sensors on some speciﬁc nodes to measure the ﬁctitious reactant concentrations to calculate the corresponding error . In the following subsections, we explain what snapshots each method captures and how to construct these Gramians correspondingly . 1) Proper Orthogonal Decomposition (POD): This method captures snapshot matrix X m that is built for speciﬁc number of steps m by concatenating the states vector into X m = [ x (0) x (1) . . . x ( m − 1)] , (19) where X ∈ R n x × m . The approximate m -step controllability Gramian W C m is deﬁned as X m X ⊤ m ∈ R n x × n x . Ne xt, we apply eigen value decomposition (ED) W C m V = V Λ and obtain V whose columns are the corresponding eigen vectors. Howe ver , in many cases applying ED for an n x × n x matrix with large n x is taxing. This can be avoided in cases of m ≪ n x by constructing f W C m = X ⊤ m X m ∈ R m × m . Accordingly , the eigen v alue decomposition procedure performing is easier and requires less computational time [ 48 ]. In this case, ED is formulated as f W C m Q = Q Λ where Λ is the diagonal matrix of eigenv alues and matrix Q is assembled with eigen vectors as columns. The transformation matrix is then calculated as V = X m Q Λ − 1 2 . For detailed step-by-step depiction of the POD method, follow Procedure 1 . This procedure is followed for both chemicals. a) Mapping and Integr ating the Nonlinearity: While applying MOR, the reason behind separating the linear term(s) and the nonlinear term(s) is to be able to capture the behavior of the latter while working in a subspace of the original system (i.e., R n r instead of R n x ). In Eq. ( 13 ), following the projection of the whole system the nonlinear term is expressed as f r = L r f ( V r x r ( t )) . Y et, the computational complexity of the nonlinear term still depends on n x ; f r = L r |{z} n r × n x f ( V r x r ( t )) | {z } n x × 1 . Henceforward, it is proposed to reduce the nonlinear term based on an approximate hyperreduction approach. The approach is to measure not the full state-space variables, but particular points and from those points we construct the nonlinear term by interpolation around these points. In our study we specify the number of these points to equal n r ; f r = L r U f r | {z } n r × n r ˆ f ( t ) |{z} n r × 1 . The goal is to project f ( V r x r ( t )) onto U f r so that f ( V r x r ( t )) ≈ U f r ˆ f ( t ) and L r U f r can be pre-computed ofﬂine. This approach is called the ”Gappy method” of Galerkin pr ojection and the Discr ete Empirical Interpolation Method (DEIM) is used to reconstruct the nonlinear vector by interpolation. W e adopt a Gr eedy sampling algorithm to construct the measurement matrix to select the entries used. W e start by stacking numerical snapshot F m only for the nonlinear term, 10 Fig. 4: An illustrative example of applying the Greedy sampling algorithm to construct the measurement matrix K for the case of n r = 5 . F m = [ f ( x (0)) f ( x (1)) . . . f ( x ( m − 1))] , (20) followed by performing a separate SVD for that snapshot, F m = U f Σ F Q ⊤ f . The next step is to deﬁne a rank- n r approximating basis U f r as the ﬁrst n r columns of U f . Next, we construct the measurement matrix K by applying the Greedy sampling algorithm as summarized in Procedure 2 . As shown in Fig. 4 , the Greedy sampling algorithm starts by choosing the index with the maximum value in the ﬁrst mode u 1 and making it the ﬁrst measurement location. In the second iteration and the subsequent ones, we compute the residual to ev aluate how the current measurement subspace projects onto the next one and decide on the next measurement point. The reason behind choosing the measurement with the maximum residual is that the modes are no longer orthogonal in the support space, hence, we calculate the residuals and locate the index with the maximum residual. Procedur e 2: Nonlinearity handling in MOR 1 Capture F m as in ( 20 ) 2 Perform SVD of F m = U f Σ F Q ⊤ f 3 Construct U f r as the ﬁrst n r columns of U f 4 Start Gr eedy sampling algorithm for selecting the indices (entries of f ) Input: U f r = [ u 1 . . . u n r ] Output: I := { i 1 . . . i n r } and K = [ e i 1 . . . e i r ] 5 [ s, i 1 ] = max {| u 1 |} 6 U f r = [ u 1 ] , K = [ e i 1 ] 7 for I = 2 : n r do 8 solve K ⊤ U f r b = K ⊤ u I for b 9 q = u I − U f r b 10 [ s, i I ] = max {| q |} 11 U f r = [ U f r , u I ] , K = [ K , e i I ] 12 end for 13 Proceed 14 Calculate ˆ f ( t ) = ( K ⊤ U f r ) − 1 f ( K ⊤ V r x r ( t )) 2) Balanced Pr oper Orthogonal Decomposition (BPOD): The advance in the BPOD method is the reﬂection of both controllability and observability in ranking the states, unlike POD. This is attained by constructing two snapshots of the 11 system, f X m which captures the impulse responses when applying impulse signal as system input (i.e., u i ( m ) = γ ( m ) ) and P m is assembled from states p ( t ) obtained from the adjoint system with impulse response in the measurements as the system’ s output. For the linearized model in ( 17 ), the adjoint system can be expressed as follows, p ( t + ∆ t ) = ˘ A ⊤ ( t )( E − 1 ( t )) ⊤ p ( t ) + C ⊤ ( t ) y ( t ) + E − 1 ( t ) Φ . (21) Next step is performing SVD to the block Hankel matrix H m = P ⊤ m f X m = U Σ Q ⊤ then specifying n r to collect the largest n r singular values in Σ and obtain the corresponding left and right singular vectors (i.e., U r and Q r ). Accordingly , V r and S r are calculated as V r = f X m Q r Σ 1 2 r , L r = Σ 1 2 r U ⊤ r P ⊤ m (22) This approach is applicable for chlorine with sensors placed to measure its lev els. For the ﬁctitious reactant representing the contaminant, matrix C 22 in ( 1b ) is a zero matrix representing a non-sensed variables in our system. T o solve this issue we assume that the contamination event is detected and the source location is determined. This is considered a valid assumption in water quality monitoring to work backward detecting, classifying, and quantify using con ventional WQ sensors [ 49 ]. This is different than the ”imaginary” sensors that are aforementioned while calculating the error to ev aluate the performance of the applied methods. Another advance of the BPOD is the ability of stabilizing it by choosing the length of the snapshots to be large enough to represent the actual Graminas shooting for inﬁnity . W e adopt an a priori stabilization method to ensure that the snapshot captures the chemicals’ ev olution from the time it is injected in the system till it is observed by the furthest sensor . This is fulﬁlled by assembling the snapshots over a period exceeding m = max l T BS ∆ t m = max l P L BS i v BS i ∆ t m where L B S i and v B S i are the length and velocity of the pipes the chemical travels through from a booster station to the furthest sensor . W ith the existence of multiple booster stations and sensors and within the simulation period, m is taken as the length corresponding to the maximum travel time T B S . Accordingly , this method is affected by the actuators’ and sensors’ locations along the network. Lastly , Procedure 3 summarizes all the steps needed for a linear(ized) WQ model. Procedur e 3: BPOD for linear(ized) WQM 1 Obtain snapshots length m = m 2 Construct snapshot f X m and P m 3 Construct the block Hankel matrix H m = P ⊤ m f X m 4 Perform SVD of H m = U Σ Q ⊤ 5 Specify n r 6 Obtain U r , Σ r , and Q r 7 Calculate V r and L r via ( 22 ) 8 Calculate E r , A r , B r , and C r I V . R E A L - T I M E R E G U L A T I O N O F M S - W Q M V I A M O D E L P R E D I C T I V E C O N T RO L A N D M C C O R M I C K R E L A X A T I O N S The water quality control problem is formulated ov er simulation period [0 , T s ] and constrained by putting standard upper and lower bounds on chlorine concentrations stated by EP A regulations [ 50 ], which are x min 1 = 0 . 2 mg/L and x max 1 = 4 mg/L. W e note that the contaminant in the system is assumed to be detected and classiﬁed. Accordingly , for some toxic or health threatening substances a constraint can be introduced to be kept lower than the allowed concentration deﬁned by EP A. These bound for both chemicals formulates the constraint x min ≤ x ( t ) ≤ x max . Additionally , the control inputs for chlorine are constrained to be non-negativ e and limited by the chlorine av ailability and capacity of each booster station. The objective of this control problem is to keep chemicals concentrations in all network’ s components within the aforementioned bounds while minimizing the cost of chlorine injections. That being said, the problem formulation is as follows, minimize x ( t ) , u 1 ( t ) J ( u 1 ( t )) = µ N p X t =1 q B ( t ) ⊤ u 1 ( t ) (23a) subject to WQM ( 1 ) , x min ≤ x ( t ) ≤ x max , u min 1 ≤ u 1 ( t ) ≤ u max 1 , (23b) where problem variables x ( t ) and u 1 ( t ) are chemicals concentrations network-wide and chlorine injections through booster stations, q B ( t ) is the ﬂo w rates at the nodes corresponding to the locations of the booster stations, µ is the unit cost of chlorine in $/mg, and WQM is the water quality model we are simulating and controlling following the representation in ( 1 ). Finally , N p is the number of time-step in the simulation period, N p = T p ∆ t . 12 (a) (b) Fig. 5: (a) Discrete MPC prediction horizon scheme, and (b) graphical representation of McCormick env elope relaxation. Nonetheless, this problem has large number of variables x ( t ) and u 1 ( t ) . This issue can be solved by transforming a constrained LP ( 23 ) to a quadratic program (QP) with fewer variables by applying real-time constrained model predictive control (WQ-MPC). The water quality control formulated in [ 39 ] is based on the linear state-space representation of the single- species WQ dynamics. In addition, same control algorithm is applied in [ 31 ] for the reduced order model of the single-species representation and it prov ed its validity and effecti veness. F or bre vity , we do not include the details and the deri vation in this paper for the case of linearized MS-WQM. Eventually , the WQC problem is formulated as quadratic program. For the nonlinear MS-WQM, the nonlinearity in the constraint can be relaxed using McCormick en velopes and integrated back to the original constrained control problem as explained in the following section. A. McCormick Relaxation The nonlinear term in the constraint is formulated as a bi-linear expression depending on the concentration of two chemicals at a speciﬁc network component. This constraint can be relaxed using McCormick Relaxation for bilinear nonlinear problems [ 51 ]. This method turns the bilinear term into two en velopes surrounded by ov erestimators and underestimators to work within, Fig. 5b . For a bilinear expression z = x 1 x 2 where x 1 and x 2 are the two chemicals concentrations under x min 1 ≤ x 1 ≤ x max 1 and x min 2 ≤ x 2 ≤ x max 2 , z is introduced as a new decision variable with the following constraints, z ≥ x min 1 x 2 + x 1 x min 2 − x min 1 x min 2 , z ≥ x max 1 x 2 + x 1 x max 2 − x max 1 x max 2 , z ≤ x max 1 x 2 + x 1 x min 2 − x max 1 x min 2 , z ≤ x min 1 x 2 + x 1 x max 2 − x min 1 x max 2 . (24) W e note that, in some cases the upper bound on x 2 is not speciﬁed or its concentration initially is lower than the maximum allowed one stated by EP A. In such cases, we specify x max 2 to be equal to this initial concentration detected to be able to tighten the overestimators env elope while having the minimum equal to zero. Eventually , the problem formulation explained for the linearized model can be adopted with these modiﬁcations. First, a new variable vector z ( t ) is introduced and it replaces f ( x 1 , x 2 , t ) in ( 1a ). Additionally , the total number of the constraints added to the optimization problem via ( 24 ) is equal to 4( n TK + P n P i =1 s L i ) as the nonlinear term is deﬁned for pipes’ segments and tanks and is the same for both chemicals at same element of the aforementioned (refer to Eq. ( 8 )). T o that end, the WQC problem described in ( 23 ) is modiﬁed as follows, minimize x ( t ) , u 1 ( t ) , z ( t ) J ( u 1 ( t )) = µ N p X t =1 q B ( t ) ⊤ u 1 ( t ) (25a) subject to WQM ( 1 ) , x min ≤ x ( t ) ≤ x max , u min 1 ≤ u 1 ( t ) ≤ u max 1 , McCormic k ( 24 ) (25b) 13 Next step is transforming ( 25 ) into a linear augmented formulation based on which the ﬁnal WQC-QP is built. First, by introducing z ( t ) into ( 1a ), the state-space representation is updates as x ( t + 1) = A ( t ) x ( t ) + B ( t ) u ( t ) + β z ( t ) . (26) where β = − k r . Then, we deﬁne the change in the states and inputs as follows ∆ x ( t + 1) = x ( t + 1) − x ( t ) , ∆ u ( t + 1) = u ( t + 1) − u ( t ) , ∆ z ( t + 1) = z ( t + 1) − z ( t ) . (27) T o concatenate these rates of change in ( 26 ), ∆ z is assembled to the vector of systems decision inputs to be optimally chosen within the en velopes deﬁned by ( 24 ). Eventually , we reach an the augmented state-space representation in ( 28 ).   ∆ x ( t + 1) y ( t + 1)   | {z } x a ( t +1) =   A ( t ) 0 C ( t ) A ( t ) I   | {z } Φ a   ∆ x ( t ) y ( t )   | {z } x a ( t ) +   B ( t ) β C ( t ) B ( t ) β C ( t )   | {z } Γ a   ∆ u ( t ) ∆ z ( t )   | {z } ∆ u a ( t ) (28) This augmented representation can be abstractly rewritten as, x a ( t + 1) = Φ a x a ( t + 1) + Γ a ∆ u a ( t ) . T o av oid redundancy , integrating this equality into WQC-MPC formation follows the same approach of [ 31 ] reaching the ﬁnal QP [ 31 , Eq. (38)]. On another note, the added constraints expressed in ( 24 ) are incorporated in the constraints on the optimization variables. B. Generalized Comprehensive W ater Quality Modeling and Control F rame work In our study , we have covered model order reduction and control for multi-species water quality dynamics where chlorine is reacting with another source of contamination in form of a bi-linear expression—refer to Section ” Multi-species Reaction and Decay Model ”. Howev er, there are other formulations for single-species and multi-species chlorine bulk decay and reaction dynamics as listed in [ 41 ]. W e include a short list of these formulations in T ab . I , ne vertheless for more details and descriptions refer to the aforementioned study . The following generalized framework described in Algorithm 1 maps out the methods adopted in this study to be applied on the different decay and reaction models. T ab. I: Chlorine bulk decay and reaction models expressions M# Model Model formulation #States L/NL a M-1 First-order dc dt = − k c ( t ) n x L M-2 First-order with stable component dc dt = − k ( c ( t ) − c L ) n x L M-3 Parallel ﬁrst-order dc 1 dt    fast = − k fast c 1 ( t ) dc 2 dt    slow = − k slow c 2 ( t ) c t ( t ) = c 1 ( t ) + c 2 ( t ) 2 n x L M-4 Parallel second-order dc F dt    fast = − k fast c ( t ) c F ( t ) dc S dt    slow = − k slow c ( t ) c S ( t ) dc dt = dc F dt + dc S dt 2 n x NL M-5 n th -order dc dt = − k c n ( t ) n x NL M-6 n th -order with stable component dc dt = − k ( c ( t ) − c L ) c ( n − 1) n x NL M-7 Second-order with ﬁctitious component dc dt = − kc ( t ) ˜ c ( t ) d ˜ c dt = − kc ( t ) ˜ c ( t ) 2 n x NL M-8 Second-order with multiple components dc i dt = − k i c ( t )˜ c i ( t ) d ˜ c i dt = − k i c ( t )˜ c i ( t ) dc dt = I X i dc i dt I n x NL a L: Linear or NL: Nonlinear model expression. For the ﬁrst-order , ﬁrst-order with stable component, and parallel ﬁrst-order (M-1&M-2&M-3) models, the dynamics are linear and accordingly follow the procedure of the linearized model represented in our study . Whilst, the second-order with 14 multiple components (M-8) is considered to be the same formula as the second-order with ﬁctitious component (M-7) we cov er in this paper except for the number of states which gets multiplied by the number of reactants in the system. That is, model order reduction for M-8 model becomes more demanded. On the other hand, the parallel second-order model (M-4) is a special form of the second-order with ﬁctitious component. Lastly , the n th -order without and with stable component models are higher order models which can be reduced as a nonlinear models or be transferred into quadratic approximation and apply piecewise linear relaxation. Algorithm 1: Generalized water quality modeling and control framew ork Input: WDN topology , components’ characteristics, and hydraulics parameters Output: Real-time water quality states x ( t ) and control inputs u ( t ) at time t of a simulation period of T s 1 Initialization 2 Deﬁne ∆ t , number of segments s i for each pipe and accordingly n x 3 Formulate WQ state-space representation ( 1 ) as explained in Section ” State-space Multi-species W ater Quality Model ” and according to the reaction dynamics in T ab. I . 4 Proceed 5 if Applying M-1/M-2/M-3 reaction model then 6 Follo w Procedure 3 to obtain ROM 7 Apply constrained real-time WQ-MPC on ( 23 ) 8 else if Applying M-4/M-7/M-8 r eaction model then 9 if F ollowing Pr ocedur e 1 then 10 Apply McCormick relaxation via ( 24 ) 11 Apply constrained real-time WQ-MPC on ( 25 ) 12 else 13 Linearizing and following Procedur e 3 then 14 Apply constrained real-time WQ-MPC on ( 23 ) 15 end if 16 else 17 Applying M-5/M-6 then 18 Follo w Procedure 1 to obtain ROM 19 T ransforming into quadratic app./Apply piecewise linear relaxation 20 Apply constrained real-time WQ-MPC on ( 23 ) 21 end if T o recapitulate, this paper is an extension of our previous work in [ 31 ], [ 41 ]. That is, some methods/aspects included in this section hav e already been covered in these studies. Y et, we hav e decided to reintroduce this material into our study in a more concise way to ensure that the reader can effecti vely follow and comprehend the new information being presented. In the following bullet points, we highlight the nov elty in our work in comparison to these studies. • Adopting a nonlinear multi-species water quality dynamics rather than the linear single-species dynamics in a control framew ork. These dynamics enable a more heightened level of realism in the system dynamics representation. • Following two different paths where different MOR methods and control algorithms are applied on the original nonlinear and a linearized forms of the model. For the linearized model, we implement the same MOR techniques (speciﬁcally , POD and BPOD) described in the paper by [ 31 ]. On the other hand, for the nonlinear model, we introduce the Gappy method, which employs a greedy algorithm to effecti vely handle the nonlinearity and reduce the model dimension. • Expanding the implementation of these MOR techniques for the case of non-zero initial conditions by developing a closed formulation that preserves the original nonlinear formulation of the model. • Likewise, for the linearized model we implement the MPC algorithm explain in [ 31 ] to control and regulate chlorine lev els under the multi-species dynamics. In contrast, we extend the MPC algorithm to incorporate the McCormick Relaxation technique, which is speciﬁcally tailored for the nonlinear model. • While the same methods and algorithms are employed for the linearized model as in the linear single-species model described in the paper by [ 31 ], special consideration is required when implementing these techniques for the linearized model. This is primarily due to the duplication of state numbers and the distinct construction of representation matrices for the ﬁctitious reactant. These factors necessitate a speciﬁc approach to ensure accurate and reliable results during the implementation of these methods in the context of the linearized model. • Inv estigating and e valuating using an explicit vs. implicit discretization schemes from a control-theoretic standpoint. Speciﬁ- cally , we apply Upwind schemes, which offer a more accurate representation of the advection-reaction 1-D partial differential equations (PDEs) compared to the Lax-W endroff scheme utilized in the work by [ 31 ]. The advance in implementing an Upwind scheme is prov ed and demonstrated in [ 41 ]. In addition, Lax-W endroff scheme is an explicit scheme and the prior study does not extensi vely in vestigate the use of an implicit scheme in this context. • A novel component of our study is in vestigating the water quality control framework performance under different system hydraulic settings. These settings directly impact the water quality dynamics and their progression over time within the same 15 network. This exploration adds a unique dimension to our study , shedding light on the interplay between system hydraulics and water quality dynamics for enhanced understanding and improved control strategies. The v alidity of these techniques and the performance of the framework are ev aluated and substantiated in the subsequent section through a series of numerical case studies. Fig. 6: Case studies’ layouts and their components count: (a) Three-node network, (b) Net1, and (c) FFCL-1 with the zone we control framed. V . C A S E S T U D I E S This section demonstrates the proposed framework for model order reduction and control of MS-WQM. Particularly , we attempt to answer the following questions: ▷ Q1: How does the number of operating points impact agreement between the linearized model and the nonlinear MS- WQM? ▷ Q2: How effecti ve are the proposed MOR producers in terms of accuracy and computational time when applied on the MS-WQM? ▷ Q3: How sensitiv e is the performance of MOR and control algorithm to the discretization methods and system’ s hydraulics? ▷ Q4: How reliable and rob ust is model predictive control when applied to control chlorine le vels under multi-species dynamics? Numerical studies in this section are performed on three different networks, three-node, Net1, and FFCL-1 networks [ 42 ]. As shown in Fig. 6 , each of the networks has different topologies and scales. The three-node network is a self-designed network to help provide simple illustrations for different approaches throughout our framework implementation. Net1 includes different types of network components and has a looped layout. The FFCL-1 network is based on the Fairﬁeld, CA, USA water distribution system on which we test the scalability of our framework and its performance with scattered dead-ends. Also, Fig. 6 illustrates and lists each of the networks’ components. In addition to the listed components for each of the test networks in Fig 6 , each network has a different number of sensors and booster stations. The three-node network has one booster station at Junction J1 and one sensor at T ank TK1. Net1 has two stations at Junctions 1 and 6 and sensors at Junctions J4 and J9. Lastly , the controlled region of the FFCL-1 network has two sensors at Junctions J56 and J67, and one rechlorination station at J89. It is worth mentioning that for any WDN, the system dimension depends on the hydraulics parameters and water quality simulation time-step which accordingly deﬁne the number of segments for each pipe (i.e., pipes state v ariables). Further, 16 0 360 720 1080 1440 Time (mins) 0 0.25 0.5 0.75 1 CL Conc. (mg/L) (a) 0 360 720 1080 1440 Time (mins) 0 0.1 0.2 0.3 0.4 FR Conc. (mg/L) 0 360 720 1080 0 5 10 -3 (b) Fig. 7: Nonlinear vs linearized models results for (a) chlorine and (b) ﬁctitious reactant at TK1 of the three-node network. changing the velocities and ﬂows from one scenario to another results in distinct chemical concentrations across the network components for each scenario. With that in mind, in some of our case studies, we feature the effect of changing the hydraulics for the same network. In some other case studies, we ﬁx the hydraulics setting in the system to inv estigate/test a technique or an approach under discussion. In addition, for all studies performed in this section we use the Implicit Upwind scheme except for Section ” Implicit vs. Explicit Discretization Schemes under Control-theoretic Perspective ” where we compare its performance with the Explicit Upwind scheme from a control-theoretic perspectiv e. A. Nonlinear vs. Linearized Models Studies [ 52 ]–[ 54 ] state that applying a linear MOR algorithm on a linearized system gives satisfactory performance when the linearized system is close to the original nonlinear one or operating within/near its linear regime. In these studies, linearization is performed around one operating point for the whole simulation horizon. W e apply the same approach by linearizing around two operating points, (0,0) and (0.2,0.05) mg/L for chlorine and ﬁctitious reactant respectively at T ank TK1 of the three-node network. In this scenario, a constant demand is dra wn from J1 and sources of 2 mg/L of chlorine and 0.5 mg/L of the ﬁctitious reactant are provided at R1 and zero initial conditions for other network components. As demonstrated in Fig. 7 , linearization around the operating point of (0,0) results in higher concentrations compared to the nonlinear model (based on the NDE ( 1 )) for both chemicals due to the fact it drops out the nonlinear term and neglects the mutual reaction. On the other hand, linearizing the model around one random operating point as (0.2,0.05) mg/L results in relatively closer values for chlorine concentrations but not as close for the ﬁctitious reactant. Furthermore, unlike this scenario, in real-time water networks hydraulics are not ﬁxed and demands are time-variant resulting in chemical ev olution with different schemes for which ﬁxing the operating point for all elements is not actively ef ﬁcient. That is, we inv estigate next taking different operating points for each network component along the simulation window every speciﬁc number of time-steps. The choice of the operating points we linearize around is critical. The narrowest the recurrent window of choosing the operating points, the closest the results to the original model. Ho we ver , if we choose to update the operating points each water quality time-step then matrices ˘ A 11 ( t ) , ˘ A 12 ( t ) , ˘ A 21 ( t ) and ˘ A 22 ( t ) in Eq. ( 10 ) should be updated that frequent instead of being updated each hydraulic time-step. Hydraulic time-step is acceptable to be within an hourly scale to reﬂect the change in demand, while the range for water quality is between minutes and seconds to allow a stable numerical simulation [ 43 ], [ 55 ]. Consequently , updating the aforementioned matrices every WQ time-step adds more computational burden to the simulation which negates the main reason for implementing linearization and model order reduction. On the contrary , widening the window to be more than the hydraulic time-step especially in cases with signiﬁcant demand change gives inaccurate approximation of the system’ s behavior . Over and above that, it is important to consider falling within the control algorithm prediction and control horizon to be able to adjust accordingly with the controller input. W ith the hydraulic setting of a patterned demand at J1 changing every 1 hr (Fig. 8c ), the model is linearized around operating points that are taken every 1 hr for each of the network elements. The same sources of chemicals are provided at R1 with zero initial conditions for the other components. Results, shown in Fig. 8a and 8b for chlorine concentrations at TK1 and P1, exhibit that updating operating points every 1 hr results in accurate representation in comparison to the original model, except for the ﬁrst hour during which operating points are taken to be the initial concentrations at those elements. T o mitigate this issue, operating points are updated after 1-10 minutes from the simulation start. The same approach is followed in scenarios where chemical dosages are increasing locally at some node for elements downstream of this node. B. MS-WQ Model Or der Reduction P erformance In this section, we assess and compare the performance of each of the proposed model order reduction procedures for multi- species water quality dynamics in terms of accuracy compared to the original full-order model, and computational time. For 17 0 360 720 1080 1440 Time (mins) 0 0.5 1 1.5 2 CL Conc. (mg/L) P1 NDE LDE-1hr 15 270 1.8 1.85 1.9 1.95 2 (a) 0 360 720 1080 1440 Time (mins) 0 0.2 0.4 0.6 0.8 1 CL Conc. (mg/L) TK1 10 140 270 0.1 0.3 0.5 0.7 0.9 (b) 0 360 720 1080 1440 Time (mins) 0 0.2 0.4 0.6 0.8 1 Multiplier Demand Pattern (c) Fig. 8: Chlorine concentrations at (a) P1 and (b) TK1 of the three-node network with (c) patterned demand at J1. Results are for the nonlinear and linearized models—linearization operating points are updated ev ery 1 hr for all network components. each network, we apply POD and BPOD on the linearized model and extended POD for the nonlinear model. W e refer to these procedures as LPOD, LBPOD, and NLPOD, respectiv ely . W e note that we record the computational time needed for assembling the snapshots, obtaining the transformation matrices, and calculating the RMSE between the original and reduced-order model for a speciﬁc simulation under the same conditions. First, we apply the three MOR methods on the three-node and Net1 networks under zero and non-zero initial conditions and static hydraulic proﬁles. The results shown in Figures 9 and 10 validate that all methods are able to reduce the model dimensions with relativ ely low RMSEs for different n r values. These RMSEs get lower with increasing the n r values and are lower for the scenario of zero initial conditions compared to the case of non-zero initial conditions. F or the scenario with non-zero initial conditions, initial chlorine concentrations are 0.5 mg/L network-wide; the initial ﬁctitious reactant concentrations at TK1 in the three-node network and T ank 11 in Net1 are 0.05 mg/L. Fig. 10 shows the chlorine and ﬁctitious reactant concentrations for both scenarios of initial conditions at TK1 of the three-node network and T ank 11 of Net1 for the full-order model and the reduced-order models using all three MOR producers. It is observed that the reduced-order models giv e almost identical results to the full-order one for the step response at TK1 and for a regular node along the network for T ank 1 for the two scenarios of zero and non-zero initial conditions. On the contrary , these results differ from [ 31 ], where the POD method was found to have higher errors for the scenario of non-zero initial concentrations under single-species dynamics as the input-output relationship is not correctly captured when the initial values are treated as inputs into the system. In our study , this effect is mitigated by building the ofﬂine snapshot with a higher impulse signal by the booster stations which results in fav oring the actual locations of booster stations. Meanwhile, it is worth mentioning that MOR methods’ performance is signiﬁcantly impacted by the locations of the sensors and actuators and their reﬂection on network-wide observability and controllability . This leads to inaccurate or unstable results in some cases and in some other scenarios. Howe ver , the allocation of these sensors and actuators for each network is out of this paper’ s scope and we solve assuming the predetermination of their locations. C. Model Order Reduction Sensitivity to System Hydraulics The construction of the transformation matrices V r and L r for both methods POD and BPOD is sensitive to the snapshots (i.e., X m and P m ) constructed of ﬂine. These snapshots need to be long enough and representati ve of the actual reaction 18 Three-node Network 1e-12 1e-10 1e-8 1e-6 1e-4 n r = 30 n r = 50 n r = 100 n r = 150 LPOD BPOD NLPOD Net1 1e-12 1e-10 1e-8 1e-6 1e-4 n r = 30 n r = 50 n r = 100 n r = 150 (a) Three-node Network 1e-12 1e-10 1e-8 1e-6 1e-4 n r = 30 n r = 50 n r = 100 n r = 150 Net1 1e-12 1e-10 1e-8 1e-6 1e-4 n r = 30 n r = 50 n r = 100 n r = 150 (b) Fig. 9: RMSEs for the three proposed MOR methods for the three-node ( n x = 204 ) and Net1 ( n x = 482 ) networks with (a) zero and (b) non-zero initial conditions for different n r values. (a) (b) Fig. 10: Chlorine and ﬁctitious reactant concentrations evolution at (a) TK1 of three-node network ( n x = 204 ) and (b) T ank 11 of Net1 ( n x = 482 ) under zero (in blue) and non-zero (in red) initial conditions simulated by full- and reduced-order models with n r = 30 for both networks. between states, inputs, and outputs. That leads to being sensitive to the hydraulic settings of the system while capturing these snapshots and also while applying the desired model reduction. Dynamic hydraulic states in a network reﬂect the consumers’ patterned consumption which can be recorded for a speciﬁc network during a speciﬁc season [ 56 ]. After validating the reliability of the three MOR methods under zero and non-zero initial conditions, we in vestigate the 19 Fig. 11: Chlorine and ﬁctitious reactant concentrations evolution at J11, J56, J76, and J107 of FFCL-1 network simulated by full- and reduced-order models. Full-order model results for at each of the junctions are in solid lines, while the LPOD method results are in dashed lines. Number of states for the full-order model is n x = 10356 and reduced to n r = 200 states. case of dynamic hydraulic demands for a bigger network; the FFCL-1 network. In Fig. 11 , the ev olution of chlorine and the ﬁctitious reactant at J11, J56, J76, and J107 of the FFCL-1 network simulated by full-order model and LPOD-based reduced- order model is presented. Note that, only LPOD is shown, which is representativ e of the behavior of all other approaches. In this scenario, an input of 0.3 mg/L for the ﬁctitious reactant is inserted at the start of the network (i.e., at the T ank) depicting an early intrusion eve nt. As demonstrated, the LPOD-based ROM is able to trace the concentrations of the chemicals at different junctions, including dead-ends and junctions that are connecting looped pipes. Nonetheless, an oscillatory effect is detected for the ﬁctitious reactant concentrations in the framed zone. This oscillation is formulated as the ﬁctitious reactant being completely consumed by the chlorine at these junctions or at pipes ﬂowing into them (e.g., J76), howe ver , the operating points around which the system is linearized force the ﬁctitious reactant to have false concentrations. Therefore, this effect is illuminated by applying NLPOD and is reduced by updating the operating points more frequently . Lastly , the computational time recorded for each of the MOR methods implementations on the three tested networks is illustrated in Fig. 12 . F or all networks, the NLPOD method requires more computational time as a result of handling the nonlinearity term separately and performing the greedy sampling algorithm. Ho we ver , the maximum increase in time is around 95 seconds compared to BPOD for the FFCL-1 network, which is considered an acceptable computational time for a network of n x = 10356 states. Three-node 0.57 0.63 0.94 LPOD BPOD NLPOD 0 0.2 0.4 0.6 0.8 1 Computational time (sec) Net1 3.1 4.8 6.98 LPOD BPOD NLPOD 0 1 2 3 4 5 6 7 FFCL-1 223.518 250.43 344.8 LPOD BPOD NLPOD 0 50 100 150 200 250 300 350 Fig. 12: Computational time to implement the three MOR methods for three tested networks. T otal number of states n x is 10356 for FFCL-1, 482 for Net1, and 204 for the three-node network. D. Implicit vs. Explicit Discretization Schemes under Contr ol-theoretic P erspective As stated in Section ” T ransport and Reaction in Pipes ”, the 1-D AR equation can be discretized by implementing either Explicit or Implicit Upwind schemes. The explicit scheme needs to be performed under satisﬁed CFL condition to ensure 20 stability which requires a small time-step in many cases and hence a higher system dimension. While the implicit scheme is unconditionally stable b ut requires more complicated mathematical calculations that add to the computational work. Therefore, it has been a pressing question that needs to be answered, ” Which is better: Implicit or Explicit discr etization schemes? ”. This has been prov en to not be an easy question with an easy answer . In our study , we reduce our system’ s dimensions while applying either of these methods. Nonetheless, while transformation matrices are calculated ofﬂine, some system matrices are updated ev ery hydraulic time-step. This adds more computational load with matrices multiplication which is higher with matrix in verse in the case of the implicit scheme. One more point to highlight, although the implicit scheme allows a bigger simulation time-step, a smaller one is more efﬁcient to be able to update the control inputs more frequently . So our question can be formulated as follows, ” F r om a control-theor etic perspective, which is better: Implicit or Explicit discretization schemes? ” 50 100 150 200 320 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 RMSE Explicit Implicit (a) 0 360 720 1080 1440 Time (mins) 0 0.5 1 1.5 2 2.5 3 Concentration (mg/L) J4 Exp-Full Exp-Red Imp-Full Imp-Red 5 10 2 2.2 2.4 2.6 2.8 3 1 4 0.5 1 1.5 2 2.5 (b) Fig. 13: (a) RMSEs for the reduced-order Explicit and Implicit Upwind schemes-based models while applying LPOD on Net1 with n x = 482 for different n r values, (b) chlorine concentrations at J4 simulated by the full- and reduced-order models (with n r = 100 ) using both schemes. As model order reduction is a prior step to applying control to our model, we test both discretization methods’ performance while applying the LPOD method for Net1. As demonstrated in Fig. 13a , the RMSEs are lower for the implicit than the explicit scheme. In addition, the change in RMSEs by increasing n r more than 150 is insigniﬁcant as the states which get retained do not hold high energy compared to the previously selected ones. On the other hand, the error values for the explicit scheme do not go lo wer than 0.003 with increasing n r —explained through the following example. Although the CFL condition is satisﬁed for the explicit scheme, it formulates sharp fronts at points with a relatively signiﬁcant change in the chemical concentrations as sho wn in Fig. 13b . Fig. 13b illustrates chlorine concentrations at J4 for both implicit and explicit schemes and the corresponding reduced models with n r = 120 of a full model with n x = 482 . It is noticeable that the reduced model based on the explicit scheme is consequently affected and shows instability behavior that dampens reaching equilibrium. This performance is recorded under a low Courant number with the network’ s pipes. T o mitigate that, the water quality time-step is required to be reduced to reach higher CN—near but less than 1. Such behavior is av oided when applying the unconditionally stable implicit scheme. T o that end, the Implicit Upwind scheme giv es more accurate results that lead to a more robust control algorithm. The computational burden of this scheme can be lowered using sparse matrix multiplication. The computational time to perform the simulation on Net1 shown in Fig. 13a is 32.9 and 43.4 seconds for the Explicit and Implicit Upwind schemes, respectiv ely for the same water quality time-step. Howe ver , the implicit scheme retains high accuracy under a higher WQ time-step while requiring lo wer computational runtime. Therefore, the Implicit Upwind scheme giv es more ﬂexibility in choosing the time-step needed to retain real-time control windows while maintaining high accuracy . E. Real-time Contr ol Implementation of MS-WQM MOR-Based MPC The main objectiv e of this paper and the prior inv estigation of the MOR methods is to integrate them into and apply a real- time control algorithm of chlorine concentrations using the booster stations distributed along the WDN under MS-WQM. W e apply the MPC algorithm on the linearized and nonlinearized MS-WQ ROM as explained in Section ” Real-T ime Regulation of MS-WQM via Model Predictiv e Control and McCormick Relaxations ”. As both LPOD and BPOD can reduce the MS-WQM effecti vely , For the linearized model we apply the BPOD method. On the other hand, we apply the NLPOD for the nonlinear model to obtain the ROM. For multi-species water quality control and regulation, while applying the McCormick relaxation the en velopes rely on the limits for both chemicals. For the network’ s components near the location of the second chemical intrusion, these en velopes put 21 0 30 60 90 120 Time (mins) 1 1.2 1.4 1.6 1.8 2 Control Input (mg/min) # 10 4 SS-LMPC MS-LMPC MS-RMPC (a) 0 30 60 90 120 Time (mins) 0 0.5 1 1.5 CL Concentration (mg/L) P1 J1 (b) Fig. 14: (a) Control action u 1 during 2 hrs of simulation on the three-node network by applying, SS-LMPC: linear single-species-based MPC, MS-LMPC: linearized multi-species-based MPC, and MS-RMPC: nonlinear multi-species-based relaxed MPC, (b) chlorine concentrations at J1 and P1 under another chemical intrusion at J1 for the ﬁrst hour . tight boundaries on the chosen v alue for z by the control problem as the x 2 is close to x max 2 . On the other hand, for components downstream from this location with lower concentrations for both chemicals the relaxation allows higher and lower values for z which leads to choosing a value of z to be as close to the underestimators so that the control inputs are lower and the cost of chlorine injections is reduced. Additionally , for higher values of the mutual reaction coefﬁcient k r , the effect of relaxation on the chosen control input increases. That is, the proposed relaxed MPC may result in overlooking/underestimating the mutual reaction and, therefore, we lower the upper bound for the ﬁctitious reactant as a procedure integrated into the looped control algorithm repeated each time-step. As explained in Section ” Nonlinear vs. Linearized Models ”, it is proposed to update the operating points around which the system is linearized ev ery signiﬁcant change window (e.g., hydraulic states change). Updating these operating points adds to the computational time by recalculating the matrices, yet it yields a more accurate representation. Therefore, we put a threshold according to which we judge changing these points after applying ev ery control input. By adopting these approaches, we start with applying the MS-WQC MPC-based method on the three-node network under a static hydraulic proﬁle and with a reduced number of states of n r = 30 states out of n x = 204 . The water quality time-step is chosen to be 5 seconds and the control horizon is 10 minutes. The ﬁctitious reactant is discharged into the system at J1 (same location of the booster station) at a concentration of 0.1 mg/L for the ﬁrst 1 hour of the simulation duration. Fig. 14 demonstrates the control actions and the corresponding control response in J1 and P1 under the multi-species linearized R OM, nonlinear ROM, and the single-species ROM that neglects the existence of the other chemical in the system for the ﬁrst consecutive 2 hours of simulation. For all scenarios chlorine concentrations at J1 and P1 are zeros at the start of the simulation. That is, MPC starts by injecting high chlorine dosage of 21284 mg/min for the case of multi-species dynamics and 20838 mg/min for the single-species model. The control input needed drops to 19158 mg/min and 17596 mg/min for multi-species and single-species dynamics, respectiv ely . After the ﬁrst hour of simulation, MPC results in the same control actions for both models as the intrusion ev ent is contained. Note that, the second substance’ s initial concentration for P1 is zero, which leads to the peak control action at the start of the simulation being relatively close as the second substance has not traveled into P1. Comprehensiv ely , this highlights the importance and effecti veness of the adopted MS-WQM and control framew ork. Furthermore, this difference between the two models’ results (i.e., chlorine concentration dynamics and optimal chlorine inputs) increases for more reactiv e components with chlorine and initial intrusion concentrations which may cause operational issues with limited chlorine av ailability and/or budget. Additionally , the linearized MS-WQC problem and relaxed one produce the same performance as illustrated in Fig. 14a . W e note that for the linearized model, the operating points are updated at the start of the simulation, with applying the peak control action, and by the end of the contamination e vent. For the relaxed MS-WQC problem, all elements are directly affected by the event resulting in tight en velops and approximating the mutual reaction near its actual value. Howe ver , the number of control variables for this procedure is higher for the ﬁrst hour . T o that end, the computational time needed for each of the two control procedures is case-oriented. For this case study , the linearized-based MPC method has computational time of 78 seconds, while it is 93 seconds for the second method. Next, we apply the proposed MS-WQC approach on Net1 under a dynamic hydraulic proﬁle deﬁned by the patterned demand at Junction 1, Fig. 15c . The FOM has 482 states which are reduced to 50 states. As both control procedures hav e proved their ability to regulate chlorine concentrations network-wide, we showcase the results from the relaxed MPC procedure only to point out case studies that can take place. In this case study , the water quality time-step is 5 seconds while the control horizon is 10 minutes and the simulation period is 24 hours. The initial concentrations of all chemicals are zeros. The ﬁctitious reactant 22 0 360 720 1080 1440 Time (mins) 10 0 10 1 10 2 10 3 10 4 10 5 Control Action (mg/min) J1 J6 (a) 0 360 720 1080 1440 Time (mins) 0 0.5 1 1.5 2 CL Concentration (mg/L) J1 J5 J6 J8 (b) 0 360 720 1080 1440 Time (mins) 0 0.25 0.5 0.75 1 Multiplier Demand Pattern (c) Fig. 15: (a) MPC control action at Junctions 1 and 6 of Net1, and the corresponding chlorine concentrations at these junctions and Junctions 5 and 8 under (c) patterned demand at Junction 1. is set to intrude the system at Junction 6 with a concentration of 0.3 mg/L mid-day . Additionally , chlorine concentrations are limited to 1.2 mg/L for cost reasons. for this case study , we introduce two types of disturbance to the system: a sudden drop in chlorine concentration at Junction 6 to 0.15 mg/L at the 12th hour of the day and a sudden increase to 2 mg/L at the 18th hour . Fig. 15a shows the control action at Junctions 1 and 6, while Fig. 15b demonstrates the corresponding chlorine concentrations at these Junctions and Junctions 5 and 8. For Junction 1, the control action is higher and almost constant at 1 . 9 × 10 4 mg/min as the junction is located at the very start of the network and all the do wnstream elements are affected by its input. On the other hand, the booster station at Junction 6 acts on the disturbances and the changes at the downstream nodes effecti vely . Results validate the performance of the control algorithm and how adaptiv e it is under these disturbances. The run-time recorded for applying the control algorithm for this case study is 278 seconds. Likewise, chlorine concentrations are regulated through the FFCL-1 network as Fig. 16 exhibits. The total number of states of the original model is 10356, while the reduced model has 200 states. W ater quality time-step, control horizon, and simulation period are same as the previous case study . T wo different ﬁctitious reactants are assumed to be detected, the ﬁrst one at J76 with initial concentration of 0.3 mg/L while the second one at J89 of 0.2 mg/L. Control actions illustrated in Fig. 16a are under the condition of hydraulic proﬁle that results in changing ﬂow directions. Y et, the control algorithm recovers effecti vely and maintain chlorine concentrations within the desired range. In short, the R OMs-based control algorithms guarantee the bounds deﬁned for the inputs and outputs while being tractable for larger networks. V I . C O N C L U S I O N , P A P E R ’ S L I M I TA T I O N S , A N D R E C O M M E N D A T I O N S F O R F U T U R E W O R K Relying on the results from the numerical case studies in Section ” Case Studies ”, we answer the research questions posed: ▷ A1: The multi-species water quality model can be effecti vely linearized around operating points updated ev ery speciﬁc moving window according to the hydraulic proﬁle, instantaneous changes, initial conditions, and control actions. Howe ver , to achieve the desired accuracy this window is reduced and accordingly , the computational time increases. ▷ A2: The presented MOR methods yield high accuracy in estimating output concentrations for both chlorine and the ﬁctitious substance in the system. The three MOR procedures: LPOD, BPOD, and NLPOD are able to handle non- 23 0 360 720 1080 1440 Time (mins) 0 1 2 3 4 5 6 Control Action (mg/min) # 10 4 J89 (a) 0 360 720 1080 1440 Time (mins) 0 0.5 1 1.5 2 2.5 CL Concentration (mg/L) J56 J67 (b) Fig. 16: (a) Control action at J89 of the FFCL-1 network, (b) the corresponding chlorine concentrations at J56 and J67. zero initial conditions by fav oring the control actuators’ inputs while building the ofﬂine snapshots. Additionally , the NLPOD method requires more computational time to handle and interpolate the nonlinearity in the system, yet, it is still computationally tractable, same for LPOD and BPOD. ▷ A3: MPC’ s behavior depends on the underlying model and its accuracy . Accordingly , the Implicit Upwind scheme is preferred over the Explicit Upwind scheme because of its ability to provide highly-accurate simulation for the full and reduced-order MS-WQM. Moreover , numerical case studies sho w that the three MOR producers are rob ust to dynamically changing hydraulic proﬁles. ▷ A4: MPC shows robustness and high ﬂexibility in re gulating chlorine le vels in WDNs under dif ferent scenarios of contamination ev ents and hydraulic proﬁles by applying feedback control on the reduced order model while maintaining affordable computational requirements. Both proposed control procedures, the linearized model- and the relax ed nonlinear model-based show reliable performance while applying adaptive approaches according to the case study considered. These approaches lead to a different lev el of complexity and computational burden for each of the procedures which results in fa voring one procedure over the other according to the case study . Our study is not limitations-free. W e highlight these limitations next along with the authors’ future work to be in vestigated. First, this work used pre-assigned ﬁxed-location booster and sensor locations. Given that these locations impact performance, future work will include optimizing these locations from a control-theoretic perspectiv e. Second, additional approaches to model linearization should be explored to potentially exploit ofﬂine pre-computed FOM trajectories. Lastly , further work is needed to improve the relaxation method because we expect opportunities to further improve computational performance compared to the linearized model. On the other hand, this study is considered a computational study that is based on a model that has been veriﬁed, howe ver , real-time experimental study to verify the considered model and our framew ork performance under various scenarios is recommended. A. Data A vailability Statement All data, models, and codes that support the ﬁndings of this study are av ailable from the corresponding author upon reasonable request. B. Acknowledgments This work is partially supported by National Science Foundation under grants 1728629, 2015603, 2015671, 2151392, and 2015658. 24 R E F E R E N C E S [1] K. N. Palansooriya, Y . Y ang, Y . F . Tsang, B. Sarkar, D. Hou, X. Cao, E. Meers, J. Rinklebe, K.-H. Kim, and Y . S. Ok, “Occurrence of contaminants in drinking water sources and the potential of biochar for water quality improvement: A review , ” Critical Reviews in Envir onmental Science and T echnology , vol. 50, no. 6, pp. 549–611, 2020. [2] B. Moore, “Principal component analysis in linear systems: Controllability, observability , and model reduction, ” IEEE T ransactions on Automatic Control , vol. 26, no. 1, pp. 17–32, Feb. 1981. [3] U. Baur, P . Benner , and L. Feng, “Model Order Reduction for Linear and Nonlinear Systems: A System-Theoretic Perspective, ” Archives of Computational Methods in Engineering , vol. 21, no. 4, pp. 331–358, Dec. 2014. [4] L. Montier, T . Henneron, B. Goursaud, and S. Clenet, “Balanced proper orthogonal decomposition applied to magnetoquasi-static problems through a stabilization methodology , ” IEEE transactions on magnetics , vol. 53, no. 7, pp. 1–10, 2017. [5] J. Rutzmoser , “Model order reduction for nonlinear structural dynamics, ” Ph.D. dissertation, T echnische Universit ¨ at M ¨ unchen, 2018. [6] C. W . Rowley , “MODEL REDUCTION FOR FLUIDS, USING B ALANCED PROPER OR THOGONAL DECOMPOSITION, ” p. 17, 2005. [7] K. W illcox and J. Peraire, “Balanced model reduction via the proper orthogonal decomposition, ” AIAA journal , vol. 40, no. 11, pp. 2323–2330, 2002. [8] E. J. Grimme, Krylov pr ojection methods for model reduction . Uni versity of Illinois at Urbana-Champaign, 1997. [9] C. A. Beattie and S. Gugercin, “Interpolation theory for structure-preserving model reduction, ” in 2008 47th ieee conference on decision and contr ol . IEEE, 2008, pp. 4204–4208. [10] S. Gugercin, “ An iterativ e svd-krylov based method for model reduction of large-scale dynamical systems, ” Linear Algebra and its Applications , vol. 428, no. 8-9, pp. 1964–1986, 2008. [11] A. C. Antoulas, D. C. Sorensen, and S. Gugercin, “ A surve y of model reduction methods for large-scale systems, ” T ech. Rep., 2000. [12] R. Kumar and D. Ezhilarasi, “ A state-of-the-art survey of model order reduction techniques for large-scale coupled dynamical systems, ” International Journal of Dynamics and Contr ol , pp. 1–17, 2022. [13] S. Lall, J. E. Marsden, and S. Glava ˇ ski, “ A subspace approach to balanced truncation for model reduction of nonlinear control systems, ” International Journal of Robust and Nonlinear Contr ol: IF AC-Afﬁliated J ournal , vol. 12, no. 6, pp. 519–535, 2002. [14] S. Barrachina, P . Benner, E. S. Quintana-Ort ´ ı, and G. Quintana-Ort ´ ı, “Parallel algorithms for balanced truncation of large-scale unstable systems, ” in Pr oceedings of the 44th IEEE Conference on Decision and Contr ol . IEEE, 2005, pp. 2248–2253. [15] K. Zhou, G. Salomon, and E. Wu, “Balanced realization and model reduction for unstable systems, ” International Journal of Robust and Nonlinear Contr ol: IF A C-Afﬁliated Journal , vol. 9, no. 3, pp. 183–198, 1999. [16] L. Sirovich, “Turb ulence and the dynamics of coherent structures. i. coherent structures, ” Quarterly of applied mathematics , vol. 45, no. 3, pp. 561–571, 1987. [17] K. Willcox and J. Peraire, “Balanced Model Reduction via the Proper Orthogonal Decomposition, ” AIAA Journal , vol. 40, no. 11, pp. 2323–2330, Nov . 2002. [18] V . B. Nguyen, S. B. Q. Tran, S. A. Khan, J. Rong, and J. Lou, “Pod-deim model order reduction technique for model predicti ve control in continuous chemical processing, ” Computers & Chemical Engineering , vol. 133, p. 106638, 2020. [19] D. Galbally , K. Fidkowski, K. W illcox, and O. Ghattas, “Non-linear model reduction for uncertainty quantiﬁcation in large-scale inverse problems, ” International journal for numerical methods in engineering , vol. 81, no. 12, pp. 1581–1608, 2010. [20] N. Akkari, F . Casenave, and D. Ryckelynck, “ A novel gappy reduced order method to capture non-parameterized geometrical variation in ﬂuid dynamics problems, ” 2019. [21] K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem, “The gnat method for nonlinear model reduction: effecti ve implementation and application to computational ﬂuid dynamics and turbulent ﬂows, ” Journal of Computational Physics , vol. 242, pp. 623–647, 2013. [22] P . Benner, S. Gugercin, and K. Willcox, “ A survey of projection-based model reduction methods for parametric dynamical systems, ” SIAM revie w , vol. 57, no. 4, pp. 483–531, 2015. [23] B. Ulanicki, A. Zehnpfund, and F . Martinez, “Simpliﬁcation of water distribution network models, ” in Pr oc., 2nd Int. Conf. on Hydroinformatics . Balkema Rotterdam, Netherlands, 1996, pp. 493–500. [24] F . Mart ´ ınez Alzamora, B. Ulanicki, and E. Salomons, “Fast and practical method for model reduction of large-scale water-distribution networks, ” Journal of W ater Resources Planning and Management , vol. 140, no. 4, pp. 444–456, 2014. [25] U. Shamir and E. Salomons, “Optimal Real-T ime Operation of Urban W ater Distribution Systems Using Reduced Models, ” Journal of W ater Resour ces Planning and Manag ement , vol. 134, no. 2, pp. 181–185, Mar . 2008. [26] A. Preis, A. J. Whittle, A. Ostfeld, and L. Perelman, “Efﬁcient hydraulic state estimation technique using reduced models of urban water networks, ” Journal of W ater Resour ces Planning and Management , vol. 137, no. 4, pp. 343–351, 2011. [27] L. Perelman and A. Ostfeld, “W ater Distribution System Aggregation for W ater Quality Analysis, ” Journal of W ater Resources Planning and Management , vol. 134, no. 3, pp. 303–309, May 2008. [28] A. Elkhashap and D. Abel, “Model Order Reduction of The Time-Dependent Adv ection-Diffusion-Reaction Equation with Time-V arying Coefﬁcients: Application to Real-Time W ater Quality Monitoring, ” in 2022 European Control Confer ence (ECC) . London, United Kingdom: IEEE, Jul. 2022, pp. 333–338. [29] Z. Li, S. G. Buchberger , and V . Tzatchkov , “Importance of dispersion in network water quality modeling, ” in Impacts of Global Climate Change , 2005, pp. 1–12. [30] F . Shang, H. W oo, J. B. Burkhardt, and R. Murray , “Lagrangian method to model advection-dispersion-reaction transport in drinking water pipe networks, ” Journal of W ater Resour ces Planning and Management , vol. 147, no. 9, p. 04021057, 2021. [31] S. W ang, A. F . T aha, A. Chakrabarty , L. Sela, and A. A. Abokifa, “Model order reduction for water quality dynamics, ” W ater Resour ces Researc h , vol. 58, no. 4, p. e2021WR029856, 2022. [32] B. Fu, J. S. Horsbur gh, A. J. Jakeman, C. Gualtieri, T . Arnold, L. Marshall, T . R. Green, N. W . T . Quinn, M. V olk, R. J. Hunt, L. V ezzaro, B. F . W . Croke, J. D. Jakeman, V . Snow , and B. Rashleigh, “Modeling W ater Quality in W atersheds: From Here to the Next Generation, ” W ater Resources Research , vol. 56, no. 11, Nov . 2020. [33] T . Lassila, A. Manzoni, A. Quarteroni, and G. Rozza, “Model Order Reduction in Fluid Dynamics: Challenges and Perspectives, ” in Reduced Or der Methods for Modeling and Computational Reduction , A. Quarteroni and G. Rozza, Eds. Cham: Springer International Publishing, 2014, pp. 235–273. [34] J. He and L. J. Durlofsky , “Reduced-order modeling for compositional simulation by use of trajectory piecewise linearization, ” SPE Journal , vol. 19, no. 05, pp. 858–872, 2014. [35] Z. Ohar and A. Ostfeld, “Optimal design and operation of booster chlorination stations layout in water distrib ution systems, ” W ater r esear ch , vol. 58, pp. 209–220, 2014. [36] A. Ostfeld and E. Salomons, “Conjunctive optimal scheduling of pumping and booster chlorine injections in water distrib ution systems, ” Engineering optimization , vol. 38, no. 03, pp. 337–352, 2006. [37] G. Munav alli and M. M. Kumar , “Optimal scheduling of multiple chlorine sources in water distrib ution systems, ” Journal of water r esources planning and management , vol. 129, no. 6, pp. 493–504, 2003. [38] I. Fisher, G. Kastl, F . Shang, and A. Sathasivan, “Frame work for optimizing chlorine and byproduct concentrations in drinking water distribution systems. ” Journal: American W ater W orks Association , vol. 110, no. 11, 2018. 25 [39] S. W ang, A. F . T aha, and A. A. Abokifa, “How effecti ve is model predictive control in real-time water quality regulation? state-space modeling and scalable control, ” W ater Resources Resear ch , vol. 57, no. 5, p. e2020WR027771, 2021. [40] S. W ang, A. Chakrabarty , and A. F . T aha, “Data-driv en identiﬁcation of dynamic quality models in drinking water networks, ” Journal of W ater Resour ces Planning and Manag ement , vol. 149, no. 4, p. 04023008, 2023. [Online]. A vailable: https://ascelibrary .org/doi/abs/10.1061/JWRMD5.WRENG- 5431 [41] S. M. Elsherif, S. W ang, A. F . T aha, L. Sela, M. H. Giacomoni, and A. A. Abokifa, “Control-theoretic modeling of multi-species water quality dynamics in drinking water networks: Surve y, methods, and test cases, ” Annual Reviews in Contr ol , 2022. [42] L. Rossman, H. W oo, M. Tryby , F . Shang, R. Janke, and T . Haxton, “Epanet 2.2 user manual water infrastructure division, center for en vironmental solutions and emer gency response, ” US En vironmental Pr otection Agency , Cincinnati, Ohio , 2020. [43] F . Shang, L. A. Rossman, and J. G. Uber , “Epanet-msx 2.0 user manual, ” US Envir onmental Pr otection Agency , Cincinnati, Ohio, EP A/600/R-22/199 , 2023. [44] C. Hirsch, “Numerical computation of internal and external ﬂows. vol. 2-computational methods for inviscid and viscous ﬂo ws, ” Chichester , 1990. [45] T . M. Apostol, Calculus, V olume 1 . John W iley & Sons, 1991. [46] S. Lall, J. E. Marsden, and S. Glava ˇ ski, “Empirical model reduction of controlled nonlinear systems, ” IF AC Proceedings V olumes , vol. 32, no. 2, pp. 2598–2603, 1999. [47] K. Diao and W . Rauch, “Controllability analysis as a pre-selection method for sensor placement in water distribution systems, ” W ater resear ch , vol. 47, no. 16, pp. 6097–6108, 2013. [48] D. M. Luchtenb urg and C. W . Rowley , “Model reduction using snapshot-based realizations, ” Bulletin of the American Physical Society , vol. 56, no. 18, 2011. [49] Y . J. Y ang, R. C. Haught, and J. A. Goodrich, “Real-time contaminant detection and classiﬁcation in a drinking water pipe using con ventional water quality sensors: T echniques and experimental results, ” Journal of en vironmental management , vol. 90, no. 8, pp. 2494–2506, 2009. [50] T . Acrylamide, “National primary drinking water regulations, ” EP A, Kidney , vol. 2, no. 4-D, pp. 0–07. [51] G. P . McCormick, “Computability of global solutions to factorable nonconv ex programs: Part i—conve x underestimating problems, ” Mathematical pr ogramming , vol. 10, no. 1, pp. 147–175, 1976. [52] Y . Chen, “Model order reduction for nonlinear systems, ” Ph.D. dissertation, Massachusetts Institute of T echnology , 1999. [53] W . H. A. Schilders, H. A. van der V orst, J. Rommes, H.-G. Bock, F . de Hoog, A. Friedman, A. Gupta, H. Neunzert, W . R. Pulleyblank, T . Rusten, F . Santosa, A.-K. T ornberg, L. L. Bonilla, R. Mattheij, and O. Scherzer, Eds., Model Order Reduction: Theory, Resear ch Aspects and Applications , ser. Mathematics in Industry. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, vol. 13. [54] Y .-Y . Liu and A.-L. Barab ´ asi, “Control principles of complex systems, ” Revie ws of Modern Physics , vol. 88, no. 3, p. 035006, 2016. [55] A. G. Seyoum and T . T . T anyimboh, “Integration of hydraulic and water quality modelling in distribution networks: Epanet-pmx, ” W ater Resour ces Management , vol. 31, no. 14, pp. 4485–4503, 2017. [56] F . Mazzoni, S. Alvisi, M. Blokker , S. G. Buchberger , A. Castelletti, A. Cominola, M.-P . Gross, H. E. Jacobs, P . Mayer , D. B. Steffelbauer et al. , “In vestigating the characteristics of residential end uses of water: a worldwide revie w , ” W ater Researc h , p. 119500, 2022.

Comprehensive Framework for Controlling Nonlinear Multi-Species Water Quality Dynamics

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment