Parabolic optimal control problems with combinatorial switching constraints, part III: branch-and-bound algorithm

Abstract

We present a branch-and-bound algorithm for globally solving parabolic optimal control problems with binary switches that have bounded variation and possibly need to satisfy further combinatorial constraints. More precisely, for a given tolerance  $\varepsilon >0$, we show how to compute in finite time an $\varepsilon$-optimal solution in function space, independently of any prior discretization. The main ingredients in our approach are an appropriate branching strategy in infinite dimension, an a posteriori error estimation in order to obtain safe dual bounds, and an adaptive refinement strategy in order to allow arbitrary switching points in the limit. The performance of our approach is demonstrated by extensive experimental results.

1 Introduction

Optimal control problems with discrete switches have recently become an increasing focus of research. Most approaches presented in the literature, however, produce only heuristic solutions without any quality guarantee. The well-known Sum-Up Rounding approach [28, 33] computes binary switching patterns by first solving a convex relaxation of the problem and then approximating the resulting continuous switching by a binary one. This approach often requires a large number of switchings when trying to come close to the optimal continuous solution. In particular, it cannot deal with an explicit bound on the number of switchings, let alone with more complex combinatorial constraints. If such constraints need to be satisfied, the Combinatorial Integral Approximation approach [32] can be applied. Since the latter again tries to approximate a given continuous control by a feasible binary one, this approach does not lead to optimal solutions to the original problem in general, even when a best-possible approximation can be computed. Other approaches aim at optimizing the switching times of the discrete switches directly [13, 14, 18, 25, 30, 31, 37], or use non-smooth penalty techniques, partly in combination with convexification, to impose the switching structure, see, e.g., [9,10,11,12, 40] and the references therein. However, both strategies in general lead to non-convex problems with potentially multiple local minima and a convexifcation of the arising problems may destroy the switching structure of the optimal solution.

In this paper, we present a branch-and-bound approach for solving parabolic optimal control problems with combinatorial switching constraints to global optimality. More precisely, we consider problems of the form

Here $T > 0$ is a given final time and $\Omega \subset {\mathbb {R}}^d$, $d\in {\mathbb {N}}$, is a bounded domain, i.e., a bounded, open, and connected set, with Lipschitz boundary $\partial \Omega$ in the sense of [24, Def. 1.2.2.1]. The form function $\psi \in H^{-1}(\Omega )$ and the initial state $y_0 \in L^2(\Omega )$ are given. Moreover, $y_{d } \in L^2(Q)$ is a given desired state and $\alpha > 0$ is a Tikhonov parameter weighting the deviation from $\tfrac{1}{2}$. Finally,

$$\begin{aligned} D \subseteq BV(0,T;\{0,1\}):=\big \{ u \in BV(0,T):u(t) \in \{0,1\} \text { f.a.a.\ } t \in (0,T)\big \} \end{aligned}$$

denotes the set of feasible switching controls and is supposed to satisfy the following assumptions:

Here BV(0, T) denotes the set of all functions in $L^1(0,T)$ with bounded variation, equipped with the norm $\Vert u\Vert _{BV(0,T)}:= \Vert u\Vert _{L^1(0,T)} + |u|_{BV(0,T)}$; see e.g., [1] for details on the space of functions with bounded variation. For simplicity, we restrict ourselves to the case of one binary switch in (P), but our main results are easily extended to the case of multiple switches, i.e., to $D\subset BV(0,T;\{0,1\}^k)$ for some $k\in {\mathbb {N}}$.

Assumption (D1) is crucial in our context. Without this condition, it would be possible to approximate any control u with $u(t)\in [0,1]\text { f.a.a.\ } t \in (0,T)$ arbitrarily well by binary switch es using an increasing number of switchings, e.g., by applying the Sum-Up Rounding approach mentioned above. Indeed, for a given control $u\in L^2(0,T;[0,1])$, there exists a sequence $\{u^n\}_{n\in {\mathbb {N}}}$ of binary controls in $BV(0,T;\{0,1\})$ such that $u^n$ converges weakly in $L^2(0,T)$ to u, and this weak convergence is sufficient to pass to the limit in the state equation. Moreover, in our branch-and-bound algorithm, the existence of a uniform bound on the variation is essential in order to ensure the effectiveness of fixings.

The Tikhonov term $\Vert u-\tfrac{1}{2}\Vert ^2_{L^2(0,T)}$ in the objective function of (P) is constant on D by construction, but not on the relaxations of D considered in the course of the algorithm. In [6, 7], we present tailored convexifications of (P) and an outer approximation approach to solve the resulting relaxations. The core of the approach is the generation of linear cutting planes describing the closed convex hull of D in $L^p(0,T)$. While the overall approach is very general, the specific shape of the cutting planes is problem-dependent. In particular, we devise results for the case of bounded total variation (without further constraints) and for the case where the switching points of the control u must satisfy given linear constraints. The latter case comprises the so-called mininum dwell time contraints [42].

Building on the results of [6, 7], our aim is thus to determine globally optimal solutions for problems of type (P) that are independent of any prior discretization, using a branch-and-bound approach. For this, we start from the convex relaxations of (P) studied in [6]. These convex relaxations correspond to the root nodes in our branch-and-bound algorithm. In order to extend this to a full branch-and-bound algorithm for computing globally optimal solutions (at least in the limit), we have to overcome several obstacles:

• Since we optimize in function space, fixing the value of the switch in finitely many points (as is common in finite-dimensional branch-and-bound algorithms) has no direct effect, or is not even well-defined. We thus have to take the bounded variation into account in order to obtain implicit restrictions on the set of admissible controls in the nodes of the branch-and-bound tree; see Sect. 3.

• The fixing of the switch at certain points in time leads to a non-closed set of admissible controls in the nodes. Moreover, the closed convex hulls of these sets are structurally different from the admissible controls arising in the root node. We study the most important classes of these sets in Sect. 4.

• In [7], we devised a semi-smooth Newton method to solve the root relaxation. However, with an increasing number of fixings, this method became less stable, so that we now propose to solve all subproblems by the alternating direction method of multipliers; see Sect. 5.

• In order to obtain globally optimal solutions, all dual bounds computed in the nodes of the branch-and-bound tree must be safe. In particular, they need to take discretization errors into account. In case the time-mesh independent dual bound is too weak to cut off a node, we may either have to branch or to refine the temporal grid, depending on the relation between the current primal bound and the time-mesh dependent dual bound. The sophisticated interplay between branching, error analysis, and adaptive refinement is at the core of our proposed approach, it is discussed in Sect. 6.

The main contribution of this paper is to present solutions for the challenges listed above. An extensive experimental evaluation presented in Sect. 7 shows that an effective and stable implementation of the resulting branch-and-bound approach is possible.

2 Preliminaries

We first collect some definitions and observations that are needed in the following sections, concerning the solution mapping for the PDE in (P) as well as functions of bounded variation and special switching constraints.

2.1 Solution mapping

The assumptions on the optimal control problem (P) listed above guarantee that, for every control $u\in D\subset L^2(0,T)$, the PDE in (P) admits a unique weak solution

$$\begin{aligned} y\in W(0,T):= H^1(0,T;H^{-1}(\Omega )) \cap L^2(0,T; H^1_0(\Omega )); \end{aligned}$$

see [38, Chapter 3]. To specify the associated solution operator

$$\begin{aligned} S:L^2(0,T)\ni u \mapsto y \in W(0,T), \end{aligned}$$

we introduce the linear and continuous (and thus Fréchet differentiable) operator

$$\begin{aligned} \Psi :L^2(0,T) \rightarrow L^2(0,T;H^{-1}(\Omega )), \quad (\Psi u) (t,x) = u(t) \psi (x) \end{aligned}$$

as well as the solution operator $\Sigma : L^2(0,T;H^{-1}(\Omega )) \rightarrow W(0,T)$ of the heat equation with homogeneous initial condition, i.e., given $w\in L^2(0,T;H^{-1}(\Omega ))$, $y = \Sigma (w)$ solves

$$\begin{aligned} \partial _t y - \Delta y = w \quad \text {in } L^2(0,T;H^{-1}(\Omega )), \quad y(0) = 0 \quad \text {in } L^2(\Omega ). \end{aligned}$$

Moreover, we introduce the function $\zeta \in W(0,T)$ as the solution for

$$\begin{aligned} \partial _t \zeta - \Delta \zeta = 0 \quad \text {in } L^2(0,T;H^{-1}(\Omega )), \quad \zeta (0) = y_0 \quad \text {in } L^2(\Omega ). \end{aligned}$$

Then the solution operator S is given by $S = \Sigma \circ \Psi + \zeta$. In particular, it is affine and continuous. Using this solution operator, the problem (P) can be written as

2.2 Functions of bounded variation

Functions of bounded variation are of central importance in the following. In order to deal with such functions, first recall that each function $u \in BV(0,T)$ admits a right-continuous representative given by ${{\hat{u}}}(t) = c+\mu ((0,t])$, $t\in (0,T)$, where $\mu$ is the regular Borel measure on [0, T] associated with the distributional derivative of u and $c\in {\mathbb {R}}$ is a constant. Note that ${{\hat{u}}}$ is unique on (0, T). Here and in the following, with a slight abuse of notation, we denote this function by the same symbol as the equivalence class in BV(0, T) and, when it comes to pointwise evaluations, we always refer to this representative function. In particular, we will often write constraints in the form $u(t)=b$ for $b\in {\mathbb {R}}$. For $t\in (0,T)$, this is well-defined by the above reasoning, it then means ${{\hat{u}}}(t)=b$, while for $t=0$, we use the same notation as shorthand for $\lim _{t\searrow 0}{{\hat{u}}}(t)=b$.

In this paper, we will mostly deal with binary controls $u\in BV(0,T;\{0,1\})$. In this case, the representative ${{\hat{u}}}$ can be parameterized through its switching points $0\le t_1\le \cdots \le t_\sigma$, where $\sigma \le {|u|}_{BV(0,T)}$. More formally, if one already counts ${{\hat{u}}}(0)=\lim _{t\searrow 0} {\hat{u}}(t)=1$ as one switching from 0 to 1, then the representative can be written in the form

