Finite Element Discretization of Elliptic Optimal Control

Overview¶

The previous lectures developed the continuous theory for linear-quadratic elliptic optimal control, derived the adjoint-based gradient formula, and established the first-order optimality system in both unconstrained and box-constrained settings.

Before introducing finite elements, we place the elliptic control problem in a general analytical framework:

1. a general variational formulation for elliptic state equations;

2. well-posedness of the control-to-state map;

3. passage from infinite-dimensional analysis to finite-dimensional approximation.

A General Elliptic Boundary Value Problem¶

Let $\Omega\subset\mathbb R^d$ be a bounded Lipschitz domain with boundary $\Gamma$. We consider the second-order linear elliptic operator

with coefficients satisfying standard assumptions:

• $\mathbf A\in L^\infty(\Omega;\mathbb R^{d\times d})$ is uniformly elliptic, i.e. there exist constants $0<\alpha\le M<\infty$ such that

  $$\alpha |\xi|^2 \le \xi^T \mathbf A(x)\xi \le M |\xi|^2 \qquad \forall \xi\in\mathbb R^d,\ \text{for q.o. }x\in\Omega;$$

  (2)

• $\mathbf b\in L^\infty(\Omega)^d$ and

  $$\|\mathbf b\|_{L^\infty(\Omega)}\le \beta;$$

  (3)

• $c\in L^\infty(\Omega)$ and

  $$\|c\|_{L^\infty(\Omega)}\le \gamma;$$

  (4)

• $\operatorname{div}\mathbf b\in L^\infty(\Omega)$;

• the lower-order part is not too negative, in the sense that there exists $\alpha_0\ge 0$ such that

  $$\frac12 \operatorname{div}\mathbf b(x) + c(x)\ge -\alpha_0 \qquad \text{for q.o. } x\in\Omega,$$

  (5)

  and $\alpha_0$ is small enough compared with the ellipticity constant $\alpha$, namely

  $$\alpha_0 < \frac{\alpha}{C_P^2},$$

  (6)

  where $C_P$ is the Poincare constant of $\Omega$.

We impose homogeneous Dirichlet boundary conditions

We work in the standard Sobolev space

The weak formulation of the state equation is: find $y\in V$ such that

where the bilinear form is

Special case:

The model studied in the previous lectures is the simplest representative of this class.

Trace Theorem and Boundary Data¶

Since $\Omega$ is Lipschitz, the trace operator

is continuous.

Consequences:

• boundary conditions are meaningful in the Sobolev setting;

• boundary control problems can also be formulated rigorously, even though in this lecture we focus on distributed controls.

For $s=1$:

The trace operator admits a continuous right inverse: there exists a bounded lifting operator

such that

This yields the standard lifting of nonhomogeneous Dirichlet data. For homogeneous Dirichlet data, the variational space is $H_0^1(\Omega)$.

Continuity and Coercivity of the Bilinear Form¶

Continuity¶

Since $\mathbf A$, $\mathbf b$, and $c$ are bounded, each term is continuous.

Transport term:

Zero-order term:

For $y,v\in H_0^1(\Omega)$, Poincare’s inequality yields

Hence

On $H_0^1(\Omega)$, the norm $\|w\|_{H^1(\Omega)}$ is equivalent to $\|\nabla w\|_{L^2(\Omega)}$. Therefore

with the explicit choice

Coercivity¶

Let $v\in V=H_0^1(\Omega)$. Then

Uniform ellipticity gives

Integrating the first-order term by parts:

The boundary term vanishes since $v=0$ on $\Gamma$. Hence

Under the assumption

Poincare’s inequality gives

Hence

If

the right-hand side is strictly positive. Since on $H_0^1(\Omega)$ the seminorm $\|\nabla v\|_{L^2(\Omega)}$ is equivalent to the full norm,

for some $\alpha_0>0$.

Nonhomogeneous Dirichlet and Neumann Conditions¶

Consider

with

If $g_D\neq 0$, one works in an affine space with prescribed trace, or reduces to the homogeneous case by lifting. If Neumann data are present, the boundary term appears on the right-hand side of the weak formulation.

For pure Dirichlet data on $\Gamma$, let

