---
jupytext:
  main_language: python
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.2
---

```{eval-rst}
.. include:: ../../../global.rst
```

(demo_poisson)=
# Distributed Control of a Poisson Problem

## Problem Formulation

In this demo we investigate the basics of cashocs for
optimal control problems. To do so, we investigate the "mother
problem" of PDE constrained optimization, i.e.,

$$
\begin{align}
    &\min\; J(y,u) = \frac{1}{2} \int_{\Omega} \left( y - y_d \right)^2
    \text{ d}x + \frac{\alpha}{2} \int_{\Omega} u^2 \text{ d}x \\
    &\text{ subject to } \qquad
    \begin{alignedat}[t]{2}
        -\Delta y &= u \quad &&\text{ in } \Omega,\\
        y &= 0 \quad &&\text{ on } \Gamma.
    \end{alignedat}
\end{align}
$$

(see, e.g., [Tröltzsch - Optimal Control of Partial Differential Equations](
https://doi.org/10.1090/gsm/112)
or [Hinze, Pinnau, Ulbrich, and Ulbrich - Optimization with PDE constraints](
https://doi.org/10.1007/978-1-4020-8839-1).

For this first example, we do not consider control constraints,
but search for an optimal control u in the entire space $L^2(\Omega)$,
for the sake of simplicitiy. For the domain under consideration, we use the unit
square $\Omega = (0, 1)^2$, since this is built into cashocs.

In the following, we will describe how to solve this problem
using cashocs. Moreover, we also detail alternative / equivalent FEniCS code which
could be used to define the problem instead.

## Implementation

The complete python code can be found in the file
{download}`demo_poisson.py
</../../demos/documented/optimal_control/poisson/demo_poisson.py>`
and the corresponding config can be found in
{download}`config.ini </../../demos/documented/optimal_control/poisson/config.ini>`.

### Initialization

We begin by importing FEniCS and cashocs. For the sake of
better readability we use a wildcard import for FEniCS

```python
from fenics import *

import cashocs
```

Afterwards, we can specify the so-called log level of cashocs. This is done in the
line

```python
cashocs.set_log_level(cashocs.LogLevel.INFO)
```

:::{hint}
There are a total of five levels of verbosity, given by
{py:class}`cashocs.LogLevel.DEBUG`, {py:class}`cashocs.LogLevel.INFO`,
{py:class}`cashocs.LogLevel.WARNING`, {py:class}`cashocs.LogLevel.ERROR`,
and {py:class}`cashocs.LogLevel.CRITICAL`. The default value is {python}`INFO`, which
would also be selected if the {py:func}`cashocs.set_log_level` method would not have
been called.
:::

Next, we have to load the config file which loads the user's input parameters into the
script. For a detailed documentation of the config files and the parameters within, we
refer to {ref}`config_optimal_control`. Note that the corresponding file is
{download}`config.ini </../../demos/documented/optimal_control/poisson/config.ini>`.
The config is then loaded via

```python
config = cashocs.load_config("config.ini")
```

Next up, we have to define the state equation. This mostly works with usual FEniCS
syntax. In cashocs, we can quickly generate meshes for squares and cubes, as well as
import meshes generated by Gmsh, for more complex geometries. In this example we take
a built-in unit square as example. This is generated via

```python
mesh, subdomains, boundaries, dx, ds, dS = cashocs.regular_mesh(25)
```

The input for regular_mesh determines the number of elements that
are placed along each axis of the square. Note, that the mesh could be
further manipulated with additional, optional arguments, and we
refer to {py:func}`regular_mesh <cashocs.regular_mesh>` for more infos. Note,
that the {python}`subdomains` object is a (empty) {py:class}`MeshFunction`, and that
{python}`boundaries` is a {py:class}`MeshFunction` that contains markers for the following
boundaries

- The left side of the square is marked by 1
- The right side is marked by 2
- The bottom is marked by 3
- The top is marked by 4,

as defined in {py:func}`regular_mesh <cashocs.regular_mesh>`.

With the geometry defined, we create a function space with the classical
FEniCS syntax

```python
V = FunctionSpace(mesh, "CG", 1)
```

which creates a function space of continuous, linear Lagrange elements.

### Definition of the state equation

To describe the state system in cashocs, we use (almost) standard
FEniCS syntax, and the differences will be highlighted in the
following. First, we define a {py:class}`fenics.Function` {python}`y` that models our
state variable $y$, and a {py:class}`fenics.Function` {python}`p` that models
the corresponding adjoint variable $p$ via

```python
y = Function(V)
p = Function(V)
```

Next up, we analogously define the control variable as {py:class}`fenics.Function`
{python}`u`

```python
u = Function(V)
```

This enables us to define the weak form of the state equation,
which is tested not with a {py:class}`fenics.TestFunction` but with the adjoint
variable {python}`p` via the classical FEniCS / UFL syntax

```python
e = inner(grad(y), grad(p)) * dx - u * p * dx
```

::::{note}
For the clasical definition of this weak form with FEniCS
one would write the following code

:::{code-block} python
y = TrialFunction(V)
p = TestFunction(V)
u = Function(V)
a = inner(grad(y), grad(p)) * dx
L = u * p * dx
:::

as this is a linear problem. However, to have greater flexibility
we have to treat the problems as being potentially nonlinear.
In this case, the classical FEniCS formulation for this as
nonlinear problem would be

:::{code-block} python
y = Function(V)
p = TestFunction(V)
u = Function(V)
F = inner(grad(y), grad(p)) * dx -u * p * dx
:::

which could then be solved via the {py:func}`fenics.solve` interface. This
formulation, which comes more naturally for nonlinear
variational problems (see the
[FEniCS examples](https://fenicsproject.org/docs/dolfin/latest/python/demos.html))
is closer to the one in cashocs. However, for the use with cashocs, the state variable
$y$ **must not** be a {py:class}`fenics.TrialFunction`, and the adjoint variable
$p$ **must not** be a {py:class}`fenics.TestFunction`. They **have to** be
defined as regular {py:class}`fenics.Function` objects, otherwise the code will not
work properly.
::::

After defining the weak form of the state equation, we now
specify the corresponding (homogeneous) Dirichlet boundary
conditions via

```python
bcs = cashocs.create_dirichlet_bcs(V, Constant(0), boundaries, [1, 2, 3, 4])
```

This creates Dirichlet boundary conditions with value 0 at the
boundaries 1,2,3, and 4, i.e., everywhere.

::::{hint}
Classically, these boundary conditions could also be defined via

:::{code-block} python
def boundary(x, on_bdry):
    return on_boundary
bc = DirichletBC(V, Constant(0), boundary)
:::

which would yield a single DirichletBC object, instead of the list returned by
{py:func}`create_dirichlet_bcs <cashocs.create_dirichlet_bcs>`. Any of the many
methods for defining the boundary conditions works here, as long as it is valid input
for the {py:func}`fenics.solve` function.
::::

With the above description, we see that defining the state system
for cashocs is nearly identical to defining it with FEniCS,
the only major difference lies in the definition of the state
and adjoint variables as {py:class}`fenics.Function` objects, instead of
{py:class}`fenics.TrialFunction` and {py:class}`fenics.TestFunction`.

### Definition of the cost functional

Now, we have to define the optimal control problem which we do
by first specifying the cost functional. To do so, we define the
desired state $y_d$ as an {py:class}`fenics.Expression` {python}`y_d`, i.e.,

```python
y_d = Expression("sin(2*pi*x[0])*sin(2*pi*x[1])", degree=1)
```

Alternatively, {python}`y_d` could also be a {py:class}`fenics.Function` or any other
object that is usable in an UFL form (e.g. generated with
{py:func}`fenics.SpatialCoordinate`).

Then, we define the regularization parameter $\alpha$ and the tracking-type
cost functional via the commands

```python
alpha = 1e-6
J = cashocs.IntegralFunctional(
    Constant(0.5) * (y - y_d) * (y - y_d) * dx + Constant(0.5 * alpha) * u * u * dx
)
```

The cost functional is defined via the {py:class}`cashocs.IntegralFunctional`, which
means that we only have to define the cost functional's integrand, which will then
further be treated by cashocs as needed for the optimization. These definitions are
also classical in the sense that they would have to be performed in this (or a
similar) way in FEniCS when one would want to evaluate the (reduced) cost functional,
so that we have only very little overhead.

### Definition of the optimization problem and its solution

Finally, we set up an
{py:class}`OptimalControlProblem <cashocs.OptimalControlProblem>` {python}`ocp` and
then directly solve it with the method {py:meth}`ocp.solve()
<cashocs.OptimalControlProblem.solve>`

```python
ocp = cashocs.OptimalControlProblem(e, bcs, J, y, u, p, config=config)
ocp.solve()
```

::::{hint}
Note that the {py:meth}`solve <cashocs.OptimalControlProblem.solve>` command without
any additional keyword arguments leads to cashocs using the settings defined in the
config file. However, there are some options that can be directly set with keyword
arguments for the {py:meth}`solve <cashocs.OptimalControlProblem.solve>` call. These
are
- {python}`algorithm` : Specifies which solution algorithm shall be used.
- {python}`rtol` : The relative tolerance for the optimization algorithm.
- {python}`atol` : The absolute tolerance for the optimization algorithm.
- {python}`max_iter` : The maximum amount of iterations that can be carried out.

Hence, we could also use the command

:::{code-block} python
ocp.solve('lbfgs', 1e-3, 0.0, 100)
:::

to solve the optimization problem with the L-BFGS method, a relative tolerance of
{python}`1e-3`, no absolute tolerance, and a maximum of 100 iterations.

The possible values for these arguments are the same as {ref}`the corresponding ones
in the config file <config_ocp_optimization_routine>`. This just allows for some
shortcuts, e.g., when one wants to quickly use a different solver.

Note, that it is not strictly necessary to supply config files in cashocs. In this
case, the user has to follow the above example and specify at least the solution
algorithm via the {py:meth}`solve <cashocs.OptimalControlProblem.solve>` method.
However, it is very strongly recommended to use config files with cashocs as
they allow a detailed tuning of its behavior.
::::

Finally, we visualize the results using matplotlib and the following code ::

```python
import matplotlib.pyplot as plt

plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
fig = plot(u)
plt.colorbar(fig, fraction=0.046, pad=0.04)
plt.title("Control variable u")

plt.subplot(1, 3, 2)
fig = plot(y)
plt.colorbar(fig, fraction=0.046, pad=0.04)
plt.title("State variable y")

plt.subplot(1, 3, 3)
fig = plot(y_d, mesh=mesh)
plt.colorbar(fig, fraction=0.046, pad=0.04)
plt.title("Desired state y_d")

plt.tight_layout()
# plt.savefig('./img_poisson.png', dpi=150, bbox_inches='tight')
```

The output should look like this
:::{image} /../../demos/documented/optimal_control/poisson/img_poisson.png
:::