Journal of Environmental Management 155 (2015) 193e203

Contents lists available at ScienceDirect

Journal of Environmental Management journal homepage: www.elsevier.com/locate/jenvman

A generic methodology for the optimisation of sewer systems using stochastic programming and self-optimizing control Miguel Mauricio-Iglesias a, 1, Ignacio Montero-Castro a, Ane L. Mollerup a, b, Gürkan Sin a, * a

CAPEC-PROCESS, Department of Chemical and Biochemical Engineering, Technical University of Denmark, Building 229, Søltofts Plads, 2800 Lyngby, Denmark b HOFOR, Orestads Boulevard, 35, 2300 Copenhagen, Denmark

a r t i c l e i n f o

a b s t r a c t

Article history: Received 21 February 2014 Received in revised form 5 February 2015 Accepted 22 March 2015 Available online 31 March 2015

The design of sewer system control is a complex task given the large size of the sewer networks, the transient dynamics of the water flow and the stochastic nature of rainfall. This contribution presents a generic methodology for the design of a self-optimising controller in sewer systems. Such controller is aimed at keeping the system close to the optimal performance, thanks to an optimal selection of controlled variables. The definition of an optimal performance was carried out by a two-stage optimisation (stochastic and deterministic) to take into account both the overflow during the current rain event as well as the expected overflow given the probability of a future rain event. The methodology is successfully applied to design an optimising control strategy for a subcatchment area in Copenhagen. The results are promising and expected to contribute to the advance of the operation and control problem of sewer systems. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Sewer systems Combined sewers overflow Self-optimising control Control design Optimisation

1. Introduction Sewer systems are an essential part of urban water management which collect the sewage (both rain water as well as domestic wastewater) and transport it to centralized wastewater treatment plants for purification, before discharging it into surface waters. Control and operation of sewer systems is a challenging problem characterized by stochastic disturbances (rain events) and transient dynamics. In addition, the EU Water Framework Directive (European Commission, 2000) requires a reduction in the combined sewer overflows (CSOs) generated in urban environments. Improvement in the system operation performance appears as the suitable means to comply with the regulations. Major design modifications of the sewer systems are unlikely, owing to the high capital costs and organizational problems related to civil works in densely populated areas. Current research in sewer system control focuses mainly on optimisation of the control action (Duchesne et al., 2001; Mollerup et al., 2013; Ocampo-Martinez, 2010). Contributions range from

* Corresponding author. E-mail address: [email protected] (G. Sin). 1 Present address: Department of Chemical Engineering, University of Santiago de Compostela, Campus Sur, E-15782 Santiago de Compostela, Spain. http://dx.doi.org/10.1016/j.jenvman.2015.03.034 0301-4797/© 2015 Elsevier Ltd. All rights reserved.

efficient model development of the sewer systems (Breinholt et al., 2011), simplified models that can be quickly solved online (Vanrolleghem et al., 2005) and on their application by predictive control and online optimisers (Darsono and Labadie, 2007; Pleau et al., 2005). In addition to the aforementioned studies which focus on coordinating the control actions within the sewers network, several authors have widened the objective of sewer systems control by including the wastewater treatment plant and the water receiving body in a complex, generally multiobjective optimisation (Fu et al., 2008; Vanrolleghem et al., 2005). In contrast, the design of the regulatory layer in sewer system control has received little or no attention to the best of our knowledge. In this contribution, we adapt and apply the ideas of self-optimising control strategies to the design of a control structure (understood as the selection of controlled variables and their relation with manipulated variables) that will keep the system close to the control objective. Self-optimising control, first formulated by Skogestad (2000) and developed during the last 15 years by his team, is based on the selection of controlled variables (CVs) which, kept with constant setpoints, lead to an acceptable operation given a defined objective. In practice, this is achieved defining a loss function that depends on the CVs of the control structure and determines the effect of the disturbances on the control objectives. The self-

194

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

