remarkRemark hypothesisHypothesis claimClaim Numerical procedure for hybrid
systemsR. Pytlak and D. Suski

# Numerical procedure for optimal control of hybrid systems with sliding
modes, Part I

RadosŁaw Pytlak Faculty of Mathematics and Information Science, Warsaw
University of Technology, 00-665 Warsaw, Poland<EMAIL_ADDRESS>Damian Suski Institute of Automatic Control and Robotics, Warsaw University of
Technology, 02-525 Warsaw, Poland<EMAIL_ADDRESS>
###### Abstract

This paper concerns the numerical procedure for solving hybrid optimal control
problems with sliding modes. The proposed procedure has several features which
distinguishes it from the other procedures for the problem. First of all a
sliding mode is coped with differential-algebraic equations (DAEs) and that
guarantees accurate tracking of the sliding motion surface. The second
important feature is the calculation of cost and constraints functions
gradients with the help of adjoint equations. The adjoint equations presented
in the paper take into account sliding motion and exhibit jump conditions at
transition instants. The procedure uses the discretization of system equations
by Radau IIA Runge–Kutta scheme and the evaluation of optimization functions
gradients with the help of the adjoint equations stated for discretized system
equations. In the first part of the paper we demonstrate the correspondence
between the discrete adjoint equations and the discretized version of the
continuous adjoint equations in the case of system equations described by
ODEs. We show that the discrete adjoint state trajectories converge to their
continuous counterparts in the case of ODEs.

###### keywords:

optimal control problems, hybrid systems, sliding modes, implicit Runge–Kutta
method, adjoint equations

65L05, 49M37, 65K10

##  1 Introduction

The systems with sliding modes considered in our paper state a special class
of so called hybrid systems ([18]). Many technical and non-technical systems
exhibit the hybrid discrete-continuous behavior, therefore optimal control
problems for hybrid systems were considered in many papers ([30], [2], [25],
[19], [20], [21], [26], [27], [29], [10]). Some of them discussed also
algorithms for these problems ([19], [20], [21], [26], [29]). These algorithms
refer to the strong version of the maximum principle for hybrid systems. It
means that their applicability is limited since the class of problems which
can be solved by algorithms using locally the strong maximum principle is much
smaller than the class of control problems which could be tackled by methods
which invoke locally the weak maximum principle condition.

The algorithm based on weak variations and aimed at solving optimal control
problems with hybrid systems was stated in ([16]). The numerical procedure
discussed therein calls the implicit Runge–Kutta method to evaluate system
equations and calculates gradients of the optimization problem functionals
with the help of adjoint equations. The adjoint equations are figured out on
the basis of discrete time system equations which are the result of the
numerical integration of system equations. It means that the discretization of
the adjoint equations is consistent with the discretization of the system
equations. The paper presents the results of employing the method to solving
some optimal control problems described by hybrid systems expressed by higher
index DAEs. However, the paper does not present the convergence analysis of
the procedure.

The first step towards the missing convergence analysis is done in the
accompanying paper ([14]). The optimal control problem discussed in the paper
includes, as the special case, the problem analyzed in ([16]) provided that
the hybrid system evaluates according to ordinary differential equations
(ODEs). On the other hand, the hybrid system examined in ([14]) is more
general since it includes the sliding motion. The paper ([14]) states
necessary optimality conditions in the form of the weak maximum principle for
optimal control problems with hybrid systems treated in such a general
setting. It also proposes a globally convergent algorithm for the considered
optimal control problems.

The aim of this paper is to show that a conceptual algorithm stated in ([14])
can be a basis for building its efficient implementable version. It can be
achieved by following the approach stated in ([16]). Since the system movement
in a sliding mode is represented by DAEs with index 2 the numerical procedure
from ([16]) is applicable to these hybrid systems. Therefore, the gradients of
the functions defining the optimal control problem can be evaluated using the
discrete adjoint equations consistent with the discretized system equations.
However, in contrast to the results of ([16]) we show that the gradients
elicited in that way converge to the gradients of functionals defining the
continuous time version of the optimal control problem provided that the step
sizes in the integration procedure approach zero.

The computational approach to hybrid systems we advocate in this paper (and in
[16],[14],[15]) assumes that for a given control function we attempt to follow
true systems trajectories as close as possible by using integration procedures
with a high order of convergence. It is in contrast to the approach which uses
approximations to the discontinuous right-hand side for the differential
equations or inclusions by a smooth right–hand side ([24],[22],[9]). The
smoothing approach demonstrates that it can be used to find approximate
solutions to optimal control problems with systems which exhibit sliding modes
however there is still a lack of evidence that it is capable of efficiently
finding optimal controls and state trajectories with high accuracy by choosing
proper values of smoothing parameters.

The organization of our paper is as follows. In the next section we describe
the systems with sliding motion. Then, we state continuous adjoint equations
for these systems. In Section 5 we formulate the discretized system equations
and the discrete adjoint equations and show the correspondence between the
discrete adjoint equations and discretization of the continuous adjoint
equations. On that basis, we show that the discrete adjoint trajectories
converge to their continuous counterparts for step sizes approaching zero. In
the first part of the paper we give an estimate for the order of convergence
of discrete adjoint equations assuming that a hybrid system is described by
ODEs. In the second part of the paper ([15]) we carry out analogical analysis
for the case when system equations correspond to index 2 DAEs. The second part
of the paper complements the result of the first part by paying attention to
the jumps of adjoint variables which occur when the hybrid system changes its
discrete state. We show that the jumps in the discrete adjoint equations
converge to jumps of their continuous counterparts provided that the
integration step sizes approach zero. In the second part of the paper we also
present the numerical results of applying our numerical procedure to several
optimal control problems.

##  2 Sliding mode

The paper concerns a numerical procedure for optimal control problems
described by hybrid systems with sliding modes. In this section we present the
concept of the sliding motion. Let us consider a hybrid system, which evolves
according to ODEs

(1) $x^{\prime}=f^{1}(x,u),$

if the condition

(2) $g(x)<0$

is satisfied. The system evolves according to ODEs

(3) $\displaystyle x^{\prime}=f^{2}(x,u)$

if the opposite condition

(4) $g(x)>0$

holds. Suppose that the system first evolves according to (1) and at the
moment $t_{t}$ reaches the switching surface

(5) $\Sigma=\\{x\in\mathbb{R}^{n}:g(x)=0\\}.$

If at $t_{t}$ the following conditions are satisfied

(6) $\displaystyle g_{x}(x(t_{t}))f^{1}(x(t_{t}),u(t_{t}))$ $\displaystyle>$
$\displaystyle 0$ (7) $\displaystyle g_{x}(x(t_{t}))f^{2}(x(t_{t}),u(t_{t}))$
$\displaystyle>$ $\displaystyle 0.$

the system starts to evolve according to ODEs (3). However if at $t_{t}$ the
condition

(8) $\displaystyle g_{x}(x(t_{t}))f^{2}(x(t_{t}),u(t_{t}))$ $\displaystyle<$
$\displaystyle 0$

is satisfied instead of (7), then both vector fields $f^{1}(x,u)$ and
$f^{2}(x,u)$ point towards the surface $\Sigma$ and we face the sliding motion
phenomenon ([3]). The system state remains on the surface $\Sigma$ and we say
that the system is in the sliding mode, which can be considered as the
separate discrete state of the system.

The system motion at the sliding mode can be handled by by applying the
concept of Filippov solutions ([3]). We say that a continuous state trajectory
is a Filippov solution of the considered hybrid system if

