Model Order Reduction for Water Quality Dynamics
Shen Wang, Ahmad F. Taha, Ankush Chakrabarty, Lina Sela, Ahmed Abokifa
11 Model Order Reduction for Water Quality Dynamics
Shen Wang † , Ahmad F. Taha † , ∗ , Ankush Chakrabarty ‡ , Lina Sela ♠ , and Ahmed Abokifa (cid:5) Abstract —A state-space representation of water quality dy-namics describing disinfectant (e.g., chlorine) transport dynam-ics in drinking water distribution networks has been recentlyproposed. Such representation is a byproduct of space- andtime-discretization of the PDE modeling transport dynamics.This results in a large state-space dimension even for smallnetworks with tens of nodes. Although such a state-space modelprovides a model-driven approach to predict water qualitydynamics, incorporating it into model-based control algorithmsor state estimators for large networks is challenging and attimes intractable. To that end, this paper investigates modelorder reduction (MOR) methods for water quality dynamics.The presented investigation focuses on reducing state-dimensionby orders of magnitude, the stability of the MOR methods, andthe application of these methods to model predictive control.
Index Terms —Water distribution network, water qualitymodel, state-space representation, model order reduction, modelpredictive control.
I. I
NTRODUCTION AND P APER C ONTRIBUTIONS W ater quality models can be expressed in the form of astate-space representation [1] representing the numeri-cal solution of the PDEs describing the spatiotemporal evolutionof disinfectant concentration (e.g., chlorine) in various networkelements in drinking water distribution networks (WDN). Un-fortunately, even for small-to-midsize networks, the dimension(or order) of the state-space models can reach or , dueto ensuring high fidelity in model predictions.This is due to discretization methods used for solving thepartial differential equations (PDE) that describe the chlorinetransport in pipes. Very high resolution in discretization yieldsimproved accuracy in the control-oriented model but resultsin high state-space dimension, leading to heavy computationalburden during simulation and difficulty for use in controllerdesign [2]. This implies that high-accuracy models, thougheffective for predicting system dynamics, are not amenableto controller or estimator design for large-scale networks—especially in the presence of state or input constraints. To thatend, model order reduction (MOR) is necessary to derive acompact model for fast simulation and efficient synthesis ofcontrollers and state estimators. † Department of Electrical and Computer Engineering, The University ofTexas at San Antonio, TX 78249. Emails: { shen.wang,ahmad.taha } @utsa.edu. ‡ Research Scientist, Mitsubishi Electric Research Laboratories, Cambridge,MA, USA. Email: [email protected] ♠ Department of Civil, Architecture, and Environmental Engineering, TheUniversity of Texas at Austin. Email: [email protected] (cid:5)
Department of Civil, Materials, and Environmental Engineering, TheUniversity of Illinois Chicago. Email: [email protected] ∗ Corresponding author.This work is partially supported by National Science Foundation undergrants 1728629, 2015603, and 2015671. A full-order model is referred to as an exact or near-exactmodel to describe water quality dynamics and usually is highin dimension. The motivation of MOR is to reduce the full-order model to a reduced-order model that has a much smallernumber of states or order without significantly decreasing modelaccuracy while maintaining input-output relationships and re-taining certain properties of the system such as controllabilityand observability. To this end, the high-dimensional state vectorfrom the full-order model is projected onto a lower-dimensionalsubspace within which the model accuracy and certain proper-ties are not significantly compromised. The MOR method hasbeen widely applied in computational fluid dynamics [3]–[5],and other fields such as computational electromagnetics, as wellas in micro- and nano-electro-mechanical systems design [6].However, reduced-order models have never been explored formodeling water quality dynamics—a gap that is filled in thispaper. A. Literature review of MOR algorithms
Various MOR algorithms have been developed in the richliterature of dynamic systems. The majority of these algorithmscan be divided into two categories: singular value decom-position (SVD) based methods [3], [5], [7]–[10] and Krylovsubspace methods or moment matching based methods [11]–[15]. There are also some studies combining SVD and Krylovmethods [16]. In general, methods based on Krylov subspacesdo not preserve important properties of the original system suchas stability and passivity [4]. Hence, we are interested in SVD-based methods in which several methods can guarantee stabilityor controller designed related properties. Moreover, SVD-basedapproaches can be further divided into subcategories, for exam-ple, the methods using Hankel operators concepts [7], [8], thebalanced realization theory [3], [9], [17], and proper orthogonaldecomposition (POD) based methods [3], [5], [10], [18], [19].This paper focuses on seeking proper SVD-based model reduc-tion methods for the water quality dynamics represented lineardiscrete-time systems with a state-space presentation.Moore [9] proposes the balanced truncation (BT) methodwhich considers or balances controllability and observabilityGramians. Unfortunately, the BT method becomes computa-tionally intractable for large-scale systems (e.g., statesor more). The POD method proposed by Sirovich [10] istractable, but not as accurate as of the BT. To compute balancingtransformations for high-order systems, Willcox [5] proposes atechnique that combines the POD and concepts from balancedrealization theory. However, when the number of outputs of asystem is large, this method is intractable, and it has the riskof truncating states that are poorly observable yet very stronglycontrollable. As a result, the balanced POD (BPOD) [3] is pro-posed to address the same issues, and it combines the advantages a r X i v : . [ ee ss . S Y ] F e b f BT and POD methods successfully. That is, BPOD balancesboth controllability and observability and is tractable for a large-scale system, which results in the computational cost similar toPOD and the accuracy similar to BT.In the water systems research community, the reduced modelsare applied in simplifying networks in the context of hydraulicsand the techniques adopted are different from the ones afore-mentioned in dynamic systems. That is, the reduced modelachieves a high-fidelity representation for the complete networkhydraulic, but greatly simplifies the computation. For example,a skeletonization method [20] removes nodes from a linearizednetwork model through Gauss elimination, and reduces thenetwork model while preserving the nonlinear characteristicsof the original network model. Salomons [21] obtain a reducedmodel by the genetic algorithm. To meet the computationalefficiency requirements of real-time hydraulic state estimation,Preis et al. [22] introduce a reduced model using a watersystem-aggregation technique. A methodology and applicationof a conjunctive hydraulic and water quality model for waterdistribution systems aggregation are presented in study [23]. Incontrast, this paper focuses on a different research problem thatreduces the order of water quality dynamics model, which hasnot been studied in the literature. B. Paper contributions
The major objective of the paper is to investigate theperformance of water quality MOR algorithms, address theirtheoretical limitations, assess their computational performance,and evaluate their potential when applied within a model-basedpredictive control framework. The below list outlines in detailthe contributions of the paper.1) We present the first attempt to identify reduced-ordermodels for water quality simulation with theoretical anal-ysis: our case studies show that the reduced-order modelproduces nearly identical water quality simulations in com-parison with the full-order models. The MOR algorithmsare tested for different scales of networks in sense ofaccuracy, low computational burden, robustness to initialconditions, and potential to be obtained as a stable reduced-order model.2) Classic MOR algorithms cannot ensure the stability of areduced-order model, despite the stability of the underlyingwater quality dynamics. Thus, we propose two methodsthat guarantee that the reduced-order model is stable. Thefirst method performs a stabilization process based on thestandard MOR procedure while trying to maintain accu-racy, and the second method simply adjusts a parameter inthe standard MOR (i.e., BPOD) procedure. Moreover, theparameter adjustment has a clear physical interpretation forwater quality dynamics, which is the largest travel time-step from booster stations to sensors.3) The obtained stable reduced-order models are then utilizedas predictive models in MPC for a WDN. Our case studyshows that the application of MPC to the reduced-ordermodel produces similar control effects compared with theapplication of MPC to the full-order model, while thecomputational burden is significantly reduced. (a) (b)
Fig. 1. (a) exemplar topology of a WDN (multi-input and multi-outputsystem with two boosters (inputs) and one sensor (output) to inject and detectchlorine), and P is split into three segments according to L-W scheme forillustration purpose; (b) demand profiles for all junctions during 24 hours. The rest of the paper is organized as follows. Section II brieflydescribes control-oriented water quality modeling, then presentsthe full-order model in state-space form. The principle of MORfor the specific water quality modeling is given at the begin-ning of Section III, followed by the fundamental proceduresof various MOR methods and the corresponding discussionsin Section IV. Section V presents two methods to obtain astable reduced-order model by adjusting the standard MORprocedures. Section VI presents case studies that corroborateour proposed approach. Limitations of our methods and futureresearch directions are given in Section VII.
Paper’s Notation.
Boldface characters represent matrices andvectors: a is a scalar, a is a vector, and A is a matrix. Matrix I denotes an identity matrix, whereas m × n denotes a zeromatrix with size m -by- n . The notation R denotes the set of realnumbers, and the notations R n and R m × n denote a columnvector with n elements and an m -by- n matrix in R . For anyvector x ∈ R n , x (cid:62) is its transpose. For any two matrices A and B with same number of columns, the notation { A , B } denotes [ A (cid:62) B (cid:62) ] (cid:62) . The norm of A is (cid:107) A (cid:107) . The i -th singular valueof A is denoted by σ i ; a vector σ ( A ) stands for all the singularvalues of A ; the k -largest singular values of A is denoted byvector σ ↓ k ( A ) . The i -th eigenvalue of A is denoted by λ i ; σ ( A ) stands for all the eigenvalues of A ; the k -largest eigenvalues of A is expressed by vector λ ↓ k ( A ) .II. S TATE -S PACE W ATER Q UALITY M ODEL
In this section, we present the necessary background for thestate-space model of water quality dynamics.
A. Water quality dynamics
We model WDN by a directed graph G = ( N , L ) . The set N defines the nodes and is partitioned as N = J ∪ T ∪ R where J , T , and R are collection of junctions, tanks, and reservoirs.Let L ⊆ N × N be the set of links, and define the partition L = P ∪ M ∪ V , where P , M , and V represent the collectionof pipes, pumps, and valves.An example of a WDN graph is shown in Fig. 1a, and adisinfectant (such as chlorine) is added into the network witha proper mass and concentration by booster stations installedat the reservoir source R and the junction J to disinfect theWDN. Each junction in a WDN is a type of node that connectstwo links. They may consume water, that is, they have water2emands. During a specific period, demands from all junctionsform a demand profile (see Fig. 1b) which is assumed as apriori for water quality simulations in water research. Moreover,all components are designed to meet a certain threshold ofwater quantity and quality. For example, the pump P M isdesigned to maintain the water pressure for all users (junctions).The concentration sensors are installed at specific junctionsfor measurement purposes (see the sensor installed at J inFig. 1a), in this way the water system operators can obtain theconcentrations around the network.In the end, the water quality model can be used to representthe movement of all chemical and/or microbial species (con-taminant, disinfectants, DBPs, metals, etc.) within a WDN asthey traverse various components of the network. Specifically,the single-species interaction and dynamics of chlorine areconsidered in our paper, and it can be expressed using astate-space model. We first introduce the concepts and thecorresponding methods that form the state-space model next.In a WDN, the pipe model is represented by chlorinetransport in differential pipe lengths by advection in additionto its decay due to reactions. Mathematically, for any pipe,the 1-D advection-reaction equation is given by a PDE. Al-though the PDE has no closed-form solution, the Lax-Wendroffscheme [24] can be adopted to discretize and solve it numer-ically. For example, pipe P in Fig. 1a is split into threesegments for illustrative purpose according to the scheme. Inpractice, the number of segments for a pipe (denoted by s L ) isdecided by the flow velocity and the length of that pipe. After thediscretization for all pipes, we define a vector c P ( t ) to collect allsegments in all pipes at time t . Similarly, the concentrations inpumps and valves at time t are denoted by c M ( t ) and c V ( t ) . The c L ( t ) (cid:44) { c P ( t ) , c M ( t ) , c V ( t ) } is a defined vector collecting allconcentrations from links (i.e., pipes, pumps, and valves). Forall nodes, we define c N ( t ) = { c J ( t ) , c R ( t ) , c T ( t ) } to collectsthe concentration from all junctions, reservoirs, and tanks. B. State-space representation
In general, the number of junctions, reservoirs, tanks, pipes,pumps, and valves is denoted by n J , n R , n TK , n P , n M , and n V . Hence, the number of states in nodes and links are n N = n J + n R + n TK and n L = n P × s L + n M + n V . Therefore, c N ( t ) ∈ R n N and c L ( t ) ∈ R n L .Summarily, a water quality state-vector x defining the con-centrations of the disinfectant (chlorine) in the network at time t is denoted as x ( t ) (cid:44) { c N ( t ) , c L ( t ) } ∈ R n x and n x = n N + n L .After applying the principle of conservation of mass, we canrewrite all component models (i.e., water quality model) usinga discrete-time state-space model of the form x ( k + 1) = A ( k ) x ( k ) + B ( k ) u ( k ) (1a) y ( k ) = C ( k ) x ( k ) + D ( k ) u ( k ) , (1b)where x ( k ) ∈ R n x is the state vector, the vectors u ( k ) ∈ R n u and y ( k ) ∈ R n y are the system inputs and outputs (i.e., n u booster stations and n y concentration sensors), and A ( k ) , B ( k ) , C ( k ) , and D ( k ) are time-varying system matrices. Notethat state matrix A ( k ) stands for the pure system dynamic(without any boosters); the element in input matrix B ( k ) rep-resents locations of booster stations and corresponding injected flow rates that are both assumed as known; output matrix C ( k ) indicates where the sensors are installed; feed-forward matrix D ( k ) describes the relationship or path directly from inputs tooutputs, for example, some booster stations are equipped with abuilt-in concentration sensor (i.e., both the booster station andthe sensor are at the same location) to calculate the injection ofchlorine accurately.Note that these system matrices are usually updated fromhalf hour to several hours according to the changing demands indemand profiles (e.g., demand changes every 2 hours in Fig. 1b).That is, during a specific period in the demand profile (demandsat all junctions do not change), the system matrices remainthe same, and the discrete-time linear time-varying (DT-LTV)system (1) degenerates into a discrete-time linear time-invariant(DT-LTI) system, that is x ( k + 1) = Ax ( k ) + Bu ( k ) (2a) y ( k ) = Cx ( k ) + Du ( k ) . (2b)Let us denote a generic n x -th order DT-LTI system (2) asthe tuple ( A , B , C , D ) for brevity. In this way, the DT-LTVsystem (1) comprises different tuples (DT-LTI systems). Forexample, DT-LTV models for 24 hours of the WDN in Fig. 1acan be represented by 12 different DT-LTI systems since thedemand changes each two hours.It is worthwhile to note that (i) the water quality model (1)is available in our recent published work [1], and the model isstable; (ii) after testing, each pipe is required to be split intoat least 100 segments on average to guarantee the accuracyof discretizing PDEs. That is, the order n x can easily reach level for a mid-sized network with hundreds of pipes; (iii) all system matrices in this model are sparse, and the sparsityrate could be from to . or higher. Moreover, thechallenging MOR problem for a DT-LTV system is equivalentto solving multiple relatively easy MOR problems for DT-LTIsystems, and the detail is presented next.III. M ODEL O RDER R EDUCTION
This section provides an overview of MOR principles for DT-LTI systems. For each specific MOR method, the correspondingprocedure is summarized.
A. Principle of MOR for DT-LTI system
The objective of MOR is to find a reduced-order model oforder n r (cid:28) n x , given by x r ( k + 1) = A r x r ( k ) + B r u ( k ) (3a) ˆ y ( k ) = C r x r ( k ) + D r u ( k ) , (3b)where the new state vector x r ∈ R n r ; the system input andoutput vectors u and ˆ y retain the same dimension. From a sim-ulation perspective, the estimated output ˆ y ( k ) from the reduced-order model is compared to the output from the full-ordersystem y ( k ) to test the performance of the MOR algorithms.The general idea of the aforementioned reduction process isto find an invertible transformation matrix T that maps state x into another space state z . Specifically, x = T z , where thereduced-order model state z ∈ R n x and transformation matrix T ∈ R n x × n x . Moreover, matrix T is constructed based on3ome specific principles so that the elements in z are orderedaccording to their importance ; more on that in Section III-B.After substituting x = T z into (2), we obtain z ( k + 1) = A z z ( k ) + B z u ( k ) (4a) y z ( k ) = C z z ( k ) + D z u ( k ) , (4b)where system matrices A z = T − AT , B z = T − B , C z = CT , and D z = D .However, the order of transformed system ( A z , B z , C z , D z ) is still n x . To obtain the reduced-order model (3), only the first n r important elements in z are kept and denoted by x r . That is, x can be replaced by T r x r while the system properties remainthe same, where T r ∈ R n x × n r is the first n r columns of T .Finally, Equation (4) turns into reduced-order model (3) where A r = S r AT r , B r = S r B , C r = CT r , D r = D . (5)Note that S r in (5) is neither the inverse or pseudo inverse of T r , it is simply the first n r row of T − .In short, MOR is achieved by defining a certain subspacewithin the original space and transforming the system dynam-ics. A question that remains unanswered is which reduced-order models are most beneficial for achieving our goal: toreduce the difficulty of designing control/estimation algorithmssuch as Kalman filter, model predictive control (MPC), orlinear–quadratic regulator (LQR) for water quality regulation.Since these control algorithms are related to the concepts ofcontrollability and observability, one could argue that furtheringour understanding of the controllable and observable subspacescould prove most useful for controller design.Controllability describes the ability of an external input u to drive system state x ( k ) from any initial state x (0) to anyother final state x ( k f ) in a finite time interval [0 , k f ] [25]. Aquantitative metric for controllability in a DT-LTI system is the controllability Gramian defined as W C = (cid:80) ∞ m =0 A m BB (cid:62) ( A (cid:62) ) m , (6)which is the solution of Lyapunov equation W C − AW C A (cid:62) = BB (cid:62) . Similarly, observability describes the ability of recon-structing the initial unknown state x (0) in finite k f steps fromthe knowledge of output y ( k ) [25]. The corresponding metricfor observability is the observability Gramian defined as W O = (cid:80) ∞ m =0 ( A (cid:62) ) m C (cid:62) CA m , (7)and it is the solution of Lyapunov equation W O − A (cid:62) W O A = C (cid:62) C .The physical meaning of MOR in terms of using control-lable and observable subspaces are that the uncontrollable (notaffected by the input u ) and unobservable (does not affect theoutput y ) parts of the system without significantly affectinginput-output relations are discarded. The methods consideringcontrollability or observability Gramian are also referred to as Gramian based MOR , and when both Gramians are consideredin a method, it is then named as cross-Gramian based MOR [4],[26]. Next, we present several specific MOR algorithms thatconsider the Gramians which are suitable for water qualityapplication.
B. Description of common MOR algorithms
Herein, we describe three reduced-order modeling algo-rithms: balanced truncation (BT), proper orthogonal decompo-sition (POD), and balanced POD (BPOD).
1) Balanced Truncation:
Balanced trunction computes aspecific T that transforms the coordinate system x into anew coordinate system z (see (4)) in which the controllabilityand observability Gramians (denoted W Cz and W Oz , respec-tively) are diagonal and identical [9]. We denote this particularGramian by Σ z . Mathematically, Σ z = W Cz = W Oz , whichimplies that Σ z = W Oz W Cz , where W Cz = (cid:80) ∞ m =0 A mz B z B (cid:62) z ( A (cid:62) z ) m = T − W C T − (cid:62) , and W Oz = T (cid:62) W O T .Note that W C and W O are the Gramians of the full-ordermodel (2) before transformation. Hence, Σ z = W Oz W Cz = T − W O W C T ⇔ W O W C T = T Σ z , which indicates the transformation T can be obtained byfinding the eigenvector matrix of W O W C (also known as crossGramian matrix W X [27]), and the Σ z is just the correspondingdiagonal eigenvalue matrix. The final T r can be obtained simplyby selecting the first n r columns of T .Furthermore, the diagonal element σ i , i = 1 , . . . , n x in Σ z measuring controllability and observability simultaneously arethe Hankel singular values (HSVs) [28] of W X . Each singleHSV σ i = √ λ i provides a measure of “energy” for each statein a system. The model order reduction procedure provided byBT can also be viewed as retaining the n r − largest HSVs (i.e., σ ↓ n r ( W X ) ) or eigenvalues (i.e., λ ↓ n r ( W X ) ) with higher energyand discarding the least important ones with lower energy. The n r can be either specified as a fixed number directly or solved byspecified energy level defined by energy = (cid:80) n r j =1 σ j / (cid:80) n x i =1 σ i .The simplified procedure of BT approach is summarized asProcedure 1. Procedure 1:
Classical BT Form cross Gramian matrix W X = W O W C after solvingtwo Lyapunov equations Find the eigenvector matrix T of W X Specify n r as needed or obtain n r via setting an energy level Extract T r from the first n r columns of T Extract S r from the first n r rows of T − Obtain reduced-order model (3) by (5)
This BT method considers and balances the controllabilityand observability simultaneously with high accuracy. Besidesthat, it preserves stability, which is a necessary property for thereduced-order system. However, it requires the computation of W O and W C or equally solving two Lyapunov equations andeigenvalue decomposition that are intractable for a large-scalenetwork.
2) Proper Orthogonal Decomposition:
The POD algo-rithm [10], [18], which is also known as principal componentanalysis or the Karhunen-Lo`eve expansion, focuses on seekinga transformation T r such that a set of given data x ( t ) ∈ R n x ,with time t ∈ [0 , m ] is projected into the subspace of fixeddimension n r as x r = T r x ∈ R n r while minimizing the total4rror (cid:80) m − t =0 (cid:107) x ( t ) − x r ( t ) (cid:107) . To solve this problem, the datais assembled into an n x × m matrix X m = (cid:2) x (0) x (1) . . . x ( m − (cid:3) , (8)which is called snapshot . Define matrix W Cm as X m X (cid:62) m ∈ R n x × n x . Matrix W Cm is controlability Gramian-like matrixcontaining the first m items when calculating W C , and itis an approximation of W C . We refer to W Cm as the m -step controllability Gramian . The optimal transformation T r in POD method can be obtained by eigenvalue decomposition W Cm T r = T r Λ . Each eigenvector or the column vector in T r is named as POD mode [3], and the eigenvalues of W Cm arethe diagonal elements of Λ .Note that a snapshot is given by (8) which, if A is sparse, canbe easily obtained even with n x ≥ since the computationaltask involves sparse matrix multiplication. However, finding theeigenvectors of an n x × n x matrix W Cm is still challenging.Sirovich [10] solves this problem by defining (cid:102) W Cm = X (cid:62) m X m which is in R m × m instead of R n x × n x with m (cid:28) n x . Thematrix form of eigenvector and eigenvalue of this new (cid:102) W Cm are denoted as U and Λ . That is, the n x × n x eigenvalue problemis solved thorough constructing an m × m eigenvalue problem.Finally, the transformation matrices T r and S r can be obtained,and the detail is in Procedure 2. Note that when n x (cid:28) m , theabove step is not necessary, and computing W Cm directly iseasier and preferred. Procedure 2:
Classical POD Construct snapshot X m using (8) Perform eigenvalue decomposition for (cid:102) W Cm = X (cid:62) m X m byfinding U and Λ via (cid:102) W Cm U = U Λ Specify n r as needed or obtain n r via setting an energy level Extract transformation T r from the first n r column of X m U Λ − Extract transformation S r from the first n r rows of Λ − U X (cid:62) m Obtain reduced-order model (3) by (5)
Since the eigenvalues of AB are the same as those of BA (see [29, Theorem 1.3.22]), i.e., λ ( AB ) = λ ( BA ) ,the (cid:102) W Cm shares the same eigenvalues with m -step control-lability Gramian W Cm . That is, POD, not as the BT methodthat exploits the controllability Gramian (summations of in-finity items), uses only the first m items of the controllabil-ity Gramian. This makes POD tractable for even large-scalesystems but less accurate compared to BT. Thus, POD onlyconsiders controllability Gramian and selects the n r modeswith high-energy states (i.e., λ ↓ n r ( W Cm ) ), and discards the low-energy states.
3) Balanced Proper Orthogonal Decomposition:
BPODaims to combine the advantages of BT and POD together, thatis, to obtain an approximation to BT that is computationallytractable for large-scale systems [3]. Different from the snapshot X m (8) defined in POD method, BPOD defines X m as X m = (cid:2) x (0) . . . x ( m − . . . x n u (0) . . . x n u ( m − (cid:3) , (9)where the x i ( m ) is the impulse response (when u i ( k ) = δ ( k ) is an impulse signal). That is, x i ( m ) = A m b i and b i is the i -thcolumn of B . Moreover, BPOD denotes Y m as the snapshot of the dualsystem of system ( A , B , C , D ) or equally (2). That is, Y m isthe snapshot of system ( A (cid:62) , C (cid:62) , B (cid:62) , D (cid:62) ) with system statedenoted as z , and Y m can be expressed as Y m = (cid:2) z (0) . . . z ( m − . . . z n y (0) . . . z n y ( m − (cid:3) , (10)where the z i ( m ) is the impulse response, that is, z i ( m ) =( A (cid:62) ) m c (cid:62) i and c (cid:62) i is the i -th column of C (cid:62) .Then, the balancing modes are computed by forming thesingular value decomposition of the block Hankel matrix H m : H m = Y (cid:62) m X m = U Σ V (cid:62) = (cid:2) U r (cid:3) (cid:20) Σ r
00 0 (cid:21) (cid:20) V (cid:62) r (cid:21) , (11)where the diagonal element in Σ r ∈ R n r × n r collects thelargest n r singular values of H m , and U r and V (cid:62) r are thecorrepsonding left- and right-singular vectors. Hence, the finaltransformation matrices are T r = X m V r Σ − r and S r = Σ − r U (cid:62) r Y (cid:62) m ; see Procedure 3 for details. Procedure 3:
Classical BPOD Form snapshots X m (9) and Y m (10) of BPOD Perform singular value decomposition of H m = Y (cid:62) m X m (11), and obtain U r , Σ r , V r Specify n r as needed or obtain n r via setting an energy level Calculate T r = X m V r Σ − r and S r = Σ − r U (cid:62) r Y (cid:62) m Obtain reduced-order model (3) by (5)
From the above procedures, it is clear that BPOD considersboth controllability and observability Gramians like BT method,successfully avoids solving the two Lyapunov equations, form-ing cross Gramian matrix W X , and finding the eigenvaluesof W X to form a traceable algorithm simply by adopting theconcept of the snapshot.IV. T HOERETICAL C OMPARISONS OF
MOR A
LGORITHMSFOR W ATER Q UALITY D YNAMICS
The advantages and disadvantages from the aspects oftractability for large-scale systems, relative accuracy, and sta-bility preserving properties are summarized in this section. Sub-sequently, we present a comparison of these MOR algorithmsto determine which properties of these algorithms fit best forour application— the water quality model.First, the complexity of solving Lyapunov equations [31],eigenvalue decomposition [32], and SVD [33] are roughlyestimated by O ( n ) for an n × n matrix. Although all BT,POD, and BPOD procedures include at least one of the abovetime-consuming computation, the matrix size appeared in theseprocedures and the number of executions of the procedures aredifferent; see Tab. I. For example, BT solves two Lyapunovequations ( n x × n x ) and performs eigenvalue decompositionfor W X once, that is, it performs algorithm with O ( m ) threetimes; POD performs eigenvalue decomposition for (cid:102) W Cm ,the complexity is roughly O ( m ) ; BPOD performs truncatedSVD (rank n r ; see (11)) for H m , the complexity is roughly O ( m n y n u n r ) . Due to m (cid:28) n x , it is clear that BT is theslowest algorithm. The computational load of POD and BPOD5 ab. IC OMPARISONS AMONG THREE CLASSICAL
MOR
METHODS . Matrixinvolved Dimension Decomposi-tion method † Moderetained Complexity(Computational load) Tractability forlarge-scale systems Relativeaccuracy (cid:5)
Stability ∗ preservingBT W X R n x × n x ED λ ↓ n r ( W X ) O ( n x ) (High) No High YesPOD ‡ (cid:102) W Cm R m × m ED λ ↓ n r ( (cid:102) W Cm ) O ( m ) (Low) Yes Low NoBPOD H m R mn y × mn u SVD σ ↓ n r ( H m ) O ( m n y n u n r ) (Low) Yes Medium Depends ♠† Eigenvalue decomposition (ED); singular value decomposition (SVD); ‡ use (cid:102) W Cm when m (cid:28) n x , otherwise, use W Cm instead. (cid:5) Under the same number of modes; ∗ asymptotically stable; ♠ Rowley [3] conjectures that BPOD can preserve stability via theoretical analysiswhile other studies [6], [30] denied this assertion based on test results. algorithms depends on the value of m and n y n u n r . When thesevalues are of the same order of magnitude, the complexity ofPOD and BPOD are similar and also tractable for large-scalesystems; see Tab. I for the final computational load.Second, a useful property of the BT method is that the errorbounds provided are close to the lower bound achieved byany reduced-order model, and the upper bound of the error isguaranteed less than (cid:80) n x i = n r +1 σ i , where σ i is the i -th HSV of W X in Section III-B1 [3]. Furthermore, BT considers W X thatbalances controllability and observability. Hence, the relativeerror in terms of preserving input-output relations of the BTmethod is theoretically the best among these three methods.As for POD, it only takes advantage of m -step controllabilityGramian W Cm , and has the risk of discarding the states that arehighly observable but less controllable. This defect leads to theinaccuracy of POD or even the failure to reduce system order.In fact, BPOD provides a good approximation of BT sincethey retain the same states/modes when performing MOR. Thecorresponding n r -largest eigenvalues (singular values) of allmodes/states retained from three procedures are summarizedin Tab. I. That is, the states retained are those with λ ↓ n r ( W X ) ( σ ↓ n r ( H m ) ) for BT (BPOD). Suppose λ i ( σ i ) is the i -th elementof λ ↓ n r ( W X ) ( σ ↓ n r ( H m ) ), and λ i (cid:117) σ i always holds truewhen m is large enough. This implies that BPOD, like BT,uses the same information (i.e., cross Gramians W X ) to obtainthe transformation T r to reduce system order. The differencelies in the way each algorithm deals with such information.For example, BT requires an eigenvalue decomposition of W X while BPOD involves singular value decomposition to H m embedded in W X . From this perspective, the transformation T obtained from the BPOD method is the balanced one thatsimultaneously diagonalizes W C and W O . Hence, BPOD hasthe same ability (when m is large enough) ensure the accuracyas BT does.The BT method guarantees stability theoretically afterMOR [9], while POD fails [4], [34] due to the projection ituses. As for the BPOD method, different studies report variousresults. Rowley [3] claims that the stability-preserving fromBPOD is guaranteed for linear systems when snapshots arelarge enough. However, there is no guideline to find what isthe proper size of snapshots. We suspect that this might bewhy the results from several studies [6], [30] do not supportthe stability-preserving property of BPOD. In other words,stability-preserving of BPOD depends on the parameter (i.e.,length of snapshots) set. We would present an approach to find the proper parameter that ensures stability in Section V.In short, POD and/or BPOD methods are still preferred sincethey are tractable for a large-scale system, even if their stability-preserving properties are not sufficient for water quality dynam-ics MOR. Toward that goal, the stability-preserving model orderreduction approaches for water quality modeling are exploredand proposed next.V. S TABILIZING
MOR A
LGORITHMS
In this section, the brief literature of stability-preservingMOR is presented first followed by proposing two methodssuitable for the stable and discrete-time system.The approaches developing stability-preserving reduced-order model involve a priori frameworks [35]–[39], and aposteriori frameworks [40]–[43].
A priori stability-preservingframework can be difficult to implement, and/or requires thesolution of PDEs or special system structure. For example,the Galerkin projection (i.e. a congruence transform) needssign-definite matrices in the full-order system. Conversely, aposteriori methods minimally modify the reduced-order systemmatrix so that the accuracy is not significantly changed. Themore detailed literature can be found in [41]. Herein, we discussstability-preservation with both these frameworks.In short, we propose a novel posteriori method to stabilizethe reduced-order model after performing the POD or BPODprocedure. Comparing to existing methods such as [42], theproposed posteriori method adjusts the final transformation A r directly instead of optimizing/adjusting the projection U r slightly while fixing the projection V (cid:62) r in (11). Furthermore,published studies focus more on continuous systems while theproposed one is for discrete systems and uses a special trick torelax the non-convex optimization problem.Moreover, a priori stabilization is proposed specifically forthe BPOD procedure. To the best of our knowledge, this isthe first priori stabilization method for BPOD. Comparing toexisting studies such as [35] focusing on finding a special andcomplex projection for POD, the proposed method is simpleand intuitive since it only requires modifying a parameter inBPOD. A. A posteriori stabilization for POD-based approach
After performing the POD procedure, the stability of thereduced-order model can not be guaranteed. That is, A r isunstable. As for the stability-preserving from BPOD, it dependson the length of snapshots m . When m is not proper, A r from6POD tends to be unstable. The motivation behind this posteriormethod is to adjust A r by ∆ A r such that the adjusted systemis stable while ensuring that the norm (cid:107) ∆ A r (cid:107) is as small aspossible to maintain the accuracy of the full-order system.From Lyapunov stability conditions [44], we know that a DT-LTI system ( A , B , C , D ) is stable when there exist a positivedefinite matrix P such that A (cid:62) P A − P is negative definite.Thus, in our case the stabilization problem can be expressed asan optimization problem below. min P , ∆ A r ,α α (12a) subject to P (cid:31) , α > (12b) (cid:107) ∆ A r (cid:107) ≤ α (12c) ( A r + ∆ A r ) (cid:62) P ( A r + ∆ A r ) − P ≺ , (12d)where P , ∆ A r in R n r × n r and α ∈ R are the optimizationvariables. However, this problem (12) is nonconvex. We cantake the following steps to tranform it into an semidefiniteprogramming (SDP) that is convex.Since P (cid:31) is invertible, let P = P P − P , and takingSchur complements, we have(12d) ⇔ ( A r + ∆ A r ) (cid:62) P P − P ( A r + ∆ A r ) − P ≺ ⇔ (cid:20) − P ( P A r + ∆ Y ) (cid:62) P A r + ∆ Y − P (cid:21) (cid:31) , (13)where ∆ Y = P ∆ A r . Furthermore, the Cauchy-Schwarzinequality yields (cid:107) ∆ Y (cid:107) = (cid:107) P ∆ A r (cid:107) ≤ (cid:107) P (cid:107) (cid:107) ∆ A r (cid:107) . Invoking (12c) yields (cid:107) ∆ Y (cid:107) ≤ α (cid:107) P (cid:107) . Taking Schur com-plements yields(12c) ⇔ (cid:107) ∆ Y (cid:107) ≤ α y ⇔ (cid:20) I ∗ ∆ Y α y I (cid:21) (cid:23) , where α y = α (cid:107) P (cid:107) . However, the α y is still nonlinear and canbe relaxed as (cid:20) I ∗ ∆ Y α y I (cid:21) (cid:23) , α y > . (14)Finally, a convex SDP stabilizing the reduced-order model isformed as min P , ∆ Y ,α y α y (15a) subject to P (cid:31) , (13) , (14) . (15b)Note that the dimension of optimization variables such as P and ∆ Y depends on the size of A r in reduced-order model (3)instead of the size of A in full-order model (2). This indicatesthe formulation (15) can be applied to solve the stabilizationproblem even for a large-scale system, provided that n r is notlarge after reduction procedures. B. A priori stabilization for BPOD approach
The a priori stabilization method introduced in this sectionis different from the posterior one since it directly modifiesthe length of snapshots m in the standard BPOD Procedure 3while generating a stable reduced-order model. Ideally, whenparameter m is set as ∞ , the procedure then has all systeminformation as the cross-Gramian does in the BT method, andconsequently BPOD reaches its best performance. However, Fig. 2. Chlorine travel paths from the input location (booster station) todifferent output locations (Sensors 1 and 2). parameter m cannot be too large (or equal to infinity) due tothe prohibitive computational cost. Hence, selecting a lowerbound for m denoted by m is important. To the best of ourknowledge, this is the first time attempt to produce such a valueor condition that ensures the stability of a DT system usingBPOD method. In particular, we propose two methods to find m . The first method is generated from insights from transportdynamics in water networks. The second method is via a controltheoretic approach. In the case studies section, we compare theperformance of these two methods.
1) First method:
The minimum number of steps m is theequivalent time-step for a chlorine parcel traveling from itsinput location to the furthest output location. Mathematically,it can be expressed as m = (cid:24) T travel ∆ t (cid:25) = (cid:24)(cid:88) L ij v ij ∆ t (cid:25) , ij ∈ L path , (16)where ∆ t is determined by the L-W scheme stability condition;the total travel time is denoted by T travel ; L path representsthe links in the travel path, and the corresponding length andvelocity of each pipe ij in that path are L ij and v ij .To illustrate the application of this method, we use thenetwork comprised of two output sensors and one boosterstation controller; see Fig. 2. After a chlorine parcel with acertain mass and concentration is injected at the booster station,it travels to the sensors with various velocities in different links.Note that both the mass and concentration of the chlorine parcelreduce gradually due to consumption by the nodes and decay ofchlorine concentration. Furthermore, parcels (the rectangle inFig. 2) with different masses and concentrations (the areas andcolors of the rectangles) arrive at the sensors at various times.The number of parcels depends on the number of paths froman input location (booster) to output one (sensor). For example,there are three paths in Fig. 2 and the furthest location from theinput is Sensor 2, there are three travel times (three paths) forit from the booster, and the longest travel time should be usedwhen calculating (16).
2) Second method:
From a control theory perspective, thelower bound m can also be estimated by the settling time ofwater quality dynamic model (2). The settling time T s of adynamic system can be roughly computed by its dominant pole p which is equal to the eigenvalue of matrix A with the largestmagnitude. In particular, T s can be computed as T s = − t ln | p | , ∆ t is the sampling time of DT-LTI system (2); thedominant pole p located within a unit disk indicates that ln | p | is always negative. Consequently, the lower bound of m is m = (cid:24) T s ∆ t (cid:25) = (cid:24) − | p | (cid:25) . (17)Note that the equivalent travel time-step in (16) has a physicalmeaning for water quality dynamics which can be computedeasily. Meanwhile, the second method is estimated by the empir-ical formula (17). Section VI-B2 investigates the performanceof these two methods.Finally, the overall procedure of general, stability-preservingMOR procedure is summarized in Algorithm 1. In this algo-rithm, we use BPOD but roughly the same approach can beapplied to POD. Algorithm 1:
Stability-Preserving BPOD (SBPOD).
Input:
Full-order model ( A , B , C , D ) (2) and control input u ( k ) Output:
Stable reduced-order model ( A r , B r , C r , D r ) (3),outputs y ( k ) and ˆ y ( k ) if a posterior stabilization method then Standard BPOD approach given in Procedure 3 Solve a convex SDP problem (15) to stabilize thereduced-order model ( A r , B r , C r , D r ) end if if a priori stabilization method then Obtain parameter m via (16) or (17), and set length ofsnapshot as m = m Execute BPOD procedure to obtain reduced-ordermodels ( A r , B r , C r , D r ) (3) end if Generate y ( k ) and ˆ y ( k ) from the full- and reduced-ordermodels VI. C
ASE S TUDIES
We present three examples (three-node, Net1, and Net3networks [45]) to illustrate the performance of three differentMOR methods in terms of accuracy, computational load/time,and stability-preserving properties. Then, we test the two sta-bilization methods proposed in Section V. Lastly, we showcasethe application of MOR algorithms when MPC is used to controlwater quality dynamics for the full-order and reduced-ordersystems. All codes, parameters and tested networks simulatedvia EPANET Matlab Toolkit [46] are available on Github [47].The testing environment is a Win 10 Precision 7920 Tower withInter(R) Xeon(R) Gold 5218 CPU @2.3GHz and 64G memory.
A. Network settings, practical conditions, and validation1) Settings:
The basic information of the tested networksare in Tab. II. The three-node network in Fig. 3a is a self-designed network for illustrative purpose that includes onejunction J2, one pipe P23, one pump PM12, one tank TK3, andone reservoir R1. A booster station and a chlorine concentrationsensor are installed at J2 and TK3, which indicates it is a single-input single-output (SISO) system ( n u = n y = 1 ). The pipeP23 connecting J2 and TK3 is split into s L = 150 segmentsaccording to L-W scheme in Fig. 1a. Hence the dimension ofthe full-order system n x = 154 ( n N = n J + n R + n TK = 3 , Tab. IIB
ASIC INFORMATION OF THREE TESTED NETWORKS . Networks ∗ Full order n x Booster stationlocations ( n u ) Sensorslocation ( n y )three-nodenetwork { }
154 J2(1) TK3(1)
Net1 { } Net3 { } ∗ Number of each component in WDN: { n J , n R , n TK , n P , n M , n V } .Fig. 3. Three tested networks: (a) three-node network; (b) Net1 network; (c)Net3 network (zoomed-in part is the controlled area). and n L = n P × s L + n M = 151 ). The dimensions and otherparameters of the other networks are listed in Tab. II.
2) Non-zero initial conditions:
In a water quality sim-ulation, initial chlorine concentrations are typically non-zero(i.e., x (0) (cid:54) = ). With that in mind, some MOR algorithmsare designed with zero initial conditions. To tailor the MORalgorithms for water quality dynamics with non-zero initialconditions, we adopt the method from [48] as follows. Supposewe have a full-order model ( A , B , C , D ) with x o = x (0) (cid:54) = ,let ˜ x = x − x o , then ˜ x ( k + 1) = ˜ A ˜ x ( k ) + ˜ B ˜ u ( k ) (18a) ˜ y ( k ) = ˜ C ˜ x ( k ) + ˜ D ˜ u ( k ) , (18b)where ˜ A = A , ˜ B = (cid:2) B Ax o (cid:3) , ˜ C = C , ˜ D = (cid:2) D Cx o (cid:3) ,and ˜ u = (cid:20) u (cid:21) . That is, the non-zero initial condition is viewedas an input for the new augmented full-order model (18) with ˜ x o = such that the introduced MOR algorithms can be applieddirectly without any modification.
3) Validation:
To quantify the performance of MOR meth-ods such as accuracy, the comparisons of step responses for the8 a) (b)(c)
Fig. 4. Step responses (with amplitude of 50) of full-/reduced- order models for the three-node network under: (a) zero-initial conditions; (b) non-zero initialconditions. (c) The step response errors between full-/reduced- order models under zero-initial conditions as n r increases from 30 to 120. full-order model (2) and reduced-order model (3) are necessary.This is customary in MOR studies. The core idea of the designedexperiments is applying the same input signal u ( k ) to boththe full- and reduced-order systems and then comparing outputdifferences (i.e., step response errors between y ( k ) and ˆ y ( k ) ).To quantify the accuracy of the MOR algorithms, we use theroot-mean-square error (RMSE) defined as RMSE = (cid:118)(cid:117)(cid:117)(cid:116) m m (cid:88) k =1 || y ( k ) − ˆ y ( k ) || to quantify the error between full-order model output y ( k ) andreduced-order model ˆ y ( k ) of a time horizon of m time-steps. B. Accuracy and stability of MOR methods1) The three-node network:
The accuracy of MOR methodsof the three-node network under zero initial conditions (i.e., x o = ) is presented in Fig. 4a(left). The input is a step signalwith amplitude of 50, that is, we inject mg chlorine persampling time ( ∆ t = 20 seconds) at J2 in Fig. 3a. The blueline is the step response of the full-order model (2) ( n x = 154 ),and the dotted lines are the step responses of the reduced-ordermodels (3) obtained from BT, POD, and SBPOD ( n r = 30 , , and ). The data is generated from Procedures 1, 2, 3, andAlgorithm 1. Fig. 4a(right) depicts the magnified step responseswhen time is limited in [36000 , seconds. It can be seenthat (i) the reduced-order models from BT, POD, and SBPOD( n r = 30 , , and ) successfully maintain the input-outputrelationship with smaller orders ( n r < n x = 154 ), and (ii) atthe end of the test ( [36000 , sec), POD starts to show vasterror while the results of BT and BPOD are still close to thefull-order system.The accuracy test under non-zero initial conditions (i.e., x o (cid:54) = ) is also evaluated for the three-node network. All settingsincluding the input signal remain the same except that the initialchlorine concentrations at R1, J2, and TK3 are . mg/L, . mg/L, and . mg/L; the initial chlorine concentration in PM12and P23 are . mg/L, and . mg/L. The corresponding resultsare shown in Fig. 4b. Note that the output of POD is not shown due to it producing unreasonable results with huge RMSE =23 . . This is because the augmented system (18) views non-zero initial conditions as inputs. In this way, the input matrixof (18) ˜ B includes two parts: the actual chlorine injection B and virtual inputs from the initial condition Ax o . If the majority ofelements in vector Ax o are non-zeros, then this system has toomany inputs making the POD hard to capture the input-outputrelationship. With that in mind, the non-zero initial conditionsare not problematic for BT and BPOD since matrices W X or H m have the term CB that cancels such effect. It indicatesPOD is unsuitable for water quality dynamics with non-zeroinitial conditions.To further test the accuracy of all MOR methods (i.e.,classical BT, POD, and BPOD and proposed SPOD) and theimpact of parameter n r , we show the norm of step responseerrors (denoted by (cid:107) Error (cid:107) = (cid:107) y − ˆ y (cid:107) ) between full-order model (2) and reduced-order models (3) under zero-initialconditions in Fig. 4c. As n r increases from to , the stepresponse errors of all MOR methods decrease. We note thefollowing. First, the amplitude of step response is in level(see Fig. 4c) while the error is in − or lower level for BTand SBPOD. This indicates that the error is tiny—it can beneglected and hence the reduced models are accurate. Second,POD does not perform well even n r reaches to ( n x = 154 ).This is due to POD only considering the controllability Gramianin Section IV. Third, classical BPOD (with m = m = 160 )performs worse than SBPOD ( m = 400 > m = 160 );this indicates classical BPOD can not ensure the stability oraccuracy and is not suitable for water quality dynamics. Thevalue of parameter m = 160 for this simple network can beobtained easily by Equation (16), that is, it equals the time-stepfrom J2 to TK3. We show details of how to find m or m usinga more complex Net1 network in the next section. All these testresults corroborate the analyses performed in Section IV andTab. I.Furthermore, the RMSEs of all MOR methods (classical BTand POD, proposed SBPOD) for the three-node network underzero or non-zero initial conditions are shown in Tab. III. The9 ab. IIIA
CCURACY WITH ( I ) STEP SIGNAL AS INPUTS , AND ( II ) NON - ZERO ( ZERO ) INITIAL CONDITIONS AND CORRESPONDING COMPUTATIONAL TIME ( INSECOND ) AMONG THREE CLASSICAL
MOR
METHODS FOR THREE TESTED NETWORKS . RMSEs when applying step signalswith x o = ( x o (cid:54) = ) Computational time(in second)three-node network Net1 Net3 three-node network Net1 Net3BT . × − ( . × − ) . × − ( . × − ) NA ∗ (NA ∗ ) 0.245 17.2 Intractable ∗ POD . × − (23.51) . × − ( . × ) . × − (0.49) 0.219 6.4 245.8 SBPOD . × − ( . × − ) . × − ( . × − ) . × − ( . × − ) 0.241 6.6 150.4 ∗ NA is “Not Applicable” as BT fails to produce the reduced-order model after 12 hours of simulation. classical BPOD cannot ensure stability under either zero or non-zero initial conditions (see the further test in Section VI-B2for more details and discussions). Hence, the corresponding
RMSEs are not included in the table. Based on these small
RMSEs , we conclude that MOR methods perform well underzero initial conditions while BT and SBPOD perform the bestunder the realistic non-zero initial conditions of water qualitydynamics.
2) Net1 network:
The accuracy tests of classical BT,classical POD, and SBPOD methods for the Net1 networkunder zero/non-zero initial conditions are performed. Insteadof presenting all overlapped lines in figures, we adopt
RMSEs in Tab. III to compare accuracy among all MOR methods. Thefull order of Net1 is n x = 1293 , and BT and SBPOD reducethe order to n r = 113 and n r = 117 under non-zero initialconditions. Note that POD fails to provide accurate reduced-order models from the corresponding RMSEs shown in Tab. III.For zero initial conditions, MOR methods (i.e., classical BT andPOD and SBPOD) successfully reduce the system orders, and n r of BT, POD, and SBPOD are , , and . From thesesmall RMSEs in Tab. III, it is evident that the reduced-ordersystems are accurate.Next, we show how parameter m is obtained for Net1, andtest the performance of BPOD (with an m less than m ), SBPODwith the posterior stabilization method (same m with BPOD),and SBPOD with the priori stabilization method ( m larger thanthe m ). Computing lower bounds for m . The first step in generatingstability-preserving reduced-order models is to compute thelower bound for m for the SBPOD method. When the lengthof snapshots m of BPOD is large enough for water quality dy-namics, BPOD can preserve stability; otherwise, it cannot. Thisis the stability-preserving condition for BPOD as introduced inSection V-B. The two approaches to finding m are tested forthe Net1 network.The first method given in (16) is used to obtain the proper m which requires computing the travel time of the last chlorineparcel injected at the booster J10 to the furthest sensor J23.The first (second) sensor installed at J22 (J23) receives two(three) step responses from the booster installed at J10, andthe corresponding time to reach the peak is marked in the thirdcolumn of Fig. 5a. The peak time of the last step response at J22(J23) is t = 16550 ( t = 45360 ), and the arrival or travel time of the last chlorine parcel (step response) is slightly less thanthe peak time. That is, 16000 for J22 and 39750 for J23 (we donot mark it in Fig. 5). Hence, parameter m can be obtained bythe longest time 39750 divided by the sampling time ∆ t = 15 seconds. That is, parameter m is set to via the first method.We also choose Equation (17), that is the second method, toobtain the proper m which requires computing the dominantpole p . For Net1, p = 0 . ± . i , and T s = − /ln | p | =2350 . resulting in m = 2351 for the second method.Although parameter m obtained by two methods are differ-ent, both are reasonable. After testing, we notice that when m ∈ [2351 , , the reduced-order model for Net1 via BPODshows instability as an oscillation tail is observed at the end ofthe step response; see the first column of figures in Fig. 5a. Computing and simulating reduced-order models.
Afterfinding the appropriate lower bounds for m , Algorithm 1 isimplemented to compute a stable reduced-order model forNet1. The step responses of the posterior stabilization method(stabilizing unstable reduced-order model from classical BPODvia solving an SDP (15)) are shown in the second columnof Fig. 5a. Comparing with the unstable step responses fromclassical BPOD with m = 2650 (i.e., the first column ofFig. 5a), the yielding responses after posterior stabilizationare not oscillatory but are unfortunately inaccurate. The stepresponses of the priori method ( m is set to that is largerthan found by stability-preserving condition) are shownin the third column of Fig. 5a. Compared with results of theposterior stabilization method, the responses are stable and themethod results in accurate estimation of the step response (seethe zoom-in areas).
3) Net3 network:
We present the results for a larger network,namely Net3 shown in Fig. 3c, under zero ( x o = 0 mg/L)or non-zero initial chlorine concentrations ( x o = 0 . mg/L).Two booster locations (J237 and J247) in Fig. 3c inject chlorinewith a rate of mg and mg per sampling time ( ∆ t =0 . second), that are two step signals with amplitudes of 50and 40. Fig. 5b presents step responses of full-order model andreduced-model from SBPOD with non-zero initial conditionsat three sensors (J255, J241, and J249). The full order of Net3is n x = 29347 , and the reduced-order n r from SBPOD is .Note that the intractable BT method fails to produce a reduced-order model after 12 hours of simulation; thus, the results of BTare not included in Fig. 5b. Furthermore, POD method is not10 a)(b) Fig. 5. (a) step responses from BPOD and SBPOD with two stabilization methods for Net1 under non-zero initial conditions; (b) step responses fromdifferent output locations of reduced-order models via BPOD methods for Net3 under on-zero initial chlorine concentrations. suitable for non-zero initial condition case, and its results arenot shown in Fig. 5b either. For zero initial condition case, bothPOD and SBPOD reduce system order successfully, and thereduced-orders n r are and . Instead of presenting stepresponse figures, we showcase the RMSEs of MOR methodsin Tab. III.
C. Computational time
Tab. III presents the computational time for computing thereduced-order models for the three algorithms for different net-works. We note that the computational time of classical BPODand SBPOD are nearly identical (seeing that the computationof m requires little to no time), and only results of SBPODare shown for simplicity. The full-order n x and the length ofsnapshots m for the three-node network are set as and when testing computational time. For Net1, they are and ; for Net3, these values are and .From Tab. III, it can be seen that BT is intractable when n x is large since it fails to give a reduced-order model for Net3after running the algorithm for 12 hours. The computationaltime for SBPOD is less than the one by POD for Net3 whichinterestingly contradicts the theoretical complexity from Tab. I.This happens because we implement a technique to avoidmatrix-matrix multiplication by exploiting the structure of theblock Hankel matrix in (11). Consequently, this results in areduction of the computational time of the SBPOD method. TheGithub codes [47] also include a detailed description of how theproblem structure is exploited to reduce the computational time.The details are omitted in this paper for brevity. D. SBPOD-Based MPC vs. Full-Order MPC
One of the main objectives of investigating MOR algorithmsfor water quality dynamics, besides the ability to performforward state-simulations more efficiently, is to showcase thealgorithms’ performance in the context of real-time feedback control (i.e., controlling chlorine concentrations x ( k ) throughthe installed booster station controllers u ( k ) ) or state observa-tion and estimation (i.e., estimating the states x ( k ) via noisysensor measurements y ( k ) ). To that end, the objective of thissection is to test the applicability of the presented SBPOD-generated, reduced-order model within a model predictive con-trol (MPC) framework that regulates chlorine concentrationsthrough controlling dosages of injected chlorine at the boosterstations. We specifically focus on SBPOD herein as it producesthe best performance (in terms of accuracy under differentinitial conditions, computational time, and stability-preserving)among the other tested methods.The motivation for using the reduced-order model is asfollows: the regulation of chlorine concentrations with strictstate constraints (e.g., as the ones mandated by the US EPA)in a water network can be posed as a constrained MPCoptimization—as we have shown in our recent work [1]. Wehave also showcased that the control routine for the full-ordersystem becomes computationally intractable when consideringMPC with state- and input-constraints and a large controlhorizon (e.g., 30 minutes) even for a mid-size network (e.g,Net1). To solve this issue, we have relaxed the strict constraintsand converted the intractable constrained MPC into an uncon-strained and tractable one. However, this leads to another issue.That is, the unconstrained MPC cannot guarantee the state andinput bounds, which is one of the limitations of our previouswork [1]. In this paper, we investigate the following researchquestion in this section: Can the water system operator utilize the SBPOD-generated, reduced-order model with constrainedMPC to determine the desired dosages of injectedchlorine? How would such an application compareto using a similar control routine for the full-ordermodel? Would utilizing the reduced-order modelresult in a significant increase in the operation costof controlling chlorine concentrations? How does b)(a) Fig. 6. The comparisons of (a) control actions and (b) subsequent chlorineconcentrations solved by the full- and reduced-order model (i.e., SBPOD) forthe constrained MPC problems. that impact the controlled chlorine consecrations?
To answer these questions, we use the Net1 network as anillustration. Similar performance is exhibited for other net-works. Note that we set a small control horizon parameter,that is 5 minutes, for Net1 to make optimization problemswith full-order models tractable so that the control performanceof full-/reduced- order models based constrained MPC can betested and compared. First, the full-order models are computedfor a time horizon of 24 hours. This is then followed bycomputing the SBPOD-generated reduced-order models via thepriori stabilization method. Second, the optimal control signalsof chlorine dosages for both the full- and reduced-order modelsare computed via the constrained MPC from study [1]. Thereference signal for the MPC algorithm to track is set as 1.5mg/L. That is, the goal of the controller is to maintain thechlorine concentration at around 1.5 mg/L at J22 and J23 (thesensor nodes in Net1) by controlling the injected mass rate ofchlorine u ( k ) at J10 (the only control node in Net1). Finally, thecontrol signals are then applied to the water quality simulation inthe Net1 network. For brevity, we do not reproduce the detailshere; the interested reader can examine the provided Githubcodes for the parameters and settings [47].To assess the performance of the SBPOD-based chlorinedosage constrained model predictive controller, in comparisonwith the similar MPC that utilizes the more computationallyexpensive full-order one, we test for three metrics: the compu-tational time of solving optimization problems, the difference inthe cost function of running the controllers (i.e., the operationalcost of running the controller), and the evolution of the chlorineconcentration. The three metrics are assessed for the twoapproaches: (i) the SBPOD-based constrained MPC and (ii) thestandard constrained MPC for the full-order model. Fig. 6a andFig. 6b show the injected mass rate u ( t ) and the correspondingchlorine concentrations at J22 and J23 for (i) and (ii) for theentire time horizon for both the reduced- and full-order models.The RMSEs for control input u and output y between thefull-order model and reduced-order one are . × − and . × − . The computational time of solving constrainedMPC problems with full- and reduced-order models are . seconds and . seconds. Both operational costs from thesetwo different controller are . × .By examining the results, we observe that the reduced-ordermodel produces nearly identical performance (in terms of chlo-rine concentration dynamics and the optimal chlorine dosages) when compared with the full-order one, while incurring smallercomputational time. We note here that solving constrainedMPC with full-order models for Net3, a larger network thanNet1, is computationally intractable. After a running time of12 hours, solutions are not obtained since this optimizationproblem has at least n x = 29374 optimization variables andmillions of constraints. However, with the reduced-order model,the constrained optimization problem only has n r = 453 optimization variables and at around ten thousand constraintsmaking it computationally tractable. In short, the number ofoptimization variables and constraints of constrained MPC withreduced-order models is small regardless of the size of thenetwork/full-order model since parameter n r is usually withina thousand after MOR. The reduced-order model enables watersystem operators to use constrained MPC that guarantees theinput and output bounds while remaining tractable for largernetworks.VII. W HAT S HOULD THE O PERATOR U SE ? C ONCLUSIONSAND F UTURE W ORK
The presented research in this paper explores the potentialof model order reduction for water quality dynamics in drink-ing water networks. The presented stability-preserving BPOD(SBPOD) algorithm is more usable than its counterparts dueto its accuracy, low computational burden, tolerance to initialconditions, and stability-insuring property. Compared with theclassical BT method, SBPOD is computationally tractable whileyielding accurate estimates of output measurements of chlorineconcentrations. Furthermore, compared to POD or BPOD,SBPOD handles a variety of initial conditions while ensuringstability and accuracy of estimated outputs from the reduced-order models. Finally, we demonstrate that when using theproposed SBPOD within an MPC framework to control chlorineconcentrations, the water system operator can use the SBPOD-generated reduced-order model instead of the full-order onewithout incurring any significant losses in the operational costof controlling the network, while saving orders of magnitudein computational time. This makes control algorithms such asMPC amenable to real-time implementation—as a result ofutilizing the reduced-order models. In short, and through thecomprehensive testing and theoretical discussions, the watersystem operator can utilize the SBPOD to perform accuratesimulations and constrained, feedback control for water qualitydynamics.Future research directions includes investigating model orderreduction algorithms for water quality dynamics that involvemulti-species interactions, in comparison with models of fo-cusing on single-species ones. A state-space model for waterquality dynamics in a multi-species framework would likelyyield a nonlinear state-space description that then necessitatedesigning MOR algorithms for nonlinear dynamics. To thatend, future work will focus on (i) thorough modeling of multi-species reaction dynamics and (ii) reducing such dynamics tolower order models.12
EFERENCES[1] S. Wang, A. F. Taha, and A. A. Abokifa, “How effective is modelpredictive control in real-time water quality regulation? state-spacemodeling and scalable control,”
Water Resources Research , 2020, inpress.[2] P. Boulos, K. Lansey, and B. Karney,
Comprehensive Water DistributionSystems Analysis Handbook for Engineers and Planners . MWH Soft,Incorporated, 2006.[3] C. W. Rowley, “Model reduction for fluids, using balanced properorthogonal decomposition,”
International Journal of Bifurcation andChaos , vol. 15, no. 03, pp. 997–1013, 2005.[4] U. Baur, P. Benner, and L. Feng, “Model Order Reduction for Linearand Nonlinear Systems: A System-Theoretic Perspective,”
Archives ofComputational Methods in Engineering , vol. 21, no. 4, pp. 331–358,2014.[5] K. Willcox and J. Peraire, “Balanced model reduction via the properorthogonal decomposition,”
AIAA Journal , vol. 40, no. 11, pp. 2323–2330, 2002.[6] L. Montier, T. Henneron, B. Goursaud, and S. Clenet, “Balanced ProperOrthogonal Decomposition Applied to Magnetoquasi-Static ProblemsThrough a Stabilization Methodology,”
IEEE Transactions on Magnet-ics , vol. 53, no. 7, 2017.[7] V. M. Adamjan, D. Z. Arov, and M. Kre˘ın, “Analytic properties ofSchmidt pairs for a Hankel operator and the generalized Schur-Takagiproblem,”
Mathematics of the USSR-Sbornik , vol. 15, no. 1, p. 31, 1971.[8] K. Glover, “All optimal hankel-norm approximations of linear multi-variable systems and their L ∞ error bounds,” International journal ofcontrol , vol. 39, no. 6, pp. 1115–1193, 1984.[9] B. Moore, “Principal component analysis in linear systems: Controllabil-ity, observability, and model reduction,”
IEEE transactions on automaticcontrol , vol. 26, no. 1, pp. 17–32, 1981.[10] L. Sirovich, “Turbulence and the dynamics of coherent structures. i.coherent structures,”
Quarterly of applied mathematics , vol. 45, no. 3,pp. 561–571, 1987.[11] E. Grimme, “Krylov projection methods for model reduction,” Ph.D.dissertation, 1997.[12] A. Antoulas, J. Ball, J. Kang, and J. Willems, “On the solution ofthe minimal rational interpolation problem,”
Linear Algebra and itsApplications , vol. 137, pp. 511–573, 1990.[13] C. A. Beattie and S. Gugercin, “Interpolation theory for structure-preserving model reduction,” in . IEEE, 2008, pp. 4204–4208.[14] Z. Bai, “Krylov subspace techniques for reduced-order modelingof large-scale dynamical systems,”
Applied Numerical Mathematics ,vol. 43, no. 1-2, pp. 9–44, 2002.[15] K. Gallivan, A. Vandendorpe, and P. Van Dooren, “Model reduction andthe solution of sylvester equations,”
MTNS, Kyoto , vol. 50, 2006.[16] S. Gugercin, “An iterative svd-krylov based method for model reductionof large-scale dynamical systems,”
Linear Algebra and its Applications ,vol. 428, no. 8-9, pp. 1964–1986, 2008.[17] S. Lall and C. Beck, “Error-bounds for balanced model-reduction oflinear time-varying systems,”
IEEE Transactions on Automatic Control ,vol. 48, no. 6, pp. 946–956, 2003.[18] J. L. Lumley,
Stochastic tools in turbulence . Courier Corporation, 2007.[19] P. Astrid, S. Weiland, K. Willcox, and T. Backx, “Missing pointestimation in models described by proper orthogonal decomposition,”
IEEE Transactions on Automatic Control , vol. 53, no. 10, pp. 2237–2251, 2008.[20] B. Ulanicki, A. Zehnpfund, and F. Martinez, “Simplification of waterdistribution network models,” in
Proc., 2nd Int. Conf. on Hydroinfor-matics . Balkema Rotterdam, Netherlands, 1996, pp. 493–500.[21] U. Shamir and E. Salomons, “Optimal real-time operation of urban waterdistribution systems using reduced models,”
Journal of Water ResourcesPlanning and Management , vol. 134, no. 2, pp. 181–185, 2008.[22] A. Preis, A. J. Whittle, A. Ostfeld, and L. Perelman, “Efficient hy-draulic state estimation technique using reduced models of urban waternetworks,”
Journal of Water Resources Planning and Management , vol.137, no. 4, pp. 343–351, 2011.[23] L. Perelman and A. Ostfeld, “Water distribution system aggregationfor water quality analysis,”
Journal of water resources planning andmanagement , vol. 134, no. 3, pp. 303–309, 2008.[24] P. D. Lax and B. Wendroff, “Difference schemes for hyperbolic equa-tions with high order of accuracy,”
Communications on pure and appliedmathematics , vol. 17, no. 3, pp. 381–398, 1964. [25] C. Chen,
Linear System Theory and Design , ser. The OxfordSeries in Electrical and Computer Engineering. Oxford UniversityPress, Incorporated, 2013. [Online]. Available: https://books.google.com/books?id=XyPAoAEACAAJ[26] C. Himpe, “emgr-The empirical Gramian framework,”
Algorithms ,vol. 11, no. 7, pp. 1–27, 2018.[27] ——, “Comparing (Empirical-Gramian-Based) Model Order ReductionAlgorithms,” 2020. [Online]. Available: http://arxiv.org/abs/2002.12226[28] S. Lall, J. E. Marsden, and S. Glavaˇski, “Empirical model reductionof controlled nonlinear systems,”
IFAC Proceedings Volumes , vol. 32,no. 2, pp. 2598–2603, 1999.[29] R. A. Horn and C. R. Johnson,
Matrix analysis . Cambridge universitypress, 2012.[30] D. Amsallem and C. Farhat, “Stabilization of projection-based reduced-order models,”
International Journal for Numerical Methods in Engi-neering , vol. 91, no. 4, pp. 358–377, 2012.[31] J.-R. Li and J. White, “Low rank solution of lyapunov equations,”
SIAMJournal on Matrix Analysis and Applications , vol. 24, no. 1, pp. 260–280, 2002.[32] V. Y. Pan and Z. Q. Chen, “The complexity of the matrix eigenproblem,”in
Proceedings of the thirty-first annual ACM symposium on Theory ofcomputing , 1999, pp. 507–516.[33] X. Feng, W. Yu, and Y. Li, “Faster matrix completion using random-ized svd,” in . IEEE, 2018, pp. 608–615.[34] S. Prajna, “Pod model reduction with stability guarantee,” in , vol. 5. IEEE, 2003, pp. 5254–5258.[35] C. W. Rowley, T. Colonius, and R. M. Murray, “Model reduction forcompressible flows using POD and Galerkin projection,”
Physica D:Nonlinear Phenomena , 2004.[36] M. F. Barone, I. Kalashnikova, D. J. Segalman, and H. K. Thornquist,“Stable Galerkin reduced order models for linearized compressibleflow,”
Journal of Computational Physics , 2009.[37] I. Kalashnikova and M. Barone, “On the stability and convergence of agalerkin reduced order model (rom) of compressible flow with solid walland far-field boundary treatment,”
International journal for numericalmethods in engineering , vol. 83, no. 10, pp. 1345–1375, 2010.[38] I. Kalashnikova and S. Arunajatesan, “A stable galerkin reduced ordermodel for compressible flow,” in , 2012.[39] S. Sirisup and G. E. Karniadakis, “A spectral viscosity method for cor-recting the long-term behavior of pod models,”
Journal of ComputationalPhysics , vol. 194, no. 1, pp. 92–116, 2004.[40] Z. Wang, I. Akhtar, J. Borggaard, and T. Iliescu, “Proper orthogonaldecomposition closure models for turbulent flows: a numerical compar-ison,”
Computer Methods in Applied Mechanics and Engineering , vol.237, pp. 10–26, 2012.[41] I. Kalashnikova, B. van Bloemen Waanders, S. Arunajatesan, andM. Barone, “Stabilization of projection-based reduced order models forlinear time-invariant systems via optimization-based eigenvalue reassign-ment,”
Computer Methods in Applied Mechanics and Engineering , vol.272, pp. 251–270, 2014.[42] B. N. Bond and L. Daniel, “Guaranteed stable projection-based modelreduction for indefinite and unstable linear systems,” in
IEEE/ACM In-ternational Conference on Computer-Aided Design, Digest of TechnicalPapers, ICCAD , 2008.[43] M. Benosman, J. Borggaard, and B. Kramer, “Robust pod modelstabilization for the 3d boussinesq equations based on lyapunov theoryand extremum seeking,” in .IEEE, 2017, pp. 1827–1832.[44] A. M. Lyapunov, “The general problem of the stability of motion,”
International journal of control , vol. 55, no. 3, pp. 531–534, 1992.[45] L. A. Rossman et al. , “EPANET 2: users manual,” 2000.[46] D. G. Eliades, M. Kyriakou, S. Vrachimis, and M. M. Polycarpou,“Epanet-matlab toolkit: An open-source software for interfacing epanetwith matlab,” in
Proc. 14th International Conference on Computing andControl for the Water Industry (CCWI) , The Netherlands, Nov 2016, p. 8.[47] ShenWang9202, “ShenWang9202/MOR4WQM: MOR FOR Waterquality dynamics,” Feb. 2021. [Online]. Available: https://doi.org/10.5281/zenodo.4554174[48] M. Heinkenschloss, T. Reis, and A. C. Antoulas, “Balanced truncationmodel reduction for systems with inhomogeneous initial conditions,”
Automatica , vol. 47, no. 3, pp. 559–564, 2011. [Online]. Available:http://dx.doi.org/10.1016/j.automatica.2010.12.002, vol. 47, no. 3, pp. 559–564, 2011. [Online]. Available:http://dx.doi.org/10.1016/j.automatica.2010.12.002