optimising structure corresponds to the set of CVs that minimises the loss function. Furthermore, the selected CVs can be measurements or linear combinations of measurements, which leads to a more robust control structure towards disturbances (Halvorsen et al., 2003). Hence, this method is particularly relevant for sewer systems. As the number of candidates for CV (i.e. the system measurements) generally exceeds the number of actuators, several measurements can be combined into a single CV providing robustness and more information to the controller. Self-optimising control is applied here to the design of sewer system control. Specific challenges in control of sewer systems are addressed, namely: i) as the future disturbances (rainfall) are stochastic and in general unknown for a long time horizon (>2 h), the optimal state of the system at a given time is difficult to define; ii) the transient dynamics, especially prevalence of unsteady state in sewer systems, prevent from defining one or several nominal operating points that can be used as a reference; iii) key parts of the sewer systems are the storage basins which, as integrators, cannot be evaluated at steady state (or very low input frequencies) and hinders the use of steady state design methods; and iv) the availability of too many measurements as alternative candidate for CVs, leads to a large design problem that, if tackled with optimisation tools, results into a large problem of combinatorial nature. The aim of this contribution is to present a generic methodology for designing self-optimising control structures in sewer systems. First, the methodology and the generic tools used are presented. Then, the methodology is illustrated through the application to a real case study: the design of a control system in a subcatchment area in the Copenhagen sewer system. All the algorithms and computer code used in this work are available from the authors upon request. 2. Methodology This methodology (Fig. 1) consists of 5 generic steps and considers that an appropriate state-space model of the sewer system is available as a starting point. The availability of a model is, in any case, the usual practice when designing a control and operation strategy. The methodology employs two important concepts, namely, a two-stage optimisation (see the explanation below in Section 2.2) and self-optimising control. The concept of selfoptimising control concept relies on defining a loss-function (Eq. (1)) that represents the difference between the optimum at nominal conditions and during the operation, i.e. subjected to disturbances, sensor noise and implementation errors.

L ¼ J  Jopt

(1)

The goal of self-optimising control is to design a control structure that minimises the loss function (L). It can be proved (Halvorsen et al., 2003) that minimising the worst case loss is equivalent to minimise the maximum singular value of the matrix M, which is defined as follows:

minL ¼

1 sðMÞ2 2

(2)

where M ¼ ðMd Mn Þ  1  1 2 Juu with Md ¼ Juu Jud  G1 Gd Wd and Mn ¼ Juu 2 G1 Wn 1

(3)

where the operator s returns the maximum singular value of a

matrix. Hence the methodology is articulated on two sections: definition of the optimal operation through a two-stage stochastic optimisation (system optimisation) and determination of the elements needed to define M (control structure design). These elements are the plant gain matrix G and the disturbance gain matrix Gd which are the ratios of the CVs w.r.t. the MVs and the disturbances respectively; the uncertainty matrix Wd which represents the maximum expected magnitude of each disturbance and Wn the implementation error for all the available measurements; Juu and Jud which are the second derivatives of the cost function with respect to the MV (u) and the disturbances (d). 2.1. Step 1: system analysis This step consists in the definition of the main characteristics of the system, the objectives and constraints. The analysis starts by identifying the main variables, namely: 1) The available measurements (Y), which may refer to those where a sensor is actually implemented (i.e. retrofitting) but also all those where a sensor can be installed (i.e. new design). 2) The manipulated variables (MV) which represent the control degrees of freedom of the system. In sewer systems, these are often the pumps. 3) The disturbances (d) including the rainfall (in most cases the main disturbance) and the dry weather flows. The user also sets the objective of the control system through the definition of a cost function (J). In sewer systems, minimising the total amount of overflow is often the main objective. Given the difficulty of forecasting rainfall, we define a cost function based both on the current and the expected overflow (Eq. (4) using a twostage optimisation problem formulation (Birge and Louveaux, 2011). Any other objective function can indeed be defined but it is likely that both current and expected terms will have to be considered and will have the same generic expression.



n X

w OFtot ðti Þ þ ð1  wÞ OFexp

(4)

i¼1

P The first term of the cost function ni¼1 OFtot ðti Þ is the total value of overflows generated in a simulation of the system for a given rain event. The second term OFexp is an estimation of the overflow caused by future rain events as defined in the next step. Finally, and for large systems, it may also be necessary to restrict the number of measurements to be used in control. In effect, if the number of available variables is too large, the size of the problem can become untreatable given the combinatorial nature of the problem. Besides, the resulting optimal control structure can be too complex, e.g. featuring many sensors. Hence, selecting the number of measurements available for control (nSubset) is a trade-off between the complexity of the controller (including sensor maintenance) and its performance. 2.2. Step 2: two stage stochastic optimisation The goal of this second step is to establish an optimal operation policy for the sewer system according to the defined objective function and a given rainfall scenario. Designed storms are useful tools to define a certain rainfall scenario according to the historical rain series and for a given return period, e.g. by the Chicago method (Keifer and Chu, 1957). After the definition of the rainfall scenario, the definition of the optimal operational policy consists of determining the value of the manipulated variables (MVs) that minimises the objective function. Minimisation of the objective function

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

195

Fig. 1. Workflow for the methodology structured in three tasks and six steps.

(Eq. (4)) represents a trade-off between the water storage to reduce overflows during the current rain event and keeping enough free capacity to handle future rain events. Hence, this optimisation is carried out as a two-stage optimisation (Birge and Louveaux, 2011) and formulated as follows:

2.2.2. Stochastic stage (offline)

2.2.1. Deterministic stage

arg min J ¼ u

tP sim

arg min OF ¼

w OFtot ðx; u; tÞ þ ð1  wÞ OFexp ðxÞ

u

t¼0

dx ¼ f ðx; u; tÞ dt y ¼ gðx; u; tÞ

s:t

(5)

n P

OFj ðti Þ

i¼1

dx ¼ f ðx; u; tÞ dt y ¼ gðx; u; tÞ

s:t

(7)

hðx; u; tÞ  0

hðx; u; tÞ  0 where x represents the states of the system, u represents the manipulated variables (pumping rates), w is the weight of each of the terms (in this work both terms are equally weighted and w ¼ 0.5). The functions f, g and h represent the evolution of the states, the relation between the outputs and the process, and the possible constraints in x,u or t. OFexp is defined by:

OFexp ðxÞ ¼ VaRa ðminOF ðx; s; uÞÞ

where s is a random variable representing the rainfall and VaRa is the value at risk at a certain probability level (i.e. the value of expected overflow that presents a probability a of being exceeded). Hence, each OF is individual solution of the following problem:

(6)

where j represents each of the potential sources of overflow. The solution of this problem follows these steps: i) the stochastic stage is solved offline determining the OF for a large number of different rain events, hence giving place to a probability distribution of overflow; ii) for a given probability threshold (e.g. at 0.95 probability level), the expected overflow becomes a function of the states of the system; and iii) the deterministic stage can be solved considering both the overflow during the current rain event as well as the probability of a future overflow. These steps are explained with detail below.

