Bruttó reakciók felbontása diszkrét matematikai eszközökkel

Papp Dávid

2003. október 27.

Tartalomjegyzék

1. Bevezetés 2

2. Bruttó és elemi reakciók 4

3. Lineáris diofantoszi egyenletrendszerek 6 3.1. Naiv rekurzív megoldás . . . . . . . . . . . . . . . . . . . . . . 6 3.2. A Contejean–Devie-algoritmus . . . . . . . . . . . . . . . . . . 9 3.3. LP relaxáció . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.4. Az algoritmusok összehasonlítása . . . . . . . . . . . . . . . . 14

4. Elemi reakciók előállítása 18 4.1. Példák . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5. A felbontások előállítása 21

5.1. A keresési tér szűkítése . . . . . . . . . . . . . . . . . . . . . . 21 5.1.1. Biztosan végbemenő reakciók . . . . . . . . . . . . . . 21 5.1.2. Nem végbemenő reakciók . . . . . . . . . . . . . . . . . 22

5.2. A legrövidebb felbontások előállítása . . . . . . . . . . . . . . 25 5.3. Az összes felbontás előállítása . . . . . . . . . . . . . . . . . . 26 5.4. Példák . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

6. Összefoglalás, további tervek 29

A. Az algoritmusok implementációja 30 A.1. A naiv rekurzív algoritmus . . . . . . . . . . . . . . . . . . . . 30 A.2. A Contejean–Devie-algoritmus . . . . . . . . . . . . . . . . . . 32 A.3. Az LP relaxációs felbontóalgoritmus . . . . . . . . . . . . . . . 33 A.4. Felbontásban nem szereplő reakciók kiszűrése LP relaxációval . 36 A.5. Volpert-indexelés . . . . . . . . . . . . . . . . . . . . . . . . . 37

1. Bevezetés

Bruttó reakcióknak nevezzük azokat a kémiai folyamatokat leíró egyenleteket, amik megmutatják, hogy a folyamat során milyen anyagokból milyen más anyagok keletkeznek. Ez a formalizmus a folyamat magas szintű leírását adja csak, nem tartalmazza azonban azt az információt, hogy a vizsgált reakció hogyan – milyen elemi lépéseken keresztül – megy végbe valójában. A bruttó reakciók felbontásának célja az, hogy megadjuk, a reakció milyen eleminek tekinthető lépések sorozataként állhat elő.

A felbontások meghatározása sok nagy jelentőségű reakció esetén régóta kutatási terület [1, 2], ennek ellenére igen kevés azoknak a reakcióknak a száma, amiket sikerült teljesen feltérképezni. (Jellemző, hogy a legjobban megismert reakciócsoport a szénhidrogének égési folyamatai [2], amelyek petrolkémiai alkalmazása közismert.)

A probléma több, önmagában is érdekes részfeladatra osztható. Meg kell határoznunk a reakcióban résztvevő anyagokat, valamint a belőlük felépíthető, eleminek tekinthető reakciókat. Ezeket korábban egyértelműen kémiai feladatnak tekintették, azaz a vizsgált bruttó reakció mellett (gyakran csupán kémiai intuícióra hagyatkozva) adottnak vették az elemi lépések listáját is. A dolgozat 4. fejezetében megvizsgáljuk, miképpen támogatható ez a feladat matematikai, illetve számítástechnikai eszközökkel.

A feladat következő lépése a bruttó reakció felbontása elemi lépésekre az összes lehetséges módon. A feladat több részből áll, melyek közül a legnehezebb egy olyan lineáris diofantoszi egyenletrendszer megoldása a természetes számok halmazán, amelynél a megoldások száma tipikusan végtelen, ezek azonban véges sok megoldással reprezentálhatók. Sok esetben ennek a véges megoldáshalmaznak a meghatározása is szinte reménytelen – annak ellenére, hogy a megkapható megoldások egy része kémiailag értelmetlen. Az ilyen reakciók kiszűrése egy másik fontos részfeladat, aminek egyik fajtája egy gráfindexelési problémára vezethető vissza. A lineáris diofantoszi egyenletrendszerekkel a 3. fejezetben, a felbontások előállításával és kémiai alapon végezhető szűrésével pedig az 5. fejezetben foglalkozunk részletesen.

A felbontások keresése során felmerülő algoritmikus problémák nagy része NP-nehéz (illetve NP-teljes), miközben a megoldások száma is rendkívül nagy, ezért a feladat megoldására nem várható hatékony algoritmus. A lineáris diofantoszi egyenletrendszerek megoldásával számos szerző foglalkozik, a felhasznált módszerek skálája igen széles, a legelemibb heurisztikus módszerektől a mély algebrai ismeretek felhasználásáig terjed [3, 4, 5, 6]. Az ezekben a munkákban leírt algoritmusok egy része csak a megoldhatóság eldöntésére használható (és nem az összes megoldás megkeresésére.) Ráadásul általában nem vizsgált feltétel a megoldások nemnegativitása, ami lényegesen nehezíti a feladatot.

Általánosan használható programcsomag a kémiai reakciók összes felbontásának előállítására mindeddig nem készült, a meglévő eredmények többnyire ad hoc megoldásokra épülnek, és felhasználják a (jellemzően csak szerves reakciók esetén) rendelkezésre álló termodinamikai adatokat. Célunk olyan diszkrét matematikai algoritmusok kidolgozása, vizsgálata és implementálása volt, amelyek általánosan használhatók kevés rendelkezésre álló kémiai adat esetén is, ugyanakkor lehetőséget adnak a felhasználónak kiegészítő feltételek (vagy akár teljes algoritmusok) beépítésére is a felbontás során.

Az 5.3. szakaszban bemutatjuk néhány szervetlen reakció felbontását is.

2. Bruttó és elemi reakciók

Elemi reakciónak azokat a kémiai reakciókat nevezzük, amelyekben legfeljebb két molekula (atom, ion) vesz részt. Ahhoz, hogy a vizsgált bruttó reakciót ilyenekre bontsuk, meg kell határoznunk, hogy a folyamat során milyen elemi reakciók mehetnek végbe. A szerzők rendszerint nem foglalkoznak ezzel a kérdéssel, a vizsgálódások során feltételezik, hogy bizonyos kémiai ismeretek alapján (illetve gyakran intuícióra hagyatkozva) a lehetséges elemi reakciók listája könnyen meghatározható. Ennek gyakran az az eredménye, hogy a folyamat modellje a valóságosnál sokkal egyszerűbb, a felhasznált elemi lépések száma olykor nem vagy alig több a reakcióban résztvevő anyagok számánál (lásd [7] és hivatkozásai, [12]) A következőkben megvizsgáljuk, hogyan lehet (legalábbis részlegesen) automatizálni ezt a lépést annak érdekében, hogy valóságosabb modellel dolgozhassunk.

Az első probléma annak meghatározása, hogy milyen anyagok vehetnek részt a reakciókban. Nyilvánvalóan ezek között van a felbontandó reakcióban szereplő összes anyag, de ezeken kívül további anyagokkal (pl. köztes lépések termékeivel) is számolnunk kell. Az előforduló elemek és vegyértékeik ismeretében algoritmikusan is előállítható egy lista azokról a molekulákról és ionokról, amelyből kiválaszthatók a szóba jövő anyagok. Ezek száma azonban rendszerint igen kicsi (legfeljebb néhány tucat), ráadásul az előbbiek helyett tisztán kémiai megfontolások és mérések segítségével is meghatározhatók; ezzel tehát a továbbiakban nem foglalkozunk. Feltesszük, hogy a résztvevő anyagok listája adott, a kérdés csak az, hogy ezekből az anyagokból milyen elemi reakciókat építhetünk fel.

A reakciókra vonatkozó egyetlen általánosan (azaz a felbontandó reakciótól függetlenül) használható megszorítás az atom- és töltésmegmaradás, ami a következőképpen formalizálható. Legyen k a vizsgált anyagok felépítésében résztvevő kémiai elemek száma. Rendeljünk minden anyaghoz egy k + 1 hosszúságú vektort, amelynek i-edik komponense 1 ≤ i ≤ k esetén megadja az anyagban található i-edik elem mennyiségét, az utolsó komponens pedig az anyag töltése. Ekkor az atom- és töltésmegmaradásnak megfelelő reakciók két oldalán szereplő anyagokhoz rendelt vektorok (reakcióbeli) együtthatóikkal súlyozott összege egyenlő. Az elemi reakciók előbbi definíciójával összevetve a feladatunk az összes felhasznált anyag vektorának, valamint a vektorpárok összegének előállítása a többi vektor segítségével. Ez lineáris diofantoszi egyenletrendszerek megoldását jelenti a természetes számok felett. Bizonyos esetekben egyéb – a tömegmegmaradásnak nem eleget tevő – reakciókat is figyelembe kell vennünk, pl. ha bizonyos reakciók során a rendszeren kívülről anyagok kerülhetnek a folyamatba, illetve távozhatnak nyom nélkül. Könnyen meggondolható, hogy ez a nullvektor előállítására vezethető vissza. A dolgozat későbbi fejezeteiben erre nem utalunk – az ott írottak módosítás nélkül alkalmazhatóak ilyen helyzetekben is.

Ha rendelkezésre áll az összes elemi reakció, hozzáláthatunk a bruttó reakció felbontásához. Az előbbiekhez hasonlóan járhatunk el: rendeljünk minden reakcióhoz egy m hosszú vektort, ahol m a résztvevő anyagok száma. Egy reakcióhoz rendelt vektor i-edik (1 ≤ i ≤ m) komponensének abszolút értéke legyen az i-edik anyagfajta együtthatója a reakcióban. Ha az i-edik anyag a reakció bal oldalán található, akkor a megfelelő komponens legyen negatív, ellenkező esetben pozitív. Egy ugyanilyen vektor a bruttó reakcióhoz is hozzárendelhető. A feladatunk ekkor a bruttó reakcióhoz rendelt vektor előállítása az elemi reakciók vektorának nemnegatív együtthatós lineáris kombinációjaként az összes lehetséges módon. Ez is egy lineáris diofantoszi egyenletrendszer megoldását jelenti a természetes számok felett. Az előzőekhez képest azonban több nehézség lép fel: ugyan csupán egyetlen, de sokkal nagyobb egyenletrendszert kell megoldanunk, másrészt az utóbbi egyenlet mátrixában tetszőleges helyen állhatnak pozitív és negatív elemek (míg az elemi reakciók előállításának feladatában csak a töltéseket jelző együtthatók lehetnek negatívak), végül az utóbbi egyenletrendszernek általában végtelen sok megoldása van.

A következő fejezetben részletesen megvizsgáljuk, hogyan kezelhető az ilyen kémiai feladatok során előálló egyenletek megoldása.

