# Nonlinear Control for Bioprocesses with Model Uncertainties and External Disturbances

Authors: Cecilia Fernández; Nadia Pantano; Leandro Rodriguez; Gustavo Scaglia
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a new alternative for profiles tracking control considering additive uncertainties is proposed. Based on a previously presented work about a nonlinear and multivariable controller design for a fed-batch bioethanol production, parametric uncertainty and process disturbance are taken into account to find a more reliable control strategy for a successful industrial implementation. To decrease the uncertainties effect, an approach based on the error estimation using Newton’s backward interpolation is included in the design equations. The proposed modification assures the error convergence to zero (demonstration is shown) despite the uncertainties, which is one of the main contributions of this work. A comparison between the new, the original proposal, and another methodology is exposed.


---

## Body

## 1. Introduction

Biological processes are becoming more frequent nowadays due to the wide variety of products obtained from them [1] and their possibility of making some processes environmentally friendly while high standard products are obtained [2, 3]. The fed-batch operation mode is which most interest awake for its main benefits [4, 5]. One of the most interesting advantages is the substrate concentration regulation in the cultivation medium by a suitable feed rate profile [6], obtaining better production yields and minimizing the production costs [5, 7]. However, bioprocesses control is required to follow a certain feed flow rate and get stability and the best productivity [8]. Furthermore, the mathematical representation of the process is the key to achieve good results.A mathematical model provides a map from inputs to responses. The quality of a model depends on how closely its responses match those of the true plant. However, a model set that includes the true physical plant can never be constructed [9]. Generally, a bioprocess modeling presents particular difficulty in their parameters determination caused by the poorly understood microorganism’s dynamics (multivariable and highly nonlinear dynamics), the strongly coupled variables, and the presence of numerous external disturbances, which leads to having many modeling uncertainties [10]. Furthermore, sometimes those parameters are determined without a previous model analysis, or their values are not informed with their respective confidence intervals [11, 12]. Moreover, time-varying parameters are usually assumed as constants [13]. Also, the uncertainties related to the processing technologies parameters are rarely considered [14]. The dismissal of all these aspects leads to a poor real-life representation and, consequently, to a bad performance with severe risks [15, 16]. Therefore, the main task to guarantee the bioprocess quality implies finding a way to control these distortions [17–19]. For this reason, the development of new control schemes that reduce the effect of uncertainties in the tracking error has become an attention focus for the scientific community [11, 15].One of the problems to take into account is the feasibility of the occurrence of events in systems with uncertainty, for this, the uncertainty theory was introduced [20]. On the other hand, a very interesting way to describe uncertainty is using the uncertain fractional differential equation (UFDE), which allows keeping a record of some properties to be considered a posteriori, which is fundamental in the evolution of uncertainties [21, 22]. Furthermore, the contemplation of uncertainties is approached with different strategies, depending on the nature of the system under study [22–24]. One of the most used strategies for model parameter identification and/or estimation involves an offline optimization using a nominal model of the process [25–30]; the main disadvantage of this methodology is that the variability of microorganisms decreases the possibility of batch-to-batch repeatability. To improve the results of nominal optimization, a methodology called “run to run” optimization appears, which uses previous runs information to optimize the operation of subsequent ones [31–37]. Another strategy is the online optimization of the model parameters [38–45]. This kind of optimization is difficult to perform since the available models might only be locally valid and thus inappropriate for predicting final concentrations [46].On the other hand, many scientists have developed some feedback control strategies to deal with bioprocess uncertainties. Optimal control, nonlinear model predictive control, hybrid control, adaptive control, fuzzy control, tracking control, and neural network are examples of them [47, 48]. But, due to the online implementation difficulty, the high computational cost, the imprecise mathematical models, and online solutions, their use for bioprocesses is limited [7] and is still a research topic [49].For the specific case of ethanol production defined by Hunag et al. [36], Fernández et al. [50] presented a controller design focused on looking for control actions, to track predefined concentration profiles. As the controller structure comes from the mathematical model of the process, it can be implemented in many systems. The procedure is characterized by its simplicity, versatility, and precision. Besides, only basic knowledge of numerical methods and linear algebra is needed to implement it. One of the main contributions of that work was to achieve the tracking error convergence to zero. Also, the technique was tested against different disturbances and compared with a typical PID controller.This manuscript aims to improve the control strategy presented in [50]. In this sense, predefined trajectories can be tracked while estimating the difference between the model and the real plant in each sampling time (error estimation). To reach this objective, another term is incorporated into the controller structure, which symbolizes model uncertainties and outside instabilities that are evaded with Newton’s backward interpolation. Moreover, the higher the interpolation order, the better the estimation and the smaller the tracking error. Besides, this change guarantees uniformity in the signal and progressive reduction of tracking error, achieving improvements of up to 98% in some cases.In [51], the problem of optimal profiles tracking control under uncertainties for a fed-batch bioprocess with two control actions is addressed with excellent results. In that manuscript, the authors add tracking error integrators in the control action calculation to reduce the additive uncertainties effect. However, the strategy presented in this paper has the advantage that adding the additive uncertainty term does not increase the order of the system, making mathematical development even simpler. Moreover, the ethanol bioprocess is an underactuated system with only one control action, so the control challenge is even greater than in [51]. Thus, a solution to the real trouble of multivariable and nonlinear tracking control in the presence of additive uncertainties is proposed, without increasing the system order.The manuscript is presented as follows: first, a summary of the process under study and the original control technique is described to contextualize the problem. Second, the contribution of this work is detailed, including the relevant demonstrations. Third, algorithms are tested and compared. Lastly, conclusions are shown.


## 2. Process and Control Description

### 2.1. Mathematical Model of the Process