196

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

2.2.3. Solution of the stochastic stage The estimation of the expected overflow consists of deriving a relation between the state of the system and the probability of an overflow caused by a future rain event. The following steps are followed: 1) The states of the system are discretised in a grid in order to map the whole range of operation. The grid is a discretization of the model states and therefore on state scale For sewer systems, the degree of filling of each of the tanks is likely to be the most important variable. 2) A node of the grid of states is selected. 3) A large number of rain events are simulated. These rain events are sampled from a long series of historic rain for the catchment studied in order to reproduce the expected variability of every geographical zone. 4) The minimum overflow caused by each of these rain events is recorded. The determination of the minimum overflow is found by optimising the MVs for each sampled rain event. The observed overflow values for each rain event are recorded and used to estimate the empirical cumulative distribution function (see Section 3.2). 5) The selected node (i.e. a representation of a system state) is linked to a cumulative probability distribution of overflow and, for a given probability value, to an expected overflow (OFexp) see Section 3.2. Tasks 2e5 are repeated until all the states of the system are mapped

2.2.4. Solution of the deterministic stage The deterministic stage (Eq. (5)) consists in finding the value of the manipulated variables that minimises the overflow for the selected rain event, while leaving a certain free capacity in order to handle possible future rain events (i.e the expected overflow). Since both terms of the equation represent the same phenomenon (accumulated volume of overflow), both can be added. However, a user may decide to weight differently any of the terms, depending on the characteristics of the network or operational priorities. 2.3. Step 3: definition of operation phases and process gains In contrast with the original formulation of self-optimising control, aimed at continuous chemical processes, sewers systems are characterised by transient dynamics lacking a relevant steady state. Furthermore, the relation between the variables in the system is highly nonlinear, leading to different behaviours depending on the state of the system. Hence, to define completely the evaluation scenario it is necessary to characterise the range of operation of the system as follows: 1) The operating points are selected to characterise the whole range of operation of the sewer systems such as the dry weather, filling or emptying phases. A simulation of the system with the optimised MVs is run in order to define the system variables which define each of the system phases. 2) Basins, the key component in sewer systems control, present a transfer function with integrating modes. As a consequence, the steady state gain cannot be defined. Instead, a range of frequency (representing the speed of changing of the disturbance dynamics) must be defined that covers the range of operation of the control system, thereby leading to the process gains (i.e. the relations between disturbances, manipulated and controlled variables) For each of the operating phases (and corresponding states) defined in task 1, the system must be linearized. The linearization

provides the expressions to evaluate the plant gain matrix (Gy ) and the disturbance gain matrix (Gyd ) at the range of frequency previy ously defined. Besides, Gy and Gd are scaled, commonly using the maximum measurement values or the range of variation expected of each variable (Skogestad and Postlethwaite, 2005). Finally, other magnitudes appearing in Eq. (3) have to be defined. The Wn and Wd matrices are positive diagonal matrices in which each element is related to the uncertainty of a certain measurement or disturbance of the system. Therefore, both Wn and Wd will be squared matrices with size defined by the total number of measurements and disturbances considered. The cost function derivatives (Juu, Jud) are calculated using second order finitedifference to the objective function previously defined (Eq. (8)). Juu is a squared matrix with size equal to the number of manipulated variables of the system. Juu is a matrix with as many rows as manipulated variables and as many columns as disturbances considered.

Juu ¼

d2 J d2 J ; Jud ¼ 2 dudd du

(8)