(9) $x^{\prime}=\left\\{\begin{array}[]{ll}f^{1}(x,u)&\text{if }g(x)<0\\\
f_{F}(x,u)&\text{if }g(x)=0\\\ f^{2}(x,u)&\text{if }g(x)>0\end{array}\right.$

where $f_{F}(x,u)$ is a convex combination of $f^{1}(x,u)$ and $f^{2}(x,u)$

(10) $f_{F}(x,u)=(1-\alpha)f^{1}(x,u)+\alpha f^{2}(x,u),\ \alpha\in[0,1].$

Since the coefficient $\alpha$ should guarantee that the system stays on the
surface $\Sigma$, the condition

(11) $g_{x}(x)f_{F}(x,u)=0.$

must be satisfied. Therefore, $\alpha$ is a function of $x$ and $u$ and can be
evaluated as

(12)
$\alpha(x,u)=\frac{g_{x}(x)f^{1}(x,u)}{g_{x}(x)\left(f^{1}(x,u)-f^{2}(x,u)\right)}.$

The hybrid system leaves the sliding mode when one of the vector fields
$f^{1}(x,u)$ and $f^{2}(x,u)$ points outside the surface $\Sigma$ at some time
$t_{t}$. This implies that a time $t_{t}$ we either have

(13) $\displaystyle\alpha(x(t_{t}),u(t_{t}))$ $\displaystyle=$ $\displaystyle
0$ (14) $\displaystyle g_{x}(x(t_{t}))f^{1}(x(t_{t}),u(t_{t}))$
$\displaystyle<$ $\displaystyle 0$ (15) $\displaystyle
g_{x}(x(t_{t}))f^{2}(x(t_{t}),u(t_{t}))$ $\displaystyle<$ $\displaystyle 0$

and the system evolves further according to (1) or

(16) $\displaystyle\alpha(x(t_{t}),u(t_{t}))-1$ $\displaystyle=$
$\displaystyle 0$ (17) $\displaystyle g_{x}(x(t_{t}))f^{1}(x(t_{t}),u(t_{t}))$
$\displaystyle>$ $\displaystyle 0$ (18) $\displaystyle
g_{x}(x(t_{t}))f^{2}(x(t_{t}),u(t_{t}))$ $\displaystyle>$ $\displaystyle 0$

and the system evolves further according to (3).

At the sliding mode the state trajectory should stay in a set $\Sigma$. During
numerical integration we will fulfill that if we integrate the differential-
algebraic equations (DAEs) (see [1] for the justification of the use of these
equations)

(19) $\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle
f_{F}(x,u)+g_{x}(x)^{T}z$ (20) $\displaystyle 0$ $\displaystyle=$
$\displaystyle g(x)$

instead of ODEs

(21) $x^{\prime}=f_{F}(x,u).$

One can show that under the condition $g_{x}(x)\neq 0$ at each $x\in\Sigma$
the equations (19)-(20) are index two DAEs ([14]).

##  3 Optimal control problem with sliding modes

The simple (due to the minimal number of discrete states) optimal control
problem with a hybrid system exhibiting a sliding mode—(P) can be stated as
follows:

(22) $\displaystyle\min_{u}\phi(x(t_{f}))$

subject to the constraints

$\displaystyle\begin{array}[]{ll}x^{\prime}=f_{1}(x,u)&\ \ \ {\rm if\ }\
g(x)<0\\\ x^{\prime}=f_{2}(x,u)&\ \ \ {\rm if\ }\ g(x)>0\\\
\left\\{\begin{array}[]{lll}x^{\prime}&=&f_{F}(x,u)+g_{x}(x)^{T}z\\\
0&=&g(x)\end{array}\right.&\ \ \ {\rm if\ }\ {\rm at\ a\ sliding\
mode}\end{array}$

and the terminal constraints

(27) $\displaystyle g^{1}_{i}(x(t_{f}))=0\ \forall i\in E$ (28) $\displaystyle
g^{2}_{j}(x(t_{f}))\leq 0\ \forall j\in I.$

Furthermore, due to the fact that $\alpha(x,u)$ depends on control functions
$u$ the admissible control is a function $u:[t_{0},t_{f}]\rightarrow{\cal U}$,
where

$\displaystyle{\displaystyle{\cal U}=\left\\{u\in{\cal L}^{2}_{m}[T]:\
u(t)=\sum_{j=1}^{N}{\bf 1}_{j}(t)u^{N}(j),\ u^{N}(j)\in
U\subset\mathbb{R}^{m}\right\\}}$

Here,

$\displaystyle{\displaystyle{\bf 1}_{j}(t)=\left\\{\begin{array}[]{ll}1,&{\rm
if}\ t\in[t_{j-1},t_{j}]\\\ 0,&{\rm if}\
t\not\in[t_{j-1},t_{j}]\end{array}\right.,}$

$t_{j}=jt_{f}/N$, $j\in\\{0,1,\ldots,N\\}$, and $U$ is a closed convex set.
The definition of the set of admissible control implies that we are dealing
with piecewise constant control functions.

Since our control functions are not continuous it is useful to introduce the
notation:

$\displaystyle{\displaystyle\lim_{t\rightarrow
t_{t},t<t_{t})}u(t)=u(t_{t}^{-}),\ \ \lim_{t\rightarrow
t_{t},t>t_{t})}u(t)=u(t_{t}^{+}),}$

which will be also applied to other functions, for example we will write
$x(t_{t}^{-})$, $\lambda(t_{t}^{+})$. The convention is also applied to
integrals, for example

$\displaystyle{\displaystyle\int_{t_{0}}^{t_{t}^{-}}p(x(t),u(t))dt=\lim_{t\rightarrow
t_{t},t<t_{t}}\int_{t_{0}}^{t}p(x(\tau),u(\tau))d\tau.}$

We can show, under some standard assumptions stated in [14], that for given
$x(t_{0})$ state trajectory $x$ is uniquely defined by $u$ (so we can write
$x^{u}$ instead of $u$) thus the considered control problem can be expressed
as an optimization problem over the set of control functions with the aid of
the functions $\bar{F}_{0}:{\cal U}\rightarrow{\cal R}$,
$\bar{g}^{1}_{i}:{\cal U}\rightarrow{\cal R}$ for $i\in E$,
$\bar{g}^{2}_{j}:{\cal U}\rightarrow{\cal R}$ for $j\in I$:

$\displaystyle\bar{F}_{0}(u)$ $\displaystyle=$
$\displaystyle\phi(x^{u}(t_{f}))$ $\displaystyle\bar{g}^{1}_{i}(u)$
$\displaystyle=$ $\displaystyle g^{1}_{i}(x^{u}(t_{f}))\ \forall i\in E$
$\displaystyle\bar{g}^{2}_{j}(u)$ $\displaystyle=$ $\displaystyle
g^{2}_{j}(x^{u}(t_{f}))\ \forall j\in I.$

Eventually, the optimal control problem we want to solve can be stated as

(30) $\displaystyle{\displaystyle\min_{u\in{\cal U}}\bar{F}_{0}(u)}$

subject to

(31) $\displaystyle\bar{g}^{1}_{i}(u)$ $\displaystyle=$ $\displaystyle 0\
\forall i\in E$ (32) $\displaystyle\bar{g}^{2}_{j}(u)$ $\displaystyle\leq$
$\displaystyle 0\ \forall j\in I.$

##  4 Calculating the optimal control

In [14] we propose the method for solving the problem ${\bf(P)}$. This is an
exact penalty algorithm discussed in [12] which was used to solve several
optimal control problems including those described by DAEs ([13],[17]). Since
the optimization method we advocate for solving for solving the problem
${\bf(P)}$ is based on an exact penalty function, instead of solving the
problem ${\bf(P)}$ we solve the problem ${\bf P_{c}}$

(33) $\displaystyle{\displaystyle\min_{u\in{\cal U}}\bar{F}_{c}(u)}$

in which the exact penalty function $\\{\bar{F}_{c}(u)$ is defined as follows

(34)
$\displaystyle{\displaystyle\bar{F}_{c}(u)=\bar{F}(u)+c\max\left[0,\max_{i\in
E}\left|\bar{g}^{1}_{i}(u)\right|,\max_{j\in I}\bar{g}^{2}_{j}(u)\right]}$

Let $H\in{\cal R}^{mN\times mN}$ be a symmetric matrix such that

(35) $\displaystyle{\displaystyle\nu_{1}\|d\|^{2}\leq
d^{T}Hd\leq\nu_{2}\|d\|^{2}\ \forall d\in{\cal R}^{mN}}$

for some positive $\nu_{1}$ and $\nu_{2}$.

For fixed $c$, $H$ satisfying (35) and $u$ the direction finding subproblem,
${\bf P_{c}(u)}$, for the problem ${\bf(P_{c})}$ is:

$\displaystyle{\displaystyle\min_{d\in{\cal U}-u,\beta\in{\cal
R}}\left[\left\langle\nabla\bar{F}_{0}(u),d\right\rangle+c\beta+1/2d^{T}Hd\right]}$

subject to

$\displaystyle\left|\bar{g}^{1}_{i}(u)+\left\langle\nabla\bar{g}^{1}_{i}(u),d\right\rangle\right|$
$\displaystyle\leq$ $\displaystyle\beta\ \ \forall i\in E$
$\displaystyle\bar{g}^{2}_{j}(u)+\left\langle\nabla\bar{g}^{2}_{j}(u),d\right\rangle$
$\displaystyle\leq$ $\displaystyle\beta\ \ \forall j\in I.$

The subproblem can be reformulated as an optimization problem over the space
${\cal R}^{mN}$ with the objective function which is strictly convex.

We define descent function $\sigma_{c}^{H}(u)$ and penalty test function
$t_{c}(u)$ (see [12] for details), to be used to test optimality of a control
$u$ and to adjust $c$, respectively, as

$\displaystyle{\displaystyle\sigma_{c}^{H}(u)=\left\langle\nabla\bar{F}_{0}(u),\bar{d}\right\rangle+c\left[\bar{\beta}-M(u)\right]}$

and

(36) $\displaystyle{\displaystyle t_{c}(u)=\sigma_{c}^{H}(u)+M(u)/c}$

for given $c>0$, $H$ and $u\in{\cal U}$. Here,

$\displaystyle{\displaystyle M(u)=\max\left[0,\max_{i\in
E}\left|\bar{g}^{1}_{i}(u)\right|,\max_{j\in I}\bar{g}^{2}_{j}(u)\right],}$

A general algorithm is as follows.
Algorithm Fix parameters: $\varepsilon>0,\ \gamma,\ \eta\in(0,1)$, $c^{0}>0$,
$\kappa>1$, matrix $H$.

  1. 1.

Choose the initial control $u_{0}\in{\cal U}$. Set $k=0$, $c_{-1}=c^{0}$.

  2. 2.

Let $c_{k}$ be the smallest number chosen from $\\{c_{k-1},\kappa
c_{k-1},\kappa^{2}c_{k-1},\ldots\\}$ such that the solution
$(d_{k},\beta_{k})$ to the direction finding subproblem ${\bf
P_{c_{k}}(u_{k})}$ satisfies

(37) $\displaystyle{\displaystyle t_{c_{k}}(u_{k})\leq 0.}$

If $\sigma_{c_{k}}^{H}(u_{k})=0$ then STOP.

  3. 3.

Let $\alpha_{k}$ be the largest number chosen from the set
$\\{1,\eta,\eta^{2},\ldots,\\}$ such that

$\displaystyle{\displaystyle u_{k+1}=u_{k}+\alpha_{k}d_{k}}$

satisfies the relation

$\displaystyle\bar{F}_{c_{k}}(u_{k+1})-\bar{F}_{c_{k}}(u_{k})$
$\displaystyle\leq$ $\displaystyle\gamma\alpha_{k}\sigma_{c_{k}}^{H}(u_{k}).$

Increase k by one. Go to Step 2.

We can show that the descent function $\sigma_{c_{k}}^{H}(u_{k})$ is
nonpositive valued at each iteration ([12]).

The main computational effort is associated with the evaluation of function
values $\bar{F}_{0}$, $\bar{g}_{i}^{1}$, $i\in E$, $\bar{g}_{j}^{2}$, $j\in I$
and their gradients $\nabla\bar{F}_{0}$, $\nabla\bar{g}_{i}^{1}$, $i\in E$,
$\nabla\bar{g}_{j}^{2}$, $j\in I$. The number of decision variables is equal
to $mN$ and for the fine discretization of control variables it can be counted
in hundreds. Therefore, it can be more efficient to use adjoint equations
instead of relying on sensitivity equations for evaluating gradients (see
Chapter 6 in [12]).

The cost of evaluating functions values can be significant since it requires
the integration of ODEs and DAEs if sliding modes occur. Higher index DAEs
usually are integrated by implicit methods (BDF, or Runge–Kutta—[1], [5]).
Since we approximate controls by piecewise constant functions we can expect
that multistep integration procedure would have to be often restarted which
would deteriorate efficiency of numerical integration. Runge–Kutta methods are
one–step methods so they are more appropriate in optimal control applications
([23],[12],[4],[11]). However, we are constrained by implicit methods and that
means that the cost of numerical integration is significant due to the
necessity of solving nonlinear algebraic equations (and possibly due to
evaluating equations functions Jacobians and their factorizations) at each
integration step.

Furthermore, adjoint equations to higher index DAEs are also higher index
DAEs, so the mentioned above efficiency issues apply also to evaluating
adjoint equations. Since each function gradient requires solution of unique
adjoint equations (due to different terminal conditions) at each optimization
step we have to solve $|E|+|I|+1$ adjoint equations.

Therefore, it is tempting to use as much as possible information acquired
during system equations integration in numerical evaluation of adjoint
equations. It can be achieved, as we will show in the next section, if adjoint
equations are built for discrete time system equations resulting from the
numerical integration. However, if we do that we have to show whether these
discrete time adjoint equations will give the solution to their continuous
counterparts provided that integration steps approach zero. This is the main
subject of the next section.

The ’reference’ continuous time adjoint equations are established in our
accompanying paper [14]. Let us recall only one case, when before the
switching time $t_{t}$ the system evolves according to ODEs

(38) $x^{\prime}=f^{1}(x,u),$

at $t_{t}$ system enters the sliding mode and then the state trajectory is the
solution to DAEs

(39) $\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle
f_{F}(x,u)+g_{x}^{T}(x)z$ (40) $\displaystyle 0$ $\displaystyle=$
$\displaystyle g(x)$

up to an ending time $t_{f}$. At a transition time $t_{t}$ the continuous
state trajectory meets the switching surface such that the following condition
holds

(41) $g(x(t_{t}))=0.$

The adjoint equations for the system (39)–(40) are as follows ([14])

$\displaystyle(\lambda_{f}^{T})^{\prime}(t)$ $\displaystyle=$
$\displaystyle-\lambda_{f}^{T}(t)(f_{F})_{x}(x(t),u(t))$
$\displaystyle-\lambda_{f}^{T}(t)(g_{x}^{T}(x(t))z(t))_{x}+\lambda_{g}^{T}(t)g_{x}(x(t))$
(43) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\lambda_{f}^{T}(t)g_{x}^{T}(x(t)),\ t\in[t_{t}^{+},t_{f}].$

and the consistent terminal values of the adjoint variables
$\lambda_{f}(t_{f})$ and $\lambda_{g}(t_{f})$ can be found by solving the
following set of equations with respect to the variables
$\lambda_{f}(t_{f}),\lambda_{g}(t_{f}),\nu_{1}$

(44) $\displaystyle\phi_{x}^{T}(x(t_{f}))+\lambda_{f}(t_{f})$ $\displaystyle=$
$\displaystyle\nu_{1}g_{x}^{T}(x(t_{f}))$ (45) $\displaystyle 0$
$\displaystyle=$ $\displaystyle g_{x}(x(t_{f}))\lambda_{f}(t_{f})$
$\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(g_{x}(x(t_{f}))\right)^{\prime}\lambda_{f}(t_{f})-$
$\displaystyle
g_{x}(x(t_{f}))\left(f_{F}\right)_{x}^{T}(x(t_{f}),u(t_{f}))\lambda_{f}(t_{f})-$
$\displaystyle
g_{x}(x(t_{f}))\left(g_{x}^{T}(x(t_{f}))z(t_{f})\right)_{x}^{T}\lambda_{f}(t_{f})+$
$\displaystyle g_{x}(x(t_{f}))g_{x}^{T}(x(t_{f}))\lambda_{g}(t_{f})$

where $\nu_{1}$ is some real number.

At the transition time $t_{t}$ the adjoint variable $\lambda_{f}$ undergoes a
jump. To calculate the value of $\lambda_{f}(t_{t}^{-})$ the following system
of equations have to be solved for the variables $\lambda_{f}(t_{t}^{-}),\
\pi$

(47) $\displaystyle\lambda_{f}(t_{t}^{-})$ $\displaystyle=$
$\displaystyle\lambda_{f}(t_{t}^{+})-\pi g_{x}^{T}(x(t_{t}))$
$\displaystyle\lambda_{f}^{T}(t_{t}^{-})f_{1}(x(t_{t}^{-}),u(t_{t}^{-}))$
$\displaystyle=$
$\displaystyle\lambda_{f}^{T}(t_{t}^{+})f_{F}(x(t_{t}^{+}),u(t_{t}^{+}))+$
$\displaystyle\lambda_{f}^{T}(t_{t}^{+})g_{x}^{T}(x(t_{t}^{+}))z(t_{t}^{+})-\lambda_{g}^{T}(t_{t}^{+})g(x(t_{t}^{+})).$

Eventually we solve adjoint ODEs

(49)
$(\lambda_{f}^{T})^{\prime}(t)=-\lambda_{f}^{T}(t)(f_{1})_{x}(x(t),u(t)),\
t\in[t_{0},t_{t}^{-}).$

Now the directional derivative $\langle\nabla\bar{F}_{0}(u),d\rangle$ of the
cost functional $\bar{F}_{0}$ takes the form

$\displaystyle\langle\nabla\bar{F}_{0}(u),d\rangle$ $\displaystyle=$
$\displaystyle-\int_{t_{0}}^{t_{t}^{-}}\lambda_{f}^{T}(t)(f_{1})_{u}(x(t),u(t))d(t)dt$
$\displaystyle-\int_{t_{t}^{+}}^{t_{f}}\lambda_{f}^{T}(t)(f_{F})_{u}(x(t),u(t))d(t)dt.$

In the next section we will examine the discrete time versions of equations
(49) and the approximation to the first term in (4). The approximation to the
second term and the analysis of discrete time versions of equations (4)–(43)
is provided in the second part of the paper [15].

##  5 Numerical calculation of reduced gradients

In this chapter we present our approach to the numerical calculation of the
reduced gradients of some endpoint state functional against the controls and
discuss the convergence of the proposed numerical scheme. We focus on the
numerical integration of continuous equations, as the numerical calculations
at transition and end points reduces to solving the system of algebraic
equations. The calculation of reduced gradients requires numerical integration
of system equations, adjoint equations and reduced gradient formula. We
discuss separately the ODEs case and DAEs case (which is presented in [15]).

We start our considerations from a control system described by ODEs

(51) $x^{\prime}(t)=f(x(t),u(t)),\ x(t_{0})=x_{0}.$

The system equations (51) are numerically integrated over the time interval
$[t_{0},t_{f}]$ using the Runge-Kutta scheme ([6])

(52) $x_{i}(k+1)=x(k)+h(k)\sum_{j=1}^{s}a_{ij}f\left(x_{j}(k+1),\
u(t(k)+c_{i}h(k))\right)$

for $i=1,\dots,s$ and

(53) $x(k+1)=x(k)+h(k)\sum_{i=1}^{s}b_{i}f\left(x_{i}(k+1),\
u(t(k)+c_{i}h(k))\right).$

where $t(k)$ is time at the actual step, $h(k)$ is the actual step size and
$x(k)\simeq x(t(k))$ is the numerical approximation of the state at $t(k)$.
The constant coefficients $a_{ij}$, $b_{i}$, $c_{i}$ for $i,j=1,\dots,s$
define the Runge-Kutta scheme. At each step of the Runge-Kutta scheme the
nonlinear system (52) is first solved for variables $x_{i}(k+1),\ i=1,\dots,s$
and then (53) is used to calculate $x(k+1)$. The scheme is repeated for steps
$k=0,\dots,K-1$, where $x(0)=x_{0}$, $t(0)=t_{0}$ and $t(K)=t_{f}$. In our
code we assume the control variable is constant along the single integration
step, so we have

(54) $u(t(k)+c_{i}h(k))=u(k),\ i=1,\ldots,s.$

In the next step we must calculate numerically the adjoint equations

(55) $\lambda^{\prime}(t)=-f_{x}^{T}(x(t),u(t))\lambda(t),$

The adjoint equations are solved backwards in time with an endpoint condition
given at time $t_{f}$. The adjoint equations are ODEs, so the Runge-Kutta
scheme can be used to numerically integrate them. The Runge-Kutta scheme for
adjoint equations is

$\displaystyle\lambda_{i}(k)$ $\displaystyle=$
$\displaystyle\lambda(k+1)-\bar{h}(k+1)\sum_{j=1}^{\bar{s}}\bar{a}_{ij}\left(-f_{x}^{T}\left(x\left(\bar{t}(k+1)-\bar{c}_{j}\bar{h}(k+1)\right),\right.\right.$
$\displaystyle\left.\left.\bar{u}(k+1)\right)\lambda_{j}(k)\right)$
$\displaystyle=$
$\displaystyle\lambda(k+1)+\bar{h}(k+1)\sum_{j=1}^{\bar{s}}\bar{a}_{ij}f_{x}^{T}\left(x\left(\bar{t}(k+1)-\bar{c}_{j}\bar{h}(k+1)\right),\right.$
$\displaystyle\left.\bar{u}(k+1)\right)\lambda_{j}(k)$

for $i=1,\dots,\bar{s}$ and

$\displaystyle\lambda(k)$ $\displaystyle=$
$\displaystyle\lambda(k+1)-\bar{h}(k+1)\sum_{i=1}^{\bar{s}}\bar{b}_{i}\left(-f_{x}^{T}\left(x\left(\bar{t}(k+1)-\bar{c}_{i}\bar{h}(k+1)\right),\right.\right.$
$\displaystyle\left.\left.\bar{u}(k+1)\right)\lambda_{i}(k)\right)$
$\displaystyle=$
$\displaystyle\lambda(k+1)+\bar{h}(k+1)\sum_{i=1}^{\bar{s}}\bar{b}_{i}f_{x}^{T}\left(x\left(\bar{t}(k+1)-\bar{c}_{i}\bar{h}(k+1)\right),\right.$
$\displaystyle\left.\bar{u}(k+1)\right)\lambda_{i}(k),$

where $\bar{t}(k)$ is time at the step $k$, $\bar{h}(k)>0$ is the step size at
step $k$ and $\lambda(k)\simeq\lambda(t(k))$ is the approximation of the
adjoint state at $t(k)$. At each step of the Runge-Kutta scheme, the system of
equations (LABEL:discrEq:rkSchemeIntAdjODE) is first solved for variables
$\lambda_{i}(k),\ i=1,\dots,\bar{s}$ and then
(LABEL:discrEq:rkSchemeStepAdjODE) is used to calculate $\lambda(k)$. The
scheme is repeated for steps $k=\bar{K}-1,\dots,0$, where $t(0)=t_{0}$ and
$t(\bar{K})=t_{f}$. In general the number of steps, time moments at particular
steps, step sizes and control values may be different for numerical
integration of system and adjoint equations. It has been denoted by adding the
bar symbol over the appropriate variables ($\bar{K}$, $\bar{t}(k)$,
$\bar{h}(k)$, $\bar{u}(k)$) in
(LABEL:discrEq:rkSchemeIntAdjODE)-(LABEL:discrEq:rkSchemeStepAdjODE). What is
more, the Runge-Kutta schemes used for numerical integration of system and
adjoint equations can also be different, what has been emphasized by
introducing the coefficients $\bar{a}_{ij}$, $\bar{b}_{i}$, $\bar{c}_{i}$ for
$i,j=1,\dots,\bar{s}$ in the Runge-Kutta scheme for adjoint equations. The
additional difficulty in implementation of the scheme
(LABEL:discrEq:rkSchemeIntAdjODE)-(LABEL:discrEq:rkSchemeStepAdjODE), is that
a subprocedure for the reliable approximation of the state $x(t)$ for
$t\in\left[\bar{t}(k),\bar{t}(k+1)\right]$ is required, such that terms
$x\left(\bar{t}(k+1)-\bar{c}_{i}\bar{h}(k+1)\right)$ are properly
approximated.

To derive the numerical scheme for calculation of reduced gradients, let us
consider the differential of the endpoint state functional (corresponding to a
perturbation $d$) of the form

(58)
$d\phi(x(t_{f}))=\int_{t_{0}}^{t_{f}}d(t)^{T}f_{u}^{T}(x(t),u(t))\lambda(t)dt.$

The control signal is approximated as a piecewise constant function

(59) $u(t)=u_{n}\ for\ t\in\left(t_{n-1},t_{n}\right],\ n=1,\dots,N,$

where

(60) $t_{n}=t_{0}+n\frac{t_{f}-t_{0}}{N}.$

Now (58) can be written as

(61)
$d\phi(x(t_{f}))=\sum_{l=1}^{N}d^{T}_{n}\int_{t_{n-1}}^{t_{n}}f_{u}^{T}(x(t),u_{n})\lambda(t)dt,$

where $d_{n}$ is a differential of the control signal $u_{n}$. The derivative
of $\phi(x(t_{f}))$ with respect to $u_{n}$ is therefore

(62)
$\frac{d\phi(x(t_{f}))}{du_{n}}=\int_{t_{n-1}}^{t_{n}}f_{u}^{T}(x(t),u_{n})\lambda(t)dt.$

To approximate numerically the above integral, we choose times

(63) $\tilde{t}(0)<\dots<\tilde{t}(k)<\dots<\tilde{t}(\tilde{K}),$

such that $\tilde{t}(0)=t_{0}$, $\tilde{t}(\tilde{K})=t_{f}$ and for each $n$
there is an index $k_{n}$, such that

(64) $\tilde{t}(k_{n})=t_{n}.$

The values $\tilde{t}(0),\dots,\tilde{t}(k),\dots,\tilde{t}(\tilde{K})$
defines the division of the interval $[t_{0},t_{f}]$ into the subintervals.
Now (62) can be written as

(65)
$\frac{d\phi(x(t_{f}))}{du_{n}}=\sum_{k=k_{n-1}}^{k_{n}-1}\int_{\tilde{t}(k)}^{\tilde{t}(k+1)}f_{u}^{T}(x(t),u_{n})\lambda(t)dt.$

The integrals in (65) can be approximated by the quadrature formula

$\displaystyle{\displaystyle\int_{\tilde{t}(k)}^{\tilde{t}(k+1)}f_{u}^{T}(x(t),u_{n})\lambda(t)dt\simeq\tilde{h}(k)\sum_{i=1}^{\tilde{s}}\tilde{b}_{i}f_{u}^{T}\left(x\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right),u_{n}\right)\times}$
(66)
$\displaystyle{\displaystyle\lambda\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right),}$