Hunag et al. [36] proposed a mathematical model for a fed-batch ethanol production, using Saccharomyces diastaticus yeast to carry out the fermentation. The temperature was fixed at 35.8°C, airflow at 1.5vvm, and pH at 5.0. The only system input is the feed rate (U), which is a 50% glucose and 50% fructose combination. The state variables are biomass (X), ethanol (P1), glycerol (P2), glucose (S1), and fructose (S2) concentration inside the bioreactor:(1)X˙t=µ1+µ2X−UVX,S˙1t=−qS1/P1YP1/S1+qS1/P2YP2/S1X+UVλSf−S1,S˙2t=−qS2/P1YP1/S2+qS2/P2YP2/S2X+UV1−λSf−S2,P˙1t=qS1/P1+qS2/P1X−UVP1,P˙2t=qS1/P2+qS2/P2X−UVP2,where(2)V˙t=U,qS1/P1=νS1P1S1KS1P1+S1kS1P1kS1P1+P1,μ1=μm1S1KS1+S1+S12/KS1IKP1KP1+P1+P12/KP1I,qS1/P2=νS2P2S2KS2P2+S2kS2P2kS2P2+P2,μ2=μm2S2KS2+S2+S22/KS2IKP2KP2+P1+P12/KP2I,qS1/P2=νS1P2S1KS1P2+S1kS1P2kS1P2+P2,qS2/P1=νS2P1S2KS2P1+S2kS2P1kS2P1+P1.V is the volume of culture medium; µ1 and µ2 are the specific yeast growth rate, qS1/P1 and qS2/P1 are the specific ethanol production rate, and qS1/P2 and qS2/P2 are the specific glycerol production rate, in all cases from glucose and fructose, respectively. U is the control action. Initial conditions for ethanol fermentation are shown in Table 1, while the nomenclature, description, and values of parameters are in Table 2.Table 1
Initial conditions for ethanol fermentation.
VariableInitial valueX (g/L)1.5P1 (g/L)5.3P2 (g/L)0.0001S1 (g/L)8.6S2 (g/L)8.6V (L)1.35Table 2
Nomenclature, description, and values of parameters.
ParameterDescriptionValueµm1Maximum specific growth rate coefﬁcient for yeast on glucose (h−1)1.8823µm2Maximum specific growth rate coefﬁcient for yeast on fructose (h−1)1.7098YP1/S1Yield coefﬁcient for ethanol from glucose0.5085YP2/S1Yield coefﬁcient for glycerol from glucose0.5331YP1/S2Yield coefﬁcient for ethanol from fructose0.5098YP2/S2Yield coefﬁcient for glycerol from fructose0.4462KS1Saturation coefﬁcient for cell growth on glucose (g/L)159.75KS1IInhibition coefﬁcient for cell growth on glucose (g/L)94.233KP1Saturation coefﬁcient for cell growth on ethanol (g/L)238.39KP1IInhibition coefﬁcient for cell growth on ethanol (g/l)2.7378KS2Saturation coefﬁcient for cell growth on fructose (g/L)0.0726KS2IInhibition coefﬁcient for cell growth on fructose (g/L)9.0048KP2Saturation coefﬁcient for cell growth on glycerol (g/L)35.958KP2IInhibition coefﬁcient for cell growth on glycerol (g/L)9.9722KS1P1Saturation coefﬁcient for ethanol production on glucose (g/L)1.3409kS1P1Inhibition coefﬁcient for ethanol production on glucose (g/L)18.612KS2P1Saturation coefﬁcient for ethanol production on fructose (g/L)0.9129kS2P1Inhibition coefﬁcient for ethanol production on fructose (g/L)1000KS1P2Saturation coefﬁcient for glycerol production on glucose (g/L)6.7116kS1P2Inhibition coefﬁcient for glycerol production on glucose (g/L)0.5863KS2P2Saturation coefﬁcient for glycerol production on fructose (g/L)0.4310kS2P2Inhibition coefﬁcient for glycerol production on fructose (g/L)1.150νS1P1Coefﬁcient of maximum speciﬁc ethanol production rate for yeast on glucose (h−1)1.5051νS2P1Coefﬁcient of maximum speciﬁc ethanol production rate for yeast on fructose (h−1)0.3321νS1P2Coefﬁcient of maximum speciﬁc glycerol production rate for yeast on glucose (h−1)0.0023νS2P2Coefﬁcient of maximum speciﬁc glycerol production rate for yeast on fructose (h−1)0.1609λProportion of glucose and fructose0.5SfSugar total feed concentration (g/L)300


### 2.2. Controller Structure Design

In [50], a technique that finds U to make the system track predefined profiles (references) is proposed. For this methodology, design is supposed that both the references and the states are known at each sampling instant. This last assertion is far from reality, so system states were estimated with neural network state estimators previously designed and published [52]. In Figure 1, reference concentration profiles and U are shown. The following is a brief description of the technique described in [50]:Figure 1
Cells, ethanol, glycerol, and fructose reference concentrations along the process. Reference feed flow rate.Firstly, the system equation (1) is integrated using numerical methods. Euler is used for its simplicity:(3)dσdt=σn+1−σnTS.In (3), σ symbolizes states variables, σn is the current value of σ measured online, and σn+1 is the σ value in the next measurement instant. TS is the sampling time (0.1 h) [53]. The process total time is 15.7 h (Tf).Then,σn+1 are approached with(4)σref,n+1−σn+1︸errorn+1=kσσref,n−σn︸errorn⟶σn+1=σref,n+1−kσσref,n−σn,where σref are the reference state variables and kσ is the controller parameter for the variable σ. For this process, the controller parameters are kX, kP1, kP2, kS1, and kS2. Then, substituting (4) in (3),(5)dσdt=σref,n+1−kσσref,n−σn︷σn+1−σnTS=Δσn.Replacing (5) in (1),(6)ΔXn=µ1+µ2Xn−UnVnXn,ΔS1,n=−qS1/P1YP1/S1+qS1/P2YP2/S1Xn+UnVnλSf−S1,n,ΔS2,n=−qS2/P1YP1/S2+qS2/P2YP2/S2Xn+UnVn1−λSf−S2,n,ΔP1,n=qS1/P1+qS2/P1Xn−UnVnP1,n,ΔP2,n=qS1/P2+qS2/P2Xn−UnVnP2,n.Stating (6) as a matrix,(7)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/Vn︷AnUn=ΔXn−µ1+µ2XnΔS1,n+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2,n+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1,n−qS1/P1+qS2/P1XnΔP2,n−qS1/P2+qS2/P2Xn︸bn.System equation (7) must have an exact solution to find Un. Consequently, bn have to be a linear combination of An [54]; that is to say, An and bn must be parallel and with the same sense. One way to accomplish this is(8)cosAn,bn=An,bnAn∗bn|=1.In (8), the operation between < > and ||.|| represent the inner product and the vectors norm in Rn space, respectively. The angle between A and b is θ = 0°; this implies a positive Un value. If θ = 180°, Un would be negative, which does not make physical sense because Un is a flow.As stated in the Introduction, bioprocess systems generally involve several control objectives that may be conflicting, to balance these objectives during the design of the presented controller, the selection of a “sacrificed variable” is required. It is denoted asS1ez, which guarantees that (7) has an exact solution and that the references are followed. For more details on its selection and calculation, see [50]. Finally, Un is obtained using least squares [54]:(9)Un=AnTAn−1AnTbn.


### 2.3. Controller Tuning

