Rational Neural Network Controllers
Neural networks have shown great success in many machine learning related
tasks, due to their ability to act as general function approximators. Recent
work has demonstrated the effectiveness of neural networks in control systems
(known as neural feedback loops), most notably by using a neural network as a
controller. However, one of the big challenges of this approach is that neural
networks have been shown to be sensitive to adversarial attacks. This means
that, unless they are designed properly, they are not an ideal candidate for
controllers due to issues with robustness and uncertainty, which are pivotal
aspects of control systems. There has been initial work on robustness to both
analyse and design dynamical systems with neural network controllers. However,
one prominent issue with these methods is that they use existing neural network
architectures tailored for traditional machine learning tasks. These structures
may not be appropriate for neural network controllers and it is important to
consider alternative architectures. This paper considers rational neural
networks and presents novel rational activation functions, which can be used
effectively in robustness problems for neural feedback loops. Rational
activation functions are replaced by a general rational neural network
structure, which is convex in the neural network's parameters. A method is
proposed to recover a stabilising controller from a Sum of Squares feasibility
test. This approach is then applied to a refined rational neural network which
is more compatible with Sum of Squares programming. Numerical examples show
that this method can successfully recover stabilising rational neural network
controllers for neural feedback loops with non-linear plants with noise and
parametric uncertainty.
Matthew Newton, Antonis Papachristodoulou
# Rational Neural Network Controllers

###### Abstract

Neural networks have shown great success in many machine learning related tasks, due to their ability to act as general function approximators. Recent work has demonstrated the effectiveness of neural networks in control systems (known as neural feedback loops), most notably by using a neural network as a controller. However, one of the big challenges of this approach is that neural networks have been shown to be sensitive to adversarial attacks. This means that, unless they are designed properly, they are not an ideal candidate for controllers due to issues with robustness and uncertainty, which are pivotal aspects of control systems. There has been initial work on robustness to both analyse and design dynamical systems with neural network controllers. However, one prominent issue with these methods is that they use existing neural network architectures tailored for traditional machine learning tasks. These structures may not be appropriate for neural network controllers and it is important to consider alternative architectures. This paper considers rational neural networks and presents novel rational activation functions, which can be used effectively in robustness problems for neural feedback loops. Rational activation functions are replaced by a general rational neural network structure, which is convex in the neural network's parameters. A method is proposed to recover a stabilising controller from a Sum of Squares feasibility test. This approach is then applied to a refined rational neural network which is more compatible with Sum of Squares programming. Numerical examples show that this method can successfully recover stabilising rational neural network controllers for neural feedback loops with non-linear plants with noise and parametric uncertainty.

## 1 Introduction

Neural networks (NNs) have shown to be highly effective in numerous machine learning tasks. Examples of these include but are not limited to: image recognition, weather prediction, natural language processing, autonomous vehicle technology, medical imaging and social media algorithms [1, 2, 3, 4, 5, 6, 7]. There have been numerous advancements that have contributed to their success such as the development of modern NN architectures [8, 9], the increase in computational power available [10, 11] and the availability of big data.

More recently, there has been an increased interest in using NNs in control systems. One reason for this is the emergence of the parallel field of reinforcement learning. By harnessing the power of NNs, deep reinforcement learning has been used to create a decision-making agent to greatly outperform humans in many complex tasks. Such examples include the board game Go [12] and the video game Dota 2 [13]. Despite their success, there are many significant issues with theses methods. These learnt policies can perform poorly when the learnt environment is different from the real environment [14, 15, 16, 17]. Additionally, bounds to quantify their safety do not sufficiently describe the performance of the algorithm and can be overly conservative [18]. However, with new advancements in robust control and the success of NNs in reinforcement learning, there is a strong motivation for work at their intersection.

We refer to control systems that contain NNs as controllers as neural feedback loops (NFLs). Most research completed in this area has addressed the robustness analysis of NFLs, where the NN's parameters are given and the task is to quantify the system's safety or robustness properties. However, the more challenging task is to obtain the NN controller's parameters, whilst enforcing robustness guarantees. NNs are useful since they can be used as general
function approximators [19, 20, 21, 22], but contain a large number of parameters. This means that optimising over all of the parameters is often computationally expensive.

One method to design the NN controller is by learning an expert control law using input-output data and then checking the robustness guarantees using analysis methods, [23, 24, 25]. However, this may lead to poor robustness certificates because no relevant objective is being optimised while training the NN. It is also possible to use reinforcement learning methods. This often involves training the controller by simulating the system trajectories and then updating the controller's parameters subject to maximising a reward function. However, there are significant challenges with this process; it is very computationally intensive, requiring a large amount of hyperparameter tuning and sometimes leads to undesirable behaviour [26]. The parameterised NN controller can then be analysed to obtain robustness certificates in a defined operating region, however these guarantees can be poor. It can be difficult for these approaches to outperform traditional control laws. Furthermore, adding additional non-linearities into the model through the NN's activation functions can increase the complexity of the closed-loop system.

Despite these drawbacks, recent methods have focused on addressing these issues by designing NN controllers whilst ensuring robustness guarantees in the process. Methods that focus on developing reinforcement learning algorithms such as [27, 28] are able to create an NN policy which can be combined with robust control guarantees. However, these approaches still suffer from other reinforcement learning issues such as requiring significant computational time and hyperparameter tuning. These drawbacks can be mitigated by instead trying to obtain an NN controller by learning from an optimal model predictive control law. This allows the NN controller to be trained offline and when implemented it can be significantly less computationally expensive than the full model predictive control law. To achieve this an SDP framework that incorporates integral quadratic constraints is presented in [29]. An iterative algorithm is used to alternate between optimising the NN parameters to fit the control law and maximising the region of attraction. However, these approaches require a known model predictive control law to optimise over the controller's input-output data, which may not be easy to obtain. This approach relies on a loop transformation and Schur complement to ensure the optimisation problem is convex. Similar approaches have been used for different NN architectures such as recurrent NNs [30] and recurrent equilibrium network controllers [31]. These methods also include a projected policy gradient algorithm and reinforcement learning to synthesise the controller, instead of requiring an expert control law.

### Our Contribution

One shortfall of many recent approaches such as [29, 30, 31] is that large convex approximations must be made, as the original problem is non-convex. For example, the ReLU and tanh activation functions must be sector bounded in the region \([0,1]\) with no constraints obtained from pre-processing bounds such as Interval Bound Propagation [32]. These sector constraints are shown in Figure 1. This can lead to the acquired robustness certificates being conservative. Another issue with these methods is that they rely on iterative approaches, where the algorithm alternates between optimising the performance of the controller and the stability guarantees. This can make the problem computationally expensive and intractable for large scale systems. They also require that the plant model is linear subject to sector bounded non-linearities.

A root cause of the problems with recent approaches is that they use traditional NN architectures with ReLU, sigmoid and tanh activation functions. These structures have shown to be very effective in many machine learning tasks, however they may not be preferable for use in control systems. By considering a class of NNs that are better aligned with control system techniques, these large convex approximations could be mitigated. Since Sum of Squares (SOS) programming

