---
abstract: 'Model-free control based on the idea of Reinforcement Learning is a promising approach that has recently gained extensive attention. However, Reinforcement-Learning-based control methods solely focus on the regulation problem or learn to track a reference that is generated by a time-invariant exo-system. In the latter case, controllers are only able to track the time-invariant reference dynamics which they have been trained on and need to be re-trained each time the reference dynamics change. Consequently, these methods fail in a number of applications which obviously rely on a trajectory not being generated by an exo-system. One prominent example is autonomous driving. This paper provides for the first time an adaptive optimal control method capable to track reference trajectories not being generated by a time-invariant exo-system. The main innovation is a novel Q-function that directly incorporates a given reference trajectory on a moving horizon. This new Q-function exhibits a particular structure which allows the design of an efficient, iterative, provably convergent Reinforcement Learning algorithm that enables optimal tracking. Two real-world examples demonstrate the effectiveness of our new method.'
address:
- 'Institute of Control Systems, Karlsruhe Institute of Technology (KIT), 76131 Karlsruhe, Germany'
- 'Intelligent Sensor-Actuator-Systems Laboratory, Institute for Anthropomatics and Robotics, Karlsruhe Institute of Technology, Karlsruhe, Germany'
author:
- Florian Köpf
- Johannes Westermann
- Michael Flad
- Sören Hohmann
bibliography:
- 'mybib.bib'
title: 'Adaptive Optimal Control for Reference Tracking Independent of Exo-System Dynamics'
---

Adaptive Dynamic Programming ,Optimal Tracking ,Reinforcement Learning ,Learning Systems ,Machine Learning ,Artificial Intelligence ,Optimal Control

Introduction
============

Reinforcement Learning (RL) has recently gained extensive attention as a model-free adaptive optimal control method which learns optimal behavior from interaction with the environment and observations of resulting states and rewards (see e.g. [@Sutton.2018; @Wang.2017; @Luo.2017] and references therein).

The existing control-oriented approaches can be classified in two groups: Some methods focus on the regulation problem, i.e. steering the system towards an equilibrium point (typically ${\boldsymbol{0}}$), e.g. [@Kamalapurkar.2016; @Luo.2017; @Vamvoudakis.2017; @Wang.2016; @Wei.2017]. The second group is devoted to a specific tracking case, where the reference ${\boldsymbol{r}}_k$ to be tracked is generated by an exo-system, i.e. it is assumed that there exists ${\boldsymbol{r}}_{k+1}={\boldsymbol{f}}_{\text{ref}}({\boldsymbol{r}}_k)$.

It is noteworthy that a controller that has learned to solve the regulation problem is not directly applicable to the tracking case as the reference trajectory influences the associated rewards that the RL agent is facing and therefore affects the value function which has to be learned.

In order to extend learning-based controllers to the tracking case, various RL-based tracking controllers belonging to the second group of approaches have been proposed recently. For unknown internal system dynamics but a known control coefficient matrix, adaptive tracking controllers approximating the value function are proposed for discrete-time [@Dierks.2009] and continuous-time [@Modares.2014; @Zhang.2017] systems. Notable critic-only Q-learning methods for completely unknown dynamics are proposed in [@Luo.2016] and [@Kiumarsi.2014], where [@Kiumarsi.2014] focuses on the linear-quadratic setting and [@Luo.2016] allows nonlinear system dynamics. Other works combine system identification and adaptive schemes including past information [@Lv.2017], use filter-based goal representation Heuristic Dynamic Programming for tracking [@Mu.2017b], focus on learning a tracking controller from input-output data rather than assuming full state information [@Kiumarsi.2015b], consider the nonzero-sum game case [@Kopf.2018], utilize the exo-system of a prescribed robot impedance model in order to learn a model-following behavior for assistive human-robot interaction [@Modares.2016], focus on systems with matched uncertainties [@Mu.2017] or consider tracking on an infinite horizon with unbounded cost [@Bernhard.2018].

All of these methods [@Dierks.2009; @Kiumarsi.2014; @Kiumarsi.2015b; @Kopf.2018; @Luo.2016; @Lv.2017; @Modares.2014; @Modares.2016; @Mu.2017; @Mu.2017b; @Zhang.2017; @Bernhard.2018] consider the case where the reference trajectory is generated by a time-invariant exo-system ${\boldsymbol{f}}_{\text{ref}}$, i.e. rely on the assumption that the reference ${\boldsymbol{r}}_k$ follows ${\boldsymbol{r}}_{k+1}={\boldsymbol{f}}_{\text{ref}}({\boldsymbol{r}}_k)$ (respectively $\dot{{\boldsymbol{r}}}(t)={\boldsymbol{f}}_{\text{ref}}({\boldsymbol{r}}(t))$ in the continuous-time case).

However, for various applications, this assumption does not hold (e.g. vehicles that follow a road, robots in sophisticated human-machine collaboration or specific time-varying sequences in process engineering). Consequently, present exo-system-based methods are not suited for these types of applications. This is because the learned parameters (i.e. Q-function or value function approximations and associated control laws) are corresponding to the system dynamics and reference dynamics ${\boldsymbol{f}}_{\text{ref}}$ and thus need to be re-learned as soon as any other reference $\tilde{{\boldsymbol{f}}}_{\text{ref}}\neq{\boldsymbol{f}}_{\text{ref}}$ should be followed. A method that tries to cope with this challenge is the multiple-model approach in [@Kiumarsi.2015], where an adaptive self-organizing map detects changes and switches between various learned models. However, in this approach new sub-models need to be trained whenever the reference dynamics ${\boldsymbol{f}}_{\text{ref}}$ changes.

An RL method that tracks reference trajectories not necessarily resulting from a time-invariant exo-system ${\boldsymbol{f}}_{\text{ref}}$ (or multiple switched models [@Kiumarsi.2015]) has not yet been proposed. Therefore, in contrast to [@Dierks.2009; @Kiumarsi.2014; @Kiumarsi.2015b; @Kopf.2018; @Luo.2016; @Lv.2017; @Modares.2014; @Modares.2016; @Mu.2017; @Mu.2017b; @Zhang.2017; @Bernhard.2018; @Kiumarsi.2015], our work provides for the first time an RL method allowing to track an arbitrary reference trajectory on a moving horizon.

Our general idea is to explicitly incorporate the reference trajectory on a moving horizon in our new Q-function that is learned without requiring the system dynamics. Due to considering not only the current tracking error but rather the reference trajectory, this Q-function yields a controller which does not only achieve reactive but predictive behavior.

In particular, we provide:

-   A novel moving horizon tracking Q-function whose minimizing control also minimizes the tracking costs.

-   Derivation and proof of the analytical solution of the novel Q-function given the cost function and the system dynamics for the linear-quadratic tracking case.

-   Convergence proofs for our learning algorithm that learns optimal tracking from data without knowledge of the system dynamics. Thus, the estimated Q-function parameters as well as the associated control law converge to the optimal solution.

The rest of this paper is organized as follows. In Section \[problem\_def\], we define the optimal tracking problem with unknown system dynamics. The novel Q-function for arbitrary reference tracking is defined and analyzed in Section \[qfunc\_section\]. In Section \[section\_learning\], we introduce our learning algorithm that is based on the previously defined Q-function and provide convergence proofs. Simulation results and a comparison of our method with an adaptive tracking controller that is trained on a reference generated by an exo-system ${\boldsymbol{f}}_{\text{ref}}$ are presented and discussed in Section \[sectionsimulation\] before the paper is concluded in Section \[conclusion\].

