Combining Learning and Model Based Control via Discrete-Time Chen-Fliess Series

TL;DR

This paper introduces a control method that combines learning and model-based approaches for nonlinear systems using discrete-time Chen-Fliess series, effectively integrating physical models with data-driven control.

Contribution

It develops a novel control framework leveraging Chen-Fliess series for nonlinear plants, incorporating algebraic structures for multivariable systems, demonstrated on a Lotka-Volterra example.

Findings

Successful control of a two-input, two-output Lotka-Volterra system

Integration of physical model knowledge with learning control

Extension of Chen-Fliess series to multivariable discrete-time systems

Abstract

A learning control system is presented suitable for control affine nonlinear plants based on discrete-time Chen-Fliess series and capable of incorporating knowledge of a given physical model. The underlying noncommutative algebraic and combinatorial structures needed to realize the multivariable case are also described. The method is demonstrated using a two-input, two-output Lotka-Volterra system.

Tables

$L$	$J$	$T$	$\Delta$	$ϵ$
100	3	6	0.06	0.05

$\Delta {\alpha }_{12}$	$\Delta {\alpha }_{21}$	$\delta y_1$	$\delta y_2$	${\Vert \hat{u}\Vert }_{\infty }$
exact model		8.66$\times {10}^{-9}$	1.25$\times {10}^{-8}$	2
-5	0	0.012	0.007	2
0	-5	0.020	0.016	2
5	0	0.004	0.006	2
0	5	0.018	0.015	2
-10	0	0.016	0.012	1.5
0	-10	0.056	0.041	1.5
10	0	0.010	0.009	1.5
0	10	0.037	0.025	1.5
-20	0	0.023	0.024	0.5
0	-20	0.144	0.113	0.5
20	0	0.012	0.016	0.5
0	20	0.071	0.047	0.5
-50	0	0.092	0.096	0.5
50	0	0.010	0.028	0.5
0	50	0.062	0.095	0.5
model free		0.191	0.897	1

$\Delta {\alpha }_{12}$	$\Delta {\alpha }_{21}$	$\Delta {\beta }_2$	$\delta y_1$	$\delta y_2$	${\Vert \hat{u}\Vert }_{\infty }$
exact model			1.547$\times {10}^{-5}$	1.165$\times {10}^{-4}$	1.4
-5	0	0	0.071	0.003	1.4
0	-5	0	0.024	0.118	1.4
0	0	-5	0.010	0.157	1.4
model free			0.602	0.330	2

$\Delta {\alpha }_{12}$	$\Delta {\alpha }_{21}$	$\Delta {\beta }_2$	$\delta y_1$	$\delta y_2$	${\Vert \hat{u}\Vert }_{\infty }$
exact model			0.007	1.93$\times {10}^{-7}$	1.4
-5	0	0	0.208	0.005	1.4
0	-5	0	0.162	0.016	1.4
0	0	-5	0.544	0.010	1.4
model free			1.970	1.680	2

$\Delta {\alpha }_{12}$	$\Delta {\alpha }_{21}$	$\Delta {\beta }_2$	$\delta y_1$	$\delta y_2$	${\Vert \hat{u}\Vert }_{\infty }$
exact model			8.73$\times {10}^{-5}$	1.599$\times {10}^{-4}$	2
-5	0	0	0.009	0.002	1.4
0	-5	0	0.094	0.118	1.4
0	0	-5	0.071	0.089	1.4
model free			0.167	0.055	1

Equations

$$l_{\infty }^{m+1}(N_0):=\{\hat{u}=(\hat{u}(N_0),\hat{u}(N_0+1),\dots ):{\Vert \hat{u}\Vert }_{\infty }<\infty \},$$

$${\hat{F}}_c[\hat{u}](N)=\sum _{\eta \in X^{\ast }}(c,\eta )S_{\eta }[\hat{u}](N)$$

$$S_{x_i\eta }[\hat{u}](N)=\sum _{k=N_0}^N{\hat{u}}_i(k)S_{\eta }[\hat{u}](k)$$

$${\hat{u}}_i(N)={\int }_{(N-1)\Delta }^{N\Delta }{u}_i(t)dt,i=0,1,\dots ,m,$$

$$\hat{y}(N)={\hat{F}}_c^J[\hat{u}](N):=\sum _{\eta \in X^{\le J}}(c,\eta )S_{\eta }[\hat{u}](N),$$

$$\dot{z}(t)$$

$$y_j(t)$$

$$\hat{y}(N)={\varphi }^T(N){\theta }_0,N\ge 1,$$

$$\varphi (N)$$

$${\theta }_0$$

$${\hat{y}}_p(N):={\varphi }^T(N)\hat{\theta }(N-1).$$

$$\hat{\theta }(N)$$

$$e(N)$$

$$g(N-1)$$

$$P(N-1)$$

$$\frac{P(N-2)\varphi (N){\varphi }^T(N)P(N-2)}{1+{\varphi }^T(N)P(N-2)\varphi (N)}$$

$${\mathbf{J}}_{\bar{N}}(\theta ):=\sum _{N=1}^{\bar{N}}[y(N\Delta )-{\varphi }^T(N)\theta {]}^2+\frac{1}{2}{\Vert \theta -\hat{\theta }(0)\Vert }^2$$

$$S[\hat{u}](N,N_0)=\sum _{\eta \in X^{\ast }}\eta S_{\eta }[\hat{u}](N,N_0),$$

$$S_{x_i\eta }[\hat{u}](N,N_0)=\sum _{k=N_0}^N{\hat{u}}_i(k)S_{\eta }[\hat{u}](k,N_0)$$

