Using graphs to find economically optimal safety targets for multiple lines of flood defences

lines of flood defences Egidius Johanna Cassianus Dupuits1, Ferdinand Lennaert Machiel Diermanse2, and Matthijs Kok1 1Delft University of Technology, Faculty of Civil Engineering and Geosciences, P.O. Box 5048, 2600 GA Delft, Netherlands 2Deltares Unit Inland Water Systems, department of Flood Risk Management, P.O. Box 177, 2600 MH Delft, Netherlands Correspondence to: E.J.C. Dupuits (e.j.c.dupuits@tudelft.nl)


Introduction
Concerns regarding the safety of people and assets in flood prone areas has led to the construction of flood defence systems all around the world.Some flood prone areas, for example a large part of the Netherlands, face huge potential loss of life and economic value in case heavy flooding occurs.This has led to extensive research regarding estimating the flood risk of flood prone areas.Coupled to this quantification of the flood risk, is the question of 'how safe' a flood prone area should be.An often used approach to help answer this question is a cost-benefit analysis.
Economic optimisation of flood defences, as applied in the Netherlands, is a cost-benefit analysis of the flood risk cost reduction balanced against the investment costs for flood defences.This type of cost-benefit analysis was already developed in the 1950's by Van Dantzig (1956), and is still used and discussed to this day (Eijgenraam, 2006;Kind, 2014).The basic principle behind the economic optimisation of flood defences is finding the minimum of the total costs and is illustrated in Figure 1.The total costs (T C, Eq. 1) are the sum of the annual risk costs ( R) and investment costs ( I) for a given time period.An annual risk cost (R) is in Eq. 2 defined as the annual probability of flooding (P f lood ) times the annual expected loss incurred due to flooding (D f lood ).An alternative term for the annual risk cost is the Annual Expected Damage, or AED.
Generally speaking, a larger investment will lead to a lower AED and this is where the economic optimisation tries to find an optimal solution (i.e. the lowest total cost).

