# Usage

This section provides details on how to instantiate and solve a constrained least-squares problem with `Enlsip.jl`

As a reminder from Home, problems to solve are of the following form:

\[\begin{aligned} \min_{x \in \mathbb{R}^n} \quad & \dfrac{1}{2} \|r(x)\|^2 \\ \text{s.t.} \quad & c_i(x) = 0, \quad i \in \mathcal{E} \\ & c_i(x) \geq 0, \quad i \in \mathcal{I}, \\ \end{aligned}\]

with:

- the residuals $r_i:\mathbb{R}^n\rightarrow\mathbb{R}$ and the constraints $c_i:\mathbb{R}^n\rightarrow\mathbb{R}$ assumed to be $\mathcal{C}^2$ functions;
- norm $\|\cdot\|$ denoting the Euclidean norm.

Note that with this formulation, bounds constraints are not distinguished from general inequality constraints. Though, for ease of use, they can be provided directly as vectors of lower and/or upper bounds (see Instantiate).

It should be borne in mind, however, that the method implemented in `Enlsip.jl`

has been conceived for nonlinear problems, as there is no other assumption made about the nature of the residuals and constraints functions, apart from being two-time continuously differentiable. The algorithm can still be used to solve linear least squares subject to linear constraints but it will not be as effective as other software where those aspects are taken into account in the design of the optimization method.

## Instantiate a model

Solving a problem with Enlsip is organized in two steps.

First, a model of type `CnlsModel`

must be instantiated.

The `CnlsModel`

constructor requires the evaluation functions of residuals, constraints, their associated jacobian matrices and dimensions of the problem.

Although the package enables one to create linear unconstrained least squares, it is recommended to use it on problems with nonlinear residuals and general constraints.

The three following positional arguments are mandatory to create a model:

`residuals`

: function that computes the vector of residuals`nb_parameters`

: number of variables`nb_residuals`

: number of residuals

The following keyword arguments are optional and deal with constraints and Jacobian matrices. If the Jacobian matrix functions are not provided, they are computed numerically by forward differences using automatic differenciation^{[Backend]}.

Argument | Details |
---|---|

`starting_point` | initial solution (can be an infeasbile point) |

`jacobian_residuals` | function computing the Jacobian matrix of the residuals |

`eq_constraints` | function computing the equality constraints |

`jacobian_eqcons` | function computing the Jacobian matrix of the equality constraints |

`nb_eqcons` | number of equality constraints |

`ineq_constraints` | function computing the inequality constraints |

`jacobian_ineqcons` | function computing the Jacobian matrix of the inequality constraints |

`nb_ineqcons` | number of inequality constraints |

`x_low` | vector of lower bounds |

`x_upp` | vector of upper bounds |

It is assumed that the the different functions passed as arguments of the `CnlsModel`

constructor are called `f(x)`

, where `x`

is a vector of `nb_parameters`

elements and `f`

is one of the functions `residuals`

, `eq_constraints`

, `jacobian_eqcons`

etc.

## Solving

Then, the `Enlsip`

solver can be used by calling the `solve!`

function on an instantiated model. One can pass keywords arguments to modify some options of the algorithm or adjust the tolerances. For the expression of the termination criteria, see Method page (Stopping).

`Enlsip.solve!`

— Function`solve!(model{T})`

Once a `CnlsModel`

has been instantiated, this function solves the optimzation problem associated by using the method implemented in the `Enlsip`

solver.

Options:

`silent::Bool`

Set to

`false`

if one wants the algorithm to print details about the iterations and termination of the solverDefault is

`true`

, i.e. by default, there is no output. If one wants to print those information afert solving, the`print_cnls_model`

method

can be called

`max_iter::Int`

Maximum number of iterations allowed

Default is

`100`

`scaling::Bool`

Set to

`true`

if one wants the algorithm to work with a constraints jacobian matrix whose rows are scaled (i.e. all constraints gradients vectors are scaled)Default is

`false`

`time_limit::T`

Maximum elapsed time (i.e. wall time)

Default is

`1000`

Tolerances:

`abs_tol::T`

Absolute tolerance for small residuals

Default is

`eps(T)`

`rel_tol::T`

Relative tolerance used to measure first order criticality and consistency