$$S_{x_i\eta }[\hat{u}](N_0,N_0)={\hat{u}}_{x_i}(N_0)S_{\eta }[\hat{u}](N_0,N_0)$$

$$S[\hat{u}](N_0,N_0)=\sum _{k=0}^{\infty }({\hat{u}}_1(N_0)x_1{)}^k=:(1-{\hat{u}}_1(N_0)x_1{)}^{-1}.$$

$$S[\hat{u}](N+1,N_0)=c_u(N+1)S[\hat{u}](N,N_0)$$

$$S[\hat{u}](N,N_0)=\overleftarrow{\prod _{i=N_0}^N}{c}_u(i),$$

$$S[\hat{u}](N,N_0)=(1-{\hat{u}}_1(N)x_1{)}^{-1}\cdots (1-{\hat{u}}_1(N_0)x_1{)}^{-1}.$$

$$S[\hat{u}](1,0)$$

$$=1+({\hat{u}}_1(1)+{\hat{u}}_1(0))x_1+({\hat{u}}_1^2(1)+$$

$${\hat{u}}_1(1){\hat{u}}_1(0)+{\hat{u}}_1^2(0))x_1^2+({\hat{u}}_1^3(1)+$$

$${\hat{u}}_1^2(1){\hat{u}}_1(0)+{\hat{u}}_1(1){\hat{u}}_1^2(0)+{\hat{u}}_1^3(0))x_1^3+\cdots$$

