# Calling Mathematica into Julia to symbolically solve a system of non-linear equations
**URL:** https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518
**Category:** Specific Domains
**Tags:** question, nonlinear, symbolics, mathematica
**Created:** [May 5, 2022, 1:51am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518 "2022-05-05T01:51:26Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@MAC](https://discourse.julialang.org/u/MAC)
#### Post date: [May 5, 2022, 1:51am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/1 "2022-05-05T01:51:26Z")
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:
```julia
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:
```julia
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:
```julia
{{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?
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [May 5, 2022, 1:58am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/2 "2022-05-05T01:58:49Z")
Seems to me you are simply solving a system of linear equations here. You can use [`JuMP`](https://jump.dev) for that.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 2:22am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/3 "2022-05-05T02:22:15Z")
> [@amrods](#):
>
> Seems to me you are simply solving a system of linear equations here.
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](https://www.juliahomotopycontinuation.org/) 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”?).
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 5, 2022, 2:51am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/4 "2022-05-05T02:51:53Z")
> [@MAC](#):
>
> 10 equations with 8 variables
This is over-determined, so may not have a solution, but clearly it does since you provide one.
> [@amrods](#):
>
> You can use [`JuMP`](https://jump.dev) for that.
Here is how you might do it with JuMP and the [Gurobi](https://github.com/jump-dev/Gurobi.jl) solver:
> **JuMP Example**
>
> ```julia
> 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)
>
> ```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 2:57am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/5 "2022-05-05T02:57:35Z")
> [@jd-foster](#):
>
> This is over-determined, so may not have a solution, but clearly it does since you provide one.
~~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
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.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 5, 2022, 3:01am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/6 "2022-05-05T03:01:14Z")
> [@stevengj](#):
>
> There are 8 equality constraints in 8 variables.
> [@stevengj](#):
>
> `System of length 10`
?
HomotopyContinuation.jl is great too.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 3:03am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/7 "2022-05-05T03:03:40Z")
> [@stevengj](#):
>
> System of length 10
Whoops, right, I miscounted commas when scanning the equations by eye. Luckily, Julia’s parser is better at counting than me.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 3:05am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/8 "2022-05-05T03:05:33Z")
> [@jd-foster](#):
>
> Here is how you might do it with JuMP and the [Gurobi](https://github.com/jump-dev/Gurobi.jl) solver:
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.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 5, 2022, 3:14am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/9 "2022-05-05T03:14:46Z")
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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 3:20am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/10 "2022-05-05T03:20:25Z")
> [@jd-foster](#):
>
> This is finding a (slightly better polished) point that matches the given one.
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?)
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 5, 2022, 3:22am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/11 "2022-05-05T03:22:41Z")
> [@stevengj](#):
>
> it’s taking an approximate solution to start with
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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 3:26am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/12 "2022-05-05T03:26:57Z")
Since this is nonconvex, is Gurobi guaranteed to find the solutions if they exist?
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 5, 2022, 3:29am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/13 "2022-05-05T03:29:58Z")
> [@jd-foster](#):
>
> `JuMP.set_optimizer_attribute(model, "NonConvex", 2)`
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.
---
### Author:  [@MAC](https://discourse.julialang.org/u/MAC)
#### Post date: [May 5, 2022, 11:56am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/14 "2022-05-05T11:56:21Z")
@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.
```julia
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]
```
```julia
{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}}}
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 5, 2022, 12:15pm UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/15 "2022-05-05T12:15:24Z")
> [@MAC](#):
>
> The benchmarking is helpful too as I will run the system over a very large parameter set.
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.
---
### Author:  [@MAC](https://discourse.julialang.org/u/MAC)
#### Post date: [May 10, 2022, 11:43pm UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/16 "2022-05-10T23:43:31Z")
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).
```julia
#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;
```
```julia
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))
```
```julia
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))
```
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 11, 2022, 12:54am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/17 "2022-05-11T00:54:02Z")
> [@MAC](#):
>
> ```julia
> model = JuMP.Model(Gurobi.Optimizer)
> JuMP.set_optimizer_attribute(model, "NonConvex", 2)
>
> ```
Hi Megan,
Add this after the above:
```julia
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:**
>
> 1.
>
> ```julia
> 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
>
> ```
>
> 1.
>
> ```julia
> 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
>
> ```
>
> 1.
>
> ```julia
> X10 => 0.199952
> X22 => 0.158478
> X11 => 0.16627
> X20 => 0.00999969
> M2 => 3.16965
> X21 => 0.831523
> X12 => 0.633778
> M1 => 8.31549
>
> ```
---
### Author:  [@MAC](https://discourse.julialang.org/u/MAC)
#### Post date: [May 11, 2022, 12:59am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/18 "2022-05-11T00:59:17Z")
Thank you for this.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [May 11, 2022, 1:14am UTC](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/19 "2022-05-11T01:14:29Z")
For future reference, here is the full example:
```julia
#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
```