# Multivariable Tracking Control of a Bioethanol Process under Uncertainties

Authors: M. Cecilia Fernández; M. Nadia Pantano; Emanuel Serrano; Gustavo Scaglia
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Bioethanol is one of the most studied alternative fuels nowadays. Due to its production process complexity and the low quality of the mathematical models that describe it, a reliable controller is needed to maximize the fuel production and minimize its environmental impact, even in the presence of uncertainty. Here, a controller for tracking optimal profiles considering model errors and external perturbations is proposed. This work is an improvement of a previously presented technique. To reduce the earlier mentioned uncertainties’ effect during the fermentation, some tracking error integrators are added in the control action calculation. This simple modification ensures the tracking error convergence to zero, even in the presence of uncertainties (demonstration available). Different tests are carried out and a performance comparison with the original controller is shown to highlight improvements in the tracking error of up to 98% when integrators are incorporated. Furthermore, a classical PI controller is contrasted with the proposed techniques.


---

## Body

## 1. Introduction

Excessive consumption of fossil fuels has led to a depletion of this natural resource and an environmental contamination caused by the greenhouse gases generated in their combustion [1, 2]. Some of the most notorious consequences are the global warming, acid rain, and urban smog [3]. For this reason, renewable and nonpolluting alternative fuels have become an investigation focus. One of the most important research lines in this topic is the bioethanol production, which is less toxic, readily biodegradable, and produces less air-borne pollutants than petroleum fuel [4]. Its production process involves a fermentation that transforms carbohydrates into alcohol. Moreover, it can be obtained from different raw materials, which classify it into three categories [4]: (i) first generation, involves feedstock like sugar cane, sugar beet, wheat, sorghum, fruits, corn, potato, rice, sweet potato, or barley [5, 6]; (ii) second generation, uses lignocellulosic biomass from wood, straw, grasses, olive mill [7], or pineapple waste [8], among others [9]; and (iii) third generation comes from macro- and microalgae biomass [10]. Furthermore, there are new investigations oriented to obtain bioethanol from fecal waste [1, 11]. However, the second generation is the most studied nowadays because it does not generate a competition with the food industry as the first one and is an alternative of final disposal for many types of waste.The bioethanol production process consists of a series of stages; some of them depend on the raw material used:(i)
Pretreatment: it allows altering the lignin and hemicellulose structure of lignocellulosic materials to leave cellulose exposed; this stage facilitates the next one and increases its efficiency.(ii)
Hydrolysis: with this treatment, long chains of carbohydrates are reduced into monomeric sugars.(iii)
Fermentation: biological process ensures the transformation of the monomeric sugars into alcohol molecules by some microorganism (yeast or bacteria).(iv)
Distillation: the alcohol is separated from the culture medium. This method requires a significant amount of energy, which makes the process more expensive.Many researchers have devoted their efforts to improve each one of the stages mentioned. Herrero, et al. [7] presented the optimization results of the sulfuric acid pretreatment variables applied to improve sugars bioavailability. Aimaretti et al. [12], studied two different enzymatic hydrolysis strategies to increase fermentable sugar concentration in the must. Kumar et al. [13] demonstrated how to increase the efficiency of sugar conversion and the ethanol yield of a fermentation applying gas stripping at high temperature. In Tgarguifa et al.’s study [14], the model and a distillation column optimization to produce ethanol from wine are proposed.The fermentation step can be improved with many strategies, such as the specific microorganism selection or genetic modification of those microorganisms [3, 4]. Fed-batch bioprocesses are intensively studied nowadays because of their main advantages [15]. For example, the medium substrate concentration can be regulated by an appropriate feed rate profile [16], obtaining better production yields and minimizing the production costs [15, 17]. However, bioprocesses control is required to follow a certain feed flow rate and get stability and the best productivity [18]. Moreover, a real-time bioprocess detailed control is a complex target but is necessary to ensure raw materials optimal use, water and energy-saving, final product consistent quality, a reduction in wastes and process cycle time, replacement of costly and slow laboratory testing, and continuous learning, which opens up the possibility of bioprocess innovation. Furthermore, the mathematical representation of the process is the key to achieve good results.A mathematical model provides a map from inputs to responses. A model quality depends on how closely its responses match those of the true plant. However, a model set which includes the true physical plant can never be constructed [19]. Generally, a bioprocess modelling 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 many modelling uncertainties [20]. Furthermore, sometimes those parameters are determined without a previous model analysis, or their values are not informed with their respective confidence intervals [21, 22]. Moreover, the time-varying parameters are usually assumed constant [23]. Also, uncertainties associated with processing technologies’ parameters are rarely considered [24]. Dismissal of aspects leads to a reality misrepresentation and, as a result, a deficient performance with serious and imminent risks [25, 26]. Hence, the main task to assure the bioprocess quality implies finding a way to control these flaws [27–29]. 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 [21].The problem of model parameters determination has motivated the development of many strategies. Donoso-Bravo et al. [30] discussed existing methodologies for parameter estimation and model validation in the field of anaerobic digestion processes. Balsa-Canto and Banga [31] proposed an iterative procedure to improve practical identifiability and to help in the design of informative experiments. Vilas, et al. [32] presented the theoretical background of a parameter identification protocol intended to deal with those mentioned challenges. To reduce the complexity of tuning a complex anaerobic digestion model, a particle swarm optimization-based smart algorithm was developed to estimate all parameters [33]. Elenchezhiyan and Prakash [34] formulated state estimation schemes for nonlinear hybrid dynamic systems subjected to stochastic state disturbances and random errors in measurements using interacting multiple-model algorithms. In Pantano et al.’s study [35], the problem of optimal profiles tracking control under uncertainties for a fed-batch bioprocess with two control actions is addressed. In Lara-Cisneros et al.’s study [36], an extreme-seeking scheme based on an approach to variable structure control for fed-batch bioreactors, which deals with uncertainty on the specific growth rate, is proposed. For some other examples, see [37] and [38].On the other hand, many scientists have developed several feedback control strategies to deal with bioprocess uncertainties, such as optimal control, nonlinear model predictive control, fuzzy control, hybrid control, adaptive control, tracking control, and neural network [37, 39–43]. Due to the online implementation difficulty, the high computational cost, the imprecise mathematical models, and online solutions, their use for bioprocesses is limited [17] and is still a research topic [44].Generally, for disturbances treatment, strategies to reject or limit their effect are used. Among the most frequent approaches can be found prefeeding, disturbance observer, disturbance and uncertainty estimation, or control for active rejection of disturbances [45–49]. Disturbance estimation needs a filter in their structure design (to ensure infinite norm less than one), which increases the mathematical complexity.For the specific case of ethanol production defined by Hunag, et al. [50], Fernández, et al. [51] presented a controller design focused on looking for control actions, to track predefined concentration profiles. As the controller structure comes from the process mathematical model, it can be implemented in many systems. This procedure is simple, versatile, and precise. Besides, only basic knowledge of numerical methods and linear algebra are needed to implement it. The technique was tested against different disturbances, showing an excellent performance and the tracking error convergence to zero. However, its weakness is handling model uncertainties and external perturbations.In this paper, the control strategy presented in Fernández, et al. [51] is improved. This new alternative allows tracking predefined optimal profiles considering additive uncertainty in an underactuated system with only one control action. To achieve this goal, a new term is added to the original controller design equation. This term represents the model uncertainties, which are avoided using tracking error integrators. It is important to highlight that as the number of integrators increases, the approximation improves and the tracking error decreases. Furthermore, this adjustment ensures signals uniformity and the tracking error convergence to zero. In this way, a positive answer to the tracking control challenging problem in multivariable nonlinear systems with additive uncertainty is provided.The paper is organized in four sections. Section2 summarizes the ethanol bioprocess and the previously presented controller technique. Section 3 explains the new controller design, which takes into account the uncertainties. Section 4 shows simulations testing the new controller performance and a comparison with other control techniques. Finally, conclusions are exposed.