where $\tilde{h}(k)=\tilde{t}(k+1)-\tilde{t}(k)$ is the step size at step $k$
and $\tilde{b}_{i},\tilde{c}_{i},\ i=1,\cdots,\tilde{s}$ are the quadrature
coefficients. Again, the number of steps, time moments at particular steps and
step sizes used for numerical integration of reduced gradients may be
different from the respective values used for numerical integration of system
or adjoint equations. It has been denoted by adding the tilde symbol over the
appropriate variables ($\tilde{K}$, $\tilde{t}(k)$, $\tilde{h}(k)$) in the
above formulas. The additional difficulty in implementation of the quadrature
formula (66), is that a subprocedure for the reliable approximation of the
state $x(t)$ and adjoint state $\lambda(t)$ for
$t\in\left[\tilde{t}(k),\tilde{t}(k+1)\right]$ is required, such that terms
$x\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right)$ and
$\lambda\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right)$ are properly
approximated.

In our code we utilize different approach to numerical calculation of adjoint
equations and reduced gradients. The proposed approach is based on the concept
of discrete state equations, discrete adjoint equations and discrete reduced
gradients. Using that approach, the subprocedure for the approximation of the
state $x(t)$ and adjoint state $\lambda(t)$ for
$t\in\left[\tilde{t}(k),\tilde{t}(k+1)\right]$ is not required.

