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

Specific Domains
, , ,
1

I’m currently using Mathematica to symbolically solve a system of non-linear ODEs (10 equations with 8 variables) as I believe this is not yet possible using Symbolics.jl. Is anyone able to comment on whether it would be possible to call Mathematica into Julia using MathLink.jl to solve the below system. The documentation provided a simple example that I wasn’t able to extrapolate to my Mathematica code (below) successfully. Thanks very much.

The system is governed by 14 parameters:

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

and is solved in the following way in Mathematica:

Solve[{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 
           && M1 >= 0 
           && M2 >= 0} , 
           {X10, M1, X11, M2, X12, X20, X22, X21}, Reals]

generating the below fixed points for variables M1 and M2:

{{X10 -> 0.000999999, M1 -> 9.93007, X11 -> 0.993006, M2 -> 29.9701, 
  X12 -> 0.00599402, X20 -> 0.000333333, X22 -> 0.999005, 
  X21 -> 0.000662005}}

How could I obtain this result in Julia?

1 Like
2

Seems to me you are simply solving a system of linear equations here. You can use JuMP for that.

3

It’s nonlinear, because there are terms like M1*X10 that are products of unknowns. However, I don’t see how it is an “ODE” as claimed by @MAC … I guess it comes from the fixed points of an ODE? — it looks like a system of polynomial equations, which can be solved very efficiently by the HomotopyContinuation.jl package if you want a numerical answer (which is what you are getting from Mathematica anyway, it seems … in what way are you solving it “symbolically”?).

4 Likes
4

This is over-determined, so may not have a solution, but clearly it does since you provide one.

Here is how you might do it with JuMP and the Gurobi solver:

JuMP Example 
using JuMP
using Gurobi

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

## Model def.:
model = JuMP.Model(Gurobi.Optimizer)
JuMP.set_optimizer_attribute(model, "NonConvex", 2)

## Variables:
V = @variables(model, begin
    X10
    M1
    X11
    M2
    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 
    M1 >= 0 
    M2 >= 0
end)

optimize!(model)

result_dict = Dict(V .=> value.(V))
primal_feasibility_report(model, result_dict)

# Provided solution:
feas_pt = Dict(X10 => 0.000999999, M1 => 9.93007, X11 => 0.993006, M2 => 29.9701, 
X12 => 0.00599402, X20 => 0.000333333, X22 => 0.999005, 
X21 => 0.000662005)

primal_feasibility_report(model, feas_pt)
4 Likes
5

There are 8 equality constraints in 8 variables. Whoops, I miscounted — 10 equality constraints in 8 variables. There are two inequalities which filter out some (in this case all but one) of the roots.

With HomotopyContinuation.jl I can reproduce the Mathematica solution:

julia> using HomotopyContinuation

julia> kb1 = 0.01;
julia> K1  = 10.0;
julia> d1  = 1.0;
julia> l11 = 1.0;
julia> m11 = 0.01;
julia> l12 = 1.0;
julia> m12 = 5.0;
julia> kb2 = 0.01;
julia> K2  = 30.0;
julia> d2  = 1.0;
julia> l22 = 1.0;
julia> m22 = 0.01;
julia> l21 = 1.0;
julia> m21 = 5.0;

julia> @var X10 M1 X11 M2 X12 X20 X22 X21
(X10, M1, X11, M2, X12, X20, X22, X21)

julia> F = System([m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10, 
                  K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20, 
                  l11*M1*X10 - m11*X11, 
                  m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20, 
                  l12*M2*X10 - m12*X12, 
                  m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20, 
                  l22*M2*X20 - m22*X22 == 0, l21*M1*X20 - m21*X21, 
                  X10 + X12 + X11 - 1, 
                  X20 + X21 + X22 - 1])