## 2. Ethanol Bioprocess Model and Controller Design

### 2.1. Ethanol Bioprocess Description

The following mathematical model describes a fed-batch ethanol fermentation. It was originally proposed by Hunag, et al. [50]. Here, the yeast Saccharomyces diastaticus is used. The system input, feed rate (U), is a 50–50 glucose and fructose mixture, while the states are cells (X), ethanol (P1), glycerol (P2), glucose (S1), and fructose (S2) concentration inside the reactor, whose profiles variation over time is wanted to be tracked:(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,μ1=μm1S1KS1+S1+S12/KS1IKP1KP1+P1+P12/KP1I,μ2=μm2S2KS2+S2+S22/KS2IKP2KP2+P1+P12/KP2I,qS1/P1=νS1P1S1KS1P1+S1kS1P1kS1P1+P1,qS2/P2=νS2P2S2KS2P2+S2kS2P2kS2P2+P2,qS1/P2=νS1P2S1KS1P2+S1kS1P2kS1P2+P2,qS2/P1=νS2P1S2KS2P1+S2kS2P1kS2P1+P1.where V is the culture volume, µ1 and µ2 are the specific yeast cells growth rates, qS1/P1 and qS2/P1 are the specific ethanol production rates, qS1/P2 and qS2/P2 are the specific glycerol production rates, in all cases from glucose and fructose, respectively. U is the control variable.Hunag, et al. [50] applied a run-to-run optimization procedure in order to estimate kinetic model parameters, fed glucose concentration, feed strategy, and fermentation time to maximize the ethanol production rate. Those results are used as references in this paper. A reactor with a total working volume of 5 L was used. Temperature was maintained at 35.8°C by controlling the circulation of the cooling water. The airflow was fixed at 1.5 vvm, and the pH was kept at 5.0 by adding 1N NaOH. The biomass concentration was determined spectrophotometrically at 540 nm, and the corresponding dry weight of cells was obtained from a calibration curve. The concentrations of glucose, fructose, ethanol, and glycerol were analyzed with a high-performance liquid chromatograph [50]. Table 1 shows system’s initial conditions, while parameters’ nomenclature, description, and values are exposed in Table 2.Table 1
Initial conditions for ethanol fermentation.
Variable
Initial value
X (g/L)
1.5
P
1 (g/L)
5.3
P
2 (g/L)
0.0001
S
1 (g/L)
8.6
S
2 (g/L)
8.6
V (L)
1.35Table 2
Nomenclature, description, and values of parameters.
Parameter
Description
Value
µ
m1
Maximum specific growth rate coefficient for yeast on glucose (h−1)
1.8823
µ
m2
Maximum specific growth rate coefficient for yeast on fructose (h−1)
1.7098
Y
P1/S1
Yield coefficient for ethanol from glucose
0.5085
Y
P2/S1
Yield coefficient for glycerol from glucose
0.5331
Y
P1/S2
Yield coefficient for ethanol from fructose
0.5098
Y
P2/S2
Yield coefficient for glycerol from fructose
0.4462
K
S1
Saturation coefficient for cell growth on glucose (g/L)
159.75
K
S1I
Inhibition coefficient for cell growth on glucose (g/L)
94.233
K
P1
Saturation coefficient for cell growth on ethanol (g/L)
238.39
K
P1I
Inhibition coefficient for cell growth on ethanol (g/l)
2.7378
K
S2
Saturation coefficient for cell growth on fructose (g/L)
0.0726
K
S2I
Inhibition coefficient for cell growth on fructose (g/L)
9.0048
K
P2
Saturation coefficient for cell growth on glycerol (g/L)
35.958
K
P2I
Inhibition coefficient for cell growth on glycerol (g/L)
9.9722
K
S1P1
Saturation coefficient for ethanol production on glucose (g/L)
1.3409
k
S1P1
Inhibition coefficient for ethanol production on glucose (g/L)
18.612
K
S2P1
Saturation coefficient for ethanol production on fructose (g/L)
0.9129
k
S2P1
Inhibition coefficient for ethanol production on fructose (g/L)
1000
K
S1P2
Saturation coefficient for glycerol production on glucose (g/L)
6.7116
k
S1P2
Inhibition coefficient for glycerol production on glucose (g/L)
0.5863
K
S2P2
Saturation coefficient for glycerol production on fructose (g/L)
0.4310
k
S2P2
Inhibition coefficient for glycerol production on fructose (g/L)
1.150
ν
S1P1
Coefficient of maximum specific ethanol production rate for yeast on glucose (h−1)
1.5051
ν
S2P1
Coefficient of maximum specific ethanol production rate for yeast on fructose (h−1)
0.3321
ν
S1P2
Coefficient of maximum specific glycerol production rate for yeast on glucose (h−1)
0.0023
ν
S2P2
Coefficient of maximum specific glycerol production rate for yeast on fructose (h−1)
0.1609
Λ
Proportion of glucose and fructose
0.5
S
f
Sugar total feed concentration (g/L)
300


### 2.2. Controller Design

The methodology described in [51] aims to find U in order to make the system follow preestablished state variables profiles (references) with minimum error. For the controller design is assumed that reference profiles and states at each sampling instant are known. Figure 1 shows the reference concentration profiles and the feed flow rate. Next, the mentioned methodology is summarized.Figure 1
Cells, ethanol, glycerol, and fructose reference concentrations along the process. Reference feed flow rate.First, equation (1) is discretized using numerical methods. Here, Euler is applied due to its simplicity and good results [51].(3)dσdt=σn+1−σnTS,where σ symbolizes each state variable, σn is the current σ value measured from the reactor (online), σn+1 is σ value in the next measurement instant, and TS is the sampling time (0.1 h) [52]. The process lasts 15.7 h (Tf).The state variables inn + 1 are approximated as follows:(4)σref,n+1−σn+1︸errorn+1=kσσref,n−σn︸errorn⟶σn+1=σrefn+1−kσσref,n−σn.where σref is the reference state variables and kσ represents the controller parameter for the variable σ. For this system, the controller parameters are kX, kP1, kP2, kS1, and kS2. As expected is the error reduction in each sampling time, it must be fulfilled that 0 ≤ kσ < 1. Then, replacing equation (4) in (3),(5)dσdt=σrefn+1−kσσrefn−σn︷σn+1−σnTS=Δσ.Applying equation (5) in (1),(6)ΔX=µ1+µ2Xn−UnVnXn,ΔS1=−qS1/P1YP1/S1+qS1/P2YP2/S1Xn+UnVnλSf−S1,n,ΔS2=−qS2/P1YP1/S2+qS2/P2YP2/S2Xn+UnVn1−λSf−S2,n,ΔP1=qS1/P1+qS2/P1Xn−UnVnP1,n,ΔP2=qS1/P2+qS2/P2Xn−UnVnP2,n.Expressing (6) in a matrix form,(7)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/Vn︷AUn︷u=ΔX−µ1+µ2XnΔS1+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1−qS1/P1+qS2/P1XnΔP2−qS1/P2+qS2/P2Xn︸b.To findU, the equation system (7) must have an exact solution. So, b has to be a linear combination of A columns [53], in other words, A and b must be parallel. One way to satisfy this condition is as follows:(8)cosA,b=A,bA∗b=1,where the operation between < > and ||.|| represents the inner product and the vectors norm in Rn space, respectively.At this point, a “sacrificed variable” selection is necessary. This variable, defined asS1ez, ensures that (7) has an exact solution. For more details on its selection and calculation, see [51]. Finally, Un is obtained using least squares [53].(9)Un=ATA−1ATb.


### 2.3. Controller Tuning

In [51], the Monte Carlo algorithm was used to find the kσ combination that minimizes the accumulated error. This methodology consists in simulating the process N times using random kσ [54]:(10)N≥log1/δlog1/1−ε,where δ is the confidence and ε the accuracy.Then, two new terms are introduced, “tracking error” and “total error”, equations (11) and (12), respectively:(11)en=Xref,n−Xn/maxXref,n2+P1ref,n−P1,n/maxP1ref,n2+P2ref,n−P2,n/maxP2ref,n2+S2ref,n−S2,n/maxS2ref,n2,(12)Ep=∑nTS=1Jen,where p represents the simulation in progress, p=1,2,…,N and nTS the measurement instant, nTS = 1, 2, …, J.Equation (12) is the function cost to be minimized with the Monte Carlo Algorithm. Then, kσ that minimize the total error are selected.Sequence of the Monte Carlo algorithm is as follows [55]:(1)
Define the controller’s parameters to be optimized.(2)
Determine the number of simulations to be performed (N). The values of δ and ε are chosen depending on the desired accuracy, for example, δ = 0.01 and ε = 0.005, resulting in N = 1000.(3)
A value is randomly assigned for each parameter of the controller.(4)
The process is simulated and the cost index (Ep) is calculated.(5)
Repeat steps 3 and 4 until completing theN iterations.(6)
Finally, the controller parameters that minimizeEp are selected.Theorem.
If the discrete system is given by equations (1), the control action is calculated with equations (9), and kσ takes values between zero and one (0 < kσ < 1), then, the tracking error convergence to zero when n tends to infinity is achieved.Demonstration. Replacing the sacrificed variable in (7) and expressing the matrix system generically:(13)a1a2a3a4a5︸AUn=b1b2b3b4b5︸b.Solving (13) with least squares,(14)Un=ATA−1ATb=a1b1+a2b2+a3b3+a4b4+a5b5a12+a22+a32+a42+a52.From (13),(15)a1a2=b1b2⟶b2=a2a1b1,a1a3=b1b3⟶b3=a3a1b1,a1a4=b1b4⟶b4=a4a1b1,a1a5=b1b5⟶b5=a5a1b1.Replacing (15) in (14),(16)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.Substituting (16) in (7),(17)Xn+1=Xref,n+1−kXXref,n−Xn+TSμ1S1ez,n,P1,n−μ1S1,n,P1,nXn.Then, the tracking error forX is defined as(18)eX,n+1=Xref,n+1−Xn+1.Introducing (17) in (18),(19)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n−μ1S1,n,P1,nXn.Theµ1 (S1n, P1n) Taylor approximation in the desired value µ1 (S1ez n, P1n) is(20)μ1S1,n,P1,n=μ1S1ez,n,P1,n+dμ1S1,P1,ndS1S1n−S1ez,nS1=Sez,n+θS1,n−S1ez,n=Sθ,where⟶0<θ<1.Substituting (20) in (19),(21)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n+dμ1S1,P1,ndS1SθS1n−S1ez,n︸−eS1,n−μ1S1ez,n,P1,nXn,eX,n+1=kXXref,n−Xn+TSdμ1S1,P1,ndS1SθeS1,nXn.Proceeding in a similar way for the other variables and joining the final expressions,(22)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1Sθ00−dqS1/P1S1,P1,ndS1Sφ−dqS1/P2S1,P2,ndS1Sα︸NLeS1,n.In equation (22), L is a linear system and NL is a bounded nonlinearity [51]. Note that if kσ = 0, the reference is reached in one step. So, if 0 < kσ < 1, the tracking error tends to zero when n ⟶ ∞ [51, 56].


### 2.4. Controller Implementation

The control system diagram is shown in Figure2, and the mathematical procedure is summarized in the following steps:Step 1: DefineTs and σref. Read σn.Step 2: Approximate the differential equations with numerical methods, equations (3).Step 3: Approximate the state variables inn + 1 using equation (4).Step 4: Define and calculate the sacrificed variable.Step 5: CalculateUn using least squares, equation (9).Figure 2
Control system diagram.


## 2.1. Ethanol Bioprocess Description

The following mathematical model describes a fed-batch ethanol fermentation. It was originally proposed by Hunag, et al. [50]. Here, the yeast Saccharomyces diastaticus is used. The system input, feed rate (U), is a 50–50 glucose and fructose mixture, while the states are cells (X), ethanol (P1), glycerol (P2), glucose (S1), and fructose (S2) concentration inside the reactor, whose profiles variation over time is wanted to be tracked:(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,μ1=μm1S1KS1+S1+S12/KS1IKP1KP1+P1+P12/KP1I,μ2=μm2S2KS2+S2+S22/KS2IKP2KP2+P1+P12/KP2I,qS1/P1=νS1P1S1KS1P1+S1kS1P1kS1P1+P1,qS2/P2=νS2P2S2KS2P2+S2kS2P2kS2P2+P2,qS1/P2=νS1P2S1KS1P2+S1kS1P2kS1P2+P2,qS2/P1=νS2P1S2KS2P1+S2kS2P1kS2P1+P1.where V is the culture volume, µ1 and µ2 are the specific yeast cells growth rates, qS1/P1 and qS2/P1 are the specific ethanol production rates, qS1/P2 and qS2/P2 are the specific glycerol production rates, in all cases from glucose and fructose, respectively. U is the control variable.Hunag, et al. [50] applied a run-to-run optimization procedure in order to estimate kinetic model parameters, fed glucose concentration, feed strategy, and fermentation time to maximize the ethanol production rate. Those results are used as references in this paper. A reactor with a total working volume of 5 L was used. Temperature was maintained at 35.8°C by controlling the circulation of the cooling water. The airflow was fixed at 1.5 vvm, and the pH was kept at 5.0 by adding 1N NaOH. The biomass concentration was determined spectrophotometrically at 540 nm, and the corresponding dry weight of cells was obtained from a calibration curve. The concentrations of glucose, fructose, ethanol, and glycerol were analyzed with a high-performance liquid chromatograph [50]. Table 1 shows system’s initial conditions, while parameters’ nomenclature, description, and values are exposed in Table 2.Table 1
Initial conditions for ethanol fermentation.
Variable
Initial value
X (g/L)
1.5
P
1 (g/L)
5.3
P
2 (g/L)
0.0001
S
1 (g/L)
8.6
S
2 (g/L)
8.6
V (L)
1.35Table 2
Nomenclature, description, and values of parameters.
Parameter
Description
Value
µ
m1
Maximum specific growth rate coefficient for yeast on glucose (h−1)
1.8823
µ
m2
Maximum specific growth rate coefficient for yeast on fructose (h−1)
1.7098
Y
P1/S1
Yield coefficient for ethanol from glucose
0.5085
Y
P2/S1
Yield coefficient for glycerol from glucose
0.5331
Y
P1/S2
Yield coefficient for ethanol from fructose
0.5098
Y
P2/S2
Yield coefficient for glycerol from fructose
0.4462
K
S1
Saturation coefficient for cell growth on glucose (g/L)
159.75
K
S1I
Inhibition coefficient for cell growth on glucose (g/L)
94.233
K
P1
Saturation coefficient for cell growth on ethanol (g/L)
238.39
K
P1I
Inhibition coefficient for cell growth on ethanol (g/l)
2.7378
K
S2
Saturation coefficient for cell growth on fructose (g/L)
0.0726
K
S2I
Inhibition coefficient for cell growth on fructose (g/L)
9.0048
K
P2
Saturation coefficient for cell growth on glycerol (g/L)
35.958
K
P2I
Inhibition coefficient for cell growth on glycerol (g/L)
9.9722
K
S1P1
Saturation coefficient for ethanol production on glucose (g/L)
1.3409
k
S1P1
Inhibition coefficient for ethanol production on glucose (g/L)
18.612
K
S2P1
Saturation coefficient for ethanol production on fructose (g/L)
0.9129
k
S2P1
Inhibition coefficient for ethanol production on fructose (g/L)
1000
K
S1P2
Saturation coefficient for glycerol production on glucose (g/L)
6.7116
k
S1P2
Inhibition coefficient for glycerol production on glucose (g/L)
0.5863
K
S2P2
Saturation coefficient for glycerol production on fructose (g/L)
0.4310
k
S2P2
Inhibition coefficient for glycerol production on fructose (g/L)
1.150
ν
S1P1
Coefficient of maximum specific ethanol production rate for yeast on glucose (h−1)
1.5051
ν
S2P1
Coefficient of maximum specific ethanol production rate for yeast on fructose (h−1)
0.3321
ν
S1P2
Coefficient of maximum specific glycerol production rate for yeast on glucose (h−1)
0.0023
ν
S2P2
Coefficient of maximum specific glycerol production rate for yeast on fructose (h−1)
0.1609
Λ
Proportion of glucose and fructose
0.5
S
f
Sugar total feed concentration (g/L)
300


## 2.2. Controller Design

The methodology described in [51] aims to find U in order to make the system follow preestablished state variables profiles (references) with minimum error. For the controller design is assumed that reference profiles and states at each sampling instant are known. Figure 1 shows the reference concentration profiles and the feed flow rate. Next, the mentioned methodology is summarized.Figure 1
Cells, ethanol, glycerol, and fructose reference concentrations along the process. Reference feed flow rate.First, equation (1) is discretized using numerical methods. Here, Euler is applied due to its simplicity and good results [51].(3)dσdt=σn+1−σnTS,where σ symbolizes each state variable, σn is the current σ value measured from the reactor (online), σn+1 is σ value in the next measurement instant, and TS is the sampling time (0.1 h) [52]. The process lasts 15.7 h (Tf).The state variables inn + 1 are approximated as follows:(4)σref,n+1−σn+1︸errorn+1=kσσref,n−σn︸errorn⟶σn+1=σrefn+1−kσσref,n−σn.where σref is the reference state variables and kσ represents the controller parameter for the variable σ. For this system, the controller parameters are kX, kP1, kP2, kS1, and kS2. As expected is the error reduction in each sampling time, it must be fulfilled that 0 ≤ kσ < 1. Then, replacing equation (4) in (3),(5)dσdt=σrefn+1−kσσrefn−σn︷σn+1−σnTS=Δσ.Applying equation (5) in (1),(6)ΔX=µ1+µ2Xn−UnVnXn,ΔS1=−qS1/P1YP1/S1+qS1/P2YP2/S1Xn+UnVnλSf−S1,n,ΔS2=−qS2/P1YP1/S2+qS2/P2YP2/S2Xn+UnVn1−λSf−S2,n,ΔP1=qS1/P1+qS2/P1Xn−UnVnP1,n,ΔP2=qS1/P2+qS2/P2Xn−UnVnP2,n.Expressing (6) in a matrix form,(7)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/Vn︷AUn︷u=ΔX−µ1+µ2XnΔS1+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1−qS1/P1+qS2/P1XnΔP2−qS1/P2+qS2/P2Xn︸b.To findU, the equation system (7) must have an exact solution. So, b has to be a linear combination of A columns [53], in other words, A and b must be parallel. One way to satisfy this condition is as follows:(8)cosA,b=A,bA∗b=1,where the operation between < > and ||.|| represents the inner product and the vectors norm in Rn space, respectively.At this point, a “sacrificed variable” selection is necessary. This variable, defined asS1ez, ensures that (7) has an exact solution. For more details on its selection and calculation, see [51]. Finally, Un is obtained using least squares [53].(9)Un=ATA−1ATb.


## 2.3. Controller Tuning

In [51], the Monte Carlo algorithm was used to find the kσ combination that minimizes the accumulated error. This methodology consists in simulating the process N times using random kσ [54]:(10)N≥log1/δlog1/1−ε,where δ is the confidence and ε the accuracy.Then, two new terms are introduced, “tracking error” and “total error”, equations (11) and (12), respectively:(11)en=Xref,n−Xn/maxXref,n2+P1ref,n−P1,n/maxP1ref,n2+P2ref,n−P2,n/maxP2ref,n2+S2ref,n−S2,n/maxS2ref,n2,(12)Ep=∑nTS=1Jen,where p represents the simulation in progress, p=1,2,…,N and nTS the measurement instant, nTS = 1, 2, …, J.Equation (12) is the function cost to be minimized with the Monte Carlo Algorithm. Then, kσ that minimize the total error are selected.Sequence of the Monte Carlo algorithm is as follows [55]:(1)
Define the controller’s parameters to be optimized.(2)
Determine the number of simulations to be performed (N). The values of δ and ε are chosen depending on the desired accuracy, for example, δ = 0.01 and ε = 0.005, resulting in N = 1000.(3)
A value is randomly assigned for each parameter of the controller.(4)
The process is simulated and the cost index (Ep) is calculated.(5)
Repeat steps 3 and 4 until completing theN iterations.(6)
Finally, the controller parameters that minimizeEp are selected.Theorem.
If the discrete system is given by equations (1), the control action is calculated with equations (9), and kσ takes values between zero and one (0 < kσ < 1), then, the tracking error convergence to zero when n tends to infinity is achieved.Demonstration. Replacing the sacrificed variable in (7) and expressing the matrix system generically:(13)a1a2a3a4a5︸AUn=b1b2b3b4b5︸b.Solving (13) with least squares,(14)Un=ATA−1ATb=a1b1+a2b2+a3b3+a4b4+a5b5a12+a22+a32+a42+a52.From (13),(15)a1a2=b1b2⟶b2=a2a1b1,a1a3=b1b3⟶b3=a3a1b1,a1a4=b1b4⟶b4=a4a1b1,a1a5=b1b5⟶b5=a5a1b1.Replacing (15) in (14),(16)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.Substituting (16) in (7),(17)Xn+1=Xref,n+1−kXXref,n−Xn+TSμ1S1ez,n,P1,n−μ1S1,n,P1,nXn.Then, the tracking error forX is defined as(18)eX,n+1=Xref,n+1−Xn+1.Introducing (17) in (18),(19)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n−μ1S1,n,P1,nXn.Theµ1 (S1n, P1n) Taylor approximation in the desired value µ1 (S1ez n, P1n) is(20)μ1S1,n,P1,n=μ1S1ez,n,P1,n+dμ1S1,P1,ndS1S1n−S1ez,nS1=Sez,n+θS1,n−S1ez,n=Sθ,where⟶0<θ<1.Substituting (20) in (19),(21)eX,n+1=kXXref,n−Xn−TSμ1S1ez,n,P1,n+dμ1S1,P1,ndS1SθS1n−S1ez,n︸−eS1,n−μ1S1ez,n,P1,nXn,eX,n+1=kXXref,n−Xn+TSdμ1S1,P1,ndS1SθeS1,nXn.Proceeding in a similar way for the other variables and joining the final expressions,(22)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1Sθ00−dqS1/P1S1,P1,ndS1Sφ−dqS1/P2S1,P2,ndS1Sα︸NLeS1,n.In equation (22), L is a linear system and NL is a bounded nonlinearity [51]. Note that if kσ = 0, the reference is reached in one step. So, if 0 < kσ < 1, the tracking error tends to zero when n ⟶ ∞ [51, 56].


## 2.4. Controller Implementation

The control system diagram is shown in Figure2, and the mathematical procedure is summarized in the following steps:Step 1: DefineTs and σref. Read σn.Step 2: Approximate the differential equations with numerical methods, equations (3).Step 3: Approximate the state variables inn + 1 using equation (4).Step 4: Define and calculate the sacrificed variable.Step 5: CalculateUn using least squares, equation (9).Figure 2
Control system diagram.


## 3. Controller Design under Uncertainties

As a main contribution, the following procedure suggests an important improvement for Fernández et al.’s [51] controller performance.


### 3.1. Controller Design under Uncertainty

In this section, the uncertainties effect in the tracking error is evaluated. With the aim of representing the additive uncertainty effect in the process mathematical model, the termsEσ,n are added in equation (7):(23)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/VnUn=ΔX−µ1+µ2XnΔS1+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1−qS1/P1+qS2/P1XnΔP2−qS1/P2+qS2/P2Xn−1TSEX,nES1,nES2,nEP1,nEP2,n︸En.The additive uncertainty (En) models several kinds of mismatches as well as perturbed systems. En might depend on the states variables and the system input. Moreover, considering a real plant like zn+1 = f (zn,un), 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 [57, 58], then Eσ,n can be modelled as a bounded uncertainty [59, 60].Taking into account this uncertainty and following the same procedure of Section2.3, the next final expression is obtained:(24)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1Sθ00−dqS1/P1S1,P1,ndS1Sφ−dqS1/P2S1,P2,ndS1SαeS1,n︸NL−EX,nES1,nES2,nEP1,nEP2,n︸En.Comparing (24) with (22), the error nonconvergence to zero is noticed due to En presence. Therefore, the next step is En approximation to reduce their effect in the tracking error, achieving the convergence to zero.


### 3.2. Integral Action

Here, a novel approach is presented to reduce the effect ofEn in the tracking error by adding some integrators in the control action calculus equation. The advantage is that new terms do not add extra mathematical complexity and maintain the simplicity of the scheme presented by Fernández et al. [51].As it was explained,En represents the additive uncertainty, typical of mathematical model formulation errors and the influence of external perturbation. Those uncertainties can not be measured, that is to say, they are unknown, but they can be assumed as a polynomial so that all kinds of possible uncertainties can be considered, which can be constant, linear, quadratic, or other. In this way, by increasing the order of the polynomial, the uncertainty approximation is improved and there is a tendency to reduce the tracking error more and more. This statement is mathematically demonstrated below.E
n differences are defined as(25)δEn=En+1−En,δ2En=δδEn=δEn+1−En=En+2−2En+1+En,δqEn=δδq−1En.


#### 3.2.1. Constant Uncertainty

IfEn is assumed as a constant, then δEn = 0. Denoting generally the error as eσ,n,(26)eσ,n=σref,n−σn.An integrator could force the tracking error convergence to zero. This integrator is defined as follows:(27)Uσ,n+1=Uσ,n+∫nTSn+1TSeσ,ntdt≅Uσ,n+eσ,nTS.Consequently, equation (4) can be rewritten as(28)σn+1=σrefn+1−kσσrefn−σn+lσUσ,n+1,where kσ and lσ are the proportional and integral controller parameters, respectively.Demonstration of error convergence to zero when uncertainties are considered is given below.Expressing equation (24) in the following summarized way,(29)eσ,n+1=kσeσ,n+NLn+En.Adding the integral,(30)eσ,n+1=kσeσ,n−lσUσ,n+1+NLn+En.Replacing equation (27) in (30),(31)eσ,n+1=kσeσ,n−lσUσ,n+eσ,nT0+NLn+En.Then, combining equation (30) with (31),(32)eσ,n+2=kσeσ,n+1−lσUσ,n+1+eσ,n+1T0+NLn+1+En+1.Otherwise,(33)eσ,n+2=kσeσ,n+1−NLn+En+kσeσ,n−eσ,n+1−lσeσ,n+1T0+NLn+1+En+1,eσ,n+2=eσ,n+1kσ−lσT0+1−kσeσ,n−NLn+NLn+1−En+En+1︸δEn=0.Remark. If En = cte, δEn = 0. Looking equation (33), it is evidenced that the tracking error tends to zero with the progress of the process.


#### 3.2.2. Linear Uncertainty

IfEn is supposed as a linear function, then δ2En = 0. Following the same reasoning as in 3.2.1, now two integrators are needed to compensate the added uncertainty. Those integrators are(34)Uσ1,n+1=Uσ1,n+∫nTSn+1TSeσntdt≅Uσ1,n+eσnTS,Uσ2,n+1=Uσ2,n+∫nTSn+1TSUσ1,n+1tdt≅Uσ2,n+Uσ1,n+1TS.Accordingly, equation (4) can be rewritten as(35)σn+1=σref,n+1−kσσref,n−σn+lσ1Uσ1,n+1+lσ2Uσ2,n+1,where kσ and lσ1 and lσ2 are the proportional and both integral controller parameters, respectively.After the same steps of3.2.1 demonstration, the error convergence to zero is achieved.


#### 3.2.3. Polynomial Uncertainty

In this case,En is supposed as an m-order polynomial function. So, if q > m, δqEn = 0. Using the same steps before mentioned, it should be added m integrators in order to compensate the uncertainties considered.Figure3 shows the block diagram of the new control system.Figure 3
New control system block diagram.


### 3.3. Controller Parameters Selection

For the error to decrease progressively and remain at low values, the linear term of equation (24), L, must tend to zero. In Fernández et al.’s study [51], this state was possible choosing kσ values between zero and one. However, when integrators are considered, two types of parameters are involved, proportional and integral, what changes the parameters selection condition to achieve the error convergence to zero.In this way, starting from the characteristic equation of the linear term of equation (33),(36)eσ,n+2+−kσ+lσTS−1︸a1eσ,n+1+kσ︸a2eσ,n=−NLn+NLn+1,r2+a1r+a2=0.If the characteristic equation roots (r1 and r2) are between zero and one, it is guaranteed that the linear term tends to zero when n tends to infinity. Hence, the root values must be determined to calculate parameters.As many authors recommend, the Monte Carlo algorithm is used to find the best roots values [61–63]. The steps to follow are similar to those described in 2.3. Considering that only one integrator is added, the process must be simulated N times (equation (10)), using random values of r1 and r2 (between zero and one). Then, kσ and lσ are obtained with(37)kσ=r1r2,lσ=−r1−r2+r1r2+1TS.For each simulation, the cost function (equation (12)) is calculated. Finally, the parameters corresponding to the simulation with minimal total error are selected.When more than one integrator are added, the process is repeated, but now the characteristic equation order increases, and consequently, the expressions for parameters calculation change. For example, using two integrators, the parameters are calculated as follows:(38)kσ=r1r2r3,l1σ=2r1r2r3+1−r1r2−r1r3−r2r3TS,l2σ=−r1−r2−r3−r1r2r3+1+r1r2+r1r3+r2r3TS2.Table3 shows the results of the three controllers tuning, the original controller [51] (C1), a controller with one integrator (C2), and another with two integrators (C3).Table 3
Optimal controller parameters determined with Monte Carlo algorithm.
Variable
C1
C2
C3
k
X
0.4946
0.9289
0.7264
k
S1
0.3009
0.5411
0.3824
k
S2
0.1126
0.3110
0.2918
k
P1
0.4219
0.9953
0.9697
k
P2
0.2876
0.9694
0.8557
L
1X
—
5.8101
5.0548
l
1S1
—
0.2746
0.1422
L
1S2
—
1.8725
1.7765
l
1P1
—
3.0267
2.9858
L
1P2
—
2.9409
3.9991
L
2X
—
—
8.8199
l
2S1
—
—
0.6212
L
2S2
—
—
0.0003
l
2P1
—
—
1.0641
L
2P2
—
—
254.47


## 3.1. Controller Design under Uncertainty

In this section, the uncertainties effect in the tracking error is evaluated. With the aim of representing the additive uncertainty effect in the process mathematical model, the termsEσ,n are added in equation (7):(23)−Xn/VnλSf−S1,n/Vn1−λSf−S2,n/Vn−P1,n/Vn−P2,n/VnUn=ΔX−µ1+µ2XnΔS1+qS1/P1YP1/S1+qS1/P2YP2/S1XnΔS2+qS2/P1YP1/S2+qS2/P2YP2/S2XnΔP1−qS1/P1+qS2/P1XnΔP2−qS1/P2+qS2/P2Xn−1TSEX,nES1,nES2,nEP1,nEP2,n︸En.The additive uncertainty (En) models several kinds of mismatches as well as perturbed systems. En might depend on the states variables and the system input. Moreover, considering a real plant like zn+1 = f (zn,un), 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 [57, 58], then Eσ,n can be modelled as a bounded uncertainty [59, 60].Taking into account this uncertainty and following the same procedure of Section2.3, the next final expression is obtained:(24)eX,n+1eS1,n+1eS2,n+1eP1,n+1eP2,n+1=kX00000kS100000kS200000kP100000kP2eX,neS1,neS2,neP1,neP2,n︸L+TSXndμ1S1,P1,ndS1Sθ00−dqS1/P1S1,P1,ndS1Sφ−dqS1/P2S1,P2,ndS1SαeS1,n︸NL−EX,nES1,nES2,nEP1,nEP2,n︸En.Comparing (24) with (22), the error nonconvergence to zero is noticed due to En presence. Therefore, the next step is En approximation to reduce their effect in the tracking error, achieving the convergence to zero.


## 3.2. Integral Action

Here, a novel approach is presented to reduce the effect ofEn in the tracking error by adding some integrators in the control action calculus equation. The advantage is that new terms do not add extra mathematical complexity and maintain the simplicity of the scheme presented by Fernández et al. [51].As it was explained,En represents the additive uncertainty, typical of mathematical model formulation errors and the influence of external perturbation. Those uncertainties can not be measured, that is to say, they are unknown, but they can be assumed as a polynomial so that all kinds of possible uncertainties can be considered, which can be constant, linear, quadratic, or other. In this way, by increasing the order of the polynomial, the uncertainty approximation is improved and there is a tendency to reduce the tracking error more and more. This statement is mathematically demonstrated below.E
n differences are defined as(25)δEn=En+1−En,δ2En=δδEn=δEn+1−En=En+2−2En+1+En,δqEn=δδq−1En.


### 3.2.1. Constant Uncertainty

IfEn is assumed as a constant, then δEn = 0. Denoting generally the error as eσ,n,(26)eσ,n=σref,n−σn.An integrator could force the tracking error convergence to zero. This integrator is defined as follows:(27)Uσ,n+1=Uσ,n+∫nTSn+1TSeσ,ntdt≅Uσ,n+eσ,nTS.Consequently, equation (4) can be rewritten as(28)σn+1=σrefn+1−kσσrefn−σn+lσUσ,n+1,where kσ and lσ are the proportional and integral controller parameters, respectively.Demonstration of error convergence to zero when uncertainties are considered is given below.Expressing equation (24) in the following summarized way,(29)eσ,n+1=kσeσ,n+NLn+En.Adding the integral,(30)eσ,n+1=kσeσ,n−lσUσ,n+1+NLn+En.Replacing equation (27) in (30),(31)eσ,n+1=kσeσ,n−lσUσ,n+eσ,nT0+NLn+En.Then, combining equation (30) with (31),(32)eσ,n+2=kσeσ,n+1−lσUσ,n+1+eσ,n+1T0+NLn+1+En+1.Otherwise,(33)eσ,n+2=kσeσ,n+1−NLn+En+kσeσ,n−eσ,n+1−lσeσ,n+1T0+NLn+1+En+1,eσ,n+2=eσ,n+1kσ−lσT0+1−kσeσ,n−NLn+NLn+1−En+En+1︸δEn=0.Remark. If En = cte, δEn = 0. Looking equation (33), it is evidenced that the tracking error tends to zero with the progress of the process.


### 3.2.2. Linear Uncertainty

IfEn is supposed as a linear function, then δ2En = 0. Following the same reasoning as in 3.2.1, now two integrators are needed to compensate the added uncertainty. Those integrators are(34)Uσ1,n+1=Uσ1,n+∫nTSn+1TSeσntdt≅Uσ1,n+eσnTS,Uσ2,n+1=Uσ2,n+∫nTSn+1TSUσ1,n+1tdt≅Uσ2,n+Uσ1,n+1TS.Accordingly, equation (4) can be rewritten as(35)σn+1=σref,n+1−kσσref,n−σn+lσ1Uσ1,n+1+lσ2Uσ2,n+1,where kσ and lσ1 and lσ2 are the proportional and both integral controller parameters, respectively.After the same steps of3.2.1 demonstration, the error convergence to zero is achieved.


### 3.2.3. Polynomial Uncertainty

In this case,En is supposed as an m-order polynomial function. So, if q > m, δqEn = 0. Using the same steps before mentioned, it should be added m integrators in order to compensate the uncertainties considered.Figure3 shows the block diagram of the new control system.Figure 3
New control system block diagram.


## 3.2.1. Constant Uncertainty

IfEn is assumed as a constant, then δEn = 0. Denoting generally the error as eσ,n,(26)eσ,n=σref,n−σn.An integrator could force the tracking error convergence to zero. This integrator is defined as follows:(27)Uσ,n+1=Uσ,n+∫nTSn+1TSeσ,ntdt≅Uσ,n+eσ,nTS.Consequently, equation (4) can be rewritten as(28)σn+1=σrefn+1−kσσrefn−σn+lσUσ,n+1,where kσ and lσ are the proportional and integral controller parameters, respectively.Demonstration of error convergence to zero when uncertainties are considered is given below.Expressing equation (24) in the following summarized way,(29)eσ,n+1=kσeσ,n+NLn+En.Adding the integral,(30)eσ,n+1=kσeσ,n−lσUσ,n+1+NLn+En.Replacing equation (27) in (30),(31)eσ,n+1=kσeσ,n−lσUσ,n+eσ,nT0+NLn+En.Then, combining equation (30) with (31),(32)eσ,n+2=kσeσ,n+1−lσUσ,n+1+eσ,n+1T0+NLn+1+En+1.Otherwise,(33)eσ,n+2=kσeσ,n+1−NLn+En+kσeσ,n−eσ,n+1−lσeσ,n+1T0+NLn+1+En+1,eσ,n+2=eσ,n+1kσ−lσT0+1−kσeσ,n−NLn+NLn+1−En+En+1︸δEn=0.Remark. If En = cte, δEn = 0. Looking equation (33), it is evidenced that the tracking error tends to zero with the progress of the process.


## 3.2.2. Linear Uncertainty

IfEn is supposed as a linear function, then δ2En = 0. Following the same reasoning as in 3.2.1, now two integrators are needed to compensate the added uncertainty. Those integrators are(34)Uσ1,n+1=Uσ1,n+∫nTSn+1TSeσntdt≅Uσ1,n+eσnTS,Uσ2,n+1=Uσ2,n+∫nTSn+1TSUσ1,n+1tdt≅Uσ2,n+Uσ1,n+1TS.Accordingly, equation (4) can be rewritten as(35)σn+1=σref,n+1−kσσref,n−σn+lσ1Uσ1,n+1+lσ2Uσ2,n+1,where kσ and lσ1 and lσ2 are the proportional and both integral controller parameters, respectively.After the same steps of3.2.1 demonstration, the error convergence to zero is achieved.


## 3.2.3. Polynomial Uncertainty

In this case,En is supposed as an m-order polynomial function. So, if q > m, δqEn = 0. Using the same steps before mentioned, it should be added m integrators in order to compensate the uncertainties considered.Figure3 shows the block diagram of the new control system.Figure 3
New control system block diagram.


## 3.3. Controller Parameters Selection

For the error to decrease progressively and remain at low values, the linear term of equation (24), L, must tend to zero. In Fernández et al.’s study [51], this state was possible choosing kσ values between zero and one. However, when integrators are considered, two types of parameters are involved, proportional and integral, what changes the parameters selection condition to achieve the error convergence to zero.In this way, starting from the characteristic equation of the linear term of equation (33),(36)eσ,n+2+−kσ+lσTS−1︸a1eσ,n+1+kσ︸a2eσ,n=−NLn+NLn+1,r2+a1r+a2=0.If the characteristic equation roots (r1 and r2) are between zero and one, it is guaranteed that the linear term tends to zero when n tends to infinity. Hence, the root values must be determined to calculate parameters.As many authors recommend, the Monte Carlo algorithm is used to find the best roots values [61–63]. The steps to follow are similar to those described in 2.3. Considering that only one integrator is added, the process must be simulated N times (equation (10)), using random values of r1 and r2 (between zero and one). Then, kσ and lσ are obtained with(37)kσ=r1r2,lσ=−r1−r2+r1r2+1TS.For each simulation, the cost function (equation (12)) is calculated. Finally, the parameters corresponding to the simulation with minimal total error are selected.When more than one integrator are added, the process is repeated, but now the characteristic equation order increases, and consequently, the expressions for parameters calculation change. For example, using two integrators, the parameters are calculated as follows:(38)kσ=r1r2r3,l1σ=2r1r2r3+1−r1r2−r1r3−r2r3TS,l2σ=−r1−r2−r3−r1r2r3+1+r1r2+r1r3+r2r3TS2.Table3 shows the results of the three controllers tuning, the original controller [51] (C1), a controller with one integrator (C2), and another with two integrators (C3).Table 3
Optimal controller parameters determined with Monte Carlo algorithm.
Variable
C1
C2
C3
k
X
0.4946
0.9289
0.7264
k
S1
0.3009
0.5411
0.3824
k
S2
0.1126
0.3110
0.2918
k
P1
0.4219
0.9953
0.9697
k
P2
0.2876
0.9694
0.8557
L
1X
—
5.8101
5.0548
l
1S1
—
0.2746
0.1422
L
1S2
—
1.8725
1.7765
l
1P1
—
3.0267
2.9858
L
1P2
—
2.9409
3.9991
L
2X
—
—
8.8199
l
2S1
—
—
0.6212
L
2S2
—
—
0.0003
l
2P1
—
—
1.0641
L
2P2
—
—
254.47


## 4. Results and Discussion

The next section demonstrates the effectiveness of the improvement proposed through four important tests. In the first one, the operation of the three controllers is evaluated with a step perturbation in the control action. In the second one, a step and a ramp disturbance are added in the control action. In the third, the controllers are evaluated under parametric uncertainty using the Monte Carlo Algorithm. Finally, the controllers are tested against all earlier mentioned perturbations together. In each situation, the better performance of controllers C2 and C3 is evidenced. Furthermore, a comparison with a typical industry controller is shown.


### 4.1. Simulation Adding a Step Perturbation in the Control Action

In this test, a hypothetical situation that may produce an unexpected variation in the production is simulated. According to this, a −30% step perturbation in the bioreactor feed rate is added to evaluate the response of the controllers. Table4 summarizes simulations results. Figure 4 shows 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, it is reduced by 89.91% just adding one integrator and by 94.55% with two integrators.Table 4
Performance with a step perturbation in the control action.
C1
C2
C3
Total Error
1.6151
0.1630
0.0881Figure 4
(a) Controllers response to a step perturbation addition in the control action. (b) Total Error comparison. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)


### 4.2. Simulation Adding Step and Ramp Perturbations in the Control Action

In this test, a ramp disturbance is added to the step perturbation previously presented. Under those conditions, controllers C1, C2, and C3 are evaluated. In Table5, the Total Errors obtained in the simulations are exposed. Figure 5 shows how the feed rate varies in each case. Besides, the decrease of the Total Error is represented in the bar graph, demonstrating the improvement of the controller when one (91.4%) and two (97.1%) integrators are incorporated.Table 5
Performance with step and ramp perturbations in the control action.
C1
C2
C3
Total Error
6.8438
0.5887
0.2032Figure 5
(a) Controllers response to a step and ramp perturbation addition in the control action. (b) Total Error comparison. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)


