Hello,
I have existing code that uses the DifferentialEquations package that was working fine yesterday, but when run today (23/5/2024) it is giving a message from the solver code:
ERROR: FunctionWrappersWrappers.NoFunctionWrapperFoundError() Stacktrace: [1] _call(::Tuple{}, arg::Tuple{…}, fww::FunctionWrappersWrappers.FunctionWrappersWrapper{…}) @ FunctionWrappersWrappers ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:23 etc, etc
I have tried several different unrelated programs which were previously working, and they all give the same error message.
I am using Julia v1.10.3, DifferentialEquations v7.13.0.
Does anyone have any idea what is going on?
One code example appears below. It is a simple DC machine control program that uses the differential equation ODE solver and a callback to calculate the control.
UPDATE May 24: DifferentialEquations v7.11.0 works OK, so something has happened in v7.12.0 and v7.13.0 to cause this error.
#=
This file simulates a DC machine equations in ODE format
and applies a current and speed controller around the equations.
=#
using Plots: plot, plot!, plotlyjs, display
using DifferentialEquations
plotlyjs()
#=
--- PARAMETER SETUP SECTION ---
=#
# Now set up the control equations and times for the
# control to occur
const START_TIME = 2.0
const END_TIME = 10.0
const Δt = 200e-6
const NUM_SETPOINTS = 2 # Number of ωᵣ setpoints
const CONTROL_DELAY = false # true if a control is delayed
ctrl_t = [i for i in 0.0 : Δt : END_TIME]
# Set up parameters for DC machine differential equations
const Rf = 240.0
const ra = 0.6
const Laa = 0.012
const Lff = 120.0
const Laf = 1.8
const J = 1.0
const Bm = 0.0
const vf = 240.0
const va = 0.0
const Tl = 0.0
const p = [Rf, ra, Laa, Lff, Laf, J, Bm, vf, va, Tl]
const Te_max = 29.2 # Nm
const Ia_max = Te_max / (Laf * vf / Rf) # Amps
const Va_max = 240.0 # Volts
#=
The matrix below is used to store the time, reference
pair for the ωᵣ set points.
=#
setP_ωᵣ = zeros(Float64, NUM_SETPOINTS, 2)
setP_ωᵣ = [START_TIME 100.0;
6.0 25.0]
#=
--- GLOBALS USED IN control! ---
=#
# The following array is used to store the set points
# for later plotting
ωᵣ_sp_arr = zeros(length(ctrl_t)) # The ωᵣ setpoint value array
ia_sp_arr = zeros(length(ctrl_t)) # the ia set point array
ctrl_va_arr = zeros(length(ctrl_t)) # The control voltage array
# Index used to store the setpoint values
setP_array_idx = 0
# The current global ωᵣ setpoint value
setpt_ωᵣ = 0.0
# For the integral error of the controller
integErr_ωᵣ = 0.0
#= ωᵣ control array
• ctrl[1] - the desired torque control value
• ctrl[2] - integral error
• ctrl[3] - ωᵣ reference error
=#
ctrl_ωᵣ = zeros(Float64, 3)
# For the integral error of the controller
integErr_ia = 0.0
#= ia control array
• ctrl[1] - desired voltage value
• ctrl[2] - integral error
• ctrl[3] - ia reference error
=#
ctrl_ia = zeros(Float64, 3)
# Global variable used to store the previous control if
# control delay of one time period is being used.
va_km1 = 0.0
#=
--- END OF control! GLOBALS ---
=#
#=
The parameters for the ωᵣ controller
p_ωᵣ - a parameter list containing
• [1] Δt - the control time step
• [2] Kp_ωᵣ - the proportional gain
• [3] Ti_ωᵣ - integrator time constant
• [4] High saturation level for the torque
• [5] Low saturation level for the torque
=#
ϕₘ_ωᵣ = 0.698 # the desired phase margin in rad
ωc_ωᵣ = (π/2.0 - ϕₘ_ωᵣ) / Δt # approximate cross-over freq rad/sec
Kp_ωᵣ = ωc_ωᵣ * J / Te_max
Ti_ωᵣ = 10.0 / ωc_ωᵣ
p_ωᵣ = zeros(Float64, 5)
p_ωᵣ = (Δt, Kp_ωᵣ, Ti_ωᵣ, Te_max, -Te_max)
#=
The parameters of the ia controller
p_ia - a parameter list containing
• [1] Δt - the control time step
• [2] Kp_ia - the proportional gain
• [3] Ti_ia - integrator time constant
• [4] High saturation level for ia
• [5] Low saturation level for ia
=#
ϕₘ_ia = 0.698 # the desired phase margin in rad
ωc_ia = (π/2.0 - ϕₘ_ia) / Δt # approximate cross-over freq rad/sec
Kp_ia = ωc_ia * Laa
Ti_ia = 10.0 / ωc_ia
p_ia = zeros(Float64, 5)
p_ia = (Δt, Kp_ia, Ti_ia, Va_max, -Va_max)
#=
--- FUNCTION DECLARATIONS ---
=#
"""
dc_mach
@brief Standard form of a differential equation representation for
the DifferentialEquations.jl package for Julia.
The u vector has the following values:
• u[1] - i_f the field current
• u[2] - i_a the armature current
• u[3] - ωᵣ the rotor angular velocity
@params • du - the derivative of u
• u - the state vector defined above
• p - the parameter vector
• t - the current value of time
@returns • du array
"""
function dc_mach!(du, u, p, t)
Rf, ra, Laa, Lff, Laf, J, Bm, vf, va, Tl = p
du[1] = -Rf/Lff * u[1] + 1/Lff * vf
du[2] = -ra/Laa * u[2] - Laf/Laa * u[1] * u[3] + 1/Laa * va
du[3] = -Bm/J * u[3] + Laf/J * u[1] * u[2] - Tl
nothing
end
"""
pi!
@brief This function implements an anti-windup PI
control algorithm.
The first parameter is an array, and is therefore
effectively passed by reference. This allows the
function to be called by multiple controllers.
@params • u - array containing:
• [1] - the control value
• [2] - the integral error
• [3] - the present time error
• ymeas - the measurement to be controlled.
• yref - the desired value for the measurement
• feedFwd - the feedforward term
• p - a parameter list containing
• [1] Δt - the control time step
• [2] Kp - the proportional gain
• [3] Ti - integrator time constant
• [4] SatH - high saturation level
• [5] SatL - low saturation level
@returns • u array:
• [1] - the calculated control value
• [2] - the integral error
• [3] - the present time error
"""
function pi!(u::Vector{Float64}, ymeas::Float64, yref::Float64, feedFwd::Float64, p::Tuple{Vararg{Float64,5}})
# function pi!(u, ymeas, yref, feedFwd, p)
Δt, Kp, Ti, SatH, SatL = p
u[3] = yref - ymeas
integErrTmp = u[2] + Δt / Ti * u[3]
u[1] = Kp * (u[3] + integErrTmp) + feedFwd
# Now check on saturation
if u[1] > SatH
u[1] = SatH
# If the error does no increase saturation then
# update the integral error
if u[3] < 0.0
u[2] = integErrTmp
end
elseif u[1] < SatL
u[1] = SatL
# If the error does no increase saturation then
# update the integral error
if u[3] > 0.0
u[2] = integErrTmp
end
else
u[2] = integErrTmp
end
end
"""
control!
@brief This function is called upon a discrete time event
from the integrator. These events occur upon the
control time period, and result in this function
being called to executed the control strategy.
This function accesses global variables for the
set points of the relevant controllers as well as
the control parameters returned from the pi!
controller algorithm.
@params The integrator object
@returns Changes to the integrator object.
"""
function control!(integrator)
global setP_array_idx, setP_ωᵣ, ωᵣ_sp_arr, ia_sp_arr
global ctrl_va_arr, ctrl_ωᵣ, ctrl_ia, integErr_ωᵣ, integErr_ia
global p_ωᵣ, p_ia, setpt_ωᵣ, START_TIME
global va_km1
setP_array_idx += 1
if integrator.t >= START_TIME
integrator.p[10] = 10.0
# The first controller is the angular velocity controller
# Find the reference value
setp_idx = findall(x -> x == integrator.t, setP_ωᵣ[:,1])
if setp_idx != []
setpt_ωᵣ = setP_ωᵣ[setp_idx[1], 2]
end
ωᵣ_sp_arr[setP_array_idx] = setpt_ωᵣ # save the setpoint
# Now call the ωᵣ controller
feedFwd = 0.0
ctrl_ωᵣ[2] = integErr_ωᵣ
pi!(ctrl_ωᵣ, integrator.u[3], setpt_ωᵣ, feedFwd, p_ωᵣ)
# setpt_ia = ctrl_ωᵣ[1] / (integrator.p[5] * integrator.u[1])
setpt_ia = ctrl_ωᵣ[1]
integErr_ωᵣ = ctrl_ωᵣ[2]
# Now limit setpt_ia - could be larger than Ia_max if
# field current is at incorrect level.
if setpt_ia > Ia_max
setpt_ia = Ia_max
elseif setpt_ia < -Ia_max
setpt_ia = -Ia_max
end
ia_sp_arr[setP_array_idx] = setpt_ia # save the set point
# Now call the ia current controller
# Feed forward the back-emf of the machine
feedFwd = integrator.p[5] * integrator.u[1] * integrator.u[3]
ctrl_ia[2] = integErr_ia
pi!(ctrl_ia, integrator.u[2], setpt_ia, feedFwd, p_ia)
# Save the control integral error
integErr_ia = ctrl_ia[2]
# now limit the control voltage
if ctrl_ia[1] > Va_max
va = Va_max
elseif ctrl_ia[1] < -Va_max
va = -Va_max
else
va = ctrl_ia[1]
end
va_km1 = va
# Put the control out with a delay if selected
@static if CONTROL_DELAY == true
# save the control voltage
ctrl_va_arr[setP_array_idx] = va_km1
# apply the control voltage to the machine
integrator.p[9] = va_km1
end
@static if CONTROL_DELAY == false
# There is no control delay so output the current
# value of the control.
# save the control voltage for later plotting
ctrl_va_arr[setP_array_idx] = va
# apply the control voltage to the machine
integrator.p[9] = va
end
else
ωᵣ_sp_arr[setP_array_idx] = 0.0
ia_sp_arr[setP_array_idx] = 0.0
ctrl_va_arr[setP_array_idx] = 0.0
end
nothing
end
#=
--- MAIN PROGRAM ---
=#
# Initial conditions for DE solution
u0 = [0.0, 0.0, 0.0]
tspan = (0.0, 10.0)
prob = ODEProblem(dc_mach!, u0, tspan, p)
# condition(u, t, integrator) = t ∈ ctrl_t
# cb = DiscreteCallback(condition, control!)
cb = PresetTimeCallback(ctrl_t, control!)
# Start the solution and control problem
sol = solve(prob, callback = cb, tstops = ctrl_t)
# Plot out the results
plot1 = plot(sol, title="DC Machine Main Performance Parameters", linewidth = 2, label=["if" "ia" "ωᵣ"])
display(plot1)
plot2 = plot(sol.t, [sol.u[i][1] for i in 1:length(sol.t)], linewidth = 2, label = "i_f")
display(plot2)
plot3 = plot(sol.t, [sol.u[i][2] for i in 1:length(sol.t)], linewidth=2, label = "iₐ", xlim = (0.0, 10.0))
plot!(ctrl_t, ia_sp_arr, linestyle = :dash, label = "iₐ Set Points", linewidth = 2,)
display(plot3)
plot4 = plot(sol.t, [sol.u[i][3] for i in 1:length(sol.t)], linewidth = 2, label = "ωᵣ")
plot!(ctrl_t, ωᵣ_sp_arr, label = "ωᵣ Set Points", linestyle = :dash, linewidth = 2,)
display(plot4)
plot5 = plot(ctrl_t, ctrl_va_arr, label = "Applied vₐ", linewidth = 2,)
display(plot5)