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

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

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. A sliding mode is coped with differential-
algebraic equations (DAEs) and that guarantees accurate tracking of the
sliding motion surface. In the second 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 DAEs. We show that the discrete adjoint state trajectories
converge to their continuous counterparts. Next, we describe the application
of the proposed procedure to three optimal control problems. The first problem
concerns optimal control of a simple mechanical system with dry friction. The
second problem is related to the planning of a haemodialysis process. The
third problem concerns the optimal steering of a racing car.

###### keywords:

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

65L80, 49M37, 65K10

##  1 Introduction

The second part of the paper continues the description of the numerical
procedure for optimal control problems for hybrid system exhibiting sliding
motion behavior. The procedure is the implementation of the algorithm for
these problems introduced in [9]. In the first part of the paper we
concentrated on describing the efficient way of evaluating reduced gradients
of the optimization problem. We paid the attention to the case the hybrid
system does not enter the sliding mode and its dynamics is described by ODEs.
In the second part we supplement our considerations by analyzing the sliding
motion case.

Let us consider the hybrid system evolving according to the following scheme.
For $t\in\left[t_{0},t_{t}\right]$ the system evolves according to ODEs

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

at a switching time $t_{t}$ system state reaches the switching surface

(2) $g(x(t_{t}))=0$

and enters the sliding mode. For $t\in\left[t_{t},t_{f}\right]$ the state
trajectory is the solution of DAEs

(3) $\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle
f_{F}(x,u)+g_{x}^{T}(x)z=f^{2}(x,z,u)$ (4) $\displaystyle 0$ $\displaystyle=$
$\displaystyle g(x)$

We now aim at finding the reduced gradient formula for the endpoint functional
([10])

(5) $\bar{F}_{0}(u)=\phi(x^{u}(t_{f}))$

with respect to controls. This formula can be found with the help of adjoint
equations. We now recall the appropriate formulation of adjoint equations in
the considered case. 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}$ ([9])

(6) $\displaystyle\phi_{x}^{T}(x(t_{f}))+\lambda_{f}(t_{f})$ $\displaystyle=$
$\displaystyle\nu_{1}g_{x}^{T}(x(t_{f}))$ (7) $\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.

Having $\lambda_{f}(t_{f})$ and $\lambda_{g}(t_{f})$ we then solve DAEs (the
meaning of $t_{t}^{-}$ and $t_{t}^{+}$ is explained in [10])

$\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))$
(10) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\lambda_{f}^{T}(t)g_{x}^{T}(x(t)),\ t\in[t_{t}^{+},t_{f}].$

backwards in time. 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$

(11) $\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

(13)
$(\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 of the functional $\bar{F}_{0}(u)$ with respect
to the control variation $d$ can be calculated from

$\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 (Section 2) we will analyze the discrete time versions of
equations (1)–(10) and the approximation to the second term in (1). The
present analysis complements the analysis given in the first part of the paper
[10] which concerned equations (13) and the first term in (1). In addition in
Section 3 we will examine how the discretization of adjoint equations
influences the jump conditions (1). Section 4 presents the efficiency of the
proposed numerical procedure in solving three optimal control problems.

##  2 Numerical calculation of reduced gradients

Let us consider the control system described by higher index differential-
algebraic equations, we focus on index 2 DAEs in Hessenberg form ([3])

(15) $\displaystyle x^{\prime}(t)$ $\displaystyle=$ $\displaystyle
f(x(t),z(t),u(t)),$ (16) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
g(x(t)),$

where $x\in\mathbb{R}^{n_{x}}$ is a differential state and
$z\in\mathbb{R}^{n_{z}}$ is an algebraic state. The differentiation index of a
system (15)-(16) is 2 provided that the matrix $g_{x}(x)f_{z}(x,z,u)$ is
nonsingular ([3]). To integrate DAEs, the following Runge-Kutta scheme can be
used ([3])

(17) $\displaystyle x_{i}^{\prime}(k+1)$ $\displaystyle=$ $\displaystyle
f\left(x_{i}(k+1),\ z_{i}(k+1),u(k)\right),$ (18) $\displaystyle 0$
$\displaystyle=$ $\displaystyle g\left(x_{i}(k+1)\right),$ (19) $\displaystyle
x_{i}(k+1)$ $\displaystyle=$ $\displaystyle
x(k)+h(k)\sum_{j=1}^{s}a_{ij}x_{j}^{\prime}(k+1),$ (20) $\displaystyle
z_{i}(k+1)$ $\displaystyle=$ $\displaystyle
z(k)+h(k)\sum_{j=1}^{s}a_{ij}z_{j}^{\prime}(k+1)$

for $i=1,...,s$ and

(21) $\displaystyle x(k+1)$ $\displaystyle=$ $\displaystyle
x(k)+h(k)\sum_{i=1}^{s}b_{i}x_{i}^{\prime}(k+1),$ (22) $\displaystyle z(k+1)$
$\displaystyle=$ $\displaystyle
z(k)+h(k)\sum_{i=1}^{s}b_{i}z_{i}^{\prime}(k+1).$

It is possible to get rid of variables $x_{i}^{\prime}(k+1)$ and
$z_{i}^{\prime}(k+1)$ from (17)-(22), provided that the Runge-Kutta matrix
$A=(a_{ij})$ is invertible. In such case, from (20) we have

(23)
$z_{i}^{\prime}(k+1)=\frac{1}{h(k)}\sum_{j=1}^{s}a_{ij}^{-}(z_{j}(k+1)-z(k)),\
i=1,\dots,s,$

where $a_{ij}^{-}$ are the coefficients of the matrix $A^{-1}$. Now the system
(17)-(22) can be rewritten as

(24) $\displaystyle x_{i}(k+1)$ $\displaystyle=$ $\displaystyle
x(k)+h(k)\sum_{j=1}^{s}a_{ij}f\left(x_{j}(k+1),\ z_{j}(k+1),u(k)\right),$ (25)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle g\left(x_{i}(k+1)\right),$

for $i=1,...,s$ and

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

where

(28) $b_{i}^{-}=\sum_{i=1}^{s}b_{j}a^{-}_{ji},\ i=1,\dots,s.$

At each step of the Runge–Kutta scheme the nonlinear system (24)-(25) is first
solved for variables $x_{i}(k+1),z_{i}(k+1)\ i=1,\dots,s$ and then (26)-(2)
are used to calculate $x(k+1)$ and $z(k+1)$. Let us notice that $z(k)$ is
required to calculate $z(k+1)$, but not for solving (24)-(25).

The adjoint equations for the control system (15)-(16) are ([9]):

(29) $\displaystyle\lambda_{f}^{\prime}(t)$ $\displaystyle=$ $\displaystyle-
f_{x}^{T}(x(t),z(t),u(t))\lambda_{f}(t)-g_{x}^{T}(x(t))\lambda_{g}(t),$ (30)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle-
f_{z}^{T}(x(t),z(t),u(t))\lambda_{f}(t).$

For the adjoint DAEs (29)-(30) to be properly defined, the term
$f_{z}^{T}(x,z,u)\lambda_{f}$ should be differentiable. The algebraic state
$z(t)$ and the control $u(t)$ are not continuous, so to keep the
differentiability we assume that $f_{z}^{T}(x(t),z(t),u(t))$ does not depend
neither on $z(t)$ nor on $u(t)$. Let us notice that the adjoint equations
formulated for the sliding motion satisfy that condition. From now we will
consider the adjoint equations of the form

(31) $\displaystyle\lambda_{f}^{\prime}(t)$ $\displaystyle=$ $\displaystyle-
f_{x}^{T}(x(t),z(t),u(t))\lambda_{f}(t)-g_{x}^{T}(x(t))\lambda_{g}(t),$ (32)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle-
f_{z}^{T}(x(t))\lambda_{f}(t).$

The adjoint equations are index 2 DAEs ([9]), so the Runge-Kutta scheme can be
used to numerically integrate them ([3]). The Runge-Kutta scheme for adjoint
equations (31)-(32) is

$\displaystyle\lambda_{fi}(k)=$
$\displaystyle\lambda_{f}(k+1)-\bar{h}(k+1)\sum_{j=1}^{\bar{s}}\bar{a}_{ij}\left[-f_{x}^{T}\left(x\left(\bar{t}_{j}(k)\right),z\left(\bar{t}_{j}(k)\right),\bar{u}(k+1)\right)\lambda_{fj}(k)\right.$
(33) $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.-g_{x}^{T}\left(x\left(\bar{t}_{j}(k)\right)\right)\lambda_{gj}(k)\right]=$
$\displaystyle\lambda_{f}(k+1)+\bar{h}(k+1)\sum_{j=1}^{\bar{s}}\bar{a}_{ij}\left[f_{x}^{T}\left(x\left(\bar{t}_{j}(k)\right),z\left(\bar{t}_{j}(k)\right),\bar{u}(k+1)\right)\lambda_{fj}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x\left(\bar{t}_{j}(k)\right)\right)\lambda_{gj}(k)\right],$
(34) $\displaystyle
0=-f^{T}_{z}\left(x\left(\bar{t}_{i}(k)\right)\right)\lambda_{fi}(k)$

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

$\displaystyle\lambda_{f}(k)=$
$\displaystyle\lambda_{f}(k+1)-\bar{h}(k+1)\sum_{i=1}^{\bar{s}}\bar{b}_{i}\left[-f_{x}^{T}\left(x\left(\bar{t}_{i}(k)\right),z\left(\bar{t}_{i}(k)\right),\bar{u}(k+1)\right)\lambda_{fi}(k)\right.$
(35) $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.-g_{x}^{T}\left(x\left(\bar{t}_{i}(k)\right)\right)\lambda_{gi}(k)\right]=$
$\displaystyle\lambda_{f}(k+1)+\bar{h}(k+1)\sum_{i=1}^{\bar{s}}\bar{b}_{i}\left[f_{x}^{T}\left(x\left(\bar{t}_{i}(k)\right),z\left(\bar{t}_{i}(k)\right),\bar{u}(k+1)\right)\lambda_{fi}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x\left(\bar{t}_{i}(k)\right)\right)\lambda_{gi}(k)\right],$
$\displaystyle\lambda_{g}(k)=$ (36)
$\displaystyle\lambda_{g}(k+1)+\sum_{i=1}^{\bar{s}}\sum_{j=1}^{\bar{s}}\bar{b}_{i}\bar{a}_{ij}^{-}(\lambda_{gj}(k)-\lambda_{g}(k+1))=$
$\displaystyle\lambda_{g}(k+1)+\sum_{i=1}^{\bar{s}}\bar{b}_{i}^{-}(\lambda_{gi}(k)-\lambda_{g}(k+1)),$

where for the sake of the shorter notation we denote

(37) $\bar{t}_{i}(k)=t(k+1)-\bar{c}_{i}\bar{h}(k+1),\ i=1,\dots,s.$

At each step of the Runge-Kutta scheme, the system of equations (33)-(34) is
first solved for variables $\lambda_{fi}(k),\lambda_{gi}(k),\
i=1,\dots,\bar{s}$ and then (35)-(36) are used to calculate $\lambda_{f}(k)$
and $\lambda_{g}(k)$.