### 4.3. Parametric Uncertainties

In all bioprocess, model parameters may change unexpectedly [64]. If this variation is not foreseen, there could be significant instabilities in the process [26]. Therefore, to achieve good performance and efficiency, bioreactors require advanced regulation procedures [29].As Monte Carlo Randomized Algorithm has been used for uncertainty quantification in many applications [65, 66], it is used for this test, following 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 [67–69].In 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 [70]. Hunag et al. [50] specified, in Table 1 of their work, the parameters’ confidence interval for the ethanol process under study. For this test, the information of the mentioned table is used for control purposes. Using C1, 1000 simulations are carried out. In each simulation, all the system parameters are randomly changed by ±10% of their original range value (Table 1 of Hunag et al.’s study [50]). Then, the Total Error is calculated with equation (12), and the system parameters for the worst situation are selected. Then, C2 and C3 are tested with those system parameters. The results are then compared considering the Total Error (Table 6). Note how, even considering the most problematic situation of the system, the Total Error decreases for each controller, and note the performance improvement for the controllers C2 (51.27%) and C3 (56.56%) over C1. Moreover, Figure 6 illustrates the accumulated error for each simulation and evidences the improvement with a bar graph.Table 6
Performance with system parametric uncertainty.
C1
C2
C3
Total Error
0.4820
0.2349
0.2094Figure 6
(a) Accumulated Error under parametric uncertainty. (b) Total Error comparison. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)Considering Theorem1 of Tempo and Ishii [54], as the tracking error (equation (11)) for the 1000 simulations remain bounded, the operation of the controllers C2 and C3 will be satisfactory with 99% of probability while the parameters vary within a range of ±10%.