3. Lineáris diofantoszi egyenletrendszerek

Formálisan a következő a feladatunk.

Probléma. Legyen adott egy egészekből álló b ∈ ℤ^{d} vektor (a felbontandó vektor), továbbá egész vektorok egy {v_{1}, v_{2}, …, v_{n}} halmaza (a particionáló vektorok, v_{i} ∈ ℤ^{d}). Keressük mindazokat az (α_{1}, α_{2}, …, α_{n}) ∈ ℕ^{n} természetes szám n-eseket, amelyekre

∑_{i=1}^{n} α_{i}v_{i} = b.  (1)

Mielőtt hozzáfognánk a megoldáshoz, megemlítjük, hogy a probléma NPnehéz. Ennek belátásához elegendő észrevenni, hogy ha a szóban forgó vektorok egydimenziósak (d = 1), akkor annak eldöntése, hogy létezik-e megoldás, az NP-teljes részhalmaz-összeg probléma [8, 8.7 szakasz]. Így a megoldhatóság eldöntése a mi feladatunkban is nyilvánvalóan NP-teljes.

Bár ennek a feladatnak dolgozatunk témáján túlmenően is számos alkalmazása van, a diofantoszi egyenletekkel foglalkozó irodalom többnyire csak a magasabbfokú egyenletek megoldásának segédeszközként fogalkozik a lineáris egyenletekkel, mint „egyszerű esettel”, az összes nemnegatív megoldás előállításának hatékony módszereiről a nagyobb összefoglaló munkákban ([4, 5]) sem olvashatunk. Az elterjedt matematikai programcsomagok, mint a Mathematica és a Maple sem tartalmaznak olyan függvényt, amivel a megoldások előállíthatók.

Elsőként egy egyszerű algoritmust mutatunk be a megoldások felsorolására.

3.1. Naiv rekurzív megoldás

Ha az (1) egyenletben szereplő vektorok egydimenziósak, akkor a klasszikus számparticionálási feladathoz jutunk [9]. Az összes felbontás előállítására ebben az esetben a Mathematica Partitions függvénye használható, ami egy rendkívül egyszerű rekurzív algoritmust implementál. Ezt a következőképpen általánosíthatjuk vektorokra.

1. [teszt] Ha nincs particionáló vektorunk, akkor ha b a nullvektor, akkor a megoldás egy „üres” együtthatólista, különben nincs megoldás. Mindkét esetben megállunk.

2. [elágazás, rekurzió] Ha az első particionáló vektornak (legyen ez v) nincs olyan pozitív (ill. negatív) komponense, amely nagyobb abszolút értékű a b megfelelő komponensénél, vagy ha a többi particionáló vektor között van olyan, amelyiknek ugyanezen komponense negatív (ill. pozitív), akkor particionáljuk az összes lehetséges módon b − v-t, majd vegyük hozzá a felbontásokhoz a v vektort.

3. [rekurzió] Dobjuk el a particionáló vektorok közül v-t, majd particionáljuk az összes lehetséges módon b-t.

4. [befejezés] Az előző két pontban kapott megoldáshalmaz uniója a megoldás.

A fenti eljárás lényegében azt a nemdetermininisztikus algoritmust hajtja végre, amely minden vektort sorra véve elágazik aszerint, hogy hányszor használjuk fel azt a felbontásban. Hogy az eljárás biztosan véget érjen, minden vektorról meg kell tudnunk állapítani, hogy hányszor használhatjuk még a felbontás során, ellenkező esetben valamelyik particionáló vektornál az eljárás végtelen ciklusba kerül.^{1} (Ezért is használtuk az eljárás szót az algoritmus helyett.) A szükséges ellenőrzést a 2. pontba kell beépítenünk. Ha a v vektornak nincs olyan komponense, amelyik alapján biztosítható, hogy a vektort véges sok lépés után eldobjuk, akkor az algoritmus hibaüzenettel leáll. Ennek a módosításnak az az ára, hogy az algoritmus nem használható minden esetben, így például akkor sem, ha minden 1 ≤ i ≤ d-hez létezik egy particionáló vektor, amelynek i-edik komponense negatív és egy másik, amelynek i-edik komponense pozitív. Szerencsére az elemi reakciók előállításánál ez a probléma nem lép fel, hiszen kizárólag az anyagok töltését mutató komponensek lehetnek negatívak.

Vegyük észre, hogy problémás particionáló vektorok esetében lehetséges, hogy pusztán a vektorok permutálásával elérhetjük, hogy az eljárás termináljon. Legyen például n = d = 2, v_{1} = (1, 0)^{T}, v_{2} = (−1, 1)^{T}. Ebben a sorrendben semmilyen b vektor esetén nem tudjuk közvetlenül alkalmazni az algoritmust, míg a vektorokat megcserélve minden jobb oldal esetén megkapjuk az összes megoldást. A terminálás szükséges és elégséges feltétele, hogy minden vektornak legyen egy pozitív (ill. negatív) komponense úgy, hogy az utána következő vektorok megfelelő komponensei mind nemnegatívak (ill. nempozitívak.) Nevezzük a {v_{i}}_{i=1}^{n} sorozatot jól rendezettnek, ha minden vektornak van ilyen (jó) komponense. A következő segédtétel lehet a kulcsa egy gyors, jól rendező algoritmusnak.

1 Az egészek particionálásánál ilyen probléma természetesen nem lép fel, a gondot a vektorokban megengedett különböző előjelű komponensek okozzák.

Lemma.

1. Ha minden k-ra a sorozat tagjainak k-adik komponensei között pozitív és negatív is található, vagy ha a sorozat valamelyik tagja nullvektor, akkor a sorozat nem jól rendezhető.

2. Legyen k egy olyan index, amire minden vektor k-adik komponense nemnegatív (ill. nempozitív.) A sorozat pontosan akkor jól rendezhető, ha elhagyva mindazokat a vektorokat, amelyek k-adik komponense nem nulla, egy jól rendezhető sorozatot kapunk. (Az elem nélküli sorozatot jól rendezettnek tekintjük.)

Bizonyítás.

1. Az első helyre kerülő vektornak, illetve a nullvektornak semmilyen permutáció esetén nincs jó komponense.

2. Az „akkor” irány bizonyítása: ha a redukált sorozat rendezhető, rendezzük azt, majd tegyük vissza az elhagyott vektorokat a sorozat elejére. Ez az eredeti sorozat egy jó rendezése, mert az el nem hagyott vektorok rendezésének jóságán nem változtatnak az eléjük tett vektorok, az elhagyott (majd visszatett) vektoroknak pedig van jó komponensük – az a nemnulla komponens, ami alapján kiválasztottuk őket. A másik irányú igazoláshoz indirekten tegyük fel, hogy az eredeti sorozat jól rendezhető, míg a redukált sorozat nem. Ez azt jelenti, hogy van olyan vektor, aminek egy jó komponense megszűnik jónak lenni azáltal, hogy más vektorokat törlünk a sorozatból. Ez nyilvánvalóan lehetetlen. ■

A segédtételből egy polinomidejű rendezőalgoritmus is kibontakozik: ha a vektoraink nem jól rendezettek, ellenőrizzük, hogy teljesül-e az előbbi lemma első pontjának feltétele, s ha nem, akkor a második pont alapján próbáljuk rekurzívan rendezni a sorozat elemeit.

A naiv felbontó algoritmus és a jó rendezéshez szükséges függvények Mathematica implementációját lásd az A.1. szakaszban.

Ez az algoritmus csak akkor működőképes, ha a megoldások száma véges. A következő szakaszban mutatunk egy, az irodalomban gyakran hivatkozott megoldást, amellyel kezelhetők a végtelen sok megoldással rendelkező egyenletek is.

3.2. A Contejean–Devie-algoritmus

Contejean és Devie [3] cikke (elsősorban) homogén lineáris egyenletrendszerek megoldásával foglalkozik a természetes számok felett, ahol ha van megoldás, szükségképpen végtelen sok van. Hogy ilyen esetben is legyen értelme az összes megoldás előállításáról beszélni, a szerzők bevezetik a minimális megoldás fogalmát.

Definíció. Az (1) egyenlet egy (α_{1}, α_{2}, …, α_{n}) megoldását minimálisnak nevezzük, ha nem található hozzá (β_{1}, β_{2}, …, β_{n}) megoldás, hogy ∀ i : α_{i} ≥ β_{i}, ∃ i : α_{i} ≠ β_{i} és ∃ j : β_{j} > 0 (jelben: (α_{1}, α_{2}, …, α_{n}) ⋗ (β_{1}, β_{2}, …, β_{n})) legyen.

Az említett cikkben nem szerepel, nekünk azonban szükségünk lesz a következő tételre.

Tétel. Az (1) egyenlet megoldásainak száma akkor és csak akkor véges, ha minden megoldás minimális.

Bizonyítás. Ha van nem minimális megoldás, azaz a, b megoldások, és a ⋗ b, akkor b + k(a − b), k ∈ ℕ végtelen sok különböző megoldás. A másik irányban be kell látnunk, hogy minimális megoldásból csak véges sok lehet. Tegyük fel indirekte, hogy (megszámlálhatóan) végtelen sok van. Ekkor a megoldások sorozatából kiválasztható végtelen sok, amelyek első koordinátája monoton növekedő sorozatot alkot. Ezek közül a vektorok közül kiválasztható végtelen sok, amelyek második koordinátája (is) monoton növekedő sorozatot alkot. Ezt folytatva beláthatjuk, hogy a vektorok közül kiválasztható végtelen sok, amelyek minden koordinátája monoton növekedő sorozatot alkot, ezek nem lehetnek mind minimális megoldások. ■

A tétel jelentősége az, hogy rámutat: minden esetben értelmes az (1) egyenlet összes megoldása helyett csak a minimális megoldásokat keresni, hiszen véges sok megoldás esetén a két halmaz megegyezik, ha pedig végtelen sok megoldás van, úgy ezeket véges módon reprezentálja a minimális megoldások M és a megfelelő homogén egyenletrendszer minimális megoldásainak M_{0} halmaza, hiszen minden megoldás az M_{0} ∪ M halmaz elemeinek nemnegatív együtthatós lineáris kombinációja.

Az algoritmus egy, az (1) egyenlet minimális megoldásait szolgáltató algoritmus ötletes gyorsításán alapul. A homogén egyenletrendszer megoldására szolgáló alapvető algoritmus a következő [3]:

1. [inicializálás] Legyen az S halmaz az n dimenziós nullvektort tartalmazó, egyelemű halmaz.

2. [ciklus] Minden s ∈ S vektor minden komponensét növeljük (különkülön) 1-gyel, az így kapott vektorokat tegyük S-be, az s vektort vegyük ki onnan.

