# Output Feedback Reinforcement Learning with Parameter Optimisation for
Temperature Control in a Material Extrusion Additive Manufacturing system

Eleni Zavrakli*1,2,3 , Andrew Parnell 1,2,3 and Subhrakanti Dey4
*Corresponding author<EMAIL_ADDRESS>of Mathematics and
Statistics, Maynooth University, Co. Kildare, Ireland2I-Form Advanced
Manufacturing Research Centre, Ireland3Hamilton Institute, Maynooth
University, Co. Kildare, Irealand4Division of Signals and Systems, Department
of Electrical Engineering, Uppsala University, Sweden

###### Abstract

With the rapid development of Additive Manufacturing (AM) comes an urgent need
for advanced monitoring and control of the process. Many aspects of the AM
process play a significant role in the efficiency, accuracy and repeatability
of the process, with temperature regulation being one of the most important
ones. In this work, we solve the problem of optimal tracking control for a
state space temperature model of a Big Area Additive Manufacturing (BAAM)
system. In particular, we address the problem of designing a Linear Quadratic
Tracking (LQT) controller when access to the exact system state is not
possible, except in the form of measurements. We initially solve the problem
with a model-based approach based on reinforcement learning concepts, with
state estimation through an observer. We then design a model-free
reinforcement-learning based controller with an internal state estimation step
and demonstrate its performance through a simulator of the systems’ behaviour.
Our results showcase the possibility of achieving comparable results while
learning optimal policies directly from process data, without the need for an
accurate, intricate model of the process. We consider this outcome to be a
significant stride towards autonomous intelligent manufacturing.

Keywords: Output Feedback; Reinforcement Learning; Intelligent Manufacturing;
Optimal tracking;

Acknowledgements: This publication has emanated from research supported by a
research grant from Science Foundation Ireland (SFI) under Grant Number
16/RC/3872. For the purpose of Open Access, the author has applied a CC BY
public copyright licence to any Author Accepted Manuscript version arising
from this submission.

Competing interests: The authors have no competing interests to declare that
are relevant to the content of this article.

Code availability: The code used for our simulations was written in R and can
be found on: https://github.com/elenizavrakli/Output-feedback-RL-temperature-
control

##  1 Introduction

Additive Manufacturing (AM) is becoming one of the most popular manufacturing
techniques, with its applications ranging from aerospace to bio-medicine to
household items [1, 2, 3]. Its appeal lies in the ability to build very
complicated structures without the need for a mold and without the waste of a
lot of materials. There is an array of different AM techniques, using
different materials, and targeting different problems. For example, in the
industrial scale, metal-based AM is being adopted at an increasing scale due
to its cost-effective and energy-saving nature [4, 5, 6]. Material extrusion
(MEX) on the other hand, a technique using different polymers such as PLA,
allows consumers to design and manufacture their own prototypes at home,
making it more appealing to hobbyists and small business owners, neither of
which would normally have access to traditional manufacturing techniques and
equipment. As AM becomes increasingly prevalent, the need naturally arises for
the process to be optimised in terms of safety, repeatability and cost
effectiveness. This objective can be addressed by efficient closed-loop
control algorithms, which are currently lacking in the AM space [7]. With many
different aspects of AM processes needing to be monitored and controlled,
developing effective feedback controllers to target them is a very active
field of study [8].

Control Theory [9, 10] offers a wide variety of methods for designing
efficient controllers, normally when a mathematical model describing the
process is available. In the AM context, proportional integral derivative
(PID) controllers have been widely used to control various aspects of the
process [11, 12, 13, 14, 15]. In AM, similarly to many other physical systems,
obtaining accurate models of the process is very challenging, which is an
active research field on its own [16]. With that in mind, the development of
controllers tuned through information learned from process data is being
increasingly investigated. Reinforcement Learning (RL) [17, 18] commonly
refers to a category of data-centric Machine Learning techniques that learn
optimal strategies within intricate and unpredictable settings by maximizing
rewards (or minimizing costs). In the field of Control Theory, decision-making
problems expressed as optimising strategies are formulated using Dynamic
Programming [19, 20, 21], which shares several commonalities with RL but
employs distinct terminology. While certain RL methods have been designed
deterministically, assuming a priori knowledge of the system dynamics, the
majority of RL applications are rooted in data-driven, model-free approaches.
In recent years, as RL has gained popularity, researchers have employed it to
address various issues within AM [22, 23, 24, 25, 26, 27]. However, none of
these papers use RL as a tool to design optimal feedback controllers. For more
references we refer to our previous work [28], where we investigated the
problem of temperature control in an AM system, when the exact system state is
directly available. The problem was solved using both standard Control Theory
approaches and RL techniques, both in a model-based and a data-driven manner.

In this paper, we address the problem of temperature control for a MEX AM
system without access to the exact system state, but with available sensor
measurements only. Specifically, we control the temperatures within the
extruder of a Big Area Additive Manufacturing (BAAM) system. We design a
controller that solves the Linear Quadratic Tracking (LQT) problem, derived
directly from the process model. However, we ultimately aim to have a
controller that can learn optimal policies directly from process data. We
utilise Reinforcement Learning techniques to design a data-driven LQT
controller [29]. In the model-based case, the state can be estimated with the
help of an observer [30]. In the case of the data-driven controller, we
incorporate a data-driven state-estimation step in the controller, using a
record of past inputs and outputs from the system only, inspired by the work
done in [31]. Furthermore, we optimise the parameters used to define the LQT
problem using Bayesian Optimisation [32, 33]. As there is no systematic way to
choose these parameters, Bayesian Optimisation has been explored in the
literature due to its ability to handle optimising complex objectives when
trial-and-error approaches do not suffice [34, 35, 36]. Finally, while setting
up our simulations, we provide detailed guidelines on how to generate
appropriate data to train the data-driven controller. These data generation
tools can be useful in many settings, in particular when one needs to train a
learning algorithm for dynamical systems.

The structure of our paper is as follows: In Section 2, we introduce the BAAM
system we wish to control and define the optimisation problem to be solved.
Section 3 deals with the design of the model-based state observer and
controller, while in Section 4 we present the design of the RL-based data-
driven controller which includes state estimation through input-output data.
Section 5 introduces the main concepts of Bayesian Optimisation and Section 6
presents the results of our simulations, using both controllers, before and
after Bayesian Optimisation. Finally, in Section 7 we provide a thorough
discussion of our findings.

##  2 Linear Quadratic Output Feedback Control for a state-space system

Faulty prototypes often occur in AM processes due to issues related to
temperature regulation in various subprocesses in AM. More specifically, the
cause of these faults may be the insufficient or excessive melting of the
material, the unsuccessful binding between layers or the non-uniform heating
of the build-plate which can stop the print from adhering properly.

In this work, we study the heating system of a Material Extrusion (MEX)
system. In particular, we are focusing on the temperatures within the extruder
of a Big Area Additive Manufacturing (BAAM) system. The extruder consists of 5
main parts, the hopper, the barrel, the hose, the nozzle and the AC motor. The
hopper is where the material is inserted and gets fed to the barrel containing
the screw which gets rotated by the motor. The screw pushes the material
towards the hose and finally to the nozzle. There are four heaters placed
along the barrel, one heater in the hose and one in the nozzle. The reason
behind the use of multiple heaters in BAAM systems is to address the large
volume of material that gets processed by the system. The areas where the heat
is applied can be viewed as thermal cells, each influencing their neighbouring
cells and all getting influenced by the motor. We adopt a discrete time,
linear state-space representation of the system inspired by the system model
obtained in [37] through system identification methods using real process
data. We make the further assumption that the temperatures within the extruder
are not directly available, instead a set of measurements can be obtained
through a number of sensor measurements of the process.

$\displaystyle x(t+1)$ $\displaystyle=$ $\displaystyle Ax(t)+Bu(t)$
$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle Cx(t)$ (1)

where $t\geq t_{0}$ is the discrete time step, $x(t)\in\mathbb{R}^{n}$ is the
system state at time $t$, $u(t)\in\mathbb{R}^{m}$ is the system input at $t$
and $y(t)\in\mathbb{R}^{p}$ is the system output obtained as measurements.
Matrices $A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{n\times m}$ and
$C\in\mathbb{R}^{p\times n}$ are the state, input and measurement matrices
respectively. For our system, there are $m=7$ inputs, the first 6 representing
the load applied to each thermal cell by the corresponding heater and the last
relaying the input applied by the motor. The system state is a vector of
length $n=6$ representing the temperature within each thermal cell. We assume
that the sensor system obtains $p=5$ measurements of the internal states of
the system. A visual representation of the system setup can be seen in Figure
1.

Figure 1: Heating system for a BAAM extruder. The inputs include the heat
applied in each thermal cell and the input provided by the AC motor rotating
the screw. The system states represent the temperatures within the extruder,
which are measured by the sensors.

The temperature control problem is defined as the infinite horizon optimal
tracking problem, where the goal is to obtain the optimal input (or control)
$u^{*}$ such that the system output accurately tracks the reference signal
$r(t)\in\mathbb{R}^{p}$, generated according to matrix
$F\in\mathbb{R}^{p\times p}$

$r(t+1)=Fr(t).$ (2)

In most applications $r(t)$ is chosen to be constant, namely $F=I_{p}$ or
$r(t)=r\in\mathbb{R}$ for all $t$. However, the general, time-varying
reference signal (2) can also be adopted in our formulation. The control
objective can be formally expressed as the minimisation of the performance
index defined as the total sum of discounted costs over an infinite horizon,
namely:

$\displaystyle\min_{\mathbf{u}}J(x(0),\mathbf{u})$ $\displaystyle=$
$\displaystyle\min_{\mathbf{u}}\sum_{t=0}^{\infty}\gamma^{t}\left\\{[y(t)-r(t)]^{T}Q[y(t)-r(t)]+u^{T}(t)Ru(t)\right\\}$
$\displaystyle=$
$\displaystyle\min_{\mathbf{u}}\sum_{t=0}^{\infty}\gamma^{t}\left\\{[Cx(t)-r(t)]^{T}Q[Cx(t)-r(t)]+u^{T}(t)Ru(t)\right\\}$