In the case of the sliding motion the differential of the endpoint functional
$\phi(t_{f})$ is of the form

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

Two important remarks need to be emphasized. First we assume that $f_{u}(x,u)$
does not depend on $z$. That is another limitation lied on the form of the
function $f(x,z,u)$, but it is satisfied by the sliding motion equations.
Second, the integral depends only on the adjoint variable $\lambda_{f}$, but
not on $\lambda_{g}$. Under that conditions the reduced gradients for DAEs are
calculated the same way as it was for ODEs (see [10])

(39)
$\frac{d\phi(x(t_{f}))}{du_{n}}\simeq\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)\lambda_{f}\left(\tilde{t}(k)+\tilde{c}_{i}\tilde{h}(k)\right),$

where for the discrete steps $k=k_{n-1},\ldots,k_{n}-1$ the control is
$u_{n}$.

Similarly to the ODEs case presented in [10], we utilize discrete adjoint
equations and discrete reduced gradients to avoid the necessity of calculation
of states and adjoint variables at arbitrary time moments. The Runge-Kutta
scheme (24)-(26) can be rewritten to a vector form (discrete step argument
omitted)

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

If we now define the augmented state vector $X(k)$ as

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

then (40) can be presented in a form of the implicit discrete time state
equation [8] (see also [12])

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

The partial derivatives matrices of the discrete state equation are

(43)
$F_{X^{+}}(k)=\left(\begin{array}[]{ccccccc}I-ha_{11}f_{x1}&-ha_{11}f_{z1}&\ldots&-ha_{1s}f_{xs}&-ha_{1s}f_{zs}&0\\\
-g_{x1}&0&\ldots&0&0&0\\\ \vdots&\vdots&&\vdots&\vdots&\vdots\\\
-ha_{s1}f_{x1}&-ha_{s1}f_{z1}&\ldots&I-ha_{ss}f_{xs}&-ha_{ss}f_{zs}&0\\\
0&0&\ldots&-g_{xs}&0&0\\\
-hb_{1}f_{x1}&-hb_{1}f_{z1}&\ldots&-hb_{s}f_{xs}&-hb_{s}f_{zs}&I\end{array}\right)$

and

(44) $F_{X}(k)=\left(\begin{array}[]{ccccccc}0&0&\ldots&0&0&-I\\\
0&0&\ldots&0&0&0\\\ \vdots&\vdots&&\vdots&\vdots&\vdots\\\
0&0&\ldots&0&0&-I\\\ 0&0&\ldots&0&0&0\\\ 0&0&\ldots&0&0&-I\end{array}\right),$

where for the sake of a shorter notation we denote

(45) $f_{xi}=f_{x}(x_{i},z_{i},u),f_{zi}=f_{z}(x_{i}),g_{xi}=g_{x}(x_{i}).$

Let us denote

$\displaystyle\Lambda(k)$ $\displaystyle=$
$\displaystyle\left(l_{f1}(k)^{T},l_{g1}(k)^{T},\ldots,l_{fs}(k)^{T},l_{gs}(k)^{T},\lambda_{f}(k)^{T}\right)^{T},$
$\displaystyle\Lambda(k+1)$ $\displaystyle=$
$\displaystyle\left(l_{f1}^{+}(k)^{T},l_{g1}^{+}(k)^{T},\ldots,l_{fs}^{+}(k)^{T},l_{gs}^{+}(k)^{T},\lambda_{f}^{+}(k)^{T}\right)^{T},$
$\displaystyle R(k)$ $\displaystyle=$
$\displaystyle\left(r_{f1}(k)^{T},r_{g1}(k)^{T},\ldots,r_{fs}(k)^{T},r_{gs}(k)^{T},r_{f}(k)^{T}\right)^{T}.$

Now the discrete adjoint equations ([8])

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

at a discrete time step $k$ take the form

(48)
$\left(\begin{array}[]{cccccc}I-ha_{11}f^{T}_{x1}&-g^{T}_{x1}&\ldots&-ha_{s1}f^{T}_{x1}&0&-hb_{1}f_{x1}^{T}\\\
-ha_{11}f^{T}_{z1}&0&\ldots&-ha_{s1}f^{T}_{z1}&0&-hb_{1}f_{z1}^{T}\\\
\vdots&\vdots&&\vdots&\vdots&\vdots\\\
-ha_{1s}f^{T}_{xs}&0&\ldots&I-ha_{ss}f^{T}_{xs}&-g^{T}_{xs}&-hb_{s}f_{xs}^{T}\\\
-ha_{1s}f^{T}_{zs}&0&\ldots&-ha_{ss}f^{T}_{zs}&0&-hb_{s}f_{zs}^{T}\\\
0&0&\ldots&0&0&I\\\ \end{array}\right)\left(\begin{array}[]{c}r_{f1}\\\
r_{g1}\\\ \vdots\\\ r_{fs}\\\ r_{gs}\\\
r_{f}\end{array}\right)=\left(\begin{array}[]{c}l^{+}_{f1}\\\ l^{+}_{g1}\\\
\vdots\\\ l^{+}_{fs}\\\ l^{+}_{gs}\\\ \lambda^{+}_{f}\end{array}\right)$

(49) $\left(\begin{array}[]{c}l_{f1}\\\ l_{g1}\\\ \vdots\\\ l_{fs}\\\
l_{gs}\\\
\lambda_{f}\end{array}\right)=-\left(\begin{array}[]{cccccc}0&0&\ldots&0&0&0\\\
0&0&\ldots&0&0&0\\\ \vdots&\vdots&&\vdots&\vdots&\vdots\\\ 0&0&\ldots&0&0&0\\\
0&0&\ldots&0&0&0\\\
-I&0&\ldots&-I&0&-I\end{array}\right)\left(\begin{array}[]{c}r_{f1}\\\
r_{g1}\\\ \vdots\\\ r_{fs}\\\ r_{gs}\\\ r_{f}\end{array}\right)$

If we rewrite (48) as a system of equations we obtain

(50) $\displaystyle r_{fi}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{s}ha_{ji}f_{xi}^{T}r_{fj}+g_{xi}^{T}r_{gi}+hb_{i}f_{xi}^{T}r_{f}+l^{+}_{fi}$
(51) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{s}ha_{ji}f_{zi}^{T}r_{fj}+hb_{i}f_{zi}^{T}r_{f}+l_{gi}^{+}$

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

(52) $r_{f}=\lambda_{f}^{+}$

Using (52), (50)-(51) can be written as

(53) $\displaystyle r_{fi}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{s}ha_{ji}f_{xi}^{T}r_{fj}+g_{xi}^{T}r_{gi}+hb_{i}f_{xi}^{T}\lambda_{f}^{+}+l^{+}_{fi}$
(54) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{s}ha_{ji}f_{zi}^{T}r_{fj}+hb_{i}f_{zi}^{T}\lambda_{f}^{+}+l_{gi}^{+}$

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

(55) $\displaystyle l_{fi}$ $\displaystyle=$ $\displaystyle 0$ (56)
$\displaystyle l_{gi}$ $\displaystyle=$ $\displaystyle 0$

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

(57) $\displaystyle\lambda_{f}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{s}r_{fi}+r_{f}.$

Using (52), (57) gives

(58) $\displaystyle\lambda_{f}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{s}r_{fi}+\lambda_{f}^{+}.$

From (55)-(56) we have that
$l_{f1}(k)=0,l_{g1}(k)=0,\ldots,l_{fs}(k)=0,l_{gs}(k)=0$ for steps
$k=0,...,K-1$. During the analysis we also assume that
$l_{f1}(K)=0,l_{g1}(K)=0,\ldots,l_{fs}(K)=0,l_{gs}(K)=0$. Under that
assumption (53)-(54) is equivalent to

(59) $\displaystyle r_{fi}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{s}ha_{ji}f_{xi}^{T}r_{fj}+g_{xi}^{T}r_{gi}+hb_{i}f_{xi}^{T}\lambda_{f}^{+}$
(60) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{s}ha_{ji}f_{zi}^{T}r_{fj}+hb_{i}f_{zi}^{T}\lambda_{f}^{+}$

Let us now introduce the auxiliary variables

(61) $\displaystyle\lambda_{fi}$ $\displaystyle=$
$\displaystyle\lambda_{f}^{+}+\sum_{j=1}^{s}\frac{a_{ji}}{b_{i}}r_{fj},$ (62)
$\displaystyle\lambda_{gi}$ $\displaystyle=$
$\displaystyle\frac{r_{gi}}{hb_{i}}$

for $i=1,\ldots,s$. (59)-(60) can then be transformed to

(63) $\displaystyle r_{fi}$ $\displaystyle=$ $\displaystyle
hb_{i}\left(f_{xi}^{T}\lambda_{fi}+g_{xi}^{T}\lambda_{gi}\right),$ (64)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
hb_{i}f_{zi}^{T}\lambda_{fi}.$

If we put (63) into (61) we obtain

(65) $\displaystyle\lambda_{fi}$ $\displaystyle=$
$\displaystyle\lambda_{f}^{+}+h\sum_{j=1}^{s}\frac{a_{ji}b_{j}}{b_{i}}\left(f_{xj}^{T}\lambda_{fj}+g_{xj}^{T}\lambda_{gj}\right).$

Dividing (64) by $hb_{i}$ results in

(66) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
f_{zi}^{T}\lambda_{fi}.$

If we put (63) into (58) we obtain

(67) $\displaystyle\lambda_{f}$ $\displaystyle=$
$\displaystyle\lambda_{f}^{+}+h\sum_{i=1}^{s}b_{i}\left(f_{xi}^{T}\lambda_{fi}+g_{xi}^{T}\lambda_{gi}\right).$

Let us rewrite equations (65)-(67) by introducing the discrete step argument

(68) $\displaystyle\lambda_{fi}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(k+1)+h(k)\sum_{j=1}^{s}\frac{a_{ji}b_{j}}{b_{i}}\left[f_{x}^{T}\left(x_{j}\left(k+1\right),z_{j}\left(k+1\right),u(k)\right)\lambda_{fj}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x_{i}\left(k+1\right)\right)\lambda_{gj}(k)\right],$
(69) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
f^{T}_{z}\left(x_{i}\left(k+1\right)\right)\lambda_{fi}(k)$

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

(70) $\displaystyle\lambda_{f}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(k+1)+h(k)\sum_{i=1}^{s}b_{i}\left[f_{x}^{T}\left(x_{j}\left(k+1\right),z_{i}\left(k+1\right),u(k)\right)\lambda_{fi}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x_{i}\left(k+1\right)\right)\lambda_{gi}(k)\right]$

and augment them by the equation for $\lambda_{g}(k)$

(71) $\displaystyle\lambda_{g}(k)$ $\displaystyle=$
$\displaystyle\lambda_{g}(k+1)+\sum_{i=1}^{s}b_{i}^{-}(\lambda_{gi}(k)-\lambda_{g}(k+1)).$

We now consider the Runge-Kutta scheme (33)-(36) under the assumption that the
discrete steps are the same as in the forward scheme and the Runge-Kutta
scheme coefficients satisfy
$\bar{a}_{ij}=\frac{a_{ji}b_{j}}{b_{i}},\bar{b}_{i}=b_{i},\bar{c}_{i}=1-c_{i}$