3. [teszt] Ha S-ben találunk megoldásvektorokat, azokat vegyük ki a halmazból. Ugyanígy töröljük azokat a vektorokat, amik (a ⋗ reláció szerint) nagyobbak valamelyik eddig megtalált megoldásnál.

4. [ciklus] Ha S üres, álljunk meg, minden megoldást megtaláltunk. Különben ugorjunk a második pontra.

Az algoritmus módosítható úgy, hogy inhomogén egyenletre is alkalmazható legyen. Ehhez vegyük a felbontandó b vektor ellentettjét a particionáló vektorok közé (pl. utolsónak), az első (inicializáló) lépésben a kiindulási vektor a nullvektor helyett legyen az n + 1 dimenziós (\underbrace{0, …, 0}_{n}, 1) vektor, majd

folytassuk az algoritmust a megadott módon, de a vektorok utolsó komponensét soha nem növelve. Végül az S halmaz minden elemének utolsó komponensét „vágjuk le.”

Ennek a [3]-beli módosítása a következő: a második pontban az s = (s_{1}, s_{2}, …, s_{n}) vektor k-adik komponensét csak abban az esetben növeljük meg, ha (az (1) egyenletben használt jelölésekkel) a ⟨∑_{i=1}^{n} s_{i}v_{i}, v_{k}⟩ skalárszorzat negatív. Igazolható, hogy az algoritmus ezzel a módosítással is megtalál minden minimális megoldást, miközben a keresési fa elágazási tényezője átlagosan legalább a felére csökken.

Az algoritmus ezzel még mindig exponenciális méretű tárat használ, és tekintettel arra, hogy a megoldások száma exponenciálisan nagy lehet, ennél sokkal jobbat nem várhatunk. Az azonban elvárható (és el is érhető), hogy a megtalált minimális megoldások mentésén kívül felhasznált tárterület kicsi legyen. A [3] cikk közli az algoritmus egy rendkívül tárhatékony változatának pszeudo-kódját – ez az új megoldások kereséséhez csupán lineáris méretű tárat (egy, a particionáló vektorok számával megegyező mélységű stack-et) használ.

Mutatunk egy további módosítást, amivel tovább redukálható az algoritmus lépéseinek száma, ha felső becslést tudunk adni a minimális megoldások méretére (a megoldás méretén az együtthatók összegét értve.)

Lemma. Ha a homogén egyenlet legfeljebb m méretű megoldásait keressük, akkor a második pontban az s = (s_{1}, s_{2}, …, s_{n}) vektor k-adik komponensét elegendő csak abban az esetben megnövelnünk, ha (m > ∑_{i=1}^{n} s_{i}, és)

⟨∑_{i=1}^{n} s_{i}v_{i}, v_{k}⟩ ≤ −\frac{∥ ∑_{i=1}^{n} s_{i}v_{i}∥^{2}}{m − ∑_{i=1}^{n} s_{i}}  (2)

Bizonyítás. Legyen s^{∗} = (s_{1}^{∗}, s_{2}^{∗}, …, s_{n}^{∗}) egy minimális megoldás, s ⋖ s^{∗}, és legyen s_{i}^{∗} = s_{i} + r_{i}. Ekkor

0 = ∥ ∑_{i} s_{i}^{∗}v_{i}∥^{2} = ∥ ∑_{i} s_{i}v_{i}∥^{2} + ∥ ∑_{i} r_{i}v_{i}∥^{2} + 2⟨∑_{i} s_{i}v_{i}, ∑_{i} r_{i}v_{i}⟩
= 2∥ ∑_{i} s_{i}v_{i}∥^{2} + 2⟨∑_{i} s_{i}v_{i}, ∑_{i} r_{i}v_{i}⟩

Ebből

−∥ ∑_{i} s_{i}v_{i}∥^{2} = ⟨∑_{i} s_{i}v_{i}, ∑_{i} r_{i}v_{i}⟩
= ∑_{k} r_{k}⟨∑_{i} s_{i}v_{i}, v_{k}⟩

A jobb oldali összeg egy ∑ r_{k} < m tagú összeg, amelynek minden tagja ⟨∑_{i=1}^{n} s_{i}v_{i}, v_{k}⟩ alakú. Ez csak akkor lehet egyenlő a bal oldalon álló kifejezéssel, ha valamelyik tagja nem nagyobb a bal oldal (1/ ∑ r_{k})-szeresénél. Ez a tag éppen kielégíti a feltételt. Tehát az algoritmus nem akad el, ha a második pontban az s vektornak csak a feltételnek eleget tevő komponenseit növeljük. Másrészt így egyetlen megoldás sem vész el, hiszen a vektorok összeadása asszociatív és kommutatív, így teljesen mindegy, hogy az s vektorból milyen lépések sorozatával jutunk el az s^{∗}-hoz. ■

Megjegyezzük, hogy az előbbiek szerint a tétel nem csak a homogén egyenletek megoldását gyorsítja, hiszen az inhomogén egyenlet megoldását egy homogén egyenlet megoldására vezettük vissza.

A kapott eredmény m → ∞ esetén éppen az idézett, [3]-beli feltétel, és első látásra nem is sokkal erősebb annál. A 3.4. szakaszban látni fogjuk, hogy valójában igen sok lépést megspórolhatunk ennek az erősebb korlátnak az alkalmazásával, ha (kémiai vagy matematikai megfontolásokkal) becsülni tudjuk a minimális megoldások méretét.

Az algoritmus erősebb korláttal dolgozó, tárhatékony Mathematica implementációját lásd az A.2. szakaszban.

3.3. LP relaxáció

A harmadik vizsgált algoritmus egy egészértékű programozási és rokon feladatoknál gyakran alkalmazott [10] ötleten, az egészértékű feladat valós relaxációján alapul. Az alapelv a következő: meghatározzuk az (1) egyenletrendszer valós számok fölötti általános megoldását, mint egy partikuláris megoldás és a megfelelő homogén egyenletrendszer általános megoldásának összegét, majd a kapott általános megoldásvektor minden komponensére felírjuk a nemnegativitási feltételt. Az így kapott egyenlőtlenségrendszer megoldásai közül kiválasztjuk azokat, amelyekre a megoldásvektor komponensei egészek. A megoldáshoz felhasználjuk az alábbi, egyszerűen igazolható állítást.

Lemma. A homogén egyenletrendszer megoldásainak létezik olyan {b_{1}, b_{2}, …, b_{k}} bázisa, amelyre teljesül, hogy minden i-re (1 ≤ i ≤ k) b_{i}-nek létezik olyan komponense, amelyik a kiválasztott b_{i} esetén pozitív egész, míg az összes többi bázisvektor ugyanezen komponense nulla.

Bizonyítás. Feltéve, hogy már meghatároztunk egy racionális vektorokból álló bázist, az első feltétel biztosítható úgy, hogy a bázisvektorokból (mint sorvektorokból) álló mátrix sorait megfelelően Gauss-elimináljuk. Ez az elemek racionalitásán nem változtat. A második feltétel teljesítése ezek után triviális, a nem megfelelő sorvektorokat −1-gyel szorozzuk. Végül a racionális bázisból megkapható az egész vektorokból álló bázis a nevezők legkisebb közös többszörösével történő szorzással. ■

Tapasztalatunk szerint a Mathematica NullSpace függvénye épp egy ilyen bázist szolgáltat.

Ha az inhomogén egyenlet egy partikuláris megoldása p = (p_{1}, p_{2}, …, p_{n}), (p_{i} ∈ ℚ), akkor az inhomogén egyenlet megoldásának általános alakja az előbbi lemmában definiált bázisvektorok segítségével így írható:

v = p + b_{1}x_{1} + ⋯ + b_{k}x_{k} (x_{i} ∈ ℝ).

Meg kell keresnünk azokat az x = (x_{1}, x_{2}, …, x_{k}) vektorokat, amelyekre v minden komponense pozitív egész. Most használhatjuk a b_{i}-kre vonatkozó előbbi feltételünket. Ha ugyanis β_{i} ∈ ℤ^{+} jelöli a b_{i} megfelelő komponensét (és ez az l-edik komponens), akkor a v l-edik komponense p_{l} + β_{i}x_{i} egész, tehát x_{i} szükségképpen egy \frac{t−p_{l}}{β_{i}} (t ∈ ℤ) alakú racionális szám. Ezek után ha találunk alsó és felső korlátot az x komponenseire, akkor az összes (immáron véges sok) esetet megvizsgálva a határok között megkereshetjük az összes megoldást.

A határokat lineáris programozással kereshetjük meg. A v ≥ 0 egyenlőtlenségrendszer, mint korlátozó feltételek mellett keressük az egyes x_{i}-k minimumait és maximumait. (Ez 2k külön programozási feladat.) A legegyszerűbb, ha ezen határok között végigszaladunk minden lehetséges (x_{1}, …, x_{k}) kombináción, ellenőrizve, hogy csakugyan minden kombináció nemnegatív egészekből álló v-t ad-e. Ez azonban nagyon időigényes lehet, ha a lineáris programozással kapott keresési tér nagy.

Ezért kezdetben csupán az egyik – például az x_{1} – változóra számítjuk ki a korlátokat, majd x_{1} minden szóbajövő értékére külön-külön meghatározzuk az x_{2}-re adható alsó és felső korlátokat, és így tovább. E módszer előnye, hogy így csak olyan szám n-eseket kapunk, amelyeket az általános megoldásba helyettesítve nemnegatív megoldásokat kapunk. Mivel pedig a határok között csak a megfelelő alakú racionális számokat vizsgáljuk meg, ezek egészek is – végső soron tehát nincs szükség semmiféle ellenőrzésre, csupán a tartománybeli szám n-esek visszahelyettesítésére, és a megoldások felsorolására. Látszólag nagy többletköltséget jelenthet, hogy az előbbinél lényegesen több lineáris programot kell megoldanunk – valójában ez a megoldás nagyságrendekkel gyorsabbnak bizonyult, ami azzal magyarázható, hogy a viszonylag gyors lineáris programozási rutin minden meghívásával exponenciálisan sok rossz (negatív számokat eredményező) „megoldást” metszhetünk le a keresési térből.

Fontos megemlítenünk, hogy az algoritmus mindhárom lényeges része, a partikuláris megoldás megkeresése, a homogén egyenlet megoldása és a lineáris programozási feladatok megoldása egyaránt tartalmaz numerikus buktatókat. A Mathematica felhasznált rutinjai azonban mind futtathatók egész típusú számokkal, amivel kicsivel lassabb, de mindig pontos (egész, illetve racionális típusú) eredményeket adó eljárásokat kapunk. (Az eljárások dupla pontosságú valós számokként megadott egészekkel nehezebb kémiai feladatok esetén numerikus hibák miatt hibás megoldásokat adtak.)

