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.