Implementacija

Da bismo realizovali predloženi model prigušenih oscilacija u jednoj dimenziji opisan običnom diferencijalnom jednačinom u sekciji Oscilator, umesto biblioteke SCIANN koristićemo nešto noviju biblioteku DeepXDE Lu et al. [LMMK21]. Ideja je da ovim materijalom pokrijemo sve značajne softverske okvire koji podržavaju metodologiju NMPFZ u trenutku pisanja (2022-2023). Kako autori kažu, DeepXDE biblioteka je dizajnirana da služi i kao obrazovno sredstvo koje će se koristiti u visokom školstvu i kao istraživački alat za rešavanje problema u kompjuterskim naukama i inženjerstvu. Konkretno, DeepXDE može da rešava probleme za unapred date početne i granične uslove, kao i inverzne probleme uz neka dodatna merenja. Podržava domene složene geometrije i omogućava da korisnički kod bude kompaktan, veoma sličan matematičkoj formulaciji.

U odnosu na SCIANN, pristup DeepXDE podrazumeva nešto viši stepen apstrakcije, pa time i lakoće upotrebe. Na primer, komplikacije oko postavljanja kolokacionih tačaka u kojima važe granični i početni uslovi, kao i veličina batch-a, ne potpadaju pod brigu korisnika, već se sama biblioteka stara o tome da specificirani broj tačaka podleže graničnim uslovima. Taj nivo apstrakcije u nekim specifičnim slučajevima može predstavljati prepreku, ali u velikoj većini slučajeva doprinosi jasnijem definisanju problema.

Na sledećem listingu dati su značajni delovi implementacije:

Listing 4 Rešenje problema prigušenih oscilacija u 1D korišćenjem DeepXDE biblioteke

import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt

m = 1
mu = 0.1
k = 2

# Pocetni uslovi
x0, v0 = 0, 2

delta = mu / (2*m)
w0 = np.sqrt(k/m)

# Da li je tacka blizu t=0 (provera pocetnog uslova)
def boundary_l(t, on_initial):
    return on_initial and np.isclose(t[0], 0)

# Jednacina ODE
def Ode(t,x):
        dxdt = dde.grad.jacobian(x,t)
        dxdtt = dde.grad.hessian(x,t)
        return m * dxdtt + mu * dxdt + k*x

# x(0)=x0
def bc_func1(inputs, outputs, X):
    return outputs - x0

# x'(0)=v0
def bc_func2(inputs, outputs, X):
    return dde.grad.jacobian(outputs, inputs, i=0,j=None) - v0

# Resava se na domenu t=(0,10)
interval = dde.geometry.TimeDomain(0, 10)

# Pocetni uslovi
ic1 = dde.icbc.OperatorBC(interval, bc_func1, boundary_l)
ic2 = dde.icbc.OperatorBC(interval, bc_func2, boundary_l)

# Definisanje problema, granicnih uslova, broja kolokacionih tacaka
data = dde.data.TimePDE(interval, Ode, [ic1, ic2], 100, 20, solution=func, num_test=100)

layers = [1] + [30] * 2 + [1]
activation = "tanh"
init = "Glorot uniform"
net = dde.nn.FNN(layers, activation, init)

model = dde.Model(data, net)

model.compile("adam", lr=.001, loss_weights=[0.01, 1, 1], metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=30000)

T = np.linspace(0, 10, 100).reshape(100,1)
x_pred = model.predict(T)

Na početku se definišu importi i konstante problema, kao i dodatne konstante delta i w0 koje smo koristili prilikom izvođenja analitičkog rešenja. Funkciju

def boundary_l(t, on_initial):
    return on_initial and np.isclose(t[0], 0)

ćemo iskoristiti za testiranje da li je data kolokaciona tačka blizu tačke t=0, tj. da li za datu tačku važi početni uslov. Upotrebićemo je pri formiranju početnih uslova za poziciju i brzinu. Dalje, kao što naziv sugeriše, naredna funkcija