Eddig nem foglalkoztunk azzal, mi történjen akkor, ha valamelyik x_{i} változóra nem kapunk felső korlátot (alsót nyilvánvalóan kapunk, ez az előbb felhasznált lemmából azonnal következik.) Ez akkor következik be, ha az egyenlet (valós) megoldásának valamelyik komponense tetszőlegesen nagy lehet. Meggondolható, hogy ilyenkor a megfelelő komponens az egészek halmazán is tetszőlegesen nagy lehet. Ez a probléma tehát csak akkor léphet fel, ha végtelen sok megoldás van, melyek felsorolása amúgy is reménytelen. Ilyenkor a Contejean–Devie-algoritmushoz hasonlóan kereshetnénk a minimális megoldásokat, erre ez az algoritmus azonban közvetlenül nem használható. Később megvizsgáljuk, hogy hogyan használható mégis ez az algoritmus, ha a megoldandó kémiai feladat ilyen egyenletre vezet.

Az eljárások Mathematica forráskódja megtalálható az A.3. szakaszban.

3.4. Az algoritmusok összehasonlítása

Ebben a pontban elméleti megfontolások és számítási példák segítségével megvizsgáljuk, mire képesek az eddigiek során ismertetett algoritmusok.

A megoldható feladatok. A három algoritmus között a legszembetűnőbb különbség, hogy különböző feladatok megoldására használhatók. A legáltalánosabban használható megoldás a Contejean–Devie-algoritmus, ami tetszőleges egyenlet összes minimális megoldását megadja, ami – mint láttuk – azt is jelenti, hogy az összes megoldást megadja, ha véges sok van.

Az LP relaxációs algoritmus minden megoldást megad, ha véges sok van, azonban közvetlenül nem használható, ha végtelen sok megoldás van. A kémiai reakciók felbontásának meghatározásakor ez azonban nem jelentős hátrány, ezt részletesebben a következő két fejezetben vizsgáljuk.

Az első, rekurzív algoritmus akkor sem mindig alkalmazható, ha a megoldások száma véges. Példaként legyen n = d = 2, v_{1} = (−1, 2)^{T}, v_{2} = (2, −1)^{T} és b = (1, 1)^{T}. Egyetlen megoldás van, de a két particionáló vektor nem jól rendezhető, ezért az algoritmusunk nem használható.

Amint azt a 3.1. szakaszban is említettük, ez a probléma nem lép fel, ha valamelyik komponens minden particionáló vektorban egyforma előjelű, vagy ha van néhány olyan komponens, amelyekre igaz, hogy mindegyik partíconáló vektorban legalább az egyik pozitív és egyik sem negatív.

Idő- és tárhatékonyság. Korábban igazoltuk, hogy az (1) egyenlet megoldhatóságának eldöntése NP-teljes, ráadásul a megoldások felsorolása – azok nagy száma miatt – rendkívül idő- és memóriaigényes lehet, így hagyományos értelemben hatékony (polinomidejű) algoritmus a megoldásra nem várható. Meghatározott feladatosztályokra azonban gyors algoritmusokat kereshetünk. Az első algoritmus várhatóan mindaddig gyors, amíg a keresés során viszonylag kevésszer folytatódik rekurzív hívással a második pontban. Különösen jellemző ez az olyan feladatokra, amelyeknél a felbontandó vektor nem sokkal nagyobb a particionáló vektoroknál. Ellenkező esetben azonban nagyon nagy válaszidőre számíthatunk az alkalmazásakor.

A Contejean–Devie-algoritmus ismertetésénél említettük, hogy implementálható rendkívül tárhatékonyan (eltekintve persze a memóriában tárolt megoldások tárigényétől), azonban a futási ideje exponenciálisan nő a megoldások nagyságának növekedtével. Erre az algoritmusra tehát általában csak addig számíthatunk, amíg a megoldandó egyenletrendszer megoldásai között kis méretűek is vannak.

A harmadik algoritmus hatékonysága előre nehezen megjósolható. Az egyetlen, ami arra utal, hogy igen gyors lehet, az az a tény, hogy a lineáris programok megoldásával körülhatárolt térben minden meghatározott alakú racionális együtthatós vektor megoldás, tehát a megoldásokat közvetlenül, részleges felbontások nélkül kapjuk meg. Ez a másik két algoritmusnál nyilvánvalóan nem teljesül. A keresési fázis előtti szakasz – az általános megoldás meghatározása – szintén rendkívül gyors. Ugyanakkor a nagy számú lineáris program megoldásának szükségessége első látásra elrettentőnek tűnik.

A programok hatékonyságát valódi kémiai feladatokon később mutatjuk be, ebben a fejezetben csupán a Contejean–Devie-algoritmus eredeti és dolgozatunkban javított változatát hasonlítjuk össze egy – nem kémiai ihletésű – feladaton.

Példa. A feladat eredetileg [3]-ban jelent meg, de több más cikkben is idézik. Legyen a particionáló vektorokból (oszloponként) képzett mátrix

( \begin{matrix}0 & 0 & 0 & 0 & −2 & 1 & 1 \\ 0 & 0 & −4 & 1 & 1 & 1 & 1 \\ −6 & 1 & 1 & 1 & 1 & 1 & 1 \end{matrix} ) ,

és keressük a homogén egyenlet minimális megoldásait. Erre a feladatra a három bemutatott eljárás közül közvetlenül csak a Contejean–Devie-algoritmus alkalmazható. A függelékbeli implementáció a ciklus 14049-dik lefutása után ér véget, és adja vissza a 95 megoldást, melyek legnagyobbika 35 nagyságú. Ha előzetesen valamilyen módszerrel megbecsüljük a megoldások maximális méretét (esetleg nem vagyunk kiváncsiak egy adott méretűnél nagyobb megoldásokra), akkor sok lépést megspórolhatunk a módszer általunk javasolt módosításával. A legszerencsésebb esetben éppen 35-re állítjuk a felső korlátot, így a megoldásokat 8672 lépésben kapjuk meg, ami közel 40%-os nyereség. Tapasztalataink szerint mindaddig levághatjuk a keresési fa néhány ágát, amíg a felső korlát 86 alatt van. Ezt kvalitatíve azzal magyarázhatjuk, hogy bár a (2) egyenlőtlenség csak az m-hez közeli méretű részeredmények esetén javít érdemben az eredeti algoritmuson, sok lépést megtakaríthatunk pusztán az eredeti eljárás által bejárt keresési fa alsó szintjein végzett vágásokkal.

A feltételezés, hogy felső korlátot adhatunk a minimális megoldások méretére, nem olyan irreális, mint első pillantásra tűnhet. Kémiai szempontok is indokolhatnak egy felső korlátot (például nem keresünk száznál több elemi lépésből álló felbontásokat), másrészt ha a megoldások száma véges (így minden megoldás minimális), akkor az összes megoldás méretére lineáris programozással adhatunk becslést, hasonlóan ahhoz, ahogyan azt az LP relaxáción alapuló megoldásban is tettük.

Megjegyezzük, hogy (az 5. fejezetben ismertetendő módon) a másik két algoritmus is alkalmas a legfeljebb N méretű megoldások megkeresésére. (N tetszőleges.) Az első algoritmus ezen változata lényegesen gyengébben teljesít a javított Contejean–Devie-algoritmusnál, míg az LP relaxáció alkalmazásával gyorsabban megkaphatjuk a 100-nál nem nagyobb megoldásokat is. (Minimálisokat és nem minimálisokat egyaránt.) Az abszolút időeredményeket tekintve az általunk közölt implementációk nem összehasonlíthatók az eredetiekkel, mivel mi asztali számítógépet, míg [3] szerzői 9 éve többprocesszoros munkaállomásokat használtak. A cikkben közölt összes feladat esetén lassabbnak bizonyult a mi PC-s Contejean–Devie implementációnk, de az LP relaxációs megoldásunk minden [3]-beli feladat esetén lényegesen gyorsabb.

A dolgozatunk szempontjából igazán érdekes, reakciók felbontásából kapható példákkal a 4.1 és az 5.4. szakaszban foglalkozunk.

Extrém méretű példák esetén mindhárom algoritmusnál szükséges lehet, hogy a megoldásokat futás közben mentsük, hogy ne foglalják le az összes memóriát. Ez a módosítás az első és a harmadik algoritmus esetében nagyon egyszerűen elvégezhető, a merőben technikai részletek ismertetésétől eltekintünk. Megjegyezzük azonban, hogy a második algoritmusnál a megoldások keresése közben a részmegoldásokat lépésenként összevetjük a megtalált megoldásokkal. Itt a megoldások mentése és kisöprése a memóriából közvetlenül nem járható út.

4. Elemi reakciók előállítása

Ebben a fejezetben megvizsgáljuk, hogyan alkalmazhatók az előző pontban ismertetett algoritmusok az elemi reakciók előállítására, majd megvizsgálunk néhány gyakorlati problémát is.

Az elemi reakciók előállításakor fellépő egyenletek az (1) egyenlet speciális esetei. A megoldások kémiai jelentéséből is azonnal kiviláglik, hogy csak olyan egyenletek lépnek fel, amelyeknek garantáltan véges sok megoldása van, ezért nem kerülnek hátrányba azok a módszerek, amik végtelen sok megoldással rendelkező egyenletekkel általában nem boldogulnak. Emellett a particionáló vektoroknak csupán az utolsó komponense lehet negatív, és minden vektornak van pozitív komponense az első d − 1 komponens között. (Minden anyagban van legalább egy atom, de egyikben sem lehet negatív semelyik elem mennyisége.) Már láttuk, hogy ez elégséges ahhoz, hogy akár az első, rekurzív algoritmust is használhassuk a megoldások megkereséséhez. Tehát az elemi reakciók felbontására mindhárom algoritmus közvetlenül alkalmazható.

4.1. Példák

Permanganát–oxálsav reakció [11]. Ezt a reakciót több mint 100 éve vizsgálják, teljes felbontása mind a mai napig nem ismert. A reakcióban 19 anyagfajta vehet részt, amelyek összesen négy elemet: mangánt, hidrogént, szenet és oxigént tartalmaznak. A felbontandó reakció:

2MnO_{4}^{−} + 6H^{+} + 5H_{2}C_{2}O_{4} → 2Mn^{2+} + 8H_{2}O + 10CO_{2}

