PINN for temperature prediction with external disturbances ; timedependent problen

[mariaade26]

Hello everyone, I'm really new to python and PINA and, by following the tutorials and examples, I was trying to implement my own PINN ( time dependent problem) to predict the temperature inside a building, according to the physics ODE, which depends on external disturbances, in particular a temperature vector. I wrote a code that works when this external disturbance it's a single value, but it doesn't work when I want to use a time-series for it. The problems started when I introduced an " init() " method. I will attach to this issue both the working and non-working code and I'll be glad if someone could help me because I didn't understardwhat is wrong and why. Thank you in advance.

WORKING CODE

from pina.problem import TimeDependentProblem
from pina import Plotter, Trainer, Condition
from pina.geometry import CartesianDomain
from pina.solvers import PINN
from pina.model import FeedForward
from pina.equation import Equation,FixedValue
from pina.operators import grad
import torch
import numpy as np
from scipy.integrate import odeint

from torch.nn import Softplus

class TempPred(TimeDependentProblem):
    output_variables = ['T']
    temporal_domain = CartesianDomain({'t': [0, 10]})

    def temp_equation(input_,output_):

        dTdt = grad(output_, input_, components = ['T'], d = ['t'])
        R= 0.005
        C = 27800
        Testerna = 10
        force_term = 1/(R*C) * (Testerna - output_.extract(['T']))
        return dTdt - force_term

    def initial_condition(input_, output_):
        T0 = 15
        return output_.extract(['T'])- T0

    conditions = {
        't0': Condition(location=CartesianDomain({'t': 0}), equation=Equation(initial_condition)),
        'D': Condition(location=CartesianDomain({'t': [0, 10]}), equation=Equation(temp_equation))
    }

problem=TempPred()
problem.discretise_domain(n=100, mode = 'random', locations=['D', 't0'])

#creazione rete neurale
model = FeedForward(
    layers = [10,10],
    func = torch.nn.Tanh,
    output_dimensions = len(problem.output_variables),
    input_dimensions = len(problem.input_variables)
)
pinn = PINN(problem=problem, model=model, optimizer_kwargs={'lr': 0.006})
trainer = Trainer(solver =pinn,accelerator='cpu', max_epochs=1000)
trainer.train()