2.4. Step 4: screening of measurements In large systems, the number of measurements can be very high, increasing the calculation load of the methodology without providing any additional improvement in the results or making the problem untreatable if it is very large. Therefore, a fraction of the measurements in the system can be ruled out without affecting the final solution. The size of the problem is thus reduced and the subsequent optimisation becomes manageable. First, a measurement is discarded when the manipulated variable does not have any influence in the measurement for the evaluation time analysed, reflected as a row of zeros in Gy . After this first screening, combinations of the resulting Gy with the dimensions fixed by nSubset, are generated. These sets of combinations are ranked following a certain standard. In this case, the minimum condition number was used as a controllability criterion (Skogestad and Postlethwaite, 2005). Under this criterion, gain matrices with high condition numbers, defined as the ratio of maximum and minimum singular values of the matrix (Eq. (9)), are a sign of poor controllability caused by interaction between loops. Therefore, only the measurements that form the set with the minimum condition number are kept for the following steps of the methodology.

CNðGÞ ¼ gðGÞ ¼

sðGÞ sðGÞ

(9)

2.5. Step 5: determination of controlled variables The matrix H defines the controlled variables (CVs) as linear combinations of the measurements as: y

G ¼ HGy ; Gd ¼ HGd

(10)

Hence, H relates the measurements (y) to the controlled variables.As an illustration, for a plant with three measurements and a single controlled variable, the latter would be obtained as:

0

1 y11 CV1 ¼ ðh11 h12 h13 Þ@ y21 A ¼ h11 y11 þ h12 y21þ h13 y31 y31 where the scalar hij are the coefficients of the linear combination of

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

measurements yji into controlled variables CVi. Hence, the determination of the self-optimising CVs is carried out finding H so that the loss function is minimised, which is equivalent to minimising the maximum singular value of the M matrix (Eq. (11)).

Hopt ≡ arg min s ðMÞ

(11)

H

This optimisation procedure is done for all the possible measurement combinations, selecting the one that provides the minimum loss value. The same procedure is followed for the different states of the system (emptying, filling …) thereby defining a control structure linked with a phase of the system operation. 2.6. Step 6: close loop evaluation and implementation A closed loop evaluation is needed to validate the results obtained previously. The users can define their own validation tests depending on the specific control structure of the sewer system. Two closed loop evaluations are proposed here in order to globally evaluate the results obtained by following this methodology: the performance during short and intense rain events and for long and mild ones. Given that different control structures are determined for each of the operation phases, two different implementations are possible. For systems where the range of operation is narrow, one control structure may be enough to handle the rain events. On the contrary, for systems where the range of operation is wide and characterised by strong nonlinearities, switching between the different control structures is more suitable although the controller becomes more complex. 3. Case study: a subcatchment area in Copenhagen, Denmark The developed methodology is now applied to a real case study: a subcatchment area of the Copenhagen sewer system (135 000 population equivalent). This system is modelled by the virtual tank (VT) approach (Ocampo-Martinez, 2010), which represents a subcatchment by a given surface area (receiving rainfall) and pipe volume. Hence, the pipe volume is lumped into a single virtual tank and the mass balance can be written as:

dV ¼ qin  qout ¼ qin  bV dt

(12)

where b is the volume/flow coefficient. Other elements include basins, represented by their level (L), pumps (P) and weirs (Fig. 2). The model of the subcatchment area was validated against a mechanistic complex model developed by the utility company HOFOR using MIKE Urban software (Mollerup et al., 2015). The model was implemented in Matlab Simulink and linearisations were performed using Matlab Simulink for different scenarios. All the calculations in the self-optimising control methodology were also implemented in Matlab. 3.1. Step 1: system analysis For this case study, all the measurements given by the model (Fig. 2) are considered. Although not all of them are actually implemented in the real sewer system, the economic cost associated to add extra sensors is not considered for this case study in order to evaluate the full potential of the methodology. Out of the two possible disturbances in this system, dry weather flow and the rainfall, only the rainfall will be assessed since the dry weather flow is negligible compared to the rainfall (when the objective considered is the CSO). The identified MVs are the pumping ratios for the three pumping stations with respect to their maximum capacities

197

