1 Introduction

We consider mixed-integer optimal control problems (MIOCPs) on the fixed and finite time horizon \({\mathscr {T}}:=[t_0,t_f]\subset \mathbb {R}\) of the following form

$$\begin{aligned}&\inf _{\mathbf {x,v}}\varPhi (\mathbf {x}(t_f))\nonumber \\&\mathrm {s.\,t.} \qquad \dot{\mathbf {x}}(t) = \ \mathbf {f}(\mathbf {x}(t), \mathbf {v}(t)) \quad \quad \text { for } t\in \mathscr {T} , \end{aligned}$$
(1.1)
$$\begin{aligned}&\qquad \quad \mathbf {x}(t_0) = \mathbf {x}_0, \quad \end{aligned}$$
(1.2)
$$\begin{aligned}&\qquad \quad \; { \mathbf {v} \in V_{\mathrm {MDT}} \subset V. } \end{aligned}$$
(1.3)

The differential states \(\mathbf {x}\in W^{1,\infty }({\mathscr {T}},R^{n_{\text {x}}})\) with fixed initial values \(\mathbf {x}_0 \in \mathbb {R}^{n_{\text {x}}}\) are governed by the right-hand side ordinary differential equation (ODE) function \( \mathbf {f}: \mathbb {R}^{{n_{\text {x}}}}\times \{v_1,\ldots ,v_{{n_\omega }}\} \rightarrow \mathbb {R}^{n_{\text {x}}}\), which is assumed to be continuous in the first argument. We further assume that there exists a solution \(\mathbf {x}\) for the above problem. Let \(V=L^{\infty }({\mathscr {T}},\{v_1,\ldots ,v_{{n_\omega }}\})\) and \(\{v_1,\ldots ,v_{{n_\omega }}\} \subset \mathbb {R}^{n_v}\) so that the discrete control function \(\mathbf {v}: {\mathscr {T}} \rightarrow \{v_1,\ldots ,v_{{n_\omega }}\}\) is assumed to be measurable and is further restricted by minimum dwell time (MDT) constraints represented by the subset \( V_{\mathrm {MDT}}\). We stress that this function takes values out of a finite set with cardinality \({n_\omega }\in \mathbb {N}\) and exclude the trivial case \({n_\omega }= 1\). We minimize \(\varPhi \in C^0(\mathbb {R}^{n_{\text {x}}},\mathbb {R})\) over the end state, which in turn depends on the discrete control function \(\mathbf {v}\).

1.1 Related work

Optimal control problems with integer control choices have been investigated in many research articles during the last few decades [12, 14, 31, 40] since they naturally arise in a range of applications. In order to apply the obtained control policy in practice, additional switching constraints are usually needed, such as minimum dwell time requirements that describe the necessity of activating an integer control for at least a given minimal duration if at all and analogously for deactivation periods. Recent case studies of MIOCPs with MDT considerations can be found, e.g., for pesticide scheduling in agriculture [2], electric transmission lines [17], solar thermal climate systems [8] and hybrid electric vehicles [30]. MDT constraints attracted also a lot of attention as part of mixed-integer linear programs (MILPs), see [29] for a study of unit-commitment problems and [25] for a corresponding polytope investigation. For a recent work about model predictive control under MDT constraints see [10].

There are MIOCP related problems and approaches that have been discussed in the literature. For example, the optimal control community has been successfully solving so-called bang-bang problems for decades. In contrast to the above problem formulation, the considered linear control problems are required to have the bang-bang property, i.e., the derivative of the Hamiltonian is strictly positive or strictly negative almost everywhere. The main challenge consists in guessing the correct switching order and numerically detecting the switching points. This approach does not work for problems that involve singular or path-constrained arcs. This can be overcome by using a discrete global maximum principle (see [15] for further references). Still, the indirect first-optimize-then-discretize approaches have some drawbacks compared to direct methods and it is not clear how combinatorial constraints such as MDT constraints could be incorporated. As part of direct methods, one common approach to tackle MIOCPs under MDT constraints is to consider two separate levels of optimization. At the upper level, a mode insertion gradient is usually evaluated in order to fix a sequence of active system modes with promising cost function value. At the lower level, the algorithm aims at minimizing the cost functional with respect to the switching times and continuous control input, if available. Such approaches can be found in [1, 4, 12, 14, 40] and are usually referred to as transition-time optimization or variable time transformation method. We remark such switched systems have also been intensively investigated with respect to stability under average dwell time, see [20]. Another approach to include MDT constraints into MIOCPs is to apply dynamic programming, see [9], which is, however, computationally expensive.

1.2 Results for the CIA decomposition

Similar to the two-stage optimization method, the combinatorial integral approximation (CIA) decomposition [31, 35] that uses the same approximation theory as the embedding transformation technique [5], is based on the idea to solve the problem in several optimization steps. The first idea is to use the partial outer convexification method [31, 34] that allows us to reformulate a problem with integer controls into an equivalent problem with a \(\{0,1\}^{{n_\omega }}\)- valued control function \(\varvec{\omega }\), which indicates that exactly one discrete realization \(v_i\in \{v_1,\ldots ,v_{{n_\omega }}\}\) is active for each time point. For this, we replace the Constraint (1.1) by

$$\begin{aligned} \dot{\mathbf {x}}(t) = \sum \limits _{i=1}^{{n_\omega }}\omega _i(t) \mathbf {f}(\mathbf {x}(t), \mathbf {v}_i) \quad \text { and } \quad \sum \limits _{i=1}^{{n_\omega }}\omega _i(t) = 1, \qquad \text { for } t\in {\mathscr {T}}. \end{aligned}$$
(1.4)

Second, discretizing the MIOCP in the spirit of first-discretize-then-optimize methods with, e.g., Direct Multiple Shooting [6] or Direct Collocation [37] results in a mixed-integer nonlinear program (MINLP). We introduce a time discretization by the ordered set \(\mathscr {G}_N:=\{t_0< \cdots < t_N=t_f \}\) denoting a grid with N intervals and lengths \(\varDelta _j:=t_{j}-t_{j-1}, \, \bar{\varDelta }:=\max _j \varDelta _j\) as well as \(\underline{\varDelta }:=\min _j \varDelta _j\) for \(j\in [N]\). The binary control functions \(\omega _i(\cdot )\) are assumed be piecewise constant, changing values only on these grid points, so that \(\omega (\cdot )\) can be uniquely represented by a matrix \(\mathbf {w}\in \{0,1\}^{{n_\omega }\times N}\). This MINLP becomes a nonlinear program (NLP), if the binary controls are relaxed. We write \(\mathbf {a}\in [0,1]^{{n_\omega }\times N}\) for denoting this relaxed value. After solving this NLP, the resulting \(\mathbf {a}^\star \) is approximated in the CIA step - which is an MILP - with binary control values \(\mathbf {w}\). The main idea of this CIA problem is to minimize the integrality gap \(\theta (\mathbf {w})\), which is the accumulated control deviation between \(\mathbf {a}^\star \) and \(\mathbf {w}\), i.e., \(\min _{\mathbf {w}} \theta (\mathbf {w})\), where

$$\begin{aligned} \theta (\mathbf {w}) :=\max \limits _{\begin{array}{c} i=1,\ldots ,{n_\omega }\\ k = 1,\ldots ,N \end{array}} \left| \sum \limits _{j=1}^k (a_{i,j}^\star -w_{i,j})\varDelta _j \right| . \end{aligned}$$

This work is based on the CIA decomposition due to its advantages:

  1. 1.

    MINLPs fall generally into the class of NP-hard problems so that using approaches that bypass the direct solution of such problems is computationally favorable.

  2. 2.

    Convergence results have been proven for MIOCPs without MDT constraints in the sense that, under mild assumptions, the obtained solution was shown to be arbitrarily close to the optimal solution with grid length \(\bar{\varDelta }\) going to zero [26, 31, 34]. Moreover, the solution of the NLP represents a useful a priori lower bound on the objective, if solved to global optimality.Footnote 1

  3. 3.

    An MILP enables the option to include intuitively a large variety of combinatorial constraints. Numerical case studies showed that the resulting feasible solution is close to the relaxed solution in case the applied combinatorial constraints are not too restrictive [7, 8].

For a generalization of the CIA decomposition to more optimization steps and MILP variants see [41] as well as to the PDE constraint case see [18, 19, 27]. Recently, an extension of the algorithm with the penalty alternating direction method has been made [16], which can be regarded as a feasibility pump [13] variant for MIOCPs.

Instead of minimizing \(\theta (\mathbf {w})\) by means of solving an MILP, fast rounding algorithms can be applied, such as sum-up rounding (SUR) [31] and next forced rounding (NFR) [21], which generate binary solutions that still converge to the relaxed binary control solution with vanishing grid length. The SUR scheme computes for \(j=1,\ldots ,N\)

