# Neural Network-Based State Estimation for a Closed-Loop Control Strategy Applied to a Fed-Batch Bioreactor

Authors: Santiago Rómoli; Mario Serrano; Francisco Rossomando; Jorge Vega; Oscar Ortiz; Gustavo Scaglia
Journal: Complexity
 (2017)

---

## Abstract

The lack of online information on some bioprocess variables and the presence of model and parametric uncertainties pose significant challenges to the design of efficient closed-loop control strategies. To address this issue, this work proposes an online state estimator based on a Radial Basis Function (RBF) neural network that operates in closed loop together with a control law derived on a linear algebra-based design strategy. The proposed methodology is applied to a class of nonlinear systems with three types of uncertainties: (i) time-varying parameters, (ii) uncertain nonlinearities, and (iii) unmodeled dynamics. To reduce the effect of uncertainties on the bioreactor, some integrators of the tracking error are introduced, which in turn allow the derivation of the proper control actions. This new control scheme guarantees that all signals are uniformly and ultimately bounded, and the tracking error converges to small values. The effectiveness of the proposed approach is illustrated on the basis of simulated experiments on a fed-batch bioreactor, and its performance is compared with two controllers available in the literature.


---

## Body

## 1. Introduction

In the last decade, several control strategies applied to different bioprocesses have been extensively investigated (e.g., [1–5]). In the particular case of fed-batch bioreactors, some nonlinear control strategies were developed to improve the process efficiency and the “batch-to-batch” reproducibility [6], as well as to deal with the unavoidable uncertainties typically present in bioprocesses [7].Many articles have investigated the effects of model uncertainties on the performance of model-based nonlinear control strategies typically applied to bioprocesses. For example, the “yield coefficient of cell to glucose” in a bioprocess is a time-varying parameter often assumed as a constant. Also, the “specific growth rate” can be represented through the Monod-Haldane model [8], but it is possible that the model becomes inaccurate because of yeast growth inhibition due to an eventual presence of ethanol. In principle, three main types of uncertainties can be envisaged [9–13]: (1) uncertain and/or time-varying model parameters, (2) uncertain nonlinearities, and (3) unmodeled dynamics. Additionally, nonmodeled external perturbations can disturb the process. In general, model uncertainties and external perturbations can cause severe risks in bioreactors, such as leading to undesirable operating conditions [14–16]. For instance, inadequate modeling of the inhibitory effect of a substrate can lead to a structural instability in the dynamical behavior of a bioprocess [15].Adaptive control strategies are a valuable choice to control bioreactors due to their ability to compensate for parametric uncertainties. The Nussbaum gain method has been applied to control a bioprocess of unknown and time-varying control gain and parameters [17]. A fed-batch bioreactor control with simultaneous dynamic identification of unknown state and parameters has been investigated through an algorithm based on the auxiliary model principle with adaptive controllers [18]. A control strategy was proposed to track a desired profile for the biomass of a yeast culture in a high-density fed-batch bioreactor [19]. This last strategy was based on a model that included an uncertainty term estimated through an adaptive technique. Also, an adaptive nonlinear feedback controller has been designed to automatically estimate uncertain time-varying parameters in order to guarantee the closed-loop control of a waste water plant [20].Fuzzy control systems have been successfully applied to several bioreactors. For instance, [12] considered an adaptive fuzzy-decentralized control for a system with unknown nonlinear uncertainties and dynamical uncertainties. The fuzzy system was used to estimate the nonlinear uncertainties. Also, a fuzzy controller for the alcoholic fermentation withSaccharomyces cerevisiaehas been proposed [21]. It combined the capability of handling uncertainties with the ability of predictive control to improve the plant performance, making use of a neural network model of the nonlinear system. Optimal control strategies have also been applied to bioprocesses [22, 23]. In order to compensate for process uncertainties, this type of controllers incorporates an additional term to the bioprocess model (typically, in the specific growth rate), which is estimated through an optimization technique. Alternative strategies for controlling bioreactors under model uncertainties are trajectory-based control [24], hybrid control [25], and model predictive control [26]. Additionally, due to the inherent nature of a bioprocess, a control scheme based on a stochastic model seems to be adequate. The specific growth rate can be modeled as affected by a white noise representing environment uncertainties, and the mathematical model consists of a set of stochastic differential equations [27] (some applications can be consulted in [28, 29]). The application of this approach onto a real plant is quite difficult due to its complex design and implementation.Based on a linear algebra methodology, a control strategy has been applied for trajectory tracking in a bioprocess [30]. The technique was used to calculate the control actions that ensure the tracking of independently determined optimal trajectories of some process variables. To this effect, a discrete approximation of the original process model must first be derived, and the control actions are then calculated through a least square problem. This methodology has been applied to a fed-batch bioprocess designed to produce a secreted protein, and it was originally modeled by Park and Ramirez [31]. Optimal profiles of such bioprocess have also been published [32].Many nonlinear control techniques of bioprocesses assume that the feedback variables can be measured online [33, 34]. However, such hypothesis is often unrealistic and therefore some process variables must be estimated through virtual sensors (or soft sensors). There are several kinds of virtual sensors published in the literature which exhibit different abilities to deal with model mismatches and perturbations. In this work, an artificial intelligence-based observer has been chosen because this type of estimator is especially recommended for highly nonlinear systems in the presence of uncertainties [35].This paper aims at designing a controller for the multivariable tracking of optimal profiles when the process is affected by uncertainties in the dynamic model and/or in the model parameters. Unmeasured process variables are estimated through a RBF neural network. The proposed approach is evaluated onto a bioprocess originally modeled by [31], with the optimal profiles determined by [32] in absence of model errors. Brief descriptions of the fed-batch bioreactor and the controller derived in [30] are presented in Appendix A.In comparison with [30], the proposed approach prevents the polynomial uncertainties from affecting the tracking control. To achieve this aim, an integrative term is incorporated in the linear algebra methodology originally proposed. Then, the conditions in order for a system of linear equations to have exact solution are analyzed. The control action is obtained by solving the linear system, even though the original process model is nonlinear.The main contribution is that the original method proposed herein ensures the convergence to zero of the tracking errors under additive uncertainties. Another important contribution is the inclusion of RBF neural network estimators: the control issue is addressed considering more realistic conditions with nonmeasurable variables. Further, the stability analysis is made using Lyapunov theory in discrete time.The article is organized as follows. In Section2, the controller design methodology and the convergence to zero of the tracking errors are presented. Section 3 shows briefly the description of the RBF estimator. A method to tune the controller is presented in Section 4, along with four testing cases designed to check the effectiveness of the new control law. Main conclusions of the work are summarized in Section 5.


