A convex penalty for switching control of partial differential equations

TL;DR

This paper introduces a convex penalty method for switching controls in PDEs, enabling efficient numerical solutions and control structure analysis through a semismooth Newton approach.

Contribution

It proposes a novel convex penalty for switching controls in PDEs and develops an efficient solution method using Moreau-Yosida approximation and semismooth Newton techniques.

Findings

Efficient numerical solution demonstrated via examples.

Controls exhibit the desired switching structure.

The method is computationally effective for PDE control problems.

Abstract

A convex penalty for promoting switching controls for partial differential equations is introduced; such controls consist of an arbitrary number of components of which at most one should be simultaneously active. Using a Moreau-Yosida approximation, a family of approximating problems is obtained that is amenable to solution by a semismooth Newton method. The efficiency of this approach and the structure of the obtained controls are demonstrated by numerical examples.

Tables

Table 1. Table 1 : Dependence on α 𝛼 \alpha for N = 7 𝑁 7 N=7 ( γ ¯ ¯ 𝛾 \bar{\gamma} is the last successful regularization parameter, #SSN is the number of semismooth Newton steps and #CG is the number of CG steps in the last Newton step for γ ¯ ¯ 𝛾 \bar{\gamma} )

$\alpha$	${10}^{-1}$	${10}^{-2}$	${10}^{-3}$	${10}^{-4}$	${10}^{-5}$	${10}^{-6}$
${\tau }_1$	$200$	$200$	$199$	$199$	$17$	$2$
${\tau }_2$	$0$	$0$	$1$	$1$	$140$	$39$
${\tau }_3$	$0$	$0$	$0$	$0$	$43$	$146$
$\overline{\gamma }$	${10}^{-12}$	${10}^{-12}$	${10}^{-12}$	${10}^{-12}$	${10}^{-6}$	${10}^{-7}$
#SSN	$1$	$1$	$4$	$2$	$8$	$15$
#CG	$3$	$3$	$1$	$1$	$28$	$50$

Table 2. Table 2 : Convergence history of the semismooth Newton method for N = 7 𝑁 7 N=7 , α = 10 − 2 𝛼 superscript 10 2 \alpha=10^{-2} and γ = 10 − 7 𝛾 superscript 10 7 \gamma=10^{-7} , where F k subscript 𝐹 𝑘 F_{k} denotes the residual norm and s k subscript 𝑠 𝑘 s_{k} denotes the number of control intervals where the active control components change

$k$	$0$	$1$	$2$	$3$	$4$
$F_k$	$3.133\cdot {10}^{-2}$	$3.131\cdot {10}^{-2}$	$1.902\cdot {10}^{-3}$	$8.285\cdot {10}^{-6}$	$3.463\cdot {10}^{-12}$
$s_k$	N/A	$200$	$152$	$4$	$0$

Table 3. Table 3 : Dependence on γ 𝛾 \gamma for N = 7 𝑁 7 N=7 , α = 10 − 3 𝛼 superscript 10 3 \alpha=10^{-3} (#SSN is the number of semismooth Newton steps, #CG is the number of CG steps in the last Newton step)

$\gamma$	${10}^{-2}$	${10}^{-3}$	${10}^{-4}$	${10}^{-5}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$	${10}^{-12}$
${\tau }_1$	$0$	$0$	$101$	$189$	$198$	$199$	$199$	$199$
${\tau }_2$	$0$	$45$	$91$	$11$	$2$	$1$	$1$	$1$
${\tau }_3$	$0$	$101$	$8$	$0$	$0$	$0$	$0$	$0$
#SSN	$3$	$4$	$3$	$3$	$3$	$4$	$2$	$4$
#CG	$4$	$5$	$5$	$4$	$5$	$4$	$1$	$1$

Table 4. Table 4 : Dependence on γ 𝛾 \gamma for N = 7 𝑁 7 N=7 , α = 5 ⋅ 10 − 5 𝛼 ⋅ 5 superscript 10 5 \alpha=5\cdot 10^{-5} (#SSN is the number of semismooth Newton steps, #CG is the number of CG steps in the last Newton step)

$\gamma$	${10}^{-2}$	${10}^{-3}$	${10}^{-4}$	${10}^{-5}$	${10}^{-6}$	${10}^{-7}$	${10}^{-8}$
${\tau }_1$	$0$	$0$	$0$	$49$	$169$	$186$	$188$
${\tau }_2$	$0$	$0$	$0$	$126$	$28$	$14$	$12$
${\tau }_3$	$0$	$0$	$70$	$19$	$3$	$0$	$0$
#SSN	$3$	$3$	$3$	$3$	$7$	$20$	$21$
#CG	$4$	$6$	$8$	$10$	$11$	$15$	$24$

Equations

$$(Bu)(t,x)=\sum _{i=1}^N{\chi }_{{\omega }_i}(x)u_i(t),$$

$$\{\begin{array}{rl}\underset{u\in L^2(0,T;{\mathbb{R}}^N)}{\min }& \frac{1}{2}\Vert y-y^d{\Vert }_{L^2(0,T;L^2({\omega }_{\mathrm{obs}}))}^2+\frac{\alpha }{2}{\int }_0^T|u(t){|}_2^{2}dt,\\ \text{s.t.}& Ly=Bu,y(0)=y_0,\end{array}$$

$$\beta {\int }_0^T\sum _{\stackrel{i,j=1}{i<j}}^N|u_i(t)u_j(t)|dt$$

$$\{\begin{array}{rl}\underset{u\in L^2(0,T;{\mathbb{R}}^N)}{\min }& \frac{1}{2}\Vert y-y^d{\Vert }_{L^2(0,T;L^2({\omega }_{\mathrm{obs}}))}^2+\frac{\alpha }{2}{\int }_0^T|u(t){|}_1^{2}dt,\\ \text{s.t.}& Ly=Bu,y(0)=y_0,\end{array}$$

$$\underset{u}{\min }\mathcal{F}(u)+\mathcal{G}(u),$$

$$\mathcal{F}(u)=\frac{1}{2}\Vert Su-y^d{\Vert }_{L^2(0,T;L^2({\omega }_{\mathrm{obs}}))}^2,\mathcal{G}(u)=\frac{\alpha }{2}{\int }_0^T|u(t){|}_1^{2}dt,$$

$$W(0,T):=L^2(0,T;H_0^1(\Omega ))\cap W^{1,2}(0,T;H^{-1}(\Omega ))\hookrightarrow C([0,T];L^2(\Omega )).$$

$$\{\begin{array}{rl}-\bar{p}& ={\mathcal{F}}^{\prime }(\bar{u}),\\ \bar{u}& \in \partial {\mathcal{G}}^{\ast }(\bar{p}),\end{array}$$

$$\partial {\mathcal{G}}_{\gamma }^{\ast }(p):=(\partial {\mathcal{G}}^{\ast }{)}_{\gamma }(p):=\frac{1}{\gamma }\left(p-{\mathrm{prox}}_{\gamma {\mathcal{G}}^{\ast }}(p)\right),$$

$${\mathrm{prox}}_{\gamma {\mathcal{G}}^{\ast }}(v):=\underset{w\in L^2(0,T;{\mathbb{R}}^N)}{\mathrm{argmin}}\frac{1}{2\gamma }\Vert w-v{\Vert }_{L^2(0,T;{\mathbb{R}}^N)}^2+{\mathcal{G}}^{\ast }(w)={\left(\mathrm{Id}+\gamma \partial {\mathcal{G}}^{\ast }\right)}^{-1}(v)$$

$$\{\begin{array}{rl}-p_{\gamma }& ={\mathcal{F}}^{\prime }(u_{\gamma }),\\ u_{\gamma }& =\partial {\mathcal{G}}_{\gamma }^{\ast }(p_{\gamma }).\end{array}$$

$$\underset{u}{\min }\mathcal{F}(u)+{\left({\mathcal{G}}_{\gamma }^{\ast }\right)}^{\ast }(u).$$

$$0\in \partial \mathcal{F}(u_{\gamma })+\partial {\left({\mathcal{G}}_{\gamma }^{\ast }\right)}^{\ast }(u_{\gamma }).$$

$${\mathcal{G}}_{\gamma }^{\ast }={\mathcal{G}}^{\ast }\square \frac{1}{2\gamma }\Vert \cdot {\Vert }_{L^2(0,T;{\mathbb{R}}^N)}^2,$$

$$({\mathcal{G}}_{\gamma }^{\ast }{)}^{\ast }$$

$$=\mathcal{G}+\frac{\gamma }{2}\Vert \cdot {\Vert }_{L^2(0,T;{\mathbb{R}}^N)}^2.$$

$$\mathcal{G}(u)={\int }_0^{T}g(u(t))dt,$$

$$g:{\mathbb{R}}^N\to \mathbb{R},g(v)=\frac{\alpha }{2}|v{|}_1^2.$$

$$g^{\ast }:{\mathbb{R}}^N\to \mathbb{R},g^{\ast }(q)=\frac{1}{2\alpha }|q{|}_{\infty }^2=\underset{1\le i\le N}{\max }\frac{1}{2\alpha }{q}_i^2.$$

$$\partial g^{\ast }(q)=\mathrm{co}\left(\bigcup _{\{i:g^{\ast }(q)=g_i^{\ast }(q)\}}\left\{(g_i^{\ast }{)}^{\prime }(q)\right\}\right),$$

$$\partial g^{\ast }(q)=\{\begin{array}{ll}\left\{\left(\frac{1}{\alpha }{q}_1,0\right)\right\}& \text{if }|q_1|>|q_2|,\\ \left\{\left(0,\frac{1}{\alpha }{q}_2\right)\right\}& \text{if }|q_1|<|q_2|,\\ \left\{t\left(\frac{1}{\alpha }{q}_1,0\right)+(1-t)\left(0,\frac{1}{\alpha }{q}_2\right):t\in [0,1]\right\}& \text{if }|q_1|=|q_2|.\end{array}$$

$$[\partial g^{\ast }(q){]}_j=\{\begin{array}{ll}\left\{\frac{1}{\alpha }{q}_j\right\}& \text{if }|q_j|={\max }_i|q_i|\text{ and }|q_j|>|q_i|\text{ for }i\ne j,\\ \left\{0\right\}& \text{if }|q_j|<{\max }_i|q_i|,\\ \left\{\frac{s_j}{\alpha }{q}_j:s_j\ge 0,{\sum }_{i\in \mathcal{A}}{s}_i=1\right\}& \text{if }|q_j|={\max }_i|q_i|\text{ and }\exists i\ne j:|q_j|=|q_i|,\end{array}$$

$$\partial {\mathcal{G}}^{\ast }(p)=\left\{u\in L^2(0,T;{\mathbb{R}}^N):u(t)\in \partial g^{\ast }(p(t))\text{ for a.e. }t\in (0,T)\right\}.$$

$${\bar{u}}_j(t)\in \{\begin{array}{ll}\left\{\frac{1}{\alpha }{\bar{p}}_j(t)\right\}& \text{if }|{\bar{p}}_j(t)|={\max }_i|{\bar{p}}_i(t)|,|{\bar{p}}_j(t)|>|{\bar{p}}_i(t)|\text{ for }i\ne j,\\ \left\{0\right\}& \text{if }|{\bar{p}}_j(t)|<{\max }_i|{\bar{p}}_i(t)|,\\ \left\{\frac{s_j}{\alpha }{\bar{p}}_j(t):s_j\ge 0,{\sum }_{i\in \mathcal{A}(t)}{s}_i=1\right\}& \text{if }|{\bar{p}}_j(t)|={\max }_i|{\bar{p}}_i(t)|,\exists i\ne j:|{\bar{p}}_j(t)|=|{\bar{p}}_i(t)|,\end{array}$$

$$\mathcal{S}:=\left\{t\in (0,T):|\mathcal{A}(t)|>1\right\},$$

$$|v_{d+1}|<\frac{\alpha }{d\alpha +\gamma }\sum _{i=1}^d|v_i|$$

$$[{\mathrm{prox}}_{\gamma g^{\ast }}(v){]}_j=\{\begin{array}{ll}\mathrm{sign}(v_j)\frac{\alpha }{d\alpha +\gamma }{\sum }_{i=1}^d|v_i|& \text{if }j\le d,\\ v_j& \text{if }j>d.\end{array}$$

$$v\in (\mathrm{Id}+\gamma \partial g^{\ast })(w)=\{w\}+\gamma \partial g^{\ast }(w).$$

$$w_i=\{\begin{array}{ll}\frac{\alpha }{\alpha +\gamma }{v}_i& \text{if }i=1,\\ v_i& \text{if }i>1.\end{array}$$

$$|v_1|>\left(1+\frac{\gamma }{\alpha }\right)|v_2|.$$

$$\mathrm{sign}(v_1)|v_1|$$

Full text

1 Introduction

Switching control refers to time-dependent optimal control problems with a vector-valued control of which at most one component should be active at every point in time. We focus here on optimal tracking control for a linear diffusion equation $Ly=Bu$ on $\Omega_{T}:=(0,T]\times\Omega$, $y(0)=y_{0}$ on $\Omega$, where $L=\partial_{t}-A$ for an elliptic operator $A$ defined on $\Omega\subset\mathbb{R}^{n}$ carrying suitable boundary conditions. The control operator $B$ is defined by

[TABLE]

where $\chi_{\omega_{i}}$ is the characteristic function of the given control domain $\omega_{i}\subset\Omega$ of positive measure. Furthermore, let $\omega_{\mathrm{obs}}\subset\Omega$ denote the observation domain and let $y^{d}\in L^{2}(0,T;L^{2}(\omega_{\mathrm{obs}}))$ denote the target. We consider the standard optimal control problem

[TABLE]

where $|v|_{2}^{2}=\sum_{j=1}^{N}v_{j}^{2}$ denotes the squared $\ell^{2}$-norm on $\mathbb{R}^{N}$. To promote the switching structure of the optimal control $\bar{u}\in L^{2}(0,T;\mathbb{R}^{N})$, we suggest adding an additional penalty term

[TABLE]

with $\beta>0$ to the objective, which can be interpreted as an $L^{1}$-penalization of the switching constraint $u_{i}(t)u_{j}(t)=0$ for $i\neq j$ and $t\in[0,T]$. For the choice $\beta=\alpha$, the sum of the control cost and the penalty can be simplified to yield the problem

[TABLE]

where $|v|_{1}=\sum_{j=1}^{N}|v_{j}|$ denotes the $\ell^{1}$-norm on $\mathbb{R}^{N}$. This is a convex optimization problem, for which we derive first-order optimality conditions in primal-dual form whose Moreau–Yosida regularization (to be introduced below) can be solved using a superlinearly convergent semismooth Newton method. The effect of other choices for $\beta$ will be discussed in Section 3.3.

The approach which we follow here is related to both the switching control problem in [5] and the distributed parabolic sparse control problem in [4]. In [5] a nonconvex formulation in the case where $N=2$ was considered; we compare its convex relaxation to the present approach below. One advantage of the approach presented in this paper over that in [5] is given by the fact that there is no significant additional technical burden when considering switching between $N>2$ controls. In [4] the $L^{2}$ norm in time of the measure norm in space was used to promote temporally varying sparsity in space. While the choice of the nonsmooth functional involving the controls in (P) is motivated by sparsity considerations, one can arrive at this functional also from controllability–observability considerations. In fact, it was shown in [23, Theorem 4.1] that – provided an appropriately defined controllability Gramian has full rank – exact null controls with perfect switching have minimal $L^{2}(0,T;\mathbb{R}^{2})$ norm, where $\mathbb{R}^{2}$ is endowed with the $\ell^{1}$ norm. In contrast, [9] follows a different approach where binary or integer decision variables are sought within a relaxation technique combined with a suitable rounding strategy.

Let us comment on further related work. While our work here aims at formulating optimal controls problems with switching controls in a way that allows an efficient numerical treatment, the larger body of work focuses on the stabilization of switching systems. For ordinary differential equations we refer to e.g., [3, 14, 19]. For partial differential equations, this problem has received comparatively little attention. In both cases, one should distinguish switching control in the sense defined above from the control of switched systems. For the latter in the context of PDEs, we refer to [8, 20]. In [10], converse Lyapunov theorems for abstract switched systems are developed. Lyapunov techniques are also used to study switches in hyperbolic systems in [17], and existence results for optimal control of switching systems modelling the use of bacteria for pollution removal are obtained in [18]. Exact null controls with switching structure for the heat and wave equation were treated in [23, 15, 22] and [7], respectively.

This work is organized as follows. Section 2 discusses the existence and first-order optimality conditions for solutions to (P) as well as its regularization within an abstract convex analysis framework. Explicit pointwise characterizations of the switching relations arising from the optimality system and for its regularization are given in Section 3. Here we also discuss the relation of the proposed switching functional in (P) to other possible choices of the penalty term. Section 4 is concerned with the numerical solution of the regularized optimality system by a semismooth Newton method. A numerical example for switching control of a two-dimensional linear heat equation is computed in Section 5.

2 Convex analysis approach

We recall the convex analysis approach for (nonconvex) switching controls for partial differential equations from [5], which is also applicable to the convex penalty considered here. For this purpose, we consider Problem (P) in the reduced form

[TABLE]

with $\mathcal{F}:L^{2}(0,T;\mathbb{R}^{N})\to\mathbb{R}$ and $\mathcal{G}:L^{2}(0,T;\mathbb{R}^{N})\to\mathbb{R}$ given by

[TABLE]

where the continuous affine solution operator $S:u\mapsto y$ assigns to any control $u\in L^{2}(0,T;\mathbb{R}^{N})$ the unique state $y\in L^{2}(\Omega_{T})$ satisfying the state equation $Ly=Bu$ with initial condition $y(0)=y_{0}$ subject to appropriate boundary conditions. Here we assume that the coefficients of $A$, the boundary and initial conditions as well as the domain $\Omega$ are sufficiently regular that the range of $S$ is contained in

[TABLE]

Since $S$ is affine, $\mathcal{F}$ is proper, convex and lower semicontinuous. Furthermore, since the squared norm $|\cdot|_{1}^{2}$ is convex, $\mathcal{G}$ is proper, convex, lower semicontinuous, and, in addition, radially unbounded. Existence of a solution thus follows from standard arguments, e.g., Tonelli’s direct method.

Proposition 2.1**.**

There exists a minimizer $\bar{u}$ for Problem (P).

We next derive first-order optimality conditions in primal-dual form. Throughout, for any proper convex function $\mathcal{H}$, we denote by $\mathcal{H}^{*}$ its Fenchel conjugate and by $\partial\mathcal{H}$ its subdifferential; see, e.g., [1, 6] for their definitions. The following proposition is a direct consequence of the sum rule and inversion formula for convex subdifferentials (see, e.g., [1, Corollary 16.24] for the latter) as well as the Fréchet-differentiability of $\mathcal{F}$.

Proposition 2.2**.**

The control $\bar{u}\in L^{2}(0,T;\mathbb{R}^{N})$ is a minimizer for (P) if and only if there exists a $\bar{p}\in L^{2}(0,T;\mathbb{R}^{N})$ such that

[TABLE]

holds.

Since $\mathcal{F}$ is a standard quadratic tracking term, the first relation in (OS) can be expressed in a straightforward manner in terms of the solution operator $S=L^{-1}B$ and its adjoint $S^{*}=B^{*}L^{-*}$ (with homogeneous boundary and initial conditions), i.e., $\bar{p}=-S^{*}(S\bar{u}-y^{d})$. For later use, we point out that due to (6) and the specific choice of $S$ there holds $\bar{p}\in V:=B^{*}(W(0,T))\hookrightarrow L^{r}(0,T;\mathbb{R}^{N})$ for any $r>2$.

The second relation is responsible for the switching structure of the optimal control $\bar{u}$, and we will give a pointwise characterization in Proposition 3.1 below.

Our numerical approach is based on the Moreau–Yosida regularization of (OS). Specifically, we replace $\partial\mathcal{G}^{*}$ for $\gamma>0$ by

[TABLE]

where

[TABLE]

is the proximal mapping of $\mathcal{G}^{*}$, which in Hilbert spaces coincides with the resolvent of $\partial\mathcal{G}^{*}$. Note that the proximal mapping and thus the Moreau–Yosida regularization of a proper and convex functional is always single-valued and Lipschitz continuous; see, e.g., [1, Corollary 23.10].

We then consider the regularized system

[TABLE]

Again, we will give an explicit formulation of (OSγ) in the next section.

Proposition 2.3**.**

For each $\gamma>0$ system (OSγ) admits a unique solution $(u_{\gamma},p_{\gamma})$.

Proof 2.4**.**

Using convex analysis techniques (see, e.g., [1, Chapter 12]), we obtain that (OSγ) is the necessary optimality condition for

[TABLE]

In fact, since $\mathcal{F}$ is globally defined and continuous, we have the necessary optimality condition

[TABLE]

Setting $p_{\gamma}:=-\mathcal{F}^{\prime}(u_{\gamma})$ we recall again the subdifferential inversion formula from [1, Corollary 16.24]. Additionally, we note that the Yosida-regularization $\partial\mathcal{G}^{*}_{\gamma}$ of the subdifferential $\partial\mathcal{G}^{*}$ coincides with the Fréchet derivative of the Moreau-envelope $\mathcal{G}^{*}_{\gamma}$ of $\mathcal{G}^{*}$; see, e.g., [1, Proposition 12.29]. Together, these yield (OSγ).

From [1, Remark 12.24], we have the alternative characterization of the Moreau-envelope via the infimal convolution

[TABLE]

where $\mathbb{R}^{N}$ is endowed with the Euclidean norm. Furthermore, from [1, Theorem 15.3], we obtain

[TABLE]

This implies that $(\mathcal{G}^{*}_{\gamma})^{*}$ is strictly convex. Moreover, $\mathcal{G}^{*}_{\gamma}\leq\mathcal{G}^{*}$ and hence $0\leq\mathcal{G}^{**}=\mathcal{G}\leq(\mathcal{G}^{*}_{\gamma})^{*}$. Therefore a minimizing sequence for (Pγ) is necessarily bounded, and by a weak subsequential limit argument, existence of a unique solution $u_{\gamma}$ to (Pγ) follows. Together with $p_{\gamma}=-\mathcal{F}^{\prime}(u_{\gamma})$, we thus obtain a unique solution $(u_{\gamma},p_{\gamma})$ to (OSγ).

Finally, convergence as $\gamma\to 0$ can be shown as in [5, Proposition 2.5]. This requires showing that $\{\partial\mathcal{G}^{*}(p_{\gamma})\}_{\gamma>0}$ is uniformly bounded in $\gamma$ provided that $\{p_{\gamma}\}_{\gamma>0}$ is bounded. But this follows from the explicit characterization in (17) below.

Proposition 2.5**.**

For any sequence $\gamma_{n}\to 0$, the sequence $\{(u_{\gamma_{n}},p_{\gamma_{n}})\}_{n\in\mathbb{N}}$ converges weakly to a solution $(\bar{u},\bar{p})$ to (OS).

3 Switching penalty

Here we compute an explicit, pointwise, characterization of $\partial\mathcal{G}^{*}$ and $\partial\mathcal{G}^{*}_{\gamma}$ by exploiting the fact that for integral functionals of the form

[TABLE]

the Fenchel conjugate and convex subdifferential can be computed pointwise via $g$; see, e.g., [6, Props. IV.1.2, IX.2.1], [1, Prop. 16.50]. We thus focus on

[TABLE]

Because the Fenchel conjugate of half the squared norm is half the squared dual norm (see, e.g., [2, Example 3.27]), a scaling argument (e.g., [1, Proposition 13.20 (i)]) yields

[TABLE]

3.1 Subdifferential

Since $g^{*}$ is the maximum of finitely many convex and differentiable functions, its convex subdifferential is given by

[TABLE]

where ${\mathrm{co}}$ denotes the closed convex hull and $g_{i}^{*}(q)=\frac{1}{2\alpha}q_{i}^{2}$; see, e.g., [11, Corollary 4.3.2]. It is instructive to first consider the case $N=2$, for which only the three cases $g^{*}(q)=g_{1}(q)\neq g_{2}(q)$, $g^{*}(q)=g_{2}(q)\neq g_{1}(q)$ and $g^{*}(q)=g_{1}(q)=g_{2}(q)$ need to be considered. The corresponding convex hulls are then given by

[TABLE]

For the general case, we proceed similarly by considering all possible cases. It is then straightforward to verify that the subdifferential can be given componentwise for $1\leq j\leq N$ as

[TABLE]

where $\mathcal{A}:=\left\{j:|q_{j}|=\max_{i}|q_{i}|\right\}$. We thus obtain for $u,p\in L^{2}(0,T;\mathbb{R}^{N})$ the pointwise characterization

[TABLE]

In particular, this yields a pointwise characterization of the second relation in (OS).

Proposition 3.1**.**

The minimizer $\bar{u}\in L^{2}(0,T;\mathbb{R}^{N})$ of (P) and $\bar{p}:=-\mathcal{F}^{\prime}(\bar{u})\in L^{2}(0,T;\mathbb{R}^{N})$ satisfy for almost every $t\in(0,T)$ and $1\leq j\leq N$

[TABLE]

where $\mathcal{A}(t):=\left\{j:|\bar{p}_{j}(t)|=\max_{i}|\bar{p}_{i}(t)|\right\}$.

From this proposition, the desired switching property follows: Outside of the singular arc

[TABLE]

we have $\bar{u}_{j_{1}}\bar{u}_{j_{2}}=0$ for all $j_{1}\neq j_{2}$. In particular, if $\mathcal{S}$ has Lebesgue measure zero, $\bar{u}$ exhibits perfect switching, i.e., $\bar{u}_{j_{1}}\bar{u}_{j_{2}}=0$ almost everywhere. Under the stronger assumption that $|\bar{p}_{j_{1}}(t)|\neq|\bar{p}_{j_{2}}(t)|$ for all $j_{1}\neq j_{2}$ and almost all $t\in(0,T)$, a similar relation for exact null controls was given in [23, equations (2.36–2.38)].

Remark 3.2**.**

It is difficult to give general conditions for the optimal control to be perfectly switching; instead one would have to exploit properties of the specific (adjoint) state equation. For example, in the case of the one-dimensional heat equation, one could use the fact that solutions are real-analytic with respect to time to argue that $|\bar{p}_{i}(t)|=|\bar{p}_{j}(t)|$ for all $t\in O$ for an open set $O\subset(0,T)$ is only possible if $\bar{p}_{i}(t)=\bar{p}_{j}(t)=0$ for all $t\in O$; cf. [23, Section 6.1].

3.2 Proximal mapping

We next characterize the proximal mapping $\mathrm{prox}_{\gamma g^{*}}(p)$ via the resolvent $(\operatorname{\mathrm{Id}}+\gamma\partial g^{*})^{-1}(p)$ by proceeding similarly to [5, Section 3.3]; see also [13].

Proposition 3.3**.**

Let the components of $v\in\mathbb{R}^{N}$ be sorted by decreasing magnitude and let $d$ be the smallest index for which

[TABLE]

holds; if no such index exists, let $d=N$. The explicit form of the proximal mapping is then given componentwise for $1\leq j\leq N$ by

[TABLE]

Proof 3.4**.**

Let $w:=(\mathrm{Id}+\gamma\partial g^{*})^{-1}(v)$. This is characterized by the subdifferential inclusion

[TABLE]

First, we observe that (23) implies that $\operatorname{\mathrm{sign}}(v_{i})=\operatorname{\mathrm{sign}}(w_{i})$. We next follow the case discrimination in the characterization of the subdifferential, where we assume that the components of $w$ are ordered by magnitude, i.e., $|w_{1}|\geq|w_{2}|\geq\dots\geq|w_{N}|$ (the remaining cases following by permutation).

1. (i)

$|w_{1}|>|w_{2}|$: In this case, all components of the subdifferentials are single-valued; solving for each $w_{i}$ in (23) yields

[TABLE]

The assumption $|w_{1}|>|w_{2}|$ is then equivalent to the condition

[TABLE] 2. (ii)

$|w_{1}|=|w_{2}|>|w_{3}|$: Using that $\operatorname{\mathrm{sign}}(w_{i})=\operatorname{\mathrm{sign}}(v_{i})$, this implies that $w_{1}=\operatorname{\mathrm{sign}}(v_{1})\bar{w}$ and $w_{2}=\operatorname{\mathrm{sign}}(v_{2})\bar{w}$ for some $\bar{w}>0$. Inserting this into the set-valued case of the subdifferential, we deduce that for some $t\in[0,1]$,

[TABLE]

Dividing by the $\operatorname{\mathrm{sign}}(v_{i})$ on both sides and adding yields

[TABLE]

from which we obtain

[TABLE]

the single-valued cases for $i>2$ following as before.

Inserting these values into $|w_{2}|>|w_{3}|$ then yields the condition

[TABLE]

Note that in addition,

[TABLE]

must hold since due to case (i), the converse is equivalent to $|w_{1}|>|w_{2}|$, in contradiction to the assumption. 3. (iii)

$|w_{1}|=|w_{2}|=|w_{3}|>|w_{4}|$: As before, the chain of equalities is equivalent to

[TABLE]

for some $t_{1},t_{2}$ with $t_{1}+t_{2}\in(0,1)$. Adding these equations yields that

[TABLE]

as well as the condition

[TABLE]

Again,

[TABLE]

must hold since otherwise we obtain a contradiction in the case distinction. 4. (iv)

The remaining cases $|w_{1}|=\dots=|w_{j}|>|w_{j+1}|$ for $j=4,\dots,N$ follow analogously, establishing the desired characterization of the proximal mapping.

It is now straight-forward to give a pointwise characterization of the Moreau–Yosida regularization $\partial\mathcal{G}^{*}_{\gamma}$ via

[TABLE]

We point out that for every $t\in(0,T)$, the number $d$ in Proposition 3.3 is the number of nonzero components in $u_{\gamma}(t)=[\partial\mathcal{G}^{*}_{\gamma}(p)](t)$. In particular for the case $d=1$, the explicit form of $\partial g^{*}_{\gamma}$ is given componentwise as

[TABLE]

(Note that the missing case $|q_{j}|\leq\left(1+\tfrac{\gamma}{\alpha}\right)|q_{i}|$ is excluded by the assumption that $d=1$.) Thus if for some $\gamma>0$, the solution $(u_{\gamma},p_{\gamma})$ to (OSγ) is such that $d=1$ for almost every $t\in(0,T)$, the regularized control is perfectly switching.

3.3 Relation to other approaches for switching control

In this section, we address the relation of the proposed approach to alternative $\ell^{1}$-type penalizations as well as to the $\ell^{0}$ formulation of [5]. For this purpose, we recall the original introduction of $g$ as a combination of a quadratic control cost and an $\ell^{1}$-penalization of the switching constraint, i.e., as

[TABLE]

We now compare this formulation with other approaches, setting $N=2$ for the sake of simplicity.

First, we address the choice $\beta\neq\alpha$. For $\beta<\alpha$, we can introduce $\gamma:=\alpha-\beta>0$ and rewrite $g_{\beta}$ as

[TABLE]

(with $\beta$ in place of $\alpha$ in the definition of $g$). The optimality system in this case is therefore equivalent to (OSγ). In case $\beta>\alpha$, we can reformulate again to

[TABLE]

which is obviously nonconvex. Computing the Fenchel conjugate leads to

[TABLE]

which coincides with (14). We thus have that $g_{\beta}^{*}=g^{*}$, from which it follows that

[TABLE]

i.e., $g$ is the lower convex envelope of the nonconvex function $g_{\beta}$.

We can also compare to the nonconvex $\ell^{0}$-switching formulation from [5],

[TABLE]

(where $|s|_{0}=0$ if $s=0$ and $1$ otherwise), whose Fenchel conjugate is given by

[TABLE]

For fixed $q$ and $\beta$ sufficiently large, this coincides with (16), since the third case (corresponding to the free arc, where both control components are active) is never attained. Therefore, the functional $g$ can be interpreted as the convex relaxation of $g_{0}$ in the limit as $\beta\to\infty$, in the following sense: Taking the formal limit $\beta\to\infty$ in (45), we arrive at the (still nonconvex) constrained functional

[TABLE]

where the switching property $v_{i}v_{j}=0$ for all $i\neq j$ is included as an explicit constraint. Proceeding as in [5, § 3.2], we find that the Fenchel conjugate of $g_{\infty}$ is given by the first two cases of (46). As in the case of $g_{\beta}$, this implies that $g$ is the lower convex envelope of $g_{\infty}$.

We conclude that our choice of $g$ is a natural convex formulation for promoting switching controls. Convexification naturally entails that switching cannot be guaranteed for arbitrarily small choices of $\alpha$ but is crucial for an efficient numerical solution as discussed in the next section.

4 Numerical solution

We apply a semismooth Newton method (see, e.g., [13, 21]) to (OSγ). By eliminating $u_{\gamma}$ via the second relation of (OSγ) and introducing the state $y_{\gamma}:=S(u_{\gamma})\in Y:=L^{2}(\Omega_{T})$, we arrive at the equivalent system

[TABLE]

where we set $H_{\gamma}:=\partial\mathcal{G}_{\gamma}^{*}$. We similarly set $h_{\gamma}:=\partial g_{\gamma}^{*}$ such that

[TABLE]

Exactly as in [5, § 5], one argues that $h_{\gamma}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is semismooth, and a Newton derivative $D_{N}h_{\gamma}(q)\in\mathbb{R}^{N\times N}$ is obtained by choosing any element of the Clarke derivative $\partial_{C}h_{\gamma}(q)$ (which for the piecewise differentiable function $h_{\gamma}$ is given by the convex hull of the piecewise derivatives). Here we take, assuming again for simplicity that the components of $q$ are sorted descending by magnitude, componentwise for $1\leq i,j\leq N$

[TABLE]

where $d$ is as in Proposition 3.3. The Newton derivative for arbitrary $q\in\mathbb{R}^{N}$ can be obtained from this by appropriate permutation of rows and columns. Note that $D_{N}h_{\gamma}(q)$ is symmetric. Furthermore, we point out that evaluation of the proximal mapping amounts to sorting for each $t\in(0,T)$ the vector $p(t)\in\mathbb{R}^{N}$, which can be carried out in $\mathcal{O}(N\log N)$ operations, and testing the $N-1$ conditions in (21) for the sorted vector, which requires $\mathcal{O}(N)$ operations. Hence, evaluating $h_{\gamma}$ and assembling the Newton derivative $D_{N}h_{\gamma}$ can be performed in $\mathcal{O}(N\log N)$ operations, thus avoiding the exponential complexity involved in computing switching points between any two (or more) control components; see, e.g., [12]. A different relaxation approach for mixed-integer optimal control that also avoids exponential complexity was presented in [9].

Since $h_{\gamma}$ is semismooth from $\mathbb{R}^{N}$ to $\mathbb{R}^{N}$, the corresponding superposition operator $H_{\gamma}$ is semismooth from $V\hookrightarrow L^{r}(0,T;\mathbb{R}^{N})$ for $r>2$ to $L^{2}(0,T;\mathbb{R}^{N})$, with Newton derivative at $p$ in direction $\delta p$ given by

[TABLE]

Considering the system (48) as an operator equation from $Y\times V$ to $Y\times V$, a semismooth Newton step for its solution thus consists of computing $(\delta y,\delta p)\in Y\times V$ for given $(y^{k},p^{k})\in Y\times V$ such that

[TABLE]

and setting $y^{k+1}=y^{k}+\delta y$ and $p^{k+1}=p^{k}+\delta p$. Here, $S_{0}$ denotes the solution operator of the state equation with homogeneous initial and boundary conditions. As in [5, Proposition 5.2], one shows uniformly bounded invertibility of each Newton step, from which locally superlinear convergence of the semismooth Newton method follows.

The disadvantage of (52) is the forward–backward structure in time, which prevents using time-stepping methods for the application of $S_{0}$ and $S^{*}$. For the practical implementation, we therefore eliminate $\delta y$ via the first relation of (52) as well as $y^{k}$ via the first relation of (48) to obtain the reduced semismooth Newton step

[TABLE]

which allows for time-stepping and matrix-free Krylov methods. The linear operator on the left-hand side becomes self-adjoint using the inner product generated by $D_{N}H_{\gamma}(p^{k})$, and hence we apply a matrix-free CG method for the solution of (53) using this inner product; see [16, Chapter 3.4.1 and Remark 3.9 ii)]. We point out that monitoring the convergence of the CG method is still based on the standard Euclidean norm of the residual.

An alternative to (53) is a reduction to $\delta y$, which would lead to a reduced Newton matrix that is already self-adjoint with respect to the standard inner product. We do not follow this approach, because it would lead to a significantly higher dimension of the Newton system (full state space versus control space). In addition, the iterates would be infeasible, since in this case the state equation is only satisfied for a solution of the reduced optimality system.

To account for the local convergence of Newton methods, we embed the semismooth Newton method within a homotopy strategy for $\gamma$, where we start with a large $\gamma$ which is successively reduced, taking the previous solution as starting point. Furthermore, we include a backtracking line search based on the residual norm to improve robustness.

5 Numerical examples

We now illustrate the behavior of the proposed approach and the structure of the resulting controls using a family of numerical examples. We consider the heat equation in two dimensions, i.e., $y=S(u)$ solves

[TABLE]

with $T=10$ and $\Omega=(-1,1)^{2}$, and zero initial and Neumann boundary conditions. We choose circular control domains $\omega_{i}=B_{0.1}(x_{i})$ with centers $x_{i}=(\cos(\varphi_{i}),\sin(\varphi_{i}))/\sqrt{2}$ regularly distributed along a circle, where $\varphi_{i}=\pi/4+2\pi(i-1)/N$, while the observation domain is $\omega_{\mathrm{obs}}=B_{0.5}(0)$; see Fig. 1.

The desired state $y^{d}\in L^{2}(0,T;L^{2}(\omega_{\mathrm{obs}}))$ is set to

[TABLE]

The corresponding optimal controls are computed with the semismooth Newton method (53), stopping at a relative tolerance of $10^{-6}$ in the residual $F_{k}=\|F(p^{k})\|_{L^{2}(0,T;\mathbb{R}^{N})}$, where $F(p)$ denotes the right-hand side of (53). The backtracking line search starts with a step size $1$ and uses a reduction factor of $1/2$. The CG method is stopped with a relative tolerance of $10^{-6}$ or after a maximal number of $50$ iterations is reached. We initialize with $p^{0}=0$. The homotopy loop starts from $\gamma=10^{-2}$, reducing $\gamma$ by a factor of $10$ as long as the Newton method converges in a prescribed number of maximal $30$ iterations; the minimal allowed $\gamma$ is $10^{-12}$. We denote by $\bar{\gamma}$ the resulting smallest $\gamma$ for which the Newton method converged. Below, we also report on the number $\tau_{1}$ of control intervals that exhibit perfect switching, i.e., in which at most one control of $N$ is active. Similarly, we define $\tau_{j}$ for $j>1$ as the number of time points with $d=j$ in the proximal mapping.

For the spatial discretization, piecewise linear finite elements on an unstructured mesh of $725$ triangular elements (maximal diameter $0.1$) are chosen. The time discretization is carried out with the cG($1$) Petrov–Galerkin method (which corresponds to the Crank–Nicolson method) and $201$ equidistant time points. The discrete control is chosen as piecewise constant, i.e., there are $200$ degrees of freedom per control component. An implementation of the presented approach in Matlab can be downloaded from https://github.com/clason/multiswitchingcontrol.

Figure 2(a), Fig. 2(b) and Fig. 3(a) show the optimal controls for $\alpha=10^{-1}$ and $N=3$, $N=5$ and $N=7$, respectively (where here and below only active (i.e., nonzero) control components are shown; components that are never active are not listed in the legend, i.e., $u_{2}$, $u_{4}$, $u_{7}$ in these cases, respectively). All three configurations terminated with $\bar{\gamma}=10^{-12}$, and the optimal controls exhibit perfect switching. Note that while all control components are discontinuous, the norm $\|u(t)\|_{l^{1}(\mathbb{R}^{N})}$ (which corresponds to the absolute value of the curve plotted in the figure) is continuous. This is expected since the $L^{2}$ norm of this quantity appears as a penalty in Problem (P).

Controls for $N=7$ and different values of $\alpha$ are depicted in Fig. 3. Since $\alpha$ determines the strength of the switching penalty, it is not surprising that the controls are perfectly switching only for sufficiently large $\alpha$, as seen in Fig. 3(a) and Fig. 3(b). For moderate values of $\alpha$ such as $\alpha=10^{-3}$, there is a single isolated control point where two components are active at the same time ($u_{3}$ and $u_{7}$ in Fig. 3(c)). For much smaller $\alpha=10^{-5}$, perfect switching is mostly lost ($\tau_{1}=17$); however, the control still exhibits a switching property, since in most intervals only two and never more than three components are active ($\tau_{2}=140$, $\tau_{3}=43$). Due to the penalization of the $L^{2}$ norm in time of the $\ell_{1}$ norm of the controls, we also observe that the magnitude of the plotted control envelope increases with decreasing $\alpha$ while also becoming less regular with respect to time. For those controls exhibiting perfect switching, the number of switching points (i.e., of time points $\tau\in(0,T)$ where $\arg\max_{j}|\bar{p}_{j}(s)|\neq\arg\max_{j}|\bar{p}_{j}(t)|$ for $s<\tau<t$ close to $\tau$) is relatively independent of $\alpha$, rising only slightly from $10$ points for $\alpha=10^{-1}$ to $11$ and $12$ points for $\alpha=10^{-2}$ and $\alpha=10^{-3}$, respectively.

The dependence of the switching structure on $\alpha$ is illustrated quantitatively in Table 1, where we also report on the effects on the numerical algorithm. As can be seen, for moderate values of $\alpha$, the homotopy terminates successfully with $\bar{\gamma}=10^{-12}$, and both the number of Newton and CG iterations in the final homotopy and Newton step, respectively, increase only slightly. For smaller $\alpha$ and $\gamma$, the Newton system becomes increasingly difficult to solve, and the inner CG and Newton iterations terminate due to the maximum number of iterations having been reached, which leads to early termination of the homotopy strategy.

To illustrate the convergence rate of the semismooth Newton method, we report in Table 2 on the convergence history for a single run of the Newton method for $N=7$, $\alpha=0.01$ and $\gamma=10^{-7}$ without homotopy. Both the norm of the residual in the reduced first order optimality condition (53) and the number of control intervals where the active control components change (denoted by $F_{k}$ and $s_{k}$, respectively) show superlinear decay. It can also be observed that the reduced optimality conditions are satisfied close to machine precision as soon as the active components no longer change. This finite termination of semismooth Newton methods is typical for problems that are linear apart from the case distinction in $H_{\gamma}$ and shows the semismooth Newton method’s relation to active set methods. The described behavior is typical for large $\alpha$ or $\gamma$ where no line search is necessary. For moderate values, the iteration first requires some reduced steps before switching to full steps with superlinear convergence and finite termination. For very small $\alpha$ and $\gamma$, the prescribed maximal number of iterations is insufficient to reach this convergence region, and the method terminates.

The behavior of the homotopy strategy for moderate $\alpha=10^{-3}$ is illustrated in Table 3, where both the number of control intervals $\tau_{1},\tau_{2},\tau_{3}$ with one, two or three active control components, respectively, and the number of Newton iterations together with the number of CG iterations in the last Newton step are shown. For large $\gamma$, the regularization naturally prevents strong switching, which starts to appear for $\gamma<10^{-4}$; for $\gamma<10^{-8}$, the structure of active components no longer changes. It can also be observed that the number of Newton and CG steps stays relatively constant throughout the whole iteration. The situation is different for small $\alpha$; see Table 4. Since the regularized subdifferential depends on the relation between $\alpha$ and $\gamma$, the switching structure starts to appear later. Depending on the problem, the limit case for small $\alpha,\gamma$ might have a solution that is not perfectly switching, in which case the unregularized subdifferential is set-valued. In this example, we observe that for small $\alpha$ and $\gamma$, the Newton steps indeed become increasingly difficult to solve, as shown by the growing number of steps necessary before the iteration terminates.

6 Conclusion

It is possible to promote optimal controls for parabolic differential equations with switching structure using a convex penalty, avoiding the need to introduce additional decision variables or explicit switching points. A combination of Moreau–Yosida regularization and a semismooth Newton method allows the efficient numerical computation of controls for an arbitrary number of control components. Since the proposed switching penalty is independent of the state equation, our approach can be adapted to hyperbolic and nonlinear problems. Furthermore, additional sparsity or control constraints can be incorporated into the presented convex analysis framework. At the same time, the convexification is unable to guarantee perfect switching for arbitrary choice of $\alpha$. Analytical issues including the dependency of the optimal solution on $\alpha$, the number of switches for the optimal control, and sufficient conditions ensuring perfect switching of the optimal control depend on the specific structure of the problem setting and are of interest for future research.

Acknowledgment

This work was supported in part by the Austrian Science Fund (FWF) under grant SFB F32 (SFB “Mathematical Optimization and Applications in Biomedical Sciences”).