Data-Driven Robust Control Using Reinforcement Learning

Abstract

This paper proposes a robust control design method using reinforcement learning for controlling partially-unknown dynamical systems under uncertain conditions. The method extends the optimal reinforcement learning algorithm with a new learning technique based on the robust control theory. By learning from the data, the algorithm proposes actions that guarantee the stability of the closed-loop system within the uncertainties estimated also from the data. Control policies are calculated by solving a set of linear matrix inequalities. The controller was evaluated using simulations on a blood glucose model for patients with Type 1 diabetes. Simulation results show that the proposed methodology is capable of safely regulating the blood glucose within a healthy level under the influence of measurement and process noises. The controller has also significantly reduced the post-meal fluctuation of the blood glucose. A comparison between the proposed algorithm and the existing optimal reinforcement learning algorithm shows the improved robustness of the closed-loop system using our method.

1. Introduction

Control of unknown dynamic systems with uncertainties is a challenge because exact mathematical models are often required. Since many processes are complicated, nonlinear, and varying with time, a control algorithm that does not depend on a mathematical model and can adapt to time-varying conditions is required. A popular approach is to develop a universal approximator for predicting the output of unknown systems [1]. Control algorithms can then be designed based on the parameters of the approximator. Based on this approach, many control techniques have been proposed using machine learning models such as neural networks and fuzzy logic. For example, Goyal et al. [2] proposed a robust sliding mode controller which can be designed from Chebyshev neural networks. Chadli and Guerra [3] introduced a robust static output feedback controller for Takagi Sugeno fuzzy models. Ngo and Shin [4] proposed a method to model unstructured uncertainties and a new Takagi Sugeno fuzzy controller using type-2 fuzzy neural networks.

However, obtaining a good approximator requires a significant amount of training data, especially for a complicated model with high-dimensional state spaces or with many inputs and outputs. The data-driven model must also be updated frequently for time-varying systems. In addition, many control design techniques assume uncertainties as functions of system parameters. However, in many cases, the causes of uncertainties are unknown and unstructured. With the development of data science and machine learning, model-free approaches such as reinforcement learning (RL) have emerged as an effective method to control unknown nonlinear systems [5,6,7,8]. The principle of RL is based on the interaction between a decision-making agent and its environment [9], and the actor–critic method is often used as the RL framework for many control algorithms. In the actor–critic framework, the critic agent uses current state information of the environment in order to update the value or action value function. The actor agent then uses the value or action value function to calculate the optimal action.

It can be seen that many data-driven algorithms lack stability analysis of the closed-loop systems. Among recent techniques focusing on the robustness of control algorithms, Yang et al. [10] presented an off-policy reinforcement learning (RL) solution to solve robust control problems for a certain class of unknown systems with structured uncertainties. In [11], a robust data-driven controller was proposed based on the frequency response of multivariable systems and convex optimization. Based on data-driven tuning, Takabe et al. [12] introduced a detection algorithm suitable for massive overloaded multiple-input multiple-output systems. In more recent works, Na et al. [13] proposed an approach to address the output-feedback robust control for continuous-time uncertain systems using online data-driven learning, while Makarem et al. [14] used data-driven techniques for iterative feedback tuning of a proportional-integral-derivative controller’s parameters. However, in many cases, stability can only be ensured for specific systems where uncertainties are structured. In addition, the value function must be estimated accurately, which is difficult to achieve, especially at the beginning of the control process when the agent has just started interacting with the environment. Additionally, in many applications, the state space is either continuous or high-dimensional. In these cases, the value function approximation is often inaccurate, potentially leading to instability. Therefore, new RL approaches for which stability can be guaranteed under uncertain conditions are essential if algorithms are to be used in critical and safety-demanding systems.

Type 1 diabetes is a disease caused by the lack of insulin secretion. The condition results in uncontrolled increase of blood glucose level if the patients are not provided with insulin doses. High blood glucose level can lead to both acute and chronic complications, and eventually result in failure of various organs. One of the major challenges in controlling the blood glucose is that the biochemical and physiologic kinetics of insulin and glucose is complicated, nonlinear, and only approximately known [15]. Additionally, the stability of the control system is essential in this case since unstable control effort will lead to life-threatening condition for the patients.