## 2. Controller Design

Derivation of the control strategy is based on the mathematical model presented in AppendixA. The effect of an additive uncertainty (En) is modeled by extending (A.4) as follows: (1)zn+1=zn+cnΔt-dnunΔt+En,where zn is the state vector, cn and dn are the input vectors, Δt is the sample time, un is the control action, and En is an additive uncertainty. In principle, En can be used to model either perturbed systems or a wide class of model mismatches. In what follows, a control law is designed to compensate for the tracking errors introduced by En.Taking into account the fact that the mismatch might depend on both system states and inputs, we here consider a real plant described byzn+1=f(zn,un). The additive uncertainty can be expressed as En=fzn,un-f^zn,un, where f^zn,un is the discrete-time nonlinear system model. Ifz andu are assumed to be bounded andf is Lipschitz, then En can be modeled as a bounded uncertainty [36–38].With the additive uncertainty of (1), the fed-batch bioreactor model of (A.3) becomes(2)PnVnTnVnXnVnSF-SnVnun∗=χnTn-Pn-Pn+1-PnΔtψnXn-Tn+1-TnΔtμnXn-Xn+1-XnΔt7.3μnXn+Sn+1-SnΔt+EP,nET,nEX,nES,n.The control action for the secreted protein concentration, P, given by (A.4) of [30], is(3)un,P=χTn-Pn-Pref,n+1-kpPref,n-Pn-Pn/ΔtPn/Vn.The discrete value of P, Pn+1, is calculated by replacing (3) into the first row of (2):(4)Pn+1=Pref,n+1-kPPref,n-Pn+ΔtTn-Pnχn-χ+ΔtEp,n,where χn and χ depend on Sn and Sez,n, respectively. Thus, the tracking error of Pn is(5)eP,n+1=Pref,n+1-Pn+1.By combining (4) and (5), one obtains(6)eP,n+1=kPeP,n-ΔtTn-Pnχn-χ-ΔtEp,n.The Taylor interpolation of χn at its desired value χ is(7)χn=χ+dχdSS=Sez,n+λSn-Sez,n=SλSn-Sez,n;0<λ<1.By replacing (7) into (6), the tracking error forPn becomes(8)eP,n+1=kPeP,n-ΔtTn-Pnχ+dχdSS=SλSn-Sez,n-χ-ΔtEp,n,eP,n+1=kPeP,n+ΔtTn-PndχdSS=SλSez,n-Sn-ΔtEp,n.If a similar procedure is followed for T, X, and S, then the error equation becomes [30](9a)en+1=eP,n+1eT,n+1eX,n+1eS,n+1=kP0000kT0000kX0000kSeP,neT,neX,neS,n+ΔtTn-PndχdSS=SλXndψdSS=SφXndμdSS=Sδ0eS,n-EP,nET,nEX,nES,n.Alternatively, in vector form,(9b)en+1=Ken+hn-En,where hn is a nonlinearity that tends to zero when eS,n tends to zero. Equations (9a) and (9b) suggest a direct effect of the additive uncertainty on the tracking error. Therefore, the presence of an uncertainty avoids a convergence of tracking errors to zero (see Theorem  1 of [30]).


### 2.1. Single Integral Action