Az anyagok adatai az 1. táblázatban láthatók. Az összes elemi reakciót megkapjuk, ha minden anyag vektorát, mindegyik kétszeresét, végül minden vektorpár összegét az összes lehetséges módon particionáljuk a többi segítségével. (Összesen 209 egyenletrendszer.) Az előző fejezetbeli megfontolások alapján az első algoritmust alkalmaztuk az első 19 egyenletben (itt a felbontandó vektor és a megoldások még nagyon kicsik), a többi esetben LP relaxációt használtunk. A Contejean–Devie-algoritmus a legkisebb egyenleteknél versenyképesnek bizonyult az egyszerű rekurzív algoritmussal, de nagyobb feladatoknál alulmaradt a relaxációs módszerrel szemben. Összesen 1022 lehetséges elemi reakció van, ezeket a két saját algoritmus alkalmazásával 40 másodperc alatt kaptuk meg.^{2}

2 AMD Athlon 1667 MHz processzoros PC, 256 MB DDR RAM, Mathematica 5.0, Windows XP

C H Mn O t
H_{2}C_{2}O_{4}  2 2 0 4 0
HC_{2}O_{4}^{−}  2 1 0 4 −1
H^{+}  0 1 0 0 1
C_{2}O_{4}^{2−}  2 0 0 4 −2
Mn^{2+}  0 0 1 0 2
MnC_{2}O_{4}  2 0 1 4 0
MnO_{4}^{−}  0 0 1 4 −1
MnO_{2}  0 0 1 2 0
Mn^{3+}  0 0 1 0 3
CO_{2}  1 0 0 2 0
H_{2}O  0 2 0 1 0
[MnO_{2}, H_{2}C_{2}O_{4}]  2 2 1 6 0
CO_{2}^{−}  1 0 0 2 −1
[Mn(C_{2}O_{4})]^{+}  2 0 1 4 1
[Mn(C_{2}O_{4})_{2}]^{−}  4 0 1 8 −1
[MnC_{2}O_{4}, MnO_{4}^{−}, H^{+}]  2 1 2 8 0
[MnC_{2}O_{4}^{2+}, MnO_{3}^{−}]^{+}  2 0 2 7 1
[MnC_{2}O_{4}^{2+}, MnO_{3}^{−}, H^{+}]^{2+} 2 1 2 7 2
[H^{+}, MnO_{2}, H_{2}C_{2}O_{4}]^{+}  2 3 1 6 1

1. táblázat. A permanganát–oxálsav reakcióban résztvevő anyagok

Levegőkémiai folyamatok. Egy, a leeds-i egyetem Üzemanyag és Energia Tanszékén folyó levegőkémiai kutatás számára készítettük el a 2. táblázatban felsorolt anyagok között végbemenő összes olyan elemi reakció listáját, amelyben szerepel ként tartalmazó anyag.

A feltételből 1668, egyenként 78 illetve 79 változós egyenletrendszer írható fel, a megoldáshoz itt is az előbb említett algoritmusokat használtuk. Mivel azonban a megoldások négyszázezret meghaladó száma kezelhetetlenül nagynak bizonyult a további vizsgálódások szempontjából, az összes elemi reakció közül ki kellett szűrnünk azokat, amiknek háromnál több terméke van.

A továbbfejlesztett Contejean–Devie-algoritmussal közvetlenül (az összes reakció meghatározása nélkül) megoldható az összes, legfeljebb három termékű minimális reakció előállítása (ezek az algoritmus nemdeterminisztikus változatával három lépésen belül elérhető megoldások), de amint arra már utaltunk, a másik két algoritmust is használhatjuk, hogy minden (nem feltétlenül minimális) reakciólépést megadjunk. (Ld. az 5. fejezetben.) Így gyorsabban megkaphatjuk a megoldásokat, mint a Contejean–Devie-módszerrel,

H_{2}  CH_{4}  C_{2}H_{2}  C_{2}H_{4} C_{2}H_{6} C_{3}H_{4} C_{3}H_{6}  C_{4}H_{2}
O_{2}  H_{2}O  H_{2}O_{2}  CO  CO_{2} CH_{2}O CH_{2}CO C
H  CH  CH_{2}  CH_{2}(S) CH_{3} C_{2}H C_{2}H_{3}  C_{2}H_{5}
C_{3}H_{2}  H_{2}CCCH H_{2}CCCCH O  OH HO_{2} HCO  CH_{3}O
CH_{2}OH HCCO  CH_{2}HCO N_{2}  CN HCN N  NH
NO  HNO  NH_{2}  H_{2}NO NCO N_{2}O NO_{2}  N_{2}H_{2}
HOCN H_{2}CN  NNH  NH_{3}  N_{2}H_{3} C_{2}N_{2} HNCO NO_{3}
S  SH  H_{2}S  SO  SO_{2} SO_{3}  HSO_{2} HOSO
HOSO_{2} SN  S_{2}  CS  COS HSNO HSO  HOS
HSOH H_{2}SO  HOSHO  HS_{2}  H_{2}S_{2} CS_{2}  HSOO H_{2}SO_{4}

2. táblázat. A második példában szereplő anyagok

ráadásul nem kell azzal foglalkoznunk, hogy a minimális megoldásokból előállítsuk a feltételeket kielégítő nem minimális megoldásokat is. A kiegészített LP relaxációs algoritmus két óra alatt meghatározta az 56677 megoldást.

Oxalát–perszulfát–ezüst oszcillátor [12]. Ebben a mechanizmusban 16 anyagfajta szerepel, ezek a 3. táblázatban láthatók. A felírható 152 egyenletrendszernek összesen 105 megoldása van. Az egyszerű rekurzív algoritmus 2,8, a Contejean–Devie-algoritmus 100, az LP relaxációs módszer 4,5 másodperc alatt találja meg a megoldásokat.

Mindhárom példa jól mutatja, hogy a Contejean–Devie módszer lényegesen kevésbé hatékony módja az elemi reakciók meghatározásának, mint a két másik ismertetett eljárás. Az egyszerű rekurzív eljárás a legjobb, amíg a felbontandó vektor és a megoldások száma kicsiny, ekkor azonban még az LP relaxációs eljárás is elég gyors. A nagyobb jobboldalú és sok megoldással rendelkező egyenletek esetén az LP relaxáció alkalmazása a legcélravezetőbb.

Ag^{+} Ag^{2+} H^{+}  SO_{4}^{−}
SO_{4}^{2−} S_{2}O_{8}^{2−} C_{2}O_{4}^{2−} Ag(C_{2}O_{4})^{−}
OH H_{2}O CO_{2}^{−}  O_{2}
HO_{2} H_{2}O_{2} O_{2}CO_{2}^{−} CO_{2}

3. táblázat. Az oxalát–perszulfát–ezüst oszcillátorban szereplő anyagok (sic!)

5. A felbontások előállítása

5.1. A keresési tér szűkítése

Nagyobb bruttó reakciók esetén – mint láttuk – a felbontásban résztvevő elemi reakciók száma igen nagy lehet, ilyenkor a felbontások keresése egy hatalmas térben történik, ezért nagyon időigényes lehet. Ebben a fejezetben megvizsgáljuk, hogyan egyszerűsíthetjük a keresést a potenciális megoldások terének szűkítésével.

5.1.1. Biztosan végbemenő reakciók

Bizonyos felbontó módszerek esetén gyorsíthatja a keresést, ha egyes reakciókról előre meg tudjuk állapítani, hogy legalább egyszer bekövetkeznek a bruttó reakció lefutása közben. Ezeket a lépéseket azonnal levonhatjuk a bruttó reakcióból. Az így kapott egyenletek megoldásai kisebbek az eredeti egyenleténél, ez különösen az első két algoritmus számára igen kedvező.

Természetesen annak eldöntése, hogy egy elemi reakció biztosan végbemegy-e, azaz egy particionáló vektor minden felbontásban szerepel-e, és ha igen, hányszor, ugyanúgy NP-nehéz, mint az eredeti feladat. (Könnyen végiggondolható, hogy ez és az eredeti feladat kölcsönösen visszavezethetők egymásra.) Ezért nem a pontos megoldásra törekszünk, csupán arra, hogy gyorsan alulról becsüljük az egyes reakciók multiplicitását az összes felbontásban. Ennek egy lehetséges módja a pontos megoldás meghatározását lehetővé tevő egészértékű programozási (ILP) feladat LP relaxációja, amely a következőképpen fogalmazható meg.

Jelöljük ismét a particionáló vektorokat v_{1}, v_{2}, …, v_{n}-nel, és legyenek ezek (a konkrét felbontás ismeretének hiányában ismeretlen) multiplicitásai x_{1}, …, x_{n}. A felbontandó vektor legyen ismét b. Ekkor a megoldandó lineáris program a k-adik reakció vizsgálatakor:

min x_{k}; ∑_{i} x_{i}v_{i} = b, ∀ i: x_{i} ≥ 0

A k-adik reakcióra kapott minimum nyilvánvalóan alsó becslés a reakció multiplicitására. A számításokat elvégző Mathematica program a beépített LinearProgramming függvénnyel könnyen elkészíthető.

A módszer legfőbb előnye, hogy gyors (mind elvi, mint gyakorlati szempontból.) Meg kell jegyeznünk azonban, hogy semmilyen javulást nem eredményezhet, ha olyan felbontó algoritmust használunk, amely maga is LP relaxációt használ, hiszen ezek ugyanezt az optimalizálást maguk is elvégzik a keresés első szakaszában. Az előző fejezetben ismertetett algoritmusok közül az első kettő alkalmazásakor azonban nagy hasznát vehetjük ennek az egyszerűsítésnek.

5.1.2. Nem végbemenő reakciók

Minden felbontó módszer esetén javulást eredményezhet egy olyan algoritmus, ami előzetesen elhagyja a particionáló vektorok közül azokat, amelyek igazolhatóan nem használhatók egyetlen felbontásban sem, ezzel csökkentve a változók számát.

A 3. fejezetben leírt algoritmusokkal minden olyan elemi reakciót megkaptunk, amelyek teljesítik a tömegmegmaradás törvényét, egyéb feltételt azonban nem szabtunk. Legelőször tehát érdemes mindazokat a reakciókat kihúznunk, amelyek kémiailag (például termodinamikai szempontból) értelmetlenek. Ebben bizonyos (tipikusan szerves) reakciók esetén számos mérési adat segít, máskor kevés az algoritmizálható kémiai elv. Akár arra is szükség lehet, hogy kézzel válogassuk ki az értelmetlen reakciókat. A Mathematica programozási környezete lehetővé teszi, hogy tetszőleges – a szóban forgó bruttó reakció esetében használható – algoritmussal kiegészítsük az általános használatra szánt rutinokat. A következőkben csak azzal foglalkozunk, milyen megoldásokat alkalmazhatunk akkor, ha a felbontandó bruttó reakció egyenletén, az (esetleg egy kémiai szűrésen már átesett) elemi reakciók listáján és a kezdeti anyagok listáján kívül semmilyen adatot nem tudunk (vagy szeretnénk) felhasználni.