Problem Definition {#problem_def}
==================

Consider a discrete-time linear system $$\begin{aligned}
{\boldsymbol{x}}_{k+1} = {\boldsymbol{A}}{\boldsymbol{x}}_{k} + {\boldsymbol{B}}{\boldsymbol{u}}_{k}, \label{eq:System}\end{aligned}$$ where $k \in \mathbb{N}_0$ is the discrete time step, $\bm{x}_{k} \in \mathbb{R}^{n}$ the state vector and $\bm{u}_k \in \mathbb{R}^{m}$ the input vector. Both the system matrix ${\boldsymbol{A}}\in\mathbb{R}^{n\times n}$ and the input matrix ${\boldsymbol{B}}\in\mathbb{R}^{n\times m}$ are unknown, albeit the system $({\boldsymbol{A}}, {\boldsymbol{B}})$ is assumed to be controllable.

Assume that at each time step k, an arbitrary reference $\tilde{{\boldsymbol{r}}}_i^{(k)}\in\mathbb{R}^n$ is given on a moving horizon of length $N$, where $i$ is the time index starting at $k$. Beyond the horizon $N$, let the reference be ${\boldsymbol{0}}$. Thus, $$\begin{aligned}
\label{eq:ref}
{\boldsymbol{r}}_i^{(k)}=\begin{cases}\tilde{{\boldsymbol{r}}}_i^{(k)}, &\text{for } i=k,\dots,k+N \\ {\boldsymbol{0}}, &\text{for } i=k+N+1,\dots,\infty\end{cases}\end{aligned}$$ follows.

Our aim is to learn to track ${\boldsymbol{r}}_i^{(k)}$ optimally w.r.t. the quadratic cost $$\begin{aligned}
J_k &= \sum \limits_{i=k}^{\infty}\gamma^{i-k}\underbrace{\frac{1}{2}\left({\boldsymbol{e}}_i^\intercal{\boldsymbol{Q}}{\boldsymbol{e}}_i + {\boldsymbol{u}}_i^\intercal{\boldsymbol{R}}{\boldsymbol{u}}_i\right)}_{c\left({\boldsymbol{x}}_i,{\boldsymbol{u}}_i, {\boldsymbol{r}}_i\right)=c_i} \label{eq:Kostena}\end{aligned}$$ that has to be minimized. Note that we write ${\boldsymbol{r}}_i={\boldsymbol{r}}_i^{(k)}$ for brevity of notation in the following as $J_k$ indicates that the reference ${\boldsymbol{r}}_i^{(k)}$ starting at $k$ is used according to . Here, ${\boldsymbol{e}}_i={\boldsymbol{x}}_i-{\boldsymbol{r}}_i$ is the deviation of the system state ${\boldsymbol{x}}_i$ from the reference ${\boldsymbol{r}}_i$ at time step $i$. ${\boldsymbol{Q}}\in\mathbb{R}^{n\times n}$, ${\boldsymbol{Q}}={\boldsymbol{Q}}^\intercal\succeq{\boldsymbol{0}}$ is a symmetric, positive semidefinite matrix penalizing deviations of the state ${\boldsymbol{x}}_i$ from the reference ${\boldsymbol{r}}_i$, ${\boldsymbol{R}}\in\mathbb{R}^{m\times m}$ and ${\boldsymbol{R}}={\boldsymbol{R}}^\intercal\succ{\boldsymbol{0}}$ is a symmetric, positive definite matrix penalizing the control effort. Furthermore, [${\gamma \in\left[0,1\right)}$]{} a discount factor and $c_i \in \mathbb{R}$ denotes the one-step cost.

The choice of our cost function is motivated as follows: On one hand, the reference to be tracked is usually known only on a finite horizon $N$ (e.g. a road course). On the other hand, the Q-learning method that we use in Section \[section\_learning\] relies on an infinite horizon cost function which allows the efficient formulation of a Bellman-like equation. Thus, this cost function accounts for a given reference of finite horizon $N$ and enables stability and convergence of the associated RL algorithm due to the infinite horizon.

Therefore, our problem definition can be summarized in Problem \[problem1\].

\[problem1\] At each time step $k$, find the optimal control input ${\boldsymbol{u}}_k^*$ and apply it to system , where ${\boldsymbol{u}}_k^*, {\boldsymbol{u}}_{k+1}^*, \dots$ is the control sequence minimizing the discounted cost subject to the system dynamics , where ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ are unknown, given ${\boldsymbol{x}}_k$ and ${\boldsymbol{r}}_k, {\boldsymbol{r}}_{k+1}, \dots, {\boldsymbol{r}}_{k+N}$.

Extended Q-Function for Reference Tracking {#qfunc_section}
==========================================

Our idea is to define a new reference-dependent Q-function in contrast to the commonly used Q-function (i.e. we define a state-action-*reference* function rather than a state-action function). This Q-function is constructed such that its minimizing control input constitutes a solution to Problem \[problem1\].

In this section, we introduce the reference-dependent Q-function and derive its analytical solution. This provides important insights in how to parametrize the Q-function in Section \[section\_learning\], which is the basic ingredient for a convergent reinforcement learning algorithm capable of learning the optimal tracking solution.

We begin with the general case of a *finite* optimization horizon $K$, later we let $K\rightarrow \infty$. The notation seems to be a little bit clumsy but is of high importance, since optimization horizon $K$ and moving horizon of the reference $N$ should not be confused. Furthermore, the notation which control input is plugged into the one-step cost and Q-functions will be useful and we need to distinguish between the actual time step $k$ which the system is in and the time step $\kappa$ on the optimization horizon $K$.

\[def:Q\] Our proposed reference-dependent Q-function is defined as $$\begin{aligned}
    \tensor*[^K]{\!Q}{_{\kappa}}&=c_{k+\kappa}+\gamma\sum_{i=k+\kappa+1}^{k+K}\gamma^{i-(k+\kappa+1)}\left.c_i\right|_{{\boldsymbol{u}}_i^*} \nonumber \\ &= c_{k+\kappa}+\gamma \left.\tensor*[^K]{\!Q}{_{\kappa+1}}\right|_{{\boldsymbol{u}}^*_{k+\kappa+1}}\label{Q_def1}
    \end{aligned}$$ where $$\begin{aligned}
            c_i&=c\left({\boldsymbol{x}}_{i}, {\boldsymbol{u}}_{i}, {\boldsymbol{r}}_{i}\right),\label{eq:c}\\
        \tensor*[^K]{\!Q}{_{\kappa}}&=Q_{K-\kappa}\left({\boldsymbol{x}}_{k+\kappa}, {\boldsymbol{u}}_{k+\kappa}, {\boldsymbol{r}}_{k+\kappa}, \dots, {\boldsymbol{r}}_{k+K}\right),\label{eq:Q}\\
                \tensor*[^K]{\!Q}{_{K}}&=c\left({\boldsymbol{x}}_{k+K},{\boldsymbol{u}}_{k+K},{\boldsymbol{r}}_{k+K}\right)\label{eq:QNN}.
    \end{aligned}$$ Here, $\kappa\in\mathbb{N}_0$, $\kappa < K$ denotes the time step on the current optimization horizon of length $K$ starting at $k$ and ${\boldsymbol{r}}_i={\boldsymbol{0}}$, $\forall i>k+N$ (see ). The notation $\left.c_i\right|_{{\boldsymbol{u}}_i^*}$ indicates that the optimal control ${\boldsymbol{u}}_i^*$ is applied in and $\left.\tensor*[^K]{\!Q}{_{\kappa+1}}\right|_{{\boldsymbol{u}}^*_{k+\kappa+1}}$ denotes that ${\boldsymbol{u}}_{k+\kappa+1}={\boldsymbol{u}}^*_{k+\kappa+1}$ in $\tensor*[^K]{\!Q}{_{\kappa+1}}$ (cf. ).

Therefore, $\tensor*[^K]{\!Q}{_{\kappa}}$ is the accumulated discounted cost from time step $k+\kappa$ to $k+K$ if the control ${\boldsymbol{u}}_{k+\kappa}$ is applied at time step $k+\kappa$ and the optimal controls ${\boldsymbol{u}}^*_{k+\kappa+1}, \dots, {\boldsymbol{u}}^*_{k+K}$ minimizing the cost-to-go are applied thereafter while the reference is provided on a moving horizon of length $N$. With the finite horizon cost function defined as $$\begin{aligned}
\label{eq:JkK}
\tensor*[^K]{\!J}{_k} = \sum_{i=k}^{k+K}\gamma^{i-k}\frac{1}{2}\left({\boldsymbol{e}}_i^\intercal{\boldsymbol{Q}}{\boldsymbol{e}}_i + {\boldsymbol{u}}_i^\intercal{\boldsymbol{R}}{\boldsymbol{u}}_i\right)\end{aligned}$$ the subsequent Lemma \[lemmauQisuJ\] follows.

\[lemmauQisuJ\] The control ${\boldsymbol{u}}_k$ minimizing the reference-dependent Q-function $\tensor*[^K]{\!Q}{_{0}}$ is a solution for ${\boldsymbol{u}}_k^*$ minimizing $\tensor*[^K]{\!J}{_k}$.

With and $$\begin{aligned}
\label{eq:minuQ}
    \min_{{\boldsymbol{u}}_k}\tensor*[^K]{\!Q}{_{0}}&=\left.c_k\right|_{{\boldsymbol{u}}^*_k}+\gamma \left.\tensor*[^K]{\!Q}{_{1}}\right|_{{\boldsymbol{u}}^*_{k+1}}\nonumber \\ &=\left.\tensor*[^K]{\!Q}{_{0}}\right|_{{\boldsymbol{u}}^*_k}=\sum_{i=k}^{k+K}\gamma^{i-k}\left.c_i\right|_{{\boldsymbol{u}}^*_i}
    \end{aligned}$$ follows $$\begin{aligned}
    \min_{{\boldsymbol{u}}_{k}, \dots, {\boldsymbol{u}}_{k+K}}\tensor*[^K]{\!J}{_k}&=\sum_{i=k}^{k+K}\gamma^{i-k}\left.c_i\right|_{{\boldsymbol{u}}^*_i}=\left.\tensor*[^K]{\!Q}{_{0}}\right|_{{\boldsymbol{u}}^*_k}.
    \end{aligned}$$

Taking the limit $$\begin{aligned}
\lim_{K\rightarrow\infty}\tensor*[^K]{\!J}{_k}=J_k\end{aligned}$$ yields $J_k$ (cf. ). With the Q-function $\tensor*[^K]{\!Q}{_{0}}$ defined according to Definition \[def:Q\] and as a result of Lemma \[lemmauQisuJ\], Problem \[problem1\] is equivalent to the following Problem \[problem2\].

\[problem2\] Given ${\boldsymbol{r}}_k, {\boldsymbol{r}}_{k+1}, \dots, {\boldsymbol{r}}_{k+N}$, in each state ${\boldsymbol{x}}_k$ at time step $k$, find the control ${\boldsymbol{u}}_k^*$ minimizing the reference-dependent Q-function $Q_{0}$, where $$\begin{aligned}
    Q_{0} = \lim_{K\rightarrow\infty}\tensor*[^K]{\!Q}{_{0}}
    \end{aligned}$$ and apply it to the system whose matrices ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ are unknown.

We will proceed in two steps. First, it is assumed that the system matrices ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ are known. Later on, this assumption will be dropped and an iterative solution based on a temporal difference error will be introduced. We will now propose the analytical solution of $\tensor*[^K]{\!Q}{_{0}}$ in Theorem \[Theo:Q\].

\[Theo:Q\] For $K\geq N$, the Q-function $\tensor*[^K]{\!Q}{_{0}}$ (cf. Definition \[def:Q\]) with the objective function is given by $$\begin{aligned}
    \tensor*[^K]{\!Q}{_{0}} = \frac{1}{2}{\begin{bmatrix*}[r]{\boldsymbol{x}}_k^\intercal & {\boldsymbol{u}}_k^\intercal & {\boldsymbol{r}}_k^\intercal & \cdots & {\boldsymbol{r}}_{k+N}^\intercal & {\boldsymbol{0}}_{\vphantom{k+N}}^\intercal
\end{bmatrix*}}{\boldsymbol{H}}_K{\begin{bmatrix*}[c]{\boldsymbol{x}}_k\\{\boldsymbol{u}}_k\\{\boldsymbol{r}}_k\\ \vdots \\ {\boldsymbol{r}}_{k+N}\\ {\boldsymbol{0}}
\end{bmatrix*}},\label{eq:Q0_mit_H}
    \end{aligned}$$ where ${\boldsymbol{H}}_K={\boldsymbol{H}}_K^\intercal\in\mathbb{R}^{\left((K+2)n+m\right)\times\left((K+2)n+m\right)}$ with $$\begin{aligned}
    {\boldsymbol{H}}_K = {\begin{bmatrix*}[l]{\boldsymbol{h}}_{xx}      & {\boldsymbol{h}}_{xu}        & {\boldsymbol{h}}_{xr_0}      & {\boldsymbol{h}}_{xr_1}      & {\boldsymbol{h}}_{xr_2}      & \cdots & {\boldsymbol{h}}_{xr_{K}}   \\ \label{eq:FormDerQ-Fkt_2}
        {\boldsymbol{h}}_{ux}          & {\boldsymbol{h}}_{uu}        & {\boldsymbol{0}}                 & {\boldsymbol{h}}_{ur_1}      & {\boldsymbol{h}}_{ur_2}      & \cdots & {\boldsymbol{h}}_{ur_{K}}   \\
        {\boldsymbol{h}}_{r_0x}        & {\boldsymbol{0}}                 & {\boldsymbol{h}}_{r_0r_0}    & {\boldsymbol{0}}                 & {\boldsymbol{0}}                 & \cdots & {\boldsymbol{0}}                \\
        {\boldsymbol{h}}_{r_1x}        & {\boldsymbol{h}}_{r_1u}      & {\boldsymbol{0}}                 & {\boldsymbol{h}}_{r_1r_1}    & {\boldsymbol{0}}                 & \cdots & {\boldsymbol{0}}                \\
        {\boldsymbol{h}}_{r_2x}        & {\boldsymbol{h}}_{r_2u}      & {\boldsymbol{0}}                 & {\boldsymbol{0}}                 & {\boldsymbol{h}}_{r_2r_2}    & \cdots & {\boldsymbol{h}}_{r_{2}r_{K}} \\ 
        {\boldsymbol{h}}_{r_3x}        & {\boldsymbol{h}}_{r_3u}      & {\boldsymbol{0}}                 & {\boldsymbol{0}}                 & {\boldsymbol{h}}_{r_3r_2}    & \cdots & {\boldsymbol{h}}_{r_3r_{K}}     \\
        \,\vdots                & \,\vdots          & \,\vdots          & \,\vdots          & \,\vdots          & \ddots & \,\vdots             \\ 
        {\boldsymbol{h}}_{r_{K}x}  & {\boldsymbol{h}}_{r_{K}u}    &{\boldsymbol{0}}                  & {\boldsymbol{0}}                 & {\boldsymbol{h}}_{r_{K}r_2}& \cdots & {\boldsymbol{h}}_{r_{K}r_{K}} 
\end{bmatrix*}}.
    \end{aligned}$$ The exact values of ${\boldsymbol{H}}_K$ follow from the subsequent proof.

The proof is of rather technical nature and given in \[append:notations\].

For $K\rightarrow \infty$ let ${\boldsymbol{H}}$ be the northwestern $\left((N+2)n+m\right)\times\left((N+2)n+m\right)$-submatrix of ${\boldsymbol{H}}_K$. Then, $$\begin{aligned}
Q_{0} &= \lim_{K\rightarrow\infty}\tensor*[^K]{\!Q}{_{0}}=\frac{1}{2}{\boldsymbol{z}}_k^\intercal{\boldsymbol{H}}{\boldsymbol{z}}_k\label{eq:H_N_steps}\end{aligned}$$ follows, where ${\boldsymbol{z}}_k={\begin{bmatrix*}[r]{\boldsymbol{x}}_k^\intercal & {\boldsymbol{u}}_k^\intercal & {\boldsymbol{r}}_k^\intercal & \cdots & {\boldsymbol{r}}_{k+N}^\intercal
\end{bmatrix*}}^\intercal$, as ${{\boldsymbol{r}}_i = {\boldsymbol{0}}}$, $\forall i>k+N$.

Thus, the Q-function for the LQ tracking problem is quadratic w.r.t. the state ${\boldsymbol{x}}_k$, the control input ${\boldsymbol{u}}_k$ and the reference ${\boldsymbol{r}}_k, \dots, {\boldsymbol{r}}_{k+N}$ and furthermore completely parametrized by ${\boldsymbol{H}}$. This obviously is a generalization of the well known structure for reference free RL-procedures [@Bradtke.1994; @Landelius.1997]. In addition, ${\boldsymbol{H}}$ is not only quadratic but has a specific structure, which is a new result and allows an efficient parametrization for the Q-learning based algorithm in the following.

Q-Learning Based Tracking {#section_learning}
=========================

We use the new Q-function for reference tracking on a moving horizon where the system matrices ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ are unknown. Thus, our aim is to determine ${\boldsymbol{H}}$ (cf. ) by observations of states and rewards.

The optimal control ${\boldsymbol{u}}_k^*$, which is equivalent to for $\kappa=0$ and $K\rightarrow\infty$, can be expressed directly by means of ${\boldsymbol{H}}$ and is given by Corollary \[corollary:uoptH\]. Here, ensures that ${\boldsymbol{u}}_k^*$ in fact minimizes the Q-function $Q_{0}$.

\[corollary:uoptH\] With Lemma \[lemmauQisuJ\] and Theorem \[Theo:Q\], the optimal control at time step $k$ is given by $$\begin{aligned}
\label{eq:uoptH}
    {\boldsymbol{u}}_{k}^* = -{\boldsymbol{h}}_{uu}^{-1}{\begin{bmatrix*}[r]{\boldsymbol{h}}_{ux} & {\boldsymbol{h}}_{ur_1}  & \cdots & {\boldsymbol{h}}_{ur_{N}}
\end{bmatrix*}} {\begin{bmatrix*}[c]{\boldsymbol{x}}_{k} \\ {\boldsymbol{r}}_{k+1} \\ \vdots \\ {\boldsymbol{r}}_{k+N}
\end{bmatrix*}}.
    \end{aligned}$$

Hence, if the Q-function $Q_{0}$ is known by means of the matrix ${\boldsymbol{H}}$, the optimal control directly results from .

In contrast to usual controllers learned by Reinforcement Learning, our control law explicitly depends on the reference values ${\boldsymbol{r}}_{k+1},\dots,{\boldsymbol{r}}_{k+N}$ and is therefore able to generalize to arbitrary reference trajectories on the horizon $N$.

Parametrization of the Reference-dependent Q-function
-----------------------------------------------------

In order to learn the Q-function, we parametrize $Q_{0}$ and perform a value iteration on the resulting squared Bellman-like temporal difference (TD) error in order to estimate the Q-function parameters as well as the corresponding optimal control law. Let the estimated Q-function be parametrized by means of a sum of weighted basis functions: $$\begin{aligned}
\label{eq:Qhat}
\hat{Q}_{0}&={\boldsymbol{w}}^\intercal{\boldsymbol{\phi}}\left({\boldsymbol{x}}_k, {\boldsymbol{u}}_k, {\boldsymbol{r}}_k, \dots, {\boldsymbol{r}}_{k+N}\right)={\boldsymbol{w}}^\intercal{\boldsymbol{\phi}}({\boldsymbol{z}}_k).\end{aligned}$$ Here, ${\boldsymbol{w}}\in\mathbb{R}^{{L}}$ is a weight vector and ${\boldsymbol{\phi}}:\mathbb{R}^{(N+2)n+m}\rightarrow\mathbb{R}^{{L}}$ is a vector of basis functions. Note that, in contrast to usual Q-function approximations, $\hat{Q}_{0}$ in explicitly incorporates the reference trajectory ${\boldsymbol{r}}_k, {\boldsymbol{r}}_{k+1}, \dots, {\boldsymbol{r}}_{k+N}$.

\[lemma:Qexakt\_long\] With $$\begin{aligned}
\label{eq:L}
    L=&~\frac{1}{2}\left((N+2)n+m\right)\left((N+2)n+m+1\right)\nonumber\\&~-\left(n^2(2N-1)+mn\right)
    \end{aligned}$$ quadratic basis functions ${{\boldsymbol{\phi}}={\begin{bmatrix*}[r]\phi_1 & \cdots & \phi_L
\end{bmatrix*}}^\intercal}$, there exists a weight vector ${\boldsymbol{w}}={\boldsymbol{w}}^*$ such that $\hat{Q}_{0}=Q_{0}$.

Due to the symmetry and the zeros in ${\boldsymbol{H}}_K$ (cf. ) and therefore also in ${\boldsymbol{H}}$, there are $L$ non-redundant elements in ${\boldsymbol{H}}$. Define quadratic basis functions $\phi_l$, $l=1,\dots,L$ of the form $$\begin{aligned}
    \phi_l= \begin{cases}\left\{{\boldsymbol{z}}_k\right\}_i\left\{{\boldsymbol{z}}_k\right\}_j, &\text{for } i\neq j \\ \frac{1}{2}\left\{{\boldsymbol{z}}_k\right\}_i^2, &\text{for } i=j,\end{cases}
    \end{aligned}$$ where $i, j$ indicate the corresponding non-redundant elements of ${\boldsymbol{H}}$, $\left\{{\boldsymbol{\cdot}}\right\}_i$ denotes the $i$-th element of a vector and ${\boldsymbol{z}}_k$ is defined as in . Thus, $Q_{0}$ in is equivalent to $\hat{Q}_{0}$ in if the weights $\left\{{\boldsymbol{w}}\right\}_i$, $i=1, \dots, L$, are equal to the corresponding non-redundant elements in ${\boldsymbol{H}}$ which we denote by ${\boldsymbol{w}}={\boldsymbol{w}}^*$.

Although $L$ in Lemma \[lemma:Qexakt\_long\] gives the maximum number of weights needed in order to parametrize $Q_{0}$ exactly, in the frequently occurring case of a sparse weighting matrix ${\boldsymbol{Q}}$ in the cost functional , the exact knowledge of the structure of ${\boldsymbol{H}}$ can be exploited further in order to drastically reduce the weights ${\boldsymbol{w}}$ that are necessary.

\[lemma:Qexakt\] If the $l$-th row and $l$-th column of ${\boldsymbol{Q}}$ equals zero, then the

-   $l$-th column of ${\boldsymbol{h}}_{xr_i} \forall i \in \left\{0,\dots,N\right\}$,

-   $l$-th column of ${\boldsymbol{h}}_{ur_i} \forall i \in \left\{0,\dots,N\right\}$ and

-   $l$-th row and $l$-th column of ${\boldsymbol{h}}_{r_ir_j} \forall i,j \in \left\{0,\dots,N\right\}$

are all equal to zero. Thus, the number of non-redundant weights $L$ in ${\boldsymbol{H}}$ (corresponding to the weight vector ${\boldsymbol{w}}$) reduces to $$\begin{aligned}
\label{eq:L_sparse}
    L &= (n-q)(n-q+1)\left(\frac{N}{2}+1\right)+\frac{1}{2}n(n+1)\\ \nonumber&+\left(m+N(n-q)\right)(m+n)+\frac{1}{2}(N-2)(N-1)(n-q)^2,
    \end{aligned}$$ where ${\boldsymbol{h}}_{xr_0}={\boldsymbol{Q}}$ and ${\boldsymbol{h}}_{r_0r_0}=-{\boldsymbol{Q}}$ has been considered and $q$ denotes the number of rows and columns of ${\boldsymbol{Q}}$ that are both zero.

The rather technical proof directly follows from considering –.

Although Lemma \[lemma:Qexakt\] is of technical nature, sparsity of ${\boldsymbol{H}}$ is essential in order to implement efficient learning controllers as will be discussed in Section \[sectionsimulation\].

Based on the specific structure of the analytical solution of $Q_{0}$ derived in Theorem \[Theo:Q\], according to Lemma \[lemma:Qexakt\], $L$ quadratic basis functions are sufficient to parametrize $Q_{0}$ exactly.

Equation and Lemma \[lemma:Qexakt\] show that our proposed Q-function generalizes over reference trajectories as the weight vector ${\boldsymbol{w}}$ does not depend on a specific reference we intend to follow.

Online Learning Algorithm
-------------------------

In this section, we propose our value iteration based algorithm in order to learn the weights ${\boldsymbol{w}}$ online by minimizing the squared temporal difference error. Let $$\begin{aligned}
\label{eq:Qhat1}
\hat{Q}_{1}&={\boldsymbol{w}}^\intercal{\boldsymbol{\phi}}\left({\boldsymbol{x}}_{k+1}, {\boldsymbol{u}}_{k+1}, {\boldsymbol{r}}_{k+1}, \dots, {\boldsymbol{r}}_{k+N}, {\boldsymbol{0}}\right)\nonumber \\&={\boldsymbol{w}}^\intercal{\boldsymbol{\phi}}\left({\boldsymbol{z}}_{k+1}\right),\end{aligned}$$ where ${\boldsymbol{z}}_{k+1}={\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+1}^\intercal & {\boldsymbol{u}}_{k+1}^\intercal & {\boldsymbol{r}}_{k+1}^\intercal & \cdots & {\boldsymbol{r}}_{k+N}^\intercal & {\boldsymbol{0}}_{\vphantom{k+N}}^\intercal
\end{bmatrix*}}^\intercal$.

The estimated optimal control input $\hat{{\boldsymbol{u}}}_{k+\kappa}^*$ based on the estimated Q-function $\left.\hat{Q}_{\kappa}\right|_{{\boldsymbol{u}}_{k+\kappa}}$ is defined as $$\begin{aligned}
\label{eq:opt_u_estimated}
\hat{{\boldsymbol{u}}}_{k+\kappa}^*&=\operatorname*{arg\,min}_{{\boldsymbol{u}}_{k+\kappa}}\left.\hat{Q}_{\kappa}\right|_{{\boldsymbol{u}}_{k+\kappa}}\nonumber \\ &={\boldsymbol{L}}{\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa}^\intercal&{\boldsymbol{r}}_{k+\kappa+1}^\intercal & \dots & {\boldsymbol{r}}_{k+\kappa+N}^\intercal
\end{bmatrix*}}^\intercal,\end{aligned}$$ where ${\boldsymbol{L}}=-\hat{{\boldsymbol{h}}}_{uu}^{-1}{\begin{bmatrix*}[r]\hat{{\boldsymbol{h}}}_{ux} & \hat{{\boldsymbol{h}}}_{ur_1}  & \cdots & \hat{{\boldsymbol{h}}}_{ur_{N}}
\end{bmatrix*}}$ (cf. ). Here, let $\hat{{\boldsymbol{H}}}$ be the estimated matrix ${\boldsymbol{H}}$ based on ${\boldsymbol{w}}$ and $\hat{{\boldsymbol{h}}}$ denote the submatrices of $\hat{{\boldsymbol{H}}}$ as in . Furthermore, $\hat{{\boldsymbol{u}}}_{k+\kappa}^*={\boldsymbol{u}}_{k+\kappa}^*$ results if ${\boldsymbol{w}}={\boldsymbol{w}}^*$ according to Lemma \[lemma:Qexakt\] and the optimal control law ${\boldsymbol{L}}^*$ results if $\hat{{\boldsymbol{H}}}={\boldsymbol{H}}$.

Then, based on the Bellman-like equation , the TD error $\epsilon_k$ [@Sutton.1988] is given in Definition \[def:TD\]. If ${{\boldsymbol{w}}={\boldsymbol{w}}^*}$, $\epsilon_k$ would vanish.

\[def:TD\] The temporal difference error $\epsilon_k$, i.e. the approximation error of the Bellman-like equation due to the deviation of the weight estimate ${\boldsymbol{w}}$ from ${\boldsymbol{w}}^*$ is defined as $$\begin{aligned}
\label{eq:td}
    \epsilon_k&=c_k+\gamma\left.\hat{Q}_{1}\right|_{\hat{{\boldsymbol{u}}}_{k+1}^*}-\hat{Q}_{0}\nonumber \\ &=c_k+\gamma{\boldsymbol{w}}^\intercal \phi\left({\boldsymbol{z}}_{k+1}^*\right)-{\boldsymbol{w}}^\intercal \phi\left({\boldsymbol{z}}_{k}\right).
    \end{aligned}$$ Here, ${\boldsymbol{z}}_{k+1}^* = \left.{\boldsymbol{z}}_{k+1}\right|_{{\boldsymbol{u}}_{k+1}=\hat{{\boldsymbol{u}}}_{k+1}^*}$.

In order to improve the estimated reference-dependent Q-function $\hat{Q}_{0}$ as well as the resulting estimated optimal control $\hat{{\boldsymbol{u}}}_k^*$, we employ a value iteration procedure. This iteration consists of a *policy evaluation* which updates the weight estimate ${\boldsymbol{w}}^{(i)}$ representing the Q-function and a *policy improvement* step, where, based on the updated Q-function weight ${\boldsymbol{w}}^{(i)}$ and the corresponding matrix ${\boldsymbol{H}}^{(i)}$, the control law ${\boldsymbol{L}}^{(i)}$ is adapted according to .

To evaluate $\epsilon_k$ in , $\hat{{\boldsymbol{u}}}_{k+1}^*$ is required. However, as the optimal weight ${\boldsymbol{w}}^*$ is unknown a priori, we initialize ${\boldsymbol{w}}^{(0)}={\boldsymbol{0}}$ and the estimated optimal control by $\hat{{\boldsymbol{u}}}_{k+1}^{*(0)} = {\boldsymbol{0}}$ which is achieved by setting ${\boldsymbol{L}}^{(0)}={\boldsymbol{0}}$.

In the *policy evaluation* step, the aim is to find an updated ${\boldsymbol{w}}^{(i+1)}$ such that $\left({\epsilon}_k^{(i)}\right)^2$ is minimized, where $$\begin{aligned}
\label{eq:td_VI}
{\epsilon}_k^{(i)} &= c_k + \gamma{{\boldsymbol{w}}^{(i)}}^\intercal{\boldsymbol{\phi}}\left({\boldsymbol{z}}_{k+1}^{*{(i)}}\right)-{{\boldsymbol{w}}^{(i+1)}}^\intercal{\boldsymbol{\phi}}\left({\boldsymbol{z}}_{k}\right)\end{aligned}$$ in analogy to . In accordance with Lemma \[lemma:Qexakt\], ${\boldsymbol{w}}\in\mathbb{R}^L$ follows. Thus, ${\epsilon}_k^{(i)}$ needs to be considered at $M\geq L$ time steps in order to perform a least-squares update, where $M$ is the number of samples used for the policy evaluation. Then, ${\boldsymbol{w}}^{(i+1)}$ results from $$\begin{aligned}
\label{eq:PE}
{\boldsymbol{w}}^{(i+1)}=\operatorname*{arg\,min}_{{\boldsymbol{w}}^{(i+1)}}\sum_{j=k-M+1}^{k}\left({\epsilon}_j^{(i)}\right)^2.\end{aligned}$$ Now, we define $$\begin{aligned}
\label{eq:Phimatrix}
{\boldsymbol{\Phi}}&={\begin{bmatrix*}[c]{\boldsymbol{\phi}}\left({\boldsymbol{z}}_{k-M+1}\right)& \dots & {\boldsymbol{\phi}}\left({\boldsymbol{z}}_{k}\right)
\end{bmatrix*}}^\intercal\end{aligned}$$ and $$\begin{aligned}
\label{eq:cmatrix}
{\boldsymbol{c}}={\begin{bmatrix*}[c]c_{k-M+1}+\gamma{{\boldsymbol{w}}^{(i)}}^\intercal\phi\left({\boldsymbol{z}}_{k-M+2}^{*(i)}\right)\\ \vdots \\ c_k+\gamma{{\boldsymbol{w}}^{(i)}}^\intercal\phi\left({\boldsymbol{z}}_{k+1}^{*(i)}\right)
\end{bmatrix*}}.\end{aligned}$$ If the excitation condition $$\begin{aligned}
\label{rankcondition}
\operatorname{rank}\left({{\boldsymbol{\Phi}}^\intercal}{\boldsymbol{\Phi}}\right) = L\end{aligned}$$ is satisfied, ${\boldsymbol{w}}^{(i+1)}$ minimizing exists, is unique and given by $$\begin{aligned}
\label{eq:LS}
{\boldsymbol{w}}^{(i+1)}&=\left({{\boldsymbol{\Phi}}^\intercal}{\boldsymbol{\Phi}}\right)^{-1}{{\boldsymbol{\Phi}}^\intercal}{\boldsymbol{c}}\end{aligned}$$ (cf. [@Astrom.1995 Theorem 2.1]).

Then, the *policy improvement* step is based on the new weight estimate ${\boldsymbol{w}}^{(i+1)}$ and its corresponding ${\boldsymbol{H}}^{(i+1)}$ and results in $$\begin{aligned}
\label{eq:policy_improvement}
{\boldsymbol{L}}^{(i+1)}=-{{\boldsymbol{h}}_{uu}^{(i+1)}}^{-1}{\begin{bmatrix*}[r]{\boldsymbol{h}}_{ux}^{(i+1)} & {\boldsymbol{h}}_{ur_1}^{(i+1)}  & \cdots & {\boldsymbol{h}}_{ur_{N}^{(i+1)}}
\end{bmatrix*}}\end{aligned}$$ (cf. Corollary \[corollary:uoptH\]).

With the time step still being fixed at $k$, this iteration is performed until the change in ${\boldsymbol{w}}$ stays below a given threshold $e_{{\boldsymbol{w}}}$. Although we would like to evaluate the Q-function corresponding to the target policy $\hat{{\boldsymbol{u}}}_k^*$ (note that ${\boldsymbol{z}}_{k+1}^{*(i)}$ is used in ), we have not yet discussed how to choose the behavior policy, i.e. the control ${\boldsymbol{u}}_k$ that is applied to the system and appears in $c_k$ and ${\boldsymbol{z}}_k$ (cf. ). During the learning process is active, let $$\begin{aligned}
\label{utilde}
{\boldsymbol{u}}_k=\tilde{{\boldsymbol{u}}}_k^*=\hat{{\boldsymbol{u}}}_k^*+{\boldsymbol{\xi}}.\end{aligned}$$ Here, the Gaussian noise ${\boldsymbol{\xi}} \sim \mathcal{N}_{m}\left({\boldsymbol{0}},{\boldsymbol{\Sigma}}\right)$ serves as exploration noise as persistent excitation is required for convergence (cf. , [@Narendra.1987; @Astrom.1995; @Lewis.2009]). Furthermore, if the reference to track is smooth, additional excitation noise should be added to the reference in order to satisfy condition . The complete algorithm is shown in Algorithm \[algorithm\] and learns to track an arbitrary reference that is given on a moving horizon of length $N$ without knowledge of the system matrices ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$. If the system dynamics is not expected to change over time, learning might be stopped after the first complete value iteration based on $M$ data tuples has been performed. This can be done by stopping Algorithm \[algorithm\] after line 11. In practice, learning might be enabled again whenever the Bellman error increases which is an indicator for suboptimal weights ${\boldsymbol{w}}$.

apply $\tilde{{\boldsymbol{u}}}_k^*$ to the system $i=0, {\boldsymbol{w}}^{(i)}={\boldsymbol{w}}$ policy evaluation: ${\boldsymbol{w}}^{(i+1)}$ policy improvement: ${\boldsymbol{L}}^{(i+1)}$ $i = i+1$ ${\boldsymbol{w}}={\boldsymbol{w}}^{(i)}$

The iterative procedure of policy evaluation and policy improvement in Algorithm \[algorithm\] is a value iteration (cf. [@Sutton.2018]). This is due to the definition of ${\epsilon}_k$ in , where $\hat{Q}_{1}$ relies on ${\boldsymbol{w}}^{(i)}$ and $\hat{Q}_{0}$ on ${\boldsymbol{w}}^{(i+1)}$ (cf. [@Lewis.2009]).

Just as regular Q-learning [@Watkins.1992], our algorithm belongs to the off-policy RL methods [@Sutton.2018] as the behavior policy $\tilde{{\boldsymbol{u}}}_k^*=\hat{{\boldsymbol{u}}}_k^*+{\boldsymbol{\xi}}$ is followed while the agent learns the Q-function belonging to the target policy $\hat{{\boldsymbol{u}}}_k^*$.

Convergence Analysis of the Learning Algorithm {#convergence}
----------------------------------------------

In this section, we will provide convergence proofs, i.e. show that the estimated reference-dependent Q-function $\hat{Q}_{0}$ converges to the underlying Q-function $Q_{0}$ and that ${{\boldsymbol{w}}^{(i)}\rightarrow {\boldsymbol{w}}^*}$, i.e. ${\boldsymbol{H}}^{(i)}\rightarrow {\boldsymbol{H}}$ as $i\rightarrow \infty$. This also implies that ${\boldsymbol{L}}^{(i)}\rightarrow {\boldsymbol{L}}^*$, i.e. the value iteration converges to the optimal control law.

Our convergence analysis is structured as follows. First, we prove that the value iteration (i.e. iterating between and ) is equivalent to a matrix sequence on ${\boldsymbol{H}}^{(i)}$. Second, we prove that this sequence of matrices is upper bounded in the sense that ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq {\boldsymbol{Y}}$ while ${\boldsymbol{H}}^{(i)}$ is monotonically increasing, i.e. ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{H}}^{(i+1)}$. Hence, the sequence converges. Finally, we show that the converged sequence fulfills the Bellman equation and the corresponding control law is optimal.

The following Lemma \[lemma:VI\_matrix\_iteration\] extends [@AlTamimi.2007 Lemma 1] to the tracking case and shows that our proposed value iteration is equivalent to a matrix sequence on ${\boldsymbol{H}}^{(i)}$.

\[lemma:VI\_matrix\_iteration\] Let ${\boldsymbol{H}}^{(0)} = {\boldsymbol{0}}$, ${\boldsymbol{R}}^\intercal = {\boldsymbol{R}} \succ {\boldsymbol{0}}$, ${\boldsymbol{Q}}^\intercal = {\boldsymbol{Q}} \succeq {\boldsymbol{0}}$ and $({\boldsymbol{A}}, {\boldsymbol{B}})$ controllable. The value iteration described by and is equivalent to the iteration $$\begin{aligned}
\label{eq:iteration_Hi1}
    {\boldsymbol{H}}^{(i+1)}={\boldsymbol{G}}+\gamma{\boldsymbol{M}}\!\left({\boldsymbol{L}}^{(i)}\right)^\intercal {\boldsymbol{H}}^{(i)}{\boldsymbol{M}}\!\left({\boldsymbol{L}}^{(i)}\right),
    \end{aligned}$$ where $$\begin{aligned}
    {\boldsymbol{G}}={\begin{bmatrix*}[r]{\boldsymbol{Q}}&{\boldsymbol{0}}&-{\boldsymbol{Q}}&{\boldsymbol{0}}\\ {\boldsymbol{0}}&{\boldsymbol{R}}&{\boldsymbol{0}}&{\boldsymbol{0}}\\ -{\boldsymbol{Q}}&{\boldsymbol{0}}&{\boldsymbol{Q}}&{\boldsymbol{0}}\\ {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}
\end{bmatrix*}}
    \end{aligned}$$ and $$\begin{aligned}
\label{eq:M}
    {\boldsymbol{M}}\!\left({\boldsymbol{L}}^{(i)}\right)&={\begin{bmatrix*}[c]{\boldsymbol{A}}&{\boldsymbol{B}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&\cdots&{\boldsymbol{0}}\\
        {\boldsymbol{L}}_x^{(i)}{\boldsymbol{A}}&{\boldsymbol{L}}_x^{(i)}{\boldsymbol{B}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{L}}_1^{(i)}&\cdots&{\boldsymbol{L}}_{N-1}^{(i)}\\
        {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{I}}_n&{\boldsymbol{0}}&\cdots&{\boldsymbol{0}}\\
        {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{I}}_n&\cdots&{\boldsymbol{0}}\\
        \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\
        {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&\cdots&{\boldsymbol{I}}_n\\
        {\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&{\boldsymbol{0}}&\cdots&{\boldsymbol{0}}
\end{bmatrix*}}
    \end{aligned}$$ with ${\boldsymbol{L}}^{(i)}={\begin{bmatrix*}[r]{\boldsymbol{L}}^{(i)}_x&{\boldsymbol{L}}^{(i)}_1&\cdots&{\boldsymbol{L}}^{(i)}_N
\end{bmatrix*}}={\boldsymbol{L}}\!\left({\boldsymbol{H}}^{(i)}\right)$ depending on ${\boldsymbol{H}}^{(i)}$ as in .

See \[append:VI\_matrix\_iteration\].

The following technical Lemma \[lemma:hilfslemma\] is required later for the proof of Lemma \[lemma:monotonically\_increasing\].

\[lemma:hilfslemma\] For ${\boldsymbol{H}}^{(0)} = {\boldsymbol{0}}$, ${\boldsymbol{R}}^\intercal = {\boldsymbol{R}} \succ {\boldsymbol{0}}$ and ${\boldsymbol{Q}}^\intercal = {\boldsymbol{Q}} \succeq {\boldsymbol{0}}$, $\forall i > 0$: ${\boldsymbol{L}}^{(i)}{\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa}^\intercal&{\boldsymbol{r}}_{k+\kappa+1}^\intercal & \dots & {\boldsymbol{r}}_{k+\kappa+N}^\intercal
\end{bmatrix*}}$ is the unique minimizer of $$\begin{aligned}
    \hat{Q}_{\kappa}^{(i)}=\frac{1}{2}{\boldsymbol{z}}_{k+\kappa}^\intercal{\boldsymbol{H}}^{(i)}{\boldsymbol{z}}_{k+\kappa}.
    \end{aligned}$$

Due to ${\boldsymbol{H}}^{(0)}={\boldsymbol{0}}\succeq{\boldsymbol{0}}$, ${\boldsymbol{R}} \succ {\boldsymbol{0}}$ and ${\boldsymbol{Q}} \succeq {\boldsymbol{0}}$, yields ${\boldsymbol{H}}^{(i)}\succeq{\boldsymbol{0}}$. Furthermore, due to ${\boldsymbol{R}} \succ {\boldsymbol{0}}$ it is obvious that ${\boldsymbol{h}}_{uu}^{(i)}\succ{\boldsymbol{0}}$, $\forall i>0$. Therefore, $\frac{\partial \hat{Q}_{\kappa}^{(i)}}{\partial {\boldsymbol{u}}_{k+\kappa}}={\boldsymbol{0}}$ yields ${\boldsymbol{L}}^{(i)}$ and $\frac{\partial ^2 \hat{Q}_{\kappa}^{(i)}}{\partial {\boldsymbol{u}}_{k+\kappa}^2}\succ{\boldsymbol{0}}$ follows due to ${\boldsymbol{h}}_{uu}^{(i)}\succ{\boldsymbol{0}}$ and completes the proof.

Define the operator $$\begin{aligned}
\label{eq:operator}
F\left({\boldsymbol{\Omega}}^{(i)}, {\boldsymbol{\Gamma}}^{(i)}\right)&={\boldsymbol{G}}+\gamma{\boldsymbol{M}}\!\left({\boldsymbol{\Gamma}}^{(i)}\right)^\intercal {\boldsymbol{\Omega}}^{(i)}{\boldsymbol{M}}\!\left({\boldsymbol{\Gamma}}^{(i)}\right),\end{aligned}$$ i.e. $F\left({\boldsymbol{H}}^{(i)}, {\boldsymbol{L}}^{(i)}\right)={\boldsymbol{H}}^{(i+1)}$ according to .

In order to prove that ${\boldsymbol{H}}^{(i)}$ given according to Lemma \[lemma:VI\_matrix\_iteration\] is upper bounded, the following technical Lemma \[lemma:monotonically\_increasing\] is required first, which generalizes [@Landelius.1997 Lemma B.1.1] to cope with the reference-dependent Q-function. Note that knowledge of the exact structure of ${\boldsymbol{H}}$ and therefore the analytical solution by means of Theorem \[Theo:Q\] plays a crucial role for the extension to the tracking case.

\[lemma:monotonically\_increasing\] Let ${\boldsymbol{W}}^{(i)}$ be an arbitrary matrix sequence and ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(0)}\preceq{\boldsymbol{Z}}^{(0)}$. Then, given the sequences ${\boldsymbol{Z}}^{(i+1)}=F\left({\boldsymbol{Z}}^{(i)},{\boldsymbol{W}}^{(i)}\right)$ and ${\boldsymbol{H}}^{(i+1)}=F\left({\boldsymbol{H}}^{(i)},{\boldsymbol{L}}\left({\boldsymbol{H}}^{(i)}\right)\right)$ with it follows that ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i+1)}\preceq{\boldsymbol{Z}}^{(i+1)}$.