In order to reduce the effect ofEn, some integrators for the tracking errors will be introduced, depending on the time variation hypothesis of En. We will assume that En is unknown and each of its components is an m-order polynomial.Remark 1.
The first-order difference ofEn is defined as δEn=En+1-En, the second-order difference is defined as δ2En=δδEn=δEn+1-En=En+2-2En+1+En, and, as a rule, the q-th order difference is defined as δqEn=δδq-1En.Remark 2.
Theq-th difference of a q-1-order polynomial is zero.Consider first a constant uncertaintyEn=const. Then, δEn=En+1-En=0. In this case, an integrator for each state variable will force the error to converge to zero. Call e(t) the continuous time error in the state vector. Then, the integral of the tracking error is defined as(10)U1,n+1=U1,n+∫nΔtn+1Δtetdt≅U1,n+enΔt.Hence, the control action (A.7) will be computed by assuming a new term in (A.5); that is,(11)zn+1=zref,n+1-Kzref,n-zn+k1U1,n+1,k1=kP10000kT10000kX10000kS1,whereK and k1 are the proportional and integral design parameters, respectively.By adding one integrator(k1U1,n+1) to (9b) and using (10), one obtains(12)en+1+k1U1,n+1+enΔt=Ken+hn-En.Rewriting (12) for the next sample time,(13)en+2+k1U1,n+2+en+1Δt=Ken+1+hn+1-En+1.Combining (10) and (12) and substituting them into (13),(14)en+2+k1en+1Δt+hn-En-en+1+k1en=Ken+1+hn+1-En+1.Alternatively,(15)en+2+-K-Ien+1+k1en+1Δt+Ken=-hn+1-hn-En+1-En.Therefore,K and k1 are chosen in order to ensure the stability of the linear system represented in the left-hand side of (15). To achieve this condition, the roots of such polynomial (ri) must lie inside the unit circle. Then, en→0 when n→∞; that is, the error tends to zero as long as uncertainties are constant.Once the integral of the error has been added, the control action is calculated following the same design procedure based on linear algebra presented in [30]. Therefore, the new control law un can be obtained through the least squares procedure as (16)u1n=PnVnχTn-Pn-Pref,n+1-kPPref,n-Pn-Pn+kP1UP1,n+1/Δt+TnVnψXn-Tref,n+1-kTTref,n-Tn-Tn+kT1UT1,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+XnVnμXn-Xref,n+1-kXXref,n-Xn-Xn+kX1UX1,n+1/Δt+20-SnVn7.3μXn+Sez,n+1-kSSez,n-Sn-Sn+kS1US1,n+1/ΔtPn2+Tn2+Xn2+20-Sn2.


### 2.2. Multiple Integral Actions

Assume that the uncertainty can be modeled by a function with a second-order difference equal to zero; that is,δ2En=δδEn=δEn+1-En=En+2-2En+1+En=0. Similar to (10), a double integrator is now introduced, and new variables U1 and U2 are defined:(17)U2,n+1=U2,n+∫nΔtn+1ΔtU1tdt≅U2,n+ΔtU1,n+1.In this case, the control action (A.7) will include an additional term in (A.5) as follows:(18)zn+1=zref,n+1-Kzref,n-zn+k1U1,n+1+k2U2,n+1,k2=kP20000kT20000kX20000kS2,whereK,k1, andk2 are the proportional, integral, and double-integral design parameters, respectively. With a similar procedure to that used for deriving (15) and after adding two integrators (k1U1,n+1 and k2U2,n+1), the error dynamics can be expressed as (19)en+3+-K-2Ien+2+k1+k2Δten+2Δt+2K+Ien+1-k1en+1Δt-Ken=hP,n+2-2hP,n+1+hP,nhT,n+2-2hT,n+1+hT,nhX,n+2-2hX,n+1+hX,n0-δ2En.As suggested by (19), under the assumption of a constant or linear varying uncertainty (δ2En=0), the uncertainty has no influence on the error dynamics. As in the case of the single integrator, the controller parametersK,k1, andk2 should be chosen in order to ensure the stability of the left-hand side of (19).The previous derivation can be generalized for uncertainties that can be approximated with aq-1-order polynomial; and, therefore, the influence of En upon en would be cancelled off by addition of q integrators.The linear algebra-based procedure enables the calculation of the following control law:(20)u2n=PnVnχTn-Pn-Pref,n+1-kPPref,n-Pn-Pn+kP1UP1,n+1+kP2UP2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+TnVnψXn-Tref,n+1-kTTref,n-Tn-Tn+kT1UT1,n+1+kT2UT2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+XnVnμXn-Xref,n+1-kXXref,n-Xn-Xn+kX1UX1,n+1+kX2UX2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+20-SnVn7.3μXn+Sez,n+1-kSSez,n-Sn-Sn+kS1US1,n+1+kS2US2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2.


## 2.1. Single Integral Action