LP-alapú megoldás. A minimális multiplicitás meghatározására alkalmazott lineáris programozás itt is alkalmazható, de másképpen, mint az előbb. A következőképpen járhatunk el [11].

1. [inicializálás] Legyen minden particionáló vektor jelöletlen.

2. [teszt] Keresünk egy „racionális felbontást”, azaz az (1) egyenlet egy megoldását a nemnegatív racionális számok felett, amelyben legalább egy jelöletlen vektor együtthatója pozitív. Ha nincs ilyen felbontás, megállunk.

3. [ciklus] A kapott felbontásban pozitív együtthatójú vektorokat megjelöljük, majd a második lépésre ugrunk.

Nyilvánvaló, hogy az eljárás véges sok lépésben véget ér, hiszen minden iterációban, amikor találunk megoldást (ezért nem állunk meg), a jelöletlen vektorok száma csökken. Az algoritmus végrehajtása után a jelölt vektorok azok, amelyek szerepelhetnek felbontásban.

Megjegyezzük, hogy a 2. pontban alkalmazott LP relaxáció (a megoldások egészértékűségi feltételének gyengítése) után kapott feladat kémiailag közelebb áll az eredeti feladathoz, mint matematikailag a kapott LP feladat megoldása az eredeti ILP-éhez. A racionális együtthatójú felbontásokból ugyanis a bruttó reakció valamely többszörösének egy felbontása adódik, ami az eredeti reakció többszöri lefutásának felel meg, így a megoldás kémiailag értelmes és hasznosítható.

Azt, hogy az újabb felbontások mindig tartalmazzanak jelöletlen vektort, úgy érhetjük el, hogy egy újabb egyenlőtlenséget veszünk a lineáris programhoz, ami ekkor a következő:

min c^{T}x; ∑_{i=1}^{n} x_{i}v_{i} = b, c^{T}x ≥ 1, ∀ i: x_{i} ≥ 0,  (3)

ahol az n dimenziós c vektor i-edik komponense 1, ha az i-edik vektor jelöletlen, és 0, ha jelölt.

Elvi szempontból teljesen mindegy, hogy a fenti lineáris programban mit optimalizálunk (akár egy konstans legkisebb értékét is kereshetnénk), a lényeges csupán az, hogy a feltételrendszer kielégíthető-e, ha kielégíthető, akkor az optimum létezzen, és hogy mutassunk egy felbontást, amely kielégíti azt. (A [11] cikk ezzel nem is foglalkozik.) A minimalizált kifejezésnek arra azonban van hatása, hogy az esetleg számos megoldás közül melyiket kapjuk meg. Ennek segítségével irányíthatjuk a keresést úgy, hogy (várhatóan) minél kevesebb iterációval érjünk az eljárás végére. Egy ennek megfelelő „mohó” heurisztikával kiegészített algoritmus a fenti helyett a jelölt vektorok multiplicitásának összegét minimalizálná abban a reményben, hogy így a kapott felbontásban sok jelöletlen vektor szerepel majd. Vegyük észre, hogy a fenti módszer éppen ez ellen dolgozik. (Ez legfeljebb abban az esetben lehet hasznos, ha sok felbontást is szeretnénk kapni, ezért kifejezetten sokat szeretnénk iterálni.) A különbség nagyobb egyenletek esetén a gyakorlatban is tettenérhető. (Ld. az 5.4. szakaszban.) Az algoritmus első, azaz a (3) lineáris programmal dolgozó változatának Mathematica kódja az A.4. szakaszban olvasható. Mást minimalizáló változatai ebből könnyen megkaphatók.

Volpert-indexelés. A fenti algoritmus kizárólag matematikai megfontolásokat használ, nem veszi figyelembe, hogy a kapott felbontások kémiailag értelmesek-e. Kémiailag értelmetlen egy felbontás akkor, ha valamelyik benne szereplő anyag reakciója a bruttó reakció lezajlása alatt mindvégig zérus. Az anyagok koncentrációjának időbeli változását a felbontás alapján felírható differenciálegyenlet-rendszer határozza meg. Volpert tétele szerint mindazon anyagok koncentrációja zérus, amelyek a következő indexelési algoritmus során végtelen indexet kapnak [1]:

1. [inicializálás] A kezdeti anyagok indexe legyen zérus.

2. [reakciók indexelése] Mindazon indexeletlen reakciókat, amelyeknek a kiindulási anyagai indexeltek, indexeljük a legnagyobb felhasznált indexszel.

3. [anyagok indexelése] Mindazon indexeletlen anyagokat, amelyek valamely indexelt reakció termékei, indexeljük a legnagyobb felhasznált indexnél eggyel nagyobb indexszel.

4. [ciklus] Ha van még a 2. pontban indexelhető reakció, lépjünk a második pontra.

5. [befejezés] Minden indexeletlen anyag és reakció indexe legyen végtelen.

A Volpert-indexelés tehát azt a páros gráfot indexeli, aminek a két színosztálya a felbontásban résztvevő reakciók és az anyagok. Az eljárást egyebek között arra használják (ld. például [1, 11]), hogy a talált felbontások kémiai helyességét ellenőrizzék. Némi módosítással azonban ugyanezt az eljárást a felbontások keresése előtti szűrésre is használhatjuk: a kezdeti anyagok listája legyen ugyanaz, mint előbb, de most vegyük figyelembe az összes elemi reakciót, és az így kapott gráfot indexeljük. Azok a reakciók, amelyek indexe végtelen, nem szerepelhetnek egyetlen kémiailag értelmes felbontásban sem. Hasonlóan: a végtelen indexet kapott anyagok semmilyen kémiailag értelmes felbontás során nem állnak elő. Ezeket a reakciókat és anyagokat a későbbi vizsgálat során figyelmen kívül hagyhatjuk.

A Volpert-indexelés egy igen gyors Mathematica implementációja olvasható az A.5. szakaszban. Érdemes megemlíteni, hogy az előbbiek alapján semmi szükség arra, hogy az indexelés pozitív egész számokkal és végtelennel történjen, a véges indexű anyagokat és reakciókat egyszerűen jelöltnek, míg a végtelen indexűeket jelöletlennek nevezve ugyanazt kapnánk, mint az LP relaxációs megoldásnál: a jelölt anyagok – és csakis azok – szerepelhetnek felbontásban. Mint azt már említettük, a Volpert-indexelést más feladatok során is használják, ezek némelyike igényli a jelölt anyagok ilyen megkülönböztetését (a már idézett [1] cikk is így közli az algoritmust), ezért választottuk mi is ezt a leírást.

A két módszer összehasonlítása. A két módszer hatékonyságát nehéz összevetni, és nem is feltétlenül érdemes, ugyanis csak a céljuk közös – a keresési tér méretének csökkentése a felesleges particionáló vektorok elhagyásával. A két algoritmus ugyanakkor teljesen különböző elven alapszik: míg az egyik tisztán matematikai eszközökkel válogat a particionáló vektorok között, a másik olyan kémiai megfontolásokat használ, amelyek az elsőbe nem építhetők be.

Ez nem jelenti azt, hogy a relaxációs megoldás nem kezelhetne kémiai feltételeket. Valójában tetszőleges olyan feltétel felhasználására képessé tehető, ami lineáris egyenlet vagy egyenlőtlenség formájában felírható. Ilyen feltételek például: bizonyos reakciók csak korlátos sokszor hajtódhatnak végre, az összes lépések száma ne lépjen túl egy meghatározott értéket, egyes reakciók kizárják egymást stb. A Volpert-indexelés eredeti céljánál többet nem érhet el, gyakorlati példák esetén azonban gyakorta több reakciót szűr ki, mint a (kiegészítő kémiai feltételek nélkül használt) relaxációs algoritmus.

A Volpert-indexelés rendkívül gyors (lineáris időben megvalósítható), de a relaxációs algoritmus is elég hatékony, a jól implementált lineáris programozási rutinnak köszönhetően. „Mellékhatásként” ráadásul a bruttó reakció számos felbontását is előállítja, ami a későbbiekben (a bruttó reakció felbontásakor) hasznos lehet. Végeredményben tehát a két algoritmust érdemes egyszerre használni. Alkalmazási példákat az 5.4. szakaszban mutatunk.

5.2. A legrövidebb felbontások előállítása

Miután az elemi reakciók listáját a lehető legjobban megkurtítottuk, következhet a bruttó reakció felbontásainak megkeresése. Elsőként arra keresünk választ, hogy hány lépésből állnak a legegyszerűbb (legkevesebb lépést tartalmazó) felbontások. Ez a feladat még akkor is könnyen megoldható lehet, ha az összes (a ⋗ reláció szerint) minimális felbontás előállítása reménytelen.

A legegyszerűbb, ha először a legfeljebb egy lépésből álló megoldásokat keressük, majd a legfeljebb két lépésből állókat, és így tovább, mígnem találunk egyet. Sok felesleges munkától szabadulhatunk meg, ha lineáris programozással becsüljük a megoldásvektorok minimális súlyát – ehhez az ismeretlenek összegét kell minializálni az (1) egyenlet, mint kényszerfeltétel mellett. Ez szerencsés esetben akár egy minimális megoldást is adhat.

A legfeljebb N lépésből álló megoldásokat a Contejean–Devie-algoritmus 3.2. szakaszbeli kiterjesztésével kereshetjük a legkönnyebben, hiszen annak opcionális paramétere a keresett megoldások maximális mérete. Azonban nem biztos, hogy ez a módszer a leghatékonyabb, különösen, ha a legrövidebb megoldások nagyon nagyok. A másik két algoritmus is használható a következő módon: vegyünk hozzá a particionáló vektorokhoz egy-egy további (d + 1-edik) komponenst, amelynek értéke legyen 1, és egészítsük ki a particionáló vektorokat egy (\underbrace{0, …, 0}_{d}, 1)^{T} vektorral. Végül a felbontandó vektort

is toldjuk meg egy további, N értékű komponenssel. Könnyen látható, hogy az így kapott particionálási feladat megoldásai (elhagyva az utolsó vektor együtthatóját) éppen az eredeti feladat azon megoldásai, amelyek legfeljebb N méretűek. Ezzel a kiegészítéssel a legrövidebb felbontások kiszámítására mindhárom algoritmus közvetlenül alkalmazható.

5.3. Az összes felbontás előállítása

Az összes (minimális) felbontás meghatározása hasonló az előző feladathoz. A Contejean–Devie-algoritmus közvetlenül alkalmazható, azonban nagy bruttó reakciók esetén nem várható tőle gyors eredmény. A másik két algoritmus az előzőekhez hasonló kiterjesztéssel alkalmazható. Ha az előbbi konstrukcióból elhagyjuk a (\underbrace{0, …, 0}_{d}, 1)^{T} vektort, akkor azokat a megoldásokat kaphatjuk