Let us start from formulation of the discrete equations. Let us rewrite the
numerical integration scheme (52)-(53) in a vector form

(67) $\left(\begin{array}[]{c}x_{1}-x-h\sum_{j=1}^{s}a_{1j}f(x_{j},u)\\\
\vdots\\\ x_{s}-x-h\sum_{j=1}^{s}a_{sj}f(x_{j},u)\\\
x^{+}-x-h\sum_{i=1}^{s}b_{i}f(x_{i},u)\\\
\end{array}\right)=\left(\begin{array}[]{c}0\\\ \vdots\\\ 0\\\
0\end{array}\right).$

For the sake of a more compact notation we omitted the discrete step argument
and introduced the symbols $x_{i}=x_{i}(k+1)$, $x=x(k)$, $x^{+}=x(k+1)$,
$h=h(k),u=u(k)$. If we now define the augmented state vector $X(k)$ as

(68) $X(k)=\left(x_{1}(k)^{T},\ldots,x_{s}(k)^{T},x(k)^{T}\right)^{T},$

then (67) can be presented in a form of the implicit discrete time state
equation

(69) $F\left(X(k+1),X(k),u(k)\right)=0.$

Having the discrete state equation (69) defined we may consider the discrete
time adjoint equations ([13])

(70) $\Lambda(k)=-F_{X}^{T}(k)\left(F_{X^{+}}^{T}(k)\right)^{-1}\Lambda(k+1),$

where $\Lambda(k)$ is the discrete adjoint variable at a discrete step $k$ and

(71) $\displaystyle F_{X^{+}}(k)$ $\displaystyle=$
$\displaystyle\frac{\partial F(X(k+1),X(k),u(k))}{\partial X(k+1)},$ (72)
$\displaystyle F_{X}(k)$ $\displaystyle=$ $\displaystyle\frac{\partial
F(X(k+1),X(k),u(k))}{\partial X(k)}.$

The equation (70) is equivalent to a system of equations

(73) $\displaystyle F_{X^{+}}^{T}(k)R(k)$ $\displaystyle=$
$\displaystyle\Lambda(k+1),$ (74) $\displaystyle\Lambda(k)$ $\displaystyle=$
$\displaystyle-F_{X}^{T}(k)R(k),$

where $R(k+1)$ is an auxiliary variable introduced to avoid calculating the
inverse of the matrix $F_{X^{+}}^{T}(k)$.

Now we are going to rewrite the discrete adjoint equations in an element-wise
manner. The partial derivatives matrices are

(75)
$F_{X^{+}}(k)=\left(\begin{array}[]{cccc}I-ha_{11}f_{x1}&\ldots&-ha_{1s}f_{xs}&0\\\
\vdots&&\vdots&\vdots\\\ -ha_{s1}f_{x1}&\ldots&I-ha_{ss}f_{xs}&0\\\
-hb_{1}f_{x1}&\ldots&-hb_{s}f_{xs}&I\end{array}\right),$ (76)
$F_{X}(k)=\left(\begin{array}[]{cccc}0&\ldots&0&-I\\\ \vdots&&\vdots&\vdots\\\
0&\ldots&0&-I\\\ 0&\ldots&0&-I\end{array}\right),$

where for the sake of a shorter notation we denote

(77) $\displaystyle f_{xi}$ $\displaystyle=$ $\displaystyle f_{x}(x_{i},u).$

Let us denote

$\displaystyle\Lambda(k)$ $\displaystyle=$
$\displaystyle\left(l_{1}(k)^{T},\ldots,l_{s}(k)^{T},\lambda(k)^{T}\right)^{T},$
$\displaystyle R(k)$ $\displaystyle=$
$\displaystyle\left(r_{1}(k)^{T},\ldots,r_{s}(k)^{T},r(k)^{T}\right)^{T}.$

Now (73) at a discrete time step $k$ takes the form

(78)
$\left(\begin{array}[]{cccc}I-ha_{11}f^{T}_{x1}&\ldots&-ha_{s1}f^{T}_{x1}&-hb_{1}f^{T}_{x1}\\\
\vdots&&\vdots&\vdots\\\
-ha_{1s}f^{T}_{xs}&\ldots&I-ha_{ss}f^{T}_{xs}&-hb_{s}f^{T}_{xs}\\\
0&\ldots&0&I\\\ \end{array}\right)\left(\begin{array}[]{c}r_{1}\\\ \vdots\\\
r_{s}\\\ r\end{array}\right)=\left(\begin{array}[]{c}l^{+}_{1}\\\ \vdots\\\
l^{+}_{s}\\\ \lambda^{+}\end{array}\right)$