In [50], the Monte Carlo algorithm was used to find the controller parameters that make the accumulated error minimum. It consists of simulating the bioprocess N times with random kσ [55]:(10)N≥log1/δlog1/1−ε.In (10), δ is the confidence and ε is the accuracy.However, many other strategies can be used to tune this controller. In [56], a genetic algorithm and a hybrid one are proposed.Next, the “tracking error (||en||)” and “total error (Ep)” concepts are introduced:(11)en=Xref,n−XnmaxXref,n2+P1ref,n−P1,nmaxP1ref,n2+P2ref,n−P2,nmaxP2ref,n2+S2ref,n−S2,nmaxS2ref,n2,(12)Ep,1=Ts∑n=1Jen,(13)Ep,2=Ts∑n=1JennTs,where p represents the simulation in progress, p = 1, 2, ..., N; subscripts 1 and 2 differentiate between one index and another; n is the sample instant, n = 1, 2, ..., J; Tf = J TS.Inkσ selection, (12) is the function cost to be minimized.Theorem 1.
If the discrete system is given by equation (1), the control action is calculated with equation (9), and kσ take values between zero and one (0 < kσ< 1); then, the tracking error convergence to zero when n tends to infinity is achieved.Demonstration [50]:Substituting the sacrificed variable in (7) and expressing the matrix system generically,(14)a1a2a3a4a5︸AUn=b1b2b3b4b5︸b.Appling least squares to (14), (15)Un=ATA−1ATb=a1b1+a2b2+a3b3+a4b4+a5b5a12+a22+a32+a42+a52.From (14),(16)a1a2=b1b2⟶b2=a2a1b1,a1a3=b1b3⟶b3=a3a1b1,a1a4=b1b4⟶b4=a4a1b1,a1a5=b1b5⟶b5=a5a1b1.Placing (16) in (15),(17)Un=a1b1+a22b1/a1+a32b1/a1+a42b1/a1+a52b1/a1a12+a22+a32+a42+a52=b1/a1a12+a22+a32+a42+a52a12+a22+a32+a42+a52=b1a1=VnXref,n+1−kXXref,n−Xn−Xn/TS−μ1S1ez,n,P1,n+μ2S2,n,P1,nXn−Xn.Replacing (17) in (7),(18)Xn+1=Xref,n+1−kXXref,n−Xn+TSμ1S1ez,n,P1,n−μ1S1,n,P1,nXn.Before, the tracking error forX is defined as(19)eX,n+1=Xref,n+1−Xn+1.Substituting (18) in (19),(20)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n−μ1S1,n,P1,nX.nTheµ1(S1n, P1n) Taylor approximation in the desired value µ1(S1ez n, P1n) is(21)μ1S1,n,P1,n=μ1S1ez,n,P1,n+dμ1S1,P1,ndS1|S1n−S1ez,nS1=Sez,n+θS1,n−S1ez,n=Sθ,where⟶0<θ<1Placing (21) in (20),(22)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n+dμ1S1,P1,ndS1|SθS1n−S1ez,n︸−eS1,n−μ1S1ez,n,P1,nXn,eX,n+1=kXXref,n−Xn+TSdμ1S1,P1,ndS1|SθeS1,nXn.Following the same procedure for the other variables and joining the final expressions,(23)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|Sα︸NLeS1,n.In equation (23), L is a linear system and NL is a bounded nonlinearity [50]. Note that if kσ = 0, the reference is reached in one step. Thus, if 0 < kσ< 1, the tracking error tends to zero when n ⟶ ∞ [50, 57].


### 2.4. Steps to Implement the Controller

Step 1.
DefineTs, σref, and σn.Step 2.
Discretize differential equations using some numerical methods, equation (3).Step 3.
Obtain the state variables inn+1 with equation (4).Step 4.
Define and calculate the sacrificed variable.Step 5.
DetermineUn with least squares, equation (9).
Figure2 outlines the control diagram.Figure 2
Control system diagram.


## 2.1. Mathematical Model of the Process

Hunag et al. [36] proposed a mathematical model for a fed-batch ethanol production, using Saccharomyces diastaticus yeast to carry out the fermentation. The temperature was fixed at 35.8°C, airflow at 1.5vvm, and pH at 5.0. The only system input is the feed rate (U), which is a 50% glucose and 50% fructose combination. The state variables are biomass (X), ethanol (P1), glycerol (P2), glucose (S1), and fructose (S2) concentration inside the bioreactor:(1)X˙t=µ1+µ2X−UVX,S˙1t=−qS1/P1YP1/S1+qS1/P2YP2/S1X+UVλSf−S1,S˙2t=−qS2/P1YP1/S2+qS2/P2YP2/S2X+UV1−λSf−S2,P˙1t=qS1/P1+qS2/P1X−UVP1,P˙2t=qS1/P2+qS2/P2X−UVP2,where(2)V˙t=U,qS1/P1=νS1P1S1KS1P1+S1kS1P1kS1P1+P1,μ1=μm1S1KS1+S1+S12/KS1IKP1KP1+P1+P12/KP1I,qS1/P2=νS2P2S2KS2P2+S2kS2P2kS2P2+P2,μ2=μm2S2KS2+S2+S22/KS2IKP2KP2+P1+P12/KP2I,qS1/P2=νS1P2S1KS1P2+S1kS1P2kS1P2+P2,qS2/P1=νS2P1S2KS2P1+S2kS2P1kS2P1+P1.V is the volume of culture medium; µ1 and µ2 are the specific yeast growth rate, qS1/P1 and qS2/P1 are the specific ethanol production rate, and qS1/P2 and qS2/P2 are the specific glycerol production rate, in all cases from glucose and fructose, respectively. U is the control action. Initial conditions for ethanol fermentation are shown in Table 1, while the nomenclature, description, and values of parameters are in Table 2.Table 1
Initial conditions for ethanol fermentation.
VariableInitial valueX (g/L)1.5P1 (g/L)5.3P2 (g/L)0.0001S1 (g/L)8.6S2 (g/L)8.6V (L)1.35Table 2
Nomenclature, description, and values of parameters.
ParameterDescriptionValueµm1Maximum specific growth rate coefﬁcient for yeast on glucose (h−1)1.8823µm2Maximum specific growth rate coefﬁcient for yeast on fructose (h−1)1.7098YP1/S1Yield coefﬁcient for ethanol from glucose0.5085YP2/S1Yield coefﬁcient for glycerol from glucose0.5331YP1/S2Yield coefﬁcient for ethanol from fructose0.5098YP2/S2Yield coefﬁcient for glycerol from fructose0.4462KS1Saturation coefﬁcient for cell growth on glucose (g/L)159.75KS1IInhibition coefﬁcient for cell growth on glucose (g/L)94.233KP1Saturation coefﬁcient for cell growth on ethanol (g/L)238.39KP1IInhibition coefﬁcient for cell growth on ethanol (g/l)2.7378KS2Saturation coefﬁcient for cell growth on fructose (g/L)0.0726KS2IInhibition coefﬁcient for cell growth on fructose (g/L)9.0048KP2Saturation coefﬁcient for cell growth on glycerol (g/L)35.958KP2IInhibition coefﬁcient for cell growth on glycerol (g/L)9.9722KS1P1Saturation coefﬁcient for ethanol production on glucose (g/L)1.3409kS1P1Inhibition coefﬁcient for ethanol production on glucose (g/L)18.612KS2P1Saturation coefﬁcient for ethanol production on fructose (g/L)0.9129kS2P1Inhibition coefﬁcient for ethanol production on fructose (g/L)1000KS1P2Saturation coefﬁcient for glycerol production on glucose (g/L)6.7116kS1P2Inhibition coefﬁcient for glycerol production on glucose (g/L)0.5863KS2P2Saturation coefﬁcient for glycerol production on fructose (g/L)0.4310kS2P2Inhibition coefﬁcient for glycerol production on fructose (g/L)1.150νS1P1Coefﬁcient of maximum speciﬁc ethanol production rate for yeast on glucose (h−1)1.5051νS2P1Coefﬁcient of maximum speciﬁc ethanol production rate for yeast on fructose (h−1)0.3321νS1P2Coefﬁcient of maximum speciﬁc glycerol production rate for yeast on glucose (h−1)0.0023νS2P2Coefﬁcient of maximum speciﬁc glycerol production rate for yeast on fructose (h−1)0.1609λProportion of glucose and fructose0.5SfSugar total feed concentration (g/L)300