(72) $\displaystyle\lambda_{fi}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(k+1)+h(k)\sum_{j=1}^{s}\frac{a_{ji}b_{j}}{b_{i}}\left[f_{x}^{T}\left(x\left(t_{j}(k)\right),z\left(t_{j}(k)\right),u(k)\right)\lambda_{fj}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x\left(t_{j}(k)\right)\right)\lambda_{gj}(k)\right],$
(73) $\displaystyle 0$ $\displaystyle=$
$\displaystyle-f^{T}_{z}\left(x\left(t_{i}(k)\right)\right)\lambda_{fi}(k)$

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

(74) $\displaystyle\lambda_{f}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(k+1)+h(k)\sum_{i=1}^{s}b_{i}\left[f_{x}^{T}\left(x\left(t_{i}(k)\right),z\left(t_{i}(k)\right),u(k)\right)\lambda_{fi}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x\left(t_{i}(k)\right)\right)\lambda_{gi}(k)\right],$
(75) $\displaystyle\lambda_{g}(k)$ $\displaystyle=$
$\displaystyle\lambda_{g}(k+1)+\sum_{i=1}^{s}b_{i}^{-}(\lambda_{gi}(k)-\lambda_{g}(k+1)),$

where

(76) $t_{i}(k)=t(k)+c_{i}h(k).$

Now the Runge-Kutta scheme (72)-(75) for adjoint equations is almost identical
to discrete adjoint equations (68)-(71). The only difference is that in
discrete equations $x_{i}(k+1)$ and $z_{i}(k+1)$ are used instead of
$x(t(k)+c_{i}h(k))$ and $z(t(k)+c_{i}h(k))$. The convergence of (72)-(75) 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)$ and $z_{i}(k+1)$ instead of
$x(t(k)+c_{i}h(k))$ and $z(t(k)+c_{i}h(k))$ does not destroy the convergence
of (68)-(71) solutions to continuous adjoint trajectory.

The partial derivative $F_{u}(k)$ is

(77) $F_{u}(k)=\left(\begin{array}[]{c}-h\sum_{j=1}^{s}a_{1j}f_{uj}\\\ 0\\\
\vdots\\\ -h\sum_{j=1}^{s}a_{sj}f_{uj}\\\ 0\\\
-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

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

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

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

From (52) and (63) we get

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

From (65) we get that

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

The discrete reduced gradient now takes the form

(94)
$\frac{\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_{fi}(k).$

Under the assumption that the discrete steps are the same as in the forward
scheme and the quadrature coefficients satisfy
$\tilde{b}_{i}=b_{i},\tilde{c}_{i}=c_{i}$, the reduced gradients formula (39)
is

(95)
$\frac{d\phi(x(t_{f}))}{du_{n}}\simeq\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_{f}\left(t(k)+c_{i}h(k)\right).$

Now the quadrature scheme (95) for continuous reduced gradients is almost
identical to discrete reduced gradients (94). The only difference is that in
(94) $x_{i}(k+1)$ and $\lambda_{fi}(k)$ are used instead of
$x(t(k)+c_{i}h(k))$ and $\lambda_{f}(t(k)+c_{i}h(k))$. The convergence of (95)
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_{fi}(k)$ instead of $x(t(k)+c_{i}h(k))$ and
$\lambda_{f}(t(k)+c_{i}h(k))$ does not destroy the convergence of discrete
reduced gradients (94) to continuous reduced gradients.

Now we want to justify the convergence of discrete adjoint trajectories and
discrete reduced gradients to their continuous counterparts. To obtain the
right error order estimated we always assume that the system functions
$f(x,z,u)$ and $g(x)$ and their partial derivatives are Lipschitz continuous
functions. To achieve that goal we use the theorems presented in [3]. That
theorems are formulated under the constant step size assumption,

(96) $h(k)=h,\ k=1,\ldots,K,$

and we also carry out our analysis under that assumption. The variable step
size case can be cumbersome to analyze, and in this section we will make a
short note on that problem.

We assume that the system equations are integrated using RADAU IIA scheme for
$s=3$. The coefficients that appear in discrete adjoint equations
$\bar{a}_{ij}=\frac{a_{ji}b_{j}}{b_{i}},\bar{b}_{i}=b_{i},\bar{c}_{i}=1-c_{i}$
define the RADAU IA scheme (see [10]). The global error of the differential
state $x(t)$ is ([3], Theorem 4.4 p. 36 and Theorem 5.9 p. 67)

(97) $x(k)-x(t(k))=O(h^{p}).$

and the global error of the algebraic state $z(t)$ is ([3], Theorem 4.6 p. 40)

(98) $z(k)-z(t(k))=O(h^{q}).$

where for the RADAU IIA scheme and $s=3$ we have [3]

(99) $p=5,\ q=3.$

We emphasize that in the DAEs case the invertibility of the matrix of
coefficients $A=(a_{ij})$ is essential for convergence of the scheme. Also the
value of so called radius of stability

(100) $R(\infty)=1-b^{T}A^{-1}\mathbf{1}.$

plays an important role. For RADAU IA and RADAU IIA schemes $A$ matrices are
invertible and $R(\infty)=0$, so the assumptions of the appropriate theorems
from [3] are satisfied.

To derive the subsequent results we consider the unperturbed Runge – Kutta
scheme

(101) $\displaystyle x_{i}^{n}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h\sum_{j=1}^{s}a_{ij}f\left(x_{j}^{n}(k+1),\
z_{j}^{n}(k+1),u(k)\right),$ (102) $\displaystyle 0$ $\displaystyle=$
$\displaystyle g\left(x_{i}^{n}(k+1)\right),$

for $i=1,...,s$ and

(103) $\displaystyle x^{n}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h\sum_{i=1}^{s}b_{i}f\left(x_{i}^{n}(k+1),\
z_{i}^{n}(k+1),u(k)\right)$ (104) $\displaystyle z^{n}(k+1)$ $\displaystyle=$
$\displaystyle z^{n}(k)+\sum_{i=1}^{s}b_{i}^{-}(z_{i}^{n}(k+1)-z^{n}(k)).$

and the perturbed Runge-Kutta scheme

(105) $\displaystyle x_{i}^{p}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h\sum_{j=1}^{s}a_{ij}f\left(x_{j}^{p}(k+1),\
z_{j}^{p}(k+1),u(k)\right)+\delta_{i},$ (106) $\displaystyle 0$
$\displaystyle=$ $\displaystyle g\left(x_{i}^{p}(k+1)\right),$

for $i=1,...,s$ and

(107) $\displaystyle x^{p}(k+1)$ $\displaystyle=$ $\displaystyle
x^{n}(k)+h\sum_{i=1}^{s}b_{i}f\left(x_{i}^{p}(k+1),\
z_{i}^{p}(k+1),u(k)\right)+\delta_{s+1}.$ (108) $\displaystyle z^{p}(k+1)$
$\displaystyle=$ $\displaystyle
z^{n}(k)+\sum_{i=1}^{s}b_{i}^{-}(z_{i}^{p}(k+1)-z^{n}(k))+\delta_{s+2}.$

From ([3], Theorem 4.2 p. 33) we have

(109) $\displaystyle x_{i}^{p}(k+1)-x_{i}^{n}(k+1)$ $\displaystyle=$
$\displaystyle O(\delta)$ (110) $\displaystyle z_{i}^{p}(k+1)-z_{i}^{n}(k+1)$
$\displaystyle=$ $\displaystyle\frac{1}{h}O(\delta)$

where $\delta=max\\{\delta_{1},\dots,\delta_{s}\\}$. Contrary to the ODEs
case, using the variable step size for integration of DAEs may be cumbersome,
because of the factor $\frac{1}{h}$ that occurs in (110). The estimate (110)
indicates that using very small step sizes may lead to big errors. On the
other hand, the results for the constant step size can be extended to the
variable step size case if we assume that the following condition holds

(111) $h_{m}\geq h(k)\geq\mu h_{m},\ k=1,\ldots,K$

where $h_{m}$ is the maximum step size and $\mu\in(0,1)$ is a constant
independent from $h_{m}$.

The main step errors can be estimated as

(112) $\displaystyle x^{p}(k+1)-x^{n}(k+1)$ $\displaystyle=$ $\displaystyle
O(\delta)+O(\delta_{s+1})$ (113) $\displaystyle z^{p}(k+1)-z^{n}(k+1)$
$\displaystyle=$ $\displaystyle\frac{1}{h}O(\delta)+O(\delta_{s+2})$

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

(114) $\displaystyle x_{i}^{n}(k+1)$ $\displaystyle=$ $\displaystyle
x(t(k))+h\sum_{j=1}^{s}a_{ij}f\left(x_{j}^{n}(k+1),\
z_{j}^{n}(k+1),u(k)\right),$ (115) $\displaystyle 0$ $\displaystyle=$
$\displaystyle g\left(x_{i}^{n}(k+1)\right),$

for $i=1,...,s$. The following local error estimates are valid for RADAU IIA
scheme ([3], Lemma 4.3 p. 34)

(116) $\displaystyle x^{n}_{i}(k+1)-x(t(k)+c_{i}h)$ $\displaystyle=$
$\displaystyle O(h^{q+1})$ (117) $\displaystyle z^{n}_{i}(k+1)-z(t(k)+c_{i}h)$
$\displaystyle=$ $\displaystyle O(h^{q})$

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)=$ (118) $\displaystyle
x(k)+h\sum_{j=1}^{s}a_{ij}f\left(x_{j}(k+1),\ z_{j}(k+1),u(k)\right)=$
$\displaystyle x(t(k))+h\sum_{j=1}^{s}a_{ij}f\left(x_{j}(k+1),\
z_{j}(k+1),u(k)\right)+(x(k)-x(t(k)),$ (119) $\displaystyle
0=g\left(x_{i}(k+1)\right),$

for $i=1,...,s$, so the perturbations are

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

From (97) we have $\delta=x(k)-x(t(k))=O(h^{p})$. From (109) and (110) we
obtain

(121) $\displaystyle x_{i}(k+1)-x^{n}_{i}(k+1)$ $\displaystyle=$
$\displaystyle O(h^{p}).$ (122) $\displaystyle z_{i}(k+1)-z^{n}_{i}(k+1)$
$\displaystyle=$ $\displaystyle O(h^{p-1}).$

By combining (116)-(117) and (121)-(122) we obtain the required global error
estimates

(123) $\displaystyle x_{i}(k+1)-x(t(k)+c_{i}h)$ $\displaystyle=$
$\displaystyle O(h^{q+1})+O(h^{p})=O(h^{min\\{p,q+1\\}})$ $\displaystyle=$
$\displaystyle O(h^{q_{x}^{g}})$ (124) $\displaystyle
z_{i}(k+1)-z(t(k)+c_{i}h)$ $\displaystyle=$ $\displaystyle
O(h^{q})+O(h^{p-1})$ $\displaystyle=$ $\displaystyle O(h^{q_{z}^{g}})$

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

(125) $q_{x}^{g}=min\\{5,3+1\\}=4,\ q_{z}^{g}=min\\{5-1,3\\}=3.\\\ $

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)$

(126) $\displaystyle\lambda_{fi}^{n}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(t(k+1))+h\sum_{j=1}^{s}\bar{a}_{ij}\left[f_{x}^{T}\left(x\left(t_{j}(k)\right),z\left(t_{j}(k)\right),u(k)\right)\lambda_{fj}^{n}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x\left(t_{j}(k)\right)\right)\lambda_{gj}^{n}(k)\right],$
(127) $\displaystyle 0$ $\displaystyle=$
$\displaystyle-f^{T}_{z}\left(x\left(t_{i}(k)\right)\right)\lambda_{fi}^{n}(k)$

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

