Control Problem with Sum of Squares Cost

[aarraf]

Hi everyone,

I want to solve the linear system resulting form an (linear elliptic) PDE constraint optimization problem using the "All-at-once" approach. The optimality system, see [1], reads as follows (y =state, p=adjoint and u=control):

I'm using no control constraints. However, my cost function includes a Sum of Squares Cost over some points, say $X_m$. This results in an Adjoint PDE with point sources located at $X_m$ as rhs. I know that point sources can be implemented in the dual space. However I'm not sure how to formulate the problem, that it can be solved with a LinearVariationalSolver, more specifically how project the current state y on the points $X_m$.

My first attempt was to use the interpolation matrix $C$ (see picture) as decried in the discussion: #4045. However, I was wondering if there is a more elegant way (with VertexOnlyMesh) which directly uses UFL Forms to do that.

Thanks in advance.

[1] Tröltzsch, "Optimal control of partial differential equations", 2010

[dham]

I think this demo has what you need: https://www.firedrakeproject.org/demos/full_waveform_inversion.py.html . In particular, look at how the point sources for the wave equation are implemented.

[aarraf]

Thanks for your fast reply.

As far as I understood, the approach in the demo works only for iterative solution of the optimality system, since for the adjoint PDE, the state is a known function. However, I am interested how to directly solve the discretized optimality system [1]:

by a single linear solve, e.g. with GMRES.

I hope this makes my problem more clear.

Thanks again for your help.

[1] Rees, Preconditioning Iterative Methods for PDE Constrained Optimization, - 2010

[dham]

I think we may be talking somewhat a cross purposes here. Above I thought you were asking how to use VertexOnlyMesh to force a PDE with point forces, which is covered in the demo I linked to. That demo does use the adjoint to go on and solve an inverse problem, but that's rather independent of the point forcing, which is there in the original forward problem.

For the more general question of how to build and solve a KKT system in Firedrake, the thing to note is that Firedrake represents a block matrix system by declaring the solution variables to be in a mixed space. So you make a MixedFunctionSpace from whichever spaces you chose for u, y, and p and then formulate the KKT system there just as if you were forming any other system in multiple variables (e.g. Navier Stokes).

An even more automated and elegant system is also possible: use Lagrange multipliers to create a UFL expression for an unconstrained quantitity of interest to optimise, once again using a mixed system to represent the solution variable(s) and the lagrange multipliers. You then use derivative to take the derivative of the quantity of interest with respect to the combined solution and lagrange multiplier variable. This gets you the residual of the KKT system and you can just use that in NonlinearVariationalProblem to solve the whole system (this approach works even if the original system is linear: in that case you can, if you wish, set the snes_type to ksponly in order to ensure that only one linear solve happens).

[aarraf]

Thanks for your reply and sorry for the confusion. I will add some more background to hopefully make it more clear.

I want to solve an Optimal Control Problem with the steady Advection-Diffusion PDE as constraint and the Sum of Squares Cost functional
$J(u) = \int_\Omega \sum_i^N (y^i - y_m^i)^2 \delta(\mathbf{x}-\mathbf{x}_m^i)$
with state $y^i$ and measurement $y^i_m$ at the measurement point $\mathbf{x}_m^i$ .

This results in an adjoint PDE of the form
$-\epsilon \nabla^2 p - \mathbf{w} \cdot \nabla p = \sum_i^N (y^i - y_m^i) \delta(\mathbf{x}-\mathbf{x}_m^i) $
where the RHS includes the residuals between the states $y^i$ and measurement $y_m^i$ evaluated at measurement point $\mathbf{x}_m^i$. Thus a sum of point sources. Since the resulting KKT system is steady and linear, I wanted to solve it directly with a single linear solver.

Regarding the implementation I can express the KKT system as mixed ufl form (named $F$) expected the RHS of the adjoint PDE. The measurement $y_m^i$ can be assembled (to the cofunction y_m) beforehand and simply added to the linear solver
problem = LinearVariationalProblem(lhs(F), rhs(F) - y_m, w, bcs=bcs)

My problem comes from the term $\sum_i^N y^i \delta(\mathbf{x}-\mathbf{x}^i) $, since it includes to evaluation of the state y at the measurement points. I tried
interpolate(y, P0DG) * q * dx
where y is the state trial function, q is the adjoint test function and P0DG point space. This raises the error
AttributeError: 'MixedMap' object has no attribute '_offset'. Did you mean: 'offset'?

So my question is how to formulate the point evaluation of the state in a way, that leads to a KKT system which ca be solved with a single call to a linear solver?

A possible solutions seems to be to use the point evaluation matrix $C$ (see e.g., #4045) which would result in this adjoint PDE
$-\epsilon \nabla^2 p - \mathbf{w} \cdot \nabla p = C^T( C y - \mathbf{y}_m ) $
where point evaluation matrix $C$ has dimension (num_measurements x num_dofs) and $\mathbf{y}_m$ is a vector of dimension num_measurements. However this requires more interaction with PETSc matrices and vectors which I'm not very familiar with.

Sorry for the amount of text and thanks in advance.

[dham]

I think the derivation of your adjoint source function isn't quite right. I note that you are writing all of your maths in strong form, which is I think the issue. The steps in your cost functional are to interpolate the solution to point space, then compute a point space integral. The adjoint to this is to compute a point space cofunction and then adjoint interpolate this cofunction into the dual to the solution space (this being the space where the adjoint PDE is posed). The code for doing this sort of point forcing in the demo is that below. using Cofunction.interpolate is the adjoint interpolation which is equivalent to $C^T$ in your maths above. Obviously this isn't the same integrand as the one you need, but the process is the same.

As noted, this is the hard way of doing this. I think you can probably get Firedrake to do this for you by writing down the full unconstrained cost function and taking its derivative.

To compute the cofunction q_s, we first construct the source mesh over the source locations, for the source number source_number:

source_mesh = VertexOnlyMesh(mesh, [source_locations[source_number]])

Next, we define a function space accordingly:

V_s = FunctionSpace(source_mesh, "DG", 0)

The point source value is coded as:

d_s = Function(V_s)
d_s.interpolate(1.0)

We then interpolate a cofunction in
onto to then have :

source_cofunction = assemble(d_s * TestFunction(V_s) * dx)
q_s = Cofunction(V.dual()).interpolate(source_cofunction)

[aarraf]

Thank you for the hint regarding the strong formulation. Your right, this caused a lot of confusion since I thought I had to somehow assemble the RHS to a vector. Sorry for asking so many question but I trying to understand how it works. I will give it another try. So...

Multiplying the adjoint RHS with the test function $q$ leads to the weak form:
$\int_\Omega \sum_i^N (y^i - y_m^i) \delta(\mathbf{x}-\mathbf{x}_m^i) q = \sum_i^N (y^i(x_m^i) - y_m^i) q(x_m^i) $

If I understood correctly, I could than define the cofunction $q(x_m^i)$ (corresponding to source_number=i) with the code you provided. But how do I than define the residual $(y^i(x_m^i) - y_m^i)$, since at the moment where the forms are defined, this depends on the unknown state y (trial function)? <\del>

Update: I realized that I did not understood it correctly. I need to rethink how dual spaces work.

Thanks for your help and patience.