Reinforcement Learning for Batch Bioprocess Optimization

TL;DR

This paper introduces a reinforcement learning approach using policy gradients and recurrent neural networks to optimize complex batch bioprocesses, addressing challenges like stochasticity and plant-model mismatch.

Contribution

It presents a novel RL-based optimization strategy that updates control policies from preliminary models using real plant data, outperforming traditional methods in complex scenarios.

Findings

Effective policy updates from preliminary models to real plant data

Successful application on three diverse bioprocess case studies

Advantages over nonlinear model predictive control demonstrated

Abstract

Bioprocesses have received a lot of attention to produce clean and sustainable alternatives to fossil-based materials. However, they are generally difficult to optimize due to their unsteady-state operation modes and stochastic behaviours. Furthermore, biological systems are highly complex, therefore plant-model mismatch is often present. To address the aforementioned challenges we propose a Reinforcement learning based optimization strategy for batch processes. In this work, we applied the Policy Gradient method from batch-to-batch to update a control policy parametrized by a recurrent neural network. We assume that a preliminary process model is available, which is exploited to obtain a preliminary optimal control policy. Subsequently, this policy is updatedbased on measurements from thetrueplant. The capabilities of our proposed approach were tested on three case studies (one of which is nonsmooth) using a more complex process model for thetruesystemembedded with adequate process disturbance. Lastly, we discussed the advantages and disadvantages of this strategy compared against current existing approaches such as nonlinear model predictive control.

Figures

Tables

Table 1. Table 1 : Parameter values for physical system ( 22 - 24 )

Parameter Values
$u_m$	0.0572	$h^{-1}$
$u_d$	0.0	$h^{-1}$
$K_N$	393.1	$mg/L$
$Y_{NX}$	504.1	$mg/g$
$k_m$	0.00016	$mg/g/h$
$k_d$	0.281	$h^{-1}$
$k_s$	178.9	$\mu mol/m^2/s$
$k_i$	447.1	$\mu mol/m^2/s$
$k_{sq}$	23.51	$\mu mol/m^2/s$
$k_{iq}$	800	$\mu mol/m^2/s$
$K_{NP}$	16.89	$mg/L$

Equations

$${\text{x}}_{t+1}\sim p({\text{x}}_{t+1}|{\text{x}}_t,{\text{u}}_t)$$

$$\begin{array}{rl}& {\text{x}}_{t+1}=f({\text{x}}_t,{\text{u}}_t,{\text{d}}_t)\end{array}$$

$$\mathcal{P}(\pi (\cdot )):=\{\begin{array}{rl}& \underset{\pi (\cdot )}{\max }\mathbb{E}[J({\text{x}}_t^k,{\text{u}}_t^k)]\\ & \text{s.t.}\\ & {\text{x}}_0^k={\text{x}}^k(0)\\ & {\text{x}}_{t+1}^k=f({\text{x}}_t^k,{\text{u}}_t^k,{\text{d}}_t^k)\forall t\in \left\{1,...,T-1\right\}\\ & {\text{u}}_t\sim \pi ({\text{x}}_t^k)\\ & \text{u}\in \mathbb{U}\\ & \text{given}\\ & {\text{x}}_t^j\forall j\in \left\{0,...,k-1\right\}\forall t\in \left\{1,...,T\right\}\end{array}$$

$$J(\mathbf{\tau })=\sum _{t=0}^T{\gamma }^t{R}_t({\text{u}}_t,{\text{x}}_t)$$

$$p(\mathbf{\tau }|\theta )=\hat{\mu }({\text{x}}_0)\prod _{t=0}^{T-1}\left[\pi ({\text{u}}_t|{\text{x}}_t,\theta )p({\text{x}}_{t+1}|{\text{x}}_t,{\text{u}}_t)\right]$$

$$\underset{\pi (\cdot )}{\max }{\mathbb{E}}_{\mathbf{\tau }\sim p(\mathbf{\tau }|\theta )}[J(\mathbf{\tau })]$$

$$\text{u}\sim {\pi }_{\theta }(\text{u}|\text{x})=\pi (\text{u}|\text{x},\theta )=p({\text{u}}_t=\text{u}|{\text{x}}_t=\text{x},{\theta }_t=\theta ).$$

$${\theta }_{m+1}={\theta }_m+{\alpha }_m{\nabla }_{\theta }{\mathbb{E}}_{\mathbf{\tau }\sim p(\mathbf{\tau }|\theta )}[J(\mathbf{\tau })]$$

$$\hat{J}(\theta )={\nabla }_{\theta }{\mathbb{E}}_{\mathbf{\tau }\sim p(\mathbf{\tau }|\theta )}[J(\mathbf{\tau })]$$

$$\text{leavevmode}=\text{leavevmode}\int {\nabla }_{\theta }p(\mathbf{\tau }|\theta )\text{leavevmode}J(\mathbf{\tau })\text{d}\mathbf{\tau }$$

$$\text{leavevmode}=\text{leavevmode}\int p(\mathbf{\tau }|\theta )\text{leavevmode}\frac{{\nabla }_{\theta }p(\mathbf{\tau }|\theta )}{p(\mathbf{\tau }|\theta )\text{leavevmode}}\text{leavevmode}J(\mathbf{\tau })\text{d}\mathbf{\tau }$$

$$\text{leavevmode}=\text{leavevmode}\int p(\mathbf{\tau }|\theta )\text{leavevmode}{\nabla }_{\theta }\text{log}\left(p(\mathbf{\tau }|\theta )\right)J(\mathbf{\tau })\text{d}\mathbf{\tau }$$

$$\text{leavevmode}=\text{leavevmode}{\mathbb{E}}_{\mathbf{\tau }}\left[J(\mathbf{\tau }){\nabla }_{\theta }\text{log}\left(p(\mathbf{\tau }|\theta )\right)\right]$$

$${\theta }_{m+1}={\theta }_m+{\alpha }_m{\mathbb{E}}_{\mathbf{\tau }}\left[J(\mathbf{\tau }){\nabla }_{\theta }\text{log}\left(p(\mathbf{\tau }|\theta )\right)\right]$$

$${\nabla }_{\theta }\text{log}\left(p(\mathbf{\tau }|\theta )\right)={\nabla }_{\theta }\sum _{t=0}^{T-1}\text{log}\left(\pi ({\text{u}}_t|{\text{x}}_t,\theta )\right)$$

$${\nabla }_{\theta }{\mathbb{E}}_{\mathbf{\tau }}[J(\mathbf{\tau })]={\mathbb{E}}_{\mathbf{\tau }}\left[J(\mathbf{\tau }){\nabla }_{\theta }\sum _{t=0}^{T-1}\text{log}\left(\pi ({\text{u}}_t|{\text{x}}_t,\theta )\right)\right]$$

$$\begin{array}{rl}{\mathbf{\nu }}_t& =h_{\theta }({\hat{\text{x}}}_t,{\text{u}}_{t-1})\\ {\text{u}}_t& =G_t({\mathbf{\nu }}_t)\end{array}$$

$${\nabla }_{\theta }{\mathbb{E}}_{\mathbf{\tau }}\approx \frac{1}{K}\sum _{k=1}^K\left[J({\mathbf{\tau }}^k){\nabla }_{\theta }\sum _{t=0}^{T-1}\text{log}\left(\pi ({\text{u}}_t^k|{\hat{\text{x}}}_t^k,\theta )\right)\right]$$

$$b=\hat{J}(\theta )\approx \frac{1}{K}\sum _{k=1}^{K}J({\mathbf{\tau }}^{(k)}),$$

$${\nabla }_{\theta }\hat{J}(\theta )\approx \frac{1}{K}\sum _{k=1}^K\left[(J({\mathbf{\tau }}^k)-b){\nabla }_{\theta }\sum _{t=0}^{T-1}\text{log}\left(\pi ({\text{u}}_t^k|{\hat{\text{x}}}_t^k,\theta )\right)\right]$$

$$\frac{d{y}_1}{dt}=-(u_1+0.5\text{leavevmode}{u}_1^2)y_1+0.5\frac{u_2{y}_2}{(y_1+y_2)}$$

$$\frac{d{y}_2}{dt}=u_1\text{leavevmode}{y}_1-0.7u_2\text{leavevmode}{y}_1$$

$$\frac{d{y}_1}{dt}=-(u_1+0.5\text{leavevmode}{u}_1^2)y_1+u_2$$

$$\frac{d{y}_2}{dt}=u_1\text{leavevmode}{y}_1-u_2\text{leavevmode}{y}_1$$

$$\begin{array}{rl}{R}_t& =0,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}t\in \{0,T-1\}\\ R_T& =y_2(T).\end{array}$$

