1 UVOD
Spremljanje procesov skupaj z zaznavanjem napak in
diagnozo je v zadnjem času postalo eno popularnejših
področij procesne avtomatike. Za spremljanje procesa
in odkrivanje napak se pogosto uporabljajo tehnike na
podlagi modela procesa, ekspertni sistemi in sistemi
za razpoznavanje vzorcev [1], [2]. Hitri senzorji in
oprema za zajem podatkov nam omogočajo zajem ve-
like količine podatkov, ki nam kažejo, kaj se dogaja s
procesom. Za obdelavo in analizo zajetih podatkov je
bilo razvitih veliko statističnih metod [3], [4]. S pomočjo
Prejet 22. junij, 2011
Odobren 25. julij, 2011
meritev lahko ves čas spremljamo procese in zaznavamo
napake ter temu primerno prožimo alarme. V našem
primeru obravnavamo zaznavanje napake na senzorju
simuliranega procesa čiščenja odpadne vode [5], [6], [7].
Ti procesi so zaradi dnevnih, mesečnih in sezonskih vpli-
vov izpostavljeni nenehnemu spreminjanju dinamike.
Nanje vplivjo sprememba temperature, količina padavin
in različne obremenitve naprave. Teoretično modeliranje
je pri čiščenju odpadnih voda zelo zapleten postopek,
katerega rezultati so dvomljivi. Zato se za spremljanje
takih procesov uporabljajo statistične metode, ki teme-
ljijo na rudarjenju s podatki. Čiščenje odpadnih voda je
zelo nelinearen proces, zato je treba njegovo nelinearnost
upoštevati v modelu. Za predprocesiranje podatkov smo
v našem primeru uporabili algoritem za mehko rojenje
podatkov. S tem smo se izognili napačnim alarmom,
ki lahko nastanejo zaradi nelinearnosti procesa (če bi
uporabili linearni model za zaznavanje), saj z mehkim
modelom lahko poljubno natančno aproksimiramo neli-
nearnosti.
1.1 Mehki model z Gustafson-Kesselsovim rojenjem
V tem poglavju bomo opisali algoritme in metode,
ki so bile uporabljene za analizo podatkov. Razložen
bo Gustafson-Kesselov algoritem rojenja in izpeljan
ter identificiran Takagi-Sugeno mehki model procesa
čiščenja odpadnih voda.
1.1.1 Gustafson-Kesselov algoritem mehkega rojenja:
Gustafson-Kesselov algoritem rojenja nam omogoča,
da definiramo roje različnih oblik, kar nam pri takih
procesih, kot je obravnavani, pride prav.
ZAZNAVANJE NAPAK IN SPREMLJANJE ČIŠČENJA ODPADNIH VODA NA PODLAGI MEHKEGA MODELA 143
Matrika vhodnih podatkov je podana kot
X ∈ Rn×p. (1)
Vektor vhodnih podatkov v trenutku k je definiran kot
xk = [xk1, . . . , xkp] , xk ∈ Rp. (2)
Nabor n meritev označimo kot
X = {xk | k = 1, 2, . . . , n} (3)
in je predstavljen kot matrika n× p:
X =


x11 x12 . . . x1p
x21 x22 . . . x2p
...
...
...
...
xn1 xn2 . . . xnp

 . (4)