$$\begin{aligned} u_{t_1,\dots ,t_\sigma }(t):= {\left\{ \begin{array}{ll} 0, & \text {if }\,|\{i \in \{1, \ldots , \sigma \} :\, t_i \le t \}|\text { is even},\\ 1, & \text {if\,} |\{i \in \{1, \ldots , \sigma \} :\, t_i \le t \}|\text { is odd}; \end{array}\right. } \end{aligned}$$

see [6] for more details. In the following, we will always regard $u_{t_1,\dots ,t_\sigma }(t)$ as a function in BV(0, T).

2.3 Examples of switching constraints

In [6, 7], we present tailored convexifications of (P’) and an outer approximation algorithm to solve the resulting relaxations. In particular, we elaborate the details of this approach for the following two relevant classes of constraints D. By Assumption (D1), the total number of switchings is bounded, and the first type of constraint arises when this is the only restriction. More specifically, we restrict the total variation of the single switch from above by $\sigma >0$, so that the set of feasible controls is

$$\begin{aligned} D(\sigma ):=\{u\in BV(0,T):u(t)\in \{0,1\} \text { f.a.a. }t\in (0,T),\ |u|_{BV(0,T)}\le \sigma \}. \end{aligned}$$

(1)

The second type of constraint imposes affine linear relations between the positions of the switching points of u. More precisely, for a given polytope $P\subseteq {\mathbb {R}}_+^\sigma$, we define

$$\begin{aligned} D(P):= \{ u_{t_1,\dots ,t_\sigma } :(t_1, \ldots , t_\sigma ) \in P,~ 0\le t_1\le \cdots \le t_\sigma < \infty \}. \end{aligned}$$

(2)

An important special case are the so-called minimum dwell time constraints, defined as

$$\begin{aligned} D(s):= \big \{ u_{t_1,\dots ,t_\sigma }:t_{i}-t_{i-1}\ge s ~ \forall \,i = 2, \dots ,\sigma ,~ t_1,\ldots ,t_\sigma \ge 0\big \} \end{aligned}$$

(3)

for some given $s>0$. In words, a minimum time span s is required between two consecutive switchings of u.

All sets $D(\sigma )$ and D(P) defined here, and hence also D(s), satisfy the general assumptions (D1) and (D2); see [6, Lemma 3.7 and Lemma 3.10]. We will investigate these sets under fixings in Sect. 4 and use $D(\sigma )$ for our experiments presented in Sect. 7.

3 Branch-and-bound algorithm

In finite-dimensional optimization, branch-and-bound is the standard approach for solving non-convex optimization problems to global optimality. First, a dual bound is computed for the original problem, corresponding to the root node of the branch-and-bound tree. Often, this is done by solving a convex relaxation of the problem. In case the optimal solution for the latter is infeasible for the original problem, a branching is applied. In the most abstract form, this means that the set of feasible solutions is subdivided into two (or more) subsets, corresponding to the child nodes of the root node. Recursive application of the branching leads to the so-called branch-and-bound tree. The bounding is now applied in order to reduce the number of nodes in this tree, which leads to a finite algorithm in many cases: one first needs to obtain so-called primal solutions, i.e., feasible solutions of the original problem. Each such solution yields (in case of a minimization problem) a global upper bound on the optimal value of the original problem. Now, if the dual bound obtained in some branch-and-bound node is larger than the best known upper bound, it follows that this node cannot contain any optimal solution, so that it can be pruned, i.e., the entire subtree rooted at this node can be ignored in the enumeration.

The most natural branching strategy for finite-dimensional binary optimization problems consists of picking a binary variable having a fractional value in the optimal solution for the convex relaxation used for computing the dual bound, and then fixing this variable to zero in the first child node and to one in the other. However, in the infinite-dimensional setting considered here, the situation is more complicated: we need to deal with infinitely many binary variables, suggesting that an infinite number of function values has to be fixed in order to uniquely determine a solution for (P). In fact, fixing a pointwise value of u has no direct effect (or is not even well-defined) in the function space $L^p(0,T)$. At this point, we can exploit Assumption (D1), which yields a finite bound on the total number of switching points. The relevant restrictions in a given node of the branch-and-bound tree are now a joint consequence of the finitely many fixing decisions taken so far and of the constraint $u\in D$.

The main challenge is now to describe these resulting restrictions. Assume that our branching strategy always picks appropriate time points $\tau \in (0,T)$ and fixes $u(\tau )=0$ in the first subproblem and $u(\tau )=1$ in the second. Then all our subproblems, corresponding to the nodes in the branch-and-bound tree, are problems in BV(0, T) of the form

with $(\tau _j,c_j) \in [0,T)\times \{0,1\}$ for $1\le j \le N$; see Sect. 2.2 for the precise meaning of the fixing constraints. In the following, we denote the feasible set of (SP) by

$$\begin{aligned} D_ SP :=\{u\in D: u(\tau _j) = c_j \ \forall j=1,\ldots ,N\}, \end{aligned}$$

where we always assume $\tau _1<\cdots <\tau _N$. Note that the set $D_ SP$ is not closed in $L^p(0,T)$ in general, and hence the subproblem (SP) does not necessarily admit a global minimizer. However, this is no problem since we are only interested in the optimal value of (SP) in our branch-and-bound framework. In fact, our approach will produce a series of dual bounds by convexifying (SP) and these co nvexifications will provide the same (primal) optimal value of (SP) in the limit; see Theorem 3.1 below. We consider the convexification

of the suproblem (SP) in $L^p(0,T)$. Here and in the following, ${\overline{\operatorname {conv}}}$ always denotes the closed convex hull in $L^p(0,T)$.

In a reasonable branching strategy, one may expect that an increasing number of fixing decisions, taken along a path in the branch-and-bound tree starting at the root node, leads to a unique solution in the limit. In particular, the dual bounds obtained in the nodes and the optimal values subject to the corresponding fixings should converge to each other. The next result shows that this is guaranteed in our infinite-dimensional setting if the fixing positions are sufficiently well-distributed.

Theorem 3.1

For $N\in {\mathbb {N}}$, let $0\le \tau _1^{N}<\dots<\tau _N^N< T$ and define

$$\begin{aligned} \Delta \tau ^{N}:= \max _{j=1,\ldots ,N+1}|\tau _j^{N}-\tau _{j-1}^{N}|, \end{aligned}$$

where $\tau _0^N:=0$ and $\tau _{N+1}^N:=T$. If $\Delta \tau ^{N} \rightarrow 0$ for $N\rightarrow \infty$, then

• (i) the diameters of the feasible sets of (SPC) and (SP) in $L^2(0,T)$ vanish and

• (ii) the optimal values of (SPC) and (SP) converge to each other.

Proof

Let $D_ SP ^{N}:=\{u\in D: u(\tau _j^{N}) = c_j^{N} \ \forall j=1,\ldots ,N\}$ denote the feasible set of (SP) for $N\in {\mathbb {N}}$. Without loss of generality, we may assume $D_ SP ^N\ne \emptyset$ for $N\in {\mathbb {N}}$, since otherwise the feasible set of (SPC) is also empty and thus both optimal values agree. We first claim that two controls $u_1,u_2\in D_ SP ^{N}$ can only differ in at most $\sigma$ of the intervals $(\tau _{j-1}^{N},\tau _j^{N})$ for $2\le j\le N$, where $\sigma$ denotes the upper bound on the total number of switchings guaranteed by Assumption (D1). Indeed, assume that $u_1$ and $u_2$ differ between $\tau _{j-1}^{N}$ and $\tau _j^{N}$. Since the values of $u_1$ and $u_2$ agree at $\tau _{j-1}^{N}$ and $\tau _j^{N}$, either one of the two functions has to switch at least twice in $(\tau _{j-1}^{N},\tau _j^{N})$, if $c_{j-1}^{N}=c_j^{N}$, or both functions have to switch at least once, if $c_{j-1}^{N}\ne c_j^{N}$. Hence, for each interval where $u_1$ and $u_2$ differ, both functions together have at least two switchings, but the total number of their switchings is bounded by $2\sigma$.

Taking into account also the intervals $(0,\tau _1^{N})$ and $(\tau _N^{N},T)$ and using that $u_1,u_2\in [0,1]$ a.e. in (0, T), we thus obtain

$$\begin{aligned} \sup _{u_1,u_2\in D_ SP ^{N}} \Vert u_1-u_2\Vert ^2_{L^2(0,T)} \le (\sigma +2)\Delta \tau ^{N} \end{aligned}$$

and consequently, for $N\rightarrow \infty$, we get

$$\begin{aligned} \sup _{u_1,u_2\in {\overline{\operatorname {conv}}}(D_ SP ^{N})} \Vert u_1-u_2\Vert _{L^2(0,T)} = \sup _{u_1,u_2\in D_ SP ^{N}} \Vert u_1-u_2\Vert _{L^2(0,T)} \rightarrow 0, \end{aligned}$$

(4)

which shows assertion (i).

Since the solution operator S in (SP) is linear and bounded, the objective $J:u \mapsto J(Su,u)$ is continuous from $L^2(0,T)$ to ${\mathbb {R}}$. Together with (4), this implies that the maximal difference of all objective values of feasible controls in (SPC) vanishes for $N\rightarrow \infty$. Since (SPC) is a relaxation of (SP), we obtain (ii). $\square$

As a consequence of Theorem 3.1 and its proof, we immediately obtain the following.

Corollary 3.2

For each $\varepsilon >0$ there exist $N\in {\mathbb {N}}$ and fixings $(\tau _j,c_j) \in (0,T)\times \{0,1\}$, $j=1,\ldots ,N$, such that the optimal value of (SPC) differs by at most $\varepsilon$ from the optimal value of the original problem (P).

In other words, up to an arbitrary desired precision $\varepsilon >0$, the optimal solution of (P) can be approximated by (SPC) using a finite number of fixings. This is crucial for the branch-and-bound algorithm we are going to present in the following. Clearly, the number of necessary fixings depends on $\varepsilon$.

To solve the subproblems in the branch-and-bound algorithm, we will use the outer approximation approach presented in [6, 7]. For this purpose, we need to discuss how to deal with the resulting projections under fixings (see Sect. 4) and how to adapt the outer approximation algorithm (see Sect. 5). For both tasks, first note that

$$\begin{aligned} {\overline{\operatorname {conv}}}(D_ SP )={\overline{\operatorname {conv}}}(\overline{D_ SP }), \end{aligned}$$

again with all closures taken in $L^p(0,T)$. While the inclusion “$\subseteq$” is obvious, it suffices for the other direction to show that $\overline{D_\text {SP}}$ is a subset of $\overline{\text {conv}}(D_\text {SP})$, as the latter is closed and convex by definition. Now if u is the limit of a sequence $\{u^m\}_{m\in {\mathbb {N}}}$ in $D_\text {SP}$, then $u^m\in {\text {conv}}(D_\text {SP})$ for all $m\in {\mathbb {N}}$ and hence $u\in \overline{\text {conv}}(D_\text {SP})$. Problem (SPC) is thus very similar to the problem without fixings addressed in [6, 7], except that D is now replaced by the more complex set $\overline{D_ SP }$. In an outer approximation approach, the impact of the fixings is then implicitly modeled by the cutting planes describing ${\overline{\operatorname {conv}}}(D_ SP )$.

However, the fixings may also directly determine significant parts of the switching pattern in such a way that u must be constant on some intervals $[\tau _{j-1},\tau _j)$, i.e., $u|_{[\tau _{j-1},\tau _j)}\equiv c_{j-1}$ for all controls $u\in D_ SP$. Indeed, as shown by the proof of Theorem 3.1, the non-fixed part of the time horizon vanishes under the assumptions of Theorem 3.1 when $N\rightarrow \infty$. In other words, the share of the time horizon (0, T) on which u is completely determined by the branching decisions converges to one. In our branch-and-bound algorithm, it is much more efficient to deal with these constraints explicity, instead of modeling them by cutting planes describing ${\overline{\operatorname {conv}}}(D_ SP )$.

Finally, note that it is also possible that the given fixings are inconsistent with the constraint D, i.e., that the feasible set of (SP) is empty, which is easy to detect for most choices of D. In this case, the subproblem is infeasible and the corresponding node in the branch-and-bound tree can be pruned.

Example 3.3

Consider the set $D(\sigma )$ defined in (1). Let $N'=|\{j\in \{2,\dots ,N\}:c_{j-1}\ne c_{j}\}|$. If $N'>\sigma$, we have $D(\sigma )_ SP =\emptyset$, since even the number of switchings enforced by the fixing is too large for a feasible solution. The subproblem can thus be pruned. If $\sigma -1\le N'\le \sigma$, we can fix all intervals $[\tau _{j-1},\tau _j)$ with $c_{j-1}=c_j$ to the value $c_{j-1}$, since any other value in this interval would increase the number of switchings by two. If $\sigma -1\le N'\le \sigma$ and $c_{j-1}\ne c_j$, then no value of u in $(\tau _{j-1},\tau _j)$ is fixed, but u has to be monotone in $[\tau _{j-1},\tau _j]$, which is modeled implicitly by cutting planes. The same is true for all further restrictions resulting from the fixings. $\square$

Example 3.4

For the minimum dwell time constraints D(s) defined in (3), we can fix an interval $[\tau _{j-1},\tau _j)$ with $c_{j-1}=c_j$ to the value $c_{j-1}$ if and only if $\tau _j-\tau _{j-1}\le s$. Otherwise, no direct fixing is possible, but the number of allowed switchings within the interval $(\tau _{j-1},\tau _j)$ reduces to $\lceil \nicefrac {\tau _j-\tau _{j-1}}{s}\rceil$. An infeasible subproblem arises whenever u is fixed to the same value at two time points having a distance of at most s, but fixed to the other value at some point in between. $\square$

4 Convex hull under fixings

As already indicated, our aim is to fully describe the convex hull of feasible switching patterns, i.e., the feasible set of (SPC), by cutting planes derived from finite-dimensional projections, extending the approach proposed in [6] for the case without fixings. For this, we project the set $\overline{D_ SP }$ to the finite-dimensional space ${\mathbb {R}}^M$, by means of local averaging

$$\begin{aligned} \Pi :BV(0,T) \ni u \mapsto \Big (\tfrac{1}{\lambda (I_i)}\int _{I_i} u(t) \,\text {d}t\Big )_{i=1}^M\in {\mathbb {R}}^{M}, \end{aligned}$$

(5)

where $I_i\subseteq (0,T)$ for $i=1,\ldots ,M$ are suitably chosen subintervals. Each projection $\Pi$ then gives rise to a relaxation

$$\begin{aligned} {\overline{\operatorname {conv}}}(\overline{D_ SP }) \subseteq \{ v \in L^p(0,T):\Pi (v) \in C_{\overline{D_ SP },\Pi }\} \end{aligned}$$

of the feasible region [6, Lemma 3.4], where $C_{\overline{D_ SP },\Pi }:= \operatorname {conv}\{\Pi (u):u \in \overline{D_ SP }\}$. By a suitable construction of projections $\Pi _k$, with increasing dimension $M_k$, a complete outer description of the finite-dimensional convex hulls $C_{\overline{D_ SP },\Pi }$ also yields a complete outer description of the convex hull of $\overline{D_ SP }$ in function space [6, Thm. 3.5], i.e.,

$$\begin{aligned} {\overline{\operatorname {conv}}}(\overline{D_ SP }) =\bigcap _{k\in {\mathbb {N}}} \{ v \in L^p(0,T):\Pi _k(v) \in C_{\overline{D_ SP },\Pi _k}\}. \end{aligned}$$

In order to solve the convexified subproblem (SPC) by means of the outer approximation algorithm presented in [7, Alg. 1], it is particularly desirable that the sets $C_{\overline{D_ SP },\Pi }$ are polyhedra for which the separation problem is tractable, in order to efficiently generate cuts of the form $a^\top \Pi (u)\le b$ for $u\in L^p(0,T)$, where $a^\top w\le b$, $a\in {\mathbb {R}}^M$, $b\in {\mathbb {R}}$ represents a valid inequality for $C_{\overline{D_ SP },\Pi }$. For prominent examples of D, it is shown in [6] that this is the case for the sets $C_{D,\Pi }$, i.e., when no fixings are considered. However, it can be shown that the fixings may destroy this property in general.

For the two classes of constraints D already discussed in [6] and defined in Sect. 2.3, one can show with a similar reasoning as in [6] that the corresponding projection sets $C_{\overline{D(\sigma )_ SP },\Pi }$ and $C_{\overline{D(P)_ SP },\Pi }$ are still polyhedra for a fixed projection $\Pi$. For convenience of the reader, we include full proofs of both results in Appendix A. In the remainder of this section, we focus on the tractability of the projection sets, since the latter is crucial for the efficient computation of dual bounds by means of outer approximation within our branch-and-bound algorithm. For this, we consider a fixed projection $\Pi$ and assume that the intervals $I_i$, $i=1,\dots ,M$, are pairwise disjoint. Moreover, without loss of generality, we may assume that the fixing points $0\le \tau _1< \cdots< \tau _N < T$ satisfy $\tau _j\notin I_i$ for all $j=1,\ldots ,N$ and $i=1,\ldots ,M$, since otherwise one can refine the projection intervals and thus generate stronger cutting planes [7, Thm. 2.2].

4.1 Restricted total variation

In this section, we show that the separation problem for $C_{\overline{D(\sigma )_ SP },\Pi }$ can be solved in polynomial time. Without any fixings, i.e., when $D(\sigma )_ SP =D(\sigma )$, it is shown in [6, Thm. 3.8] that the separation problem for $C_{D(\sigma ),\Pi }$ reduces to the separation problem of

$$\begin{aligned} \operatorname {conv}\Big \{v\in \{0,1\}^M: \sum _{l=2}^M |v_l-v_{l-1}|\le \sigma \Big \}. \end{aligned}$$

For the slightly different setting where $v_1$ is fixed to zero, it is shown in [8] that the separation problem for the above set can be solved in polynomial time. It is easy to see that the separation problem remains tractable also without this fixing, i.e., for $C_{D(\sigma ),\Pi }$. Our aim is now to efficiently reduce the separation problem for $C_{\overline{D(\sigma )_ SP },\Pi }$ to the separation problem for $C_{D(\sigma ),\Pi }$. To this end, we extend the vector v by the fixings $c_1,\ldots ,c_N$. More precisely, for all $j\in \{1,\dots ,N\}$, let $i_j\in \{0,\ldots ,M\}$ be the index such that $b_{i_{j-1}}\le \tau _j\le a_{i_j}$ holds, where $b_0:=0$ and $I_i=(a_i,b_i)$ for $i=1,\ldots ,M$. In addition, define $E:{\mathbb {R}}^M \rightarrow {\mathbb {R}}^{M+N}$ by

$$\begin{aligned} Ev=(v_1,\ldots ,v_{i_1},c_1,v_{i_1+1},\ldots ,v_{i_2},c_2,v_{i_2+1},\ldots ,v_{i_N},c_N,v_{i_N+1},\ldots ,v_M)^\top . \end{aligned}$$

The desired reduction is now based on the fact that, analogously to [6, Thm. 3.8], one can again show that $C_{\overline{D(\sigma )_ SP },\Pi }$ is the convex hull of all projection vectors resulting from feasible controls that are constant almost everywhere on each of the intervals $I_1,\ldots ,I_N$, i.e., $C_{\overline{D(\sigma )_ SP },\Pi }=\operatorname {conv}(K)$ with

$$\begin{aligned} \begin{aligned} K:=\{\Pi (u) :&u \in \overline{D(\sigma )_ SP } \text { and for all }i=1,\ldots ,M \text { there exists } w_i\in \{0,1\}\\ &\text {with } u(t)\equiv w_i \text { f.a.a.\ } t\in I_i\}, \end{aligned} \end{aligned}$$

see Theorem A.1, and on the following:

Lemma 4.1

A vector $v\in {\mathbb {R}}^M$ belongs to K if and only if Ev belongs to

$$\begin{aligned} {{{\mathcal {C}}}}:=\Big \{w\in \{0,1\}^{M+N}: \sum _{l=2}^{M+N}|w_l-w_{l-1}|\le \sigma \Big \}. \end{aligned}$$

Proof

For the first direction, let $v=\Pi (u)\in K$ for some $u\in \overline{D(\sigma )_ SP }$ being constant almost everywhere on each projection interval. Then there exists a sequence $\{u^m\}_{m\in {\mathbb {N}}}\subseteq D(\sigma )_ SP$ with $u^m\rightarrow u$ in $L^p(0,T)$ for $m\rightarrow \infty$. For every $m\in {\mathbb {N}}$, the control $u^m$ has at most $\sigma$ switchings and satisfies $u^m(\tau _j)=c_j$ for $j=1,\ldots ,N$, so that we have

$$\begin{aligned} \sum _{l=2}^{M+N} |E\Pi (u^m)_l-E\Pi (u^m)_{l-1}|\le \sigma . \end{aligned}$$

The continuity of $\Pi$ in $L^p(0,T)$ yields $v=\Pi (u)=\lim _{m\rightarrow \infty } \Pi (u^m)$ and hence

$$\begin{aligned} \sum _{l=2}^{M+N} |Ev_l-Ev_{l-1}|\le \lim _{m\rightarrow \infty }\sum _{l=2}^{M+N} |E\Pi (u^m)_l-E\Pi (u^m)_{l-1}| \le \sigma , \end{aligned}$$

i.e., we have $Ev\in {{{\mathcal {C}}}}$ as desired.

We next show the opposite direction. So, let $Ev\in {{{\mathcal {C}}}}$ for some vector $v\in {\mathbb {R}}^M$. In addition, let $0=z_0<z_1<\cdots <z_r=T$ include all endpoints of the intervals $I_1\ldots ,I_M$ and the fixed positions $\tau _1,\ldots ,\tau _N$. Construct functions $u^m$ for $m\in {\mathbb {N}}$ such that

$$\begin{aligned} \begin{array}{ll} u^m(t)=v_{i} & \ \text{ for } t\in [a_i+\tfrac{\lambda (I_i)}{2m},b_i-\tfrac{\lambda (I_i)}{2m}) \text{ and } i=1,\ldots ,N,\\ u^m(t)=c_j & \ \text{ for } t\in [\tau _j,\tau _j+\tfrac{\varepsilon _j}{2m}) \text{ and } j=1,\ldots ,N, \end{array} \end{aligned}$$

where $\varepsilon _j=\min \{|z_i-\tau _j|: i\in \{1,\ldots ,r\}, z_i\ne \tau _j\}>0$. For points in (0, T) not covered by the above intervals, we copy the value of the left neighboring interval. The construction is illustrated in Fig. 1a. We have $u^m(\tau _j)=c_j$ for every $m\in {\mathbb {N}}$ and $j=1,\ldots ,N$, hence all fixings are respected. Moreover, $Ev\in {{{\mathcal {C}}}}$ guarantees that $u^m$ switches at most $\sigma$ times, i.e., we get $u^m\in D(\sigma )_ SP$. By copying always the value of the left neighboring interval, we guarantee that the control functions $u^m$ converge in $L^p(0,T)$ to some function u; see Fig. 1b.

Fig. 1

Illustration of the second part of the proof of Lemma 4.1

Moreover, by construction, the limit u is $v_1,\ldots ,v_M$ almost everywhere on the projection intervals $I_1,\ldots ,I_M$, respectively, and due to $\{u^m\}\subseteq D(\sigma )_ SP$, we have $u\in \overline{D(\sigma )_ SP }$. Therefore $v=\Pi (u)\in K$. $\square$

Theorem 4.2

The separation problem for $C_{\overline{D(\sigma )_ SP },\Pi }$ can be solved in polynomial time.

Proof

Using Lemma 4.1 and $C_{\overline{D(\sigma )_ SP },\Pi }=\operatorname {conv}(K)$, we obtain that $v\in C_{\overline{D(\sigma )_ SP },\Pi }$ if and only if $Ev\in \operatorname {conv}(\mathcal{C})$. The separation problem for $C_{\overline{D(\sigma )_ SP },\Pi }$ thus reduces to the separation problem for $\operatorname {conv}({{{\mathcal {C}}}})$. $\square$

The separation algorithm used in the outer approximation approach devised in [6] even computes the most violated cutting plane. The same can be done when considering fixings: our aim is thus to find the most violated cutting plane $({\bar{a}},{\bar{b}})\in {\mathbb {R}}^{M+1}$ in the set

$$\begin{aligned} H_{\overline{D(\sigma )_ SP },\Pi }=\{(a,b)\in [-1,1]^M\times {\mathbb {R}}: a^\top w \le b \ \forall w\in C_{\overline{D(\sigma )_ SP },\Pi }\} \end{aligned}$$

of all valid inequalities for $C_{\overline{D(\sigma )_ SP },\Pi }$, where $a\in [-1,1]^M$ can be assumed without loss of generality by scaling. It is easily verified that this aim can be achieved by first computing the most violated cutting plane $(a,b)\in {\mathbb {R}}^{M+N+1}$ for Ev in

$$\begin{aligned} H_{{{{\mathcal {C}}}}}=\{(a,b)\in [-1,1]^{M+N}\times {\mathbb {R}}: a^\top w \le b \ \forall w\in {{{\mathcal {C}}}}\}\ \end{aligned}$$

and then replacing the $(i_j+j)$ th variable by the constant $c_j$ for all $j=1,\dots ,N$, i.e., ${\bar{a}}$ results from a by deleting the $(i_j+j)$ th variables and ${\bar{b}}=b-\sum _{j=1}^N a_{i_j+j}\,c_j$. The first task can again be reduced to the case without fixings.

4.2 Dwell time constraints

We now focus on the special case of dwell time constraints, as defined in (3). Here, a minimum time span $s>0$ between two switchings is required. For the case without fixings, it is stated in [6, Thm. 3.13] that there exists a separation algorithm with polynomial running time in M and in the implicit bound $\sigma =\lceil \tfrac{T}{s}\rceil$ on the number of allowed switchings. In the presence of fixings $(\tau _j,c_j)$, $1\le j\le N$, we show in the following that the separation problem for $C_{\overline{D(s)_ SP },\Pi }$ is still tractable. More precisely, we claim that there exists a separation algorithm with polynomial time in M, $\sigma$, and the number N of fixings.

We thus consider the set

$$\begin{aligned} \begin{aligned} D(s)_ SP := \big \{ u_{t_1,\dots ,t_\sigma }:&t_{i}-t_{i-1}\ge s\; \forall \,i = 2, \dots ,\sigma ,~t_1,\ldots ,t_\sigma \ge 0,\\&u_{t_1,\dots ,t_\sigma }(\tau _j)=c_j\; \forall \,j = 1, \dots ,N\big \} \end{aligned} \end{aligned}$$

and first argue, similar to [6, Sect. 3.2], that it suffices to consider as switching points the finitely many points in the set

$$\begin{aligned} S:=[0,T]\cap \Big ({\mathbb {Z}}s+\big (\{0,T\}\cup \{a_i,b_i:i=1,\dots ,M\}\cup \{\tau _j:j=1,\dots ,N\}\big )\Big ) \end{aligned}$$

where $I_i=(a_i,b_i)$ for $i=1,\dots ,M$. The set S thus contains, as before, all end points of the intervals $I_1,\dots ,I_M$ and [0, T] shifted by arbitrary integer multiples of s, as long as they are included in [0, T]. In addition, we now need to consider all fixing points $\tau _1,\ldots ,\tau _N$ and their corresponding shiftings. Clearly, we can compute S in $O((M+N)\sigma )$ time.

Lemma 4.3

Let v be a vertex of $C_{\overline{D(s)_ SP },\Pi }$. Then there exists $u\in \overline{D(s)_ SP }$ with $\Pi (u)=v$ such that u switches only in S.

Proof

The proof is similar to that of [6, Lemma 3.12], but one needs to pay attention to the fixings when shifting switching points outside of S. The full proof can be found in Appendix A.2. $\square$

We next show our main result that there exists an efficient separation algorithm for $C_{\overline{D(s)_ SP },\Pi }$ by specifying an efficient optimization algorithm over $C_{\overline{D(s)_ SP },\Pi }$. Let $\omega _1\dots ,\omega _{|S|}$ be the elements of S sorted in ascending order.

Theorem 4.4

One can optimize over $C_{\overline{D(s)_ SP },\Pi }$ (and hence also separate from $C_{\overline{D(s)_ SP },\Pi }$) in time polynomial in M, $\sigma$, and N.

Proof

By Lemma 4.3, it suffices to optimize over the projections of all $u\in \overline{D(s)_ SP }$ with switchings only in S. This can be done by a dynamic programming approach similar to the one presented in [6, Thm. 3.13]; we mainly need to change the recursion formula for the fixing points $\tau _1,\ldots ,\tau _N$. So assume that an arbitrary objective function $c\in {\mathbb {R}}^M$ is given. Then we compute the optimal value

$$\begin{aligned}c^*(t,b):=\min \;c^\top \Pi (u\cdot \chi _{[0,t]})\; \text { s.t. }u\in \overline{D(s)_ SP },\;u(t)=b\text { if }t<T\end{aligned}$$

recursively for all $t\in S$ as follows. Starting with $c^*(\omega _1,b)=0$ if $\tau _1\ne 0$ and

$$\begin{aligned} c^*(\omega _1,b)={\left\{ \begin{array}{ll} \infty , & \text { if } c_1=1 \text { and } b=0 \\ 0, & \text { otherwise} \end{array}\right. } \end{aligned}$$

otherwise, we obtain for $k\in \{2,\ldots ,|S|\}$ with $\omega _k\in S{\setminus }\{\tau _1,\ldots ,\tau _N\}$ that

where for $b\in \{0,1\}$ we define $\tau (b):=\{\tau _j: c_j=b, j=1,\ldots ,N\}$. For $k\in \{2,\ldots ,|S|\}$ with $\tau _j=\omega _k\in S$, $i\in \{1,\ldots ,N\}$, we get

and

The desired optimal value is $\min \{c^*(T,0),c^*(T,1)\}$ then, and a corresponding optimal solution can be derived if this value is finite. Otherwise, the problem is infeasible due to the fixings, i.e., the polytope $C_{\overline{D(s)_ SP },\Pi }$ is empty. $\square$

In the proof of Theorem 4.4, the recursion formula of the dynamic optimization approach over $C_{\overline{D(s)_ SP },\Pi }$ is the same for the fixing points $\tau _1,\ldots ,\tau _N$ as for the points in $S{\setminus } \{\tau _1,\ldots ,\tau _N\}$, as long as the fixings are respected. This is not surprising, since in this case we do not know whether the control is already constantly $c_j$ before or after $\tau _j$. However, if the fixing is not respected, then it is clear that the control has to be constantly $c_j$ on $[\tau _j-s,\tau _j)$ and one has to check whether this is possible taking the other fixings and the start value zero into account. In particular, if $\tau _1=0$ and $c_1=1$, all controls in $\overline{D(s)_ SP }$ have to respect the fixing due to the start value zero, so that we have $c^*(0,0)=\infty$ in this case.

5 Computation of primal and dual bounds

The main task in every branch-and-bound algorithm is the fast computation of primal and dual bounds. While primal bounds are often obtained by applying rather straightforward heuristics to the original problem (P), see Sect. 5.2, the computation of dual bounds is a more complex task, see Sect. 5.1.

5.1 Dual bounds

Our goal is to obtain strong dual bounds by solving the convexified subproblems (SPC); see Sect. 3. To this end, we can use the outer approximation algorithm developed in [7], since ${\overline{\operatorname {conv}}}(D_ SP )={\overline{\operatorname {conv}}}(\overline{D_ SP })$ as already noted above. This approach is applicable whenever we have a separation algorithm for ${\overline{\operatorname {conv}}}(\overline{D_ SP })$ at hand; see Sect. 4. Within the outer approximation algorithm, we thus need to repeatedly solve problems of the form

where $G:L^2(0,T)\rightarrow {\mathbb {R}}^k$ with $(Gu)_\ell =a_\ell ^\top \Pi _\ell (u)$ for $\ell =1,\ldots ,k$. The latter constraints represent the cutting planes for the sets $C_{\overline{D_ SP },\Pi }$ that have been generated so far.

As discussed in Sect. 3, our branching strategy will implicitly fix the control u on certain subintervals of the time horizon [0, T]; see Example 3.3 and Example 3.4. Let ${\mathcal {A}}$ be the union of all such fixed intervals and set ${\mathcal {I}}:=[0,T]\setminus {\mathcal {A}}$. Denote the restrictions to ${\mathcal {A}}$ and ${\mathcal {I}}$ by $\chi _{{\mathcal {A}}}:L^2(0,T) \rightarrow L^2({\mathcal {A}})$ and $\chi _{{\mathcal {I}}}:L^2(0,T) \rightarrow L^2({\mathcal {I}})$, respectively, and let $\chi _{\mathcal {A}}^*$ and $\chi _{{\mathcal {I}}}^*$ be the respective extension-by-zero operators mapping from $L^2({\mathcal {A}})$ and $L^2({\mathcal {I}})$ to $L^2(0,T)$, respectively. Then we can restrict ($\textrm{SPC }_k$) to the unfixed control $u|_{{\mathcal {I}}}$, which leads to

$$\begin{aligned} \left\{ \quad \begin{aligned} \text{ min }&J(S(\chi _{{\mathcal {I}}}^*u|_{ {\mathcal {I}}}+\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}}),\chi _{{\mathcal {I}}}^*u|_{ {\mathcal {I}}}+\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}})=:f(u|_{ {\mathcal {I}}}) \\ \text {s.t. }~&u|_{ {\mathcal {I}}} \in [0,1] \quad \text {a.e.~in } {\mathcal {I}}, \\&G(\chi _{{\mathcal {I}}}^*u|_{ {\mathcal {I}}}) \le b-G(\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}}), \end{aligned} \qquad \right. \end{aligned}$$