$$\begin{aligned} w_{i,j} :={\left\{ \begin{array}{ll} 1, &{} \text { if } i=\mathop {{\text {argmax}}}\limits _{k=1,\ldots ,{n_\omega }} \left\{ \sum \limits _{l=1}^j a_{k,l}^\star \varDelta _l - \sum \limits _{l=1}^{j-1} w_{k,l}\varDelta _l \right\} \qquad \text {(break ties arbitrarily),} \\ 0, &{} \text { else, } \end{array}\right. } \qquad \text {for } i=1,\ldots ,{n_\omega }. \end{aligned}$$

For defining the NFR algorithm, one needs for all \(i=1,\ldots ,{n_\omega }\) and iteratively for \(j=1,\ldots ,N\) the quantity

$$\begin{aligned} \mathscr {N}_j(i):=\left\{ \begin{array}{ll} \mathop {{\text {argmin}}}\limits _{k=j,\ldots ,N} \{ \sum _{l=1}^{k} a_{i,l}^\star \varDelta _l - \sum _{l=1}^{j-1} w_{i,l}\varDelta _l> \bar{\varDelta } \}, \qquad &{} \quad \text {if } \quad \sum _{l=1}^{N} a_{i,l}^\star \varDelta _l - \sum _{l=1}^{j-1} w_{i,l}\varDelta _l > \bar{\varDelta },\\ \infty , &{} \quad \text {else.} \end{array} \right. \end{aligned}$$
(1.5)

A control with index \(i \in [{n_\omega }]\) on interval j is defined to be next forced, if and only if

$$\begin{aligned} \mathscr {N}_j(i) = \min \limits _{k\in [{n_\omega }]}\mathscr {N}_j(k) \quad \text {and} \quad \mathscr {N}_j(i)<\infty . \end{aligned}$$
(1.6)

Then, the NFR algorithm sets iteratively for \(j=1,\ldots ,N\) the next forced control equal to one (break ties arbitrarily) and if there is no such control, the active control is chosen according to the SUR scheme. We summarize established integrality gap bounds for these schemes in the form of \(\theta (\mathbf {w})\le C({n_\omega }) \bar{\varDelta }\) in Table 1.

Table 1 So far known integrality gap bounds for binary control approximation algorithms
Full size table

1.3 Contribution

We conduct a theoretical analysis of the integrality gap in the presence of MDT constraints. In particular, we prove the upper bound

$$\begin{aligned} \min _{\mathbf {w}} \theta (\mathbf {w}) \le \tfrac{2{n_\omega }-3}{2{n_\omega }-2}(C_{UD}+\bar{\varDelta }), \end{aligned}$$

where \(C_{UD}\) denotes the maximum of the minimum up (MU) and minimum down (MD) duration of all controls. We show that this bound is tight for MU time constraints. As a consequence of this result, the tightest possible bound for the CIA problem is \( C({n_\omega })=\frac{2{n_\omega }-3}{2{n_\omega }-2}.\) The proof is constructive, as we introduce a generalization of the NFR scheme to the MDT setting. We present a rounding modification so that MD times can be included explicitly allowing us to deduce an improved integrality gap bound for this case. In the same spirit, we modify SUR so that MDT requirements are satisfied by the obtained binary solution. We test our algorithms on the three tank problem from a benchmark library [33] and evaluate with increasingly restrictive MDT constraints how large the gap between the constructed integer feasible and relaxed solution becomes.

1.4 Outline

We give a problem definition of the MIOCP in Sect. 2 and describe the proposed CIA decomposition algorithm with (CIA) as subproblem in detail in Sect. 3. A generalization of NFR to the MDT setting is presented in Sect. 4. It provides the tools to derive upper bounds on the integrality gap in Sect. 5. The commonly used SUR scheme can be extended to include also MDT constraints, which we present in Sect. 6 together with a discussion of the integrality gap. Section 7 provides a numerical example and we finish this article with conclusions in Sect. 8.

1.5 Notations

Let \([n] :=\{1, \ldots , n\}, [n]_0 :=\{0\}\cup [n],\) for \( n \in \mathbb {N}\). We use Gauss’ bracket notation, i.e., \(\lfloor x \rfloor :=\max \{k\in \mathbb {Z}\, \mid \, k\le x\}, \, x\in \mathbb {R}\), and analogously for \(\lceil x \rceil \).

2 Mixed-integer optimal control problem

We introduce our problem at hand in the discretized and (partial) outer convexified setting, but refer to [22, 31] and, more recently, to [23, 24] for extensive descriptions of partial outer convexification, the continuous MIOCP setting and its relation to the discretized problem.

2.1 Definition of binary and relaxed controls

We base the following definitions on the grid \(\mathscr {G}_N\) we already introduced in the previous section.

Definition 1

(Convex combination constraint (Conv) and Matrix sets W, A) We express the requirement that the columnwise entries of a matrix \((m_{i,j})\in [0,1]^{{n_\omega }\times N}\) sum up to one by

$$\begin{aligned} \sum \limits _{i\in [{n_\omega }]} m_{i,j} =1, \qquad \text { for } j\in [N], \end{aligned}$$
(Conv)

and call it convex combination constraint (Conv) in the sequel. Based on this constraint, we define the sets

$$\begin{aligned} W:=\left\{ \mathbf {w} \in \{0,1\}^{{n_\omega }\times N} \, \mid \, \mathbf {w} \text { fulfills } (Conv) \right\} , \quad A:=Conv( W), \end{aligned}$$

where we denote by Conv(W) the convex hull of W.

We notice the geometric nature of W and A. They are the vertices respectively the set of faces of the N-fold iterated standard simplex without the origin and spanned by the \({n_\omega }\) unit vectors.

Definition 2

(Binary \(\varvec{\omega }\) and relaxed control functions \(\varvec{\alpha }\)) Let the vector of binary control functions \(\varvec{\omega }\) and its corresponding vector of relaxed control functions \(\varvec{\alpha }\) be defined by their function space domains

$$\begin{aligned} \varOmega :=\{\varvec{\omega }: {\mathscr {T}}\rightarrow \{0,1\}^{n_\omega }\, \mid \, \omega (t)= & {} \mathbf {w}_{\cdot ,j}, \,\, \text { for } t\in [t_{j-1},t_j),\,\, t_{j-1},t_j\in \mathscr {G}_N, \,\, j\in [ N ], \,\, \mathbf {w} \in W \}, \\ \mathscr {A} :=\{ \varvec{\alpha }: {\mathscr {T}}\rightarrow [0,1]^{n_\omega }\, \mid \, \alpha (t)= & {} \mathbf {a}_{\cdot ,j}, \,\, \text { for } t\in [t_{j-1},t_j), \,\, t_{j-1},t_j\in \mathscr {G}_N,\,\, j\in [ N ],\,\, \mathbf {a} \in A \}, \end{aligned}$$

where \(\mathbf {w}_{\cdot ,j}\), respectively \(\mathbf {a}_{\cdot ,j}\), denotes the jth column of \(\mathbf {w}\), respectively of \(\mathbf {a}\).Footnote 2

2.2 Optimal control problem class

We take interest in the discretized binary control problem (DBCP) with MDT constraints and its naturally relaxed problem, the discretized relaxed control problem (DRCP), which arises from replacing \(\omega \) by \(\alpha \).

Definition 3

(Problems (DBCP) and (DRCP)) Let a fixed discretization grid \(\mathscr {G}_N\) be given together with an MU time \(C_U\ge 0\) and an MD time \(C_D\ge 0\). We refer to the following general problem class as (DBCP)

figure a

where we denote the intervals affected by the MDT \(C_1=C_U,C_D\) from interval \(k\in [N]\) on with the set

$$\begin{aligned} {\mathscr {J}}_{k}(C_1):=\{k\} \cup \{j \ \mid \ t_{j-1}\in \mathscr {G}_N \cap [t_{k-1} ,t_{k-1}+ C_1) \}. \end{aligned}$$

We define a control \(i\in [{n_\omega }]\) to be active on the interval starting with \(t_{j-1}\in \mathscr {G}_N\), if and only if \(\omega _i(t)=1\) for \(t\in [t_{j-1},t_j)\) is true and the other way around for inactive controls. If a binary control is active after a switch on \(t_j\), it has to stay active for a time period of at least \(C_U\) as required by (2.1), whereas (2.2) enforces the analogous case for deactivating a control. Finally, we define (DRCP) as the canonical relaxation of problem (DBCP) where we optimize over \(\varvec{\alpha }\in \mathscr {A}\) instead of \(\varvec{\omega }\in \varOmega \).

We remark that our study assumes no mode specific MU times \(C_{i,U}\) or MD times \(C_{i,D}\), but may include them by setting \(C_U=\max _{i\in [{n_\omega }]}C_{i,U}\) and accordingly \(C_D=\max _{i\in [{n_\omega }]}C_{i,D}\), even though this simplification may result in suboptimal solutions.

Without loss of generality, we omit in our problem definition further constraints and continuous controls \(\mathbf {u}\in L^\infty ({\mathscr {T}},\mathbb {R}^{n_u})\). See [32] for a discussion of extensions and [27] for PDE constraints.

3 Combinatorial integral approximation decomposition

We recapitulate the CIA decomposition algorithm to solve MIOCPs in Fig. 1. This approach is justified by a convergence result in the situation without MDT constraints establishing that the objective integer gap of (DBCP) compared to (DRCP) depends linearly on the integrated control deviation under certain regularity assumptions [26, 31].

Fig. 1
figure 1

Schematic representation of the CIA decomposition. MIOCPs can be equivalently reformulated into their partially outer convexified counterpart problem that is thereafter transformed into (DBCP) via a temporal discretization. Next, allowing a convex combination \(\mathbf {a}\) in (Conv) yields the relaxed problem (DRCP). After solving this problem we obtain \(\mathbf {a}^\star \), which are then approximated with binary values \(\mathbf {w}^\star \) in the rounding problem (CIA). Finally, we use \(\mathbf {w}^\star \) as fixed variables for solving (DBCP) as a continuous variable problem

Full size image

Thus, we are interested in an admissible \(\omega \) represented by its corresponding value \(\mathbf {w}\in W\) with integrality gap bounded by the grid length in the sense of

$$\begin{aligned} \min \limits _{\mathbf {w}\in W} \max \limits _{i\in [{n_\omega }],\, k \in [N]} \left| \sum \limits _{j\in [k]} ( a_{i,j}^\star -w_{i,j})\varDelta _j \right| \le C({n_\omega }) \bar{\varDelta }, \end{aligned}$$
(3.1)

as already mentioned in the introduction, so that the solution of (DRCP) can be arbitrarily well approximated by an integer solution through refining the discretization grid [34]. The SUR scheme constructs solutions \(\omega \) with this property; in fact, Kirches et al.  [24] showed \(C\in \mathscr {O}(\log ({n_\omega }))\) as conjectured in [34]. Hence, the integrality gap obtained by this rounding scheme is going to be arbitrarily small with arbitrarily small grid length, which is why SUR is often applied in the CIA decomposition as rounding step after solving (DRCP). Rather than using SUR, it has been proposed [35] to formulate the problem (3.1) as an MILP for an improved approximation and to be able to consider also different norms (such as the Manhattan norm) and combinatorial constraints. Therefore, we recall its definition and state its MDT variants.

Definition 4

(Problems (CIA), (CIA-U), (CIA-D), and (CIA-UD)) Let \(\mathbf {a}^\star \in A\) be the (local) optimal solution of (DRCP) and assumed to be given. Then, we define the problem (CIA) to be

$$\begin{aligned} \theta ^\star :=\underset{\theta \in \mathbb {R}_{\ge 0},\, \mathbf {w}\in W}{\mathrm {min}}&\quad \theta \\ \mathrm {s.\,t.}&\quad \theta \ge \pm \sum \limits _{l\in [j]} (a_{i,l}^\star -w_{i,l})\varDelta _l, \qquad \text { for } i\in [{n_\omega }], \, j\in [N]. \end{aligned}$$

The (CIA) problem with added MU time constraint (2.1) from Definition 3 is in the remainder referred to as (CIA-U). Similarly, let us (CIA) with added MD time constraint (2.2) call (CIA-D) and (CIA) with both (2.1) and (2.2) (CIA-UD).

Clearly, (CIA) is a reformulation of minimizing (3.1). In Sect. 5 we are going to elaborate upper bounds for (3.1) in the presence of MDT constraints by investigating its (CIA) variants. Since the lower bound on its objective is trivially zero (and can be reached), bounds always refer to upper bounds in this article.

With respect to the constraints (2.1) and (2.2) we stress that there are other, in many cases computationally more efficient, formulations of MDT constraints, e.g. in the spirit of extended formulations [25]. The latter may lead to relaxations that are less likely to deliver fractional solutions [25] and thus can be beneficial for including the constraints (2.1) and (2.2) into the NLP solving procedure. But, as this study does not focus on the NLP formulation, i.e., how \(\mathbf {a}^\star \) is achieved, and since we propose to solve (CIA-UD) by means of tailored rounding heuristics or a branch and bound algorithm and not with a standard MILP solver exploiting extended formulations, we would benefit neither numerically nor theoretically from alternative MILP formulations and we therefore skip the presentation of these.

4 Dwell time next forced rounding

We introduced in (1.5)–(1.6) the NFR scheme as an algorithm that can compute approximations to solutions of (CIA) in \(\mathscr {O}({n_\omega }N^2)\) [21] and that constructs feasible solutions of (CIA) with objective no larger than \(\bar{\varDelta }\). In this section, we introduce dwell time next forced rounding (DNFR) as a generalization, aiming for a scheme that constructs MDT constraint feasible solutions and from which we derive bounds for the (CIA) objective and its MDT variants. Several definitions are needed for DNFR and we begin with a definition of sequences of intervals that are grouped into blocks in the presence of MDT constraints.

Definition 5

(Dwell time block interval sets) Let an MDT \(C_1 \ge 0\) be given. We define iteratively the dwell time invoked interval sets \({\mathscr {J}}_b\) and their last indices \(l_{b}\) for \(b=1,\ldots ,n_b\) and with \(l_0:=0\):

$$\begin{aligned} {\mathscr {J}}_{b}&:=\{l_{b-1}+1\} \cup \{j \ \mid \ t_{j-1}\in \mathscr {G}_N \cap [t_{l_{b-1}} ,t_{l_{b-1}}+ C_1) \}, \\ l_{b}&:=\max \{j \ | \ j \in {\mathscr {J}}_{b}\}, \end{aligned}$$

where \(n_b :=\min \{ b \ \mid \ l_{b}=N\}\) represents the number of interval blocks.

In the following we will sometimes write loosely block instead of dwell time block for shortening our language. Next, we establish the lengths of dwell time blocks.

Definition 6

(Dwell time block length) Let a family of dwell time block interval sets \(\{{\mathscr {J}}_b\}_{b\in [n_b]}\) be given. We denote by \(\mathscr {L}_b\) the length of dwell time block \(b\in [n_b]\) and name the maximum, respectively minimum, length of all dwell time blocks \(\overline{\mathscr {L}}\), respectively \(\underline{\mathscr {L}}\), i.e.,

$$\begin{aligned} \mathscr {L}_b&:=t_{l_b}-t_{l_{b-1}}, \qquad b\in [n_b],\\ \overline{\mathscr {L}}&:=\max \limits _{b\in [n_b]} \mathscr {L}_b, \qquad \underline{\mathscr {L}} :=\min \limits _{b\in [n_b]} \mathscr {L}_b. \end{aligned}$$

By the definition of dwell time blocks, we see that \(\mathscr {L}_b\) depends both on the time discretization \(\mathscr {G}_N\) and on \(C_1\). If \(C_1 \le \underline{\varDelta }\), then, the blocks are in fact the grid intervals, i.e., \(\mathscr {L}_j=\varDelta _j, \ j\in [N]\) and \(n_b=N\). As soon as \(C_1 > \underline{\varDelta }\) holds, there is at least one block b with length of two subsequent intervals \(\mathscr {L}_b=\varDelta _j+\varDelta _{j+1}, \ j\in [N-1]\). Overall, one recognizes that \(\overline{\mathscr {L}}\) increases monotonically with increasing \(C_1\) and obviously stops as soon as \(C_1 > t_f-t_0\). The DNFR scheme relies crucially on the block dependent accumulated control deviation, which is why we introduce it as auxiliary variable in the next definition.

Definition 7

(Accumulated control deviation \(\theta _{i,j},\varTheta _{i,b}, \gamma _{i,j}, \varGamma _{i,b}\)) Let \(\mathbf {a}\in A\) and \(\mathbf {w}\in W\). For controls \(i\in [{n_\omega }]\) we define the accumulated control deviation on interval \(j\in [N]\) as

$$\begin{aligned} \theta _{i,j} :=\sum \limits _{l=1}^j (a_{i,l}-w_{i,l}) \varDelta _l, \qquad \gamma _{i,j} :=\sum \limits _{l=1}^j a_{i,l}\varDelta _l - \sum \limits _{l=1}^{j-1}w_{i,l} \varDelta _l, \end{aligned}$$

and further define \(\theta _{i,0} :=0\) for all \(i\in [{n_\omega }]\). We introduce for blocks \(b\in [n_b]\) the analogous variables

$$\begin{aligned} \varTheta _{i,b} :=\theta _{i,l_{b}}, \qquad \varGamma _{i,b}:=\varTheta _{i,b-1} + \sum \limits _{j\in {\mathscr {J}}_b}a_{i,j}\varDelta _j. \end{aligned}$$

In the sequel, we sometimes write forward control deviation for control i in order to distinguish \(\gamma _{i,j}, \varGamma _{i,b}\) from the (accumulated) control deviation \(\theta _{i,j}, \varTheta _{i,b}\).

The defined quantities of Definitions 5–7 are illustrated in Fig. 2.

Remark 1

If a small MDT \(C_1 \le \underline{\varDelta }\) is given, \(\theta _{i,j}\) equals trivially \(\varTheta _{i,j}\) and the same holds for \(\gamma _{i,j}\) and \(\varGamma _{i,j}\). Nevertheless, with an MDT of \(C_1 > \underline{\varDelta }\) one needs interval and block related variables to be able to clearly distinguish between both values.

Remark 2

Let \(\mathbf {w}\in W\) and we denote by \(\theta (\mathbf {w})\) its (CIA) objective value. With the above definition, we conclude \(\theta (\mathbf {w})=\max \nolimits _{i\in [{n_\omega }], j\in [N]} |\theta _{i,j}|\). Generally, one notices that the maximum of the \(|\theta _{i,j}|\) values must be assumed at an interval before a switch happens (i.e., \(\mathbf {w}_{\cdot ,j}\ne \mathbf {w}_{\cdot ,j+1}\)) or on the last interval, since \(|\theta _{i,j}|\) increases monotonically with constant \(\mathbf {w}_{\cdot ,j}\) and increasing j. Hence, with constant \(\mathbf {w}_{\cdot ,j}\) on the dwell time blocks, we also have that \(\theta (\mathbf {w})=\max \nolimits _{i\in [{n_\omega }], b\in [n_b]} |\varTheta _{i,b}|\).

Fig. 2
figure 2

Left: example binary and relaxed control values on nine intervals \(\varDelta _j\), respectively three blocks \(\mathscr {L}_b\). Right: corresponding accumulated control deviation. The forward control deviation with respect to intervals \(\gamma _i\) is greater or equal to \(\theta _i\) with equality if the control i is inactive. Also, \(\varGamma _{i,b}\) is greater or equal to \(\varTheta _{i,b}\), where the latter equals via definition the last \(\theta _{i}\) before block b begins. Notice that \(\varGamma _{i,2}\) sums up the weighted \(\alpha \) values for the intervals 3-6, which results in a large value

Full size image

We introduced in (1.5)–(1.6) the concept of a next-forced control that depends on the maximum grid length \(\bar{\varDelta }\). We generalize this concept by using blocks and a generic rounding threshold factor \(C_2>0\) instead of using always \(C_2=1\) as in the NFR scheme. To this end we present a definition of different types of control variable activations that depend on prior variables choices and on \(\mathbf {a}\in A\).

Definition 8

(Admissible, forced, and future forced activation) Let the rounding threshold factor \(C_2>0\) and \(\mathbf {a}\in A\) be given. The choice \(w_{i,j}=1\) for \(i\in [n_{\omega }], \ j \in {\mathscr {J}}_b, \ b=1,\ldots , n_b\) is admissible, if

$$\begin{aligned} \varGamma _{i,b} \ge - C_2 \overline{\mathscr {L}} + \mathscr {L}_b \end{aligned}$$

holds. Denote by \(W^b_a\) the set of admissible controls for block b. Similarly, the choice \(w_{i,j}=1\) for \(i\in [n_{\omega }], \ j \in {\mathscr {J}}_b, \ b=1,\ldots , n_b\) is forced, if

$$\begin{aligned} \varGamma _{i,b} > C_2 \overline{\mathscr {L}} \end{aligned}$$

holds. We define a control \(i\in W^b_a\) on block b to be l-future forced, if

$$\begin{aligned} \varTheta _{i,b-1} + \sum \limits _{k=b}^{l} \sum \limits _{j\in {\mathscr {J}}_k} a_{i,j}\varDelta _j > C_2 \overline{\mathscr {L}} \end{aligned}$$

holds, with the special case \(l=b\) meaning i is actually forced on b. If the above inequality holds for any \(l\le n_b\), we call the control \(i\in W^b_a\) on block b to be future forced, and group these controls into the set \(W^b_{f}\).

Definition 9

(Minimum down time forbidden control) We introduce the constant \(\chi _D\in \{0,1\}\). If the CIA problem involves an MD time constraint with \(C_D > \underline{\varDelta }\), we set \(\chi _D=1\) and otherwise \(\chi _D=0\). We define \(i_b^D, \ b=3,\ldots ,{n_\omega },\) as the index of the control that has been activated on block \(b-2\) and deactivated on block \(b-1\) - if such a control exists:

$$\begin{aligned} \exists i \in [{n_\omega }] \ : \ w_{i,j} =1, j\in {\mathscr {J}}_{b-2} \ \ \wedge \ \ w_{i,j} =0, j\in {\mathscr {J}}_{b-1} \quad \Rightarrow \quad i_b^D :=i. \end{aligned}$$

Then, let \(\mathscr {I}_b^D\) denote the \(\chi _D\) dependent set of the minimum down time forbidden control

$$\begin{aligned} \mathscr {I}_b^D :={\left\{ \begin{array}{ll} \{i_b^D\} ,\qquad &{} \qquad \text {if } \chi _D=1, \text { and } b\ge 3,\\ \emptyset , &{} \qquad \text {otherwise.} \end{array}\right. } \end{aligned}$$

Note that \(\mathscr {I}_b^D\) is either the empty set or contains exactly one control index. It may seem unintuitive to declare only one control per block as minimum down time forbidden, because a sufficiently large chosen MD time can comprise more than two intervals and therefore more than one control could be minimum down forbidden on certain blocks. However, in such a situation, where several controls are minimum down forbidden, it could happen that only one control may be allowed to be active, resulting in a large control deviation. Consequently, a fine granular definition would be critical for deriving bounds for (CIA-UD) using the DNFR scheme. We will specify such a worst case later in Example 1 as part of “Appendix B” and argue thereby why we tolerate at most one minimum down time forbidden control per block. We illustrate the different control activation types of Definitions 8 and 9 in Fig. 3.

Fig. 3
figure 3

Exemplary visualization of the defined quantities. Left: binary and relaxed control values on four blocks. Right: corresponding block accumulated control deviation. Control i is admissible on block \(b_{j}\), not admissible on block \(b_{j+1}\), down time forbidden and \(b_{j+3}\)-future forced on block \(b_{j+2}\) as well as forced on block \(b_{j+3}\)

Full size image

Finally, we use these control properties to declare the DNFR scheme in Algorithm 1. In contrast to the original NFR, we do not iterate over all intervals, but over all dwell time \(C_1\) invoked blocks (line 2) and check on each block whether there is a forced control and activate it in this case (line 3–4). Otherwise, we test if there is an earliest future forced control and if so, we set it to be active (line 5–8). Else, the algorithm selects the control with the maximum forward control deviation (line 9–13), which represents a fallback to the classical SUR scheme. In case the MD mode is turned on by setting \(\chi _D=1\), we exclude the set \(\mathscr {I}_b^D\) from our control selection task (line 3, 5, 11). This consideration of minimum down time forbidden controls is a further extension of the original NFR scheme.

figure b

5 Tight bounds on the integrality gap for (CIA) with dwell time constraints

We now consider the question of how large the objective function value \(\theta ^\star \) of (CIA) and its MDT variants can become. In other words we examine

$$\begin{aligned} \theta ^{\max }:=\max \limits _{\mathbf {a}\in A} \min \limits _{\mathbf {w} \in W} \max \limits _{i\in [n_{\omega }],j\in [N]} | \theta _{i,j} | \quad \text {s.t.} \quad \text {{MDT} constraints } (2.1), (2.2). \end{aligned}$$

It turns out that the DNFR scheme is particularly suitable for deriving these integrality gap bounds. We state approximation results for (CIA) by means of DNFR constructed solutions. These results are presented as two theorems in Sect. 5.1 with specified parameter choices for \(C_2\) and \(\chi _D\). Afterwards, we deduce specific bounds for (CIA-U), (CIA-D) and (CIA-UD) and evaluate how tight they are in the Sects. 5.2–5.3. We begin this section with the trivial upper bound

$$\begin{aligned} \theta ^{\max } \le \sum \limits _{j\in [N]} \varDelta _j = t_f - t_0 \end{aligned}$$

and another weak result in the following remark.

Remark 3

Neglecting for a moment MDT constraints, it is known from Minimax theory [39] that

$$\begin{aligned} \max \limits _{\mathbf {a}\in A} \min \limits _{\mathbf {w} \in W} \max \limits _{i\in [n_{\omega }],j\in [N]} | \theta _{i,j} | \ \ \le \ \ \min \limits _{\mathbf {w} \in W} \max \limits _{\mathbf {a}\in A} \max \limits _{i\in [n_{\omega }],j\in [N]}| \theta _{i,j} | \end{aligned}$$

holds. In the right hand side, we maximize over \(\mathbf {a}\) for given \(\mathbf {w}\) and check, which one of the latter leads to an overall minimum objective. In this way \(\mathbf {a}\) manipulates the control deviation to be as large as possible. That means with given \(\mathbf {w}\) it is possible to set \(a_{i^{\min },j}=1, j\in [N],\) where \(i^{\min }\) is the control with smallest total accumulation \(\sum \nolimits _{j\in [N]} \mathbf {w}_{i,j} \varDelta _j\) so that we obtain the (CIA) objective value \(\theta ^\star =\sum \nolimits _{j\in [N]}(1 - \mathbf {w}_{i^{\min },j})\varDelta _j\). With these arguments one may derive

$$\begin{aligned} \theta ^{\max } \le \left( N- \left\lfloor \frac{N}{{n_\omega }}\right\rfloor \right) \bar{\varDelta }. \ \end{aligned}$$

We omit the exact proof since this bound is generally weak as we will see later in this section.

5.1 Integrality gap results through dwell time next forced rounding

We examine in Theorem 1 how large the control deviation can become as part of the DNFR algorithm during an MD time phase. Based on this result we derive in Theorem 2 that DNFR constructs (CIA) feasible solutions with objective bounds depending on the rounding threshold \(C_2\) and whether down time forbidden controls are allowed, i.e., \(\chi _D=1\). The proofs and the corresponding lemmata are moved to “Appendix A” to enhance readability for readers interested in the results and algorithms.

Theorem 1

(Control accumulation of a minimum down time forbidden control) Let \(\mathbf {a}\in A\), \((C_2,\chi _D)=\left( \tfrac{3}{2},1 \right) \) and \(C_1\ge 0\) be given and assume there is a minimum down time forbidden control \(i_{D}\in \mathscr {I}_b^D\) on block \(b\ge 3\) after Algorithm 1 was executed. Then, the forward control deviation satisfies

$$\begin{aligned} \varGamma _{i_{D},b} \le \tfrac{3}{2} \overline{\mathscr {L}}. \end{aligned}$$

Proof

See “Appendix A.2”. \(\square \)

Theorem 2

(Rounding gap of solution constructed by DNFR) Let \(\mathbf {a}\in A\) and the following parameter settings be given:

  1. I.

    \((C_2,\chi _D)=\left( \frac{2{n_\omega }-3}{2 {n_\omega }- 2}, 0 \right) \),

  2. II.

    \((C_2,\chi _D)=\left( \frac{3}{2},1 \right) \),

and \(C_1\ge 0\). Then, \(\mathbf {w}^{\text {DNFR}}\) obtained by Algorithm 1 is a feasible solution of (CIA) for both cases with approximation quality

$$\begin{aligned} \theta (\mathbf {w}^{\text {DNFR}}) \ \le \ C_2 \overline{\mathscr {L}}. \end{aligned}$$

Proof

See “Appendix A.3”. \(\square \)

5.2 Implications for the Objectives of (CIA-U) and (CIA)

Theorem 2 states only generic approximation results for (CIA) with an MDT parameter \(C_1\). We are going to assess the consequences for (CIA-U) by specifying \(C_1\) and proving tightness of the resulting upper bound. Clearly, (CIA) is a special case of (CIA-U), where \(C_U=0\), so that results for (CIA-U) are inherited by (CIA).

Proposition 1

(Upper bound for (CIA-U)) Let any MU time \(C_U\ge 0\), grid \(\mathscr {G}_N\) and \(\mathbf {a}\in A\) be given. Then, for (CIA-U) holds:

$$\begin{aligned} \theta ^\star \le \frac{2{n_\omega }- 3}{2{n_\omega }- 2} \left( C_U + \bar{\varDelta } \right) . \end{aligned}$$

Proof

We consider the DNFR scheme with \((C_1,C_2,\chi _D)=\left( C_U,\frac{2{n_\omega }- 3}{2{n_\omega }- 2},0\right) \). Then, \(\mathbf {w}^{\text {DNFR}}\) is a feasible solution by Theorem 2 and by the property of DNFR to activate dwell time blocks of intervals with length at least \(C_1=C_U\). From the definition of block lengths we conclude \(\overline{\mathscr {L}} < C_U + \bar{\varDelta }\) so that the assertion follows directly from Theorem 2. \(\square \)

In order to analyze the obtained bound with respect to tightness, we introduce a grid where the MDT \(C_1\) overlaps the grid points by a small \(\varepsilon >0\). We determine the length of the resulting blocks in the following lemma.

Definition 1

(Minimal \(C_1\)-overlapping grid) Let us consider a non-degenerate MDT length, i.e., \(C_1>0\), and let \(\varepsilon \) be \(C_1 \gg \varepsilon > 0\). Let further a time horizon \([t_0,t_f]\) be given with length at least \(4 C_1\), i.e.,

$$\begin{aligned} t_f\ge t_0+4 C_1 . \end{aligned}$$

We define a minimal \(C_1\)-overlapping grid \(\mathscr {G}_N\) recursively as follows

$$\begin{aligned} t_1&:=t_0 + C_1 - \varepsilon , \\ t_2&:=t_1 + C_1, \\ t_j&:={\left\{ \begin{array}{ll} t_{j-1} + C_1 - \varepsilon ,\qquad &{} \text {if } j \text { odd,}\\ t_{j-1} + C_1 ,\qquad &{} \text {if } j \text { even,} \end{array}\right. } \qquad \text { for } j=3,\ldots , N-1, \end{aligned}$$

where we set \(N-1:=\max \{ j \, |\, t_j < t_f\}\), so that \(\mathscr {G}_N\) consists of N intervals (Fig. 4).

Fig. 4
figure 4

Visualization of a minimal \(C_1\)-overlapping grid

Full size image

Lemma 1

(Length of blocks of a minimal \(C_1\)-overlapping grid) The dwell time invoked blocks \({\mathscr {J}}_b, \ b\in [n_b]\) of a minimal \(C_1\)-overlapping grid as introduced in Definition 1 have the length

$$\begin{aligned} \mathscr {L}_b&= 2 C_1 - \varepsilon , \qquad \qquad b\in [n_b-1],\\ \mathscr {L}_{n_b}&= t_f - \left( t_0 + (n_b-1) (2 C_1 - \varepsilon )\right) . \end{aligned}$$

Moreover, we have that

$$\begin{aligned} \bar{\varDelta } = C_1, \qquad \overline{\mathscr {L}}=2 C_1 - \varepsilon . \end{aligned}$$

Proof

We keep in mind Definition 5 from which we deduce \({\mathscr {J}}_1=\{1,2\}\), because \(t_0+C_1\) (minimally) overlaps \(t_1\). The next dwell time block begins at \(t_2= t_0 + 2 C_1 - \varepsilon \) and consists again of two intervals. This argumentation can be extended to the first \(n_b-1\) blocks and by the definition of block lengths we conclude \(\mathscr {L}_b= 2 C_1 - \varepsilon \). The length of the last block \(\mathscr {L}_{n_b}\) is directly computed by the definition of \(N-1\) to be the last index of the grid point recursion before \(t_f\). Finally, the definition of a minimal \(C_1\)-overlapping grid and the obtained block lengths implies

$$\begin{aligned} \bar{\varDelta } = C_1, \qquad \overline{\mathscr {L}}=2 C_1 - \varepsilon . \end{aligned}$$

\(\square \)

With this preliminary work, we show that the deduced MU time bound is tight.

Proposition 2

(Tightness of the bound for (CIA-U)) Let an MU time \(C_U\ge 0\) and a grid \(\mathscr {G}_N\) be given with

$$\begin{aligned} t_f - t_0 \ge 2\, C_U({n_\omega }-1). \end{aligned}$$

Then, the objective bound for (CIA-U) mentioned in Proposition 1 is the tightest possible bound.

Proof

Let us first consider \(C_U>0\) and construct an example with the desired objective value by means of a minimal \(C_1\)-overlapping grid, where we set \(C_1=C_U\). The proposition assumes a time horizon length of at least \(2 C_1 ({n_\omega }-1 )\) so that the grid induced by Lemma 1 consists of at least \(n_b \ge n_{\omega }-1\) blocks. Let \(\mathbf {a}\in A\) be given as

Consequently, all controls \(i\in [{n_\omega }]\), \(i\ne 1\), assume in \(\mathbf {a}\) the same values on each interval. After the first block, we set control \(i=1\) to zero, while all other variables are set to \(\frac{1}{n_{\omega }-1}\) for the remaining intervals, i.e. blocks. Next, we discuss how the optimal solution of (CIA-U) on the first \({n_\omega }-1\) blocks may be chosen. Let us calculate the control deviation if we were to activate a control \(i= 2 \ldots n_{\omega }\) on the first block:

$$\begin{aligned} \varTheta _{i,1}&= \left| \sum \limits _{j\in {\mathscr {J}}_1}\frac{1}{2n_{\omega }-2}\varDelta _j - \mathscr {L}_1 \right| = \frac{2n_{\omega }-3}{2n_{\omega }-2} \mathscr {L}_1 = \frac{2n_{\omega }-3}{2n_{\omega }-2}\, (2 C_U - \varepsilon ) \\&= \frac{2n_{\omega }-3}{2n_{\omega }-2}\, ( C_U +\bar{\varDelta } - \varepsilon ) . \end{aligned}$$

In the second and third equality we used Lemma 1. The values of the relaxed controls \(\mathbf {a}\) for \(i=2\ldots n_{\omega }\) and blocks \(1,\ldots ,{n_\omega }-1\) sum up to

$$\begin{aligned} \sum \limits _{b\in [{n_\omega }-1]} \sum \limits _{j\in {\mathscr {J}}_b}a_{i,j} \varDelta _j&= \frac{1}{2n_{\omega }-2} \mathscr {L}_1 + \sum \limits _{b=2,\ldots ,{n_\omega }-1} \frac{1}{n_{\omega }-1} \mathscr {L}_b\\&= \frac{1}{2n_{\omega }-2} ( C_U +\bar{\varDelta } - \varepsilon ) + \sum \limits _{b=2,\ldots ,{n_\omega }-1} \frac{1}{n_{\omega }-1} ( C_U +\bar{\varDelta } - \varepsilon )\\&=\frac{2n_{\omega }-3}{2n_{\omega }-2} ( C_U +\bar{\varDelta } - \varepsilon ). \end{aligned}$$

Thus, there are \({n_\omega }-1\) controls with this control accumulation on \({n_\omega }- 1\) blocks, but activating any of these controls on the first block yields already the same control deviation. Hence, the objective of (CIA-U) with this \(\mathbf {a}\) is at least \(\frac{2n_{\omega }-3}{2n_{\omega }-2} ( C_U +\bar{\varDelta } - \varepsilon )\), where \(\varepsilon \) is arbitrarily small. If we combine this result with Proposition 1, we get that \(\frac{2n_{\omega }-3}{2n_{\omega }-2} ( C_U +\bar{\varDelta } )\) is the tightest possible bound. Last, we argue for the degenerate case, \(C_U=0\), that we can create an example with length of all blocks set to \(\bar{\varDelta }\) and obtain the same tight bound. \(\square \)

Corollary 1

(Tight bound on the rounding gap for (CIA)) Consider \(\mathscr {G}_N\) and \(\mathbf {a}\in A\). The objective of (CIA) is bounded by

$$\begin{aligned} \theta ^\star \ \le \ \frac{2n_{\omega }-3}{2n_{\omega }-2} \ \bar{\varDelta }. \end{aligned}$$

If \(N\ge {n_\omega }-1\) holds, then this bound is tight.

Proof

The bound follows from Proposition 1 with \(C_U=0\) and if \(N\ge {n_\omega }-1\), we are able to construct the same worst-case example as in the proof of Proposition 2, with intervals applied instead of blocks. \(\square \)

5.3 Implications for the Objectives of (CIA-D) and (CIA-UD)

The bound obtained for (CIA-U) can be transferred in a straightforward manner to (CIA-D) by using \(C_1=C_D\) as MDT in the DNFR scheme. However, we notice the increased number of degrees of freedom when dealing with MD times rather than MU times: only the down time forbidden control is fixed for a certain time duration in comparison with the MU time constraint situation where all controls are fixed due to the fixed active control. With this observation, we introduced in the DNFR scheme the min down mode \(\chi _D=1\) and are going to deduce in the sequel an alternative upper bound compared to the one obtained by DNFR with \(\chi _D=0\). As we are going to detect, this alternative bound is independent of \({n_\omega }\) but not always an improvement, so that we declare the minimum of both bounds as the upper bound in the following proposition.

Proposition 3

(Bounds on the objective of (CIA-D) and (CIA-UD)) Consider any grid \(\mathscr {G}_N\) and \(\mathbf {a}\in A\). Let MU and MD times \(C_U, C_D \ge 0\) be given. Then

  1. 1.

    (CIA-D) is bounded by

    $$\begin{aligned} \theta ^\star \le \min \left\{ \frac{3}{4} C_D + \frac{3}{2} \bar{\varDelta }, \frac{2{n_\omega }- 3}{2{n_\omega }- 2} \left( C_D + \bar{\varDelta } \right) \right\} . \end{aligned}$$
  2. 2.

    (CIA-UD) is bounded by

    $$\begin{aligned} \theta ^\star \le \ {\left\{ \begin{array}{ll} \frac{2{n_\omega }- 3}{2{n_\omega }- 2} \left( C_U + \bar{\varDelta } \right) ,\qquad &{} \text {if } C_U\ge C_D,\\ \min \left\{ \frac{3}{2} C_U + \frac{3}{2} \bar{\varDelta }, \frac{2{n_\omega }- 3}{2{n_\omega }- 2} \left( C_D + \bar{\varDelta } \right) \right\} ,\qquad &{} \text {if } C_D>C_U> C_D/2,\\ \min \left\{ \frac{3}{4} C_D + \frac{3}{2} \bar{\varDelta }, \frac{2{n_\omega }- 3}{2{n_\omega }- 2} \left( C_D + \bar{\varDelta } \right) \right\} , &{} \text {if } C_D/2 \ge C_U. \end{array}\right. } \end{aligned}$$

Proof

Generally, if \(C_D> C_U\) or if MU constraints are absent, we may apply the DNFR scheme with \((C_1,C_2,\chi _D)=\left( C_D,\frac{2{n_\omega }- 3}{2{n_\omega }-2},0 \right) \), which constructs feasible solutions for (CIA-D) respectively (CIA-UD) with objective bound \(\frac{2{n_\omega }- 3}{2{n_\omega }- 2} (C_D + \bar{\varDelta })\). We are left with the other case \(\chi _D=1\):

  1. 1.

    If we set \(C_1=\frac{ 1}{ 2} C_D\), we have \(\overline{\mathscr {L}} < \frac{ 1}{ 2} C_D + \bar{\varDelta }\). With this MDT and the choice \(\chi _D=1\) the DNFR scheme constructs a feasible solution for (CIA-D). Then, by virtue of Theorem 2, case II., we deduce with \(C_2=\frac{3}{2}\) the bound \(\frac{3}{4} C_D + \frac{3}{2} \bar{\varDelta }\).

  2. 2.
    1. (a)

      If \(C_U\ge C_D\) is given, we can set \(C_1=C_U\) and all block lengths are at least as large as those of the MD time \(C_D\). Therefore, the solution constructed by DNFR with \(\chi _D=0\) and \(C_2=\frac{2{n_\omega }- 3}{2{n_\omega }- 2}\) fulfills both the MU and MD time constraint.

    2. (b)

      We set \(\chi _D=1\) and \(C_1= C_U,\ C_2=\frac{3}{ 2}\), when \(C_D>C_U> C_D/2\) is given. By this choice, the solution of DNFR fulfills a MD time of \(2 C_U\) because of

      $$\begin{aligned} 2 \underline{\mathscr {L}}> 2 C_U > 2 C_D/2 = C_D. \end{aligned}$$

      Furthermore, by setting \(C_1= C_U\) it is clear that \(\mathbf {w}^{\text {DNFR}}\) does not violate the MU time.

    3. (c)

      \(C_D/2\ge C_U\): DNFR with down time configuration \(\chi _D=1\) and \(C_1= C_D/2 \ge C_U, \ C_2=\frac{3}{ 2}\) can be executed without violating the MU time constraint. \(\square \)

Since tightness results for the problems (CIA-D) and (CIA-UD) are not as straightforward obtained as for the problem (CIA-U), we move the discussion on the quality of the bounds obtained in Proposition 3 to the “Appendix B”.

6 Sum-up rounding in the Dwell time context

SUR is computationally less expensive (\(\mathscr {O}(n_{\omega } N)\)) than the DNFR scheme executed on intervals. But, on the other hand, the last section showed DNFR constructs solutions for (CIA) with a maximum integrality gap that is less than the one of solutions constructed by SUR for \(n_{\omega }\ge 3\) (equality for the case \(n_{\omega }=2\)). Since the SUR scheme is very often used for finding approximative solutions of (CIA), but does not necessarily fulfill MDT constraints, we discuss in this section a canonical extension of the algorithm to this setting.

6.1 Dwell time sum-up rounding (DSUR)

We introduce the concept of a currently activated control and dwell time blocks that depend on the initial interval and the MDT duration \(C_1\).

Definition 10

(Initial interval dwell time block index sets) Let an MDT \(C_1 \ge 0\) be given. We define for all intervals \(k\in [N]\) the initial interval dependent dwell time index sets to be

$$\begin{aligned} {\mathscr {J}}_{k}^{\text {SUR}}(C_1):=\{k\} \cup \{j \ \mid \ t_{j-1}\in \mathscr {G}_N \cap [t_{k-1} ,t_{k-1}+ C_1) \}. \end{aligned}$$

Definition 11

(Currently activated control) We call a control index i currently activated at interval \(j=2,\ldots ,N\), if

$$\begin{aligned} w_{i,j-1}=1 \end{aligned}$$

holds. Otherwise, or if \(j=1\), we call the binary control i currently deactivated.

In contrast to the DNFR scheme we are now interested in considering intervals individually and work independently of the constant \(\chi _D\). Hence, grouping of minimum down time forbidden controls for each interval into sets \(\mathscr {I}_j^{\text {SUR}}\) is necessary.

Definition 12

(SUR minimum down time forbidden control set) Let a MDT \(C_D\ge 0\) be given. We define the set of down time forbidden controls \(\mathscr {I}_j^{\text {SUR}} \subset [{n_\omega }] \) on interval \(j\in [N]\) as follows:

$$\begin{aligned} \mathscr {I}_j^{\text {SUR}} :=\{ i \in [{n_\omega }] \, \mid \, \exists \, k < j: \ t_{j-1} \le t_{k-1}+C_D, \, t_{k-1}\in \mathscr {G}_N \ \wedge \ w_{i,k} = 1 \}. \end{aligned}$$

We say \(i \in [{n_\omega }]\) is MD time admissible on \(j\in [N]\), if \(i \notin \mathscr {I}_j^{\text {SUR}}\) holds.

Note that the above definition assumes implicitly fixed control variables for previous intervals \([N] \ni k<j\). We have \(\mathscr {I}_1^{\text {SUR}}=\emptyset \), because there are no down time forbidden controls on the first interval. Moreover, the set \(\mathscr {I}_j^{\text {SUR}}\) may contain several controls, but at most \({n_\omega }-1\).

Next, we give a definition of the DSUR scheme in Algorithm 2. It iterates over all intervals \(j\in [N]\) and selects initially the interval representing the beginning of the time horizon, where a currently activated control does not yet exist. The control dependent MDT \(C_{i}\) is updated in line 3 for each iteration inside the while loop so that \(C_{i}\) equals the maximum of MU time \(C_U\) and MD time \(C_D\) for a currently activated control, and otherwise it is set to the MU time \(C_U\). The algorithm sets \(C_{i}=C_U\) for all controls in the first while iteration. Next, in line 4, one searches for the MD time admissible control \(i^\star \) with maximum forward control deviation on the upcoming intervals covering the dwell time \(C_{i^\star }\). If it is the currently activated one, we fix this control to be active also on the current selected interval j and increase the interval index (line 5–6). Else, the control is activated on the whole dwell time block represented by its interval indices \({\mathscr {J}}_{j}^{\text {SUR}}(C_{i^\star })\) and the interval index is updated accordingly (line 8–9). Lastly, DSUR updates the set of down time forbidden controls for the next iteration in line 11. The algorithm stops as soon as the control choice has been made for the last interval N.

figure c

Clearly, \(\mathbf {w}^{\text {DSUR}}\) is a feasible solution for (CIA), because exactly one control is active per interval. It is also feasible for (CIA-U), because whenever a currently deactivated control is activated, it dwells on active for at least the duration \(C_U\) (line 8–9). The solution also satisfies MD time constraints by the definition of \(\mathscr {I}_j^{\text {SUR}}\) and altogether \(\mathbf {w}^{\text {DSUR}}\) is a feasible solution for (CIA-UD).

We transferred the main idea from the original SUR scheme to the problem setting with MDT constraints by selecting in each iteration the control with maximum forward control deviation. In the presence of MU time requirements, we need to calculate this forward accumulation for the set of next intervals with total length at least \(C_U\). For a given MD time larger than the MU time, Algorithm 2 compares the forward accumulation with length at least \(C_D\) of the currently activated control with the ones of other controls with length at least \(C_U\). The idea behind this approach is to prevent a situation in which a control gets deactivated, but is going to accumulate a large control deviation during its down time forbidden period.

Remark 4

(Run time of DSUR) Algorithm 2 is in \(\mathscr {O}(n_{\omega } N^2)\), since it sums up, on each interval and for all controls, the relaxed control values \(\mathbf {a}\) on the next dwell time induced intervals.

6.2 Rounding gap bounds for Dwell time sum-up rounding

Kirches et al.  [24] have proven the tightest possible bound on the integrality gap for the original SUR. From this we can derive implications for DSUR in the absence of MD conditions.

Theorem 3

(Tight bound for SUR integrality gap, cf. Theorem 7.1, [24]) Let \(\mathbf {w}^{\text {SUR}}\) be constructed from \(\mathbf {a}\in A\) by means of SUR for an equidistant discretization of \([t_0 , t_f ]\) and let denote by \(\theta (\mathbf {w}^{\text {SUR}})\) its (CIA) objective value. Then, the rounding gap is bounded by

$$\begin{aligned} \theta (\mathbf {w}^{\text {SUR}}) \le \bar{\varDelta }\sum \limits _{i=2}^{n_\omega }\frac{1}{i}, \end{aligned}$$

which is the tightest possible upper bound.

Corollary 2

(Tight bound for DSUR integrality gap without MD times) Let \(C_D < \underline{\varDelta }\) and an MU time \(C_U >0 \) be given. Let the time horizon \([t_0 , t_f ]\) with a minimal \(C_U\)-overlapping grid be discretized and let \(\mathbf {w}^{\text {DSUR}}\) be constructed from \(\mathbf {a}\in A\) by means of DSUR. Then, the rounding gap \(\theta (\mathbf {w}^{\text {DSUR}})\) of its (CIA) objective value is bounded by

$$\begin{aligned} \theta (\mathbf {w}^{\text {DSUR}}) \le (C_U+\bar{\varDelta }) \sum \limits _{i=2}^{n_\omega }\frac{1}{i}, \end{aligned}$$

which is the tightest possible upper bound.

Proof

As in the proof of Theorem 3 in [24] a dynamic programming argument can be applied, here with equidistant block length of \((C_U+ \bar{\varDelta } - \varepsilon )\) as derived in Lemma 1. With a time horizon length of \( {n_\omega }(C_U+\bar{\varDelta }-\varepsilon )\) we may analogously to the proof of Theorem 3 construct an example indicating that the bound is tight as follows.

$$\begin{aligned} a_{i,j} = {\left\{ \begin{array}{ll} 0, \qquad &{} \text {if } 2i+1 \le j \le N, \\ 1/({n_\omega }+ 1 - j/2),\qquad &{} \text {if } j \text { is even},\\ 1/({n_\omega }+ 1 - (j+1)/2), \qquad &{} \text {if } j \text { is odd}, \end{array}\right. } \qquad 1\le j \le N = 2{n_\omega }. \end{aligned}$$

The DSUR scheme constructs for this example a solution that switches directly after each block with length \((C_U+ \bar{\varDelta } - \varepsilon )\). Moreover, the controls \(i\in [{n_\omega }-1]\) are active each on block i so that the last control \({n_\omega }\) accumulates the asserted rounding gap until the end of block \({n_\omega }-1\). \(\square \)

Remark 5

(DSUR as generalization of SUR) The last corollary implicitly states that DSUR can be seen as a generalization of the original SUR algorithm, since it reduces to the latter for a negligible MDT \(C_U,C_D \le \underline{\varDelta }\).

Theorem 3 allows no direct conclusion for the case with absent MU times and an active MD time \(C_D> \underline{\varDelta }\). At least, it is possible to provide worst-case examples for \(\mathbf {a}\) in order to get a clue of how large the upper bound can be for the DSUR rounding gap without MU times. This expresses the following proposition.

Proposition 4

(Rounding gap for DSUR without MU times) Consider an inactive MU time constraint, i.e., \(C_U \le \underline{\varDelta }\) and an equidistant grid \(\mathscr {G}_N\). We assume for the MD time

$$\begin{aligned} C_D > (2({n_\omega }-1)-1)\bar{\varDelta }. \end{aligned}$$
(6.1)

Let for the grid hold

$$\begin{aligned} N\ge ({n_\omega }-1)(1+M_D) + \lceil M_D/2 \rceil -1, \end{aligned}$$
(6.2)

where \(M_D\) denotes the number of minimum down time intervals constructed by \(C_D\), i.e. \(M_D :=\lceil C_D/\bar{\varDelta } \rceil \). Then, there is an \(\mathbf {a}\in A\) that yields a (CIA-D) objective value \(\theta (\mathbf {w}^{\text {DSUR}})\) of \(\mathbf {w}^{\text {DSUR}}\) constructed by DSUR with

$$\begin{aligned} \theta (\mathbf {w}^{\text {DSUR}}) \ge \left( \frac{M_D}{2} + ({n_\omega }-2)\right) \bar{\varDelta }. \end{aligned}$$
(6.3)

Proof

See “Appendix C”. \(\square \)

Remark 6

(Rounding gap for DSUR with MU and MD constraints) Generally, when the problem setting involves both MU and MD time constraints, i.e., \(C_D,C_U>\underline{\varDelta }\), the worst-case rounding gap constructed by the DSUR scheme is at least the maximum of the bounds obtained in Corollary 2 and in Proposition 4.

7 Computational Experiments

We consider a three tank flow system problem with three controlling modes in order to evaluate the integrality gap in practice and to test the proposed rounding methods. It models the dynamics of an upper, middle and lower level tank, connected to each other with pipes. The goal is to minimize the deviation of certain fluid levels \(k_2\), \(k_4\) in the middle, respectively lower, level tank. This problem type was discussed in a variety of publications for the optimal control of constrained switched systems [11, 36] and is taken from the benchmark https://mintOC.de library [33]. The problem reads

$$\begin{aligned} \begin{aligned} \min \limits _{\mathbf {x},\, \omega }&\int _0^T k_1(x_2(s)-k_2)^2+k_3(x_3(s)-k_4)^2 \, \mathop {}\!\mathrm {d}s \\&\begin{aligned} \text {s.t.}\ \dot{x}_1(t)&= -\sqrt{x_1(t)}+c_1\omega _1(t)+c_2\omega _2(t)-\omega _3(t)\sqrt{c_3x_1(t)}, \quad&\text {for a.e.}\ t \in [0,T], \\ \dot{x}_2(t)&= \sqrt{x_1(t)} - \sqrt{x_2(t)}, \qquad&\text {for a.e.}\ t \in [0,T], \\ \dot{x}_3(t)&= \sqrt{x_2(t)} - \sqrt{x_3(t)} +\omega _3(t)\sqrt{c_3x_1(t)}, \qquad&\text {for a.e.}\ t \in [0,T], \\ \mathbf {x}(0)&= \mathbf {x_0},&\\ 1&= \sum _{i=1}^{3}\omega _i(t), \ \ \omega (t)\in \{0,1\}^3, \qquad&\text {for a.e.}\ t \in [0,T]. \end{aligned} \end{aligned} \end{aligned}$$
(P)

The additional parameters are

$$\begin{aligned} k:=(2,3,1,3)^T, \qquad c:=(1,2,0.8)^T, \qquad \mathbf {x_0}:=(2,2,2)^T, \qquad T:=12. \end{aligned}$$

Furthermore, we add the MU and MD time constraints (2.1)–(2.2) to the three tank problem with varying \(C_U\) and \(C_D\) parameters. We leave open the question of the regularity of the differential states \(\mathbf {x}\), but we assume that there exists a unique solution that is continuous.

We solve this test problem with the CIA decomposition. We applied Direct Multiple Shooting for temporal discretization with varying number of grid intervals N together with a fourth order Runge-Kutta scheme for obtaining the differential state’s evolution and thus the objective value.Footnote 3 We used CasADi v3.4.5 [3] to parse the NLP with efficient derivative calculation of Jacobians and Hessians to the solver IPOPT 3.12.3 [38]. Here, the helper classes Opti stack are useful, as they allow a compact syntax for NLP modeling. For finding the optimal solution of the resulting (CIA) problem and its MDT variants we used the branch and bound solver of the open source software package pycombinaFootnote 4 [8]. We published the python source code of solving (P) via the CIA decomposition online.Footnote 5

We stress that the obtained feasible solutions for (P) via the CIA decomposition are in general not global optimal solutions. In fact, the Problem (P) appears to be nonconvex so that IPOPT may construct a local solution, just like the rounding via (CIA) may do. Nevertheless, finding a global optimal solution is computational expensive as argued in the introduction.

Figure 5 depicts the state and control trajectories with relaxed (DRCP) and binary (DBCP) control values and a required MU time of \(C_U = 0.3\). We remark that the binary solution’s objective value under MU time constraints is about 1.3% larger than the one obtained by the relaxed solution and about 1.2% larger than the binary solution’s objective value without MU time constraints.

Fig. 5
figure 5

Differential state and control trajectories for the test problem (P) with relaxed binary controls, i.e. for problem (DRCP), on the left and approximated binary controls, i.e. for problem (DBCP), on the right with MU time \(C_U=0.3\) and a temporal discretization with \(N=1280\) intervals. The objective value for (DRCP) amounts to \(\varPhi =8.776\), while the one of (DBCP) is \(\varPhi =8.888\)

Full size image

We illustrate in Fig. 6 the effect of an increasing MU time on the objective values of (CIA-U) and (DBCP). As expected, the finer the discretization grid and the smaller the required MDT time, the better the objective values of both problems. A small MU time results in a weak restriction for (DBCP) so that its objective value is close to the one of (DRCP), which is \(\varPhi =8.776\). But, from \(C_U=0.7\) on, a refinement of the grid cannot compensate the MU time restriction anymore and the (DBCP) objective value is about 25% larger than (DRCP) in this case. Interestingly, this objective value increases hardly in \(C_U\ge 0.7\), but it even decreases slightly after \(C_U=0.7\) before increasing again and staying constant from \(C_U\approx 2.0\) on. We observe a few outlier instances, e.g., \(N=20\) with \(C_U=1.2\) or \(N=40\) with \(C_U\in \{0.4,0.5\}\), where the objective value appears to be unexpectedly large. This can be explained by the coarse grid choices, which in turn stresses the importance of a fine time discretization regarding the stability of the obtained solution for (DBCP).

On the other hand, (CIA-U) objective’s value increases roughly linearly in \(C_U\) on fine grids until it reaches \(\theta ^\star \approx 0.87\), which seems to be the maximum value for P in this setting. Thus, while small values of the (CIA-U)’s objective correspond to promising objective values of (DBCP), the relationship of (CIA-U) and (DBCP) appears to be quite uncorrelated from \(C_U\ge 0.7\) on. We computed similar results for (P) with MD time constraints (not shown). We also tested whether including the relaxed MDT constraints into the NLP or not has a significant impact on the solution - this was not the case.

Fig. 6
figure 6

Objective values of (CIA-U) and (DBCP) based on the test problem (P) and different control discretizations N and MU time durations \(C_U\)

Full size image
Fig. 7
figure 7

CIA objective function values \(\theta \) for test problem (P) with time discretization \(N=1280\) and varying MU time \(C_U\) (left), respectively varying MD time \(C_D\) (right). The optimal solutions denoted with (CIA-U), respectively (CIA-D), are obtained via pycombina’s branch and bound algorithm and compared with the solutions constructed by DNFR and DSUR. We also show the upper bound (UB) for (CIA-U) respectively (CIA-D) from Propositions 1–3 and the bounds derived for DSUR from Corollary 2 and Proposition 4. We note that although Proposition 4 derives only a lower bound for the upper bound of DSUR with MD time constraints, this bound is not violated by the computational results

Full size image

We analyze the performance of DNFR and DSUR for both MU and MD time constraints and with respect to \(\theta ^\star \) in Fig. 7. The obtained solutions are compared with the global minima for (CIA-UD) from pycombina. We observe that DNFR seems to perform better for MU time constraints, while DSUR performs better for the instances with MD time requirements. We plotted also the theoretical upper bound (UB) from Propositions 1 and 3, which are here \(\frac{3}{4}(C_1+\bar{\varDelta })\), \(C_1=C_U,C_D\). As already observed for Fig. 6, the minima of (CIA-U) and (CIA-D) do hardly increase, if at all, for large MDTs and therefore diverge compared with the theoretical upper bound. We explain this behavior by the problem specific given relaxed values, which induce here an objective value of \(\theta ^\star \approx 0.87\) for (CIA-U) and (CIA-D) even if no switches are used in the binary solution.

We show also the upper bound derived for DSUR with \({\hbox {MU}}\) time constraints from Corollary 2, i.e. \(\frac{5}{6}(C_U+\bar{\varDelta })\), and the lower bound for the upper bound for DSUR with \({\hbox {MD}}\) time constraints by Proposition 4, i.e. \(\frac{1}{2}(C_D+\bar{\varDelta })\). While the solutions constructed by DSUR may violate the upper bounds for (CIA-U) and (CIA-D), as happening for the MU time case, the bounds for DSUR are not violated. We observe that the (CIA-D) objective values are not only by far smaller than their upper bound, but even smaller or equal to the DSUR bound by Proposition 4.

We implemented DNFR and DSUR in Python 3.6 as additional solvers for pycombina. The execution of the heuristics took for all instances not more than 0.02 seconds on a workstation with 4 Intel i5-4210U CPUs (1.7 GHz) and 7.7 GB RAM. We conclude that the heuristics can be used to quickly generate robust solutions in terms of competitive objective values. They might also be useful for initializing pycombina’s branch and bound algorithm with a good upper bound. However, a numerical study is needed for verifying the added benefit, which could be elaborated in future work.

8 Conclusions

In this article, we have derived tight bounds for the integrality gap of the CIA decomposition applied to MIOCPs under MDT constraints. The presented proofs are constructive and take advantage of the introduced DNFR scheme. Numerical experiments show that the CIA decomposition performs notably well in terms of objective quality for problems with small MDT requirements, since the deviation from its relaxed solution is negligible. For more restrictive dwell time constraints, the algorithm may provide feasible solutions that differ a lot from the relaxed solution. Hence, the constructed solution quality might be low compared to the optimal solution. Nevertheless, when considering the runtime of MINLP solvers, the CIA decomposition solution is computationally inexpensive and can be easily assessed with the relaxed solution for its quality. We have extended the SUR scheme so that it constructs dwell time feasible solutions and tested the proposed algorithm on a benchmark problem. The resulting (CIA) problem solutions of this and the DNFR scheme are close to the optimal ones. Hence, we propose that the DNFR or DSUR scheme may be beneficial in the context of huge (CIA) problems or in the setting of model predictive control, where a branch and bound algorithm struggles to find the optimal solution efficiently.