and (74) takes the form

(79) $\left(\begin{array}[]{c}l_{1}\\\ \vdots\\\ l_{s}\\\
\lambda\end{array}\right)=-\left(\begin{array}[]{cccc}0&\ldots&0&0\\\
\vdots&&\vdots&\vdots\\\ 0&\ldots&0&0\\\ -I&\ldots&-I&-I\\\
\end{array}\right)\left(\begin{array}[]{c}r_{1}\\\ \vdots\\\ r_{s}\\\
r\end{array}\right).$

For the sake of a more compact notation we omitted the discrete step argument
and introduced the symbols $r_{i}=r_{i}(k)$, $r=r(k)$, $l^{+}_{i}=l_{i}(k+1)$,
$\lambda^{+}=\lambda(k+1)$, $l_{i}=l_{i}(k)$, $\lambda=\lambda(k)$. If we
rewrite (78) as a system of equations we obtain

(80) $r_{i}=\sum_{j=1}^{s}ha_{ji}f_{xi}^{T}r_{j}+hb_{i}f_{xi}^{T}r+l^{+}_{i}$

for $i=1,\ldots,s$ and

(81) $r=\lambda^{+}.\\\ $

Using (81), (80) can be written as

(82)
$r_{i}=\sum_{j=1}^{s}ha_{ji}f_{xi}^{T}r_{j}+hb_{i}f_{xi}^{T}\lambda^{+}+l^{+}_{i}.$

If we rewrite (79) as a system of equations we get

(83) $l_{i}=0$

for $i=1,\ldots,s$ and

(84) $\lambda=\sum_{i=1}^{s}r_{i}+r.$

Using (81), (84) becomes

(85) $\lambda=\sum_{i=1}^{s}r_{i}+\lambda^{+}.$

From (83) we have that $l_{1}(k)=0,\ldots,l_{s}(k)=0$ for steps $k=0,...,K-1$.
During the analysis we also assume that $l_{1}(K)=0,\ldots,l_{s}(K)=0$ . Under
that assumption (82) is equivalent to

(86) $r_{i}=\sum_{j=1}^{s}ha_{ji}f_{xi}^{T}r_{j}+hb_{i}f_{xi}^{T}\lambda^{+}$

for $i=1,\dots,s$.

To show the resemblance of discrete adjoint equations and the Runge-Kutta
scheme for continuous adjoint equations, we introduce the auxiliary variables
in a manner proposed in ([4])

(87) $\lambda_{i}=\sum_{j=1}^{s}\frac{a_{ji}}{b_{i}}r_{j}+\lambda^{+}$

for $i=1,\ldots,s$. We have then

(88) $r_{i}=hb_{i}f_{xi}^{T}\lambda_{i}$

and (87) can be transformed to the following form

(89)
$\lambda_{i}=h\sum_{j=1}^{s}\frac{a_{ji}b_{j}}{b_{i}}f_{xj}^{T}\lambda_{j}+\lambda^{+}$

and (85) can be transformed to

(90) $\lambda=h\sum_{i=1}^{s}b_{i}f_{xi}^{T}\lambda_{i}+\lambda^{+}.$

Let us rewrite (89) and (90) by inserting the discrete step argument

(91)
$\lambda_{i}(k)=\lambda(k+1)+h(k)\sum_{j=1}^{s}\frac{a_{ji}b_{j}}{b_{i}}f_{x}^{T}(x_{j}(k+1),u(k))\lambda_{j}(k)$

for $i=1,\dots,s$ and

(92)
$\lambda(k)=\lambda(k+1)+h(k)\sum_{i=1}^{s}b_{i}f_{x}^{T}(x_{i}(k+1),u(k))\lambda_{i}(k).$

For the convenience of reading we recall here the Runge-Kutta scheme for
adjoint equations
(LABEL:discrEq:rkSchemeIntAdjODE)-(LABEL:discrEq:rkSchemeStepAdjODE)

$\lambda_{i}(k)=\lambda(k+1)+\bar{h}(k+1)\sum_{j=1}^{\bar{s}}\bar{a}_{ij}f_{x}^{T}\left(x\left(\bar{t}(k+1)-\bar{c}_{j}\bar{h}(k+1)\right),\bar{u}(k+1)\right)\lambda_{j}(k)$

for $i=1,\dots,s$ and

$\lambda(k)=\lambda(k+1)+\bar{h}(k+1)\sum_{i=1}^{\bar{s}}\bar{b}_{i}f_{x}^{T}\left(x\left(\bar{t}(k+1)-\bar{c}_{i}\bar{h}(k+1)\right),\bar{u}(k+1)\right)\lambda_{i}(k).$

Let us suppose, that the adjoint equations are integrated such that discrete
times are the same as during the numerical integration of system equations so

(93) $\bar{t}(k)=t(k).$

In such case we obtain

(94) $\displaystyle\bar{h}(k+1)$ $\displaystyle=$ $\displaystyle h(k),$ (95)
$\displaystyle\bar{u}(k+1)$ $\displaystyle=$ $\displaystyle u(k)$

and

(96)
$\bar{t}(k+1)-\bar{c}_{i}\bar{h}(k+1)=t(k+1)-\bar{c}_{i}h(k)=t(k)+(1-\bar{c}_{i})h(k).$

By introducing the above relations into
(LABEL:discrEq:rkSchemeIntAdjODE)-(LABEL:discrEq:rkSchemeStepAdjODE) we obtain

(97)
$\lambda_{i}(k)=\lambda(k+1)+h(k)\sum_{j=1}^{\bar{s}}\bar{a}_{ij}f_{x}^{T}\left(x\left(t(k)+(1-\bar{c}_{j})h(k)\right),u(k)\right)\lambda_{j}(k)$

for $i=1,\dots,s$ and

(98)
$\lambda(k)=\lambda(k+1)+h(k)\sum_{i=1}^{\bar{s}}\bar{b}_{i}f_{x}^{T}\left(x\left(t(k)+(1-\bar{c}_{i})h(k)\right),u(k)\right)\lambda_{i}(k).$

If we now assume that

(99) $\bar{s}=s,\ \bar{a}_{ij}=\frac{a_{ji}b_{j}}{b_{i}},\ \bar{b}_{i}=b_{i},\
\bar{c}_{i}=1-c_{i},\ i,j=1,\dots,s$

(97)-(98) can be rewritten as

(100)
$\lambda_{i}(k)=\lambda(k+1)+h(k)\sum_{j=1}^{s}\frac{a_{ji}b_{j}}{b_{i}}f_{x}^{T}\left(x\left(t(k)+c_{j}h(k)\right),u(k)\right)\lambda_{j}(k)$

for $i=1,\dots,s$ and

(101)
$\lambda(k)=\lambda(k+1)+h(k)\sum_{i=1}^{s}b_{i}f_{x}^{T}\left(x\left(t(k)+c_{i}h(k)\right),u(k)\right)\lambda_{i}(k).$

Now the Runge-Kutta scheme (100)-(101) for adjoint equations is almost
identical to discrete adjoint equations (91)-(92). The only difference is that
in discrete equations $x_{i}(k+1)$ is used instead of $x(t(k)+c_{i}h(k))$. The
convergence of (100)-(101) is guaranteed if only the Runge-Kutta scheme
$\bar{a}_{ij},\ \bar{b}_{i},\ \bar{c}_{i},\ i,j=1,\dots,s$ satisfies
appropriate conditions. We will also justify that the usage of $x_{i}(k+1)$
instead of $x(t(k)+c_{i}h(k))$ does not break the convergence of (91)-(92) to
continuous adjoint trajectory.

The discrete reduced gradients are defined as ([13])

(102)
$\frac{d\phi(X(K))}{du_{n}}=\sum_{k=k_{n-1}}^{k_{n}-1}-F_{u}^{T}(k)\left(F_{X^{+}}^{T}(k)\right)^{-1}\Lambda(k+1)=\sum_{k=k_{n-1}}^{k_{n}-1}-F_{u}^{T}(k)R(k+1),$

where

(103) $F_{u}(k)=\frac{\partial F(X(k+1),X(k),u(k))}{\partial u(k)}.$

The above partial derivative rewritten in the element-wise manner is

(104) $F_{u}(k)=\left(\begin{array}[]{c}-h\sum_{j=1}^{s}a_{1j}f_{uj}\\\
\vdots\\\ -h\sum_{j=1}^{s}a_{sj}f_{uj}\\\
-h\sum_{i=1}^{s}b_{i}f_{ui}\end{array}\right),$

where for the sake of the shorter notation we omitted the discrete step
argument and introduced

(105) $f_{ui}=f_{u}(x_{i},u).$

Let us now inspect $-F_{u}^{T}(k)R(k+1)$

$\displaystyle-F_{u}^{T}(k)R(k+1)$ $\displaystyle=$
$\displaystyle-\left(\begin{array}[]{c}-h\sum_{j=1}^{s}a_{1j}f_{uj}\\\
\vdots\\\ -h\sum_{j=1}^{s}a_{sj}f_{uj}\\\
-h\sum_{i=1}^{s}b_{i}f_{ui}\end{array}\right)^{T}\left(\begin{array}[]{c}r_{1}\\\
\vdots\\\ r_{s}\\\ r\end{array}\right)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{s}h\left(\sum_{j=1}^{s}a_{ij}f^{T}_{uj}\right)r_{i}+h\sum_{i=1}^{s}b_{i}f^{T}_{ui}r.$

From (81) and (88) we get

$\displaystyle-F_{u}^{T}(k)R(k+1)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{s}h\left(\sum_{j=1}^{s}a_{ij}f^{T}_{uj}\right)hb_{i}f_{xi}^{T}\lambda_{i}+h\sum_{i=1}^{s}b_{i}f^{T}_{ui}\lambda^{+}$
$\displaystyle=$ $\displaystyle
h\sum_{j=1}^{s}f^{T}_{uj}\left(h\sum_{i=1}^{s}a_{ij}b_{i}f_{xi}^{T}\lambda_{i}\right)+h\sum_{i=1}^{s}b_{i}f^{T}_{ui}\lambda^{+}$
$\displaystyle=$ $\displaystyle
h\sum_{j=1}^{s}f^{T}_{uj}b_{j}\left(h\sum_{i=1}^{s}\frac{a_{ij}b_{i}}{b_{j}}f_{xi}^{T}\lambda_{i}\right)+h\sum_{i=1}^{s}b_{i}f^{T}_{ui}\lambda^{+}.$