This paper proposes a novel method to capture uncertainty in estimating the value function in reinforcement learning based on observation data. Using the uncertainty information, the paper also presents a new technique to improve the policy while guaranteeing the stability of the closed-loop system under uncertainty conditions for partially-unknown dynamical systems. The proposed methodology is applied to a blood glucose model for testing its effectiveness in controlling the blood glucose level in patients with Type 1 diabetes.

Structure of Paper

The content of the paper is organized as follows. Section 2 describes the proposed robust RL algorithm. Section 3 shows the simulation results of the methodology. The conclusions are given in Section 4.

2. Materials and Methods

In this section we present the robust RL method and the simulation setup used for evaluation of the algorithm.

2.1. Robust Control Using Reinforcement Learning

In this paper, a class of dynamical systems is considered, which can be described by the following linear state-space equation:

$$\dot{x}\left(t\right)=Ax\left(t\right)+Bu\left(t\right),$$

(1)

where $x\in {\mathbb{R}}^{\mathit{n}}$ is the vector of n state variables, $u\in {\mathbb{R}}^{\mathit{m}}$ is the vector of m control inputs, $A\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}$ is the state matrix, and $B\in {\mathbb{R}}^{\mathit{n}\times \mathit{m}}$ is the input matrix. It is assumed that matrix A is a squared $n\times n$ unknown matrix and the system (A and B) is stabilized. Our target is to derive a control algorithm $u\left(t\right)$ that can regulate the state variables contained in $x\left(t\right)$ based on input and output data without knowing matrix A.

As an RL framework, the proposed robust control algorithm consists of an agent that takes actions and learns the consequences of its actions in an unknown environment. The environment is defined by a state vector $x\left(t\right)$ that describes its states at time t. The action at time t is represented by $u\left(t\right)$. As a consequence of the action, a cost $r\left(t\right)$ is incurred and accumulated. The cost function $r\left(t\right)$ is assumed to be known and predefined as a function of the current state and action. The objective of the learning process is to minimize the total cost accumulation in the future.

At each decision time point, the agent receives information about the state of the environment and chooses an action. The environment reacts to this action and transitions to a new state, which determines whether the agent receives a positive or negative reinforcement. Current RL techniques propose optimal actions by minimizing the predicted cost accumulation. However, uncertainties due to noises in the data or inaccurate estimation of the cost accumulation can lead to suboptimal actions and even unstable responses. Our target is to provide the agent with a robust and safe action that can guarantee the reduction of the future cost accumulation in the presence of uncertainties. The action calculated by the proposed algorithm may not be the optimal action that reduces the cost in the fastest way, but it can always guarantee the stability of the system, which is imperative in many critical applications.

2.1.1. Estimation of the Value Function by the Critics

In the RL context, the accumulation of cost over time, when starting in the state $x\left(t\right)$ and following policy $\pi$, is defined as the value function of policy $\pi$, i.e.,

$$V^{\pi }\left(x\left(t\right)\right)=E_{\pi }\left\{{\int }_t^{\infty }{\gamma }^{\tau -t}r\left(\tau \right)d\tau \right\},$$

(2)

where $\gamma$ is the discount factor. The cost $r\left(t\right)$ is assumed to be a quadratic function of the states:

$$r\left(t\right)=x^T\left(t\right)Qx\left(t\right),$$

(3)

where the positive definite matrix $Q\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}$ is symmetric, positive semidefinite (since the cost is assumed to be non-negative), and contains the weighting factors of the variables that are minimized.

In order to facilitate the formulation of the stability condition in the form of linear matrix inequalities (LMI), the value function $V\left(x\right(t\left)\right)$ is approximated by a quadratic function of the states:

$$V^{\pi }\left(x\left(t\right)\right)\approx x^T\left(t\right)Px\left(t\right),$$

(4)

where the kernel matrix $P\in {\mathbb{R}}^{\mathit{n}\times \mathit{n}}$ is symmetric and positive semidefinite (since matrix Q in the cost function is symmetric and positive semidefinite).

By using the Kronecker operation, the approximated value function can be expressed as a linear combination of the basis function $\varphi \left(x\right(t\left)\right)=\left(x\right(t)\otimes x(t\left)\right)$:

$$\begin{array}{cc}\hfill V^{\pi }\left(x\left(t\right)\right)& \approx x^T\left(t\right)Px\left(t\right)=vec{\left(P\right)}^T(x\left(t\right)\otimes x\left(t\right))\hfill \\ \hfill & =w^T(x\left(t\right)\otimes x\left(t\right))=w^T\varphi \left(x\left(t\right)\right),\hfill \end{array}$$