In order to reduce the effect ofEn, some integrators for the tracking errors will be introduced, depending on the time variation hypothesis of En. We will assume that En is unknown and each of its components is an m-order polynomial.Remark 1.
The first-order difference ofEn is defined as δEn=En+1-En, the second-order difference is defined as δ2En=δδEn=δEn+1-En=En+2-2En+1+En, and, as a rule, the q-th order difference is defined as δqEn=δδq-1En.Remark 2.
Theq-th difference of a q-1-order polynomial is zero.Consider first a constant uncertaintyEn=const. Then, δEn=En+1-En=0. In this case, an integrator for each state variable will force the error to converge to zero. Call e(t) the continuous time error in the state vector. Then, the integral of the tracking error is defined as(10)U1,n+1=U1,n+∫nΔtn+1Δtetdt≅U1,n+enΔt.Hence, the control action (A.7) will be computed by assuming a new term in (A.5); that is,(11)zn+1=zref,n+1-Kzref,n-zn+k1U1,n+1,k1=kP10000kT10000kX10000kS1,whereK and k1 are the proportional and integral design parameters, respectively.By adding one integrator(k1U1,n+1) to (9b) and using (10), one obtains(12)en+1+k1U1,n+1+enΔt=Ken+hn-En.Rewriting (12) for the next sample time,(13)en+2+k1U1,n+2+en+1Δt=Ken+1+hn+1-En+1.Combining (10) and (12) and substituting them into (13),(14)en+2+k1en+1Δt+hn-En-en+1+k1en=Ken+1+hn+1-En+1.Alternatively,(15)en+2+-K-Ien+1+k1en+1Δt+Ken=-hn+1-hn-En+1-En.Therefore,K and k1 are chosen in order to ensure the stability of the linear system represented in the left-hand side of (15). To achieve this condition, the roots of such polynomial (ri) must lie inside the unit circle. Then, en→0 when n→∞; that is, the error tends to zero as long as uncertainties are constant.Once the integral of the error has been added, the control action is calculated following the same design procedure based on linear algebra presented in [30]. Therefore, the new control law un can be obtained through the least squares procedure as (16)u1n=PnVnχTn-Pn-Pref,n+1-kPPref,n-Pn-Pn+kP1UP1,n+1/Δt+TnVnψXn-Tref,n+1-kTTref,n-Tn-Tn+kT1UT1,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+XnVnμXn-Xref,n+1-kXXref,n-Xn-Xn+kX1UX1,n+1/Δt+20-SnVn7.3μXn+Sez,n+1-kSSez,n-Sn-Sn+kS1US1,n+1/ΔtPn2+Tn2+Xn2+20-Sn2.


## 2.2. Multiple Integral Actions

Assume that the uncertainty can be modeled by a function with a second-order difference equal to zero; that is,δ2En=δδEn=δEn+1-En=En+2-2En+1+En=0. Similar to (10), a double integrator is now introduced, and new variables U1 and U2 are defined:(17)U2,n+1=U2,n+∫nΔtn+1ΔtU1tdt≅U2,n+ΔtU1,n+1.In this case, the control action (A.7) will include an additional term in (A.5) as follows:(18)zn+1=zref,n+1-Kzref,n-zn+k1U1,n+1+k2U2,n+1,k2=kP20000kT20000kX20000kS2,whereK,k1, andk2 are the proportional, integral, and double-integral design parameters, respectively. With a similar procedure to that used for deriving (15) and after adding two integrators (k1U1,n+1 and k2U2,n+1), the error dynamics can be expressed as (19)en+3+-K-2Ien+2+k1+k2Δten+2Δt+2K+Ien+1-k1en+1Δt-Ken=hP,n+2-2hP,n+1+hP,nhT,n+2-2hT,n+1+hT,nhX,n+2-2hX,n+1+hX,n0-δ2En.As suggested by (19), under the assumption of a constant or linear varying uncertainty (δ2En=0), the uncertainty has no influence on the error dynamics. As in the case of the single integrator, the controller parametersK,k1, andk2 should be chosen in order to ensure the stability of the left-hand side of (19).The previous derivation can be generalized for uncertainties that can be approximated with aq-1-order polynomial; and, therefore, the influence of En upon en would be cancelled off by addition of q integrators.The linear algebra-based procedure enables the calculation of the following control law:(20)u2n=PnVnχTn-Pn-Pref,n+1-kPPref,n-Pn-Pn+kP1UP1,n+1+kP2UP2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+TnVnψXn-Tref,n+1-kTTref,n-Tn-Tn+kT1UT1,n+1+kT2UT2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+XnVnμXn-Xref,n+1-kXXref,n-Xn-Xn+kX1UX1,n+1+kX2UX2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2+⋯+20-SnVn7.3μXn+Sez,n+1-kSSez,n-Sn-Sn+kS1US1,n+1+kS2US2,n+1/ΔtPn2+Tn2+Xn2+20-Sn2.


## 3. Neural State Estimator