### 4.4. Simulation Adding a Step and a Ramp Perturbation in the Control Action and Parametric Uncertainty

This last test aims to demonstrate how the controllers can fix an important perturbation. In this way, the perturbations described in4.1, 4.2, and 4.3 are simultaneously considered. Table 7 compares the resulting Total Errors. Also, Figure 7 shows the feed rate evolution in the three simulations, the accumulated error obtained with each controller, and a visible comparison of the performance improvement by a bar graph, when controllers C2 and C3 are used. In this case, the performance error can be reduced by 87.87% just adding one integrator and by 98.07% with two integrators.Table 7
Performance with a simultaneous step and ramp perturbation in the control action under parametric uncertainties.
C1
C2
C3
Total Error
6.8371
0.8297
0.1322Figure 7
Controllers response comparison to a simultaneous step and ramp perturbation in the control action under parametric uncertainties. (a) Control action. (b) Accumulated Error. (c) Total Error. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)
(c)To summarize, the addition of integrators in the controller design has advantages, such as performance improvement, no significant mathematical complexity increase, and additive uncertainty consideration, although it has a disadvantage as integrators are added, the equations order increases, making mathematical development less simpler.


### 4.5. Comparison with a Typical Controller

For their simple structure and easy parameter adjustment, PID or PI controllers are still the most widely used in factories [71], although their weakness with tracking variable set points is well known [72–74]. Figure 8 shows a performance comparison between the proposed controllers and a traditional PI in the same ethanol production system. Here, a step and a ramp perturbation in the control action are added. The best PI parameters were selected with a Monte Carlo algorithm. PI accumulated error is much higher than that of C1, due to the increase in error over time.Figure 8
Total Error comparison when (a) a step and (b) a ramp perturbation in the control action are added. (PI) Proportional-integral controller; (C1) original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)It is important to note that a rigorous comparison between C1 and PI is detailed in [51].