From (89) we get that

$\displaystyle-F_{u}^{T}(k)R(k+1)$ $\displaystyle=$ $\displaystyle
h\sum_{j=1}^{s}f^{T}_{uj}b_{j}\left(\lambda_{j}-\lambda^{+}\right)+h\sum_{i=1}^{s}b_{i}f^{T}_{ui}\lambda^{+}$
$\displaystyle=$ $\displaystyle
h\sum_{i=1}^{s}f^{T}_{ui}b_{i}\left(\lambda_{i}-\lambda^{+}\right)+h\sum_{i=1}^{s}b_{i}f^{T}_{ui}\lambda^{+}$
$\displaystyle=$ $\displaystyle h\sum_{i=1}^{s}b_{i}f^{T}_{ui}\lambda_{i}.$

The discrete reduced gradient now takes the form

(117)
$\frac{d\phi(X(K))}{du_{n}}=\sum_{k=k_{n-1}}^{k_{n}-1}h(k)\sum_{i=1}^{s}b_{i}f_{u}^{T}(x_{i}(k+1),u_{n})\lambda_{i}(k).$

The quadrature formula for the continuous reduced gradient calculation is

(118) $\displaystyle\left(\frac{d\phi(x(t_{f}))}{du_{n}}\right)_{Q}$
$\displaystyle=$
$\displaystyle\sum_{k=k_{n-1}}^{k_{n}-1}\tilde{h}(k)\sum_{i=1}^{\tilde{s}}\tilde{b}_{i}f_{u}^{T}\left(x\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right),u_{n}\right)\times$
$\displaystyle\lambda\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right).$

Let us now suppose, that the quadrature is calculated such that the discrete
times are the same as during the numerical integration of system equations so

(119) $\tilde{t}(k)=t(k).$

In such case we obtain $\tilde{h}(k)=h(k)$ and

(120)
$\left(\frac{d\phi(x(t_{f}))}{du_{n}}\right)_{Q}=\sum_{k=k_{n-1}}^{k_{n}-1}h(k)\sum_{i=1}^{\tilde{s}}\tilde{b}_{i}f_{u}^{T}\left(x\left(t(k)+\tilde{c}_{i}h(k)\right),u_{n}\right)\lambda\left(t(k)+\tilde{c}_{i}h(k)\right).$

Let us also assume that

(121) $\tilde{s}=s,\ \tilde{b}_{i}=b_{i},\ \tilde{c}_{i}=c_{i},\ i=1,\dots,s.$

In such case (120) becomes

(122)
$\left(\frac{d\phi(x(t_{f}))}{du_{n}}\right)_{Q}=\sum_{k=k_{n-1}}^{k_{n}-1}h(k)\sum_{i=1}^{s}b_{i}f_{u}^{T}\left(x\left(t(k)+c_{i}h(k)\right),u_{n}\right)\lambda\left(t(k)+c_{i}h(k)\right).$

Now the quadrature scheme (122) for continuous reduced gradients is almost
identical to discrete reduced gradients (117). The only difference is that in
(117) $x_{i}(k+1)$ and $\lambda_{i}(k)$ are used instead of
$x(t(k)+c_{i}h(k))$ and $\lambda(t(k)+c_{i}h(k))$. The convergence of (122) is
guaranteed if only the quadrature with coefficients $b_{i},\ c_{i},\
i=1,\dots,s$ satisfies appropriate conditions. We will also justify that the
usage of $x_{i}(k+1)$ and $\lambda_{i}(k)$ instead of $x(t(k)+c_{i}h(k))$ and
$\lambda(t(k)+c_{i}h(k))$ does not break the convergence of (117) to
continuous reduced gradients.

After showing the correspondence between discrete equations and their
continuous counterparts, we want to justify the convergence of the discrete
equations. To obtain the right error order estimated we always assume that the
system function $f(x,u)$ and its partial derivatives are Lipschitz continuous
functions. The convergence analysis for any Runge-Kutta scheme is a
challenging task, so we limit ourselves to the analysis of the RADAU IIA
integration scheme for $s=3$. This scheme has been chosen for the numerical
integration of the system equations because of its advantageous properties,
especially for the integration of DAEs ([5]). The coefficients of the RADAU
IIA scheme for $s=3$ collected in matrices are ([7], p. 74)

(123)
$A=(a_{ij})=\left(\begin{array}[]{ccc}\frac{11}{45}-\frac{7\sqrt{6}}{360}&\frac{37}{225}-\frac{169\sqrt{6}}{1800}&-\frac{2}{225}+\frac{\sqrt{6}}{75}\\\
&&\\\
\frac{37}{225}+\frac{169\sqrt{6}}{1800}&\frac{11}{45}+\frac{7\sqrt{6}}{360}&-\frac{2}{225}-\frac{\sqrt{6}}{75}\\\
&&\\\
\frac{4}{9}-\frac{\sqrt{6}}{36}&\frac{4}{9}+\frac{\sqrt{6}}{36}&\frac{1}{9}\end{array}\right),$
$b=(b_{i})=\left(\begin{array}[]{ccc}\frac{4}{9}-\frac{\sqrt{6}}{36},&\frac{4}{9}+\frac{\sqrt{6}}{36},&\frac{1}{9}\end{array}\right),\
c=(c_{i})=\left(\begin{array}[]{ccc}\frac{2}{5}-\frac{\sqrt{6}}{10},&\frac{2}{5}+\frac{\sqrt{6}}{10},&1\end{array}\right).$

To calculate the convergence order $p$ of the Runge-Kutta scheme the following
conditions imposed on the Runge-Kutta coefficients are checked ([6], p. 208)

(124) $\displaystyle B(p)$ $\displaystyle:$ $\displaystyle\
\forall_{l=1,\dots,p}\ \sum_{i=1}^{s}b_{i}c_{i}^{l-1}=\frac{1}{l},$ (125)
$\displaystyle C(q)$ $\displaystyle:$ $\displaystyle\ \forall_{l=1,\dots,q}\
\forall_{i=1,\dots,s}\ \sum_{j=1}^{s}a_{ij}c_{j}^{l-1}=\frac{c_{i}^{l}}{l},$
(126) $\displaystyle D(r)$ $\displaystyle:$ $\displaystyle\
\forall_{l=1,\dots,r}\ \forall_{j=1,\dots,s}\
\sum_{i=1}^{s}b_{i}c_{i}^{l-1}a_{ij}=\frac{b_{j}}{l}\left(1-c_{j}^{l}\right),$

If the conditions $B(p),C(q),D(r)$ are satisfied for $p\leq 2q+2$ and $p\leq
q+r+1$ then the order of convergence of the Runge-Kutta scheme is $p$ ([6],
Theorem 7.4 p. 208), i.e. the difference between approximated state $x(k)$ and
exact state at time $t(k)$ satisfies

(127) $x(k)-x(t(k))=O(h_{m}^{p}),$

where

(128) $h_{m}=max\\{h(0),\dots,h(K-1)\\}.$

The RADAU IIA scheme satisfies conditions $B(p=5),C(q=3),D(r=2)$ ([7], p. 77)
and so the order of the method is $p=5$.

To justify the convergence of the discrete equations we consider the
unperturbed (nominal) Runge-Kutta scheme

(129) $\displaystyle x_{i}^{n}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h(k)\sum_{j=1}^{s}a_{ij}f\left(x_{j}^{n}(k+1),u(k)\right),\
i=1,\dots,s$ (130) $\displaystyle x^{n}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h(k)\sum_{i=1}^{s}b_{i}f\left(x_{i}^{n}(k),u(k)\right)$

and the perturbed Runge-Kutta scheme

(131) $\displaystyle x_{i}^{p}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h(k)\sum_{j=1}^{s}a_{ij}f\left(x_{j}^{p}(k+1),u(k)\right)+\delta_{i},\
i=1,\dots,s$ (132) $\displaystyle x^{p}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h(k)\sum_{i=1}^{s}b_{i}f\left(x_{i}^{p}(k),u(k)\right)+\delta_{s+1}.$

where $\delta_{i},\ i=1,\dots,s+1$ are perturbations. The estimates of the
errors in the perturbed system are ([7], Theorem 14.3 p. 219)

(133) $\displaystyle x^{p}_{i}(k+1)-x^{n}_{i}(k+1)$ $\displaystyle=$
$\displaystyle O(\delta),\ i=1,\dots,s$ (134) $\displaystyle
x^{p}(k+1)-x^{n}(k+1)$ $\displaystyle=$ $\displaystyle O(\delta),$

where $\delta=max\\{\delta_{1},\dots,\delta_{s+1}\\}$. It should be noted that
the error estimates has been formulated for stiff differential equations ([7],
Theorem 14.3 p. 219). For that reason some additional conditions has been
imposed on the Runge-Kutta scheme. Those conditions concern the invertibility
and coercivity of the matrix $A=(a_{ij})$. It is possible that those
conditions could be relaxed for non-stiff systems, but as the considered
schemes (RADAU IA and RADAU IIA) satisfy them ([7], Theorems 14.3 and 14.5) we
do not discuss this issue further. The error estimates (133)-(134) play a
crucial role for the proposed convergence analysis. To use these estimates
nominal and perturbed schemes must be formulated in an appropriate way in each
case.

Let us now derive the estimate of the global error
$x_{i}(k+1)-x(t(k)+c_{i}h(k))$. The nominal Runge-Kutta scheme is formulated
assuming that the exact state value is known at $t(k)$

$\displaystyle x^{n}_{i}(k+1)$ $\displaystyle=$ $\displaystyle
x(t(k))+h(k)\sum_{j=1}^{s}a_{ij}f\left(x^{n}_{j}(k+1),\ u(k)\right),\
i=1,\dots,s.$