def Ode(t,x):
        dxdt = dde.grad.jacobian(x,t)
        dxdtt = dde.grad.hessian(x,t)
        return m * dxdtt + mu * dxdt + k*x

predstavlja postavku problema u svom izvornom obliku jednačine (20). Ovde je očigledna jedna od glavnih prednosti NMPFZ, tj. da neke dodatne transformacije u integracionim tačkama nisu potrebne, već samo postavka diferencijalne jednačine u vidu funkcije gubitka i graničnih uslova u istom obliku. Uslužne metode dde.grad.jacobian i dde.grad.hessian vraćaju prve, odnosno druge izvode po ulaznim varijablama primenjujući tzv. automatsku diferencijaciju. Podrazumevano se u pozadini koristi Tensorflow za tenzorske operacije niskog nivoa.

Postavka dva početna uslova u formi funkcije gubitka data je u sledeće dve metode, za koordinatu i brzinu respektivno:

def bc_func1(inputs, outputs, X):
    return outputs - x0

def bc_func2(inputs, outputs, X):
    return dde.grad.jacobian(outputs, inputs, i=0,j=None) - v0

Nakon postavke jednodimenzionog vremenskog domena u kome se problem rešava:

interval = dde.geometry.TimeDomain(0, 10)

možemo da formiramo i objekte graničnih uslova kombinujući funkcije gubitka sa funkcijom lokacije boundary_l():

ic1 = dde.icbc.OperatorBC(interval, bc_func1, boundary_l)
ic2 = dde.icbc.OperatorBC(interval, bc_func2, boundary_l)

Sada imamo sve elemente da formiramo objekat problema koji rešavamo. Ovde ćemo to učiniti metodom dde.data.TimePDE za vremenski zavisne probleme:

data = dde.data.TimePDE(interval, Ode, [ic1, ic2], 100, 20, solution=func, num_test=100)

Specificiramo redom računski domen, osnovnu jednačinu, listu graničnih uslova, broj kolokacionih tačaka za osnovni domen (100), broj kolokacionih tačaka za granične uslove (20), egzaktno rešenje (ako postoji) i broj testnih tačaka (za poređenje sa egzaktnim rešenjem). U ovom primeru ćemo ignorisati egzaktno rešenje. Odmah se vidi razlika u postavci u odnosu na SCIANN pristup u načinu navođenja kolokacionih tačaka. Naime, kod DeepXDE kolokacione tačke se ne generišu manuelno, već se prepušta biblioteci da to uradi za nas, što ukazuje na jedan viši nivo apstrakcije.

Naredne linije koda konstruišu neuronsku mrežu koja će se koristiti kao aproksimacija problema, sa svim svojim hiper-parametrima:

layers = [1] + [30] * 2 + [1]
activation = "tanh"
init = "Glorot uniform"
net = dde.nn.FNN(layers, activation, init)

Naša mreža ima jedan ulaz, jedan izlaz i dva skrivena sloja od po 30 neurona, sa aktivacijom skrivenih slojeva u vidu tanh funkcije i odgovarajućom inicijalizacijom. Na kraju, možemo da krenemo u obučavanje, kada spojivši problem i generisanu neuronsku mrežu formiramo model:

model = dde.Model(data, net)
model.compile("adam", lr=.001, loss_weights=[0.01, 1, 1], metrics=["l2 relative error"])
losshistory, train_state = model.train(epochs=30000)

Ovo su standardne metode koje se široko koriste u oblasti dubokog učenja, pa je dovoljno samo pomenuti da se navodi algoritam optimizacije, brzina učenja i način proračuna greške koja upravlja ovim procesom. Specifičnost za NMPFZ je što ovde listom loss_weights možemo i da „ponderišemo” težine osnovne diferencijalne jednačine, prvog i drugog graničnog uslova, respektivno. U narednoj sekciji Rešenja direktnog problema ćemo razmotriti rešenja za sva tri slučaja prigušenog oscilovanja.