## 4.1. Simulation Adding a Step Perturbation in the Control Action

In this test, a hypothetical situation that may produce an unexpected variation in the production is simulated. According to this, a −30% step perturbation in the bioreactor feed rate is added to evaluate the response of the controllers. Table4 summarizes simulations results. Figure 4 shows 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, it is reduced by 89.91% just adding one integrator and by 94.55% with two integrators.Table 4
Performance with a step perturbation in the control action.
C1
C2
C3
Total Error
1.6151
0.1630
0.0881Figure 4
(a) Controllers response to a step perturbation addition in the control action. (b) Total Error comparison. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)


## 4.2. Simulation Adding Step and Ramp Perturbations in the Control Action

In this test, a ramp disturbance is added to the step perturbation previously presented. Under those conditions, controllers C1, C2, and C3 are evaluated. In Table5, the Total Errors obtained in the simulations are exposed. Figure 5 shows how the feed rate varies in each case. Besides, the decrease of the Total Error is represented in the bar graph, demonstrating the improvement of the controller when one (91.4%) and two (97.1%) integrators are incorporated.Table 5
Performance with step and ramp perturbations in the control action.
C1
C2
C3
Total Error
6.8438
0.5887
0.2032Figure 5
(a) Controllers response to a step and ramp perturbation addition in the control action. (b) Total Error comparison. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)


