Abstract
This paper considers a linear-quadratic optimal control problem where the control function appears linearly and takes values in a hypercube. It is assumed that the optimal controls are of purely bang–bang type and that the switching function, associated with the problem, exhibits a suitable growth around its zeros. The authors introduce a scheme for the discretization of the problem that doubles the rate of convergence of the Euler’s scheme. The proof of the accuracy estimate employs some recently obtained results concerning the stability of the optimal solutions with respect to disturbances.
Similar content being viewed by others
1 Introduction
Discretization schemes for optimal control problems have been largely investigated in the last 60 years (see, e.g., [6,7,8,9, 17, 21], and the more recent paper [5] and the references therein). In the aforementioned papers and in most of the literature, the optimal controls are typically assumed to be sufficiently smooth (at least Lipschitz continuous) and results are usually based on second-order optimality conditions. On the other hand, whenever the control appears linearly in the system, the lack of coercivity typically leads to discontinuities of the optimal controls.
Recently, new second-order optimality conditions for systems that are linear with respect to the control have been developed. We refer to [3, 22] for analysis of second-order necessary conditions for bang–bang and singular-bang controls, respectively. Results on the stability of solutions with respect to disturbances were also recently obtained, see [4, 12, 14, 25] and the bibliography therein. Based on these results, error estimates for the accuracy of the Euler discretization scheme applied to various classes of affine optimal control problems were obtained in [1, 2, 13, 18, 26, 27]. The error estimates are at most of first order with respect to the discretization step, which is natural in view of the discontinuity of the optimal control. For the same reason, using higher order Runge–Kutta discretization schemes on a fixed grid does not help to improve the order of accuracy. Seemingly, the first paper that addresses the issue of accuracy of discrete approximations for affine problems is [30], where a higher order Runge–Kutta scheme is applied to a linear system, but the error estimate is of first order or less.
A new type of discretization scheme was recently presented in [23] for Mayer’s problems for linear systems. The idea behind this scheme goes back to [15, 20, 28] and is based on a truncated Volterra–Fliess-type expansion of the state and adjoint equations. The analysis of the convergence and of the error estimate makes use of the strong Hölder metric sub-regularity of the map associated with the Pontryagin maximum principle, proved in [25].
The goal of the present paper is to extend this discretization scheme and the pertaining error analysis to affine linear-quadratic problems. This extension is not a routine task, due to the appearance of the state function in the associated with the problem switching function, and of both the state and the control, in the adjoint equation.
More precisely, we consider the problem
subject to
Here [0, T] is a fixed time horizon, the state x is n-dimensional, the initial state \(x_0\) is given, the control u is m-dimensional, \(Q\in {\mathbb {R}}^{n\times n}, q\in {\mathbb {R}}^n\), and the matrix functions \(A,\, W:[0,T]\rightarrow {\mathbb {R}}^{n\times n}\) and \(S,\, B:[0,T]\rightarrow {\mathbb {R}}^{n\times m}\) are given data; the superscript “\({}^\top \)” denotes transposition. Admissible controls are all measurable functions with values in the set U for a.e. \(t \in [0,T]\). Linear terms are not included in the integrand in (1.1), since they can be shifted in a standard way into the differential equation (1.2).
The optimal controls in the problem (1.1)–(1.3) are typically concatenations of bang–bang and singular paths. In this paper, we assume that the optimal controls are of strictly bang–bang type with a finite number of switches, and the components of the switching function have a certain growth rate at their zeros, characterized by a number \(\kappa \ge 1\). This number appears in the error estimate obtained in this paper for the proposed discretization scheme. “Generically”, \(\kappa = 1\), and in this case the error estimate is of second order.
The paper is organized as follows. In Sect. 2 we recall some notations and formulate the assumptions. In Sect. 3 we introduce our discretization scheme and present the main result—the error estimate. Section 4 contains the proof. Section 5 presents an error estimate in case of inexact solution of the discretized problem. A numerical experiment confirming the theoretical findings is given in Sect. 6. Concluding remarks and further perspectives are discussed in Sect. 7.
2 Preliminaries
First of all we pose some assumptions, which are standard in the context of problem (1.1)–(1.3).
Assumption (A1)
The matrix functions A(t), B(t), W(t) and S(t), \(t \in [0,T]\), have Lipschitz continuous first derivatives, Q and W(t) are symmetric.
The set of all admissible controls will be denoted by \({\mathcal {U}}\subset L^\infty \). A pair (x, u) formed by an admissible control u and the corresponding solution x of (1.2) is referred to as an admissible process, and the set of all admissible processes is denoted by \({\mathcal {F}}\). The set \({\mathcal {F}}\) will be considered as a subset of the space \(W^{1,1}_{x_0} \times L^{1}\), where \(W^{1,1}_{x_0} = W^{1,1}_{x_0}(0,T)\) is the (affine) space of all absolutely continuous functions \(x:[0,T]\rightarrow {\mathbb {R}}^n\) with \(x(0) = x_0\), and \(L^{1}= L^{1}(0,T)\) has the usual meaning.Footnote 1
Due to the compactness and the convexity of U, the set \({\mathcal {F}}\) is compact with respect to the \(L^2\)-weak topology for u and the uniform norm for x. Thus a minimizer \((\hat{x},\hat{u})\) does exist in the space \(W^{1,1}_{x_0}\times L^{1}\) (in fact, also in \(W^{1,\infty }_{x_0} \times L^{\infty }\)). The following assumption requires a sort of “directional convexity” of the objective functional J at \((\hat{x}, \hat{u})\).
Assumption (A2)
Let \((\hat{x},\hat{u})\) be an optimal process in problem (1.1)–(1.3). According to the Pontryagin (minimum) principle, there exists \(\hat{p}\in W^{1,\infty }\) such that \((\hat{x},\hat{u},\hat{p})\) satisfies the following system of generalized equations: for a.e. \(t\in [0,T]\),
where \(N_U(u)\) is the normal cone to U at u:
System (2.1)–(2.4) can be shortly rewritten as
where F is the set-valued map defined by
The mapping F is considered as acting from the space \({\mathcal {X}}\) to the space \({\mathcal {Y}}\), where
The norms in these spaces are defined as usual: for \((x,p,u) \in {\mathcal {X}}\) and \((\xi ,\pi ,\rho ,\nu ) \in {\mathcal {Y}}\),
and
where \(\Vert x\Vert _{1,1}\) abbreviates \(\Vert x\Vert _{W^{1,1}}\) and \(\Vert \cdot \Vert _{s}\) is the norm in \(L^s\), \(s \in \{1, \infty \}\). Notice that the normal cone \(N_{\mathcal {U}}(u)\) to the set \({\mathcal {U}}\subset L^1\) at \(u(\cdot ) \in {\mathcal {U}}\) has the point-wise representation \(\lbrace \xi \in L^\infty : \xi (t) \in N_U(u(t)) \text{ a.e. } \text{ on } [0,T]\rbrace \).
We also define the distance
in the space \(L^\infty \). We mention that \({\mathcal {U}}\subset L^\infty \) is a complete metric space with respect to this metric ([11, Lemma 7.2]).
Observe that the inclusion (2.3) is equivalent to \(u(t)\in \mathop {\text{ Argmin }}\limits _{w\in U} \sigma (t)^{\top }w\), where \(\sigma :[0,T]\rightarrow {\mathbb {R}}^m\) is the so-called switching function, defined for all \(t \in [0,T]\) as
Thus, for \(j=1,\ldots , m\),
where \(\sigma ^j\) and \(\hat{u}^j\) stay for the j-component of \(\sigma \) and \(\hat{u}\), respectively.
Assumption (B)
(Strict bang–bang property) There exist positive real numbers \(\kappa \ge 1\), \(m_0\) and \(\delta \) such that for every \(j=1,\ldots ,m\), and for every zero \(\tau \in [0,T]\) of \(\sigma ^j\), the inequality \(|\sigma ^j (t)| \ge m_0 | t-\tau |^{\kappa }\) holds for all \(t\in [\tau -\delta ,\tau +\delta ] \cap [0,T]\).
Remark 2.1
Clearly, Assumption (B) implies that each component \(\sigma ^j\) has a finite number of zeros in [0, T] and then each component of \(\hat{u}\) is piecewise constant with values \(-1\) and 1.
The following theorem plays a crucial role in the error analysis of the discretization scheme presented below. It is a modification (under weaker conditions and a different space setting) of [4, Theorem 8], and is proved in [24].
Theorem 2.2
Let Assumption (A1) and (A2) be fulfilled. Let \((\hat{x},\hat{p},\hat{u})\) be a solution of generalized equation (2.5) and let Assumption (B) be fulfilled with some real number \(\kappa \ge 1\). Then for any \(b > 0\) there exists \(c>0\) such that for every \(y := (\xi ,\pi ,\rho ,\nu ) \in \mathcal {Y} \) with \(\Vert y \Vert \le b\), there exists a triple \((x,p,u)\in \mathcal {X}\) solving \(y \in F(x,p,u)\) and every such triple satisfies the inequality
We mention that the above property of the mapping F and the reference point \((\hat{x},\hat{u},\hat{p}) \in {\mathcal {X}}\) and \(0 \in {\mathcal {Y}}\) is a stronger version (non-local with respect to (x, p, u)) of the so-called metric sub-regularity [10].
3 Discretization scheme
In this section we propose a discretization scheme for problem (1.1)–(1.3) which has a higher accuracy than the Euler scheme without a substantial increase of the numerical complexity of the discretized problem. We recall that the Euler method has already been profoundly investigated in the case where bang–bang controls appear (e.g. [1, 2, 13, 18, 26]). As mentioned in the introduction, in doing this we use an idea that originates from [15, 28] and was implemented in [23] in the case of Mayer’s problems. The approach uses second order truncated Volterra–Fliess series, as described in the next subsection.
3.1 Truncated Volterra–Fliess series
Given a natural number N, denote \(h = T/N\), and \(t_{i}= i h\) for \(i=0,\ldots ,N\). Let u be an admissible control on \([t_i, t_{i+1}]\). The solution x of (1.2) on the interval \([t_i, t_{i+1}]\) can be represented as (see [23, Section 3])
where for shortness we skip the fixed argument \(t_i\) in the appearing functions, that is, \(A:=A(t_i)\), \(B:=B(t_i)\), etc. As usual, we denote by \(O(\varepsilon )\), \(\varepsilon > 0\), any function that satisfies \(|O(\varepsilon )|/\varepsilon \le c\), where c is a “generic” constant, that is, depending only on the data of the problem (thus, independent of i and t, although \(O(h^3)\) may depend on i and t in the above context). The second-order expansion of the solution of the adjoint equation (2.2) differs from that in [23, Section 3] due to the presence of an integral term in the objective functional (1.1), therefore we derive it below. For all \(t\in [t_i,t_{i+1}]\),
Applying the first order Taylor expansion for A, W, S at \(t_i\), the representations (3.1) of x(s) and (3.2) of p(s), and skipping all third order terms we obtain that
Hence, we obtain the following truncated Volterra–Fliess expansion for the adjoint function:
Now we shall derive a second order approximation to the integral term of the objective functional J(x, u) on \([t_i,t_{i+1}]\) in term of the first and second order momentum of the controls.
Concerning the quadratic term in the x-variable we make use of the first order part of representation (3.1) and of the Taylor expansion at \(t_i\) for W. Remembering that W is a symmetric-matrix-valued function, we obtain that
Note that an easy calculation implies that
Now we consider the mixed term in the integral in (1.1). A calculation of the same fashion as the previous one yields
Let us focus on the last term:
Integrating by parts we obtain the relation
Following [28], in order to obtain a second-order expansion expressed in term of the first and second-order momentum of u, we assume the following.
Assumption (I)
The matrix \(B^{\top }(t)S(t)\) is symmetric for all \(t\in [0,T]\).
Indeed, using Assumption (I) we obtain from the last exposed equality the expression
which can be substituted in (3.5).
Notice that Assumption (I) is always fulfilled if \(m=1\). The above obtained second order approximations will be used in the next subsection to define an appropriate discrete-time approximation of problem (1.1)–(1.3).
3.2 The numerical scheme
First of all, observe that the representation (3.1) of x(t) for \(t \in [t_i,t_{i+1}]\) depends on the control u only through the integrals \(\int _{t_i}^{t_{i+1}} u(t) {\mathrm{\,d}}t\) and \(\int _{t_i}^{t_{i+1}} (t-t_i) u(t) {\mathrm{\,d}}t\). The same applies to the approximations of the integral terms of the objective functional obtained in the last subsection. By changing the variable \(t = t_i + h s\), this pair of integrals can be represented in the form \(h z_1\) and \(h^2 z_2\), respectively, where
and \(\varphi (s) = u(t_i+hs) \) is a measurable function with values in \([-1,1]\). By varying u, hence \(\varphi \), in the set of all admissible controls on [0, T], the couple \((z_1,z_2) \in {\mathbb {R}}^{2m}\) generates a strictly-convex and compact set \(Z^m\subset {\mathbb {R}}^{2m}\). Note that \(Z^m\) can be expressed as the Cartesian product \(\prod _{1}^m Z\), where Z is the Aumann integral
As pointed out in [23], it is a matter of standard calculation to represent the set Z in the more convenient way as
where \(\phi _1(\alpha ):=\frac{1}{4}\left( -1+ 2 \alpha +\alpha ^2\right) \) and \(\phi _2(\alpha ):=\frac{1}{4}\left( 1+ 2 \alpha -\alpha ^2\right) \).
Following the hint provided by the representation (3.1), we introduce the notations
and replace the differential equation (1.2) with the discrete-time controlled dynamics
Taking into account the approximations of the objective functional in the previous subsection, we introduce its discrete-time counterpart: for \(x=(x_0,\ldots ,x_N)\), \(u=(u_0,\ldots ,u_{N-1})\), \(v=(v_0,\ldots ,v_{N-1})\),
We denote by (\(P^h\)) the discrete problem of minimizing (3.9) subject to (3.7)–(3.8). The Karush–Kuhn–Tucker theorem gives the following necessary conditions for optimality of \((x_0, \ldots , x_N)\), \((w_0,\ldots , w_{N-1})\), with \(w_i := (u_i, v_i) \in Z^m\): there is an (adjoint) sequence \((p_0,\ldots , p_N)\) such that
In order to obtain (3.12) we use again Assumption (I).
3.3 Construction of continuous-time controls and order of convergence
Let \(\{(x_i,u_i,v_i,p_i)\}\) be any solution of system (3.10)–(3.13). Based on the sequence \(\{(u_i,v_i)\}_{i=0}^{N-1}\) we shall construct a continuous-time admissible control u such that
The construction is by idea similar to that in [23] with the essential difference that now u takes values only in the set \(\{-1,1\}\) and the construction is simpler.
For \((\alpha , \beta ) \in Z\), with \(\alpha \not = -1\) (that is, \(\alpha \in (-1, 1]\)) and \(\beta \in [\phi _1(\alpha ), \phi _2(\alpha )]\) define
For \(\alpha = -1\) we set \(\tau = \theta =0\). (Given that \(\beta \in [\phi _1(\alpha ), \phi _2(\alpha )]\), this is, in fact, an extension by continuity.) Clearly, \(\tau \le \theta \), while \(\tau \ge 0\) is implied by \(\beta \ge \phi _1(\alpha )\) and \(\theta \le 1\) is implied by \(\beta \le \phi _2(\alpha )\). Then define the admissible control u component-wise as follows: for \(j = 1, \ldots , m\) and \(i = 0, \ldots , N-1\) set \(\tau _i^j = \tau (u^j_i,v^j_i)\), \(\theta _i^j = \theta (u^j_i,v^j_i)\), and
The functions \(\tau (\cdot )\) and \(\theta (\cdot )\) are defined in such a way that the relations (3.14) are fulfilled. To show this, it is enough to substitute the above defined \(u(\cdot )\) in (3.14) and calculate the integrals. We skip this trivial but cumbersome calculation.
We mention that in our framework the pairs \((u_i^j, v_i^j)\), \(i=0, \ldots , N-1\), \(j = 1, \ldots , M\), typically belong to the boundary of the set Z. In such a case every component of the control u defined in (3.15) has at most one switching point per mesh interval \([t_i,t_{i+1}]\) and we can distinguish the following possibilities:
-
(i)
if \(u_i^j = -1\) or \(u_i^j = 1\), then \(u^j(t) = -1\), respectively \(u^j(t) = 1\), in \([t_i,t_{i+1}]\);
-
(ii-a)
if \(u_i^j \in (-1,1)\) and \(v_i^j = \phi _1(u_i^j)\) then \(\tau _i^j = 0\), \(\theta _i^j = (1+u_i^j)/2\), thus
$$\begin{aligned} u^{j}(t) := \left\{ \begin{array}{ll} \, \;\;1 &{}\quad \text{ if } t \in [t_i, t_i + h \theta ^j_i],\\ -1 &{} \quad \text{ if } t \in (t_i + h \theta ^j_i, t_{i+1}]; \end{array} \right. \end{aligned}$$(3.16) -
(ii-b)
if \(u_i^j \in (-1,1)\) and \(v_i^j = \phi _2(u_i^j)\) then \(\tau _i^j = (1-u_i^j)/2\), \(\theta _i^j = 1\), thus
$$\begin{aligned} u^{j}(t) := \left\{ \begin{array}{ll} -1 &{}\quad \text{ if } t \in [t_i, t_i + h \tau ^j_i],\\ \,\;\;1 &{} \quad \text{ if } t \in (t_i + h \tau ^j_i, t_{i+1}]. \end{array} \right. \end{aligned}$$A third possibility is that \((u_i^j, v_i^j)\) happens to belong to the interior of Z:
-
(iii)
if \(u_i^j \in (-1,1)\) and \(v_i^j \in ( \phi _1(u_i^j), \, \phi _2(u_i^j) )\) then formula (3.15) has to be used to define \(u(\cdot )\).
In fact, formula (3.15) gives a unified description of all the above cases, where some of the three subintervals intervals in (3.15) degenerate in the (typical) cases (i) and (ii).
Theorem 3.1
Let Assumption (A1), (A2) and (I) be fulfilled. Let \((\hat{x}, \hat{u})\) be a solution of problem (1.1)–(1.3) for which Assumption (B) is fulfilled with some \(\kappa \ge 1\). Let \(\hat{p}\) be the corresponding adjoint function (so that \((\hat{x}, \hat{p}, \hat{u})\) satisfies the Pontryagin system (2.1)–(2.4)). Then for every natural number N system (3.10)–(3.13) has a solution \(\{(x_i,u_i,v_i,p_i)\}\). Moreover, for the continuous embedding of \((u_i,v_i)\) defined in (3.15), it holds that
We mention that, for time-invariant problems without singular arcs, Assumption (B) is typically fulfilled with a number \(\kappa \in \{1, 2, 3 , \ldots \}\), corresponding to the multiplicity of the zeros of the switching function. As argued in [23], the case \(\kappa = 1\) is in a certain sense “generic” and the error estimate (3.17) is of second order in this case. Also in the case \(\kappa > 1\) the order of accuracy is doubled in comparison with that proved in [26] for the Euler scheme. Utilization of higher order Runge–Kutta schemes on a fixed mesh could not improve the accuracy of the Euler scheme due to the discontinuity of the optimal control.
A solution \(\{(x_i,u_i,v_i,p_i)\}\) of system (3.10)–(3.13) can be obtained by any method for solving the discrete problem (3.7)–(3.9). The adjoint variables \(p_i\) do not need to be directly involved. For example, we use for numerical computations a version of the steepest descent method, where the adjoint equation is only indirectly involved for calculation of the derivative of the objective function (3.9) with respect to the variables \((u_i,v_i)\). In any case, the adjoint functions \(\hat{p}(\cdot )\) and \(\{p_i\}\) are well defined and the error estimate (3.17) is valid.
As we argue in Sect. 5, the solution \(\{(x_i,u_i,v_i,p_i)\}\) can be inexact, which leads to a modification of the error estimate as stated there.
4 Proof of Theorem 3.1
1. Preliminaries.
Let \(\{(x_i,u_i,v_i,p_i)\}\) be a solution of the discrete system (3.10)–(3.13) and let \(u(\cdot )\) be the continuous embedding of \(\{(u_i,v_i)\}\) defined in (3.15).
We embed the sequences \(\lbrace x_i \rbrace \) and \(\lbrace p_i\rbrace \) into the spaces \(W^{1,1}_{x_0}\) and \(W^{1,1}\) using the hint provided by the expansions developed in Sect. 3.1. Namely, for \(t \in [t_i,t_{i+1})\), we define
and
We show below that \(p\in W^{1,1}\). From (3.4) and (3.14) it follows that
Since
the right-hand side of (4.3) is equal to the expression of \(p_i\) given by (3.11). Thus, p is continuous at \(t_i\), and hence \(p\in W^{1,1}\). The proof that \(x\in W^{1,1}_{x_0}\) is analogous and can be found in [23, Section 5].
By Theorem 2.2, for every \(b > 0\) there exist a number c such that if \(\Vert y \Vert \le b\) then
where \(y = (\xi ,\pi ,\nu ,\rho )\) is the residual that (x, p, u) gives in (2.1)–(2.4), that is, \(y \in F(x,p,u)\). Thus we have to estimate the norm of this residual.
The estimate \(\Vert \xi \Vert _1 \le O(h^2)\) of the first residual is obtained in [23, Section 4], where the primal differential equation is the same as in the present paper. We shall analyze below the residual in the remaining Eqs. (2.2)–(2.4).
2. Residual in (2.2) and (2.4).
First, we differentiate the expression in (4.2) for \(t\in [t_i,t_{i+1})\):
Here and below \(O(t;h^2) \le C h^2\) for some constant \(C>0\) which is independent of \(t \in [0,T]\). Using (4.1) and then (4.2) we obtain that
Hence, we deduce that \(\Vert \pi \Vert _{\infty } \le O(h^2)\). Notice that the Eq. (3.12) gives that the residual in (2.4) is zero, that is, \(\nu =0\).
3. Residual in (2.3).
First of all, we derive a second order expansion of the term \(B^{\top }p+S^{\top }x\) appearing in (2.3). By (4.1), (4.2) and the Taylor expansion for B and r we have
Then, by using the definition of \(B_i\) and \(C_i\) we obtain that
Since Assumption (I) means \(B^{\top }S=S^{\top }B\),
Our goal is now to estimate the norm \(\parallel \cdot \parallel _{\infty }\) of the residual in (2.3). Since \(N_U (u)=\prod _{j=1,\ldots ,m} N_{[-1,1]}(u^j)\), we analyze a single component j of (2.3). Moreover, (2.3) is a point-wise relation, therefore we consider it on an arbitrarily fixed interval \([t_i, t_{i+1}]\). We also mention that the set Z is the area surrounded by the two parabolas \(\beta = \phi _1(\alpha )\) and \(\beta = \phi _2(\alpha )\), where \(\phi _1(\alpha ) \le \phi _2(\alpha )\). Thus the normal cone to Z is easy to calculate and the following expression is provided in [23, Section 4]:
where \(\zeta = \mathrm{sgn}(\alpha - 2\beta )\).
We consider separately each of the cases (i), (ii) and (iii) for construction of continuous time control that appear in Sect. 3.3.
Case (i) \(u_i^j \in \lbrace -1,1\rbrace \). To be specific, let us assume that \(u_i^j = -1\), hence \(v_i^j = \phi _1(-1) = \phi _2(-1) = -1/2\). (The case \(u_i^j = 1\) is similar). The normal cone to Z at the point \((-1, - 1/2)^\top \) is [see the second line in (4.6)]
Then due to (3.12) for every \(j= 1, \ldots , m\) there exist \(\mu \ge 0\) and \(\lambda \ge 0\) such that
Observe that, for \(t\in [t_i,t_{i+1}]\),
thus the quantity above is non-negative for all \(t\in [t_i,t_i+h)\). Thus, adding up (4.7) and (4.8), the latter multiplied by \((t-t_i)/h\), we obtain that
By (3.14) the quantity above is identical to the j-th component of the right-hand side of (4.5), modulo \(O(t;h^2)\). By the fact that \(u^j(t)=-1\) in case (i), we thus obtain
Case (ii) \(u_i^j \in (-1,1)\), \(v_i^j \in \{ \phi _1(u_i^j), \, \phi _2(u_i^j)\}\). We consider the case \(v_i^j = \phi _1(u_i^j)\); the case \(v_i^j = \phi _2(u_i^j)\) can be treated similarly. The continuous-time control \(u^j(\cdot )\) is defined by (3.16), where the jump point \(\theta _i^j = (1-u_i^j)/2\). The normal cone to Z at \((u_i^j,\phi (u_i^j))\) is (see the third line in (4.6))
By (3.12), there exists \(\mu \ge 0\) such that
Observe that, from the definition of \(\theta ^j_i\), it follows that the quantity
is non-positive whenever \(t\in [t_i,t_i+h \theta _i^j)\), and non-negative whenever \(t\in [t_i+h \theta _i^j , t_{i+1}]\). Thus, adding up (4.11) and (4.12), the latter multiplied by \((t-t_i)/h\), and using the definition of \(u^j\) we obtain that
By (3.14) and (4.5), the left-hand side term of the relation above is equal to \(( B^{\top }p(t)+S(t)^{\top }x(t)+r(t))^j +O(h^2)\). This proves (4.10) in the case (ii).
Case (iii) By (3.12) and the fact that \((u_k^j,v_k^j) \in \text{ Int } Z\), we have
Then,
This and (4.5) yield (4.10) in the case (iii). We can finally conclude that \(\Vert \rho \Vert _{\infty } = O(h^2)\).
Summarizing, we have obtained that \(\Vert y\Vert \le c_1 h^2\), where \(c_1\) is independent of N. Since \(c_1 h^2 \le c_1 T^2 =: b\), Theorem 2.2 implies existence of c such that for every natural N
We know that \(x(t_i) = x_i\) and \(p(t_i) = p_i\), hence
Now we focus on the last term in the left-hand side. Since \(\hat{u}\) and u take only values \(\pm 1\), as already pointed out for instance in [25, Section 4], we have that \(d^{\#}(u,\hat{u})\le c_3 \Vert u - \hat{u} \Vert _1\) for some \(c_3 \ge 0\), and this concludes the proof.
5 Error estimate in case of inexact solution of the discrete problem
The estimation (3.17) in Theorem 3.1 is valid on the assumption that the discrete-time problem (3.7)–(3.9) is exactly solved. In the present section we incorporate in the error estimation possible inaccuracy in solving the discrete-time problem. The basic argument for that is identical with the one for Mayer’s problems, presented in [23, Section 5], therefore we only sketch it.
We assume that as a result of a numerical procedure for solving the mathematical programming problem (3.7)–(3.9) we have obtained an approximate solution \((\{\tilde{x}_i \}\), \(\{\tilde{p}_i \}\), \(\{\tilde{w}_i\})\) of the first order optimality (Karush–Kuhn–Tucker) system (3.10)–(3.13). This means that the relations (3.10)–(3.13) are satisfied by the sequences \((\{\tilde{x}_i \}\), \(\{\tilde{p}_i \}\), \(\{\tilde{w}_i\})\) with some residual \((\xi , \pi , \rho , \nu ) = \left( \{\xi _i\}_0^{N-1}, \{\pi _i\}_0^{N-1}, \{\rho _i\}_0^{N-1}, \nu \right) \). We measure the size or the residual by the number
Using the approximate solution \(\{\tilde{w}_i\}\), one can define an approximation, \(\tilde{u}(\cdot )\), of the optimal control \(\hat{u}\) in the same way as described in Sect. 3.3. Then the estimation (3.17) in Theorem 3.1 takes the form
The proof of this statement is not straightforward, but the argument is identical with that in [23, Section 5], therefore we do not repeat it here.
Clearly, in order to make the error arising in the proposed discretization scheme and the error in solving the discretized problem consistent, one has to solve the latter with accuracy \(\varepsilon \) proportional to \(h^2\).
6 A numerical experiment
Example 6.1
Let us consider the following optimal control problem on the plane:
subject to
with control constraint \(u \in [-1,1]\) and for \(a > 1/2\), \(b > 0\).
Here (for appropriate values of a and b) there is a unique optimal solution with a switch from \(u=-1\) to \(u=1\) at time \(\tau \) which is a solution of the equation
Moreover, \(\tau \) is a simple zero of the switching function, thus \(\kappa =1\). Taking, for example, \(a = 1\) and \(b=0.1\), the equation above becomes
and the single real solution of this equation in [0, 1] is \(\tau = 0.492487520 \) with all digits being correct. We solved system (3.10)–(3.13) for various values of the discretization step \(h = T/N\) (N is a natural number) by using a version of the gradient projection method. The computation of the approximate solutions \((\{\tilde{x}_i \}\), \(\{\tilde{p}_i \}\), \(\{\tilde{w}_i\})\) is done with accuracy (measured by the residual \(\varepsilon \)—see Sect. 5) higher than \(h^2\), so that the theoretical error estimate (5.1) is \(c h^2\). As before, \((\hat{x}, \hat{p}, \hat{u})\) is the exact solution of the considered problem, and \(\tilde{u}\) is the continuous time control constructed as described in Sect. 3.3 for the computed \(\{\tilde{w}_i\} = \{(\tilde{u}_i, \tilde{v}_i)\}\) (note that \(\tilde{u}\) depends on size of the mesh N, therefore further we use the notation \(\tilde{u}^N\) instead of \(\tilde{u}\)).
In Table 1 we report numerical results, focusing on the most critical error \(e_N =d^{\#}(\hat{u} - \tilde{u}^N)\). We also calculate the value \(e_N/h^2\), which according to (5.1) should be bounded. This is confirmed by the results.
7 Concluding remarks
In this paper we extend the analysis of the discretization scheme introduced in [23] for Mayer’s problems to embrace convex linear-quadratic problems, affine with respect to the control. The optimal controls are assumed to be purely bang–bang, with an additional assumption involving the associated “switching function”. Our discretization approach opts for solving a discrete-time optimization problem involving an additional control variable. This yields an order of convergence which doubles that of the Euler’s method. The price for that is that the discrete problem involves quadratic control constraints (instead of the box-type constraints in the original problem).
It is worth noting that the components of the optimal controls could be, in general, concatenations of bang–bang and singular arcs. This challenging case will be a subject of further investigation. It requires, among other things, a deeper analysis of the metric sub-regularity of the system of first order optimality conditions under perturbations (a step in this direction is made in [13]). Another challenging issue is to avoid Assumption (I). Assumption of this kind is present also in [29], as well as in [16], in a nonlinear context, where it requires the Lie brackets of the involved controlled vector fields to vanish. An idea in this direction is presented in [19], but it does not seem to be efficient for numerical implementation.
Notes
To avoid a confusion, we mention that the admissible controls, hence the derivative of the state function, \(\dot{x}\), as well as the derivative of the adjoint function, \(\dot{p}\), appearing below, belong to \(L^\infty \). However, we use the \(L^1\)-norms of these derivatives in most of the considerations.
References
Alt, W., Baier, R., Gerdts, M., Lempio, F.: Error bounds for Euler approximations of linear-quadratic control problems with bang–bang solutions. Numer. Algebra Control Optim. 2(3), 547–570 (2012)
Alt, W., Baier, R., Lempio, F., Gerdts, M.: Approximations of linear control problems with bang–bang solutions. Optimization 62(1), 9–32 (2013)
Aronna, M.S., Bonnans, J.F., Dmitruk, A.V., Lotito, P.A.: Quadratic order conditions for bang-singular extremals. Numer. Algebra Control Optim. 2(3), 511–546 (2012)
Alt, W., Schneider, C., Seydenschwanz, M.: Regularization and implicit Euler discretization of linear-quadratic optimal control problems with bang–bang solutions. Appl. Math. Comput. 287/288, 104–124 (2016)
Bonnans, J.F., Festa, A.: Error estimates for the Euler discretization of an optimal control problem with first-order state constraints. SIAM J. Numer. Anal. 55(2), 445–471 (2017)
Dontchev, A.L.: An a priori estimate for discrete approximations in nonlinear optimal control. SIAM J. Control Optim. 34, 1315–1328 (1996)
Dontchev, A.L., Hager, W.W.: Lipschitzian stability in nonlinear control and optimization. SIAM J. Control Optim. 31, 569–603 (1993)
Dontchev, A.L., Hager, W.W., Malanowski, K.: Error bounds for Euler approximation of a state and control constrained optimal control problem. Numer. Func. Anal. Optim. 21, 653–682 (2000)
Dontchev, A.L., Hager, W.W., Veliov, V.M.: Second-order Runge–Kutta approximations in control constrained optimal control. SIAM J. Numer. Anal. 38, 202–226 (2000)
Dontchev, A.L., Rockafellar, R.T.: Implicit Functions and Solution Mappings: A View from Variational Analysis, 2nd edn. Springer, New York (2014)
Ekeland, I.: On the variational principle. J. Math. Anal. Appl. 47, 324–353 (1974)
Felgenhauer, U.: On stability of bang–bang type controls. SIAM J. Control Optim. 41(6), 1843–1867 (2003)
Felgenhauer, U.: Discretization of semilinear bang-singular-bang control problems. Comput. Optim. Appl. 64, 295–326 (2016). doi:10.1007/s10589-015-9800-2
Felgenhauer, U., Poggolini, L., Stefani, G.: Optimality and stability result for bang–bang optimal controls with simple and double switch behavior. Control Cybern. 38(4B), 1305–1325 (2009)
Ferretti, R.: High-order approximations of linear control systems via Runge–Kutta schemes. Computing 58(4), 351–364 (1997)
Grüne, L., Kloeden, P.E.: Higher order numerical schemes for affinely controlled nonlinear systems. Numer. Math. 89, 669–690 (2001)
Hager, W.W.: Runge–Kutta methods in optimal control and the transformed adjoint system. Numer. Math. 87, 247–282 (2000)
Haunschmied, J., Pietrus, A., Veliov, V.M.: The Euler method for linear control systems revisited. In: Lirkov, I., Margenov, S., Wasniewski J. (eds.) Large-Scale Scientific Computing, Lecture Notes in Computer Science, vol. 8353, pp. 90–97. Springer, Berlin (2014)
Krastanov, M.I., Veliov, V.M.: High-order approximations of nonholonomic affine control systems. In: Lirkov, I., Margenov, S., Wasniewski, J. (eds.) Large-Scale Scientific Computing, Lecture Notes in Computer Science, vol. 5910, pp. 294–301. Springer, Berlin (2010)
Lempio, F., Veliov, V.M.: Discrete approximations of differential inclusion. Bayreuther Mathematische Schriften 54, 149–232 (1998)
Malanowski, K.: Convergence of approximations vs. regularity of solutions for convex control-constrained optimal control problems. Appl. Math. Optim. 8, 65–95 (1981)
Osmolovskii, N.P., Maurer, H.: Applications to Regular and Bang–Bang Control: Second-Order Necessary and Sufficient Conditions in Calculus of Variations and Optimal Control. SIAM, Philadelphia (2012)
Pietrus, A., Scarinci, T., Veliov, V.: High order discrete approximations to Mayer’s problems for linear system. SIAM J. Control Optim. http://orcos.tuwien.ac.at/fileadmin/t/orcos/Research_Reports/2016-04.pdf
Preininger, J., Scarinci, T., Veliov, V.M.: Metric regularity properties in bang–bang type linear-quadratic optimal control problems. Preprint. http://orcos.tuwien.ac.at/fileadmin/t/orcos/Research_Reports/2017-04.pdf
Quincampoix, M., Veliov, V.M.: Metric regularity and stability of optimal control problems for linear systems. SIAM J. Control. Optim. 51(5), 4118–4137 (2013)
Seydenschwanz, M.: Convergence results for the discrete regularization of linear-quadratic control problems with bang–bang solutions. Comput. Optim. Appl. 61(3), 731–760 (2015)
Schneider, C., Wachsmuth, G.: Regularization and discretization error estimates for optimal control of ODEs with group sparsity. ESAIM: Control Optim. Calc. Var. (2017). doi:10.1051/cocv/2017049
Veliov, V.M.: Approximations of differential inclusions by discrete inclusions. IIASA Working Paper WP-89-017 (1989)
Veliov, V.M.: On the time-discretization of control systems. SIAM J. Control Optim. 35(5), 1470–1486 (1997)
Veliov, V.M.: Error analysis of discrete approximation to bang–bang optimal control problems: the linear case. Control Cybern. 34(3), 967–982 (2005)
Acknowledgements
Open access funding provided by Austrian Science Fund (FWF).
Author information
Authors and Affiliations
Corresponding author
Additional information
T. Scarinci and V. M. Veliov were supported by the Austrian Science Foundation (FWF) under Grant No. P26640-N25. T. Scarinci was also supported by the Doctoral Programme Vienna Graduate School on Computational Optimization, funded by Austrian Science Foundation under Project No. W1260-N35.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Scarinci, T., Veliov, V.M. Higher-order numerical scheme for linear quadratic problems with bang–bang controls. Comput Optim Appl 69, 403–422 (2018). https://doi.org/10.1007/s10589-017-9948-z
Received:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10589-017-9948-z
Keywords
- Optimal control
- Numerical methods
- Bang–bang control
- Linear-quadratic optimal control problems
- Time-discretization methods