$$(\hat{v}{\#}_M\hat{u})(N)$$

$$\displaystyle=\left\{\begin{array}[]{rcl}\hat{u}(N)&:&N_{a}\leq N\leq M\\ \hat{v}((N-M)+N_{c})&:&M<N\leq M+(N_{d}-N_{c})\end{array}\right.$$

$$l_{\infty ,e}^{m+1}(0):=l_{\infty }^{m+1}(0)\cup \{\hat{0}\},$$

$$S[\hat{v}]((N-M)+N_c,N_c)S[\hat{u}](M,N_a)=S[\hat{v}{\#}_M\hat{u}](N,N_a).$$

$$S[\hat{v}](N-M)S[\hat{u}](M)=S[\hat{v}{\#}_M\hat{u}](N).$$

$${\mathcal{M}}_C(X)$$

$$0\le N\le N_f<\infty \}.$$

$${\hat{y}}_p(N+1)$$

Full text

††thanks: Corresponding author

1 Introduction

The central attraction of applying learning/adaptive data science to control is the ability to learn and generalize plant dynamics from partial input-output data in order to react properly to new situations. This can be done off-line via a training phase and/or online during closed-loop operation. The most widely used methods at present are based on artificial neural networks (ANNs), reinforcement learning control, and local adaptive control.

Most of the modern ANN approaches to learning control have their origins in the work of McCulloch and Pitts in the 1950’s and is based on the computational capabilities of networks of individual units called neurons (McCulloch & Pitts, 1943). The approach was further developed by Rosenblatt to produce what is now known as perceptron multilayer feedforward nets (Rosenblatt, 1962). At its core, an ANN in learning control realizes a static map parameterized by a structured set of real parameters that operates on an input signal (Hunt, et al., 1992). These parameters are adjusted based on some learning strategy such as backpropagation. The family of mappings in the controller is assumed to be sufficiently rich to represent a wide variety of potential control laws. But usually only simulation based justifications are possible. A long standing criticism is that there are few theoretical results strongly linking the properties of ANNs to the learning/adaptive control problem (Polycarpou & Ioannou, 1992). Adding dynamics to ANNs to form recurrent neural networks (RNNs) was a natural step in the development of learning methodologies in control. The hope here was that they would better approximate dynamic input-output behaviors (Baldi & Hornik, 1996; Jin, et al., 1995; Kambhampati, et al., 2000; Shaefer & Zimmermann, 2006; Xiao-Dong, et al., 2005). They also permit online learning as the network’s state evolves over time in response to an applied input much like the state of the plant. But again a major drawback concerning RNNs is that they lack theoretical support for the intended application of control. One attempt to address this issue is the use of input convex neural networks for optimal control (Chen, et al., 2019). They take advantage of recent progress in deep learning optimization, but have significant computational overhead. Hence, they are not suitable for every control application.

Reinforcement learning control is based on the classical Lyapunov/Bellman function. Simply stated, the aim is to find the proper control action with respect to an overall long term objective. A dynamic programming paradigm is used to minimize at each time instant a quantity that measures the overall goodness of a given state. This is classically known in optimal control as the Bellman value function or a Lyapunov energy function. The actor-critic approach relies on the principle of assigning a cost to a given state and/or a proposed action in a way that this cost becomes a good predictor of eventual long term outcomes (Lewis & Vrabie, 2009; Lewis, et al., 2012a, b; Vrabie & Lewis, 2009; Vrabie, et al., 2009). It originated from specializing the work of Barto, who first introduced the so called adaptive critic (Barto, 1990). Lewis et al. provided a more systematic version of the idea using a more suitably posed goal. This work is also related to that described in Mendel & Fu, (1970); Narendra & Thathachar (1989). Specifically, a parametric form of a Lyapunov function (the critic system) is given and one attempts to fit parameters for the Lyapunov function and the proposed feedback law simultaneously. This is done by adjusting the parameters after a training event (occurring on a different time scale) via a steepest descent step.

The local adaptive control methodology is based on operating the system about a set of predetermined operating points for which robust controllers are designed using a corresponding set of linear models (Åström & Wittenmark, 1994; Slotine & Li, 1991). The learning consists of building an association between the current state and the appropriate controller. This is conceptually a variation of a gain scheduled controller combined with a pattern recognition device to choose the most suitable gain. While a very effective approach in some applications, the overall design is mainly verified by simulation.

The main goal of this paper is to present a type of learning control system for control affine nonlinear systems based on a discretization of the Chen-Fliess functional series or Fliess operator (Fliess, 1981, 1983; Isidori, 1995). It is well known that any analytic control affine state space system in continuous-time has an input-output map with a Fliess operator representation. Therefore, the structure of the proposed learning system contains learning units that are known a priori to be capable of approximating the input-output behavior of the plant to an arbitrary desired accuracy (Gray, et al., 2017a). The learning system is also capable of incorporating a given physical model or can be used to provide purely data-driven control. The approach is distinct from model based adaptive control in that the plant is not made to track the output of a reference model, and there is no adaptation of the model. It is also distinct from the local adaptive control approach in that there is no need to linearize models or develop a gain scheduling strategy. The method is demonstrated using a two-input, two-output Lotka-Volterra system. Some of these results have appeared in preliminary form (often without proof) in Gray, et al. (2017b, 2019a, 2019b); Venkatesh, et al. (2019).

The paper is organized as follows. In the next section, a brief summary of the key concepts concerning discrete-time Fliess operators is given. A type of learning unit is then described in Section 3 based on discrete-time Fliess operators along with a purely inductive implementation. In the subsequent section, it is shown how to combine this type of learning with model based control. The main conclusions of the paper are given in the last section, as well as directions for future research.

2 Discrete-Time Fliess Operators

In this section, a brief review of discrete-time Fliess operators is presented. For additional details, see Duffaut Espinosa, et al. (2018); Gray, et al. (2017a).

An alphabet $X=\{x_{0},x_{1},$ $\ldots,x_{m}\}$ is any nonempty and finite set of noncommuting symbols referred to as letters. A word $\eta=x_{i_{1}}\cdots x_{i_{k}}$ is a finite sequence of letters from $X$. The number of letters in a word $\eta$, written as $\left|\eta\right|$, is called its length. The empty word, $\emptyset$, is taken to have length zero. The collection of all words having length $k$ is denoted by $X^{k}$. Define $X^{\ast}=\bigcup_{k\geq 0}X^{k}$ and $X^{\leq J}=\bigcup_{k=0}^{J}X^{k}$. The former is a monoid under the concatenation product. Any mapping $c:X^{\ast}\rightarrow{\mathbb{R}}^{\ell}$ is called a formal power series. Often $c$ is written as the formal sum $c=\sum_{\eta\in X^{\ast}}(c,\eta)\eta$, where the coefficient $(c,\eta)$ is the image of $\eta\in X^{\ast}$ under $c$. The set of all noncommutative formal power series over the alphabet $X$ is denoted by ${\mathbb{R}}^{\ell}\langle\langle X\rangle\rangle$. It forms an associative ${\mathbb{R}}$-algebra under the Cauchy product.

Inputs are assumed to be sequences of vectors from the normed linear space

[TABLE]

where $\hat{u}(N):=[\hat{u}_{0}(N),\hat{u}_{1}(N),\ldots,\hat{u}_{m}(N)]^{T}$, $N\geq N_{0}$ with $\hat{u}_{i}(N)\in{\mathbb{R}}$, $\left|\hat{u}(N)\right|:=\max_{i\in\{0,1,\ldots,m\}}$ $\left|\hat{u}_{i}(N)\right|$, and $\left\|\hat{u}\right\|_{\infty}:=\sup_{N\geq N_{0}}\left|\hat{u}(N)\right|$. The subspace of finite sequences over $[N_{0},N_{f}]$ is denoted by $l_{\infty}^{m+1}[N_{0},N_{f}]$.

Definition 1**.**

Given a generating series c\in\mbox{{\mathbb{R}}^{\ell}\langle\langle X\rangle\rangle}, the corresponding discrete-time Fliess operator is defined as

[TABLE]

for any $N\geq N_{0}$, where

[TABLE]

with $x_{i}\in X$, $\eta\in X^{\ast}$, and $\hat{u}\in l_{\infty}^{m+1}[N_{0}]$. By assumption, $S_{\emptyset}[\hat{u}](N):=1$.

Following Grune & Kloeden (2001), select some fixed $u\in L_{1}^{m}[0,T]$ with $T>0$ finite. Choose an integer $L\geq 1$, let $\Delta:=T/L$, and define the sequence of real numbers

[TABLE]

where $N\in\{1,2,\ldots,L\}$. Assume $u_{0}=1$ so that $\hat{u}_{0}(N)=\Delta$. A truncated version of $\hat{F}_{c}$ will be useful,

[TABLE]

since numerically only finite sums can be computed. The main assertion proved in (Gray, et al., 2017a, Theorems 6 and 7) is that the class of truncated, discrete-time Fliess operators acts as a set of universal approximators with computable error bounds for their continuous-time counterparts described in Fliess (1981, 1983); Isidori (1995). In which case, they can be used to approximate any input-output system corresponding to an analytic control affine state space realization

[TABLE]

with increasing accuracy as $L$ and $J$ increase. This fact is exploited in the next subsection to create a type of learning unit for data generated by such dynamical systems.

3 Learning Unit Based on Discrete-Time Fliess Operator

The main objective of this section is to introduce a learning unit based on a truncated discrete-time Fliess operator whose coefficients are identified via a standard least-squares algorithm. In general, there is one learning unit per output channel. So without loss of generality it is assumed that $\ell=1$. First the basic architecture of the learning unit is described, and then an inductive implementation of the underlying learning algorithm is developed.

3.1 Learning Unit Architecture

The first step is to write (1) as an inner product

[TABLE]

where

[TABLE]

with $l={\rm card}(X^{\leq J})=\sum_{k=0}^{J}(m+1)^{k}=((m+1)^{J+1}-1)/m$ and assuming some ordering $(\eta_{1},\eta_{2},\ldots)$ has been imposed on the words in $X^{\ast}$. If some estimate of $\theta_{0}$ is available at time $N-1$, say $\hat{\theta}(N-1)$, then (3) gives a corresponding estimate of $\hat{y}(N)$:

[TABLE]

The following least-squares algorithm is used to update the series coefficients:

[TABLE]

for any $N\geq 1$ with the initial estimate $\hat{\theta}(0)$ given, and $P(-1)$ is any positive definite matrix $P_{0}$ (Goodwin & Sin, 2009, p. 65). Covariance resetting is done periodically to enhance convergence. The corresponding learning unit is shown in Figure 1. Here input-output data $(u,y)$ from some unknown continuous-time plant (or the error system between the plant and an assumed model) is fed into the unit. The only assumption is that the data came from a system which has a Fliess operator representation, for example, any system modeled by (2).

In general, the learning unit has no a priori knowledge of the system, so $\hat{\theta}(0)$ is initialized to zero. Setting $P_{0}=I$, it is known that this algorithm minimizes the performance index

[TABLE]

with respect to the parameter $\theta$. It should be stated that since the model class consists of truncated versions of $\hat{F}_{c}$, there is no reason to expect the parameter vector $\hat{\theta}(N)$ to converge to $c$ in any fashion as $N$ increases. But this is not a problem since the only objective is to ensure that the underlying continuous-time input-output map $F_{c}$ is well approximated by $\hat{F}^{J}_{\hat{\theta}(N)}$. On the other hand, the approximation theory presented in Gray, et al. (2017a) guarantees that if the underlying system has such a Fliess operator representation then the true generating series $c$ is a feasible limit point for the sequence $\hat{\theta}(N)$, $N\geq 0$.

3.2 Inductive Implementation

To devise an inductive implementation of the learning unit, it is necessary to identify the algebraic structure underling the iterated sums in the definition of the discrete-time Fliess operator. The starting point for this is the following concept.

Definition 2**.**

Gray, et al. (2017a)** Given any $N\geq N_{0}$ and $\hat{u}\in l_{\infty}^{m+1}(N_{0})$, a discrete-time Chen series is defined as

[TABLE]

where

[TABLE]

with $x_{i}\in X$, $\eta\in X^{\ast}$, and $S_{\emptyset}[\hat{u}](N,N_{0}):=1$. If $N_{0}=0$ then $S[\hat{u}](N,0)$ is abbreviated as $S[\hat{u}](N)$.

Let $X$ be arbitrary and define $\hat{u}_{\eta}(N)=\hat{u}_{i_{k}}(N)\cdots\hat{u}_{i_{1}}(N)$ for any $\eta=x_{i_{k}}\cdots x_{i_{1}}\in X^{\ast}$ and $N\geq N_{0}$ with $\hat{u}_{\emptyset}(N):=1$. In addition, $c_{u}(N):=\sum_{\eta\in X^{\ast}}\hat{u}_{\eta}(N)\eta$. Then

[TABLE]

so that $S_{\eta}[\hat{u}](N_{0},N_{0})=\hat{u}_{\eta}(N_{0})$, and thus, $S[\hat{u}](N_{0},N_{0})=c_{u}(N_{0})$.

Example 3**.**

If $X=\{x_{1}\}$ and $\hat{u}_{x_{1}}(N_{0})=\hat{u}_{1}(N_{0})$, then

[TABLE]

A key observation is that a discrete-time Chen series $S[\hat{u}](N,N_{0})$ satisfies a difference equation as described next and proved in Appendix A

Theorem 4**.**

(Gray, et al., 2019a)* For any $\hat{u}\in l_{\infty}^{m+1}(N_{0})$ and $N\geq N_{0}$*

[TABLE]

with $S[\hat{u}](N_{0},N_{0})=c_{u}(N_{0})$ so that

[TABLE]

where $\overleftarrow{\prod}$ denotes a directed product from right to left.

Example 5**.**

Consider the case in Example 3 where $X=\{x_{1}\}$ and $\hat{u}_{x_{1}}(i)=\hat{u}_{1}(i)$ for all $i\geq N_{0}$. Then $c_{u}(i)=\sum_{k\geq 0}(\hat{u}_{1}(i)x_{1})^{k}=(1-\hat{u}_{1}(i)x_{1})^{-1}$ and

[TABLE]

For instance,

[TABLE]

In this case, $S[\hat{u}](N,N_{0})$ is always a rational series (Berstel & Reutenauer, 1988).

Consider next two input sequences $(\hat{u},\hat{v})\in l_{\infty}^{m+1}[N_{a},N_{b}]\times l_{\infty}^{m+1}[N_{c},N_{d}]$ with $N_{b}>N_{a}$ and $N_{d}>N_{c}$. The concatenation of $\hat{u}$ and $\hat{v}$ at $M\in[N_{a},N_{b}]$ is taken to be

[TABLE]

as shown in Figure 2. Define the set of sequences

[TABLE]

where $\hat{\mathbf{0}}$ denotes the empty sequence with duration zero so that formally $\hat{v}\#_{M}\hat{\mathbf{0}}=\hat{\mathbf{0}}\#_{M}\hat{v}:=\hat{v}$ for all $\hat{v}\in l_{\infty,e}^{m+1}(0)$. In which case, $l_{\infty,e}^{m+1}(0)$ is a monoid under this input concatenation operator. Define $S[\hat{\mathbf{0}}]=1$. The following is a straightforward generalization of Theorem 4.

Theorem 6**.**

(Gray, et al., 2019a)* (Discrete-time Chen’s identity) Given $(\hat{u},\hat{v})\in l_{\infty}^{m+1}[N_{a},N_{b}]\times l_{\infty}^{m+1}[N_{c},N_{d}]$, $M\in[N_{a},N_{b}]$, and $N\in[M,M+(N_{d}-N_{c})]$ it follows that*

[TABLE]

In particular, when $N_{a}=N_{c}=0$ then

[TABLE]

Define the set of discrete-time Chen series

[TABLE]

Theorem 7**.**

(Gray, et al., 2019a)* ${\mathcal{M}}_{C}(X)$ is a monoid under the Cauchy product. In addition, $S:l^{m+1}_{\infty,e}(0)\rightarrow{\mathcal{M}}_{C}(X)$ is a monoid homomorphism.*

{pf}

The results follow directly from (8).

Let ${\rm End}({\mathbb{R}}^{\infty})$ be the set of endomorphisms on the ${\mathbb{R}}$-vector space of real right-sided infinite sequences. This set can be viewed as the monoid of doubly infinite matrices with well defined matrix products and unit $I={\rm diag}(1,1,\ldots)$. A monoid $M$ is said to have an infinite dimensional real representation, $\Pi$, if the mapping $\Pi:M\rightarrow{\rm End}({\mathbb{R}}^{\infty})$ is a monoid homomorphism. The representation is faithful if $\Pi$ is injective.

Theorem 8**.**

(Gray, et al., 2019a)* The monoid ${\mathcal{M}}_{C}(X)$ has a faithful infinite dimensional real representation $\Pi$ given by $\Pi(S[\hat{u}](N))=\overleftarrow{\prod_{i=0}^{N}}{\mathcal{S}}(i)$, where ${\mathcal{S}}(i)$ is any matrix representation of the ${\mathbb{R}}$-linear map on ${\mathbb{R}}\langle\langle X\rangle\rangle$ given by the left concatenation map ${\mathcal{C}}:d\mapsto c_{u}(i)d$.*

{pf}

The representation claim follows from (7). To see that $\Pi$ is injective, assume a fixed ordering of the words in $X^{\ast}$, say $\{\eta_{1},\eta_{2},\ldots\}$. Then define the matrix $[{\mathcal{S}}(i)]_{jk}=(c_{u}(i)\eta_{k},\eta_{j})=\hat{u}_{\xi}(i)$, where $\xi\eta_{k}=\eta_{j}$. Thus, ${\mathcal{S}}(i)$ is a lower triangular matrix with ones along the diagonal since $u_{\emptyset}(i)=1$, $i\geq 0$. The first column is comprised of the coefficients of $c_{u}(i)$ in the order given to $X^{\ast}$. Hence, the map $\Pi$ on the monoid ${\mathcal{M}}_{C}(X)$ is injective since $c_{u}(i)$ can be uniquely identified from ${\mathcal{S}}(i)=\Pi(S[\hat{u}](i,i))$.

Note that the above theorem implies that (4) can be written in the form

[TABLE]

for $N\geq N_{0}$, where ${\mathcal{S}}(N+1)$ and $S[\hat{u}](N)$ have been suitable truncated, and $e_{1}:=[1\,0\,0\cdots 0]^{T}\in{\mathbb{R}}^{l}$. Equation (9) can also be written in the form $\hat{y}_{p}(N+1)=Q(\hat{u}(N+1))$, where $Q$ is a polynomial in the components of $\hat{u}(N+1)$ with maximum degree $l-1$.

Example 9**.**

Suppose $X=\{x_{1}\}$ as in Example 5. Assuming the ordering on $X^{\ast}$ to be $\{\emptyset,x_{1},x_{1}^{2},\ldots\}$. Then for all $i\geq 0$

[TABLE]

and $c_{u}(i)=\sum_{k\geq 0}\hat{u}^{k}_{1}(i)x_{1}^{k}$. In addition,

[TABLE]

As expected, the first column coincides with the coefficients of $S[\hat{u}](1)$ in Example 5. Setting $J=3$ so that $l={\rm card}(X^{\leq J})=4$ gives the truncated versions

[TABLE]

where $S_{x_{1}^{k}}(N):=(S[\hat{u}](N),x_{1}^{k})$. Therefore, the output

[TABLE]

where the coefficients $q_{i}(N)$ are functions of $(c,x_{1}^{k})$ and $S_{x_{1}^{k}}(N)$, $k=0,1,2,3$.

As a specific example, consider a plant modeled by the Fliess operator

[TABLE]

where the generating series is $c=\sum_{k\geq 0}x_{1}^{k}$, and

[TABLE]

with $E_{\emptyset}[u]:=1$. The system has the state space realization

[TABLE]

since for all $t\geq 0$

[TABLE]

The output $y$ shown in Figure 3 is computed from a numerical simulation of the state space model (11) when the input $u(t)=2e^{-t/3}\sin(2\pi t)$ is applied. The output of the learning unit $\hat{y}_{p}(N)$, $N\geq 0$ as implemented using (5), (9), and (10) is also shown in the figure. As the learning unit processes more data, its estimate of the output $y$ improves asymptotically.

The more challenging problem is systematically building a real representation of ${\mathcal{M}}_{C}(X)$ when $X$ has more than one letter, as in the multivariable case or when the drift letter $x_{0}$ is present. A partial ordering $\preceq$ is first defined on all words in $X^{\ast}$. For all \zeta,\eta\in\mbox{X^{\ast}}, let $\zeta\preceq\eta$ if and only if there exists a \gamma\in\mbox{X^{\ast}} such that $\gamma^{-1}(\eta)=\zeta$, where $\gamma^{-1}$ denotes the left-shift operator. The following theorem is proved in Appendix B.

Theorem 10**.**

(Venkatesh, et al., 2019)* The pair \left(\mbox{X^{\ast}},\preceq\right) is a partially ordered set.*

The partial order (\mbox{X^{\ast}},\preceq) can be graphically represented by a Hasse diagram. Starting with $\emptyset$ at the root, the Hasse diagram of (\mbox{X^{\ast}},\preceq) when $X=\{x_{0},x_{1},\ldots,x_{m}\}$ forms a ($m$+$1$)-ary infinitely branching tree. Define an injective map $R:X\longrightarrow{\mathcal{C}}$, where ${\mathcal{C}}$ is a set of colors. Color the edge between the nodes $\eta$ and $x_{i}\eta$ with the color $R(x_{i})$ in the tree. As an illustration, the tree for the case when $m=2$ is shown in Figure 4, where $R(x_{0})=\text{black}$, $R(x_{1})=\text{red}$, and $R(x_{2})=\text{blue}$.

In (9), the underlying discrete-time Fliess operator $\hat{F}_{c}\left[\hat{u}\right]$ has been truncated up to words of length $J$. Therefore, the tree is pruned at the $J$-th level. Next, a depth-first search (DFS) algorithm is employed to traverse the graph and generate words. The corresponding vector of words, $\chi^{J}(X)$, is called the order vector of degree $J$ and is given by $\chi^{0}(X)=[\emptyset]$ and

[TABLE]

for $J\geq 0$.

Example 11**.**

The tree for words $\eta\in X^{\leq 2}$ when $X=\{x_{0},x_{1},x_{2}\}$ is given by

The DFS algorithm gives the order vector

[TABLE]

Let ${\mathcal{S}}^{J}(N+1)$ denote the matrix ${\mathcal{S}}(N+1)$ truncated for words up to length $J$, i.e., ${\mathcal{S}}^{J}(N+1)\in{\mathbb{R}}^{l\times l}$ with $l={\rm card}(X^{\leq J})$. An inductive algorithm to build such matrices is developed next.

Definition 12**.**

Define $C_{i}$ as the colored tree of the Hasse diagram of (\mbox{X^{\ast}},\preceq) up to the $i$-th level, that is, the $(m+1)$-ary tree with $\emptyset$ as the root and $\eta\in X^{i}$ as leaves of the Hasse diagram. Let $C:=\{C_{i}:i\in\mathbb{N}_{0}\}$ be the set of colored trees given by (\mbox{X^{\ast}},\preceq) of all levels.

It is useful to define a product $\dagger$ on $C$ as follows: $C_{i}\dagger C_{j}:=${tree with each leaf node $\beta\in X^{i}$ replaced by the tree $C_{j}$, where all the nodes of $C_{j}$ are right concatenated with $\beta$}. The following theorem is proved in Appendix C.

Theorem 13**.**

$(C,\dagger)$* is a commutative monoid isomorphic to the additive monoid $({\mathbb{N}}_{0},+)$. Specifically, $C_{i}\dagger C_{j}$ $=C_{i+j}$ for all $C_{i},C_{j}\in C$.*

Example 14**.**

Let $X=\{x_{0},x_{1}\}$ and define the color map $R$ as: $R(x_{0})=\text{red}$, $R(x_{1})=\text{blue}$. Observe that

[TABLE]

[TABLE]

[TABLE]

The above tree can be expanded as

{forest}

[$\emptyset$ [$x_{0}$, edge = red [$x_{0}^{2}$, edge = red[$x_{0}^{3}$, edge = red] [$x_{1}x_{0}^{2}$, edge = blue]] [$x_{1}x_{0}$, edge = blue [$x_{0}x_{1}x_{0}$, edge = red] [$x_{1}^{2}x_{0}$, edge = blue]] ] [$x_{1}$, edge = blue [$x_{0}x_{1}$, edge = red[$x_{0}^{2}x_{1}$, edge = red] [$x_{1}x_{0}x_{1}$, edge = blue]] [$x_{1}^{2}$, edge = blue [$x_{0}x_{1}^{2}$, edge = red] [$x_{1}^{3}$,edge = blue]] ] ]

This final tree is identified as $C_{3}$ so that $C_{1}\dagger C_{2}=C_{3}$.

Now assume that each color $R(x_{i})$, $x_{i}\in X$, is given the weight $\hat{u}_{i}(N+1)$ at the discrete time instant $N+1$. Then it follows for any $\eta_{j},\eta_{k}\in X^{\ast}$ that

[TABLE]

By Theorem 13, in the case where $X=\{x_{0},x_{1}\}$ with color map $R$ defined as in Example 14, $C_{J+1}=C_{1}\dagger C_{J}$. That is,

[TABLE]

Hence, from the structure of the order vector in (12) and the above tree recursion, one can deduce for any $m\geq 1$ that the block structure of the matrix ${\mathcal{S}}^{J+1}(N+1)$ can be written inductively in terms of ${\mathcal{S}}^{J}(N+1)$ as:

[TABLE]

where ‘$\otimes$’ denotes the Kronecker matrix product, and the block diagonal matrix is comprised on $m+1$ blocks.

Example 15**.**

Let $X=\{x_{0},x_{1}\}$. For $J=2$, the words are indexed by $\chi^{2}(X)=[\emptyset\;x_{0}\;x_{0}^{2}\;x_{1}x_{0}\;x_{1}\;x_{0}x_{1}\;x_{1}^{2}]$. ${\mathcal{S}}^{2}(N+1)$ can be computed directly from $C_{2}$ to be

[TABLE]

(For brevity, the argument $(N+1)$ is suppressed in the elements of the matrix.) But the same matrix can also be computed inductively from (17). For the base case, ${\mathcal{S}}^{0}(N+1)=1$ so that

[TABLE]

Applying (17) once more gives (18).

4 Combining Learning and Model Based Control

Suppose $y_{d}$ is a desired output known to be in the range of a given plant with an underlying but unknown Fliess operator representation $F_{c}$. It is most likely in applications that $y_{d}$ was designed using an assumed model (2) with perhaps the aid of some expert knowledge. Already this implies that the model is not too poor an approximation of the plant, otherwise $y_{d}$ may not be in its range of the true plant. When both the plant and model are given the same input, a modeling error $e_{i}=y_{i}-\hat{y}_{i}$ is generated for the $i$-th channel of the plant’s output as shown in Figure 5. This signal and the applied input are then fed to a learning unit of the type presented in the previous section in order to learn the input-output behavior of each error map $u\mapsto e_{i}$, $i=1,\ldots,m$, which in this case must also have a Fliess operator representation. At any time instant the output of the plant is approximated by $\hat{y}+\hat{e}_{p}=\hat{y}+Q(\hat{u})$, where $Q$ was defined using (9). A suitable input $u$ for tracking $y_{d}$ can be approximated by a piecewise constant function taking values for $N\in\{1,2,\ldots,L\}$ equivalent to

[TABLE]

for some fixed bound $\bar{u}>0$ and where

[TABLE]

with $W\in{\mathbb{R}}^{\ell\times\ell}$ being a fixed symmetric positive semi-definite weighting matrix. The MatLab command fmincon can be used to compute these local minima over the interval $[-\bar{u},\bar{u}]$. In summary then equations (5), (9), and (19) provide a fully inductive implementation of a one step ahead predictive controller with learning. If the model is omitted from this set up, the resulting controller is still viable and can be viewed as a type of data-driven/model free closed-loop system as first proposed for SISO systems in Gray, et al. (2017b).

As an example, consider the classical Lotka-Volterra model

[TABLE]

used to describe the population dynamics of $n$ species in competition (Chauvet, et al., 2002; May & Leonard, 1975; Smale, 1976). Here $z_{i}$ is the biomass of the $i$-th species, $\beta_{i}$ represents the growth rate of the $i$-th species, and the parameter $\alpha_{ij}$ describes the influence of the $j$-th species on the $i$-th species. More recently in Jafarian, et al. (2018) it was shown that a power network, where each node voltage $z_{i}$ is regulated by a quadratic droop controller, has dynamics described by a Lotka-Volterra model. In general, this model can exhibit a wide range of behaviors including the presence of multiple stable equilibria, stable limit cycles, and even chaotic behavior.

Consider the case where a subset of system parameters $\beta_{i_{j}}$, ${j=1,\ldots,m}$ in (20) can be actuated and thus viewed as inputs $u_{i}$, $i=1,\ldots,m$. Assume some set of output functions is given

[TABLE]

Since the inputs enter the dynamics linearly, it is clear that (20)-(21) constitute a control affine analytic state space system. In which case, the input-output map $u\mapsto y$ has an underlying Fliess operator representation $F_{c}$ with generating series c\in\mbox{{\mathbb{R}}^{\ell}\langle\langle X\rangle\rangle} computable directly from (20)-(21) and a given initial condition $z_{0}$ (Fliess, 1981, 1983; Isidori, 1995). Of particular interest here is the special case of a predator-prey system, which is a two dimensional Lotka-Volterra system

[TABLE]

where $y_{1}=z_{1}$ and $y_{2}=z_{2}$ are taken to be the population of prey and predator species respectively, and (20) has been re-parameterized so that $\beta_{i},\alpha_{ij}>0$. This positive system has precisely two equilibria when all the parameters are fixed, namely, a saddle point equilibrium at the origin and a center at $z_{e}=(\beta_{2}/\alpha_{21},\beta_{1}/\alpha_{12})$ corresponding to periodic solutions.

Taking the system inputs in (22) and (23) to be $u_{1}=\beta_{1}$ and $u_{2}=\beta_{2}$, the orbit transfer problem as shown in Figure 6 is to determine an input to drive the system from some initial orbit to within an $\epsilon$ neighborhood of a final orbit using a given orbit transfer trajectory. The proposed controller was tested in simulation assuming all the plant’s parameters are set to unity. The discretization parameters were selected as in Table 1, and the sampled input was bounded by $\|\hat{u}\|_{\infty}$ as given in Table 2 (the positivity constraint on the input was not enforced). The tracking performance for various choices of model parameter errors is shown in Table 2. For each plant parameter $\lambda$, $\Delta\lambda:=(\lambda_{\text{model}}-\lambda_{\text{plant}})\times 100\%$. In addition, $\delta y_{i}$ for $i=1,2$ is the RMS error per sample normalized by the sample value of desired trajectory for the given output channel. First the control system was tested assuming the exact plant model is available. In which case, the modeling error is zero and the learning units are inactive. The closed-loop performance is therefore determined solely by the predictive controller, which is quite accurate as shown in Table 2. Next a variety of parametric errors were introduced in the model. For all such cases, the weighting matrix $W$ was set to the identity matrix. As an example, the simulation results for the case of $+20\%$ error in $\alpha_{21}$ are shown in Figures 7-9. The performance of the system for the case where there is $-20\%$ error in $\alpha_{12}$ was very similar. The case where $\Delta\alpha_{21}=-20\%$ is an extreme scenario as decreasing $\alpha_{21}$ any further resulted in the plant’s response being oscillatory. The simulation results pertaining to this case are shown in Figures 10-12. Finally, the model free case was also simulated as shown in Figures 13-15. Here it was necessary to select a nontrivial weighting matrix, in this case $W=[1\;0.25;\;0.25\;1]$, in order for the optimizer to compensate for the cross coupling between the input-output channels, something that was done automatically when a model was present. Note that tracking was achieved, but the performance was about an order of magnitude worse than most cases employing a model. While in practice input bounds are dictated by the physical application, it was observed here that the larger the modeling error, the more conservative the input bounds needed to be in order to avoid instabilities. On the other hand, if the bounds were too conservative, then there was not enough actuation energy available to follow the desired trajectory. For the sake of comparison to earlier work reported in Gray, et al. (2019b), where only a single learning unit was employed and thus only SISO and SIMO control (with $W=[1\;0.25;\;0.25\;2]$) was possible, the simulation results are summarized in Tables 3-5. As a general rule, the lack of a second control input rapidly reduced performance when parameter errors exceeded five percent. But of course controller complexity was also significantly reduced in this case.

5 Conclusions and Future Work

A learning system for nonlinear control was presented based on discrete-time Chen-Fliess series and capable of incorporating a given physical model. A fully inductive implementation in the multivariable case required one to exploit the underlying noncommutative algebraic and combinatorial structures in order to identify a convenient basis to represent the learning dynamics. The method was demonstrated using a two-input, two-output Lotka-Volterra system.

Future work will include the introduction of measurement noise in the system, exercising the method on more complex engineering plants, and identifying conditions under which closed-loop stability can be guaranteed.

{ack}

This research was supported by the National Science Foundation under grants CMMI-1839378 and CMMI-1839387.

Appendix A Proof of Theorem 4

The first identity is addressed by proving that

[TABLE]

via induction on the length of $\eta$. When $\eta=\emptyset$ then trivially $S_{\emptyset}[\hat{u}](N+1,N_{0})=1=\hat{u}_{\emptyset}(N+1)S_{\emptyset}[\hat{u}](N,N_{0})$. If $\eta=x_{i}\in X$ then from (6)

[TABLE]

Finally, assume the identity holds for all words up to some fixed length $n\geq 0$. Then for any $\eta\in X^{n}$ and $x_{i}\in X$ it follows that

[TABLE]

which proves the claim for all $\eta\in X^{\ast}$. The second identity in the theorem follows directly from the first.

Appendix B Proof of Theorem 10

Let \eta,\zeta,\gamma,\alpha,\beta\in\mbox{X^{\ast}}. Reflexivity is trivial since $\emptyset^{-1}(\eta)=\eta$ if and only if $\eta\preceq\eta$. To prove transitivity, first observe that

[TABLE]

so that

[TABLE]

Therefore,

[TABLE]

To prove anti-symmetry, note that

[TABLE]

and therefore,

[TABLE]

Hence, (\mbox{X^{\ast}},\preceq) is a partially ordered set.

Appendix C Proof of Theorem 13

A preliminary lemma is needed first. For any fixed $\eta\in X^{\ast}$ define the right concatenation map as

[TABLE]

Lemma 16**.**

Every right concatenation map is an order embedding map on $(X^{\ast},\preceq)$. That is, $\zeta\preceq\gamma$ if and only if $\zeta\eta\preceq\gamma\eta$ for all $\zeta,\gamma,\eta\in X^{\ast}$.

{pf}

From the definition of $\preceq$ it follows that

[TABLE]

Applying ${\mathcal{R}}_{\eta}$ to both sides of the equality above gives

[TABLE]

PROOF (Theorem 13). First the identity $C_{i}\dagger C_{j}=C_{i+j}$ is proved. From the definition of the tree $C_{j}$ and Lemma 16 it is clear that $C_{j}\eta$ has a Hasse diagram with $\eta$ as the root and $X^{j}\eta$ as the set of leaf nodes. By the definition of the dagger product, every leaf node $\eta$ of $C_{i}$ is replaced by $\eta C_{j}$ since

[TABLE]

Therefore, $C_{i}\dagger C_{j}$ has a Hasse diagram with $\emptyset$ as the root and $X^{i+j}$ as the set of leaf nodes, that is, $C_{i}\dagger C_{j}=C_{i+j}$. It is now easily checked using this identity that $(C,\dagger)$ is associative, commutative, and has $C_{0}$ as the unit. Hence, $(C,\dagger)$ forms a commutative monoid. The monoid isomorphism between ${\mathbb{N}}_{0}$ and $C$ is given by the bijection $i\mapsto C_{i}$ for all $i\in{\mathbb{N}}_{0}$.