## 2.2. Controller Structure Design

In [50], a technique that finds U to make the system track predefined profiles (references) is proposed. For this methodology, design is supposed that both the references and the states are known at each sampling instant. This last assertion is far from reality, so system states were estimated with neural network state estimators previously designed and published [52]. In Figure 1, reference concentration profiles and U are shown. The following is a brief description of the technique described in [50]:Figure 1
Cells, ethanol, glycerol, and fructose reference concentrations along the process. Reference feed flow rate.Firstly, the system equation (1) is integrated using numerical methods. Euler is used for its simplicity:(3)dσdt=σn+1−σnTS.In (3), σ symbolizes states variables, σn is the current value of σ measured online, and σn+1 is the σ value in the next measurement instant. TS is the sampling time (0.1 h) [53]. The process total time is 15.7 h (Tf).Then,σn+1 are approached with(4)σref,n+1−σn+1︸errorn+1=kσσref,n−σn︸errorn⟶σn+1=σref,n+1−kσσref,n−σn,where σref are the reference state variables and kσ is the controller parameter for the variable σ. For this process, the controller parameters are kX, kP1, kP2, kS1, and kS2. Then, substituting (4) in (3),(5)dσdt=σref,n+1−kσσref,n−σn︷σn+1−σnTS=Δσn.Replacing (5) in (1),(6)ΔXn=µ1+µ2Xn−UnVnXn,ΔS1,n=−qS1/P1YP1/S1+qS1/P2YP2/S1Xn+UnVnλSf−S1,n,ΔS2,n=−qS2/P1YP1/S2+qS2/P2YP2/S2Xn+UnVn1−λSf−S2,n,ΔP1,n=qS1/P1+qS2/P1Xn−UnVnP1,n,ΔP2,n=qS1/P2+qS2/P2Xn−UnVnP2,n.Stating (6) as a matrix,(7)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/Vn︷AnUn=ΔXn−µ1+µ2XnΔS1,n+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2,n+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1,n−qS1/P1+qS2/P1XnΔP2,n−qS1/P2+qS2/P2Xn︸bn.System equation (7) must have an exact solution to find Un. Consequently, bn have to be a linear combination of An [54]; that is to say, An and bn must be parallel and with the same sense. One way to accomplish this is(8)cosAn,bn=An,bnAn∗bn|=1.In (8), the operation between < > and ||.|| represent the inner product and the vectors norm in Rn space, respectively. The angle between A and b is θ = 0°; this implies a positive Un value. If θ = 180°, Un would be negative, which does not make physical sense because Un is a flow.As stated in the Introduction, bioprocess systems generally involve several control objectives that may be conflicting, to balance these objectives during the design of the presented controller, the selection of a “sacrificed variable” is required. It is denoted asS1ez, which guarantees that (7) has an exact solution and that the references are followed. For more details on its selection and calculation, see [50]. Finally, Un is obtained using least squares [54]:(9)Un=AnTAn−1AnTbn.


## 2.3. Controller Tuning

In [50], the Monte Carlo algorithm was used to find the controller parameters that make the accumulated error minimum. It consists of simulating the bioprocess N times with random kσ [55]:(10)N≥log1/δlog1/1−ε.In (10), δ is the confidence and ε is the accuracy.However, many other strategies can be used to tune this controller. In [56], a genetic algorithm and a hybrid one are proposed.Next, the “tracking error (||en||)” and “total error (Ep)” concepts are introduced:(11)en=Xref,n−XnmaxXref,n2+P1ref,n−P1,nmaxP1ref,n2+P2ref,n−P2,nmaxP2ref,n2+S2ref,n−S2,nmaxS2ref,n2,(12)Ep,1=Ts∑n=1Jen,(13)Ep,2=Ts∑n=1JennTs,where p represents the simulation in progress, p = 1, 2, ..., N; subscripts 1 and 2 differentiate between one index and another; n is the sample instant, n = 1, 2, ..., J; Tf = J TS.Inkσ selection, (12) is the function cost to be minimized.Theorem 1.
If the discrete system is given by equation (1), the control action is calculated with equation (9), and kσ take values between zero and one (0 < kσ< 1); then, the tracking error convergence to zero when n tends to infinity is achieved.Demonstration [50]:Substituting the sacrificed variable in (7) and expressing the matrix system generically,(14)a1a2a3a4a5︸AUn=b1b2b3b4b5︸b.Appling least squares to (14), (15)Un=ATA−1ATb=a1b1+a2b2+a3b3+a4b4+a5b5a12+a22+a32+a42+a52.From (14),(16)a1a2=b1b2⟶b2=a2a1b1,a1a3=b1b3⟶b3=a3a1b1,a1a4=b1b4⟶b4=a4a1b1,a1a5=b1b5⟶b5=a5a1b1.Placing (16) in (15),(17)Un=a1b1+a22b1/a1+a32b1/a1+a42b1/a1+a52b1/a1a12+a22+a32+a42+a52=b1/a1a12+a22+a32+a42+a52a12+a22+a32+a42+a52=b1a1=VnXref,n+1−kXXref,n−Xn−Xn/TS−μ1S1ez,n,P1,n+μ2S2,n,P1,nXn−Xn.Replacing (17) in (7),(18)Xn+1=Xref,n+1−kXXref,n−Xn+TSμ1S1ez,n,P1,n−μ1S1,n,P1,nXn.Before, the tracking error forX is defined as(19)eX,n+1=Xref,n+1−Xn+1.Substituting (18) in (19),(20)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n−μ1S1,n,P1,nX.nTheµ1(S1n, P1n) Taylor approximation in the desired value µ1(S1ez n, P1n) is(21)μ1S1,n,P1,n=μ1S1ez,n,P1,n+dμ1S1,P1,ndS1|S1n−S1ez,nS1=Sez,n+θS1,n−S1ez,n=Sθ,where⟶0<θ<1Placing (21) in (20),(22)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n+dμ1S1,P1,ndS1|SθS1n−S1ez,n︸−eS1,n−μ1S1ez,n,P1,nXn,eX,n+1=kXXref,n−Xn+TSdμ1S1,P1,ndS1|SθeS1,nXn.Following the same procedure for the other variables and joining the final expressions,(23)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|Sα︸NLeS1,n.In equation (23), L is a linear system and NL is a bounded nonlinearity [50]. Note that if kσ = 0, the reference is reached in one step. Thus, if 0 < kσ< 1, the tracking error tends to zero when n ⟶ ∞ [50, 57].


