Difficulty in using JuMP with Complex Functions

Thank you for your reply. So this is a reduced version of the code that is reproducible. I will explain why I am writing the code in this way.

I had a very complicated symbolic expression that needed to be integrated. I couldn’t get good accuracy with Julia packages, so I used Sympy to calculate the integration. Then, I saved the Julia code generated by sympy as M_julia2.jl M_julia2.jl (954.1 KB), which is also attached to this reply. Julia parser had a lot of problem parsing its own generated code Parsing Problem In Julia. That is why I am loading the file in this way.

I’m working on an optimization problem, where I need to enforce the dynamic constraint of my system as equality constraints + some other constraint. In order to do that I need to calculate the inertia matrix + some other terms. The M function here calculates only an element of inertia matrix [the first 3-3 elements of a 9-9 matrix], which is a function of state of the system or decision variables which I defined here as qq.

I started learning Julia and Jump with the hope that I can use it to speed up my optimization and leverage some of the solvers that developed in Julia such as complementarity here. But unfortunately, I am stuck generating my terms to use them in the optimization. I am continuing my work in MATLAB again, but I am hoping that by solving this problem I can use Julia and JuMP.

Regarding eval, I agree that it might not be the best option, I tried to use build_function too without evaluating it, and then use jump but still cannot make the code output anything when a symbolic vector is passed.

using Symbolics
Symbolics.@variables xi r L0 q[1:3]

M_julia = read("M_julia2.jl", String)
function parsed_data(expr)
    expr_fixed = replace(expr, r"\.([/*+^])" => s"\1")
    expr_parsed = Meta.parse(expr_fixed)
    return expr_parsed
end
M = parsed_data(M_julia)
@eval @fastmath const M_func(xi,r,L0,q) = $M

function TermsReduced(q)
    L0 = 0.15
    xi = 1
    r = 24e-3
    M = M_func(xi, r, L0, q)
    return M
end

q1 = reshape([0.05, 0.03, 0.01], (3, 1))
M = TermsReduced(q1)
@time TermsReduced(q1)
println("numerical calculation is finished")

using Complementarity
using JuMP
m = MCPModel()
@variable(m, qq[1:3,1:10], start = 0.0)
println("Using JuMP variables in the function")
@time TermsReduced(qq[:,1])