The proposed feedback control technique requires reliable estimates of several state variables. The state estimation error,x~n, is defined as(21)x~n=xn-x^n,where x^n is the vector of estimated state variables and xn=[P^n,T^n,X^n] is the vector of online measured states. The following neural estimator represents the nonlinear dynamic of the bioreactor described in (7):(22)xn+1=θ∗Tξxn,υn+εn,x0=x0,where υn=[un,Sn,Vn] is the input vector to the neural estimator, θ∗ is the optimal weight vector, ε is the neural approximation error, ξi is the RBF that represents each neuron in the hidden layers, subindex i indicates the neuron number of the RBF, and m is the maximum number of neurons. Every neuron of the hidden layer is modeled as(23)ξixn,υn=exp-12σi2xnυn-ci2,i=1,2,…,m,where ci and σi are two parameter vectors that represent the centers and widths of the RBF, respectively.In the bioreactor under study, biomass(X) cannot be measured online and the protein concentrations (P,T) cannot be measured at all; then a state estimator function based on (22) is determined as follows:(24)x^n+1=θ^nTξxn,υn.From the difference between (22) and (24), the neural identification error can be described by(25)x~n+1=xn+1-x^n+1=θ∗Tξxn,υn-θnTξxn,υn+εn=θ~nTξxn,υn+εn,where θ~n represents the neural weight vector estimation error and it is defined as(26)θ~n=θ∗-θn.In order to train the neural network for identification, an offline data set of (x,u) will be used. This training data will be obtained after following the same procedure presented in [39] and then applied in [40, 41]. A similar procedure was also used in [42]. The learning rule utilized to train the neural network is demonstrated in the next theorem.Theorem 3.
The model of the fed-batch bioreactor (see (A.1) and (2)) can be approximated by the neural network (24) using a neuronal adjustment law defined by(27)Δθ~n=-γξxn,υnx~nTΛ.Proof.
The demonstration of this theorem is presented in AppendixB.As recommended by [35], the estimator was evaluated to test its performance and avoid biased estimates. Table 1 summarizes the implemented tests, and their results were evaluated through the average values of the mean squared errors (MSE) corresponding to P, T, and X. The obtained results show an acceptable performance of the RBF neural network for estimating the nonmeasurable states. Note that none of the MSE values exceeds the thresholds proposed by [43], 1×10-3. Figure 1 shows the proposed closed-loop architecture. Note that the controller requires information on both measured and estimated variables.Table 1
Performance evaluation of the RBF neural network.
Test
M
S
E
T
M
S
E
P
M
S
E
X
Normal operating conditions
1.58 × 10−6
2.14 × 10−6
7.37 × 10−4
Continuous variations in the substrate concentration (15%)
2.00 × 10−6
2.47 × 10−6
8.17 × 10−4
20% of error in model parameters (SF and YS/X)
1.98 × 10−6
2.70 × 10−6
9.51 × 10−4
20% lesser in the initial nominal operating conditions inV, S, and X
1.83 × 10−6
2.25 × 10−6
3.80 × 10−4
Normally distributed noise with zero mean and variance of 0.2 inu and S
1.93 × 10−6
2.35 × 10−6
7.51 × 10−4Figure 1
The proposed closed-loop control architecture: the linear algebra-based controller and the Radial Basis Function (RBF) neural network estimator.


## 4. Results and Discussion

All simulation examples are based on the bioprocess originally modeled by [31]. The available optimal profiles for such process are also adopted [32]. Equation (A.1) represents the fed-batch bioreactor model. The effectiveness of the new control law is tested through some simulation examples. First, the optimal controller parameters are tuned through a Monte Carlo experiment. Then, a sinusoidal perturbation is added to some model parameters. Also, the controller performance under uncertain time-varying parameters is verified through a Monte Carlo Randomized Algorithm. In other examples, the functionality of the specific growth rate is changed. Finally, some parameters are contaminated by additive white noise.As previously discussed, the behavior of the controlled system depends on parametersk,K1,…,Kp, where p is the number of integrators. For p=0,1,2, the controller parameters are determined through a Monte Carlo method, where the number of simulations is determined by application of Theorem 2 of [30]. All simulations are implemented in Matlab™ (MathWorks, Inc.), with a sample time Δt=0.01h. For each controller (i.e., for each p), 1,000 simulations are run; and the controller parameters are randomly chosen, such that the roots (ri) of the linear systems defined in the right-hand side of (9a) and (9b) and left-hand side of (15) and (19) verify ri∈(0,1). More specifically, for p=0, the controller parameters are directly taken from [30]. For p = 1, 1,000 new simulations are performed, and each controller parameter is chosen according to(28)ki=r1r2,K1i=-r1-r2+r1r2+1Δt,where the superscript i refers to each state variable and r1 and r2 are the roots of the linear system in the left-hand side of (15). For p = 2, 1,000 simulations are performed, and each controller parameter is chosen according to(29)ki=r1r2r3,K1i=2r1r2r3+1-r1r2-r1r3-r2r3Δt,K2i=-r1-r2-r3-r1r2r3+1+r1r2+r1r3+r2r3Δt2,where r1, r2, and r3 are the roots of the linear system in (19).Table2 summarizes the results obtained with each controller. For comparison purposes, the performance of the proposed strategy is compared with that provided by a Parallel Two-Degree-of-Freedom PID Controller (2-DOF PID), which computes a weighted difference signal for each of the proportional, integral, and derivative actions according to the set point weights specified by the user [44]. In this case, the culture cell density has been chosen as the controlled variable. The 2-DOF PID was tuned with the automatic tuning of the Simulink™ Control Design software that provides the proportional, integral, and derivative gains, as well as the filter coefficient.Table 2
Main simulation results in the Monte Carlo experiments.
Controller
Controller parameters
p
=
0
K
P = 0.23
K
T = 0.14
K
X = 0.52
K
s = 0.23
p
=
1
K
P = 0.7758
k
P
1 = 0.043
K
T = 0.6181
k
T
1 = 0.083
K
X = 0.8173
k
X
1 = 0.008
K
s = 0.0356
k
s
1 = 0.013
p
=
2
K
P = 0.7758
k
P
1 = 0.053
k
P
2 = 0.037
K
T = 0.6181
k
T
1 = 0.084
k
T
2 = 0.052
K
X = 0.8173
k
X
1 = 0.097
k
X
2 = 0.065
K
s = 0.0356
k
S
1 = 0.099
k
s
2 = 0.015


