1 Introduction↩︎

Optimal controllers minimize a given cost function subject to a dynamic model. This design method has received a lot of attention, especially due to its potential in applications such as autonomous vehicle navigation and robotics, economics and management, as well as energy optimization [1], to cite a few. Traditionally, optimal control design relies on well-established techniques like dynamic programming, calculus of variations, and Pontryagin’s maximum principle [2]. While these methods have proven effective in many cases, they often face challenges when dealing with large-scale, nonlinear, or uncertain systems [[3]][4]. To address this issue, one often approximates a nonlinear system by a linear model around an operating point, thus enabling the design of a linear quadratic regulator (LQR) [5]. One of the main advantages of LQR [5] is that it has a closed-form solution [[6]][7] as a function of a matrix that is the positive definite solution of a Riccati equation [2]. An iterative technique for solving the Riccati equation is discussed in article [8], which requires an accurate linear model of the system. However, there are instances where complex nonlinear systems cannot be adequately approximated by linear models or where the system’s dynamics remain unknown. Therefore, conventional techniques for solving optimal control problems may not be applicable. As a result, researchers have explored data-driven methods such as RL, imitation learning, and Gaussian process regression as alternatives for solving optimal control problems [[9]][10][11]. RL is a powerful data-driven approach that enables an agent to learn optimal control policies by interacting with the environment and adjusting its actions or control policies to minimize a cumulative cost [12]. The motivation to use RL as an approxiation to optimal control designs stems from several key factors, including adaptability, the ability to handle complex dynamics, learning from real-world experience, and overcoming inaccuracies in system models [[13]][14]. The majority of research dedicated to RL in optimal control investigates discrete-time optimal control using value-based RL algorithms [15]. In value-based RL, the core concept is to estimate value functions, such as the Q-function [13], which can then be used to derive optimal policies. A common approach to estimate value functions is to use the temporal difference (TD) equation [16]. In such approach the Bellman error [17] is reduced by iteratively updating the value function approximation. This is called adaptive dynamic programming (ADP) [15]. ADP is widely favored due to its link with dynamic programming and the Bellman equation [18]. The application of ADP methods such as the policy iteration (PI) and the value-iteration (VI) to feedback control systems are discussed in references [[15]][19]. For optimal control, PI is preferred over VI due to its stability arising from the policy improvement guarantee, which leads to reliable convergence to the optimal policy [20]. Additionally, Q-learning, an ADP algorithm that approximates the Q-function, can be used to achieve optimal control without relying on a system model [21].

It is common to use a neural network (NN) to approximate the value function during the policy evaluation step of a PI algorithm [[22]][23]. However, the drawbacks of using NNs as the value function approximator (VFA) are (i) the optimization of the neural network’s weights is not convex and training the neural network yields only locally optimal weights, (ii) selecting the appropriate architecture for the NN typically involves a trial and error procedure, (iii) the input-output mapping of the NN lacks an analytical expression and as a consequence a second neural network is needed to perform policy improvement, (iv) providing a comprehensive proof of convergence to the optimal control is difficult or even unattainable. To address the mentioned issues, a two-layer quadratic neural network (QNN) trained by a convex optimization introduced in reference [24] will be chosen as the VFA. The advantages of using the two-layer QNN as the VFA compared to other neural networks are: (i) Two-layer QNNs are trained by solving a convex optimization and therefore the global optimal weights are found [24], (ii) the optimal number of neurons in the hidden layer is obtained ad a by-product [24], (iii) the input-output mapping of the QNN is a quadratic form [24], [25] and therefore one can analytically minimize the value-function with respect to the control policy. Reference [25] addresses applications of QNNs to system identification and control of dynamical systems.

LQR design problems have been addressed using various RL approaches as they can serve as benchmarks for evaluating the performance of RL algorithms [26], [27] because their solution is known in closed-form. This fact has motivated the use of ADP to solve discrete-time LQR in reference [28], discrete-time LQ Gaussian in references [[29]][30], and discrete-time LQ tracking in references [[31]][32] without using the system model. Moreover, in the Q-learning scheme proposed in reference [28] one can obtain the desired LQR controller by solving a least-squares optimization to derive the Q-function provided the system is controllable and a stabilizing linear state feedback controller is known.