Default is

`sqrt(eps(T))`

`c_tol::T`

Tolerance used to measure feasability of the constraints

Default is

`sqrt(eps(T))`

`x_tol::T`

Tolerance used to measure the distance between two consecutive iterates

Default is

`sqrt(eps(T))`

Diagnosis of the conduct of the algorithm can be printed by either setting the `silent`

keyword argument of the function `solve!`

to `false`

or by calling `print_cnls_model`

after solving. Here are some details on how to read and understand the different columns of the output:

Column | Description |
---|---|

`iter` | iteration number |

`objective` | value of the sum of squared residuals (i.e. objective function) at current point |

$\vert\vert$ `active_constraints` $\vert\vert^2$ | value of the sum of squared active constraints at current point |

$\vert\vert$ `p` $\vert\vert$ | norm of the search direction computed at current iteration |

$\alpha$ | value of the steplength computed at current iteration |

`reduction` | reduction in the objective function performed after moving to the next iterate |

One can get additional info about termination of the algorithm by calling one of the following functions:

Name |
---|

`solution` |

`status` |

`constraints_values` |

`sum_sq_residuals` |

`Enlsip.solution`

— Function`solution(model)`

Once the given `model`

has been solved, this function returns the optimal solution, or last solution obtained if no convergence, as a `Vector`

of approriate dimension.

`Enlsip.status`

— Function`status(model)`

This functions returns a `Symbol`

that gives brief information on the solving status of `model`

.

If a model has been instantiated but the solver has not been called yet, it will return `:unsolved`

.

Once the solver has been called and if a first order stationary point satisfying the convergence criteria has been computed, it will return `:found_first_order_stationary_point`

.

If the algorithm met an abnormall termination criteria, it will return one of the following:

`:failed`

: the algorithm encoutered a numerical error that triggered termination`:maximum_iterations_exceeded`

: a solution could not be reached within the maximum number of iterations`:time_limit_exceeded`

: the algorithm stopped because solving time exceeded the time limit

`Enlsip.constraints_values`

— Function`constraints_values(model)`

Computes values of all the constraints in `model`

at the solution.

The vector returned is the concatenation of equalities, inequalities and box constraints (in that order).

For instance, let `xₛ`

be the solution found. If functions `h`

, `g`

compute equality and inequality constraints and `xₗ, xᵤ`

are vectors of lower and lower bounds, it will return `[h(xₛ); g(xₛ); xₛ-xₗ; xᵤ-xₛ]`

.

If one wants to compute each type of constraints seperately, see `equality_constraints_values`

, `inequality_constraints_values`

and `bounds_constraints_values`

.

`Enlsip.sum_sq_residuals`

— Function`sum_sq_residuals(model)`

Once the given `model`

has been solved, returns the value of the objective function, i.e. sum of squared residuals functions, computed at the optimal solution. If no convergence, this value is computed at the last solution obtained.

## Examples

### Problem 65 from Hock and Schittkowski collection^{}[HS80]

We show how to implement and solve the following problem:

\[\begin{aligned} \min_{x_1, x_2, x_3} \quad & (x_1-x_2)^2 + \dfrac{(x_1+x_2-10)^2}{9}+(x_3-5)^2 \\ \text{s.t.} \quad & 48-x_1^2-x_2^2-x_3^2 \geq 0\\ & -4.5\leq x_i \leq 4.5, \quad i=1,2\\ & -5 \leq x_3 \leq 5. \end{aligned}.\]

The expected optimal solution is $(3.650461821, 3.65046168, 4.6204170507)$.

Associated value of objective function equals $0.9535288567$.

First, we provide the dimensions of the problems.

```
# Dimensions of the problem
n = 3 # number of parameters
m = 3 # number of residuals
```

Then, we define the functions required to compute the residuals, constraints, their respective jacobian matrices and a starting point. In this example, we use the starting point given in the reference^{[HS80]}, i.e. $(-5, 5, 0)$

```
# 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] # evaluation function for the equality constraint
jac_c(x::Vector) = [ -2x[1] -2x[2] -2x[3]] # Jacobian matrix of the equality constraint
x_l = [-4.5, -4.5, -5.0] # lower bounds
x_u = [4.5, 4.5, 5.0] # upper bounds
# Starting point
x0 = [-5.0, 5.0, 0.0]
```