See \[append:monotonically\_increasing\].

In the next step, upper boundedness of ${\boldsymbol{H}}^{(i)}$ in the sense that ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq {\boldsymbol{Y}}$ is shown. For the regulation case, boundedness was shown in [@Landelius.1997 Lemma B.1.2]. In contrast to that, we consider the specific structure of the iteration in the tracking case (cf. ) and prove boundedness of ${\boldsymbol{H}}^{(i)}$ for the more generalized tracking formulation. For reasons of self-consistency, the complete proof is given.

\[lemma:upper\_bounded\] Let $\left({\boldsymbol{A}}, {\boldsymbol{B}}\right)$ be controllable and ${\boldsymbol{H}}^{(i)}$ be the sequence with ${\boldsymbol{H}}^{(0)}={\boldsymbol{0}}$. Then, there exists ${\boldsymbol{Y}}$ such that ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{Y}}$.

See \[append:upper\_bounded\].

We will now get to the main result of our convergence analysis and show that the proposed value iteration converges to the optimal weight vector ${\boldsymbol{w}}^*$ and the optimal control law ${\boldsymbol{L}}^*$ in the tracking case described by Problem \[problem1\].

Let ${\boldsymbol{R}}^\intercal = {\boldsymbol{R}} \succ {\boldsymbol{0}}$, ${\boldsymbol{Q}}^\intercal = {\boldsymbol{Q}} \succeq {\boldsymbol{0}}$, $({\boldsymbol{A}}, {\boldsymbol{B}})$ controllable and ${\boldsymbol{w}}^{(0)}={\boldsymbol{0}}$, i.e. ${\boldsymbol{H}}^{(0)} = {\boldsymbol{0}}$. Then, iterating between and yields ${\boldsymbol{H}}^{(i)}\rightarrow{\boldsymbol{H}}$, i.e. ${\boldsymbol{w}}^{(i)}\rightarrow{\boldsymbol{w}}^*$ as well as ${\boldsymbol{L}}^{(i)}\rightarrow{\boldsymbol{L}}^*$.