This paper designs a discrete-time LQR controller using Q-learning. More specifically, we propose to use a two-layer QNN as the VFA in Q-learning for a LQR problem, for which the value function is known to be quadratic. To the best of our knowledge, this represents the first result with a convex optimization-trained QNN used as a VFA. Additionally, the designed controller is proven to converge to the LQR controller provided the system is controllable and an initial stabilizing policy is known. Simulations conducted in MATLAB illustrate the convergence of the learning algorithm to the optimal controller for different initial stabilizing policies.

This paper is organized as follows. Section 2 presents the problem statement. Section 3 reviews two-layer QNNs. Section 4 presents the general form of the Q-learning algorithm. Section 5 contains the proposed approach with the proof of convergence. Section 6 focuses on examples, which is followed by the conclusions.

2 Problem statement↩︎

An unknown controllable linear system model is written as \[\begin{gather} x_{k+1}=Ax_k+Bu_k \label{linear95dynamics} \end{gather}\tag{1}\]

where \(x_k \in \mathbb{R}^{n_x}\) is the state vector, \(u_k \in \mathbb{R}^{n_u}\) is the input vector and \(A, \: B\) are unknown matrices constrained to be such that the pair \((A,B)\) is controllable. Define the policy \(\pi (.)\) as a linear mapping from the state vector to the input vector as \[\begin{gather} u_{k}=\pi (x_k) = -K^\pi x_k \end{gather}\]

where \(K^\pi\) is a matrix to be determined by the designer. Define the local cost as \[\begin{gather} c(x_k,u_k)= x_k^T Q x_k + u_k^T R u_k \end{gather}\]

where \(Q\), and \(R\) are both positive definite. The value function when one follows the control policy \(\pi(.)\) is defined as \[V^\pi(x_k) = \sum_{i=k}^{\infty} \gamma^{i-k} \left( x_i^T Q x_i + u_i^T R u_i \right)\]

where \(0<\gamma \leq 1\) is a discount factor.

The objective is to obtain the optimal policy \(\pi^* (.)\) that minimizes \(V^{\pi}(x_k)\) for all states \(x_k\) subject to the unknown dynamics of the system. This is done using QNNs as a VFA.

3 Two-layer QNNs↩︎

Consider the neural network in Fig 1 with one hidden layer, one output, and a degree two polynomial activation function, where \(\mathcal{X}_i \in \mathbb{R}^n\) is the \(i\)-th input data vector, \(\mathcal{\hat{Y}}_i \in \mathbb{R}\) is the corresponding output, \(\mathcal{Y}_i \in \mathbb{R}\) is the output label corresponding to the input \(\mathcal{X}_i\), \(L\) is the number of hidden neurons, \(f(z)=az^2+bz+c\) is the polynomial activation function, and \(a \neq 0\), \(b\), and \(c\) are pre-defined constant coefficients. The notation \(w_{kj}\) represents the weight from the \(k\)-th input-neuron to the \(j\)-th hidden-neuron, and \(\rho_{j}\) represents the weight from the j-th hidden-neuron to the output. The input-output mapping is

\[\mathcal{\hat{Y}}_i=\sum_{j=1}^{L} f(\mathcal{X}_i^TW_j)\rho_{j}\] where \(W_j = \begin{bmatrix} w_{1j} & w_{2j} & \ldots & w_{nj} \end{bmatrix}^T\).

Reference [24] proposes the training optimization \[\begin{align} \label{eq:32changed32neural32net32optimizaion} \min_{W_k,\rho_k} \: \: & l(\mathcal{\hat{Y}}-\mathcal{Y}) + \beta \sum_{j=1}^{L} | \rho_j |\\ s.t. \: \: & \mathcal{\hat{Y}}_i=\sum_{j=1}^{L} f(\mathcal{X}_i^TW_j)\rho_{j}\\ & \Vert W_k \Vert_2 =1, \: \: \: \: k=1,2,...,L, \: \: \: \:i=1,2,...,N \end{align}\tag{2}\]