(128) $\displaystyle\lambda_{f}^{n}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(t(k+1))+h\sum_{i=1}^{s}\bar{b}_{i}\left[f_{x}^{T}\left(x\left(t_{i}(k)\right),z\left(t_{i}(k)\right),u(k)\right)\lambda_{fi}^{n}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x\left(t_{i}(k)\right)\right)\lambda_{gi}^{n}(k)\right],$
(129) $\displaystyle\lambda_{g}^{n}(k)$ $\displaystyle=$
$\displaystyle\lambda_{g}(t(k+1))+\sum_{i=1}^{s}\bar{b}_{i}^{-}(\lambda_{gi}^{n}(k)-\lambda_{g}(t(k+1))),$

The coefficients $\bar{a}_{ij},\bar{b}_{i},\bar{c}_{i}$ constitutes RADAU IA
scheme ([10]), so the following local error estimates are valid ([3], Theorem
4.3 p. 43)

(130) $\displaystyle\lambda^{n}_{fi}(k)-\lambda_{f}(t(k)+c_{i}h)$
$\displaystyle=$ $\displaystyle O(h^{\bar{q}+1}),\ i=1,\dots,s\ ,$ (131)
$\displaystyle\lambda^{n}_{gi}(k)-\lambda_{g}(t(k)+c_{i}h)$ $\displaystyle=$
$\displaystyle O(h^{\bar{q}}),\ i=1,\dots,s\ ,$ (132)
$\displaystyle\lambda^{n}_{f}(k)-\lambda_{f}(t(k))$ $\displaystyle=$
$\displaystyle O(h^{\bar{q}+1}),$ (133)
$\displaystyle\lambda^{n}_{g}(k)-\lambda_{g}(t(k))$ $\displaystyle=$
$\displaystyle O(h^{\bar{q}}).$

where for the RADAU IA scheme with $s=3$ we have ([3])

(134) $\bar{p}=5,\ \bar{q}=2.$

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

(135) $\displaystyle\lambda_{fi}^{p}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(t(k+1))+h\sum_{j=1}^{s}\bar{a}_{ij}\left[f_{x}^{T}\left(x_{j}(k+1),z_{j}(k+1),u(k)\right)\lambda_{fj}^{p}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x_{j}(k+1)\right)\lambda_{gj}^{p}(k)\right],$ (136)
$\displaystyle 0$ $\displaystyle=$
$\displaystyle-f^{T}_{z}\left(x_{i}(k+1)\right)\lambda_{fi}^{p}(k)$

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

(137) $\displaystyle\lambda_{f}^{p}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(t(k+1))+h\sum_{i=1}^{s}\bar{b}_{i}\left[f_{x}^{T}\left(x_{i}(k+1),z_{i}(k+1),u(k)\right)\lambda_{fi}^{p}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x_{i}(k+1)\right)\lambda_{gi}^{p}(k)\right],$ (138)
$\displaystyle\lambda_{g}^{p}(k)$ $\displaystyle=$
$\displaystyle\lambda_{g}(t(k+1))+\sum_{i=1}^{s}\bar{b}_{i}^{-}(\lambda_{gi}^{p}(k)-\lambda_{g}(t(k+1))),$

The perturbed system can be rewritten as

(139) $\displaystyle\lambda_{fi}^{p}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(t(k+1))+h\sum_{j=1}^{s}\bar{a}_{ij}\left[f_{x}^{T}\left(x(t_{j}(k)),z(t_{j}(k)),u(k)\right)\lambda_{fj}^{p}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x(t_{j}(k))\right)\lambda_{gj}^{p}(k)\right]+\delta_{i},$
(140) $\displaystyle 0$ $\displaystyle=$
$\displaystyle-f^{T}_{z}\left(x(t_{i}(k))\right)\lambda_{fi}^{p}(k)$

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

(141) $\displaystyle\lambda_{f}^{p}(k)$ $\displaystyle=$
$\displaystyle\lambda_{f}(t(k+1))+h\sum_{i=1}^{s}\bar{b}_{i}\left[f_{x}^{T}\left(x(t_{i}(k)),z(t_{i}(k)),u(k)\right)\lambda_{fi}^{p}(k)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\ \
\left.+g_{x}^{T}\left(x(t_{i}(k))\right)\lambda_{gi}^{p}(k)\right]+\delta_{s+1},$
(142) $\displaystyle\lambda_{g}^{p}(k)$ $\displaystyle=$
$\displaystyle\lambda_{g}(t(k+1))+\sum_{i=1}^{s}\bar{b}_{i}^{-}(\lambda_{gi}^{p}(k)-\lambda_{g}(t(k+1))),$

where the perturbations result from differences
$x_{i}(k+1)-x\left(t(k)+c_{i}h(k)\right)=O(h^{q_{x}^{g}})=O(h^{4})$ and
$z_{i}(k+1)-z\left(t(k)+c_{i}h(k)\right)=O(h^{q_{z}^{g}})=O(h^{3})$. The
perturbations satisfy

(143) $\delta_{i}=O(h^{q_{z}^{g}}),\ i=1,\dots,s+1.$

From (109)-(113) we obtain the error estimates

(144) $\displaystyle\lambda^{p}_{fi}(k)-\lambda^{n}_{fi}(k)$ $\displaystyle=$
$\displaystyle O(h^{q_{z}^{g}}),$ (145)
$\displaystyle\lambda^{p}_{gi}(k)-\lambda^{n}_{gi}(k)$ $\displaystyle=$
$\displaystyle O(h^{q_{z}^{g}-1}),$ (146)
$\displaystyle\lambda^{p}_{f}(k)-\lambda^{n}_{f}(k)$ $\displaystyle=$
$\displaystyle O(h^{q_{z}^{g}}),$ (147)
$\displaystyle\lambda^{p}_{g}(k)-\lambda^{n}_{g}(k)$ $\displaystyle=$
$\displaystyle O(h^{q_{z}^{g}-1}).$

By combining (130)-(133) with (144)-(147) we obtain the following local error
estimates

(148) $\displaystyle\lambda^{p}_{fi}(k)-\lambda_{f}(t(k)+c_{i}h)$
$\displaystyle=$ $\displaystyle
O(h^{min\\{q_{z}^{g},\bar{q}+1\\}})=O(h^{\bar{q}_{f}^{l}}),\ i=1,\dots,s\ ,$
(149) $\displaystyle\lambda^{p}_{gi}(k)-\lambda_{g}(t(k)+c_{i}h)$
$\displaystyle=$ $\displaystyle
O(h^{min\\{q_{z}^{g}-1,\bar{q}\\}})=O(h^{\bar{q}_{g}^{l}}),\ i=1,\dots,s\ ,$
(150) $\displaystyle\lambda^{p}_{f}(k)-\lambda_{f}(t(k))$ $\displaystyle=$
$\displaystyle O(h^{min\\{q_{z}^{g},\bar{q}+1\\}})=O(h^{\bar{p}_{f}^{l}}),$
(151) $\displaystyle\lambda^{p}_{g}(k)-\lambda_{g}(t(k))$ $\displaystyle=$
$\displaystyle O(h^{min\\{q_{z}^{g}-1,\bar{q}\\}})=O(h^{\bar{p}_{g}^{l}}).$

For the RADAU IA scheme with $s=3$ we have

(152) $\displaystyle\bar{q}_{f}^{l}=\bar{p}_{f}^{l}$ $\displaystyle=$
$\displaystyle min\\{q_{z}^{g},\bar{q}+1\\}=min\\{3,2+1\\}=3$ (153)
$\displaystyle\bar{q}_{g}^{l}=\bar{p}_{g}^{l}$ $\displaystyle=$ $\displaystyle
min\\{q_{z}^{g}-1,\bar{q}\\}=min\\{3-1,2\\}=2$

The global error of $\lambda^{p}_{f}(k)$ can be obtained from ([3], Theorem
4.4 p. 36)

(154)
$\lambda_{f}(k)-\lambda_{f}(t(k))=O(h^{\bar{p}_{f}^{l}-1})=O(h^{\bar{p}_{f}^{g}}).$

For RADAU IIA scheme with $s=3$ we obtain
$\bar{p}_{f}^{g}=min\\{q_{z}^{g},\bar{q}+1\\}-1=min\\{3,2+1\\}-1=2$. Let us
now derive the global error estimate of
$\lambda_{fi}(k)-\lambda_{f}(t(k)+c_{i}h(k))$. Using the global error estimate
(154) and the local error estimate (144) we can repeat the reasoning for
discrete system equations and derive the global error estimate

(155)
$\lambda_{fi}(k)-\lambda_{f}(t(k)+c_{i}h(k))=O(h^{min\\{\bar{p}_{f}^{g},\bar{q}_{f}^{l}\\}})=O(h^{\bar{q}_{f}^{g}}).$

For RADAU IIA scheme with $s=3$ we obtain $\bar{q}_{f}^{g}=$
$min\\{\bar{p}_{f}^{g},\bar{q}_{f}^{l}\\}$ $=min\\{2,3\\}=2$.

Having the global error estimates (123) for $x_{i}(k+1)-x(t_{i}(k))$ and (155)
for $\lambda_{fi}(k)-\lambda_{f}(t_{i}(k))$, we can repeat the derivation of
the reduced gradients calculation order presented in the first part of the
paper (for ODEs case) to obtain

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

where $\tilde{p}_{d}=min\\{\tilde{p},q_{x}^{g},\bar{q}_{f}^{g}\\}$. For RADAU
IIA scheme we obtain

(157) $\tilde{p}_{d}=min\\{5,4,2\\}=2.$

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

In Table 1 we summarize the proven minimal integration orders for system and
adjoint equations derived in this paper and in [10]. It follows that if a
state trajectory includes sections with sliding modes then the reduced
gradient can be approximated with the accuracy at least $O(h^{2})$, otherwise
the accuracy is at least $O(h^{3})$. This accuracy estimate takes into account
the fact that if a discrete state changes then we have to start the
integration of system and adjoint equations with the perturbed initial values
of system and adjoint variables respectively—the influence of that
perturbation on the accuracy of system and adjoint variables determination
follows from Theorem 4.3 in [4] (cf. (120)–(121) and (139) in [10]), or from
Theorem 4.2 in [3] (cf. (112)–(113) and (144)–(147)).

Table 1: Integration orders of convergence if RADAU IIA is applied to systems equations Equations | Orders $(x,z)$ | Orders $(\lambda_{f},\lambda_{g})$ | gradient orders
---|---|---|---
$x^{\prime}=f(x,u)$ | $p=5$ | $\bar{q}_{f}=4$ | $p_{d}=3$
$x^{\prime}=f(x,u)+g_{x}(x)z$ |  |  |
$0=g(x)$ |  $p=5$, $q=3$ |  $\bar{q}_{f}=2$, $\bar{q}_{g}=2$ | $p_{d}=2$