Cilj rojenja je razdeliti nabor podatkov X v c pod-
množic, ki jih imenujemo roji. Mehko razdeljen nabor
podatkov X je skupek mehkih podmnožic (rojev) {Ai |
1 ≤ i ≤ c}. Roji so definirani s pripadnostnimi funkci-
jami, ki so implicitno vsebovane v pripadnostni matriki
U = [µik] ∈ Rc×n. i-ta vrsta matrike U vsebuje vre-
dnosti pripadnostne funkcije i-tega roja Ai iz množice
podatkov X . Matrika U izpolnjuje naslednje pogoje:
pripadnostne stopnje so realna števila iz intervala od
nič do ena (µik ∈ [0, 1] , 1 ≤ i ≤ c, 1 ≤ k ≤ n;),
pripadnost vzorca xk k uniji rojev je ena (
∑c
i=1 µik =
1, 1 ≤ k ≤ n;), noben roj ni prazen oziroma ne vsebuje
vseh podatkov (0 <
∑n
k=1 µik < n, 1 ≤ i ≤ c). Zgoraj
našteto pomeni, da matrika U pripada množici, ki je
definirana kot:
M = {U ∈ Rc×n | µik ∈ [0, 1] ,∀i, k;
c∑
i=1
µik = 1,∀k; 0 <
n∑
k=1
µik < n, ∀i}.
(5)
V našem primeru je matrika pripadnostni (U ) dobljena
z mehkim ”c-means” algoritmom, ki temelji na Ma-
halanobisovi normi. Algoritem dobimo z minimizacijo
kriterijske funkcije ob upoštevanju omejitev iz En. 5:
J(X,U, V, λ) =
c∑
i=1
n∑
k=1
µmikd
2
ik + λ
c∑
i=1
n∑
k=1
(µik − 1) ,
(6)
kjer je U matrika pripadnostni množice podatkov X , V
je vektor centrov rojev
V = [v1, v2, . . . , vc] , vi ∈ Rp, (7)
d2ik pa je kvadrat razdalje
d2ik = (xk − vi)
T
Ai (xk − vi) .
Matrika Ai je definirana kot:
Ai = (ρidet (Ci))
1/p
C−1i ,
kjer je ρi = 1, i = 1, ..., c in p enak številu merjenih
spremenljivk. Ci je mehka kovariančna matrika i-tega
roja, definirana kot
Ci =
∑n
k=1 µ
m
ik (xk − vi) (xk − vi)
T∑n
k=1 µ
m
ik
.
Slednje nam omogoča detekcijo hiperelipsoidnega roja
v porazdelitvi podatkov. Če so podatki porazdeljeni ob
nelinearni hiperpovršini, algoritem najde roje, ki so lo-
kalni linearni aproksimatorji hiperpovršine. Prekrivanje
rojev je definirano s parametrom m ∈ [1,∞).
Za določanje števila rojev lahko uporabimo funk-
cije za merjenje veljavnosti rojev ali pa jih iterativno
dodajamo oziroma odstranjujemo glede na odstopanje
modela. Parameter prekrivanja m vpliva na mehkost
roja; od trdega (m = 1) do čisto mehkega (m→∞). V
našem primeru smo parameter nastavili m = 2, kar je
tudi običajna vrednost tega parametra.
1.1.2 Koraki Gustafson-Kesselovega algoritma: Po-
stopek Gustafson-Kesselovega algoritma za rojenje na-
bora podatkov X lahko opišemo z naslednjimi koraki:
• Inicializacija Izberemo število rojev c, faktor pre-
krivanja m (v našem primeru m = 2) in napako,
pri kateri se algoritem konča �end > 0 ( v našem
primeru �end = 0.001). Inicializiramo matriko
pripadnosti: U ∈M (naključno) in epoh r = 0.
• Zanka
r = r + 1
izračun centrov rojev:
v
(r)
i =
∑n
k=1
(
µ
(r)
ik
)m
xk∑n
k=1
(
µ
(r)
ik
)m , 1 ≤ i ≤ c. (8)
izračun kovariančne matrike in matrike Ai:
Ci =
∑n
k=1 µ
m
ik (xk − vi) (xk − vi)
T∑n
k=1 µ
m
ik
, (9)
Ai = (ρidet (Ci))
1/p
C−1i , 1 ≤ i ≤ c (10)
izračun razdalje vzorcev do centrov rojev
d2ik =
(
xk − v
(r)
i
)T
Ai
(
xk − v
(r)
i
)
,
1 ≤ i ≤ c, 1 ≤ k ≤ n
(11)
osvežitev matrike pripadnosti:
če je dik > 0, µ
(r)
ik =
1∑c
j=1
(
dik
djk
) 2
m−1
(12)
• dokler ni ||U (r) − U (r−1)|| < �end
1.1.3 Mehki model Takagi-Sugeno : Mehki TS-
modeli aproksimirajo nelinearen sistem z interpolacijo
lokalnih linearnih modelov. Vsak lokalni linearni model
prispeva svoj delež k izhodu globalnega modela. Če
imamo nabor vhodnih vektorjev
X = [x1, x2, . . . , xn]
T (13)
in nabor pripadajočih izhodov
Y = [y1, y2, . . . , yn]
T
, (14)
144 DOVŽAN, LOGAR, HVALA, ŠKRJANC
potem je mehki model podan v obliki pravil Ri:
Ri :
če je xk iz Ai potem ŷk = φi(xk), i = 1, . . . , c
(15)
Vektor xk je vektor vhodnih podatkov, ŷk je izhod
lokalnega linearnega modela v trenutku k. Vektor xk
pripada mehkim množicam (A1, . . . , Ac) z določeno
pripadnostjo µAi(xk) oziroma µik : R → [0, 1]. Funk-
cije φi(·) so lahko poljubne zvezne funkcije, čeprav so
največkrat uporabljene linearne afine funkcije.
Izhod globalnega modela izračunamo po enačbi:
ŷk =
∑c
i=1 µikφi(xk)∑c
i=1 µik
. (16)
En. (16) poenostavimo tako, da definiramo funkcijo:
βi(xk) =
µik∑c
i=1 µik
, i = 1, . . . , c (17)
Funkcija βi(xk) nam podaja normirano pripadnost
vzorca k posameznemu roju. S tem nam pove tudi,
kakšna je veljavnost pravila. Vsota funkcije po rojih
je ena (
∑c
i=1 βi(xk) = 1) ne glede na xk, če je
imenovalec βi(xk) različen od nič. To pa je lahko
zagotoviti s pravilno definicijo pripadnostnih funkcij nad
roji. Z združitvijo enačb (16) in (17) pridemo do slednje
poenostavljene enačbe za izhod globalnega modela:
ŷk =
c∑
i=1
βi(xk)φi(xk), k = 1, . . . , n (18)
Izhod lokalnega modela je običajno definiran kot line-
arna kombinacija elementov podatkovnega vektorja:
φi(xk) = xkθi, i = 1, . . . , c,
θTi =
[
θi1, . . . , θi(p+q)
]
. (19)
Vektor zmehčanih vhodov v trenutku k zapišemo kot:
ψk = [β1(xk)xk, . . . , βc(xk)xk] , k = 1, . . . , n, (20)
matriko zmehčanih podatkov pa kot:
ΨT =
[
ψT1 , ψ
T
2 , . . . , ψ
T
n
]
. (21)
Če zapišemo matriko koeficientov celotnega nabora pra-
vil kot
ΘT =
[
θT1 , ..., θ
T
c
]
, (22)
potem lahko izhod globalnega modela (En. (18))
zapišemo v matrični obliki:
ŷk = ψkΘ. (23)
Relacija izhoda modela in vhodnih podatkov za celoten
nabor zapišemo kot:
Ŷ = ΨΘ, (24)
kjer je Ŷ vektor izhodov modela ŷk (k = 1, . . . , n)
Ŷ = [ŷ1, ŷ2, . . . , ŷn]
T
. (25)
Mehki model, podan z enačbo En. (23), imenujemo
afini model Takagi-Sugeno. S tem modelom lahko s
poljubno natančnostjo aproksimiramo poljubno funkcijo
[8], [9], [10]. Splošnost modela lahko dokažemo s Stone-
Weierstrassovim teoremom [11]. Ta pravi, da lahko
vsako zvezno funkcijo aproksimiramo z razširitvijo
osnovnih mehkih funkcij [12].
1.1.4 Ocena parametrov mehkega modela: Za oceno
parametrov mehkega modela bo uporabljena metoda naj-
manjših kvadratov. Meritve zadoščajo nelinearni enačbi
sistema:
yi = g(xi), i = 1, . . . , n (26)
Po Stone-Weierstrassevm teoremu obstaja za realno zve-
zno funkcijo g, definirano nad naborom U c ⊂ Rp, mehki
sistem f tako, da je
max
xi∈X
|f(xi)− g(xi)| < δ, ∀i, (27)
kjer je δ > 0 poljubno majhna konstanta. Pri aproksima-
ciji zveznih funkcij z mehkimi funkcijami iz razreda Fp,
ki so definirane z En. (23), je treba opozoriti, da manjše
vrednosti δ povečajo število rojev c, da bi zadovoljili
En. (27). Pri aproksimaciji funkcije z mehkim modelom
je napaka med meritvami in vrednostmi izhoda mehkega
modela definirana kot:
ei = yi − f(xi) = yi − ŷi, i = 1, ..., n, (28)
kjer yi pomeni merjene izhode sistema, ŷk pa izhode
mehkega modela v trenutku k. Parametre mehke funkcije
(Θ) dobimo z minimizacijo vsote kvadratov napak:
E =
n∑
i=1
e2i =
= (Y − Ŷ )T (Y − Ŷ ) == (Y −ΨΘ)T (Y −ΨΘ).
(29)
Parameter Θ dobimo kot parcialni odvod ∂E
∂Θ
= 0:
Θ =
(
ΨT Ψ
)−1
ΨTY.
2 BIOLOŠKA ČISTILNA NAPRAVA ZA
ČIŠČENJE ODPADNIH VODA
Procesi čiščenja odpadnih voda so veliki nelinearni
sistemi, ki so izpostavljeni velikim motnjam predvsem v
toku in obremenitvi. Hkrati pa se spreminja tudi sestava
odpadne vode, ki pride v čistilno napravo. Za nepristran-
sko primerjavo različnih sistemov oziroma strategij vo-
denja in sistemov spremljanja je bil razvit enoten simu-
lacijski model. Sestoji iz petih zaporedno vezanih reak-
torjev, ki se nato iztekajo v desetplastni sekundarni tank,
kjer se tvorijo usedline. Model procesa, pripadajoče
enačbe in podroben opis najdemo na spletni strani http :
//www.ensic.inpl − nancy.fr/COSTWWTP/. V
našem primeru smo spremljali del sistema, kjer je odpa-
dana voda čiščena (tanki). Shema procesa je prikazana
na sliki 1.
ZAZNAVANJE NAPAK IN SPREMLJANJE ČIŠČENJA ODPADNIH VODA NA PODLAGI MEHKEGA MODELA 145
ANOXIC TANKS
AERATION TANKS
Qin
Qair
Qw
Qout
tanki
Slika 1: Shema simuliranega procesa.
Za izračun mehkega modela obravnavanega dela pro-
cesa smo uporabili naslednje spremenljivke: koncentra-
cija amonijaka v vhodnem toku v proces Qin, ki je
označena kot CNH4Nin , koncentracija stopljenega kisika
v prvem tanku C1O2 , koncentracija stopljenega kisika v
drugem tanku C2O2 in koncentracija amonijaka v drugem
tanku CNH4Nout . Mehki model je bil zgrajen tako, da
aproksimira odvisost koncentracije amonijaka v drugem
tanku od preostalih merjenih spremenljivk:
CNH4Nout(k) = G
(
CNH4Nin(k), C
1
O2
(k), C2O2(k)
)
,
(30)
kjer G označuje nelinearno povezavo med merje-
nimi spremenljivkami. Prvih 15000 vzorcev (s časom
vzorčenja Ts = 120s) je bilo uporabljenih za gra-
dnjo (identifikacijo) mehkega modela Takagi-Sugeno.
Ob 17000. vzorcu začne počasi naraščati napaka na
senzorju, ki je odpravljena ob 18000. vzorcu. Senzor
za merjenje koncentracije amonijaka v drugem tanku
CNH4Nout se je torej pokvaril. K nominalnemu signalu
senzorja smo dodali eksponentni signal za simulacijo
napake. Celoten set meritev je prikazan na sliki 2.
0 0.5 1 1.5 2 2.5 3
x 10
4
0
50
100
Sample
C
N
H
4
N
in
0 0.5 1 1.5 2 2.5 3
x 10
4
0
5
10
Sample
C
1 O
2
0 0.5 1 1.5 2 2.5 3
x 10
4
0
5
10
Sample
C
2 O
2
0 0.5 1 1.5 2 2.5 3
x 10
4
0
5
10
Sample
C
N
H
4
N
o
u
t
Slika 2: Celoten nabor meritev. Koncentracija amonijaka v
pritoku v sistem CNH4Nin , koncentracija kisika v prvem tanku
C1O2 ter koncentracija kisika C
2
O2
in amonijaka CNH4Nout v
drugem tanku.
Mehki model smo zgradili na prvih 15000 vzorcih.
Izhod mehkega modela ĈNH4Nout in izhod procesa
CNH4Nout sta prikazana na sliki 3. Indeks za detekcijo
0 5000 10000 15000
0
1
2
3
4
5
6
7
8
9
10
Sample
C
(-
),
C
(-
-)
N
H
4
N
o
u
t
N
H
4
N
o
u
t
Slika 3: Verifikacija mehkega modela, kjer je ĈNH4Nout izhod
modela, CNH4Nout pa izhod procesa.
napake je definiran kot:
f =
(
CNH4Nout − ĈNH4Nout
ĈNH4Nout
)2
. (31)
Ko indeks detekcije preseže indeks tolerance napake, se
sproži alarm. Indeks tolerance je defeniran kot največji
indeks napake, pomnožen s konstanto: ftol = γmax f .
Največji indeks napake smo zaznali pri učenju modela
na podatkih, kjer ni bilo simulirane napake. V našem
primeru smo pomnožili maksimalno vrednost indeksa
napake s konstanto γ = 1, 5. To pomeni, da smo za
tolerančni indeks dobili vrednost ftol = 0, 15. Slika 4
prikazuje tolerančni indeks, indeks detekcije napake in
signal za alarm. Napako, ki nastane pri vzorcu 17000, je
predlagani sistem za spremljanje procesa detektiral pri
17556. vzorcu. Detekcija napake je zakasnjena, kar je
pri napakah, ki naraščajo počasi običajno.
0 0.5 1 1.5 2 2.5 3
x 10
4
0
0.2
0.4
0.6
0.8
1
Sample
f,
f
to
l
0 0.5 1 1.5 2 2.5 3
x 10
4
0
0.5
1
1.5
2
Sample
fa
u
lt
(
--
),
f
a
u
lt
d
e
te
c
te
d
(
-)
Slika 4: Indeks detekcije napake, tolerančni indeks ftol dejan-
ska in detektirana napaka.
146 DOVŽAN, LOGAR, HVALA, ŠKRJANC
3 SKLEP
V članku smo obravnavali detekcijo napake na sen-
zorju in spremljanje čiščenja odpadnih voda. Sistem za
detekcijo je bil realiziran z mehkim modelom Takagi-
Sugeno. Model smo identificirali s pomočjo Gustafson-
Kesslovega algorita rojenja, parametri pa so bili dobljeni
s pomočjo najmanjših kvadratov. Razvit sistem je bil te-
stiran na modelu procesa, kjer je bila simulirana napaka
na enem od senzorjev. Za gradnjo mehkega modela smo
uporabili naslednje meritve: vhodno koncentracijo amo-
nijaka in koncentracijo kisika v prvem tanku ter tempe-
raturo, koncentracijo kisika in koncentracijo amonijaka
v drugem. Zaradi narave napake, simulirane na senzorju
koncentracije amonijaka na drugem tanku, je bila le-ta
detektirana brez lažnih alarmov in z majhno zakasni-
tvijo. Glede na to, da je dinamika sistema izpostavljena
spreminjanju zaradi različnih obremenitev, padavin itd.,
bomo poskušali sistem za detekcijo napake implemen-
tirati z rekruzivno mehko identifikacijo, ki bo sposobna
prilagajati dinamiko sistema tem spremambam, hkrati pa
bo sistem sposoben detektirati tudi napake.