# Convert SymPy.jl equation system to JuMP.jl Model

Hello together,

I am new to Julia and ran into a problem that I am stuck with for quite some while now.

I want to convert equation systems from SymPy into JuMP Models. However, the symbolic variables and equations are dynamic so variables cannot be hardcoded.

I have found some discussion regarding SymPy and JuMP exchange here. However, I was not able to find a solution for my proble with that.

Here is a small example:

``````using SymPy
using JuMP
using Ipopt

@syms a, b
eq = [a*b-0.5, b - 1]

model = Model(Ipopt.Optimizer)
@variable(model, x[1:2])

var_map = Dict(a => x[1], b => x[2])

expr1 = :(\$(var_map[a]) * \$(var_map[b]) - 0.5)
expr2 = :(\$(var_map[b]) - 1)

optimize!(model)
solution = value.(x)
``````

How can I create expr1 and expr2 automatically from eq. Is there a way?
Also, I am not using NLsolve, because I found that incorporating additional constraints appears to be difficult.

Thanks a lot!

Hi @lufu, welcome to the forum.

You can use:

``````julia> using JuMP

julia> import Ipopt

julia> import SymPy

julia> SymPy.@syms(a, b)
(a, b)

julia> eq = [a * b - 0.5, b - 1]
2-element Vector{Sym{PyCall.PyObject}}:
a⋅b - 0.5
b - 1

julia> model = Model(Ipopt.Optimizer);

julia> set_silent(model)

julia> @variable(model, x[1:2]);

julia> for expr in eq
@constraint(model, SymPy.lambdify(expr, [a, b])(x...) == 0)
end

julia> optimize!(model)

julia> solution = value.(x)
2-element Vector{Float64}:
0.5000000000188031
1.0
``````

But a better question might be: what are you using SymPy for? Can it be done using JuMP instead?

Hi @odow,

Following your example, I was able to dynamically add equations as constraints. However, I experience convergence problems which are probably related to the use of @constraint and my non-linear equation system. If I use @NLconstraint, the convergence problems are gone as you can see in the code below, where I manually replaced my variables with x[1:16]. Unfortunately, `@NLconstraint(model, SymPy.lambdify(expr, [a, b])(x...) == 0)` does not appear to be possible. The numbers like `19076.7467277171` come to be, because some variables have been dynamically replaced with fixed values.

``````using JuMP
using Ipopt

model = Model(Ipopt.Optimizer)
@variable(model, x[1:16] >= 0)

x_init = [
5000, 5000, 5000, 4700, 4000, 4000, 5000, 5000,
5000, 300, 300, 2200, 2200, 11e6, 11e6, 11e6
]

for i in 1:16
set_start_value(x[i], x_init[i])
end

@NLconstraint(model, -x[1] - x[5] - x[6] - x[9] + 19076.7467277171 == 0)
@NLconstraint(model, 4062.12165255267 - x[1] == 0)
@NLconstraint(model, x[2] + x[11] - x[9] == 0)
@NLconstraint(model, 6728.17631205133 - x[2] == 0)
@NLconstraint(model, x[3] + x[10] - x[8] == 0)
@NLconstraint(model, 3255.58343367359 - x[3] == 0)
@NLconstraint(model, 5030.86532943948 - x[7] == 0)
@NLconstraint(model, -x[5] - x[11] - x[6] - x[10] + x[8] + x[7] == 0)
@NLconstraint(model, -x[4] - x[10] + x[7] == 0)
@NLconstraint(model, 2260 * x[2] - x[14] == 0)
@NLconstraint(model, 2260 * x[3] - x[15] == 0)
@NLconstraint(model, 11369755.6445332 - x[16] == 0)
@NLconstraint(model, 2260 * x[1] - x[14] + 826349.638741093 == 0)
@NLconstraint(model, x[12] * x[2] - x[15] + 600594.104706997 == 0)
@NLconstraint(model, x[13] * x[3] - x[16] + 424914.455821467 == 0)

optimize!(model)

solution = value.(x)

println("Solution: ", solution)
``````

Regarding your question: what are you using SymPy for?
I am currently working on a project that requires to dynamically:

• Solve non-linear equation systems analytically for variables
• Replace variable with fixed values or fix variables at a specific value
• Numerically solve non-linear equation systems given initial values and constraints.

SymPy appeared to be the easiest way to achieve the first two points.

Is there something equivalent to SymPy.solve in JuMP to solve non-linear equation systems analytically for variables? I have though a bit about `Symbolics.solve_for` but I am not so sure about the support for non-linear equation systems. If the exchange from JuMP variable to SymPy variable is easier I could also programm with JuMP and then only convert the variables to SymPy for the analytical solver.