$$d{y}_1=\left[-(u_1+0.5\text{leavevmode}{u}_1^2)y_1+0.5\frac{u_2{y}_2}{(y_1+y_2)}\right]dt$$

$$d{y}_2=\left[u_1\text{leavevmode}{y}_1-0.7u_2\text{leavevmode}{y}_1\right]dt+\left[0.1\sqrt{y_1}\right]dW$$

$$\frac{d{c}_x}{dt}$$

$$\frac{d{c}_N}{dt}$$

$$\frac{d{c}_q}{dt}$$

$$w(t)$$

$${\sigma }_d$$

$${\sigma }_n$$

$$noise(t)\sim \mathcal{N}(\text{0},diag(4\times 10^{-4},.1,1\times 10^{-8})).$$

$$\left[\begin{array}{ccc}{c}_x(0)& c_N(0)& c_q(0)\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{ccc}1.& 150.& 0.\end{array}\right],diag(1\times 10^{-3},22.5,0.)).$$

$$\begin{array}{rl}{R}_t& =-\Delta {\text{u}}_t^{T}diag(3.125\times 10^{-8},3.125\times 10^{-6})\Delta {\text{u}}_t^T,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}t\in \{0,T-1\}\\ R_T& =c_q(T),\end{array}$$

Full text

In this work we applied the Policy Gradient method from batch-to-batch to update a control policy parametrized by a recurrent neural network. We assume that a preliminary process model is available, which is exploited to obtain a preliminary optimal control policy. Subsequently, this policy is updated based on measurements from the true plant. The capabilities of our proposed approach were tested on three case studies (one of which is nonsmooth) using a more complex process model for the true system embedded with adequate process disturbance. Lastly, we discussed advantages and disadvantages of this strategy compared against current existing approaches such as nonlinear model predictive control.

1 Introduction

There has been a global interest in using sustainable bio-production systems to produce a broad range of chemicals to substitute fossil derived synthetic routes [11, 21]. Bioprocesses exploit microorganisms to synthesize platform chemicals and high-value products by utilizing different means of resources [22]. Compared to a traditional chemical process, a biochemical process is highly complex due to the intricate relationships between metabolic reaction networks and culture fluid dynamics [16]. As a result, it is difficult to construct accurate physics-based models to simulate general large-scale biosystems, and plant-model mismatch is inevitable. Furthermore, the behaviour of bioprocesses is often stochastic in the macro-scale, because biological metabolic pathways are sensitive even to mild changes of process operating conditions [58, 51]. Therefore, developing control and optimization strategies for bioprocesses remains an open challenge. Given these critical limitations in mechanistic models, in this work we propose a data-driven approach to address this issue.