where $u|_{ {\mathcal {A}}}$ is fixed and implicitly given through the fixings. As a first attempt to solve this problem, we applied the semi-smooth Newton method described in [7], but, as the branching implicitly fixed larger parts of the switching structure, i.e., ${\mathcal {A}}$ got larger, the semi-smooth Newton method matrix became singular. To overcome these numerical issues, we decided to replace the semi-smooth Newton method by the alternating direction method of multipliers (ADMM), which was first mentioned in [21] for nonlinear elliptic problems and is widely applied to elliptic control problems [2, 4, 27]. Its convergence for convex optimization problems is well-studied, see, e.g., [15, 16, 20, 23]. Recently, [22] also addressed linear parabolic problems with state constraints by the ADMM and proved its convergence without any assumptions on the existence and regularity of the Lagrange multiplier.

We first need to rewrite our problem in the form

where ${\bar{b}}:= b-G(\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}})$ and

$$\begin{aligned} I_{(-\infty ,{\bar{b}}]}(v)={\left\{ \begin{array}{ll}0, & v\le {\bar{b}} \\ \infty , & \text{ otherwise }, \end{array}\right. }\quad \text{ and } \quad I_{[0,1]}(w)={\left\{ \begin{array}{ll}0, & w(t)\in [0,1] \text{ f.a.a. } t \\ \infty , & \text{ otherwise }. \end{array}\right. } \end{aligned}$$

Note that ($\textrm{SPC}{'_k}$) is still a convex optimization problem, but no longer strictly convex. The first-order algorithm ADMM is an alternating minimization scheme for computing a saddle point of the augmented Lagrangian

$$\begin{aligned} \begin{aligned} L_{\rho ,\beta }(u|_{{\mathcal {I}}},v,w,\lambda ,\mu ) =~&f(u|_{\mathcal {I}}) +I_{(-\infty ,{\bar{b}}]}(v)+I_{[0,1]}(w)\\&+ \lambda ^\top (G(\chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}})-v)+\langle \mu , u|_{\mathcal {I}}-w\rangle _{L^2({\mathcal {I}}),L^2({\mathcal {I}})}\\&+ \tfrac{\rho }{2}\Vert G(\chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}})-v\Vert ^2+\tfrac{\beta }{2}\Vert u|_{{\mathcal {I}}}-w\Vert ^2_{L^2({\mathcal {I}})}, \end{aligned} \end{aligned}$$

which differs from the Lagrangian by the penalty terms $\tfrac{\rho }{2}\Vert G(\chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}})-v\Vert ^2$ for the cutting planes and $\tfrac{\beta }{2}\Vert u|_{{\mathcal {I}}}-w\Vert ^2_{L^2({\mathcal {I}})}$ for the box constraints, but has the same saddle points as the Lagrangian [15]. First, the augmented Lagrangian is minimized with respect to the unfixed control variables

$$\begin{aligned} u|_{{\mathcal {I}}}=\arg \min _{u|_{{\mathcal {I}}}} L_{\rho ,\beta }(u|_{{\mathcal {I}}},v,w,\lambda ,\mu ), \end{aligned}$$

then with respect to v and w, i.e.,

$$\begin{aligned} \begin{aligned} v&=\arg \min _{v} L_{\rho ,\beta }(u|_{{\mathcal {I}}},v,w,\lambda ,\mu ),\\ w&=\arg \min _w L_{\rho ,\beta }(u|_{{\mathcal {I}}},v,w,\lambda ,\mu ), \end{aligned} \end{aligned}$$

and finally, the dual variables $\lambda$ and $\mu$ are updated by a gradient step as follows

$$\begin{aligned} \begin{aligned} \lambda&=\lambda + \gamma _\rho \rho \,\partial _\lambda L_{\rho ,\beta }(u|_{{\mathcal {I}}},v,w,\lambda ,\mu ), \\ \mu&=\mu +\gamma _\beta \beta \,\partial _\mu L_{\rho ,\beta }(u|_{{\mathcal {I}}},v,w,\lambda ,\mu ). \end{aligned} \end{aligned}$$

For $\gamma _\rho ,\gamma _\beta \in (0,\tfrac{1+\sqrt{5}}{2})$, the convergence of ADMM is guaranteed [19], but these parameters and the penalty parameters influence the convergence performance and numerical stability of the algorithm. For instance, the penalty parameter $\beta$ should be chosen close to $\alpha$ in order to balance the Tikhonov term $\tfrac{\alpha }{2}\Vert \chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}}+\chi _{{\mathcal {A}}}^*u|_{{\mathcal {A}}}-\tfrac{1}{2}\Vert _{L^2(0,T)}$ and the penalty term of the box constraints in the augmented Lagrangian. Moreover, the best choice for $\gamma _\rho$ and $\gamma _\beta$ generally seems to be one [19]. We thus use $\gamma _\rho =\gamma _\beta =1$ in the following.