##  3 Calculating adjoint variables jumps

In this section we want to discuss the correspondence between jump conditions
for discrete and continuous adjoint equations at transition times. Let us
consider the hybrid system trajectory described by (1)-(4). As a result of the
numerical integration we obtain the discrete state equations (42).

In Section 1 we have stated that at a transition time $t_{t}$ adjoint
variables undergo jumps, the extent of the jump depends on the sequence of
discrete variables before and after the jump. The analysis which follows
concerns the case of the equations (1)-(4), in that case the extent of the
jump can be determined according to the equations (11)–(1)

These equations can be solved with respect to $\pi_{t}$ and
$\lambda_{f}(t_{t}^{-})$ giving:

$\displaystyle{\displaystyle\pi_{t}=-\frac{\lambda(t_{t}^{+})^{T}\left(f^{2}(x(t_{t}^{+}),z(t_{t}^{+}),u(t_{t}^{+}))-f^{1}(x(t_{t}^{-}),u(t_{t}^{-}))\right)}{g_{x}(x(t_{t}))f^{1}(x(t_{t}^{-}),u(t_{t}^{-}))}}$

(we have taken into account that $g(x(t_{t}^{+}))=0$ ).

Between transitions, a system of differential-algebraic equations is
integrated with the help of an appropriate numerical integration scheme. The
numerical integration scheme is represented by the equation

(158) $F(X(k+1),X(k),u(k),h(k))=0$

in which the discrete step $k$ corresponds to a time instant $t(k)$.

During the numerical integration of a hybrid trajectory, a possible violation
of invariant set conditions has to be monitored. This task is realized by
checking the sign changes of

$g(x(k))$

in subsequent steps, where $x(k)=x(t(k))$. When a sign change of $g(x(k))$
between discrete steps $k$ and $k+1$ is detected, the following problem is
solved

(159) ${\rm find}\ \ t_{t}\in[t(k),t(k+1)],\ {\rm s.t.}\ \ \hat{g}(t_{t})=0$

where $\hat{g}(\cdot)$ is a function, which approximates $g(\cdot))$ on a time
interval $[t(k),t(k+1)]$. When a transition time $t_{t}$ is found, the actual
iteration of numerical integration is repeated but with a step-size
$h(k)=t_{t}-t(k)$ instead of $h(k)=t(k+1)-t(k)$. The discrete step at which
the transition takes place we denote by $k_{t}$.

We assume that at discrete times $0,\ldots,k_{t}$ the system evolves according
to the equations:

(160) $F^{1}(X^{1}(k+1),X^{1}(k),u(k),h(k))=0,$

and at times $k_{t},\ldots,K-1$ by the equations

(161) $F^{2}(X^{2}(k+1),X^{2}(k),u(k),h(k))=0$

The optimal control problem with that system was investigated in [11].
Therein, the adjoint equations for the functional $\phi(x(t_{f}))$ were
established, herein parts of these equations are presented to expose jumps in
adjoint variables.

$\displaystyle{\rm for}\ k=N-1,\ldots,k_{t}+1$ (162a)
$\displaystyle\Lambda^{2}(k)$
$\displaystyle=-F^{2}_{X}(k_{t})^{T}\left[F^{2}_{X^{+}}(k_{t})\right]^{-T}\Lambda^{2}(k+1)$
(162b) $\displaystyle\Lambda^{2+}(k_{t})$
$\displaystyle=-F^{2}_{X}(k_{t})^{T}\left[F^{2}_{X^{+}}(k_{t})\right]^{-T}\Lambda^{2}(k_{t}+1)$
(162c)
$\displaystyle\tilde{\Lambda}^{1-}(k_{t})+\pi(k_{t})\left(g_{x}(k_{t})\right)^{T}$
$\displaystyle=\tilde{\Lambda}^{2+}(k_{t})$ (162d)
$\displaystyle\Lambda^{1-}(k_{t})^{T}\left[F^{1}_{X^{+}}(k_{t}-1)\right]^{-1}F^{1}_{h}(k_{t}-1)$
$\displaystyle=\Lambda^{2}(k_{t}+1)^{T}\left[F^{2}_{X^{+}}(k_{t})\right]^{-1}F^{2}_{h}(k_{t}),$
(162e) $\displaystyle\Lambda^{1}(k_{t}-1)$
$\displaystyle=-F^{1}_{X}(k_{t})^{T}\left[F^{1}_{X^{+}}(k_{t})\right]^{-T}\Lambda^{1-}(k_{t})$
$\displaystyle{\rm for}\ k=k_{t}-2,\ldots,1$ (162f)
$\displaystyle\Lambda^{1}(k)$
$\displaystyle=-F^{1}_{X}(k_{t})^{T}\left[F^{1}_{X^{+}}(k_{t})\right]^{-T}\Lambda^{1}(k+1).$

Here, $\tilde{\Lambda}^{1-}(k_{t})$ and $\tilde{\Lambda}^{2+}(k_{t})$ are
parts of vectors $\Lambda^{1}(k_{t})$ and $\Lambda^{2}(k_{t})$ respectively,
defined in such a way to be able to extract the essential part of Eq. (32d) in
[11]. Notice that

$\displaystyle\Lambda^{2}(k)$ $\displaystyle=$
$\displaystyle\left(\left(l_{f1}^{2}(k)\right)^{T},\left(l_{g1}^{2}(k)\right)^{T},\ldots,\left(l_{fs}^{2}(k)\right)^{T},\left(l_{gs}^{2}(k)\right)^{T},\left(\lambda_{f}^{2}(k)\right)^{T}\right)^{T},$

for $k=k_{t}+1,\ldots,K$, but

(163) $\displaystyle l_{fi}^{2}(k)$ $\displaystyle=$ $\displaystyle 0$ (164)
$\displaystyle l_{gi}^{2}(k)$ $\displaystyle=$ $\displaystyle 0$

for $i=1,\ldots,s$ and $k=k_{t}+1,\ldots,K$ (the justification of that is
given in Section 2), and similarly

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

for $k=0,\ldots,k_{t}-1$ and

(165) $\displaystyle l_{i}^{1}(k)$ $\displaystyle=$ $\displaystyle 0$

for $i=1,\ldots,s$ and steps $k=0,...,k_{t}$ (the justification for that is
given by Eq. (78) of Section 5 in [10]).

Therefore, we can take

$\displaystyle\tilde{\Lambda}^{1-}(k_{t})$ $\displaystyle=$
$\displaystyle\lambda^{1-}(k_{t})$ $\displaystyle\tilde{\Lambda}^{2+}(k_{t})$
$\displaystyle=$ $\displaystyle\lambda_{f}^{2+}(k_{t})$

and then (162c) becomes

(166)
$\displaystyle\lambda^{1-}(k_{t})+\pi(k_{t})\left(g_{x}(k_{t})\right)^{T}$
$\displaystyle=$ $\displaystyle\lambda_{f}^{2+}(k_{t})$

Our aim is to determine $\pi(k_{t})$ on the basis of equations (166) and
(162d). To this end we need analytical formula for $F^{1}_{X}(k_{t})$,
$F^{1}_{X^{+}}(k_{t})$, $F^{2}_{X}(k_{t})$, $F^{2}_{X^{+}}(k_{t})$,
$F^{1}_{h}(k_{t})$, $F^{2}_{h}(k_{t})$. The mappings $F^{1}$ and $F^{2}$ are
stated in Eq. (62) (in [10]) and (40) respectively. Furthermore, matrices
$F^{1}_{X+}$, $F^{1}_{X}$ are given by Eqns (70)–(71) in [10] and matrices
$F^{2}_{X+}$, $F^{2}_{X}$ by (43)–(44). It remains to provide formula for
$F^{1}_{h}(k_{t})$ and $F^{2}_{h}(k_{t})$. According to Eq. (62) (in [10]) and
(40) we have

(171) $\displaystyle F^{1}_{h}(k)=\frac{\partial
F^{1}(X^{1}(k+1),X^{1}(k),u(k),h(k))}{\partial h(k)}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}-\sum_{j=1}^{s}a_{1j}f^{1}(x_{j},u)\\\
\vdots\\\ -\sum_{j=1}^{s}a_{sj}f^{1}(x_{j},u)\\\
-\sum_{i=1}^{s}b_{i}f^{1}(x_{i},u)\end{array}\right)$ (178) $\displaystyle
F^{2}_{h}(k)=\frac{\partial F^{2}(X^{2}(k+1),X^{2}(k),u(k),h(k))}{\partial
h(k)}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}-\sum_{j=1}^{s}a_{1j}f^{2}(x_{j},z_{j},u)\\\
0\\\ \vdots\\\ -\sum_{j=1}^{s}a_{sj}f^{2}(x_{j},z_{j},u)\\\ 0\\\
-\sum_{i=1}^{s}b_{i}f^{2}(x_{i},z_{i},u)\end{array}\right).$

In order to derive analytical formula for $\pi(k_{t})$ we have to figure out
two vectors:
$w(k_{t})=\left[F^{2}_{X^{+}}(k_{t})\right]^{-1}F^{2}_{h}(k_{t})$,
$w(k_{t}-1)=\left[F^{1}_{X^{+}}(k_{t}-1)\right]^{-1}F^{1}_{h}(k_{t}-1)$ — both
these vectors are solutions to linear equations:

(191)
$\displaystyle{\displaystyle\left(\begin{array}[]{ccccccc}I-ha_{11}f_{x1}^{2}&-ha_{11}f_{y1}^{2}&\ldots&-ha_{1s}f_{xs}^{2}&-ha_{1s}f_{ys}^{2}&0\\\
-g_{x1}^{2}&0&\ldots&0&0&0\\\ \vdots&\vdots&&\vdots&\vdots&\vdots\\\
-ha_{s1}f_{x1}^{2}&-ha_{s1}f_{y1}^{2}&\ldots&I-ha_{ss}f_{xs}^{2}&-ha_{ss}f_{ys}^{2}&0\\\
0&0&\ldots&-g_{xs}^{2}&0&0\\\
-hb_{1}f_{x1}^{2}&-hb_{1}f_{y1}^{2}&\ldots&-hb_{s}f_{xs}^{2}&-hb_{s}f_{ys}^{2}&I\end{array}\right)\left(\begin{array}[]{c}w_{f}^{1}\\\
w_{g}^{1}\\\ \vdots\\\ w_{f}^{s}\\\ w_{g}^{s}\\\ w_{f}\end{array}\right)=}$
(198)
$\displaystyle{\displaystyle\left(\begin{array}[]{c}-\sum_{j=1}^{s}a_{1j}f^{2}(x_{j},z_{j},u)\\\
0\\\ \vdots\\\ -\sum_{j=1}^{s}a_{sj}f^{2}(x_{j},z_{j},u)\\\ 0\\\
-\sum_{i=1}^{s}b_{i}f^{2}(x_{i},z_{i},u)\end{array}\right),}$
$\left(\begin{array}[]{cccc}I-ha_{11}f_{x1}^{1}&\ldots&-ha_{1s}f_{xs}^{1}&0\\\
\vdots&&\vdots&\vdots\\\ -ha_{s1}f_{x1}^{1}&\ldots&I-ha_{ss}^{1}f_{xs}&0\\\
-hb_{1}f_{x1}^{1}&\ldots&-hb_{s}f_{xs}^{1}&I\end{array}\right)\left(\begin{array}[]{c}w^{1}\\\
\vdots\\\ w^{s}\\\
w\end{array}\right)=\left(\begin{array}[]{c}-\sum_{j=1}^{s}a_{1j}f^{1}(x_{j},u)\\\
\vdots\\\ -\sum_{j=1}^{s}a_{sj}f^{1}(x_{j},u)\\\
-\sum_{i=1}^{s}b_{i}f^{1}(x_{i},u)\end{array}\right).$