### 4.1. Time-Varying Parameters and Input Disturbances

As in [18, 45], the following expressions are adopted for simulating uncertainties in the parameters and the feed flow rate:(30)Y^S/X=7.3+0.05sin⁡t,S^F=20+0.05sin⁡t,u^=u+0.05sin⁡t.The tracking error (en) is defined as(31)en2=Pref,n-Pn2+Tref,n-Tn2+Xref,n-Xn2+Sez,n-Sn2.In what follows, we shall call C0, C1, and C2 the controllers that use p = 0, p = 1, and p = 2 integrators, respectively, and C3 the 2-DOF PID controller. Particularly note that C0 corresponds to the controller designed in [30]. In spite of the assumed time-varying parameters, Figure 2(a) shows that C2 produces the minimum tracking error. According to Figure 2(b), both C2 and C1 are effective to reach the desired protein concentration profile, while C0 and C3 produce large errors.Figure 2
(a) Comparison of the tracking errors (see (31)) (left ordinate axis: C0 and C3; right axis: C2 and C1). (b) Total protein concentration profiles: C1 and C2 successfully reach the maximum production of protein in the presence of time-varying parameters.
(a)
(b)In several research fields, probabilistic methods have been found to be useful for dealing with problems related to robustness of systems affected by uncertainties [46]. Particularly, the Monte Carlo Randomized Algorithm has been used in many applications such as radioactive decay, power system generation, and traffic on roads [47]. In the control area, Monte Carlo methods are effective tools for the analysis of probabilistically robust control schemes [30, 46]. A Monte Carlo simulation is here performed to demonstrate the effectiveness of C2 from a statistical point of view [48–50], when uncertainties are introduced by time-varying parameter. To this effect, 1,000 random values of SF and YS/X are chosen in a range of ±25% of their nominal values. Figures 3(a), 3(c), 3(e), and 3(g) show that only C2 and C1 were able to adequately track the reference profile. In fact, both controller responses are quite similar, thus suggesting that the time-varying parameter has not affected their performances. For comparison, the tracking error was evaluated through the following performance index J:(32)J=0.25Δt∑n=1NrefPref,n-Pn2+Tref,n-Tn2+Xref,n-Xn2+Sez,n-Sn2Nref,where Nref represents the number of points of the reference profile. Figures 3(b), 3(d), 3(f), and 3(h) confirm a better performance of C2 and C1 over C0 and C3. In fact, J is quite small for C2 and C1, even though 25% of uncertainty in the main bioreactor parameters was adopted. Since 1,000 simulations with random values of SF and YS/X were carried out and the tracking errors remained bounded, the performance of C2 and C1 will be acceptable with 99% of probability whatever the parameter values are in a range of ±25% [46].Figure 3
Total protein concentration ((a), (c), (e), and (g)) and performance indexJ ((b), (d), (f), and (h)) for 1,000 random values of SF and YS/X (controllers are indicated in each figure).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


### 4.2. Uncertain Nonlinearities

Uncertain nonlinearities are often present in many bioprocesses. For instance, the specific growth rate,µ, can suffer from important changes from batch to batch, and unfortunately such variability is impossible to be accurately modeled. There are a number of models in the literature for μ [8], for example, the third expression of (A.2), for the bioprocess currently under study [31]. To demonstrate the controller performance, an additional term can be included in the calculation of μ as follows [18, 51]:(33)μ^=21.87+0.28sin⁡tSS+0.4S+62.5.Figure 4(a) compares the performances reached by C2, C1, C0, and C3 under the simulated uncertain nonlinearities in μ. The lower tracking errors correspond to C2. Figure 4(b) shows that only C2 and C1 are able to reach the maximum concentration of the total protein.Figure 4
(a) Comparison of the tracking errors (see (31)) under uncertain nonlinearities in μ (left ordinate axis: C0, C1, and C3; right axis: C2). (b) Total protein concentration profiles.
(a)
(b)


### 4.3. Unmodeled Dynamics