With the solution mapping $S=\Sigma \circ \Psi +\zeta$, as defined in Sect. 2, the reduced objective in ($\textrm{SPC}{'_k}$) reads

$$\begin{aligned} f(u|_{{\mathcal {I}}})=\tfrac{1}{2}\, \Vert \Sigma \Psi (\chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}}+\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}}) + \zeta - y_{d }\Vert _{L^2(Q)}^2 + \tfrac{\alpha }{2} \,\Vert \chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}}+\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}} - \tfrac{1}{2}\Vert _{L^2(0,T)}^2 \end{aligned}$$

such that, by the chain rule, its Fréchet derivative at $u|_{ {\mathcal {I}}} \in L^2({\mathcal {I}})$ is given by

$$\begin{aligned} f'(u|_{ {\mathcal {I}}}) = \chi _{{\mathcal {I}}}\Psi ^* \Sigma ^* (\Sigma \Psi (\chi _{{\mathcal {I}}}^*u|_{{\mathcal {I}}}+\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}}) + \zeta - y_{d }) + \alpha \big (u|_{ {\mathcal {I}}} - \tfrac{1}{2}\big ) \in L^2({\mathcal {I}}), \end{aligned}$$

where we identified $L^2({\mathcal {I}})$ with its dual using the Riesz representation theorem. For the penalty term associated with the cutting planes, the Fréchet derivative at $u|_{ {\mathcal {I}}} \in L^2({\mathcal {I}})$ is

$$\begin{aligned} \rho \, \chi _{{\mathcal {I}}} G^*\big (G(\chi _{{\mathcal {I}}}^*u|_{ {\mathcal {I}}})-v\big ). \end{aligned}$$

