using Turing
using DifferentialEquations
using Plots
# Define the differential equation model for the microgrid system
function microgrid_model(dx, x, p, t)
P_pv, P_wt, P_dg, P_load = x
A_pv, R_pv, A_wt, R_wt, A_dg, R_dg, A_load, R_load = p
dx[1] = A_pv * G_t[t] - P_pv / R_pv
dx[2] = A_wt * wind_speed[t] - P_wt / R_wt
dx[3] = A_dg - P_dg / R_dg
dx[4] = A_load - P_load / R_load
end
# Define the prior distributions for the model parameters
A_pv = Normal(5, 2)
R_pv = Gamma(2, 0.5)
A_wt = Normal(5, 2)
R_wt = Gamma(2, 0.5)
A_dg = Normal(5, 2)
R_dg = Gamma(2, 0.5)
A_load = Normal(5, 2)
R_load = Gamma(2, 0.5)
p=[5, 2, 5, 2, 5, 2, 5, 2]
G_t = [0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3] # Solar radiation data in W/m^2
wind_speed = [2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0] # Wind speed data in m/s
# Use the `solve` function from the DifferentialEquations.jl package to solve the differential equations
# and get the solution at each point in time
tspan = (0.0, 10.0)
x0 = [1.0, 2.0, 3.0, 4.0]
prob = ODEProblem(microgrid_model, x0, tspan, p)
sol = solve(prob, Tsit5())
# Generate some synthetic data from the solution of the differential equation model
data = [sol(t) + randn() * 0.1 for t in tspan[1]:0.1:tspan[2]]
# Perform inference
microgrid_trace = sample(microgrid_model, data, HMC(1000, 0.1), p)
You will certainly get help from multiple users shortly But just to help everyone, it would be easier to read your code if you format it as code by using ``` backquotes. Please read: make it easier to help you
2 Likes
In your microgrid_model
function, you have G_t[t]
and windspeed[t]
. Indexing expects t
to be an Int
, but inside solve
, t
is a Float64
.
I don’t know what the intention of the code is, maybe you want something like G_t[Int(ceil(t))]
. More likely, you want something such as
using Interpolations
Ginterp = linear_interpolation(ts, G_t,extrapolation_bc=Line())
# then inside microgrid_model replace G_t[t] with Ginterp(t)
2 Likes
You cannot index arrays with a floating point number; in your code you do it twice: wind_speed[t] and G_t[t]
2 Likes
i will try thanks
yes i see you are right thanks