The presence of white noise inµmax leads to a mathematical model governed by a set of stochastic differential equations. In order to properly deal with this type of uncertainties, a kind of stochastic control approach should be used. As an alternative to the complex theory of stochastic control, a simpler deterministic strategy is here proposed. To this effect, note first that all parameters in (A.1) and (A.2) are truly affected by environmental noise, and consequently they always fluctuate at long times around some average values. For simplicity, we will consider that only µmax is affected by randomness and changes according to the following expression [52, 53]:(34)dμmax=μmax+αdBt,where dB(t) is a standard Brownian motion defined on a complete probability space and α ≥ 0 represents the intensity of the noisy signal.A white noise can be approximated through a band-limited white noise (BLWN), which is simply implemented through a proper block of Simulink. Theoretically, a white noise is uncorrelated; it exhibits a flat power spectral density, and therefore it is characterized by an infinite energy. The BLWN block can simulate the effect of white noise by using a random sequence with a correlation time much smaller than the shortest time constant of the system. The correlation time of the noise is the sample rate of the block [54] and can be estimated as(35)tc≈2π1001fmax=2π10010.5/Δtmax,where tc is the correlation time, fmax is the bandwidth of the system, and Δtmax is the maximum sample time for which the system is stable. For the feedback control presented in this work, Δtmax = 0.0145 h. Then, from (35), tc = 0.0018 h.Many researchers have utilized this methodology to prove the robustness of the bioreactors controllers [19, 22–24, 55]. In the current simulations, a BLWN with a noise power of 0.9 and α = 0.5 is introduced along with the feedback signal. Thus, the maximum specific growth rate affected with BLWN is calculated as(36)μ~max=μmax+0.05BLWN.Figure 5(a) indicates that C2 and C1 adequately work in the presence of unmodeled dynamics. In fact, the tracking errors are quite small in comparison with C0 and C3. Figure 5(b) shows the profile tracking for each controller. Again, C2 and C1 are able to follow the reference profile, while C0 and C3 fail to reach the maximum protein concentration.Figure 5
(a) Comparison of the tracking errors (see (31)) (right ordinate axis: C2; left axis: C0, C1, and C3). (b) Total protein concentration profiles.
(a)
(b)


## 4.1. Time-Varying Parameters and Input Disturbances

As in [18, 45], the following expressions are adopted for simulating uncertainties in the parameters and the feed flow rate:(30)Y^S/X=7.3+0.05sin⁡t,S^F=20+0.05sin⁡t,u^=u+0.05sin⁡t.The tracking error (en) is defined as(31)en2=Pref,n-Pn2+Tref,n-Tn2+Xref,n-Xn2+Sez,n-Sn2.In what follows, we shall call C0, C1, and C2 the controllers that use p = 0, p = 1, and p = 2 integrators, respectively, and C3 the 2-DOF PID controller. Particularly note that C0 corresponds to the controller designed in [30]. In spite of the assumed time-varying parameters, Figure 2(a) shows that C2 produces the minimum tracking error. According to Figure 2(b), both C2 and C1 are effective to reach the desired protein concentration profile, while C0 and C3 produce large errors.Figure 2
(a) Comparison of the tracking errors (see (31)) (left ordinate axis: C0 and C3; right axis: C2 and C1). (b) Total protein concentration profiles: C1 and C2 successfully reach the maximum production of protein in the presence of time-varying parameters.
(a)
(b)In several research fields, probabilistic methods have been found to be useful for dealing with problems related to robustness of systems affected by uncertainties [46]. Particularly, the Monte Carlo Randomized Algorithm has been used in many applications such as radioactive decay, power system generation, and traffic on roads [47]. In the control area, Monte Carlo methods are effective tools for the analysis of probabilistically robust control schemes [30, 46]. A Monte Carlo simulation is here performed to demonstrate the effectiveness of C2 from a statistical point of view [48–50], when uncertainties are introduced by time-varying parameter. To this effect, 1,000 random values of SF and YS/X are chosen in a range of ±25% of their nominal values. Figures 3(a), 3(c), 3(e), and 3(g) show that only C2 and C1 were able to adequately track the reference profile. In fact, both controller responses are quite similar, thus suggesting that the time-varying parameter has not affected their performances. For comparison, the tracking error was evaluated through the following performance index J:(32)J=0.25Δt∑n=1NrefPref,n-Pn2+Tref,n-Tn2+Xref,n-Xn2+Sez,n-Sn2Nref,where Nref represents the number of points of the reference profile. Figures 3(b), 3(d), 3(f), and 3(h) confirm a better performance of C2 and C1 over C0 and C3. In fact, J is quite small for C2 and C1, even though 25% of uncertainty in the main bioreactor parameters was adopted. Since 1,000 simulations with random values of SF and YS/X were carried out and the tracking errors remained bounded, the performance of C2 and C1 will be acceptable with 99% of probability whatever the parameter values are in a range of ±25% [46].Figure 3
Total protein concentration ((a), (c), (e), and (g)) and performance indexJ ((b), (d), (f), and (h)) for 1,000 random values of SF and YS/X (controllers are indicated in each figure).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 4.2. Uncertain Nonlinearities

