[slack] <Adedotun Agbemuko> #diffeq-bridged #sciml
This issue is based on the DifferentialEquations package:
To give a proper context of my problem, I have a system of equations i want to solve in time in the form
dx = Ax + Bu
where:
x = [x(1), x(2), ...x(24)] # state vector
u = [u(1), u(2), ...u(24)] # control reference input
I am attempting to implement a couple of hybrid callbacks. By this i mean one a discrete call back (e.g. to make a step change in one or more
of inputs u at arbitrary, but known times) and another a continuous call back for a continuous check of certain states x(13) to x(24).
In principle, the values of the solutions of x(13) to x(24) must be saturated if they exceed a value on the low or high side.
For example, if any of x(13):x(24) at any time is > 1.2 then it is saturated to 1.2. If it is less than -1.2 then the value is saturated at
-1.2. This saturation is important to the credibility of the solutions. It has to do with equipment rating/tolerance and actuator capability
The discrete call back alone was implemented and this worked perfectly as expected. However, when i try to combine with the continuous callback,
i get too many errors to understand or make sense. What am i doing wrong? below is a snippet of what am trying to do:
function time2step_disco(x,t,integrator)
t==10.0
end
function step_change_disco!(integrator)
matrix_size = size(integrator.p)
integrator.p[6, matrix_size(2)]-= integrator.p[6, matrix_size(2)]
end
basically making driving input 6 to zero value from the original value by subtracting the original value from itself.
function limit_check(I_limits,x,t,integrator)
I_limits[1] = x[13:24] > 1.2
I_limits[2] = x[13:24] < -1.2
end
function current_limiting!(integrator,idx)
if idx==1
integrator.u[13:24] = 1.2
elseif idx==2
integrator.u[13:24] = -1.2
end
end
Then in the call to the solver
nodal_disconxn = DiscreteCallback(time2step_disco,step_change_disco!)
current_saturation = VectorContinuousCallback(limit_check,current_limiting!,2) # 2 implies there are only 2 checks made e.g. the high and low sides of the saturation curve
event_set = CallbackSet(nodal_disconxn,current_saturation)
so event_set is passed on to the ODEproblem.
[slack] <Dan Burns> I’d like to solve a simple ODE with units from Unitful
, but I get an error related to initialization. MWE:
```using DifferentialEquations, Unitful
tspan = (0.0u"s", 10.0u"s")
x0 = [0.0u"m", 0.0u"m/s"]
F(t) = 1
A = [0u"s^-1" 1; 0u"s^-2" 0u"s^-1"]
B = [0u"m/s"; 1u"m/s^2"]
di(x,u,t) = Ax + Bu(t)
prob = ODEProblem(di, x0, tspan, F)
sol = solve(prob, Tsit5())Which yields:
ERROR: ArgumentError: zero(Quantity{Float64,D,U} where U where D) not defined.
Stacktrace:
[1] zero(::Type{Quantity{Float64,D,U} where U where D}) at /Users/dan/.julia/packages/Unitful/m6pR3/src/quantities.jl:374``
This is the same error one gets when entering
zero(x0)` , which makes me suspect initialization. Any suggestions? Thanks in advance!
x0 = ArrayPartition([0.0u"m", 0.0u"m/s"])
ERROR: LoadError: MethodError: Cannot `convert` an object of type Array{Quantity{Float64,D,U} where U where D,1} to an object of type ArrayPartition{Quantity{Float64,D,U} where U where D,Tuple{Array{Quantity{Float64,D,U} where U where D,1}}}
x0 = ArrayPartition(0.0u"m", 0.0u"m/s")
, but the https://diffeq.sciml.ai/stable/features/diffeq_arrays/#ArrayPartitions indicate that this is the correct use case, so thanks for the pointer!
x0 = ComponentArray(pos=0.0u"m", vel=0.0u"m/s")
. https://jonniedie.github.io/ComponentArrays.jl/stable/examples/adaptive_control/ is an example of using ComponentArrays for more deeply nested state-space problems too. It doesn't use Unitful.jl in this example, but there's no reason you wouldn't be able to do this with units as well.