where $0<\gamma\leq 1$ is the discount factor, putting greater emphasis on
short term costs. $Q\in\mathbb{R}^{p\times p}$ and $R\in\mathbb{R}^{m\times
m}$ are positive definite symmetric matrices applying weight to the tracking
error and the input cost respectively. $x(0)$ gives rise to all future states
$x(t)$ and we define $\mathbf{u}$ as the infinite sequence containing all
inputs $\\{u(0),u(1),u(2),\dots\\}$, also known as a policy. This means that
the minimisation problem is equivalent to the search for an optimal policy
$\mathbf{u}^{*}$. In the Reinforcement Learning literature [19, 20, 18], the
minimisation problem is addressed through the principle of optimality:

######  Definition 1 (Principle of optimality)

[19] An optimal policy has the property that whatever the initial state and
initial input are, the remaining inputs must constitute an optimal policy with
regard to the state resulting from the first decision.

To utilise the principle of optimality, we define the Value Function as the
cost-to-go at time $t$, namely the sum of discounted future costs from time
$t$ onwards.

$\displaystyle V(x(t),\mathbf{u}_{t})$ $\displaystyle=$
$\displaystyle\sum_{k=t}^{\infty}\gamma^{k-t}\left\\{[y(k)-r(k)]^{T}Q[y(k)-r(k)]+u^{T}(k)Ru(k)\right\\}$
(4) $\displaystyle=$
$\displaystyle\sum_{k=t}^{\infty}\gamma^{k-t}\left\\{[Cx(k)-r(k)]^{T}Q[Cx(k)-r(k)]+u^{T}(k)Ru(k)\right\\}$

where $\mathbf{u}_{t}=\\{u(t),u(t+1),u(t+2),\dots\\}$. The Value function
allows for the principle of optimality to be expressed mathematically as the
Bellman optimality equation [20, 18]. Let
$c(x(t),u(t))=[Cx(t)-r(t)]^{T}Q[Cx(t)-r(t)]+u^{T}(t)Ru(t)$ be the cost at time
$t$. Then for the optimal policy $\mathbf{u}^{*}$, the principle of optimality
yields

$V(x(t),\mathbf{u}^{*}_{t})=c(x(t),u^{*}(t))+\gamma
V(x(t+1),\mathbf{u}^{*}_{t+1})$ (5)

for all $t$. In the following sections we explore two different approaches to
solving the optimal tracking problem. The first is based on the design of a
state observer while the second utilises a record of past inputs and outputs
to estimate the internal state. While the observer design requires knowledge
of the system dynamics, the second approach can be adapted to be completely
model free and data-driven.

##  3 State Observer based Output Feedback

An important issue in the design of efficient controllers is the lack of
access to the exact system state $x(t)$. Assuming knowledge of the system
dynamics (2), a state observer can be implemented to provide an informed
estimate of the state, that can then be used for control design. This observer
is also known as a closed loop observer or a Luenberger observer [30, 10].

In practice, the Luenberger observer is another state space system, made up of
the same state, input and output matrices as the original system, where the
state equation has an additional term associated with the error between the
system output and the observer output.

$\displaystyle\hat{x}(t+1)$ $\displaystyle=$ $\displaystyle
A\hat{x}(t)+Bu(t)+L\left[y(t)-\hat{y}(t)\right]$ $\displaystyle\hat{y}(t)$
$\displaystyle=$ $\displaystyle C\hat{x}(t)$ (6)

The matrix $L\in\mathbb{R}^{n\times p}$ is known as the observer gain. The
goal is to have an observer that emulates the behaviour of the original system
as closely as possible. In control theory terms, we need to find a gain matrix
$L$ such that the matrix $A-LC$ is stable [10]. This problem can be approached
as a pole placement problem, where the objective is to place the poles of
$A-LC$ inside the unit circle. For that to be possible, the system needs to be
completely state observable [10]. The observability condition requires for the
column rank of the matrix

$\mathcal{O}=\begin{bmatrix}C\\\ CA\\\ CA^{2}\\\ \vdots\\\
CA^{n-1}\end{bmatrix}$ (7)

to be $n$. If this condition is satisfied, the system state can be fully
determined from the system inputs and outputs. One approach to solve the pole
placement problem is to choose a small positive number $\tau$ and calculate
the observer gain through:

$L=AC^{T}\Psi^{-1}\text{ where }\Psi=CC^{T}+\tau I_{p}$ (8)

Once $L$ is determined, the state estimate $\hat{x}$ can be used in place of
$x$ to design a controller. This means that the value function (4) can be
rewritten as follows:

$V_{obs}(x(t),\mathbf{u}_{t})=\sum_{k=t}^{\infty}\gamma^{k-t}\left\\{[C\hat{x}(k)-r(k)]^{T}Q[C\hat{x}(k)-r(k)]+u^{T}(k)Ru(k)\right\\}$
(9)

For the optimal control problem to be successfully solved, the system needs to
be controllable [10]. System (2) is said to be controllable when the matrix

$\mathcal{C}=\begin{bmatrix}B&AB&A^{2}B&\cdots&A^{n-1}B\end{bmatrix}$ (10)

known as the controllability matrix, has full row rank $n$. The
controllability condition implies that the system can be controlled to have
the desired behaviour within a finite number of steps, through an appropriate
choice of control function $u$.

For the infinite horizon LQT problem (9), a controller can be designed by
adopting the approach in [29]. Consider the augmented state

$X_{obs}(t)=\begin{bmatrix}\hat{x}(t)\\\ r(t)\end{bmatrix}$ (11)

and the corresponding system dynamics

$X_{obs}(t+1)=TX_{obs}(t)+B_{1}u(t).$ (12)

where $T=\begin{bmatrix}A&0\\\ 0&F\end{bmatrix}$ and
$B_{1}=\begin{bmatrix}B\\\ 0\end{bmatrix}$. Using (11), (9) can be rewritten
as

$V_{obs}(X(t),\mathbf{u}_{t})=\sum_{k=t}^{\infty}\gamma^{k-t}\\{X_{obs}^{T}(k)Q_{1}X_{obs}(k)+u^{T}(k)Ru(k)\\}$
(13)

where

$Q_{1}=\begin{bmatrix}C^{T}QC&-C^{T}Q\\\ -QC&Q\end{bmatrix}.$ (14)

The optimal policy is in feedback form of the augmented system state

$u^{*}(t)=-KX_{obs}(t).$ (15)

where $K$ is the optimal gain. $K$ is commonly obtained through the solution
of an Algebraic Riccati Equation (ARE) which can be a challenging task and
computationally expensive. An alternative approach to determining $K$,
particularly useful for large scale models, is through Reinforcement Learning
(RL) techniques, namely, by utilising the Bellman equation [21, 18, 17]. In
order to form the Bellman equation, the following lemma is needed.

###### Lemma 1

[29] For the optimal tracking problem with value function of the form (9) and
reference generated by $r_{obs}(t+1)=F_{obs}x_{obs}(t)$, then for any
stabilising policy of the form (15), the value function can be written in
quadratic form as

$V_{obs}(x(t),\mathbf{u}_{t})=V_{obs}(X(t))=\frac{1}{2}X_{obs}^{T}(t)PX(_{obs}t)$
(16)

for some matrix $P=P^{T}>0$.

Utilising the quadratic form of the value function, it can be rewritten as the
Bellman equation

$X_{obs}^{T}(t)PX_{obs}(t)=X_{obs}^{T}(t)Q_{1}X_{obs}(t)+u^{T}(t)Ru(t)+\gamma
X_{obs}^{T}(t+1)PX_{obs}(t+1).$ (17)

Consider a stabilising policy $\hat{K}^{*}$ in (17). A matrix $K$ is said to
be stabilising when the matrix $A-BK$ is stable, namely has eigenvalues with
strictly negative real parts [9]. Using the augmented system dynamics (12), we
obtain the Lyapunov equation

$P=Q_{1}+(\hat{K}^{*})^{T}R\hat{K}^{*}+\gamma(T+B_{1}\hat{K}^{*})^{T}P(T+B_{1}\hat{K}^{*}).$
(18)

Iterative solutions of the Lyapunov equation while updating the optimal
control estimate $\hat{u}^{*}=-\hat{K^{*}}X_{obs}$ can effectively lead to the
solution to the Optimal Tracking Problem. A step by step description of the RL
based solution to the LQT problem using state estimation through a state
observer is given in Algorithm 1.

Result: State estimate $\hat{x}(t)$, Optimal control gain $K^{*}$, optimal
controller $u^{*}(t)$

Input: State input and output matrices $A,B,C$, Weighting matrices for
performance index $Q,R$, Initial state estimate $\hat{x}(t_{0})$, small
positive number $\tau$, Initial reference signal $r(t_{0})$, Discount factor
$\gamma$, Initial control policy $K^{0}$, Error threshold $\epsilon$

  1. 1.

Create augmented system matrices $T=\begin{bmatrix}A&0\\\
0&F\end{bmatrix},B_{1}=\begin{bmatrix}B\\\ 0\end{bmatrix}$.

  2. 2.

Augment weighting matrix $Q_{1}=\begin{bmatrix}C^{T}QC&-C^{T}Q\\\
-QC&Q\end{bmatrix}.$

  3. 3.

Determine optimal control gain $K^{*}$:

    1. (i)

Set $j=0$.

    2. (ii)

Solve the LQT Lyapunov equation
$P^{j+1}=Q_{1}+(K^{j})^{T}RK^{j}+\gamma(T-B_{1}K^{j})^{T}P^{j+1}(T-B_{1}K^{j})$.

    3. (iii)

