Solve Large-Scale Problem with DCISolver

In this tutorial we use dci 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

Definition of the domain and discretization

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

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])

Integration machinery

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

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

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)

Definition of the initial guess

npde = Gridap.FESpaces.num_free_dofs(Ypde)
ncon = Gridap.FESpaces.num_free_dofs(Ycon)
x0 = zeros(npde + ncon);
nothing #hide

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)

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 dci) uses by default nlp.meta.x0 as an initial guess.

using JSOSolvers, NLPModelsModifiers

nls = FeasibilityResidual(nlp)
stats_trunk = trunk(nls)

We check the solution from the stats returned by trunk:

norm(cons(nlp, stats_trunk.solution))

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)

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

stats_ipopt.counters

Reinitialize the counters before re-solving.

reset!(nlp);
nothing #hide

NullLogger avoids printing iteration information.

using DCISolver, Logging

stats_dci = with_logger(NullLogger()) do
  dci(nlp, stats_trunk.solution, atol = 1e-5, rtol = 0.0)
end

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

stats_dci.counters

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

stats_ipopt.elapsed_time, stats_dci.elapsed_time

and also check objective value, feasibility and dual feasibility of ipopt and dci.

(stats_ipopt.objective, stats_ipopt.primal_feas, stats_ipopt.dual_feas),
(stats_dci.objective, stats_dci.primal_feas, stats_dci.dual_feas)

Overall DCISolver is doing great for solving large-scale optimization problems!

This page was generated using Literate.jl.