(Table 1). The objective function is defined as depending on the current overflow OFtot, obtained through model simulation and an expected overflow OFexp, based on the state of the system (Eq. (4)). It is considered that the capacity of the system to handle a rain event depends on the degree of filling of tanks 2 and 3 (tank 1 has negligible volume in comparison). The objective function, then, provides a trade-off between the minimisation of the current overflow and the availability of storage in the tanks. The total number of measurements available for optimisation nSubset is fixed at 6, i.e. the maximum number of measurements that will be used for control. This choice depends on the user, on the desired complexity of the final control structure and on the possibility of sensor implementation. The number of measurements (nY) used to define the controlled variables will vary between 2 and 6 as it will be seen next. Depending on the number of nY considered, different results will be obtained, generating a range of control structures. 3.2. Step 2: two-stage stochastic optimisation Obviously the optimal policy depends on the scenario defined. Thus, a synthetic rain event was created using Chicago method which captures the relevant statistic characteristics of historical rain, in particular the maximum average intensities and the volume. Hence, it can be ensured that the system behaviour is well represented and displays the whole range of operation of the sewer system. In this particular case study, the rain event used is characterised by the following parameters: a return period of 5 years, a duration of 240 min and an average intensity of 630 mm/year, obtaining the resulting rain event (Fig. 3). The goal of the control system is to drive the system to the optimal policy according to a given objective (here minimise the current and expected overflow). Hence, defining the optimal policy for the selected rain event is equivalent to solve the optimisation set in Eq. (5), i.e. find the pumping rates that minimise the combination of overflow. First the stochastic stage (Eq. (7)) was solved offline. The degrees of filling (DF) of tanks 2 and 3 were taken as the most relevant states (tank 1 has negligible volume in comparison). They are defined as the ratio between the total volume of water in the tank and the total volume of the vessel (Eq. (13)).

DF ¼

VT;water VT;tank

(13)

Then, solving the stochastic stage consisted of relating the degrees of filling with the probability of overflow in the future. 1) The space of degrees of filling (DF2 and DF3) was discretized and sampled in 52 nodes by Latin Hypercube Sampling (LHS). 2) A node, i.e. a value of DF2 and DF3 from the LHS grid is selected. 3) 250 rain events are randomly sampled from the historical rain series in Copenhagen for the years 2011 and 2012. Dry weather is also included. Each of the rain events lasts 18 h which is the time it takes to the whole catchment area to discharge all the volume by the outflow (global emptying time). The rationale for selecting 18 h as the rain event time window is the following: for any rain event that starts 18 h after the previous one, the reservoirs would be empty, therefore, it is considered as independent of the previous one. The number of rain events was set to 250 as this sampling number was found sufficient to estimate the empirical distribution function for the expected overflows. 4) The overflow caused by each of these rain events is recorded for values of the pumping ratios (vP2 and vP3) that minimise the overflow generated each rain event (Eq. (14))

198

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

½vP2; vP3opt ¼ arg min OFexp ¼ vP2; vP3

n X ½UO17ðti Þ þ UO32ðti Þ þ UO38ðti Þ þ UO42ðti Þ þ UO44ðti Þ i¼1

dx ¼ f ðx; u; tÞ dt y ¼ gðx; u; tÞ

(14)

s:t

u≡vP2; vP3 2½0; 1

and ti 2½0; 18 h

The third manipulated variable, vP1, is fixed as a constant value of 0.717. This value is estimated as the fraction of the outflow from the catchment (Flout) that corresponds to the precipitation in the surface of the subcatchment area located upstream the pump number 1. This strategy weights evenly the flow delivered by the three pumps in relation to the downstream bottleneck (Flout in Fig. 2). Taking into account the maximum flow of discharge (Floutmax), the total area of the subcatchment area (Atot, syst) and the area that must be discharged by vP1 (AP1), the pump ratio is

calculated as:



L L 1100 Floutmax s ¼ 3:08 s ¼ 356:6 ha ha Atot;syst

L Fl3;std ¼ AP1 $r ¼ 645 /vP1 ¼ 0:717 s

Fig. 2. Schematic representation of the virtual tank model diagram for the Copenhagen sewer system used as case study.

(15)

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203 Table 1 Possible measurements in the subcatchment area. Type

Variable identification

Overflows Flows Add. flows Levels Pumping ratios

UO17, UO32, UO38, UO42, UO44 q1, q2, q3, q4, q5, q6, q7, ip1, ip2, ip3 Fl1, Fl2, Fl3, Fl4, F5, Flout L1, L2, L3 vP1, vP2, vP3

5) The 250 values of overflow recorded are used to build an empirical cumulative probability function (CDF, estimated using Matlab ecdf command) overflow for each of the nodes of the LHS grid (Fig. 4 shows the procedure for one node in the grid). The reproducibility of the empirical CDFs was checked by subsampling, i.e. re-estimate the CDF using a smaller number of LHS samples and perform a statistical test on the mean of two distributions to verify there is no statistical difference. For a probability level (here chosen as 0.05), the expected value of the overflow is determined from the empirical CDF. This value

199

