# Modelling heating curves with a PINN

**URL:** <https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521>\
**Category:** Machine Learning\
**Tags:** question, flux, ode\
**Created:** [May 28, 2023, 4:54pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521 "2023-05-28T16:54:14Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 28, 2023, 4:54pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/1 "2023-05-28T16:54:14Z")

</div>

I am Trying to model some temperature curves with a physically informed NN (PINN).

I the physical model I am trying to apply is an [FOPDT](https://apmonitor.com/pdc/index.php/Main/FirstOrderOptimization).

```julia
NN = Chain(
                Dense(1 => 32, tanh),
		Dense(32 => 32, tanh),
		Dense(32 => 1) |> gpu)
	
function lossₚ(i, y) 

        run = floor(Int, i / 130)
	i₀ = run * 130
	iᵣ = i - i₀
	U = (iᵣ < 60) ? 100 : 0

	# currently from scipy.optimize fsolve
	Κ = 66.3
	τ = 145.7
	θ = 73.4

	yₑ = grad.(y)
	yₘ = (.-y .+ Κ .* U .*(i₀ .- θ)) ./ τ
		
	return sum((yₑ .- yₘ) .^2)
end

lossₘ(y_hat, y) = sum((y_hat .- y).^2)
	
loss_ratio = 1.0 / 5.0
lossₜ(y_hat, y) = (1 - loss_ratio) * lossₘ(y_hat, y)+ loss_ratio * lossₚ(enumerate(y))

	
train_mask = 1: trunc(Int, length(static_df.T1)//2) 

X_train = static_df[train_mask, :Q1]
y_train = static_df[train_mask, :T1]
train_set = [(x, y) for x in X_train, y in y_train]

η = 0.001
opt = Flux.setup(Adam(η), NN)
	
for epoch in 1:100 
   Flux.train!(NN, train_set, opt) do m, x, y
       lossₜ(m.(x) .+ y_train[1], y)
   end
end

```

beside the normal loss of the MSE I am trying to add a term for the ODE. But I am not so sure how to add the time depened term, that is why I used enumerate. In the training I also add the start condition to the NN.  
Current problem is something with the type of the tanh activation function, but I don’t get it.

> MethodError: no method matching (::Flux.Dense{typeof(tanh), Matrix{Float32}, Vector{Float32}})(::Float64)

So my detailed questions is how to solve this bug, but I addition I want whether I am god on track, or whether you have any other tip.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 28, 2023, 9:22pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/2 "2023-05-28T21:22:38Z")

</div>

> [@Ordens\_Ritter](#):
>
> MethodError: no method matching (::Flux.Dense{typeof(tanh), Matrix{Float32}, Vector{Float32}})(::Float64)

It’s saying that your neural network is taking in a `Float64` which isn’t allowed, since it’s supposed to be a vector or matrix.

> [@Ordens\_Ritter](#):
>
> `lossₜ(m.(x) .+ y_train[1], y)`

you don’t need to broadcast the model evaluation `m`, just put the different batches in different columns, so use `x'`.

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 28, 2023, 10:51pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/3 "2023-05-28T22:51:42Z")

</div>

Sorry. I did not really understand you answer. Which might be the reason of the problem?

What I now did substitute the real data with a muck up and used a simpler loss function without the casting and changed the definition of train\_set.  
I commented out the older code and added the new one the line below.  
I can only say something is computed slowly.

```julia

# X_train = static_df[train_mask, :Q1]
X_train = 1:130
# y_train = static_df[train_mask, :T1]
y_train = [x^2+21.453 for x in 1:130]
# train_set = [(x, y) for x in X_train, y in y_train]
train_set = [([x,], y) for x in X_train, y in y_train]

η = 0.001
opt = Flux.setup(Adam(η), NN)
	
for epoch in 1:100
  Flux.train!(NN, train_set, opt) do m, x, y
    # lossₜ(m(x) .+ y_train[1], y)
    sum((m(x) .- y).^2)
  end
end

```

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [May 30, 2023, 8:37pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/4 "2023-05-30T20:37:27Z")

</div>

I have tried to go from the toy example from the toy example to an the real data.  
My net is. With 21.543 C being the starting condition.

```julia
	NN = Chain(
		Dense(1 => 16, tanh_fast),
		Dense(16 => 16, tanh_fast),
		Dense(16 => 1))
	

	g(t) = NN(t) .+ 21.543

```

Which I have adapted from @ChrisRackauckas [course](https://book.sciml.ai/notes/03-Introduction_to_Scientific_Machine_Learning_through_Physics-Informed_Neural_Networks/).

Well, and it uses only a single core. Hence, tried

- gpu usage
- mini batching

For batching I am struggling with the documentation to go from:

```julia
train_set = [([x,], y) for x in X_train, y in y_train]

```

to use the DataLoader. My best guess was:

```julia
begin
	train_mask = 1: Int.(length(static_df.T1)//2)
		
	x_train = Float32.(static_df[train_mask, :Q1])

	y_train = Float32.(static_df[train_mask, :T1])

	loader = Flux.DataLoader((data=X_train, label=y_train), batchsize=64, shuffle=false)
	
	η = 0.001
	opt = Flux.setup(Adam(η), NN)
	
	η = 0.001
	opt = Flux.setup(Adam(η), NN)
	
	for epoch ∈ 1:5
		for mini_batch ∈ loader
		  	train!(g, mini_batch, opt) do m, x, y
			    # lossₜ(m(x) .+ y_train[1], y)
				ŷ = m(x)
				sum((ŷ .- y).^2)
			 end
	  end
	end
end

```

When I execute; I get the error:

> MethodError: no method matching (::var"#57#58"{typeof(-), typeof(^), typeof(sum)})(::typeof(Main.var"workspace#30".g), ::Vector{Float32})  
> I understand it, as I have to convert my vector X into a Matrix? But I thought I already had.

About the GPU: I understood, that when I want to train on the GPU I only have to transfer the net onto the GPU. But

```julia

NN = Chain(
		Dense(1 => 16, tanh_fast),
		Dense(16 => 16, tanh_fast),
		Dense(16 => 1)) |> gpu

```

But this results in the error:

> no method matching unsafe\_convert(::Type{Ptr{Float32}}, ::CUDA.CuPtr{Float32})

I wanted to include the model: as second term for the loss function. I hope, I have correctly written the formula.

function lossₚ((i, x)::Tuple{Int, Float32})

```
	run = trunc(Int, i / 130)
	i₀ = run * 130
	iᵣ = i - i₀
	U = (iᵣ < 60) ? 100 : 0 # Input heat is either fully on, or off

	Κ = 0.997887
	τ = 1439.8
	θ = -661.1

	ŷ = g(x)
	Δyₗ = ŷ - g(x - 1)
	
	yᵣ = (.-ŷ .+ Κ .* U .*(i₀ .- θ)) ./ τ
	
	return sum((Δyₗ .- yᵣ) .^2)
end

lossₘ(ŷ, y) = sum((ŷ .- y).^2)

loss_ratio = 1.0 / 5.0
lossₜ(ŷ, y, x) = (1 - loss_ratio) * lossₘ(ŷ, y)+ loss_ratio * lossₚ.(enumerate(y))

```

But the code breaks. Probably at `ŷ = g.(x)` with the error: MethodError: no method matching (::Flux.Dense{typeof(NNlib.tanh\_fast), Matrix{Float32}, Vector{Float32}})(::Float32)

1. I do not understand, why m(x) works in the training loop than.
2. I am also not sure, whether en enumerate is the correct thing, because I need the index of the train data set, as it marks the point in time. Or should I include it in the training dataset?

To be honest, I find the documentation of flux confusion, because I it does mix several styles. For instance the training is either done with train in a loop or functional? and with a withgradient function. I would be nice to have some paragraphs about theses styles. And how I can the several same things in these different ways.  
For instance saving and print the loss.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 9, 2023, 10:10pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/5 "2023-06-09T22:10:57Z")

</div>

@xtalax can you take a look at this one? It’s been sitting in my list too long.

Big thing though I can see is that it’s mixing Float32 and Float64.

> [@Ordens\_Ritter](#):
>
> `g(t) = NN(t) .+ 21.543`

Should be `g(t) = NN(t) .+ Float32(21.543)` and there’s a few other places with that. @Ordens_Ritter did you already correct that? I know someone asked a similar question on Slack with a bunch of Float64 values mixed into a Float32 computation.

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [June 11, 2023, 10:37am UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/6 "2023-06-11T10:37:41Z")

</div>

Hi @ChrisRackauckas ,

I did a total overhaul.  
I changed the input and am able to get the PINN running without the physical loss term, with the caveat that I not able to use mini batches,  
or not able to use the gpu.

The network is:

```julia
NN = Chain(
	Dense(2 => 16, relu),
	Dense(16 => 1)
) |> f64
	

```

Streaming the network to the gpu… would be nice, but without the physical loss everything more than just is fast enough.

The X is the heating either 0 or 100 and the previous temperature and y the resulting updated temperature (in Kelvin). I get the data from a DataFrame.

```julia
train_mask = 1: Int.(length(static_df.T1)//2)
	
train_df = static_df[train_mask, :]

x_train_df = train_df[:,[:Q1, :T1prev]]
x_train_m = x_train_df|> Matrix
x_train_m = [x for x in eachrow(x_train_m)]

y_train = train_df[:, :T1]

train_set = [(x, y) for (x, y) in zip(x_train_m, y_train)]
	
η = 1e-5
opt = Flux.setup(Adam(η), NN)
epochs = 20

losses = []
	
for i in 1:epochs
	local_loss = []
	for (x, y) in train_set
		loss, grads = Flux.withgradient(NN) do m
			ŷ = m(x)
			lossₜ(ŷ, y)
		end
		Flux.update!(opt, NN, grads[1])
		push!(local_loss, loss)
	end
	push!(losses, mean(local_loss))
end

```

I try to take the physical loss not from the training data set but from a set of interesting points. Therefore I selected data points, where the first and second derivative are zero.  
For each Data point I compute:

```julia
function loss_physical(x)

	Κ = 0.997887
	τ = 1439.8
	θ = -661.1

	# Q1 in df is index 4 in vector, T1prev is index 10
	ŷ = [x[4], x[10]] |> NN |> first
	prev_index = convert(Int, x[1]) - 1
	xₚᵣₑᵥ = static_m[prev_index, :]

	ŷₚᵣₑᵥ = [xₚᵣₑᵥ[4], xₚᵣₑᵥ[10]] |> NN |> first
	Δyₗ = ŷ - ŷₚᵣₑᵥ
	U = x[4]
	i₀ = x[8] # duration
	
	yᵣ = (.-ŷ .+ Κ .* U .*(i₀ .- θ)) ./ τ
	
	return sum((Δyₗ .- yᵣ) .^2)
end

```

Perhaps there is a better version of getting:  
τ\_p \frac{dy(t)}{(dt)}=−y(t)+K\_pu(t−θ\_p)

The overall loss is:

```julia
loss_physical() = eachrow(sample_time_m) .|> loss_physical |> sum

loss_ratio = 1 / size(sample_time_df)[1]
	
loss_total(ŷ, y) = (1 - loss_ratio) * mae(ŷ, y)+ loss_ratio * loss_physical()

```

The current problem from my understanding is, that I compute the complete physical loss for every training point. Therefore I would need to get batching working.  
In fact it is so slow, that I stopped after halve an hour.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 11, 2023, 1:32pm UTC](https://discourse.julialang.org/t/modelling-heating-curves-with-a-pinn/99521/7 "2023-06-11T13:32:25Z")

</div>

Instead of doing it point-by-point via broadcast, make the input to the NN a matrix where each column is a sample and GPU it. It should be a rather quick computation like that. You’ll need to adjust `[x[4], x[10]] |> NN |> first` to accommodate that form.