where \(\beta\ge 0\) is a pre-defined regularization parameter, \(l(.)\) is a convex loss function, N is the number of data points, \(\mathcal{\hat{Y}} = \begin{bmatrix} \mathcal{\hat{Y}}_1 & \mathcal{\hat{Y}}_2 & \ldots & \mathcal{\hat{Y}}_N \end{bmatrix}^T\), and \(\mathcal{Y} = \begin{bmatrix} \mathcal{Y}_1 & \mathcal{Y}_2 & \ldots & \mathcal{Y}_N \end{bmatrix}^T\).

The optimization problem 2 can be equivalently solved by the dual convex optimization [24] \[\begin{align} \label{convex95formula} \min_{Z^+ , Z^-} \: \: & l(\mathcal{\hat{Y}} - \mathcal{Y}) + \beta \: Trace(Z_{1}^+ + Z_{1}^-) \\ s.t. \: \: & \mathcal{\hat{Y}}_{i}=a\mathcal{X}_{i}^T(Z_{1}^+ - Z_{1}^-)\mathcal{X}_{i} + b\mathcal{X}_{i}^T(Z_{2}^+ - Z_{2}^-) + \\ & \: \: \: \: \: \: \: \: \: \: c \: Trace(Z_{1}^+ - Z_{1}^-), \\ & Z^+=\begin{bmatrix} Z_{1}^+ & Z_{2}^+ \\(Z_{2}^+)^T &Trace(Z_{1}^+) \end{bmatrix} \geq 0 , \\ & Z^-=\begin{bmatrix} Z_{1}^- & Z_{2}^- \\(Z_{2}^-)^T &Trace(Z_{1}^-) \end{bmatrix} \geq 0, \\ & i=1, 2, \ldots N \end{align}\tag{3}\] where \(Z_1^+\), \(Z_2^+\), \(Z_1^-\), \(Z_2^-\) are optimization parameters [24]. After training the neural network and obtaining \(Z^+\), \(Z^-\) from 3 , the quadratic input-output mapping is [24], [25]

\[\label{quadratic95mapping95QNN} \mathcal{\hat{Y}}_i=\begin{bmatrix} \mathcal{X}_i \\ 1 \end{bmatrix}^T H \begin{bmatrix} \mathcal{X}_i \\ 1 \end{bmatrix}\tag{4}\] where \[\begin{align} H=\begin{bmatrix} a(Z_{1}^+ - Z_{1}^-) & 0.5b(Z_{2}^+ - Z_{2}^-) \\0.5b(Z_{2}^+ - Z_{2}^-)^T & c \left [ Trace \left (Z_{1}^+ - Z_{1}^- \right ) \right ] \end{bmatrix} \end{align}\]

4 policy iteration-based Q-learning↩︎

This section presents the Q-learning algorithm used to calculate the optimal policy \(\pi^* (.)\). The Q-function [33] is \[Q^{\pi}(x_k,u_k) = c(x_k,u_k) + \gamma V^{\pi}(x_{k+1}) \label{eqn:32Q-function95def}\tag{5}\]

Using Bellman’s principle of optimality [15] and the definition of Q-function, equation 5 can be rewritten as \[\begin{gather} Q^{\pi}(x_k,u_k) = c(x_k,u_k) + \gamma Q^{\pi}(x_{k+1},\pi(x_{k+1})) \label{eqn:32Bellman-eqn} \end{gather}\tag{6}\] Obtaining the Q-function by equation 6 is known as the policy evaluation step [21]. After the policy evaluation step, the optimal stabilizing policy \(\pi^\prime (.)\) for \(Q^\pi\) is given by \[\begin{gather} \pi^\prime (x_k) = \arg \min_{u_k} Q^{\pi}(x_k,u_k) \label{eq:32policy32improvement32vanila} \end{gather}\tag{7}\] where \(V^{\pi^\prime}(x_k) < V^{\pi}(x_k)\) for all \(x_k\) [23]. Solving 7 is known as the policy improvement step. The optimal policy \(\pi^* (.)\) can be attained from any initial stabilizing policy by repeatedly performing policy evaluation and policy improvement steps until convergence, starting from an initial stabilizing policy \(\pi_0(.)\). Algorithm 2 is the PI-based Q-learning algorithm.

