Error from gradient of a constant

[wupeifan]

Similar to what I posted here #209. #209 was closed though.
It seems that there are still some issues when taking the gradient of an MTK constant.

using ModelingToolkit
using ForwardDiff2: DI
using ChainRules

@scalar_rule(^(x, y), (y * x^(y - 1), Ω * log(x)))

function equations(vars, theta, n_exo, n_endo)
    # X: (k, z)
    # Y: (c)
    XP = vars[1:n_exo]
    X = vars[(n_exo + 1):(2 * n_exo)]
    YP = vars[(2 * n_exo + 1):(2 * n_exo + n_endo)]
    Y = vars[(2 * n_exo + n_endo + 1):(2 * (n_exo + n_endo))]
    beta, alpha, delta, rho, sigma = theta
    equ = [ 1 / Y[1] - beta / YP[1] * (alpha * exp(XP[2]) * XP[1] ^ (alpha - 1) + 1 - delta),
            Y[1] + XP[1] - (1 - delta) * X[1] - exp(X[2]) * X[1] ^ alpha,
            XP[2] - rho * X[2]]
    return equ
end

@variables beta, alpha, delta, rho, sigma
@variables k, z, c, k_p, z_p, c_p

vars_sym = [k_p, z_p, k, z, c_p, c]
theta_sym = [beta, alpha, delta, rho, sigma]

# this one fails
julia> d2HdθdX = DI(vars_arg -> DI(theta_arg -> equations(vars_arg, theta_arg, 2, 1))(theta_sym))(vars_sym)
ERROR: MethodError: no method matching AbstractFloat(::ModelingToolkit.Constant)
Closest candidates are:
  AbstractFloat(::T) where T<:Number at boot.jl:718
  AbstractFloat(::Bool) at float.jl:258
  AbstractFloat(::Int8) at float.jl:259

# but the following works
julia> d2HdθdX = DI(theta_arg -> DI(vars_arg -> equations(vars_arg, theta_arg, 2, 1))(vars_sym))(theta_sym)
18×5 Array{Operation,2}: (omitted)

[wupeifan]

I think it is because MTK.constant is a struct <:Expression. The value of it is a number...