Choose a lifting

Set

Then $z\in H_0^1(\Omega)$ and

Define

By continuity of $a(\cdot,\cdot)$,

Hence

Lax-Milgram gives

and therefore

Well-Posedness of the State Equation¶

By Lax-Milgram, for every $u\in L^2(\Omega)$ and every $F\in V'$, there exists a unique state $y\in V$ such that

The solution depends continuously on the data. This defines the control-to-state operator

Cea and Bramble-Hilbert Lemmas¶

Let $V_h\subset V$ be a conforming finite-dimensional space and let $y_h\in V_h$ be the Galerkin approximation of the state:

Cea’s lemma gives the quasi-optimality estimate

The finite element error is controlled, up to the factor $C/\alpha_0$, by the best approximation error in the chosen discrete space.

The approximation term

is estimated by the Bramble-Hilbert lemma. For standard finite element spaces of polynomial degree $k$ on a shape-regular mesh, if the exact solution is sufficiently smooth, for example $y\in H^{k+1}(\Omega)$, then the interpolation error satisfies

Combining this with Cea’s lemma gives

• the weak formulation is posed in Sobolev spaces;

• Galerkin approximations inherit stability from coercivity;

• approximation theory converts best approximation into concrete mesh-dependent error bounds.

Linear-Quadratic Optimal Control Problem¶

Let

Consider

subject to

The state equation defines the control-to-state map

The reduced functional is

The reduced problem is

Continuous Optimality System¶

State equation:

Adjoint equation:

Gradient condition in the unconstrained case:

Equivalently, in dual notation,

If control constraints are present, this last equation is replaced by the corresponding variational inequality or projection formula.

Optimize-Then-Discretize and Discretize-Then-Optimize¶

Two routes are available.

Optimize-Then-Discretize¶

1. derive the continuous state equation;

2. derive the continuous adjoint equation;

3. derive the continuous optimality system;

4. discretize the resulting system.

Discretize-Then-Optimize¶

1. discretize the state equation;

2. define the discrete cost functional;

3. formulate the finite-dimensional optimization problem;

4. derive its KKT conditions.

In linear-quadratic elliptic settings these routes often lead to the same algebraic system. In more general settings they may differ.

Choice of Discrete Spaces¶

Let $\mathcal T_h$ be a finite element mesh of $\Omega$. We choose a finite-dimensional subspace

for the state and adjoint variables.

Typical choice:

• $V_h$ = continuous piecewise polinomial finite elements.

For the control, there are several possible choices. The two most common are:

1. choose a discrete control space

   $$Q_h\subset Q=L^2(\Omega),$$

   (62)

   for example piecewise polynomials or continuous piecewise polynomials;

2. do not discretize the control explicitly at first, and only discretize the state and adjoint.

In this lecture we focus first on the standard fully discrete setting

Why are these choices important?

• the state equation must be well posed in the discrete space;

• the control space determines the number and meaning of the optimization variables;

• the discrete gradient and the discrete KKT system depend on the pairing between $V_h$ and $Q_h$.

The passage from the continuous problem to a numerical method is not only the replacement of functions by vectors: it also requires a choice of finite-dimensional spaces.

Discrete State Equation¶

Given $u_h\in Q_h$, the finite element state equation is: find $y_h\in V_h$ such that

This is the Galerkin approximation of the continuous weak problem.

Choose bases

Expand

Then the discrete state equation becomes the linear system

where

Interpretation of the matrices:

• $A$ is the stiffness matrix of the PDE;

• $B$ transfers the control unknowns into the state equation;

• $f$ is the load vector.

If $Q_h=V_h$ and the same basis is used, then $B$ is simply the mass matrix. If $Q_h$ is piecewise constant, then $B$ is a rectangular coupling matrix.

The important structural point is: for each discrete control vector $u$, there is a unique discrete state vector $y$, provided the stiffness matrix $A$ is invertible.

The discrete control-to-state map is

Discrete Cost Functional¶

Once the state has been discretized, the cost functional becomes

In matrix form, this reads

where

The fully discrete constrained problem is

The conceptual path of Lecture 1 reappears in finite element form:

1. choose a control vector $u$;