Solving equation [eq:32PEA1] can be challenging as it is a Lyapunov function [15]. However, exploiting its fixed-point nature, we can obtain \(Q^\pi(.)\) for the stabilizing policy \(\pi(.)\) as \[\label{eq:32iterative32Bellman32van} Q^{\pi}_{i+1}(x_{k},u_k) = c(x_k,u_k) + \gamma Q^{\pi}_{i}(x_{k+1},\pi(x_{k+1}))\tag{8}\] Under certain conditions that will be detailed later, starting with any initial value \(Q^\pi_0\) and using the iteration (8 ) will lead to \(Q^\pi_i \rightarrow Q^\pi\). Consequently, equation [eq:32PEA1] will be replaced by the iterative equation 8 .

5 The proposed approach↩︎

This section presents how to perform policy evaluation and policy improvement steps by designing a two-layer QNN. According to reference [28], the Q-function is a quadratic form \[\begin{gather} \label{Q95function95linear95systems} Q^{\pi}(x_k,u_k) = \begin{bmatrix} x_k \\ u_k \end{bmatrix}^T H^\pi \begin{bmatrix} x_k \\ u_k \end{bmatrix} \end{gather}\tag{9}\] where \(H^\pi \in \mathbb{R}^{(n_x + n_u) \times (n_x + n_u)}\). As a result, a QNN with coefficients \(b=c=0\), \(a=1\) is the perfect candidate to approximate the Q-function. The policy evaluation step then obtains \(H^\pi\) for the stabilizing policy \(\pi(.)\) by solving, \[\begin{gather} \label{eq:32policy95evaluation95vanilla95LQR950} \begin{bmatrix} x_k \\ u_k \end{bmatrix}^T {H}^\pi \begin{bmatrix} x_k \\ u_k \end{bmatrix} = x_k^T Q x_k + u_k ^T R u_k + \\ \gamma \begin{bmatrix} x_{k+1} \\ u_{k+1} \end{bmatrix}^T {H}^\pi \begin{bmatrix} x_{k+1} \\ u_{k+1} \end{bmatrix} \end{gather}\tag{10}\]

We propose to obtain \(H^\pi\) in the policy evaluation step by employing the iterative equation 11 . The proof of convergence is presented in Lemma 1.

Due to the implications of Lemma 1, the problem of policy evaluation transforms into the task of computing \(H^\pi_{i+1}\) from equation 11 given \(H^\pi_{i}\). This can be done by training a two-layer QNN. Since it is assumed that one has access to the state \(x_k\), one can calculate \(\mathcal{Y}_k\) defined as \[\begin{gather} \mathcal{Y}_k = x_k^T Q x_k + u_k^T R u_k + \gamma \begin{bmatrix} x_{k+1} \\ u_{k+1} \end{bmatrix}^T {H}_{i}^\pi \begin{bmatrix} x_{k+1} \\ u_{k+1} \end{bmatrix}. \end{gather}\] Therefore, from (11 ), \[\mathcal{X}_k^T {H}_{i+1}^\pi \mathcal{X}_k = \mathcal{Y}_k\] where \(\mathcal{X}_k^T = \begin{bmatrix} x_k^T & u_{k}^T \end{bmatrix}\). Thus, \(H^\pi_{i+1}\) can be obtained from the training of a QNN as the solution of the convex optimization 3 using a set of input data points \(\mathcal{X}_k\) and the corresponding labels \(\mathcal{Y}_k\), as well as the coeffficients \(a=1, \: b=0, \: c=0\).

5.1 Policy improvement:↩︎