pl = Plotter()
pl.plot(solver=pinn)`

NON WORKING CODE

import numpy as np
from pina.problem import TimeDependentProblem
from pina import Plotter, Trainer, Condition
from pina.geometry import CartesianDomain
from pina.solvers import PINN
from pina.model import FeedForward
from pina.equation.equation import Equation
from pina.operators import grad
import torch
from torch.nn import Softplus
from scipy.integrate import odeint
from pina.callbacks import MetricTracker
import matplotlib.pyplot as plt

########################
#EQUAZIONE ORIGINALE

tspan = np.linspace(0, 500, 500)
T0 = 15.0
Text = []

with open("Testerna.csv", 'r', newline = '') as file:
    for line in file:
        Text.append(float(line.strip()))
       # print(line.strip())
T0 = 15
p= [0.005, 27800, 1.5]
#tspan = np.linspace(0,52560,52560)
disturbance = Text[0:len(tspan)]

def model(u,t,tspan,p,disturbance_fun):

    R = p[0]
    C = p[1]
    a = p[2]

    disturbance_interp = np.interp(t,tspan, disturbance_fun)
    dudt = 1/(R*C) * (disturbance_interp - u)

    return dudt

sol= odeint(model, T0, tspan, args = (tspan,p, disturbance))

plt.plot(tspan, sol, 'b', label = "Tinterna", linewidth = 1)
plt.grid()
plt.show()

##############################     PINN

class TempPred(TimeDependentProblem):

    output_variables = ['T']
    temporal_domain = CartesianDomain({'t': [0.0, 10.0]})

    def __init__(self, disturbance):
        super().__init__()
        self.disturbance = disturbance

    def temp_equation(self, input_, output_):

        dTdt = grad(output_, input_, components = ['T'], d = ['t'])

        R= 0.005
        C = 27800.0
        #Testerna = 10.0
        disturbance_interp = np.interp(input_.extract(['t']), tspan, self.disturbance)

        force_term = 1/(R*C) * (disturbance_interp - output_.extract(['T']))
        return dTdt - force_term

    def initial_condition(input_, output_):
        T0 = 15.0
        return output_.extract(['T'])- T0

    conditions = {
        't0': Condition(location=CartesianDomain({'t': 0.0}), equation=Equation(initial_condition)),
        'D': Condition(location=CartesianDomain({'t': [0.0, 500.0]}), equation=Equation(temp_equation))
    }

problem=TempPred(disturbance)
problem.discretise_domain(n=500, mode = 'random', locations=['D', 't0'])

#creazione rete neurale
model = FeedForward(
    layers = [10,10],
    func = torch.nn.Tanh,
    output_dimensions = len(problem.output_variables),
    input_dimensions = len(problem.input_variables)
)
pinn = PINN(problem=problem, model=model, optimizer_kwargs={'lr': 0.006})
trainer = Trainer(solver = pinn,callbacks=[MetricTracker()] , accelerator='cpu',enable_model_summary=False, max_epochs=1000)
trainer.train()
#trainer.logged_metrics

pl = Plotter()
pl.plot(solver=pinn)

#pl.plot_loss(trainer=trainer, label = 'mean_loss', logy=True)

[dario-coscia]

Ciao @mariaade26 👋🏻

It's great to see PINNs used for scientific applications, and I hope the tutorials helped you start with PINA!

Regarding your problem, the issue is not related to the definition of an init but to the fact that you have passed the self argument into a function (temp_equation) where the residual must be computed. Those functions in pina expect only input, output, and (optionally for inverse problem modelling) parameters params.

You can solve your problem very easily by using static variables, which are variables shared among all instances of a class, rather than being unique to each class instance. You can define a new static variable by simply calling TempPred.disturbance=... before or after initialization. To call a static variable you simply type TempPred.disturbance. Here is an example following a snippet of your code.

class TempPred(TimeDependentProblem):

    output_variables = ['T']
    temporal_domain = CartesianDomain({'t': [0.0, 10.0]})
    # # ADDING disturbance (INSIDE the class)
    # TempPred.disturbance = disturbance

    def temp_equation(input_, output_):

        dTdt = grad(output_, input_, components = ['T'], d = ['t'])

        R= 0.005
        C = 27800.0
        #Testerna = 10.0
        print(TempPred.disturbance)
        disturbance_interp = np.interp(input_.extract(['t']), tspan, TempPred.disturbance)

        force_term = 1/(R*C) * (disturbance_interp - output_.extract(['T']))
        return dTdt - force_term

    def initial_condition(input_, output_):
        T0 = 15.0
        return output_.extract(['T'])- T0

    conditions = {
        't0': Condition(location=CartesianDomain({'t': 0.0}), equation=Equation(initial_condition)),
        'D': Condition(location=CartesianDomain({'t': [0.0, 500.0]}), equation=Equation(temp_equation))
    }

# OR ADDING disturbance (OUTSIDE the class)
TempPred.disturbance = disturbance
problem=TempPred()

# code ...

Notice that static methods can also be updated during training. For example, suppose you want to change the disturbance at the end of every epoch before the optimization step is called. This can be simply achieved by creating a Callback:

# update
from pytorch_lightning.callbacks import Callback
class Update(Callback):
    def on_train_epoch_end(self, trainer, __):
        trainer.solver.problem.__class__.disturbance = ...

For more on callbacks and where to put them have a look at the lightning documentation.

Let me know how it goes and if you like the package leave us a star ⭐️ which helps us grow the community!😄

[mariaade26]

Ciao! I added the disturbance outside the class, as you suggested and ,after fixing some errors related to the Tensors definitions, the code runs.

In addition, I would like to ask you how can I plot, at the end of the training, both the true solution and the predicted one, because right now it only plots the predicted one (which is still really bad😅). Here is the actual code!

import numpy as np
from pina.problem import TimeDependentProblem
from pina import Plotter, Trainer, Condition
from pina.geometry import CartesianDomain
from pina.solvers import PINN
from pina.model import FeedForward
from pina.equation.equation import Equation
from pina.operators import grad
import torch
from torch.nn import Softplus
from scipy.integrate import odeint
from pina.callbacks import MetricTracker
import matplotlib.pyplot as plt

########################
#EQUAZIONE ORIGINALE

tspan = np.linspace(0, 1000, 1000)
T0 = 15.0
Text = []

with open("Testerna.csv", 'r', newline = '') as file:
    for line in file:
        Text.append(float(line.strip()))
       # print(line.strip())
T0 = 15
p= [0.005, 27800, 1.5]
#tspan = np.linspace(0,52560,52560)
disturbance = Text[0:len(tspan)]

def model(u,t,tspan,p,disturbance_fun):

    R = p[0]
    C = p[1]
    a = p[2]

    disturbance_interp = np.interp(t,tspan, disturbance_fun)
    dudt = 1/(R*C) * (disturbance_interp - u)

    return dudt

sol= odeint(model, T0, tspan, args = (tspan,p, disturbance))

plt.plot(tspan, sol, 'b', label = "Tinterna", linewidth = 1)
plt.grid()
plt.show()

##############################     PINN

class TempPred(TimeDependentProblem):

    output_variables = ['T']
    temporal_domain = CartesianDomain({'t': [0.0, 1000.0]})

    #bisogna aggiungere la disturbance nella classe o fuori
    #TempPred.disturbance = disturbance

    def temp_equation(input_, output_):

        dTdt = grad(output_, input_, components = ['T'], d = ['t'])

        R= 0.005
        C = 27800.0
        #print(TempPred.disturbance)
        disturbance_interp = np.interp(input_.extract(['t']).detach().numpy(), tspan, TempPred.disturbance)
        disturbance_interp = torch.Tensor(disturbance_interp)

        force_term = 1/(R*C) * (disturbance_interp - output_.extract(['T']))
        return dTdt - force_term

    def initial_condition(input_, output_):
        T0 = 15.0
        return output_.extract(['T']) - T0

    conditions = {
        't0': Condition(location=CartesianDomain({'t': 0.0}), equation=Equation(initial_condition)),
        'D': Condition(location=CartesianDomain({'t': [0.0, 1000.0]}), equation=Equation(temp_equation))
    }

TempPred.disturbance = disturbance
problem=TempPred()
problem.discretise_domain(n=1000, mode = 'random', locations=['D', 't0'])

#creazione rete neurale
model = FeedForward(
    layers = [10,10],
    func = torch.nn.Sigmoid,
    output_dimensions = len(problem.output_variables),
    input_dimensions = len(problem.input_variables)
)
pinn = PINN(problem=problem, model=model, optimizer_kwargs={'lr': 0.006})
trainer = Trainer(solver = pinn,callbacks=[MetricTracker()] , accelerator='cpu',enable_model_summary=False, max_epochs=1000)
trainer.train()

pl = Plotter()
pl.plot(solver=pinn)

[dario-coscia]

Hi, I am glad you managed to have a working code! In the Parametric Poisson Tutorial you can see how to plot also the original solution. You have to define a static function, something like:

def my_solution(self, pts):
        # pts are your times
        return .....   
 truth_solution = my_solution

We do not have support for non-analytical solution plotting at the moment. But you can always achieve it using matplotlib, a PINA dependency. Something like:

import matplotlib.pyplot as plt
# evaluation points
pts = ...
# true solution
nn_ = solver.neural_net(pts).detach()
true_ = compute_true(pts) # here you compute your true solution
# plot
plt.subplots(1, 3, 1) # true
plt.plot(pts, true_)
plt.subplots(1, 3, 2) # neural net
plt.plot(pts, nn_)
plt.subplots(1, 3, 3) # absolute difference
plt.plot(pts, (true_-nn_).abs())
plt.show()

Hope this is helpful!

[dario-coscia]

Closing this issue, as completed! @mariaade26 If you need further help, feel free to open new issues😄👋🏻