(5)

where w is the parameter vector, $\varphi \left(x\right(t\left)\right)$ is the vector of basis functions, and ⊗ is the Kronecker product. The transformation between w and P can be performed as follows:

$$w=vec\left(P\right)={[P_{11},P_{21},\dots ,P_{n1},P_{12},\dots ,P_{1n},P_{nn}]}^T,$$

(6)

where $P_{i,j}$ is the element of matrix P in the $i\mathrm{th}$ row and $j\mathrm{th}$ column. With T as the interval time for data sampling, the integral RL Bellman equation can be used to update the value function [8]:

$$V^{\pi }\left(x\left(t\right)\right)={\int }_t^{t+T}{\gamma }^{\tau -t}r\left(\tau \right)d\tau +V^{\pi }\left(x(t+T)\right).$$

(7)

By using the quadratic cost function (Equation (3)) and the approximated value function (Equation (5)), the integral RL Bellman equation can be written as follows:

$$x^T\left(t\right)Px\left(t\right)={\int }_t^{t+T}x{\left(\tau \right)}^{T}Qx\left(\tau \right)d\tau +x^T(t+T)Px(t+T)$$

(8)

or

$$\begin{array}{cc}\hfill w^T\varphi \left(x\left(t\right)\right)={\int }_t^{t+T}& x{\left(\tau \right)}^{T}Qx\left(\tau \right)d\tau +w^T\varphi \left(x(t+T)\right).\hfill \end{array}$$

(9)

At each iteration, n samples along the state trajectory are collected ($x^{\left(1\right)}\left(t\right),x^{\left(2\right)}\left(t\right),$$\dots ,x^{\left(n\right)}\left(t\right)$). The mean value of w can be obtained by using least-square technique:

$$\widehat{w}=\left(X{X}^T\right)XY,$$

(10)

where

$$X={[{\varphi }_{\Delta }^1{\varphi }_{\Delta }^2\dots {\varphi }_{\Delta }^N]}^T,$$

(11)

$${\varphi }_{\Delta }^i=\varphi \left(x^i\left(t\right)\right)-\varphi \left(x^i(t+T)\right),$$

(12)

$$Y={[d\left(x^1\left(t\right)\right)d\left(x^2\left(t\right)\right)\dots d\left(x^n\left(t\right)\right)]}^T$$

(13)

and

$$d\left(x^i\left(t\right)\right)={\int }_t^{t+T}{x}^i{\left(\tau \right)}^{T}Q{x}^i\left(\tau \right)d\tau$$

(14)

with $i=1,2,\dots ,N$.

The confidence interval for the coefficient $w^{\left(j\right)}$ is given by

$$w^{\left(j\right)}\in [{\widehat{w}}^{\left(j\right)}-q_{1-\frac{\theta }{2}}\sqrt{{\tau }_j{\widehat{\sigma }}^2},{\widehat{w}}^{\left(j\right)}+q_{1-\frac{\theta }{2}}\sqrt{{\tau }_j{\widehat{\sigma }}^2}],$$

(15)

where $1-\theta$ is the confidence level, $q_{1-\frac{\theta }{2}}$ is the quantile function of standard normal distribution, ${\tau }_j$ is the jth element on the diagonal of ${\left(X{X}^T\right)}^{-1}$, and ${\widehat{\sigma }}^2=\frac{{\widehat{ϵ}}^T\widehat{ϵ}}{n-p}$, with $ϵ=Y-\widehat{w}X$. From that, the uncertainty $\Delta w$ is defined as the deviation interval around the nominal value:

$$\Delta w=\left[-q_{1-\frac{\theta }{2}}\sqrt{{\tau }_j{\widehat{\sigma }}^2},-q_{1-\frac{\theta }{2}}\sqrt{{\tau }_j{\widehat{\sigma }}^2}\right].$$

(16)

Matrices $\widehat{P}$ and $\Delta P$ can be obtained by placing elements of $\widehat{w}$ and $\Delta w$ into columns.

2.1.2. Policy Improvement by the Actor

Linear feedback controllers have been widely used as a stabilization tool for nonlinear systems where dynamic behavior is considered approximately linear around the operating condition [16,17,18]. Hence, in this paper, we use linear functions of the states with gain $K_i$ as the control policy at iteration i:

