Convert SymPy.jl equation system to JuMP.jl Model

Specific Domains Optimization (Mathematical)
1

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)

add_nonlinear_constraint(model, :($(expr1) == 0))
add_nonlinear_constraint(model, :($(expr2) == 0))

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!

2

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?

3

Hi @odow,

thank you for the fast and helpful answer.

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.

Regarding your question: Can it be done using JuMP instead?
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.

Thank you in advance!

4

Why not just include the nonlinear equality constraints directly in your model?

Note that you should use @constraint instead of @NLconstraint. See ANN: JuMP v1.15 is released and JuMP 1.15.0 is released | JuMP

5

Thank you for the links regarding @constraint and @NLconstraint, I will use @constraint from now on.
I have trouble understanding what you mean by:

include the nonlinear equality constraints directly in your model

Do you mean while setting the variables in @variable(model, x[1:16] >= 0)?

Thanks!

6

No, I meant just write then like:

@constraint(model, x[1] * x[2] - 0.5 == 0)
@constraint(model, x[2] - 1 == 0)

But it depends on how you are using SymPy.