Figure 1: Sector constraints that are used for common activation functions when designing neural network controllers.
is effective for systems with polynomial and rational functions, we investigate what happens when the NN is built upon them.

* Motivated by the success of rational neural networks in machine learning tasks, we propose novel rational activation functions to approximate the traditional sigmoid and tanh activation functions. We then show their effectiveness in analysing neural feedback loops that contain these rational activation functions.
* We note that we could be falling short of the potential expressivity of the neural network by using fixed activation functions. Therefore we consider a general neural network structure which is built upon equations that are convex in the design parameters. We show that a neural network of this form can be expressed in a similar way to that of a feed forward neural network with rational activation functions.
* We consider the Lyapunov stability criteria for constrained dynamics systems. We then propose a novel convex Sum of Squares procedure to recover a stabilising controller for a non-linear system through solving a Sum of Squares feasibility test.
* This procedure is extended to the generalised rational neural network architecture, to recover a stabilising rational neural network in a convex way. We then adapt the rational neural network architecture to make it more compatible with Sum of Squares programming and allow the stabilising controller to be recovered.
* We show for numerous examples that our proposed procedure and neural network architecture are able to effectively recover stabilising controllers for unstable systems with non-linear plants, noise and parametric uncertainty.

## 2 Preliminaries

### Sum of Squares Programming and The Positivstellensatz

In this section we outline how to formulate and solve SOS optimisation problems by solving an equivalent Semidefinite Program (SDP). For more details on SOS programming the reader is referred to [33, 34, 35, 36].

The fundamental idea behind SOS is to replace a polynomial positivity condition with a condition that enforces that the polynomial is a Sum of Squares.

**Definition 1**.: _A polynomial \(p(x)\) is said to be a Sum of Squares (SOS) polynomial if it can be expressed as_

\[p(x)=\sum_{i=1}r_{i}^{2}(x),\]

_where \(r_{i}(x)\in\mathbb{R}[x]\). We denote the set of polynomials that admit this decomposition by \(\Sigma[x]\) and say '\(p(x)\) is SOS'. Note that \(\mathbb{R}[x_{1},\ldots,x_{n}]\) is defined as the set of polynomials in \(x_{1},\ldots,x_{n}\) with real coefficients and we denote \(x=(x_{1},\ldots,x_{n})\) for simplicity._

**Definition 2**.: _A monomial in \(x=(x_{1},\ldots x_{n})\) is_

\[x^{\beta}=x_{1}^{\beta_{1}}x_{2}^{\beta_{2}}\ldots x_{n}^{\beta_{n}},\]

_where the exponent and degree are denoted as \(\beta=(\beta_{1},\ldots,\beta_{n})\in\mathbb{N}^{n}\) and \(|\beta|=\beta_{1}+\cdots+\beta_{n}\) respectively._

The column vector of monomials with only certain exponents is expressed as \(x^{\mathbb{B}}=(x^{\beta})_{\beta\in\mathbb{B}}\), where \(\mathbb{B}\subset\mathbb{N}^{n}\) is the set of exponents that are used in the monomials. The summation operation on \(\mathbb{B}\) is defined as

\[\mathbb{B}+\mathbb{B}:=\{\beta+\gamma\;:\;\beta,\gamma\in\mathbb{B}\}.\]

A polynomial \(p(x)\) can be written as a linear combination of a set of monomials \(x^{\mathbb{B}}\) for a set of coefficients \(p_{\beta}\in\mathbb{R}\) such that

\[p=\sum_{\beta\in\mathbb{N}^{n}_{d}}p_{\beta}x^{\beta},\]

where \(\mathbb{N}^{n}_{d}=\{\beta\in\mathbb{N}^{n}\;:\;|\beta|\leq d\}\) is the set of all \(n\)-variate exponents of degree \(d\) or less.

**Theorem 1**.: _A polynomial \(p(x)\) is SOS if and only if it can be written in what is referred to as a Gram matrix representation such that_

\[p(x)=(x^{\mathbb{B}})^{T}Qx^{\mathbb{B}},\]

_where \(Q\in\mathbb{S}^{|\mathbb{B}|}_{+}\) is a positive semidefinite matrix._
The Gram matrix representation can be rewritten as

\[(x^{\mathbb{B}})^{T}Qx^{\mathbb{B}}=\langle Q,x^{\mathbb{B}}(x^{\mathbb{B}})^{T} \rangle=\sum_{\alpha=\mathbb{B}+\mathbb{B}}\langle Q,A_{\alpha}\rangle x^{\alpha},\]

where the symmetric binary matrix \(A_{\alpha}\in\mathbb{S}^{|\mathbb{B}|}\) for each exponent \(\alpha\in\mathbb{B}+\mathbb{B}\) is

\[[A_{\alpha}]_{\beta,\gamma}:=\begin{cases}1,\;\beta+\gamma=\alpha,\\ 0,\;\mathrm{otherwise}.\end{cases}\]

The following is therefore true

\[p(x)\in\Sigma[x]\Leftrightarrow\exists Q\in\mathbb{S}^{|\mathbb{B}|}_{+}\; \text{where}\;\langle Q,A_{\alpha}\rangle=p_{\alpha},\;\forall\alpha\in\mathbb{ B}+\mathbb{B}.\]

This means that checking the SOS condition is equivalent to solving an SDP, which can be achieved with parsers such as SOSTOOLS [37] in MATLAB.

A central theorem of real algebraic geometry is known as the Positivstellensatz (Psatz) [38], which will now be briefly outlined. The Psatz provides an equivalent relation between an algebraic condition and the emptiness of a semi-algebraic set.

We express the semi-algebraic set with notation

\[S=\{x\in\mathbb{R}^{n}\mid g_{i}(x)\geq 0,\;h_{j}(x)=0,\;\forall\;i=1,\ldots,q_ {1},\;j=1,\ldots,q_{2}\}, \tag{1}\]

where \(g_{i}\) and \(h_{j}\) are polynomial functions.

**Definition 3**.: _Given \(f_{1},\ldots,f_{r}\in\mathbb{R}[x]\), the multiplicative monoid generated by the \(f_{k}\)'s is the set of all finite products of \(f_{k}\)'s, including 1 (i.e. the empty product). It is denoted as \(\mathcal{M}(f_{1},\ldots,f_{r})\)._

**Definition 4**.: _Given \(g_{1},\ldots,g_{q_{1}}\in\mathbb{R}[x]\), the cone generated by the \(g_{i}\)'s is_

\[\mathrm{cone}\{g_{1},\ldots,g_{q_{1}}\}=\Bigg{\{}s_{0}+\sum_{i=1}^{q_{1}}s_{i} G_{i}\mid s_{i}\in\Sigma[x],\;G_{i}\in\mathcal{M}(g_{1},\ldots,g_{q_{1}}) \Bigg{\}}. \tag{2}\]

**Definition 5**.: _Given \(h_{1},\ldots,h_{q_{2}}\in\mathbb{R}[x]\), the ideal generated by the \(h_{k}\)'s is_

