# Solve Large-Scale Problem with FletcherPenaltySolver.jl

In this tutorial we use `fps_solve`

to solve a large-scale optimization problem resulting from the discretization of a PDE-constrained optimization problem and compare the solve with Ipopt.

## Problem Statement

Let Ω = (-1,1)², we solve the following distributed Poisson control problem with Dirichlet boundary:

\[ \left\lbrace \begin{aligned} \min_{y \in H^1_0, u \in H^1} \quad & \frac{1}{2} \int_\Omega |y(x) - y_d(x)|^2dx + \frac{\alpha}{2} \int_\Omega |u|^2dx \\ \text{s.t.} & -\Delta y = h + u, \quad x \in \Omega, \\ & y = 0, \quad x \in \partial \Omega, \end{aligned} \right.\]

where yd(x) = -x₁² and α = 1e-2. The force term is h(x₁, x₂) = - sin(ω x₁)sin(ω x₂) with ω = π - 1/8.

We refer to Gridap.jl for more details on modeling PDEs and PDENLPModels.jl for PDE-constrained optimization problems.

`using Gridap, PDENLPModels`

```
WARNING: method definition for SparseMatrixCSR at /juliateam/.julia/packages/Gridap/EZQEK/src/Algebra/SparseMatrixCSR.jl:33 declares type variable Ti but does not use it.
WARNING: method definition for SparseMatrixCSR at /juliateam/.julia/packages/Gridap/EZQEK/src/Algebra/SparseMatrixCSR.jl:33 declares type variable Tv but does not use it.
WARNING: method definition for push_coo! at /juliateam/.julia/packages/Gridap/EZQEK/src/Algebra/SparseMatrixCSR.jl:166 declares type variable Bi but does not use it.
WARNING: could not import Base._rangestyle into ArrayLayouts
WARNING: could not import LinearAlgebra.Abuf into ArrayLayouts
WARNING: could not import LinearAlgebra.Bbuf into ArrayLayouts
WARNING: could not import LinearAlgebra.Cbuf into ArrayLayouts
WARNING: method definition for unsafe_convert at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/ArrayLayouts.jl:90 declares type variable P but does not use it.
WARNING: method definition for unsafe_convert at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/ArrayLayouts.jl:90 declares type variable N but does not use it.
WARNING: method definition for similar at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/mul.jl:79 declares type variable N but does not use it.
WARNING: method definition for similar at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/muladd.jl:41 declares type variable N but does not use it.
WARNING: method definition for similar at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/lmul.jl:22 declares type variable N but does not use it.
WARNING: method definition for similar at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/lmul.jl:22 declares type variable N but does not use it.
WARNING: method definition for materialize! at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/triangular.jl:186 declares type variable UPLO but does not use it.
WARNING: method definition for materialize! at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/triangular.jl:192 declares type variable UPLO but does not use it.
WARNING: method definition for MemoryLayout at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/factorizations.jl:97 declares type variable T but does not use it.
WARNING: method definition for MemoryLayout at /juliateam/.julia/packages/ArrayLayouts/9PQtL/src/factorizations.jl:99 declares type variable T but does not use it.
WARNING: could not import Base._maybetail into BlockArrays
WARNING: method definition for _BlockArray at /juliateam/.julia/packages/BlockArrays/cOFZs/src/blockarray.jl:74 declares type variable T but does not use it.
WARNING: method definition for _BlockArray at /juliateam/.julia/packages/BlockArrays/cOFZs/src/blockarray.jl:76 declares type variable T but does not use it.
WARNING: method definition for reshape at /juliateam/.julia/packages/BlockArrays/cOFZs/src/pseudo_blockarray.jl:300 declares type variable N but does not use it.
WARNING: method definition for to_indices at /juliateam/.julia/packages/BlockArrays/cOFZs/src/views.jl:26 declares type variable R but does not use it.
WARNING: method definition for materialize! at /juliateam/.julia/packages/BlockArrays/cOFZs/src/blockbroadcast.jl:146 declares type variable Style but does not use it.
WARNING: method definition for similar at /juliateam/.julia/packages/BlockArrays/cOFZs/src/blocklinalg.jl:47 declares type variable N but does not use it.
WARNING: method definition for SymTensorValue at /juliateam/.julia/packages/Gridap/EZQEK/src/TensorValues/SymTensorValueTypes.jl:23 declares type variable T but does not use it.
WARNING: method definition for SymFourthOrderTensorValue at /juliateam/.julia/packages/Gridap/EZQEK/src/TensorValues/SymFourthOrderTensorValueTypes.jl:23 declares type variable T but does not use it.
WARNING: method definition for isless at /juliateam/.julia/packages/Gridap/EZQEK/src/TensorValues/Operations.jl:32 declares type variable T but does not use it.
WARNING: method definition for isless at /juliateam/.julia/packages/Gridap/EZQEK/src/TensorValues/Operations.jl:32 declares type variable D but does not use it.
WARNING: method definition for get_polytopes at /juliateam/.julia/packages/Gridap/EZQEK/src/Geometry/GridTopologyMocks.jl:81 declares type variable d but does not use it.
WARNING: method definition for _compute_hess_structure at /juliateam/.julia/packages/PDENLPModels/pW0Iv/src/hessian_struct_nnzh_functions.jl:70 declares type variable T but does not use it.
WARNING: method definition for _compute_hess_structure at /juliateam/.julia/packages/PDENLPModels/pW0Iv/src/hessian_struct_nnzh_functions.jl:74 declares type variable T but does not use it.
```

Definition of the domain and discretization

```
n = 20
domain = (-1, 1, -1, 1)
partition = (n, n)
model = CartesianDiscreteModel(domain, partition)
```