Eventually, to transform equation (162d) to its useful form we take into
account (163)–(164) and (165), then the equation (162d) becomes

$\displaystyle{\displaystyle\lambda^{1-}(k_{t})^{T}\left(h(k_{t}-1)\sum_{i=1}^{s}b_{i}f^{1}_{xi}(k_{t}-1)w^{i}(k_{t}-1)\right.}$
$\displaystyle{\displaystyle\left.-\sum_{i=1}^{s}b_{i}f^{1}(x_{i}(k_{t}),u(k_{t}-1))\right)}$
$\displaystyle{\displaystyle=\lambda_{f}^{2}(k_{t}+1)^{T}\left(h(k_{t})\sum_{i=1}^{s}b_{i}f^{2}_{xi}(k_{t})w^{i}_{f}(k_{t})+h(k_{t})\sum_{i=1}^{s}b_{i}f^{2}_{yi}(k_{t})w^{i}_{g}(k_{t})\right.}$
(199)
$\displaystyle{\displaystyle\left.-\sum_{i=1}^{s}b_{i}f^{2}(x_{i}(k_{t}+1),z_{i}(k_{t}+1),u(k_{t}))\right)}$

By plugging $\lambda^{1-}(t_{k})$ from equations (199) into equation (166) one
will get

$\displaystyle{\displaystyle\lambda_{f}^{2+}(k_{t})^{T}\left(h(k_{t}-1)\sum_{i=1}^{s}b_{i}f^{1}_{xi}(k_{t}-1)w^{i}(k_{t}-1)\right.}$
$\displaystyle{\displaystyle\left.-\sum_{i=1}^{s}b_{i}f^{1}(x_{i}(k_{t}),u(k_{t}-1))\right)}$
$\displaystyle{\displaystyle-\pi
g_{x}(k_{t})\left(h(k_{t}-1)\sum_{i=1}^{s}b_{i}f^{1}_{xi}(k_{t}-1)w^{i}(k_{t}-1)\right.}$
$\displaystyle{\displaystyle\left.-\sum_{i=1}^{s}b_{i}f^{1}(x_{i}(k_{t}),u(k_{t}-1))\right)}$
$\displaystyle{\displaystyle=\lambda^{2}_{f}(k_{t}+1)^{T}\left(h(k_{t})\sum_{i=1}^{s}b_{i}f^{2}_{xi}(k_{t})w^{i}_{f}(k_{t})+h(k_{t})\sum_{i=1}^{s}b_{i}f^{2}_{yi}(k_{t})w^{i}_{g}(k_{t})-\right.}$
(200)
$\displaystyle{\displaystyle\left.\sum_{i=1}^{s}b_{i}f^{2}(x_{i}(k_{t}+1),z_{i}(k_{t}),u(k_{t}))\right)}$

and eventually

$\displaystyle{\displaystyle\pi(k_{t})=\left[\lambda^{2}_{f}(k_{t}+1)^{T}\left(-h(k_{t})\sum_{i=1}^{s}b_{i}f^{2}_{xi}(k_{t})w^{i}_{f}(k_{t})-h(k_{t})\sum_{i=1}^{s}b_{i}f^{2}_{yi}(k_{t})w^{i}_{g}(k_{t})\right.\right.}$
$\displaystyle{\displaystyle\left.+\sum_{i=1}^{s}b_{i}f^{2}(x_{i}(k_{t}+1),z_{i}(k_{t}),u(k_{t}))\right)-\lambda_{f}^{2+}(k_{t})^{T}\times}$
$\displaystyle{\displaystyle\left.\left(-h(k_{t}-1)\sum_{i=1}^{s}b_{i}f^{1}_{xi}(k_{t}-1)w_{i}(k_{t}-1)+\sum_{i=1}^{s}b_{i}f^{1}(x_{i}(k_{t}),u(k_{t}-1))\right)\right]\bigg{/}}$
$\displaystyle{\displaystyle
g_{x}(k_{t})\left(h(k_{t}-1)\sum_{i=1}^{s}b_{i}f^{1}_{xi}(k_{t}-1)w_{i}(k_{t}-1)-\sum_{i=1}^{s}b_{i}f^{1}(x_{i}(k_{t}),u(k_{t}-1))\right).}$
(201)

We assume that all functions values and their partial derivatives: $f^{1}$,
$f^{2}$, $f^{1}_{xi}$, $f^{2}_{xi}$, $f^{2}_{yi}$ are uniformly bounded in a
neighborhood of the considered solution (cf. assumptions (H1) and (H2) in
[9]). Furthermore, in the neighborhood, matrices $F^{1}_{X+}$, $F^{2}_{X+}$
are invertible and their elements uniformly bounded. Therefore, there exists
$0<L<+\infty$ such that $\|z\|\leq L$ and $\|w\|\leq L$.

We consider RADAU IIA scheme which has the property $\sum_{i=1}^{s}b_{i}=1$.
That together with the boundedness of $z$ and $w$ show that

(202) $\displaystyle{\displaystyle\pi(t_{k})\rightarrow\pi_{t}}$

provided that $h(k)\rightarrow 0$.

However that property can be not sufficient to guarantee high accuracy of
solutions to a continuous optimal control problem with hybrid system when a
Runge–Kutta method with variable step sizes is applied. For that method we
prefer a scheme for adjoint variables jumps which is determined by the formula
$\pi(k_{t})$ with the property

(203) $\displaystyle{\displaystyle|\pi(k_{t})-\pi_{t}|\leq
O(h^{\bar{q}_{\pi}}),}$

where the order of convergence $\bar{q}$ can be estimated on the basis of the
integration orders for state and adjoint variables.

That can be achieved for many new formula for $\pi(k_{t})$ which does not
dispatch far away from the formula (201). For example, consider the formula

$\displaystyle{\displaystyle\hat{\pi}(k_{t})=\frac{\lambda_{f}^{2+}(k_{t})^{T}\left(f^{2}(x(k_{t}),z(k_{t}),u(k_{t}-1))-f^{1}(x(k_{t}),u(k_{t}-1))\right)}{g_{x}(k_{t})f^{1}(x(k_{t}),u(k_{t}-1))}.}$
(204)

Then, we will have

$\displaystyle{\displaystyle\hat{\pi}(t_{k})-\pi_{t}=}$
$\displaystyle{\displaystyle
g_{x}(x(t_{t}))f^{1}(x(t_{t}^{-}),u(t_{t}^{-}))\lambda_{f}^{2+}(k_{t})^{T}\left(f^{2}(x(k_{t}),z(k_{t}),u(k_{t}-1))-\right.}$
$\displaystyle{\displaystyle\left.f^{1}(x(k_{t}),u(k_{t}-1))\right)\bigg{/}g_{x}(k_{t})f^{1}(x(k_{t}),u(k_{t}-1)){g_{x}(x(t_{t}))f^{1}(x(t_{t}^{-}),u(t_{t}^{-}))}-}$
$\displaystyle{\displaystyle
g_{x}(k_{t})f^{1}(x(k_{t}),u(k_{t}-1))\lambda(t_{t}^{+})^{T}\left(f^{2}(x(t_{t}^{+}),z(t_{t}^{+}),u(t_{t}^{+}))-\right.}$
$\displaystyle{\displaystyle\left.f^{1}(x(t_{t}^{-}),u(t_{t}^{-}))\right)\bigg{/}g_{x}(k_{t})f^{1}(x(k_{t}),u(k_{t}-1))g_{x}(x(t_{t}))f^{1}(x(t_{t}^{-}),u(t_{t}^{-})).}$

Since we have to assume that in the neighborhood of $x(t_{t})$ the following
regularity assumption must hold (cf. the assumption (H3) in [9])

(205) $\displaystyle{\displaystyle\left|g_{x}(x)f^{1}(x,u(t_{t}))\right|\geq
M>0,}$

and the product of Lipschitz functions is a Lipschitz function on the bounded
domain, after some transformations we will arrive at the relation

$\displaystyle\left|\hat{\pi}(t_{k})-\pi_{t}\right|$ $\displaystyle\leq$
$\displaystyle
L\left(\|g_{x}(x(k_{t}))-g_{x}(x(t_{t}))\|+\|\lambda_{f}^{2+}(k_{t})-\lambda(t_{t}^{+})\|+\right.$
$\displaystyle\|f^{1}(x(k_{t}),u(k_{t}))-f^{1}(x(t_{t}),u(t_{t}))\|+$
$\displaystyle\left.\|f^{2}(x(k_{t}),z(k_{t}),u(k_{t}))-f^{2}(x(t_{t}),z(t_{t}),u(t_{t}))\|\right)$

Since

$\displaystyle\lambda_{f}^{2+}(k_{t})-\lambda(t_{t}^{+})$ $\displaystyle=$
$\displaystyle O(h^{\bar{q}_{f}})$ $\displaystyle x(k_{t})-x(t_{t})$
$\displaystyle=$ $\displaystyle O(h^{p})$ $\displaystyle z(k_{t})-z(t_{t})$
$\displaystyle=$ $\displaystyle O(h^{q})$

($\bar{q}_{f}$, $p$ and $q$ are the appropriate integration orders) and if we
assume that functions $g_{x}$, $f^{1}$, $f^{2}$ are Lipschitz continuous then
(203) holds with $\bar{q}_{\pi}=min\\{\bar{q}_{f},p,q\\}$. From (166) the
adjoint variable $\lambda^{1-}(k_{t})$ is calculated with the error
$O(h^{\bar{q}_{f}})$. $\lambda^{1-}(k_{t})$ states an initial condition for
the calculation of the adjoint variables in the non-sliding phase. The order
of the global error of the adjoint variables in the non-sliding phase is
therefore reduced relative to the estimates presented in Table 1. The order of
convergence for reduced gradients is therefore equal to $\tilde{p}_{d}=2$ for
both sliding and non sliding phase.

##  4 Numerical results

In this section we present results of solving three optimal control problems
with hybrid systems by using the methods discussed in this paper and in the
papers accompanying it ([9],[10]). The results have been obtained by our
preliminary software based on its two core subroutines: the first one aimed at
solving optimal control problems with piecewise constant approximations to
control functions ([7]); the second one for evaluating solutions to
differential–algebraic equations and their corresponding adjoint equations.