\[\mathrm{ideal}\{h_{1},\ldots,h_{q_{2}}\}=\Bigg{\{}\sum_{j=1}^{q_{2}}t_{j}h_{j} \mid t_{j}\in\mathbb{R}[x]\Bigg{\}}. \tag{3}\]

**Theorem 2**.: _(Positivstellensatz, [38]) Given the semi-algebraic set \(S\) in (1), the following are equivalent:_

1. _The set_ \(S\) _is empty._
2. _There exist_ \(s_{i}\in\Sigma[x]\) _in (_2_) and_ \(t_{j}\in\mathbb{R}[x]\) _in (_3_) such that_ \[\mathrm{cone}\{g_{1},\ldots,g_{q_{1}}\}+\mathrm{ideal}\{h_{1},\ldots,h_{q_{2}} \}=0.\]

Based on Theorem 2 one can create a convex test for computational purposes by using a representation of the function \(p(x)\).

**Proposition 1**.: _Consider the set \(S\) in (1). If_

\[p=1+\sum_{j=1}^{q_{2}}t_{j}h_{j}+s_{0}+\sum_{i=1}^{q_{1}}s_{i}g_{i}, \tag{4}\]

_where \(s_{i}\in\Sigma[x]\) and \(t_{j}\in\mathbb{R}[x]\), then \(p(x)>0,\;\forall x\in S\)._

To test if the polynomial \(p(x)\geq 0\), \(\forall x\in S\) using SOS programming we can rewrite (4) as

\[p+\sum_{j=1}^{q_{2}}t_{j}h_{j}-\sum_{i=1}^{q_{1}}s_{i}g_{i}\in\Sigma[x],\]

where \(s_{i}\in\Sigma[x]\) and \(t_{j}\in\mathbb{R}[x]\). By selecting higher degree multipliers \(s_{i}\), \(t_{j}\) etc. we can obtain a series of set emptiness tests with increasing complexity and non-decreasing accuracy.
### Stability of Neural Feedback Loops using Sum of Squares

In this section we outline the methods proposed in [25], which presents a method to determine if the equilibrium a closed loop NFL is stable and then also use this to compute an inner approximation of the region of attraction.

Consider a continuous-time system

\[\dot{z}(t)=f(z(t),u(t)), \tag{5}\]

where \(f\) is the plant model \(z(t)\in\mathbb{R}^{n_{z}}\) and \(u(t)\in\mathbb{R}^{n_{u}}\) are the system states and inputs respectively. \(n_{z}\) and \(n_{u}\) are the number of system states and inputs respectively. This system is a continuous time NFL if the controller \(u\) is an NN. Consider a state feedback controller \(u(t)=\pi(z(t)):\mathbb{R}^{n_{z}}\rightarrow\mathbb{R}^{n_{u}}\) as a feed-forward fully connected NN such that

\[\begin{split} x^{0}(t)&=z(t),\\ v^{k}(t)&=W^{k}x^{k}(t)+b^{k},\;\mathrm{for}\;k=0, \ldots,\ell-1,\\ x^{k+1}(t)&=\phi(v^{k}(t)),\,\mathrm{for}\;k=0, \ldots,\ell-1,\\ \pi(z(t))&=W^{\ell}x^{\ell}(t)+b^{\ell},\end{split} \tag{6}\]

where \(W^{k}\in\mathbb{R}^{n_{k+1}\times n_{k}}\), \(b^{k}\in\mathbb{R}^{n_{k+1}}\) are the weights matrix and biases of the \((k+1)^{\text{th}}\) layer respectively and \(z(t)=x^{0}(t)\in\mathbb{R}^{n_{z}}\) is the input into the NN. The activation function \(\phi\) is applied element-wise to the \(v^{k}(t)\) terms. The number of neurons in the \(k^{\text{th}}\) layer is denoted by \(n_{k}\).

As shown in [25], an NFL can be viewed as a dynamical system with equality and inequality constraints that arise from the input-output description of the NN. These constraints can be described as the semi-algebraic set in (1). Further constraints can be added to the semi-algebraic set when only local asymptotic stability is being verified. The region

\[D^{z}=\left\{z\in\mathbb{R}^{n_{z}}\,|\,d_{k}(z)\geq 0,\;k=1,\ldots,n_{d}\right\}, \tag{7}\]

is considered where the stability conditions will need to be satisfied.

**Proposition 2**.: _([25]) Consider System (5) in feedback with an NN controller given by (6). Suppose the input-output properties of the NN are described by (1), and consider the region given by (7). Suppose there exists a polynomial function \(V(z)\) satisfying the following conditions_

\[\begin{split} V(z)-\rho(z)&\in\Sigma[z],\\ \rho(z)&>0,\\ -\frac{\partial V}{\partial z}(z)f(z,\pi(z))-\sum_{k=1}^{n_{d}} p_{k}(X)d_{k}(z)-\sum_{j=1}^{q_{2}}t_{j}(X)h_{j}(x)-\sum_{i=1}^{q_{1}}s_{i}(X)g _{i}(x)\in\Sigma[X],\\ p_{k}(X)&\in\Sigma[X],\;\forall k=1,\ldots,n_{d},\\ s_{i}(X)&\in\Sigma[X],\;\forall i=1,\ldots,q_{1},\\ t_{j}(X)&\in\mathbb{R}[X],\;\forall j=1,\ldots,q_{2}, \end{split} \tag{8}\]

_where \(X\) is a vector of all the system and NN states, i.e. \(X=(x,z)\). Then the equilibrium of the NFL is stable._

If a Lyapunov function is constructed using Proposition 2 then the region of attraction can be approximated. To achieve this we find the largest level set of the Lyapunov function \(V(z)\) that is contained within the region that the Lyapunov conditions are satisfied. This can be cast as an SOS program. Consider a Lyapunov function \(V(z)\), if the SOS optimisation problem

\[\begin{split}|z|^{k}(V(z)-\gamma)+p(z)d(z)&\in \Sigma[z],\\ p(z)&\in\Sigma[z],\end{split} \tag{9}\]

is feasible, where \(\gamma\) is a variable to be maximised and \(k\) is a positive integer, then \(V(z)\leq\gamma\) is an estimate of the region of attraction.

## 3 Rational Neural Networks

There has been little prior work investigating NNs with rational expressions for control. [20] showed the expressive power of using rational activation functions and how they can be used to approximate commonly used activation functions such as ReLU. The effectiveness of rational NNs in machine learning tasks was shown in [39]. Polynomial
activation functions have also been used in previous work [40, 41, 42, 43, 44], however they have not seen much investigation within the machine learning research community due to vanishing and exploding gradients that can be exhibited when used with the backpropagation algorithm [45]. Another class of NNs that has seen recent interest are quadratic NNs. These can take many different forms, which are summarised in [46]. Quadratic NNs have been shown to be beneficial in many aspects [47, 22, 48], such as being general universal function approximators. Recent work from [49] has explored the use of a two-layer NN controller that has quadratic activation functions and how a convex formulation can be achieved using this structure.

One of the benefits of using polynomial or rational activation functions in the NN is that they can contain trainable parameters, which can improve the performance of the NN and reduce the number of neurons in the network. For example, the rational activation function used in [39] is expressed as