## 4.3. Parametric Uncertainties

In all bioprocess, model parameters may change unexpectedly [64]. If this variation is not foreseen, there could be significant instabilities in the process [26]. Therefore, to achieve good performance and efficiency, bioreactors require advanced regulation procedures [29].As Monte Carlo Randomized Algorithm has been used for uncertainty quantification in many applications [65, 66], it is used for this test, following 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 [67–69].In 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 [70]. Hunag et al. [50] specified, in Table 1 of their work, the parameters’ confidence interval for the ethanol process under study. For this test, the information of the mentioned table is used for control purposes. Using C1, 1000 simulations are carried out. In each simulation, all the system parameters are randomly changed by ±10% of their original range value (Table 1 of Hunag et al.’s study [50]). Then, the Total Error is calculated with equation (12), and the system parameters for the worst situation are selected. Then, C2 and C3 are tested with those system parameters. The results are then compared considering the Total Error (Table 6). Note how, even considering the most problematic situation of the system, the Total Error decreases for each controller, and note the performance improvement for the controllers C2 (51.27%) and C3 (56.56%) over C1. Moreover, Figure 6 illustrates the accumulated error for each simulation and evidences the improvement with a bar graph.Table 6
Performance with system parametric uncertainty.
C1
C2
C3
Total Error
0.4820
0.2349
0.2094Figure 6
(a) Accumulated Error under parametric uncertainty. (b) Total Error comparison. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)Considering Theorem1 of Tempo and Ishii [54], as the tracking error (equation (11)) for the 1000 simulations remain bounded, the operation of the controllers C2 and C3 will be satisfactory with 99% of probability while the parameters vary within a range of ±10%.