## 2.4. Steps to Implement the Controller

Step 1.
DefineTs, σref, and σn.Step 2.
Discretize differential equations using some numerical methods, equation (3).Step 3.
Obtain the state variables inn+1 with equation (4).Step 4.
Define and calculate the sacrificed variable.Step 5.
DetermineUn with least squares, equation (9).
Figure2 outlines the control diagram.Figure 2
Control system diagram.


## 3. Control Structure Design under Additive Uncertainties

### 3.1. Controller Design under Uncertainty

Next, uncertainties effect in the tracking error are considered by adding the termsEσ,n in (7):(24)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/VnUn=ΔXn−µ1+µ2XnΔS1,n+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2,n+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1,n−qS1/P1+qS2/P1XnΔP2,n−qS1/P2+qS2/P2Xn−1TSEX,nES1,nES2,nEP1,nEP2,n︸En.Additive uncertainty (En) can be used to model several kinds of uncertainties as well as external perturbations (measurement errors are not considered). It might depend on the state variables and the system input. Moreover, considering a real plant like zn+1 = f (zn,un), therefore the additive uncertainty can be expressed as En = f (zn,un) –f^(zn,un), where f^ (zn,un) is the discrete-time nonlinear system model. If z and u are assumed to be bounded and f is Lipschitz [58], then Eσ,n can be modeled as a bounded uncertainty [59, 60].The uncertainties terms in (24) affect the error convergence to zero of the tracking error [50]. This can be observed following the same procedure as in Section 2.3 Theorem demonstration:(25)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|SαeS1,n︸NL−EX,nES1,nES2,nEP1,nEP2,n︸En.Comparing (25) with (23), the error nonconvergence to zero is noticed due to En presence.Therefore, the next step isEσ,n estimation to reduce their effect on the tracking error, achieving the convergence to zero.


### 3.2. Uncertainty Estimation

The following procedure, as the main contribution of this work, develops a control strategy based on linear algebra that reduces the effect of uncertainty in tracking error by its estimation using Newton’s backward interpolation [61]. The advantage is the estimation development with an easy-to-understand numerical procedure, which does not add further complexity to the control methodology previously presented. Thus, En is estimated with Ên. Then, Ên is added in the control action calculation:(26)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/VnnUn=ΔXn−µ1+µ2XnΔS1,n+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2,n+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1,n−qS1/P1+qS2/P1XnΔP2,n−qS1/P2+qS2/P2Xn−1TSE^X,nE^S1,nE^S2,nE^P1,nE^P2,n.Taking into account thatEσ,n is unknown, but it is assumed as a polynomial, its differences can be defined as(27)δEσ,n=Eσ,n+1−Eσ,n,δ2Eσ,n=δδEσ,n=δEσ,n+1−Eσ,n=Eσ,n+2−2Eσ,n+1+Eσ,n,δqEσ,n=δδq−1Eσ,n.


#### 3.2.1. Constant Uncertainty

IfEσ,n is assumed as a constant, δEσ,n = 0. The uncertainty estimation is represented by(28)Eσ,n=σn+1−σ^n+1.where σ is the real state variable and σ^ represents the state variable calculated with equation (3). Then, the uncertainty is estimated with a zero-order approximation:(29)E^σ,n=Eσ,n.Demonstration of error convergence to zero when uncertainties are considered:Replacing the control action given by equation (26) in (24) and applying the same steps described in Section 2.3 demonstration,(30)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|SαeS1,n︸NL+EX,nES1,nES2,nEP1,nEP2,n︸En−E^X,nE^S1,nE^S2,nE^P1,nE^P2,n︸E^n.Remark 1.
If Eσ,n = c, where c is a constant value, and Êσ,n is calculated with equation (29), then δEσ,n = 0.
Thus, ifEn is unknown and each component is an m-order polynomial, En = Ên, then En and Ên can be simplified from (30). Therefore, it is demonstrated that the tracking error tends to zero with the process progress.


#### 3.2.2. Linear Uncertainty

IfEσ,n is supposed as a linear function, then δ2Eσ,n = 0. Following the same reasoning as in 3.2.1, Êσ,n is defined as(31)E^σ,n=Eσ,n+Eσ,n−Eσ,n−1.Following the same steps of3.2.1, the error convergence to zero is demonstrated.


#### 3.2.3. Polynomial Uncertainty

In this case,Eσ,n is supposed as an m-order polynomial function. Then, if q > m, δqEσ,n = 0. Using the same procedure explained before, the representation of Êσ,n is generically expressed as(32)E^σ,n=∑j=0m∑i=0jji−1iEσ,n−i−1j!.


### 3.3. Controller Parameter Selection

To select the bestkσ values (between zero and one for stability guarantee), as many authors recommend [62–64], the Monte Carlo algorithm is used following the procedure described in 2.3. Therefore, 1000 simulations were done to find the parameters for three different controllers, the original, described in [50] (C1), other with a zero-order estimator (C2), and the last with a first-order estimator (C3). Figure 3 shows the difference between the three controllers. Note the error improvement evidenced with both indexes ((12) and (13)), with C2, the total error Ep, 1 decreases by 22.22% and Ep, 2 by 19.31% compared to C1, while with C3, the total error is reduced by 52.38% and 48.73%, respectively, concerning C1. Due to the results similarity with both indexes, from now on only Ep,1 will be used to evaluate the algorithms.Figure 3
Controllers comparison. Total error decrease: (a) Ep 1 (12); (b) Ep 2 (13).
(a)(b)


## 3.1. Controller Design under Uncertainty

