# Solve Constrained Nonlinear System (?)

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 `kph`simultaneously 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.

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

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

1 Like

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.

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

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
b = x
Tf = x
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})
``````