Risk cost Investment Total cost Minimum
Figure 1.Schematic view of an economic cost-benefit analysis for a flood defence.The total costs are the sum of the risk and investment costs, and the optimum can be found at the minimum of the total costs.
Recent publications regarding economically optimal safety targets, for The Netherlands, can be found in Eijgenraam (2006); Brekelmans et al. (2012); Zwaneveld and Verweij (2014b, a).In Eijgenraam (2006) & Eijgenraam et al. (2016), a set of equations was derived which describe the economically optimal safety target for a single homogeneous flood defence system (dike ring).Because of the incorporated influence of time-dependent parameters such as economic growth or climate model parameters, these equations also describe the quantity of (repeated) investments, as well as the optimal time between these investments.These repeated investments are necessary to 'repair' the effect of for example economic growth (i.e. a higher expected losses in case of a flood) or subsidence (i.e. a higher flood probability).A schematic view of the result of such an economic optimisation with time-dependent parameters is shown in Figure 2.This figure shows that the safety level goes down over time, and recurring investments are needed to repair the effects over time of time-dependent parameters such as economic growth or climate model parameters.
The equations described in Eijgenraam et al. (2016) are analytically solvable, and the method resulted in a global minimum of the total costs for relatively simple homogeneous systems.However, dike rings in the Netherlands often consist of mutually different, non-homogeneous sections.This meant the homogeneous case needed to be extended.In Brekelmans et al. (2012), a possible, heuristic solution is given by modelling the problem as a mixed-integer nonlinear programming (MINLP) problem.Zwaneveld and Verweij (2014b) improved on this method by developing a graph based modelling approach to solve the nonhomogeneous case to proven optimality.Eijgenraam (2006); Brekelmans et al. (2012); Zwaneveld and Verweij (2014b)  As the AED changes, the economic optimisation will also be affected.Therefore, it makes sense to explicitly integrate the effect of multiple lines of defence on the AED in the economic optimisation routines.A method to provide a modelling approach to the economic optimisation of a flood defence system with multiple dependent and independent dikes was first presented in 10 Zwaneveld and Verweij (2014a).In Zwaneveld and Verweij (2014a) (in Dutch), a graph based modelling approach is used to obtain economically optimal safety norms and heights for multiple lines of flood defences.Furthermore, they mentioned that the economic optimisation problem can be formulated in the form of a minimal cost flow graph or a shortest path problem.Three approaches to solve economic optimal safety problems for multiple flood defences were identified by Zwaneveld andVerweij (2014a) &Verweij (2014): a heuristic approach based on closed form formulas, a dynamic programming/shortest path approach (as also used in Eijgenraam et al. (2016)) and a branch-and-cut/ILP approach.In Zwaneveld and Verweij (2014a) the branch-and-cut/ILP-approach is preferred and applied.An English description of of the model in Zwaneveld and Verweij (2014a) can be found in (Yüceoglu, 2015, Chapter 5).
However, a consequence of using an ILP approach in Zwaneveld and Verweij (2014a) is that, prior to starting the optimisation routine, all AED estimates for each and every possible combination of flood defences in time need to be computed.Generally speaking, finding AED estimates for a number of these combinations will not be necessary.For example, it is unlikely that is economically optimal to keep all flood defences at their lowest level for the next 300 years.Calculating these AED estimates can be costly, especially if hydrodynamic interactions are included: acquiring a single AED estimate can take hours (De Bruijn et al., 2014) or even days (Courage et al., 2013).In these cases, computationally efficiency will be largely determined by the time it takes to compute AED estimates.
Furthermore, the method of Zwaneveld and Verweij (2014a) is modelled in GAMS and solved using CPLEX.While the method of Zwaneveld and Verweij (2014a) is in principle applicable to an arbitrary number of lines of defence, in practice this means manually extending the model code with new equations that implement any additional lines of defence.While these extensions are trivial for anyone with experience in integer programming and GAMS, we believe that by automating these steps the threshold for using and applying these models can be lowered.
Reducing the number of AED estimates that will be computed can be done based on the principle of 'lazy evaluation', which delays calculations until they are actually required.However, 'lazy evaluation' requires a tight coupling between the AED estimation and the economic optimisation routine.This tight coupling needs to be technically and organisationally possible.
Organisationally, it is possible for projects which are carried out by a single team combining all relevant disciplines; in the remainder of this research we assume that the organisational requirement is fulfilled.Technically, the economic optimisation routine needs to be able to dynamically call the AED estimation function during its optimisation process.However, optimisation routines typically expect a pre-calculated set of data, which means an optimisation routine will need to be modified in order to support 'lazy evaluation'.Therefore, we intend to further investigate the shortest path based approach in order to solve the problem of an economic optimisation for multiple lines of flood defences.The aim of this paper is to find a generic, computationally efficient approach for finding the economically optimal configuration of a flood defence system with an arbitrary number of defence lines.These multiple lines of defence can be dependent on each other (i.e.influence each other's AED).The reliability and performance (in terms of number of AED calculations) of finding economically optimal targets will be tested by comparing the results of the proposed method with a number of benchmark problems.We will accomplish this aim with the following approach: -Computational efficiency will be primarily obtained by minimising the number of (time-consuming) Annual Expected Damage (AED) computations in the algorithm until they are actually required (i.e.'lazy evaluation') -A generically applicable, flexible representation of the problem space will be presented which is able to use an arbitrary number of defences.Specifically, this entails generating a graph in an automated way based on an arbitrary number of lines of defence Section 2 starts with a description of the application and a description of the applied algorithm.Implementation details, focusing on the computational efficiency of the algorithm, are discussed in Section 3, as well as a list of potential future improvements to the algorithm.Next, the proposed approach is applied to some simplified case studies in Section 4, and is followed by a discussion (Section 5) regarding the relevance of the proposed approach.Finally, the results and experiences are concluded in Section 6.
2 An algorithm for flood defence systems with multiple lines of defence

Programmatic representation of the solution space
A common choice to present optimisation problems is to use graph algorithms (Cormen, 2009).Regarding the economic optimisation of flood defences, this choice was also made in Zwaneveld and Verweij (2014b).An example of a graph for a single flood defence is shown in Figure 4.The graph shows the possible investments over time for a single flood defence.In this graph the vertices (dots) are the possible heights the flood defence can have at a certain point in time.In order to go the next point in time, edges are drawn which connect a vertex to all the possible vertices in the next point of time.These points in time are not fixed; the amount and position can be altered to the needs of a particular problem.In practice, these points in time can be related to the (political) decision process of a particular problem: if the relevant flood defences are reviewed and (if necessary) reinforced every five years, it would make sense to have a graph that corresponds to these points in time.
Generally speaking, edges in a graph can be directed or undirected.However, steps backwards in time do not make sense for investment schemes.Therefore, only edges directed forward in time are used.The edge cost (or weight) of an edge is the total cost (AED plus investment cost) of moving between the connected vertices.Furthermore, it is assumed that flood defences will not be intentionally decreased to a lower level, which is why, for example, there are no edges running from h 1 to h 0 .The starting point of the graph is denoted with start in Figure 4 at time t start at a height equal to the current height (h 0 ).
In case of multiple lines of defence, our method takes into account that flood defences can be interdependent and interact with each other hydrodynamically.This means that the AED of the system of defences can potentially be influenced by each defence, which also means that each combination of flood defence levels has to be considered relevant.For a graph with multiple interdependent flood defences, these combinations replace the height of a single flood defence on the y-axis in Figure 4.These combinations of heights for multiple flood defences can be obtained by computing the Cartesian product of the flood defence levels of all the involved defences.For n flood defences, the Cartesian product equation for determining the combinations is shown in Eq. 3: where X i is a vector containing all the flood defence levels of flood defence i, and x i is a realisation of vector X i (i.e. a flood defence level for flood defence i).If all vectors X i are of the same length y, the total number of combinations will be y n .
The number of relevant system combinations reduces significantly if each flood defence can be optimised independently of the other flood defences in the system.The assumption of independence can be made if none of the flood defences in a system have a (significant) influence on the AED estimates of the other flood defences.The total number of system configurations under the independence assumption is n • y, as each flood defence could then be optimised separately (e.g. using a graph per flood defence similar to the graph shown in Figure 4.If only some defences are independent from other flood defences in the system, this is considered a special case of our approach.In that case, our method of using the Cartesian product is still valid and applicable, although it will result in larger than necessary graph.In case of some independent elements, a possible approach to reduce the size of the graph is discussed in Section 3.4.are only drawn to vertices containing sets of heights equal or greater than the set of heights in the vertex at the origin of the edge.However, because Figure 5 has two defences instead of one, the outgoing edges are slightly different when compared to Figure 4.For example, the height combination h A0 , h B1 is never connected to h A1 , h B0 (since that would correspond to a reduction in height for defence B).
. Graph with vertices (dots) and edges (arrows) for two defences (A and B).Each defence has two possible heights.

Implementation of a graph algorithm
In general terms, a graph algorithm will iterate over vertices in a graph in an effort to find the path with the lowest costs between a given start and end vertex.However, in the graphs of Figure 4 and Figure 5 t end contains a number of possible end points, which means that the algorithm will need to find as many optimal paths as there are end points in the graph.In order to only have to run the algorithm once, a stop vertex is added; the graph of Figure 5 with an additional stop vertex is shown in Figure 6.

5
The edges running towards this stop vertex are all given a weight of zero.Now, the algorithm only has to find a single optimal path between t start and t stop .Why this is an efficient contribution is illustrated in Section 2.4.The graph as shown in Figure 6 is a graph with directed non-negative edges.For this kind of graph, a number of algorithms can be used to find the shortest (optimal) path in a graph, for example: the Dijkstra algorithm (Dijkstra, 1959), the A* algorithm (Hart et al., 1968), and the Uniform Cost Search (e.g.(Verwer et al., 1989)).All three can be considered to be part of the family 10 of best-first search algorithms, where both the Dijkstra and the Uniform Cost Search (UCS) algorithms can be seen as a special case of the A* algorithm.
Typically, the best-first search algorithms are implemented with a min-priority queue.A min-priority queue holds a sorted list of vertices, where the sorting is based on the cost of reaching that vertex from the start vertex; the vertex with the lowest cost is at the top of the queue.This list of vertices in the priority queue constitutes of, depending on the implementation, either all vertices in the graph (Dijkstra as implemented in Cormen ( 2009)), or only the vertices already visited by the graph algorithm (UCS).A comparison between the two algorithms can be found in Felner (2011), where the priority queue as implemented by UCS was found to be faster and using less memory.We consider this a relevant advantage, as the number of vertices can be large when using the Cartesian product of flood defence levels (Section 2.1).For this reason, we chose to implement the UCS algorithm.
In Eijgenraam et al. (2016), a dynamic programming approach was used, which is related to the shortest-path algorithms discussed thus far.However, the Dijkstra algorithm (and by extension the UCS and A* algorithms) are seen in Cormen (2009) as a part of the greedy shortest-path algorithms family, which in Cormen (2009) is clearly defined as a different type of algorithm than dynamic programming.Greedy algorithms are typically much faster than dynamic programming approaches, at the expense of not always finding the optimal solution (because less possible solutions are considered).The optimality condition is further discussed in Section 2.4.Nevertheless, because less possible solutions are considered in a greedy algorithm, this can also lead to a part of the graph never being visited by a greedy algorithm.Combined with 'lazy evaluation', this can lead to a significant reduction in the number of AED calculations which are actually executed; see also Section 3.3.
Applying the UCS algorithm to a graph such as shown in Figure 6 begins with creating a priority queue which only contains the start vertex.After this initialization, the iteration process is started.Each iteration starts with taking out the vertex with the lowest cost known thus far from the priority queue (which is the top entry in the queue).Taking out means the optimal route (lowest cost) from the start vertex to this vertex now known.The vertex that has just been taken out of the priority queue is then queried in the graph to find all the connecting vertices in the next time step.Each connecting vertex is added to the priority queue if the vertex is not already in the queue.If the vertex already exists in the queue, the weight is only updated if the newly proposed cost is lower than the known cost so-far.Iteration continues until at the start of an iteration the stop vertex is the top entry in the priority queue.An actual example of the application of this algorithm will be elaborated in Section 2.3.

Example application of the algorithm in an economic optimisation
This section shows a simple example of an economic optimisation for a single flood defence.While this example uses a single flood defence for simplicity, the same principles apply for multiple flood defences.Regarding the investment costs and AED estimates, if a vertex at t 1 is connected to another vertex with a larger height at t 2 , it is assumed that the actual heightening occurs at t 1 .This leads to a slightly different graph than the conceptual implementation shown in Section 2.1 & 2.2, and is emphasized by drawing the edges of the figures in this example (i.e.Figures 7, 8 & 9) in a way which is visually more consistent with the timing of the investment decision.
The result of the first two iterations is shown in Figure 7, where the start vertex is labelled with the number 1.In this example, the start vertex is associated with a height of 4.25 meter and starts at t = 0, identical to vertex 2. Because the path to vertex 2 is the only possible path, vertex 2 is the only addition to the priority queue.In the next iteration, vertex 2 is taken out of the priority queue as it is the vertex with the lowest total cost.The total cost to reach vertex 2 is 0, because there was no heightening (height remains at 4.25 meter) and no time expired (t start = 0); the AED is zero because time needs to expire for risk to occur.
From vertex 2, the number of possible next steps and associated total costs are computed and added to the priority queue, as illustrated in Figure 7.Note that the total costs to reach for example vertex 22 consists of the total cost from t start to t = 100, not just the cost from t = 50 to t = 100.
The algorithm will continue for a while, until the situation of Figure 8 is reached where vertex 24 is taken out of the priority queue.The new found total costs for vertex 29, 30 and 31 are not lower than the total cost for vertex 25, which means the algorithm takes a step back and continues from vertex 25.From vertex 25, vertex 30 and 31 are re-evaluated in Figure 9, where only vertex 30 results in lower costs than the existing options.This means that only vertex 30 is updated with the new, lower, total cost in the priority queue.Additionally, if hypothetically vertex 31 is the vertex with the lowest cost, the optimal path would revert back to using vertex 24 instead of vertex 25 (because the path from vertex 25 to 31 has higher costs than the path from vertex 24 to 31). which have total costs above their respective vertices, and the first three entries are shown in the column on the right.Note that because the choice was made to connect the start vertex to vertex 2, vertices 3 -6 will not be visited.

Global optimal solution
The UCS algorithm finds the shortest path in a graph, see for example Felner (2011) for a recent elaboration regarding the 'correctness' of the UCS algorithm or Gelperin (1977) for a proof regarding A* (UCS can be considered a special case of A*).
What remains is whether the additional stop vertex of Section 2.2 leads to a potential heuristic solution or still to the optimal path.However, assuming that the optimal path towards the stop vertex is found, whichever vertex at t end is part of that optimal path has to be the optimal choice.Otherwise, the path towards the end vertex is not optimal, which contradicts the earlier mentioned proofs.In order to further test the performance of the proposed method, Section 4 will compare numeric results from our proposed method to other approaches.These approaches are known to give global optimal results.

Overview of the approach
A general overview of the approach discussed in the previous sections is shown in Figure 10.The method is composed of four steps: input, pre-processing, processing and post-processing.Of these steps, user interaction is only required at the input step.
The rest of the steps run automatically.Specifically, the user needs to supply vectors of flood defence levels per flood defence, a time vector and a function which can calculate the cost of an edge in the graph.In the following steps, the graph is created (pre-process), the optimal path is found (process), and the optimal path is shown (post-process).

Efficiency improvements
The economic optimisation of multiple lines of defence, implemented in a graph using Section 2, can potentially lead to large numbers of vertices and even larger numbers of edges.For example, for eight lines of defences with six possible heights the number of vertices per time step is approximately 1.68 million (6 8 ), while the number of edges per time step is even larger at approximately 35 billion.For large problems such as these, storing all the possible vertices and edges would lead to huge data structures and a huge number of AED calculations.This requires both an efficient implementation of the graph, and an efficient evaluation of AED calculations (i.e. as few as possible).An efficient graph implementation is discussed in Sections 3.1 & 3.2, while the efficient evaluation of AED calculations is discussed in Section 3.3.Potential further efficiency improvements are discussed in Section 3.4.

Repetitiveness in lists of vertices
Even though the graphs of Section 2.1 can be classified as sparse graphs (number of edges is much smaller than the number of vertices squared, Cormen (2009)), the number of edges is still much larger than the number of vertices.Therefore, we first focused on data structures related to the edges of a graph.For sparse graphs, these are the adjacency lists: a group of vertices connected via edges stemming from a source vertex in a previous time step.In these adjacency lists, repetitiveness can be found with respect to two aspects.
The first repetitive aspect is the similarity of adjacency lists for the same combination of flood defence levels at different time steps (except for the adjacency lists at t end ).In Figure 7, vertices with the same combination of flood defence levels at different time steps are for example vertices 2, 7 and 12.The adjacency lists for these vertices are shown in Figure 8, where it is apparent that the adjacency list for the next time step can be found by adding an offset to the elements of the adjacency list of the current time step.For example, the adjacency list of vertex 2 can be turned into the adjacency list of vertex 7 by adding the total number of combinations in each time step (which is five in Figure 11) The second repetitive aspect is for adjacency lists between vertices in the same time step.Because the lowest vertex in each time step (e.g.vertices 2, 7, 12, 17, 22 and 27 in Figure 7) has outgoing edges running to each and every vertex in the next time step, higher vertices (e.g. in Figure 7, vertex 8 is 'higher' than vertex 7) contain a subset of the adjacency list of the lowest vertex.In other words, outgoing edge lists in a single time step can be generated dynamically by shrinking the adjacency list of the lowest vertex in a time step.This is shown in Figure 12.
The combination of these two repetitive characteristics results in that only a single adjacency list needs to be stored in memory (i.e. the adjacency list of the lowest vertex in the first time step).This single adjacency list can be adapted to most vertices in the graph by means of offsetting and shrinking the stored adjacency list.Notable exceptions are the adjacency lists for the vertices at t end , but the adjacency lists for these vertices are already known and only contain the stop vertex.

Conditionally removing edge connections
Besides reducing the size of the data structures associated with a graph, the adjacency list associated with a vertex can also be reduced under certain conditions.Typically, the time between improvements in flood defences is large (in the order of 50 years), due to either high (fixed and variable) costs associated with investments in flood defences, or long planning periods (Zwaneveld and Verweij, 2014b).Therefore, if one or multiple flood defences have been strengthened recently, the adjacency 5 list can be reduced to only contain vertices that keep the recently strengthened flood defence(s) at the current level(s).However, this so called 'waiting time' before new investments are considered has to be chosen with care, because the waiting time should not influence the optimal time between investments.Nevertheless, a correctly chosen waiting time can greatly improve the run time of the algorithm, because of the significant reduction in number of edges that need to be evaluated.This reduction is shown in Figure 13, where the total number of visited vertices is plotted as a function of the 'waiting time'; the underlying 10 problem that is solved by the algorithm is the same problem as shown in Section 4.3.

Reducing the number of AED calculations
In the overview of Figure 10 it is implied that the AED calculations belonging to an edge are only carried out when that edge is visited by the graph algorithm.Provided that a graph algorithm does not visit all vertices, delaying AED calculations belonging to an edge until that edge is visited leads to less AED calculations than the total number of possible AED calculations in a particular graph.In contrast, if AED (or more generally, cost) calculations are done before a graph algorithm is initialised, all possible AED calculations need to be calculated beforehand.
As an example, Figure 14 shows the number of times each vertex is visited is in the example of Section 2.2.The majority of the vertices in Figure 14 get visited once, but a significant proportion is never visited by the graph algorithm; these vertices have a zero above their indices.A small proportion of the vertices, specifically vertices 30 and 31, are visited twice; the reason for this re-visiting can be seen in Figure 9.To avoid completely re-doing cost calculations upon a revisit, parts of a calculation can be cached in order to reduce the computational penalty incurred by revisiting vertices.The total number of possible AED calculations is the number of years multiplied with the number of options on the y-axis; for Figure 14 this leads to a total number of AED calculations of 1505 (or 301 • 5).Because a number of vertices do not get visited by the algorithm, the number of actual executed AED calculations goes down to 1000, or approximately 66% of all possible AED calculations.
Furthermore, the number of AED calculations can be further reduced by using the 'waiting time' of Section 3.2.In Section 3.2, it was found that a minimum waiting time between investments will lead to less edges being evaluated by the algorithm.This also implies that less AED calculations will be executed.Using the same example as in Figure 13, the reduction in the percentage of actual executed AED calculations is given as a function of the waiting time in Figure 15.Between using a waiting of 0 years (i.e.no minimum waiting time at all) and a waiting time of 50 years the number of AED calculations goes down from approximately 60% to 40% for the example of Section 4.3.

Potential improvements and special cases
Further improvements can be made both to the graph implementation and to the implementation of the algorithm.The algorithm was implemented as a single process; a performance improvement might be found by utilizing parallel programming.The first place where parallel programming could be beneficial is the loop over an adjacency list.This is because the potentially expensive AED calculations are done as part of determining an edge weight.Therefore, parallelising the loop over an adjacency list over multiple computational nodes can lead to significant performance improvements.
Furthermore, regarding the graph implementation, a special case is a flood defence system which has independent flood defences.Section 2.1 uses the Cartesian product of flood defence options, which has the underlying notion that all flood defence lines are interdependent.If some flood defences are independent (i.e. the defences protect different, independent areas), this leads to an inefficient graph.The independency of flood defences can be used in an adapted graph representation in order to get an efficient graph.While we did not implement this, a way to solve this for the system in Figure 16 is shown in Figure 17, which uses 'subgraphs' to reduce the number of combinations.
These subgraphs are small graphs which only contain the number of strengthening options for a single defence for a single time period (e.g. in Figure 17, from t i−1 to t i ).Additionally, the subgraphs take into account what the level is of the influential defences (e.g. in Figure 17, the front defence B is the only influential defence for the rear defences).The use of subgraphs leads to a smaller number of combinations, as the Cartesian product would have resulted in a total number of combinations (or vertices per time step) of 5120 (5 • 4 5 ).With subgraphs, the number of combinations is reduced to 100 (5 • 5 • 4).
4 Results for simplified flood defence systems In order to test the performance of the proposed algorithm versus some existing approaches, three cases are investigated.For simplicity, these three cases will be based upon a common set of investment and AED relations, as well as a common set of input values.The values and symbols used in this section are largely copied from (Eijgenraam, 2006, page 34) and reproduced in Table 1, with only minimal changes.These AED and investment cost relations consist of simple formulations which were Table 1 in Eijgenraam (2006).The sum of the investment cost and AED is the total cost, which needs to be minimised in order to get economically optimal safety targets: R (t) = P 0 e −α(Ht−H0−ηt) V 0 e γt e −δt (5) where sign(u t ) is used to prevent fixed costs in case there is no heightening u t .This sign(u t ) function returns zero if the heightening u t is equal to zero, and returns one when the heightening u t is larger than zero.

Single flood defence
For a single flood defence, with the values of Table 1, an analytical solution can be found in (Eijgenraam, 2006, page 35).This solution consists out of an initial dike height increase coupled with a periodical, constant dike increase over an infinite time 5 horizon.The initial increase was found to be 236 centimetres, with a periodical increase of 129 centimetres every 73 years.
Because the approach introduced in this paper is a numerical approach, a finite time period has to be used instead of an infinite time horizon.Similar to Zwaneveld and Verweij (2014b), we choose to use a time period of 300 years with, for this application, steps of one year.The possible heights were discretized using a range starting from 425 to 1225 cm, with steps of one centimetre.Note that these step sizes (and dimensions) were deliberately chosen to be on par with the accuracy level of in this problem is 241,101 (or 801 • 301).Of these, 137,971 were actually executed by the UCS algorithm, which corresponds to using only 57% of all possible AED calculations.If a 'waiting time' of 50 years is used, the solution is unaffected but the percentage of executed AED calculations goes down to 43%.
A comparison of the results found using the algorithm and the analytical solution is shown in Figure 18.The algorithm found an initial increase of 235 cm, with three additional increases in height at 73 years apart.These three were found to be 129 cm, 130 cm, and 132 cm.The last increase is different from the analytical solution, and can be attributed to being close to the end of the time horizon.A finite time horizon implies that there is no AED beyond the time horizon, which explains why there is no investment found by the algorithm in year 292.To compensate for the lack of an investment in year 292, the investment in year 219 is slightly larger.This explanation is supported by results with a time horizon of 400 years, where the heightening in year 219 changes to an expected 129 centimetres.These deviations near to the time horizon underline that if a certain point in time is considered relevant, the used time horizon should stretch significantly beyond that point in time.However, this is a general problem with all numerical methods, because of the required finite time horizon, and not a specific issue related to the approach proposed in this paper.Furthermore, in practice this problem can be circumvented by setting the time horizon used in the algorithm to sufficiently exceed the practically required time horizon.

Two independent lines of defence
In the next example two defences are investigated using the graph algorithm, both with the same characteristics as the single flood defence in the previous section.However, the step size for the heights is increased to 20 cm in order to test the response of the algorithm to larger step sizes.Expected is that, despite the less detailed step size, the investment scheme for both defences should be identical to each other and close to the analytical solution provided in the previous section.
Indeed, the results of the algorithm, illustrated in Figure 19, show that both defences are initially increased with 240 cm, while in both year 75 and 143 the defences are increased with 120 cm, and finally in year 212 with 140 cm.Clearly, the larger step size in height leads to larger differences when compared to the analytical solution.Nevertheless, any overshoot/undershoot of the height is 'repaired' in the duration between investments, keeping the solution of the optimal path stable and close to the analytical solution.Furthermore, the total number of possible AED calculations in this problem is 24,682 (or 2 • 41 • 301).
Of these, 14,510 were actually executed by the UCS algorithm, which corresponds to using only 59% of all possible AED calculations.If a 'waiting time' of 50 years is used, the solution is unaffected but the percentage of executed AED calculations goes down to 48%.Analytical Algorithm Figure 19.The two (independent) lines of defence have an identical solution with the approach proposed in this paper, and are (even with the usage of larger step sizes) good approximations of the known analytical solution.

Two dependent lines of defence 5
The final case is similar to the case with two independent lines of defence, however the second defence is now dependent on the performance of the first defence.This dependency is illustrated in Figure 16, and is a simplified version of the case discussed in Dupuits et al. (2016).where P i is a generic formulation used for the failure probabilities P 1 , P 2|1 and P 2|1 .The probabilities and are the failure probabilities of the second defence, dependent on the failure (P 2|1 ) or non-failure (P 2|1 ) of the first defence, where the failure probability of the first defence is denoted by P 1 .Similarly, the investment equation in Eq. 6 is expanded to include different should result in two optimal configurations (one with a closure dam, the other with a storm surge barrier), which can then be compared using some metric, for example their benefit-cost ratios.

River Retention area
Figure 22.A retention area can also be optimised using the approach proposed in this paper.In this example, the surface area of the retention basin is used instead of the height of a flood defence.

Conclusions
This paper presented a generic, computationally efficient approach for finding the economically optimal configuration of a flood defence system with an arbitrary number of defence lines.The computational efficiency was achieved by delaying AED (flood risk cost) calculations until they are actually needed in an optimisation routine (i.e.'lazy evaluation'), which leads to a reduction in the number of AED calculations that need to be done.In the few examples shown in this research, the reduction in number of AED calculations was at least 40%.This is a significant and relevant reduction, as the AED calculations relevant for this approach often have a high computational cost.This is especially the case when multiple flood defences interact with each other hydrodynamically in a larger flood defence system.
The proposed approach is implemented using a best-first graph algorithm, which itself is a simple algorithm.Furthermore, the approach is flexible towards the number and type of flood defences because the graph representation shown in this paper can trivially accommodate an arbitrary number of interdependent flood defences.The proposed approach utilizes the repetitive properties of the graphs in order to efficiently store the representation of the graph.In case independent flood defences are present in a system, the proposed approach of generating a graph can be adapted to a more efficient method which makes use of the attractive properties of independence.To that end, an improvement in the generation of the graphs has been proposed.
Assuming that the graph and combinations of flood defences are portrayed correctly, the best-first graph algorithm has been proven in literature to return the shortest (or optimal) path in a graph.To corroborate this for our implementation and intended application, the method was tested on a number of benchmark problems with known solutions.The tests show that indeed the optimal path is found with the approach proposed in this paper, which shows that the implementation was done correctly.

Figure 2 .
Figure 2. Schematic view of an economic cost-benefit analysis for a flood defence, with time-dependent parameters.Because of these time-dependent parameters (e.g. economic growth or subsidence), recurring investments in safety are needed.

Figure 4 .
Figure 4. Graph where the vertices (dots) at each time step are connected via edges (arrows) to the next time step.

Figure 5
Figure 5 shows an example of the Cartesian product for two flood defences where each flood defence has two possible heights.The graph in Figure 5 resembles the graph in Figure 4 for a single flood defence.Similar to Figure 4, edges in Figure 5

Figure 6 .
Figure 6.The graph of Figure 5 with an additional stop vertex.

Figure 7 .
Figure 7.The first two iterations of the graph algorithm with a min-priority queue.The vertices available in the priority queue are those

Figure 8 .Figure 9 .
Figure8.After six iterations with the graph algorithm, vertices 29, 30 and 31 are added to the priority queue.However, in this case the algorithm makes a step back in time, because vertex 25 is the item with the lowest total cost in the priority queue.

IFFigure 10 .
Figure10.Overview of the approach using a graph and graph algorithm.In our approach, the graph algorithm is the UCS algorithm.The input column is the only part what the user should provide, the other steps run automatically.

Figure 11 .
Figure 11.The adjacency lists for vertices 7 and 12 of Figure 7 can be obtained by adding an offset to the adjacency list of vertex 2.

Figure 12 .
Figure 12.The adjacency lists for vertices 3 and 4 of Figure 7 are reduced sets of the adjacency list for vertex 2.

Figure 13 .
Figure 13.Total number of visited vertices as a function of the waiting time for the example of Section 4.3.

Figure 14 .Figure 15 .
Figure 14.Number of times each vertex is visited by the algorithm for the example in Section 2.2.The dotted area emphasizes that a part of the graph is never visited, while vertex 30 and 31 get visited twice.

Figure 18 .
Figure18.The investment scheme found using the algorithm is almost identical to the analytical solution.

Figure 20 .
Figure 20.A coastal system with two lines of defence.This figure is an adaption from an illustration found in Dupuits et al. (2016).

Table 1 .
Eijgenraam (2006)ues taken fromEijgenraam (2006)for the AED and investment equations in this section.NLG refers to the currency used in the Netherlands prior to the euro.