A `CnlsModel`

can now be instantiated.

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

Finally, the `solve!`

function can be called on our model. In this example, keyword arguments remain to default values.

`Enlsip.solve!(hs65_model)`

Once `Enlsip`

solver has been executed on our problem, a summary of the conduct of the algorithm can be printed by calling `print_cnls_model`

.

`Enlsip.print_cnls_model(hs65_model)`

```
****************************************************************
* *
* Enlsip.jl v0.9.7 *
* *
* This is the Julia version of the ENLSIP algorithm, initially *
* conceived and developed in Fortran77 by Per Lindstrom and *
* Per Ake Wedin from the Institute of Information Processing *
* (University of Umea, Sweden). *
* *
****************************************************************
Characteristics of the model
Number of parameters.................: 3
Number of residuals..................: 3
Number of equality constraints.......: 0
Number of inequality constraints.....: 1
Number of lower bounds...............: 3
Number of upper bounds...............: 3
Constraints internal scaling.........: false
Iteration steps information
iter objective ||active_constraints||² ||p|| α reduction
1 1.3611111e+02 4.50e+00 5.00e+00 4.81e-01 1.599e+01
2 1.1589067e+02 4.66e+01 6.84e+00 4.30e-01 2.657e+01
3 7.8319097e+01 6.61e-02 8.19e+00 7.58e-01 7.366e+01
4 4.6650091e+00 2.18e+01 1.44e+00 1.00e+00 5.459e+00
5 9.5479609e-01 4.29e+00 4.66e-01 1.00e+00 4.411e-01
6 9.3767054e-01 4.70e-02 6.67e-02 1.00e+00 1.766e-02
7 9.5320179e-01 1.98e-05 7.92e-03 1.00e+00 3.653e-04
8 9.5352428e-01 3.94e-09 9.97e-04 1.00e+00 5.152e-06
9 9.5352878e-01 9.89e-13 1.25e-04 1.00e+00 8.164e-08
10 9.5352886e-01 2.46e-16 1.57e-05 1.00e+00 1.287e-09
11 9.5352886e-01 6.11e-20 1.98e-06 1.00e+00 2.030e-11
12 9.5352886e-01 1.52e-23 2.48e-07 1.00e+00 3.205e-13
13 9.5352886e-01 3.65e-27 3.12e-08 1.00e+00 4.774e-15
Number of iterations...................: 13
Square sum of residuals................: 9.5352886e-01
Number of function evaluations.........: 152
Number of Jacobian matrix evaluations..: 30
Solving time (seconds).................: 0.080
Termination status.....................: found_first_order_stationary_point
```

If one just wants to know about termination of the algorithm, calling `status`

will tell if the problem has been successfully solved or not.

`Enlsip.status(hs65_model)`

`:found_first_order_stationary_point`

Then, calling `solution`

and `sum_sq_residuals`

will respectively return the optimal solution obtained and the value of objective function at that point.

`hs65_solution = Enlsip.solution(hs65_model)`

```
3-element Vector{Float64}:
3.6504617268193496
3.650461726819378
4.620417552781781
```

`hs65_objective = Enlsip.sum_sq_residuals(hs65_model)`

`0.9535288568047824`

The solution obtained is relatively close to the expected optimal solution, although it differs from more than the tolerance used.

`maximum(abs.(hs65_solution - [3.650461821, 3.65046168, 4.6204170507])) < sqrt(eps(Float64))`

`false`

However, the difference between the objective value obtained with `Enlsip`

and the expected one does not exceed the default tolerance.

`abs(hs65_objective - 0.9535288567) < sqrt(eps(Float64))`

`true`

### Chained Rosenbrock problem^{}[LV99]

We now consider the following problem:

\[\begin{aligned} \min_{x} \quad & \sum_{i=1}^{n-1} 100(x_i^2-x_{i+1})^2 + (x_i-1)^2 \\ \text{s.t.} \quad & 3x_{k+1}^3 + 2x_{k+2}+\sin(x_{k+1}-x_{k+2})\sin(x_{k+1}+x_{k+2})+4x_{k+1} - x_k\exp(x_k-x_{k+1}) - 8 = 0 \\ & k=1,\ldots, n-2, \end{aligned}\]