\[\phi(x)=\frac{P(x)}{Q(x)}=\frac{\sum_{i=0}^{r_{P}}a_{i}x^{i}}{\sum_{j=0}^{r_{ O}}b_{j}x^{j}}, \tag{10}\]

where \(r_{P}=\deg(P(x))\), \(r_{Q}=\deg(Q(x))\) and \(a_{i}\) and \(b_{j}\) are the trainable parameters within the activation function. Another proposed activation function based on rational functions is the Pade activation unit, which has shown to be useful when applied to image classification [50]. This activation function is expressed as

\[\phi(x)=\frac{P(x)}{Q(x)}=\frac{\sum_{i=0}^{r_{P}}a_{i}x^{i}}{1+|\sum_{j=0}^{r _{O}}b_{j}x^{j}|}\]

and can be used to learn and approximate commonly used activation functions, whilst training the NN. The resulting NNs can have compact representations and can perform similarly to state-of-the-art NNs.

### Rational Approximation of Tanh Activation Function

Despite their success in some machine learning applications, rational activation functions have yet to be applied to NN controllers. Most methods to train NNs in control systems use traditional activation functions such as ReLU, sigmoid and tanh. In addition, NNs that are used in control systems are often significantly smaller than those used for machine learning tasks such as image classification. The reason behind this is that methods to train NNs for control often require the use of an SDP, meaning that having a large number of parameters in the NN makes the problem become intractable. However, the requirement for the NN to only contain a small number of neurons to learn a function that can sufficiently control a system has not been well explored. One justification for this is the low number of input and output dimensions that are required for small NFLs. However, if the number of system states increases, then larger networks may be required, which may be intractable to obtain with current methods.

For NNs to be effectively used as controllers, it would be of interest to investigate the use of alternative NN structures and activation functions that may give a sufficient level of expressivity in the network, whilst being easier to compute or ensure robustness guarantees. Motivated by this we propose a novel activation function defined by a simple rational expression. We name this function 'Rtanh' as it is an approximation of the tanh activation function

\[\tanh(x)=\phi(x)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}\]

and is defined as

\[\mathrm{Rtanh}(x)=\phi(x)=\frac{4x}{x^{2}+4}.\]

This function is shown in Figure 1(a) against the tanh function and the error between these functions is shown in Figure 1(b).

The tanh function can be over-approximated by sector constraints, as demonstrated in [29, 30, 31]. These constraints can then be used to analyse the NFL by creating an optimisation problem through an SDP or SOS programming. However, these bounds are very conservative as shown in Figure 0(b). A big advantage of the 'Rtanh' function is that it can be represented by a single equality constraint such that

\[\phi(x)(x^{2}+4)-4x=0. \tag{11}\]

If this activation function were to be used in an NN, then to test its robustness properties using the Psatz, the equality constraint (11) can be used directly. This is beneficial as no conservatism is introduced when abstracting the input-output properties of the NN using a semi-algebraic set.
To demonstrate that the Rtanh activation function is useful when analysing NFLs, we take an existing NFL that uses tanh activation functions. We consider the inverted pendulum from [24], which uses a five layer NN with five nodes in each layer and tanh activation functions. The dynamics of this system are expressed as

\[\ddot{\theta}(t)=\frac{mgl\sin\theta(t)-\mu\dot{\theta}(t)+\text{sat}(u(t))}{ml^ {2}}, \tag{12}\]

where \(\text{sat}(\cdot)\) is the saturation function. The system is discretised with time step \(\Delta t=0.2\) and is parameterised by \(m=0.15\text{kg}\), \(l=0.5\text{m}\), \(\mu=0.5\text{Nms/rad}\), \(g=10\text{m/s}^{2}\), \(u_{\text{max}}=1\) where \(u_{\text{max}}\) is the saturation limit in the saturation function. As shown in [51], we can use 'ReachSparsePastz', which is an SOS optimisation technique to approximate the reachable set at each time step. Using the tanh activation function as in the original NN controller, the reachable sets can be computed and are shown in Figure 3.

We then replace the tanh function with Rtanh and observe the system behaviour. We conduct the same reachability analysis as in [51] by computing the reachable sets for six time steps. In Figure 4 we can see that this activation function behaves similarly to that of the tanh function and that using the Psatz with sparse polynomial optimisation (ReachSparsePastz) gives very tight approximations of the reachable sets.

We also compute the region of attraction using the approach outlined in Section 2.2, refereed to as 'NNSSOSStability'. The region of attraction when using the tanh and Rtanh functions are shown in Figure 5 and Figure 6 respectively. We

Figure 3: Reachable sets for the inverted pendulum in (12) with the tanh activation function. The reachable sets are computed using ReachSparsePastz and are shown in black. The true reachable sets are shown in blue.

Figure 2: Plots showing a comparison between the \(\mathrm{Rtanh}(x)\) and \(\tanh(x)\) functions.
find that the region of attraction is increased significantly when using Rtanh over the tanh activation function. The areas of the region of attraction on the phase plane are 58 and 0.94 per square unit for the Rtanh and tanh activation functions respectively.

These results show that the Rtanh function can be useful in NNs and NFLs. However, this function is not compatible with most recent methods to learn NN controllers due to it being represented by an equality that is a polynomial of degree three. Indeed, most methods require sector inequality constraints and the Schur complement, which is used to

Figure 4: Reachable sets for the inverted pendulum in (12) where the tanh activation function is replaced with the Rtanh function. The reachable sets are computed using ReachSparsePastz and are shown in black. The true reachable sets are shown in blue.

Figure 5: Diagram showing the region of attraction for the inverted pendulum in (12) with tanh activation functions in the neural network controller. The region of attraction is computed using NNSSStability (black) with a fourth order Lyapunov function. The system trajectories are shown in blue.
make the problem convex. Therefore, this function cannot be used in those formulations. This requires the use of an alternative method to obtain an NN controller with these activation functions, which we will investigate later in this paper.

### Rational Approximation of Sigmoid Activation Function

It is also possible to approximate the sigmoid function as a rational function. We define the function

\[\mathrm{Rsig}(x)=\phi(x)=\frac{(x+4)^{2}}{2(x^{2}+16)}.\]

As shown in Figure 6(a), this is a good approximation to the sigmoid function. Figure 6(b) shows the error with the sigmoid function. Rsig can be expressed as the equality constraint

\[2(x^{2}+16)\phi(x)-(x+4)^{2}=0.\]

Figure 6: Diagram showing the region of attraction for the inverted pendulum in (12) where the tanh activation function is replaced with the Rtanh function. The region of attraction is computed using NNNOSSSstability (black) with a fourth order Lyapunov function. The system trajectories are shown in blue.

Figure 7: Plots showing a comparison between the \(\mathrm{Rsig}(x)\) and \(\mathrm{sig}(x)\) functions.
### Irrational Approximation of ReLU Activation Function

Approximating the ReLU function with a rational function is more difficult. Here we consider the class of irrational activation functions; to demonstrate how they could be used we propose the following function

\[\mathrm{IReLU}(x)=\phi(x)=\sqrt{x^{2}+1}+x-1.\]

