Difficulty in using JuMP with Complex Functions

Specific Domains Optimization (Mathematical)
, ,
1

I am trying to create an optimization problem with JuMP. In order to create the constraints I have a function, which I generated its elements with Symbolics. When I pass numerical values, I can calculate the terms in tenth of a mili second, but when I pass Symbolics or JuMP variables, even after waiting two hours I couldn’t make it to work. Is there any work around that I can solve this problem? I will put the pseudo code here:

Here how each function is generated with symbolics:

@eval @fastmath const T_func(xi, r, L0, q) = $T
println("T_func_generated")
@eval @fastmath const TE1_diff_1_func(xi, r, L0, q) = $TE1_diff_1
println("TE1_diff_1_func_generated")
@eval @fastmath const TE1_diff_2_func(xi, r, L0, q) = $TE1_diff_2
println("TE1_diff_2_func_generated")
@eval @fastmath const TE1_diff_3_func(xi, r, L0, q) = $TE1_diff_3
println("TE1_diff_3_func_generated")
@eval @fastmath const p1_int_func(xi, r, L0, q) = $p1_int
println("p1_int_func_generated")
@eval @fastmath const p1_p1_transpose_int_func(xi, r, L0, q) = $p1_p1_transpose_int
println("p1_p1_transpose_int_func_generated")
@eval @fastmath const p1_diff_1_int_func(xi, r, L0, q) = $p1_diff_1_int
println("p1_diff_1_int_func_generated")
@eval @fastmath const p1_diff_2_int_func(xi, r, L0, q) = $p1_diff_2_int
println("p1_diff_2_int_func_generated")
@eval @fastmath const p1_diff_3_int_func(xi, r, L0, q) = $p1_diff_3_int
println("p1_diff_3_int_func_generated")
@eval @fastmath const p1_p1_diff_1_t_int_func(xi, r, L0, q) = $p1_p1_diff_1_t_int
println("p1_p1_diff_1_t_int_func_generated")
@eval @fastmath const p1_p1_diff_2_t_int_func(xi, r, L0, q) = $p1_p1_diff_2_t_int
println("p1_p1_diff_2_t_int_func_generated")
@eval @fastmath const p1_p1_diff_3_t_int_func(xi, r, L0, q) = $p1_p1_diff_3_t_int
println("p1_p1_diff_3_t_int_func_generated")
@eval @fastmath const M_func(xi,r,L0,q) = $M

This is the rest of the code that I wrap these elements to get my terms:

function Kin(L0,r,q)
    T = T_func(1,r,L0,q)
    TE1_diff_1 = TE1_diff_1_func(1,r,L0,q)
    TE1_diff_2 = TE1_diff_2_func(1,r,L0,q)
    TE1_diff_3 = TE1_diff_3_func(1,r,L0,q)
    p1_int = p1_int_func(1,r,L0,q)
    p1_p1_transpose_int = p1_p1_transpose_int_func(1,r,L0,q)
    p1_diff_1_int = p1_diff_1_int_func(1,r,L0,q)
    p1_diff_2_int = p1_diff_2_int_func(1,r,L0,q)
    p1_diff_3_int = p1_diff_3_int_func(1,r,L0,q)
    p1_p1_diff_1_t_int = p1_p1_diff_1_t_int_func(1,r,L0,q)
    p1_p1_diff_2_t_int = p1_p1_diff_2_t_int_func(1,r,L0,q)
    p1_p1_diff_3_t_int = p1_p1_diff_3_t_int_func(1,r,L0,q)
    M = M_func(1,r,L0,q)
    return T,TE1_diff_1,TE1_diff_2,TE1_diff_3,p1_int,p1_p1_transpose_int,p1_diff_1_int,p1_diff_2_int,p1_diff_3_int,p1_p1_diff_1_t_int,p1_p1_diff_2_t_int,p1_p1_diff_3_t_int, M
end

include("JC_func.jl")
include("DynamicsTerm.jl")
include("BiasT.jl")
include("RP.jl")