`CartesianDiscreteModel()`

Definition of the FE-spaces

```
reffe = ReferenceFE(lagrangian, Float64, 2)
Xpde = TestFESpace(model, reffe; conformity = :H1, dirichlet_tags = "boundary")
y0(x) = 0.0
Ypde = TrialFESpace(Xpde, y0)
reffe_con = ReferenceFE(lagrangian, Float64, 1)
Xcon = TestFESpace(model, reffe_con; conformity = :H1)
Ycon = TrialFESpace(Xcon)
Y = MultiFieldFESpace([Ypde, Ycon])
```

`MultiFieldFESpace()`

Integration machinery

```
trian = Triangulation(model)
degree = 1
dΩ = Measure(trian, degree)
```

`Measure()`

Objective function

```
yd(x) = -x[1]^2
α = 1e-2
function f(y, u)
∫(0.5 * (yd - y) * (yd - y) + 0.5 * α * u * u) * dΩ
end
```

`f (generic function with 1 method)`

Definition of the constraint operator

```
ω = π - 1 / 8
h(x) = -sin(ω * x[1]) * sin(ω * x[2])
function res(y, u, v)
∫(∇(v) ⊙ ∇(y) - v * u - v * h) * dΩ
end
op = FEOperator(res, Y, Xpde)
```

`FEOperatorFromWeakForm()`

Definition of the initial guess

```
npde = Gridap.FESpaces.num_free_dofs(Ypde)
ncon = Gridap.FESpaces.num_free_dofs(Ycon)
x0 = zeros(npde + ncon);
```

Overall, we built a GridapPDENLPModel, which implements the NLPModels.jl API.

```
nlp = GridapPDENLPModel(x0, f, trian, Ypde, Ycon, Xpde, Xcon, op, name = "Control elastic membrane")
(nlp.meta.nvar, nlp.meta.ncon)
```

`(1962, 1521)`

## Find a Feasible Point

Before solving the previously defined model, we will first improve our initial guess. We use `FeasibilityResidual`

from NLPModelsModifiers.jl to convert the NLPModel as an NLSModel. Then, using `trunk`

, a solver for least-squares problems implemented in JSOSolvers.jl, we find An improved guess which is close to being feasible for our large-scale problem. By default, a JSO-compliant solver such as `trunk`

(the same applies to `fps_solve`

) uses by default `nlp.meta.x0`

as an initial guess.

```
using JSOSolvers, NLPModelsModifiers
nls = FeasibilityResidual(nlp)
stats_trunk = trunk(nls)
```

`"Execution stats: first-order stationary"`

We check the solution from the stats returned by `trunk`

:

`norm(cons(nlp, stats_trunk.solution))`

`1.6058259477971256e-5`

We will use the solution found to initialize our solvers.

## Solve the Problem

Finally, we are ready to solve the PDE-constrained optimization problem with a targeted tolerance of `1e-5`

. In the following, we will use both Ipopt and DCI on our problem.

```
using NLPModelsIpopt
stats_ipopt = ipopt(nlp, x0 = stats_trunk.solution, tol = 1e-5, print_level = 0)
```

`"Execution stats: first-order stationary"`

The problem was successfully solved, and we can extract the function evaluations from the stats.

`nlp.counters`

```
Counters:
obj: ██████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 9 grad: ███████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 10 cons: ████████████⋅⋅⋅⋅⋅⋅⋅⋅ 18
cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons_nln: ████████████⋅⋅⋅⋅⋅⋅⋅⋅ 18 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ███████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 10 jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jac_nln: ███████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 10 jprod: ███████████████⋅⋅⋅⋅⋅ 22 jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod_nln: ███████████████⋅⋅⋅⋅⋅ 22 jtprod: ████████████████████ 31 jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jtprod_nln: ████████████████████ 31 hess: ██████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 8 hprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
```

Reinitialize the counters before re-solving.

`reset!(nlp);`

`NullLogger`

avoids printing iteration information.

```
using FletcherPenaltySolver, Logging
stats_fps_solve = with_logger(NullLogger()) do
fps_solve(nlp, stats_trunk.solution, atol = 1e-5, rtol = 1e-5)
end
```

`"Execution stats: first-order stationary"`

The problem was successfully solved, and we can extract the function evaluations from the stats.

`nlp.counters`

```
Counters:
obj: ██████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 3 grad: ████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 4 cons: ████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 4
cons_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 cons_nln: ████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 4 jcon: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jgrad: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jac: ████████████⋅⋅⋅⋅⋅⋅⋅⋅ 6 jac_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jac_nln: ████████████⋅⋅⋅⋅⋅⋅⋅⋅ 6 jprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jprod_nln: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jtprod: ██████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 5 jtprod_lin: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
jtprod_nln: ██████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 5 hess: ████████⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 4 hprod: ████████████████████ 10
jhess: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 jhprod: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
```

We now compare the two solvers with respect to the time spent,

`stats_ipopt.elapsed_time, stats_fps_solve.elapsed_time`

`(8.66, 105.24244213104248)`

and also check objective value, feasibility and dual feasibility of `ipopt`

and `fps_solve`

.

```
(stats_ipopt.objective, stats_ipopt.primal_feas, stats_ipopt.dual_feas),
(stats_fps_solve.objective, stats_fps_solve.primal_feas, stats_fps_solve.dual_feas)
```

`((0.005425026428348349, 2.2204460492503135e-18, 4.2724309076258577e-7), (0.005425025968573663, 2.2204460492503135e-18, 2.968652346957424e-7))`

Overall `FletcherPenaltySolver`

is doing great for solving large-scale optimization problems!

*This page was generated using Literate.jl.*