According to Lemma \[lemma:VI\_matrix\_iteration\], the value iteration is equivalent to iterating on ${\boldsymbol{H}}^{(i)}$ (cf. ). With ${\boldsymbol{Z}}^{(0)}={\boldsymbol{H}}^{(0)}$ and ${\boldsymbol{Z}}^{(i+1)}=F\left({\boldsymbol{Z}}^{(i)}, {\boldsymbol{L}}\left({\boldsymbol{H}}^{(i+1)}\right)\right)$, Lemma \[lemma:monotonically\_increasing\] yields ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{Z}}^{(i)}$. With ${\boldsymbol{H}}^{(0)}={\boldsymbol{0}}$ follows ${\boldsymbol{H}}^{(1)}={\boldsymbol{G}}\succeq{\boldsymbol{0}}$ and hence ${\boldsymbol{H}}^{(1)}-{\boldsymbol{Z}}^{(0)}\succeq{\boldsymbol{0}}$. The proof is drawn by induction: Assume the induction hypothesis ${\boldsymbol{H}}^{(i)}-{\boldsymbol{Z}}^{(i-1)}\succeq{\boldsymbol{0}}$. Then, $$\begin{aligned}
    &{\boldsymbol{H}}^{(i+1)}-{\boldsymbol{Z}}^{(i)}\nonumber \\ &=\gamma {\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{H}}^{(i)}\right)\right)^\intercal\left({\boldsymbol{H}}^{(i)}-{\boldsymbol{Z}}^{(i-1)}\right){\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{H}}^{(i)}\right)\right)\nonumber \\&\succeq{\boldsymbol{0}}.
    \end{aligned}$$ This implies ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{Z}}^{(i)}\preceq{\boldsymbol{H}}^{(i+1)}$.