This function and the error with the ReLU function are shown in Figure (a)a and Figure (b)b respectively. By making a substitution, the IReLU function can be expressed as a semi-algebraic set with two equality constraints and one inequality constraint

\[\phi(x)-y-x+1 =0,\] \[y^{2}-x^{2}-1 =0,\] \[y \geq 0.\]

### General Rational Neural Networks

The rational approximations of the tanh and sigmoid functions proposed in Section 3 can be used as activation functions in an NN. However, if we were to use these structures we could be missing the potential expressivity that can be obtained from a general class of rational functions. As in (10), rational activation functions can contain training parameters. Therefore, instead of considering a rigid structure for the activation function and fixing the parameters, we can consider a general rational expression that contains the NN parameters.

If we were to consider a feed-forward fully connected NN with predefined rational activation functions, then the NN can be written as

\[x^{0} =u, \tag{13}\] \[v^{k} =W^{k}x^{k}+b^{k},\;\mathrm{for}\;k=0,\dots,\ell-1,\] \[x^{k+1} =\frac{p(v^{k})}{q(v^{k})},\;\mathrm{for}\;k=0,\dots,\ell-1,\] \[\pi(u) =W^{\ell}x^{\ell}+b^{\ell},\]

where \(p(v^{k})\) and \(q(v^{k})\) are polynomial functions with specified coefficients. However, if we substitute the preactivation value \(v^{k}_{j}\) into the polynomial expressions, we obtain a rational expression in \(x^{k}_{j}\), where the coefficients are parameterised by the values of the weight matrices \(W^{k}\) and biases vector \(b^{k}\). Any coefficients in the rational activation function will be multiplied by the weights and bias terms. Therefore, we can instead consider an NN with no affine transformation and parameters only in the rational activation function. We can then write the rational NN as

\[x^{0} =u, \tag{14}\] \[x^{k+1}_{i} =\frac{p_{i}(x^{k})}{q_{i}(x^{k})},\;\mathrm{for}\;i=0,\dots,n_{k},\;\mathrm{for}\;k=0,\dots,\ell-1,\] \[\pi_{i}(u) =\frac{p_{i}(x^{\ell})}{q_{i}(x^{\ell})},\;\mathrm{for}\;i=0, \dots,n_{\ell},\]

Figure 8: Plots showing a comparison between the \(\mathrm{IReLU}(x)\) and \(\mathrm{ReLU}(x)\) functions.
where \(p_{i}(x^{k})\) and \(q_{i}(x^{k})\) are general polynomial functions which can be written as

\[p_{i}(x) =\sum_{\alpha\in\mathbb{N}_{q_{pi}}^{n}}\lambda_{\alpha}x^{\alpha},\] \[q_{i}(x) =\sum_{\beta\in\mathbb{N}_{q_{q_{i}}}^{n}}\gamma_{\beta}x^{\beta},\]

where \(d_{q_{i}}=\deg(q_{i})\), \(d_{p_{i}}=\deg(p_{i})\), \(\alpha\) and \(\beta\) are the exponents that are defined in Section 2.1 and \(\lambda_{\alpha}\) and \(\gamma_{\beta}\) are the coefficients of \(p_{i}(x)\) and \(q_{i}(x)\) respectively. To show the similarity between (13) and (14) we present the following simple example.

**Example 1**.: _Consider an NN with structure defined by (13) with two layers and two nodes in each layer. The rational activation functions with \(\deg(p)=\deg(q)=2\) can be written as_

\[p(v) =c_{1}v^{2}+c_{2}v+c_{3},\] \[q(v) =d_{1}v^{2}+d_{2}v+d_{3}.\]

_The first node in the second layer can be expressed as_

\[v_{1}^{1} =W_{1}^{1}x_{1}^{1}+W_{2}^{1}x_{2}^{1}+b_{1}^{1},\] \[x_{1}^{2} =\frac{p(v_{1}^{1})}{q(v_{1}^{1})}.\]

_We can substitute in the preactivation terms into the activation functions to obtain the polynomials_

\[p(v_{1}^{1}) =c_{1}(W_{1}^{1})^{2}(x_{1}^{1})^{2}+2c_{2}W_{1}^{1}W_{2}^{1}(x_{1 }^{1}x_{2}^{1})+c_{3}(W_{2}^{1})^{2}(x_{2}^{1})^{2}+\] \[(c_{1}W_{1}^{1}+2c_{2}W_{1}^{1}b_{1}^{1})(x_{1}^{1})+(c_{1}W_{2} ^{1}+2c_{2}W_{2}^{1}b_{1}^{1})(x_{2}^{1})+(c_{1}(b_{1}^{1})^{2}+c_{2}b_{1}^{1}+ c_{3}),\]

\[q(v_{1}^{1}) =d_{1}(W_{1}^{1})^{2}(x_{1}^{1})^{2}+2d_{2}W_{1}^{1}W_{2}^{1}(x_{1 }^{1}x_{2}^{1})+d_{3}(W_{2}^{1})^{2}(x_{2}^{1})^{2}+\] \[(d_{1}W_{1}^{1}+2d_{2}W_{1}^{1}b_{1}^{1})(x_{1}^{1})+(d_{1}W_{2} ^{1}+2d_{2}W_{2}^{1}b_{1}^{1})(x_{2}^{1})+(d_{1}(b_{1}^{1})^{2}+d_{2}b_{1}^{1} +d_{3}),\]

_which form the rational expression for \(x_{1}^{2}\). However, we can instead write the rational expression as_

\[x_{1}^{2}=\frac{\lambda_{1}(x_{1}^{1})^{2}+\lambda_{2}x_{1}^{1}x_{2}^{1}+ \lambda_{3}(x_{2}^{1})^{2}+\lambda_{4}x_{1}^{1}+\lambda_{5}x_{2}^{1}+\lambda_ {6}}{\gamma_{1}(x_{1}^{1})^{2}+\gamma_{2}x_{1}^{1}x_{2}^{1}+\gamma_{3}(x_{2}^ {1})^{2}+\gamma_{4}x_{1}^{1}+\gamma_{5}x_{2}^{1}+\gamma_{6}},\]

_which is in the form of (14). This can be generalised to larger NNs with higher degree polynomials in the rational functions. This approach can reduce the number of parameters in the NN and each polynomial is convex in the decision variables. This allows the parameters in the rational expression to be tuned directly instead of simultaneously tuning the weights, biases and rational activation parameters._

## 4 Recovering Stabilising Controllers using Sum of Squares

In this section, we propose a novel procedure to obtain a stabilising controller for a non-linear polynomial system using SOS programming. To do this we leverage the Psatz and exploit its structure to generate a feasibility test for a stabilising controller.

**Proposition 3**.: _Consider the polynomial \(P(x)\in\mathbb{R}[x]\) and partition \(x=[y,z]\) so that \(y\in\mathbb{R}^{n},z\in\mathbb{R}\). Consider the set_

\[S=\{y\in\mathbb{R}^{n},\;z\in\mathbb{R}\;|\;p_{i}(y,z)-q_{i}(y)z=0\;\;\forall\; i=1,\ldots,m\},\]