Compute the control gain estimate $K^{j+1}=(R+\gamma
B_{1}^{T}P^{j+1}B_{1})^{-1}$ $\gamma B_{1}^{T}P^{j+1}T$.

    4. (iv)

Set $j=j+1$.

Repeat steps $(ii)-(iv)$ while $|K^{j}-K^{j-1}|>\epsilon$.

  4. 4.

Set $K^{*}=K^{j}$.

  5. 5.

Determine and apply optimal controller to system while estimating the state
through the observer:

    1. (a)

Set $t=t_{0}$.

    2. (b)

Design augmented state $X_{obs}(t)=\begin{bmatrix}\hat{x}(t)\\\
r(t)\end{bmatrix}$

    3. (c)

Obtain system input $u^{*}(t)=-K^{*}X_{obs}(t)$

    4. (d)

Obtain output estimation $\hat{y}(t)=C\hat{x}(t)$.

    5. (e)

Observe system output $y(t)=Cx(t)$.

    6. (f)

Determine $L=AC^{T}\Psi^{-1}$ where $\Psi=CC^{T}+\tau I_{p}$.

    7. (g)

State estimate $\hat{x}(t+1)=A\hat{x}(t)+Bu^{*}(t)+L(y(t)-\hat{y}(t))$.

    8. (h)

System state update (not measurable) $x(t+1)=Ax(t)+Bu^{*}(t)$

    9. (i)

Reference signal update $r(t+1)=Fr(t)$

    10. (j)

Set $t=t+1$.

Repeat steps $(b)-(j)$ for the desired length of the experiment or simulation.

Algorithm 1 LQT control with state estimation through Luenberger observer

If the initial conditions given to the observer are the same or close in value
to the true initial conditions of the system, the observer estimates the state
very closely within a few time steps. The further away the initial values are,
the longer it takes for the observer to estimate the state correctly, so the
longer it takes for the system’s behaviour to be optimised, as the controller
is designed based on the state estimates.

While the observer is a powerful estimator of the state, it includes an
important assumption which is access to the system dynamics, namely matrices
A, B and C. This is an unrealistic assumption in many applications, where a
model of the process is either not available or not accurate enough. In the
next section, we explore another approach for control design with state
estimation.

##  4 Output Feedback using Input-Output data

Consider the performance index (4) to be minimised with $r$ generated as in
(2). Similarly to the previous section, we augment the state with the
reference and construct the augmented system dynamics:

$X(t+1)=\begin{bmatrix}x(t+1)\\\ r(t+1)\end{bmatrix}=\begin{bmatrix}A&0\\\
0&F\end{bmatrix}\begin{bmatrix}x(t)\\\ r(t)\end{bmatrix}+\begin{bmatrix}B\\\
0\end{bmatrix}u(t)=TX(t)+B_{1}u(t).$ (19)

The value function can be written in terms of $X$ as

$V(X(t),\mathbf{u}_{t})=\sum_{k=t}^{\infty}\gamma^{k-t}\\{X^{T}(k)Q_{1}X(k)+u^{T}(k)Ru(k)\\}$
(20)

where $Q_{1}$ defined as in (14). Assuming a feedback control policy, 1 allows
for 20 to be written as

$V(x(t),\mathbf{u}_{t})=V(X(t))=\frac{1}{2}X^{T}(t)P_{1}X(t)$ (21)

which gives rise to the Bellman equation

$X^{T}(t)P_{1}X(t)=X^{T}(t)Q_{1}X(t)+u^{T}(t)Ru(t)+\gamma
X^{T}(t+1)P_{1}X(t+1).$ (22)

Unlike the previous section, in addition to not having access to the internal
system state, we no longer assume knowledge of the system dynamics. Instead,
we estimate the state through a record of past inputs and outputs, under the
assumption that the system (2) is controllable and observable. Consider a time
horizon $N$ and recursively express the state $x(t)$ in terms of past inputs
$u(t-1),\dots,u(t-N)$

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle Ax(t-1)+Bu(t-1)$ (23)
$\displaystyle=$ $\displaystyle A^{2}x(t-1)+ABu(t-2)+Bu(t-1)$
$\displaystyle\dots$ $\displaystyle=$ $\displaystyle
A^{N}x(t-N)+\begin{bmatrix}B&AB&A^{2}B&\dots&A^{N-1}B\end{bmatrix}\begin{bmatrix}u(t-1)\\\
u(t-2)\\\ \vdots\\\ u(t-N)\end{bmatrix}$ $\displaystyle:=$ $\displaystyle
A^{N}x(t-N)+U_{N}\bar{u}(t-1,t-N).$

Similarly, the system output can be expressed recursively as

$y(t)=Cx(t)=CA^{N}x(t-N)+CU_{N}\bar{u}(t-1,t-N).$ (24)

Consider a sequence of N outputs

$\bar{y}(t-1,t-N)=\begin{bmatrix}y(t-1)\\\ y(t-2)\\\ y(t-3)\\\ \vdots\\\
y(t-N)\end{bmatrix}$ (25)

which can be expressed in terms of $x(t-N)$ and $\bar{u}(t-1,t-N)$ as

$\displaystyle\bar{y}(t-1,t-N)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}CA^{N-1}\\\ CA^{N-2}\\\ \vdots\\\ CA\\\
C\end{bmatrix}x(t-N)$ (26)
$\displaystyle+\begin{bmatrix}0&CB&CAB&\dots&CA^{N-2}B\\\
0&0&CB&\dots&CA^{N-3}B\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\dots&CB\\\ 0&0&0&\dots&0\end{bmatrix}\begin{bmatrix}u(t-1)\\\ u(t-2)\\\
\vdots\\\ u(t-N)\end{bmatrix}$ $\displaystyle:=$ $\displaystyle
W_{N}x(t-N)+D_{N}\bar{u}(t-1,t-N)$

Matrix $W_{N}$ has full rank for $N\geq K$ where $K$ is the rank of the
observability matrix (7). Since we assume our system to be observable, then
$K=n$. If $N$ is chosen to be equal to $n$, matrix $U_{N}$ is equal to the
controllability matrix $\mathcal{C}$ 10 and $W_{N}$ is equal to the
observability matrix $\mathcal{O}$ 7 flipped from top to bottom.

We aim to express the state $x(t)$ in terms of past inputs and outputs. Based
on 23, $x(t-N)$ needs to be replaced by an expression in terms of
$u(t-1),\dots,u(t-N)$ and $y(t-1),\dots,y(t-N)$. This can be achieved by
solving 24 for $x(t-N)$. To that end, consider the generalised Moore-Penrose
inverse [38, 39] of $W_{N}$, $W^{+}_{N}=(W^{T}_{N}W_{N})^{-1}W_{N}^{T}$. Then

$x(t-N)=W^{+}_{N}(\bar{y}(t-1,t-N)-D_{N}\bar{u}(t-1,t-N)).$ (27)

Considering that $r(t)=F^{N}r(t-N)$, the augmented state $X(t)$ can be written
as:

$\displaystyle X(t)$ $\displaystyle=$ $\displaystyle\begin{bmatrix}A^{N}&0\\\
0&F^{N}\end{bmatrix}\begin{bmatrix}x(t-N)\\\
r(t-N)\end{bmatrix}+\begin{bmatrix}U_{N}\\\ 0\end{bmatrix}\bar{u}(t-1,t-N)$
(28) $\displaystyle=$ $\displaystyle\begin{bmatrix}A^{N}&0\\\
0&F^{N}\end{bmatrix}\begin{bmatrix}W^{+}_{N}(\bar{y}(t-1,t-N)-D_{N}\bar{u}(t-1,t-N))\\\
r(t-N)\end{bmatrix}$ $\displaystyle+\begin{bmatrix}U_{N}\\\
0\end{bmatrix}\bar{u}(t-1,t-N)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}U_{N}-A^{N}W^{+}_{N}D_{N}&A^{N}W^{+}_{N}&0\\\
0&0&F^{N}\end{bmatrix}\begin{bmatrix}\bar{u}(t-1,t-N)\\\ \bar{y}(t-1,t-N)\\\
r(t-N)\end{bmatrix}$ $\displaystyle:=$ $\displaystyle
M\begin{bmatrix}\bar{u}(t-1,t-N)\\\ \bar{y}(t-1,t-N)\\\ r(t-N)\end{bmatrix}$

The above result establishes that the system state can be expressed in terms
of past inputs and outputs with the help of matrix $M$. Let
$Z(t)=\begin{bmatrix}X(t)\\\ u(t)\end{bmatrix}$ and consider the quadratic
form of the value function 21. It can be rewritten as

$V(Z(t))=\frac{1}{2}Z^{T}(t)HZ(t)$ (29)

with $H=\begin{bmatrix}Q_{1}+\gamma T^{T}P_{1}T&\gamma T^{T}P_{1}B_{1}\\\
\gamma B_{1}^{T}P_{1}T&R+\gamma B_{1}^{T}P_{1}B_{1}\end{bmatrix}$. Using 28 in
29 we obtain

$\displaystyle V(Z(t))$ $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{bmatrix}X(t)\\\
u(t)\end{bmatrix}^{T}H\begin{bmatrix}X(t)\\\ u(t)\end{bmatrix}$ (30)
$\displaystyle=$
$\displaystyle\frac{1}{2}\begin{bmatrix}M\begin{bmatrix}\bar{u}(t-1,t-N)\\\
\bar{y}(t-1,t-N)\\\ r(t-N)\end{bmatrix}\\\
u(t)\end{bmatrix}^{T}H\begin{bmatrix}M\begin{bmatrix}\bar{u}(t-1,t-N)\\\
\bar{y}(t-1,t-N)\\\ r(t-N)\end{bmatrix}\\\ u(t)\end{bmatrix}$