Next, uncertainties effect in the tracking error are considered by adding the termsEσ,n in (7):(24)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/VnUn=ΔXn−µ1+µ2XnΔS1,n+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2,n+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1,n−qS1/P1+qS2/P1XnΔP2,n−qS1/P2+qS2/P2Xn−1TSEX,nES1,nES2,nEP1,nEP2,n︸En.Additive uncertainty (En) can be used to model several kinds of uncertainties as well as external perturbations (measurement errors are not considered). It might depend on the state variables and the system input. Moreover, considering a real plant like zn+1 = f (zn,un), therefore the additive uncertainty can be expressed as En = f (zn,un) –f^(zn,un), where f^ (zn,un) is the discrete-time nonlinear system model. If z and u are assumed to be bounded and f is Lipschitz [58], then Eσ,n can be modeled as a bounded uncertainty [59, 60].The uncertainties terms in (24) affect the error convergence to zero of the tracking error [50]. This can be observed following the same procedure as in Section 2.3 Theorem demonstration:(25)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|SαeS1,n︸NL−EX,nES1,nES2,nEP1,nEP2,n︸En.Comparing (25) with (23), the error nonconvergence to zero is noticed due to En presence.Therefore, the next step isEσ,n estimation to reduce their effect on the tracking error, achieving the convergence to zero.


## 3.2. Uncertainty Estimation

The following procedure, as the main contribution of this work, develops a control strategy based on linear algebra that reduces the effect of uncertainty in tracking error by its estimation using Newton’s backward interpolation [61]. The advantage is the estimation development with an easy-to-understand numerical procedure, which does not add further complexity to the control methodology previously presented. Thus, En is estimated with Ên. Then, Ên is added in the control action calculation:(26)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/VnnUn=ΔXn−µ1+µ2XnΔS1,n+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2,n+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1,n−qS1/P1+qS2/P1XnΔP2,n−qS1/P2+qS2/P2Xn−1TSE^X,nE^S1,nE^S2,nE^P1,nE^P2,n.Taking into account thatEσ,n is unknown, but it is assumed as a polynomial, its differences can be defined as(27)δEσ,n=Eσ,n+1−Eσ,n,δ2Eσ,n=δδEσ,n=δEσ,n+1−Eσ,n=Eσ,n+2−2Eσ,n+1+Eσ,n,δqEσ,n=δδq−1Eσ,n.


### 3.2.1. Constant Uncertainty

IfEσ,n is assumed as a constant, δEσ,n = 0. The uncertainty estimation is represented by(28)Eσ,n=σn+1−σ^n+1.where σ is the real state variable and σ^ represents the state variable calculated with equation (3). Then, the uncertainty is estimated with a zero-order approximation:(29)E^σ,n=Eσ,n.Demonstration of error convergence to zero when uncertainties are considered:Replacing the control action given by equation (26) in (24) and applying the same steps described in Section 2.3 demonstration,(30)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|SαeS1,n︸NL+EX,nES1,nES2,nEP1,nEP2,n︸En−E^X,nE^S1,nE^S2,nE^P1,nE^P2,n︸E^n.Remark 1.
If Eσ,n = c, where c is a constant value, and Êσ,n is calculated with equation (29), then δEσ,n = 0.
Thus, ifEn is unknown and each component is an m-order polynomial, En = Ên, then En and Ên can be simplified from (30). Therefore, it is demonstrated that the tracking error tends to zero with the process progress.


### 3.2.2. Linear Uncertainty

IfEσ,n is supposed as a linear function, then δ2Eσ,n = 0. Following the same reasoning as in 3.2.1, Êσ,n is defined as(31)E^σ,n=Eσ,n+Eσ,n−Eσ,n−1.Following the same steps of3.2.1, the error convergence to zero is demonstrated.


### 3.2.3. Polynomial Uncertainty

In this case,Eσ,n is supposed as an m-order polynomial function. Then, if q > m, δqEσ,n = 0. Using the same procedure explained before, the representation of Êσ,n is generically expressed as(32)E^σ,n=∑j=0m∑i=0jji−1iEσ,n−i−1j!.


## 3.2.1. Constant Uncertainty

IfEσ,n is assumed as a constant, δEσ,n = 0. The uncertainty estimation is represented by(28)Eσ,n=σn+1−σ^n+1.where σ is the real state variable and σ^ represents the state variable calculated with equation (3). Then, the uncertainty is estimated with a zero-order approximation:(29)E^σ,n=Eσ,n.Demonstration of error convergence to zero when uncertainties are considered:Replacing the control action given by equation (26) in (24) and applying the same steps described in Section 2.3 demonstration,(30)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1|Sθ00−dqS1/P1S1,P1,ndS1|Sφ−dqS1/P2S1,P2,ndS1|SαeS1,n︸NL+EX,nES1,nES2,nEP1,nEP2,n︸En−E^X,nE^S1,nE^S2,nE^P1,nE^P2,n︸E^n.Remark 1.
If Eσ,n = c, where c is a constant value, and Êσ,n is calculated with equation (29), then δEσ,n = 0.
Thus, ifEn is unknown and each component is an m-order polynomial, En = Ên, then En and Ên can be simplified from (30). Therefore, it is demonstrated that the tracking error tends to zero with the process progress.


## 3.2.2. Linear Uncertainty

IfEσ,n is supposed as a linear function, then δ2Eσ,n = 0. Following the same reasoning as in 3.2.1, Êσ,n is defined as(31)E^σ,n=Eσ,n+Eσ,n−Eσ,n−1.Following the same steps of3.2.1, the error convergence to zero is demonstrated.


## 3.2.3. Polynomial Uncertainty

In this case,Eσ,n is supposed as an m-order polynomial function. Then, if q > m, δqEσ,n = 0. Using the same procedure explained before, the representation of Êσ,n is generically expressed as(32)E^σ,n=∑j=0m∑i=0jji−1iEσ,n−i−1j!.


## 3.3. Controller Parameter Selection

To select the bestkσ values (between zero and one for stability guarantee), as many authors recommend [62–64], the Monte Carlo algorithm is used following the procedure described in 2.3. Therefore, 1000 simulations were done to find the parameters for three different controllers, the original, described in [50] (C1), other with a zero-order estimator (C2), and the last with a first-order estimator (C3). Figure 3 shows the difference between the three controllers. Note the error improvement evidenced with both indexes ((12) and (13)), with C2, the total error Ep, 1 decreases by 22.22% and Ep, 2 by 19.31% compared to C1, while with C3, the total error is reduced by 52.38% and 48.73%, respectively, concerning C1. Due to the results similarity with both indexes, from now on only Ep,1 will be used to evaluate the algorithms.Figure 3
Controllers comparison. Total error decrease: (a) Ep 1 (12); (b) Ep 2 (13).
(a)(b)