_where \(p_{i}(y,z)\in\mathbb{R}[x]\), \(q_{i}(z)\in\mathbb{R}[z]\). If_

\[P(x)-\sum_{i=1}^{m}(p_{i}(y,z)-q_{i}(y)z)\in\Sigma[x],\;q_{i}(y)\neq 0,\; \forall\;i=1,\ldots,m, \tag{15}\]

_then \(P(x)\geq 0\) on the set_

\[T=\bigg{\{}y\in\mathbb{R}^{n},\;z\in\mathbb{R}\;\bigg{|}\;\frac{p_{i}(y,z)}{q _{i}(y)}-z=0,\;q_{i}(y)\neq 0\;\;\forall\;i=1,\ldots,m\bigg{\}}.\]
Proof.: Consider the set \(S\), if

\[P(x)-\sum_{i=1}^{m}t_{i}(p_{i}(y,z)-q_{i}(y)z)\in\Sigma[x],\]

where \(t_{i}=1\), since \(t_{i}\in\mathbb{R}[x]\), then \(P(x)\geq 0\) on \(S\). If we include the condition \(q_{i}(y)\neq 0,\;\forall\;i=1,\ldots,m\), then we can rewrite (15) as

\[P(x)-\sum_{i=1}^{m}q_{i}(y)\bigg{(}\frac{p_{i}(y,z)}{q_{i}(y)}-z\bigg{)}\in \Sigma[x],\;q_{i}(y)\neq 0,\;\forall\;i=1,\ldots,m,\]

since \(q_{i}(y)\in\mathbb{R}[x]\), then \(P(x)\geq 0\) on \(T\). 

Proposition 3 considers the Psatz in a particular form to obtain a positivity certificate of a function over a set of rational functions. In essence, this formulation fixes the multipliers in the Psatz to unity and then searches over the polynomials to find a constraint set that satisfies the feasibility test. This is useful if we want to find a constraint set for the nonnegativity of a function, instead of testing that function over a known constraint set. This can be leveraged to find a controller that satisfies the Lyapunov stability conditions.

**Proposition 4**.: _Consider the non-linear system_

\[\dot{z} =f(z)+g(z)u, \tag{16}\] \[u =\frac{p(z)}{q(z)},\;q(z)\neq 0,\]

_where \(z\in\mathbb{R}^{n_{x}}\) are the system states, \(u\in\mathbb{R}^{n_{u}}\) is the controller input and \(f(z)\) and \(g(z)\) are polynomials. Suppose that there exists a Lyapunov function \(V(z)\) such that \(V(z)\) is positive definite in a neighbourhood of the origin and polynomials \(p(z),\;q(z)\) satisfying_

\[-\frac{\partial V}{\partial z}(f(z)+g(z)u)-(p(z)-q(z)u)\geq 0\;\forall z,u.\]

_Then the origin of the state space is a stable equilibrium of the system._

Proof.: We consider the stability of constrained dynamical systems as in [33]. Using the same argument as in Proposition 3, if

\[-\frac{\partial V}{\partial z}(f(z)+g(z)u),\]

is nonnegative on the set

\[\{z\in\mathbb{R}^{n_{z}},\;u\in\mathbb{R}^{n_{u}}\;|\;p(z)-q(z)u=0,\;q(z)\neq 0\},\]

then it is also nonnegative on the set

\[\bigg{\{}z\in\mathbb{R}^{n_{z}},\;u\in\mathbb{R}^{n_{u}}\;\bigg{|}\;\frac{p(z )}{q(z)}-u=0,\;q(z)\neq 0\bigg{\}},\]

which defines the controller in the closed-loop system. 

To compute the Lyapunov function and polynomials that define the rational controller in Proposition 4 we can formulate an SOS program.

**Proposition 5**.: _Consider the dynamical system (16) in Proposition 4. Suppose there exists polynomial functions \(V(z),\;p(z),\;q(z)\), a positive definite function \(\rho(z)\) such that_

\[V(z)-\rho(z) \in\Sigma[z],\] \[-\frac{\partial V}{\partial z}(f(z)+g(z)u)-(p(z)-q(z)u) \in\Sigma[X],\] \[p(z) \in\mathbb{R}[X],\] \[q(z) \in\mathbb{R}[X],\] \[q(z) \neq 0,\]

_where \(X=(z,u)\) is a vector of all of the states. Then the origin of the system is a stable equilibrium._

Theorem 5 allows us to reconstruct a stabilising controller from a feasibility test using SOS programming. However, the structure of the controller is limited as it is a simple rational function. We therefore extend this approach to a more expressive class of functions through rational NNs.
### Extension to Rational Neural Network Controllers

The technique outlined in the previous section can be used to recover a stabilising controller that is a rational function of the system states. However, we can expand this approach to consider an NN architecture that contains rational functions similar to the one proposed in (14). We consider a state feedback controller \(u(t)=\pi(z(t)):\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{u}}\) as a rational NN such that

\[\begin{split} x^{0}(t)&=z(t),\\ x_{i}^{k+1}(t)&=\frac{p_{i}^{k}(x^{k}(t))}{q_{i}^{ k}(x^{k}(t))},\,\text{for}\ i=1,\ldots,n_{k},\ k=0,\ldots,\ell-1,\\ u_{i}(t)=\pi_{i}(z(t))&=\frac{p_{i}^{\ell}(x^{ \ell}(t))}{q_{i}^{\ell}(x^{\ell}(t))},\,\text{for}\ i=1,\ldots,n_{z},\end{split} \tag{17}\]

where \(p_{i}^{k}(x^{k}(t))\), \(q_{i}^{k}(x^{k}(t))\) are the polynomials that form the rational expression associated with the \(i^{\text{th}}\) node in the \((k+1)^{\text{th}}\) layer. The number of neurons in the \(k^{\text{th}}\) layer is denoted by \(n_{k}\). We will drop the time dependence notation throughout the rest of this paper for simplicity.

We can apply Proposition 5 and the theory of constrained dynamical systems as in [33] due to the well-defined structure of this controller. The following proposition shows how Lyapunov stability over constrained dynamical systems can be used to recover a controller of this form.

**Proposition 6**.: _Consider the non-linear system_

\[\begin{split}\dot{z}&=f(z)+g(z)u,\\ u&=\pi(z),\end{split}\]

_where \(z\in\mathbb{R}^{n_{x}}\) are the system states, \(u\in\mathbb{R}^{n_{u}}\) is the controller input and \(f(z)\) and \(g(z)\) are polynomials. Consider the controller structure \(\pi(z)\) defined in (17) and the region given by (7). Suppose there exist polynomial functions \(V(z)\), \(p_{i}^{k}(x^{k})\ \forall i=1,\ldots,n_{k+1},\ k=0,\ldots,\ell\), and \(q_{i}^{k}(x^{k})\ \forall i=1,\ldots,n_{k+1},\ k=0,\ldots,\ell\) satisfying the following conditions_