represents the expected overflow that has a probability of 0.05 of being exceeded. Tasks 2e5 were then repeated for all the 52 nodes of the LHS grid. Basically the solution of this grid provides us with the offline solution of the stochastic stage of the problem. This means that for a given state in the system (defined by a 2-dimensional space), we can now infer about the expected probability of overflow OFexp. (at 95% probability). This expected overflow is then added to the objective function together with the value obtained from the deterministic part. For on-line implementation, an interpolation was used to obtain a continuous representation of OFexp (using Matlab function TriScatteredInterp). The results of the interpolation are displayed in Fig. 5. The deterministic stage consists now of solving Eq. (5) for u ≡ vP1, vP2, vP3 2 [0,1] and with the function OFexp determined by interpolation. Solving this equation led to optimal pumping rates of vP1 ¼ 0.467, vP2 ¼ 0.996 and vP3 ¼ 0.703. These are the nominal operating conditions for the system during a rain event. Likewise, a simulation of the system for these pumping rates provides the optimal trajectory for the design rain event. 3.3. Step 3: definition of process gains and operation phases

Fig. 3. Designed rain event for the Copenhagen sewer system subcatchment area used as disturbance for the control structure design stage.

As the optimal trajectories have been defined, the next step is to design the self-optimising controller that will drive the system close to the optimum. In order to do so, process gains relating all possible inputs and outputs are needed to define the matrices appearing in Eq. (3). However, the sewer system does not present a steady state operating point that can be used to evaluate G and Gd; instead, different operation phases must be defined in order to capture all the dynamics of the sewer system. In particular, some variables are inherently difficult to capture, as, e.g., the overflows, which are zero until the basins are full and then can grow very quickly. Five operating phases were defined in this system to model the whole range of operation: dry weather, filling, overflow from structures with little or no volume, overflow from structures with large volumes and emptying. A simulation of the rain episode (Fig. 6) was used to record the time at which each of these events took place, namely: 1) dry weather (t1 ¼ 10 min), 2) filling (t2 ¼ 50 min), 3) overflow from structures with little or no volume (t3 ¼ 120 min), 4) overflow from structures with large volumes

Fig. 4. Illustration of workflow to determinte the function OFexp(x). The CDF is determined from a large number of Monte Carlo simulation of rain events for each node in the LHS grid. The value at risk for a ¼ 0.05 is displayed as an example.

200

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

Fig. 5. Representation of the DFs space and the overflow values according to the stochastic optimisation stage for a value at risk of 0.05.

Fig. 6. Evolution of the main variables in the subcatchment area during simulation of the rain event. Naming of the overflow streams and levels is explained in Fig. 2.

(t4 ¼ 170 min) and 5) emptying (t5 ¼ 330 min). As the gains of integrating processes are not defined at steady state, a range of frequency was conservatively estimated as relevant for the operation. The lowest frequency was taken as the inverse of the expected period of a rain event (determined by time-frequency analysis of the series of two years of rain in the Copenhagen area) and corresponded to a frequency u2 ¼ 0.046 rad/min. The highest frequency corresponds to the fastest variations in the disturbances (rainfall) that will be propagated through the system. Since the virtual tanks are first-order systems, they act as low pass filters depending on their time constant. Hence, the highest bound is

taken as the highest break frequency of the virtual tanks, in this case VT 2 (the one with the lowest b as defined in Eq. (12)). This frequency is equal to u1 ¼ 0.59 rad/min. The gains of G and Gd remain relatively constant in this range of frequency; hence, only the slowest frequency is taken for the rest of the analysis. The expressions of G and Gd for the five operation phases were obtained using the linearisation routine from the control design tool (Matlab-Simulink). The linearisation was done as snapshots at the times recorded earlier (t1, … ,t5). The process gains were obtained evaluating G and Gd at u2 ¼ 0.046 rad/s. The uncertainty matrix Wn elements were defined as the 10% of

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203 Table 2 Table of measurements selected through the screening step. Operation phase

Second screening

1. 2. 3. 4. 5.

ip2, Fl2, Fl3, L1, L2, L3 ip2, Fl1, Fl2, Fl3, L2, L3 ip2, Fl1, Fl2, Fl3, L2, L3 UO44, ip1, ip2, ip3, Fl1, Fl4 ip2, Fl1, Fl2, Fl3, L2, L3

Dry weather Filling Overflow from structures with little or no volume Overflow from structures with large volume Emptying

the nominal value of each measured variable at the evaluation time. It represents the measurement error of the sensors. Wd matrix, which represents the expected magnitude of the disturbance, is a scalar value of 0.001 m3/s, which serves to scale the rainfall flow to a value close to unity. Finally, the cost function derivatives were obtained by numerical differentiation (second order central finite difference).

0

Juu

1:64 $104 ¼@ 294 1:14 $103

Jud

1 6:54 $108 ¼ @ 1:97$103 A 3:34 $107

294 3:91 $106 6:75 $107

1 1:14 $103 6:75 $107 A 2:46$103

0

3.4. Step 4: screening of measurements A screening of measurements was carried out to reduce the size of the problem. First, all the measurements that are not influenced by all the manipulated variables are ruled out (reflected as a row of zeros in Gy). This brought the number of measurements down from 24 to 10. Secondly, the remaining measurements are combined in sets of 6 (since nSubset ¼ 6) giving place to 210 sets of 6 variables