To solve the reported problems we applied the SQP code described in [7],
however instead of using range–space active set method for solving QP
subproblems the new version of the SQP code is based on the implementation of
the interior point method as described in [2]. The matrices $H_{k}$ used in
the direction finding subproblems ${\bf P_{c}(u)}$ ([10]) were evaluated
according to BFGS updates with the Powell’s modifications ([6])—see the
discussion in Section 4 of Chapter 5 of [7].

We have used RADAU5 subroutine, which is the Fortran implementation of the
RADAU IIA scheme. The subroutine does not have the facility for locating
switching points and then restarting the procedure with a new description of
differential–algebraic equations. In order to enhance RADAU5 procedure
applicability to DAEs with hybrid description we incorporated into it the
subroutine ROOTS (for finding a root of nonlinear algebraic equations with a
secant method) from the SUNDIALS package. Since the ROOTS subroutine was
implemented with multistep integration methods in mind (in which polynomial
approximations of state variables $(x,y)$ between mesh points are provided by
methods themselves) we had to work out an approach to the interpolation of
state variables between mesh points which would be suitable for Runge–Kutta
methods.

In our implementation switching points $t_{t}$ are evaluated with the help of
$x_{i}(k)$, $i=1,\ldots,s$ determined at the intermediate points $t(k)$,
$t(k)+c_{i}h(k)$, $i=1,\ldots,s$. Having $c_{i}$ and $x_{i}(k)$ we evaluate
derivatives at these points:

$\displaystyle{\displaystyle
x_{i}^{{}^{\prime}}(k)=f^{1}(x_{i}(k),u(k-1),t(k)+c_{i}h(k)),\ i=1,\ldots,s}$

and then construct the Hermite polynomials by taking into accounts points
$x_{i}(k)$ and their derivatives $x_{i}^{{}^{\prime}}(k)$, $i=1,\ldots,s$. In
our current implementation we use two vectors (and their corresponding
derivatives) to build the polynomial approximation to vector $x$ on the time
interval $[t(k),t(k)+h(k)]$ ([1]):

(206) $\displaystyle x(t)$ $\displaystyle=$
$\displaystyle(2\tau^{3}-3\tau^{2}+1)x(k)+h(k)(\tau^{3}-2\tau^{2}+\tau)x^{{}^{\prime}}(k)+$
$\displaystyle(-2\tau^{3}+3\tau^{2})x_{s}(k)+h(k)(\tau^{3}-\tau^{2})x_{s}^{{}^{\prime}}(k),$

That interpolating scheme guarantees the accuracy $O(h(k)^{4})$ ([1]). We
notice, that the order of the interpolating scheme is lower than the numerical
integration order $p=5$ and affects the integration order after the discrete
transition. Nevertheless, the estimated orders for reduced gradients
calculation presented in Table 1 remains valid if we assume integration order
$p=4$.

The use of interpolating polynomials in locating switching times can influence
also the accuracy with which these times are determined. Suppose that on an
open interval $(a,b)$, such that $t_{t}\in(a,b)$, the state trajectory $x(t)$
is perturbed by $\delta x(t)$. To that state perturbation corresponds the
perturbation of the switching time, $t_{t}+\delta t_{t}$. We will have

(207) $\displaystyle{\displaystyle g(x(t_{t}+\delta t_{t})+\delta
x(t_{t}+\delta t_{t}))=0.}$

Evaluating $g$ around $(t_{t},0)$ will result in

(208) $\displaystyle 0$ $\displaystyle=$ $\displaystyle g(x(t_{t}+\delta
t_{t})+\delta x(t_{t}+\delta t_{t}))=g_{x}(x(t_{t}))\delta x(t_{t})+$
$\displaystyle g_{x}(x(t_{t}))f^{1}(x(t_{t}),u(t_{t}))\delta t_{t}+o(\delta
t_{t},\delta x).$

Eventually, we have

(209) $\displaystyle{\displaystyle\delta
t_{t}=-\frac{1}{g_{x}(x(t_{t}))f^{1}(x(t_{t}),u(t_{t}))}g_{x}(x(t_{t}))\delta
x(t_{t})+o(\delta t_{t},\delta x).}$

where $\tau=(t-t(k))/h(k)$.

This together with (205) imply that to the perturbation of $x$ by $O(h^{p})$
corresponds the perturbation of the switching time $t_{t}$ of the same order
$O(h^{p})$. It is straightforward to show that this perturbation will have the
same effect on order of convergence of integrating procedures as the
perturbation of initial conditions of the same order.

The next issue concerning the implementation of adjoint equations evaluation
scheme is that related to the jump formula. To obtain the numerical results
reported in the paper we have applied the formula (201). Since it guarantees
convergence of the adjoint equations Runge–Kutta scheme to their continuous
counterpart under the assumption $h(k)\rightarrow 0$ we set absolute and
relative tolerances to $10^{-9}$ in RADAU5 procedure. Consequently our
tolerances for satisfying optimality conditions had to be essentially less
demanding. Our optimization procedure stopped iterating when the following set
of conditions was satisfied:

$\displaystyle\sigma_{c_{k}}^{H_{k}}(u_{k})$ $\displaystyle\geq$
$\displaystyle-\varepsilon,$ $\displaystyle\left|g^{1}_{i}(u_{k})\right|$
$\displaystyle\leq$ $\displaystyle\varepsilon,\ i\in E,$ $\displaystyle
g^{2}_{j}(u_{k})$ $\displaystyle\leq$ $\displaystyle\varepsilon,\ j\in I.$

We set $\varepsilon=10^{-6}$ in our calculations.

We show the results of application of our numerical procedures to three
optimal control problems with nonlinear differential equations. In the first
two examples functions $g(x)$ defining switching surface are linear, in the
third example it is nonlinear.

Example 1 The first example concerns the Coulomb–Stribeck friction model. A
mass $m$ is attached to inertial space with a spring $k$. The mass is riding
on a belt, which itself is moving with a constant velocity $v_{dr}$ (see
Figure 6.4a in [5]).The relative velocity of the mass with respect to the belt
is equal to $v_{rel}=v-v_{dr}$. Between the mass and the belt there is the dry
friction with a friction force $F_{T}$. In the slip phase it is the function
of $v_{rel}$ and is given by the relation

(210) $\displaystyle{\displaystyle
F_{T}=-\frac{\mu_{s}}{1+\delta|v_{rel}|}F_{N}{\rm sign}(v_{rel}).}$

Here, $F_{N}=mg$. Furthermore, in the stick phase the friction force is
limited by the relation $|F_{T}|\leq F_{s}=\mu_{s}mg$. The functions $f_{1}$,
$f_{2}$ and $g$ of the hybrid system are

(213) $\displaystyle f^{1}(x,u)=\left[\begin{array}[]{c}x_{2}\\\
-\frac{k}{m}x_{1}+\frac{1}{m}\frac{F_{s}}{1+\delta|x_{2}-v_{dr}|}+x_{3}\\\
u\end{array}\right],$ (216) $\displaystyle
f^{2}(x,u)=\left[\begin{array}[]{c}x_{2}\\\
-\frac{k}{m}x_{1}-\frac{1}{m}\frac{F_{s}}{1+\delta|x_{2}-v_{dr}|}+x_{3}\\\
u\end{array}\right]$

and

(217) $\displaystyle{\displaystyle g(x)=x_{2}-v_{dr}.}$

Here, $x_{1}$ corresponds to the mass position, $x_{2}$ to its velocity and
$x_{3}$ influences the mass movement through the control $u$.

The optimal control problem is as follows

(218) $\displaystyle\min_{u\in{\mathcal{U}}}x_{2}(t_{f})$ (219)
$\displaystyle{\rm s.\ t.}$ (220) $\displaystyle\ x(t_{0})=x_{0}$ (221)
$\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle f^{1}(x,u),\
\text{if }g(x)<0$ (222) $\displaystyle x^{\prime}$ $\displaystyle=$
$\displaystyle f^{2}(x,u),\ \text{if }g(x)>0$ (223) $\displaystyle x^{\prime}$
$\displaystyle=$ $\displaystyle f_{F}(x,u)+g_{x}^{T}(x)z,\ \text{if sliding
mode occurs}$ $\displaystyle 0$ $\displaystyle=$ $\displaystyle g(x)$

and

(224) $\displaystyle{\displaystyle x_{1}(t_{f})-0.6=0,}$

with the additional constraint on the control signal value $u$

$\displaystyle{\displaystyle-2.5\leq u(t)\leq 2.5\ t\in[t_{0},t_{f}],}$

and the time interval endpoints $t_{0}=0$ $t_{f}=1$.

The optimal trajectories are shown in Fig. 1(a) and the optimal control is
shown in Fig. 1(b). We used $N=100$ in piecewise constant approximations of
control functions. One can observe that from the time $t\approx 0.8$ the
system is in the sliding mode. The program needed $4$ iterations to find the
solution.

(a) Optimal state trajectories.

(b) Optimal control.

Figure 1: Mass–spring example optimal solution

Example 2 The example concerns the application of our approach to the problem
of planning a haemodialysis process. The problem is fully described in [13].
In this paper we present some results obtained for a variant of the problem.

The system equations that determine the concentrations of urea and phosphorus
in intracellular fluid – $C_{IC}^{urea}$, $C_{IC}^{PO_{4}}$, urea and
phosphorus concentrations in extracellular fluid – $C_{EC}^{urea}$,
$C_{EC}^{PO_{4}}$ and ultrafiltration volume – $UFR$ are as follows