As the matrix sequence is upper bounded by ${\boldsymbol{Y}}$ according to Lemma \[lemma:upper\_bounded\] and ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{H}}^{(i+1)}$, the limit ${\boldsymbol{H}}^{(\infty)}$ exists, i.e. the value iteration converges to $$\begin{aligned}
    {\boldsymbol{H}}^{(\infty)}={\boldsymbol{G}}+\gamma{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{H}}^{(\infty)}\right)\right)^\intercal{\boldsymbol{H}}^{(\infty)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{H}}^{(\infty)}\right)\right).
    \end{aligned}$$ Furthermore, ${\boldsymbol{L}}\left({\boldsymbol{H}}^{(\infty)}\right)$ minimizes $\hat{Q}_{0}^{(\infty)}$ according to Lemma \[lemma:hilfslemma\]. Thus, $$\begin{aligned}
    \lim_{i\rightarrow \infty}\epsilon_k^{(i)}&=c_k+\gamma{{\boldsymbol{z}}_{k+1}^{*}}^\intercal{\boldsymbol{H}}^{(\infty)}{\boldsymbol{z}}_{k+1}^{*}\nonumber \\ &-{\boldsymbol{z}}_k^\intercal \, F\left({\boldsymbol{H}}^{(\infty)}, {\boldsymbol{L}}\left({\boldsymbol{H}}^{(\infty)}\right)\right){\boldsymbol{z}}_k = 0.
    \end{aligned}$$ Therefore, ${\boldsymbol{H}}^{(i)}\rightarrow{\boldsymbol{H}}^{(\infty)}={\boldsymbol{H}}$, i.e. ${\boldsymbol{w}}^{(i)}\rightarrow{\boldsymbol{w}}^*$ and ${\boldsymbol{L}}^{(i)}\rightarrow{\boldsymbol{L}}^*$ which completes the proof.