$$u\left(t\right)=\pi \left(x\left(t\right)\right)=-K_{i}x\left(t\right),$$

(17)

and the level of uncertainty is constant during the controlling process. The task of the actor is to robustly improve the current policy such that the value function is guaranteed to be reduced during the next policy implementation. If the following differential inequality is satisfied:

$${\dot{V}}_i\left(x\left(t\right)\right)+\alpha V_i\left(x\left(t\right)\right)\le 0$$

(18)

with some positive constant $\alpha$, then by using the comparison lemma (Lemma 3.4 in [19]), the derivative of function ${\dot{V}}_i\left(x\left(t\right)\right)$ can be bounded by

$${\dot{V}}_i\left(x\left(t\right)\right)\le V_i\left(x\left(t_0\right)\right)e^{-\alpha (t-t0)}.$$

(19)

Therefore, maximizing the rate $\alpha$ will ensure a maximum exponential decrease in the value of ${\dot{V}}_i\left(x\left(t\right)\right)$.

The following part shows the main results of the paper, which describe how the policy gain can be improved during the learning process. Derivations of the results are provided in the stability analysis (Section 2.1.3).

Definition 1.

Assume A is a square matrix with dimension $n\times n$ and x is a vector with dimension $n\times 1$. The maximize operation on matrix A and vector x is defined as follows:

$$\mathrm{maximize}(A,x)=C,$$

(20)

where