Set $\bar{Z}(t)=\begin{bmatrix}\bar{u}(t-1,t-N)\\\ \bar{y}(t-1,t-N)\\\
r(t-N)\\\ u(t)\end{bmatrix}$ and $\bar{H}=\begin{bmatrix}M&0\\\
0&I_{m}\end{bmatrix}^{T}H\begin{bmatrix}M&0\\\ 0&I_{m}\end{bmatrix}$. Then we
can rewrite 29 as

$\displaystyle V(\bar{Z}(t))$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bar{Z}^{T}(t)\bar{H}\bar{Z}(t)$ $\displaystyle=$
$\displaystyle{\small\frac{1}{2}\begin{bmatrix}\bar{u}(t-1,t-N+1)\\\
\bar{y}(t,t-N+1)\\\ r(t-N+1)\\\
u(t)\end{bmatrix}^{T}\begin{bmatrix}H_{\bar{u}\bar{u}}&H_{\bar{u}\bar{y}}&H_{\bar{u}r}&H_{\bar{u}u}\\\
H_{\bar{y}\bar{u}}&H_{\bar{y}\bar{y}}&H_{\bar{y}r}&H_{\bar{y}u}\\\
H_{r\bar{u}}&H_{r\bar{y}}&H_{rr}&H_{ru}\\\
H_{u\bar{u}}&H_{u\bar{y}}&H_{ur}&H_{uu}\end{bmatrix}\begin{bmatrix}\bar{u}(t-1,t-N+1)\\\
\bar{y}(t,t-N+1)\\\ r(t-N+1)\\\ u(t)\end{bmatrix}}$

The decomposition of $\bar{H}$ in blocks of the form $H_{ij}$ is according to
the dimensions of each set of $i,j$. For example block $H_{\bar{u}\bar{u}}$
will be a $Nm\times Nm$ matrix and block $H_{u\bar{y}}$ will be a $m\times Np$
matrix.

Due to the quadratic form of the value function, the optimality condition
$\frac{\partial V}{\partial u}=0$ is sufficient to ensure its minimisation.
Solving for $u(t)$ yields

$u(t)=-H_{uu}^{-1}(H_{u\bar{u}}\bar{u}(t-1,t-N)+H_{u\bar{y}}\bar{y}(t-N,t-N)+H_{ur}r(t-N)).$
(32)

If the system model is available, 32 can be used to directly determine the
optimal LQT controller using input-output information. The kernel matrix
$\bar{H}$ can be directly determined using

$\displaystyle\bar{H}$ $\displaystyle=\begin{bmatrix}M&0\\\
0&I_{m}\end{bmatrix}^{T}H\begin{bmatrix}M&0\\\
0&I_{m}\end{bmatrix}=\begin{bmatrix}U_{N}-A^{N}W^{+}_{N}D_{N}&A^{N}W^{+}_{N}&0&0\\\
0&0&F^{N}&0\\\ 0&0&0&I_{m}\end{bmatrix}^{T}$
$\displaystyle\begin{bmatrix}Q_{1}+\gamma T^{T}P_{1}T&\gamma
T^{T}P_{1}B_{1}\\\ \gamma B_{1}^{T}P_{1}T&R+\gamma
B_{1}^{T}P_{1}B_{1}\end{bmatrix}\begin{bmatrix}U_{N}-A^{N}W^{+}_{N}D_{N}&A^{N}W^{+}_{N}&0&0\\\
0&0&F^{N}&0\\\ 0&0&0&I_{m}\end{bmatrix}$

where $P_{1}$ is obtained by solving the standard LQT problem as in 18.

We are interested in the case when a reliable model of the system is not
available, meaning matrices $A,B$ and $C$ are unknown. In this case, the
kernel matrix $\bar{H}$ needs to be estimated strictly from data. To that end,
consider the Bellman equation derived from the quadratic form 29

$\displaystyle\bar{Z}(t)^{T}\bar{H}\bar{Z}(t)$ $\displaystyle=$
$\displaystyle(y(t)-r(t))^{T}Q(y(t)-r(t))+u^{T}(t)Ru(t)$ (34)
$\displaystyle+\gamma\bar{Z}(t+1)^{T}\bar{H}\bar{Z}(t+1).$

To obtain an approximation formula for $\bar{H}$, it needs to be isolated from
the quadratic form. This can be achieved by firstly vectorising 34:

$\displaystyle vec\left(\bar{Z}^{T}(t)\bar{H}\bar{Z}(t)\right)$
$\displaystyle=$ $\displaystyle(y(t)-r(t))^{T}Q(y(t)-r(t))+u^{T}(t)Ru(t)$ (35)
$\displaystyle+\gamma vec\left(\bar{Z}^{T}(t+1)\bar{H}\bar{Z}(t+1)\right)$

where the $vec$ operator stacks the columns of a matrix to create a single
column. We then apply the "vector trick" associated with the Kronecker product

$\displaystyle\left(\bar{Z}^{T}(t)\otimes\bar{Z}^{T}(t)\right)vec(\bar{H})$
$\displaystyle=$ $\displaystyle(y(t)-r(t))^{T}Q(y(t)-r(t))+u^{T}(t)Ru(t)$ (36)
$\displaystyle+\gamma\left(\bar{Z}^{T}(t+1)\otimes\bar{Z}^{T}(t+1)\right)vec(\bar{H}).$

To solve 36, we may employ the Policy Iteration (PI) or Value Iteration (VI)
algorithm [18]. In the standard Reinforcement Learning framework, where the
problem is modelled as a Markov Decision Process with discrete state and
action spaces, there are two main differences between the algorithms: Firstly,
the PI assumes that the initial policy is an admissible one while VI has no
such requirement. The second difference lies in the way the optimal policy is
built. VI focuses on optimising the value function over all possible actions
and extracting the optimal policy at the end. PI updates the value function
based on the current policy estimate and extracts a new, better policy
afterwards, iteratively until the optimal is reached. In the case of
continuous state and action spaces, such as our control problem, we will focus
on the VI algorithm, mainly due to the lack of need for an initial stabilising
controller (admissible policy). When the system dynamics are not known, the
design of such a controller becomes very challenging, though some work has
been done in recent literature on obtaining stabilising controllers through
measured data [40]. The algorithm iterates between a Policy Evaluation and a
Policy Improvement step. We start with an arbitrarily chosen initial kernel
matrix $\bar{H}^{0}$ that obtains an initial policy $u^{0}(t)$. We then
iterate between the following two steps:

  1. 1.

Policy Evaluation

$\displaystyle\left(\bar{Z}^{T}(t)\otimes\bar{Z}^{T}(t)\right)vec(\bar{H}^{i+1})$
$\displaystyle=$ $\displaystyle(y(t)-r(t))^{T}Q(y(t)-r(t))+u^{T}(t)Ru(t)$
$\displaystyle+\gamma\left(\hat{\bar{Z}}^{T}(t+1)\otimes\hat{\bar{Z}}^{T}(t+1)\right)vec(\bar{H^{i}})$

  2. 2.

Policy Improvement

$\displaystyle u^{i+1}(t)$ $\displaystyle=$
$\displaystyle-(H_{uu}^{i+1})^{-1}(H_{u\bar{u}}\bar{u}^{i+1}(t-1,t-N)+H_{u\bar{y}}^{i+1}\bar{y}(t-1,t-N)$
$\displaystyle+H_{ur}^{i+1}r(t-N))$

The Policy evaluation step can be solved through the Least Squares (LS)
algorithm using measured data $\bar{Z}(t)$ and calculating the corresponding
cost term $(y(t)-r(t))^{T}Q(y(t)-r(t))+u^{T}(t)Ru(t)$. It is worth noting that
for the discounted future cost term
$\gamma\left(\hat{\bar{Z}}^{T}(t+1)\otimes\hat{\bar{Z}}^{T}(t+1)\right)vec(\bar{H^{i}})$,
we are not using the measured input $u(t+1)$ but instead the estimated optimal
input $u^{i}(t+1)$, namely
$\hat{\bar{Z}}(t+1)=\begin{bmatrix}\bar{u}(t,t-N+1)&\bar{y}(t,t-N+1)&r(t-N+1)&u^{i}(t+1)\end{bmatrix}^{T}$.
In RL terms, this yields the expected future reward based on the current
policy. Matrix $\bar{H}$ is an $((N+1)m+(N+1)p)\times((N+1)m+(N+1)p)$
symmetric matrix which means that its determination is a problem that requires
a minimum of $((N+1)m+(N+1)p)\times((N+1)m+(N+1)p)/2$ data points. In practice
many more data points are usually needed, especially when dealing with more
complicated systems with larger state and action spaces. When the data are
highly correlated, which is often the case with sequential input output data
from a control system, the inversion step in the LS algorithm becomes
challenging. This can be rectified with a regularized LS approach, which
includes an appropriate regularisation parameter $\mu$. A step by step
description of the method introduced in this section is given in Algorithm 2.

Result: Optimal kernel matrix estimate $H^{*}$, optimal controller estimate
$u^{*}(t)$

Input: Weighting matrices for performance index $Q_{1},R$, regularisation
parameter $\mu$, Reference signal $r(t)$, Discount factor $\gamma$, Initial
kernel matrix $H^{0}$, Error threshold $\epsilon_{rl}$, time horizon $N$

  1. 1.

Collect $M\gg((N+1)m+(N+1)p)\times((N+1)m+(N+1)p)/2$ consecutive data tuples
$\left\\{y(t),u(t)\right\\}$.

  2. 2.

Create vectors $\bar{u}(t-1,t-N),\bar{y}(t-1,t-N)$ and
$\bar{Z}(t)=\begin{bmatrix}\bar{u}(t-1,t-N)&\bar{y}(t-1,t-N)&r(t-N)&u(t)\end{bmatrix}^{T}$.

  3. 3.

