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!