$$C_{ij}=\left\{\begin{array}{c}\max \left(A_{ij}\right)\mathrm{if}{x}_i{x}_j\ge 0\\ \min \left(A_{ij}\right)\mathrm{if}{x}_i{x}_j<0\end{array}\right.\mathrm{with}i,j=1\dots n.$$

(21)

Assuming that the sign of all state variables cannot be changed between each policy update interval, the improved policy $K_{i+1}$ can be obtained by minimizing $\alpha$ subject to

$$\begin{array}{c}\hfill \left[\begin{array}{cc}V& K_{i+1}^T{B}^T\\ B{K}_{i+1}& -\gamma I\end{array}\right]\le 0\end{array}$$

(22)

and

$$\left[\begin{array}{cc}\zeta & K_{i+1}\\ K_{i+1}^T& -I\end{array}\right]\le 0,$$

(23)

where:

$$V=M+\Delta P_{i,\max }^T\Delta P_{i,\max }\gamma -{\widehat{P}}_{i}B{K}_{i+1}-K_{i+1}^T{B}^T{\widehat{P}}_i+\alpha ({\widehat{P}}_i+\frac{1}{2}\Delta P_{i,\max }^T\Delta P_{i,\max }+I)$$

(24)

and

$$\begin{array}{cc}\hfill M=& -Q-K_i^{T}R{K}_i+{\widehat{P}}_{i}B{K}_i+K_i^T{B}^T{\widehat{P}}_i+H_i\hfill \end{array}$$

(25)

with $\Delta P_{i,\max }=\mathrm{maximize}(\Delta P_i,x)$ and $H_i=\mathrm{maximize}(\Delta P_{i}B{K}_i+K_i^T{B}^T\Delta P_i,x)$, utilizing the maximize operation defined in Definition 1. Inequality (22) provides the stable condition and its derivation is provided in Section 2.1.3. Inequality (23) provides the upper bound for the updated gain $K_{i+1}$ through the user-defined parameter $\zeta$. The value of $-\zeta$ limits the maximum $L^2$ gain of $K_{i+1}$ since inequality (23) is equivalent to $K_{i+1}{K}_{i+1}^T\le -\zeta$.

2.1.3. Stability Analysis

With the control policy as described in Equation (17), the equation for the closed-loop system can be derived as follows:

$$\dot{x}\left(t\right)=Ax\left(t\right)-BKx\left(t\right)=(A-BK)x\left(t\right).$$

(26)

Lemma 1.

Assuming that the closed-loop system described by Equation (26) is stable, solving for P in Equation (8) is equivalent to finding the solution of the underlying Lyapunov equation [8]:

$$P(A-BK)+{(A-BK)}^{T}P=-Q.$$

(27)

Proof of Lemma 1. 

We start with Equation (27) and try to prove that matrix P is also the solution of Equation (8). Consider $V\left(x\left(t\right)\right)=x^T\left(t\right)Px\left(t\right)$, where P is the solution of Equation (27):

$$\begin{array}{cc}\hfill \dot{V}\left(x\left(t\right)\right)& =\frac{d\left(x^T\left(t\right)Px\left(t\right)\right)}{dt}\hfill \\ \hfill & ={\dot{x}}^T\left(t\right)Px\left(t\right)+x^T\left(t\right)P\dot{x}\left(t\right)\hfill \\ \hfill & =x^T\left(t\right)\left[{(A-BK)}^{T}P+P(A-BK)\right]x\left(t\right)\hfill \\ \hfill & =-x^T\left(t\right)Qx\left(t\right)(\mathrm{using}\mathrm{Equation}\left(27\right)).\hfill \end{array}$$

(28)

Since the closed-loop system is stable, the Lyapunov Equation (27) has a unique solution, $P_i>0$. From (28), this solution will satisfy

$$\frac{d\left(x^T\left(t\right)P_{i}x\left(t\right)\right)}{dt}=-x^T\left(t\right)Qx\left(t\right),$$

(29)

which is equivalent to

$$x^T(t+T)Px(t+T)-x^T\left(t\right)Px\left(t\right)={\int }_t^{t+T}-x^T\left(\tau \right)Qx\left(\tau \right)d\tau .$$

(30)

Therefore, P is also the solution of Equation (8). □

Lemma 2.

Given matrices E and F with appropriate dimensions, the following LMI can be obtained:

$$E{F}^T+F{E}^T\le E{E}^T+F{F}^T.$$

(31)

Proof of Lemma 2. 

From the properties of matrix norm, we have

$$(E-F){(E-F)}^T\ge 0,$$

(32)

which is equivalent to

$$E{E}^T+F{F}^T-E{F}^T-F{E}^T\ge 0$$

(33)

or

$$E{F}^T+F{E}^T\le E{E}^T+F{F}^T.$$

(34)

□

Lemma 3.

Given A as a square matrix with dimension $n\times n$ and x as a vector with dimension $n\times 1$, the following LMI can be obtained:

$$x^{T}Ax\le x^{T}Cx,$$

(35)

where $C=\mathrm{maximize}(A,x)$ as in Definition 1.

Proof of Lemma 3. 

We have

$$\begin{array}{cc}\hfill x^{T}Ax& =\sum _{i,j=1,2\dots n}{a}_{ij}{x}_i{x}_j\le \sum _{i,j=1,2\dots n}\left|a_{ij}{x}_i{x}_j\right|\hfill \\ \hfill & =\sum _{i,j=1,2\dots n}{c}_{ij}{x}_i{x}_j=x^{T}Cx,\hfill \end{array}$$

(36)

where $c_{ij}=\left\{\begin{array}{c}\max \left(a_{ij}\right)\mathrm{if}{x}_i{x}_j\ge 0\\ \min \left(a_{ij}\right)\mathrm{if}{x}_i{x}_j<0\end{array}\right.$ with $i,j=1\dots n$. □

Theorem 1.

Consider a dynamic system that can be represented by Equation (1) with the state matrix A unknown. Assume that the sign of all state variables cannot be changed between each policy update interval and the estimated value function at iteration i is $V_i\left(x\left(t\right)\right)=x^T\left(t\right)P_{i}x\left(t\right)$ with $P_i={\widehat{P}}_i+\Delta P_i$. If

• The current control policy $u\left(t\right)={\pi }_i\left(x\left(t\right)\right)=-K_{i}x\left(t\right)$ is stabilizing;
• The LMI given in (22) is satisfied with some positive constant γ;

then the closed-loop system with the control policy $u\left(t\right)=-K_{i+1}x\left(t\right)$ is quadratic stable with convergence rate α.

Proof of Theorem 1. 

Since the current control policy is stable, the estimated parameter matrix $P_i$ is positive definite. Hence, $V_i\left(x\left(t\right)\right)=x_t^T{P}_i{x}_t>0$. Here, $V_i\left(x\left(t\right)\right)$ is used as the Lyapunov function for the updated control policy $u\left(t\right)={\pi }_{i+1}\left(x\left(t\right)\right)=-K_{i+1}x\left(t\right)$. For notation convenience, the state vector $x\left(t\right)$ and input vector $u\left(t\right)$ are denoted as $x_t$ and $u_t$, respectively. By using Equation (27) in Lemma 1 and the representation $P_i={\widehat{P}}_i+\Delta P_i$, we can calculate the left side of Equation (18) as follows:

$$\begin{array}{cc}\hfill {\dot{V}}_i& \left(x\left(t\right)\right)+\alpha V_i\left(x\left(t\right)\right)\hfill \\ \hfill =& {\dot{x}}_t^T{P}_i{x}_t+x_t{P}_i{\dot{x}}_t^T+\alpha x_t^T{P}_i{x}_t\hfill \\ \hfill =& {(A{x}_t+B{u}_t)}^T{P}_i{x}_t+x_t{P}_i{(A{x}_t+B{u}_t)}^T+\alpha x_t^T{P}_i{x}_t\hfill \\ \hfill =& x_t^T[P_i(A-B{K}_{i+1})+{(A-B{K}_{i+1})}^T{P}_i+\alpha P_i]x_t\hfill \\ \hfill =& x_t^T[P_i(A-B{K}_i)+{(A-B{K}_i)}^T{P}_i+\alpha P_i]x_t+x_t^T[P_{i}B(K_i-K_{i+1})+{(K_i-K_{i+1})}^T{B}^T{P}_i+\alpha P_i]x_t\hfill \\ \hfill =& -x_t^T[Q+K_i^{T}R{K}_i]x_t+x_t^T[({\overline{P}}_i+\Delta P_i)B(K_i-K_{i+1})+{(K_i-K_{i+1})}^T{B}^T({\overline{P}}_i+\Delta P_i)+\alpha {\overline{P}}_i+\alpha \Delta P_i]x_t\hfill \\ \hfill =& x_t^T[-Q-K_i^{T}R{K}_i+{\overline{P}}_{i}B{K}_i+K_i^T{B}^T{\overline{P}}_i+\alpha {\overline{P}}_i+\Delta P_{i}B{K}_i+K_i^T{B}^T\Delta P_i-\Delta P_{i}B{K}_{i+1}\hfill \\ \hfill -& K_{i+1}^T{B}^T\Delta P_i-{\overline{P}}_{i}B{K}_{i+1}-K_{i+1}^T{B}^T{\overline{P}}_i+\alpha \Delta P_i]x_t.\hfill \end{array}$$

By using Lemma 3, we have the following inequality:

$$\Delta P_{i}B{K}_i+K_i^T{B}^T\Delta P_i\le H_i,$$

(37)

and the following inequality can be obtained by Lemma 2:

$$\begin{array}{cc}\hfill -\Delta P_{i}B{K}_{i+1}-K_{i+1}^T{B}^T\Delta P_i& \le \gamma \Delta P_i\Delta P_i^T+\frac{1}{\gamma }{\left(B{K}_{i+1}\right)}^T\left(B{K}_{i+1}\right)\hfill \\ & \le \gamma \Delta P_{i,\max }\Delta P_{i,\max }^T+\frac{1}{\gamma }{K}_{i+1}^T{B}^{T}B{K}_{i+1}\hfill \end{array}$$

(38)

Additionally,

$$\alpha \Delta P_i\le \alpha \left(\frac{1}{2}\Delta P_i\Delta P_i^T+I\right)\le \alpha \left(\frac{1}{2}\Delta P_{i,\max }\Delta P_{i,\max }^T+I\right),$$

(39)

where $H_i=\mathrm{maximize}(\Delta P_{i}B{K}_i+K_i^T{B}^T\Delta P_i,x)$, and $\Delta P_{i,\max }=\mathrm{maximize}(\Delta P_i,x)$, utilizing the maximize operator defined in Definition 1.

Hence, ${\dot{V}}_i\left(x\left(t\right)\right)+\alpha V_i\left(x\left(t\right)\right)$ can be bounded by

$$\begin{array}{c}\hfill {\dot{V}}_i(x\left(t\right))+\alpha V_i\left(x\left(t\right)\right)\le x_t^T[-Q-K_i^{T}R{K}_i+{\overline{P}}_{i}B{K}_i+K_i^T{B}^T{\overline{P}}_i+\alpha \left({\overline{P}}_i+\frac{1}{2}\Delta P_{i,\max }\Delta P_{i,\max }^T\right)\\ \hspace{1em}\hspace{1em}-{\overline{P}}_{i}B{K}_{i+1}-K_{i+1}^T{B}^T{\overline{P}}_i+\gamma \Delta P_{i,\max }\Delta P_{i,\max }^T+\frac{1}{\gamma }{K}_{i+1}^T{B}^{T}B{K}_{i+1}]x_t.\hfill \end{array}$$

(40)

Using the Lyapunov theory, the system will be quadratic stable with the convergent rate $\alpha$ if ${\dot{V}}_i\left(x\left(t\right)\right)\le -\alpha V_i\left(x\left(t\right)\right)$. This condition is satisfied if

$$\begin{array}{cc}\hfill & x_t^T[-Q-K_i^{T}R{K}_i+{\overline{P}}_{i}B{K}_i+K_i^T{B}^T{\overline{P}}_i+\alpha \left({\overline{P}}_i+\frac{1}{2}\Delta P_{i,\max }\Delta P_{i,\max }^T\right)-{\overline{P}}_{i}B{K}_{i+1}-K_{i+1}^T{B}^T{\overline{P}}_i\hfill \\ \hfill & +\gamma \Delta P_{i,\max }\Delta P_{i,\max }^T+\frac{1}{\gamma }{K}_{i+1}^T{B}^{T}B{K}_{i+1}]x_t\le 0.\hfill \end{array}$$

The above condition can be written in the matrix form, as shown in Theorem 1. □

By using Theorem 1, it can be seen that with the proposed improved policy, the closed-loop system will be asymptotically stable. It is also noted that Theorem 1 is also applicable for unknown nonlinear systems if they can be approximated by a linear state-space equation (Equation (1)) and if their nonlinearity is within the uncertainty bound $\Delta P$ calculated from $\Delta w$ in Equation (16).

2.1.4. Robust Reinforcement Learning Algorithm

The robust RL algorithm for controlling partially unknown dynamically systems includes the following steps:

Initialization

(Step $i=0$)

• Select an initial policy $u\left(t\right)=-K_{0}x\left(t\right)$.

Estimation of the Value Function

(Step $i=1,2,\dots$)

• Apply the control action $u\left(t\right)$ based on the current policy $u\left(t\right)=-K_{i}x\left(t\right)$.
• At time $t+T$, collect and compute the dataset $(X,Y)$, which are defined in Equations (11) and (13).
• Update vector w by using the batch least-square method (Equation (10)).

Control Policy Update

• Transform vector w into the kernel matrix P using the Kronecker transformation.
• Update the policy by solving the LMI in Theorem 1.

Figure 1 shows the simplified diagram of the above algorithm. It is noted that the estimation of the value function is an on-policy learning since it updates $V^{\pi }\left(x\left(t\right)\right)$ using the V-value of the next state and the current policy’s action.

2.1.5. Simulation Setup

A simulation study of the proposed robust RL controller was conducted on a glucose kinetics model, which can be described by [20,21,22,23]:

$$\frac{d{D}_1\left(t\right)}{dt}=A_{G}D\left(t\right)-\frac{D_1\left(t\right)}{{\tau }_D},$$

(41)

$$\frac{d{D}_2\left(t\right)}{dt}=\frac{D_1\left(t\right)}{{\tau }_D}-\frac{D_2\left(t\right)}{{\tau }_D},$$

(42)

$$\frac{dg\left(t\right)}{dt}=-p_{1}g\left(t\right)-\chi \left(t\right)g\left(t\right)+\frac{D_2\left(t\right)}{{\tau }_D}+w\left(t\right)$$

(43)

and

$$\frac{d\chi \left(t\right)}{dt}=-p_2\chi \left(t\right)+p_{3}V\left(i\left(t\right)-i_b\left(t\right)\right).$$

(44)

In this model, parameter and variable descriptions can be found in

Table 1. Glucose kinetics model parameters.

Parameter	Description	Unit
$p_1$	Glucose effectiveness	${\min }^{-1}$
$p_2$	Insulin sensitivity	${\min }^{-1}$
$p_3$	Insulin rate of clearance	${\min }^{-1}$
$A_G$	Carbohydrate bioavailability	${\min }^{-1}$
${\tau }_D$	Glucose absorption constant	min
V	Plasma volume	mL
$i_b\left(t\right)$	Initial basal rate	$\mathsf{\mu }$IU/(mL·min)

and

Table 2. Variables of the glucose kinetics model.

Variable	Description	Unit
D	Amount of carbohydrate intake	mmol/min
$D_1$	Glucose in compartment 1	mmol
$D_2$	Glucose in compartment 2	mmol
$g\left(t\right)$	Plasma glucose concentration	mmol/L
$\chi \left(t\right)$	Interstitial insulin activity	${\min }^{-1}$
$i\left(t\right)$	Plasma insulin concentration	$\mathsf{\mu }\mathrm{IU}/\mathrm{mL}$

, respectively. The values of the parameters are selected based on [20,21]. Variable $w\left(t\right)$ in Equation (43) is the process noise. The measured blood glucose value is affected by a random noise $v\left(t\right)$:

$$\widehat{g}\left(t\right)=g\left(t\right)+v\left(t\right).$$

(45)

The inputs of the model are the amount of carbohydrate intake D and the insulin concentration i. The value of $i\left(t\right)-i_b\left(t\right)$ must be non-negative:

$$i\left(t\right)-i_b\left(t\right)\ge 0.$$

(46)

3. Results and Discussion

In order to evaluate the performance of the robust RL controller, we implemented the controller on the glucose kinetics model as described in the previous section under a daily scenario of patients with Type 1 diabetes. In order to make the scenario realistic, three different levels of uncertainties were used in the model. Uncertainties include process noise ($w\left(t\right)$) and measurement noise ($v\left(t\right)$). It is assumed that the noises are Gaussian distributions with standard deviations for each case as shown in

Table 3. Standard deviations of process and measurement noises.

Uncertainty Case	Process Noise ($\mathit{w}\left(\mathit{t}\right)$)	Measurement Noise ($\mathit{v}\left(\mathit{t}\right)$)
1	0	0
2	0	0.002
3	0.1	0.1
4	0.1	1

.

3.1. Without Meal Intake

This part describes the simulation results during the fasting period (without meal intake). The purpose of the simulation is to compare the performances of the robust RL algorithm with the conventional optimal RL algorithm [24] in the nominal condition (uncertainty case 1). The initial blood glucose for both scenario was set at 290 mg/dL and the target blood glucose is 90 mg/dL. The initial policy at the beginning of the simulation was chosen as follows:

$$u\left(t\right)=-K_{0}x\left(t\right)=-0.27g\left(t\right)+266.00\chi \left(t\right).$$

(47)

Figure 2 shows the comparison in blood glucose level between the robust RL and the optimal RL algorithm in the nominal condition. From the results, it can be seen that the robust RL successfully reduces the blood glucose level while the optimal RL becomes unstable when the blood glucose approaches the desired value. The instablity of the optimal RL in this case can be explained by the nonlinearity of the system (due to the coupling term $\chi \left(t\right)g\left(t\right)$ in Equation (43)), the saturation of the insulin concentration (Equation (46)), and the lack of perturbed data when the blood glucose approaches the steady-state value. The insulin concentration during the simulation can be found in Figure 3. In this figure, the dotted blue line indicates the unstable insulin profile.

Figure 4 shows the blood glucose responses from the robust RL in different uncertain conditions without meal intake. The results show similar and stable responses in all the uncertain conditions with settling time to the desired blood glucose level of approximately 45 min. The insulin concentration and the update of controller gains can be found in Figure 5 and Figure 6.

3.2. With Meal Intake

In this part, the performance of the robust RL controller was tested under conditions for which the system is subjected to meal intakes with the carbohydrate profile as shown in Figure 7.

During the simulation period with meal intakes, blood glucose responses throughout the day of the robust RL control systems under four uncertain cases are shown in Figure 8. The insulin concentration during the process can also be found in Figure 9. The results show that the controller provides the most aggressive action under case 1 (no uncertainty) and the least aggressive action under case 4 (with highest level of measurement and process noises). This leads to the largest and smallest reduction of postprandial blood glucose in case 1 and case 4, respectively. Most importantly, the robust RL algorithm kept the system in stable condition and there is no hypoglycemia event during the simulation for all four cases under different level of uncertainties.

4. Conclusions

The paper proposes a robust RL algorithm for dynamical systems with uncertainties. The uncertainties can be approximated by the critic and represented in the value function. LMI techniques were used to improve the controller gain. The algorithm was simulated on a blood glucose model for patients with Type 1 diabetes. The objective of the simulation is to control and maintain a healthy blood glucose level. The comparison between the robust RL algorithm and the optimal RL algorithm shows a significant improvement in the robustness of the proposed algorithm. Simulation results show that the algorithm successfully regulated the blood glucose and kept the system stable under different levels of uncertainty.