We must seek a strategy that can optimize our process and handle both the system’s stochastic behaviours (e.g. process disturbances) and plant-model mismatches. It is here that we have opted to use Reinforcement learning and more specifically, Policy Gradients, as an alternative to current existing methods. Reinforcement learning (RL) was developed to address nonlinear and stochastic optimal control problems [6]. Two main branches have been established on how to solve dynamic optimization problems via RL.

The first branch is based on Dynamic Programming (DP), hence termed Approximate Dynamic Programming (ADP). DP relies on the Hamilton-Jacobi-Bellman equation (HJBE), the solution of which becomes intractable for small size problems with nonlinear dynamics and continuous state and control actions. Because of this, past research has relied on using ADP techniques to find approximate solutions to these problems  [48].

The second branch is to use Policy Gradients, which directly obtains a policy by maximizing a desired performance index. This approach is particularly well suited to deal with problems where both the state space and the control space are continuous. Given the advantages that Policy Gradients can offer when confronted with bioprocess optimization, we have adopted this approach in the current work. Policy Gradient methods, along with their benefits, are further explained in Section  2.2.

1.1 Related Work

Given that chemical engineers have always dealt with complex and uncertain systems there have been several approaches that address specific instances of the aforementioned problems, we highlight some related previous work in the sections below.

One example to track stochastic batch-to-batch systems is Iterative Learning Control (ILC) which was initially introduced for robot manipulators [2], and later implemented by the process control community [56]. ILC deals with the problem of tracking the control performance in batch processes given a reference trajectory for runs that last a fixed time, and where the process state is reset to the same value at the start of each run. An overview of ILC strategies in process control can be found in [28].

Real-time optimization (RTO) is another method that deals with uncertain processes. The main idea is to represent the process dynamics by a nonlinear input/output mapping where the disturbances are explicitly accounted for. This mapping is then used to optimize some desired performance index. For the interested reader, further details can be found in [8] and [12]. A recent review on this topic can be found in [32]. For the dynamic systems, these methodologies are usually referred to as Dynamic Real-time optimization (DRTO) which is closely related to NMPC. More details can be found in [41, 42].

Another technique that deals with stochastic systems is model predictive control (MPC), and its extension to nonlinear systems, NMPC. NMPC has a vast variety of methods that can incorporate uncertainty or maintain properties under the presence of stochastic environments. The most common paradigms are the stochastic NMPC [33] and the Robust NMPC [4], where the former incorporates the uncertainty by minimizing (usually) the expectation of the objective function, whilst the latter approach solves a min-max optimization by minimizing the worst case scenario of the uncertainty. Both of these approaches require knowledge regarding the nature of the uncertainty in order to proceed.

There are different approaches that have been proposed for NMPC frameworks, including scenario [5] based multi-stage schemes for nonlinear systems [31, 25], where stochastic programming is utilized and future information is incorporated in an adaptive manner. Another approach is the use of Gaussian processes [9, 10] or using (generalized) polynomial chaos expansions [23] to model effectively the uncertainties of the process. In the case where no proper information for the uncertainty is available, e.g. there is not enough data to conduct uncertainty quantification, optimal control is explored using the nominal linear or nonlinear available model. In terms of solution procedures for the dynamic optimization problem, it is common to use a direct approach after parametrizing and discretizing the control inputs [52] or the system dynamics [7] resulting in a nonlinear programming problem. Although much less common, indirect approaches can also be used, where the necessary conditions of optimality are solved explicitly [3]. If no information on structural information is known, conservative assumptions can be made in order to establish stability conditions [19, 39, 38].

Reinforcement Learning (in an Approximate Dynamic Programming philosophy), has lately caught significant attention for chemical process control. For example, in [30] a model-based strategy and a model-free strategy for control of nonlinear processes were proposed, in [37] ADP strategies were used to address fed-batch reactor optimization, in [27] mixed-integer decision problems were addressed with applications to scheduling. In [50] with the inclusion of distributed optimization techniques, an input-constrained optimal control problem solution technique was presented, among other works (e.g. [13], [45]). All these approaches rely on the (approximate) solution of the HJBE, and have been shown to be reliable and robust for several problem instances.

In this paper, we present another take on RL, that of using Policy Gradients. Policy Gradient methods directly estimate the control policy, without the need of a model, or the solution of the HJBE, its advantages are highlighted in the following section.

In addition to the above, for recent reviews of Machine Learning and Artificial Intelligence applied to chemical engineering the reader is referred to [29] and [53]. A shorter review focused on modelling bioprocesses with ML tools can be found in [17].

1.2 Motivation

The process systems engineering community has been dealing with stochastic batch-to-batch systems for a long time. For example, nonlinear dynamic optimization and particularly NMPC are a powerful methodology to address uncertain dynamic systems, however there are several properties that make its application less attractive. All the approaches in NMPC require the knowledge of a detailed model that describe the system dynamics, and stochastic NMPC additionally requires an assumption for the uncertainty quantification/propagation. Furthermore, the online computational time may be a bottleneck for real time applications since a (possibly) nonlinear optimization problem has to be solved.

In contrast, RL directly accounts for the effect of future uncertainty and its feedback in a proper ‘closed-loop’ manner, whereas conventional NMPC assumes open-loop control actions at future time points in the prediction, which can lead to overly conservative control actions [30]. In addition, policy gradients can establish a policy in a model-free fashion and excel at on-line computational time. This is because the online computations require only evaluation of a policy, since all the computational cost is shifted off-line.

As mentioned previously, Real-time optimization (RTO) has been used to address many instances of batch-to-batch problems. Interestingly, some recent approaches have suggested a hybrid modeling strategy, where function approximates are used in conjunction with a pre-existing model [14, 20]. From some perspectives these recent algorithms could be thought of as model-based Reinforcement learning approaches. However, there is not yet a clear consensus on how to address problems with plant-model mismatch, measurement noise, and disturbances in an RTO framework.

