UDK/UDC: 359.42:628.144
Prejeto/Received: 18.05.2017 Sprejeto/Accepted: 02.11.2017
Izvirni znanstveni članek - Original scientific paper
Uporaba metode spektralnega razbitja grafov za zasnovo merilnih
OBMOČIJ V VODOVODNEM OMREŽJU
An APPLICATION OF SPECTRAL GRAPH PARTITION FOR DESIGNING DISTRICT METERED AREAS IN WATER SUPPLY NETWORKS
Daniel Kozelj1'*, Marjan Gorjup2, Marjeta Kramar Fijavž1
'Fakulteta za gradbeništvo in geodezijo, Univerza v Ljubljani, Jamova 2, 1000 Ljubljana 2Petrol d.d., Dunajska cesta 50, 1000 Ljubljana, Slovenija
Izvleček
Pri upravljanju vodovodnih sistemov se spopadamo z velikim deležem neprodane vode, pri čemer največji delež odpade na dejanske vodne izgube. Za zmanjšanje količine neprodane vode lahko v vodovodnem sistemu oblikuj emo merilna območj a, kjer se na vtoku merita parametra tlak in pretok vode. Zasnova merilnih območij upravljalcem omogoča hitrejše ugotavljanje vodnih izgub in učinkovitejši nadzor nad delovanjem celotnega sistema. Prikazana je metoda za razdelitev kompleksnih vodovodnih omrežij na posamezna območja, ki uporablja algoritem za spektralno razbitje grafov. Za določitev ustreznih realnih rešitev smo vključili tudi hidravlične vhodne podatke. Opisano metodo smo testirali na realnem primeru.
Ključne besede: vodovodni sistem, merilna območja, algoritem, teorija grafov, spektralno razbitje.
One of the major challenges in managing water supply systems is the high percentage of non-revenue water, wherein the largest contributor is the actual loss of water. In order to decrease the volume of the non-revenue water, district metered areas (DMAs) are introduced, where hydraulic parameters such as pressure and flow are measured. Dividing water distribution systems into different DMAs allows the water management to identify water losses and effective control over the whole system more quickly. We present a method for partitioning complex networks, which uses the algorithm for spectral partitioning of a graph. In order to generate hydraulically suitable solutions also hydraulic input data was also used. The proposed method was tested on a real case study.
Keywords: water supply system, district metered areas, algorithm, graph theory, spectral partition.
* Stik / Correspondence: daniel.kozelj@fgg.uni-lj.si
© Kozelj D. et al.; Vsebina tega članka se sme uporabljati v skladu s pogoji licence Creative Commons Priznanje avtorstva -Nekomercialno - Deljenje pod enakimi pogoji 4.0.
© Kozelj D. et al.; This is an open-access article distributed under the terms of the Creative Commons Attribution - Non Commercial - ShareAlike 4.0 Licence.
Abstract
1 Uvod
Merilno območje (angl. District Metered Area, DMA) je zaključen del vodovodnega sistema, za katerega lahko izmerimo količino vode, ki vanj vteka oziroma izteka (Morrison et al., 2007). Koncept uvedbe merilnih območij je ponazorjen na sliki 1. Na shemi so prikazana na novo ustvarjena merilna območja, lokacije merilnih mest (M) na vtokih v območja ter zaprte čevi (rdeče prečrtane), ki dovoljujejo pretok vode le preko merjenih čevnih odsekov.
Na podlagi meritev izračunamo skupno bi-lančo vode v opazovanem merilnem območju v izbranem časovnem obdobju (Lambert, 2003). Z vpeljavo merilnih območij v vodovodni sistem enostavneje razumemo delovanje čelotnega sistema. S pomočjo evidentiranih vrednosti tlakov in pretokov hitreje prepoznamo novo nastale vodne izgube zaradi poskodb čevi. Ker so posamezna merilna območja ločena od preostalega dela vodovodnega sistema, jih lahko enostavneje upravljamo in nadziramo (Morrison, 2004). Običajno na merilnih mestih merimo pretok in tlak, mozno pa je tudi merjenje parametrov, povezanih s kakovostjo vode.
Poleg vzpostavitve merilnih območij lahko vodovodne sisteme upravljamo in nadziramo tudi s pomočjo numeričnih modelov, t. i. hidravličnih modelov vodovodnih sistemov. Hidravlični model vodovodnega sistema je sestavljen iz numeričnega opisa objektov, ki jih delimo na linijske in točkovne elemente. Med linijske elemente pristevamo čevi, ventile in črpalke, medtem ko so točkovni elementi vo-zlisča, vodni viri in vodohrani. Taksna struktura gradnikov hidravličnega modela ustreza matematičnemu objektu, imenovanem graf, pri čemer linijski elementi predstavljajo povezave in točkovni elementi točke ali vozlišča grafa. Matematično področje, ki se ukvarja s takimi objekti, se imenuje teorija grafov (Wilson & Watkins, 1997). Namen nasega dela je poiskati in primerno prirediti algoritme iz teorije grafov
za učinkovito razdelitev vodovodnega omrežja na merilna območja, pri čemer bi radi dodatno upoštevali tudi hidravlične kriterije.
Mnogi raziskovalci so si pri zasnovi merilnih območij ze pomagali s teorijo grafov. Dobljene metode so bile sprva bolj robustne in so podale manj ustrezne resitve končno zasnovanih merilnih območij, tudi hidravlični kriteriji so bili manj upostevani. Prvi so se s tem problemom ukvarjali Tzačhkov, Alčočer-Yamanaka in Bourguett-Ortiz (Tzatčhkov et al., 2006). Za določitev neodvisnih območij so uporabili t. i. načelo LIFO (angl. last in first out). Sempewo, Pathirana in Vairavamoorthy (Sempewo et al., 2009) so razvili orodje za oblikovanje območij, ki so ga poimenovali METIS. To orodje je ze imelo moznost vključitve nekaterih hidravličnih kriterijev v obliki utezi na povezavah grafa.
Di Nardo in soavtorji (Di Nardo et al., 2014) so uporabili Dijkstrov algoritem iskanja naj-krajsih poti, za utezi so vzeli energijo v sistemu. Gomes, Sa Marques in Sousa (Gomes et al., 2012) so uporabili Floyd-Warshallov algoritem. Pri novih zasnovah so preverjali spremenjene vrednosti hitrosti vode in tlaka. Upostevana je tudi investičijska vrednost zaradi spremenjenega sistema. Herrera in soavtorji (Herrera et al., 2012) so izbrali algoritem kopičenja, ki pa je omejen po stevilu merilnih območij.
Alvisi in Frančhini (Alvisi & Frančhini, 2014) sta poskusila z algoritmom iskanja v sirino. V metodi z algoritmom kumulativno sestevata porabo vode in hkrati isčeta naj-krajso pot od vodnega vira do preostalih vo-zlisč. Diao, Zhou in Raučh (Diao et al., 2013) so vozlisča pripisali posameznih skupnostim, ki so kasneje zdruzene, če ne krsijo pogoja o velikosti.
Ferrari, Savič in Beččiu (Ferrari et al., 2014) z rekurzivno uporabo bisekčije razpola-vljajo območja, dokler obe novi podobmočji ne ustrezata predpisanim mejam glede porabe
Slika 1: Shema enovitega in na merilna območja razdeljenega vodovodnega sistema.
Figure 1: Scheme of a uniform and DMA water supply system.
vode. Tudi tu je uporabljen algoritem iskanja v Širino. Hajebi in soavtorji (Hajebi et al., 2014) so predlagali metodo, v kateri je upoštevana večparametrska optimizacija in medsebojna neodvisnost območij. Galdiero (Galdiero, 2015) je ponovno uporabil Floyd-Warshallov algoritem. Pri metodi je uposteval več hidravličnih in ekonomskih kriterijev ter izmed dobljenih resitev nato izbral optimalno.
V nasem delu bomo prikazali uporabo algoritma, ki ga je razvil Hespanha (He-spanha, 2004) in temelji na spektralnem razbitju grafa. V naslednjem razdelku najprej razložimo osnovne enačbe hidravlike, ki nastopajo pri izračunu hidravličnih spremenljivk v vodovodnem sistemu. Nato predstavimo algoritem spektralnega razbitja grafov ter opisemo, kako ga primerno prilagodimo za uporabo na vodovodnih sistemih. Na konču opisano metodo uporabimo na konkretnem primeru in komentiramo rezultate.
2 Hidravlično ozadje
2.1 Hidravlično modeliranje cevovodnih sistemov
Hidravlični modeli so pomembno orodje v pro-česu odločanja, ki omogočajo numerične simu-lačije delovanja realnega vodovodnega sistema. Na podlagi matematičnih formulačij fizikalnih zakonitosti toka vode je mogoče s podatki o fizikalnih karakteristikah sistema, obteznih primerih in robnih pogojih simulirati delovanje vodovodnih sistemov, ki odraza stanje sistema v naravi. Za vzpostavitev hidravličnega modela se zahteva zajem fizikalnih in nefizikalnih lastnosti posameznih gradnikov vodovodnega sistema.
Matematične formulačije hidravličnih modelov so razdeljene v stiri skupine, ki se delijo po zajemu časovne spremenljivosti pretoka: stalni tok, razsirjena časovna simulačija (kvazi) stalnega toka, nestalni nestisljivi tok in nestalni stisljivi tok oziroma analiza vodnega udara (Boulos et al., 2006). Najpogosteje uporabljeni simulačiji za analizo obseznih vodovodnih sistemov sta statična simulačija stalnega toka in razsirjena časovna simulačija stalnega toka. V
nadaljevanju bodo podani splošni izrazi za zapis zakona o ohranitvi mase za točkovne gradnike in zakona o ohranitvi energije za linijske gradnike, ki jih obravnavamo kot "dolge cevovode". Navedene enačbe tvorijo sistem nelinearnih enačb, katerega resitev predstavlja tudi resitev hidravličnega modela vodovodnega sistema.
Zakon o ohranitvi mase pove, da je za vsako vozlisče v dovedeni oziroma odvzeti pretok qv enak razliki vsote pretokov Qe v ceveh e, kjer tekočina v vozlisče v vteka, in vsote pretokov Qe' v čeveh e', kjer tekočina iz vozlisča v izteka (Boulos et al., 2006),
X/Qe - X/Qe' = •
(1)
Zakon o ohranitvi energije pravi, da je razlika med energijama na odseku med dvema vo-zlisčema v in v' enaka sestevku energijskih izgub v zaporednih čevovodih e\,...,ei obravnavanega odseka med v in v' (Boulos et al., 2006):
Hv — Hv' = 2J KeiQei |Qe-i|
i= 1
m— 1
(2)
V enačbi (2) člen Kei zdruzuje koefičiente, ki predstavljajo karakteristike posameznega čevo-voda za izračun energijskih izgub, eksponent m pa določimo glede na uporabljeni model za izračun energijskih izgub (npr. za Darčy-Weisbačhov model je m = 2). Karakteristika čevovoda Kei je v Darčy-Weisbačhovem modelu linijskih energijskih izgub podana z izrazom:
Ke = —
1 fiLi
2g diSf
(3)
kjer je fi brezdimenzijski koefičient trenja [—], Li dolzina čevovoda [m], d premer čevovoda [m], g gravitačijska konstanta [m/s2] in Si prerez čevovoda [m2], i = 1,..., l.
Za vzporedne čevovode, ki tvorijo zanko, je sestevek energijskih izgub čelotnega odseka zanke enak nič (Boulos et al., 2006):
Hv - Hv' =J2 KeiQei |Qe-i|m-1 = 0 (4) i=1
Nabor enačb (1)-(4) za vsa vozlisča in vse odseke čevovodov je osnova za določitev nelinearnega sistema hidravličnega modela čevovo-dnega sistema, katerega resitev daje informa-čije o neznanih pretokih Qe v čevovodih in neznanih količinah energije Hv v vozlisčih. Skupno stevilo neznank v sistemu je vsota vseh neznanih pretokov in neznanih energij v vo-zlisčih. Za resevanje tega nelinearnega sistema enačb je bilo razvitih več metod. Najbolj raz sirjena med njimi je gradientna metoda, ki stajo leta 1987 predlagala Todini in Pilati (Todini & Pilati, 1987) in je implementirana v hidravlični simulačijski model EPANET 2 (Ros-sman, 2000).
2.2 Todinijev odpornostni indeks Ir
Formulačijo odpornostnega indeksa Ir je podal Todini (Todini, 2000). Odpornostni indeks Ir prikaze, kolik sno stopnjo varnosti oskrbe s pitno vodo je mogoče se zagotoviti ob nastopu motenj (npr. pozar ali lom čevi) v čelotnem vodovodnem sistemu. Zaradi povečanega odjema oziroma iztekanja iz čevovoda se povečajo pretoki v čeveh in zmanj sa tlak v vodovodnem omrezju. Vrednost indeksa Ir določa stopnjo varnosti oskrbe. Pove nam torej, ali je na odje-mnih mestih na razpolago dovolj tlaka za zagotovitev nemotene oskrbe s pitno vodo. Brezdimenzijski indeks Ir je definiran s pomočjo razmerij različnih moči. Todini navaja tri različne primere moči: celotno moč Pa, porabljeno moč PD in razpoločljivo moč PN. Enačba (5) podaja zvezo med njimi:
Pa = PD + PN •
(5)
Celotna moč Pa (6) je določena z energijskim potencialom na vodnem viru oziroma na vtoku v vodovodni sistem. Definirana je s produktom specifične teZe vode y [N/m3] in vsote produkta izdatnosti vodnih virov Q? [l/s] ter hidravlične visine H? [m] (angl. head) vseh ns vodnih virov:
preostalo moč. Minimalna hidravlična visina je v tem primeru sestevek nadmorske visine vozlisča in minimalno potrebnega tlaka v tem vozlisču za nemoteno oskrbo s pitno vodo. Z razliko čelotne moči Pa in minimalne preostale moči v vozlisčih dobimo maksimalno vrednost porabljene moči:
pa = Y £ QHS.
(6)
i=1
Pa = Y £ QS H + £ Pp.
i=1
i=1
Pd = Y £ QfHf.
(8)
i=1
Razliko čelotne moči Pa in porabljene moči PD predstavlja razpolozljiva moč PN (9), ki ostane na razpolago posameznim vozlisčem. Izračunamo jo kot produkt spečifične teze vode Y [N/m3] in vsote produkta porabe vode qjV [l/s] ter hidravlične visine HV [m] v vseh nv vozli sčih sistema:
P,
N
nv
Y £ HV.
i=1
PDmax = Pa - Y £ H
(10)
i=1
Ker vodovodni sistemi vključujejo tudi črpalke, to dodatno moč Pp vključimo v enačbo (6) ter jo s tem preoblikujemo v:
Todinijev odpornostni indeks Ir končno izračunamo po enačbi:
Ir = 1 -
Pd
(7)
P Dm
(11)
Moč PD (8) predstavlja porabljeno moč v vodovodnem sistemu zaradi linijskih energijskih izgub, ki nastanejo zaradi trenja vodnega toka z ostenjem čevovoda. Skupno porabljeno moč PD dobimo tako, da za vsako čev posebej določimo linijske energijske izgube Hte [m], nato sestejemo produkte linijskih izgub in pretokov Qe [l/s] v vseh ne čeveh v vodovodnem omrezju ter pomnozimo s spečifično gostoto vode 7 [N/m3]:
Vrednost odpornostnega indeksa Ir je vedno manj sa od 1. Večja kot je vrednost odpornostnega indeksa Ir, večja je zanesljivost oskrbe s pitno vodo končnih porabnikov, ki se oskrbujejo iz tega vodovodnega sistema (Todini, 2000).
3 Metoda
3.1 Algoritem spektralnega razbitja grafov
Denimo, da imamo povezan graf G = (V, E), kjer je V mnoziča n točk (vozli sč) grafa in E mnoziča povezav med temi n točkami, za nek n G N. Graf zelimo razdeliti na k < n podgra-fov Gi, G2,..., Gk tako, da bo vsaka točka iz originalnega grafa pripadala natančno enemu izmed dobljenih podgrafov. Ce z Vi označimo mnozičo točk podgrafa Gj, torej velja:
V = V1 U V2 u ■ ■ ■ U Vk, kjer
vi n v j =
i =j.
(12)
(9)
Ce v vsakem vozliSču upoštevamo le minimalno hidravlično višino H^ [m], dobimo minimalno
Pri tem je vsak Gj = (Vj, Ej), i = 1,..., k, indučiran podgraf grafa G, kar pomeni, da prevzame vse povezave grafa G med točkami znotraj Vj. Povezava originalnega grafa, ki ima
kraji sci v različnih podmnožicah iz (12), ni vsebovana v nobenem od dobljenih podgrafov G i. Opisano razbitje grafa označimo s P:
P := {Gi,G2,...,Gk}.
(13)
Graf G je utezen: cvv/ > 0 označuje utez na povezavi vv' £ E. Te utezi se prenesejo na povezave v induciranih podgrafih Gi. Graf z utezmi podamo z utezeno matriko sosednosti A = (Aij) z elementi:
A	I Cvivj,
j = 0,
ViVj G E, sicer,
(14)
i, j = 1,...,n.
Ker povezave niso usmerjene, velja cvv/ = cv/v, torej je matrika A simetrična. Privzeli bomo, da so vse utezi na povezavah pozitivne in normirane tako, da je matrika A tudi (dvojno) stohastična:
vev
cvv/ = 1 za vse v' e V.
C (P ):= E E
Cvv' .
(16)
i=j vv' eE veVi,v' eVj
i =
S tem zagotovimo, da posamezni podgrafi niso preveliki. Ker zelimo tudi, da imajo vsi pod-grafi nasega razbitja priblizno enako stevilo točk, bomo iskali razbitje, kjer je \Vj\ & £ za vse i = 1,...,k.
Imamo torej dve zahtevi, ki jima skusamo zadostiti:
•	minimalna vrednost funkčije C (P),
•	stevilo točk v vsakem podgrafu razbitja je kar se da enako in navzgor omejeno z £.
V kolikor teh dveh zahtev ni mogoče izpolniti, zelimo poiskati enolično re sitev, ki se jima najbolj pribliza. Opisimo postopek, po katerem bomo to re sitev dobili.
Vsakemu razbitju mnoziče točk (12) lahko priredimo matriko M := (Mj) velikosti n x k, katere elementi so:
(15)
Mij : =
1, Vi G Vj, 0, sicer,
(18)
Ce je vrednost utezi na povezavi enaka 0, je to enako, kot da ta povezava v grafu ni prisotna.
Razbitju P priredimo vrednostno funkcijo C(P), ki je definirana z naslednjo enačbo:
Funkčija C (P) sesteje utezi na povezavah, ki v G povezujejo različne podgrafe iz razbitja P. Nas čilj je poiskati razbitje, za katero je vrednost te funkčije čim manj sa.
Naslednji pogoj, ki ga za razbitje zahtevamo, je maksimalno dovoljeno stevilo točk £ v posameznem podgrafu. Stevilo £ izračunamo kot navzgor zaokrozeni čeli del kvočienta med n, stevilom vseh vozli sč v grafu G, in k, stevilom podgrafov v razbitju P:
i = 1,..., n, j = 1,..., k. Matriko M imenujemo matrika k-razbitja. Opazimo, da velja tudi obratno: vsaka matrika k-razbitja poda neko razbitje točk grafa oblike (12).
Za matriko k-razbitja M lahko enostavno preverimo naslednje lastnosti:
(i)	V vsaki vrstiči matrike M je samo en element neničeln (to je, enak 1). To pomeni, da vsaka točka grafa pripada natančno enemu podgrafu v razbitju.
(ii)	Vsota elementov (oziroma stevilo neničelnih elementov) v j-tem stolpču matrike M je enaka stevilu točk v j-tem podgrafu razbitja, to je v mnoziči Vj.
(iii)	Ml k = 1 n, kjer z tm označimo vektor z m kooordinatami, ki so vse enake 1.
(17) (iv) MTM = diag(|Vi|, |— iVk|).
(v) Za Frobeniusovo normo matrike M velja
\
n k
EE M
i= i j=i
||M||f :=
= ^/sl(M TM) = Vn
(tu sl označuje sled matrike, tj. vsoto diagonalnih elementov (Meyer, 2000)).
Z upostevanjem teh lastnosti lahko vrednostno funkcijo C (P) iz (16) zapisemo s pomočjo matrike sosednosti A in matrike k-razbitja M kot:
C (P) = E(M TAM )ij
i=j
= 1 TmTAM1 k - E(MTAM)ii i
= 1TAln - sl(MTAM).	(19)
Torej, če zelimo, da je vrednost funkcije C (P) čim manjsa, mora biti sled produkta matrik MTAM čim večja. Poleg tega zahtevamo se
MTM < £Ik,
Pri vsakem k-razbitju potrebujemo le k največjih lastnih vrednosti matrike A, tako da se velikosti matrik D in U spreminjata s Številom k. Vedno velja:
U TU = Ik in AU = UD, (22)
kjer smo z Ik označili enotsko matriko velikosti k x k.
Denimo, da znamo poiskati tako k x k matriko Z, da je: M = UZ. Iz enačb (22) in (23) sledi:
M tM = Z TU TUZ = ZT Z
in:
sl(M 1 AM) = sl(Z 1 U 1 AUZ) = sl(Z 1 DZ).
Vrednost tega izraza je največja, če je D = Ik, zahteva (20) pa se prepise v:
Z TZ < nk.
(23)
(20)
da zadostimo tudi pogoju o čimbolj enakomernem razbitju.
Ker so vse lastne vrednosti simetrične sto-hastične matrike A realne in je največja med njimi enaka 1 (Meyer, 2000), velja
sl(MTAM) < sl(MTM) = n. (21)
Označimo k največjih lastnih vrednosti A:
Ai = 1 > A2 > ■ > Ak,
in njim pripadajoče ortonormirane lastne vektorje: u1 ,u2,...,uk. Te lastne vrednosti za-pi semo po diagonali v diagonalno matriko D velikosti k x k, lastne vektorje pa po stolpčih v matriko U velikosti n x k:
D := diag (A1, A2,..., Ak) in
U := [ui U2 ■ ■ ■ Uk] .
To seveda ni vedno dosegljivo. Hespanha (Hespanha, 2004) pokaze, da se tema dvema pogojema zadovoljivo priblizamo, če matriko k-razbitja M poi sčemo z optimizačijsko enačbo:
M = argmin ||M' - UZ||F, (24)
M'
kjer M' predstavlja vse mozne matrike k-razbitja (tj., preverimo vsa mozna k-razbitja grafa) in je Z taka k x k matrika, da zanjo (priblizno) velja (23) in je UZ dovolj blizu neki matriki k-razbitja.
Za poznano matriko Z dobimo matriko M iz enačbe (24) z metodo projekcije. S pomočjo lastnosti matrik k-razbitja in linearnosti operatorja sledi namreč dobimo
||M' - UZ||F = sl ((M' - UZ)t(M' - UZ)) = n + ||Z||f - 2 ■ sl(M'TUZ).
Posebej, če označimo M := UZ, je
k n
sl (M'tm) = EE Mj Mj j=1i=i
n
= E ^' i= 1
kjer smo upoštevali lastnost (i) matrike k-razbitja M' in z ji označili njen edini neničelni element v i-ti vrstici (ki je enak 1). Vrednost izraza ||M' — UZ||f bo torej najmanjsa, če v matriki M' enice postavimo na mesta (i, ji), kjer indekse določimo tako, da je Mj največji element v i-ti vrstiči matrike M = UZ.
Zdaj moramo povedati le se, kako poiskati ustrezno matriko Z. Denimo, da lahko n vr-stič n x k matrike U zdruzimo v skupine oziroma kopiče (angl. člusters) okoli k ortonor-miranih vektorjev z1, z2,..., zk. To pomeni, da je vsaka od vrstič blizu natanko enemu od vektorjev z j. Razdaljo definiramo s skalarnim produktom tako, da je za vrstičo v j-ti skupini skalarni produkt z zj najblizje 1 in z ostalimi vektorji zi,l = j, blizu 0. Ce definiramo matriko
Z := [zi Z2 ■ ■ ■ Zk] ,
bo i-ta vrstiča matrike MMf = UZ, ki je sestavljena iz skalarnih produktov i-te vrstiče matrike U z vektorji z j, blizu nekemu standardnemu enotskemu vektorju v Rk. To pomeni, da je matrika M = UZ blizu neki matriki k-razbitja.
Stolpče matrike Z torej lahko poisčemo z znano metodo kopičenja z voditelji (angl. k-means člustering), ki je implementirana v različna programska orodja, mi smo uporabili MATLAB 2015b (MATLAB, 2015). Algoritem lahko prilagodimo tako, da so velikosti razredov navzgor omejene z t (Hespanha, 2004; Zhong & Ghosh, 2003).
Povzemimo zdaj osnovne korake algoritma spektralnega razbitja:
1.	poisči k največjih lastnih vrednosti matrike A in pripadajoče ortonormirane lastne vektorje, iz katerih sestavi matriko U;
2.	z metodo kopičenja z voditelji s pomočjo matrike U določi ustrezno matriko Z;
3.	s projekčijo na mesta največjih elementov po vrstičah matrike U Z zapisi matriko k-razbitja M.
Opisani algoritem za izvajanje potrebuje dva vhodna podatka:
•	utezeno matriko sosednosti A,
•	stevilo podgrafov v razbitju: k.
Matrika sosednosti A je nespremenljiv vhodni podatek, ki je definiran s topologijo vodovodnega omrezja. Število podgrafov k je spremenljiv vhodni podatek, ki potem v končni fazi predstavlja s tevilo merilnih območij.
3.2 Implementacija algoritma spektralnega razbitja grafov v algoritem za zasnovo merilnih območij
Algoritem za zasnovo merilnih območij smo implementirali v programskem orodju MA-TLAB 2015b (MATLAB, 2015). Poleg algoritma za spektralno razbitje grafov smo v algoritmu uporabili tudi hidravlični simulačijski program EPANET 2.0 (Rossman, 2000). Izvedba hidravličnih simulačij znotraj programskega okolja MATLAB je mogoča s programskim orodjem EPANET-MATLAB (Eliades et al., 2016). To programsko orodje uporablja pristop objektnega programiranja z določitvijo razredov v MATLAB-u, ki zagotavlja standardiziran način upravljanja podatkovne strukture hidravličnega modela. Programsko orodje EPANET-MATLAB ponuja obsezen nabor funkčij, ki uporabniku omogočajo pridobivanje in spremljanje podatkov hidravličnega modela ter simulačije hidravličnih stanj in dinamike kakovosti vode v vodovodnem sistemu.
Slika 2: Shema algoritma zasnove merilnih območij. Figure 2: Algorithm of the design of district metering areas.
Prvi korak algoritma za zasnovo merilnih območjih je branje vhodne datoteke hidravličnega modela v tekstovnem formatu datoteke *.inp, ki je običajno ustvarjena v samostojnem programskem okolju EPANET 2.0. V tem izhodisčem koraku se ustvari objekt in ini-čializirajo spremenljivke objekta glede na podatke iz EPANET-ove vhodne datoteke *.inp. V naslednjem koraku se izvede hidravlična simulačija izhodisčne topologije vodovodnega sistema, ki najpogosteje predstavlja obstoječe ali predvideno stanje obravnavanega vodovodnega sistema. Rezultati hidravlične simula-čije, kot npr. tlaki v vozli sčih p, pretoki po če-veh Q in vozli sčna poraba q, se shranijo in uporabijo v nadaljevanju algoritma. Kot prvo so uporabljeni za izračun odpornostnega indeksa Ir (11).
Za uporabo algoritma spektralnega razbitja grafov je treba oblikovati utezeno matriko so-sednosti A (14) danega grafa, ki ustreza topologiji obravnavanega vodovodnega omrezja. Utezi na povezavah cViVj morajo odrazati ustrezne hidravlične lastnosti vodovodnega sistema. Paziti moramo tudi, da so utezi pozitivne in normirane tako, da je matrika A tudi (dvojno) stohastična (15).
Poleg utezene matrike sosednosti A algoritem spektralnega razbitja potrebuje se podatek o stevilu podgrafov v razbitju k, ki ustreza stevilu merilnih območij. Za iskanje optimalnega razbitja zelimo preveriti resitve pri različnih vrednostih k, zato je treba določiti spodnjo in zgornjo mejo, kmin in kmax. Vedno formiramo čimbolj enakomerna merilna območja, za kar poskrbi pogoj o velikosti podgrafov (20). Z večanjem s tevila merilnih območij se velikost posameznega merilnega območja manj sa. Zato stevilo k spreminjamo, dokler so resitve se smiselne oziroma se formirajo dovolj velika merilna območja.
Algoritem spektralnega razbitja grafov nam pri izbranem stevilu k vrne k podgrafov in označi pripadnost vsakega vozlisča posame-
znemu podgrafu. Ker algoritem ne določi tudi ustrezne pripadnosti povezav, moramo to storiti naknadno. Z analizo pripadnosti za vse pare vozli sč Vi,Vj £ V preverimo, ali pripadajo istemu podgrafu Gi v razbitju (13). Ce je to izpolnjeno, je v podgrafu Gi posledično vsebovana tudi povezava Vj v j. V nasprotnem primeru pa imamo mejno povezavo oziroma mejno cev. Identifikacija mejnih povezav med posameznimi podgrafi Gi je pri oblikovanju merilnih območij ključna. Med sosednjimi podgrafi zelimo največ eno mejno povezavo, na kateri se spremljajo hidravlične količine, tj. pretok Q, ostale mejne povezave pa naj bi se zaprle oziroma ukinile.
Iz nabora vseh mejnih povezav m zelimo ohraniti odprtih k — 1 povezav, kar zagotavlja povezanost čelotnega grafa in prisotnost minimalnega vpetega drevesa (Wilson & Watkins, 1997). Vzpostavimo mnozičo vseh moznih kombinačij odprtih povezav C med k podgrafi, ki jih je v tem primeru (fcmJ. Z uporabo funk-čij programskega orodja EPANET-MATLAB je mogoče čevi, ki so v posamezni kombina-čiji C £ C predvidene za zaprtje, tudi v hidravličnem modelu zapreti in s tem zasnovati novo topologijo vodovodnega sistema. Po posodobitvi topologije se izvede hidravlični preračun, s katerim preverjamo ustreznost zasnov merilnih območij glede na kriterij minimalnega preostalega tlaka v vodovodnem sistemu. C e so izračunani tlaki pri izbrani kombinačiji odprtih mejnih povezav nizji od minimalno zahtevanih, to kombinačijo odstranimo.
V primeru ustreznih tlačnih razmer pri izračunu hidravlične simulačije za dano kom-binačijo odprtih čevi izračunamo odpornostni indeks Ir. Nato vse kombinačije mejnih povezav C £ C, ki ustrezajo hidravlični preverbi minimalnih tlakov, med seboj primerjamo glede na odpornostni indeks Ir. Za posamezni k kot najustreznej so zasnovo merilnih območij izberemo tisto resitev, ki izkazuje najvi sjo vrednost odpornostnega indeksa Ir.
Algoritem zasnove merilnih območij nato preveri rešitve, ki ustrezajo k + 1 merilnim območjem, in zaključi z izvajanjem, ko je k = kmax. Slika 2 prikazuje shematski potek celotnega algoritma zasnove merilnih območij z vključenim algoritmom spektralnega razbitja grafov.
4 Delovni primer
Algoritem za zasnovo merilnih območij je bil preizkusen na podsistemu vodovodnega sistema Kranj. Ozje področje obravnave je območje Primskovega in severnega dela Planine. Vzpostavljeni hidravlični model je sestavljen iz 476 vozlisč in 500 čevi. Območje, na katerem smo uporabili algoritem zasnove merilnih območij, je na sliki 3 označeno s črtkano črto.
Obravnavani del vodovodnega sistema Kranj ni hidravlično ločen od ostalega vodovodnega sistema, zato je bilo v hidravlično analizo zajeto sirse območje vodovodnega sistema, ki pa ni neposredno vključeno v analizo določevanja merilnih območij. Vključitev sirsega območja vodovodnega sistema je bila potrebna zaradi robnih pogojev hidravličnega izračuna, saj spreminjanje topologije čevovodnega omrezja lahko bistveno vpliva na robne pogoje hidravličnega modela. Z zajemom sirsega območja so bili vplivi na robne pogoje bistveno zmanjsani, kar daje večje zaupanje v hidravlične simulačije, izračune odpornostnih indeksov Ir in pridobljene zasnove merilnih območij.
Določitev utezene matrike sosednosti A zahteva oblikovanje utezi, ki naj čim bolje vključuje čilje in vodila pri zasnovi merilnih območij. V ta namen sta bila izbrana dva parametra vodovodnega sistema, to sta notranji premer čevi d in hidravlična visina H v vo-zlisčih. Oba parametra v utezi povezave posredno predstavljata mero hidravlične prevodno-
Slika 3: Delovni primer Planine (smeri vtokov in iztokov so označene s puščicami).
Figure 3: Case study Planina (directions of inflows and outflows are denoted by arrows).
sti vodovodnega sistema in zajemata vplive, ki so zajeti v izračunu odpornostnega indeksa Ir. Medtem ko je hidravlična visina neposredni parameter pri izračunu odpornostnega indeksa, je premer čevovoda d nadomestni parameter za pretok Q po čevi. Ker ima vsaka povezava dve vozlisči, smo za izračun utezi izbrali nizjo vrednost hidravlične visine, saj se tako ohranja in-formačija o usmeritvi pretokov v vodovodnem sistemu. Oba parametra sta normirana in v enakem delezu prispevata k vrednosti posameznih utezi na povezavah grafa G.
Dodatno k izbiri utezi je v algoritmu za zasnovo merilnih območij treba nastaviti se nekatere vhodne podatke. Pri algoritmu za spek-
tralno razbitje grafov so bile analizirane zasnove za različno S tevilo merilnih območij. Za minimalno Stevilo merilnih območij smo vzeli kmin = 2, za maksimalno pa kmax = 8. Hidravlično preverjanje zasnov merilnih območij smo opravili z mejo minimalnega preostalega tlaka, ki je znasala 2,5 bar oziroma 25 m vodnega stolpča. Obtezbeni primer, ki smo ga uporabili za preverjanje hidravlične ustreznosti različnih zasnov merilnih območij, je bila srednja poraba v vodovodnem sistemu. Pri tem smo uporabili hidravlične simulačije stalnega toka.
5 Rezultati
Analiza rezultatov algoritma zasnove merilnih območij z uporabo algoritma spektralnega razbitja grafov je pokazala, da je med vrednostmi od kmin do kmax najustreznejsa res itev pri k = 4. Ključna mera, na podlagi katere smo izbrali najustreznej so re sitev, je vrednost odpor-nostnega indeksa Ir. Slika 4 prikazuje maksimalne vrednosti odpornostnega indeksa Ir gleda na k podgrafov.
Odpornostni indeks Iri za obstoječe stanje je enak 0,684 (glej pregledničo 1). Po vzpostavitvi merilnih območij se vrednost Ir zniza na 0, 626 pri k = 5. Kot najbolj sa je bila izbrana resitev pri k = 4, kjer je odpornostni indeks dosegel vrednost 0, 646. Iz slike 4 je razvidno, da pride pri k = 4 do povisanja vrednosti odpornostnega indeksa Ir, kar nakazuje na dobro ravnovesje med izbranimi vo-zlisči za posamezne podgrafe Gi, ki določajo merilna območja, in mejnimi povezavami oziroma čevmi, ki omogočajo boljso pretočnost in s tem boljse hidravlične razmere tako zasnovanega vodovodnega sistema.
Primerjava podatkov o stevilu mejnih povezav in stevilu merilnih območij je prikazana na sliki 5. Zaradi uporabe končepta minimalnega vpetega drevesa je stevilo pretočnih oziroma
Slika 4: Odpornostni indeks Ir v odvisnosti od števila merilnih območij (Gorjup, 2016).
Figure 4: Resilience index Ir depending on the number of DMAs (Gorjup, 2016).
odprtih mejnih povezav enako vrednosti k — 1 in s tem linearno narasčajoče s stevilom merilnih območij k. Ob pogledu na podatke o stevilu vseh mejnih povezav je razvidno, da se te povečujejo do k = 5, nato pa se skupno stevilo mejnih povezav ustali. To je posle-diča definičije vrednostne funkčije C (P) (16), ki lahko vpliva tudi na stevilo mejnih povezav. Razlika med stevilom mejnih povezav in stevilom odprtih čevi, na katerih se merijo hidravlične količine, predstavlja stevilo zaprtih mejnih čevi (slika 5). Tudi v tem pogledu je razvidno, da se stevilo zaprtih povezav povečuje do k = 5, pri čemer nato stevilo upada.
Interpretačije podatkov iz slike 5 imajo velik vpliv na sprejemanje odločitev o zasnovi merilnih območij na vodovodnem sistemu. Razvidno je, da z večanjem stevila merilnih območij
16
0	2	4	6	8
Št. merilnih območij
■	Št. odprtih mejnih povezav
■	Št. zaprtih mejnih povezav p Št. vseh mejnih povezav
Slika 5: Mejne povezave pri različnih vrednostih k (Gorjup, 2016).
Figure 5: Boundary links for different k values (Gorjup, 2016).
prihaja do razčlenjenosti omrežja in s tem tudi večjega Števila ukrepov, ki jih je treba izvesti na cevovodnem omrežju. Zagotavljanje povezanosti celotnega vodovodnega sistema in uporaba koncepta minimalnega vpetega drevesa omogoča generiranje resitev, ki imajo minimalno stevilo odprtih povezav med merilnimi območji in s tem minimalno stevilo merilnih jaskov. Zagotovitev merilnih jaskov lahko predstavlja finančno intenziven ukrep (Farley, 2001), (De Paola et al., 2013). Ravno tako pa tudi stevilo zaprtih mejnih povezav zahteva ukrepe na čevovodnem omrezju, pri čemer pa ti nimajo tako močnih finančnih posledič. Pri sprejemanju odločitev o stevilu merilnih območij je tako gotovo potrebno večkriterijsko odločanje, kjer se upostevajo različni in včasih nasprotujoči si tehnični, finančni in prostorski vidiki.
Pregledniča 1: Povzetek ključnih vrednosti za k = 4.
Table 1: Summary of results for k = 4.
Količina	Vrednost
k	4
Iri	0,684
Ir4	0,646
St. resitev	20 (120)
St. mejnih čevi	10
St. zaprtih čevi	7
St. merilnih mest	3
Na podlagi analiziranih rezultatov o odpor-nostnih indeksih Ir in mejnih povezavah pri različnem Številu merilnih območij, smo reSitev s k = 4 merilnimi območji izbrali za najustreznejšo z vidika zagotavljanja hidravličnih sposobnosti in z vidika zahtevanih ukrepov vzpostavitve merilnih območij.
Pregledniča 1 in slika 6 podajata podatkovni in grafični prikaz resitve pri k = 4. Pri k = 4 je bilo vseh kombinačij mejnih povezav 120, a so nekatere resitve izkazovale neustrezno hidravlično zasnovo. Stevilo resitev, pri katerih so se pojavili negativni tlaki, je bilo kar 57 od skupnih 120. Dodatno je bilo iz preostalega nabora kombinačij 43 resitev taksnih, ki so imele vsaj v enem vozlisču tlak manjsi od meje minimalnega preostalega tlaka, ki je znasala 25 m vodnega stolpča. Tako je bilo le 20 resitev pri k = 4 hidravlično ustreznih.
Resitev za k = 4 pri maksimalnem odporno-stnem indeksu Ir je prikazana na sliki 6, medtem ko je povzetek lastnosti ustvarjenih stirih merilnih območij podan v pregledniči 2.
6 Zaključek
V prispevku obravnavamo zasnovo merilnih območij vodovodnega sistema, za katero smo uporabili algoritem spektralnega razbitja gra-
Slika 6: Zasnova merilnih območij pri k = 4. Figure 6: Designed DMA at k = 4.
fov. Za oblikovanje utežene matrike sosedno-sti smo upoštevali dva normirana hidravlična kriterija, notranji premer cevi d in hidravlično višino vozlišč H, ki sta v enakih deležih prispevala k oblikovanju uteži na povezavah celotnega grafa. Spektralno razbitje grafa smo implementirali v programsko okolje MATLAB, skupaj ž evalvacijskim algoritmom hidravličnih razmer v povezavi s kombinacijami odprtosti mejnih povezav med razlicnimi merilnimi obmocji.
Predstavljena metoda zasnove merilnih mest omogoca ucinkovito in uspesno razbitje na podgrafe, ki ponazarjajo merilna obmocja v vodovodnem sistemu. Implementirani algoritem s spektralnim razbitjem na grafe smo preizkusili na realnem delovnem primeru.
Med analiziranimi zasnovami za k merilnih obmocij so bile izbrane tiste, ki so izkazovale najvisje vrednosti odpornostnega indeksa Ir. Rezultati kazejo na obstoj resitev, ki imajo glede na hidravlicne karakteristike in obseg ukrepov za vzpostavitev merilnih obmocij prednosti pred drugimi zasnovami. Podrobna
Preglednica 2: Osnovne značlnosti merilnih območij pri k = 4 (Gorjup, 2016).
Table 2: Main properties of DMAs at k = 4 (Gorjup, 2016).
	St. vozlisc	St. cevi	Skupna poraba	Skupna dolžina cevi	Povprecna nadmorska visina
	[-]	[-]	[l/s]	[m]	[m]
DMA1	92	94	7,55	1743	383,63
DMA2	83	85	1,51	3410	384,09
DMA3	107	108	5,53	3282	384,00
DMA4	170	179	14,58	4662	387,36
hidravlična in topoloska analiza kaže na dejstvo, da rešitvam, ki uporabljajo enake transportne poti vode, kot so prisotne v obstoječem vodovodnem sistemu, ustrezajo visje vrednosti odpornostnega indeksa. Na primeru je bilo ugotovljeno, da imajo te resitve topoloske lastnosti, ki so blizu lastnostim obstoječega celovitega vodovodnega sistema.
Viri
Alvisi, S., Franchini, M. (2014). Heuristic procedure for the automatic creation of district metered areas in water distribution systems, Urban Water Journal 11 (2), 137-159.
Boulos, P. F., Lansey, K. E., Karney, B. W. (2006). Comprehensive water distribution systems analysis handbook for engineers and planners. MWH Soft, Inc., Pasadena, California, USA.
De Paola, F., Fontana, N., Galdiero, E., Giu-gni, M., Sorgenti degli Uberti, G., Vitaletti, M. (2013). Optimal design of district me-
tered areas in water distribution networks, "12th International Conference on Computing and Control for the Water Industry, CCWI2013", Procedia Engineering 70, 449457.
Di Nardo, A., Di Natale, M., Guida, M., Mu-smarra, D. (2013). Water Network Protection from Intentional Contamination by Sec-torization, Water Resources Management 27 (6), 1837-1850.
Diao, K., Zhou, Y., Rauch, W. (2013). Automated creation of district metered areas boundaries in water distribution systems, Journal of Water Resources Planning and Management 139 (2), 184-190.
Eliades, D. G., Kyriakou, M., Stelios Vra-chimis, S., and Polycarpou, M. M. (2016). EPANET-MATLAB Toolkit: An Open-Source Software for Interfacing EPANET with MATLAB, "14th International Conference on Computing and Control for the Water Industry, CCWI 2016", Amsterdam, Netherlands, str. 8.
Farley, M. (2001). Leakage management and control. A best practice training manual, World Health Organization, Geneva.
Ferrari, G., Savic, D., and Becciu, G. (2014). Graph-Theoretic Approach and Sound Engineering Principles for Design of District Metered Areas,J. Water Resources Planning and Management 140 (12), 040140361-04014036-13.
Galdiero, E. (2015). Multi-Objective Design of District Metered Areas in Water Distribution Networks, PhD Thesis, Department of Civil, Architectural and Environmental Engineering, University of Naples Federico II.
Gorjup, M. (2016). Uporaba teorije grafov za zasnovo merilnih obmocij v vodovodnih
omrezjih. Magistrsko delo, Univerza v Ljubljani, Fakulteta za gradbenistvo in geodezijo, 68 str.
Gomes, R., Sa Marques, A., Sousa, J. (2012). Identification of the optimal entry points at District Metered Areas and implementation of pressure management, Urban Water 9 (6), 365-384.
Hajebi, S., Temate, S., Barrett, A., Clarke, A., Clarke, S. (2014). Distribution Network Sec-torization Using Structure Partitioning and Multi-Objective Optimization, "16th Water Distribution System Analysis Conference, WDSA2014 - Urban Water Hydroinforma-tics and Strategic Planning", Procedia Engineering 89, 1144-1151.
Herrera, M., Izquierdo, J., Perez-Garcia, R., Montalvo, I. (2012). Multi-agent adaptive boosting on semi-supervised water supply clusters, Advances in Engineering Software 50,131-136.
Hespanha, J. (2004). An efficient MATLAB Algorithm for Graph Partitioning, Technical Report, University of California.
MATLAB 2015b, The MathWorks, Inc., Na-tick, Massachusetts, United States.
Meyer, C. (2000). Matrix analysis and applied linear algebra, SIAM, Philadelphia.
Morrison, J. (2004). Managing leakage by District Metered Areas: A practical approach, Water 21 (2), 44-46.
Morrison, J., Rogers, D., Tooms, S. (2007). District Metered Areas Guidance Notes, Water Loss Task Force, IWA Publication.
Lambert, A (2003). Assessing non-revenue water and its components: A practical approach, Water 21, IWA, Issue 5.4, 50-51.
Rossman, L. A. (2000). Epanet2 Users Manual. "U. S. Environmental Protection Agency, National Risk Management Research Laboratory", Cincinnati, Ohio, USA.
Sempewo, J., Pathirana, A., Vairavamoorthy, K. (2009). Spatial analysis tool for development of leakage control zones from the analogy of distributed computing, "Proc. 10th Annual Water Distribution Systems Analysis Conference (WDSA2008)", Kruger National Park, South Africa, 1-15.
Todini, E., Pilati, S. (1987). A gradient method for the analysis of pipe networks, International Conference on Computer Applications for Water Supply and Distribution, Leicester Polytechnic, UK.
Todini, E. (2000). Looped water distribution networks design using a resilience index based heuristic approach, Urban Water 2, 115122.
Tzatchkov, V. G., Alcocer-Yamanaka, V. H., Bourguett-Ortíz, V.J. (2006). Graph theory based algorithms for water distribution network sectorization projects, "Proc. 8th Annual Water Distribution Systems Analysis Symposium, WDSA 2006", Cincinnati, Ohio, USA.
Wilson, R.J., Watkins, J.J. (1997). Uvod v teorijo grafov. DMFA Slovenije.
Zhong, S., Ghosh, J. (2003). Scalable balanced model-based clustering, "Proc.3rd SIAM Int. Conf. Data Minning", San Francisco, California, USA, 71-82.