# Nonlinear System of Equations with Bounds/Constraints on Unknowns

**URL:** <https://discourse.julialang.org/t/nonlinear-system-of-equations-with-bounds-constraints-on-unknowns/37257>\
**Category:** Optimization (Mathematical)\
**Created:** [April 9, 2020, 12:11am UTC](https://discourse.julialang.org/t/nonlinear-system-of-equations-with-bounds-constraints-on-unknowns/37257 "2020-04-09T00:11:34Z")\
**Posts on this page:** 1\
**Showing post:** 16

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [April 11, 2020, 12:35am UTC](https://discourse.julialang.org/t/nonlinear-system-of-equations-with-bounds-constraints-on-unknowns/37257/16 "2020-04-11T00:35:55Z")

</div>

I apologize for the long answer. Below are 3 experiments.

Here’s a solution based on solving a feasibility problem (minimizing a constant objective function subject to your system as constraints). This allows you to include bounds on your variables. However, it seems to me that your problem is over-constrained: 6 variables and 7 equality constraints if you include w\_1 + w\_2 = 1. I left this last constraint out and it appears to be (approximately) satisfied at the solution found by IPOPT:

```julia
(@v1.4) pkg> activate test-feas

(test-feas) pkg> add NLPModels NLPModelsIpopt

julia> using NLPModels, NLPModelsIpopt

julia> function c(x)
           M = [2.250, 9.675, 57.263, 427.219, 3836.109, 40234.852]
           F = [M[1] - x[5]*x[3]*x[1] - x[6]*x[4]*x[2] ;
                 M[2] - x[5]*x[3]^2*(x[1]+1)*x[1] - x[6]*x[4]^2*(x[2]+1)*x[2] ;
                 M[3] - x[5]*x[3]^3*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^3*(x[2]+2)*(x[2]+1)*x[2] ;
                 M[4] - x[5]*x[3]^4*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^4*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2] ;
                 M[5] - x[5]*x[3]^5*(x[1]+4)*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^5*(x[2]+4)*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2] ;
                 M[6] - x[5]*x[3]^6*(x[1]+5)*(x[1]+4)*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^6*(x[2]+5)*(x[2]+4)*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2] ]
                 #x[5] + x[6] - 1] # left this out
           return F
       end

julia> x0 = [0.5; 0.5; 0.5; 0.5; 0.5; 0.5];

julia> # define a model with derivatives computed by ForwardDiff
julia> # bounds on the variables are included
julia> model = ADNLPModel(x -> 0.0, x0; c=c, lcon=zeros(6), ucon=zeros(6), lvar=zeros(6))

julia> stats = ipopt(model) # solution found after 481 iterations

julia> x = stats.solution
6-element Array{Float64,1}:
 2.0064087361563248
 0.46019054128577414
 1.4990976942166683
 1.0902146315580712
 0.6949587215544324
 0.3183161594144738

julia> x[5] + x[6]
1.0132748809689063

```

Choosing a different starting point might give you the solution you’re looking for.

As an experiment, I tried formulating your problem as an over-determined nonlinear least-squares problem with bounds and include the term w\_1 + w\_2 - 1 in the residual. IPOPT doesn’t manage to solve it but KNITRO does (KNITRO is a commercial solver but you can download a demo version from the Artelys website):

```julia
(test-feas) pkg> add NLPModelsKnitro

julia> using NLPModelsKnitro

julia> function c(x)
           M = [2.250, 9.675, 57.263, 427.219, 3836.109, 40234.852]
           F = [M[1] - x[5]*x[3]*x[1] - x[6]*x[4]*x[2] ;
                 M[2] - x[5]*x[3]^2*(x[1]+1)*x[1] - x[6]*x[4]^2*(x[2]+1)*x[2] ;
                 M[3] - x[5]*x[3]^3*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^3*(x[2]+2)*(x[2]+1)*x[2] ;
                 M[4] - x[5]*x[3]^4*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^4*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2] ;
                 M[5] - x[5]*x[3]^5*(x[1]+4)*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^5*(x[2]+4)*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2] ;
                 M[6] - x[5]*x[3]^6*(x[1]+5)*(x[1]+4)*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1] - x[6]*x[4]^6*(x[2]+5)*(x[2]+4)*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2] ;
                 x[5] + x[6] - 1]
           return F
       end

julia> model = ADNLSModel(c, x0, 7, lvar=zeros(6)) # NSL = nonlinear least squares

julia> stats = knitro(model) # solved in 154 iterations

julia> x = stats.solution
6-element Array{Float64,1}:
 2.0014732680262566
 0.4965607369186837
 1.4997796318745336
 1.0165903327466042
 0.6988962139789188
 0.30110695603014176

```

Finally, you can try one of our pure Julia solvers that handles bound constraints:

```julia
(test-feas) pkg> add JSOSolvers
julia> using JSOSolvers

julia> stats = tron(model) # solved in 1006 iterations

julia> stats.solution
6-element Array{Float64,1}:
 1.9441520958770582
 1.9441520958770582
 1.5104044760203716
 1.5104044760203716
 0.36776275007834164
 0.36776275007834164

```

---

_[View the full topic](https://discourse.julialang.org/t/nonlinear-system-of-equations-with-bounds-constraints-on-unknowns/37257)._