With the above Fréchet derivatives at hand, we are able to write down the ADMM method for ($\textrm{SPC}{'_k}$). Algorithm 1 shows the procedure, where m is the iteration counter.

Algorithm 1

ADMM method for ($\textrm{SPC}{'_k}$)

The primal and dual residuals

$$\begin{aligned} r_P^{m}=\begin{pmatrix} G(\chi _{{\mathcal {I}}}^*u_{\ | {\mathcal {I}}}^{m})-v^{m} \\ u_{\ | {\mathcal {I}}}^{m}-w^{m} \end{pmatrix}, \quad r_D^{m}=\rho \, \chi _{\mathcal {I}}G^*(v^{m-1}-v^{m})+\beta (w^{m-1}-w^{m}) \end{aligned}$$

of the optimality conditions for ($\textrm{SPC}{'_k}$) can be used to bound the primal objective sub-optimality [5], i.e., $f(u_{\ |{\mathcal {I}}}^{m})-f(u^\star )$. More precisely, [5] derived sub-optimality estimates for problems in ${\mathbb {R}}^n$ based on their primal and dual residuals, but the arguments readily carry over to our setting. We thus have

$$\begin{aligned} f(u_{\ |{\mathcal {I}}}^{m})-f(u^\star ) \le - (r_P^{m})^\top \begin{pmatrix} \lambda ^{m} \\ \mu ^{m} \end{pmatrix} + (u_{\ |{\mathcal {I}}}^{m}-u_{\ |{\mathcal {I}}}^\star , r_D^{m})_{L^2({\mathcal {I}})} \end{aligned}$$

so that we can estimate

$$\begin{aligned} f(u_{\ |{\mathcal {I}}}^{m})-f(u^\star ) \le - (r_P^{m})^\top \begin{pmatrix} \lambda ^{m} \\ \mu ^{m} \end{pmatrix} + \sqrt{T}\, \Vert r_D^{m}\Vert _{L^2({\mathcal {I}})}=:e^{m}, \end{aligned}$$

(6)

since $u_{\ |{\mathcal {I}}}^{m},\ u_{\ |{\mathcal {I}}}^\star \in \{0,1\}$ a.e. in ${\mathcal {I}}\subset [0,T]$. As a reasonable stopping criterion, we choose that the primal and dual residual must be small, as well as the primal objective sub-optimality. As tolerances for the residuals, we may use an absolute and relative criterion, such as

$$\begin{aligned} \begin{aligned} \Vert r_P^m\Vert&\le (\sqrt{k}+1) \varepsilon ^\text {abs} + \varepsilon ^\text {rel} \max \{\Vert G(\chi _{{\mathcal {I}}}^*u_{\ | {\mathcal {I}}}^{m})\Vert _2+\Vert u^m\Vert _{L^2({\mathcal {I}})},\Vert v^m\Vert _2+\Vert w^m\Vert _{L^2({\mathcal {I}})}\}\, \\ \Vert r_D^m\Vert&\le \varepsilon ^\text {abs} + \varepsilon ^\text {rel} \Vert \chi _{\mathcal {I}}G^\star \lambda ^m+\mu ^m\Vert _{L^2({\mathcal {I}})}\, \end{aligned} \end{aligned}$$

where $\varepsilon ^\text {abs}>0$ is an absolute tolerance, whose scale depends on the scale of the variable values, and $\varepsilon ^\text {rel} > 0$ is a relative tolerance, which might be $\varepsilon ^\text {rel}=10^{-3}$ or $\varepsilon ^\text {rel}=10^{-4}$. The factor $\sqrt{k}$ accounts for the fact that ($\textrm{SPC}{'_k}$) contains k cutting plane constraints. In addition, the absolute error $e^{m}$ in the primal objective should be less than a chosen tolerance $\varepsilon ^{\text {pr}}>0$.

When the algorithm stops, we obtain $f(u_{\ |{\mathcal {I}}}^{m})-e^{m}$ as a dual bound for the subproblem (SP) of the branch-and-bound algorithm, and we can either proceed by calling the separation algorithm again, in order to generate another violated cutting plane, if possible, or by stopping the outer approximation algorithm. When proceeding with the cutting plane algorithm, one has to solve another parabolic optimal control problem of the form ($\textrm{SPC}_{k}$) with an additional cutting plane $k+1$ by Algorithm 1. The performance of the algorithm can be improved by choosing the prior solution $(u,v,\lambda ,w,\mu )$ as initialization in Step 1, and setting the auxiliary variable to $v_{k+1}=b-G(\chi _{{\mathcal {A}}}^*u|_{ {\mathcal {A}}})$ as well as the dual variable to $\lambda _{k+1}=0$ for the new cutting plane, since the latter is violated by u for sure.

5.2 Primal bounds

Another crucial ingredient in the branch-and-bound framework are primal heuristics, i.e., algorithms for computing good feasible solutions of the original problem (P), which yield tight primal bounds. It is common to call such primal heuristics in each subproblem, where the heuristic is often guided by the optimal solution for the convexified problem being solved in this subproblem for obtaining a dual bound. In our case, we can apply problem-specific rounding strategies from the literature to the solution of ($\textrm{SPC}{'_k}$) found by the ADMM method, e.g., the Dwell Time Sum-up Rounding and Dwell Time Next Force Rounding algorithms [42] for the case of a minimum time span between two switchings, and the Adaptive Maximum Dwell Rounding strategy [35] for the case of an upper bound on the total number of switchings.

Moreover, it is often possible to efficiently optimize a linear objective function over the set $C_{D,\Pi }$, as shown in [6]. We can benefit from this as follows. First, we define an appropriate objective function based on the solution u of ($\textrm{SPC}{'_k}$). Second, we can use the resulting minimizer $v^\star \in C_{D,\Pi }$ and construct a control $u'\in D$ with $\Pi (u')=v^\star$. For the first task, one can consider the distance of u to $\tfrac{1}{2}$ over the intervals $I_i$ defining the local averaging operators of the projection $\Pi$ and define the i-th objective coefficient as

$$\begin{aligned} \int _{I_i} (\tfrac{1}{2}-u) \,\text {d}t = \lambda (I_i)(\tfrac{1}{2}-\Pi (u)_i). \end{aligned}$$

(7)

The intuition in this definition is that a bigger objective coefficient, i.e., a smaller average value of u on $I_i$, will promote a smaller entry $v^\star _i$ in the minimizer $v^\star$, and vice versa. The minimizer $v^\star$ will thus agree with $\Pi (u)$ as much as possible while guaranteeing $v^\star \in C_{D,\Pi }$. In fact, if $C_{D,\Pi }$ is a 0/1-polytope, then the minimization problem

$$\begin{aligned} \min _{v\in C_{D,\Pi }} ~~ \sum _{i=1}^M \lambda (I_i)~|v_i-\Pi (u)_i| \end{aligned}$$

(8)

can be reformulated as a linear optimization problem over $C_{D,\Pi }$, which is equivalent to the one with the objective coefficients given in (7). Moreover, if the intervals $I_1,\dots ,I_M$ agree with the given discretization, the minimization problem (8) is equivalent to the CIA problem addressed in [26, 34], which tracks the average of the relaxed solution over the given temporal grid of the discretization while respecting the considered switching constraints.

Example 5.1

For $D(\sigma )$, the set $C_{D(\sigma ),\Pi }$ is a 0/1-polytope by [6, Thm. 3.8], and any linear objective function can be optimized in linear time over $C_{D(\sigma ),\Pi }$ [8]. The minimizer $v^\star$ can thus be guaranteed to be binary and it can be computed very efficiently, which even allows to choose as intervals $I_1,\ldots ,I_M$ for the projection exactly the ones given by the currently used discretization in time. In this case, the minimizer $v^\star$ solves the CIA problem over $D(\sigma )$ and it is trivial to find a control $u'$ with $\Pi (u')=v^\star$. Indeed, on each interval $I_i$, we can set $u'$ constantly to $v^\star _i$. $\square$

Example 5.2

The set $C_{D(s),\Pi }$ of the minimum dwell time constraints is not necessarily a 0/1-polytope, but one can optimize over $C_{D(s),\Pi }$ in $O(M\sigma )$ time, and, by backtracking, one can reconstruct the corresponding solution $u'\in D(s)$ in $O(M\sigma )$ time; see [6, Thm. 3.13]. $\square$

The implicit fixings of the control in a subproblem of the branch-and-bound algorithm can also be considered explicitly in the optimization over $C_{D,\Pi }$ by setting the corresponding objective coefficients in (7) to $\infty$ and $-\infty$, respectively. More precisely, one may use sufficiently large/small objective coefficients in this case.

In the above examples, a feasible control $u\in D$ can be computed quickly. Nevertheless, in order to obtain the corresponding primal bound, one needs to first calculate the resulting state $y=S(u)$ and then to evaluate the objective function.

6 Discretization error and adaptive refinement

The dual bounds computed by the outer approximation algorithm described in the previous section are safe bounds for ($\textrm{SPC}_{k}$), as long as we do not take discretization errors into account. However, our objective is to solve (P) in function space. This implies that we need to (a) estimate the discretization error contained in these bounds and (b) devise a method to deal with situations where the discretization-dependent dual bound allows to prune a subproblem but the discretization-independent dual bound does not, i.e., where the current primal bound lies between the two dual bounds. In the latter case, the only way out is the refinement of the discretization.

In order to address the first task, we will estimate the a posteriori error of the discretization with respect to the cost functional. We use the dual weighted residual (DWR) method, which has already achieved good results in practice, and combine the results from [29] and [39] to obtain an error analysis for the suproblems ($\textrm{SPC}_{k}$) arising in our branch-and-bound tree. First, we describe the finite element discretization of the optimal control problems arising in the branch-and-bound algorithm (Sect. 6.1). Then we discuss how to compute safe dual bounds (Sect. 6.2) as well as safe primal bounds (Sect. 6.3). Finally, we describe our adaptive refinement strategy (Sect. 6.4).

6.1 Finite element discretization

To solve problems of the form ($\textrm{SPC}_{k}$) in practice, we need to discretize the PDE constraint given as

$$\begin{aligned} & \langle \partial _t y, \varphi \rangle _{L^2(0,T;H^{-1}(\Omega )),L^2(0,T;H_0^1(\Omega ))}+(\nabla y,\nabla \varphi )_{L^2(0,T;L^2(\Omega ))}+(y(0),\varphi (0))_{L^2(\Omega )}\nonumber \\ & \qquad =(\Psi (u),\varphi )_{L^2(0,T;L^2(\Omega ))}+(y_0,\varphi (0))_{L^2(\Omega )}\qquad \forall \varphi \in W(0,T) \end{aligned}$$

(9)

in its weak formulation, as well as the control function, so that we implicitly discretize the Lagrangian $L:W(0,T)\times L^2(0,T)\times W(0,T) \times L^2(0,T)\times L^2(0,T)\times {\mathbb {R}}^k$ corresponding to ($\textrm{SPC}_{k}$) given as

$$\begin{aligned} \begin{aligned}L(y,u,p,\mu ^+,\mu ^-,\lambda )&= ~ J(y,u)-\langle \partial _t y, \varphi \rangle _{L^2(0,T;H^{-1}(\Omega )),L^2(0,T;H_0^1(\Omega ))}\\&-(\nabla y,\nabla p)_{L^2(0,T;L^2(\Omega ))}\\&-(y(0)-y_0,p(0))_{L^2(\Omega )}+(\Psi (u),p)_{L^2(0,T;L^2(\Omega ))} \\&+(\mu ^+,u-1)_{L^2(0,T)}- (\mu ^-,u)_{L^2(0,T)} + \lambda ^\top (Gu-b). \end{aligned} \end{aligned}$$

By calculating the derivative of L w.r.t. y in arbitrary direction $\varphi \in W(0,T)$, as well as applying interval-wise integration by parts to the equation, we get the adjoint equation

$$\begin{aligned} \begin{aligned} -&\langle y, \partial _t p\rangle _{L^2(0,T;H_0^1(\Omega )),L^2(0,T;H^{-1}(\Omega ))}{+}(\nabla \varphi ,\nabla p)_{L^2(0,T;L^2(\Omega ))}{+}(\varphi (T),p(T))_{L^2(\Omega )} \\&\quad =(\varphi ,y-y_d )_{L^2(0,T;L^2(\Omega ))} \qquad \forall \varphi \in W(0,T). \end{aligned}\nonumber \\ \end{aligned}$$

(10)

We use a discontinuous Galerkin element method for the time discretization of the PDE constraint with piecewise constant functions. Let

$$\begin{aligned} {\bar{J}}=\{0\}\cup J_1 \cup \cdots \cup J_{L-1}\cup J_L \end{aligned}$$

be a partition of [0, T] with half open subintervals $J_l=(t_{l-1},t_l]$ of size $s_l=t_l-t_{l-1}$ with time points $0=t_0<t_1<\cdots<t_{L-1}<t_L=T$. Define $s:=\max _{l=1,\ldots ,N} s_l$ as the maximal length of a subinterval. The spatial discretization of the state equation uses a standard Galerkin method with piecewise linear and continuous functions, where the domain $\Omega$ is partitioned into disjoint subsets $K_i$ of diameter $h_i:=\max _{p,q\in K_i}\Vert p-q\Vert _2$ for $i=1,\ldots , R$, i.e., $\overline{\Omega }=\cup _{i=1}^R \overline{K_i}$. For the one-dimensional domain $\Omega$ used in our experiments in Sect. 7, this means that we subdivide $\Omega$ into R disjoint intervals of length $h_i$. Set $h:=\max _{i=1,\ldots ,R} h_i$ and ${\mathcal {K}}_h=K_1\cup \cdots \cup K_R$. We define the finite element space

$$\begin{aligned} V_{h}:=\{v\in C({\bar{\Omega }})\cap H_0^1(\Omega ):v|_K\in P_1(K),\ K\in {\mathcal {K}}_h\} \end{aligned}$$

and associate with each time point $t_l$ a partition ${\mathcal {K}}_h^l$ of $\Omega$ and a corresponding finite element space $V_h^l\subset H_0^1(\Omega )$ which is used as spatial trial and test space in the time interval $J_l$. Denote by $P_0(J_l,V_h^l)$ the space of constant functions on $J_l$ with values in $V_h^l$. Then we use as a trial and test space for the state equation in (P) the space

$$\begin{aligned} X_{s,h}=\{v_{sh}\in L^2(I,L^2(\Omega )):v_{sh}|_{J_l}\in P_0(J_l,V_h^l),\ l=1,\ldots ,L\}. \end{aligned}$$

By introducing the notation

$$\begin{aligned} \begin{aligned}&y_{sh,l}^+=\lim \limits _{t\searrow 0}y_{sh}(t_l+t), \\&y_{sh,l}^-=\lim \limits _{t\searrow 0} y_{sh}(t_l-t)=y_{sh}(t_l),\; \text{ and } \\&[y_{sh}]_l:=y_{sh,l}^+ - y_{sh,l}^-\end{aligned} \end{aligned}$$

for the discontinuities of functions $y_{sh}\in X_{sh}$ in time, we obtain the following fully discretized state equation: Find for a given $u_{sh}\in L^2(0,T)$ a state $y_{sh} \in X_{s,h}$ such that

$$\begin{aligned} \begin{aligned}\sum _{l=1}^{L} \langle \partial _t y_{sh}, \varphi \rangle _{J_l}+\sum _{l=1}^{L}&(\nabla y_{sh}, \nabla \varphi )_{J_l}+\sum _{l=1}^{L-1} ([y_{sh}]_l,\varphi _l^+)+(y_{sh,0}^+,\varphi _0^+)\\ &=\sum _{l=1}^{L} (\Psi (u_{sh}),\varphi )_{J_l}+ (y_0,\varphi _0^+)\quad \forall \varphi \in X_{s,h},\end{aligned}\end{aligned}$$

(11)

where $\langle \cdot ,\cdot \rangle _{J_l}:=\langle \cdot ,\cdot \rangle _{L^2(J_l;H^{-1}(\Omega )),L^2(J_l;H_0^1(\Omega ))}$, $(\cdot ,\cdot )_{J_l}:=(\cdot ,\cdot )_{L^2(J_l;\Omega )}$, and $(\cdot ,\cdot ):=(\cdot ,\cdot )_{L^2(\Omega )}$. Note that, for piecewise constant states $y_{sh}\in X_{s,h}$, the term $\langle \partial _t y_{sh}, \varphi \rangle _{J_l}$ in (11) is zero for all $l=1,\ldots ,L$. We denote the discrete solution operator by $S_{sh}:L^2(0,T)\rightarrow X_{s,h}$, i.e., $y_{sh}=S_{sh}(u_{sh})$ satisfies the discrete state equation (11) for $u_{sh}\in L^2(0,T)$. Finally, we use piecewise constant functions for the temporal discretization of the control function on the same temporal grid as for the state equation, i.e., we use the space

$$\begin{aligned} Q_\rho =\{w\in L^2(0,T): w|_{J_l}=w_l \text{ for } \text{ all } l=1,\ldots ,L\}. \end{aligned}$$

Altogether, the discretization of ($\textrm{SPC}_{k}$) is given as

Moreover, the Lagrangian ${\tilde{L}}:X_{s,h}\times Q_\rho \times X_{s,h}\times Q_\rho \times Q_\rho \times {\mathbb {R}}^k \rightarrow {\mathbb {R}}$ associated with ($\textrm{SPC}{_{k\rho }}$) results as

$$\begin{aligned} \begin{aligned}&{\tilde{L}}(y_\rho ,u_\rho ,p_\rho , \mu _\rho ^+,\mu _\rho ^-,\lambda _\rho )\\&\quad =\,J(y_\rho ,u_\rho )-\sum _{l=1}^{L} \langle \partial _t y_\rho ,p_\rho \rangle _{J_l}-\sum _{l=1}^{L} (\nabla y_\rho ,\nabla p_\rho )_{J_l}\\&\qquad -\sum _{l=1}^{L-1} ([y_\rho ]_l,p_{\rho ,l}^+)-(y_{\rho ,0}^+-y_0,p_{\rho ,0}^+)+\sum _{l=1}^{L} (\Psi (u_\rho ),p_\rho )_{J_l} \\&\qquad + \sum _{l=1}^{L} \lambda (J_l)(\mu _\rho ^+|_{J_l})^\top (u_\rho |_{J_l}-1)- \sum _{l=1}^{L} \lambda (J_l)(\mu _\rho ^-|_{J_l})^\top u_\rho |_{J_l} + \lambda _\rho ^\top (Gu_\rho -b). \end{aligned} \end{aligned}$$

Based on this, we will devise a posteriori error estimates for both primal and dual bounds in the next sections.

6.2 A posteriori discretization error of dual bounds

Following the ideas of [29, 39], we now derive an a posteriori estimate for the error term $J(y,u)-J(y_\rho ,u_\rho )$, where $(y,u)\in W(0,T)\times L^2(0,T)$ denotes the optimizer of ($\textrm{SPC}_{k}$) and $(y_\rho ,u_\rho )\in X_{s,h}\times Q_\rho$ the one of ($\textrm{SPC}{_{k\rho }}$). For this, let us write down the first-order optimality conditions of ($\textrm{SPC}_{k}$) and ($\textrm{SPC}{_{k\rho }}$) by means of the Lagrangian L and ${\tilde{L}}$, respectively. If $(y,u)\in W(0,T)\times L^2(0,T)$ is optimal for ($\textrm{SPC}_{k}$), then there exist multipliers $p\in W(0,T)$, $\mu ^+\in L^2(0,T)$, $\mu ^-\in L^2(0,T)$ and $\lambda \in {\mathbb {R}}^k$ such that for $\chi :=(y,u,p,\mu ^+,\mu ^-,\lambda )$, we have

$$\begin{aligned} & L^\prime (\chi )(\delta y,\delta u,\delta p)=0 \quad \forall (\delta y,\delta u,\delta p) \in W(0,T)\times L^2(0,T)\times W(0,T) \end{aligned}$$

(12a)

$$\begin{aligned} & \quad \mu ^+\ge 0,\quad \mu ^+(u-1)=0,\quad u\le 1 \quad \text{ a.e. } \text{ in } (0,T) \end{aligned}$$

(12b)

$$\begin{aligned} & \quad \mu ^-\ge 0,\quad \quad \mu ^-u=0,\quad u\ge 0 \quad \text{ a.e. } \text{ in } (0,T) \end{aligned}$$

(12c)

$$\begin{aligned} & \quad \lambda \ge 0,\quad \lambda ^\top (Gu-b)=0, \quad Gu\le b \end{aligned}$$

(12d)

Analogously, if $(y_\rho ,u_\rho )\in X_{s,h}\times Q_\rho$ is optimal for ($\textrm{SPC}{_{k\rho }}$), then there exist $p_\rho \in X_{s,h}$, $\mu _\rho ^+\in Q_\rho$, $\mu _\rho ^-\in Q_\rho$ and $\lambda _\rho \in {\mathbb {R}}^k$ such that for $\chi _\rho :=(y_\rho ,u_\rho ,p_\rho ,\mu _\rho ^+,\mu _\rho ^-,\lambda _\rho )$ we have

$$\begin{aligned} & {\tilde{L}}^\prime (\chi _\rho )(\delta y,\delta u,\delta p)=0 \quad \forall (\delta y,\delta u,\delta p) \in X_{s,h}\times Q_\rho \times X_{s,h} \end{aligned}$$

(13a)

$$\begin{aligned} & \quad \mu _\rho ^+|_{J_l}\ge 0,\quad \mu _\rho ^+|_{J_l}(u_\rho |_{J_l}-1)=0,\quad u_\rho |_{J_l}\le 1 \quad \forall l=1,\ldots ,L \end{aligned}$$

(13b)

$$\begin{aligned} & \quad \mu _\rho ^-|_{J_l}\ge 0,\quad \quad \mu _\rho ^-|_{J_l}u_\rho |_{J_l}=0,\quad u_\rho |_{J_l} \ge 0 \quad \forall l=1,\ldots ,L \end{aligned}$$

(13c)

$$\begin{aligned} & \quad \lambda _\rho \ge 0,\quad \lambda _\rho ^\top (Gu_\rho -b)=0,\quad Gu_\rho \le b \end{aligned}$$

(13d)

Using the shorthand notation

$$\begin{aligned} \begin{aligned} {\mathcal {Y}}&=W(0,T)\times L^2(0,T)\times W(0,T)\times L^2(0,T) \times L^2(0,T) \times {\mathbb {R}}^k \text{ and } \\ {\mathcal {Y}}_\rho&= X_{s,h}\times Q_\rho \times X_{s,h}\times Q_\rho \times Q_\rho \times {\mathbb {R}}^k, \end{aligned} \end{aligned}$$

we have everything at hand to combine the results from [29] and [39] to obtain the following a posteriori discretization error estimation.

Theorem 6.1

Let $\chi =(y,u,p,\mu ^+,\mu ^-,\lambda )\in {\mathcal {Y}}$ fulfill the first-order necessary optimality conditions (12a)–(12d) for ($\textrm{SPC}_{k}$) and $\chi _\rho =(y_\rho ,u_\rho ,p_\rho ,\mu _\rho ^+,\mu _\rho ^-,\lambda _\rho )\in ~{\mathcal {Y}}_\rho$ the first-order necessary optimality conditions (13a)–(13d) for the discretized problem ($\textrm{SPC}{_{k\rho }}$). Then

$$\begin{aligned} J(y,u)-J(y_\rho ,u_\rho )&=\tfrac{1}{2}{\tilde{L}}^\prime (\chi )(\chi -\chi _\rho )+\tfrac{1}{2}{\tilde{L}}^\prime (\chi _\rho )(\chi -\chi _\rho ) \\&=\tfrac{1}{2}\Big (\, {\tilde{L}}_y^\prime (\chi _\rho )(y-y_\rho )+{\tilde{L}}_p^\prime (\chi _\rho )(p-p_\rho )+{\tilde{L}}_{u}^\prime (\chi _\rho )(u-u_\rho )\\&~\quad +{\tilde{L}}_{\mu ^+}^\prime (\chi )(\mu ^+-\mu _\rho ^+)+{\tilde{L}}_{\mu ^-}^\prime (\chi )(\mu ^--\mu _\rho ^-)+{\tilde{L}}_{\lambda }^\prime (\chi )(\lambda -\lambda _\rho )\\&~\quad +{\tilde{L}}_{\mu ^+}^\prime (\chi _\rho )(\mu ^+-\mu _\rho ^+)+{\tilde{L}}_{\mu ^-}^\prime (\chi _\rho )(\mu ^--\mu _\rho ^-)+{\tilde{L}}_{\lambda }^\prime (\chi _\rho )(\lambda -\lambda _\rho ) \,\Big )\;. \end{aligned}$$

Proof

The main arguments of the following proof are taken from the proofs of [29, Thm. 4.1] and [39, Thm. 4.2]. From the first-order optimality system (12a)–(12d) of $\chi \in {{\mathcal {Y}}}$ for ($\textrm{SPC}_{k}$) we obtain $J(y,u)=L(\chi )$. Analogously, the first-order conditions (13a)–(13d) of $\chi _\rho \in {\mathcal {Y}}_\rho$ for ( $\textrm{SPC}{_{k\rho }}$) lead to $J(y_\rho ,u_\rho )={\tilde{L}}(\chi _\rho )$. Moreover, it holds $L(\chi )={\tilde{L}}(\chi )$ since the continuous embedding $W(0,T)\hookrightarrow C([0,T];L^2(\Omega ))$ [41, Prop. 23.23] guarantees $y\in W(0,T)$ to be continuous in time such that the additional jump terms in ${\tilde{L}}$ compared to L vanish. We thus obtain

$$\begin{aligned} J(y,u)-J(y_\rho ,u_\rho )={\tilde{L}}(\chi )-{\tilde{L}}(\chi _\rho )=\int _{0}^1 {\tilde{L}}^\prime (\chi _\rho +s(\chi -\chi _\rho ))(\chi -\chi _\rho ) \,\text {d}s. \end{aligned}$$

Evaluation of the integral by the trapezoidal rule leads to

$$\begin{aligned} {\tilde{L}}(\chi )-{\tilde{L}}(\chi _\rho )=\tfrac{1}{2}{\tilde{L}}^\prime (\chi )(\chi -\chi _\rho )+\tfrac{1}{2}{\tilde{L}}^\prime (\chi _\rho )(\chi -\chi _\rho )+R \end{aligned}$$

(14)

with the residual

$$\begin{aligned} R= \tfrac{1}{2} \int _{0}^1 {\tilde{L}}^{\prime \prime \prime }(\chi + \zeta (\chi - \chi _\rho ))(\chi - \chi _\rho ,\chi - \chi _\rho ,\chi - \chi _\rho )\zeta (\zeta - 1)\,\text {d}\zeta . \end{aligned}$$

Since the PDE contained in ($\textrm{SPC}_{k}$) as well as the control constraints in u are linear, and the objective is quadratic in y and u, respectively, we have $R=0$.

We now have a closer look at the different error terms arising in (14). First, we have

$$\begin{aligned} {\tilde{L}}^\prime (\chi )(\chi -\chi _\rho )={\tilde{L}}_{\mu ^+}^\prime (\chi )(\mu ^+-\mu _\rho ^+)+{\tilde{L}}_{\mu ^-}^\prime (\chi )(\mu ^--\mu _\rho ^-)+{\tilde{L}}_{\lambda }^\prime (\chi )(\lambda -\lambda _\rho ), \end{aligned}$$

because the other terms are zero thanks to the condition (12a), which can be seen as follows: since $y\in W(0,T)$ is continuous in time due to $W(0,T)\hookrightarrow C([0,T];L^2(\Omega ))$ by [41, Prop. 23.23], the additional terms in ${\tilde{L}}_y^\prime$ compared to ${L}_y^\prime$ and ${\tilde{L}}_p^\prime$ compared to ${L}_p^\prime$, respectively, vanish, so that (12a) immediately yields ${\tilde{L}}_y^\prime (\chi )(y)=0$ and ${\tilde{L}}_p^\prime (\chi )(p)=0$. Moreover, the continuity of y in time implies that ${\tilde{L}}_p^\prime (\chi )(p_\rho )=0$ can equivalently be expressed as

$$\begin{aligned} \sum _{l=1}^{L} \langle \partial _t y, p_\rho \rangle _{J_l}+\sum _{l=1}^{L} (\nabla y, \nabla p_\rho )_{J_l}+(y_0^+,p_{\rho ,0}^+)=\sum _{l=1}^{L} (\Psi (u),p_\rho )_{J_l}+(y_0,p_{\rho ,0}^+). \end{aligned}$$

For the continuous state y, the state equation (9) implies that $(\varphi ,y(0))=(\varphi ,y_0)$ holds for all $\varphi \in L^2(\Omega )$, so that the term $(y_0^+,p_{\rho ,0}^+)$ containing $y(0)=y_0^+$ cancels out with $(y_0,p_{\rho ,0}^+)$, as $p_{\rho ,0}^+\in L^2(\Omega ),$ and it remains to ensure

$$\begin{aligned} \begin{aligned} \langle \partial _t y, p_\rho \rangle _{L^2(0,T;H^{-1}(\Omega )),L^2(J_l;H_0^1(\Omega ))}+(\nabla y,\nabla p_\rho )_{L^2(0,T;L^2(\Omega ))}=(\Psi (u),p_\rho )_{L^2(0,T;L^2(\Omega ))}. \end{aligned} \end{aligned}$$

Again from the continuous state equation (9), the latter equation is satisfied by y such that we obtain ${\tilde{L}}_p^\prime (\chi )(p_\rho )=0$, as desired. It remains to prove ${\tilde{L}}_y^\prime (\chi )(y_\rho )=0$. To this end, note that $p\in W(0,T)$ is continuous with respect to time by [41, Prop. 23.23], so that we can rewrite ${\tilde{L}}_y^\prime (\chi )(y_\rho )=0$ after interval-wise integration by parts in W(0, T) [17] as

$$\begin{aligned} \begin{aligned} -\sum _{l=1}^L \langle \partial _t p, y_\rho \rangle _{J_l}+\sum _{l=1}^L(\nabla y_\rho ,\nabla p)_{J_l}+(y_{\rho ,L}^-,p_{L}^-) =\sum _{l=1}^L (y_\rho ,y-y_d )_{J_l}. \end{aligned} \end{aligned}$$

Using $p_L^-=p(T)=0$ for the adjoint $p\in W(0,T)$, the above equation becomes

$$\begin{aligned} \begin{aligned} -\sum _{l=1}^L \langle \partial _t p,y_\rho \rangle _{J_l}+\sum _{l=1}^L(\nabla y_\rho ,\nabla p)_{J_l} =\sum _{l=1}^L (y_\rho ,y-y_d )_{J_l}. \end{aligned} \end{aligned}$$

By the adjoint equation (10) and the density of W(0, T) in $L^2(0,T;H_0^1(\Omega ))$, the equation is satisfied by $y_\rho \in L^2(0,T;H_0^1(\Omega ))$. We thus get ${\tilde{L}}_y^\prime (\chi )(y_\rho )=0$. Finally, (12a) directly yields ${\tilde{L}}_u^\prime (\chi )(u-u_\rho )=0$ because of $(u-u_\rho )\in L^2(0,T)$. The second term in (14) is given as

$$\begin{aligned} \begin{aligned} {\tilde{L}}^\prime (\chi _\rho )(\chi -\chi _\rho )=~&{\tilde{L}}_y^\prime (\chi _\rho )(y-y_\rho )+{\tilde{L}}_p^\prime (\chi _\rho )(p-p_\rho )+{\tilde{L}}_{u}^\prime (\chi _\rho )(u-u_\rho )\\ &+{\tilde{L}}_{\mu ^+}^\prime (\chi _\rho )(\mu ^+-\mu _\rho ^+)+{\tilde{L}}_{\mu ^-}^\prime (\chi _\rho )(\mu ^--\mu _\rho ^-)+{\tilde{L}}_{\lambda }^\prime (\chi _\rho )(\lambda -\lambda _\rho ), \end{aligned} \end{aligned}$$

which completes the proof. $\square$

We need to further specify the estimation of the a posteriori error given in Theorem 6.1, since it contains the unknown solution $\chi \in {\mathcal {Y}}$. A common approach in the context of the DWR method is to use higher-order approximations, which work satisfactorily in practice; see, e.g., [3]. Since our control function can only vary over time and the novelty of our approach lies primarily in the determination of the finitely many switching points, we assume for simplicity that there is no error caused by the spatial discretization of the state equation to keep the discussion concise. Thus, we only use a higher-order interpolation in time. For that, we introduce the piecewise linear interpolation operator $I_s^{(1)}$ in time and map the computed solutions to the approximations of the interpolation errors

$$\begin{aligned} y-y_\rho \approx I_s^{(1)} y_\rho -y_\rho \text{ and } p-p_\rho \approx I_s^{(1)} p_\rho -p_\rho . \end{aligned}$$

Then we obtain the approximations

$$\begin{aligned} \begin{aligned} {\tilde{L}}_y^\prime (\chi _\rho )(y-y_\rho )&\approx {\tilde{L}}_y^\prime (\chi _\rho )(I_s^{(1)} y_\rho -y_\rho ), \\ {\tilde{L}}_p^\prime (\chi _\rho )(p-p_\rho )&\approx {\tilde{L}}_p^\prime (\chi _\rho )(I_s^{(1)} p_\rho -p_\rho ). \end{aligned} \end{aligned}$$

Since the space of the Lagrange multiplier $\lambda$ of the cutting planes is finite-dimensional and thus not implicitly discretized by the discretization of the control space, we may choose $\lambda _\rho$ as higher-order interpolating and consequently neglect the error terms in $\lambda$, i.e.,

$$\begin{aligned} {\tilde{L}}_{\lambda }^\prime (\chi )(\lambda -\lambda _\rho ) +{\tilde{L}}_{\lambda }^\prime (\chi _\rho )(\lambda -\lambda _\rho ) \approx 0. \end{aligned}$$

Finally, as mentioned in [39], the control u typically does not possess sufficient smoothness, due to the box and cutting plane constraints. We thus suggest, as in [39], based on the gradient equation

$$\begin{aligned} L_u^\prime (\chi )=\alpha (u-\tfrac{1}{2})+\Psi ^\star p + \mu ^+-\mu ^- +G^\star \lambda =0 \end{aligned}$$

and the resulting projection formula

$$\begin{aligned} u=\min \{\max \{-\tfrac{1}{\alpha }(\Psi ^\star p+ G^\star \lambda )+\tfrac{1}{2},0\},1\}, \end{aligned}$$

the choice of

$$\begin{aligned} {\tilde{u}}=\min \{\max \{-\tfrac{1}{\alpha }(\Psi ^\star I_s^{(1)}p_\rho + G^\star \lambda _\rho ) +\tfrac{1}{2},0\},1\} \end{aligned}$$

and

$$\begin{aligned} {\tilde{\mu }}=-\alpha ({\tilde{u}}-\tfrac{1}{2})-\Psi ^\star I_s^{(1)}p_\rho -G^\star \lambda _\rho =:{\tilde{\mu }}^+-{\tilde{\mu }}^- \end{aligned}$$

with ${\tilde{\mu }}^+,{\tilde{\mu }}^-\ge 0$ a.e. on (0, T). The computable error estimate is thus given as

with ${\tilde{\chi }}=(I_s^{(1)}y_\rho ,{\tilde{u}},I_s^{(1)}p_\rho ,{\tilde{\mu }}^+,{\tilde{\mu }}^-,\lambda _\rho )$.

As in [29], one could split the error $J(y,u)-J(y_\rho ,u_\rho )$ into (a) the error caused by the semi-discretization of the state equation in time, (b) the error caused by the additional spatial discretization of the state equation, which we would consider as zero again, and (c) the error caused by the control space discretization. This would allow to choose different time grids for the state equation and the control space, where the former has to be at least as fine as the latter [29]. Since we are mostly interested in the combinatorial switching constraints, so that our focus is on the controls, we decided not to split the error and thus not to consider a finer temporal grid for the state.

As discussed in Sect. 3, the given fixings may determine parts of the switching pattern of u in ($\textrm{SPC}_{k}$). In this case, we need to calculate the a posteriori error ($\textrm{E}_{\eta }$) only on the unfixed control variables $u|_{\mathcal {I}}$, as well as on the Lagrange multipliers $\mu ^+,\mu ^-\in L^2({\mathcal {I}})$ corresponding to the box constraints, since we explicitly eliminated the fixed control variables from the problem ($\textrm{SPC}_{k}$). Then, it is clear that the terms ${\tilde{L}}_{u}^\prime (\chi _\rho )({\tilde{u}}-u_\rho )$, ${\tilde{L}}_{\mu ^+}^\prime ({\tilde{\chi }})({\tilde{\mu }}^+-\mu _\rho ^+)$, ${\tilde{L}}_{\mu ^-}^\prime ({\tilde{\chi }})({\tilde{\mu }}^--\mu _\rho ^-)$, ${\tilde{L}}_{\mu ^+}^\prime (\chi _\rho )({\tilde{\mu }}^+-\mu _\rho ^+)$, and ${\tilde{L}}_{\mu ^-}^\prime (\chi _\rho )({\tilde{\mu }}^--\mu _\rho ^-)$ in the error estimator ($\textrm{E}_{\eta }$) tend to zero for an increasing number of fixings satisfying the assumptions of Theorem 3.1, since the non-fixed part of the time horizon vanishes in this case. On the other hand, the error terms ${\tilde{L}}_y^\prime (\chi _\rho )(I_s^{(1)} y_\rho -y_\rho )$ and ${\tilde{L}}_p^\prime (\chi _\rho )(I_s^{(1)} p_\rho -p_\rho )$ reflect the error $J(Su_\rho ,u_\rho )-J(y_\rho ,u_\rho )$ in the cost functional caused by calculating the discretized state $y_\rho =S_{sh}(u_\rho )$ rather than $Su_\rho$. This error is also taken into account in the primal bounds throughout our branch-and-bound scheme; see Sect. 6.3 below.

In summary, in order to numerically compute a safe dual bound for the subproblem (SP), we first calculate a solution $u_\rho$ of the fully discretized problem ($\textrm{SPC}{_{k\rho }}$) with objective value $J(y_\rho ,u_\rho )$ by means of the ADMM method, as described in Sect. 5.1. Second, we use $J(y_\rho ,u_\rho )-e+\eta$ as a dual bound, where e denotes the absolute error in the primal objective caused by the ADMM algorithm, see (6), and $\eta$ the a posteriori error of the discretization of ($\textrm{SPC}_{k}$); compare ($\textrm{E}_{\eta }$).

6.3 A posteriori discretization error of primal bounds

Every feasible solution $u\in D$, e.g., obtained by applying primal heuristics as described in Sect. 5.2, leads to a primal bound J(Su, u) for the original problem (P). However, this bound is again subject to discretization errors. To estimate the latter, we first need to solve the fully discretized equation (11) to get a state $y_{sh}=S_{sh}(u)$ and then to estimate the a posteriori error $\nu :=J(Su,u)-J(S_{sh}u,u)$ in the cost functional. For the latter, we can again use the DWR method, which was originally invented to estimate the error in the cost function caused by the discretization of the state equation, see, e.g., [3]. We may directly apply [3, Prop. 2.4] to get the approximation

$$\begin{aligned} \begin{aligned} \nu \approx \ p_y(y_{sh},u,p_{sh})(p-p_{sh}):=&-\sum _{l=1}^L (\nabla y_{sh},p-p_{sh})_{J_l}-\sum _{l=1}^{L-1} ([y_{sh}]_l,p_l^+-p_{sh,l}^+) \\&-(y_{sh,0}^+-y_0,p_0^+-p_{sh,0}^+)+\sum _{l=1}^L (\Psi (u),p-p_{sh})_{J_l}, \end{aligned} \end{aligned}$$

with $\langle \partial _t y_{sh}, p_{sh}\rangle _{J_l}=0$ for $l=1,\ldots ,L$, where $p=S^*(y)$ and $p_{sh}=S_{sh}^*(y_{sh})$ denotes the adjoint corresponding to the state $y=S(u)$ and $y_{sh}=S_{sh}(u)$, respectively. Assuming again that there is no error caused by the spatial discretization, we may use the piecewise linear interpolation $I_s^{(1)}p_{sh}$ of $p_{sh}$ in time to obtain the computable a posteriori error

$$\begin{aligned} \nu \approx p_y(y_{sh},u,p_{sh})(I_s^{(1)}p_{sh}-p_{sh}). \end{aligned}$$

Then $J(S_{sh}u,u)+\nu$ is a safe primal bound for (P).

6.4 Adaptive refinement strategy

The central feature of our branch-and-bound algorithm is the approximate computation of an optimal solution for (P) in function space. In the limit, this solution does not depend on any predetermined discretization of the time horizon. However, in practice, we need to discretize our subproblems (SP) in order to numerically compute dual bounds, as described in Sect. 6.2. The main idea of our approach is to use a coarse temporal grid at the beginning, when the branchings have not yet determined a significant part of the switching structure, and then to refine the subintervals (only) if necessary.

More specifically, as long as the time-mesh dependent dual bound $J(y_\rho ,u_\rho )-e$ for (SP) is below the best known primal bound, we proceed with the given discretization. Otherwise, we cannot find a better solution for (SP) for the given discretization. We then must decide whether better solutions for (SP) may potentially exist when using a finer temporal grid. This is the case if and only if the time-mesh independent bound $J(y_\rho ,u_\rho )-e+\eta$ is still below the primal bound PB. We thus have to refine the grid whenever

$$\begin{aligned} J(y_\rho ,u_\rho )-e+\eta \le PB < J(y_\rho ,u_\rho )-e. \end{aligned}$$

If even $J(y_\rho ,u_\rho )-e+\eta$ exceeds the primal bound, we can prune the subproblem. Indeed, in this case we cannot find better solutions for the subproblem even in function space.

The adaptive refinement of the temporal grid is guided by the a posteriori error estimation of the discretization proposed in Sect. 6.2. The error estimator ($\textrm{E}_{\eta }$) can be easily split into its contribution on each subinterval $J_l$, i.e.,

$$\begin{aligned} \eta =\sum _{l=1}^L \eta _l, \end{aligned}$$

with the local error contributions $\eta _l$ on $J_l$ for $l=1,\ldots ,L.$ Note that this splitting is directly possible since we assumed that there is no error caused by the spatial discretization of the state equation, and thus no further localization on each spatial mesh is needed. A popular strategy for mesh adaptation is to order the subintervals according to the absolute values of their error indicators in descending order, i.e., to find a permutation $\varrho$ of $\{1,\ldots ,L\}$ such that $|\eta _{\varrho (1)}| \ge \cdots \ge |\eta _{\varrho (L)}|$, and then to refine the subintervals which make up a certain percentage $\gamma >0$ of the total absolute error, i.e., the subintervals $J_{\varrho (1)},\ldots ,J_{\varrho (L_\gamma )}$ with

$$\begin{aligned} L_\gamma :=\min \Big \{j\in \{1,\ldots ,L\}:\sum _{l=1}^j |\eta _{\varrho (l)}| > \gamma \sum _{l=1}^L |\eta _l|\Big \}. \end{aligned}$$

The resulting subproblem ($\textrm{SPC}{_{k\rho }}$) with respect to the refined discretization again has to be solved by Algorithm 1. As a reoptimization strategy, the values of the prior discretized solution $(u_\rho ,v_\rho ,\lambda _\rho ,w_\rho ,\mu _\rho )$ returned by Algorithm 1 can be used to initialize the variables in Step 1. More precisely, the values of $(u_\rho ,w_\rho ,\mu _\rho )$ can be duplicated according to the refinement of the subintervals and $(v_\rho ,\lambda _\rho )$ can be kept unchanged. In this way, we produce a primal feasible solution $(u_\rho ,v_\rho ,w_\rho )$ for the new subproblem ($\textrm{SPC}{_{k\rho }}$), but note that $(\lambda _\rho ,\mu _\rho )$ is not feasible for the corresponding dual problem.

7 Numerical experiments

We now report the results of an extensive numerical evaluation of our branch-and-bound algorithm presented in the previous sections. The overall branch-and-bound method has been implemented in C++, using the DUNE-library [36] for the discretization of the PDE. The source code can be downloaded at https://github.com/agruetering/dune-bnb. For all experiments, we discretize the problems as described in Sect. 6.1. This means that the spatial discretization uses a standard Galerkin method with continuous and piecewise linear functionals, while the temporal discretization for the control, the state, and the desired state $y_d$ uses piecewise constant functionals in time. The spatial integrals in the weak formulation of the state equation (9) and the adjoint equation (10), respectively, are approximated by a Gauss-Legendre rule with order 3. This means that all spatial integrals except for the one containing the form function $\varphi$ are calculated exactly. The discretized systems, arising from the discretization of the state and adjoint equation, are solved by a sequential conjugate gradient solver preconditioned with AMG smoothed by SSOR. All computations have been performed on a 64bit Linux system with an Intel Xeon E5-2640 CPU @ 2.5 GHz and 32 GB RAM.

7.1 Algorithmic framework

We start the branch-and-bound algorithm with an equidistant time grid with 20 nodes and, if necessary, we refine the subintervals that account for $\gamma =50\,\%$ of the total error; see Sect. 6.4. The choice of the time point $\tau$ for the branching is crucial for the practical performance of the algorithm, since the implicit restrictions on the controls are highly influenced by the branching points; see Example 3.3 and Example 3.4. Thus, the quality of the dual bounds of each node in the branch-and-bound tree strongly depends on the branching decisions. As already mentioned in Sect. 3, it is natural to take the last computed relaxed control of the outer approximation algorithm into account, which we know up to a discretization of (0, T); see Sect. 6.1. As a branching point, we choose the point of the time grid where the control has the highest deviation from 0/1, i.e., where the distance to $\{0,1\}$ multiplied by the length of the corresponding grid cell is maximal. This branching strategy corresponds to the choice of the variable with the most fractional value in finite-dimensional integer optimization. Finally, we use breadth-first search as an enumeration strategy since our computed primal bounds track the average of the relaxed solution over the given temporal grid of the discretization, i.e., solve the CIA problem over $D(\sigma )$; compare Example 5.1. In depth-first search, the shape of the computed relaxed controls for the subproblems hardly changed, so that our primal heuristic always produced the same feasible solution and good primal bounds were found late. As a result, many nodes had to be examined before pruning. This effect is avoided by breadth-first search.

The results presented in [7] suggest to add only a few cutting planes before resorting to branching, because a significant increase in the dual bound was mostly obtained in the first cutting plane iterations. Moreover, we observed that the dual bounds got better with a decreasing Tikhonov parameter $\alpha$, but the time needed to compute them increased with decreasing $\alpha$. Thus, we investigate in Sect. 7.3 whether a good quality or a quick computation of the dual bounds have a greater influence on the overall performance.

The parabolic optimal control problems arising in each iteration of the outer approximation algorithm are solved by the ADMM algorithm; see Algorithm 1 in Sect. 5.1. As tolerances for the primal and dual residuals in the ADMM algorithm, we chose $\varepsilon ^\text {rel}=\varepsilon ^\text {abs}=10^{-3}$ and required the absolute error of the discretization of ($\textrm{SPC}{'_k}$) to be less than $\varepsilon ^\text {pr}=10^{-5}$. In order to guarantee the numerical stability of the ADMM algorithm, the penalty parameter of the cutting planes was set to $\rho =\tfrac{1+\sqrt{5}}{2}$. The best choice of the penalty term $\beta$ of the box constraints depending on the Tikhonov term $\alpha$ is investigated in Sect. 7.3. The resulting linear system in Step 3 of Algorithm 1 is solved by the conjugate gradient method, preconditioned with $P=(\alpha +\beta ) I + \rho G^\star G.$

7.2 Instances

In all experiments, we focus on the case of an upper bound $\sigma$ on the number of switchings, i.e., we consider the feasible set

$$\begin{aligned} D(\sigma ) = \big \{ u \in BV(0,T):\; u(t) \in \{0,1\} \text { f.a.a.\ } t \in (0,T),\; |u|_{BV(0,T)} \le \sigma \big \} \end{aligned}$$

as defined in Sect. 2.3. However, we assume that u is fixed to zero before the time horizon, so that we already count it as one switching if u is 1 at the beginning. Notwithstanding this slight modification, the most violated cutting plane for a given vector $v\notin C_{D(\sigma )_ SP ,\Pi }$ can be computed in $O(M+N)$ time as discussed in Sect. 4.1, using the separation algorithm presented in [8]. This separation algorithm is thus fast enough to allow to choose the intervals for the projection exactly as the intervals given by the discretization in time; compare Sect. 6.1.

We created instances of (P) with $\Omega =(0,1)$, $T=1\textrm{, and }\psi (x)=\exp (x)\sin (\pi \,x)+0.5.$ In order to obtain challenging instances, we produced the desired state $y_d$ as follows: we first generated a control $u_d :[0,T] \rightarrow \{0,1\}$ with a total variation $|u_d |_{BV(0,T)}=\theta$ and chose the desired state $y_d \textrm{as }S(u_d )$, such that $u_d$ is the optimal solution for Problem (P) if we allow $\theta$ switchings. More specifically, we randomly chose $\theta$ jump points $0< t_1< \cdots<t_{\theta } < T$ on the equidistant time grid with 320 nodes. Then, we chose $u_d :[0,T]\rightarrow \{0,1\}$ as the binary control starting in zero and having the switching points $t_1,\ldots ,t_{\theta }$. In this way, we generated non-trivial instances, where the constraint $D(\sigma )$ strongly affects the optimal solution of (P) in case $\sigma \ll \theta$.

7.3 Parameter tuning

Before testing the potential of our approach, we investigate the influence of some parameters on the overall performance. We first consider the Tikhonov term $\alpha$ and the penalty term $\beta$ of the box constraints; see Sect. 5.1. Afterwards, we investigate how time-consuming it is to solve the subproblems arising in the branch-and-bound algorithm, depending on when we stop the outer approximation algorithm for each subproblem (SP). Here, we resort to branching if the relative change of the bound is less than a certain percentage (RED) in three successive iterations. Finally, we vary the allowed relative deviation (TOL) of the objective value of the returned solution from the optimal value of (P); a subproblem in the branch-and-bound node is pruned when the remaining gap between primal and dual bound falls below this relative threshold. We start with RED $=$ TOL $=$ 1 %.

For all results presented in this subsection, we have chosen the same instance with $\theta =8$ jump points and allowed $\sigma =3$ switchings, since we observed the typical behavior of the algorithm with these settings. We always report the overall number of investigated subproblems (Subs), of cutting plane iterations (Cuts), and of ADMM iterations (ADMM). Moreover, the average number of fixings ($\varnothing$ FixPoints) and the average percentage of control variables that are implicitly fixed ($\varnothing$ FixIndices) are reported, where both averages are taken over all pruned subproblems. We also provide the overall run time (Time) in CPU hours.

Table 1 Influence of the Tikhonov parameter $\alpha$ and the penalty term $\beta$ of the box constraints on the branch-and-bound algorithm

The results for different values of $\alpha$ and $\beta$ can be found in Table 1. The main message of Table 1 is that a small value of $\alpha$ is generally favorable for the branch-and-bound algorithm, since a smaller value of $\alpha$ leads to stronger dual bounds and consequently, fewer fixings are needed on average to prune a subproblem. So, as long as no numerical issues arise with the ADMM algorithm and the DWR error estimator, one should choose $\alpha =0.001$. But, with smaller value of $\alpha$ it becomes more likely that the higher-order approximation of the unknown quantities (see Sect. 6.3) is too imprecise to estimate the error in the cost functional, so that the branch-and-bound algorithm returns wrong solutions. This was also observed in our experiments: in many instances, the obtained solutions for $\alpha \in \{0.01,0.005\}$ switched three times and had very similar switching times for all values of $\beta$. In contrast, the obtained solutions for $\alpha =0.001$ frequently switched only twice and differed enormously from the others. By recalculating the objective on such a fine grid that all returned solutions are piecewise constant on it, it turned out that the solutions obtained for $\alpha \in \{0.01,0.005\}$ were indeed better than the ones for $\alpha =0.001$. Moreover, the primal heuristic even produced some of the better solutions within the branch-and-bound scheme for $\alpha =0.001$, but due to the DWR error estimator, their time-mesh independent objective values were worse. For that reason, we choose $\alpha =\beta =0.005$ in all subsequent experiments.

Table 2 Impact of the ratio between branching and cutting plane iterations on the branch-and-bound algorithm

We next investigate the interplay between branching and outer approximation. Table 2 demonstrates that a good balance is important: a stronger focus on the outer approximation leads to fewer branching decisions needed to cut off a subproblem. However, this does not necessarily imply that fewer fixings are needed to prune a subproblem, since the branching points strongly depend on the shape of the relaxed solutions. Moreover, it is more time-consuming to solve each node due to the increased number of cutting plane iterations. On the other hand, it is also not beneficial to resort to branching too early because more subproblems need to be investigated then. We thus use RED  =  2% in the following.

Finally, the impact of the relative allowed deviation from the optimal objective value on the performance of the branch-and-bound algorithm is shown in Table 3. As expected, a higher tolerance leads to an earlier pruning of the subproblems, as indicated by the number of fixings required to prune a subproblem, and the running time decreases significantly. At the same time, however, the best known primal bound (Obj) found by the algorithm obviously increases, so that ultimately the user has to decide which deviation is still acceptable. We choose TOL  =  2% in the following, which we think is a reasonable optimality tolerance.

Table 3 Influence of the relative allowed deviation (TOL) from the optimum on the branch-and-bound algorithm

Fig. 2

Complete branch-and-bound tree of an instance generated with $\theta =3$ jump points and with $\sigma =1$ allowed switchings. The path of the optimal solution is marked in bold and the branching decisions along the optimal path are listed. In the case of a single child node, the temporal discretization of the subproblem has been refined

7.4 Performance of the algorithm

Before reporting running times and other key performance indicators of our algorithm, we first illustrate the interplay between branching and adaptive refinement by an example. Figure 2 shows the complete branch-and-bound tree obtained for an instance with $\theta =3$ jump points and only one allowed switching, i.e., $\sigma =1$. Whenever a node has a single child node in the illustration, the discretization of the subproblem has been refined. The branch-and-bound tree shows that a large part of the generated subproblems can already be pruned without any refinement. Moreover, in relatively few branches the subproblems need to be refined multiple times in order to decide whether a solution of desired quality can be found in these branches. The branching decisions taken along the path leading to the returned solution illustrate that, e.g., the generated subproblem 16 was refined in order to choose the sixth fixing point as $\tau _6=0.225$. This was not possible with the previous discretization of the problem. In particular, the fourth and fifth fixing point together have limited the switching point to be in the interval (0.2, 0.25]. The last branching decision in this tree serves to determine $t=0.2375$ as the switching point of the returned solution.

Table 4 Performance of the branch-and-bound algorithm for instances generated with $\theta$ switching points, allowing $\sigma$ switchings

Table 4 shows the performance of the branch-and-bound algorithm for various instances generated with $\theta \in \{1,\ldots ,8\}$ and $\sigma \in \{1,\ldots ,4\}$ for the total number of switching points. We were able to solve problems with up to four allowed switchings, but, as could be expected, the number of generated subproblems strongly increases in $\sigma$. However, we note that the ratio between generated subproblems and total cutting plane iterations is not affected by $\sigma$. While the branch-and-bound algorithm is able to solve problems with $\sigma =3$ within 14 CPU hours, the algorithm does not terminate within 60 CPU hours for most instances with $\sigma =4$ allowed switchings. However, the results of Table 4 show that the average number of subproblems in the branch-and-bound-tree remains relatively small for all instances, showing that the dual bounds computed by our algorithm are rather tight, and that the main challenge in terms of running times is the fast computation of these dual bounds.

Moreover, the reported results show that our approach to globally solve parabolic optimal control problems with dynamic switches by means of branch-and-bound, combined with an adaptive refinement strategy, works in practice. Whenever the maximal number of refinements of a grid cell within the branch-and-bound algorithm was larger than 4 in our experiments, a grid cell was refined this often in less than $10\,\%$ of the subproblems. Here, the finest grid mesh size decreases with the number of allowed switching points. This means that, if more switchings are allowed, a finer temporal discretization is needed to detect the optimal positions of the switching points.

In summary, our proposed branch-and-bound method is an effective and robust algorithm to globally solve control problems of the form (P). A few pointwise fixings of the controls suffice to significantly truncate the set of feasible switching patterns. Moreover, thanks to the computation of tight dual bounds by means of outer approximation, relatively few subproblems need to be inspected and refined within the branch-and-bound algorithm.

Finite-dimensional convex hulls under fixings

1.1 A.1 Restricted total variation

We show that if we restrict the total variation of a single switch to be less than $\sigma >0$, i.e., the set of feasible switching patterns is given by

$$\begin{aligned} D(\sigma )=\{u\in BV(0,T):u(t)\in \{0,1\} \text { a.e. in }(0,T),\ |u|_{BV(0,T)}\le \sigma \}, \end{aligned}$$

then the convex hull $C_{\overline{D(\sigma )_ SP },\Pi }$ of the finite dimensional projections $\{\Pi (u):u\in D(\sigma )_ SP \}$ under arbitrary fixings is a 0/1 polytope.

Theorem A.1

The set $C_{\overline{D(\sigma )_ SP },\Pi }$ is a 0/1 polytope.

Proof

The proof is similar to the one of Theorem 3.8 in [6], where no fixings have been considered. We claim that $C_{\overline{D(\sigma )_ SP },\Pi }=\operatorname {conv}(K)$, where

$$\begin{aligned} \begin{aligned} K:=\{\Pi (u) :&u \in \overline{D(\sigma )_ SP } \text { and for all }i=1,\ldots ,M \text { there exists } w_i\in \{0,1\}\\&\text {with } u(t)\equiv w_i \text { f.a.a.\ } t\in I_i\}. \end{aligned} \end{aligned}$$

From this, the result follows directly, as $K\subseteq \{0,1\}^M$ holds by definition.

Since K is a subset of $\{\Pi (u) :u \in \overline{D(\sigma )_ SP }\}$, the direction “$\supseteq$” is trivial. It thus remains to show “$\subseteq$”. For this, let $u\in \overline{D(\sigma )_ SP }$. We need to prove that $\Pi (u)$ can be written as a convex combination of vectors in K. Let $l\in \{0,\dots ,M\}$ denote the number of intervals in which the switch u is switched at least once. We prove the assertion by means of complete induction over the number l. For $l=0$, we clearly have $\Pi (u)\in K\subseteq \operatorname {conv}(K)$. So let $l>0$ and choose an index $\ell \in \{1,\ldots ,M\}$ such that u switches at least once in $I_\ell$. For $k=0,1$, define the function $u_k$ as follows:

$$\begin{aligned}u_k(t):={\left\{ \begin{array}{ll} k,& \text {if }t\in I_\ell \\ u(t), & \text {otherwise.} \end{array}\right. } \end{aligned}$$

Then, by construction, $\Pi (u)_\ell = \lambda \Pi (u_1)_\ell + (1-\lambda )\Pi (u_0)_\ell$ for $\lambda :=\Pi (u)_\ell \in [0,1]$ and $u_k$ has at most as many switching as u.

We next show that the controls $u_0$ and $u_1$ belong to $\overline{D(\sigma )_ SP }$. So let $k\in \{0,1\}$ be arbitrary. Due to $u\in \overline{D(\sigma )_ SP }$, there exists a sequence $\{v^m\}_{m\in {\mathbb {N}}}\in D(\sigma )_ SP$ such that $v^m\rightarrow u$ in $L^p(0,T)$ for $m\rightarrow \infty$. In particular, there exists a subsequence, which we denote by the same symbol for simplicity, with $v^m(t)\rightarrow u(t)$ f.a.a. $t\in (0,T)$ for $m\rightarrow \infty$. Since u switches at least once in the interval $I_\ell$ and $v^m$ converges pointwise almost everywhere to u, there exists $m_0\in {\mathbb {N}}$ such that for all $m\ge m_0$ the controls $v^m$ also switch at least once in $I_\ell$. When constructing a sequence in $D(\sigma )_ SP$ converging to $u_k$ with the help of $\{v^m\}_{m\in {\mathbb {N}}}$, we need to consider that fixing points $\tau _j$ may coincide with the interval limits of $I_\ell =(a_\ell ,b_\ell )$ so that we are only able to change the values in the inner of $I_\ell$. Thus, we define

$$\begin{aligned} w_k^m(t)={\left\{ \begin{array}{ll} k, & t\in [a_\ell +\tfrac{\lambda (I_\ell )}{2m}, b_\ell -\tfrac{\lambda (I_\ell )}{2m}) \\ v^m(a_\ell ), & t\in [a_\ell ,a_\ell +\tfrac{\lambda (I_\ell )}{2m}) \\ v^m(b_\ell ), & t\in [b_\ell -\tfrac{\lambda (I_\ell )}{2m},b_\ell ) \\ v^m(t), & \text {otherwise}. \end{array}\right. } \end{aligned}$$

Due to $v^m\in \{0,1\}$ a.e. in (0, T), also $w_k^m(t)\in \{0,1\}$ holds f.a.a. $t\in (0,T)$. By our general assumption, we have $\tau _j\notin (a_\ell ,b_\ell )$ for all $j=1,\ldots ,N$, so that $w_k^m(\tau _j)=v^m(\tau _j)=c_j$ follows with $v^m\in D(\sigma )_ SP$. Furthermore, for $m\ge m_0$, $w_k^m$ has at most as many switchings as $v^m$ in total and we thus obtain $w_k^m\in D(\sigma )_ SP$ for $m\ge m_0$. It is easy to see that $w_k^m \rightarrow u_k$ in $L^p(0,T)$ for $m\rightarrow \infty$, so that we get $u_k\in \overline{D(\sigma )_ SP }$, as claimed.

By the induction hypothesis, the vectors $\Pi (u_k)$ can thus be written as convex combinations of vectors in K and consequently, also $\Pi (u)$ is a convex combination of vectors in K. $\square$

1.2 A.2 Switching point constraints

In the following, we show that the finite-dimensional convex hulls of switching point constraints

$$\begin{aligned} \begin{aligned} D(P):= \{ u_{t_1,\dots ,t_\sigma } :(t_1, \ldots , t_\sigma ) \in P,~ 0\le t_1\le \cdots \le t_\sigma < \infty \}, \end{aligned} \end{aligned}$$

where $P\subseteq {\mathbb {R}}_+^\sigma$ is a given polytope, under fixings are still polytopes. Moreover, we prove some auxiliary results for Theorem 4.4 in Sect. 4.2, stating that the corresponding separation problems are tractable in the case of $P=\{t\in {\mathbb {R}}_+^\sigma : t_i-t_{i-1}\ge s \ \forall i=1,\ldots ,\sigma \}$ for some $s>0$.

For this purpose, let $-1=z_0< z_1< \cdots z_{r-1} < z_r=\infty$ include all end points of the intervals $I_1,\ldots ,I_M$ defining $\Pi$ and the fixing points $\tau _j$, $j=1,\ldots ,N$. In addition, let ${\mathcal {Z}}$ be the set of all maps $\varphi :\{1,\ldots ,\sigma \}\rightarrow \{1,\ldots ,r\}$ which assign each switching point to an interval $[z_{\varphi (i)-1},z_{\varphi (i)}]$ such that an even number of switching points is assigned to $(\tau _{j-1},\tau _j]$ if $c_{j-1}=c_j$ and an odd number otherwise, for all $j=1,\ldots ,N$, where $\tau _0:=-1$ and $c_0:=0$. Define

$$\begin{aligned} \begin{aligned}Q_\varphi :=\{(t_1,\ldots ,t_\sigma )\in P:\ &t_1\le \cdots \le t_\sigma , \ z_{\varphi (i)-1}\le t_i \le z_{\varphi (i)}\;\forall i=1,\dots ,\sigma ,\;\\ &z_{\varphi (i)-1}<t_i \; \forall i \in J\}\end{aligned} \end{aligned}$$

for all $\varphi \in {\mathcal {Z}}$, where $J_\varphi := \{i\in \{1,\ldots ,\sigma \}:\exists j\in \{1,\ldots ,N\} \text { s.t. }z_{\varphi (i)-1}=\tau _j \}$ as above, and

$$\begin{aligned}V_\varphi :=\{u_{t_1,\ldots ,t_\sigma }\mid (t_1,\ldots ,t_\sigma )\in Q_\varphi \}.\end{aligned}$$

Then the following holds true.

Lemma 1

$$\begin{aligned} D(P)_ SP =\bigcup _{\varphi \in {\mathcal {Z}}} V_\varphi . \end{aligned}$$

Proof

Let $u=u_{t_1,\ldots ,t_\sigma }\in D(P)_ SP$ with $(t_1,\ldots ,t_\sigma )\in P$. Define ${\bar{\varphi }}: \{1,\ldots ,\sigma \}\rightarrow \{1,\ldots ,r\}$ such that $z_{{\bar{\varphi }}(i)-1}<t_i\le z_{{\bar{\varphi }}(i)}$ holds for $i=1,\ldots ,\sigma$. Due to $u_{t_1,\ldots ,t_\sigma }(\tau _1)=c_1$, the other fixings $u_{t_1,\ldots ,t_\sigma }(\tau _j)=c_j$, $2\le j\le N$, can only be satisfied if the number of switching points in $(\tau _{j-1},\tau _{j}]$ is even in the case $c_{j-1}=c_{j}$ and odd, otherwise. If $c_1=0$, then $u_{t_1,\ldots ,t_\sigma }(\tau _1)=0$ only holds if an even number of switching points is less or equal to $\tau _1$, and in the other case $c_1=1$, this number must be odd. Consequently, we obtain ${\bar{\varphi }}\in {\mathcal {Z}}$ and $u\in V_{{\bar{\varphi }}}$.

For the reverse inclusion, let $u\in V_\varphi$ for some $\varphi \in {\mathcal {Z}}$. Then there exists $(t_1,\ldots ,t_\sigma )\in Q_\varphi$ such that $u=u_{t_1,\ldots ,t_\sigma }$. With $Q_\varphi \subseteq P$ it follows that $u\in D(P)$. Since $\varphi \in {\mathcal {Z}}$, we know that the correct number of switching points is assigned between $\tau _{j-1}$ and $\tau _{j}$ in order to respect the given fixings in $D(P)_ SP$. Moreover, the last requirement in the definition of $Q_\varphi$ ensures that no switching point assigned to the right neighboring interval of $\tau _j$ is equal to $\tau _j$, so the given fixings $u_{t_1,\ldots ,t_\sigma }(\tau _j)=c_j$ are indeed satisfied for all $j\in \{1,\ldots ,N\}$, which completes the proof. $\square$

To show that the convex hull of all projection vectors from controls $u\in \overline{D(P)_ SP }$ is a polytope, we can use that $\overline{D(P)_ SP }=\bigcup _{\varphi \in {\mathcal {Z}}} \overline{V_\varphi }$ holds, thanks to Lemma 1 and the fact that ${\mathcal {Z}}$ is finite. Consequently, we essentially need that $\Pi (\overline{V_\varphi })$ is a polytope for every $\varphi$ to deduce the polyhedricity of $C_{\overline{D(P)_ SP },\Pi }$. For this, we now prove that we simply need to consider the closure of the sets $Q_\varphi$ in ${\mathbb {R}}^\sigma$ to obtain $\overline{V_\varphi }$, with the help of the following auxiliary result:

Lemma 2

If there exists a sequence $\{u_{t_1^m,\ldots ,t_\sigma ^m}\}_{m\in {\mathbb {N}}}$ with $u_{t_1^m,\ldots ,t_\sigma ^m}\rightarrow u$ in $L^p(0,T)$ for some $u\in L^p(0,T)$ and $t^m:=(t_1^m,\ldots ,t_\sigma ^m) \rightarrow {\bar{t}}$ in ${\mathbb {R}}^\sigma$, then $u=u_{{{\bar{t}}}_1,\dots ,{{\bar{t}}}_\sigma }$.

Proof

The assertion is proven in [6, Lemma 3.10] and is based on the continuity of the mapping ${\mathbb {R}}^\sigma \ni (t_1,\dots ,t_\sigma )\mapsto u_{t_1,\ldots ,t_\sigma }\in L^p(0,T)$. $\square$

Lemma 4

$$\begin{aligned} \overline{V_\varphi }= \{u_{t_1,\ldots ,t_\sigma }:(t_1,\ldots ,t_\sigma )\in \overline{Q_\varphi }\} \end{aligned}$$

Proof

First, let $u\in \overline{V_\varphi }$ and consider a sequence $\{u^m\}_{m\in {\mathbb {N}}}$ in $V_\varphi$ with $u^m =u_{t^m_1,\dots ,t^m_\sigma }\rightarrow u$ in $L^p(0,T)$, where $t^m=(t^m_1,\dots ,t^m_\sigma )\in Q_\varphi$. The strong convergence in $L^p$ implies that there is a subsequence, denoted by the same symbol for convenience, which converges pointwise almost everywhere in (0, T) to u. Furthermore, as a polytope, P is bounded by definition, so that $Q_\varphi$ is bounded as well and thus there is yet another subsequence such that $t^m$ converges to ${{\bar{t}}} \in \overline{Q_\varphi }$. With Lemma 2, we may conclude that $u_{{{\bar{t}}}_1,\dots ,{{\bar{t}}}_\sigma }=u$ and, thanks to ${{\bar{t}}} \in \overline{Q_\varphi }$, this finishes the proof of the first inclusion.

For the reverse inclusion, consider $u_{t_1,\ldots ,t_\sigma }$ with switching points $t=(t_1,\ldots ,t_\sigma )\in \overline{Q_\varphi }$. Since $t\in \overline{Q_\varphi }$, there exists a sequence $t^m=(t^m_1,\dots ,t^m_\sigma )\in Q_\varphi$ with $t^m\rightarrow t$ in ${\mathbb {R}}^\sigma$. Again thanks to the continuity of the mapping $(t_1,\dots ,t_\sigma )\mapsto u_{t_1,\ldots ,t_\sigma }$ [6, Lemma 3.10], the sequence $\{u_{t^m_1,\ldots ,t^m_\sigma }\}_{m\in {\mathbb {N}}}\subseteq V_\varphi$ converges to $u_{t_1,\ldots ,t_\sigma }$ in $L^p(0,T)$, so that the latter belongs to the closure of $V_\varphi$ in $L^p(0,T)$. $\square$

We have everything at hand to prove the polyhedricity of $C_{\overline{D(P)_ SP },\Pi }$.

Theorem A.5

The set $C_{\overline{D(P)_ SP },\Pi }$ is a polytope.

Proof

We have $\Pi (\overline{D(P)_ SP })=\bigcup _{\varphi \in {\mathcal {Z}}} \Pi (\overline{V_{\varphi }})$ due to Lemma 1 and the fact that ${\mathcal {Z}}$ is finite. Since P is a polytope, also $\overline{Q_\varphi }$ is a polytope. Moreover, analogously to Theorem 3.11 in [6], one obtains that the function $\overline{Q_\varphi }\ni (t_1,\dots ,t_\sigma )\mapsto \Pi (u_{t_1,\dots ,t_\sigma })\in {\mathbb {R}}^M$ is linear for every $\varphi \in {\mathcal {Z}}$, so that $\Pi (\overline{V_\varphi })$ is a polytope using Lemma 4. In summary, we obtain that $\Pi (\overline{D(P)_ SP })$ is a finite union of polytopes and consequently $C_{D(P)_ SP ,\Pi }$, as the convex hull of a finite union of polytopes, is a polytope as well. $\square$

Besides the fact that $C_{\overline{D(P)_ SP },\Pi }$ is a polytope, it is also crucial for our approach that there exists an efficient separation algorithm for this set. Indeed, for the special case

$$\begin{aligned} \begin{aligned} D(s)_ SP :=\big \{ u_{t_1,\dots ,t_\sigma }:&t_{i}-t_{i-1}\ge s ~ \forall \,i = 2, \dots ,\sigma ,~ t_1,\ldots ,t_\sigma \ge 0,\\ &u_{t_1,\dots ,t_\sigma }(\tau _j)=c_j\; \forall \,j = 1, \dots ,N\big \} \end{aligned} \end{aligned}$$

of dwell time constraints with fixings $(\tau _j,c_j)\in [0,T)\times \{0,1\}$, $1\le j\le N$, the separation problem is polynomially solvable in the number M of projection intervals, the number $\sigma$ of allowed switchings, and the number N of fixings, as claimed in Theorem 4.4. For the proof of the latter result, it remains to show the following lemma, using the definition of S given in Sect. 4.2.

Lemma A.6

Let v be a vertex of $C_{\overline{D(s)_ SP },\Pi }$. Then there exists $u\in \overline{D(s)_ SP }$ with $\Pi (u)=v$ such that u switches only in S.

Proof

Choose $c\in {\mathbb {R}}^M$ such that v is the unique minimizer of $c^\top v$ subject to $v\in C_{\overline{D(s)_ SP },\Pi }$. Moreover, choose any $u\in \overline{D(s)_ SP }$ with $\Pi (u)=v$ as well as a sequence $\{u^m\}_{m\in {\mathbb {N}}}\subset D(s)_ SP$ such that $u^m \rightarrow u$ in $L^p(0,T)$. Let $t^m_1,\dots ,t^m_\sigma$ be the switching points of $u^m$ for $m\in {\mathbb {N}}$, i.e., let $0\le t^m_1\le \dots \le t^m_\sigma <\infty$ such that $u_{t^m_1,\dots ,t^m_\sigma }=u^m$. Then there exists a subsequence of $t^m:= (t_1^m, \ldots , t_\sigma ^m)$ that converges to some $t\in {\mathbb {R}}^\sigma$ with $0\le t_1 \le \cdots \le t_\sigma <\infty$ and, thanks to Lemma 2, we have $u_{t_1,\ldots ,t_\sigma }=u$. For the following, for $j=1,\dots ,\sigma$ and $m\in {\mathbb {N}}\cup \{\infty \}$, we define

$$\begin{aligned} S^{m}_j:=\{t^m_\ell :\ell \in \{1,\dots ,\sigma \},\ t^m_\ell -t^m_j=s(\ell -j)\}, \end{aligned}$$

where we set $t^\infty :=t$. The set $S_j^m$ thus contains all switching points in $t^m$ that have the minimal possible distance to $t^m_j$.

Assume first that $t_j\in (a_i,b_i)\setminus S$ for some $i\in \{1,\dots ,M\}$ and some  $j\in \{1,\dots ,\sigma \}$. Due to $t^m \rightarrow t$ in ${\mathbb {R}}^\sigma$, we deduce for m sufficiently large that $t_j^m\in (t_j-\tfrac{\varepsilon }{2},t_j+\tfrac{\varepsilon }{2})$, where $\varepsilon >0$ is given by $\varepsilon :=\min _{q\in S}|t_j-q|>0$. Then $t_j^m\not \in S$ and $S^{m}_j\cap S=\emptyset$ by definition of S. Now all points in $S^{m}_j$ can be shifted by some $0<\delta <\tfrac{\varepsilon }{2}$, in both directions, maintaining feasibility with respect to $D(s)_ SP$, since none of these points is shifted to one of the fixing points $\tau _1,\ldots ,\tau _N$. Consequently, all points in $S^{\infty }_j$ can be slightly shifted simultaneously in both directions, maintaining feasibility with respect to $\overline{D(s)_ SP }$ and without any of these points leaving or entering any of the intervals $I_1,\dots ,I_M$ or [0, T]. This shifting changes the value of $c^\top \Pi (u)$ linearly, compare [6, Thm. 3.11], which is a contradiction to unique optimality of v.

We have thus shown that all switching points of u are either in S or outside of any interval $I_i$. Let $t_j\not \in S$ be any switching point of u not belonging to any interval $I_i$. Then, for sufficiently large m, we have $t_j^m\notin S$ and $t^m_j\notin I_i$ for any $i\in \{1,\ldots ,M\}$. The idea is, as in the proof of [6, Lemma 3.12], to shift the switching point $t_j^m\notin S$ for each m to the next point on the left belonging to S, but if this point belongs to $[0,T]\cap ({\mathbb {Z}}s+(\{\tau _j:j=1,\dots ,N\}))$, we can only shift $t^m_j$ arbitrarily close to the latter point in order to maintain feasibility in $D(s)_ SP$. For small enough $\delta >0$, we thus shift all switching points in $S^{m}_j$ simultaneously to the left until

$$\begin{aligned} \textrm{dist}({S^{m}_j},S):=\mathop {\min }\limits _{p\in S^{m}_j, q\in S}|p-q|=\delta , \end{aligned}$$

(15)

taking into account that the set $S^{m}_j$ may increase when $t^m_j$ decreases. Consequently, for all $\delta$, we obtain another sequence $\{u_{\delta }^{m}\}_{m\in {\mathbb {N}}}$. By construction, no switching point is moved beyond the next point in S to the left of its original position and no switching point is moved to any of the fixing points $\tau _1,\ldots ,\tau _N$, so that we conclude $u_{\delta }^{m}(\tau _j)=c_j$ for $j=1,\ldots ,N$ and thus $u_{\delta }^{m}\in D(s)_ SP$. In particular, none of the switching points being moved enters any of the intervals $I_i$, so that we derive

$$\begin{aligned} \Pi (u_{\delta }^{m})=\Pi (u^m)\rightarrow \Pi (u)=v\qquad \text{ for } m\rightarrow \infty \end{aligned}$$

(16)

by the continuity of the projection $\Pi$. We know that $\{u_{\delta }^{m}\}_{m\in {\mathbb {N}}}$ is a bounded sequence in BV(0, T) and hence by [1, Thm. 10.1.3 and Thm. 10.1.4] there exists a strongly convergent subsequence, which we again denote by $\{u_{\delta }^{m}\}_{m\in {\mathbb {N}}}$, such that $u_{\delta }^{m} \rightarrow u_{\delta }\in \overline{D(s)_ SP } \text { for } m\rightarrow \infty$. By (16) and the continuity of $\Pi$, we obtain $\Pi (u_{\delta })=v$ for $\delta >0$. Now $\{u_{\delta }:\delta >0\}\subset \overline{D(s)_ SP }$ is bounded in BV(0, T) as well, so that it contains an accumulation point $u'\in \overline{D(s)_ SP }$ and, again by the continuity of the projections, we have $\Pi (u')=v$. Thanks to (15), $u'$ then has at least one switching point more in S than u, but still satisfies $\Pi (u')=v$. By repeatedly applying the same modification, we eventually obtain a function projecting to v with all switching points in S. $\square$