Numerical solution of elliptic distributed optimal control problems with boundary value tracking

This paper presents a finite element discretization for elliptic distributed optimal control problems with boundary value tracking, reformulating the problem as a state-based variational form to derive optimal error estimates and fast solvers, which are validated through numerical experiments.

The Gist

Imagine you are a master chef trying to bake the perfect cake (the target). You have a recipe that tells you exactly how the cake should look and taste on the outside (the boundary). However, you can't just magically shape the cake; you have to control the oven's temperature distribution (the control) to get the result you want.

This paper is about solving a very specific, tricky version of that problem using math and computers. Here is the breakdown in everyday language:

1. The Problem: The "Too Hot, Too Cold" Oven

In the real world, if you want a cake to have a specific pattern on the crust, you need to adjust the heat inside the oven perfectly.

• The Goal: Make the outside of the cake (the boundary) look exactly like a picture you have in your head.
• The Catch: You can only control the heat inside the oven, not the outside directly. Also, the oven has a rule: no heat can escape through the walls (this is the Neumann boundary condition).
• The Difficulty: If you try to force the outside to look exactly like your picture, the math says you might need infinite energy or create a chaotic, unstable oven. To fix this, the authors add a "penalty" (a regularization parameter). Think of it as a rule that says, "You can change the heat, but don't go crazy with it." This balances getting a good-looking cake with keeping the oven stable.

2. The Trick: Turning the Problem Inside Out

Usually, to solve this, you'd have to guess the heat, bake the cake, check the outside, guess again, and repeat. That's slow.

The authors found a clever shortcut. They realized that because of the physics of the oven, there is a direct, one-to-one link between the "heat settings" and the "final cake shape."

• The Analogy: Instead of asking, "What heat setting gives me this cake?", they asked, "If I want this specific cake shape, what does the heat setting have to be?"
• By flipping the problem, they turned a complex two-step puzzle into a single-step math problem. This made it much easier to solve on a computer.

3. The Method: Building a Digital Lego Model

To solve this on a computer, you can't treat the oven as a smooth, continuous block. You have to break it down into tiny chunks, like Legos.

• Tensor-Product Mesh: Imagine a giant 3D grid of cubes. The authors used a very specific, neat way of stacking these cubes (like a perfect grid of bricks) rather than a messy pile.
• Why it matters: This neat structure allows them to use "fast solvers." Think of it like this: If you have a messy pile of 1,000,000 Lego bricks, finding a specific one takes forever. If they are stacked in a perfect 100x100x100 grid, you can find any brick instantly using a simple coordinate system.

4. The Speed: The "Magic Elevator" (Fast Solvers)

When you break the problem into millions of tiny Lego pieces, you end up with a massive system of equations. Solving this normally would take a supercomputer days.

• The Innovation: The authors developed a "magic elevator" (called a Schur Complement method with Preconditioning).
• How it works: Instead of trying to solve the whole 3D building at once, the math allows them to focus only on the "skin" of the building (the boundary) and use a shortcut to figure out the inside.
• The Result: No matter how many Lego bricks they use (whether 27 or 16 million), the computer solves the problem in roughly the same amount of time (about 4 to 20 steps). It's like having an elevator that takes the same amount of time to go from the 1st floor to the 10th floor as it does from the 1st to the 100th.

5. The Experiments: Testing the Recipe

The authors tested their method with three different "cakes" (targets):

1. A Smooth Cake: A gentle, wavy pattern. The math worked perfectly, and the computer got very close to the target quickly.
2. A Bumpy Cake: A pattern that was a bit rougher. The computer was still fast, but the accuracy dropped slightly, just as the math predicted.
3. A Blocky Cake: A pattern that was a sharp square (discontinuous). This is the hardest case. The computer still solved it quickly, though the accuracy was lower, which was expected.

The Bottom Line

This paper is about finding a fast, reliable, and efficient way to control a system (like an oven, a heat shield, or a medical imaging device) to match a desired shape on the surface, without wasting computing power.

They proved that by using a specific type of grid (Lego bricks) and a clever mathematical shortcut (the magic elevator), you can solve these complex control problems almost instantly, even as the details get finer and finer. This is huge for fields like medical imaging (seeing inside the body) and heat transfer (managing temperature in engines), where speed and accuracy are everything.

Technical Summary

Here is a detailed technical summary of the paper "Numerical solution of elliptic distributed optimal control problems with boundary value tracking."

1. Problem Statement

The paper addresses a boundary value tracking optimal control problem (OCP) constrained by an elliptic partial differential equation (PDE) with Neumann boundary conditions.

• Objective: Minimize a cost functional $J(y_\varrho, u_\varrho)$ that balances the tracking error of the state $y_\varrho$ against a target $y$ on the boundary $\Gamma$, and a regularization term for the control $u_\varrho$:
  $$J(y_\varrho, u_\varrho) = \frac{1}{2} \|y_\varrho - y\|^2_{L^2(\Gamma)} + \frac{1}{2} \varrho \|u_\varrho\|^2_{U}$$
• Constraints: The state is governed by the Neumann boundary value problem:
  $$-\Delta y_\varrho + y_\varrho = u_\varrho \quad \text{in } \Omega, \quad \partial_n y_\varrho = 0 \quad \text{on } \Gamma$$
• Key Distinction: Unlike previous works using $L^2(\Omega)$ regularization, this paper investigates energy regularization where the control space is chosen as $U = \tilde{H}^{-1}(\Omega) := [H^1(\Omega)]^*$. This choice ensures the state-to-control map is an isomorphism within $H^1(\Omega)$, avoiding the need for higher-order basis functions required by $L^2$ regularization.

