- Email: [email protected]
Robust parameter estimation and output prediction for reactive carrier-load nonlinear dynamic networks
13th IFAC Symposium on Large Scale Complex Systems: Theory and Applications July 7-10, 2013. Shanghai, China
Robust parameter estimation and output prediction for reactive carrier-load nonlinear dynamic networks Mietak A. Brdys ∗,∗∗ Tomasz Zubowicz ∗ Krzysztof Arminski ∗ ∗
Department of Control Systems Engineering Gdask University of Technology, ul. G. Narutowicza 11/12, 80-233 Gdask, Poland (e-mail: [email protected]; [email protected]; [email protected]) ∗∗ Department of Electronic, Electrical and Computer Engineering, College of Engineering and Physical Sciences University of Birmingham, Edgbaston, Birmingham B15 2TT, UK (e-mail: [email protected]) Abstract: In this paper an extension of on-line model simplification technique for a class of networked systems, namely reactive carrier-load nonlinear dynamic networked system (RCLNDNS), kept within point-parametric model (PPM) framework is addressed. The PPM is utilised to acquire a piece wise constant time-varying parameter linear structure for the RCLNDNS suitable for the on-line one step ahead prediction that may be applied to monitoring and model predictive control purposes. The advantageous structure of PPM offers a considerable simplification of those algorithms in terms of computational burden in comparison to the nonlinear structures. Moreover, the availability of technical tools makes the analysis of the final form algorithms straightforward. An application to water quality estimation/prediction in drinking water distribution system illustrates the technology. Keywords: networked systems, robust parameter estimation, prediction 1. INTRODUCTION Development of a modern world has contributed to enlarging the number of large scale complex systems (LSCS) which deliver fundamental goods driving the existence and sustainable development of modern societies. These are called the critical infrastructure systems (CIS). Most of the CIS are characterised by complex, highly interconnected, networked structure e.g. power grid, drinking water distribution system, telecommunication networks etc. The first two constitute a very specific type of networked systems which can be described as reactive carrier-load networked system where one can distinguish a medium carrying a load of certain quality which can be altered by some mechanisms resulting from occurring processes. These processes result from given system operation and in majority of cases are guided by a nonlinear dynamics. Henceforth, the system at hand, to be tackled with, is a reactive carrier-load nonlinear dynamic networked system (RCLNDNS). Monitoring, control and security of RCLNDNS-CIS constitute timely as well as highly important research topic having strong impact on daily life of the societies worldwide. Moreover, it requires appropriate structures and algorithms in order to achieve reliable (in terms of robustness) solutions. In this work one of the most fundamental CIS is considered as RCLNDNS application example, namely the drinking water distribution system (DWDS). 978-3-902823-39-7/2013 © IFAC
426
The DWDS has been selected since not only that it is well fitted example of RCLNDNS but also it is characterised by multiple time-scale dynamics not only in the context of transportation-reaction part but also in terms of biochemical processes constituting the reactive part of the system. Furthermore, available measurements in terms of their number and diversity are very limited resulting in poor access to process information. A robust monitoring of DWDS has been addressed in (Langowski and Brdys (2007)) by applying an interval observer technology based on cooperative systems theory. Control of RCLNDNS have for the long time now been addressed in the model predictive algorithm framework (MPC). The MPC family constitutes of several mainstream approaches: richly referenced e.g. within (Maciejowski (2002)). Although, nonlinear MPC algorithms have been studied extensively in recent years it still comes as a fact, that nonlinear optimisation task that needs to be solved to find an optimal/suboptimal solution are very computationally demanding tasks, which is why having a neat linear model is still attractive in many applications. Handling such a model in manner to account for changing operational conditions, uncertainty in both structure and parameters most certainly requires proper mechanisms. In previous works (Brdys and Chang (2002); Chang (2003); Chang et al. (2004)) authors introduced a pointparametric model (PPM) approach to deliver such a model 10.3182/20130708-3-CN-2036.00026
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
of linear structure and time-varying parameters as well as algorithms for its parameter robust estimation and in consequence robust representation of DWDS process dynamics. The structure of the PPM was derived based on the modified path-tracking algorithm (Zierolf et al. (1998); Shang et al. (2000, 2002)) in order to account for the storage facilities within the network. In this approach a bound model for the uncertainty was utilised. The PPM framework was built and verified using linear quality model in the reactive part of the system. This paper extends the concept of PPM framework from the linear networked systems to RCLNDNS. To achieve that, a more sophisticated implicit simulation model is introduced within the PPM estimator/predictor structure in order to reduce the structural error. Moreover, the verification process utilises also the enhanced, nonlinear, model for plant simulation briefly introduced in Section 3. The paper is organised as follows. The preliminaries are given in forthcoming Section 2. Section 3 introduces and characterises the structures and algorithms utilised for robust estimation and prediction of PPM parameters and outputs. In Section 4 a numerical example based on unique DWDS model to illustrate the technology has been shown. Finally, Section 5 concludes the paper. 2. PRELIMINARIES Consider a nonlinear dynamic advective-reactive transport processes in axial direction within the network environment: { vjp ∂xjp Ψijp + Ξip , bulk flow i (1) ∂t Ψj p = Ξip , link wall ( ) ∆ where: t is time; Ψjp i = Ψjp i t, xjp is the i-th concen( ) ∆ tration in the jp -th pipe; Ξip = Ξip Ψjp , vjp , sjp are the rates dependent on all process variables; sjp is the bulk flow contact surface with the pipe wall; D is the dispersion coefficient; i and jp denotes the state variables and pipes of the network respectively. The processes may be described in general as: ( ) ( ) Ξp Ψjp , vjp , sjp = Arp Ψjp , vjp , sjp (2) ( ) where: Ξip Ψjp , vjp , sjp are the reaction rates and Ξp = [ 1 2 ]T Ξp , Ξp , ..., Ξnp where n is the total number of processes and ( state variables ) as well; A is the coefficient matrix; rip (Ψ)jp , vjp , sjp is the process rate vector and rp = [ 1 2 ]T rp , rp , ..., rpz where z denotes the processes involved in the model. Note that the interaction between the bulk and wall species is integrated into the Ξip . The storage facilities may be described as:( ) d Ψil Vl = Fin Ψiinl Fout Ψil + Vl Ξil (3) dt where: Fin , Fout are the inflow and outflows respectively; Ψiinl is the concentration of the i th specie in the inflow; ∆
Vl = Vl (t) is the volume of the medium contained in ∆ the l th tank; Ξil = Ξil (Ψl ) are the quality processes dependent on all species concentration and is assumed that may be different than Ξip ; l is the tank number. 427
The definition and mechanisms related to Ξit are given analogously. It is assumed that the robust information about an arbitrary, measurement accessible, variable of the selected process may be obtained at a given point in space utilising the simplified parameter time-varying model structure: y (k) = θ (k) φ (k 1) + ε (k) (4) where: y (k) is the output; θ (k) and ε (k) represent model parameter vector and structure error and disturbance impact respectively and furthermore constitute the uncertainty within the model; φ (k 1) is the regressor vector; k 2 [k0 , k0 + Tm ]Z denotes time where k0 is current time instant and Tm is a modelling horizon length; Z yields the integer field. The internal structure of (4) results from the discrete MISO model structure: yjo (k) = +
nT ∑
i=1 nI ∑
bjo ,i (k) yjo (k ∑
i)
ajo ,jI ,i (k) ujI (k
i) + εjo (k)
(5)
jI =1 i∈I
∧jo 2 1, nM where: yjo is the jo -th monitored node output and jo 2 1, nM and nM is the total number of monitored nodes; ujI is the jI -th plant input and jI 2 1, nI and nI is the total number of inputs contributing to the output yjo ; bjo ,i and ∆
ajo ,jI ,i are time-varying model parameters; I = [il , ..., iu ] and il and iu define range of delays in the control input; and the constructions are as follows: θjo (k) = [θjo ,1 , ..., θjo ,n ] = ] [ bjo ,1 (k), ..., bjo ,nT (k), (6) = ajo ,1,il (k), ..., ajo ,1,iu (k) ..., ajo ,nI ,il (k), ..., ajo ,nI ,iu (k) ∧n = nT + nI (k ijo ,nI ,u ) ] [ yjo (k 1) , ..., yjo (k nT ) , φjo (k 1) = u1 (k ijo ,1,l ) , ..., u1 (k ijo ,1,u ) , ..., unI (k ijo ,nI ,l ) , ..., unI (k ijo ,nI ,u ) (7) Moreover, if the input u[k0 ,k] () is valued on a compact set over the time horizon [k0 , k] then there exist an upper and lower bound on the trajectories of θ (k) and ϵ (k); θu,0 (k), θl,0 (k) and ϵu,0 (k), ϵl,0 (k) respectively, which guarantee that the model output y[k0 ,k] () is equal to p p the plant output y[k (), y[k0 ,k] () = y[k () and in 0 ,k] 0 ,k] particular k = k0 + Tm . Furthermore, it is assumed that the information about the bounds on the input uncertainty umin u(k) umax are known a priori. It is also assumed that a good approximation of the tightest bounds on parameter values can be made at time t0 as in fact they are not known. Although θ (k) can change abruptly its value at the time instant kj 2 [k0 , k0 + Tm ], j 2 1, Λ, the changes may be considered slow over a shorter time horizon [kj−1 , kj ]. The time instances kj , to account for the uncertainty in the plant transportation process, are to be found in an online manner. Therefore, the modelling time horizon Tm is partitioned into Λ time-slots (9) which in fact change their numerical amount with time thus Λ should be in fact
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
replaced with Λ(k) resulting in the following definition of the j-th time-slot: { } ∆ Sj = k 2 Z : kj−1 k kj , j 2 1, Λ(k) (8) Moreover, not all parameters contained within θ (k) have non-zero values. Hence, an active parameter sub-vector { } may be defined θj = θij ij ∈Ij where Ij is a set of active parameter indexes over Sj . Furthermore, θ (k) = θj for is assumed to be constant. Therefore, a significant reduction in the bound determination problem complexity may be gained. In consequence, the orthotope being searched for is given as: ( j ) θ , ϵ (k) 2 Θ, nc nc (θ , ϵ ) 2 Rn : n n c c =[ θ φ](k 1) + ϵ (k) , ( l u l u ) ∆ yn(k) Θ θ , θ , ϵ , ϵ = θ c 2 [ θl , θu] , nc l u ϵ 2 ϵ , ϵ , ne = 1, Ne , k 2 [k0 , k0 + Tm ] (9) The final time-space point parametric piecewise model structure is as follows: y (k) = θ (k) φ (k 1) + ϵ (k) θ( (k) = θ)j ∧ k 2 Sj ∧ j 2 1, Λ θj , ϵ (k) 2 Θ
(10)
The proof verifying the existence of such a model can be found in (Chang (2003)). In order to ensure the PPM to be suitable for robust MPC purposes the uncertainty radius must be kept sufficiently small. On the other hand partitioning the modelling horizon into slots will most certainly increase the width of the estimated parameter envelopes; however the piece-wise continuity may be attractive to handle by MPC internal algorithms. Finding a suitable compromise between those two is a problem on its own and will not be addressed within this work. The robust parameter and output estimation is achieved throughout solving dedicated optimisation tasks defined in the forthcoming Section 3. 3. DWDS CASE STUDY 3.1 Implementation of PPM framework nc Let y[k () denote the model response to an input 0 ,k] nc u[k0 ,k] () over a time horizon [k0 , k + Tm ]. Then there nc exists a trajectory θnc () that results in y[k () = 0 ,k] p nc y[k0 ,k] (). Different u[k0 ,k] () scenarios will require different parameter scenarios, therefore: y nc (k) = θnc (k) φ (k 1) + ϵnc (k) (11) where: ne is the input scenario number and ne = 1, ..., Ne where Ne is the total number of scenarios.
The total number of input scenarios required to find the parameter bounding envelopes, in the direct approach, should be equal to a total number of possible input combinations that is Ne = 2N i where Ni is the total number of input time horizon samples corresponding to output observations on [k0 , k0 + Tm ]. However, in (Chang (2003)) 428
Fig. 1. PPM information structure it has been shown that a much smaller population of welldesigned experiments is sufficient to obtain satisfactory results. Thus Ne will be limited within this paper to a number of required experiments as proposed in (Chang (2003)). Structure of information exchange In order to utilise the modified PPM framework a proper information structure between the subsystems is required – see Fig. 1. The initial information about process parameters and measurements is utilised by the Calibrator subsystem in order to feed the data into the implicit model Simulator and build up the Plant information. The Plant information enriches this information with the disturbance predictions and plant structural information. Appropriate data packages are then redistributed between the internal Simulator, Explicit Model Structure Determination unit and Estimator. The Estimator sub-module utilises the explicit model simulator in order to acquire the data for robust parameter prediction based on previously defined structure (10) build up within Model Structure Determination unit. Based on the calculated model parameter envelopes the Explicit model sub-module produces the robust output prediction and a nominal one assumed to be a Chebyshev centre. The information about the parameter and output robust predictions (as well as the nominal one) is fed to the PPM estimator output. Robust parameter estimation and output prediction algorithms Following the work done in (Brdys and Chang (2002); Chang (2003)) the algorithms for the extended PPM approach were derived in the object oriented programming framework. There are three key algorithms underlying the main one - PPM estimator, that is: structure determination algorithm, parameter estimation algorithm and output prediction algorithm. Structure determination algorithm Since the plant considered within this paper is a CIS considered to have network structure, the first thing to do is the to determine the paths that contributed into building an output at current time instant. For that purpose a recursive path tracking algorithm is exploited based on (Zierolf et al. (1998); Shang et al. (2000, 2002); Chang (2003)). Parameter estimation algorithm This algorithm is used for parsimonious piece-wise robust estimation of the parameters along the modelling horizon. Within this calculation routine two optimisation tasks, announced previously
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
in Section 2, are solved. It constitutes the core optimisation task designed to obtain the information about parameter envelopes and is given by: [ l u l u ]∗ { ( )} θ ,θ ,ϵ ,ϵ = arg min J θl , θu , ϵl , ϵu l u l u [θ ,θ ,ϵ ,ϵ ] ( ) (12) s.t. (θnc , ϵnc ) 2 Θ θl , θu , ϵl , ϵu
where (for the clarity of presentation): θl , θu , ϵl , ϵu are the lower and upper bounds on parameter vector and structural uncertainty respectively; pair (θne , ϵne ) represents the coefficients that explain the ne -th ( ) experiment outputs; the cost function J θl , θu , ϵl , ϵu is set up to assess the distance between the parameter and structure error bounds as: ( ) ( ) ( ) J θl , θu , ϵl , ϵu = (θu θl) P ( θu θl) + (13) + ϵu ϵl Q ϵu ϵl where: P and Q are assumed to be diagonal positive defined matrices. The optimization algorithm is repeated within each timeslot to obtain the general PPM structure for the whole modelling horizon Tm . Since matrices P and Q weight the impact of the uncertainty contained in the parametric and structure error terms on the overall estimation and prediction that is why choosing the appropriate values of those matrices has a significant influence on the result of given optimisation task. The optimal weight determination lays outside the scope of this work. The second optimisation task is utilised in order to verify the quality of the parameter estimation for the MPC control purposes and is given by: { } W = max Ypu Ypl (14) where: W is the model uncertainty radius which serves as a measure of model uncertainty impact on the output prediction which is calculated under the assumption that the input u (t) = H (t); the bounds Ypu and Ypl are defined as: [ ] Ypu = [ypu (k + 1jk) , ..., ypu (k + Tm jk)] ; (15) Ypl = ypl (k + 1jk) , ..., ypl (k + Tm jk) ; where ypu (k + ijk) and ypl (k + ijk) are calculated according to the following scheme: ypu (k + ijk) = max fy (k + ijk)g Vp,k+i
s.t. y (k + ijk) 2 Ypi l yp (k + ijk) = min fy (k + ijk)g
(16)
Vp,k+i
s.t. y (k + ijk) 2 Ypi where the set Vp,k+i is given by: { } y (k + 1jk) , ..., y (k + Tm jk) , Vp,k+i = φ (k + 1jk) , ..., φ (k + Tm jk) , θ, ϵ (k + 1jk) , ..., ϵ (k + Tm jk)
(17)
and to complete Ypi is defined as: y (k + ijk) 2 R : y (k + ιjk) = θφ (k + ιjk) + ϵ (k + ιjk) , ∆ Ypi = θ 2 [θl , θu ] , ϵ (k + ιjk) 2 [ϵl , ϵu ] , (18) φ (k + ιjk) 2 [φl , φu ] , ι 2 [1, i]
Output prediction algorithm The robust output prediction is acquired by utilising (15)-(18). 429
Fig. 2. Model structure The full optimisation task is highly computationally demanding due to very high dimensionality of the problem obtained from wrapping the time horizon into a single point in time-space. In order to perform the optimisation a significant number of input scenarios is required. This however may be reduced as in (Chang (2003)) by a proper experiment design utilising a set of ’rich’ input signals. The same technique was utilised within this work, however due to strong nonlinear character of the implicit model an adequate handling of the model initialisation within experiments was highly required. 3.2 Simulation models Two simulation models are utilised within this paper. First is the plant model utilised to simulate the processes occurring within this DWDS. It is a unique quality model derived in the work (Arminski and Zubowicz (2011)), recently submitted as (Arminski et al. (2013)), to account for the normal and disturbance operational states. A single space-point scheme representing the reaction part of the model has been illustrated in Fig. 2. This complex DWDS quality model in general consists of two general modules first accounting for the normal and disturbance operational states (see Fig. 1). This firs module is divided into two sub-modules, modelling the chemical and biological processes respectively and a boundary biochemical layer constituting interactions between the modules. The chemical part of the model consists of the chloramine decay model and its interactions with, nitrogen, bromine, carbon and organic compounds. The later results in formation of disinfection by-products (DBPs). The biological module consists of heterotrophic bacteria and two step nitrification mechanisms conducted by nitrite and nitrate bacteria species. The interaction layer consists of the substances influencing the bacteria life-cycle, including inhibitory influence of the Chlorine species, and the products of their growth undergoing aqueous biochemical reactions within DWDS. Second model is utilised within the simulator of the PPM estimator/predictor mechanisms. Within this work a nonlinear structure of quality processes, tuned according to the previously introduced full model, has been utilised (Gallard and von Guten (2002); Clark and Sivaganesan (2002)). The reactive part yields:
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
Fig. 3. Network layout ( ) { 1 dt Cl = ΞCl = kCl Cl sDBP kDBP DBP)p1 DBP ( 2 dt DBP = ΞDBP = kDBP DBPp2 DBP Cl (19) where: Cl = Cl(t) is disinfectant concentration; DBP = DBP (t) is DBP concentration in mgCl; DBP p is DBP forming potential assumed to be constant; kcl , kDBP are the kinetic constants; sDBP is stoichiometric constant.
Fig. 4. Path delay variability 24h horizon
The advantage of introducing the nonlinear dynamics within the PPM simulator is the structure error reduction resulting from more precise disinfectant count. Since both the structure and algorithms are in place for the extended PPM estimator/predictor an illustrative numerical example can be given. 4. NUMERICAL EXMAPLE To illustrate the work-flow of the described technology a widely known EPANET DWDS case study plant Rossman (2000) has been selected ( Fig. 3). From the input – output layout and hydraulic analysis of DWDS (see Fig. 3) it is clear that one may distinguish 3 distinct paths on 24h prediction horizon, which are as follows: P 1 : 122 112; P 2 : 122 21 111; P 3 : 31 121 111. The corresponding time delays are given in Fig. 4. In Fig. 5 an example of robust parameter prediction of point parametric model was presented while in Fig. 6 a robust output prediction of chlorine concentration at a monitoring node (see the red marker in Fig. 3). The vertical red line in Fig. 5 denotes the current time instant so the upper plot represents the 24h prediction horizon calculated at time instant being the beginning of simulation experiment while the lower plot illustrates the same however at time instant being six hours after the start of the experiment. In order to obtain the results the combined MATLAB and EPANET 2.0 (Rossman (2000)) along with MSX extension (Shang et al. (2008)) simulation environments have been utilised. 5. CONCLUSION In this work an extension of the point-parametric model framework given in (Brdys and Chang (2002); Chang (2003); Chang et al. (2004)) to a class of reactive carrierload nonlinear dynamic networked system has been shown. 430
Fig. 5. Robust parameter estimations 24h horizon The goal was to derive a simplified model structure suitable for on-line monitoring and control carried out utilising the model predictive control scheme. The highly nonlinear spatially distributed structure was reintroduced by utilising the piece wise time-varying parameter linear structure. The applied point-parametric model approach resulted in obtaining robust information about the parameters and outputs. This can be utilised to perform robust monitoring (one step prediction) or for model predictive control purposes (prediction on a given horizon). A numerical example of drinking water distribution system application illustrates the capabilities of the technology. ACKNOWLEDGEMENTS This work was supported by the EU COST ACTION IC0806: Intelligent Monitoring, Control and Security of Critical Infrastructure Systems (IntelliCIS) and Polish MNiSW No. 638/N COST/09/20/2010/0: Intelligent systems for monitoring, control and security of critical infrastructure plants: methodology, structures, algorithms and applications to drinking water distribution networks (InSIK). The authors wish to express their thanks for the support.
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
Rossman, L. (2000). EPANET 2 User‘s manual. Risk Reduction Engineering Laboratory, U.S. Environmental Protection Agency, Cincinnati, OH, 45268, USA. Shang, F., Uber, J., and Polycarpou, M. (2000). Inputoutput model of water quality in water distribution systems. In A.S.C.E. 2000 Joint Conference on Water Resources Engineering and Water Resources Planning and Management. Minneapolis, Minnesota, USA. Shang, F., Uber, J., and Polycarpou, M. (2002). Particle backtracking algorithm for water distribution system analysis. Journal of Environmental Engineering, 128(5), 441–450. Shang, F., Uber, J., and Rossman, L. (2008). EPANET Multi Specie Extension Users Manual. National Risk Management Research Laboratory, National Homeland Security Research Center Office of Research and Development, U.S. Environmental Protection Agency, Cincinnati, OH, 45268, USA. Zierolf, M., Polycarpou, M., and Uber, J. (1998). Development and auto-calibration of an input-output model of chlorine transport in drinking water distribution systems. IEEE Transactions on Control Systems Technology, 6(4), 543–553.
Fig. 6. Robust output prediction 24h horizon REFERENCES Arminski, K. and Zubowicz, T. (2011). Multispecies quality model for drinking water distribution system. Technical report, Gdansk University of Technology, Gdansk, Poland. Arminski, K., Zubowicz, T., and Brdys, M. (2013). Biochemical multi-specie quality model of drinking water distribution system for simulation and design. Submitted to International Journal of Applied Mathematics and Computer Science. Brdys, M. and Chang, T. (2002). Robust model predictive control under output constraints. In IFAC 15th Triennial World Congress. Barcelona, Span. Chang, T. (2003). Robust Model Predictive Control of Water Quality in Drinking Water Distribution Systems. Ph.D. thesis, University of Birmingham. Chang, T., Duzinkiewicz, K., and Brdys, M. (2004). Bounding approach to parameter estimation without priori knowledge on model structure error. In IFAC 10th Symposium Large Scale Systems: Theory and Applications, volume 9. Osaka, Japan. Clark, R. and Sivaganesan, M. (2002). Predicting chlorine residuals in drinking water: Second order model. Water Resour. Plann. Manage., 128(2), 152–161. Gallard, H. and von Guten, U. (2002). Chlorination of natural organic master: kinetics of chlorination and of THM formation. Water Research, 36, 65–74. Langowski, R. and Brdys, M. (2007). Monitoring of chlorine concentration in drinking water distribution systems using an interval estimator. International Journal of Applied Mathematics and Computer Science, 17(2), 199–216. Maciejowski, J. (2002). Predictive Control with Constraints. Prentice Hall, London. 431
Robust parameter estimation and output prediction for reactive carrier-load nonlinear dynamic networks Mietak A. Brdys ∗,∗∗ Tomasz Zubowicz ∗ Krzysztof Arminski ∗ ∗
Department of Control Systems Engineering Gdask University of Technology, ul. G. Narutowicza 11/12, 80-233 Gdask, Poland (e-mail: [email protected]; [email protected]; [email protected]) ∗∗ Department of Electronic, Electrical and Computer Engineering, College of Engineering and Physical Sciences University of Birmingham, Edgbaston, Birmingham B15 2TT, UK (e-mail: [email protected]) Abstract: In this paper an extension of on-line model simplification technique for a class of networked systems, namely reactive carrier-load nonlinear dynamic networked system (RCLNDNS), kept within point-parametric model (PPM) framework is addressed. The PPM is utilised to acquire a piece wise constant time-varying parameter linear structure for the RCLNDNS suitable for the on-line one step ahead prediction that may be applied to monitoring and model predictive control purposes. The advantageous structure of PPM offers a considerable simplification of those algorithms in terms of computational burden in comparison to the nonlinear structures. Moreover, the availability of technical tools makes the analysis of the final form algorithms straightforward. An application to water quality estimation/prediction in drinking water distribution system illustrates the technology. Keywords: networked systems, robust parameter estimation, prediction 1. INTRODUCTION Development of a modern world has contributed to enlarging the number of large scale complex systems (LSCS) which deliver fundamental goods driving the existence and sustainable development of modern societies. These are called the critical infrastructure systems (CIS). Most of the CIS are characterised by complex, highly interconnected, networked structure e.g. power grid, drinking water distribution system, telecommunication networks etc. The first two constitute a very specific type of networked systems which can be described as reactive carrier-load networked system where one can distinguish a medium carrying a load of certain quality which can be altered by some mechanisms resulting from occurring processes. These processes result from given system operation and in majority of cases are guided by a nonlinear dynamics. Henceforth, the system at hand, to be tackled with, is a reactive carrier-load nonlinear dynamic networked system (RCLNDNS). Monitoring, control and security of RCLNDNS-CIS constitute timely as well as highly important research topic having strong impact on daily life of the societies worldwide. Moreover, it requires appropriate structures and algorithms in order to achieve reliable (in terms of robustness) solutions. In this work one of the most fundamental CIS is considered as RCLNDNS application example, namely the drinking water distribution system (DWDS). 978-3-902823-39-7/2013 © IFAC
426
The DWDS has been selected since not only that it is well fitted example of RCLNDNS but also it is characterised by multiple time-scale dynamics not only in the context of transportation-reaction part but also in terms of biochemical processes constituting the reactive part of the system. Furthermore, available measurements in terms of their number and diversity are very limited resulting in poor access to process information. A robust monitoring of DWDS has been addressed in (Langowski and Brdys (2007)) by applying an interval observer technology based on cooperative systems theory. Control of RCLNDNS have for the long time now been addressed in the model predictive algorithm framework (MPC). The MPC family constitutes of several mainstream approaches: richly referenced e.g. within (Maciejowski (2002)). Although, nonlinear MPC algorithms have been studied extensively in recent years it still comes as a fact, that nonlinear optimisation task that needs to be solved to find an optimal/suboptimal solution are very computationally demanding tasks, which is why having a neat linear model is still attractive in many applications. Handling such a model in manner to account for changing operational conditions, uncertainty in both structure and parameters most certainly requires proper mechanisms. In previous works (Brdys and Chang (2002); Chang (2003); Chang et al. (2004)) authors introduced a pointparametric model (PPM) approach to deliver such a model 10.3182/20130708-3-CN-2036.00026
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
of linear structure and time-varying parameters as well as algorithms for its parameter robust estimation and in consequence robust representation of DWDS process dynamics. The structure of the PPM was derived based on the modified path-tracking algorithm (Zierolf et al. (1998); Shang et al. (2000, 2002)) in order to account for the storage facilities within the network. In this approach a bound model for the uncertainty was utilised. The PPM framework was built and verified using linear quality model in the reactive part of the system. This paper extends the concept of PPM framework from the linear networked systems to RCLNDNS. To achieve that, a more sophisticated implicit simulation model is introduced within the PPM estimator/predictor structure in order to reduce the structural error. Moreover, the verification process utilises also the enhanced, nonlinear, model for plant simulation briefly introduced in Section 3. The paper is organised as follows. The preliminaries are given in forthcoming Section 2. Section 3 introduces and characterises the structures and algorithms utilised for robust estimation and prediction of PPM parameters and outputs. In Section 4 a numerical example based on unique DWDS model to illustrate the technology has been shown. Finally, Section 5 concludes the paper. 2. PRELIMINARIES Consider a nonlinear dynamic advective-reactive transport processes in axial direction within the network environment: { vjp ∂xjp Ψijp + Ξip , bulk flow i (1) ∂t Ψj p = Ξip , link wall ( ) ∆ where: t is time; Ψjp i = Ψjp i t, xjp is the i-th concen( ) ∆ tration in the jp -th pipe; Ξip = Ξip Ψjp , vjp , sjp are the rates dependent on all process variables; sjp is the bulk flow contact surface with the pipe wall; D is the dispersion coefficient; i and jp denotes the state variables and pipes of the network respectively. The processes may be described in general as: ( ) ( ) Ξp Ψjp , vjp , sjp = Arp Ψjp , vjp , sjp (2) ( ) where: Ξip Ψjp , vjp , sjp are the reaction rates and Ξp = [ 1 2 ]T Ξp , Ξp , ..., Ξnp where n is the total number of processes and ( state variables ) as well; A is the coefficient matrix; rip (Ψ)jp , vjp , sjp is the process rate vector and rp = [ 1 2 ]T rp , rp , ..., rpz where z denotes the processes involved in the model. Note that the interaction between the bulk and wall species is integrated into the Ξip . The storage facilities may be described as:( ) d Ψil Vl = Fin Ψiinl Fout Ψil + Vl Ξil (3) dt where: Fin , Fout are the inflow and outflows respectively; Ψiinl is the concentration of the i th specie in the inflow; ∆
Vl = Vl (t) is the volume of the medium contained in ∆ the l th tank; Ξil = Ξil (Ψl ) are the quality processes dependent on all species concentration and is assumed that may be different than Ξip ; l is the tank number. 427
The definition and mechanisms related to Ξit are given analogously. It is assumed that the robust information about an arbitrary, measurement accessible, variable of the selected process may be obtained at a given point in space utilising the simplified parameter time-varying model structure: y (k) = θ (k) φ (k 1) + ε (k) (4) where: y (k) is the output; θ (k) and ε (k) represent model parameter vector and structure error and disturbance impact respectively and furthermore constitute the uncertainty within the model; φ (k 1) is the regressor vector; k 2 [k0 , k0 + Tm ]Z denotes time where k0 is current time instant and Tm is a modelling horizon length; Z yields the integer field. The internal structure of (4) results from the discrete MISO model structure: yjo (k) = +
nT ∑
i=1 nI ∑
bjo ,i (k) yjo (k ∑
i)
ajo ,jI ,i (k) ujI (k
i) + εjo (k)
(5)
jI =1 i∈I
∧jo 2 1, nM where: yjo is the jo -th monitored node output and jo 2 1, nM and nM is the total number of monitored nodes; ujI is the jI -th plant input and jI 2 1, nI and nI is the total number of inputs contributing to the output yjo ; bjo ,i and ∆
ajo ,jI ,i are time-varying model parameters; I = [il , ..., iu ] and il and iu define range of delays in the control input; and the constructions are as follows: θjo (k) = [θjo ,1 , ..., θjo ,n ] = ] [ bjo ,1 (k), ..., bjo ,nT (k), (6) = ajo ,1,il (k), ..., ajo ,1,iu (k) ..., ajo ,nI ,il (k), ..., ajo ,nI ,iu (k) ∧n = nT + nI (k ijo ,nI ,u ) ] [ yjo (k 1) , ..., yjo (k nT ) , φjo (k 1) = u1 (k ijo ,1,l ) , ..., u1 (k ijo ,1,u ) , ..., unI (k ijo ,nI ,l ) , ..., unI (k ijo ,nI ,u ) (7) Moreover, if the input u[k0 ,k] () is valued on a compact set over the time horizon [k0 , k] then there exist an upper and lower bound on the trajectories of θ (k) and ϵ (k); θu,0 (k), θl,0 (k) and ϵu,0 (k), ϵl,0 (k) respectively, which guarantee that the model output y[k0 ,k] () is equal to p p the plant output y[k (), y[k0 ,k] () = y[k () and in 0 ,k] 0 ,k] particular k = k0 + Tm . Furthermore, it is assumed that the information about the bounds on the input uncertainty umin u(k) umax are known a priori. It is also assumed that a good approximation of the tightest bounds on parameter values can be made at time t0 as in fact they are not known. Although θ (k) can change abruptly its value at the time instant kj 2 [k0 , k0 + Tm ], j 2 1, Λ, the changes may be considered slow over a shorter time horizon [kj−1 , kj ]. The time instances kj , to account for the uncertainty in the plant transportation process, are to be found in an online manner. Therefore, the modelling time horizon Tm is partitioned into Λ time-slots (9) which in fact change their numerical amount with time thus Λ should be in fact
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
replaced with Λ(k) resulting in the following definition of the j-th time-slot: { } ∆ Sj = k 2 Z : kj−1 k kj , j 2 1, Λ(k) (8) Moreover, not all parameters contained within θ (k) have non-zero values. Hence, an active parameter sub-vector { } may be defined θj = θij ij ∈Ij where Ij is a set of active parameter indexes over Sj . Furthermore, θ (k) = θj for is assumed to be constant. Therefore, a significant reduction in the bound determination problem complexity may be gained. In consequence, the orthotope being searched for is given as: ( j ) θ , ϵ (k) 2 Θ, nc nc (θ , ϵ ) 2 Rn : n n c c =[ θ φ](k 1) + ϵ (k) , ( l u l u ) ∆ yn(k) Θ θ , θ , ϵ , ϵ = θ c 2 [ θl , θu] , nc l u ϵ 2 ϵ , ϵ , ne = 1, Ne , k 2 [k0 , k0 + Tm ] (9) The final time-space point parametric piecewise model structure is as follows: y (k) = θ (k) φ (k 1) + ϵ (k) θ( (k) = θ)j ∧ k 2 Sj ∧ j 2 1, Λ θj , ϵ (k) 2 Θ
(10)
The proof verifying the existence of such a model can be found in (Chang (2003)). In order to ensure the PPM to be suitable for robust MPC purposes the uncertainty radius must be kept sufficiently small. On the other hand partitioning the modelling horizon into slots will most certainly increase the width of the estimated parameter envelopes; however the piece-wise continuity may be attractive to handle by MPC internal algorithms. Finding a suitable compromise between those two is a problem on its own and will not be addressed within this work. The robust parameter and output estimation is achieved throughout solving dedicated optimisation tasks defined in the forthcoming Section 3. 3. DWDS CASE STUDY 3.1 Implementation of PPM framework nc Let y[k () denote the model response to an input 0 ,k] nc u[k0 ,k] () over a time horizon [k0 , k + Tm ]. Then there nc exists a trajectory θnc () that results in y[k () = 0 ,k] p nc y[k0 ,k] (). Different u[k0 ,k] () scenarios will require different parameter scenarios, therefore: y nc (k) = θnc (k) φ (k 1) + ϵnc (k) (11) where: ne is the input scenario number and ne = 1, ..., Ne where Ne is the total number of scenarios.
The total number of input scenarios required to find the parameter bounding envelopes, in the direct approach, should be equal to a total number of possible input combinations that is Ne = 2N i where Ni is the total number of input time horizon samples corresponding to output observations on [k0 , k0 + Tm ]. However, in (Chang (2003)) 428
Fig. 1. PPM information structure it has been shown that a much smaller population of welldesigned experiments is sufficient to obtain satisfactory results. Thus Ne will be limited within this paper to a number of required experiments as proposed in (Chang (2003)). Structure of information exchange In order to utilise the modified PPM framework a proper information structure between the subsystems is required – see Fig. 1. The initial information about process parameters and measurements is utilised by the Calibrator subsystem in order to feed the data into the implicit model Simulator and build up the Plant information. The Plant information enriches this information with the disturbance predictions and plant structural information. Appropriate data packages are then redistributed between the internal Simulator, Explicit Model Structure Determination unit and Estimator. The Estimator sub-module utilises the explicit model simulator in order to acquire the data for robust parameter prediction based on previously defined structure (10) build up within Model Structure Determination unit. Based on the calculated model parameter envelopes the Explicit model sub-module produces the robust output prediction and a nominal one assumed to be a Chebyshev centre. The information about the parameter and output robust predictions (as well as the nominal one) is fed to the PPM estimator output. Robust parameter estimation and output prediction algorithms Following the work done in (Brdys and Chang (2002); Chang (2003)) the algorithms for the extended PPM approach were derived in the object oriented programming framework. There are three key algorithms underlying the main one - PPM estimator, that is: structure determination algorithm, parameter estimation algorithm and output prediction algorithm. Structure determination algorithm Since the plant considered within this paper is a CIS considered to have network structure, the first thing to do is the to determine the paths that contributed into building an output at current time instant. For that purpose a recursive path tracking algorithm is exploited based on (Zierolf et al. (1998); Shang et al. (2000, 2002); Chang (2003)). Parameter estimation algorithm This algorithm is used for parsimonious piece-wise robust estimation of the parameters along the modelling horizon. Within this calculation routine two optimisation tasks, announced previously
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
in Section 2, are solved. It constitutes the core optimisation task designed to obtain the information about parameter envelopes and is given by: [ l u l u ]∗ { ( )} θ ,θ ,ϵ ,ϵ = arg min J θl , θu , ϵl , ϵu l u l u [θ ,θ ,ϵ ,ϵ ] ( ) (12) s.t. (θnc , ϵnc ) 2 Θ θl , θu , ϵl , ϵu
where (for the clarity of presentation): θl , θu , ϵl , ϵu are the lower and upper bounds on parameter vector and structural uncertainty respectively; pair (θne , ϵne ) represents the coefficients that explain the ne -th ( ) experiment outputs; the cost function J θl , θu , ϵl , ϵu is set up to assess the distance between the parameter and structure error bounds as: ( ) ( ) ( ) J θl , θu , ϵl , ϵu = (θu θl) P ( θu θl) + (13) + ϵu ϵl Q ϵu ϵl where: P and Q are assumed to be diagonal positive defined matrices. The optimization algorithm is repeated within each timeslot to obtain the general PPM structure for the whole modelling horizon Tm . Since matrices P and Q weight the impact of the uncertainty contained in the parametric and structure error terms on the overall estimation and prediction that is why choosing the appropriate values of those matrices has a significant influence on the result of given optimisation task. The optimal weight determination lays outside the scope of this work. The second optimisation task is utilised in order to verify the quality of the parameter estimation for the MPC control purposes and is given by: { } W = max Ypu Ypl (14) where: W is the model uncertainty radius which serves as a measure of model uncertainty impact on the output prediction which is calculated under the assumption that the input u (t) = H (t); the bounds Ypu and Ypl are defined as: [ ] Ypu = [ypu (k + 1jk) , ..., ypu (k + Tm jk)] ; (15) Ypl = ypl (k + 1jk) , ..., ypl (k + Tm jk) ; where ypu (k + ijk) and ypl (k + ijk) are calculated according to the following scheme: ypu (k + ijk) = max fy (k + ijk)g Vp,k+i
s.t. y (k + ijk) 2 Ypi l yp (k + ijk) = min fy (k + ijk)g
(16)
Vp,k+i
s.t. y (k + ijk) 2 Ypi where the set Vp,k+i is given by: { } y (k + 1jk) , ..., y (k + Tm jk) , Vp,k+i = φ (k + 1jk) , ..., φ (k + Tm jk) , θ, ϵ (k + 1jk) , ..., ϵ (k + Tm jk)
(17)
and to complete Ypi is defined as: y (k + ijk) 2 R : y (k + ιjk) = θφ (k + ιjk) + ϵ (k + ιjk) , ∆ Ypi = θ 2 [θl , θu ] , ϵ (k + ιjk) 2 [ϵl , ϵu ] , (18) φ (k + ιjk) 2 [φl , φu ] , ι 2 [1, i]
Output prediction algorithm The robust output prediction is acquired by utilising (15)-(18). 429
Fig. 2. Model structure The full optimisation task is highly computationally demanding due to very high dimensionality of the problem obtained from wrapping the time horizon into a single point in time-space. In order to perform the optimisation a significant number of input scenarios is required. This however may be reduced as in (Chang (2003)) by a proper experiment design utilising a set of ’rich’ input signals. The same technique was utilised within this work, however due to strong nonlinear character of the implicit model an adequate handling of the model initialisation within experiments was highly required. 3.2 Simulation models Two simulation models are utilised within this paper. First is the plant model utilised to simulate the processes occurring within this DWDS. It is a unique quality model derived in the work (Arminski and Zubowicz (2011)), recently submitted as (Arminski et al. (2013)), to account for the normal and disturbance operational states. A single space-point scheme representing the reaction part of the model has been illustrated in Fig. 2. This complex DWDS quality model in general consists of two general modules first accounting for the normal and disturbance operational states (see Fig. 1). This firs module is divided into two sub-modules, modelling the chemical and biological processes respectively and a boundary biochemical layer constituting interactions between the modules. The chemical part of the model consists of the chloramine decay model and its interactions with, nitrogen, bromine, carbon and organic compounds. The later results in formation of disinfection by-products (DBPs). The biological module consists of heterotrophic bacteria and two step nitrification mechanisms conducted by nitrite and nitrate bacteria species. The interaction layer consists of the substances influencing the bacteria life-cycle, including inhibitory influence of the Chlorine species, and the products of their growth undergoing aqueous biochemical reactions within DWDS. Second model is utilised within the simulator of the PPM estimator/predictor mechanisms. Within this work a nonlinear structure of quality processes, tuned according to the previously introduced full model, has been utilised (Gallard and von Guten (2002); Clark and Sivaganesan (2002)). The reactive part yields:
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
Fig. 3. Network layout ( ) { 1 dt Cl = ΞCl = kCl Cl sDBP kDBP DBP)p1 DBP ( 2 dt DBP = ΞDBP = kDBP DBPp2 DBP Cl (19) where: Cl = Cl(t) is disinfectant concentration; DBP = DBP (t) is DBP concentration in mgCl; DBP p is DBP forming potential assumed to be constant; kcl , kDBP are the kinetic constants; sDBP is stoichiometric constant.
Fig. 4. Path delay variability 24h horizon
The advantage of introducing the nonlinear dynamics within the PPM simulator is the structure error reduction resulting from more precise disinfectant count. Since both the structure and algorithms are in place for the extended PPM estimator/predictor an illustrative numerical example can be given. 4. NUMERICAL EXMAPLE To illustrate the work-flow of the described technology a widely known EPANET DWDS case study plant Rossman (2000) has been selected ( Fig. 3). From the input – output layout and hydraulic analysis of DWDS (see Fig. 3) it is clear that one may distinguish 3 distinct paths on 24h prediction horizon, which are as follows: P 1 : 122 112; P 2 : 122 21 111; P 3 : 31 121 111. The corresponding time delays are given in Fig. 4. In Fig. 5 an example of robust parameter prediction of point parametric model was presented while in Fig. 6 a robust output prediction of chlorine concentration at a monitoring node (see the red marker in Fig. 3). The vertical red line in Fig. 5 denotes the current time instant so the upper plot represents the 24h prediction horizon calculated at time instant being the beginning of simulation experiment while the lower plot illustrates the same however at time instant being six hours after the start of the experiment. In order to obtain the results the combined MATLAB and EPANET 2.0 (Rossman (2000)) along with MSX extension (Shang et al. (2008)) simulation environments have been utilised. 5. CONCLUSION In this work an extension of the point-parametric model framework given in (Brdys and Chang (2002); Chang (2003); Chang et al. (2004)) to a class of reactive carrierload nonlinear dynamic networked system has been shown. 430
Fig. 5. Robust parameter estimations 24h horizon The goal was to derive a simplified model structure suitable for on-line monitoring and control carried out utilising the model predictive control scheme. The highly nonlinear spatially distributed structure was reintroduced by utilising the piece wise time-varying parameter linear structure. The applied point-parametric model approach resulted in obtaining robust information about the parameters and outputs. This can be utilised to perform robust monitoring (one step prediction) or for model predictive control purposes (prediction on a given horizon). A numerical example of drinking water distribution system application illustrates the capabilities of the technology. ACKNOWLEDGEMENTS This work was supported by the EU COST ACTION IC0806: Intelligent Monitoring, Control and Security of Critical Infrastructure Systems (IntelliCIS) and Polish MNiSW No. 638/N COST/09/20/2010/0: Intelligent systems for monitoring, control and security of critical infrastructure plants: methodology, structures, algorithms and applications to drinking water distribution networks (InSIK). The authors wish to express their thanks for the support.
IFAC LSS 2013 July 7-10, 2013. Shanghai, China
Rossman, L. (2000). EPANET 2 User‘s manual. Risk Reduction Engineering Laboratory, U.S. Environmental Protection Agency, Cincinnati, OH, 45268, USA. Shang, F., Uber, J., and Polycarpou, M. (2000). Inputoutput model of water quality in water distribution systems. In A.S.C.E. 2000 Joint Conference on Water Resources Engineering and Water Resources Planning and Management. Minneapolis, Minnesota, USA. Shang, F., Uber, J., and Polycarpou, M. (2002). Particle backtracking algorithm for water distribution system analysis. Journal of Environmental Engineering, 128(5), 441–450. Shang, F., Uber, J., and Rossman, L. (2008). EPANET Multi Specie Extension Users Manual. National Risk Management Research Laboratory, National Homeland Security Research Center Office of Research and Development, U.S. Environmental Protection Agency, Cincinnati, OH, 45268, USA. Zierolf, M., Polycarpou, M., and Uber, J. (1998). Development and auto-calibration of an input-output model of chlorine transport in drinking water distribution systems. IEEE Transactions on Control Systems Technology, 6(4), 543–553.
Fig. 6. Robust output prediction 24h horizon REFERENCES Arminski, K. and Zubowicz, T. (2011). Multispecies quality model for drinking water distribution system. Technical report, Gdansk University of Technology, Gdansk, Poland. Arminski, K., Zubowicz, T., and Brdys, M. (2013). Biochemical multi-specie quality model of drinking water distribution system for simulation and design. Submitted to International Journal of Applied Mathematics and Computer Science. Brdys, M. and Chang, T. (2002). Robust model predictive control under output constraints. In IFAC 15th Triennial World Congress. Barcelona, Span. Chang, T. (2003). Robust Model Predictive Control of Water Quality in Drinking Water Distribution Systems. Ph.D. thesis, University of Birmingham. Chang, T., Duzinkiewicz, K., and Brdys, M. (2004). Bounding approach to parameter estimation without priori knowledge on model structure error. In IFAC 10th Symposium Large Scale Systems: Theory and Applications, volume 9. Osaka, Japan. Clark, R. and Sivaganesan, M. (2002). Predicting chlorine residuals in drinking water: Second order model. Water Resour. Plann. Manage., 128(2), 152–161. Gallard, H. and von Guten, U. (2002). Chlorination of natural organic master: kinetics of chlorination and of THM formation. Water Research, 36, 65–74. Langowski, R. and Brdys, M. (2007). Monitoring of chlorine concentration in drinking water distribution systems using an interval estimator. International Journal of Applied Mathematics and Computer Science, 17(2), 199–216. Maciejowski, J. (2002). Predictive Control with Constraints. Prentice Hall, London. 431