$\displaystyle\frac{dC_{EC}^{urea}}{dt}$ $\displaystyle=$
$\displaystyle\frac{K_{IE}^{urea}\cdot(C_{IC}^{urea}-C_{EC}^{urea})-C_{EC}^{urea}\cdot(K_{D}^{urea}+K_{r}^{urea}+K^{ufr})}{0.34\cdot
V(0)-UFR}$ (226) $\displaystyle\frac{dC_{IC}^{urea}}{dt}$ $\displaystyle=$
$\displaystyle\frac{K_{IE}^{urea}\cdot(C_{EC}^{urea}-C_{IC}^{urea})+G^{urea}}{0.66\cdot
V(0)}$ (227) $\displaystyle\frac{dC_{EC}^{PO_{4}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{K_{IE}^{PO_{4}}\cdot(C_{IC}^{PO_{4}}-C_{EC}^{PO_{4}})-K_{D}^{PO_{4}}\cdot
C_{EC}^{PO_{4}}}{0.34\cdot V(0)-UFR}+K_{3}^{PO_{4}}+K_{4}^{PO_{4}}$ (228)
$\displaystyle\frac{dC_{IC}^{PO_{4}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{K_{IE}^{PO_{4}}\cdot(C_{EC}^{PO_{4}}-C_{IC}^{PO_{4}})}{0.66\cdot
V(0)}$ (229) $\displaystyle\frac{dUFR}{dt}$ $\displaystyle=$ $\displaystyle
U_{ufr}$ (230) $\displaystyle K_{4}^{PO_{4}}$ $\displaystyle=$
$\displaystyle\alpha\cdot\max\left(C_{min}^{PO_{4}}-C_{IC}^{PO_{4}},0\right)$
(231) $\displaystyle K_{3}^{PO_{4}}$ $\displaystyle=$
$\displaystyle\beta\cdot\max\left(C_{max}^{PO_{4}}-C_{IC}^{PO_{4}},0\right)$

where

(232) $\displaystyle{\displaystyle
K_{D}^{urea}=\frac{e^{K_{0}A\left(Q_{D}-Q_{B}\right)/\left(Q_{B}Q_{D}\right)}-1}{e^{K_{0}A\left(Q_{D}-Q_{B}\right)/\left(Q_{B}Q_{D}\right)}-Q_{D}/Q_{B}}}.$

The model coefficients are explained in [13]—we used the model parameters as
stated for the first run of optimization. The algebraic equations are
responsible for hybrid behavior of the system equations. We have
$0<C_{min}^{PO_{4}}-C_{IC}^{PO_{4}}<C_{max}^{PO_{4}}<+\infty$ and two
switching functions:

(233) $\displaystyle g_{1}(C_{IC}^{PO_{4}})$ $\displaystyle=$ $\displaystyle
C_{IC}^{PO_{4}}-C_{min}^{PO_{4}}$ (234) $\displaystyle g_{1}(C_{IC}^{PO_{4}})$
$\displaystyle=$ $\displaystyle C_{IC}^{PO_{4}}-C_{max}^{PO_{4}}.$

Having combined kinetic models of urea and phosphorus we look for proper
concentrations of urea and phosphorus at the end of the haemodialysis process
by controlling the parameters $Q_{B}$, $Q_{D}$ and $U_{ufr}$ (see [13] for
details). In other words, by solving the optimal control problem we want to
choose a proper dialysis membrane in order to achieve final parameters of
haemodialysis.

The optimization problem is as follows:

(235) $\min_{Q_{B},Q_{D},U_{ufr}}C_{EC}^{urea}(t_{f})$

subject to the constraints (LABEL:dialysis1)–(231), the following constraints
at final time

(236) $\displaystyle C_{IC}^{urea}(t_{f})$ $\displaystyle\leq L^{urea}_{IC}$
(237) $\displaystyle L_{min}^{UFR}\leq$ $\displaystyle UFR(t_{f})$
$\displaystyle\leq L_{max}^{UFR},$

and the constraints on the control variables

(238) $\displaystyle Q_{B}^{min}\leq$ $\displaystyle Q_{B}(t)$
$\displaystyle\leq Q_{B}^{max}$ (239) $\displaystyle Q_{D}^{min}\leq$
$\displaystyle Q_{D}(t)$ $\displaystyle\leq Q_{D}^{max}$ (240) $\displaystyle
U_{ufr}^{min}\leq$ $\displaystyle U_{ufr}(t)$ $\displaystyle\leq
U_{ufr}^{max}.$

on $[t_{0},t_{f}]$. Parameters for the inequality constraints are the same as
in [13] with the exception that now we allow $U_{ufr}$ to vary by assuming:
$U_{ufr}^{min}=0$, $U_{ufr}^{max}=75$.

Some of the optimal trajectories are shown in Figs. 2–3. One can observe that
from the time $t\approx 125$ the system is in the sliding mode. Optimal
control variable $U_{ufr}$ is shown in Fig. 4, the other control variables
reached their allowed maximum values—we used $N=10$ in the piecewise constant
approximations of control variables. Optimization procedure needed $69$
iterations to find the solution with the specified accuracy.

Figure 2: State trajectories for the urea $C_{IC}^{urea}$ and $C_{EC}^{urea}$.
Figure 3: State trajectories for the phosphorus $C_{IC}^{PO_{4}}$ and
$C_{EC}^{PO_{4}}$ (rebound of $C_{EC}^{PO_{4}}$ can be seen in 125 minute of
haemodialysis).  Figure 4: Control trajectory for the $U_{ufr}$.

Example 3

The example concerns the optimal control of a of a race car discussed in [14]
and depicted in Fig.5.

$n(\theta)$$t(\theta)$$v$$Fn(\theta)$$x_{2}$$x_{1}$$\theta$ Figure 5: Race car
model – mathematical description

We introduce the following notation

  * •

$x=(x_{1},x_{2})^{T}\in\mathbb{R}^{2}$ is a position of a car in the plane

  * •

$v=(v_{1},v_{2})^{T}\in\mathbb{R}^{2}$ is a velocity of a car

  * •

$\theta\in\mathbb{R}$ is an angle of the car orientation

  * •

$t(\theta)=(\cos\theta,\sin\theta)^{T}$ is an unit vector pointing in the
direction $\theta$

  * •

$n(\theta)=(-\sin\theta,\cos\theta)^{T}$ is an unit vector normal to
$t(\theta)$

  * •

$a(t)\in\mathbb{R}$ is an acceleration force

  * •

$s(t)\in\mathbb{R}$ is a steering control

The equations of motion are [14]

(241) $\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle v,$ (242)
$\displaystyle v^{\prime}$ $\displaystyle=$ $\displaystyle
a(t)t(\theta)+Fn(\theta),$ (243) $\displaystyle\theta^{\prime}$
$\displaystyle=$ $\displaystyle s(t)(t(\theta)^{T}v),$

where $Fn(\theta)$ is a friction force vector. $F$ is given by

(244) $F=-\mu N{\rm sign}(n(\theta)^{T}v),$

where $\mu$ and $N$ are a friction coefficient and a normal contact force
respectively. The switching surface is therefore defined as a solution set for
an equation

(245) $n(\theta)^{T}v=-v_{1}\sin\theta+v_{2}\cos\theta=0.$

During the non-sliding motion the amplitude of a friction force vector is
constant and equal to $\mu N$. The friction force is normal to $t(\theta)$ and
it acts such that it decreases the component of a velocity which is normal to
a direction $t(\theta)$. In other words the friction force always attempts to
reduce the angle between $t(\theta)$ and $v$ to zero. Physically, the non-
sliding motion corresponds to a drift of a car.

We face the the sliding motion when the following conditions are satisfied

(246) $\displaystyle n(\theta)^{T}v$ $\displaystyle=$ $\displaystyle 0,$ (247)
$\displaystyle-\mu N-s(t)(t(\theta)^{T}v)^{2}$ $\displaystyle<$ $\displaystyle
0,$ (248) $\displaystyle\mu N-s(t)(t(\theta)^{T}v)^{2}$ $\displaystyle>$
$\displaystyle 0.$

During the sliding motion the car is always oriented in the direction of a car
velocity ($t(\theta)$ and $v$ are collinear), there is no drift and the
friction force prevents the motion in a direction normal to $t(\theta)$.

We solved the following optimization problem:

(249) $\displaystyle{\displaystyle\min_{a,s}x_{1}(t_{f})}$

subject to the system equations (241)–(243), the endpoint constraints

(250) $\displaystyle x_{2}(t_{f})$ $\displaystyle=$ $\displaystyle 0$ (251)
$\displaystyle-v_{1}(t_{f})$ $\displaystyle\leq$ $\displaystyle 0$ (252)
$\displaystyle v_{2}(t_{f})$ $\displaystyle=$ $\displaystyle 0$ (253)
$\displaystyle\theta(t_{f})$ $\displaystyle=$ $\displaystyle 0,$

and the constraints on controls

(254) $\displaystyle-0.3\leq$ $\displaystyle a(t)$ $\displaystyle\leq 0.3$
(255) $\displaystyle-1.0\leq$ $\displaystyle s(t)$ $\displaystyle\leq 1.0$

on an interval $[t_{0},t_{f}]$ ($t_{f}=3.0$). We assumed $\mu N=0.5$ and
initial values: $x_{1}(0)=0$, $x_{2}(0)=1$, $v_{1}(0)=1$, $v_{2}(0)=0$.

Our program needed 52 iterations to find an approximation to an optimal
solution with the required accuracy, we assumed $N=10$ in the piecewise
constant approximations of control functions. The optimal controls, positions
and velocities are presented in Fig.6, 7(a) and 7(b) respectively.

Figure 6: Optimal controls for the racing car problem.

(a) Optimal trajectories of positions.

(b) Optimal trajectories of velocities.

Figure 7: Racing car – optimal solution

The optimal orientation trajectory as well as the value of $n^{T}(\theta)v$
along the trajectory are presented in Fig.8(a).

Near the end of the time interval the condition $n^{T}(\theta)v=0$ becomes
satisfied, so the system enters the sliding mode. The relevant fragment of
system trajectories is presented in Fig.8(b).

(a) $\theta$ and $n^{T}(\theta)v$ trajectories.

(b) $\theta$ and $n^{T}(\theta)v$ near the transition to the sliding mode.

Figure 8: Racing car – optimal solution

##  5 Conclusions

The paper presents the computational approach to hybrid optimal control
problems with sliding modes. It seems to be the first method for optimal
control problems with hybrid systems which can exhibit sliding modes. Our
computational method is based on a Runge–Kutta method for integrating system
and adjoint equations. The possible improvements of our code are related to
the approximating rules for the jumps of adjoint variables (such as (204)) and
polynomial approximations of state variables used in the procedure of the
switching times $t_{t}$ localization (cf. (206)). These issues are the subject
of our current research.

## References

  * [1] K.E. Atkinson, An Introduction to Numerical Analysis, J. Wiley & Sons, (1974).
  * [2] E. M. Gertz and S. J. Wright, ”Object–Oriented Software for Quadratic Programming”, ACM Transactions on Mathematical Software, vol. 29, pp. 58–81 (2003).
  * [3] 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).
  * [4] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Springer–Verlag, Berlin, Heidelberg, (1996).
  * [5] R. Leine and H. Nimeijer, Dynamics and Bifurcations of Non–smooth Dynamical Systems, Berlin, Springer, (2004).
  * [6] M. J. D. Powell, ”The convergence of variable metric methods for nonlinearly constrained optimization calculations”, Technical Memorandum No. 315, Argonne National Laboratory, Argonne, Illinois, 1977.
  * [7] R. Pytlak, Numerical Methods for Optimal Control Problem With State Constraints, Lecture Notes in Mathematics, vol. 1707, Springer–Verlag, (1999), .
  * [8] R. Pytlak, ”Numerical procedure for optimal control of higher index DAEs”, J. Discret. Dyn. Nat. Soc. (USA), vol. 29, pp. 1–24, (2011).
  * [9] R. Pytlak and D. Suski, ”Algorithms for optimal control of hybrid systems with sliding motion”, arXiv:2101.04754, (2021).
  * [10] R. Pytlak and D. Suski, ”Numerical procedure for optimal control of hybrid systems with sliding modes, Part I”, arXiv:2101.05829, (2021).
  * [11] R. Pytlak and D. Suski, ”On solving hybrid optimal control problems with higher index DAEs”, Optimization Methods & Software vol. 32, pp. 940–962, (2017).
  * [12] R. Pytlak and T. Zawadzki, ”On solving optimal control problems with higher index DAEs”, Optimization Methods & Software vol. 29, pp. 1139–1162, (2014).
  * [13] W. Stecz, R. Pytlak, A. Rymarz, S. Niemczyk, ”Application of dynamic optimisation for planning a haemodialysis process”, BMC Nephrology, pp. 1–11, (2019).
  * [14] D. Stewart and M. Anitescu, ”Optimal control of systems with discontinuous differential equations”, Num. Math., vol. 114, pp. 653-695, (2010).