For each $\bar{Z}(t)$, calculate the Kronecker product
$K_{\bar{Z}}(t):=\bar{Z}^{T}(t)\otimes\bar{Z}^{T}(t)$.

  4. 4.

Set $i=0$.

  5. 5.

For t=N to M:

    1. (i)

Determine
$u^{i}(t+1)=-(H_{uu}^{i+1})^{-1}(H_{u\bar{u}}\bar{u}^{i+1}(t-1,t-N)+H_{u\bar{y}}^{i+1}\bar{y}(t-1,t-N)+H_{ur}^{i+1}r(t-N))$

    2. (ii)

Create the updated augmented state
$\hat{\bar{Z}}(t+1)=\begin{bmatrix}\bar{u}(t,t-N+1)&\bar{y}(t,t-N+1)&r(t-N+1)&u^{i}(t+1)\end{bmatrix}^{T}$
and the corresponding Kronecker product
$K_{\hat{\bar{Z}}}(t+1)=\hat{\bar{Z}}^{T}(t+1)\otimes\hat{\bar{Z}}^{T}(t+1)$

    3. (iii)

Compute the future expected cost term $c(t)=X^{T}(t)Q_{1}X(t)$
$+u^{T}(t)Ru(t)+\gamma K_{\hat{\bar{Z}}}(t+1)vec(H^{i})$

  6. 6.

Calculate the sums $L=\sum_{t=N+1}^{M}K^{T}_{Z}(t)K_{Z}(t)$ and
$R=\sum_{t=N+1}^{M}K^{T}_{Z}(t)c(t)$ .

  7. 7.

Regularise matrix L: $L_{r}=L+\mu I$

  8. 8.

Produce new H estimate $H^{i+1}=L_{r}^{-1}R$

  9. 9.

Set $i=i+1$

Repeat steps (5)-(9) until $\|H^{i}-H^{i+1}\|\leq\epsilon_{rl}$ or for a
predetermined number of iterations.

Algorithm 2 LQT output-feedback control using input-output measured data

From a software point of view, the implementation of Algorithm 2 can become
very computationally expensive with larger system dimensions. There are a few
ways to work around this problem. The construction of matrix $L$ in step (6)
consists of the Kronecker product terms obtained directly from the dataset.
This means that $L$ can be determined only once outside the training
iterations. Similarly for its regularised version $L_{r}$ and more importantly
for the inverse $L_{r}^{-1}$. Matrix inversion is a very expensive operation
and an approximation of the inverse might be necessary. Specifically, we chose
to calculate the inverse with the use of the Cholesky decomposition [41].

##  5 Parameter Tuning through Bayesian Optimisation

There are multiple aspects that influence the performance of a controller,
both in the model-based scenario and especially in the model-free, data-driven
scenario. In the model-based case, these aspects concern mainly the control
objective design. In the case of data-driven control, the resulting
controller’s performance is influenced as much by characteristics of the data
as by the definition of the optimisation problem. In the next section we
discuss the setup of our simulations so as to generate informative data, that
can be used to effectively train a controller. In this section, we focus on
optimising the design of the performance index, which dictates the control
objective.

Apart from the system dynamics (2), the main parameters influencing the
optimisation are the weighting matrices $Q$ and $R$ in the performance index
(2). However, there is no systematic way to choose the optimal values for
these matrices, as their direct relationship to the controller performance
cannot be expressed by a known function. Hence, they are commonly chosen
empirically by the engineer, or for simpler problems, exhaustive search may be
performed to obtain the best choice. For our system, an exhaustive search
would not be feasible, especially in the model-free case where the
calculations become vastly more computationally expensive due to the system
augmentation. Motivated by this we explore Bayesian Optimisation [32, 33, 42],
which is a powerful optimization tool in which uncertainty over the objective
function is typically represented as a Gaussian process (GP)[43].

In particular, we want to design a fitness function which quantifies the
effect of $Q$ and $R$ on the controller’s performance. We consider a
controller to have good performance if it tracks the reference signal closely,
whilst not applying costly inputs to the system. A choice for the fitness
function could be the performance index itself as defined in 2. There are two
main drawbacks to this choice. Firstly, an infinite sum is not feasible, which
means it would need to be replaced by a finite sum corresponding to the length
of the experiment. Secondly, using (2) as the fitness function to be
minimised, may imply that using very small values for $Q$ and $R$ is
sufficient to have a good performance. Such a choice however, would only
ensure that the sum of discounted costs is small, without necessarily
guaranteeing good tracking performance and low input costs. Hence, we define
the simpler fitness function

$F(\theta)=\sum_{t=0}^{l}\gamma^{t}(\|y(t)-r(t)\|^{2}+\|u(t)\|^{2})$ (37)

where $\|\cdot\|$ denotes the Euclidian norm and $l$ is the length of the
experiment. $\theta\in\mathcal{D}$ represents the parameterised choice of $Q$
and $R$ sampled from a bounded search domain $\mathcal{D}$. Considering that
$Q$ and $R$ are symmetric positive definite matrices, we can chose them to be
diagonal matrices with all positive entries, as is common practice for
simplicity. This means that $\theta$ can be written as
$\theta=\\{q_{1},\dots,q_{p},r_{1},\dots,r_{m}\\}$ and
$\mathcal{D}\subset\mathbb{R}_{+}^{p+m}$. This choice for $F$ allows us to
investigate the tracking and input costs as an indirect result of $Q$ and $R$,
with the weighting matrices being used in the Bellman equation to obtain the
optimal policy and thus the optimal trajectory.

We want to minimise the fitness function with respect to $\theta$. $F$ can be
modelled through a GP prior as

$F(\theta)\sim\mathcal{GP}(\mu(\theta),k(\theta,\theta^{\prime}))$ (38)

The GP is characterised by a mean $\mu(\theta)$ and a covariance or kernel
function $k(\theta,\theta^{\prime})$.

Consider a sample of $l$ points $\\{x_{1},\dots,x_{l}\\}$ and calculate the
covariance for each pair $\\{x_{i},x_{j}\\}$. The mean for this initial sample
may be chosen as $\mu_{0}=0$ which is a common choice for the mean when no
prior information is available. For the kernel, we choose one of the most
popular and simpler kernel function which is the Square Exponential or
Gaussian kernel

$k_{0}(\theta,\theta^{\prime})=\exp(-\|\theta-\theta^{\prime}\|^{2})\text{
where
}\|\theta-\theta^{\prime}\|^{2}=\sum_{i=1}^{p+m}(\theta_{i}-\theta_{i}^{\prime})^{2}.$
(39)

All the above yield the prior distribution

$F(x_{1:l})\sim\mathcal{N}(\mu_{0}(x_{1:l}),k_{0}(x_{1:l}))$ (40)

where $F(x_{1:l})=\\{F(x_{1}),\dots,F(x_{l})\\}$, $\mu_{0}(x_{1:l})=0_{l}$ and
$k_{0}(x_{1:l})=\\{k(x_{1},x_{1}),\dots,$
$k(x_{1},x_{l});\dots;k(x_{l},x_{1}),\dots,k(x_{l},x_{l})\\}$. Once the prior
distribution has been determined for the initial sample of points, we can
evaluate $F$ at a new point $x$ and compute the conditional distribution of
$F(x)$ using Bayes’ rule

$F(x)|F(x_{1:l})\sim\mathcal{N}(\mu_{l}(x),k_{l}(x))$ (41)

where
$\mu_{l}(x)=k_{0}(x,x_{1:l})k_{0}(x_{1:l},x_{1:l})^{-1}(F(x_{1:l})-\mu_{0}(x_{1:l}))+\mu_{0}(x)$
and
$k_{l}(x)=k_{0}(x,x)-k_{0}(x,x_{1:l})k_{0}(x_{1:l},x_{1:l})^{-1}k_{0}(x_{1:l},x)$.
Distribution (41) is also known as the posterior probability distribution. The
more points that are sampled, the more informed the posterior becomes.

Returning to the optimisation problem at hand, we are looking for a choice of
$Q$ and $R$ that minimise $F$. The tool used in Bayesian Optimisation to find
such a point is known as the acquisition function. Acquisition functions
provide a candidate for the best point to sample next, based on the current
estimate for the distribution. There are a few different choices for
acquisition functions, the most common including Expected Improvement,
Probability of Improvement and Gaussian Process Upper Confidence Bound. For
our application, we will be using the latter, defined as follows

$a_{UCB}(x;x_{1:l})=\mu_{l}(x)-\epsilon k_{l}(x)$ (42)

where $\epsilon$ is a parameter balancing exploration against exploitation.
The next point to sample $V$ is chosen as $x_{next}=\arg\min a_{UCB}$. The new
point can be used to update the posterior distribution, while resulting in a
more informed acquisition function that will in turn suggest a better point to
sample next. The process is repeated for a predefined number of steps and the
optimal is extracted at the end if $M$ points have been sampled in total, as

$x_{best}=\arg\min_{x}\mu_{M}(x).$ (43)

##  6 Simulation setup and results

In this section we perform numerical simulations to investigate the behaviour
of the extruder heating system introduced in Section 2 after a controller is
introduced to the system. We assume partial access to the system state in the
form of measurements and explore both the case of full knowledge of the system
dynamics and the case where a system model is unknown. For our simulations, we
need a model of the extruder dynamics which will be used for the design of the
model-based controller and to generate data to train the data-driven
controller. In [37], system identification methods were used to obtain the
following state and input matrices for a BAAM system:

$\displaystyle\footnotesize A$ $\displaystyle=$
$\displaystyle\begin{bmatrix}0.992&0.0018&0&0&0&0\\\
0.0023&0.9919&0.0043&0&0&0\\\ 0&-0.0042&1.0009&0.0024&0&0\\\
0&0&0.0013&0.9979&0&0\\\ 0&0&0&0&0.9972&0\\\ 0&0&0&0&0&0.9953\end{bmatrix}$
$\displaystyle B$ $\displaystyle=$
$\displaystyle\begin{bmatrix}1.0033&0&0&0&0&0&-0.2175\\\
0&1.0460&0&0&0&0&-0.0788\\\ 0&0&1.0326&0&0&0&-0.0020\\\
0&0&0&0.4798&0&0&-0.0669\\\ 0&0&0&0&0.8882&0&0.1273\\\
0&0&0&0&0&1.1699&-0.1792\end{bmatrix}$