System of length 10
 8 variables: M1, M2, X10, X11, X12, X20, X21, X22

 0.01*X11 + 5.0*X12 - 1.0*M1*X10 - 1.0*M2*X10
 -1.0*M1 + 0.01*X10 + 10.01*X11 + 5.0*X21 - 1.0*M1*X10 - 1.0*M1*X20
 -0.01*X11 + 1.0*M1*X10
 -1.0*M2 + 5.0*X12 + 0.01*X20 + 30.01*X22 - 1.0*M2*X10 - 1.0*M2*X20
 -5.0*X12 + 1.0*M2*X10
 5.0*X21 + 0.01*X22 - 1.0*M1*X20 - 1.0*M2*X20
 0
 -5.0*X21 + 1.0*M1*X20
 -1 + X10 + X11 + X12
 -1 + X20 + X21 + X22

julia> result = solve(F)
Tracking 6 paths... 100%|███████████████████████████████| Time: 0:00:16
  # paths tracked:                  6
  # non-singular solutions (real):  4 (4)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         4 (4)
Result with 4 solutions
=======================
• 6 paths tracked
• 4 non-singular solutions (4 real)
• random_seed: 0x86ff6d6b
• start_system: :polyhedral

julia> real_solutions(result) # extract the real roots
4-element Vector{Vector{Float64}}:
 [-1.0070462962057125e-5, 29.990003354585465, 0.142918523818753, -0.00014392557011470473, 0.8572254017513615, 0.00033333329628380104, -6.713646621069738e-10, 0.9996666673750833]
 [9.930069784028507, 29.97014319708603, 0.0009999989929586964, 0.9930059784038578, 0.005994022603182795, 0.0003333332962592715, 0.0006620045786384862, 0.9990046621251022]
 [9.9900100166733, -3.3366674193768253e-6, 0.0009999989990009555, 0.9990000016683313, -6.67332813376323e-10, 0.33359260855930595, 0.6665187001991278, -0.00011130875843388742]
 [-1.0009999973788885e-5, -3.3344444459124908e-6, 1.001002671228785, -0.0010020036712260715, -6.675575594819791e-7, 1.0003355590062681, -2.002671784362474e-6, -0.00033355633448375795]

julia>  filter(x -> x[1] >= 0 && x[2] >= 0, real_solutions(result)) # extract the M1 >= 0, M2 >= solutions
1-element Vector{Vector{Float64}}:
 [9.930069784028507, 29.97014319708603, 0.0009999989929586964, 0.9930059784038578, 0.005994022603182795, 0.0003333332962592715, 0.0006620045786384862, 0.9990046621251022]

Note that the final line is the solution for M1, M2, X10, X11, X12, X20, X21, X22 (== F.variables), and seems to match what you are reporting from Mathematica.

9 Likes
6

?

HomotopyContinuation.jl is great too.

7

Whoops, right, I miscounted commas when scanning the equations by eye. Luckily, Julia’s parser is better at counting than me.

1 Like
8

I’m curious how long Gurobi took to solve this? HomotopyContinuation took 16 seconds (though it was finding all the roots, and I only filtered them a posteriori). Looks like this is mostly compilation time; solving a second time is about 0.02 seconds.

9

Initial:
7.042219 seconds (19.25 M allocations: 1.061 GiB, 5.44% gc time, 99.02% compilation time)
After compiliation:
0.058017 seconds (1.80 k allocations: 96.125 KiB)

This is finding a (slightly better polished) point that matches the given one.

10

Oh, so it’s not really comparable — it’s taking an approximate solution to start with, whereas HomotopyContinuation is finding the solution from scratch?

(@MAC, I wonder how HomotopyContinuation’s performance compares to Mathematica here?)

11

Sorry if that wasn’t clear, it is starting from scratch, just the result satisfies the constraints to a smaller tolerance than the provided one.

2 Likes
12

Since this is nonconvex, is Gurobi guaranteed to find the solutions if they exist?

13

