Calling Mathematica into Julia to symbolically solve a system of non-linear equations

Specific Domains
, , ,
19

For future reference, here is the full example:

#Parameters
kb1 =  0.01;
K1  =  50.0;
d1  =  1.0;
l11 =  1.0;
m11 =  10.0;
l12 =  1.0;
m12 =  1.0;
kb2 =  0.01;
K2  =  20.0;
d2  =  1.0;
l22 =  5.0;
m22 =  1.0;
l21 =  10.0;
m21 =  1.0;

using JuMP
using Gurobi

model = JuMP.Model(Gurobi.Optimizer)
JuMP.set_optimizer_attribute(model, "NonConvex", 2)
JuMP.set_optimizer_attribute(model, "PoolSearchMode", 2)
JuMP.set_optimizer_attribute(model, "PoolSolutions", 100)

V = @variables(model, begin
    X10
    M1 >= 0
    X11
    M2 >= 0
    X12
    X20
    X22
    X21
end)

C = @constraints(model, begin
    m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0
    K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0
    l11*M1*X10 - m11*X11 == 0
    m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0
    l12*M2*X10 - m12*X12 == 0
    m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0
    l22*M2*X20 - m22*X22 == 0
    l21*M1*X20 - m21*X21 == 0
    X10 + X12 + X11 == 1
    X20 + X21 + X22 == 1
end)

optimize!(model)

#Seems to find only 1 fixed point M1=>40, M2=>0
result_dict = Dict(V .=> value.(V))

result_dicts = [Dict(V .=> value.(V; result=i)) for i in 1:result_count(model)]

"Determine if values of two `Dict`s compared on their common keys are all within `atol`"
function diff_dicts(d1::Dict, d2::Dict; atol=1e-4)

    different_keys = []

    for k in intersect(keys(d1), keys(d2))
        if abs(d1[k] - d2[k]) > atol
            push!(different_keys, k)
        end
    end

    return isempty(different_keys)
end

candidates = [result_dicts[1]]

for dc1 in result_dicts
    if dc1 in candidates
        continue
    end
    if all([!diff_dicts(dc1, dc2) for dc2 in candidates])
        push!(candidates, dc1)
    end
end

display(candidates) # gives 3 different solutions
2 Likes