meg, amelyek pontosan N lépésből állnak. Az N szám fokozatos növelésével újabb és újabb megoldásokat generálhatunk, amelyeket összevetve a kisebb megoldásokkal eldönthetjük, hogy minimálisak-e. Az egyetlen probléma, hogy nem tudhatjuk, N mekkora értékénél állhatunk meg, ha nem tudjuk becsülni a vizsgált egyenlet minimális megoldásainak a méretét.

5.4. Példák

Permanganát–oxálsav reakció [11]. Az előző fejezetben előállítottuk azt az 1022 elemi reakciót, amely esetleg végbemehet a bruttó reakció során. Ezek egy része kémiailag értelmetlen. Megfelelő adatok híján termodinamikai szempontok alapján nem tudtunk számítógép segítéségével szűrni, ehelyett egy vegyész közreműködésével húztunk ki a reakciók közül 349-et. A fennmaradó 673 termodinamikailag lehetséges reakciót már az 5.1. szakaszbeli módszerekkel vizsgáltuk.

Elsőként megállapítottuk, hogy két reakció minden felbontás során végbemegy. A H_{2}C_{2}O_{4} + MnO_{2} → [MnO_{2}, H_{2}C_{2}O_{4}] reakció legalább négyszer, a 2 CO_{2}^{−} → C_{2}O_{4}^{2−} reakció pedig legalább egyszer végrehajtódik. Ezzel a megoldások méretét öttel csökkenthetjük, ami az első két módszer esetén jelentős előny. A keresés alig 2 percet vett igénybe.

Mindkét ismertetett módszerrel kerestünk nem végbemenő reakciókat és elő nem állítható anyagokat. Az első (LP-n alapuló) algoritmus további kémiai feltételek nélkül 314 lépést talált, amelyek semmiképpen nem mehetnek végbe, kiegészítő feltételek alkalmazásával még kettőt. Megvizsgáltuk, milyen hatása van annak, ha a (3) LP feladatban megváltoztatjuk a minimalizálandó kifejezést. A várakozásoknak megfelelően kevesebb iterációra volt szükség, amikor a c^{T}x kifejezés helyett az (1 − c^{T})x kifejezést minimalizáltuk minden lépésben. (313 helyett 280 LP feladatot kell megoldanunk.) Végül 359 (illetve 357) reakció marad. Az eljárások futásideje 5–7 perc volt.

Az így kapott reakciókat beindexeltük Volpert algoritmusával, különböző kezdetianyag-kombinációk mellett. A [11]-ben szereplő öt kezdeti anyag feltételezése esetén éppúgy 297 reakció (és csak 16 anyagfajta) indexelhető be, mintha az összes nem komplex anyagot a kezdeti anyagok közé választjuk. Sőt, pontosan ugyanazokat a reakciókat kapjuk akkor is, ha az első szűrés után megmaradt 673 reakciólépést indexeljük be. Az indexelés minden esetben a másodperc törtrésze alatt végbement.

A 297 felhasználható reakcióval megkezdhettük a bruttó reakció összes felbontásának a meghatározását. A legrövidebb lépéssorozatok meghatározásához először lineáris programozással becsültük a megoldások méretét. Az egyes reakciólépések együtthatóit jelentő változók összegének minimumára (az (1) egyenlet, mint kényszerfeltétel mellett) 15 adódott, a talált optimum ráadásul minden változóhoz egész értéket rendelt, így azonnal egy minimális felbontást is megkaptunk:

4 × H_{2}C_{2}O_{4} + MnO_{2}  → [MnO_{2}, H_{2}C_{2}O_{4}]
2 × H_{2}C_{2}O_{4} + [MnO_{2}, H_{2}C_{2}O_{4}] → 2CO_{2} + 2H_{2}O + MnC_{2}O_{4}
1 × H^{+} + HC_{2}O_{4}^{−}  → H_{2}C_{2}O_{4}
1 × C_{2}O_{4}^{2−} + H^{+}  → HC_{2}O_{4}^{−}
2 × MnC_{2}O_{4} + MnO_{4}^{−}  → CO_{2}^{−} + CO_{2} + 2MnO_{2}
1 × 2CO_{2}^{−}  → C_{2}O_{4}^{2−}
2 × H^{+} + [MnO_{2}, H_{2}C_{2}O_{4}]  → [H^{+}, MnO_{2}, H_{2}C_{2}O_{4}]^{+}
2 × H^{+} + [H^{+}, MnO_{2}, H_{2}C_{2}O_{4}]^{+} → 2CO_{2} + 2H_{2}O + Mn^{2+}

Érdemes megjegyeznünk, hogy az LP relaxációs algoritmus azon kiterjesztése, amellyel a legfeljebb N lépéses megoldásokat megkereshetjük, fél másodperc alatt felismeri, hogy nincsen legfeljebb 14 lépéses megoldás. A Contejean–Devie-algoritmus tárhatékony implementációjának ehhez 25 percre volt szüksége (a ciklus többszázezerszeri lefutásával) még azután is, hogy a biztosan végbemenő öt lépést kivontuk az egyenletből.

A többi minimális felbontást is megkaphatjuk, ha a Contejean–Deviealgoritmus kiegészített változatát használjuk a legfeljebb 15 lépéses megoldások előállítására, vagy ha a másik két módszer valamelyikének módosított változatát alkalmazzuk a pontosan 15 lépéses megoldások megkeresésére. A biztosan végbemenő öt lépést itt is kivonhatjuk az egyenletből, ekkor természetesen 10 lépéses megoldásokat keresünk. 11 különböző 15 lépéses megoldás van, ezek természetesen mind minimálisak. Előállításuk az LP relaxációs algoritmussal (ami ilyen méretű feladatoknál már lényegesen gyorsabb a másik két eljárásnál) másfél percig tart. További felbontásokat állíthatunk elő, ha egyesével növelve a lépésszámot előállítjuk a többlépéses megoldásokat. Itt már elérjük az LP relaxációs módszer korlátait – a 14 lépéses felbontások előállítására nem volt elég egy óra sem.

A megoldás utolsó lépéseként a megtalált felbontások közül Volpert-indexe-léssel ki kell választanunk a kémiailag értelmeseket.

6. Összefoglalás, további tervek

Dolgozatunk célja az volt, hogy ismertessük a bruttó reakciók felbontásának kémiai feladatát, a megoldása során felmerülő matematikai problémákat, és hogy nagyobb méretű gyakorlati példák megoldása esetén is használható algoritmusokat mutassunk be, azok Mathematica implementációjával együtt.

A feladat első két része lényegében azonos gyökerű, a megoldáshoz lineáris diofantoszi egyenletrendszerek összes minimális megoldásának felsorolása szükséges. Megállapítottuk, hogy a diofantoszi egyenletekkel foglalkozó (igen bőséges) irodalomban meglepően kevés figyelmet szentelnek ennek a feladatnak. Igazoltuk, hogy a megoldhatóság eldöntése NP-teljes, így az összes minimális megoldás előállítására sem várható gyors algoritmus. Bemutattunk két saját, és egy, az irodalomban fellelhető eljárást a megoldásra, és alkalmazásukat gyakorlati kémiai feladatokon.

Megvizsgáltunk néhány, a keresési tér csökkentésére szolgáló módszert is. Ezek egyike arra irányul, hogy a kémiailag értelmetlen megoldások számát csökkentsük, míg két másik módszerrel tetszőleges – nem csak reakciófelbontásból származó – egyenletek megoldását gyorsíthatjuk a megoldások méretének csökkentése, illetve a változók számának csökkentése révén, polinomidejű preprocesszálással.

A függelékben megadtuk az algoritmusok Mathematica implementációit is. A fejlesztési tervek között szerepel az elkészített függvények egy olyan Mathematica csomagba szerkesztése, amellyel a reakciók felbontásának összes vizsgált részfeladata automatizálható.

A. Az algoritmusok implementációja

Minden ismertetett algoritmust Mathematicában implementáltunk, ügyelve arra, hogy az elkészített függvények robosztusak és nagy (de még elvileg kezelhető) méretű bemenet esetetén is elég gyorsak legyenek a lehetőségekhez képest. Helykímélés céljából a programoknak egy minimális hibakezeléssel közölt változatát mutatjuk be. A függvényeknek mindig csak egy változatát közöljük. Amint arra a dolgozat megfelelő részeiben utaltunk, bizonyos problémák esetén szükség lehet könnyen elvégezhető, apróbb módosításokra.

A közölt kódok egy része lényegében bármilyen Mathematica kernellel futtatható, néhányukhoz Mathematica 5.0 szükséges.

A.1. A naiv rekurzív algoritmus

SimplePartitions[M_?MatrixQ, b_?VectorQ] :=
    Block[{$IterationLimit = $RecursionLimit = Infinity)},
            SimplePartitions[M, b, 1]
    ]
SimplePartitions[M_?MatrixQ, b_?VectorQ, index_Integer] :=
    Which[
        Max[Abs[b]] == 0,
            Array[0&, {1, Length[M]}],
        SmallerComponentQ[b, M�index�,
                        NonNegativeColumns[Drop[M, index - 1]]],
            If[index == Length[M], {}, SimplePartitions[M, b, index + 1]],
        index == Length[M],
            IncreaseColumn[SimplePartitions[M, b - M�index�, index], -1],
        True,
            Join[
                IncreaseColumn[SimplePartitions[M, b - M�index�, index], index],
                SimplePartitions[M, b, index + 1]
            ]
    ]

Az eljárásunk három segédfüggvényt használ. Az első meghatározza, hogy egy megadott vektornak van-e a másikénál kisebb komponense (paraméterezve avval, hogy hányadik komponenseket vesszük figyelembe.)

A AZ ALGORITMUSOK IMPLEMENTÁCIÓJA 31

SmallerComponentQ[v1_?VectorQ, v2_?VectorQ, selected_?VectorQ] :=
    Max[(v2 - v1)�selected�] > 0

A második egy mátrix kijelölt oszlopának minden eleméhez egyet ad.

IncreaseColumn[M_?MatrixQ, col_Integer] :=
    Module[{n = Transpose[M]}, n�col� += 1; Transpose[n]]
IncreaseColumn[{}, _] := {}

A harmadik egy mátrix csupa nemnegatív elemeit tartalmazó oszlopainak indexét adja vissza.