In terms of previous RL approaches in chemical engineering to address process control and optimization, they have relied on action-value methods (e.g. Q-learning, solution of the HJBE). However, to address continuous nonlinear action domains Policy Gradient methods present several advantages:

In Policy Gradient methods, the approximate policy can naturally approach a deterministic policy, whereas action-value methods (that use epsilon-greedy or Boltzmann functions) select a random control action with some heuristic rule [48].

• 2.

Although it is possible to estimate the objective value of state-action pairs in continuous action spaces by function approximators, this does not help choose a control action. Therefore, on-line optimization over the action space for each time-step should be performed, which can be slow and inefficient. Policy Gradient methods work directly with policies that emit probability distributions, which is much faster and does not require an online optimization step.

• 3.

Policy Gradient methods are guaranteed to converge at least to a locally optimal policy even in high dimensional continuous state and action spaces, unlike action-value methods where convergence to local optima is not guaranteed [48].

• 4.

Policy Gradient methods enable the selection of control actions with arbitrary probabilities. In such cases, the best approximate policy may be stochastic [48].

Due to the above advantages, in this work we propose an optimization strategy that uses a Policy Gradient algorithm to optimize batch-to-batch bioprocesses. This work extends our proposed methodology in [40], the new approach presents a much faster adaptation time by implementing transfer learning for the efficient adaptation of the policies. Additional more complex case studies and a comparison against NMPC are included. Furthermore, we exemplify both approaches (NMPC and our approach) in a system described by a nonsmooth differential equation model. The difficulty for nonsmooth models is highlighted in [46].

2 Methodology

2.1 Problem Statement

In this work, we assume that the system’s dynamics are given by an (generally) unknown probability distribution, following a Markov process:

[TABLE]

This system can be approximated by the following discrete time stochastic nonlinear system represented as a state-space model:

[TABLE]

where $t$ represents the discrete time, $\textbf{x}_{t}\in\mathbb{R}^{n_{x}}$ is the vector of states, $\textbf{u}_{t}\in\mathbb{R}^{n_{u}}$ is the vector of inputs, $\textbf{d}_{t}\in\mathbb{R}^{n_{d}}$ is the the vector of process disturbances, and $f(\cdot)$ are the nonlinear dynamics of the physical system.

Our strategy seeks to find the optimal policy for a batch process under the presence of disturbances and measurement noise. Then, the problem can be written as an Optimal Control Problem (OCP):

[TABLE]

The objective is to maximize the expectation of an economic criterion $J$, where $k$ is the current batch, while $j$ refers to previous batch realizations. Additionally, the optimization problem (3) searches for a set of functions $\pi(\cdot)$ that maps the probability of $u_{t}^{k}$ given $x_{t}^{k}$. Notice that in problem (3) we make no assumptions about the nature of d. Even in the case where the dynamics are fully known, the solution of problem (3) may be intractable for medium size systems.

To overcome this limitation a novel strategy is proposed, where a policy $\pi_{\theta}(\cdot)$, parametrized by the parameters $\theta$, is constructed that maximizes the expectation of a performance index $J$. The states at the time $t+1$ are assumed to be given by the probability density $p(\textbf{x}_{t+1}|\textbf{x}_{t},\textbf{u}_{t})$. The interaction with the policy can be depicted as a closed-loop, see Fig. 1.

Let $\boldsymbol{\tau}$ denote a joint random variable of states and controls defining a trajectory with a time horizon $T$: $\boldsymbol{\tau}=(\textbf{x}_{0},\textbf{u}_{0},R_{0},...,\textbf{x}_{T-1},\textbf{u}_{T-1},R_{T-1},\textbf{x}_{T},R_{T})$, the performance index being

[TABLE]

where $\gamma\in(0,1]$ is the discount factor and $R_{t}$ a given reward at the time instance $t$ for the values of $u_{t}$, $x_{t}$. We represent the probability density of a trajectory as:

[TABLE]

where $\hat{\mu}(\textbf{x}_{0})$ is the probability density of the initial state. We can therefore state the following optimization problem:

[TABLE]

notice that the process dynamics are implicit in $\boldsymbol{\tau}$. To solve problem (6) we turn our attention to policy gradient methods.

2.2 Policy Gradient Methods

Policy gradient methods compute a policy that maximizes the expectation over some objective function ($i.e.$ problem (6) ). They rely on a parametrized policy function $\pi_{\theta}(\cdot)$ that returns an action u given a state of the system x and a set of intrinsic parameters $\theta$ of the policy. In the case of stochastic policies, the policy function returns the defining parameters of a probability distribution over all possible actions, from which the actions are sampled:

[TABLE]

In this work, a Recurrent neural network (RNN) is used as the parametrized policy, which takes (a number of past) states and control actions as inputs and returns the moments of a probability distribution. Then the next control action is drawn from the corresponding probability distribution. For example, if the control actions live in a normal distribution then a mean and a variance are computed, from these mean and variance a control action can be drawn. In this setting, the exploitation-exploration trade-off is represented explicitly by the value of the variance of the underlying distribution of the policy. Deterministic policies can be seen as a limiting case where the variance converges to zero.

To learn the optimal policy, we seek to maximize our performance metric, and hence we can follow a gradient ascent strategy:

[TABLE]

