sNewton4PDEOpt

The Python module implements a semismooth Newton method for solving finite-element discretizations of the strongly convex, linear elliptic PDE-constrained optimization problem

$\begin{align*} \min_{u \in U_{\text{ad}}}\, (1/2) \|S(u)-y_d\|_{L^2(D)}^2 + (\alpha/2)\|u\|_{L^2(D)}^2 + \beta \|u\|_{L^1(D)}, \end{align*}$

where $S(u) \in H_0^1(D)$ solves the linear elliptic PDE:

$\begin{align*} \text{find} \quad y \in H_0^1(D) \quad \text{with} \quad (\kappa \nabla y,\nabla v)_{L^2(D)^d} = (u+g,v)_{L^2(D)} \quad \text{for all} \quad v \in H_0^1(D). \end{align*}$

The feasible set is $U_{\text{ad}} = \{u \in L^2(D): \lb \leq u \leq \ub \}$$

The control space is discretized using piecewise constant functions and the state space is discretized using piecewise continuous functions. FEniCS is used to preform the discretization. Sparse linear systems are solved using npsolve of scipy.

The parameter alpha must be positive and beta be nonnegative. The domain D is (0,1)^d with d being one or two. The diffusion coefficient kappa maps from the domain to the positive reals. The lower and upper bounds lb and ub, the desired state yd, the diffusion coefficient kappa, and the source term g must be instances of either Constant, Function, or Expression.

The semismooth Newton method is applied a normal map, a reformulation of the first-order optimality conditions as a nonsmooth operator equation (see eq. (3.3) in Ref. [3]).

The module implements a globalization via a restricted residual monotonicity test (see MonotonicityTest), while NewtonStep chooses the step size equal to one. The implementation of the restricted residual monotoncitity test is based on eq. (3.32) in Ref. [1].

Installation

The code can be downloaded using git clone.

Examples

example1, example2, example3, and example4 implement Examples 1–4 in Ref. [4].
poisson1d and possion2d use an L1-heuristic to determine beta (compare with p. 199 in Ref. [2] and Lemma 3.1 in Ref. [4]).
example73 implements Example 7.3 in Ref. [5] and it is used for code verification.
random_poisson2d demonstrates how to use the module to solve sample average approximations of a simple risk-neutral problem.
cdc provides an example where globalization using MonotonicityTest requires fewer iterations than NewtonStep.

Dependencies

Some tests use the following packages:

The code was tested using python version 3.8.11, scipy version 1.7.0, numpy version 1.21.1, fenics version 2019.1.0, moola version 0.1.6, and matplotlib version 3.4.2.