## 4. Results and Discussion

In the following section, two important tests are developed to demonstrate the proposed estimation effectiveness. First, the Monte Carlo Algorithm is applied to test the controller operation under parametric uncertainty. Second, two different perturbations in the control action are added. Furthermore, both above tests are performed simultaneously and are compared with the performance of two other controllers using another methodology [65]. In this section, the original controller is noted as C1, the controller with a zero-order estimation is C2, and the controller with a first-order estimator is C3.


### 4.1. Simulation under Parametric Uncertainties

In all bioprocess, model parameters may vary in an unpredictable way [66]. This can lead to structural instability in the system dynamical behavior [15]. Therefore, a strict and efficient control system is required to deal with this problem.In several research fields, probabilistic methods are useful for dealing with problems related to the robustness of systems affected by uncertainties [55]. Particularly, Monte Carlo Randomized Algorithm has been used for uncertainty quantification in many applications [67–69]. In this paper, the Monte Carlo Randomized Algorithm is applied using the procedure described in 2.3. The number of simulations, N = 1000, is obtained with equation (10), adopting a confidence (δ) of 0.01 and an accuracy (ε) of 0.005. Therefore, the following test demonstrates the technique success from a statistical point of view [70–72].In a simulation, a way to quantify uncertainties and perturbations is to specify the parameters real range of variation instead of using a constant value with greater error [73]. Hunag et al. [36] specified, in Table 1 of their work, the parameters confidence intervals for the ethanol process under study. For this test, in each simulation, all the system parameters are randomly changed by ±10% of their original range value (Table 1 of [36]). Then, the total error is calculated with equation (12). Table 3 shows the error range for each controller. Note how the total error range considerably decreases and maintains bounded to lower values when the estimation order increases. Figure 4 shows the total error for 1000 simulations to evaluate the three-controller performance under parametric uncertainty. As can be seen, the total tracking error range is visibly reduced.Table 3
Parametric uncertainty of ±10%.
ControllerTotal error rangeC10.5267–0.0098C20.2958–0.0056C30.0812–0.0049Figure 4
Total error for 1000 simulations under parametric uncertainty (±10%). (C1) Original controller; (C2) controller with a zero-order estimation; and (C3) controller with first-order estimator.Expressing the total error for each case into a percentage can easily quantify the controller improvement when estimators are incorporated. Figure5 suggests the better performance of the controllers C2 and C3 over C1. The worst case presented for C2 has a total error of 43.83% lower than that of the worst situation of C1, while the maximum error presented with C3 is 84.58% less than that of C1.Figure 5
Controllers comparison of maximum total error.Considering Theorem 1 of [55], as the tracking error (equation (11)) for 1000 simulations remains bounded, C2 and C3 controllers operation will be satisfactory with 99% of probability while the parameters vary within a ±10% range.


### 4.2. Simulation for Uncertainties Prevention

In the following test, controllers are tuned taking into account the worst possible situation presented in 4.1 forC1, to prevent it. Figure 6 shows the total error decreases when each controller is tuned. Note how even considering the most problematic situation of the system, the performance error can be reduced by 52.55% just by adding a zero-order estimator and by 56.15% with a first-order one.Figure 6
Total error decrease when the controllers are tuned considering the worst situation of C1 (total error:C1 = 0.1334, C2 = 0.0633, C3 = 0.0585).Then, the controllers using the chosen parameters were tested under normal operation conditions, demonstrating that their performance is not negatively affected. Figure7 compares the accumulated error for the three controllers. Moreover, total errors (C1 = 0.1416, C2 = 0.0062, C3 = 0.0031) contrast is shown.Figure 7
Errors comparison for the three controllers: (C1) original controller; (C2) controller with a zero-order estimation; (C3) controller with the first-order estimator.This test demonstrates another advantage of the proposed estimation technique. Table4 compares the errors of Figure 3 with those of Figure 7. In both cases, the controllers are tested under normal conditions; however, the difference lies in the tuning conditions considered. Analyzing Table 2, it can be observed that C2 and C3 controllers present similar results, which does not happen with C1. Therefore, it can be said that the operation of C2 and C3 will have a minimum error over the entire range of possible uncertainties between the nominal system parameters values and their worst variation.Table 4
Total error comparison under normal conditions.
C1C2C3Total errora)0.00630.00490.0030b)0.14160.00620.0031a) Tuned in normal conditions; b) tuned considering the worst situation of 4.1 for C1.


### 4.3. Simulation Adding Perturbations in the Control Action

In this test, a hypothetical situation that may produce an unexpected variation in the production is simulated. Firstly, a -30% step perturbation in the bioreactor feed rate is added to evaluate the response of controllers. Secondly, a ramp disturbance is added to the step perturbation previously presented. Figures8 and 9 show the control action variation compared to the reference and the percentage error, considering C1 error as 100%. Note how the total error improves when the new algorithm is applied. In the first test, it is reduced by 90.46% with C2 and by 94.99% with C3, while in the second test, the results are improved 91.87% with C2 and 97.15% with C3.Figure 8
(a) Controllers response to a step perturbation addition in the control action. (b) Total error comparison (C1= 1.5968, C2 = 0.1523, C3 = 0.0801).
(a)(b)Figure 9
(a) Controllers response to a step and ramp perturbation addition in the control action. (b) Total error comparison (C1 = 6.9351, C2 = 0.6542, C3 = 0.1977).
(a)(b)


### 4.4. Simulation Adding Parametric Uncertainties and Perturbations in the Control Action

This latest test aims to demonstrate how the controllers can fix a large drift. In this way,4.1 and 4.3 disturbances are considered simultaneously. Besides, the methodology proposed in [55] was also implemented to this ethanol bioprocess [65] and the results obtained are compared with those of the proposed technique. Figure 10 shows the biomass and ethanol profiles obtained with each controller throughout the process and a comparison of the errors. Here, C4 refers to a controller with one integrator and C5 with two ones. An improvement can be seen with the use of estimators concerning integrators (C1 = 100%; C2 = 12.5%; C3 = 1.9%; C4 = 13.1%; C5 = 2.4%). Furthermore, the error estimation has the advantage that it does not increase the order of the system by incorporating the estimation term and does not modify the form tuning as it does the integrators methodology, making mathematical development simpler.Figure 10
(a) Controllers response to parametric uncertainty and a step and ramp perturbation addition in the control action. (b) Total error comparison.
(a)(b)


## 4.1. Simulation under Parametric Uncertainties