where $m$ is the current iteration that the parameters are updated (epoch), $\mathbb{E}_{\boldsymbol{\tau}\sim p(\boldsymbol{\tau}|\theta)}[J(\boldsymbol{\tau})]$ is the expectation of $J$ over $\boldsymbol{\tau}$ and $\alpha_{m}$ is the step size (also termed learning rate in the RL community) for the $m^{th}$ iteration. Computing $\hat{J}(\theta)=\mathbb{E}_{\boldsymbol{\tau}\sim p(\boldsymbol{\tau}|\theta)}[J(\boldsymbol{\tau})]$ directly is difficult, therefore we use the Policy Gradient Theorem [49], which shows the following:

[TABLE]

Notice from (9b) that, $p(\boldsymbol{\tau}|\theta)\leavevmode\nobreak\ J(\boldsymbol{\tau})$ is an objective function value multiplied by its probability density, therefore, integrating this over all possible values of $\boldsymbol{\tau}$ we obtain the expected value. From there we arrive at (9e), where, dropping the explicit distribution of $\boldsymbol{\tau}$, gives us an unbiased gradient estimator, (8) now becomes:

[TABLE]

Using the expression for $p(\boldsymbol{\tau}|\theta)$ in (5) and taking its logarithm, we obtain:

[TABLE]

Note that since $p(\textbf{x}_{t+1}|\textbf{x}_{t},\textbf{u}_{t})$ and $\hat{\mu}(\textbf{x}_{0})$ are independent of $\theta$ they disappear from the above expression. Then we can rewrite (9e) for a trajectory as:

[TABLE]

Notice that expression (12) does not require the knowledge of of the dynamics of the physical system. However, the above update presents two challenges: the selection of the policy $\pi(\textbf{u}_{t}|\textbf{x}_{t},\theta)$ and the computation of the expectation. To address these possible issues, in this work, recurrent neural networks are used to parametrize the policy of the policy gradient (presented in Section 2.3), while a Monte-Carlo method is utilized to approximate the expectation (presented in Section 2.4).

2.3 Recurrent Neural Network

Recurrent neural networks, (RNNs) [43], are a type of artificial neural network that has been tailored to address sequential data. RNNs produce an output at each time step and have recursive connections between hidden units. This allows them to have a ‘memory’ of previous data and hence be well suited to model time-series. In general, RNNs can be depicted as a folded computational graph as presented in Fig. 2.

Fig. 2 shows how an input $\hat{\textbf{x}}_{t}$ is presented to the network as well as the recursive state of RNN, $\textbf{u}_{t-1}$ and outputs $\boldsymbol{\nu}_{t}$. A more detailed representation of an RNN is depicted in Fig. 3, which is equivalent to a series of unfolded nodes associated with a particular time instance.

In Fig. 3 we can appreciate that each node receives two inputs $\textbf{x}_{t-i}$ and $\textbf{u}_{t-i-1}$. Generally speaking, the input $\textbf{x}_{t-i}$ corresponds to the data supplied to that node, such as in a traditional artificial neural network (also referred to as feedforward). The unfolded computational graph in this case can be represented as a dynamic system:

[TABLE]

where $\hat{\textbf{x}}_{t}$ is the vector that contains all the external variables for the RNN, $G$ the function that computes the output of the network $u$, and $h_{\theta}$ represents the layers of the neural networks. Deep structures (which means having more than one hidden layer) can be employed to enhance the performance of the network (deep neural networks) which have been combined with Reinforcement learning recently in [35, 36]. Previous realizations of the states $x$ and the controls $u$ are also used as input variables to the network, e.g. $\hat{\textbf{x}}_{t}=\left[\textbf{x}_{t}^{T},\dots,\textbf{x}_{t-N}^{T},\textbf{u}_{t-2}^{T},\dots,\textbf{u}_{t-N-1}^{T}\right]^{T}$, to model dynamic systems. RNNs have previously been applied either as a surrogate model of the process dynamics [47] or as a parametrization of the agent (control policy) [34].

