Using JuMP for solving a symbolic algebraic system

Specific Domains Optimization (Mathematical)
,
1

Hi All,

I’d like to solve a system of algebraic equations with the flexibility to specify one variable as independent and one as the dependent variable. In the toy example below, I am calculating the Reynolds number, Re but as shown in the output below, the solver did not converge.

Questions:

  1. Is JuMP the best package for this? I tried Symbolics.jl previously but I received an error indicating the solver is only compatible with integers and rational values, which is not applicable to my system.

  2. Is there a way I can see more information about the termination status?

  3. Did the model not solve because the equations are non-linear and I need to use a different type of optimizer?

Thank you!

Code

# Sample problem to solve for variables associated with Reynolds number
using JuMP, Optim

# Set up actual values
ρ_val = 1000.0;
u_val = 1.0;
L_val = 1.0;
μ_val = 1.67e-3;
Re_val = 5.98e+5; # Expected answer


model = Model(Optim.Optimizer)

@variable(model, Re >= 0)
@variable(model, ρ == ρ_val)
@variable(model, u == u_val)
@variable(model, L == L_val)
@variable(model, μ == μ_val)
@objective(model, Min, 0)
@constraint(model, c1, Re == ρ * u * L / μ)

println(model)

optimize!(model)

if !is_solved_and_feasible(model)
    error("Solver did not find an optimal solution. Termination status: $(termination_status(model))")
end

println("Solver found an optimal solution")
println("Objective value: $(objective_value(model))")
println("Re: $(value.(Re))")

Output

ERROR: Solver did not find an optimal solution. Termination status: OTHER_ERROR
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] top-level scope
   @ c:\Users\jlym\ReynoldsNumberToy\ReynoldsNumberToy.jl:28
2

FYI, I solved the issue by switching to Ipopt.jl. For my understanding, I would still be interested to learn why I received the error above from Optim.jl.

Updated Code

# Sample problem to solve for variables associated with Reynolds number
using JuMP, Ipopt

# Set up actual values
ρ_val = 1000.0;
u_val = 1.0;
L_val = 1.0;
μ_val = 1.67e-3;
Re_val = 5.98e+5; # Expected answer


model = Model(Ipopt.Optimizer)

@variable(model, Re >= 0)
@variable(model, ρ == ρ_val)
@variable(model, u == u_val)
@variable(model, L == L_val)
@variable(model, μ == μ_val)
@objective(model, Min, 0)
@constraint(model, c1, Re == ρ * u * L / μ)

println(model)

optimize!(model)

if !is_solved_and_feasible(model)
    error("Solver did not find an optimal solution. Termination status: $(termination_status(model))")
end

println("Solver found an optimal solution")
println("Objective value: $(objective_value(model))")
println("Re: $(value.(Re))")

Output

Min 0
Subject to
 c1 : Re - (((ρ*u) * L) / μ) == 0
 ρ == 1000
 u == 1
 L == 1
 μ == 0.00167
 Re >= 0

This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

Number of nonzeros in equality constraint Jacobian...:        1
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:        1
                     variables with only lower bounds:        1
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 5.99e+05 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 0.00e+00 0.00e+00  -1.0 5.99e+05    -  1.65e-08 1.00e+00f  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   0.0000000000000000e+00    0.0000000000000000e+00


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 0.026

EXIT: Optimal Solution Found.
Solver found an optimal solution
Objective value: 0.0
Re: 598802.3952095808
3

Ipopt is probably the recommended option here. There are many reasons why Optim failed; it is less robust, particularly with a bunch of fixed variables like this.

1 Like