In all bioprocess, model parameters may vary in an unpredictable way [66]. This can lead to structural instability in the system dynamical behavior [15]. Therefore, a strict and efficient control system is required to deal with this problem.In several research fields, probabilistic methods are useful for dealing with problems related to the robustness of systems affected by uncertainties [55]. Particularly, Monte Carlo Randomized Algorithm has been used for uncertainty quantification in many applications [67–69]. In this paper, the Monte Carlo Randomized Algorithm is applied using the procedure described in 2.3. The number of simulations, N = 1000, is obtained with equation (10), adopting a confidence (δ) of 0.01 and an accuracy (ε) of 0.005. Therefore, the following test demonstrates the technique success from a statistical point of view [70–72].In a simulation, a way to quantify uncertainties and perturbations is to specify the parameters real range of variation instead of using a constant value with greater error [73]. Hunag et al. [36] specified, in Table 1 of their work, the parameters confidence intervals for the ethanol process under study. For this test, in each simulation, all the system parameters are randomly changed by ±10% of their original range value (Table 1 of [36]). Then, the total error is calculated with equation (12). Table 3 shows the error range for each controller. Note how the total error range considerably decreases and maintains bounded to lower values when the estimation order increases. Figure 4 shows the total error for 1000 simulations to evaluate the three-controller performance under parametric uncertainty. As can be seen, the total tracking error range is visibly reduced.Table 3
Parametric uncertainty of ±10%.
ControllerTotal error rangeC10.5267–0.0098C20.2958–0.0056C30.0812–0.0049Figure 4
Total error for 1000 simulations under parametric uncertainty (±10%). (C1) Original controller; (C2) controller with a zero-order estimation; and (C3) controller with first-order estimator.Expressing the total error for each case into a percentage can easily quantify the controller improvement when estimators are incorporated. Figure5 suggests the better performance of the controllers C2 and C3 over C1. The worst case presented for C2 has a total error of 43.83% lower than that of the worst situation of C1, while the maximum error presented with C3 is 84.58% less than that of C1.Figure 5
Controllers comparison of maximum total error.Considering Theorem 1 of [55], as the tracking error (equation (11)) for 1000 simulations remains bounded, C2 and C3 controllers operation will be satisfactory with 99% of probability while the parameters vary within a ±10% range.


## 4.2. Simulation for Uncertainties Prevention

In the following test, controllers are tuned taking into account the worst possible situation presented in 4.1 forC1, to prevent it. Figure 6 shows the total error decreases when each controller is tuned. Note how even considering the most problematic situation of the system, the performance error can be reduced by 52.55% just by adding a zero-order estimator and by 56.15% with a first-order one.Figure 6
Total error decrease when the controllers are tuned considering the worst situation of C1 (total error:C1 = 0.1334, C2 = 0.0633, C3 = 0.0585).Then, the controllers using the chosen parameters were tested under normal operation conditions, demonstrating that their performance is not negatively affected. Figure7 compares the accumulated error for the three controllers. Moreover, total errors (C1 = 0.1416, C2 = 0.0062, C3 = 0.0031) contrast is shown.Figure 7
Errors comparison for the three controllers: (C1) original controller; (C2) controller with a zero-order estimation; (C3) controller with the first-order estimator.This test demonstrates another advantage of the proposed estimation technique. Table4 compares the errors of Figure 3 with those of Figure 7. In both cases, the controllers are tested under normal conditions; however, the difference lies in the tuning conditions considered. Analyzing Table 2, it can be observed that C2 and C3 controllers present similar results, which does not happen with C1. Therefore, it can be said that the operation of C2 and C3 will have a minimum error over the entire range of possible uncertainties between the nominal system parameters values and their worst variation.Table 4
Total error comparison under normal conditions.
C1C2C3Total errora)0.00630.00490.0030b)0.14160.00620.0031a) Tuned in normal conditions; b) tuned considering the worst situation of 4.1 for C1.


## 4.3. Simulation Adding Perturbations in the Control Action

In this test, a hypothetical situation that may produce an unexpected variation in the production is simulated. Firstly, a -30% step perturbation in the bioreactor feed rate is added to evaluate the response of controllers. Secondly, a ramp disturbance is added to the step perturbation previously presented. Figures8 and 9 show the control action variation compared to the reference and the percentage error, considering C1 error as 100%. Note how the total error improves when the new algorithm is applied. In the first test, it is reduced by 90.46% with C2 and by 94.99% with C3, while in the second test, the results are improved 91.87% with C2 and 97.15% with C3.Figure 8
(a) Controllers response to a step perturbation addition in the control action. (b) Total error comparison (C1= 1.5968, C2 = 0.1523, C3 = 0.0801).
(a)(b)Figure 9
(a) Controllers response to a step and ramp perturbation addition in the control action. (b) Total error comparison (C1 = 6.9351, C2 = 0.6542, C3 = 0.1977).
(a)(b)


## 4.4. Simulation Adding Parametric Uncertainties and Perturbations in the Control Action

This latest test aims to demonstrate how the controllers can fix a large drift. In this way,4.1 and 4.3 disturbances are considered simultaneously. Besides, the methodology proposed in [55] was also implemented to this ethanol bioprocess [65] and the results obtained are compared with those of the proposed technique. Figure 10 shows the biomass and ethanol profiles obtained with each controller throughout the process and a comparison of the errors. Here, C4 refers to a controller with one integrator and C5 with two ones. An improvement can be seen with the use of estimators concerning integrators (C1 = 100%; C2 = 12.5%; C3 = 1.9%; C4 = 13.1%; C5 = 2.4%). Furthermore, the error estimation has the advantage that it does not increase the order of the system by incorporating the estimation term and does not modify the form tuning as it does the integrators methodology, making mathematical development simpler.Figure 10
(a) Controllers response to parametric uncertainty and a step and ramp perturbation addition in the control action. (b) Total error comparison.
(a)(b)


## 5. Conclusions

This paper presents an improvement for a tracking control strategy previously published [50]. This technique lets tracking reference concentration profiles, even in the presence of model uncertainties and external perturbances. To consider those uncertainties, a new term is included in the mathematical model. The error estimation is approximated with Newton’s backward interpolation. In this manuscript, the main contribution is to decrease additive uncertainties effect on the tracking error without increasing the controller mathematical complexity. Moreover, the controller tuning is simpler than in conventional controllers, since varying the parameters between zero and one the error convergence to zero is achieved.The Monte Carlo Randomized Algorithm is used to tune the controllers and carry out the tests with parametric uncertainty. Those tests show the effectiveness of this methodology, which was demonstrated in Sections3 and 4.Comparing with other methodologies that deal with similar uncertain control problems, such as [74–79], the proposed controller presents the advantage of avoiding the stochastic modeling needed to deal with parameters under perturbation of white noise. Besides, this nonlinear control does not require a great mathematical effort and does not add significant complexity to the original controller.



---

*Source: 2865264-2021-05-18.xml*