201

per operating phase. For each of the combinations, the condition number of the matrix (Eq. (9)) is calculated. Only the 6 variables in the set with the minimum condition number are kept for each operating phase. The whole screening procedure allowed to select a set of 6 variables from 134,596 sets, hence greatly decreasing the computational time required for the optimisation. The 6 selected measurements for each operation phase appear in Table 2, which shows that the operation phases 1, 2 and 3 result in the same selection of measurements. 3.5. Step 5: determination of controlled variables Minimization of the maximum singular value of M (Eq. (3)) led to matrix Hopt, which defines the dependency between the controlled variables and the measurements. For each operation phase, H was determined for an increasing number of measurements (from 2 to 6). It was seen that the value of the loss (Lopt) dropped suddenly when using four or five measurements; therefore, these are the simplest structures that provide a limited loss of the objective function, which is the deviation from the optimal operation point found in step 2. The resulting control structures can be found in the supplementary material. 3.6. Step 6: closed loop evaluation Different rain events were selected for this validation stage. Two rain events that characterise short and intense rain events (event A) and long rain events with low intensities (event B) were chosen from real rain events data for the city of Copenhagen (Fig. 7). The other two events (event A2 and B2) are modifications of event A and B respectively such that the rainfall intensity is not homogeneous in the subcatchment area. Instead, the areas which are located downstream (virtual tanks 3, 4 and 5) receive a more intense rainfall than the upstream areas (virtual tanks 1, 2, 6 and 7). Thus, virtual tank (VT) 5 receives a rain intensity which is 80% higher than in the non-modified rain event. Likewise, the rain

Fig. 7. Rain events A (above) and B (below) used for validation. The first 6 h of event are magnified in the figure to show the rain intensity.

202

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

Table 3 Accumulated overflow (m3) after four different rain events. SOC ¼ self-optimising control, OL ¼ open loop, PI ¼ PI loop structure based on the actual plant. Event/controller

SOC

OL

PI

Event Event Event Event 2011

1.60$105 1.36$105 2.24$104 8.25$103 1.93$105

1.60$105 1.36$105 2.30$104 8.96$103 2.06$105

1.61$105 1.37$105 3.73$104 9.50$103 2.69$105

A A2 B B2

intensity is 60% higher for VT4 and 40% higher for VT3; on the other hand VT 1, 2, 6 and 7 receive a 35% lower rainfall intensity, compared to the non-modified rain events. Events A2 and B2 aim at checking whether the SOC structure enables better feedback from the process thanks to the information provided by more measurements, and therefore can improve the control of the overflow. Finally, a more realistic validation (which did not rely on synthetic rain events) was done by simulating all the rain events corresponding to the year 2011 in Copenhagen. Based on the defined Hopt matrices (supplementary material), the designed controller is benchmarked with two other strategies. The designed self-optimising controller (SOC controller) is based on a switching structure that selects the Hopt depending on the operating phase of the system. The SOC controller is compared to a strategy without control (open loop, OL) with the pumping rates fixed at optimal values, and to a standard decentralised proportional-integral control structure (PI controller). The PI control structure resembles the one which is actually implemented in this subcatchment area with respect to pairing of measurements and controlled variables (supplementary material). 4. Discussion Designing the nominal pumping flows by a two-stage stochastic optimisation was proved as a suitable means to adjust the pumping policy to the sewer system, taking into account its particular features (bottlenecks, tank volume, typical rainfall pattern …). The results demonstrated that this step is essential for a sound

regulation of the overflows from the sewer systems, even in the absence of automatic control. Furthermore, if any system of automatic control is implemented, using self-optimizing control to identify the most suitable controlled variables and, in particular, combining the information provided by the measurements allows to design a control structure adapted to the aim of the system (in most cases avoiding CSO). According to the accumulated overflow of each of the three tested controllers during the validation tests (Table 3), it is seen that the self-optimising controller (SOC) has the best performance for the rainfall events. For both A events, it can be seen that the rainfall peak is simply too large for the system to handle it and the results are basically the same for the three strategies. On the contrary, the control has a large influence on the use of the storage capacity and therefore, it plays an important role during event B and the 2011 collection of rain events. Comparing the SOC with the OL shows that the difference is quite modest (ranging between 0 and 10\% of the total volume), which points out the importance of selecting a sound pump policy, well adapted to a given subcatchment area. The SOC represents a dramatic improvement compared to the PI (67 \% reduction for event B, 15 \% for event B2, 39\% for the 2011 events). For example, in event B it can be seen (Fig. 8) that the PI control system is not able to utilise all the storage capacity, and overflow takes place earlier than for the two other systems (min 350 compared to min 300). However, it must be noted that the PI is a simplified version of the actual control system. A sounder assessment about the performance of SOC compared to a realistic control system will be carried out in a forthcoming publication. Such validation would be the previous step to replacing the current control system, a complex rule-based control structure, with a decentralised structure where the setpoints are found by optimisation and the controlled variables are selected according to self-optimizing control principles.

