Inverzni problem

Kao što smo već u više navrata naglasili, prava snaga NMPFZ dolazi do izražaja kod problema u kojima je potrebno identifikovati nepoznate parametre. Razlog tome je što NMPFZ metoda inverzne probleme tretira na isti način kao direktne. Važno je samo da se problem postavi ne dozvoljavajući višeznačnost, što ćemo demonstrirati na sledećem primeru. Zamislimo eksperiment postavljen kao na Sl. 15 u kome nam je poznat samo koeficijent trenja \(\mu=0,6\), a ne znamo ni masu kuglice \(m\) ni koeficijent elastičnosti opruge \(k\). Zadatak nam je da identifikujemo ova dva parametra tako što ćemo pustiti kuglicu da osciluje i na osnovu njenog kretanja dati NMPFZ da odredi elemente koji nedostaju.

Recimo da smo vizuelno utvrdili da je u pitanju podprigušeni slučaj. Prva strategija merenja koja bi mogla da bude realno izvodljiva je da kuglicu otklonimo za \(x_0=2\) i samo pustimo (\(v_0=0\)), a onda štopericom odredimo vremenske trenutke u kojima kuglica prolazi kroz ravnotežni položaj i da te trenutke zabeležimo:

t	1,18	3,27	5,37	7,46	9,55
x	0	0	0	0	0

Dakle, pored graničnih i početnih uslova \((x_0=2,v_0=0)\) uvešćemo i granične uslove tipa PointSet koji definišu vrednost modelovane funkcije u pojedinim tačkama. Programski kod koji implementira ovaj inverzni problem dat je na sledećem listingu.

Listing 5 Inverzni problem prigušenih oscilacija u 1D. Nepoznati parametri su m i k.

import deepxde as dde
from deepxde.backend import tf
import numpy as np
import matplotlib.pyplot as plt

m = dde.Variable(0.5)
mu = 0.6
k = dde.Variable(2.)

x0, v0 = 2, 0

delta = mu / (2*m)
w0 = tf.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)

bc_x = np.array([1.18, 3.27, 5.37, 7.46, 9.55]).reshape(6,1)
bc_y = np.array([0, 0, 0, 0, 0 ]).reshape(6,1)
ic3 = dde.icbc.PointSetBC(bc_x, bc_y, component=0)

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

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

model = dde.Model(data, net)

# Callback funkcija koja stampa varijablu na svakih 1000 epoha
variable1 = dde.callbacks.VariableValue(k, period=1000)
variable2 = dde.callbacks.VariableValue(m, period=1000)

model.compile("adam", lr=.001, loss_weights=[0.01, 1, 1, 1], metrics=["l2 relative error"], external_trainable_variables=[k,m])
losshistory, train_state = model.train(epochs=50000, callbacks=[variable1, variable2])

Objasnićemo samo delove koji se razlikuju u odnosu na direktni problem rešen u sekciji Implementacija. Za početak, tu su nepoznati parametri koje deklarišemo na sledeći način:

m = dde.Variable(0.5)
k = dde.Variable(2.)

U zagradi se daju početne vrednosti parametra. Sledeće linije definišu pomenuti dodatni PointSet granični uslov (uslove) koji važe u pojedinim tačkama unutar domena:

bc_x = np.array([1.18, 3.27, 5.37, 7.46, 9.55]).reshape(6,1)
bc_y = np.array([0, 0, 0, 0, 0 ]).reshape(6,1)
ic3 = dde.icbc.PointSetBC(bc_x, bc_y, component=0)

Kako bismo obezbedili praćenje vrednosti nepoznatih parametara tokom obuke, potrebno je da postavimo tzv. callback funkcije, koje će se pozivati na svakih 1000 epoha:

variable1 = dde.callbacks.VariableValue(k, period=1000)
variable2 = dde.callbacks.VariableValue(m, period=1000)

Prilikom postavljanja modela postavljamo eksterne varijable za treniranje k i m, dok se pri pozivu treninga navode callback funkcije:

model.compile("adam", lr=.001, loss_weights=[0.01, 1, 1, 1], metrics=["l2 relative error"], external_trainable_variables=[k,m])
losshistory, train_state = model.train(epochs=50000, callbacks=[variable1, variable2])

Nakon završenog obučavanja, dobija se, očigledno pogrešno, rešenje koje je prikazano na Sl. 23. U odnosu na analitičko rešenje koje je postavljeno koristeći vrednosti parametara m=1 i k=2,25, dobijene vrednosti m=0,287 i k=1,22 se puno razlikuju.

Sl. 23 NMPFZ rešenje inverznog problema sa nepoznatim parametrima. Pronađene vrednosti su m=0,287 i k=1,22.

Zašto smo dobili ovako loše rešenje? Odgovor se krije u loše postavljenim graničnim uslovima koji dovode do nejednoznačnosti inverznog problema. Naime, postoji više scenarija, tj. parova (m,k) koji zadovoljavaju granične uslove postavljene samo u tačkama prolaska tega kroz ravnotežni položaj. Očigledno je da moramo da dodamo još neku tačku van ravnotežnog položaja, kako bismo obezbedili jednoznačno rešenje. Zamislimo da smo izmerili i trenutak kada je teg bio na najvećoj negativnoj udaljenosti u odnosu na ravnotežni položaj i kolika je ta udaljenost bila. Dodajmo sada i tu tačku u postavljeni PointSet, koji sada izgleda ovako:

t	1,18	3,27	5,37	7,46	9,55	2,12
x	0	0	0	0	0	-1,67

Pogledom na Sl. 24 odmah se vidi da je poklapanje sa analitičkim rešenjem u ovako postavljenom problemu skoro pa idealno.

Sl. 24 NMPFZ rešenje inverznog problema sa nepoznatim parametrima. Pronađene vrednosti su m=0,98 i k=2,26.

Ispravnost rešenja dodatno potvrđuju parametri m=0,98 i k=2,26, čije su vrednosti veoma bliske onima koje su date u analitičkoj postavci m=1 i k=2,25. Na ovaj način smo pokazali da pristup rešavanju inverznog problema, iako metodološki sličan, ima specifičnosti o kojima treba voditi računa. Kod direktnih problema rešenje je uvek jednoznačno, dok kod inverznih moraju da se obezbede odgovarajući uslovi koji u dovoljnoj meri determinišu rešenje.