In addition, we design an output matrix for the system as follows:

$C=\begin{bmatrix}0.992&0.00018&0&0&-0.0001&0\\\ 0.0023&1.3&0.0043&0&0&0\\\
0&-0.0042&1.0109&0.0024&0&0.201\\\ 0&0&0.0013&0.989&0.00031&0.64\\\
0&0&0&0&0.923&0.3\end{bmatrix}$ (45)

We can easily verify that the system described by the above matrices is
controllable and observable. As highlighted in the previous section, the most
important design parameters for the optimisation are the weighting matrices Q
and R. Before performing Bayesian Optimisation to obtain the best values, we
chose as a baseline identity matrices of appropriate dimensions, namely
$Q=I_{p}$ and $R=I_{m}$.

In BAAM systems, the multiple heaters located across the barrel heat the
material in steadily increasing temperatures until the material reaches a
suitable temperature for extrusion. For our simulations, we assume the polymer
used is PLA, whose melting temperature ranges between 160 and 180 °C. With
this in mind, we set the reference signal to be
$r(t)=\\{150,160,170,175,180\\},\forall t\geq t_{0}$. This means that the
reference generating matrix is $F=I_{5}$. This reference will be used
throughout the simulations, both for the model-based controller and the data-
driven one.

###  6.1 State Observer based Output Feedback

We first focus on the approach introduced in Section 3, starting with the
design of the State Observer, for which we choose $\tau=0.002$. For the
controller design, we use a discount factor $\gamma=0.99$ and an error
threshold for convergence $\epsilon=0.01$. We arbitrarily initialise the state
estimate with $\hat{x}(t_{0})=\\{50,\dots,50\\}$ and the control policy with
random normally distributed numbers as $K^{0}=\\{k_{1},\dots,k_{n+p}\\}$ where
$k_{i}\sim\mathcal{N}(0,10)$. Finally, even though we assume that the system
state is not directly available, we need to initialise it so that it can be
updated and measured during the simulations. That starting point is chosen to
be $x(t_{0})=\\{20,\dots,20\\}$. Table 1 summarises all parameter values used
for the observer-based simulations.

Table of parameters - Observer based output feedback
---
Parameter | Symbol | Value
Tracking reference | r | $\left\\{150,160,170,175,180\right\\}$
Initial state | $x(t_{0})$ | $\left\\{20,\dots,20\right\\}$
Initial observer state | $\hat{x}(t_{0})$ | $\left\\{50,\dots,50\right\\}$
Observer pole | $\tau$ | $0.002$
Error threshold | $\epsilon$ | $0.01$
Initial control gain | $K^{0}=\\{k_{i}\\}_{i=1}^{n+p}$ | $k_{i}\sim\mathcal{N}(0,10)$
Table 1: Parameters chosen for observer-based control implementations. The
information listed includes the parameter names, the corresponding symbol and
the value that was chosen.

Figure 2 shows the output trajectories obtained as a result of applying the
model-based controller for 100 time steps, as described in Algorithm 1. The
controller manages to bring all measurements within $0.03$°C of the optimal.
The measured temperatures arrive within $0.1$ degrees of the target within
only 10 steps and remain there. The values that the measurements finally
converge to are $y^{*}=\\{149.9798,159.9932,169.9795,174.9815,179.9859\\}$.
The error between the final measurement values and the reference, calculated
through the Euclidean norm, is $\|r-y^{*}\|_{2}=0.0376$. The value of the
performance index obtained by this controller, approximated as the discounted
sum of the first 100 cost terms is 183,362.5. When calculated over 1000 steps,
it reaches 183.436.2. The performance index trajectory for 100 time steps can
be seen in Figure 3(a). The good performance of the controller is due to the
fact that the observer approximates the system state very closely, thanks to
the knowledge of the system dynamics. The estimation error, calculated as the
Euclidean norm $e(t)=\|\hat{x}(t)-x(t)\|_{2}$, is decreasing with each time
step, reaching the value $0.008$ in 100 steps. It is worth noting that the
closer the initial state estimate is to the initial state, the faster the
error decreases. This decreasing trend can be seen in Figure 3(b) where we
plot the estimation error for the first 100 time steps.

Figure 2: Observer-based Output Feedback. The system output achieves
temperatures very close to the optimal within a few steps.

(a) Sum of discounted cost terms for 100 time steps.

(b) Estimation error for the system state by the observer.

Figure 3: Performance index and state estimation error for the Observer-based
Output feedback.

We perform Bayesian Optimisation to obtain better choices for the weighting
matrices $Q$ and $R$ and study the effect they have on our results. For
simplicity, as is common practice, we assume these matrices are diagonal and
positive definite. Through some trial and error attempts, we find that the
best range of values for the diagonal entries is positive numbers $\leq 1$. To
avoid numerical instability, we assume a lower limit for these values to be
$0.01$. Additionally, we want to put a higher emphasis on the tracking of the
reference signal. This is expressed by choosing appropriate upper limits for
the entries of each matrix. Specifically, we chose $1$ as the upper limit for
$Q$ and $0.4$ as the upper limit for $R$.

We run the Bayesian Optimisation procedure as described in Section 5 with 50
initial randomly sampled points for 100 iterations. The optimal values
obtained for the matrices Q and R rounded up to 3 decimal points are

$Q^{*}=\begin{bmatrix}0.943&0&0&0&0\\\ 0&0.762&0&0&0\\\ 0&0&0.542&0&0\\\
0&0&0&0.420&0\\\ 0&0&0&0&0.514\end{bmatrix}$ (46)

and

$R^{*}=\begin{bmatrix}0.300&0&0&0&0&0&0\\\ 0&0.270&0&0&0&0&0\\\
0&0&0.281&0&0&0&0\\\ 0&0&0&0.092&0&0&0\\\ 0&0&0&0&0.054&0&0\\\
0&0&0&0&0&0.269&0\\\ 0&0&0&0&0&0&0.318\end{bmatrix}$ (47)

Using the optimised values $Q^{*}$ and $R^{*}$ we obtain the updated
trajectories seen in Figure 4(a). The values that the outputs converge to are:
$\\{149.9940,159.9946,169.9843,174.9873,179.9834\\}$ which means they all
arrive within $0.02$°C of the optimal (cf $0.03$°C for the previous version).
The error of the final values from the optimal is $\|r-y^{*}\|_{2}=0.0274$,
verifying that the tracking performance is improved. Additionally, the
performance index is significantly lower in this case, as can be seen in
Figure 4(b). Specifically, the performance index obtains the value 76,060.36
after 100 steps and converges to 76,072.6 over 1000 steps. This means that the
updated weighting matrices result in a $58.5\%$ reduction of the performance
index.

(a) System output trajectories

(b) Comparison of performance index before and after using Bayesian
Optimisation

Figure 4: Output trajectories and performance index comparison for the
Observer-based Output feedback where the weighting matrices $Q$ and $R$ were
obtained through Bayesian Optimisation. The output trajectories are similar to
the default case while the performance index is minimised more effectively
through the Bayesian Optimisation approach.

###  6.2 Data generation

When designing and training model-free, data-driven controllers, much like any
learning algorithm, it is important to have data that is rich in information.
In the case of the algorithm introduced in Section 4, the iterative algorithm
learns the kernel matrix $H$ through input-output tuples. Matrix $H$ includes
information of the system dynamics as seen in Equation 4, meaning that
learning $H$ is equivalent to learning the system model. This means that the
problem can also be viewed as a system identification problem [44].

An important advantage of working with simulated data generated from an
assumed model is that we can ensure the data is appropriate to train our
algorithms. The way to achieve that is by using a well designed input signal
to the system. Firstly, we need the system output to stay bounded and not
diverge towards infinite values, in which case the data would not be usable.
To that end, we can design the input to be in feedback form

$u(t)=-K_{data}x(t)$ (48)

where the feedback gain $K_{data}$ is stabilising. In the case where the
training data is obtained from a real process, this step is not applicable,
however, the process would be tuned from the manufacturer to avoid extreme
divergence. In the case of the system introduced in 6 and 45, we can use the
Matlab Control System Toolbox to obtain a stabilising gain. In particular, two
of the most popular approaches to stabilising a system’s behaviour are the
Pole Placement method and the Linear Quadratic Regulator (LQR). Using the LQR
approach, we obtain the gain:

$K_{data}=\begin{bmatrix}[r]0.7395&-0.0076&-0.0003&-0.0264&0.0194&-0.0170\\\
-0.0076&0.7430&0.0031&-0.0093&0.0068&-0.0060\\\
-0.0003&-0.0033&0.7599&0.0021&0.0002&-0.0002\\\
-0.0126&-0.0042&0.0016&1.0971&0.0092&-0.0079\\\
0.0171&0.0058&0.0002&0.0170&0.8179&0.0108\\\
-0.0198&-0.0067&-0.0002&-0.0193&0.0143&0.6823\\\
-0.1525&-0.0519&-0.0018&-0.1412&0.1091&-0.0977\end{bmatrix}$ (49)

Secondly, in order to ensure the output data is varied enough to efficiently
train the control algorithm, we require the input signal to not be ‘too
simple’. This requirement can be fulfilled by the addition of an appropriate
probing noise signal to the input, meaning that it will now be of the form

$u(t)=-K_{data}x(t)+\omega_{pr}.$ (50)