\[\begin{split} V(z)-\rho(z)&\in\Sigma[z],\\ \rho(z)&>0,\\ -\frac{\partial V}{\partial z}(z)(f(z)+g(z)u)&-\sum_ {k=1}^{n_{d}}s_{k}(X)d_{k}(z)\ldots\\ &-\sum_{k=1}^{\ell}\sum_{i=1}^{n_{k}}\left(p_{i}^{k}(x^{k})-q_{i} ^{k}(x^{k})x_{i}^{k+1}\right)\ldots\\ &-\sum_{i=1}^{n_{u}}\left(p_{i}^{\ell}(x^{\ell})-q_{i}^{\ell}(x^{ \ell})u_{i}\right)\in\Sigma[X],\\ s_{k}(X)&\in\Sigma[X],\,\forall k=1,\ldots,n_{d},\\ q_{i}^{k}(x^{k})&\neq 0,\,\forall\,i=1,\ldots,n_{k}, \ k=0,\ldots,\ell,\end{split}\]

_where \(X\) is a vector of all the system and NN states, i.e. \(X=(x,u,z)\). Then the equilibrium of the system is stable._

The above proposition can generate a feasible SOS program, however the rational NN controller may be difficult to recover. This is because the coefficients in the \(q_{i}^{k}(x^{k})\) terms will be set to very small values by the SOS program, due to each term cancelling with the adjacent layers. We therefore propose an alternative rational NN controller structure in the following section to mitigate this issue.

### Refined Rational Neural Network Controller

The Lyapunov condition in Proposition 6 may result in numerical issues when solving the SOS program due to the structure of the constraints, making the rational NN controller difficult to recover. To overcome this issue, we enrich the
NN structure by considering a state feedback controller \(u=\pi(z):\mathbb{R}^{n_{z}}\rightarrow\mathbb{R}^{n_{u}}\) as a rational NN such that

\[\begin{split} y^{0}&=z,\\ x_{i}^{k+1}&=\frac{\sum_{j=1}^{n_{k}}p_{i,j}^{k}(y^{k}) y_{j}^{k}}{q_{i}^{k}(z)},\,\text{for}\,i=1,\dots,n_{k+1},\,k=0,\dots,\ell-1,\\ y_{i}^{k+1}&=x_{i}^{k+1}+1,\,\text{for}\,i=1,\dots,n _{k},\,k=0,\dots,\ell-1,\\ u_{i}=\pi_{i}(z)&=\frac{\sum_{j=1}^{n_{\ell}}p_{i,j} ^{\ell}(y^{\ell})y_{j}^{\ell}}{q_{i}^{\ell}(z)}\Big{(}\sum_{m=1}^{n_{k}}z_{m}^{ 2}\Big{)},\,\text{for}\,i=1,\dots,n_{z},\end{split} \tag{18}\]

where \(p_{i,j}^{k}(x^{k})\), \(q_{i}^{k}(z)\) are the \(j^{\text{th}}\) polynomials that form the rational expression associated with the \(i^{\text{th}}\) node in the \((k+1)^{\text{th}}\) layer.

Each \(x_{j}^{k}\) node in the NN contains a rational activation function and each \(y_{j}^{k}\) node is equal to the \(x_{j}^{k}\) term with the addition of a bias term which we set to unity. The \(y_{j}^{k}\) term that appears in the numerator of the rational activation function is to ensure that all terms in the \(x_{i}^{k+1}\) node are a function of all of the nodes in the \(k^{\text{th}}\) layer and to impose more structure on the NN. The denominator \(q_{i}^{k}(z)\) is a function of the system states to ensure that all nodes in the network are tied to the system states and not just the nodes in the previous layer. This will ensure that the SOS program does not set the coefficients in the \(q_{i}^{k}(z)\) terms to very small values. The final layer contains a multiplier term which is a sum of all of the system states \(\sum_{m=1}^{n_{k}}z_{m}^{2}\) to ensure that the controller input goes to zero at the origin.

We can then adapt Proposition 6 for this modified rational NN structure, to obtain a controller that can be recoverable from the feasibility test.

**Proposition 7**.: _Consider the non-linear system_

\[\begin{split}\dot{z}&=f(z)+g(z)u,\\ u&=\pi(z),\end{split}\]

_where \(z\in\mathbb{R}^{n_{z}}\) are the system states, \(u\in\mathbb{R}^{n_{u}}\) is the controller input and \(f(z)\) and \(g(z)\) are polynomials. Consider the controller structure \(\pi(z)\) defined in (18) and the region given by (7). Suppose there exist polynomial functions \(V(z)\), \(p_{i,j}^{k}(x^{k})\;\forall i=1,\dots,n_{k+1},\;j=1,\dots,n_{k},\;k=0,\dots,\ell\) and \(q_{i}^{k}(z)\;\forall i=1,\dots,n_{k+1},\;k=0,\dots,\ell\) satisfying the following conditions_

\[\begin{split} V(z)-\rho(z)&\in\Sigma[z],\\ \rho(z)&>0,\\ -\frac{\partial V}{\partial z}(z)(f(z)+g(z)u)&- \sum_{k=1}^{n_{d}}s_{k}(X)d_{k}(z)\dots\\ &-\sum_{k=1}^{\ell}\sum_{i=1}^{n_{k}}\left(\Big{(}\sum_{j=1}^{n_ {k-1}}p_{i,j}^{k}(y^{k})y_{j}^{k}\Big{)}-q_{i}^{k}(z)x_{i}^{k+1}\right)\dots \\ &-\sum_{k=1}^{\ell}\sum_{i=1}^{n_{k}}t_{i,k}(X)(y_{i}^{k}-x_{i}^{k }-1)\dots\\ &-\sum_{i=1}^{n_{u}}\bigg{(}\sum_{j=1}^{n_{\ell}}\bigg{(}p_{i,j} ^{\ell}(y^{\ell})y_{j}^{\ell}\bigg{(}\sum_{m=1}^{n_{z}}z_{m}^{2}\bigg{)} \bigg{)}-q_{i}^{\ell}(z)u_{i}\bigg{)}\in\Sigma[X],\\ s_{k}(X)&\in\Sigma[X],\;\forall k=1,\dots,n_{d},\\ t_{i,k}(X)&\in\mathbb{R}[X],\;\forall k=1,\dots,\ell,\;i=1,\dots,n_{k},\\ q_{i}^{k}(x^{k})&\neq 0,\;\forall\;i=1,\dots,n_{k},\;k=0, \dots,\ell,\end{split}\]

_where \(X\) is a vector of all the system and NN states, i.e. \(X=(u,x,y,z)\). Then the equilibrium of the system is stable._

The SOS program in the above proposition is convex in the rational NN parameters and can hence be solved using SOS programming. Note that saturation and any uncertainty and robustness conditions can be incorporated in the same way that is demonstrated in [25]. Proposition 7 presents a method to obtain a stabilising NN controller for a non-linear polynomial system in a convex way by solving one SOS optimisation problem. As described in Section 1, other recent approaches such as [27, 28, 29, 30, 31] rely on iterative algorithms or reinforcement learning formulations that are significantly more expensive to compute.
## 5 Numerical Examples

