Enlsip.jl

Package Enlsip.jl is the Julia version of ENLSIP, an open source algorithm originally written in Fortran77 and designed to solve nonlinear least-squares problems subject to nonlinear constraints. The optimization method implemented in ENLSIP is described in

Per Lindström and Per-Åke Wedin, Gauss Newton based algorithms for constrained nonlinear least squares problems. Technical Report S-901 87, Institute of Information processing, University of Umeå, Sweden, 1988.

The source code of the Fortran77 library is available at https://plato.asu.edu/sub/nonlsq.html.

Problems that can be solved using Enlsip.jl are modeled as follows:

where:

• the residuals $r_i:\mathbb{R}^n\rightarrow\mathbb{R}$ and the constraints $c_i:\mathbb{R}^n\rightarrow\mathbb{R}$ are assumed to be $\mathcal{C}^2$ functions;
• norm $||\cdot||$ denotes the Euclidean norm;
• $l$ and $u$ are respectively vectors of lower and upper bounds.

In the formulation above, bound constraints are written seperately but they are treated as classical inequality constraints in the method implemented in ENLSIP.

To add Enlsip, use Julia's package manager by typing the following command inside the REPL:

   using Pkg
   Pkg.add("Enlsip")

Solving a problem with Enlsip is organized in two steps.

First, you need to create a model of your problem with the CnlsModel structure.

An object of type CnlsModel can be created using a constructor, whose arguments are the following:

• residuals : function that computes the vector of residuals

• nb_parameters : number of variables

• nb_residuals : number of residuals

• stating_point : initial solution

• jacobian_residuals : function that computes the jacobian matrix of the residuals. If not passed as argument, it is computed numericcaly by forward differences

• eq_constraints : function that computes the vector of equality constraints

• jacobian_eqcons : function that computes the jacobian matrix of the equality constraints. If not passed as argument, it is computed numericcaly by forward differences

• nb_eqcons : number of equality constraints

• ineq_constraints : function that computes the vector of inequality constraints

• jacobian_ineqcons : function that computes the jacobian matrix of the inequality constraints. If not passed as argument, it is computed numericcaly by forward differences

• nb_ineqcons : number of inequality constraints

• x_low and x_upp : respectively vectors of lower and upper bounds

Then, once your model is instantiated, you can call the solve! function to solve your problem.

# Import Enlsip
using Enlsip

# Dimensions of the problem

n = 3 # number of parameters
m = 3 # number of residuals

# Residuals and jacobian matrix associated
r(x::Vector) = [x[1] - x[2]; (x[1]+x[2]-10.0) / 3.0; x[3]-5.0]

jac_r(x::Vector) = [1. -1. 0;
    1/3 1/3 0.;
    0. 0. 1.]

# Constraints (one nonlinear inequality and box constraints)
c(x::Vector) = [48.0 - x[1]^2-x[2]^2-x[3]^2]
jac_c(x::Vector) = [ -2x[1] -2x[2] -2x[3]]
x_l = [-4.5, -4.5, -5.0]
x_u = [4.5, 4.5, 5.0] 

# Starting point 
x0 = [-5.0, 5.0, 0.0]

# Instantiate a model associated with the problem 
hs65_model = Enlsip.CnlsModel(r, n, m ;jacobian_residuals=jac_r, starting_point=x0,
ineq_constraints = c, jacobian_ineqcons=jac_c, nb_ineqcons = 1, x_low=x_l, x_upp=x_u)


# Call of the `solve!` function
Enlsip.solve!(hs65_model)

# Print solution and objective value

println("Algorithm termination status: ", Enlsip.status(hs65_model))
println("Optimal solution: ", Enlsip.solution(hs65_model))
println("Optimal objective value: ", Enlsip.sum_sq_residuals(hs65_model))