Solve Constrained Nonlinear System (?)

General Usage
1

I’m a little sad, because it seems like we don’t have any packages for solving nonlinear equations (not odes, not pdes) on a domain range. An equivalent function in the MATLAB world is vpasolve() in which I can feed it a vector of functions, and I can specify the range in which the solution is. The problem here is that NLSolve.jl cannot specify ranges, and the functions I feed it depend on other functions with things like square roots and when I start introducing complex numbers the solver breaks. I seemingly also can’t use JuMP.jl, as the functions I am using can’t accept the abstract datatype that @objective() imposes on the function I feed it.

Here is my example. I am trying to solve for an adiabatic flame temperature, and I have enthalpy and equilibrium constant table lookup functions. My setup is thus:


kpc_lhs(T) = kp(10,T)/kp(9,T); # Eq. 1
kph_lhs(T) = kp(4,T); # Eq. 3
kpc_rhs(a,b) = a / ( (4-a)*sqrt(5.475-0.5*(a+b)) ); # Eq. 2
kph_rhs(a,b) = b / ( (5-a)*sqrt(5.475-0.5*(a+b)) ); # Eq. 4

ΔH(a,b,T) = a*(hs("CO2",T)-94.054) + (4-a)*(hs("CO",T)-26.42) + b*(hs("H2O",T)-57.798) + (5-a)*hs("H2",T) + (5.475-0.5*a-0.5*b)*hs("O2",T) + 28.1209*hs("N2",T) + 30.065; # = 0
kpc(a,b,T) = kpc_lhs(Tf) - kpc_rhs(a,b) # Eq. 1 - Eq. 2 = 0
kph(a,b,c) = kph_lhs(Tf) - kph_rhs(a,b); # Eq. 3 - Eq. 4 = 0

where I need to solve ΔH, kpc, and kphsimultaneously for a, b, and T. In the past I tried using Roots.jl, and solving one of the k equations for a or b, but this gets complex quickly when doing this over and over. As you can see, a+b needs to remain smaller than 10.95 per the formulation in Eq. 2 and Eq. 4. I know their values are between 1 and 5, so I need my solver to stay within these bounds.

I can’t be the only one who needs constraints in problems like these, and as much as I’d like to turn towards something like JuMP, I am unable to due to my enthalpy lookup functions having no clue what to do with what it places in for T. Moreover, JuMP is not a solving library but an optimization library. I’m hoping I’ve just missed a package out there, if anyone knows of something capable of performing this.

2

If you let c^2 = 5.475 - 0.5(a+b), couldn’t you change variables in your solver from (a,b) to (a,c) by letting b = 10.95 - a - 2c^2? That way a and c can be unconstrained real numbers and you only take sqrt(c^2), which is real.

I agree that it would be nice to add box constraints to NLsolve, however. Doesn’t seem like it would be too hard to add this feature to trust-region Newton, since you can just intersect the trust region with the box?

4 Likes
3

That is T? Do you have a minimal working example?

1 Like
4

In principle this is possible using IntervalRootFinding.jl, but I guess only if hs is a Julia function.
Please make a minimal example that we can run to define the functions.

5

See my answers to a very similar question here: Nonlinear System of Equations with Bounds/Constraints on Unknowns - #16 by dpo

6

Unfortunately, this package seems to not like my functions either when I set it up as you suggest:

using NLPModels, NLPModelsIpopt
function c(x)
    a = x[1]
    b = x[2]
    Tf = x[3]
    F = [kpc_lhs(Tf) - kpc_rhs(a,b); # Eq. 1 - Eq. 2 = 0
        kph_lhs(Tf) - kph_rhs(a,b); # Eq. 3 - Eq. 4 = 0
        ΔH(a,b,Tf) ]; # = 0 
    return F
end

x0 = [3,4,2100]

model = ADNLPModel(x -> 0.0, x0; c=c, lcon=[1,1,2000], ucon=[5,5,2200])

stats = ipopt(model)

with:

MethodError: no method matching createProblem(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Array{Int64,1}, ::Array{Int64,1}, ::Int64, ::Int64, ::NLPModelsIpopt.var"#eval_f#3"{ADNLPModel}, ::NLPModelsIpopt.var"#eval_g#4"{ADNLPModel,Int64}, ::NLPModelsIpopt.var"#eval_grad_f#5"{ADNLPModel}, ::NLPModelsIpopt.var"#eval_jac_g#6"{ADNLPModel}, ::NLPModelsIpopt.var"#eval_h#7"{ADNLPModel})

Thank you for your suggestion!

7

This is an interesting approach, I’ll have to play around with it further. I should have known there was a better mathematical path. Thank you.

8

Is the issue that your x0, lcon and ucon are integer? Try making them float vectors.

9

Ah, that may have been it. It does run now, but oddly the solutions it is returning are outside the bounds given, so I will have to look at it further and see what the issue is. It also never converges, so that may be it (although I’m not surprised, everything is so nonlinear with each variable that it makes it tough. Graphical inspection with limited iteration seems to be how this is solved in the past).

10

Ipopt doesn’t guarantee that general constraints are satisfied at intermediate iterates. If that’s a requirement for you, you may want to try Knitro, which has a feasible mode. You can download an evaluation version and give it a go (see the link I posted above). Also, those solvers assume your problem data is C1. Is that the case?