If $C(q)$ holds, the following local error estimate is valid ([6], Lemma 7.5
p. 210)

(136) $x^{n}_{i}(k+1)-x(t(k)+c_{i}h(k))=O(h^{q+1}(k)),\ i=1,\dots,s.$

In this case the perturbed Runge-Kutta scheme is actually the regular Runge-
Kutta scheme, with the approximation of the state $x(k)$ at time $t(k)$ used

$\displaystyle x_{i}(k+1)$ $\displaystyle=$ $\displaystyle
x(k)+h(k)\sum_{j=1}^{s}a_{ij}f\left(x_{j}(k+1),\ u(k)\right)$ $\displaystyle=$
$\displaystyle x(t(k))+h(k)\sum_{j=1}^{s}a_{ij}f\left(x_{j}(k+1),\
u(k)\right)+(x(k)-x(t(k)))$

for $i=1,\dots,s$, so the perturbations are

(138) $\delta_{i}=\delta=x(k)-x(t(k)),\ i=1,\dots,s.$

From (127) we have $\delta=x(k)-x(t(k))=O(h_{m}^{p})$. From (133) we obtain

(139) $x_{i}(k+1)-x^{n}_{i}(k+1)=O(h_{m}^{p}).$

By combining (136) and (139) we obtain the required global error estimate

(140) $\displaystyle x_{i}(k+1)-x(t(k)+c_{i}h(k))$ $\displaystyle=$
$\displaystyle O(h^{q+1}(k))+O(h_{m}^{p})=O(h^{q+1}_{m})+O(h_{m}^{p})$
$\displaystyle=$ $\displaystyle O(h_{m}^{min\\{p,q+1\\}}).$

For the RADAU IIA scheme with $s=3$ we get

(141) $x_{i}(k+1)-x(t(k)+c_{i}h(k))=O(h_{m}^{min\\{5,3+1\\}})=O(h_{m}^{4}).$

Let us now consider the Runge–Kutta scheme for adjoint equations (100)-(101).
The coefficients $\bar{a}_{ij},\ \bar{b}_{i},\ \bar{c}_{i},\ i,j=1,\dots,s$ of
the Runge-Kutta scheme collected in matrices are

(142)
$\bar{A}=(\bar{a}_{ij})=\left(\begin{array}[]{ccc}\frac{11}{45}-\frac{7\sqrt{6}}{360}&\frac{11}{45}+\frac{43\sqrt{6}}{360}&\frac{1}{9}\\\
&&\\\
\frac{11}{45}-\frac{43\sqrt{6}}{360}&\frac{11}{45}-\frac{7\sqrt{6}}{360}&\frac{1}{9}\\\
&&\\\
-\frac{1}{18}+\frac{\sqrt{6}}{18}&-\frac{1}{18}-\frac{\sqrt{6}}{18}&\frac{1}{9}\end{array}\right),$
$\bar{b}=(\bar{b}_{i})=\left(\begin{array}[]{ccc}\frac{4}{9}-\frac{\sqrt{6}}{36},&\frac{4}{9}+\frac{\sqrt{6}}{36},&\frac{1}{9}\end{array}\right),\
\bar{c}=(\bar{c}_{i})=\left(\begin{array}[]{ccc}\frac{3}{5}+\frac{\sqrt{6}}{10},&\frac{3}{5}-\frac{\sqrt{6}}{10},&0\end{array}\right).$

It turns out that the above Runge-Kutta scheme is RADAU IA scheme ([7], p.
73). The RADAU IA scheme satisfies conditions
$B(\bar{p}=5),C(\bar{q}=2),D(\bar{r}=3)$ ([7], p. 77) and so the order of the
method is $\bar{p}=5$.

Let us now define the unperturbed Runge-Kutta scheme as the Runge-Kutta scheme
for continuous adjoint equations with the exact adjoint state known at time
$t(k+1)$

(143)
$\lambda_{i}^{n}(k)=\lambda(t(k+1))+h(k)\sum_{j=1}^{s}\bar{a}_{ij}f_{x}^{T}\left(x\left(t(k)+c_{j}h(k)\right),u(k)\right)\lambda_{j}^{n}(k)$

for $i=1,\dots,s$ and

(144)
$\lambda^{n}(k)=\lambda(t(k+1))+h(k)\sum_{i=1}^{s}\bar{b}_{i}f_{x}^{T}\left(x\left(t(k)+c_{i}h(k)\right),u(k)\right)\lambda_{i}^{n}(k).$

As $\bar{q}=2$ and $\bar{p}=5$ the following local error estimates are valid

(145) $\displaystyle\lambda^{n}_{i}(k)-\lambda(t(k)+c_{i}h(k))$
$\displaystyle=$ $\displaystyle O(h^{\bar{q}+1}(k)),\ i=1,\dots,s\ ,$ (146)
$\displaystyle\lambda^{n}(k)-\lambda(t(k))$ $\displaystyle=$ $\displaystyle
O(h^{\bar{p}+1}(k)).$

The perturbed Runge-Kutta scheme is defined as the discrete adjoint equations
with the exact adjoint state known at time $t(k+1)$

(147)
$\lambda_{i}^{p}(k)=\lambda(t(k+1))+h(k)\sum_{j=1}^{s}\bar{a}_{ij}f_{x}^{T}(x_{j}(k+1),u(k))\lambda_{j}^{p}(k)$

for $i=1,\dots,s$ and

(148)
$\lambda^{p}(k)=\lambda(t(k+1))+h(k)\sum_{i=1}^{s}\bar{b}_{i}f_{x}^{T}(x_{i}(k+1),u(k))\lambda_{i}^{p}(k).$

Using (140), (147)-(148) can be rewritten as

(149)
$\lambda_{i}^{p}(k)=\lambda(t(k+1))+h(k)\sum_{j=1}^{s}\bar{a}_{ij}f_{x}^{T}\left(x\left(t(k)+c_{j}h(k)\right),u(k)\right)\lambda_{j}^{p}(k)+\delta_{i}$

for $i=1,\dots,s$ and

(150)
$\lambda^{p}(k)=\lambda(t(k+1))+h(k)\sum_{i=1}^{s}\bar{b}_{i}f_{x}^{T}\left(x\left(t(k)+c_{i}h(k)\right),u(k)\right)\lambda_{i}^{p}(k)+\delta_{s+1},$

where the perturbations result from differences
$x_{i}(k+1)-x\left(t(k)+c_{i}h(k)\right)$ and so the perturbations satisfy

(151) $\delta_{i}=O(h(k)h^{q+1}_{m}),\ i=1,\dots,s+1.$

From (134) we obtain

(152) $\lambda^{p}(k)-\lambda^{n}(k)=O(h(k)h^{q+1}_{m}).$

By combining (146) and (152) we obtain the following local error estimate for
discrete adjoint equations

(153)
$\lambda^{p}(k)-\lambda(t(k))=O(h^{\bar{p}+1}(k))+O(h(k)h^{q+1}_{m})=O(h(k)h^{min\\{\bar{p},q+1\\}}_{m}).$

If at each step the local error of the numerical integration scheme is
$O(h(k)h^{p}_{m})$, then the global error of the scheme is $O(h^{p}_{m})$
([6], Theorem 3.6 p. 162). Having the local error estimate (153) the global
error estimate for the discrete adjoint equations (91)-(92) is therefore

(154)
$\lambda(k)-\lambda(t(k))=O(h^{min\\{\bar{p},q+1\\}}_{m})=O(h^{\bar{p}_{d}}_{m}),$

where $\bar{p}_{d}=min\\{\bar{p},q+1\\}$. For RADAU IIA scheme we obtain
$\bar{p}_{d}=min\\{5,3+1\\}=4$. Let us now derive the global error estimate of
$\lambda_{i}(k)-\lambda(t(k)+c_{i}h(k))$. The unperturbed system is (143) and
the perturbed system is (91). The perturbed system can be rewritten as

(155)
$\lambda_{i}(k)=\lambda(k)+h(k)\sum_{j=1}^{s}\bar{a}_{ij}f_{x}^{T}\left(x_{i}(k+1),u(k)\right)\lambda_{j}(k)+\delta_{i}$

for $i=1,\dots,s$. The perturbations result from differences
$\lambda(k)-\lambda(t(k+1))$ and $x_{i}(k+1)-x\left(t(k)+c_{i}h(k)\right)$ and
so the perturbations satisfy

(156)
$\displaystyle\delta_{i}=O(h^{\bar{p}_{d}}_{m})+O(h(k)h^{q+1}_{m})=O(h^{\bar{p}_{d}}_{m})+O(h^{q+2}_{m})=O(h_{m}^{min\\{\bar{p}_{d},q+2\\}}).$

From (133) we obtain

(157) $\lambda_{i}(k)-\lambda_{i}^{n}(k)=O(h_{m}^{min\\{\bar{p}_{d},q+2\\}}).$

By combining (145) and (157) we obtain

$\displaystyle\lambda_{i}(k)-\lambda(t(k)+c_{i}h(k))$ $\displaystyle=$
$\displaystyle O(h^{\bar{q}+1}(k))+O(h_{m}^{min\\{\bar{p}_{d},q+2\\}})$
$\displaystyle=$ $\displaystyle
O(h_{m}^{min\\{\bar{q}+1,\bar{p}_{d},q+2\\}})=O(h_{m}^{\bar{q}_{d}+1}),$

where $\bar{q}_{d}=min\\{\bar{q}+1,\bar{p}_{d},q+2\\}-1$. For RADAU IIA scheme
we obtain $\bar{q}_{d}=min\\{2+1,4,3+2\\}-1=2$.

Let us finally derive the order of the discrete reduced gradients (117). If
the quadrature coefficients $\tilde{b}_{i},\tilde{c}_{i},\ i=1,\dots,s$
satisfy the condition $B(\tilde{p})$, then the quadrature integrates exactly
the polynomials of a degree up to $\tilde{p}-1$ on an interval $[0,1]$ ([5],
p.15). We consider the quadrature with coefficients
$\tilde{b}_{i}=b_{i},\tilde{c}_{i}=c_{i},\ i=1,\dots,s$, which satisfy the
condition $B(\tilde{p}=5)$. The local error resulting from the usage of the
quadrature (122) over single integration step can therefore be estimated as

