Q: Symbolics.jl | How to convert symbolic multivariate polynomial to DynamicPolynomials type?

[Edit] The following test code solves the above variable substitution issue using a bit of trivial metaprogramming.

# test_symbolics_to_dynamicpolynomials.jl
using FixedPolynomials
using Symbolics

function main()

    @variables q[1:2], p[1:2], t

    # Generate the numbered variables {qi, pi | i = 1, 2,  . . .}
    numbered_cv_str = "@variables"

    for i = 1:2
        numbered_cv_str = string(numbered_cv_str, " ", "q", string(i), " ", "p", string(i))
    end

    @show eval(Meta.parse(numbered_cv_str))

    # Test expression
    expr = 
        p[1] - q[1]*t - (1//2)*p[1]*(t^2) - (2//1)*q[1]*q[2]*t -
        p[1]*q[2]*(t^2) - p[2]*q[1]*(t^2) + (1//6)*q[1]*(t^3) -
        (2//3)*p[1]*p[2]*(t^3)

    # Substitute {q[i] => qi | i = 1, 2, . . .}
    expr = 
        substitute(
            expr, 
            Dict(
                [q[i] => Num(eval(Symbol("q$(i)"))) for i = 1:2]
            )
        )

    # Substitute {p[i] => pi | i = 1, 2, . . .}
    expr = 
        substitute(
            expr, 
            Dict(
                [p[i] => Num(eval(Symbol("p$(i)"))) for i = 1:2]
            )
        )

    @show expr

    pf_expr = PolyForm(Symbolics.unwrap(expr))

    @show pf_expr.p

    @show Polynomial{Float64}(pf_expr.p)
end

begin
    main()
end

yields

pf_expr.p = p1 - q1t - 2//1q1tq2 - 1//2p1t² - q1t²p2 + 1//6q1t³ - p1t²q2 - 2//3p1t³p2

and

Polynomial{Float64}(pf_expr.p) = -0.6666666666666666p₁t³p₂-p₁t²q₂+0.16666666666666666q₁t³-q₁t²p₂-0.5p₁t²-2.0q₁tq₂-q₁t+p₁
-0.6666666666666666p₁t³p₂-p₁t²q₂+0.16666666666666666q₁t³-q₁t²p₂-0.5p₁t²-2.0q₁tq₂-q₁t+p₁

as desired.