We define the probing noise as a sum of a Gaussian noise signal and sinusoidal
signals, of varying ranges and frequencies. Specifically, to ensure great
variety in the data, we chose ranges between 10 and 100. The choice of
sampling frequencies is more particular, where an appropriate choice is the
system’s bandwidth, i.e. the part of frequency domain where most of signal’s
energy is contained. Formally, a system’s bandwidth is the frequency at which
the magnitude of the system’s response drops below 0.707 (-3dB). In the case
of multi-input multi-output (MIMO) systems, the bandwidth is the frequency
where the maximum singular value of the frequency response crosses 0.707 from
below [45]. In Matlab, this value can be easily obtained for any system using
the sigma function. The bandwidth for our system is 1.65, so we use it as the
maximum frequency from which to sample sinusoidal signals. Finally, we add a
normally distributed random vector to our noise signal, to add an extra
element of randomness. We consider the variance for this vector to be
$\sigma=1.5$. Bringing all the above elements together, we have

$\displaystyle\omega_{pr}$ $\displaystyle=$
$\displaystyle\omega_{1}+100\sin(\omega_{2}t)+90\sin(\omega_{3}t)+80\sin(\omega_{4}t)+70\sin(\omega_{5}t)$
(51)
$\displaystyle+60\sin(\omega_{6}t)+50\sin(\omega_{7}t)+40\sin(\omega_{8}t)+30\sin(\omega_{9}t)$
$\displaystyle+20\sin(\omega_{10}t)+10\sin(\omega_{11}t)$

where

$\omega_{1}\sim\mathcal{N}(0,\sqrt{\sigma}I_{7}),\;\omega_{2}=(16.5,\dots,16.5)\in\mathbb{R}^{7},\text{
and }\omega_{3},\dots\omega_{11}\sim\mathcal{U}_{7}(0,1.65).$ (52)

###  6.3 Data-driven Output Feedback

Using the input signal designed in the previous section, we can obtain data
tuples $\\{y(t),u(t)\\}$ to train the data-driven RL based controller
introduced in Section 4. In particular, we generate $M=15000$ samples and set
a maximum number of iterations to be 1000 in case the algorithm does not
converge according to the error threshold chosen as $\epsilon_{rl}=0.001$. As
the system is fully observable, we choose the time horizon $N$ for the
augmented data tuples to be $N=6$. We choose the discount factor to be
$\gamma=0.99$ and the regularisation parameter to be $\mu=0.0001$. Finally we
initialise the kernel matrix $H$ as the identity matrix of appropriate
dimensions, namely $H^{0}=I_{(N+1)m+(N+1)p}$. We run two sets of simulations,
one with the weighting matrices $Q,R$ chosen to be identity matrices of
appropriate dimensions, and one using the values for $Q$ and $R$ obtained
through Bayesian Optimisation.

Table of parameters - Model based methods
---
Parameter | Symbol | Value
Universal parameters
Tracking reference | r | $\left\\{150,160,170,175,180\right\\}$
Initial state | $x(t_{0})$ | $\left\\{50,\dots,50\right\\}$
Discount factor | $\gamma$ | $0.99$
Error threshold | $\epsilon_{rl}$ | $0.001$
Regularisation parameter | $\mu$ | $0.0001$
Time horizon | $N$ | $6$
Initial kernel matrix | $H^{0}$ | $I_{(N+1)m+(N+1)p}$
Table 2: Parameters chosen for data-driven RL based implementations. The
information listed includes the parameter names, the corresponding symbol and
the value that was chosen.

Firstly, we focus on the case of the arbitrarily chosen weighting matrices $Q$
and $R$. We chose them to be identity matrices of appropriate dimensions,
namely $Q=I_{5}$ and $R=I_{7}$. Figure 5 shows the output trajectories
obtained through the data-driven controller described in Algorithm 2, and the
corresponding performance index. The measurements converge to values within
$1$°C of the optimal, specifically to the vector
$y^{*}=\\{150.0857,160.1181,171.0027,175.7926,180.1416\\}$. The error between
these values and the optimal ones, calculated through the Euclidean norm, is
$\|r-y^{*}\|_{2}=1.2942$. The performance index, defined as the infinite sum
of discounted costs, is estimated by a large finite sum which for 1000 steps
is 1,074,985. We plot it for the first 100 steps, in which period the sum has
obtained the value 1,065,077.

(a) System output trajectories

(b) Performance Index

Figure 5: Output trajectories and the performance index the data-driven Output
feedback algorithm using input-output data, where the weighting matrices $Q$
and $R$ were chosen as identity matrices of appropriate dimensions.

We now perform Bayesian Optimisation on $Q$ and $R$ using the same upper and
lower bounds for the search space as in Section 6.1, namely we chose $0.1$ as
the lower bound for both matrices and chose $1$ as the upper bound for $Q$ and
$0.4$ as the upper bound for $R$. We generate a grid of 50 random initial
samples and then train the algorithm for 100 iterations. The optimal values
found for $Q$ and $R$ rounded to 3 decimal points are

$Q^{*}=\begin{bmatrix}0.174&0&0&0&0\\\ 0&0.056&0&0&0\\\ 0&0&0.010&0&0\\\
0&0&0&0.010&0\\\ 0&0&0&0&0.160\end{bmatrix}$ (53)

and

$R^{*}=\begin{bmatrix}0.145&0&0&0&0&0&0\\\ 0&0.389&0&0&0&0&0\\\
0&0&0.020&0&0&0&0\\\ 0&0&0&0.116&0&0&0\\\ 0&0&0&0&0.099&0&0\\\
0&0&0&0&0&0.316&0\\\ 0&0&0&0&0&0&0.010\end{bmatrix}.$ (54)

With these values for $Q$ and $R$, Algorithm 2 obtains a controller that
produces the output seen in Figure 6(a). The measurements converge to the
vector $y^{*}=\\{150.069160.154170.011175.199180.176\\}$ after 100 steps,
namely they converge within $0.2$°C of the optimal. The error associated with
these values is $\|r-y^{*}\|=0.3154$. This means that a tracking performance
is significantly improved compared to the previous case. In addition, the
overall cost is much lower, as can be seen by the comparison of the two
performance indices, before and after Bayesian Optimisation, in Figure 6(b).
In particular, the new infinite sum of discounted costs, when approximated by
the first 100 steps, obtains the value $204,635.7$, while the value reached
after 1000 steps is $206,358.8$. This means that Bayesian Optimisation
achieves an 80.8% decrease in the cost, as expressed through the performance
index.

(a) System output trajectories after Bayesian Optimisation

(b) Cost Comparison before and after Bayesian Optimisation

Figure 6:  (a) Output trajectories where the weighting matrices $Q$ and $R$
were obtained through Bayesian Optimisation. (b) A comparison of the
performance index before and after the optimisation of $Q$ and $R$. There is a
significance decrease of the total cost after Bayesian Optimisation, while
simultaneously .

Different implementations might yield slightly different numerical results,
depending on the randomness introduced into the algorithm. In particular,
randomness is introduced in the data generation step, as well as the sampling
for Bayesian Optimisation. However, as long as the learning algorithm is
designed according to Algorithm 2 and the data is generated appropriately, as
described in Section 6.2, the resulting controller will successfully track the
reference signal, and Bayesian Optimisation will consistently reduce the cost
while maintaining good tracking performance.

##  7 Conclusion

In this work, we investigated the problem of temperature LQT control within
the extruder of a Big Area Additive Manufacturing system, through tracking a
reference signal. In particular, we focused on the case where direct access to
the system state is not possible, instead some measurements of the state are
available. We considered a state-space model identified and validated based on
real data describing the system dynamics, and used it to simulate the system’s
behaviour and response to different controllers. The lack of access to the
system state introduces the need for a state estimation step, to better inform
the controller. Assuming the system model is available, the task of state
estimation can be performed through the implementation of a state observer.
The resulting state estimate is then used to inform a feedback controller,
designed using model-based Reinforcement Learning Techniques, with the Bellman
equation being the main tool.

Next, we investigated solving the Output Feedback problem in a model-free and
data-driven way. We used the state space temperature model for data generation
but no longer assumed it available for the control design. We instead used
Least Squares to train a controller directly from process data, while
incorporating a state estimation procedure in the algorithm. As with many
learning algorithms, the data quality is really important when training data-
driven controllers. For this reason, we provide an in depth discussion on the
input signals that would generate "good" data when applied to the system.
Finally, we explored the use of Bayesian Optimisation to improve the choice of
weighting parameters in the performance index.

We presented our results before and after the implementation of Bayesian
Optimisation, both for the model-based and data-driven controllers. For the
model-based scenario, the obtained controller successfully tracks the
reference signal, which was expected given the access to full information
about the system. We were able to achieve comparable results while alleviating
the need for knowledge of the system dynamics. This is a very encouraging
result, considering how challenging it is to have accurate models of AM
processes. Our results also indicate the improvement Bayesian Optimisation
offers, as expressed through the significant reduction of the performance
index, namely a 58.5% improvement in the model-based case and an 80.8%
improvement in the data-driven case.

While our approach is a valuable tool to solving the LQT problem, it presents
certain limitations. The basis of our approach is the quadratic form of the
objective function (performance index). When defining a control problem, the
objective function can take many forms, especially when introducing
constraints into the problem. For example, we can impose constraints on the
output response of the system so that the values follow steadily increasing
trajectories and the dip seen in Figures 6(a) and 5(a) is avoided. Another
important drawback of our approach stems from its mathematical and
computational complexity. In particular, the implementation of Least Squares
for the estimation of kernel matrix $H$ involves matrix operations between
very large matrices. Normally, the more complex the system to be controlled
is, the larger the state space needed to describe it is. This means that the
matrices involved in the introduced method can become significantly larger
than the ones in our case study.

