Time-Domain Decomposition for Mixed-Integer Optimal Control Problems

We consider mixed-integer optimal control problems, whose optimality conditions involve global combinatorial optimization aspects for the corresponding Hamiltonian pointwise in time. We propose a time-domain decomposition, which makes this problem class accessible for mixed-integer programming using parallel-in-time direct discretizations. The approach is based on a decomposition of the optimality system and the interpretation of the resulting subproblems as suitably chosen mixed-integer optimal control problems on subintervals in time. An iterative procedure then ensures continuity of the states at the boundaries of the subintervals via co-state information encoded in virtual controls. We prove convergence of this iterative scheme for discrete-continuous linear-quadratic problems and present numerical results both for linear-quadratic as well as nonlinear problems.


Problem Statement and Optimality Conditions
We consider optimal control problems of the form where t 0 and t f define a fixed and finite time horizon, x : [t 0 , t f ] → R n is a state function, u : [t 0 , t f ] → R m is a control function, and U is an arbitrary subset of R m . The functions χ j , ψ j : R n → R impose constraints on the initial and terminal state, respectively, and the functions ϕ 0 , ϕ f : R n → R model initial and terminal costs. The minimum is taken over all absolutely continuous functions x(·) as well as over all measurable and essentially bounded functions u(·). Such problems are numerically solved predominantly by gradient-based optimization methods. Firstly, because of their close relation to classic problems from the calculus of variations leading to multiplier-based optimality conditions such as Pontryagin's principle, which then in turn can be solved using Newton-type shooting approaches. Secondly, because of the difficulty that their discretization naturally leads to large-scale finite-dimensional problems being typically solved locally based on finite-dimensional optimality conditions. Hence, these problems are mostly studied for compact and convex constraints on the control values.
In this work, it is important to observe that we do not impose any connectivity or convexity assumptions on the control constraint set U . Thus, (1) includes integrality constraints for, e.g., modeling discrete decisions on some components of u. Such optimal control problems exhibit combinatorial aspects and received a lot of attention in recent years. Indeed, many applications require to include integer restrictions on the control value as to model, e.g., decisions such as opening or closing valves in gas network operation [21]. Similar logical restrictions occur, e.g., in autonomous driving in case of vehicles with gear shifting [13], in contact problems such as robotic multiarm transport [6], or in the control of chemical reactors in process engineering [53].
Similar to solution approaches based on hybrid extensions of the maximum principle, we consider in this work the infinite-dimensional optimality conditions of (1) for discrete U . To this end, our goal will be to find a Pontryagin-minimum of (1), i.e., an admissible state-control point (x * , u * ), such that for any constant N there exists an ε = ε(N ) > 0 with the properties that for any admissible point (x, u) satisfying see, e.g., [8,42]. Note that global minima are Pontryagin-minima.
We suppose that f : Q → R n is continuous with a continuous partial derivative f x on a set Q ⊆ R n+m and the functions χ j , ψ j , ϕ 0 , ϕ f : P → R are continuously differentiable on a set P ⊆ R n . Here, Q and P are open super-sets of the admissible state-control set and its projection onto the states, respectively.
In particular, under this assumption, the Pontryagin maximum principle together with the constraints of Problem (1) yield the following necessary optimality conditions for (x, u) being a Pontryagin-minimum: So far, this is well-known and, indeed, does not require any assumption on U in the original proof in [42] using needle-variations, the terminal cone, and the separating hyperplane theorem as, for instance, also being exploited for the so-called hybrid maximum principle [8]. It is also well-known that Problem (1) being in so-called Mayer form is not restrictive. One can, for instance, include running costs in Problem (1) by defining a new state variable y witḣ and adding y(t f ) to the objective function. Nevertheless, we note that our convergence result in Sect. 4 applies to problems in Mayer form with linear dynamics and a quadratic objective function. Hence, nonlinear terms in L(x, u) that depend on x are not covered. For running costs L that are independent of x, however, (3) is equivalent to considering . Besides this, we choose the Mayer form to keep the optimality system as simple as possible. Moreover, a free terminal time can be transformed into a fixed terminal time, etc.; see, e.g., [33]. We also note that the regularity assumptions on f and on the constraint functions can be relaxed [1]. As already indicated above, most approaches in the literature yet focused on solving (2) using further necessary optimality conditions for the pointwise maximization of the Hamiltonian. Additional properties on U are then needed. Most prominently, e.g., for convex sets U , the pointwise maximization of the Hamiltonian can be replaced using Karush-Kuhn-Tucker type optimality conditions and the two-point boundary value problem (2) can then be solved using Newton's method in the fashion of single or multiple shooting. If U is a discrete-continuous set, switching-time reparameterization can be used. This, however, leaves a gap on efficient treatments concerning the combinatorial aspects. For an overview of shooting techniques we refer to [5]. In any case, all of these approaches solve (2) to some form of local optimality with respect to the maximization of the Hamiltonian. However, the conditions in (2) hold with global optimality even for (locally optimal) Pontryagin-minima.
This motivates us to consider a solution approach for (2) based on (global) mixed-integer nonlinear programming (MINLP) techniques. This, of course, needs to overcome the difficulty that (2) is an MINLP subject to differential equations with no direct advantage over the original problem (1). Instead, (2) has even more variables for the adjoint equation. In particular, a direct discretization of (2) as, e.g., considered in [20] yields a large-scale MINLP that in many applications becomes computationally intractable.
We overcome this limitation by using time-domain decomposition methods originally developed for partial differential equations (PDEs); see, e.g., [23,24,[27][28][29][30][31][34][35][36]. The main idea of these methods is to split the time domain [t 0 , t f ] into nonoverlapping subintervals and then to iteratively decouple the optimality system (2) such that on each subinterval, subproblems are solved together with conditions at the breakpoints that couple the states over two subsequent iterations. For PDEs, this decomposition approach goes back to the so-called parareal-scheme introduced in [37]. In the following, we show that this approach can also be used to decompose (2) for arbitrary sets U . As in [29][30][31], we exploit that the decomposition yields subproblems that correspond to so-called virtual control problems on the subintervals. For discrete-continuous U , this yields an iterative scheme of mixed-integer optimal control problems on smaller time horizons. For sufficiently fine decompositions, these can eventually be solved using direct transcription methods such as collocation [22] or Runge-Kutta discretizations [49], both resulting in reasonably sized finite-dimensional mixed-integer nonlinear problems. The limit behavior of such approximations are discussed in [20]. The method can be interpreted as an alternating approach [12] or as a structure-exploiting decomposition for MINLPs rather than a generic one [41].
Our main contribution is that we prove convergence of this iterative scheme for the important case of linear-quadratic problems. We support our theoretical result with encouraging numerical results and also include experiments for nonlinear problems. In particular, we demonstrate that the iterative method can provide solutions for problems where the same solver applied to full direct discretization fails to reach an optimal solution within amply time limits.
The remainder of the paper is structured as follows. In Sect. 2, we introduce the timedomain decomposition of the optimality system and discuss the iterative idea to recover the transmission conditions at the boundaries of the subintervals. Afterward, we reinterpret the decomposed (primal-dual) optimality conditions as so-called (primal) virtual control problems in Sect. 3 and state the overall iterative procedure. For the case of linear-quadratic problems, we prove convergence of the method in Sect. 4. We present some case studies for linear as well as nonlinear problems in Sect. 5 to give a numerical proof of concept. Finally, we conclude in Sect. 6.

Time-Domain Decomposition of the Optimality System
We now consider the time-domain decomposition of the entire time horizon of Problem (1). Accordingly, we define for k = 0, . . . , K ; see Fig. 1.
The optimality conditions (2) then are equivalent to the following sets of conditions.
These conditions are completed with the conditions which ensure the continuity of the state x and of the adjoint λ at the boundaries of the sub-intervals. The main idea of the time-domain decomposition method is to decouple the transmission conditions (7) and to construct an iterative procedure that converges to points that satisfy the decomposed optimality system above. To this end, we consider an iteration index , iterates x k , λ k , and u k for x k , λ k , and u k , respectively, and a scalar parameter γ > 0, to decouple the transmission conditions (7) as with the update rules for ε ∈ (0, 1); see, e.g., [29][30][31][34][35][36]. Let us make a first observation. To this end, assume that the iterates x −1 k , λ −1 k for k = 0, . . . , K converge for → ∞ to x k and λ k , respectively. Then, substituting φ k,k+1 and φ k,k−1 and dividing by (1 − ε) yields Next, we shift the k-index in Equation (10b) by 1 and obtain Adding this to (10a) yields which is, for γ > 0, equivalent to Finally, (10a) implies Hence, we have shown that if the state x k and the adjoint λ k converge, then the iteratively decoupled transmission conditions (8) tend to the continuity conditions (7). If we combine the decomposed optimality conditions (4-6) for all k = 0, . . . , K with the iterative transmission conditions (8), we geṫ together with the update rules (9).

Virtual Control Problems and the Iterative Procedure
We will now observe that the iteratively decomposed optimality systems (11) have a primal interpretation as other optimal control problems with additional and so-called virtual controls h; see, e.g., [29][30][31][34][35][36]. For the inner time sub-intervals k, the respective primal problem reads The ODE constraint (12b) and the control constraint (12e) carry over from Problem (1) but are now restricted to the time interval [t k , t k+1 ]. The virtual controls enter in the initial conditions (12c). These controls allow the state to be bounded away from the current transmission conditions at t k . The violation of the current transmission conditions at t k and t k+1 is penalized in the objective function. Since this virtual control problem is defined for the inner time sub-intervals, there are no initial or final conditions as in Problem (1) in the constraints or in the objective function. In order to bring (12) into the form of Problem (1), one can remove (12d) and model the virtual control h k as a trivial state variable, i.e., Technically speaking, h k is not a control but a constant state. We will, however, keep referring to it as a virtual control to be in line with the pertinent literature; see, e.g., [31,34]. We use the Hamiltonian in which we omit the adjoint state corresponding to h k because it vanishes. The corresponding optimality conditions are given bẏ If we substitute h k (t k ) in (13c) with (13f) and rewrite (13g), we obtaiṅ This exactly corresponds to the iteratively decomposed optimality system (11) for k = 1, . . . , K − 1. Note that we scaled the objective function (12a) with the factor 1/(2γ ) to be in line with with the notation of System (11). However, any other positive factor would be valid as well.
For the first time sub-interval, i.e., k = 0, the respective primal problem reads in which we have no transmission conditions at t 0 but the original initial conditions instead. For the last time sub-interval, i.e., k = K , we get min in which we have no transmission conditions at t K +1 but the original final conditions instead. Now we can state the iterative time-domain decomposition method; see Algorithm 1. In Step 3 the virtual control problems (12), (14), and (15) are solved. This is equivalent to solving the iteratively decomposed optimality systems (11) for all k = 0, . . . , K . One can compute the values of the adjoint variables λ k at the transmission points t k and t k+1 by using (11i) and (11j) in Step 4. In Step 6, the update rules in (9) are used to obtain the transmission conditions for the next iteration.
To conclude this section, also note that the solution of the K + 1 problems in Step 3 of Algorithm 1 can be done in parallel.

Convergence Analysis for Linear-Quadratic Problems
In this section, we will prove the convergence of Algorithm 1 in the sense of continuity at the boundaries of the sub-intervals, i.e., we show that To this end, we restrict ourselves to the case of linear dynamics with objective functions being quadratic with respect to the state and include control costs, i.e., we consider cf.
We assume that (17) and (18) have a solution. In general, it is not to be expected that this solution is unique, which is one of the main differences to the assumptions made for similar time-domain decompositions of PDEs in, e.g., [29][30][31]. Let (x k , u k , λ k ) be a solution of (17) and let (x k , u k , λ k ) be a solution of the iteratively decomposed optimality conditions (18) for k = 0, . . . , K and ∈ N. We introduce the errors These errors satisfy the conditionsẋ together with the update rules Moreover, we setŨ To derive System (19) we use that both (17) and (18) are linear. This is the case, because Problem (16) has linear dynamics and in the objective function, the only quadratic terms in x are initial and terminal costs. There are approaches to extend this method to more general nonlinear right-hand sides f (x, u) by applying Lipschitz-type conditions on the nonlinearities and their derivatives; see, e.g., [27] where this was done for hyperbolic semilinear PDEs. This, however, exceeds the scope of this work. First, we prove the following result.
Proof From (17h) we have for all u ∈ U a. e. in (t k , t k+1 ). Analogously, using (18h) we get for all u ∈ U a. e. in (t k , t k+1 ). If we set u = u k in the last inequality, we can write a. e. in (t k , t k+1 ).

Definition 1
We define the state of all updates in iteration as Moreover, we define the mapping With this notation at hand and I denoting the identity, the update in Step 6 of Algorithm 1 is given by the relaxed mapping T ε := (1 − ε)T + ε I , i.e., Consequently, we consider the fixed point iteration Clearly, X is a fixed point of T and T ε if and only if the transmission conditions for the errorsx are satisfied. The errors fulfill these transmission conditions if and only if x k and λ k from the iteratively decomposed optimality conditions (18) do, too. This is only the case in a solution of (17). Since we assume that (17) has at least one solution, T and T ε have a fixed point.

Remark 1 For anyX
the corresponding right-hand sides of the transmission conditions (18i) and (18j) are given by where (x k , u k , λ k ) is the chosen solution of (17). If System (18) has a solution (x k ,ū k ,λ k ), then the errors (x k ,ũ k ,λ k ) = (x k ,ū k ,λ k ) − (x k , u k , λ k ) are well-defined and, thus, the mapping T : R 2K → R 2K is well-defined, too.

Definition 2
We define the energies Lemma 2 For all ∈ N, the energies E and F are non-negative, i.e., Proof It is clear that E is non-negative because it is a sum of norms. The non-negativity of F follows from Lemma 1 and of Q 0 and Q f being positive semi-definite.
Next, we show that we can use the energies to describe the state of all updates.

Lemma 3 It holds X
Proof We multiply the state equation (19a) byλ k and integrate to obtain It follows Now, we can write We get a similar result if we apply the mapping T first.

Lemma 4 It holds T X
Proof Similarly to the proof of Lemma 3 we can write

Remark 2 A direct consequence of Lemma 3 and 4 is
Since F is non-negative, it follows Hence, it is sufficient for Algorithm 1 to be well-defined that (18) has a solution for X ∈ {x ∈ R 2K : x ≤ X 1 } instead of for all vectors in R 2K .
Next, we observe that T is a non-expansive mapping.

Lemma 5 The mapping T
be the the solutions of (18) corresponding to X 1 and X 2 . The errors that solve (19) are given by We define the differences for all k = 0, . . . , K . Because of the linear nature of (19), these differences fulfill the systemẏ Next, we investigate the relation between y k (t k ), μ k (t k ) and y k (t k+1 ), μ k (t k+1 ) .
To this end, we multiply (25a) with μ k and integrate to get From (22), we have (t k , t k+1 ). Therefore, it holds Finally, we can conclude that holds because Q f and Q 0 are positive semi-definite.
We can now apply Schaefer's theorem to T .
Theorem 1 (Schaefer [48]) Let I be a convex, closed, and bounded set in a uniformly convex Banach space X and let T : I → I be a non-expansive mapping with at least one fixed point. Then, for any ε ∈ (0, 1) the mapping T ε = (1 − ε)T + ε I is asymptotically regular, i.e., This theorem allows us to prove convergence for Algorithm 1 applied to the linearquadratic problem (16). (16), the iterates (x k , λ k ) for k = 0, . . . , K converge in the sense

Theorem 2 If Algorithm 1 is applied to Problem
as → ∞. We rewrite the scalar product of T X and X as Using Lemmas 3,4,and (27) we get Because of (26), it follows for all k = 0, . . . , K − 1 as → ∞. Since this holds for the errors, it is also true for the iterates, i.e.,

Remark 3
Note that Theorem 2 states convergence of Algorithm 1 with respect to the error. In order to conclude for convergence of the iterates, further assumptions are needed. To this end, we recall that if we have continuity at the boundaries of the sub-intervals, i.e., for all k = 0, . . . , K − 1, then the optimality conditions (17) are fulfilled as well. Therefore we can conclude that the iterates (x k , λ k ) converge to (x k , λ k ) in L 2 (t 0 , t f ) for all k = 0, . . . , K if (17) has a unique solution (x k , λ k ).
Moreover, we can conclude that there is a subsequence of (x k , λ k ) that converges to a solution (x k , λ k ) of (17) if the boundary points x k (t k+1 ), x k+1 (t k+1 ), λ k (t k+1 ), λ k+1 (t k+1 ) are contained in a compact set. This is the case if (17) has more than one but still finitely many solutions.
If System (17) has infinitely many solutions, a similar conclusion can, in general, not be made. However, Algorithm 1 can still be applied in practice, where we have to impose a reasonable stopping criterion anyway such as, e.g., (28) for all k = 0, . . . , K − 1 with given tolerances δ x , δ λ > 0. This way, we can still compute an approximate solution that is arbitrarily close to a solution of (17) by choosing δ x and δ λ sufficiently small.

Case Studies
In this section, we apply Algorithm 1 to test its practical performance on some exemplary problems. We implemented the algorithm in Julia 1.5.3 [3]. 1 All computations were done on a machine with an Intel(R) Core(TM) i7-8550U CPU with 4 physical cores (and 8 logical cores), 1.8 GHz to 4.0 GHz, and 16 GB RAM. To solve the virtual control problems in Step 3, we apply a first-discretize-then-optimize approach using a Runge-Kutta method [20,49]. To this end, we equidistantly partition each time domain [t k , t k+1 ], k = 0, . . . , K , using N + 1 ∈ N time points Runge-Kutta methods approximate the solution of the ODĖ of the virtual control problems (12), (14), and (15) as Moreover, we use the initial condition (12c) to solve for the virtual control We can then substitute h k in the objective function and remove the constraint (12c). The discretized virtual control problem (12) is thus given by min The problems (14), (15) can be discretized analogously. Problem (30) is a finitedimensional optimization problem, which we model using GAMS 25.1.2 [39]. The resulting instances are solved using ANTIGONE 1.1 [40]. All virtual control problems are solved in parallel using Julia's build-in multi-threading capabilities.
In our implementation, the iteration is started with the zero transmission conditions φ 0 k,k−1 = 0 ∈ R n for k = 1, . . . , K and φ 0 k,k+1 = 0 ∈ R n for k = 0, . . . , K − 1. We initialize the algorithm with all variables set to 0 and after the first iteration, we warmstart all virtual control problems with their solution from the previous iteration. As a stopping criterion for Algorithm 1, we use (28) for all k = 0, . . . , K − 1 with tolerances δ x = δ λ = 10 −2 . In addition, we use the time limit of 1000 s for all subproblems unless stated otherwise.
After termination of the algorithm, we fix the resulting control u and compute the corresponding state x on the entire time horizon (t 0 , t f ) to obtain an accurate objective value.

A Mixed-Integer Linear-Quadratic Problem
We consider the mixed-integer linear-quadratic problem Note that this problem is of the form (16). We use the time-domain decomposition with t i = i/4 for i = 0, . . . , 4 as well as parameters γ = 1, ε = 0.5, and Δt = 1/100. Figure 2 shows the solution of Problem (31) using the direct discretization and solution approach from above on a single interval  Table 1. While these smaller tolerances obviously increase the running time of the algorithm, they have a positive impact on the objective value. This impact increases if more domains are used. Since (31) is a mixed-integer linear-quadratic problem, it is also possible to solve it with the solver Gurobi 9.1.2. However, our preliminary numerical results showed that for almost all cases, ANTIGONE performs better w.r.t. the running time.
For comparison, we include these results in the fourth part of Table 1.
To make sure that our discretization step size Δt is small enough, we used the resulting controls from the first part of Table 1 to compute the corresponding statex with the halved step size Δt/2. The subsequent errors max i=0,...,N +1 x(t i ) −x(t i ) ∞ are all smaller than 1.4 × 10 −9 which is significantly smaller than the tolerances δ x and δ λ at the interfaces.
Next, we study the impact of the parameter γ . It functions as a weighing factor of the state x compared to the adjoint λ in Step 4 of the algorithm. It is thus an important parameter regarding the errors that are made w.r.t. the continuity of the states and adjoints at the interfaces of the time domains. Figure 3 displays these errors for γ ∈ {0.2, 1, 5} and the case of 4 domains. One can clearly see that the smallest value γ = 0.2 leads to larger errors for the adjoint λ compared to the errors for the state x. For the larger value γ = 5 it is the other way around. For γ = 1, both errors have roughly the same size. Since in our case, the tolerances δ x and δ λ are the same, it takes less iterations for the algorithm to reach the stopping criterion (28) for The overall time ("Time") includes the time spent in the algorithm ("Alg. time") and the time to a posteriori compute the objective value for fixed controls γ = 1 than for the other values. More precisely, we need 35 iterations instead of 45 for γ = 0.2 or 44 for γ = 5. Figure 3 also shows the typical convergence behavior of Algorithm 1, i.e., the errors are reduced quickly during the first iterations but only decrease slowly when they are already closer to zero. This suggest that the method has an asymptotically sublinear convergence rate. However, estimates on these rates are beyond the scope of this paper. Lastly for this example, we investigate the role of the parameter ε, which weighs the boundary data of the domain k against the boundary data of its neighboring domain k−1 or k + 1. Table 2 shows the effect of different parameters ε on the number of iterations, the running time, and the objective value of Problem (31) when solved using 4 time domains. One can see that ε does not have a significant impact on the objective value but on the number of iterations and, thus, the running time. We see that we obtain a roughly monotonically increasing relation between ε and the running times of the algorithm. For this example, the fastest algorithm is obtained with ε = 0.1. We also additional sampling points for ε between 0 and 0.2, which confirm that ε = 0.1 is the best choice for this problem.

Mixed-Integer Nonlinear Problems
The convergence theory from Sect. 4 does not apply to problems with nonlinearities in the differential equation. Nevertheless, we can conclude from considering Algorithm 1 as the fixed-point iteration (23) that upon convergence, the error in the transmission conditions vanishes and a candidate for a Pontryagin minimum is found. We demonstrate this for two nonlinear examples. As a first nonlinear example we consider Fuller's initial value problem from the benchmark library mintOC of mixed-integer optimal control problems (see [45] and https://mintoc.de): Again, we choose the parameters γ = 1 and ε = 0.5 for our computations. For this example the step size for the discretization is Δt = 1/50. We solve Problem (32) with Algorithm 1 using 1, 2, 4, 8, and 16 domains. For a single domain, i.e., without using the presented time-domain decomposition approach, ANTIGONE reaches the time limit and cannot solve the problem. For two domains, our algorithm also reaches   Table 3. For 4 domains, the algorithm finds a solution in 56 s. Figure 4 shows the solution for 4 domains. For more domains, the number of iterations and the running time increases.
The objective values for the last three cases are all close to zero. This shows the potential of Algorithm 1 when applied to nonlinear problems since it can outperform applying a global MINLP solver to the discretized problem directly, i.e., on the entire time horizon. However, let us comment on that this comparison needs to be interpreted carefully since the goal of ANTIGONE as a global MINLP solver is different to the one of our method that aims to compute a Pontryagin minimum of the problem. Thus, we additionally applied ANTIGONE with a time limit of 56 s, which is the solution time of our method, and compared the feasible point found by ANTIGONE within this time limit with the point that our method computed. For the considered instance, ANTIGONE found a feasible point of comparable quality in 26 s. Note, however, that this point comes with no quality guarantee whereas we can guarantee that we computed a Pontryagin minimum. Again, we used the resulting controls from Table 3 to compute the corresponding statex with the halved step size Δt/2. The errors max i=0,...,N +1 x(t i ) −x(t i ) ∞ are all smaller than 1.9 × 10 −9 which is, again, significantly smaller than the tolerances δ x and δ λ at the interfaces.
As for the example in the linear-quadratic case we again test, which values of ε lead to the best results. Table 4 contains the results for the case of 4 domains. We skip the case of ε = 0 since we reached the time limit for this case. The algorithm also takes rather long for ε = 0.1. For ε between 0.2 and 0.7, one achieves the shortest running times. Again, if the value of ε is too large, the number of iterations and the running time increase considerably. As a second example, we consider a variation of Problem (32), namely Fuller's initial value multimode problem from mintOC (see again [45] and https://mintoc.de): s.t.ẋ 1 = x 2 a. e. in (0, 1), (33b) x 3 = x 2 1 a. e. in (0, 1), (33d) Here, we have three more binary controls and an SOS-1 constraint w.r.t. the binary controls. This makes the problem harder to solve, which is why we increased the maximal running time for this example to 2000 s. Note that u 4 does not occur in the ODE system (33b-33d). Therefore, u 4 (t) = 1 corresponds to not controlling the system at time t.
The parameters γ = 1, ε = 0.5, and Δt = 1/50 stay unchanged and we, again, solve the problem using 1, 2, 4, 8, and 16 domains. The results are displayed in Table 5 except for a single domain and the case of 16 domains. For a single domain, we again reach the time limit-which is also the case for our algorithm applied to 16 domains, where the time limit is reached in iteration 737. For this example, we achieved the fastest running time (of 209 s) for 8 domains. One can see that the number of iterations for 2 and 4 domains is lower than for 8 domains but the running times are higher. This is the case because, in a few iterations, the solution time of a virtual control problem is much higher than usual. Thus, we can see the trade-off between (i) a usually higher coordination effort of the method to obtain continuity at the interfaces of the time domains for a larger number of domains and (ii) a usually larger running time for less domains since the virtual control problems then are larger. The solution is shown in Fig. 5.
We again compared the solution that we obtain with the best feasible point computed by ANTIGONE within the time limit given by the best running time of our method, i.e., within 209 s. In this case, the feasible point found by ANTIGONE is slightly worse than the one we computed. The first feasible point found by ANTIGONE that is better than our solution needs 899 s to be computed by ANTIGONE.
For this example, we also used the resulting controls from Table 5 to compute the corresponding statex with the halved step size Δt/2. This time, the errors max i=0,...,N +1 x(t i ) −x(t i ) ∞ are all smaller than 3.1 × 10 −9 which is, again, significantly smaller than the tolerances δ x and δ λ at the interfaces.
An example for which convergence of our implementation of Algorithm 1 was not observed is a reformulation of the highly nonlinear problem "F-8 aircraft" from the mintOC library (see [45] and https://mintoc.de) to a fixed-time horizon problem, where either the time limit of 2000 s was reached or at some point one of the virtual control problems reported infeasibility. This happened for all combinations of parameters ε ∈ {0, 0.1, . . . , 0.9}, Δt ∈ {1/50, 1/100}, δ x , δ λ ∈ {10 −1 , 10 −2 } and all numbers of domains in {1, 2, 4, 8, 16} and could be caused by the presence of multiple local minima. See [15] for a successful application of a global maximum principle to the aircraft problem with numerical results.
The state x(t) ∈ R represents the position of a car with bounded acceleration u. The car starts at −7 and has an initial speed of 2. The goal is to park the car in the origin in the shortest possible time T . The optimal bang-bang solution to this problem is to accelerate (u(t) = 1) until t * = 1 and to break (u(t) = −1) after this until T = 4. To fit this in our setting we use a time transformation to reformulate the problem as where x 1 is the position, x 2 the speed, and x 3 the time T . This is a nonlinear problem with initial and terminal conditions. The optimal solution stays unchanged but the switching point between acceleration and breaking is now t * = 1/4. If we try to solve this problem with Δt = 1/100 using ANTIGONE, no feasible point is found. For 2 domains of equal length, the algorithm reaches the time limit of 1000 s in iteration 19. However, for 4 domains of equal length and the 2 domains [0, 1/4], [1/4, 1] the global optimal solution is found in 60.886 s and 182.550 s, respectively. This can be explained by both decompositions containing the switching point t * = 1/4 in an interface between two subdomains, which decreases the difficulty of the virtual control problems. Therefore, it could be an interesting topic for future research to compute the switching points of a bang-bang solution before applying Algorithm 1. This could be combined with considering constant control decisions for each subdomain, which could further decrease the difficulty of solving the virtual control problems.

Conclusion
Our results show that linear-quadratic mixed-integer optimal control problems can be solved iteratively by computing solutions to suitably chosen virtual mixed-integer optimal control problems on smaller time-horizons. In order to reach a certain accuracy of the solution, these subproblems then require less discretization variables and they can be solved in parallel. Both aspects can provide significant computational advantages compared to direct discretization of the original problem on the entire time horizon.
We exemplarily demonstrate this with our numerical results, where we also discuss the choice of additional parameters of the proposed algorithm. Our computational experiments show that this advantage also pays off for the case of nonlinear problems if the iterative procedure converges. However, a more detailed numerical study that might also include a comparison with other methods are out of scope of this paper and part of our future work. In particular, the iterative method may even yield solutions where full direct discretization fails. Guarantees for nonlinear problems, however, requires further algorithmic and theoretical investigations. The presented study is primarily motivated by the need to incorporate integer restrictions on the control values into optimal control problems as to model logical constraints. However, the findings are also interesting for classic optimal control problems with other nonconvex control constraints, because the proposed iterative procedure provides an alternative to shooting-type methods as a solution approach with its well-known limitations such as proper initialization or achieving global optimality in the Hamiltonian maximization condition.
Future work will therefore concern extensions of this approach to more general problem classes including further nonlinearities, state constraints, and also partial differential equations. Moreover, the general type of methods discussed in this paper is known to have slow convergence rates although approximate solutions of acceptable quality might be obtained rather quickly. It is thus another future research topic to consider possible crossover-approaches in which an acceptable solution found with the proposed method is handed over to another method that has favorable local convergence properties. Additionally, investigations on the convergence rates and on accelerating mechanisms would be of great interest. Finally, the method presented in this paper could be combined with techniques that try to compute the time points at which the optimal control switches, see, e.g., [9,50], and to use this information to set up the time blocks of the domain decomposition.