2. solve the discrete state equation to get $y=S_h(u)$;

3. evaluate the discrete objective;

4. apply finite-dimensional optimization ideas.

Reduced Discrete Problem¶

Since the discrete state is uniquely determined by the discrete control, we may eliminate $y$ and define the reduced discrete functional

Then the discrete optimal control problem becomes

or, with bounds,

where $U_{\mathrm{ad}}^h\subset \mathbb R^{n_u}$ is the discrete admissible set.

Remarks:

• the unknown is now a vector of control coefficients;

• each evaluation of $f_h(u)$ still requires the solution of one PDE-like linear system;

• all numerical optimization methods from Lecture 3 now apply directly to $f_h$.

From the algorithmic point of view, the reduced finite element problem is a standard optimization problem in which function and gradient evaluations require PDE solves.

Discrete Adjoint Equation¶

To compute the gradient of $f_h$, we follow the continuous theory. Given $u_h\in Q_h$, first solve the discrete state equation for $y_h$, and then define the discrete adjoint $p_h\in V_h$ by

In matrix form, this becomes

For symmetric elliptic operators such as Poisson with standard Dirichlet conditions, the stiffness matrix is symmetric, so $A^T=A$. The transpose is kept in the general derivation: the adjoint equation is governed by the adjoint operator.

Once $p_h$ has been computed, the reduced gradient is

Finite element counterpart of the continuous identity:

when the control space is identified with $L^2(\Omega)$ through the Riesz map.

Computational pattern:

1. state solve;

2. adjoint solve;

3. gradient assembly.

The adjoint principle is unchanged by discretization.

Discrete KKT System¶

Instead of reducing the problem to the control alone, we may keep state, control, and adjoint unknowns together. In the unconstrained case, the discrete first-order system is

Reordering the equations as before, this becomes the block system

This is the discrete saddle-point system associated with the optimality conditions.

Structure:

• the $(1,1)$ block comes from state tracking;

• the $(2,2)$ block comes from Tikhonov regularization;

• the off-diagonal blocks encode state and adjoint couplings;

• the zero block in the lower-right corner is the saddle-point signature.

Discretization preserves the continuous KKT structure and yields a sparse block linear system.

This algebraic form motivates:

• sparse direct solvers;

• Krylov methods;

• block preconditioners;

• all-at-once implementations.

Box Constraints in the Discrete Problem¶

Suppose now that the control is constrained componentwise:

Then the reduced discrete problem is

where $U_{\mathrm{ad}}^h$ is a box in $\mathbb R^{n_u}$.

The first-order condition becomes the finite-dimensional variational inequality

Using the discrete gradient formula,

If the control basis is chosen so that the discrete control inner product is diagonal or easily invertible, this leads to the discrete projection formula

In the simplest lumped-mass or nodal setting, this reduces componentwise to

The discrete constrained problem has the same geometry as the finite-dimensional KKT picture from Lecture 2 and the continuous projection picture from Lecture 5.

A Minimal 1D Finite Element Thought Experiment¶

To make the previous formulas concrete, imagine $\Omega=(0,1)$ and let $V_h$ be the space of continuous piecewise linear functions on a uniform mesh.

Then:

• the stiffness matrix $A$ is tridiagonal;

• the mass matrix $M_y$ is sparse and symmetric positive definite;

• the discrete state solve means solving one sparse linear system;

• the discrete adjoint solve has the same matrix and a different right-hand side.

One reduced-gradient iteration has the form:

1. choose the current control vector $u_k$;

2. solve

   $$A y_k = B u_k + f;$$

   (89)

3. solve

   $$A^T p_k = M_y(y_k-y_d);$$

   (90)

4. compute

   $$g_k=\alpha M_u u_k+B^T p_k;$$

   (91)

5. update the control with a finite-dimensional optimization step.

This is the finite element version of the state-adjoint-gradient loop of the continuous reduced formulation.

References¶

• Fredi Troeltzsch, Optimal Control of Partial Differential Equations, Chapters 3 and 4.

• Juan Carlos De Los Reyes, Numerical PDE-Constrained Optimization.

• Alfio Borzi and Volker Schulz, Computational Optimization of Systems Governed by Partial Differential Equations.