Hello all,
I’m quite new to Turing and Bayesian inference in general, and I’m trying to estimate parameters for the power coefficient Cp in a wind turbine aerodynamic model. The parametric expression for Cp has 6 parameters, c1, c2, c3, c4, c5, c6, where the 1st one is not estimated as it is just scaling c2, c3, c4 and thus it is redundant. The model has 4 states wr, Tg, β, wr_err and is excited by an exogenous input, vr, introduced in the model as a spline interpolation vr_interp(t). As the states have very different magnitudes and I wasn’t sure how that would affect the algorithm, I first compute the operating points of the system and scale each variable in the data by the respective operating point. Then, I also scale the states in the model, and I plot the data and the model simulation states for comparison.
Then I define the Turing @model, I set the prior for σ as a truncated Normal with variance 0.03, as acceptable (scaled) trajectories should not diverge much more than that with respect to the (scaled) data, and I also set priors for my parameters c2, c3, c4, c5, c6 for now. To test the algorithm, I cheat a bit and set the mean values of the Normal priors to the values that I know yield a quite accurate trajectory, which are also the values that I’ve used in the model simulation shown in the previous graphs.
Finally I simulate 2 chains, with NUTS(0.3), 1000 samples per chain. Issues:
- It is taking a really long time to simulate for a model that has just 4 states and 5 parameters, even with relatively narrow priors and NUTS(0.3), it takes around 8 hours. If I increase the NUTS acceptance rate it takes days.
- The chains do not converge, as can be seen in the attached pictures.
Any help on how to solve these issues would be greatly appreciated. My code can be seen below:
using Turing
using DifferentialEquations
using StatsPlots
using LinearAlgebra
using XLSX
using Random
Random.seed!(14);
# 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["C"] # Assuming column C contains pitch angle data (rad)
β = 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_4000.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)
all_data = hcat(wr_data, Tg_data, β_data)'
tspan = (30.0, 229.95) # 2000->129.95, 4000->229.95, 12000 -> 629.95
t_data = range(tspan[1], tspan[2], length=4000) # Time points
using DataInterpolations
vr_interp = BSplineInterpolation(vr_data[1:end], t_data[1:end], 2, :ArcLen, :Average)
# **********************************
# COMPUTE EQUILIBRIUM POINT TO SCALE VARIABLES
# **********************************
using NLsolve
# Define constants
μd = 0.05
Rr = 120.998
ρ = 1.225
Ar = π * Rr^2
Jr = 321699000
Jg = 3.777e6
c1 = 0.1661
c2 = 115.98
c3 = 1.4138
c4 = 5.6336
c5 = 16.2881
c6 = 0.0410
tau_Tg = 0.05
tau_β = 0.01
wr_rated = 0.7917
Kp_Tg = 2.2476e7
Ki_Tg = 0
Tg_off = 2.05e7
Kp_β = 30.2670
Ki_β = 7.9968
vr_eq = 14
function f(x)
[((1 - μd) * 1 / (2 * x[1]) * ρ * Ar * vr_eq^3 * (c1 * (c2 * (1 / (x[1] * Rr / vr_eq + 0.08 * x[3]) - 0.035 / (x[3]^3 + 1)) - c3 * x[3] - c4) * exp(-c5 * (1 / (x[1] * Rr / vr_eq + 0.08 * x[3]) - 0.035 / (x[3]^3 + 1))) + c6 * x[1] * Rr / vr_eq) - x[2]) / (Jr + Jg),
(-Kp_Tg * (x[1] - wr_rated) + Tg_off - x[2]) / tau_Tg,
(Kp_β * (x[1] - wr_rated) + Ki_β * x[4] - x[3]) / tau_β,
x[1] - wr_rated]
end
op = nlsolve(f, [0.79, 2.05e7, 9, 0])
wr_op = op.zero[1]
Tg_op = op.zero[2]
β_op = op.zero[3]
wr_err_op = op.zero[4]
op_points = [wr_op, Tg_op, β_op, wr_err_op] # Equilibrium points for linear system with constant exogenous input 'vr'
# Scaling parameters
wr_scale = wr_op
Tg_scale = Tg_op
β_scale = β_op
wr_err_scale = wr_err_op
wr_data_scaled = wr_data/wr_scale
Tg_data_scaled = Tg_data/Tg_scale
β_data_scaled = β_data/β_scale
all_data_scaled = hcat(wr_data_scaled, Tg_data_scaled, β_data_scaled)'
# **********************************
# WIND TURBINE MODEL
# **********************************
# Define constants
μd = 0.05
Rr = 120.998
ρ = 1.225
Ar = π * Rr^2
Jr = 321699000
Jg = 3.777e6
c1 = 0.1661
c2 = 115.98
c3 = 1.4138
c4 = 5.6336
c5 = 16.2881
c6 = 0.0410
tau_Tg = 0.05
tau_β = 0.01
wr_rated = 0.7917
Kp_Tg = 2.2476e7
Ki_Tg = 0
Tg_off = 2.05e7
Kp_β = 30.2670
Ki_β = 7.9968
# Define wind turbine model with exogenous inputs
function wind_turbine(du, u, p, t)
# Parameters
c2, c3, c4, c5, c6 = p
# States (rel)
wr = u[1]
Tg = u[2]
β = u[3]
wr_err = u[4]
# States (abs)
wr_abs = wr * wr_scale
Tg_abs = Tg * Tg_scale
β_abs = β * β_scale
wr_err_abs = wr_err * wr_err_scale
# Controller
Tg_ref_abs = -Kp_Tg * (wr_abs - wr_rated) + Tg_off
β_ref_abs = Kp_β * (wr_abs - wr_rated) + Ki_β * wr_err_abs
# Calculate tip-speed ratio λ
λ = wr_abs * Rr / vr_interp(t)
# Calculate 1/λi (numerical approx. for Cp(θ, λ))
λi_inv = 1 / (λ + 0.08 * β_abs) - 0.035 / (β_abs^3 + 1)
# Calculate Cp (numerical approx. for Cp(θ, λ))
Cp = c1 * (c2 * λi_inv - c3 * β_abs - c4) * exp(-c5 * λi_inv) + c6 * λ
# Calculate aerodynamic torque
Tr = 1 / (2 * wr_abs) * ρ * Ar * vr_interp(t)^3 * Cp # max(0, Cp)
# State derivatives (rel)
du[1] = ((1 - μd) * Tr - Tg_abs) / (Jr + Jg) / wr_scale
du[2] = (Tg_ref_abs - Tg_abs) / tau_Tg / Tg_scale
du[3] = (β_ref_abs - β_abs) / tau_β / β_scale
du[4] = (wr_abs - wr_rated) / wr_err_scale
end
# Define initial-value problem.
wr0_scaled = wr_data_scaled[1]
Tg0_scaled = Tg_data_scaled[1]
β0_scaled = β_data_scaled[1]
wr_err0_scaled = wr0_scaled
vr0 = vr_data[1]
u0_scaled = [wr0_scaled, Tg0_scaled, β0_scaled, wr_err0_scaled]
p = [c2, c3, c4, c5, c6] # Model parameters
prob_model = ODEProblem(wind_turbine, u0_scaled, tspan, p)
dt = 0.05 # Sampling time
sol_model = solve(prob_model, Tsit5(), saveat=dt)
t_sim = sol_model.t
u_sim = sol_model.u
wr_sim = [inner_vector[1] for inner_vector in u_sim]
Tg_sim = [inner_vector[2] for inner_vector in u_sim]
β_sim = [inner_vector[3] for inner_vector in u_sim]
wr_err_sim = [inner_vector[4] for inner_vector in u_sim]
# Compare data and model simulation
plot(t_sim, wr_sim, xlabel="Time", ylabel="Rotor Speed", label="Rotor speed (model)")
plot!(t_data, wr_data_scaled, label="Rotor speed (data)")
plot!(legend=:topright)
plot(t_sim, Tg_sim, xlabel="Time", ylabel="Generator torque", label="Generator torque (model)")
plot!(t_data, Tg_data_scaled, label="Generator torque (data)")
plot!(legend=:topright)
plot(t_sim, β_sim, xlabel="Time", ylabel="Pitch angle", label="Pitch angle (model)")
plot!(t_data, β_data_scaled, label="Pitch angle (data)")
plot!(legend=:topright)
@model function fitlv(data, prob)
# Priors for variances
σ ~ truncated(Normal(0, 0.03); lower=0, upper=Inf)
# Priors for parameters
c2 ~ truncated(Normal(115.98, 1); lower=113, upper=118)
c3 ~ truncated(Normal(1.4138, 0.5); lower=0, upper=3)
c4 ~ truncated(Normal(5.6336, 0.5); lower=3, upper=7)
c5 ~ truncated(Normal(16.2881, 2); lower=12, upper=20)
c6 ~ truncated(Normal(0.04, 0.02); lower=0, upper=0.08)
# Simulate wind turbine model.
p = [c2,c3,c4,c5,c6]
predicted = solve(prob_model, Tsit5(); p=p, saveat=dt, save_idxs=[1, 2, 3])
# Observations.
for i in 1:length(predicted)
data[:, i] ~ MvNormal(predicted[i], σ^2 * I)
end
return nothing
end
model = fitlv(all_data_scaled, prob_model)
chain = sample(model, NUTS(0.3), MCMCSerial(), 1000, 2; progress=true)