Uncertain nonlinearities are often present in many bioprocesses. For instance, the specific growth rate,µ, can suffer from important changes from batch to batch, and unfortunately such variability is impossible to be accurately modeled. There are a number of models in the literature for μ [8], for example, the third expression of (A.2), for the bioprocess currently under study [31]. To demonstrate the controller performance, an additional term can be included in the calculation of μ as follows [18, 51]:(33)μ^=21.87+0.28sin⁡tSS+0.4S+62.5.Figure 4(a) compares the performances reached by C2, C1, C0, and C3 under the simulated uncertain nonlinearities in μ. The lower tracking errors correspond to C2. Figure 4(b) shows that only C2 and C1 are able to reach the maximum concentration of the total protein.Figure 4
(a) Comparison of the tracking errors (see (31)) under uncertain nonlinearities in μ (left ordinate axis: C0, C1, and C3; right axis: C2). (b) Total protein concentration profiles.
(a)
(b)


## 4.3. Unmodeled Dynamics

The presence of white noise inµmax leads to a mathematical model governed by a set of stochastic differential equations. In order to properly deal with this type of uncertainties, a kind of stochastic control approach should be used. As an alternative to the complex theory of stochastic control, a simpler deterministic strategy is here proposed. To this effect, note first that all parameters in (A.1) and (A.2) are truly affected by environmental noise, and consequently they always fluctuate at long times around some average values. For simplicity, we will consider that only µmax is affected by randomness and changes according to the following expression [52, 53]:(34)dμmax=μmax+αdBt,where dB(t) is a standard Brownian motion defined on a complete probability space and α ≥ 0 represents the intensity of the noisy signal.A white noise can be approximated through a band-limited white noise (BLWN), which is simply implemented through a proper block of Simulink. Theoretically, a white noise is uncorrelated; it exhibits a flat power spectral density, and therefore it is characterized by an infinite energy. The BLWN block can simulate the effect of white noise by using a random sequence with a correlation time much smaller than the shortest time constant of the system. The correlation time of the noise is the sample rate of the block [54] and can be estimated as(35)tc≈2π1001fmax=2π10010.5/Δtmax,where tc is the correlation time, fmax is the bandwidth of the system, and Δtmax is the maximum sample time for which the system is stable. For the feedback control presented in this work, Δtmax = 0.0145 h. Then, from (35), tc = 0.0018 h.Many researchers have utilized this methodology to prove the robustness of the bioreactors controllers [19, 22–24, 55]. In the current simulations, a BLWN with a noise power of 0.9 and α = 0.5 is introduced along with the feedback signal. Thus, the maximum specific growth rate affected with BLWN is calculated as(36)μ~max=μmax+0.05BLWN.Figure 5(a) indicates that C2 and C1 adequately work in the presence of unmodeled dynamics. In fact, the tracking errors are quite small in comparison with C0 and C3. Figure 5(b) shows the profile tracking for each controller. Again, C2 and C1 are able to follow the reference profile, while C0 and C3 fail to reach the maximum protein concentration.Figure 5
(a) Comparison of the tracking errors (see (31)) (right ordinate axis: C2; left axis: C0, C1, and C3). (b) Total protein concentration profiles.
(a)
(b)


## 5. Conclusions

A new control law for tracking an optimal concentration profile in a fed-batch bioreactor has been presented. The proposed method allows controlling a nonlinear system in presence of unavoidable uncertainties in bioreactors: time-varying parameters, uncertain nonlinearities, and unmodeled dynamics.In order to deal with these uncertainties, some integrators have been added to the original controller proposed in [30], and a Radial Basis Function (RBF) neural network was employed to estimate X, P, and T. The integrators have been chosen based on the variation hypothesis of the modeling error (En): if a q-grade polynomial can approximate En, then q+1 integrators will be added in order to avoid its influence on the tracking error (en).Several tests were carried out to demonstrate the effectiveness of the proposed methodology. The controllers are tuned through a Monte Carlo experiment and under stability conditions. To evaluate the performance of the designed controllers, different uncertainties were simulated. It was proven that the proposed methodology successfully compensates the uncertainties, and the tracking errors decrease as long as the number of integrators increases. The performance of the proposed controller with two integrators (C2) is remarkable, and the complexity of the control algorithm is relatively small. Moreover, the addition of integrators does not increase significantly the calculation time.In comparison to previous strategies that deal with similar control problems (e.g., [17, 21, 22, 24, 25, 27]), the proposed controllers present the advantage of avoiding a stochastic modeling, as it is typically required to deal with systems where the parameters are contaminated by white noise. The trade-off between the response speed and the stability is solved through a simple Monte Carlo experiment. In addition, the controller design is mathematically simple. Currently, the developed methodology for the controller design is being applied to another bioprocess. Besides, the inclusion of saturation of the control signals and sliding mode control schemes in the formulation of the problem will be addressed in future contributions.



---

*Source: 9391879-2017-09-05.xml*