We demonstrate how the approach in Proposition 7 can be used to recover stabilising rational NN controllers through numerical examples. These examples were run on a four-core Intel Xeon processor @3.50GHz with 16GB of RAM. The SOS programs were implemented using MATLAB and SOSTOOLS to parse the SOS constraints into an SDP, which is solved using MOSEK [52].

### One Dimensional System

To show that this method is able to obtain a stabilising rational NN controller, we consider a very simple one dimensional linear system of the form

\[\dot{z}=z+u.\]

This system is unstable without a feedback controller. To recover a stabilising controller we set the size of the rational NN to be a small two layer network with a single node in each layer. The equations of the controller can be written as

\[x_{1} =\frac{\lambda_{1,1}z^{4}+\lambda_{1,2}z^{2}}{\gamma_{1,1}z^{2}+ \gamma_{1,2}},\] \[y_{1} =x_{1}+1,\] \[u =\frac{(\lambda_{2,1}y_{1}^{2}+\lambda_{2,2}y_{1})z}{\gamma_{2,1} z^{2}+\gamma_{2,2}},\]

where \(\gamma_{1,1}\geq 0\), \(\gamma_{1,2}>0\), \(\gamma_{2,1}\geq 0\), \(\gamma_{2,2}>0\). We also add saturation to the controller such that \(-10\leq u\leq 10\) and we enforce the local region of the state space to be \(-10\leq z\leq 10\). The SOS program is able to recover a feasible controller which can stabilise the system to the zero equilibrium. The state trajectory and controller input for this NFL over time are shown in Figure 9 and 10 respectively.

### Three Dimensional Non-linear System

We now consider a three dimensional non-linear system system given by

\[\dot{z}_{1} =-z_{1}+z_{2}-z_{3},\] \[\dot{z}_{2} =-z_{1}(z_{3}+1)-z_{2},\] \[\dot{z}_{3} =-z_{1}+u,\]

Figure 9: Plot showing the trajectories of the system state in Section 5.1 over time.
and attempt to find a stabilising rational NN controller for the system. We set the NN to have two layers and three nodes in each layer, with saturation \(-10\leq u\leq 10\). The system states are defined to operate in the region

\[1^{2}-z_{1}^{2}-z_{2}^{2}-z_{3}^{2}\geq 0.\]

Each polynomial in the rational NN is assigned to contain zeroth to fourth order terms. We define a quartic Lyapunov function and second order polynomials for the \(s_{k}\) and \(t_{i,k}\) terms. The trajectories for this system are shown in Figure 11, which shows that the controller can successfully stabilise the system.

Figure 11: Plot showing the trajectories of the system states for the example in Section 5.2.

Figure 10: Plot showing the controller input of the system in Section 5.1 over time.
### Non-linear Inverted Pendulum

We now consider the inverted pendulum proposed in [29] with dynamics given by

\[\ddot{\theta}(t)=\frac{mgl\sin{(\theta(t))}-\mu\dot{\theta}(t)+\mathrm{sat}(u(t)) }{ml^{2}}.\]

As shown in [33], we can rewrite the dynamics as a four dimensional polynomial system. We let \(z_{1}=\theta\), \(z_{2}=\dot{\theta}\) and by making the substitution \(z_{3}=\sin(z_{1})\), \(z_{4}=\cos(z_{1})\) the system can be written as

\[\dot{z}_{1} =z_{2},\] \[\dot{z}_{2} =\frac{g}{l}z_{3}-\frac{\mu}{ml^{2}}z_{2}+\frac{1}{ml^{2}}u,\] \[\dot{z}_{3} =z_{2}z_{4},\] \[\dot{z}_{4} =-z_{2}z_{3},\]

where \(m=0.15\,\mathrm{kg},\,l=0.5\,\mathrm{m},\,\mu=0.5\,\mathrm{Nmsrad}^{-1},\,g=9.81 \,\mathrm{m}\mathrm{s}^{-2}\) and the controller input is saturated such that \(-1\leq u\leq 1\). The system also requires the equality constraint

\[z_{3}^{2}+z_{4}^{2}-1=0,\]

to be enforced. We do not define any region of the state space and instead consider global stability. We include the additional robustness constraints on the length of the pendulum to be \(\pm 0.1\) its original length and additive white noise \(w\) on the angular velocity such that \(||w||_{\infty}\leq 0.1\). By making the substitution \(\delta=1/l\) the full dynamical system can be written as

\[\dot{z}_{1} =z_{2},\] \[\dot{z}_{2} =g\delta z_{3}-\frac{\mu\delta^{2}}{m}z_{2}+\frac{\delta^{2}}{m}u +w,\] \[\dot{z}_{3} =z_{2}z_{4},\] \[\dot{z}_{4} =-z_{2}z_{3},\] \[0 =z_{3}^{2}+z_{4}^{2}-1,\] \[0 \leq 1^{2}-u^{2},\] \[0 \leq 0.1^{2}-w^{2},\] \[0 \leq(1/0.4-\delta)(\delta-1/0.6).\]

The Lyapunov function must be carefully constructed due to the \(z_{4}\) state being equal to one at the origin. We therefore define the Lyapunov function to be the sum of two Lyapunov functions, the first of which is defined as

\[V_{1}(z_{3},z_{4})=a_{1}z_{3}^{2}+a_{2}z_{4}^{2}+a_{3}z_{4}+a_{4}\]

and the second is quadratic in \(z_{1}\) and \(z_{2}\). To ensure that the Lyapunov function is zero at the origin we must enforce

\[a_{2}+a_{3}+a_{4}=0.\]

To ensure that the Lyapunov function is positive definite, we define

\[\rho(z)=\epsilon_{1}z_{1}^{2}+\epsilon_{2}z_{2}^{2}+\epsilon_{3}(1-z_{4}),\]

where \(\epsilon_{1}\geq 0.1\), \(\epsilon_{2}\geq 0.1\), \(\epsilon_{3}\geq 0.1\).

The size of the rational NN controller is a two layer network with four nodes in each layer. By setting the size of the polynomials in the network to be between zeroth and fourth order, we are able to recover a controller. The trajectories for this NFL are shown in Figure 12. We can see that the controller initially drives the system states towards a manifold and then moves it towards the equilibrium. Since the control system is discontinuous at the manifold, further analysis is required to show stability of the system.

## 6 Conclusion

In this paper, we analyse the use of rational NNs in previous application areas. We present novel rational activation functions to approximate the traditional sigmoid and tanh functions and show how they can be used in robustness
problems for NFLs. We argue that rational activation functions can be replaced with a general rational NN structure where each layer is convex in the NN's parameters. We then propose a method to recover a stabilising controller from a feasibility test and then extend this approach to rational NNs. This structure is refined to make it more compatible when used in conjunction with SOS programming. Through numerous numerical examples we show how this approach can be used to recover stabilising rational NN controllers for NFLs with non-linear plants with noise and parametric uncertainty.

## Acknowledgements

This work was supported by EPSRC grants EP/L015897/1 (to M. Newton) and EP/M002454/1 (to A. Papachristodoulou) and the Tony Corner Research Fund.