Yes, this setting in the newer versions of Gurobi guarantees global lower bounds and optimality for bilinear problems. That said, this is a feasbility problem. I also ran it asking it to find all available solutions, and it returns this point only. It becomes much harder if we drop the lower bounds M1 >= 0, M2 >= 0 and filter a posteriori. I think HomotopyContinuation is more impressive here in finding multiple solutions.

3 Likes
14

@jd-foster and @stevengj: Thanks so much for all these great responses - it’s my first time posting here and it has been very helpful.

In the Mathematica code the system is solved exactly and the results are subsequently numericized - sorry for the confusion. A numerical solution will still work for me so these examples are excellent. The benchmarking is helpful too as I will run the system over a very large parameter set.

@stevengj in Mathematica this example takes on average 0.0097s to solve.

RepeatedTiming[  
Solve[{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 
&& M1 >= 0 && M2 >= 0} , 
{X10, M1, X11, M2, X12, X20, X22, X21}, Reals], 30]

{0.0097, {{X10 -> 0.000999999, M1 -> 9.93007, X11 -> 0.993006, M2 -> 29.9701, X12 -> 0.00599402, X20 -> 0.000333333, X22 -> 0.999005, X21 -> 0.000662005}}}
15

If you are varying parameters and computing roots over and over, the problem becomes much easier, because then you just need to track how solutions evolve. (This is called a “homotopy” or “numerical continuation”. HomotopyContinuation has built-in methods for this.)

Even without taking advantage of this, in HomotopyContinuations it looks like much of the 16 seconds I quote above was compilation time (or some other initialization? I don’t know how much it is caching things internally). Solving again is about 0.02 seconds.

3 Likes
16

As an aside: I ran a different parameter set that should produce 3 fixed points (2 stable and 1 unstable) as determined by Mathematica. Using JuMP with Gurobi it was only able to identify 1 fixed point while HomotopyContinuation was able to find all 3 (unless I have implemented something incorrectly).

#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 HomotopyContinuation

@var X10 M1 X11 M2 X12 X20 X22 X21

F = System([m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10),
                  K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20,
                  l11*M1*X10 - m11*X11,
                  m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20,
                  l12*M2*X10 - m12*X12,
                  m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20,
                  l22*M2*X20 - m22*X22,
                  l21*M1*X20 - m21*X21,
                  X10 + X12 + X11 - 1,
                  X20 + X21 + X22 - 1])

result = HomotopyContinuation.solve(F)

#Finds all 3 fixed points (M1=>40, M2=>0); (M1=>8.32, M2=>3.17); (M1=>0, M2=>19.8)
all_fps = filter(x -> x[1] >= 0 && x[2] >= 0, real_solutions(result))

using JuMP
using Gurobi

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

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))
17

Hi Megan,
Add this after the above:

JuMP.set_optimizer_attribute(model, "PoolSearchMode", 2)
JuMP.set_optimizer_attribute(model, "PoolSolutions", 1000)

where the 1000 in "PoolSolutions" is some arbitrary but large enough number for the number of solutions you expect to capture.
The default in Gurobi is to return just one optimal solution, and since the objective is just 0 here, it returns one feasible point.
With the above, I get 14 solutions, but these are not all different up to small tolerances.

The 'unique' solutions are: 
  X10 => 0.19999
  X22 => 4.14194e-7
  X11 => 0.800003
  X20 => 0.00249363
  M2  => 3.32202e-5
  X21 => 0.997506
  X12 => 6.6437e-6
  M1  => 40.0022

 X10 => 0.0480795
  X22 => 0.989937
  X11 => 3.04079e-6
  X20 => 0.00999995
  M2  => 19.7988
  X21 => 6.32807e-5
  X12 => 0.951917
  M1  => 0.00063283

  X10 => 0.199952
  X22 => 0.158478
  X11 => 0.16627
  X20 => 0.00999969
  M2  => 3.16965
  X21 => 0.831523
  X12 => 0.633778
  M1  => 8.31549
18

Thank you for this.

1 Like
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