2. Methodology

A. State-Based Variational Reformulation

The authors reformulate the OCP into a state-based variational problem. By utilizing duality arguments, they establish that the norm of the control in the dual space is equivalent to the $H^1(\Omega)$ norm of the state. This allows the elimination of the control variable $u_\varrho$, reducing the problem to minimizing a functional dependent solely on the state $y_\varrho$:
$$\tilde{J}(y_\varrho) = \frac{1}{2} \|y_\varrho - y\|^2_{L^2(\Gamma)} + \frac{1}{2} \varrho \|y_\varrho\|^2_{H^1(\Omega)}$$
The minimizer satisfies a gradient equation (5) involving the $L^2(\Gamma)$ and $H^1(\Omega)$ inner products.

B. Finite Element Discretization

The paper focuses on a conforming tensor-product finite element (FE) discretization on a unit hypercube domain $\Omega = (0,1)^d$.

• Basis Functions: They construct a specific FE space $Y_h \subset Y$ using modified piecewise linear basis functions. Crucially, these basis functions are designed to have vanishing normal derivatives ($\partial_n \phi = 0$) on the boundary, satisfying the Neumann condition naturally without Lagrange multipliers.
• Discrete System: The discretization leads to a linear system of the form:
  $$(\hat{M}_h + \varrho(\tilde{K}_h + \mathring{K}_h))y_h = y_h$$
  where the matrices represent boundary mass, boundary stiffness, and interior stiffness contributions.

C. Fast Solvers

To solve the resulting linear system efficiently, the authors propose a Schur Complement (SC) approach:

1. Elimination: Interior degrees of freedom are eliminated, reducing the system to a boundary Schur complement system $S_{BB}y_B = y_B$.
2. Spectral Equivalence: They prove that for $\varrho \leq h$, the Schur complement $S_{BB}$ is spectrally equivalent to the boundary mass matrix $M_{BB}$.
3. Iterative Solvers:
   • The SC system is solved using the Conjugate Gradient (CG) method.
   • Preconditioning: A lumped mass matrix is used as a diagonal preconditioner for the PCG method, ensuring optimal scaling.
   • Inner Solver: The action of the inverse of the interior stiffness matrix (required for the Schur complement) is approximated using an Algebraic MultiGrid (AMG) preconditioned CG method.

3. Key Contributions

1. Theoretical Error Estimates: The paper derives rigorous error estimates for the discretized state $y_{\varrho h}$ relative to the target $y$. The convergence rates depend on the regularity of the target:
   • If $y \in H^{1/2}(\Gamma)$, error is $O(h^{1/2})$ (with $\varrho = h$).
   • If $y \in H^1(\Gamma)$, error is $O(h)$ (with $\varrho = h$).
   • If $y \in H^2(\Gamma)$ (smooth target), error is $O(h^2)$ (with $\varrho = h^2$).
2. Optimal Solver Construction: The derivation of a fast solver based on the Schur complement with a lumped mass preconditioner. Theoretical analysis shows that the condition number is bounded independently of the mesh size $h$ and the regularization parameter $\varrho$ (provided $\varrho \leq h$).
3. Handling Neumann Conditions: The construction of a conforming FE space with basis functions explicitly satisfying homogeneous Neumann conditions, simplifying the implementation compared to methods requiring Lagrange multipliers.

4. Numerical Results

The authors validate their theory with numerical experiments in 3D ($\Omega = (0,1)^3$) using three different target functions:

• Case 1: Smooth Target ($y \in H^2(\Gamma)$)

  • Target: $y = \cos(\pi x_1)\cos(\pi x_2)\cos(\pi x_3)$.
  • Parameter choice: $\varrho = h^2$.
  • Result: Achieved second-order convergence ($O(h^2)$) in the $L^2(\Gamma)$ error.
  • Solver Performance: The number of PCG iterations for the Schur complement system remained independent of the mesh level (approx. 4–29 iterations depending on preconditioner), confirming the robustness of the solver.

• Case 2: Less Regular Target ($y \in H^{3/2-\epsilon}(\Gamma)$)

  • Target: A quadratic polynomial not satisfying homogeneous Neumann conditions.
  • Parameter choice: $\varrho = h^{3/2}$.
  • Result: Observed convergence rate of $O(h^{1.5})$, matching theoretical predictions.

• Case 3: Discontinuous Target ($y \in L^2(\Gamma)$)

  • Target: A piecewise constant function (1 inside a sub-region, 0 outside).
  • Parameter choice: $\varrho = h$.
  • Result: Observed convergence rate of $O(h^{0.5})$, consistent with the lowest regularity case.

In all cases, the AMG-PCG and Schur Complement PCG solvers demonstrated mesh-independent iteration counts, proving the efficiency of the proposed fast solvers.

5. Significance and Outlook

• Significance: This work provides a complete framework for solving boundary tracking OCPs with energy regularization. It bridges the gap between abstract variational reformulation and practical, efficient numerical implementation using tensor-product meshes. The derivation of optimal error estimates and the construction of a solver with mesh-independent convergence are critical for large-scale 3D applications.
• Applications: The methodology is relevant for inverse problems such as inverse heat transfer and photoacoustic tomography, where boundary data is used to infer internal states or sources.
• Future Work: The authors note that for general domains (non-tensor-product), the homogeneous Neumann conditions may need to be enforced via Lagrange multipliers, which is identified as a topic for future research. Additionally, the current conforming discretization is restricted to tensor-product meshes.