In this section, the policy improvement step is addressed for the stabilizing policy \(\pi(\cdot)\) using \(H^\pi\) from the policy evaluation. We first partition \(H^\pi\) as \[\begin{gather} \begin{bmatrix} x_k \\ u_k \end{bmatrix}^T {H}^\pi \begin{bmatrix} x_k \\ u_k \end{bmatrix} = \begin{bmatrix} x_k \\ u_k \\ \end{bmatrix}^T \begin{bmatrix} {H}_{xx}^\pi & {H}_{xu}^\pi \\ ({H}_{xu}^\pi)^T & {H}_{uu}^\pi \end{bmatrix} \begin{bmatrix} x_k \\ u_k \end{bmatrix} \end{gather}\]

The complete procedure for solving the LQR is presented in Algorithm 3.

6 Simulations↩︎

This sectiom considers a quadrotor flying at a constant altitude movimg in the \(x_1\) direction. It is shown that the control policy converges to the LQR solved by conventional methods. The quadrotor is modeled by the state-space \[\begin{gather} \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3\\ \dot{x}4 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & -0.1 & 10 & 0\\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4 \end{bmatrix} + \begin{bmatrix} 0 \\ 0\\ 0\\ 4.35 \end{bmatrix}u \end{gather}\] where \(x_1\) is the horizontal position, \(x_2\) is the speed, \(x_3\) is the pitch angle, \(x_4\) is the pitch rate and the initial state is \[\begin{bmatrix} x_1(0) \\ x_2(0) \\ x_3(0) \\ x_4(0) \end{bmatrix} = \begin{bmatrix} -10 \\ 0 \\ 0 \\ 0 \end{bmatrix}\] Assume that the sampling time is \(T=0.1\) seconds. The discretized model 19 is obtained using Tustin’s method in order to compare the result with the discrete-time LQR. \[\begin{gather} \begin{bmatrix} x_{1,k+1} \\ x_{2,k+1} \\ x_{3,k+1} \\ x_{4,k+1} \end{bmatrix} = \begin{bmatrix} 1 & 0.1 & 0.05 & 0.003 \\ 0 & 0.99 & 0.99 & 0.05 \\ 0 & 0 & 1 & 0.1 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_{1,k } \\ x_{2,k } \\ x_{3,k } \\ x_{4,k } \end{bmatrix} + \begin{bmatrix} 0.001 \\ 0.011 \\ 0.022 \\ 0.435 \end{bmatrix}u_k \label{eq:32linear95quadrotor95tustin} \end{gather}\tag{19}\]

Rewrite the control policy \(\pi_j(.)\) as: \[\begin{gather} \pi_j(x_k)= \begin{bmatrix} k_1^{\pi_j} & k_2^{\pi_j} & k_3^{\pi_j} & k_4^{\pi_j} \end{bmatrix} \begin{bmatrix} x_{1,k} \\ x_{2,k} \\ x_{3,k} \\ x_{4,k} \end{bmatrix} \label{eq:32shand} \end{gather}\tag{20}\] In this example, we choose the following parameter values: \[\begin{gather} \nonumber R = 100,\: N=100, \: \gamma=1, \:\beta=0.005 \\ Q = \begin{bmatrix} 0.01 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 10 \end{bmatrix} \end{gather}\]

The optimal LQR controller is \[\begin{gather} \pi^*(x_k)=- \begin{bmatrix} 0.046 & 0.464 & 4.347 & 2.014 \end{bmatrix} \begin{bmatrix} x_{1,k} \\ x_{2,k} \\ x_{3,k} \\ x_{4,k} \end{bmatrix} \end{gather}\] We now pretend that we do not know the model of the quadrotor and run the Algorithm 3. Fig 4 demonstrates the convergence of \(k_1^{\pi_j}\), \(k_2^{\pi_j}\), \(k_3^{\pi_j}\), \(k_4^{\pi_j}\) to their optimal values. To illustrate this convergence, we conducted five simulations using random initial stabilizing policies. The random initial policies are given in the Table 1.

The trajectory of the quadrotor’s position and its speed over time using the optimal policy are depicted in Fig. 5 where one sees the convergence to the origin.

References↩︎