Posing constraints on component array optimizations

Specific Domains Optimization (Mathematical)
1

I have a component array p I want to optimize, p.ps_τ and p.ps_J are parameters of neural nets, and p.ps_F1 and p.ps_F2 are two floats (I followed this example: Simultaneous Fitting of Multiple Neural Networks )

The constraints are:
0<= p.ps_F1+ p.ps_F2<=1, while keeping them both positive.

Here’s my code:

using DifferentialEquations, Lux, Optimization, OptimizationOptimJL, Random, NNlib, SciMLSensitivity

function func!(du, u, p, t)
    τ, J, F1, F2 = p
    du[1] = (0.1*F1 + 0.5*F2 - τ*J - u[1]) / τ
end

u0 = [200.0]
tspan = (1., 100.)

τ(t) = 1000/t
J(t) = 4-t^0.1

function parameterized_func!(du, u, p, t)
    p_ = [τ(t), J(t), 0.1, 0.9]
    func!(du, u, p_, t)
end

prob1 = ODEProblem(parameterized_func!, u0, tspan, saveat=1)
sol1 = solve(prob1)

rng = Random.default_rng()
Random.seed!(rng,1)

m_τ = Lux.Chain(Lux.Dense(1, 30, tanh), Lux.Dense(30, 1, relu))
ps_τ, st_τ = Lux.setup(rng, m_τ)

m_J = Lux.Chain(Lux.Dense(1, 30, tanh), Lux.Dense(30, 1, relu))
ps_J, st_J = Lux.setup(rng, m_J)

ps_F1 = 0.5
ps_F2 = 0.49

ps_τ = Lux.ComponentArray(ps_τ)
ps_J = Lux.ComponentArray(ps_J)
p = Lux.ComponentArray{eltype(ps_τ)}()
p = Lux.ComponentArray(p;ps_τ)
p = Lux.ComponentArray(p;ps_J)
p = Lux.ComponentArray(p;ps_F1)
p = Lux.ComponentArray(p;ps_F2)

function parameterized_func2!(du, u, p, t)
    p_ = [m_τ([t], p.ps_τ, st_τ)[1][1], m_J([t], p.ps_J, st_J)[1][1], p.ps_F1, p.ps_F2]
    func!(du, u, p_, t)
end

prob2 = ODEProblem(parameterized_func!, u0, tspan, p, saveat=sol1.t)

function predict_p(p)
    solve(prob2, Tsit5(), p=p, sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)))
end

function loss_p(p)
    pred = Array(predict_p(p))
    sum(abs2, Array(sol1) .- pred), pred
end

cons(res, p, x) = (res .= [p.ps_F1+p.ps_F2, p.ps_F1, p.ps_F2])
optfun = OptimizationFunction((x, p) -> loss_p(x), Optimization.AutoForwardDiff(); cons=cons)
optprob = OptimizationProblem(optfun, p, lcons=[0.0, 0.0, 0.0], ucons=[1.0, 1.0, 1.0])

res1_uode = Optimization.solve(optprob, IPNewton(), iterations=100)

But I get the following error:

MethodError: no method matching Optim.BarrierLineSearchGrad(::Vector{Float32}, ::Matrix{Float32}, ::Optim.BarrierStateVars{Float64}, ::Optim.BarrierStateVars{Float64})
Closest candidates are:
  Optim.BarrierLineSearchGrad(::Vector{T}, ::Matrix{T}, ::Optim.BarrierStateVars{T}, ::Optim.BarrierStateVars{T}) where T at ~/.julia/packages/Optim/Zq1jM/src/multivariate/solvers/constrained/ipnewton/interior.jl:153

My question is, what is the proper way to define the constraints here? Thank you for your help!