5. Conclusions A novel and generic methodology was developed and evaluated for the optimization of sewer systems. The methodology features

Fig. 8. Total overflow during rain event A (above, time magnified around the peak rain) and rain event B (below).

M. Mauricio-Iglesias et al. / Journal of Environmental Management 155 (2015) 193e203

the following novel methods and tools:  A two-stage stochastic optimisation is successfully used to find the optimal operation points for the sewer systems. The twostage optimization accounted for both the overflow during the design rain event as well as for the expected overflow given the historical rain in the area.  To reduce the computational complexity, the stochastic optimization was first solved offline using Monte Carlo simulations (LHS used for sampling from rain events). For online implementation, the expected overflow was found using interpolation of the offline solution.  Self-optimizing control principle was successfully used to identify the optimal pairing of measurements with actuators.  A fast screening of the available measurements based on controllability criteria is successfully used so that only the most influential measurements are retained as candidates for control structure screening using self optimizing control principles.  The transient dynamics of sewer systems were addressed successfully by using an operation trajectory based on characteristic operation phases, e.g. dry weather, filling, emptying, etc.  The performance of sewer system operation using stochasticoptimization as well as self-optimizing control principles was found promising compared to PI and open-loop operation of the system. Last but not least, this generic optimization methodology is expected to contribute to rational based design and better management and mitigation of overflows in sewer systems. Acknowledgements The work is part of the Storm and Wastewater Informatics (SWI) project partly financed by the Danish Agency for Science, Technology and Innovation (09-063208). Authors would like to express their gratitude to Dr. Alberto Quaglia for useful advices on stochastic optimisation.

203

Appendix A. Supplementary data Supplementary data related to this article can be found at http:// dx.doi.org/10.1016/j.jenvman.2015.03.034.

References Birge, J., Louveaux, F., 2011. Introduction to Stochastic Programming. Springer, Berlin. Breinholt, A., Thordarson, F., Moeller, J., Grum, M., Mikkelsen, P., Madsen, H., 2011. Grey-box modelling of flow in sewer systems with state-dependent diffusion. Environmetrics 22, 946e961. Darsono, S., Labadie, J., 2007. Neural-optimal control algorithm for real-time regulation of in-line storage in combined sewer systems. Environ. Model. Softw. 22 (9), 1349e1361. Duchesne, S., Mailhot, A., Dequidt, E., Villeneuve, J., 2001. Mathematical modeling of sewers under surcharge for real time control of combined sewer overflows. Urban Water 3, 241e252. European Commission, 2000. Directive 2000/60/ec of the European Parliament and of the Council Establishing a Framework for the Community Action in the Field of Water Policy (Technical report). Fu, G., Butler, D., Khu, S.-T., 2008. Multiple objective optimal control of integrated urban wastewater systems. Environ. Model. Softw. 23 (2), 225e234. Halvorsen, I., Skogestad, S., Morud, J., Alstad, V., 2003. Optimal selection of controlled variables. Industry Eng. Chem. Res. 42, 3273e3284. Keifer, J., Chu, H., 1957. Synthetic storm patterns for drainage design. J. Hydraulics Div. 83. Mollerup, A.L., Grum, M., Muschalla, D., van Velzen, E., Vanrolleghem, P., Mikkelsen, P., Sin, G., 2013. Integrated control of the wastewater system: potential and barriers. Water 21, 15, 2, 39e41. Mollerup, A., Mikkelsen, P.S., Thornberg, D., Sin, G., 2015. Regulatory control analysis and design for sewer systems. Environ. Model. Softw. 66, 153e166, 10.1016/ j.envsoft.2014.12.001. Ocampo-Martinez, C., 2010. Model Predictive Control of Wastewater Systems. Springer, Berlin. e, P., Pelletier, G., Bonin, R., 2005. Global optimal realPleau, M., Colas, H., Lavalle time control of the quebec urban drainage system. Environ. Model. Softw. 20 (4 Spec. iss.), 401e413. Skogestad, S, 2000. Plantwide control: the search for the self-optimizing control structure. J. Process Control 10, 487e507. Skogestad, S., Postlethwaite, I., 2005. Multivariable Feedback Control. Analysis and Design. John Wiley and Sons, England. Vanrolleghem, P., Benedetti, L., Meirlaen, J., 2005. Modelling and real-time control of the integrated urban wastewater system. Environ. Model. Softw. 20 (4 Spec. iss.), 427e442.

A generic methodology for the optimisation of sewer systems using stochastic programming and self-optimizing control.

The design of sewer system control is a complex task given the large size of the sewer networks, the transient dynamics of the water flow and the stoc...
2MB Sizes 0 Downloads 9 Views