$\displaystyle\int_{t(k)}^{t(k+1)}$ $\displaystyle
f_{u}^{T}(x(t),u_{n})\lambda(t)dt$ $\displaystyle=$
$\displaystyle\int_{0}^{1}f_{u}^{T}(x(t(k)+\tau h(k),u_{n})\lambda(t(k)+\tau
h(k))h(k)d\tau$ $\displaystyle=$ $\displaystyle
h(k)\int_{0}^{1}f_{u}^{T}(x(t(k)+\tau h(k),u_{n})\lambda(t(k)+\tau h(k))d\tau$
$\displaystyle=$ $\displaystyle
h(k)\sum_{i=1}^{s}b_{i}f_{u}^{T}\left(x\left(t(k)+c_{i}h(k)\right),u_{n}\right)\lambda\left(t(k)+c_{i}h(k)\right)+$
$\displaystyle O(h^{\tilde{p}+1}(k)).$

From (141) and (5) we obtain

$\displaystyle h(k)$
$\displaystyle\sum_{i=1}^{s}b_{i}f_{u}^{T}(x_{i}(k+1),u_{n})\lambda_{i}(k)$
$\displaystyle=$ $\displaystyle
h(k)\sum_{i=1}^{s}b_{i}f_{u}^{T}\left(x\left(t(k)+c_{i}h(k)\right),u_{n}\right)\lambda\left(t(k)+c_{i}h(k)\right)$
$\displaystyle+O(h(k)h_{m}^{min\\{q+1,\bar{q}_{d}+1\\}}).$

The local error of the discrete reduced gradients is therefore

$\displaystyle\int_{t(k)}^{t(k+1)}$ $\displaystyle
f_{u}^{T}(x(t),u_{n})\lambda(t)dt-h(k)\sum_{i=1}^{s}b_{i}f_{u}^{T}(x_{i}(k+1),u_{n})$
$\displaystyle=$ $\displaystyle
O(h^{\tilde{p}+1}(k))+O(h(k)h_{m}^{min\\{q+1,\bar{q}_{d}+1\\}})=$
$\displaystyle=$ $\displaystyle
O(h(k)h_{m}^{min\\{\tilde{p},q+1,\bar{q}_{d}+1\\}})=O(h(k)h_{m}^{\tilde{p}_{d}}),$

where $\tilde{p}_{d}=min\\{\tilde{p},q+1,\bar{q}_{d}+1\\}$. For RADAU IIA
scheme we obtain $\tilde{p}_{d}=min\\{5,3+1,2+1\\}=3$. If the local error of
the discrete reduced gradient is $O(h(k)h_{m}^{\tilde{p}_{d}})$ then the
global error can be estimated as ([6], Theorem 3.6 p. 162)

(162)
$\frac{d\phi(X(K))}{du_{n}}-\frac{d\phi(x(t_{f}))}{du_{n}}=O(h_{m}^{\tilde{p}_{d}})$

This result confirms that the discrete reduced gradients (117) converge to the
continuous reduced gradients (62) with an order at least $\tilde{p}_{d}=3$
(for RADAU IIA scheme). The discrete reduced gradients provide therefore an
efficient and reliable method for approximation of reduced gradients.

##  6 Conclusions

The paper presents the novel computational approach to hybrid optimal control
problems with sliding modes. Our computational method is based on a
Runge–Kutta scheme for integrating system and adjoint equations. We show that
in case of ODEs system equations, discrete adjoint equations correspond to the
discretized continuous adjoint equations and on that basis we justify the
convergence of the discrete adjoint equations to their continuous counterparts
for step sizes approaching zero, we also give an estimate of the convergence
order. In the second part of the paper we derive a similar result for index 2
DAEs describing the system motion in a sliding mode. In the second part we
also present several numerical results of the application of the proposed
computational procedure to optimal control problems with hybrid systems.

## References

  * [1] U. Ascher and L. Petzold, Computer methods for ordinary differential equations and differential-algebraic equations , Society for Industrial and Applied Mathematics, Philadelphia, (1998).
  * [2] A. Bryson and Y. Ho, Applied Optimal Control , Hemisphere, New York, (1975).
  * [3] L. Dieci and L. Lopez, ”Sliding Motion in Filippov Differential Systems: Theoretical Results and a Computational Approach”, SIAM J. Numer. Anal., vol. 47(3), pp. 2023–2051, (2009).
  * [4] W. Hager, ”Runge–Kutta methods in optimal control and the transformed adjoint equations”, Numer. Math., vol. 87, pp. 247–282, (2000).
  * [5] E. Hairer, Ch. Lubich and M. Roche, The Numerical Solution of Differential–Algebraic Equations by Runge–Kutta Methods, Lecture Notes in Mathematics, vol. 1409, Springer–Verlag, Berlin, Heidelberg, (1989).
  * [6] E. Hairer, S. Norsett and G. Wanner, Solving Ordinary Differential Equations I, Springer–Verlag, Berlin, Heidelberg, (1993).
  * [7] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Springer–Verlag, Berlin, Heidelberg, (1996).
  * [8] J. Lygeros, K. H. Johansson, S. S. Sastry and M. Egerstedt, ”On the existence of executions of hybrid automata”, in Proceed. of the 38th IEEE CDC, Phoenix, Arizona, 1999, pp. 2249–2254, (1999).
  * [9] I. Mynttinen, A. Hoffmann, E. Runge, P. Li, ”Smoothing and regularization strategies for optimization of hybrid dynamic systems”, Optim Eng., Vol. 16, pp. 541–569, (2015).
  * [10] A. Pakniyat and P.E. Caines, ”The Hybrid Minimum Principle in the Presence of Switching Costs ”, Proceedings of the 52nd IEEE Conference on Decision and Control, December 10-13, Florence, Italy, pp. 3831–3836, (2013).
  * [11] O. Pironneau, E. Polak, ”Consistent Approximations and Appropriate Functions and Gradients in Optimal Control”, SIAM J. Control and Optimization, vol. 41, pp. 487–510, (2002).
  * [12] R. Pytlak, Numerical Methods for Optimal Control Problem With State Constraints, Lecture Notes in Mathematics, vol. 1707, Springer–Verlag, (1999), .
  * [13] R. Pytlak, ”Numerical procedure for optimal control of higher index DAEs”, J. Discret. Dyn. Nat. Soc. (USA), vol. 29, pp. 1–24, (2011).
  * [14] R. Pytlak and D. Suski, ”Algorithms for optimal control of hybrid systems with sliding motion”, arXiv:2101.04754, (2021).
  * [15] R. Pytlak and D. Suski, ”Numerical procedure for optimal control of hybrid systems with sliding modes, Part 2”, submitted to arXiv, (2021).
  * [16] R. Pytlak and D. Suski, ”On solving hybrid optimal control problems with higher index DAEs”, Optimization Methods & Software vol. 32, pp. 940–962, (2017).
  * [17] R. Pytlak and T. Zawadzki, ”On solving optimal control problems with higher index DAEs”, Optimization Methods & Software vol. 29, pp. 1139–1162, (2014).
  * [18] A. van der Schaft and H. Schumacher, An Introduction to Hybrid Dynamical Systems, Springer-Verlag London, (2000).
  * [19] M. S. Shaikh, ”Optimal Control of Hybrid Systems: Theory and Algorithms”, Ph.D. dissertation, McGill University, Montreal, May 2004. [Online]. Available: http://digitool.library.mcgill.ca/R/?func=dbin-jump-full&object_id=85095&local_base=GEN01-MCG02, (2004)
  * [20] M. S. Shaikh, and P.E. Caines, ”On the Hybrid Optimal Control Problem: Theory and Algorithms”, IEEE Transactions on AC, Vol. 52, No. 9, pp. 1587–1603, (2007).
  * [21] M. S. Shaikh, and P.E. Caines, ”Correction to ’On the Hybrid Optimal Control Problem: Theory and Algorithms’”, IEEE Transactions on AC, Vol. 54, No. 6, p. 1, (2009).
  * [22] H. Schiller, M. Arnold, ”Convergence of continuous approximations for discontinuous ODEs”, Applied Numerical Mathematics, Vol. 62, pp. 1503–1514, (2012).
  * [23] A. Schwartz and E. Polak ”Consistent approximations for optimal control problems based on Runge–Kutta integration”, SIAM J. Control and Optimization, vol. 34, pp. 1235–1269, (1998).
  * [24] D. Stewart and M. Anitescu, ”Optimal control of systems with discontinuous differential equations”, Num. Math., vol. 114, pp. 653-695, (2010).
  * [25] H. Sussmann, ”A maximum principle for hybrid optimal control problems”, Proceed. of the 38th IEEE CDC, Phoenix, Arizona, 1999, pp. 425–430, (1999).
  * [26] F. Taringoo and P.E. Caines, ”Gradient-Geodesic HMP Algorithms for the Optimization of Hybrid Systems Based on the Geometry of Switching Manifolds”, Proceedings of the 49th IEEE Conference on Decision and Control, December 15-17, Atlanta, GA, USA, pp. 1534–1539, (2010).
  * [27] F. Taringoo and P.E. Caines, ”On the Extension of the Hybrid Minimum Principle to Riemannian Manifolds”, Proceedings of the 50th IEEE Conference on Decision and Control, December 12-15, Orlando, Fl, USA, pp. 3301–3306, (2011).
  * [28] F. Taringoo and P.E. Caines, ”Newton-Geodesic HMP Algorithms for the Optimization of Hybrid Systems and the Geometric Properties of Hybrid Value Functions”, Proceedings of the 51st IEEE Conference on Decision and Control, December 10-13, Maui, Hawaii, USA, pp. 4211–4216, (2012).
  * [29] F. Taringoo and P.E. Caines, ”On the Optimal Control of Hybrid Systems On Lie Groups and the Exponential Gradient HMP Algorithm”, Proceedings of the 52nd IEEE Conference on Decision and Control, December 10-13, Florence, Italy, pp. 2653–2658, (2013).
  * [30] H. Witsenhausen, ”A class of hybrid-state continuous-time dynamic systems”, IEEE Trans. on Autom. Control, vol. 11, pp. 161–167, (1966).