function Terms(q)
    L0 = 0.15
    beta1 = pi/3
    beta2 = -pi/3
    b1 = 40e-3
    b2 = 40e-3
    m1, m2, m3 = 1.17, 0.54, 0.265
    r1, r2, r3 = 46.32e-3, 32e-3, 24e-3
    g = [0,0,-9.81]

    q1 = q[1:3]
    q2 = q[4:6]
    q3 = q[7:9] 


    TE1, TE1_diff_1, TE1_diff_2, TE1_diff_3, P1_int, p1_p1_int, p1_diff_1_int, p1_diff_2_int, p1_diff_3_int, p1_p1_diff_1_int, p1_p1_diff_2_int, p1_p1_diff_3_int, JJ1 = Kin(L0, r1, q1)
    TE2, TE2_diff_1, TE2_diff_2, TE2_diff_3, P2_int, p2_p2_int, p2_diff_1_int, p2_diff_2_int, p2_diff_3_int, p2_p2_diff_1_int, p2_p2_diff_2_int, p2_p2_diff_3_int, JJ2 = Kin(L0, r2, q2)
    TE3, TE3_diff_1, TE3_diff_2, TE3_diff_3, P3_int, p3_p3_int, p3_diff_1_int, p3_diff_2_int, p3_diff_3_int, p3_p3_diff_1_int, p3_p3_diff_2_int, p3_p3_diff_3_int, JJ3 = Kin(L0, r3, q3)

    R1B, P1B = BiasT(beta1,b1)
    R2B, P2B = BiasT(beta2,b2)
    R1, P1 = RP(TE1);
    R2, P2 = RP(TE2);
    R3, P3 = RP(TE3);

    R1_diff_1, p1_diff_1 = RP(TE1_diff_1)
    R1_diff_2, p1_diff_2 = RP(TE1_diff_2)
    R1_diff_3, p1_diff_3 = RP(TE1_diff_3)

    R2_diff_1, p2_diff_1 = RP(TE2_diff_1)
    R2_diff_2, p2_diff_2 = RP(TE2_diff_2)
    R2_diff_3, p2_diff_3 = RP(TE2_diff_3)

    R3_diff_1, p3_diff_1 = RP(TE3_diff_1)
    R3_diff_2, p3_diff_2 = RP(TE3_diff_2)
    R3_diff_3, p3_diff_3 = RP(TE3_diff_3)
    (Jv1, Jw1, Jv2, Jw2, Jv3, Jw3, Jv1b, Jw1b, Jv2b, Jw2b, Jv1_base, Jw1_base, Jv2_base, Jw2_base, Jv3_base, Jw3_base, P_out, R_out) = JC(R1, R1B, R2, R2B, R3, P1, P1B, P2, P2B, P3,
    p1_diff_1, p1_diff_2, p1_diff_3, p2_diff_1, p2_diff_2, p2_diff_3, p3_diff_1, p3_diff_2, p3_diff_3,
    R1_diff_1, R1_diff_2, R1_diff_3, R2_diff_1, R2_diff_2, R2_diff_3, R3_diff_1, R3_diff_2, R3_diff_3)

    (M,G) = Dyn(JJ1, JJ2, JJ3,
    Jv1, Jw1, Jv2, Jw2, Jv3, Jw3,
    P1_int, P2_int, P3_int,
    p1_p1_int, p2_p2_int, p3_p3_int,
    p1_diff_1_int, p1_diff_2_int, p1_diff_3_int,
    p2_diff_1_int, p2_diff_2_int, p2_diff_3_int,
    p3_diff_1_int, p3_diff_2_int, p3_diff_3_int,
    p1_p1_diff_1_int, p1_p1_diff_2_int, p1_p1_diff_3_int,
    p2_p2_diff_1_int, p2_p2_diff_2_int, p2_p2_diff_3_int,
    p3_p3_diff_1_int, p3_p3_diff_2_int, p3_p3_diff_3_int,
    m1, m2, m3,
    R1, R1B, R2, R2B, R3,
    g)

    return P_out, R_out, Jv3_base, Jw3_base, M, G

end

when I do
@time Terms(numerica_variables)
I get sub,mili second time
but when I do

@variables qq[1:9]
@time Terms(qq)

It just get stuck unfortunately.

2

Can you provide a reproducible example? I don’t see anything related to JuMP in your example.

The issue is almost certainly that the expression returned by Terms is too complicated. JuMP likes many variables and many constraints that are sparse, not a single complicated expression.

As a different approach, you might be able to write this as a user-defined operator:

I will also caution you that using @eval is almost certainly the wrong thing to do. How and why are you constructing functions like this?

3

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])
4

I didn’t run your code, but consider something like:

using JuMP
model = Model()
@variable(model, qq[1:9])
@operator(model, op_Terms, 9, (x...) -> Terms(collect(x)))
@expression(model, op_Terms(qq...))

Note that JuMP now supports complementarity problems natively. You don’t need to use Complementarity.jl. See Example: mixed complementarity problems · JuMP