NonNegativeColumns[M_?MatrixQ] :=
    Flatten[Position[And @@ #& /@ NonNegative[Transpose[M]], True]]

A particionáló vektorok jól rendezettségét a következő függvényekkel érhetjük el:

RightComponents[M_?MatrixQ] :=
    Flatten[Position[Transpose[M],
        {___?NonNegative, _?Positive, ___?NonNegative} |
        {___?NonPositive, _?Negative, ___?NonPositive}
    ]]
Orderable[M_?MatrixQ] :=
    With[{rc = RightComponents[M]},
        Length[rc] > 0 && Orderable[Select[M, #�First[rc]� === 0&]]
    ]
Orderable[{}] = True;
RightOrder[M_?Orderable] :=
    With[{rc = First[RightComponents[M]]},
        Join[
            Select[M, #�rc� =!= 0&],
            RightOrder[Select[M, #�rc� === 0&]]
        ]
    ]
RightOrder[{}] = {};
RightOrder[_] := {}

A.2. A Contejean–Devie-algoritmus

Az alábbi program az algoritmus javított változatát valósítja meg, a stepno paramétert az alapértelmezés szerinti (végtelen) értékén hagyva az eredeti változatot kapjuk meg.

ContejeanDeviePartitions[M_?MatrixQ, b_?VectorQ, stepno_:Infinity] :=
Module[{vecs = Append[M, -b],
im = IdentityMatrix[Length[M] + 1],
sols = {},
stack, szorzat, mo, hatra, valhatok, jok},
stack = Array[{M�#� - b, im�#� + im�Length[vecs]�, stepno, Range[#]}&, Length[M]];
While[Length[stack] =!= 0,
{szorzat, mo, hatra, valhatok} = First[stack];
Which[
    hatra === 0 || !MinimalQ[mo, sols],
        stack = Rest[stack],
    ZeroVectorQ[szorzat],
        PrependTo[sols, mo]; stack = Rest[stack],
    True,
        jok = Select[valhatok, szorzat.vecs�#� ≤ -(szorzat.szorzat/hatra) &&
                               MinimalQ[mo + im�#�, sols]&];
        stack = Join[
            Map[Function[x, {szorzat + vecs�x�, mo + im�x�, hatra - 1,
                             Complement[valhatok, Select[jok, # > x&]]}],
                jok],
            Rest[stack]
        ]
]
];
Most /@ sols
]

A felhasznált segédfüggvények:

ZeroVectorQ[v_List] := Max[Abs[v]] == 0
MinimalQ[v_List, solutions_List] := Negative[Max[ Min[v-#]& /@ solutions ]]

A.3. Az LP relaxációs felbontóalgoritmus

Az alábbi függvényt használjuk az általános, valós megoldás meghatározására.

UnrestrictedSolutions[m_?MatrixQ, b_List] :=
Block[{$Messages, ls, ns},
Check[
    ls = LinearSolve[m, b],
    Return[{}],
    LinearSolve::"nosol"
];
If[(ns = NullSpace[m]) === {},
   ls,
   ls + Transpose[ns].x[#]& /@ Range[Length[ns]])
]
]

A particionáló rutin a következő.

LPPartitions[m_?MatrixQ, v_List] :=
Block[{generalsol, vars, sols},
If[(generalsol = UnrestrictedSolutions[Transpose[m], v]) === {}, Return[{}] ];
vars = Union[Cases[generalsol, x[_], Infinity]];
sols = FindAllSolutions[
            generalsol,
            vars,
            First[Cases[
                generalsol,
                a_. + b_.vars�1� /; NumberQ[a] → 1 - FractionalPart[a]
            ]],
            Max[Cases[
                generalsol,
                a_. + b_.vars�1� /; NumberQ[a] → 1/b
            ]]
];
generalsol /. Thread[vars → #]& /@ sols
]

A függvény előállítja az általános megoldást paraméteres formában, majd megkeresi a paraméterek összes megfelelő értékét a FindAllSolutions hívással. A FindAllSolutions függvény határozza meg a paraméterekre adható alsó és felső korlátokat lineáris programozás segítségével. Ehhez használni fogjuk a következő segédfüggvényeket.

MyFloor[x_, prec_]   := Floor[x/prec]*prec
MyCeiling[x_, prec_] := Ceiling[x/prec]*prec

Végül a FindAllSolutions kódja:

FindAllSolutions[conlist_List, var_, offset_?NumberQ, step_?NumberQ] :=
Block[{$Messages, lower, upper},
Check[
    lower =
        MyCeiling[First@Minimize[
                            {var, Thread[conlist ≥ 0]},
                            {var}] - offset, step] + offset,
        Return[{}],
        Minimize::natt
];
Check[
    upper =
        MyFloor[First@Maximize[
                            {var, Thread[conlist ≥ 0]},
                            {var}] - offset, step] + offset,
        Return[{}],
        Maximize::natt
];
If[lower > upper, {},
   Transpose[{Select[Range[lower, upper, step],
                     And @@ IntegerQ /@ (conlist /. var → #)& ]}]
]
]
FindAllSolutions[conlist_List, vars_List, offset_?NumberQ, step_?NumberQ] :=
Block[{$Messages, var = First[vars], lower, upper, i, newcons, a, b},

Check[
lower = MyCeiling[First@Minimize[
                            {var, Thread[conlist ≥ 0]},
                            vars] - offset, step] + offset,
Return[{}],
Minimize::natt
];
Check[
upper = MyFloor[First@Maximize[
                            {var, Thread[conlist ≥ 0]},
                            vars] - offset, step] + offset,
Return[{}],
Maximize::natt
];
If[lower > upper, {},
Join @@ Table[
    newcons = conlist /. var → i;
    If[! And @@ IntegerQ /@ Select[newcons, NumberQ],
       {},
       Prepend[#, i]& /@
         FindAllSolutions[
            newcons,
            Rest[vars],
            First[Cases[
                newcons,
                a_. + b_.vars�2� /; NumberQ[a] → 1 - FractionalPart[a]
            ]],
            1 / Min[Abs[Cases[
                    newcons,
                    a_. + b_.vars�2� /; NumberQ[a] → b
                ]]]
         ]
    ],
    {i, lower, upper, step}]
]
]

A.4. Felbontásban nem szereplő reakciók kiszűrése LP relaxációval

A dolgozatban említettük, hogy a program több változatban elkészíthető attól függően, hogy milyen kémiai feltételeket kívánunk figyelembe venni, illetve hogy a felbontások keresésére használt LP feladatban mit minimimalizálunk. Az alábbiakban a legegyszerűbb, semmilyen kémiai feltételt nem használó, a (3) LP feladatot megoldó változat olvasható.

Kizarhato[reakciok_?MatrixQ, brutto_?VectorQ] :=
    Block[{$Messages, egyesek, pozitivak = {}, d1 = Length[reakciok], cm, poz},
        egyesek = Array[1&, d1];
        While[True,
            If[Total[egyesek] === 0, Break[]];
            lp =
                Check[
                    LinearProgramming[egyesek,
                                      Prepend[Transpose[reakciok], egyesek],
                                      Prepend[{#, 0}& /@ brutto, {1, 1}]
                    ],
                    Break[],
                    LinearProgramming::"lpsnf"
                 ];
            poz = Complement[Range[d1], Flatten[Position[lp, 0]]];
            AppendTo[pozitivak, {poz, Select[lp, Positive]}];
            egyesek�poz� = 0;
        ];
        pozitivak
    ]

A függvény paraméterezése értelemszerű, a kimenet pedig két listából álló párok listája, ahol minden pár egy felbontást ad meg: az első tagja a pozitív együtthatójú reakciók indexlistája, második tagja pedig az együtthatólista. Így a nem kizárható reakciók indexeinek a listája közvetlenül megkapható a kimenetből. Hogy mégsem csak ezt adjuk vissza (ehhez csupán az utolsó előtti sort kellene módosítani), annak az az oka, hogy a kiszámított felbontások kémiai vizsgálódások során jól hasznosíthatók.

A.5. Volpert-indexelés

VolpertIndexing[A_?MatrixQ, B_?MatrixQ, M0_List] :=
    Module[{iterationlist =
                FixedPointList[
                    {#�1�(1 - Sign[B.(1 - #�2�)]), #�2�}& @
                                    {#�1�, #�2� Sign[#�1�.A]} &,
                    {MapAt[# - 1&, Array[1&, Length[A]], Transpose[{M0}]],
                     Array[1&, Length[First[A]]]}
                ]
           },
           MapAt[# - 1&, (Last[iterationlist] /. (1 → Infinity))
                                                  + Total[iterationlist], 2]
    ]

A program első két paramétere a reakciókat leíró mátrixok. Az A mátrix iedik sorának j-edik eleme k, ha az i-edik reakció bal oldalán a j-edik anyag k együtthatóval szerepel. Hasonlóképpen kapjuk a reakciók jobb oldalaiból a B mátrixot. Végül a harmadik paraméter a kezdeti anyagok indexeinek listája. A kimenet a reakciók és az anyagok indexeiből álló listapár.

Készült a T037491 sz. OTKA pályázat részleges támogatásával.

A szerző köszönetet mond Daniel Lichtblaunak (Wolfram Research) és Vizvári Bélának (ELTE, Operációkutatási Tanszék) értékes segítségükért.

Hivatkozások

[1] Volpert, A. I., Chudayev, S. I., Analysis in classes of discontinuous functions and the equations of mathematical physics, Nauka, Moszkva, 1975.

[2] Tomlin, A. S., Turányi, T., Pilling, M. J., Mathematical tools for the construction, investigation and reduction of combustion mechanisms, In: Comprehensive Chemical Kinetics, Ed. Compton, R. G. and Hancock, G., Vol. 35, Elsevier, 1997.

[3] Contejean, E., Devie, H., „An efficient incremental algorithm for solving systems of linear Diophantine equations”, Information and Computation, 113 (1994), 143–172.

[4] de Weger, B. M. M., Algorithms for Diophantine equations, CWI, Amsterdam, 1989.

[5] Smart, N., The algorithmic resolution of Diophantine equations, Cambridge UP, Cambridge, 1988.

[6] Lichtblau, D., „Revisiting strong Gröbner bases over Euclidean domains”, Wolfram Research, Inc., Champaign, 2001.

[7] Szalkai I., „A new general algorithmic method in reaction syntheses using linear algebra”, Journal of Mathematical Chemistry, 28 (2000)

[8] Rónyai L., Ivanyos G., Szabó R., Algoritmusok, Typotex, Budapest, 2000.

[9] Ahlgren, S., Ono, K.: „Addition and Counting: the Arithmetic of Partitions”, Notices of the AMS, 48 (9) (2001), 978–984.

[10] Vizvári B., Egészértékű programozás, Tankönyvkiadó, Budapest, 1992.

[11] Kovács K., Riedel M, Tóth J., Vizvári B., „Decomposition of the permanganate/oxalic acid overall reaction to elementary steps based on integer programming theory”, előkészületben

[12] Clarke, B. L., „Stoichiometric network analysis of the oxalate–persulfate– silver oscillator”, Journal of Chemical Physics, 97 (4) (1992), 2459–2472.