In this work, RNNs are applied to parameterize the stochastic policy. We must remark that in theory the Markov decision process does not require RNNs (due to the Markov property), however in practice the use of RNNs can improve the performance of the policy by exploiting additional memory that is provided. In the current work the RNN initially computes the mean and the variance of a multivariate normal distribution where the control actions live. Subsequently, the actual control action is drawn. Precisely, $\boldsymbol{\nu}_{t}=\left[\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t}\right]$ where the $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are the mean and variance, respectively, and $\textbf{u}_{t}=G_{t}(\boldsymbol{\nu}_{t})$ is substituted by $\textbf{u}_{t}\sim\mathcal{N}(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$ making it a stochastic policy. Under the presence of uncertainty (stochastic in nature) a deterministic policy will fail as the control action will always be the same for the same states since it learns a deterministic mapping from states to control actions at the exact same state. On the contrary, a stochastic policy draws a control action from a probability distribution which can account for stochastic environments.

In Figure 4, a schematic representation of a policy network where the stochastic policy follows a Gaussian distribution is depicted. In this figure, we can observe how the states $\hat{x}_{t}$ are used as input to the network, and how the network computes the mean $\mu_{t+1}$ and standard deviation $\sigma_{t+1}$ for the subsequent time step. Then, the control action $\textbf{u}_{t+1}$ is drawn from the distribution defined by the outputs of the network.

2.4 Reinforce Algorithm

Given that the parametrized policy used in this work is a RNN, it must be trained to adjust its weights so that the output corresponds to an optimal control action. To this end, we use the steepest ascent strategy mentioned in (8). However, computing the expectation in (12) can be an intractable problem, and this expression is needed to compute the steepest ascent update ( 10). Therefore, we propose to use the Reinforce algorithm to compute the policy gradient. The Reinforce algorithm  [55] approximates the gradient of the policy to maximize the expected reward with respect to the parameters $\theta$ without the need of a dynamic model of the process. To compute the expectation we take several sample trajectories and then approximately calculate $\mathbb{E}_{\boldsymbol{\tau}}\left[J(\boldsymbol{\tau})\nabla_{\theta}\sum_{t=0}^{T-1}\text{log}\left(\pi(\textbf{u}_{t}|\textbf{x}_{t},\theta)\right)\right]$ as an average of $K$ samples:

[TABLE]

where we denote the sample $k$ as a super-index. The variance of this estimation can be reduced with the aid of an action-independent baseline $b$, which does not introduce a bias [48]. A simple but effective baseline is the expectation of reward under the current policy, approximated by the mean over the sampled paths:

[TABLE]

which leads to:

[TABLE]

This selection increases the log likelihood of an action by comparing it to the expected reward of the current policy. (16) is the gradient that we can now incorporate into our steepest ascent strategy. The algorithm that trains the RNN and obtains the optimal policy network is the following.

The steps in the Algorithm 1 are explained below.

Initialization: The RNN policy network and its weights $\theta$ are initialized, along with the algorithm’s hyperparameters such as learning rate, number of episodes and number of epochs.

Training loop: The weights on the RNN are updated by a policy gradient scheme for a total of $N$ epochs. In Step 1 $K$ trajectories are computed, each trajectory consists of $T$ time steps, and states and control actions are collected. In Step 2 the weights of the RNN are updated based on the policy gradient framework. In Step 3, either the algorithm terminates or returns to Step 1.

2.5 Reinforcement Learning for Bioprocess Optimization under Uncertainty

The methodology presented aims to overcome plant-model mismatch in uncertain dynamic systems, a usual scenario in bioprocesses. It is common to construct simple deterministic models according to a hypothesized mechanism, however the real system is more complex and presents disturbances. We propose the following methodology to address this problem (following from Algorithm 2).

Step 0, Initialization: The algorithm is initialized by considering an initial policy network (e.g. RNN policy network) with untrained parameters $\theta_{0}$.

Step 1, Preliminary Learning (Off-line): It is assumed that a preliminary mechanistic model can be constructed from previous existing process data, hence, the policy learns this preliminary mechanistic model. This is done by running Algorithm 1 in a simulated environment by the mechanistic model. This allows the policy to incorporate previously believed knowledge about the system. The policy will therefore end with an optimal control policy for the mechanistic model. The termination criteria can be defined either by the designer or by the difference from the solution of the OCP, since the process model is known.

Given that the experiments are in silico, a large number of episodes and trajectories can be generated that corresponds to different actions from the probability distribution of $\textbf{u}_{t}$, and a specific set of parameters of the RNN, respectively. The resulting control policy is a good approximation of the optimal policy. Notice that if a stochastic preliminary model exists, this approach can immediately exploit it, contrary to traditional NMPC approaches. This finishes the in silico part of the algorithm, subsequent steps would be run in the true system. Therefore, emphasis after this step is given on sampling as least as possible, as every new sample would result in a new batch run from the real plant.

Step 2-3, Transfer Learning: The policy could directly be retrained using the true system and adapt all the weights according to the Reinforce algorithm. However, this may result in undesired effects. The control policy proposed in this work has a deep structure, as a result a large number of weights could be present. Thus, the optimization to update the policy may easily be stuck in a low-quality local optima or completely diverge. To overcome this issue the concept of transfer learning is adopted. In transfer learning, a subset of training parameters is kept constant to avoid the use of a large number of epochs and episodes, applying knowledge that has been stored in a different but related problem. This technique is originated from the task of image classification, where several examples exists, $e.g.$ in [26], [44], [18].

Using transfer learning, the current work only retrained the last hidden layers, and the policy is able to adapt to new situation without losing previously obtained knowledge, as shown in  Fig.5. Alternatively, an additional set of layers could be added on the top of the network.

Step 4, Transfer Learning Reinforce (On-line): In this step, Algorithm 1 is applied again, but now, on the true system. This algorithm aims to maximize a given reward (e.g. product concentration, economic objective)

Step 5: Terminate policy update and output $\Theta$ that defines the optimal RNN policy.

The methodology is described in Algorithm 2 and depicted in Fig. 6.

3 Computational Case Studies

In this section three case studies are presented to illustrate the effectiveness of the proposed batch-to-batch strategy. Our strategy is applied to 3 fed-batch bioreactors, where the objective is to maximize the concentration of a target product ($y_{2}$ or $c_{q}$) at the end of the batch time, using light and an inflow rate ($u_{1}$ or $I$ and $u_{2}$ or $F_{N}$) as manipulated variables.

3.1 Case Study 1 - Ordinary Differential Equations

In the first case study, the “real” photo-production system (plant) is described by the following equations plus an additional random disturbance:

[TABLE]

where ${u}_{1}$, ${u}_{2}$ and $y_{1}$, $y_{2}$ are the manipulated variables and the outlet concentrations of the reactant and product, respectively. The batch operation time course is normalized to 1. Additionally, a random disturbance is assumed, which is given by a Gaussian distribution with mean value 0 and standard deviation 0.02 on the states $y_{1}$ and $y_{2}$. We discretize the time horizon into 10 intervals of the dimensionless time, with one constant control input in each interval, resulting in a total of 20 control inputs.

The exact model is usually not known, and a simplified deterministic model is assumed according to some set of parameters. This preliminary model, given in (19 - 20), is utilized in an extensive offline training in order to construct the control policy network. As illustrated in the previous section 2.4, there is a potential to have a close approximation of the solution of the OCP using the RNN-Reinforce.

[TABLE]

The training consist of 100 epochs and 800 episodes using the simplified model to search the optimal control policy that maximizes the reward of (19 - 20) in this case the concentration of $y_{2}$ at the final instance (21).

[TABLE]

The control actions are constrained to be in the interval $[0,5]$. The control policy RNN is designed to contain 2 hidden layers, each of which comprises 20 neurons embedded by a hyperbolic tangent activation function. It was observed that 2 hidden layers are sufficient to approximate the optimal control policy, however there is the potential to use deeper structures with more layers for more complex systems. Furthermore, we employed two policy networks instead of one for simplicity. This approach assumes that the two manipulated variables are independent resulting in diagonal variance.

The algorithm is implemented in Pytorch version 0.4.1. Adam [24] is employed to compute the network’s parameter values using a step size of $10^{-2}$ with the rest of hyperparameters at their default values. After the training, using the simplified model the reward has the same value with the one computed by the optimal control problem, as expected. It should be noted that the computation cost of the control action using the policy is insignificant since it only requires the evaluation of the corresponding RNN, and does not depend directly on the complexity or the number of variables. In contrast, the solution of the OCP scale very badly with respect to both the complexity and the number of variables. Precisely, the maximum rewards for RL and OCP for both cases is 0.64. The reward for its epoch is illustrated in Fig. 8 and the process trajectories after the final update of the policy networks are shown in Fig. 7.

Fig. 8 shows that the reward has a large variance at the beginning of the training but is undetectable at the end. This can be explained as the trade-off between exploration and exploitation, where initially there is a lack of information and policy explores possible control actions, while at the end the policy exploits the available information. This policy can be considered as an initialization of the Reinforce algorithm which uses transfer learning to incorporate new knowledge gained from the true plant (Steps 3-5 in Algorithm 2). New data-sets from 25 batches are used (i.e. 25 real plant epochs) to update the true plant’s RL policy. The solution after only 4 epochs is 0.591 while the stochastic-free optimal solution identified using the unknown (complex) model of the plant is 0.583. This results show that the stochastic nature of the system can also affect the performance. The reward for each epoch is depicted in Fig. 10 and the process trajectories after the last epoch are depicted in Fig. 9. Notice that even before having any interaction with the “real” system the proposed approach has a superior performance than NMPC. This is because RL directly accounts for the effect of future uncertainty and its feedback in a proper closed-loop manner, whereas NMPC assumes open-loop control actions at future time points in the prediction.

There is a variation on the results after the last batch upon which the policy is updated. This makes sense, since the system appears to have some additive noise (i.e. Gaussian disturbance) and the policy maintains its stochastic nature.

The results are also compared with the use of NMPC using shrinking horizon. The results can be seen in Fig. 10, where 100 Monte-Carlo simulations were conducted. The optimization using our approach appears to be superior to the one given by the NMPC, showing the significance of our result. Furthermore, it should be noted that the performance of our proposed policy is better even in epoch 1, before the adaptation is started. In addition, in Fig. 11, the comparison between the control inputs of that are computed using our approach and the NMPC.

3.2 Case Study 2 - Stochastic Differential Equations

In this case study the same type of reaction is assumed to follow a stochastic differential equations:

[TABLE]

where $W$ is Wiener stochastic process. The simplified model is assumed to be the same with the previous case study. As a result the same policy that is trained off-line is used here. The purpose of this case study is to observe how the same policy can adapt in different environments. Now the model that describes the real system is not only structurally different, but also stochastic in nature.

The same hyperparmeters and networks are utilized for the policies in both stages, in order to show that the same policy can adapt to different environments successfully.

The same validation is conducted here using 100 Monte-Carlo simulations. Through comparison, our approach is found to be superior to the NMPC. In this case, our proposed algorithm adapts more rapidly to the new conditions, reducing significantly the requirement for a large number of episodes and epochs, as it can be seen in Fig. 12. This is attributed to the systematic transfer learning proposed in our algorithm. The computationally intensive part has been shifted off-line where the preliminary inaccurate model was used to train the policy. Then the (deep) recurrent neural network adapts successfully to the new environment that consists of a system of stochastic differential equations. The comparison is also depicted in Fig. 13. This result is also observed in the previous case study where the stochastic part of the physical system has a different nature. The control inputs are depicted in Fig. 14.

It should be noted that in both case studies the NMPC produced very similar control actions, with the only difference being the variance, compared to our approach which shapes the control actions to fit the needs of the different dynamics and uncertainty.

The methods used in the Reinforce algorithm usually require substantial number of episodes and epochs, therefore a good initial solution in combination with transfer learning is paramount so that Step 4 can be completed with a few batch-to-batch runs. In order to keep the problem realistic, only a small number of batches is utilized in Step 2-3 to refine the policy network.

3.3 Case Study 3 - Nonsmooth Model

The last case study in this paper focuses on the photo-production of phycocyanin synthesized by cyanobacterium Arthrospira platensis. Phycocyanin is a high-value bioproduct and its biological function is to enhance the photosynthetic efficiency of cyanobacteria and red algae. It has applications as a natural colorant to replace other toxic synthetic pigments in both food and cosmetic production. Additionally, the pharmaceutical industry considers it as beneficial because of its unique antioxidant, neuroprotective, and anti-inflammatory properties.

Both the “real” and simplified model in this case are considered to be nonmooth. Due to different growth phases, nonmooth behaviour is observed for the physical system. To accommodate this difficulty, switching functions have been proposed [15, 57]. In this work the nonsmooth behaviour is modelled using a $\text{sign}(\cdot)$ function. The “real” dynamic system consists of three nonsmooth ODEs describing the evolution of the concentration of biomass ($X$), nitrate ($N$), and product ($q$). The dynamic model is based on Monod kinetics, which describes microorganism growth in nutrient sufficient cultures, where intracellular nutrient concentration is kept constant because of the rapid replenishment. We assume a fixed volume fed-batch. The manipulated variables as in the previous examples are the light intensity ($I$) and inflow rate ($F_{N}$). The mass balance equations are

[TABLE]

where the parameters are given in Table 1. The real physical system consists of additive disturbance

[TABLE]

and measurement noise

[TABLE]

Additionally, uncertainty is assumed for the initial concentration, where

[TABLE]

The reward is additionally penalized by the change of the control actions $\textbf{u}(t)=\left[I,F_{N}\right]^{T}$. As a result the reward is given as:

[TABLE]

where $\Delta\textbf{u}_{t}=\textbf{u}_{t}-\textbf{u}_{t-1}$.

The simplified deterministic model is assumed without the noise or the additive disturbance. This preliminary model, is utilized in an extensive offline training in order to construct the control policy network. As illustrated in the previous section 3.4, there is a potential to have a close approximation of the solution of the OCP using RNN-Reinforce.

The training consists on 100 epochs and 500 episodes and the optimal control policy that maximizes the reward in equation 30. The control actions are constrained to be in the interval $0\leq F_{N}\leq 40$ and $120\leq T\leq 400$. The control policy RNN is designed to contain 4 hidden layers, each of which comprises 20 neurons embedded by a leaky rectified linear unit (ReLU) activation function. Furthermore, in this case a unified policy network with diagonal variance is utilized such that the control actions share memory and the previous states are used from the RNN (together the current measured states).

The algorithm is implemented in Pytorch with the same configurations. It should be noted that the computational cost of computing the control action online is insignificant since it only requires the evaluation of the corresponding RNN, and does not depend directly on the complexity or the number of variables. In contrast, the solution of the OCP scales badly in this case due to the presence of integer variables. The reward for each epoch is depicted in Fig 15. In this case the probability density is shown due to the nonsmoothness of the model, that may result in multiple peaks. The lines are faded out towards earlier epochs. Additionally, there is no guarantee of global optimality in the current work, as a result the reward may get stuck in other local minima. It should be noted that in this case the uncertain initial conditions are applied during this training phase.

The nominal control actions are depicted in Fig. 16, where the shaded areas are the 98% and 2% percentiles. The corresponding states are depicted in Fig.17 with their 98% and 2% percentiles. The nominal behaviour is subject to the corresponding initial conditions since no other uncertainty is taken into account in the offline procedure. The probability density of the product $c_{q}$ has clearly only one peak, this is shown in Fig. 17(c).

As in the previous case studies, the results are compared with the use of NMPC using shrinking horizon. The optimization is a mixed integer nonlinear programming problem (MINLP). Local optimization is used in order to be numerically tractable. Orthogonal collocation is implemented and integer variables have been used to model the switches. It should be noted that this MINLP takes 2-4 mins to be solved.

The results can be seen in Fig. 18 with their 98% & 2% percentile respectively, where 100 Monte-Carlo simulations were conducted. The optimization using our approach appears to be superior to the one given by the NMPC, showing the significance of our result. After the adaptation the probability densities are depicted in Fig. 19. Next, the control inputs are depicted in Fig. 20 with their 98% & 2% percentile respectively. It is clear that the NMPC control actions have large variance compare to the ones produce by our proposed methodology. This is due to the nonsmoothness of the model and the uncertainty which the NMPC struggles with.

In addition, in Fig. 20, the comparison between the control inputs of our approach and the NMPC is presented.

4 Conclusions and Future Work

In this work we propose a new methodology for batch-to-batch learning by adapting Reinforcement learning techniques to uncertain and complex bioprocesses. The results reveal that it is possible to obtain a near optimal policy for a stochastic system when the true dynamics are unknown. In real systems with the absence of a true model, it is impossible to generate highly accurate datasets to train the policy network. As a result we propose a 2-stage framework where first an approximate (possibly stochastic) model is used to train the policy network. Subsequently, this policy is implemented into the true system. In this way, there is no need for a large number of evaluations of the true system which can be costly and time consuming.

A systematic adaptation to the new environment is achieved using transfer learning. In Step 4: Transfer Learning Reinforce the policy is trained using $T<<T_{0}$ episodes conducting the Steps 1-3 of Algorithm 1. The proposed algorithm is validated using two case studies for different nature of stochastic processes. Our proposed methodology results in a policy that overcomes the performance of the NMPC, where only simple policy evaluations are needed.

The off-line CPU time is 3 hours, however the online implementation of the needs only 0.002 secs. This means that all the computational complexity is shifted offline and an efficient optimal control policy is constructed. One should also keep in mind that after the off-line training the solution to a nonlinear stochastic dynamical system is provided, in the form of a stochastic policy. This is a more complete and efficient solution as it is a closed-loop solution, rather than an open-loop optimization. Furthermore, a nonsmooth system was integrated with Casadi [1], which is more time consuming than integrating a smooth dynamic system.

For both the case studies 4 epochs and 25 batches were implemented. In this work, the training was stopped after the $4^{th}$ epoch, but the training could have been continued or stopped earlier. Here, the total number of batches is $4\times 25=100$; however, the policy for all case studies performs better from the beginning of the online implementation. This means that a smaller number of batches can be used and still outperform NMPC.

We emphasize that our considered systems contain both stochasticity and plant-model mismatch, and there is no process structure available. The optimisation of such systems is generally known to be intractable. Given the early stage of this research, there are still disadvantages of this method which must be accommodated in the future, including the robust satisfaction of constraints. In addition, there is a wide discussion regarding the safety in reinforcement learning  [54], which is also a result of the difficulty of robust satisfaction of constraints. Future work will focus on the robust satisfaction of constraints in RL methods.

The codes are available at: https://gitlab.com/Panos108/rl-with-nonsmooth

Acknowledgements

This project has received funding from the EPSRC project (EP/P016650/1).