Hello all,
I’m trying to estimate 10 parameters (c1 to c10) of an ODE representing the aerodynamic model of a wind turbine, and I’ve decided to do it via Bayesian inference using physics-informed neural networks (PINN) with the package NeuralPDE.jl. The code attached below is a modification of the one in the tutorials, Parameter Estimation with PINNs for ODEs · NeuralPDE.jl, replacing the Lotka-Volterra system by the wind turbine aerodynamic model (and considering 3 exogenous inputs, Tg, θ and vr). For now, I decided to use NNODE() instead of BNNODE() to avoid having to deal with priors.
My issue is the following: Why does the PINN estimation basically kill my dynamics (see attached pictures)? It seems that, rather than adjusting the parameters to better match the data, the PINN makes minimal modifications on them and averages out the dynamics instead. Could it be that the neural network is currently defined as a map from time to my state variable (wr), and it should be a map from time AND inputs to my state variable, or am I mistaken? How could this be done?
I was also going to try another strategy based on the Universal Differential Equations example Automatically Discover Missing Physics by Embedding Machine Learning into Differential Equations · Overview of Julia's SciML (so instead of assuming a parametric model for Cp(λ,θ) and estimate parameters, I would define it as a neural network) to see if that works better, but I’m quite sure it should be possible to obtain much better results for such a simple ODE with PINNs.
Thank you in advance!
using NeuralPDE, Lux, Plots, OrdinaryDiffEq, Distributions, Random
using XLSX, Plots
using OptimizationOptimJL, LineSearches
# Define your wind turbine model with exogenous inputs
function wind_turbine(u, p, t)
# Extract model parameters
c1, c2, c3, c4, c5, c6, c7, c8, c9, c10 = p
# Extract current state
wr = u
# Define constants
μd = 0.05
Rr = 120.998
ρ = 1.225
Ar = π * Rr^2
Jr = 321699000
Jg = 3.777e6
x = 0
# Calculate tip-speed ratio λ
λ = wr * Rr / vr_interp(t)
# Calculate 1/λi (numerical approx. for Cp(θ, λ))
λi_inv = 1 / (λ + c9 * θ_interp(t)) - c10 / (θ_interp(t)^3 + 1)
# Calculate Cp (numerical approx. for Cp(θ, λ))
# Cp = max(c1 * (c2 * λi_inv - c3 * θ_interp(t) - c4) * exp(-c5 * λi_inv) + c6 * λ)
Cp = max(0, c1 * (c2 * λi_inv - c3 * θ_interp(t) - c4 * 1/λi_inv * θ_interp(t) - c5 * θ_interp(t)^x - c6) * exp(-c7 * λi_inv) + c8 * λ)
Tr = 1 / (2 * wr) * ρ * Ar * vr_interp(t)^3 * Cp
# Evaluate differential equation for rotor speed
dwr = ((1 - μd) * Tr - Tg_interp(t)) / (Jr + Jg)
return dwr
end
# Define initial-value problem.
wr0 = 0.8246
u0 = wr0 # Initial state with exogenous inputs
p = [0.5, 116, 0.4, 0, 0, 5, 21, 0.005, 0.089, 0.035] # Model parameters
tspan = (30.0, 229.95)
# Load data from Excel sheet
function load_data(filename)
data = XLSX.readxlsx(filename)["Ark1"] # Assuming your data is in Sheet1
wr = data["A"][1:4000] # Assuming column A contains angular speed data (rad/s)
Tg = data["B"][1:4000] # Assuming column B contains generator torque data (Nm)
θ = data["F"][1:4000] # Assuming column F contains pitch angle data (degrees)
vr = data["D"][1:4000] # Assuming column D contains wind speed data (m/s)
return wr, Tg, θ, vr
end
# Load data from Excel sheet
filename = "…\\data_14ms.xlsx"
wr_data, Tg_data, θ_data, vr_data = load_data(filename)
# ‘Any’ to ‘Float64’ (needed for interpolation)
wr_data = map(Float64, wr_data)
Tg_data = map(Float64, Tg_data)
θ_data = map(Float64, θ_data)
vr_data = map(Float64, vr_data)
t_data = range(tspan[1], tspan[2], length=4000) # Time points
using DataInterpolations
# Create function by interpolating inputs
Tg_interp = BSplineApprox(Tg_data, t_data, 49, 50, :ArcLen, :Average)
θ_interp = BSplineApprox(θ_data, t_data, 49, 50, :ArcLen, :Average)
vr_interp = BSplineApprox(vr_data, t_data, 49, 50, :ArcLen, :Average)
prob_model = ODEProblem(wind_turbine, u0, tspan, p)
# Solve the ODE using stiffness detection and auto-switching algorithm with the standard Tsit5 and Rosenbrock23 for non-stiff and stiff ODE
# and specify saveat to save the solution at specific time points
dt = 0.05 # Sampling time
sol_model = solve(prob_model, AutoTsit5(Rosenbrock23()), saveat=dt)
t_ = sol_model.t
u_ = sol_model.u
#plot(sol_model, xlabel="Time", ylabel="Rotor Speed", label="Rotor speed (model)")
#plot!(t_data, wr_data, label="Rotor speed (data)")
#plot!(legend=:topright)
# Define NN architecture
rng = Random.default_rng()
Random.seed!(rng, 0)
n = 20
chain = Lux.Chain(
Lux.Dense(1, n, Lux.σ),
Lux.Dense(n, n, Lux.σ),
Lux.Dense(n, n, Lux.σ),
Lux.Dense(n, 1)
)
ps, st = Lux.setup(rng, chain) |> Lux.f64
# A function additional_loss(phi, θ) where phi are the neural network trial solutions, θ are the weights of the neural network(s)
function additional_loss(phi, θ) # should it include p?
return sum(abs2, phi(t_data, θ) .- wr_data)/size(wr_data, 1)
end
opt = LBFGS(linesearch = BackTracking())
alg_NN = NNODE(chain, opt, ps; strategy = WeightedIntervalTraining([0.7, 0.2, 0.1], 500), param_estim = true, additional_loss = additional_loss)
sol_NN = solve(prob_model, alg_NN, verbose = true, abstol = 1e-6, maxiters = 500, saveat = t_data)
# plotting solution for x,y for chain
plot(sol_NN, label = "estimated wr")
# comparing it with the original solution
plot!(sol_model, labels = "model wr")
plot!(t_data, wr_data, label = "data wr")
# Check estimated parameters
sol_NN.k.u.p
Estimated parameters (last line):
10-element view(::Vector{Float64}, 902:911) with eltype Float64:
0.500109466140297
115.99999759675735
0.40097664757801016
0.005697616093540125
-1.2050434295049353e-5
4.999987949565712
21.000018552553325
0.00011821424564469212
0.08842105778473758
0.03504212549242426
Last plot (comparison between model simulation, data and PINN estimation):
I’ll attach the data below for reproducibility.