## 4.4. Simulation Adding a Step and a Ramp Perturbation in the Control Action and Parametric Uncertainty

This last test aims to demonstrate how the controllers can fix an important perturbation. In this way, the perturbations described in4.1, 4.2, and 4.3 are simultaneously considered. Table 7 compares the resulting Total Errors. Also, Figure 7 shows the feed rate evolution in the three simulations, the accumulated error obtained with each controller, and a visible comparison of the performance improvement by a bar graph, when controllers C2 and C3 are used. In this case, the performance error can be reduced by 87.87% just adding one integrator and by 98.07% with two integrators.Table 7
Performance with a simultaneous step and ramp perturbation in the control action under parametric uncertainties.
C1
C2
C3
Total Error
6.8371
0.8297
0.1322Figure 7
Controllers response comparison to a simultaneous step and ramp perturbation in the control action under parametric uncertainties. (a) Control action. (b) Accumulated Error. (c) Total Error. (C1) Original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)
(c)To summarize, the addition of integrators in the controller design has advantages, such as performance improvement, no significant mathematical complexity increase, and additive uncertainty consideration, although it has a disadvantage as integrators are added, the equations order increases, making mathematical development less simpler.


## 4.5. Comparison with a Typical Controller

For their simple structure and easy parameter adjustment, PID or PI controllers are still the most widely used in factories [71], although their weakness with tracking variable set points is well known [72–74]. Figure 8 shows a performance comparison between the proposed controllers and a traditional PI in the same ethanol production system. Here, a step and a ramp perturbation in the control action are added. The best PI parameters were selected with a Monte Carlo algorithm. PI accumulated error is much higher than that of C1, due to the increase in error over time.Figure 8
Total Error comparison when (a) a step and (b) a ramp perturbation in the control action are added. (PI) Proportional-integral controller; (C1) original controller [51]; (C2) controller with one integrator; (C3) controller with two integrators.
(a)
(b)It is important to note that a rigorous comparison between C1 and PI is detailed in [51].


## 5. Conclusions

A previously presented control technique for tracking optimal profiles in nonlinear and multivariable systems was improved. The proposed methodology allows tracking predefined concentration profiles, even in the presence of typical bioprocess uncertainties like unknown nonlinearities or parameters uncertainties. To take into account the mentioned uncertainties, one or more integrators of the tracking error are added in the controller design. The main contribution of this paper is to reduce additive uncertainties influence on the tracking error without increasing the controller complexity. Note that C2 and C3 preserve the original controller advantages (simplicity, fast design, etc.) but manage to reduce uncertainties effect.The Monte Carlo Randomized Algorithm is used to tune the controllers and carry out the test with parametric uncertainty. The different tests carried out show the effectiveness of the proposed procedure, reaching an improvement of 98% in the worst case considered.Comparing with other previous strategies that deal with analogous uncertain control problems, the proposed controller presents the advantage of avoiding the stochastic modelling 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: 8263690-2020-01-22.xml*