for a given natural number $n\geq 3$. In this example, we use $n=1000$. Though analytic Jacobians are easy to define, they will be passed as arguments in the model instantiation, so to let the algorithm use automatic differentiation.

```
# Dimensions
n = 1000 # Number of variables
m = 2(n-1) # Number of residuals
nb_eq = n-2 # Number of equality constraints
# Residuals
function r(x::Vector)
n = length(x)
m = 2(n-1)
rx = Vector{eltype(x)}(undef,m)
rx[1:n-1] = [10(x[i]^2 - x[i+1]) for i=1:n-1]
rx[n:m] = [x[k-n+1] - 1 for k=n:m]
return rx
end
# Constraints
function c(x::Vector)
n = length(x)
cx = [3x[k+1]^3 + 2x[k+2] - 5 + sin(x[k+1]-x[k+2])*sin(x[k+1]+x[k+2]) + 4x[k+1] -
x[k]*exp(x[k]-x[k+1]) - 3 for k=1:n-2]
return cx
end
x0 = [(mod(i,2) == 1 ? -1.2 : 1.0) for i=1:n] # Starting point
# Instantiation of the model
cr_enlsip = CnlsModel(r,n,m; starting_point=x0, eq_constraints=c, nb_eqcons=nb_eq)
# Solving
Enlsip.solve!(cr_enlsip)
# Show solving status
Enlsip.status(cr_enlsip)
```

`:found_first_order_stationary_point`

To give an insight on how the performance of the algorithm scales when dimensions increase, we give a table of the solving time obtained with different values of parameter `n`

. For comparison purposes, the calculation times obtained with the Julia version of the IPOPT general solver for nonlinear programming ^{[WB06]} are also indicated. Both solvers were used with their respective default tolerances. With the same starting point, both solvers reach similar local solutions despite different stopping criteria.

We show the modeling of the problem with JuMP and its solving with `IPOPT`

for `n=1000`

.

```
cr_ipopt = Model(Ipopt.Optimizer); set_silent(cr_ipopt)
@variable(cr_ipopt, x[i=1:n], start = x0[i])
for k=1:n-2
@NLconstraint(cr_ipopt, 3x[k+1]^3 + 2x[k+2] - 5 + sin(x[k+1]-x[k+2])*sin(x[k+1]+x[k+2]) + 4x[k+1] - x[k]*exp(x[k]-x[k+1]) - 3 == 0)
end
@NLobjective(cr_ipopt, Min, sum(100(x[i]^2 - x[i+1])^2 + (x[i]-1)^2 for i=1:n-1))
optimize!(cr_ipopt)
termination_status(cr_ipopt)
```

`LOCALLY_SOLVED::TerminationStatusCode = 4`

The two solvers reach similar solutions.

`Enlsip.sum_sq_residuals(cr_enlsip) ≈ objective_value(cr_ipopt)`

`true`

`Enlsip.solution(cr_enlsip) ≈ value.(cr_ipopt[:x])`

`true`

The benchmark^{[BT]} of the solving times (in seconds) for `Enlsip.jl`

and `Ipopt.jl`

is given in the following table:

Value of `n` | `Enlsip.jl` | `IPOPT` |
---|---|---|

`10` | `3.616e-4` | `2.703e-3` |

`100` | `3.322e-2` | `4.759e-3` |

`1000` | `2.325e0` | `2.520e-2` |

`5000` | `3.172e2` | `1.490e-1` |

- Backend
`ForwardDiff.jl`

https://juliadiff.org/ForwardDiff.jl/stable/ - HS80W. Hock and K. Schittkowski.
*Test Examples for Nonlinear Programming Codes*, volume 187 of Lecture Notes in Economics and Mathematical Systems. Springer, second edition, 1980. - LV99L. Lukšan and J. Vlček.
*Sparse and Partially Separable Test Problems for Unconstrained and Equality Constrained Optimization*. Technical report 767, 1999. - BTThe values shown in the table were obtained using the
`@btime`

macro from the BenchmarkTools.jl package. - WB06A. Wächter and L. Biegler,
*On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming*, Mathematical Programming, vol. 106(1), pages 25-57, 2006.