The above concerns can be addressed through the introduction of more abstract
Reinforcement Learning algorithms, based on incorporating Deep Neural Networks
[21]. These are powerful function approximators that can handle large amounts
of data, thus being able to handle much more complicated systems and
alieviating the need for a particular form to the value function. These
approaches are widely known as Deep Reinforcement Learning methods. A very
popular one among them is the Actor-Critic framework [46], where one network
approximates the value function (critic) and another produces the new action
to be taken (actor). We leave these approaches as future work which we explore
elsewhere.

## References

  * [1] B. Berman, 3-d printing: The new industrial revolution, Business Horizons 55 (2) (2012) 155–162.  doi:https://doi.org/10.1016/j.bushor.2011.11.003.
URL https://www.sciencedirect.com/science/article/pii/S0007681311001790

  * [2] C. Buchanan, L. Gardner, Metal 3d printing in construction: A review of methods, research, applications, opportunities and challenges, Engineering Structures 180 (2019) 332–348.  doi:https://doi.org/10.1016/j.engstruct.2018.11.045.
URL https://www.sciencedirect.com/science/article/pii/S0141029618307958

  * [3] U. M. Dilberoglu, B. Gharehpapagh, U. Yaman, M. Dolen, The Role of Additive Manufacturing in the Era of Industry 4.0, Procedia Manufacturing 11 (2017) 545–554, 27th International Conference on Flexible Automation and Intelligent Manufacturing, FAIM2017, 27-30 June 2017, Modena, Italy.  doi:https://doi.org/10.1016/j.promfg.2017.07.148.
URL https://www.sciencedirect.com/science/article/pii/S2351978917303529

  * [4] W. E. Frazier, Metal additive manufacturing: a review, Journal of Materials Engineering and performance 23 (2014) 1917–1928.
  * [5] A. Uriondo, M. Esperon-Miguez, S. Perinpanayagam, The present and future of additive manufacturing in the aerospace sector: A review of important aspects, Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering 229 (11) (2015) 2132–2147.
  * [6] D. D. Gu, W. Meiners, K. Wissenbach, R. Poprawe, Laser additive manufacturing of metallic components: materials, processes and mechanisms, International materials reviews 57 (3) (2012) 133–164.
  * [7] F. J. Mercado Rivera, A. J. Rojas Arciniegas, Additive manufacturing methods: techniques, materials, and closed-loop control applications, The International Journal of Advanced Manufacturing Technology 109 (2020) 17–31.
  * [8] H. Lhachemi, A. Malik, R. Shorten, Augmented reality, cyber-physical systems, and feedback control for additive manufacturing: A review, IEEE Access 7 (2019) 50119–50135.
  * [9] B. D. Anderson, J. B. Moore, Optimal control: linear quadratic methods, Courier Corporation, 2007.
  * [10] K. Ogata, et al., Modern control engineering, Vol. 5, Prentice hall Upper Saddle River, NJ, 2010.
  * [11] F. Previdi, S. M. Savaresi, A. Panarotto, Design of a feedback control system for real-time control of flow in a single-screw extruder, Control engineering practice 14 (9) (2006) 1111–1121.
  * [12] M. H. Farshidianfar, A. Khajepour, A. Gerlich, Real-time control of microstructure in laser additive manufacturing, The International Journal of Advanced Manufacturing Technology 82 (5) (2016) 1173–1186.  doi:10.1007/s00170-015-7423-5.
URL https://doi.org/10.1007/s00170-015-7423-5

  * [13] D. Hu, R. Kovacevic, Sensing, modeling and control for laser-based additive manufacturing, International Journal of Machine Tools and Manufacture 43 (1) (2003) 51–60.  doi:https://doi.org/10.1016/S0890-6955(02)00163-3.
URL https://www.sciencedirect.com/science/article/pii/S0890695502001633

  * [14] J.-P. Kruth, P. Mercelis, J. Van Vaerenbergh, T. Craeghs, Feedback control of selective laser melting, in: Virtual and Rapid Manufacturing, Crc Press, 2007, pp. 521–528.
  * [15] P. M. Sammons, M. L. Gegel, D. A. Bristow, R. G. Landers, Repetitive Process Control of Additive Manufacturing With Application to Laser Metal Deposition, IEEE Transactions on Control Systems Technology 27 (2) (2019) 566–575.  doi:10.1109/TCST.2017.2781653.
  * [16] H. Bikas, P. Stavropoulos, G. Chryssolouris, Additive manufacturing methods and modelling approaches: a critical review, The International Journal of Advanced Manufacturing Technology 83 (2016) 389–405.
  * [17] J. Si, A. G. Barto, W. B. Powell, D. Wunsch, Handbook of learning and approximate dynamic programming, Vol. 2, John Wiley & Sons, 2004.
  * [18] R. S. Sutton, A. G. Barto, Reinforcement learning: An introduction, MIT press, 2018.
  * [19] R. Bellman, Dynamic programming, Science 153 (3731) (1966) 34–37.
  * [20] D. Bertsekas, Dynamic programming and optimal control: Volume I, Vol. 1, Athena scientific, 2012.
  * [21] D. P. Bertsekas, J. N. Tsitsiklis, Neuro-dynamic programming, Athena Scientific, 1996.
  * [22] B. Yao, F. Imani, H. Yang, Markov Decision Process for Image-Guided Additive Manufacturing, IEEE Robotics and Automation Letters 3 (4) (2018) 2792–2798.  doi:10.1109/LRA.2018.2839973.
  * [23] K. Wasmer, T. Le-Quang, B. Meylan, S. A. Shevchik, In situ quality monitoring in AM using acoustic emission: A reinforcement learning approach, Journal of Materials Engineering and Performance 28 (2019) 666–672.
  * [24] M. Alicastro, D. Ferone, P. Festa, S. Fugaro, T. Pastore, A reinforcement learning iterated local search for makespan minimization in additive manufacturing machine scheduling problems, Computers & Operations Research 131 (2021) 105272.
  * [25] S. Patrick, A. Nycz, M. Noakes, Reinforcement learning for generating toolpaths in additive manufacturing, in: 2018 International Solid Freeform Fabrication Symposium, University of Texas at Austin, 2018.
  * [26] A. G. Dharmawan, Y. Xiong, S. Foong, G. S. Soh, A model-based reinforcement learning and correction framework for process control of robotic wire arc additive manufacturing, in: 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2020, pp. 4030–4036.
  * [27] F. Ogoke, A. B. Farimani, Thermal control of laser powder bed fusion using deep reinforcement learning, Additive Manufacturing 46 (2021) 102033.
  * [28] E. Zavrakli, A. Parnell, A. Dickson, S. Dey, Data-driven Linear Quadratic Tracking based Temperature Control of a Big Area Additive Manufacturing System, arXiv preprint arXiv:2307.07039 (2023).
  * [29] B. Kiumarsi, F. L. Lewis, H. Modares, A. Karimpour, M.-B. Naghibi-Sistani, Reinforcement Q-learning for optimal tracking control of linear discrete-time systems with unknown dynamics, Automatica 50 (4) (2014) 1167–1175.
  * [30] D. Luenberger, An introduction to observers, IEEE Transactions on automatic control 16 (6) (1971) 596–602.
  * [31] C. Chen, W. Sun, G. Zhao, Y. Peng, Reinforcement Q-Learning Incorporated With Internal Model Method for Output Feedback Tracking Control of Unknown Linear Systems, IEEE Access 8 (2020) 134456–134467.  doi:10.1109/ACCESS.2020.3011194.
  * [32] P. I. Frazier, A tutorial on Bayesian optimization, arXiv preprint arXiv:1807.02811 (2018).
  * [33] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, N. De Freitas, Taking the human out of the loop: A review of Bayesian optimization, Proceedings of the IEEE 104 (1) (2015) 148–175.
  * [34] R. Calandra, A. Seyfarth, J. Peters, M. P. Deisenroth, Bayesian optimization for learning gaits under uncertainty: An experimental comparison on a dynamic bipedal walker, Annals of Mathematics and Artificial Intelligence 76 (2016) 5–23.
  * [35] A. Marco, P. Hennig, J. Bohg, S. Schaal, S. Trimpe, Automatic LQR tuning based on Gaussian process global optimization, in: 2016 IEEE international conference on robotics and automation (ICRA), IEEE, 2016, pp. 270–277.
  * [36] K. Miyamoto, J. She, D. Sato, N. Yasuo, Automatic determination of LQR weighting matrices for active structural control, Engineering Structures 174 (2018) 308–321.
  * [37] D. Gootjes, Applying feedback control to improve 3d printing quality, Master’s thesis, Delf University of Technology (2017).
  * [38] E. H. Moore, On the reciprocal of the general algebraic matrix, Bull. Am. Math. Soc. 26 (1920) 394–395.
  * [39] R. Penrose, A generalized inverse for matrices, Mathematical Proceedings of the Cambridge Philosophical Society 51 (3) (1955) 406–413.  doi:10.1017/S0305004100030401.
  * [40] V. G. Lopez, M. Alsalti, M. A. Müller, Efficient off-policy Q-learning for data-based discrete-time lqr problems, IEEE Transactions on Automatic Control (2023).
  * [41] A. Krishnamoorthy, D. Menon, Matrix inversion using Cholesky decomposition, in: 2013 Signal Processing: Algorithms, Architectures, Arrangements, and Applications (SPA), 2013, pp. 70–72.
  * [42] J. Snoek, H. Larochelle, R. P. Adams, Practical bayesian optimization of machine learning algorithms, Advances in neural information processing systems 25 (2012).
  * [43] C. E. Rasmussen, C. K. Williams, et al., Gaussian processes for machine learning, Vol. 1, Springer, 2006.
  * [44] L. Ljung, System identification, in: Signal analysis and prediction, Springer, 1998, pp. 163–173.
  * [45] S. M. Shinners, Modern control system theory and design, John Wiley & Sons, 1998\.
  * [46] V. Konda, J. Tsitsiklis, Actor-critic algorithms, Advances in neural information processing systems 12 (1999).