Minimaliseren van kosten van dijkverhoging (Minimization of costs of dike raising)
ZOË LAGERWEIJ Verslag ten behoeve van het Delft Institute for Applied Mathematics als onderdeel ter verkrijging van de graad van Bachelor of Science in Technische Wiskunde Delft, Nederland Juni 2008
Copyright © 2008 door Zoë Lagerweij. Alle rechten voorbehouden.
BSc verslag Technische Wiskunde
Minimaliseren van kosten van dijkverhoging (Minimization of costs of dike raising)
ZOË LAGERWEIJ
Technische Universiteit Delft
Begeleider Prof.dr.ir. C. Roos Overige commissieleden Prof.dr. F.M. Dekking Drs. A. Verweij / Dr. J.G. Spandaw Delft, Nederland Juni 2008
Samenvatting In 1953 vond in Nederland de watersnoodramp plaats. Als gevolg daarvan werd de Deltacommissie ingesteld. Haar doel was om de overheid te adviseren welke maatregelen getroffen moesten worden zodat een tweede ramp voorkomen werd. De adviezen van de Deltacommissie hebben bijvoorbeeld geresulteerd in de afsluiting van de Nieuwe IJssel, de Oosterschelde, de Grevelingen en het Haringvliet. In Nederland zijn echter nog veel meer beschermingsmaatregelen. Er zijn 57 dijkringen. Een dijkring is een in de wet vastgelegd binnengebied dat voor water beschermd wordt door onder andere dijken, sluizen en gemalen. Ook de dijkringen zijn onderhavig aan verzakking en waterpeilstijging. Hierdoor is het belangrijk dat vaak onderzocht wordt of de dijken en de andere waterkeringen aangepast moeten worden. De verhoging van dijken staat centraal in dit verslag. Wanneer en hoeveel moet je een dijk ophogen zodat de kosten minimaal zijn. Om dit probleem op te lossen zijn verschillende formules nodig, zoals de kans op overstroming, verwachte schade bij overstroming en de investeringskosten voor het ophogen van de dijken. Aan de hand van deze formules kan het probleem gediscretiseerd worden tot een kortste pad probleem. Dit probleem wordt opgelost aan de hand van een algoritme dat in Matlab is ge¨ımplementeerd. In de minimale oplossing zijn de ophogingen nagenoeg periodiek. Dit is niet doordat de maximale toelaatbare overstromingskans wordt overschreden, maar doordat de verwachte schade, die als kostenpost wordt genomen, kleiner wordt als de kans op overstroming ook klein is. In de formules zijn verschillende waarden die afhangen van de economische situatie in het land. Dit zijn de rentevoet en de groei van de welvaart. Verschillende waarden hiervoor zijn ingevoerd in het model. Als de rente hoger wordt door hoge inflatie, worden de kosten van de ophogingen lager. Als de welvaart groeit, worden de kosten voor dijkophogingen hoger. Dit komt doordat bij overstroming meer verlies wordt geleden in een welvarender gebied. In New Orleans ging het in 2005 ook mis. Van deze ramp kan men ook in Nederland leren: New Orleans is vergelijkbaar met Nederland doordat het ook onder zeeniveau ligt en beschermd wordt door dijken. Uit deze ramp is gebleken dat men meer voorbereid moet zijn op een ramp en alle kosten van zo’n ramp ingecalculeerd moeten worden.
Inhoudsopgave 1 Inleiding
6
2 Verschillen tussen Eijgenraam en Van Dantzig
7
3 Beschrijving van het probleem 3.1 Dijkring 10 . . . . . . . . . . . . . . 3.2 Relevante formules . . . . . . . . . 3.2.1 Overstromingskans . . . . . 3.2.2 Schade bij overstroming . . . 3.2.3 Verwachte schade . . . . . . 3.2.4 Kostenformules . . . . . . . 3.3 Actuele waarde en groei van welvaart 3.3.1 Actuele waarde . . . . . . . 3.3.2 Groeipercentage van welvaart
. . . . . . . . .
8 8 9 9 10 10 11 12 12 13
4 Achtergrondinformatie 4.1 Waar het wel mis ging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 New Orleans (2005) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Nederland (1953) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15 15 15 16
5 Werkwijze 5.1 Discretiseren en dynamisch programmeren . . . . . . . . . . . . . . . . . . . . . 5.2 Programmeren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17 17 18
6 Resultaten 6.1 Uitkomsten voor Dijkring 10 . . . . . . . 6.2 Uitkomsten voor andere dijkringen . . . . 6.3 Uitkomsten voor verschillende constanten 6.3.1 Uitkomsten voor verschillende δ . 6.3.2 Uitkomsten voor verschillende γ . 6.3.3 Uitkomsten voor verschillende γ en
19 19 22 25 25 25 26
. . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . δ
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . . . . .
. . . . . .
. . . . . .
7 Conclusie
27
A MatLab-modellen A.1 Gebruikte functies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1.1 Overstromingskans . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31 31 31
2
A.1.2 Schade . . . . . . . . . . . . . A.1.3 Verwachte schade op tijdstip t A.1.4 Kosten . . . . . . . . . . . . . A.2 Algoritme . . . . . . . . . . . . . . . A.3 Toepassingen . . . . . . . . . . . . . A.3.1 Algoritme . . . . . . . . . . . A.3.2 Rente . . . . . . . . . . . . . A.3.3 Groei van welvaart . . . . . . A.3.4 Vari¨eren δ en γ . . . . . . . .
3
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
31 31 31 32 35 35 38 38 39
Lijst van gebruikte symbolen Ht = dijkhoogte op tijdstip t
(1)
Pt = overstromingskans op tijdstip t
(2)
Vt = schade bij overstroming op tijdstip t (miljoen euro)
(3)
St = Pt Vt = verwachte schade op tijdstip t
(4)
α = parameter exponenti¨ele verdeling extreme waterpeilen (1/cm)
(5)
η = structurele verhoging van waterpeil (cm/jaar)
(6)
γ = groeipercentage van de welvaart in de dijkring (per jaar)
(7)
ζ = groei van schade per cm dijkverhoging (1/cm)
(8)
δ1 = discount factor welvaartskosten (1/jaar)
(9)
δ2 = discount factor investeringskosten (1/jaar)
(10)
T
(11)
= aantal jaren waarover het model loopt
Pmax = maximaal toegestane overstromingskans β = αη + γ β
0
= β − δ1 = αη + γ − δ1
θ = α−ζ
(12) (13) (14) (15)
Voorwoord Dit is het verslag van mijn bachelorproject. Het gaat over kosten van het ophogen van dijken in Nederland. In eerste instantie wist ik alleen in welke richting waarin ik dit project wilde doen, namelijk optimalisering. Professor Cees Roos was zelf al bezig met een onderzoek over de dijken in Nederland. Dit kon ik als onderwerp nemen. Dit leek mij erg interessant door de actualiteit. In New Orleans is het net misgegaan en door de waterpeilstijging loopt Nederland steeds meer risico. Maar ook door de grootsheid: de dijken zijn voor iedereen in Nederland even belangrijk en iedereen heeft ermee te maken. Het specifieke probleem: Minimaliseren van kosten van dijkverhogingen, is een wiskundig probleem, maar heeft een economische inval. De economie van een land heeft namelijk invloed op de schade bij overstroming. Aangezien ik in het verleden bijna geen economie heb gehad, was het een goede mogelijkheid om daar meer over te leren. Alhoewel het economische deel niet diep gaat heb ik nu toch een beter idee over hoe de economie werkt. Mede door dit onderwerp heb ik veel plezier gehad aan dit project, vooral het maken van het model, wat toch altijd een grote puzzel is. Maar als je eindelijk een oplossing hebt, geeft dat een enorme voldoening. En dat is wat dit project mij heeft gegeven. Zo¨e Lagerweij Delft, 3 juni 2008
Hoofdstuk 1
Inleiding In 1956 schreef Professor D. Van Dantzig naar aanleiding van de watersnoodramp in 1953 een rapport ([1]) over de kosten-batenanalyse voor de optimale hoogte van dijken. De formule die hij daarvoor opstelde wordt nog steeds gebruikt. Echter, volgens C.J.J. Eijgenraam(CPB) is zijn formule onjuist. In zijn rapport ([2]) geeft hij een (volgens hem) betere formule voor de kostenbatenanalyse. In hoofdstuk 2 staan deze verschillen uitgelegd. Professor C. Roos heeft samen met Professor D. Den Hertog(Universiteit van Tilburg) aan de hand van Eijgenraams formules verschillende modellen gemaakt om de tijdstippen van dijkverhogingen te berekenen wanneer de kosten minimaal zijn. Dit hebben zij met vier verschillende modellen gedaan. E´en van deze modellen was gebaseerd op dynamisch programmeren. Dynamisch programmeren herleid het probleem tot een een kortste padprobleem. In dit geval wordt eindtijdstip T genomen en kan er een aantal keer op verschillende tijdstippen worden verFiguur 1.1: Grafische weergave van kortste hoogd om een veilige hoogte te bereiken. Daarna pad probleem wordt gekeken via welke weg het goedkoopst opgehoogd kan worden (zie figuur 1.1). Om de kosten te minimaliseren is het nodig om bepaalde variabelen vast te leggen. Soms zijn dit gegevens, maar vaak is het een moeilijke afweging. In dit verslag zal ik de invloed van twee variabelen gaan onderzoeken het groeipercentage van de welvaart en de actuele waarde(Engels: discount rate), γ en δ. Dit zal ik doen door in Matlab een model te schrijven met data die ook gebruikt zijn in [3]. Vervolgens zal ik de twee waarden gaan vari¨eren. Voor het model dat ik schrijf zal ik ´e´en dataset gebruiken, namelijk de data van dijkring nummer 10, ofwel Mastenbroek, gelegen in Overijssel. Alleen voor deze dataset zal ik δ en γ vari¨eren. Om te laten zien hoe de verschillen tussen de dijkringen in het model uitkomen heb ik de data van nog drie andere dijkringen in het model gezet. In dit verslag zal ik allereerst uitleggen wat de situatie is met betrekking tot de dijken en achtergrond informatie geven. Daarna zal ik mijn oorspronkelijke idee¨en over mijn aanpak uiteenzetten en vervolgens het verloop. Daarna zal ik de resultaten uitleggen en een conclusie trekken. 6
Hoofdstuk 2
Verschillen tussen Eijgenraam en Van Dantzig In dit hoofdstuk zal ik de verschillen tussen de modellen van Eijgenraam en Van Dantzig bespreken. Tussen de twee modellen zitten grote en kleine verschillen. Ik zal de verschillen uitleggen die van invloed zijn op mijn model. Het grootste verschil is dat Van Dantzig alleen antwoord gaf op de vraag hoe duur de dijkophogingen zouden zijn. Eijgenraam geeft ook antwoord op de vraag wanneer een dijk moet worden opgehoogd ([2]). Verder gaat Van Dantzig er in zijn model van uit dat de noodzakelijke verhoging van een dijk gelijk is aan de hoogte die de dijk verloren heeft door waterspiegelstijging en verzakking (ηt). Eijgenraam neemt ook mee dat de welvaart in een gebied kan groeien, waardoor de schade bij overstroming groter wordt. Daardoor wordt de dijk iedere keer meer opgehoogd dan het waterpeil gestegen is en de dijk verzakt is. De investeringskosten zijn daardoor iedere keer gelijk op de term e−δt na. Een ander verschil is dat Van Dantzig de vaste kosten van de dijkverhoging alleen neemt bij de eerste investering. Ook heeft hij in de investeringskosten de groei van welvaart meegenomen en van de rentevoet, δ, afgetrokken. Van Dantzig heeft de grotere schade bij overstroming als de dijk hoger wordt niet opgenomen in zijn model. Dit geeft geen groot verschil tussen de modellen en geldt volgens Eijgenraam alleen bij rivieren ([2]).
7
Hoofdstuk 3
Beschrijving van het probleem In dit hoofdstuk zal ik uitleggen wat een dijkring is en wat achtergrondinformatie geven over de in dit verslag meest gebruikte dijkring. Ook zal ik de meest gebruikte formules uitleggen en uitleggen wat γ en δ zijn.
3.1
Dijkring 10
Een dijkring of dijkringgebied is een gebied dat door een stelsel van waterkeringen en hoge gronden beschermd wordt tegen buitenwater. De waterkeringen kunnen bestaan uit dijken, duinen en constructies voor het waterbeheer zoals sluizen en gemalen. Voorbeelden van hoge gronden zijn natuurlijk hooggelegen gebieden zoals de Veluwe en de Utrechtse heuvelrug. In de wet is bepaald hoe groot de kans op overstroming in iedere dijkring mag zijn. Deze voorschriften zijn bepaald door de Deltacommissie. De dijkringen zijn ook in deze wet vastgelegd ([5]). Dijkring nummer 10, de dijkring die centraal staat in dit verslag, is Mastenbroek. De polder Mastenbroek ligt in Overijssel tussen de IJssel en het Zwarte Water en tussen de steden Hasselt, Kampen en Zwolle (Zie figuur 3.1) [6]. De polder werd vanaf 1364 ingepolderd en ontgonnen ([7]). De data die voor deze polder gelden zijn: H0 T η α P0 Pmax ω
= = = = = = =
0 100 0.320 0.033027
V0 C b λ ζ γ δ1 = δ2
1 2270 P0 Pmax 1 α log( P0 )
8
= = = = = = =
1564.9 16.6939 0.6258 0.0014 0.003774 0.02 0.04
Figuur 3.1: Een kaart van Mastenbroek uit 1750 ([6])
3.2
Relevante formules
Om het model te kunnen maken zijn bepaalde formules nodig. Hier zal ik uitleggen wat die formules zijn en waarom ze nodig zijn. Deze formules staan ook in [2] en [3].
3.2.1
Overstromingskans
De kans op overstroming hangt af van verschillende factoren, om te beginnen de hoogte van een dijk. De overstromingkans van een dijk is p(h) = 1 − F (h), waarbij F (h) de cumulatieve verdeling van hoogtij is en p(h) de kans dat de hoogte h van hoogtij overschreden wordt. Deze kans is ceαh , waarbij α de parameter voor de exponenti¨ele verdeling voor extreme waterniveaus is. In 1939 benaderde P.J. Wemelsfelder deze formule op basis van data van de overschrijdingsfrequenties in Hoek van Holland ([1]). In de formule nemen we p0 = ceαH0 : de kans op overstroming bij de beginhoogte. Hierdoor wordt de formule: p(h) = p0 e−α(h−H0 ) . (3.1) 9
Deze formule hangt nog steeds alleen van de hoogte van een dijk af. De hoogte van een dijk hangt echter ook van de tijd af, rekening houdend met de stijging van het waterpeil. De parameter van de stijging van het waterpeil wordt aangegeven met η in centimeter per jaar. De formule voor de kans wordt nu: Pt = P0 eαηt e−α(Ht −H0 ) .
(3.2)
In het model zal deze waarde nodig zijn om te kunnen bepalen of de maximaal toelaatbare kans op overstroming, Pmax , overschreden wordt.
3.2.2
Schade bij overstroming
De schade bij overstroming op tijdstip t = 0 is V0 . Deze schade is gebaseerd op materi¨ele schade maar ook immateri¨ele schade. Als de welvaart van een gebied groeit, zal met de jaren ook de schade groeien. De formule van de groei van de welvaart wordt analoog aan de actuele waarde berekend. Zoals in paragraaf 3.3.1 staat, is de formule voor de toekomstige waarde, als de groei γ is: toekomstige waarde = actuele waarde × eγt . Hieruit volgt dat Vt = V0 eγt . Er is echter nog een factor die de schade be¨ınvloedt. Als een dijk verhoogd wordt en het waterpeil toch hoger komt dan de dijk is dan zal de schade groter zijn dan als de dijk lager is. Dit komt doordat aangenomen wordt dat het waterpeil in het overstroomde gebied even hoog is als het laagste punt van de dijk. Bij een hogere dijk zal het waterpeil ook hoger zijn ([2]). Dit wordt weergegeven door Vt te vermenigvuldigen met eζ(Ht −H0 ) hier is ζ de groei van schade per centimeter dijkverhoging. De uiteindelijke formule voor de schade is: Vt = V0 eγt eζ(Ht −H0 ) .
(3.3)
Deze formule is vooral nodig om de verwachte schade op tijdstip t te kunnen berekenen.
3.2.3
Verwachte schade
De verwachtingswaarde van de schade op tijdstip t is gegeven door: St = Pt Vt = P0 eαηt e−α(Ht −H0 ) V0 eγt eζ(Ht −H0 ) = P0 V0 e(αη+γ)t e(α−ζ)(Ht −H0 ) = S0 eβt e−θ(Ht −H0 )
(3.4)
met S0 = P0 V0 ,
β = αη + γ,
θ = α − ζ.
De verwachte schade op tijdstip t is nodig om de totale kosten te kunnen berekenen (zie paragraaf 3.2.4). 10
3.2.4
Kostenformules
De kostenformule is opgedeeld in drie delen: de investeringskosten, de totale verwachte scha- de (welvaartskosten) en de verwachte schade na tijdstip T . In het model is het de bedoeling om deze som te minimaliseren. De verwachte schade is ook opgenomen in deze functie omdat overwogen moet worden of het de investering waard is om voor een bepaald bedrag de dijken op te hogen. Ter illustratie: als de verwachte schade in een gebied de helft van de investeringskosten is, dan is het misschien niet verstandig om die investering te doen. De investeringskosten van het ophogen van een dijk worden gegeven door: I = (C + b(htnieuw − htoud ))eλhtnieuw e−δ2 t
(3.5)
htoud is de hoogte op tijdstip t, htnieuw is de hoogte waarnaar wordt opgehoogd. htnieuw − htoud is hoeveel de dijk wordt opgehoogd (in cm). De vaste lasten van een dijkophoging worden gegeven door C, de variabele kosten door b. De term eλhtnieuw geeft weer dat de kosten om een dijk op te hogen hoger worden als houd hoger wordt. Hoe hoger een dijk wordt, hoe breder hij moet worden. De term e−δ2 t zet de kosten op een tijdstip om in de actuele waarde(zie paragraaf 3.3.1). Als een dijk niet verhoogd wordt, zijn de investeringskosten 0. De totale verwachte schade wordt gegeven door: Z Z T E= St e−δ1 t = S0
T
0
eβ t eθht
(3.6)
0
0
met β 0 = β − δ1 = αη + γ − δ1 en ht = Ht − H0 . Dit is de totale verwachte schade over een interval [0, T ] vermenigvuldigd met e−δ1 t , zodat de kosten in de actuele waarde genomen zijn. Als we er van uitgaan dat er K updates in interval [0, T ] zijn dan kan E geschreven worden als: K Z tk+1 X 0 eβ t e−θht dt E = S0 = S0 = S0 =
S0 β0
k=0 tk K X −θht
Z
e
k=0 K X
tk+1
0
eβ t dt
tk
e−θht [
k=0 K X
1 β 0 tk+1 1 0 e − 0 eβ t k ] 0 β β 0
0
e−θht [eβ tk+1 − eβ tk ].
(3.7)
k=0
Hieruit volgt dat de verwachte schade op interval [tk , tk+1 ) S0 β 0 tk+1 0 [e − eβ tk ]e−θht 0 β is. De verwachte schade is niet 0 als de dijk niet wordt opgehoogd. E=
(3.8)
De verwachte schade na tijdstip T is gegeven door ST . De totale kosten van de ophogingen met als horizon T zijn: TK
= I + E + ST
(3.9)
= (C + b(htnieuw − htoud ))eλhtnieuw e−δ2 11
S0 0 0 + 0 [eβ tk+1 − eβ tk ]e−θHt + ST β
(3.10)
3.3 3.3.1
Actuele waarde en groei van welvaart Actuele waarde
De actuele waarde(Engels: discount rate) is een economisch begrip. Het geeft aan hoeveel een bedrag in de toekomst waard is in hedendaagse valuta. Het verschil heeft te maken met rente: als je nu een euro belegt, krijg je er rente over en is het bedrag over een jaar iets meer. De formule voor de actuele waarde wordt gegeven door: actuele waarde = discount factor × toekomstige waarde.
(3.11)
De discount factor wordt gegeven door: 1 , (3.12) 1+r met r de rente (in percentage gedeeld door honderd): de vergoeding die investeerders ontvangen als beloning voor het uitlenen voor hun geld ([12]). De actuele waarde over ´e´en jaar wordt gegeven door: discount factor = δ =
actuele waarde =
toekomstige waarde . 1+r
(3.13)
Hieruit volgt dat de actuele waarde over n jaar gegeven wordt door: actuele waarde =
toekomstige waarde = δ n × toekomstige waarde. (1 + r)n
(3.14)
Dit is het geval voor ’normale’ rente, als de rente echter niet opgevraagd wordt, zal over de rente ook weer rente geheven worden. Hierdoor stijgt de rente niet lineair maar exponentieel ([13]). Dit wordt samengestelde rente genoemd. De formule voor het berekenen van de samengestelde rente is (uit [14]): A = A0 (1 +
δ nt ) n
(3.15)
met A = waarde op tijdstip t A0 = initi¨ele waarde δ = rentepercentage × 100 n = aantal keren per jaar dat de rente wordt vastgesteld t = aantal jaren Als de rente continu wordt vastgesteld geldt dat n → ∞. De formule wordt dan: A = A0 lim (1 + n→∞ δt
δ nt ) n
A = A0 e Deze formule staat gelijk aan:
toekomstige waarde = actuele waarde × eδt 12
(3.16) (3.17)
De formule voor de actuele waarde wordt: actuele waarde = toekomstige waarde × e−δt In Nederland schommelt de rente tussen de 2 en 4 %([15]). In figuur (3.2) is te zien hoe de waarde van 1000 euro zich ontwikkelt met rentes van 3, 5 en 7 procent, de doorlopende lijn is met ’normale’ rente, de gestreepte met samengestelde rente. Ik zal de rente vari¨eren tussen de 1 en 10%.
Figuur 3.2: Actuele waarde van 1000 euro na n jaar.
3.3.2
Groeipercentage van welvaart
Welvaart kan op verschillende manieren ge¨ınterpreteerd worden: op een financi¨ele(rijkdom) en niet- financi¨ele manier [16]. Onder de financi¨ele welvaart wordt verstaan dat men zich geen zorgen hoeft te maken om inkomsten en uitgaven. Het kan ook wel ge¨ınterpreteerd worden als het hebben van meer geld dan men kan uitgeven. Niet-financi¨ele welvaart betekent dat iets in overvloed aanwezig is, zoals bijvoorbeeld grondstoffen. In [17] staan echter nog twee andere vormen van welvaart, in een meer economische context. Wederom zijn er twee verschillende vormen van welvaart: reproduceerbare en niet-reproduceerbare. Onder reproduceerbare welvaart verstaat men onder andere gebouwen, middelen en software. Onder niet-reproduceerbare welvaart verstaat men bijvoorbeeld land, dat bewerkt kan worden. In een artikel over de welvaart van verschillende nieuwe lidstaten van de Europese Unie ([18]) wordt welvaart gemeten in het Bruto Binnenlands Product per hoofd van de bevolking(BBP). Het BBP is het totaal wat geproduceerd wordt in een land aan goederen en diensten. In Nederland groeide het BBP in 2006 volgens [19] met 3.25%. In jaren waar het in Nederland slecht ging met 13
de economie, zoals in 2002, groeide het BBP met 0.1%, in 2003 groeide het BBP met 0.3%. Ter vergelijking: In Estland en Letland groeide het BBP in 2006 met meer dan 11% ([18]). Die groei is in Nederland niet re¨eel, toch zal ik dit percentage vari¨eren tot van 1 tot 10%.
14
Hoofdstuk 4
Achtergrondinformatie 4.1
Waar het wel mis ging
In deze paragraaf zal ik vertellen wat er in New Orleans in 2005 en in Nederland bij de watersnoodramp in 1953 mis ging.
4.1.1
New Orleans (2005)
In augustus 2005 loeide orkaan Katrina over New Orleans heen. De orkaan richtte een dusdanige schade aan de waterkeringen aan dat hij er in combinatie met de hoge waterpeilen verantwoordelijk voor was dat een heel groot gedeelte van New Orleans onder water kwam te staan. Dit heeft zeer ernstige gevolgen gehad voor de stad, zowel financieel als niet-financieel. Wat voor dit model interessant is aan dit probleem is dat de omstandigheden in New Orleans gelijk zijn aan de omstandigheden in Nederland. New Orleans ligt evenals Nederland voor een deel onder de zeespiegel. Ook is er een dubbele dreiging van zowel de zee als de rivieren. Verder bestaat New Orleans uit moerassen die drooggelegd zijn en ingepolderd zijn. ([8]) New Orleans werd beschermd door een stelsel van waterkeringen met een lengte van 500 km. Deze waterkeringen moesten de stad beschermen tegen golven die veroorzaakt werden door orkanen in de Figuur 4.1: New Orleans na de orkaan. Foto derde categorie. Toen Katrina aan land kwam was van NOAA ([9]) het een orkaan van de derde categorie. De overstroming is veroorzaakt door twee factoren: ten eerste veroorzaakten de orkaanwinden een wateropstuwing van een uitloper van de Golf van Mexico die resulteerde in hoge waterstanden aan de oostkant van de stad. Daarna werden de hoge waterstanden voortgeplant in de rest van het gebied. Ten tweede zorgde de orkaan er zelf voor dat het water nog meer werd opgestuwd wat resulteerde in dijkdoorbraken. ([10]) In Nederland zullen geen orkanen voorkomen, maar zware stormen kunnen in Nederland dezelfde schade toebrengen. Het is belangrijk dat er goed wordt gekeken wat er in New Orleans mis ging. 15
Aantal inwoners getroffen gebied Aantal getroffenen Aantal slachtoffers Overstroomd gebied Economische schade Aantal dijkdoorbraken
Watersnood 1953 250.000 250.000 1836 2000 km2 1,5 miljard (in guldens 1953) 140
New Orleans 2005 500.000 100.000 1100 260 km2 $ 30 miljard (in US $ 2006) 30
Tabel 4.1: Schade door de watersnoodramp van 1953 in Nederland en de overstroming van New Orleans door de orkaan Katrina (gegevens over New Orleans betreffen het overstroomde deel van de stad en de schade door overstroming) ([8]).
Gebleken is dat de schade in New Orleans groter is dan verwacht. In Nederland kan men daar van leren want in veel scenario’s zijn schadeposten als schoonmaken, redding en herstel van de waterkeringen niet voldoende ingecalculeerd in de modellen. ([10])
4.1.2
Nederland (1953)
In de nacht van 31 januari op 1 februari 1953 zijn in Zeeland door een combinatie van een sterke noordwesterwind en springtij vele dijken doorgebroken. Dit veroorzaakte een overstroming van bijna heel Zeeland en delen van Zuid-Holland en Noord-Brabant. De oorzaak van de overstroming waren de te lage en zwakke dijken. Rijkswaterstaat kende dit gevaar en was bezig met plannen voor de afsluiting van de zeearmen, maar andere projecten zoals de Afsluitdijk hadden voorrang gekregen ([11]). Als gevolg van de watersnoodramp werd de Deltacommissie in het leven geroepen. Haar taak was de overheid te adviseren zodat er niet nog een watersnoodramp zou komen. Haar adviezen hebben onder andere geleid tot de afsluiting van de Hollandse IJssel, de Oosterschelde, de Grevelingen en het Haringvliet. In tabel 4.1 staan de gegevens van de watersnoodramp en de ramp van New Orleans ter vergelijking. Deze tabel is overgenomen uit [8], pagina 35.
16
Hoofdstuk 5
Werkwijze Om te kunnen kijken wat de invloed is van δ en γ, is het nodig om in korte tijd makkelijk en snel data te kunnen genereren. Het is belangrijk een model te schrijven waar je zo min mogelijk waarden hoeft te veranderen om andere data te krijgen. Om het model te schrijven heb ik er eerst voor gekozen een kloppend model te maken voor de waarden die ook in [3] genomen zijn. Om het model te verifi¨eren vergelijk ik de resultaten met de resultaten gevonden in [3]. Ook heb ik alleen voor T = 10 een zeer eenvoudig model gemaakt dat alle toegestane wegen afging en vervolgens de weg met de laagste kosten gaf. Dit model was niet optimaal doordat het alle wegen afging. Voor grotere T duurde de berekening zodanig lang (als T = 100 en er waren 20 verschillende hoogtes toegestaan, dan waren er ongeveer 20100 iteraties), dat het niet mogelijk was om die door te laten gaan.
5.1
Discretiseren en dynamisch programmeren
Om de minimale kosten te berekenen is het makkelijk om een kortste pad probleem te maken van de dijkring. Dit wordt gedaan door het probleem te discretiseren. De y-richting van het netwerk representeert de hoogte van de dijk in centimeters, de knopen liggen ieder op gelijke afstand. De x-richting staat voor de tijd in jaren. Een tak representeert de kosten om van een knoop naar de volgende te komen. Als een volgende knoop op dezelfde hoogte ligt zijn er geen kosten om de dijk te verhogen. Er zijn geen dijkverlagingen mogelijk en als de kans op overstroming op tijdstip t + 1 en hoogte hnieuw groter is dan Pmax , de maximale overstromingskans, moet er verhoogd worden. Zie voor een grafische weergave van het kortste padprobleem figuur 1.1. De oplossingsmethode van het kortste pad probleem is niet ingewikkeld. Op tijdstip t = 0 wordt gekeken hoe duur het is om naar de verschillende hoogtes op te hogen, deze waarden worden bewaard. Vervolgens wordt vanaf iedere knoop (waarvoor geldt Pt < Pmax ) gekeken hoe duur het is om van t = 1 naar een knoop op t = 2 te verhogen. Dan wordt gekeken wat de goedkoopste weg is om van de knoop op t = 0 naar iedere knoop op t = 2 te komen. Er zijn dan evenveel waarden als er op t = 2 knopen zijn. Vervolgens wordt op dezelfde manier gekeken wat de goedkoopste weg is naar de knopen op t = 3. Dit gaat door tot t = T en dan wordt uit de uitkomsten een minimale waarde gekozen ([20]). De recursieformule die bij dit probleem hoort is: Kt+1,hnieuw =
min
hnieuw ≥houd
(Kt,houd + ct,houd ,hnieuw )
17
(5.1)
Hierbij is Kt,h de totale kosten op tijdstip t en bij hoogte h. ct,houd ,hnieuw is de kosten om op tijdstip t de dijk van hoogte houd naar hnieuw te verhogen.
5.2
Programmeren
Om dit in Matlab te zetten heb ik besloten om sommige functies in en los Matlab bestand te zetten, zodat het belangrijkste deel van de methode er overzichtelijker zou uitzien. Dit waren de kosten als functie van de oude dijkhoogte, de nieuwe hoogte, het tijdstip van verhoging,de kans op overschrijding op tijdstip t en de verwachte schade bij overstroming op tijdstip t beide als functie van tijdstip t en de dijkhoogte op tijdstip t. Verder ben ik begonnen om een klein, overzichtelijk model te maken, dat over kleine waarden gaat, namelijk T = 10 en 8 verschillende toegestane hoogtes. Toen dat eenmaal werkte en klopte heb ik waarden vervangen door parameters en die bovenaan de functie een waarde gegeven, zodat die makkelijk te veranderen zijn. Later heb ik er ook voor gezorgd dat de resultaten makkelijk af te lezen zijn en er wordt een grafiek gemaakt van Pt . Om de resultaten te verifi¨eren heb ik de resultaten vergeleken met de resultaten uit [3]. Daar staan ook de resultaten van OptimaliseRing, een programma dat het HKV(Lelystad) ontwikkeld heeft en dat het CPB gebruikt. Om vervolgens te kijken wat de invloed van δ en γ is heb ik van het model een functie gemaakt die makkelijk te gebruiken was en heb ik de kosten in een grafiek uitgezet tegen verschillende waarden van δ en γ.
18
Hoofdstuk 6
Resultaten In dit hoofdstuk zal ik allereerst de resultaten van het dynamisch programmeren met vaste waarden van γ en δ geven. Ook zal ik voor andere dijkringen de uitkomsten geven. Daarna zal ik de resultaten van het dynamisch programmeren geven als γ en δ gevari¨eerd worden.
6.1
Uitkomsten voor Dijkring 10
De resultaten van het dynamisch programmeren met vaste waarden zal ik geven en daarna vergelijken met de resultaten van Professor Roos en Professor Den Hertog. De uitkomsten heb ik gekregen met het programma uit appendix A.3.1. Ik zal voor verschillende tijdshorizonnen de oplossing van het dynamisch programmeren geven. Allereerst voor tijdshorizon T = 100 en welvaartskosten 0. In de grafieken is de rode lijn Pmax de blauwe lijn is Pt , de overstromingkans op tijdstip t, van de oplossing. Naast de grafiek staan vectoren met de tijdstippen van ophoging en hoogte van de ophoging. Daaronder staan de kosten.
0 12.8 t = 40 , u = 12.8 80 6.4 Investeringskosten = 31.2029 Welvaartskosten = 0 totale kosten = 31.2029 Te zien is dat op t = 0, t = 40 en t = 80 opgehoogd Figuur 6.1: T = 100 en geen welvaartskos- wordt. Dit is een periodieke oplossing. De totale ten. kosten van deze oplossing zijn 31,2029. Dit is gelijk aan de investeringskosten, de welvaartskosten zijn 0.
19
Voor T = 100 maar met de welvaartskosten is de oplossing: · t=
0 59
¸
· , u=
19.2 38.4
¸
Investeringskosten = 33.6596 Welvaartskosten = 20.2807 totale kosten = 53.9403
Figuur 6.2: T = 100.
Het verschil tussen deze oplossing en de vorige is dat hier wel de welvaartskosten berekend zijn. Er wordt ook nog maar ´e´en keer opgehoogd. De investeringskosten zijn wel hoger.
Voor T = 300 is de oplossing: t=
0 59 117 174 233
, u=
19.2 51.2 57.6 57.6 64
Investeringskosten = 35.2272 Welvaartskosten = 20.3077 totale kosten = 55.5348
Figuur 6.3: T = 300.
Deze oplossing is bijna periodiek met een periode van ongeveer 59 jaar. De totale kosten zijn iets hoger dan bij T = 100. De kans op overstroming gaat naar 0, dit komt doordat de welvaartskosten lager zijn als de kans op overstroming lager is.
20
De oplossing voor T = 300 zonder welvaartskosten: t=
Figuur 6.4: T = 300 zonder welvaartskosten.
0 40 80 120 159 200 239 280
, u=
12.8 12.8 12.8 12.8 12.8 12.8 12.8 6.4
Investeringskosten = 31.6598 Welvaartskosten = 0 totale kosten = 31.5698
Deze oplossing verschilt duidelijk van de vorige. Hier wordt ook duidelijk wat de welvaartskosten eigenlijk doen in de formule. In deze oplossing worden de dijken iedere keer verhoogd als het nodig is. In de vorige oplossing was dat niet zo. Dat komt doordat de welvaartskosten laag zijn als de kans op overstroming klein is. Immers St = Pt Vt . Er wordt dan niet opgehoogd omdat het nodig is (dus als Pt = Pmax ). De dijk wordt opgehoogd zodat de kans op overstroming z´o klein is dat de verwachte schade bijna 0 is. Investeren is dan beter dan afwachten tot Pmax behaald is. De oplossing voor T = 300 met een kleinere stapgrootte is: t=
0 59 117 175 233
, u=
19.20 51.84 57.60 57.60 57.60
Investeringskosten = 35.271 Welvaartskosten = 20.264 totale kosten = 55.5349 Figuur 6.5: T = 300 en een kleinere stapgrootte. Deze oplossing is bijna gelijk aan de vorige, alhoewel de stapgrootte 4 keer kleiner is.
21
De oplossing voor T = 500 is: t=
0 59 117 174 232 290 348 405 457
, u=
19.2 51.2 57.6 57.6 57.6 57.6 57.6 54.4 44.8
Figuur 6.6: T = 500. Investeringskosten = 35.2277 Welvaartskosten = 20.3079 totale kosten = 55.25356 Deze resultaten verschillen niet veel met de resultaten uit [3]. De oplossing van figuur 6.1 is zelfs helemaal gelijk. De verschillen tussen de resultaten zijn vrij klein en te verklaren doordat op t = 60 volgens mijn model Pmax is overschreden. De dijk wordt op t = 59 opgehoogd. In het model uit [3] wordt op t = 60 opgehoogd. Ik heb wel de kosten volgens mijn model maar met de tijd en stapgroottes uit [3] berekend en die zijn wel gelijk aan de kosten uit [3].
6.2
Uitkomsten voor andere dijkringen
Om verschillen tussen dijkringen te laten zien heb ik voor nog drie andere dijkringen de uitkomsten bepaald. Dit heb ik gedaan voor de dijkringen Alblasserwaard en Vijfheerenlanden, Heerewaarden, en de Gelderse vallei. Deze heb ik gekozen vanwege hun uiterste waarden. In Alblasserwaard en 1 Vijfheerenlanden is de kans op overstroming 906 , terwijl de schade bij overstroming op tijdstip t = 0 vrij groot is. In Heerewaarden is het precies andersom: de schade is klein en de kans op overstroming ook. In de Gelderse vallei, een verstedelijkt gebied, is de kans op overstroming klein en de schade bij overstroming groot. De data staan ter vergelijking in tabel 6.1. Voor alle oplossingen is T = 300 genomen.
22
De oplossing voor Alblasserwaard en Vijfheerenlanden is: t=
0 56 108 161 214 267
, u=
51.2 51.2 51.2 51.2 51.2 57.6
Investeringskosten = 885.7257 Welvaartskosten = 212.9379 Figuur 6.7: Alblasserwaard en Vijfheerenlanden.
totale kosten = 1098.6636
De oplossing voor Heerewaarden is: t=
0 60 121 186 246
, u=
25.6 25.6 70.4 76.8 76.8
Investeringskosten = 10.027 Welvaartskosten = 1.8698 totale kosten = 11.8968
Figuur 6.8: Heerewaarden.
23
dijkring 10: Mastenbroek C b λ α η ζ V0 P0 γ δ H0
dijkring 40: Heerewaarden
dijkring 45: Gelderse Vallei
16.6939 0.6258 0.0014 0.033027 0.320 0.003774 1564.9
dijkring 16: Alblasserwaard en Vijfheerenlanden 324.6287 2.1304 0.01 0.0574 0.76 0.002032 22656.5
5.8601 0.1000 0.0026 0.025321 0.422 0.003905 43.5
3.4375 0.1375 0.0069 0.033027 0.320 0.002397 10421.2
1 2270
1 906
1 2855
1 6120
0.02 0.04 0
0.02 0.04 0
0.02 0.04 0
0.02 0.04 0
Tabel 6.1: Data van de verschillende dijkringen
De oplossing voor de Gelderse vallei is: t=
0 47 95 142 191 246
, u=
64 38.4 38.4 38.4 44.8 44.8
Investeringskosten = 22.3706 Welvaartskosten = 11.3766 Figuur 6.9: Gelderse vallei.
totale kosten = 33.7473
De oplossing van Alblasserwaard en Vijfheerenlanden is heel duur doordat de kans op overstroming groot mag blijven en de verwachte schade op tijdstip t = 0 ook groot is. Hierdoor zijn zowel de inversteringskosten als de welvaartskosten zeer hoog. De dijkring van Heerewaarden is wat data betreft juist het tegenovergestelde van Alblasserwaard en Vijfheerenlanden: de kosten zijn laag en de kans op overstroming is klein. Hierdoor wordt in plaats van ´e´en keer twee keer opgehoogd omdat de kans op overstroming bijna Pmax overschrijdt. De Gelderse vallei is een groot, verstedelijkt gebied, hierdoor is de schade bij overstroming groot, maar doordat de kans op overstroming maar heel klein mag zijn, zijn de totale kosten niet heel hoog. 24
6.3
Uitkomsten voor verschillende constanten
In deze paragraaf zal ik de oplossing met verschillende δ en γ geven. De uitkomsten heb ik verkregen met de programma’s uit appendices A.3.2 en A.3.3. Ik heb bij de programma’s gekozen voor tijdshorizon T = 300 en de kleine, vaste stapgrootte die ook gebruikt is bij de oplossing in figuur 6.5.
6.3.1
Uitkomsten voor verschillende δ
Allereerst ben ik met gelijke γ δ gaan vari¨eren, doordat δ1 = δ2 hoeft er geen onderscheid gemaakt te worden tussen de δ’s. Bij een stijgende δ wordt, zoals ook in figuur 3.3.1 te zien is, de actuele waarde minder. Doordat de actuele waarde lager wordt worden de totale kosten ook lager. De resultaten zijn te zien in figuur (6.10). Met links de kosten voor δ vari¨erend tussen 1 en 10 procent, met een stapgrootte van 1 procent. Rechts de kosten voor een voor Nederland re¨ele δ. Het is duidelijk te zien dat voor een groter wordende δ de kosten minder worden. De daling is te verklaren doordat bij hoge inflatie de rente ook hoger zal zijn ([13]). Als de rente hoger wordt zal de actuele waarde op zijn beurt juist weer minder worden. Daardoor worden de kosten lager. In Japan is de rente momenteel 0,5%. Daar zullen de kosten heel hoog zijn als de rente zo laag blijft. Maar een lage rente betekent ook dat het makkelijker is om te investeren.
Figuur 6.10: Links de kosten voor δ 1-10 %, rechts voor δ 2-5% met een kleinere stapgrootte.
6.3.2
Uitkomsten voor verschillende γ
Allereerst kijken we wat een grotere γ betekent. Een grotere waarde voor γ betekent dat de welvaart meer groeit, een gebied wordt meer waard. Onder die omstandigheden rendeert het om vaker te investeren in ophogingen van dijken. Dit is inderdaad te zien in figuur 6.11. Tot 4% stijgen de kosten niet heel veel, maar daarna gaat het heel snel. Zoals in paragraaf (3.3.2) staat is een groei van meer dan 4% in Nederland echter niet re¨eel. Aangezien in bijvoorbeeld Oost-Europese landen de welvaart wel zo snel groeit, zullen daar de kosten van een dergelijke operatie erg hoog zijn. De stijging van de kosten bij de groei van 25
welvaart is economisch te verklaren doordat bij lagere economische groei de inflatie ook laag is ([21]). Hierdoor zullen de kosten laag blijven bij kleine groei.
Figuur 6.11: Links de kosten voor γ 1-10 %, rechts voor γ 1-3% met een kleinere stapgrootte.
6.3.3
Uitkomsten voor verschillende γ en δ
In figuur 6.12 zijn de kosten voor verschillende γ en δ berekend. Het verloop is als verwacht met een maximum bij kleine δ en grote γ en een minimum bij grote δ en kleine γ. Het is ook te zien dat de kosten bij groei in welvaart sneller groeien dan de kosten dalen bij de stijgende inflatie.
Figuur 6.12: De kosten voor verschillende γ en δ. 26
Hoofdstuk 7
Conclusie Allereerst zal ik het hebben over de verschillen tussen Eijgenraam en Van Dantzig, daarna over de formules en het model en ten slotte zal ik het nog over de resultaten hebben. De verschillen tussen de modellen van Eijgenraam en Van Dantzig zijn groot, vooral het verschil in de investeringskosten. Verder valt het op dat Van Dantzig slechts de verloren hoogte van de dijk aan de waterpeilstijging en dijkverzakking wil goedmaken met een verhoging. Eijgenraam let er ook op dat als een gebied rijker wordt dat de dijken dan meer opgehoogd worden zodat de verwachte schade kleiner wordt. Het belangrijkste verschil in de concepten van de modellen is dat Van Dantzig slechts uitzoekt hoeveel het kost om op te hogen en niet wanneer. Eijgenraam doet dit wel. De formules die nodig zijn om het model te maken zijn de volgende: Pt = P0 eαηt e−α(Ht −H0 ) γt ζ(Ht −H0 )
Vt = V0 e e TK
(7.1) (7.2)
= I + E + ST = (C + b(htnieuw − htoud ))eλhtnieuw e−δ2 S0 0 0 + 0 [eβ tk+1 − eβ tk ]e−θHt + ST β
(7.3)
De kans op overstroming is 7.1, de schade is 7.2 en de totale kosten zijn 7.3. Het maken van het model was vrij ingewikkeld en het duurde vrij lang voor ik een definitief, kloppend model had. De moeilijkheden zaten in het algoritme dat makkelijk is om te begrijpen en uit te leggen, maar ingewikkelder om te programmeren. Doordat ik de kostenfunctie los had geprogrammeerd van het algoritme, was het niet erg toen de kostenfunctie niet klopte, zoals in het begin het geval was. Daardoor kon de functie makkelijk gewijzigd worden en kostte dit maar weinig extra tijd. Eerst heb ik ervoor gezorgd dat het algoritme werkte daarna heb ik de formule gebruikt om een nieuw programma te maken waarin de data voor verschillende waarden makkelijk konden worden bekeken. De verschillen werden daardoor makkelijk zichtbaar. De resultaten met de waarden γ = 0.02 en δ = 0.04 waren bijna gelijk aan de resultaten uit [3]. De verschillen ontstonden doordat in het model dat ik gebruikte niet accepteerde dat op het 27
volgende tijdstip Pmax overschreden werd. Hieruit volgde dat het model klopte en ik verder kon gaan met de variatie van γ en δ. Aangezien het belangrijkste deel al werkte was dit niet lastig. De resultaten van de verschillende γ en δ kwamen uit zoals te verwachten was: de kosten stegen voor hogere γ en de kosten daalden door hogere δ. De gevolgen zijn dat als de Nederlandse welvaart groeit, worden de kosten om dijken te verhogen hoger: het is dan namelijk meer waard om op te hogen aangezien de welvaartskosten groter zijn. Als er weinig inflatie is en de rente laag is worden de kosten hoger. Deze kosten groeien minder snel dan bij groter wordende γ. De uitkomsten van het project geven aan dat het in de toekomst zeker mogelijk is om Nederland droog te houden met behulp van dijkverhogingen maar het is verstandig om uit te kijken naar andere oplossingen. Ook is het belangrijk dat de kosten bij overstroming goed worden ingeschat. Na de overstroming in New Orleans is gebleken dat veel factoren niet zijn meegenomen in bestaande modellen. Het gaat om kosten van schoonmaken, redding en herstel van de waterkeringen. Aan de hand van de gegevens uit New Orleans kan dit bepaald worden.
28
Bibliografie [1] D. van Dantzig: Economic decision problems for flood prevention,(1956) [2] C.J.J. Eijgenraam: Optimal safety standards for dike-ring areas, Centraal Planbureau (2006) [3] Dick den Hertog, Kees Roos: Computing safe dike heights at minimal costs, (2007) [4] Jarl Kind, Carlijn Bak: Optimalisatie van veiligheid tegen overstromen, (2007) [5] www.wikipedia.nl: Dijkring [6] www.wikipedia.nl: Mastenbroek [7] Polders: een theater van land en water; Polder Mastenbroek(1364), Nederlands Architectuur Instituut, http://static.nai.nl/polders/nl/index.html [8] Matthijs Kok, Rob Theunissen, Bas Jonkman, Han Vrijling: Schade door overstoming: Ervaringen uit New Orleans, TU Delft en HKV Lijn in water (2006) [9] National Oceanographic and Atmospheric Administration: www.photolib.noaa.gov [10] W. Kanning, M. Kok en J. Vrijling: New Orleans: toen de dijken bezweken, H2 O nummer 3 - 2007, HKV [11] www.wikipedia.nl; Watersnood van 1953 [12] Brealy, Myers: Principles of corporate finance, McGraw-Hill (2000), versie uit dictaat Inleiding Technische Bestuurskunde, Pieter W.G. Bots, TU Delft [13] www.wikipedia.nl: Rente [14] www.wikipedia.com: Compound interest [15] Tabel T1.1, Rentetarieven van de Europese Centrale Bank en De Nederlandsche Bank, De Nederlandsche Bank (versie 30 januari 2008) [16] www.wikipedia.nl: Rijkdom [17] P. R. Brahmananda: Growth rates of wealth and capital in the US, Business Line, Internet edition (2001) [18] Bart Meijer: Nederlanders relatief rijk, Webmagazine CBS (4 juli 2007) [19] Macro Economische Verkenning 2007, CPB (september 2006) 29
[20] Hillier, Lieberman: Introduction to operations research, McGraw-Hill (2005) [21] Martin Mellens: Economische groei en inflatie, Webmagazine CBS (21 juni 1999)
30
Bijlage A
MatLab-modellen A.1 A.1.1
Gebruikte functies Overstromingskans
%overstromingskans op tijdstip t en hoogte ht function P = P(t,Ht,u,P0,alfa,eta,H0) P = P0 * (exp(alfa*eta*t))* (exp(-alfa*(Ht-H0)*u));
A.1.2
Schade
%loss by flooding at time t function V = V(t,Ht,u,V0,gamma,dzeta,H0) V = V0 * (exp(gamma*t))* (exp(dzeta*(Ht-H0)*u));
A.1.3
Verwachte schade op tijdstip t
%verwachte schade op tijdstip t function S = S(t,Ht,u,S0,beta,theta,H0) S = S0 * (exp((beta)*t))* (exp(-theta*(Ht-H0)*u));
A.1.4
Kosten
%Kosten om een dijk te verhogen van hnieuw naar houd op tijdstip t function K = KH(t,houd,hnieuw,u,beta2,S0,theta,C,b,labda,delta2) K = zeros(3,1); E = (S0/beta2) * ((exp(beta2*(t+1))) - (exp(beta2*t))) * (exp(-theta*(u*hnieuw))); %welfare costs I = (C + b*(hnieuw-houd)*u) * (exp(labda*hnieuw*u)) * (exp(-delta2*t)); %investment costs %E = 0; K(1,1) = E + I; %total costs K(2,1) = I; K(3,1) = E; if houd == hnieuw 31
K(2,1) = 0; K(1,1) = K(3,1); end
A.2
Algoritme
%functie van het algoritme function [DPH2,DPH2R] = DPH2(T,gamma,delta1,ho,u) H0 = 0; %beginhoogte dijk %structurele groei van waterniveau(cm/jaar) eta = 0.320; %parameter exponentiele verdeling voor extreme waterniveaus alfa = 0.033027; %exceedance probability op tijdstip t = 0 P0 = 1/2270; %limiet voor exceedance probability Pmax = P0; %verlies door overstroming op tijdstip t(miljoen euros) V0 = 1564.9; C = 16.6939; b = 0.6258; labda = 0.0014; %groei van verlies per cm dijkverhoging(1/cm) dzeta = 0.003774; delta2 = delta1; S0 = P0*V0; beta = alfa*eta + gamma; beta2 = alfa*eta + gamma - delta1; theta = alfa - dzeta; TK = zeros(3,1); %voor t = 0 kijken wat de minimale waarden zijn om van 0 naar hnieuw te %verhogen t = 0; L = zeros(ho,T+2); Weg = zeros(ho,T+2); for hnieuw = 0:(ho-1) %kijken of op het nieuwe tijdstip de nieuwe hoogte de kans kleiner is %dan Pmax if P(t+1,hnieuw,u,P0,alfa,eta,H0) <= Pmax K = KH(t,0,hnieuw,u,beta2,S0,theta,C,b,labda,delta2); L(hnieuw+1,1) = K(1,1); Weg(hnieuw+1,1) = 0; % in matrix Weg aangeven wat de oude hoogte was end end
32
%voor t van 1 tot de eindgrens kijken wat de minimale kosten zijn voor %verhoging op t naar hnieuw for t = 1:T for houd = 0:(ho-1) for hnieuw = 0:(ho-1) %kijken of op het nieuwe tijdstip de nieuwe hoogte de kans kleiner is %dan Pmax if (P(t,houd,u,P0,alfa,eta,H0) <= Pmax) && (P(t+1,hnieuw,u,P0,alfa,eta,H0) <= Pmax) K = KH(t,houd,hnieuw,u,beta2,S0,theta,C,b,labda,delta2); %als er op dezelfde hoogte gebleven mag worden zijn de kosten 0 s = L(houd+1,t) + K(1,1); %kosten berekenenen if L(hnieuw+1,t+1) == 0 %als er nog geen kosten geweest zijn L(hnieuw+1,t+1) = s; %wordt de waarde in de matrix gezet Weg(hnieuw+1,t+1) = houd; end if s < L(hnieuw+1,t+1) %er wordt gekeken of er een kleinere L(hnieuw+1,t+1) = s; %waarde is dan de zojuist berekende Weg(hnieuw+1,t+1) = houd;%waarde, als dat niet zo is, wordende end %waarde en de oude hoogte opgeslagen end end end end %in de matrix Weg wordt gekeken wanneer met hoeveel wordt opgehoogd en in %matrix route opgeslagen route = ones(4,ho-1); w = ho-1; z = 1; for ht = 1:T-1 s = T - ht; if not(Weg(w+1,s) == w) route(1,z) = s-1; route(2,z) = Weg(w+1,s); route(3,z) = Weg(w+1,s+1); route(4,z) = route(3,z) - route(2,z); w = Weg(w+1,s); z = z+1; end end %uitkomsten omzetten van matrix Weg in twee vectoren U en t t = zeros(z-1,1); U = zeros(z-1,1); for s = 1:(ho-1) 33
w = (ho) - s; if not((route(2,w) == 1) && (route(3,w) == 1)) t(z-w,1) = route(1,w); U(z-w,1) = route(4,w)*u; Z = KH(route(1,w),route(2,w),route(3,w),u,beta2, S0,theta,C,b,labda,delta2); TK(2,1) = TK(2,1) + Z(2,1); if (w > 1) TK(3,1) = TK(3,1) + ((S0/beta2) * ((exp(beta2*(route(1,w-1)))) - (exp(beta2*route(1,w)))) * (exp(-theta*(u*route(3,w))))); end if (w == 1) TK(3,1) = TK(3,1) + ((S0/beta2) * ((exp(beta2*T)) - (exp(beta2*route(1,w)))) * (exp(-theta*(u*route(3,w))))); end TK(1,1) = TK(2,1) + TK(3,1); end end %grafiek van Pt maken hoog = zeros(T+1,1); PMAX = zeros(T+1,1); tijd = zeros(T+1,1); PT = zeros(T+1,1); ht = 0; w = 1; for i = 0:T if (w < z) && (i == t(w,1)) ht = ht + U(w,1); w = w + 1; end hoog(i+1,1) = ht; PMAX(i+1,1) = Pmax; tijd(i+1,1) = i; end for i = 1:T+1 PT(i,1) = P(i,hoog(i,1)/u,u,P0,alfa,eta,H0); end ST = S(T,ho-1,u,S0,beta,theta,H0); DPH2R = zeros(4,ho-1); DPH2(1,:) = t’; DPH2(2,:) = U’; DPH2(3,1) = TK(1,1); DPH2(3,2) = ST; 34
DPH2R(1,:) DPH2R(2,:) DPH2R(3,:) DPH2R(4,:)
A.3 A.3.1
= = = =
route(1,:); route(2,:); route(3,:); route(4,:);
Toepassingen Algoritme
%uitvoering van het algoritme clf; clc; clear; %function [DPH2,DPH2R] = DPH2(T,gamma,delta1,ho,u) H0 = 0; %beginhoogte dijk %structurele groei van waterniveau(cm/jaar) eta = 0.320; %parameter exponentiele verdeling voor extreme waterniveaus alfa = 0.033027; %exceedance probability op tijdstip t = 0 P0 = 1/2270; %limiet voor exceedance probability Pmax = P0; %verlies door overstroming op tijdstip t(miljoen euros) V0 = 1564.9; C = 16.6939; b = 0.6258; labda = 0.0014; %groei van verlies per cm dijkverhoging(1/cm) gamma = 0.02; delta1 = 0.04; dzeta = 0; %0.003774; delta2 = delta1; S0 = P0*V0; beta = alfa*eta + gamma; beta2 = alfa*eta + gamma - delta1; theta = alfa - dzeta; ho = 40; u = 6.4; TK = zeros(3,1); T = 300; %voor t = 0 kijken wat de minimale waarden zijn om van 0 naar hnieuw te %verhogen t = 0; L = zeros(ho,T+2); Weg = zeros(ho,T+2); 35
for hnieuw = 0:(ho-1) %kijken of op het nieuwe tijdstip de nieuwe hoogte de kans kleiner is %dan Pmax if P(t+1,hnieuw,u,P0,alfa,eta,H0) <= Pmax K = KH(t,0,hnieuw,u,beta2,S0,theta,C,b,labda,delta2); L(hnieuw+1,1) = K(1,1); Weg(hnieuw+1,1) = 0; % in matrix Weg aangeven wat de oude hoogte was end end %voor t van 1 tot de eindgrens kijken wat de minimale kosten zijn voor %verhoging op t naar hnieuw for t = 1:T for houd = 0:(ho-1) for hnieuw = 0:(ho-1) %kijken of op het nieuwe tijdstip de nieuwe hoogte de kans kleiner is %dan Pmax if (P(t,houd,u,P0,alfa,eta,H0) <= Pmax)&& (P(t+1,hnieuw,u,P0,alfa,eta,H0) <= Pmax) K = KH(t,houd,hnieuw,u,beta2,S0,theta,C,b,labda,delta2); %als er op dezelfde hoogte gebleven mag worden zijn de kosten 0 s = L(houd+1,t) + K(1,1); %kosten berekenenen if L(hnieuw+1,t+1) == 0 %als er nog geen kosten geweest zijn L(hnieuw+1,t+1) = s; %wordt de waarde in de matrix gezet Weg(hnieuw+1,t+1) = houd; end if s < L(hnieuw+1,t+1) %er wordt gekeken of er een kleinere L(hnieuw+1,t+1) = s; %waarde is dan de zojuist berekende Weg(hnieuw+1,t+1) = houd;%waarde, als dat niet zo is, worden de end %waarde en de oude hoogte opgeslagen end end end end %in de matrix Weg wordt gekeken wanneer met hoeveel wordt opgehoogd en in %matrix route opgeslagen route = ones(4,ho-1); w = ho-1; z = 1; for ht = 1:T-1 s = T - ht; if not(Weg(w+1,s) == w) route(1,z) = s-1; %tijdstip van ophoging route(2,z) = Weg(w+1,s); %houd route(3,z) = Weg(w+1,s+1); %hnieuw route(4,z) = route(3,z) - route(2,z); %hoeveel niveaus opgehoogd 36
w = Weg(w+1,s); z = z+1; end end %uitkomsten omzetten van matrix Weg in twee vectoren U en t t = zeros(z-1,1); U = zeros(z-1,1); for s = 1:(ho-1) w = (ho) - s; if not((route(2,w) == 1) && (route(3,w) == 1)) t(z-w,1) = route(1,w); U(z-w,1) = route(4,w)*u; Z = KH(route(1,w),route(2,w),route(3,w),u,beta2,S0,theta,C,b,labda,delta2); TK(2,1) = TK(2,1) + Z(2,1); %investment costs if (w > 1) TK(3,1) = TK(3,1) + ((S0/beta2) * ((exp(beta2*(route(1,w-1)))) - (exp(beta2*route(1,w)))) * (exp(-theta*(u*route(3,w))))); end if (w == 1) TK(3,1) = TK(3,1) + ((S0/beta2) * ((exp(beta2*T)) - (exp(beta2*route(1,w)))) * (exp(-theta*(u*route(3,w))))); end TK(1,1) = TK(2,1) + TK(3,1); end end %grafiek van Pt maken hoog = zeros(T+1,1); PMAX = zeros(T+1,1); tijd = zeros(T+1,1); PT = zeros(T+1,1); ht = 0; w = 1; for i = 0:T hoog(i+1,1) = ht; %vector met lengte op tijdstip i-1 PMAX(i+1,1) = Pmax; %vector met overal PMAx tijd(i+1,1) = i; %vector met de tijd if (w < z) && (i == t(w,1)) ht = ht + U(w,1); w = w + 1; end end %vector met overschrijdingskans op tijdstip met de hoogte volgens de 37
%optimale oplossing for i = 1:T+1 PT(i,1) = P(i-1,hoog(i,1)/u,u,P0,alfa,eta,H0); end plot(tijd,PT); hold on plot(tijd,Pmax,’r-’); %samenvattingsmatrix ST = S(T,ho-1,u,S0,beta,theta,H0); DPH2(1,:) = t’; DPH2(2,:) = U’; DPH2(3,1) = TK(1,1); DPH2(3,2) = ST;
A.3.2
%verwachte schade
Rente
clc; clf; clear; T = 300; gamma = 0.02; delta1 = 0.04; ho = 128; u = 1.92; grafiek = zeros(10,1); percentage = zeros(10,1); for i = 1:10 delta1 = i*0.01; %voor 1-10% i loopt van 1-10 %delta1 = 0.02 + i*0.001; %voor 2-5% i loopt van 0-30 [A,B] = DPH2(T,gamma,delta1,ho,u); grafiek(i,1) = A(3,1); percentage(i,1) = delta1; end plot(percentage,grafiek,’mo-’);
A.3.3
Groei van welvaart
clc; clear; T = 300; gamma = 0.02; delta1 = 0.04; 38
u = 1.92; ho = 128; grafiek = zeros(10,1); percentage = zeros(10,1); for i = 1:10 %gamma = 0.01 + i*0.001; %voor gamma 1-3% i loopt van 0-20 gamma = i*0.01; %voor gamma 1-10% i loopt van 1-10 [A,B] = DPH2(T,gamma,delta1,ho,u); grafiek(i,1) = A(3,1); percentage(i,1) = gamma; end for i = 1:8 V(i,1) = i; end plot(percentage,grafiek,’go-’);
A.3.4
Vari¨ eren δ en γ
clc; clf; clear; T = 300; gamma = 0.02; delta1 = 0.04; ho = 128; u = 1.92; grafiek = zeros(21,21); percg = zeros(21,1); percd = zeros(21,1); for i = 0:20 for j = 0:20 gamma = 0.01 + i*0.001; %voor gamma 1-3% i loopt van 0-20 %gamma = i*0.01; delta1 = 0.03 + j*0.001; %voor 3-5% i loopt van 0-20 %delta1 = j*0.01; [A,B] = DPH2(T,gamma,delta1,ho,u); grafiek(i+1,j+1) = A(3,1); percg(i+1,1) = gamma; percd(j+1,1) = delta1; end 39
end surf(percd,percg,grafiek);
40