Results {#simulationResults}
=======

\[sectionsimulation\] In order to validate our proposed method, we show results for two LQ-tracking problems with unknown system dynamics. We furthermore compare our results with an RL tracking method which assumes that the reference can be described by a time-invariant exo-system ${\boldsymbol{f}}_{\text{ref}}$.

For excitation purposes (cf. ), the variance of the input excitation ${\boldsymbol{\xi}}$ in is set to ${\boldsymbol{\Sigma}}=0.1$ and Gaussian noise $\xi_{\text{ref}}\sim \mathcal{N}\left({0},0.1\right)$ is added to the reference trajectory which we intend to track while learning is active. We furthermore choose the number of data tuples to $M=1.2L$, with $L$ according to , the stopping criterion to $e_{{\boldsymbol{w}}}=\SI{1e-6}{}$, the discount factor to $\gamma = 0.9$, the horizon on which the reference is known to $N=10$ and finally ${\boldsymbol{x}}_0= {\boldsymbol{0}}$ for all experiments.

Simulation Examples
-------------------

The first system is a second-order rotatory mass-spring-damper system, whereas the second system is a sixth-order linear single-track steering model. Both systems were initially given in continuous time and hence discretized by means of a Tustin approximation with a sampling time of $\SI{1}{\second}$.

### System 1

Consider a rotatory mass-spring-damper system modeled by the discrete-time second order linear state space representation $$\begin{aligned}
\label{sys_diskret}
{\boldsymbol{x}}_{k+1}={\begin{bmatrix*}[c]0.99&0.9\\-0.02&0.8
\end{bmatrix*}}{\boldsymbol{x}}_k+{\begin{bmatrix*}[c]0.01\\0.02
\end{bmatrix*}}u_k.\end{aligned}$$

The control input $u_k$ is a torque command applied to the system. Note that this system is not known to the controller and only needed for simulation as well as validation purposes. Furthermore, let $$\begin{aligned}
{\boldsymbol{Q}}=\text{diag}(100, 0) \text{ and } {\boldsymbol{R}}=1\end{aligned}$$ be the parameters of the cost function , i.e. we focus on tracking a reference angle $x_{1,\text{ref}}=\alpha_\text{ref}$ that is given on a moving horizon of length $N=10$.

### System 2

Our second example system is a sixth order linear single-track steering model $$\label{continuous_ss}
{\boldsymbol{x}}_{k+1} = {{\boldsymbol{A}}}{\boldsymbol{x}}_k + {{\boldsymbol{B}}}u_k$$ with $$\setlength{\arraycolsep}{3pt}
{{\boldsymbol{A}}} = \begin{bmatrix*}[r]
-0.741          &-0.033 &0      &0  &-\SI{2.0E-4}{} &-\SI{6.4E-3}{}\\
4.146   &-0.914 &0      &0  &\SI{4.7E-3}{}  &0.151\\
2.073                               &0.043              &1      &0  &\SI{2.4E-3}{}  &0.076\\
23.326                          &0.106                                          &20 &1  &0.022  &0.693\\
23.499      &-2.593     &0      &0  &-0.939 &-2.053\\
11.749                              &-1.297                                         &0      &0  &0.031  &-0.027
\end{bmatrix*}$$ and $$\begin{aligned}
\setlength{\arraycolsep}{3pt}
{{\boldsymbol{B}}} = \begin{bmatrix}
-\SI{3.8E-4}{} & 0.0091 & 0.0046 & 0.041 & 0.117  & 0.059
\end{bmatrix}^\intercal,\end{aligned}$$ where $u_k$ is a torque applied to the steering wheel by the controller. The parameters of the corresponding continuous-time model can be found in [@Kopf.2018]. The system state ${\boldsymbol{x}}_k$ is given by $$\begin{aligned}
 \label{eq:Zustaende Fzg}
{\boldsymbol{x}}_k = \begin{bmatrix}
\beta_k & {\psi}_{r,k} & \psi_k & y_k & {\delta}_{v,k} & \delta_k
\end{bmatrix}^\intercal,\end{aligned}$$ with sideslip angle $\beta_k$, yaw angle $\psi_k$, yaw rate ${\psi}_{r,k}$, lateral deviation from the origin $y_k$, steering wheel angle $\delta_k$ and angular velocity ${\delta}_{v,k}$. Geometric relations are depicted in Fig. \[fig:einspurmodell\], where $v$ is the constant velocity ($\SI{20}{\meter\per \second}$ in our example) and $\delta_s = 0.0625\delta$ denotes the steering angle in contrast to the steering wheel angle $\delta$.

As we desire to track the lateral position $x_{4,\text{ref}} = y_{\text{ref}}$, we choose $$\begin{aligned}
{\boldsymbol{Q}}=\text{diag}(0, 0, 0, 100, 0, 0) \text{ and } {\boldsymbol{R}}=1.\end{aligned}$$

Evaluation Method
-----------------

In order to compare our method with the class of RL tracking algorithms where the reference trajectory is assumed to be generated by a time-invariant exo-system ${\boldsymbol{f}}_{\text{ref}}({\boldsymbol{r}}_k)={\boldsymbol{F}}_{\text{ref}}{\boldsymbol{r}}_k$, let both the reference angle $\alpha_{\text{ref}}$ (system 1) and the reference lateral position $y_{\text{ref}}$ (system 2) be generated by $$\begin{aligned}
\label{eq:refe}
{\boldsymbol{r}}_{k+1} &= \underbrace{{\begin{bmatrix*}[r]0.9801&0.1987\\-0.1987&0.9801
\end{bmatrix*}}}_{{\boldsymbol{F}}_{\text{ref}}}{\boldsymbol{r}}_k, \\
\alpha_{\text{ref}} &= y_{\text{ref}} = {\begin{bmatrix*}[r]1&0
\end{bmatrix*}}{\boldsymbol{r}}_k\label{eq:refe1}\end{aligned}$$ during the training procedure. For comparison reasons, we then train both our proposed method and an algorithm as in [@Luo.2016; @Kiumarsi.2014] which assumes that the reference always follows the dynamics of an exo-system ${\boldsymbol{f}}_{\text{ref}}({\boldsymbol{r}}_k)$ (in the following termed *comparison algorithm*). The comparison algorithm is trained on the augmented system $$\begin{aligned}
\label{eq:augmented_sys}
{\boldsymbol{x}}_{\text{aug},k+1}={\begin{bmatrix*}[c]{\boldsymbol{A}}&{\boldsymbol{0}}\\{\boldsymbol{0}}&{\boldsymbol{F}}_{\text{ref}}
\end{bmatrix*}}{\boldsymbol{x}}_{\text{aug},k}+{\begin{bmatrix*}[c]{\boldsymbol{B}}\\ {\boldsymbol{0}}
\end{bmatrix*}}{\boldsymbol{u}}_k,\end{aligned}$$ where ${\boldsymbol{x}}_{\text{aug},k}={\begin{bmatrix*}[r]{\boldsymbol{x}}_k^\intercal&{\boldsymbol{r}}_k^\intercal
\end{bmatrix*}}^\intercal$. After the training, we vary $\alpha_{\text{ref}} = y_{\text{ref}}$ (arbitrary references such as different frequencies, ramps and steps) in order to show that the controller learned with our method successfully generalizes to these references. In contrast, we show the resulting behavior of the comparison algorithm which is constructed on the assumption that the reference dynamics always follows .

Our evaluation is twofold. On one hand, after the learning process has finished, we analyze the root mean square (RMS) tracking errors $\alpha_{\text{RMS}}$ and $y_{\text{RMS}}$ between the learned tracking behavior by means of the trajectory $\alpha_{\text{learned}}$ and $y_{\text{learned}}$ where the system dynamics is unknown (both for our algorithm and the comparison algorithm) and the optimal solution $\alpha_{\text{opt}}$ and $y_{\text{opt}}$ which results from Theorem \[Theo:Q\] and known system dynamics. On the other hand, we compare the learned weights ${\boldsymbol{w}}$ of our algorithm with the optimal solution ${\boldsymbol{w}}^*$ (weights corresponding to Theorem \[Theo:Q\] and Lemma \[lemma:Qexakt\]). In order to achieve comparability for different ranges of ${\boldsymbol{w}}$, we normalize the absolute error of each weight with the maximum absolute weight $\max_j {\left\lvert \left\{{\boldsymbol{w}}^*\right\}_j \right\rvert}$ and define the average of this normalized absolute error by $$\begin{aligned}
\label{eq:e1}
e_{\text{I}}&=\frac{1}{L}\sum_{i=1}^{L}\frac{{\left\lvert \left\{{\boldsymbol{w}}\right\}_i-\left\{{\boldsymbol{w}}^*\right\}_i \right\rvert}}{\max_j {\left\lvert \left\{{\boldsymbol{w}}^*\right\}_j \right\rvert}}.\end{aligned}$$ Its maximum is given by $$\begin{aligned}
\label{eq:e2}
e_{\text{II}}&=\max_i \frac{{\left\lvert \left\{{\boldsymbol{w}}\right\}_i-\left\{{\boldsymbol{w}}^*\right\}_i \right\rvert}}{\max_j {\left\lvert \left\{{\boldsymbol{w}}^*\right\}_j \right\rvert}}.\end{aligned}$$ Note that these measures are only reasonable to judge how well our method has learned the unknown optimal weights ${\boldsymbol{w}}^*$ as the comparison algorithm does not learn ${\boldsymbol{w}}^*$ corresponding to Theorem \[Theo:Q\] and Lemma \[lemma:Qexakt\] but the optimal weights corresponding to .

= 0.2mm

[lcc]{}

&\
&&\
&proposed method&comparison alg. [@Luo.2016; @Kiumarsi.2014]\
$\alpha_{\text{RMS}}\!$&$\SI{2.1E-3}{}$&$1.55$\
$e_{\text{I}}$&$\SI{6.2E-5}{}$&–\
$e_{\text{II}}$&$\SI{2.1E-3}{}$&–\
\

&\
&&\
&proposed method&comparison alg. [@Luo.2016; @Kiumarsi.2014]\
$y_{\text{RMS}}\!$&$\SI{2.9E-5}{}$&$1.01$\
$e_{\text{I}}$&$\SI{2.1E-7}{}$&–\
$e_{\text{II}}$&$\SI{1.4E-5}{}$&–\

Results of the Adaptive Tracking Controller
-------------------------------------------

The RMS tracking errors $\alpha_{\text{RMS}}$ for system 1, $y_{\text{RMS}}$ for system 2 and the corresponding weight estimation errors $e_{\text{I}}$ and $e_{\text{II}}$ are given in Table \[table1\].

Plots of the corresponding tracking performances of our proposed method and the comparison method are given in Fig. \[fig:sys1\_x\] for system 1 and Fig. \[fig:sys2\_x\] for system 2. Here, the vertical dash-dotted lines indicate the weight update of our method (i.e. when $k=M$). The reference trajectories $\alpha_{\text{ref}}$ and $y_{\text{ref}}$ are depicted in gray. The black dashed lines show the optimal solutions $\alpha_{\text{opt}}$ and $y_{\text{opt}}$ calculated using full system knowledge and the red line the learned behavior $\alpha_{\text{learned, our}}$ and $y_{\text{learned, our}}$ without knowledge of the system matrices ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ when using our proposed method. The resulting tracking behavior $\alpha_{\text{learned, comparison}}$ and $y_{\text{learned, comparison}}$ of the comparison method [@Luo.2016; @Kiumarsi.2014] is depicted in blue.

The weight estimation measures $e_{\text{I}}$ and $e_{\text{II}}$ in Table \[table1\] indicate that for both systems the optimal Q-function weights have successfully been learned. The decay of $e_{\text{I}}$ and $e_{\text{II}}$ during the value iteration (with iteration index $i$ in Algorithm \[algorithm\]) is depicted in Fig. \[fig:sys1\_e\] (system 1) and Fig. \[fig:sys2\_e\] (system 2).

Discussion {#sec:discussion}
----------

Under the excitation condition , our learning controller converges to the optimal control law according to Section \[convergence\]. Comparing $\alpha_{\text{RMS}}$ (system 1), $y_{\text{RMS}}$ (system 2) and considering Fig. \[fig:sys1\_x\], Fig. \[fig:sys2\_x\] and Fig. \[fig:sys1\_x\_detail\] it is obvious that our algorithm is successfully tracking arbitrary references that have not been seen during training. State of the art methods which assume that the reference follows a time-invariant exo-system ${\boldsymbol{f}}_{\text{ref}}$ (e.g. [@Luo.2016; @Kiumarsi.2014]) are very effective as long as this assumption holds but their tracking performance decreases as soon as the reference to track deviates from the sine described by –. This behavior is hardly surprising as these controllers were specifically designed under the assumption of time-invariant exo-system dynamics.

Due to the explicit dependency of our Q-function on the reference on a moving horizon $N$, the learned weights generalize to references that are unknown during the learning procedure such as the ramps, steps and curves in the examples. Due to the optimal behavior according to our proposed controller exhibits predictive rather than only reactive behavior as can be seen in Fig. \[fig:sys1\_x\_detail\] which depicts a more detailed view of the step in Fig.\[fig:sys1\_x\].

We would further like to point out that exact knowledge of the structure of ${\boldsymbol{H}}_K$ and thus ${\boldsymbol{H}}$ which results from Theorem \[Theo:Q\] is very beneficial if not vital for the Q-learning method to work efficiently as it helps to reduce the number of weights that have to be learned. If one would only assume ${\boldsymbol{H}}$ to be quadratic and symmetric, $L=325$ (system 1) and $L=2701$ (system 2) weights would have to be estimated. Considering Lemma \[lemma:Qexakt\_long\], these numbers reduce to $L=247$ (system 1) and $L=2011$ (system 2) and exploiting the sparsity properties of ${\boldsymbol{Q}}$ according to Lemma \[lemma:Qexakt\], $L=84$ (system 1) and $L=146$ (system 2) result which renders the complexity tractable.

The choice of the moving horizon length $N$ depends on the available information regarding the reference to track. Furthermore, the larger $N$, the more predictive the learned controller will be but the more unknown weights ${\boldsymbol{w}}$ need to be learned (see ). Thus, an appropriate choice of $N$ obviously depends on the specific application.

Conclusion
==========

In this paper, we proposed a new Reinforcement-Learning-based algorithm that is able to track an arbitrary reference trajectory which is given on a moving horizon while the system dynamics is unknown. In contrast to state-of-the-art methods that are based on RL respectively Adaptive Dynamic Programming, our method does not require the reference trajectory to be generated by an exo-system. It explicitly incorporates arbitrary reference values in a new Q-function. This Q-function, which is constructed such that its minimizing control is part of the solution of the optimal LQ tracking problem, generalizes to arbitrary reference trajectories given on a moving horizon.

We showed that the analytical solution to this Q-function has a specific structure w.r.t. the current state and control as well as the reference on the given horizon. Based thereon, sparsity properties of the resulting structure were exploited in order to reduce the Q-function weights that have to be estimated. The temporal difference error of the reference-dependent Q-function serves as a target in order to learn the optimal tracking behavior online when the system dynamics is unknown. Here, the choice of basis functions is based on the findings regarding the specific structure of the analytical solution. We proved that this iterative algorithm converges to the optimal solution.

Proof of Theorem \[Theo:Q\] {#append:notations}
===========================

Theorem \[Theo:Q\] is vital for an efficient Q-function parametrization as it yields the exact structure of the analytical solution of the novel reference-dependent Q-function which reduces the numbers of function approximation weights needed (see Lemma \[lemma:Qexakt\_long\] and Lemma \[lemma:Qexakt\] and the discussion regarding the number of weights to be learned in Section \[sec:discussion\]).

Although the notation is complex, it is required to keep the technical proof correct. This is due to the fact that clear distinction of the current time step $k$ and the time step $\kappa$ in the current optimization horizon $K$ starting at $k$ is required. The main idea of this proof is to use dynamic programming and prove the analytical solution of the Q-function by means of backwards induction. Writing down this dynamic programming procedure is cumbersome but recompensates by means of the exact structure of our new Q function. For $l=1,\dots,d$, the $l$-th submatrix of a matrix ${\boldsymbol{\Pi}}\in\mathbb{R}^{p\times nd}$ is defined as $${\boldsymbol{\Pi}}\left[l\right]={\begin{bmatrix*}[c]{\boldsymbol{\Pi}}(1,(l-1)n+1) & \cdots & {\boldsymbol{\Pi}}(1,nl)\\
    \vdots & \ddots & \vdots\\
    {\boldsymbol{\Pi}}(p,(l-1)n+1) & \cdots & {\boldsymbol{\Pi}}(p,nl)
\end{bmatrix*}}.$$

Furthermore, note that ${\boldsymbol{I}}_n$ denotes the $n\times n$ identity matrix and $\zeta\in\mathbb{N}_0$ is a placeholder. Later, $\zeta$ will be replaced by $\eta$ (respectively $\eta+1$ in the inductive step), where $\eta=K-\kappa$ denotes the remaining time steps on the horizon $K$. Further, $p$ is an index with $p\in\mathbb{N}:$ $p>1$. We define the following shorthand notations.

$$\begin{aligned}
&\begin{aligned}
{\boldsymbol{X}}_{\zeta}^{0} =& {\boldsymbol{X}}^{0}={\begin{bmatrix*}[r] {\boldsymbol{I}}_n & -{\boldsymbol{I}}_n 
\end{bmatrix*}},\\
{\boldsymbol{X}}_{\zeta}^{1} =& \sqrt{\gamma}\left(-{\boldsymbol{X}}^{0}[1]{\boldsymbol{B}}{\boldsymbol{G}}_\zeta{\begin{bmatrix*}[r]{\boldsymbol{F}}_\zeta & {\boldsymbol{L}}_\zeta^{\zeta-1} & \cdots & {\boldsymbol{L}}_\zeta^{0}
\end{bmatrix*}} \right.\\&+\left. {\begin{bmatrix*}[r]{\boldsymbol{X}}^{0}[1]{\boldsymbol{A}} & {\boldsymbol{X}}^{0}[2] & {\boldsymbol{0}} & \cdots & {\boldsymbol{0}}
\end{bmatrix*}}\right)\\
\intertext{and}
{\boldsymbol{X}}_{\zeta}^{p} =&         \sqrt{\gamma}\left(-{\boldsymbol{X}}_{\zeta}^{p-1}[1]{\boldsymbol{B}}{\boldsymbol{G}}_\zeta{\begin{bmatrix*}[r]{\boldsymbol{F}}_\zeta & {\boldsymbol{L}}_\zeta^{\zeta-1} & \cdots & {\boldsymbol{L}}_\zeta^{0}
\end{bmatrix*}}\right.\nonumber
\end{aligned}
\\
&\begin{aligned}
+\left. {\begin{bmatrix*}[r]{\boldsymbol{X}}_{\zeta}^{p-1}[1]{\boldsymbol{A}} & {\boldsymbol{0}} & {\boldsymbol{X}}_{\zeta}^{p-1}[2] & \cdots & {\boldsymbol{X}}_{\zeta}^{p-1}[\zeta-1]
\end{bmatrix*}} \right)\!,
\end{aligned}\end{aligned}$$ as well as $$\begin{aligned}
&\begin{aligned}
{\boldsymbol{U}}_{\zeta}^{1} =& -{\boldsymbol{G}}_\zeta{\begin{bmatrix*}[r]{\boldsymbol{F}}_\zeta & {\boldsymbol{L}}_\zeta^{\zeta-1} & \cdots & {\boldsymbol{L}}_\zeta^{0} 
\end{bmatrix*}}\\
\intertext{and}
{\boldsymbol{U}}_{\zeta}^{p} =& \sqrt{\gamma}\left(-{\boldsymbol{U}}_{\zeta}^{p-1}[1]{\boldsymbol{B}}{\boldsymbol{G}}_\zeta{\begin{bmatrix*}[r]{\boldsymbol{F}}_\zeta & {\boldsymbol{L}}_\zeta^{\zeta-1} & \cdots & {\boldsymbol{L}}_\zeta^{0}
\end{bmatrix*}}\right.\nonumber
\end{aligned}
\\
&\begin{aligned}
+\left.{\begin{bmatrix*}[r]{\boldsymbol{U}}_{\zeta}^{p-1}[1]{\boldsymbol{A}} & {\boldsymbol{0}} & {\boldsymbol{U}}_{\zeta}^{p-1}[2] & \cdots & {\boldsymbol{U}}_{\zeta}^{p-1}[\zeta-1]
\end{bmatrix*}}\right)\!,
\end{aligned}\end{aligned}$$ with $$\begin{aligned}
{\boldsymbol{M}}_\zeta &= \gamma{\boldsymbol{B}}^\intercal\left(\vphantom{\sum\limits_{i=2}^{\zeta}\left({\boldsymbol{X}}_{K-\zeta}^{i}[1]\right)^\intercal {\boldsymbol{Q}} {\boldsymbol{X}}_{K-\zeta}^{i}[1]}\right.\sum\limits_{i=0}^{\zeta-2}\left({\boldsymbol{X}}_{\zeta}^{i}[1]\right)^\intercal {\boldsymbol{Q}} {\boldsymbol{X}}_{\zeta}^{i}[1]\nonumber& \\   &\hphantom{=\gamma{\boldsymbol{B}}^\intercal\left(\vphantom{\sum\limits_{i=2}^{\zeta-1}{\boldsymbol{X}}[1]_{K-\zeta}^{i^\intercal} {\boldsymbol{Q}} {\boldsymbol{X}}[1]_{K-\zeta}^{i}}\right.}+\left.\sum\limits_{i=1}^{\zeta-2}\left({\boldsymbol{U}}_{\zeta}^{i}[1]\right)^\intercal {\boldsymbol{R}} {\boldsymbol{U}}_{\zeta}^{i}[1]\right), &\\
{\boldsymbol{F}}_\zeta &= {\boldsymbol{M}}_\zeta{\boldsymbol{A}},&\\
{\boldsymbol{G}}_\zeta^{-1} &= {\boldsymbol{M}}_\zeta{\boldsymbol{B}}+{\boldsymbol{R}}, &\\
{\boldsymbol{L}}_\zeta^j &= \begin{cases}\gamma{\boldsymbol{B}}^\intercal\!\big({\boldsymbol{X}}^{0}[1]\big)^\intercal {\boldsymbol{Q}} {\boldsymbol{X}}^{0}[2], \text{ for }j=\zeta-1,&\vspace{0.2cm}\\
\gamma{\boldsymbol{B}}^\intercal\!\left(\sum\limits_{i=1}^{\zeta-2}\left({\boldsymbol{X}}_{\zeta}^{i}[1]\right)^\intercal {\boldsymbol{Q}} {\boldsymbol{X}}_{\zeta}^{i}\left[\zeta-j\right] \right. \\
\hphantom{\gamma{\boldsymbol{B}}^\intercal\!\left(\vphantom{\sum\limits_{i=2}^{\zeta-1}{\boldsymbol{U}}[1]_{K-\zeta}^{i^\intercal} {\boldsymbol{R}}}\right.}  +\left. \sum\limits_{i=1}^{\zeta-2}\left({\boldsymbol{U}}_{\zeta}^{i}[1]\right)^\intercal {\boldsymbol{R}} {\boldsymbol{U}}_{\zeta}^{i}\left[\zeta-j\right]\right)\!,&\\ \hphantom{\gamma{\boldsymbol{B}}^\intercal{\boldsymbol{X}}[1]_{K-\zeta}^{j+1^\intercal} {\boldsymbol{Q}} {\boldsymbol{X}}[2]_{K-\zeta}^{j+1}, \,\quad}\text{for }j<\zeta-1,&\end{cases}\end{aligned}$$ with $j\in\mathbb{N}_0$.

Let furthermore $$\begin{aligned}
{\boldsymbol{\rho}}_0^\kappa &= {\begin{bmatrix*}[r]{\boldsymbol{x}}_{k_\kappa}^\intercal & {\boldsymbol{r}}_{k_\kappa}^\intercal
\end{bmatrix*}}\left({\boldsymbol{X}}^{0}\right)^\intercal, \label{eq_rho0}\\
{\boldsymbol{\rho}}_1^\kappa &= {\begin{bmatrix*}[r]{\boldsymbol{x}}_{k_\kappa}^\intercal & {\boldsymbol{u}}_{k_\kappa}^\intercal & {\boldsymbol{r}}_{k_\kappa+1}^\intercal
\end{bmatrix*}} {\begin{bmatrix*}[c]\left({\boldsymbol{X}}^{0}[1]{\boldsymbol{A}}\right)^\intercal  \\ \left({\boldsymbol{X}}^{0}[1]{\boldsymbol{B}}\right)^\intercal \\\left({\boldsymbol{X}}^{0}[2]\right)^\intercal
\end{bmatrix*}}\!, \\
{\boldsymbol{\mu}}_i^\kappa &= {\begin{bmatrix*}[r]{\boldsymbol{x}}_{k_\kappa}^\intercal & {\boldsymbol{u}}_{k_\kappa}^\intercal & {\boldsymbol{r}}_{k_\kappa+2}^\intercal & \cdots & {\boldsymbol{r}}_{k+K}^\intercal
\end{bmatrix*}}{\begin{bmatrix*}[c]\left({\boldsymbol{U}}_{\eta}^{\eta-i}[1]{\boldsymbol{A}}\right)^\intercal  \\ \left({\boldsymbol{U}}_{\eta}^{\eta-i}[1]{\boldsymbol{B}}\right)^\intercal \\\left({\boldsymbol{U}}_{\eta}^{\eta-i}[2]\right)^\intercal \\ \vdots \\ \left({\boldsymbol{U}}_{\eta}^{\eta-i}[\eta]\right)^\intercal
\end{bmatrix*}}\!, \\
{\boldsymbol{\chi}}_i^\kappa &= {\begin{bmatrix*}[r]{\boldsymbol{x}}_{k_\kappa}^\intercal & {\boldsymbol{u}}_{k_\kappa}^\intercal & {\boldsymbol{r}}_{k_\kappa+2}^\intercal & \cdots & {\boldsymbol{r}}_{k+K}^\intercal
\end{bmatrix*}}{\begin{bmatrix*}[c]\left({\boldsymbol{X}}_{\eta}^{\eta-i}[1]{\boldsymbol{A}}\right)^\intercal  \\ \left({\boldsymbol{X}}_{\eta}^{\eta-i}[1]{\boldsymbol{B}}\right)^\intercal \\\left({\boldsymbol{X}}_{\eta}^{\eta-i}[2]\right)^\intercal \\ \vdots \\ \left({\boldsymbol{X}}_{\eta}^{\eta-i}[\eta]\right)^\intercal
\end{bmatrix*}}\label{eq_xi}\!,\end{aligned}$$ where ${k_\kappa}=k+K-\eta=k+\kappa$ and $i\in\mathbb{N}$.

In the first step, we use backwards induction to prove that the Q-function $\tensor*[^K]{\!Q}{_{\kappa}}$ (cf. Definition \[def:Q\]) for system with the objective function is given by $$\begin{aligned}
    \tensor*[^K]{\!Q}{_{\kappa}} &= \frac{1}{2}\!\left(\vphantom{\sum_{i=1}^{N-2}\left({\boldsymbol{\chi}}_i{\boldsymbol{Q}}{\boldsymbol{\chi}}_i^\intercal+{\boldsymbol{\mu}}_i{\boldsymbol{R}}{\boldsymbol{\mu}}_i^\intercal\right)}
    {\boldsymbol{\rho}}_0^\kappa{\boldsymbol{Q}}\left({\boldsymbol{\rho}}_0^{\kappa}\right)^\intercal+{\boldsymbol{u}}_{k+\kappa}^\intercal {\boldsymbol{R}}{\boldsymbol{u}}_{k+\kappa}+\gamma{\boldsymbol{\rho}}_1^\kappa{\boldsymbol{Q}}\left({\boldsymbol{\rho}}_1^{\kappa}\right)^\intercal\right.\nonumber\\ &\hphantom{=\frac{1}{2}} \left.+ \gamma\sum_{i=1}^{K-\kappa-2}\left({\boldsymbol{\chi}}_i^\kappa{\boldsymbol{Q}}\left({\boldsymbol{\chi}}_i^{\kappa}\right)^\intercal+{\boldsymbol{\mu}}_i^\kappa{\boldsymbol{R}}\left({\boldsymbol{\mu}}_i^{\kappa}\right)^\intercal\right)\right).\label{eq:Q_kappa}
    \end{aligned}$$ Starting from $\tensor*[^K]{\!Q}{_{K}}$ (cf. ), ${\boldsymbol{u}}_{k+K}^*={\boldsymbol{0}}$ directly follows from Definition \[def:Q\] and $$\begin{aligned}
    \left.\frac{\partial \tensor*[^K]{\!Q}{_{K}}}{\partial {\boldsymbol{u}}_{k+K}}\right|_{{\boldsymbol{u}}_{k+K}^{*}}={\boldsymbol{0}}, \qquad \left.\frac{\partial^2 \tensor*[^K]{\!Q}{_{K}}}{\partial {\boldsymbol{u}}_{k+K}^2}\right|_{{\boldsymbol{u}}_{k+K}^{*}}={\boldsymbol{R}}\succ{\boldsymbol{0}}.
    \end{aligned}$$ Then, by iterating backwards in time, applying and the system dynamics , with $\eta = K-\kappa$, can be shown to hold for $\eta = 0, 1, 2$, i.e. $\kappa = K, K-1, K-2$. Furthermore, $$\begin{aligned}
    {\boldsymbol{u}}_{k+\kappa}^*=-{\boldsymbol{G}}_\eta\left({\boldsymbol{F}}_\eta{\boldsymbol{x}}_{k+\kappa}+\sum_{j=0}^{\eta-1}{\boldsymbol{L}}_\eta^j{\boldsymbol{r}}_{k+K-j}\right)\label{eq:uopt_G}
    \end{aligned}$$ minimizes because $$\begin{aligned}
\label{dQdu_d2Qdu2}
    \left.\frac{\partial \tensor*[^K]{\!Q}{_{\kappa}}}{\partial {\boldsymbol{u}}_{k+\kappa}}\right|_{{\boldsymbol{u}}_{k+\kappa}^{*}}={\boldsymbol{0}},\quad \text{where }
    \left.\frac{\partial^2 \tensor*[^K]{\!Q}{_{\kappa}}}{\partial {\boldsymbol{u}}_{k+\kappa}^2}\right|_{{\boldsymbol{u}}_{k+\kappa}^{*}}\succ{\boldsymbol{0}}
    \end{aligned}$$ is guaranteed as ${\boldsymbol{R}}\succ{\boldsymbol{0}}$ and ${\boldsymbol{Q}}\succeq{\boldsymbol{0}}$. The induction hypothesis $\tensor*[^K]{\!Q}{_{\kappa-1}}$ (see with $\kappa \rightarrow \kappa-1$) is then proven in the inductive step. This is done by representing $\tensor*[^K]{\!Q}{_{\kappa-1}}$ by means of and utilizing ${\boldsymbol{u}}_{k+\kappa}^{*}$ from . This yields $$\begin{aligned}
\label{eq:kappa-1induction}
    \tensor*[^K]{\!Q}{_{\kappa-1}}&=\frac{1}{2}{\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa-1}\\ {\boldsymbol{r}}_{k+\kappa-1}
\end{bmatrix*}}^\intercal\left({\boldsymbol{X}}^{0}\right)^\intercal{\boldsymbol{Q}}{\boldsymbol{X}}^{0}{\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa-1}\\ {\boldsymbol{r}}_{k+\kappa-1}
\end{bmatrix*}}\nonumber \\
    &\hphantom{ = }+\frac{1}{2}{\boldsymbol{u}}_{k+\kappa-1}^\intercal{\boldsymbol{R}}{\boldsymbol{u}}_{k+\kappa-1}\nonumber \\
    &\hphantom{ = }+\frac{1}{2}{\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa}\\ {\boldsymbol{r}}_{k+\kappa}
\end{bmatrix*}}^\intercal\left({\boldsymbol{X}}^{0}\right)^\intercal{\boldsymbol{Q}}{\boldsymbol{X}}^{0}{\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa}\\ {\boldsymbol{r}}_{k+\kappa}
\end{bmatrix*}}\nonumber \\
    &\hphantom{ = }+\frac{1}{2}\gamma\bar{{\boldsymbol{z}}}_k^\intercal \sum_{i=1}^{\eta-1}\left(\left({\boldsymbol{X}}_{\eta+1}^{i}\right)^\intercal{\boldsymbol{Q}}{\boldsymbol{X}}_{\eta+1}^{i}\right.\nonumber \\ &\left.\hphantom{ = +\frac{1}{2}\gamma\bar{{\boldsymbol{z}}}_k^\intercal \sum_{i=1}^{\eta-1}}+\left({\boldsymbol{U}}_{\eta+1}^{i}\right)^\intercal{\boldsymbol{R}}{\boldsymbol{U}}_{\eta+1}^{i}\right) \bar{{\boldsymbol{z}}}_k,
    \end{aligned}$$ where $\bar{{\boldsymbol{z}}}_k^\intercal={\begin{bmatrix*}[r]{\boldsymbol{x}}_{k+\kappa}^\intercal&{\boldsymbol{r}}_{k+\kappa+1}^\intercal&\dots&{\boldsymbol{r}}_{k+K}^\intercal
\end{bmatrix*}}$.

Then, replace ${\boldsymbol{x}}_{k+\kappa}={\boldsymbol{A}}{\boldsymbol{x}}_{k+\kappa-1}+{\boldsymbol{B}}{\boldsymbol{u}}_{k+\kappa-1}$ (cf. ) in which results in $$\begin{aligned}
    \tensor*[^K]{\!Q}{_{{\kappa-1}}} &= \frac{1}{2}\!\left(\vphantom{\sum_{i=1}^{N-2}\left({\boldsymbol{\chi}}_i{\boldsymbol{Q}}{\boldsymbol{\chi}}_i^\intercal+{\boldsymbol{\mu}}_i{\boldsymbol{R}}{\boldsymbol{\mu}}_i^\intercal\right)}
    {\boldsymbol{\rho}}_0^{\kappa-1}{\boldsymbol{Q}}\left({\boldsymbol{\rho}}_0^{{\kappa-1}}\right)^\intercal+{\boldsymbol{u}}_{k+{\kappa-1}}^\intercal {\boldsymbol{R}}{\boldsymbol{u}}_{k+{\kappa-1}}\right.\nonumber \\&\hphantom{=\frac{1}{2}X}+\gamma{\boldsymbol{\rho}}_1^{\kappa-1}{\boldsymbol{Q}}\left({\boldsymbol{\rho}}_1^{{\kappa-1}}\right)^\intercal\nonumber\\ &\hphantom{=\frac{1}{2}X} + \gamma\sum_{i=1}^{K-(\kappa-1)-2}\left({\boldsymbol{\chi}}_i^{\kappa-1}{\boldsymbol{Q}}\left({\boldsymbol{\chi}}_i^{{\kappa-1}}\right)^\intercal\nonumber\right. \\ &\hphantom{=\frac{1}{2}X}\left.\left.+{\boldsymbol{\mu}}_i^{\kappa-1}{\boldsymbol{R}}\left({\boldsymbol{\mu}}_i^{{\kappa-1}}\right)^\intercal\right)\vphantom{\sum_{i=1}^{N-2}\left({\boldsymbol{\chi}}_i{\boldsymbol{Q}}{\boldsymbol{\chi}}_i^\intercal+{\boldsymbol{\mu}}_i{\boldsymbol{R}}{\boldsymbol{\mu}}_i^\intercal\right)}\right)
    \end{aligned}$$ and yields the induction hypothesis ( with $\kappa\rightarrow\kappa-1$) and thus proves .

Thus, the analytical solution of [$\tensor*[^K]{\!Q}{_{0}}$]{} is quadratic w.r.t. ${\boldsymbol{\rho}}_0^\kappa$, ${\boldsymbol{u}}_{k}$, ${\boldsymbol{\rho}}_1^\kappa$, ${\boldsymbol{\chi}}_i^\kappa$ and ${\boldsymbol{\mu}}_i^\kappa$. As each of these components is linear w.r.t. ${\boldsymbol{x}}_k$, ${\boldsymbol{u}}_k$ and ${\boldsymbol{r}}_{k+1},\dots,{\boldsymbol{r}}_{k+N}$ according to –, Theorem \[Theo:Q\] follows directly for $\kappa=0$ and $K\geq N$.

Proof of Lemma \[lemma:VI\_matrix\_iteration\] {#append:VI_matrix_iteration}
==============================================

Let $v(\cdot)$ be a function that transforms a symmetrical squared matrix to a vector such that $v\!\left({\boldsymbol{H}}^{(i)}\right)={\boldsymbol{w}}^{(i)}$. With $$\begin{aligned}
    &c_k+\gamma{{\boldsymbol{w}}^{(i)}}^\intercal{\boldsymbol{\phi}}\left({\boldsymbol{z}}_{k+1}^{*(i)}\right)\nonumber\\&=\frac{1}{2}{\boldsymbol{z}}_k^\intercal\underbrace{\left({\boldsymbol{G}}+\gamma{\boldsymbol{M}}\!\left({\boldsymbol{L}}^{(i)}\right)^\intercal {\boldsymbol{H}}^{(i)}{\boldsymbol{M}}\!\left({\boldsymbol{L}}^{(i)}\right)\right)}_{={\boldsymbol{H}}^{(i+1)}}{\boldsymbol{z}}_k,
    \end{aligned}$$ , and the definition of ${\boldsymbol{\phi}}$ according to Lemma \[lemma:Qexakt\] follows that can be written as $$\begin{aligned}
    {\boldsymbol{w}}^{(i+1)}&=\left({{\boldsymbol{\Phi}}^\intercal}{\boldsymbol{\Phi}}\right)^{-1}{{\boldsymbol{\Phi}}^\intercal}{\begin{bmatrix*}[c]\frac{1}{2}{\boldsymbol{z}}_{k-M+1}^\intercal{\boldsymbol{H}}^{(i+1)}{\boldsymbol{z}}_{k-M+1}\\\vdots\\\frac{1}{2}{\boldsymbol{z}}_k^\intercal{\boldsymbol{H}}^{(i+1)}{\boldsymbol{z}}_k
\end{bmatrix*}}\nonumber \\ &=\underbrace{\left({{\boldsymbol{\Phi}}^\intercal}{\boldsymbol{\Phi}}\right)^{-1}\left({{\boldsymbol{\Phi}}^\intercal}{\boldsymbol{\Phi}}\right)}_{={\boldsymbol{I}}}v\left({\boldsymbol{H}}^{(i+1)}\right).\label{eq:lemma_iteration_I}
    \end{aligned}$$ Thus, as ${\boldsymbol{H}}^{(i+1)}$ is symmetrically constructed from ${\boldsymbol{w}}^{(i+1)}$, it follows from that iterating on ${\boldsymbol{w}}^{(i)}$ by means of the value iteration is equivalent to iterating on .

Proof of Lemma \[lemma:monotonically\_increasing\] {#append:monotonically_increasing}
==================================================

According to Lemma \[lemma:hilfslemma\], the control law ${\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)$ minimizes ${\boldsymbol{z}}_{k+1}^\intercal{\boldsymbol{Z}}^{(i)}{\boldsymbol{z}}_{k+1}$, $\forall i>0$. Thus, $$\begin{aligned}
    {\boldsymbol{z}}_k^\intercal{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{W}}^{(i)}\right)\right)^\intercal{\boldsymbol{Z}}^{(i)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{W}}^{(i)}\right)\right){\boldsymbol{z}}_k\nonumber \\
    \geq {\boldsymbol{z}}_k^\intercal{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)^\intercal{\boldsymbol{Z}}^{(i)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right){\boldsymbol{z}}_k
    \end{aligned}$$ follows. This yields $$\begin{aligned}
    {\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{W}}^{(i)}\right)\right)^\intercal{\boldsymbol{Z}}^{(i)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{W}}^{(i)}\right)\right)\nonumber \\-{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)^\intercal{\boldsymbol{Z}}^{(i)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)\succeq {\boldsymbol{0}}
    \end{aligned}$$ and hence $$\begin{aligned}
\label{eq:Z_Zdach}
    {\boldsymbol{Z}}^{(i+1)}=F\left({\boldsymbol{Z}}^{(i)},{\boldsymbol{W}}^{(i)}\right)\succeq F\left({\boldsymbol{Z}}^{(i)},{\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right) \eqqcolon \hat{{\boldsymbol{Z}}}^{(i+1)}.
    \end{aligned}$$ This also implies $$\begin{aligned}
\label{eq:this_implies}
    F\left({\boldsymbol{H}}^{(i)},{\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)\succeq F\left({\boldsymbol{H}}^{(i)},{\boldsymbol{L}}\left({\boldsymbol{H}}^{(i)}\right)\right) = {\boldsymbol{H}}^{(i+1)}.
    \end{aligned}$$ With ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(0)}\preceq{\boldsymbol{Z}}^{(0)}$ and the induction hypothesis ${\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{Z}}^{(i)}$, $$\begin{aligned}
    &F\left({\boldsymbol{H}}^{(i)},{\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)\nonumber \\
    &={\boldsymbol{G}}+\gamma{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)^\intercal{\boldsymbol{H}}^{(i)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)\nonumber \\
    &\preceq{\boldsymbol{G}}+\gamma{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)^\intercal{\boldsymbol{Z}}^{(i)}{\boldsymbol{M}}\left({\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)\nonumber \\
    &=F\left({\boldsymbol{Z}}^{(i)},{\boldsymbol{L}}\left({\boldsymbol{Z}}^{(i)}\right)\right)=\hat{{\boldsymbol{Z}}}^{(i+1)}
    \end{aligned}$$ follows. With , this yields ${\boldsymbol{H}}^{(i+1)}\preceq\hat{{\boldsymbol{Z}}}^{(i+1)}$ and incorporating completes the proof.

Proof of Lemma \[lemma:upper\_bounded\] {#append:upper_bounded}
=======================================

Let ${\boldsymbol{Z}}^{(0)}={\boldsymbol{H}}^{(0)}$, ${\boldsymbol{Z}}^{(i+1)}=F\left({\boldsymbol{Z}}^{(i)},\tilde{{\boldsymbol{L}}}\right)$, where $\tilde{{\boldsymbol{L}}}$ is chosen such that all eigenvalues of $\left({\boldsymbol{A}}+{\boldsymbol{B}}\tilde{{\boldsymbol{L}}}_x\right)$ are inside the unit circle. Note that existence of $\tilde{{\boldsymbol{L}}}$ is guaranteed due to $\left({\boldsymbol{A}}, {\boldsymbol{B}}\right)$ controllable. With ${\boldsymbol{W}}^{(i)}=\tilde{{\boldsymbol{L}}}$ in Lemma \[lemma:monotonically\_increasing\], ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{Z}}^{(i)}$ holds. With $$\begin{aligned}
    {\boldsymbol{Z}}^{(i+1)}-{\boldsymbol{Z}}^{(i)}&=F\left({\boldsymbol{Z}}^{(i)},\tilde{{\boldsymbol{L}}}\right)-F\left({\boldsymbol{Z}}^{(i-1)},\tilde{{\boldsymbol{L}}}\right)\nonumber \\
    &=\gamma{\boldsymbol{M}}\left(\tilde{{\boldsymbol{L}}}\right)^\intercal\left({\boldsymbol{Z}}^{(i)}-{\boldsymbol{Z}}^{(i-1)}\right){\boldsymbol{M}}\left(\tilde{{\boldsymbol{L}}}\right),
    \end{aligned}$$ $\text{vec}({{\boldsymbol{\cdot}}})$ stacking the columns of a matrix and $\otimes$ being the Kronecker product, $$\begin{aligned}
    &\text{vec}\left({\boldsymbol{Z}}^{(i+1)}-{\boldsymbol{Z}}{(i)}\right)\nonumber \\&=\underbrace{\gamma{\boldsymbol{M}}\left(\tilde{{\boldsymbol{L}}}\right)^\intercal\otimes {\boldsymbol{M}}\left(\tilde{{\boldsymbol{L}}}\right)}_{={\boldsymbol{E}}}\text{vec}\left({\boldsymbol{Z}}^{(i)}-{\boldsymbol{Z}}^{(i-1)}\right)
    \end{aligned}$$ follows, thus $$\begin{aligned}
\label{eq:Ehochi1}
    \text{vec}\left({\boldsymbol{Z}}^{(i)}-{\boldsymbol{Z}}^{(i-1)}\right)&={\boldsymbol{E}}^{i-1}\text{vec}\left({\boldsymbol{Z}}^{(1)}-{\boldsymbol{Z}}^{(0)}\right).
    \end{aligned}$$ If all eigenvalues of $\sqrt{\gamma}{\boldsymbol{M}}\left(\tilde{{\boldsymbol{L}}}\right)$ are inside the unit circle, this also holds for the eigenvalues of ${\boldsymbol{E}}$. Due to its specific structure (cf. ), $(N+1)n$ eigenvalues of ${\boldsymbol{M}}\left(\tilde{{\boldsymbol{L}}}\right)$ are at the origin. Therefore, consider the remaining eigenvalues, i.e. the eigenvalues of ${\boldsymbol{D}}={\begin{bmatrix*}[c]{\boldsymbol{A}}&{\boldsymbol{B}}\\ \tilde{{\boldsymbol{L}}}_x{\boldsymbol{A}}&\tilde{{\boldsymbol{L}}}_x{\boldsymbol{B}}
\end{bmatrix*}}$ (cf. [@Landelius.1997 Lemma B.1.2]). Let ${\left\lVert {\boldsymbol{\cdot}} \right\rVert}$ be the spectral norm of a matrix, respectively the euclidean norm of a vector. Then, $$\begin{aligned}
    \lim_{i\rightarrow \infty}{\left\lVert {\boldsymbol{D}}^i \right\rVert}&=\lim_{i\rightarrow \infty}{\left\lVert {\begin{bmatrix*}[c]{\boldsymbol{I}}_n\\\tilde{{\boldsymbol{L}}}_x
\end{bmatrix*}}\left({\boldsymbol{A}}+{\boldsymbol{B}}\tilde{{\boldsymbol{L}}}_x\right)^{i-1}{\begin{bmatrix*}[c]{\boldsymbol{A}}&{\boldsymbol{B}}
\end{bmatrix*}} \right\rVert}\nonumber \\
    &\leq \lim_{i\rightarrow \infty}{\left\lVert {\begin{bmatrix*}[c]{\boldsymbol{I}}_n\\\tilde{{\boldsymbol{L}}}_x
\end{bmatrix*}} \right\rVert}{\left\lVert {\begin{bmatrix*}[c]{\boldsymbol{A}}&{\boldsymbol{B}}
\end{bmatrix*}} \right\rVert}{\left\lVert \left({\boldsymbol{A}}+{\boldsymbol{B}}\tilde{{\boldsymbol{L}}}_x\right)^{i-1} \right\rVert}\nonumber \\ &=0.
    \end{aligned}$$ As $\lim_{i\rightarrow \infty}{\boldsymbol{D}}^i={\boldsymbol{0}}$, all eigenvalues of ${\boldsymbol{D}}$ are inside the unit circle. Hence, all eigenvalues of ${\boldsymbol{E}}$ are also inside the unit circle and $e={\left\lVert {\boldsymbol{E}} \right\rVert}<1$. With $$\begin{aligned}
    \text{vec}\left({\boldsymbol{Z}}^{(j)}\right)&=\text{vec}\left({\boldsymbol{Z}}^{(0)}\right)+\sum_{i=1}^{j}\text{vec}\left({\boldsymbol{Z}}^{(i)}-{\boldsymbol{Z}}^{(i-1)}\right)\nonumber \\ &\overset{\eqref{eq:Ehochi1}}=\text{vec}\left({\boldsymbol{Z}}^{(0)}\right)+\sum_{i=1}^{j}{\boldsymbol{E}}^{i-1}\text{vec}\left({\boldsymbol{Z}}^{(1)}-{\boldsymbol{Z}}^{(0)}\right)
    \end{aligned}$$ follows $$\begin{aligned}
    &{\left\lVert \text{vec}\left({\boldsymbol{Z}}^{(j)}\right) \right\rVert}\nonumber \\&\leq{\left\lVert \text{vec}\left({\boldsymbol{Z}}^{(0)}\right) \right\rVert}+\sum_{i=0}^{\infty}{\left\lVert {\boldsymbol{E}} \right\rVert}^i{\left\lVert \text{vec}\left({\boldsymbol{Z}}^{(1)}-{\boldsymbol{Z}}^{(0)}\right) \right\rVert}=e_0,
    \end{aligned}$$ where the upper bound $e_0$ is independent of $j$. As ${\left\lVert \text{vec}\left({\boldsymbol{Z}}^{(j)}\right) \right\rVert}$ is upper bounded by $e_0$, there exists $e_1$ such that ${\left\lVert {\boldsymbol{Z}}^{(j)} \right\rVert}\leq e_1$, $\forall j$. With ${\boldsymbol{Y}}=e_1{\boldsymbol{I}}_{\text{dim}({\boldsymbol{H}})}$, ${\boldsymbol{0}}\preceq{\boldsymbol{H}}^{(i)}\preceq{\boldsymbol{Z}}^{(i)}\preceq{\left\lVert {\boldsymbol{Z}}^{(i)} \right\rVert}{\boldsymbol{I}}_{\text{dim}({\boldsymbol{H}})}\preceq e_1{\boldsymbol{I}}_{\text{dim}({\boldsymbol{H}})}={\boldsymbol{Y}}$ results which completes the proof.
