STATE ESTIMATION AND TRAJECTORY TRACKING CONTROL FOR A NONLINEAR AND MULTIVARIABLE BIOETHANOL PRODUCTION SYSTEM

ABSTRACT

In this paper a controller is proposed based on linear algebra for a fed-batch bioethanol production process. It involves finding feed rate profiles (control actions obtained as a solution of a linear equations system) in order to make the system follow predefined concentration profiles. A neural network states estimation is designed in order to know those variables that cannot be measured. The controller is tuned using a Monte Carlo experiment for which a cost function that penalizes tracking errors is defined. Moreover, several tests (adding parametric uncertainty and perturbations in the control action) are carried out so as to evaluate the controller performance. A comparison with another controller is made. The demonstration of the error convergence, as well as the stability analysis of the neural network, are included.

Key words:
Fed-batch bioprocess; Nonlinear and multivariable system; Profiles tracking control; Numerical methods/linear algebra; State estimation

INTRODUCTION

In recent years, the bioprocess industry, so called white biotechnology (Heux et al., 2015), has gained an important position. This is because it has an important role in the production of high-added value products; such as recombinant proteins, vaccines and antibiotics in the pharmaceutical industry, or beer, wine, yeast in the manufacturing of agro-food goods, or biogas and compost in the treatment of urban and industrial solid organic wastes and wastewater (Mangesh and Jana, 2008), between others.

Alcoholic fermentation is an ancient practice, commonly used in the production of alcoholic beverages such as beer, wine, cider, sake and distilled drinks (Wood, 2012). Moreover, research on using ethanol as an alternative fuel has gained tremendous attention all over the world since the petroleum crisis in the 1970s (Ajbar and Ali, 2017). It is noteworthy that all these processes are carried out in bioreactors.

There are many operating modes for bioreactors: batch, fed-batch and continuous. Among the different modes, the fed-batch is preferred because of the operational flexibility that it provides (Mangesh and Jana, 2008). It consists of changing the feed rate of nutrients, inhibitors, catalysts or inducers along the process, whereas cells and products remain in the fermenter until the operation ends (Hecklau et al., 2016). Moreover, this cultivation technique has several advantages over other operation modes (Liu et al., 2013). Also, it avoids formation of unwanted products and microorganism inhibition caused by overfeeding, and prevents the microorganism starvation induced by underfeeding (Ochoa, 2016). Ethanol production in a fed-batch fermenter is common in white biotechnology (Mangesh and Jana, 2008).

One of the main benefits of fed-batch processes is that the concentration of substrate in the cultivation medium can be externally regulated with a suitable feed rate profile (Liu et al., 2013), in order to obtain better production yields (Jin et al., 2014). For this reason, many efforts have focused on bioprocess optimization and control so as to minimize the production costs while increasing the yield and productivity (Ochoa, 2016, Pantano et al., 2017a). Nevertheless, that is an arduous task because usually microorganisms have a complex dynamic behavior (nonlinear and sometimes unstable) that leads to strong modeling approximations; furthermore, there are other complications like the presence of numerous external disturbances or the difficulty of measuring representative variables. These problems avoid the possibility of using PI, PID or other classic industrial controller, which means that the development of an specific control algorithm is necessary for each bioprocess (De Battista et al., 2012).

The worldwide market requirements for high standard products, as well as safe and environmentally friendly processes, force chemical industries to look for optimization methods and controllers that can afford nonlinearities and transient behavior of chemical and biochemical processes (Fujiki et al., 2009). On the one hand, many scientists have established optimization methods to find the best feed rate profile for different processes involving fed-batch fermentations (for some examples see: Ye et al., 2014; Kookos, 2004; del Rio-Chanona et al., 2016 and Dai et al., 2014, as this issue is beyond the aim of this paper). However, while finding the optimal feed policy of a determined process is not an easy task, once it is known a bigger challenge is presented: trying to track the optimal profiles and obtain repeatability from one batch to another. Recent papers have been focused on solving this problem, such as (Chang et al., 2016), (Jin et al., 2014), (Fu and Chai, 2007), (Lehouche et al., 2012), (Hofmann et al., 2017), and (Mohammadzaheri and Chen, 2011, Hofmann et al., 2017).

Some advanced control techniques for fermentation processes like: on-line adaptive control (Guay et al., 2004), optimal control (Logist et al., 2010), fuzzy control (Karakuzu et al., 2006), neural networks (Imtiaz et al., 2013) and predictive control based on nonlinear model (Preuβ et al., 2003) have gained popularity because of their strong capability in dealing with process non-linearity, dynamics and optimization. This last (NMPC), although it has been successfully applied in numerous practical applications, encounters a lot of limitations such as: difficulties in building accurate dynamic models, complexity of online implementation, the insufficient accuracy of on-line solutions and the computational time required to find the solution, limiting its applications to bioprocesses (Jin et al., 2014). Therefore, probably the application of the other techniques can be difficult or even infeasible for bioprocess. For all this reasons, the PID controller is still the most widely used in factories (Imtiaz et al., 2014). PID controllers are chosen mainly for their simple structure and ease of adjustment of their parameters. On the other hand, a major challenge is the control set point variations, avoiding oscillations and delays, trying to reduce tracking errors at the same time. Therefore, although PID controllers are sufficient to solve the control problem of many applications in industry, for bioprocesses they do not give good results with traditional tuning procedures. Few works in this field have claimed satisfactory results tuning the parameters based on alternative approaches such as: neural networks (Andrášik et al., 2004) and metaheuristic algorithms (Roeva and Slavov, 2012) to deal with the inherent complex behavior of bioprocess.

The aim of this paper is to present a controller for application in bioethanol production. The system is represented with a nonlinear and multivariable model. The controller consists of obtaining a feed rate profile (control action) in order to track optimal concentration profiles (desired variables). The main innovation is the simultaneously excellent tracking of four time varying profiles with only one control action. The controller has a fast and easy design because it allows obtaining the control action as a solution of linear equations, even though the original system model is nonlinear. Thus, only algebra is needed to understand and apply this methodology. As the controller structure comes from the mathematical model of the process, it can be implemented in many systems. Besides, it is versatile against different changes and disturbances in the process and system. This last assertion is confirmed with different tests and, in each one, it can be observed that the error remained at low levels (the mathematical demonstration of this fact is also presented in Appendix A). Additionally, the controller parameters are selected with a Monte Carlo Randomized algorithm and a comparison with a PID controller is shown.

In order to achieve the proposed objective, it is assumed that the bioethanol process is well represented by the mathematical model, the desired profiles are known and the control action can be obtained at each instant of time. Furthermore, the state variables must be able to be measured online, but as is known, the states of some systems are immeasurable. Accordingly, the estimation of those states plays a significant role in control design and system identification (Xu et al., 2017). Therefore, a neural network states estimation is designed. Fig. 1 shows a flowchart of the proposed procedure.

Figure 1
Flowchart of the developed technique.

The paper is ordered as follows. In the next section the process and the mathematical model that governs it are described. The controller design is then described; this includes the determination of the optimal parameters of the controller, which are obtained with a Monte Carlo experiment. Section 4 shows the neural network states estimator. The subsequent section is organized in three different subsections; first, a simulation of the process under normal operation conditions is carried out; after that, as a demonstration of the good performance, two tests are executed: addition of parametric uncertainty and disturbances in the control action. Then, the subsequent section shows a comparison with a typical PID controller used in the industry. Finally, the conclusions are presented.

PROCESS DESCRIPTION

In this study, an ethanol fermentation process carried out in a fed-batch reactor is considered. The microorganism used is Saccharomyces diastaticus. This yeast has high ethanol tolerance and can utilize glucose and fructose to produce ethanol and glycerol simultaneously (Hunag et al., 2012). The mathematical modeling of this system was proposed by Hunag and coworkers (Hunag et al., 2012). They used a batch fermenter to generate time series data, which were then applied to estimate the kinetic model parameters. After that, they applied them to a fed-batch process to determine an optimal control policy. This system has one input, namely, the feed rate, and various outputs (whose profiles must follow determined variation over time), i.e., the cells, ethanol, glycerol, glucose and fructose concentrations inside the reactor. The reactor’s feed consists of a 50-50 mixture of glucose and fructose. The developed mathematical model is:

$$\left\{\begin{array}{l}\dot{X}\left(t\right)=\left({\mu }_1+{\mu }_2\right)X-\frac{U}{V}X\\ {\dot{S}}_1\left(t\right)=-\left(\frac{q_{S_1/P_1}}{Y_{P_1/S_1}}+\frac{q_{S_1/P_2}}{Y_{P_2/S_1}}\right)X+\frac{U}{V}\left(\lambda S_f-S_1\right)\\ {\dot{S}}_2\left(t\right)=-\left(\frac{q_{S_2/P_1}}{Y_{P_1/S_2}}+\frac{q_{S_2/P_2}}{Y_{P_2/S_2}}\right)X+\frac{U}{V}\left(\left(1-\lambda \right)S_f-S_2\right)\\ {\dot{P}}_1\left(t\right)=\left(q_{S_1/P_1}+q_{S_2/P_1}\right)X-\frac{U}{V}{P}_1\\ {\dot{P}}_2\left(t\right)=\left(q_{S_1/P_2}+q_{S_2/P_2}\right)X-\frac{U}{V}{P}_2\end{array}\right.$$ (1)

where:

$$\begin{array}{l}\dot{V}\left(t\right)=U{q}_{S_1/P_1}=\frac{{\nu }_{S_1{P}_1}{S}_1}{K_{S_1{P}_1}+S_1}\frac{k_{S_1{P}_1}}{k_{S_1{P}_1}+P_1}\\ {\mu }_1=\frac{{\mu }_{m_1}{S}_1}{K_{S_1}+S_1+S_1^2/K_{S_{1I}}}\frac{K_{P_1}}{K_{P_1}+P_1+P_1^2/K_{P_{1I}}}{q}_{S_2/P_2}=\frac{{\nu }_{S_2{P}_2}{S}_2}{K_{S_2{P}_2}+S_2}\frac{k_{S_2{P}_2}}{k_{S_2{P}_2}+P_2}\\ {\mu }_2=\frac{{\mu }_{m_2}{S}_2}{K_{S_2}+S_2+S_2^2/K_{S_{2I}}}\frac{K_{P_2}}{K_{P_2}+P_1+P_1^2/K_{P_{2I}}}{q}_{S_1/P_2}=\frac{{\nu }_{S_1{P}_2}{S}_1}{K_{S_1{P}_2}+S_1}\frac{k_{S_1{P}_2}}{k_{S_1{P}_2}+P_2}\\ q_{S_2/P_1}=\frac{{\nu }_{S_2{P}_1}{S}_2}{K_{S_2{P}_1}+S_2}\frac{k_{S_2{P}_1}}{k_{S_2{P}_1}+P_1}\end{array}$$ (2)

Here, the state variables are: biomass (X), glucose (S ₁ ), fructose (S ₂ ), ethanol (P ₁ ) and glycerol (P ₂ ) concentrations inside the reactor. The manipulated variable is the substrate feed rate (U). Furthermore, µ ₁ and µ ₂ represent the specific yeast cells growth rate, q _S1/P1 and q _S2/P1 the specific ethanol production rate, and q _S1/P2 and q _S2/P2 the specific glycerol production rate, in all cases from glucose and fructose, respectively. V is the culture volume. All these relations describe the metabolic activity of microorganisms.

The variables initial values are shown in Table 1. The parameter nomenclature, description and values are in Table 2.

Table 1
Initial variable values for ethanol fermentation (Hunag et al., 2012).

Variable	Initial value
X (g/L)	1.5
P1 (g/L	5.3
P2 (g/L)	0.0001
S1 (g/L)	8.6
S2 (g/L)	8.6
V (L)	1.35
U (L/h)	0.9325

Table 2
Parameters nomenclature, description and values (Hunag et al., 2012).

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
YP1/S1	Yield coefficient for ethanol from glucose	0.5085
YP2/S1	Yield coefficient for glycerol from glucose	0.5331
YP1/S2	Yield coefficient for ethanol from fructose	0.5098
YP2/S2	Yield coefficient for glycerol from fructose	0.4462
KS1	Saturation coefficient for cell growth on glucose (g/L)	159.75
KS1I	Inhibition coefficient for cell growth on glucose (g/L)	94.233
KP1	Saturation coefficient for cell growth on ethanol (g/L)	238.39
KP1I	Inhibition coefficient for cell growth on ethanol (g/l)	2.7378
KS2	Saturation coefficient for cell growth on fructose (g/L)	0.0726
KS2I	Inhibition coefficient for cell growth on fructose (g/L)	9.0048
KP2	Saturation coefficient for cell growth on glycerol (g/L)	35.958
KP2I	Inhibition coefficient for cell growth on glycerol (g/L)	9.9722
KS1P1	Saturation coefficient for ethanol production on glucose (g/L)	1.3409
kS1P1	Inhibition coefficient for ethanol production on glucose (g/L)	18.612
KS2P1	Saturation coefficient for ethanol production on fructose (g/L)	0.9129
kS2P1	Inhibition coefficient for ethanol production on fructose (g/L)	1000
KS1P2	Saturation coefficient for glycerol production on glucose (g/L)	6.7116
kS1P2	Inhibition coefficient for glycerol production on glucose (g/L)	0.5863
KS2P2	Saturation coefficient for glycerol production on fructose (g/L)	0.4310
kS2P2	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
Sf	Sugar total feed concentration (g/L)	300

CONTROLLER DESIGN

Controller Structure

Generally, the structure of a controller is constructed based on a mathematical model. This kind of model gives a scheme of the process from inputs to outputs. Obviously, its quality depends on how closely those outputs match the actual process. Furthermore, it is important to highlight that a model that accurately replicates a real process will never be created (Zhou et al., 1996). For this reason, a controller that allows monitoring predefined profiles (obtained with a specific model and determined condition) in the presence of different perturbations, with minimal error, is necessary.

Below, the main contribution of this work, a control methodology based on linear algebra for tracking predefined profiles in a bioreactor is presented. For its development the mathematical model of the process (shown in Section 2) and the reference profiles that are wanted to be followed by the system are required. Those profiles are the concentrations of cells, ethanol, glycerol and fructose that are the variables of interest in this bioprocess. To determine them, an open-loop simulation of Eq. (1) was done using the initial conditions of Table 1, the parameters of Table 2 and the optimal feed rate policy determined in (Hunag et al., 2012). Fig. 2: illustrates the obtained reference concentration profiles, the culture volume and the corresponding feed flow rate.

Figure 2
Cells, ethanol, glycerol and fructose reference concentrations along the process. Reference feed flow rate.

To begin with the development of the technique, as the mathematical model is given as a system of ordinary differential equations, a numerical method is used to integrate them. Euler is applied in this study because it is simple and yields good results.

$$\left(\frac{d\sigma }{dt}\right)=\frac{{\sigma }_{n+1}-{\sigma }_n}{T_S}$$ (3)

In Eq. (3), σ symbolizes each state variable, σ _n is the current value of σ measured from the reactor (on-line), and σ _n+1 is the value of σ in the next measurement instant. T _S is the sampling time; for this study, its value is 0.1 h, adopted taking into account the recommendations of (Griffiths and Smith, 2006), ensuring that no important event is neglected, and that the error is minimized and compensated with the computational cost. The process lasts 15.7 h.

As the state variables in n+1 are unknown, it is necessary to find some way to approximate them. In this paper, the value at n+1 time for each variable is adopted by a function of the error at n time (see Eq. (4)). This allows expressing the unknown variables (σ _n+1 ) as a function of the known ones: σ _n , references in n+1 and in n, and a constant.

$$\underset{erro{r}_{n+1}}{\underbrace{{\sigma }_{ref,}{}_{n+1}-{\sigma }_{n+1}}}=k_{\sigma }(\underset{erro{r}_n}{\underbrace{{\sigma }_{ref,}{}_n-{\sigma }_n}})\to {\sigma }_{n+1}={\sigma }_{ref}{}_{,n+1}-k_{\sigma }({\sigma }_{ref}{}_{,n}-{\sigma }_n)$$ (4)

In this equation, σ _ref are the reference state variable values at the corresponding instant of time. k _σ is a constant that represents the controller parameter for the variable σ. There are five different k _σ for this model of ethanol production, one for each variable, k _X , k _S1 , k _S2 , k _P1 and k _P2 .

Substituting Eq. (4) in Eq. (3) gives the following expression that allows the approximation of the derivatives:

$$\left(\frac{d\sigma }{dt}\right)=\frac{\stackrel{{\sigma }_{n+1}}{\overbrace{\left[{\sigma }_{ref}{}_{,n+1}-k_{\sigma }({\sigma }_{ref}{}_{,n}-{\sigma }_n)\right]}}-{\sigma }_n}{T_S}$$ (5)

Substituting Eq. (5) in the mathematical model:

$$\left\{\begin{array}{l}\frac{\left(X_{ref,n+1}-k_X\left(X_{ref,n}-X_n\right)-X_n\right)}{T_S}=\left({\mu }_{1,n}+{\mu }_{2,n}\right)X_n-\frac{U_n}{V_n}{X}_n\\ \frac{\left(S_{1ref,n+1}-k_{S_1}\left(S_{1ref,n}-S_{1,n}\right)-S_{1,n}\right)}{T_S}=-\left(\frac{{q_{S_1/P_1,}}_n}{Y_{P_1/S_1}}+\frac{q_{S_1/P_2,n}}{Y_{P_2/S_1}}\right)X_n+\frac{U_n}{V_n}\left(\lambda S_f-S_{1,n}\right)\\ \frac{\left(S_{2ref,n+1}-k_{S_2}\left(S_{2ref,n}-S_{2,n}\right)-S_{2,n}\right)}{T_S}=-\left(\frac{q_{S_2/P_1,n}}{Y_{P_1/S_2}}+\frac{q_{S_2/P_{2,}n}}{Y_{P_2/S_2}}\right)X_n+\frac{U_n}{V_n}\left(\left(1-\lambda \right)S_f-S_{2,n}\right)\\ \frac{\left(P_{1ref,n+1}-k_{P_1}\left(P_{1ref,n}-P_{1n}\right)-P_{1,n}\right)}{T_S}=\left(q_{S_1/P_1,n}+q_{S_2/P_1,n}\right)X_n-\frac{U_n}{V_n}{P}_{1,n}\\ \frac{\left(P_{2ref,n+1}-k_{P_2}\left(P_{2ref,n}-P_{2,n}\right)-P_{2,n}\right)}{T_S}=\left(q_{S_1/P_2,n}+q_{S_2/P_2,n}\right)X_n-\frac{U_n}{V_n}{P}_{2,n}\end{array}\right.$$ (6)

With this new appreciation of Eq. (1) as a system of linear equations, a simple possibility for calculating the control action is available. Besides, it could be expressed in a matrix form by placing the state variables as a function of U. In this way, it is more straightforward to clear the control action:

$${\displaystyle \stackrel{A}{\overbrace{\left[\begin{array}{c}-X_n/V_n\\ (\lambda S_f-S_{1,n})/V_n\\ ((1-\lambda )S_f-S_{2,n})/V_n\\ -P_{1,n}/V_n\\ -P_{2,n}/V_n\end{array}\right]}}\stackrel{U}{\overbrace{U_n}}=\underset{b}{\underbrace{\left[\begin{array}{l}\frac{(X_{ref,n+1}-k_X(X_{ref,n}-X_n)-X_n)}{T_S}-({\mu }_{1,n}+{\mu }_{2,n})X_n\\ \frac{(S_{1ref,n+1}-k_{S_1}(S_{1ref,n}-S_{1,n})-S_{1,n})}{T_S}+(\frac{q_{S_1/P_1,n}}{Y_{P_1/S_1}}+\frac{q_{S_2/P_2,n}}{Y_{P_2/S_1}})X_n\\ \frac{\left(S_{2ref,n+1}-k_{S_2}\left(S_{2ref,n}-S_{2,n}\right)-S_{2,n}\right)}{T_S}+\left(\frac{q_{S_2/P_1,n}}{Y_{P_1/S_2}}+\frac{q_{S_2/P_2,n}}{Y_{P_2/S_2}}\right)X_n\\ \frac{\left(P_{1ref,n+1}-k_{P_1}\left(P_{1ref,n}-P_{1,n}\right)-P_{1,n}\right)}{T_S}-\left(q_{S_1/P_1,n}+q_{S_2/P_1,n}\right)X_n\\ \frac{\left(P_{2ref,n+1}-k_{P_2}\left(P_{2ref,n}-P_{2,n}\right)-P_{2,n}\right)}{T_S}-\left(q_{S_1/P_2,n}+q_{S_2/P_2,n}\right)X_n\end{array}\right]}}}$$ (7)

To simplify the mathematical expression of the problem, Eq. (7) is expressed generically as AU=b:

$$\left[\begin{array}{l}{a}_1\\ a_2\\ a_3\\ a_4\\ a_5\end{array}\right]U=\left[\begin{array}{l}{b}_1\\ b_2\\ b_3\\ b_4\\ b_5\end{array}\right]$$ (8)

To find U, the system must have an exact solution. To accomplish this, b has to be a linear combination of A columns (Strang, 2006), that is to say, A and b must be parallel. This condition can be satisfied in different ways; one of them is that the angle between A and b must be zero:

$$\cos (A,b)=\frac{<A,b>}{\left|\right|A\left|\right|.\left|\right|b\left|\right|}=1$$ (9)

Here, the operation between < > and ||.|| represent the inner product and the norm of the vectors in the R ⁿ space, respectively. Consequently, Eq. (9) can be expressed as:

$$(a_1{b}_1+a_2{b}_2+a_3{b}_3+a_4{b}_4+a_5{b}_5)-\sqrt{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}\sqrt{b_1^2+b_2^2+b_3^2+b_4^2+b_5^2}=0$$ (10)

For Eq. (10) to have solution, the “sacrificed variable” is defined, which is denoted by the subscript “ez”. To select it, it is essential to examine and interpret the role of each variable in the process. In a bioprocess, the substrate concentration, which can be adjusted by varying the supply flow rate, directly affects the rate of substrate consumption, growth rate of cells, and formation rates of products and byproducts (Öztürk et al., 2016). As can be appreciated in Eq. (1), the kinetic model for cell growth and ethanol formation considers the interactions between two types of sugars: glucose and fructose (Hunag et al., 2012). However, few articles have addressed this issue because the influence of glucose is more significant than that of fructose. Therefore, S ₁ is chosen as the sacrificed variable.

Then, by replacing S _1ref by S _1ez in Eq. (7) and applying the parallelism condition of Eq. (10), the sacrificed variable is readily to be found. Once S _1ez,n is calculated, the control action (U _n ) can be determined at any sampling time. To accomplish this, S _1ez is replaced in Eq. (7) and U _n is obtained by using least squares (Strang, 2006).

$$U={\left(A^{T}A\right)}^{-1}{A}^{T}b=\frac{a_1{b}_1+a_2{b}_2+a_3{b}_3+a_4{b}_4+a_5{b}_5}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}$$ (11)

Replacement with component values of each corresponding matrix:

$$\begin{array}{l}{U}_n=\frac{-X_n{V}_n\left(\frac{\left(X_{ref,n+1}-k_X\left(X_{ref,n}-X_n\right)-X_n\right)}{T_S}-\left({\mu }_{1,n}+{\mu }_{2,n}\right)X_n\right)}{-X_n^2+{\left(\lambda S_f-S_{1,n}\right)}^2+{\left(\left(1-\lambda \right)S_f-S_{2,n}\right)}^2-P_{1,n}^2-P_{2,n}^2}\\ +\frac{\left(\lambda S_f-S_{1,n}\right)V_n\left(\frac{\left(S_{1ez,n+1}-k_{S_1}\left(S_{1ez,n}-S_{1,n}\right)-S_{1,n}\right)}{T_S}+\left(\frac{{q_{S_1/P_1,}}_n}{Y_{P_1/S_1}}+\frac{q_{S_1/P_2,n}}{Y_{P_2/S_1}}\right)X_n\right)}{-X_n^2+{\left(\lambda S_f-S_{1,n}\right)}^2+{\left(\left(1-\lambda \right)S_f-S_{2,n}\right)}^2-P_{1,n}^2-P_{2,n}^2}\\ +\frac{\left(\left(1-\lambda \right)S_f-S_{2,n}\right)V_n\left(\frac{\left(S_{2ref,n+1}-k_{S_2}\left(S_{2ref,n}-S_{2,n}\right)-S_{2,n}\right)}{T_S}+\left(\frac{q_{S_2/P_1,n}}{Y_{P_1/S_2}}+\frac{q_{S_2/P_2,n}}{Y_{P_2/S_2}}\right)X_n\right)}{-X_n^2+{\left(\lambda S_f-S_{1,n}\right)}^2+{\left(\left(1-\lambda \right)S_f-S_{2,n}\right)}^2-P_{1,n}^2-P_{2,n}^2}\\ +\frac{-P_{1,n}{V}_n\left(\frac{\left(P_{1ref,n+1}-k_{P_1}\left(P_{1ref,n}-P_{1,n}\right)-P_{1,n}\right)}{T_S}-\left(q_{S_1/P_1,n}+q_{S_2/P_1,n}\right)X_n\right)}{-X_n^2+{\left(\lambda S_f-S_{1,n}\right)}^2+{\left(\left(1-\lambda \right)S_f-S_{2,n}\right)}^2-P_{1,n}^2-P_{2,n}^2}\\ +\frac{-P_{2,n}{V}_n\left(\frac{\left(P_{2ref,n+1}-k_{P_2}\left(P_{2ref,n}-P_{2,n}\right)-P_{2,n}\right)}{T_S}-\left(q_{S_1/P_2,n}+q_{S_2/P_2,n}\right)X_n\right)}{-X_n^2+{\left(\lambda S_f-S_{1,n}\right)}^2+{\left(\left(1-\lambda \right)S_f-S_{2,n}\right)}^2-P_{1,n}^2-P_{2,n}^2}\end{array}$$ (12)

Equation is the mathematical representation of the control action required to achieve the concentration profiles shown in Fig. 2.

Note that one of the advantages of this technique is that it can be applied in many systems, both single-input multiple-output and multiple-inputs multiple-outputs, if the information required is available.

Controller Parameters Selection.

Here, a new term is introduced, the “tracking error”, defined as follows:

$${\displaystyle \parallel e_n\parallel =\sqrt{\begin{array}{c}((X_{ref,n}-X_n)/max{X}_{ref}{)}^2+((P_{1ref,n}-P_{1,n})/max{P}_{1ref}{)}^2\\ +((P_{2ref,n}-P_{2,n})/max{P}_{2ref}{)}^2+((S_{2ref,n}-S_{2,n})/max{S}_{2ref}{)}^2\end{array}}}$$ (13)

where max σ _ref is the maximum value of the corresponding reference variable.

As was introduced aboce, the performance of the bioreactor is directly affected by the controller parameters (k _σ ). Those parameters take values between zero and one (0 < k _σ < 1), hence the tracking error tends to zero when n tends to infinity (go to Appendix A to see the demonstration).

The following procedure aims to find the best values for these parameters, between 0 and 1, such that the tracking error is minimized. So, viewing this as a problem of searching in an unknown environment, an appropriate algorithm must be chosen in order to solve it. Several variants of search problems have been proposed (e.g., Chrobak et al., 2008). Genetic Algorithms are widely used to solve complex optimization problems (Holland, 1975, Sadatsakkak et al., 2015, Asadi et al., 2014, Yu et al., 2015, Ismail et al., 2014, Alvarez et al., 2008); this method imitates the theory of biological evolution proposed by Darwin for the resolution of problems. Ant Colony Optimization is a probabilistic technique to solve computational problems that can be reduced to looking for the best paths in graphs; it is based on the behavior of ants when they are looking for a way between the colony and a food source (Chiha et al., 2012, Castillo et al., 2015, Omar et al., 2013). The Monte Carlo Algorithm gives statistical estimates of the required solution by performing random sampling of a random variable, whose mathematical expectation is the desired solution (Dimov et al., 2015).

In this test, the Monte Carlo algorithm is applied to tune the controller. This method is chosen because it has several advantages: i) there is low probability of producing an incorrect result (Motwani and Raghavan, 1995); ii) less computational complexity in comparison with other algorithms (Tempo and Ishii, 2007); iii) often, a slight modification of the algorithm for solving systems of linear algebraic equations allows one to solve other linear algebra problems such as matrix inversion and computing the extreme eigenvalues (Dimov et al., 2015); iv) high reliability and easy application (Cheein and Scaglia, 2014, de Oliveira et al., 2012, Mohammadi et al., 2014).

The experiment consists of simulating the process a number of times while using random values of k _σ . Then, the total error is calculated for each iteration. The k _σ values that make a minimum total error are selected.

To determine the number of simulations (N), Eq. (14) is used (Tempo and Ishii, 2007). Note that in order to limit the chance of a wrong answer, appropriate confidence (δ) and accuracy (ε) must be indicated.

$$N\ge \left[\frac{\log \frac{1}{\delta }}{\log \frac{1}{1-\epsilon }}\right]$$ (14)

According to the desired precision of the results, δ and ε values are designated. The selected values are: δ=0.01 and ε=0.005. Consequently, N=1000.

The total error (E _p ) is calculated by adding tracking errors (||e _n ||, defined in Eq. (13)) at each sampling instant n = 1, 2, ..., 157:

$$E_p{\mathrm{=T}}_s{\displaystyle \sum _{n=1}^{157}‖e_n‖}$$ (15)

In the last equation, p = 1, 2, ..., N.

Steps summary:

- Define δ and ε, (Tempo and Ishii, 2007).

- Calculate N, Eq. (14).

- Compute the tracking error in each time instant, Eq. (13).

- Obtain the total error at the end of the simulation, Eq. (15).

- Simulate the process N times.

- Compare all the total errors.

- Select the k values that minimize the total error.

It is noteworthy that the invested time for this procedure is less than 1 hour. The k_σ values that minimize the tracking are listed in Table 3.

Table 3
Optimal Controller Parameters.

kX	kP1	kP2	kS1	kS2
0.9992	0.9997	0.9988	0.6771	0.0093

NEURAL STATE ESTIMATION DESIGN

State estimation techniques have a long development history focused mainly to supply the lack of system measurements (Salau et al., 2012). During a fermentation process, variables such as cell and product concentrations are determined by off-line laboratory analysis, making this set of variables of limited use for control purposes (de Assis and Filho, 2000). Although a specific sensor could be available, this kind of hardware is usually expensive. According to this, there are many studies for state estimation in chemical and biochemical processes (Bogaerts and Coutinho, 2014, Bogaerts and Wouwer, 2003, Dewasme et al., 2015, Hulhoven et al., 2006, Lara-Cisneros et al., 2016, Heidarinejad et al., 2012, Kravaris et al., 2007, Hulhoven et al., 2008, Oliveira-Esquerre et al., 2002, Pantano et al., 2017b).

To implement the proposed control technique, it is necessary to have a good variable state estimation to feedback all control variables. Consequently, a neural network state estimation is developed in order to provide cells, bioethanol and glycerol concentrations on-line. Moreover, the estimator has been trained in order to obtain good responses against different perturbations. The estimator design is detailed below.

First it is necessary to define the estimated state error ${\tilde{x}}_n$ as follow:

$${\tilde{x}}_n=x_n-{\widehat{x}}_n$$ (16)

where ${\hat{x}}_n$ represents the estimated state variables of the process, and x_n is the off-line measured state (whose components are: X_n, P_1n, P_2n ). The nonlinear dynamics of the system, described in Eq. (7), can be represented by an exact neural estimator denoted by:

$$\begin{array}{cc}{x}_{n+1}=W^{*T}\xi \left({\upsilon }_n\right)+{\epsilon }_n& x_0=x\left(0\right)\end{array}$$ (17)

where υ _n is the input vector to the neural estimator, being υ _n = [U _n , X _n , P _1,n , P _2,n ], and W* ? R ^m×l is the optimal weight vector, ε ? R ^l×1 is the neural approximation error and ξ _i is the RBF that represent each neuron in the hidden layers, sub-index i indicates the neuron number of radial based function (RBF) functions (ξ ? R ^m×1 ), l is the number of estimated variables (l=3) and m is the maximum number of neurons (m=10)

$$\begin{array}{cc}{\xi }_i\left({\upsilon }_n\right)=\exp \left(-\frac{1}{2{\tau }_i^2}{‖{\upsilon }_n-c_i‖}^2\right)& i=\mathrm{1,2,...}m\end{array}$$ (18)

where c _i and τ _i are parameter vectors of centers and widths of the RBF respectively.

Since the output state vector is non-measured or affected by disturbances, then there is a need to estimate the values. A state estimator function based on Eq. (17) is determined as follows:

$${\widehat{x}}_{n+1}={\widehat{W}}_n^T\xi \left({\upsilon }_n\right)$$ (19)

Now, from the difference between Eq. (17) and Eq. (19), the neural identification error can be described by:

$$\begin{array}{c}{\tilde{x}}_{n+1}=x_{n+1}-{\widehat{x}}_{n+1}=W^{*T}\xi \left({\upsilon }_n\right)-W_n^T\left(k\right)\xi \left({\upsilon }_n\right)+{\epsilon }_n=\mathrm{...}\\ ={\tilde{W}}_n^T\xi \left({\upsilon }_n\right)+{\epsilon }_n\end{array}$$ (20)

where ${\tilde{W}}_n$ represents the neural weight vector estimation error and can be defined as:

$${\tilde{W}}_n=W^{*}-W_n$$ (21)

To train the neural network for identification, an off-line data set of (x, υ) pairs was used. The learning rule to train the neural will be demonstrated in the next theorem:

Theorem: Considering the bioethanol process defined in Eq. (1), it can be approximated by the neural network Eq. (19) using a neuronal adjustment law defined by:

$$\Delta {\tilde{W}}_n=-\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda$$ (22)

where Λ = diag [λ ₁ , λ ₂ , λ ₃ ] is a diagonal definite positive matrix and γ is an arbitrary positive constant.

Proof: the demonstration of this theorem was added in Appendix B; the convergence of this estimator is very important, because an exact representation of the variables to be estimated (${\hat{x}}_n=\left[{\hat{X}}_n,{\hat{P}}_{1,n},{\hat{P}}_{2,n}\right]$) is necessary.

CONTROLLER PERFORMANCE EVALUATION

In this section two different tests are carried out with the intent of demonstrating the reliability of the technique. Those tests are: adding parametric uncertainty and perturbations in the control action. Moreover, the controller is shown while working under normal operation conditions.

Controller performance under normal conditions.

This simulation is made with the aim of having a point of comparison when the different tests are analyzed. To carry it out, a close-loop simulation of the process is executed by using the information from Table 1, Table 2, Table 3 and Fig. 2:. It should be noted that it is being considered that no external disturbances can affect the process.

Figure 3 shows the evolution of cells, ethanol, glycerol and fructose concentration profiles during the simulation; both results with and without estimation are shown. Those profiles are also compared with the references; note how the state variables follow them perfectly. This figure also shows the accumulated error and the tracking error (||e_n|| ). The neural network introduces a limited intrinsic error (ɛ, neuronal approach error, Eq. (17)), this error is accumulated with the progress of the process, the reason why the accumulated error is different for the results with estimation and no-estimation. As can be seen, as the process moves forward, the tracking error tends to decrease and remains limited to low values, i.e., the controller is progressively approaching to the reference in each instant of time.

Figure 3
Comparison between the real cells, ethanol, glycerol and fructose concentrations with the references under normal operation conditions. Note how the tracking error tends to decrease and remains limited to low values.

Test adding parametric uncertainty.

As is known, a particular feature of bioprocesses is the difficulty of measuring their parameters, especially because often they vary over time (Wechselberger et al., 2010). Thus, the next test intends to show how the controller responds when the values of the system parameters are not accurate or fluctuate throughout the process.

To fulfill what is proposed, the Monte Carlo algorithm is used. The number of simulations (N) is determined with Eq. (14), taking into account the same values for δ and ε. During the N simulations, all system parameters are randomly changed by +15% of their original value at the same time (see Table 2).

Fig. 4, shows the concentrations of biomass, fructose and products throughout the process taking into account the changes in the system parameters. It is important to highlight that every reference is almost perfectly followed, which means that this controller gives excellent results even in the presence of uncertainties.

Figure 4
Desired concentration profile variations under parametric uncertainty.

Controller performance with perturbations in the control action.

This test aims to simulate unforeseen events that may cause unwanted variations in the production. To carry it out, a random perturbation in the feed rate of the bioreactor is added. This perturbation affects the control action by +20% of its original value, as is revealed in Fig. 5. This can be explained as a random noise that results in non-zero-mean Gaussian disturbances (George, 2014).

Figure 5
Controller response to the addition of a disturbance in the control action.

In Fig. 5 the perturbed feed rate profile can be seen in comparison with the original one. Moreover, this figure also illustrates how the reference concentrations are followed to perfection despite the applied disturbance. Finally, note how this perturbations cause an increase of the tracking error (with respect to the one analyzed in Figure 3), however, it remains at acceptable levels.

CONTROLLER COMPARISON

As was mentioned before, despite the abundance of sophisticated control techniques for nonlinear systems, PID or PI controllers are still the most widely used in factories (Imtiaz et al., 2014) because they have a simple structure and it is easy to adjust their parameters in comparison with the complexity of online implementation, the insufficient accuracy of online solutions and the computational time required to find a solution offered by other controllers (Jin et al., 2014). Nevertheless, it is well-known that tracking variable set-points is a weakness for those controllers (Aiba, 1979).

This section shows a comparison between the proposed controller and a traditional PID controller in the same bioethanol production system. The best PID parameters were selected with a Monte Carlo algorithm. The technique developed in this study does not present the disadvantages of the aforementioned sophisticated controllers, so it can be implemented in nonlinear systems with minimum requirements. In addition, its performance is significantly better than a PID.

In Fig. 6 a comparison of both controllers’ activity under normal condition is shown. Note the oscillations of the PID response and the difference reflected in the accumulated error. In Fig. 7 and Fig. 8, the functioning of the controllers is manifested in the presence of white noise and step disturbances in the control action, respectively. It is remarkable that, in all three cases, the PID signal saturates, causing a sudden decrease in the substrate concentration which leads to an abrupt drop in cell concentration, and consequently in the product formation. This verifies that a PID controller does not work correctly in non-linear systems. Note that this effect does not occur with the linear algebraic controller, which achieves the tracking of the desired profiles with minimal error.

Figure 6
Comparison of the performance of both controllers under normal conditions.

Figure 7
Performance comparison of both controllers with non-zero-mean Gaussian disturbances in the control action.

Figure 8
Performance comparison of both controllers with -30% step perturbation in the control action.

CONCLUSION

A technique for multivariable system control has been proposed. It has been tested with a fed-batch bioprocess of ethanol production. The way used to calculate the control action allows the minimization of the tracking error obtained, which tends to zero as the process progresses.

One relevant contribution of this paper is the Monte Carlo algorithm application, which enables one to find the controller parameters (k _σ ) successfully and allows evaluating its performance by adding parametric uncertainties. Furthermore, the results of the different tests carried out demonstrate the excellence of the controller.

Moreover, some advantages of this control law are: i) the controller is independent of the operating point because it does not use the linearized model; ii) one does not need specific knowledge to implement it, one is able to use it just with some basic algebra concepts; iii) the form of this technique to represent a system of nonlinear equations allows the calculation of the control action from linear equations, which reduces the mathematical complexity; iv) as it uses discrete equations, it can be programmed in any computer; v) it can be used in several kinds of systems, such as SIMO and MIMO, besides, it cannot be utilized only for bioprocesses, it could be extended to many fields, hence it is a promising method; vi) it is a reliable controller (proved in the simulations section), and its application is inexpensive and fast.

ACKNOWLEDGMENT

The National Council of Scientific and Technological Research (CONICET) and the Chemical Engineering Institute (IIQ) from the National University of San Juan contributed with the support of this investigation.

REFERENCES

• Aiba, S. Review of process control and optimization in fermentation. Biotechnol Bioeng Symp Ser., 269-281. (1979).
• Ajbar, A., Ali, E. Study of advanced control of ethanol production through continuous fermentation. Journal of King Saud University-Engineering Sciences, 29 (1), 1-11, (2017). https://doi.org/10.1016/j.jksues.2015.10.005
  » https://doi.org/10.1016/j.jksues.2015.10.005
• Alvarez, V., Larico, R., Ianos, Y., Aznar, M. Parameter estimation for VLE calculation by global minimization: the genetic algorithm. Brazilian Journal of Chemical Engineering, 25 (2), 409-418, (2008). https://doi.org/10.1590/S0104-66322008000200018
  » https://doi.org/10.1590/S0104-66322008000200018
• Andrášik, A., Mészáros, A., De Azevedo, S. F. On-line tuning of a neural PID controller based on plant hybrid modeling. Computers & Chemical Engineering, 28 (8), 1499-1509, (2004). https://doi.org/10.1016/j.compchemeng.2003.12.002
  » https://doi.org/10.1016/j.compchemeng.2003.12.002
• Asadi, E., Da Silva, M. G., Antunes, C. H., Dias, L., Glicksman, L. Multi-objective optimization for building retrofit: A model using genetic algorithm and artificial neural network and an application. Energy and Buildings, 81 444-456, (2014). https://doi.org/10.1016/j.enbuild.2014.06.009
  » https://doi.org/10.1016/j.enbuild.2014.06.009
• Bogaerts, P., Coutinho, D. Robust nonlinear state estimation of bioreactors based on H∞ hybrid observers. Computers & Chemical Engineering , 60 315-328, (2014). https://doi.org/10.1016/j.compchemeng.2013.09.013
  » https://doi.org/10.1016/j.compchemeng.2013.09.013
• Bogaerts, P., Wouwer, A. V. Software sensors for bioprocesses. ISA transactions, 42 (4), 547-558, (2003). https://doi.org/10.1016/S0019-0578(07)60005-6
  » https://doi.org/10.1016/S0019-0578(07)60005-6
• Castillo, O., Neyoy, H., Soria, J., Melin, P., Valdez, F. A new approach for dynamic fuzzy logic parameter tuning in ant colony optimization and its application in fuzzy control of a mobile robot. Applied soft computing, 28 150-159, (2015). https://doi.org/10.1016/j.asoc.2014.12.002
  » https://doi.org/10.1016/j.asoc.2014.12.002
• Chang, L., Liu, X., Henson, M. A. Nonlinear model predictive control of fed-batch fermentations using dynamic flux balance models. Journal of Process Control, 42 137-149, (2016). https://doi.org/10.1016/j.jprocont.2016.04.012
  » https://doi.org/10.1016/j.jprocont.2016.04.012
• Cheein, F. A., Scaglia, G. Trajectory tracking controller design for unmanned vehicles: A new methodology. Journal of Field Robotics, 31 (6), 861-887, (2014). https://doi.org/10.1002/rob.21492
  » https://doi.org/10.1002/rob.21492
• Chiha, I., Liouane, N., Borne, P. Tuning PID controller using multiobjective ant colony optimization. Applied Computational Intelligence and Soft Computing, 2012 1-11, (2012). https://doi.org/10.1155/2012/536326
  » https://doi.org/10.1155/2012/536326
• Chrobak, M., Kenyon, C., Noga, J., Young, N. E. Incremental medians via online bidding. Algorithmica, 50 (4), 455-478, (2008). https://doi.org/10.1007/s00453-007-9005-x
  » https://doi.org/10.1007/s00453-007-9005-x
• Dai, W., Word, D. P., Hahn, J. Modeling and dynamic optimization of fuel-grade ethanol fermentation using fed-batch process. Control Engineering Practice, 22 231-241, (2014). https://doi.org/10.1016/j.conengprac.2013.01.005
  » https://doi.org/10.1016/j.conengprac.2013.01.005
• De Assis, A. J., Filho, R. M. Soft sensors development for on-line bioreactor state estimation. Computers & Chemical Engineering , 24 (2), 1099-1103, (2000). https://doi.org/10.1016/S0098-1354(00)00489-0
  » https://doi.org/10.1016/S0098-1354(00)00489-0
• De Battista, H., Picó, J., Picó-Marco, E. Nonlinear PI control of fed-batch processes for growth rate regulation. Journal of Process Control , 22 (4), 789-797, (2012). https://doi.org/10.1016/j.jprocont.2012.02.011
  » https://doi.org/10.1016/j.jprocont.2012.02.011
• De Oliveira, L. P., Verstraete, J. J., Kolb, M. A Monte Carlo modeling methodology for the simulation of hydrotreating processes. Chemical engineering journal, 207 94-102, (2012). https://doi.org/10.1016/j.cej.2012.05.039
  » https://doi.org/10.1016/j.cej.2012.05.039
• Del Rio-Chanona, E. A., Zhang, D., Vassiliadis, V. S. Model-based real-time optimisation of a fed-batch cyanobacterial hydrogen production process using economic model predictive control strategy. Chemical Engineering Science, 142 289-298, (2016). https://doi.org/10.1016/j.ces.2015.11.043
  » https://doi.org/10.1016/j.ces.2015.11.043
• Dewasme, L., Fernandes, S., Amribt, Z., Santos, L., Bogaerts, P., Wouwer, A. V. State estimation and predictive control of fed-batch cultures of hybridoma cells. Journal of Process Control , 30 50-57, (2015). https://doi.org/10.1016/j.jprocont.2014.12.006
  » https://doi.org/10.1016/j.jprocont.2014.12.006
• Dimov, I., Maire, S., Sellier, J. M. A new Walk on Equations Monte Carlo method for solving systems of linear algebraic equations. Applied Mathematical Modelling, 39 (15), 4494-4510, (2015). https://doi.org/10.1016/j.apm.2014.12.018
  » https://doi.org/10.1016/j.apm.2014.12.018
• Fu, Y., Chai, T. Nonlinear multivariable adaptive control using multiple models and neural networks. Automatica, 43 (6), 1101-1110, (2007). https://doi.org/10.1016/j.automatica.2006.12.010
  » https://doi.org/10.1016/j.automatica.2006.12.010
• Fujiki, T. L., Sanzovo, M. R., Leite, M. S., Silva, F. V., Fileti, A. M. Optimizing Control Action Online Using a Neural Model and the Solver of an Electronic Worksheet. Computer Aided Chemical Engineering, 27 1461-1466, (2009). https://doi.org/10.1016/S1570-7946(09)70634-0
  » https://doi.org/10.1016/S1570-7946(09)70634-0
• George, J. On adaptive loop transfer recovery using Kalman filter-based disturbance accommodating control. Control Theory & Applications, IET, 8 (4), 267-276, (2014). https://doi.org/10.1049/iet-cta.2013.0671
  » https://doi.org/10.1049/iet-cta.2013.0671
• Griffiths, D. V.,Smith, I. M. Numerical methods for engineers, CRC press. (2006)
• Guay, M., Dochain, D.,Perrier, M. Adaptive extremum seeking control of continuous stirred tank bioreactors with unknown growth kinetics. Automatica, 40 (5), 881-888, (2004). https://doi.org/10.1016/j.automatica.2004.01.002
  » https://doi.org/10.1016/j.automatica.2004.01.002
• Hecklau, C., Pering, S., Seibel, R., Schnellbaecher, A., Wehsling, M., Eichhorn, T., Von Hagen, J.,Zimmer, A. S-Sulfocysteine simplifies fed-batch processes and increases the CHO specific productivity via anti-oxidant activity. Journal of biotechnology, 218 53-63, (2016). https://doi.org/10.1016/j.jbiotec.2015.11.022
  » https://doi.org/10.1016/j.jbiotec.2015.11.022
• Heidarinejad, M., Liu, J.,Christofides, P. D. State-estimation-based economic model predictive control of nonlinear systems. Systems & Control Letters, 61 (9), 926-935, (2012). https://doi.org/10.1016/j.sysconle.2012.06.007
  » https://doi.org/10.1016/j.sysconle.2012.06.007
• Heux, S., Meynial-Salles, I., O’donohue, M.,Dumon, C. White biotechnology: State of the art strategies for the development of biocatalysts for biorefining. Biotechnology advances, 33 (8), 1653-1670, (2015). https://doi.org/10.1016/j.biotechadv.2015.08.004
  » https://doi.org/10.1016/j.biotechadv.2015.08.004
• Hofmann, S., Bajcinca, N., Raisch, J.,Sundmacher, K. Optimal control of univariate and multivariate population balance systems involving external fines removal. Chemical Engineering Science , 168 101-123, (2017). https://doi.org/10.1016/j.ces.2016.12.032
  » https://doi.org/10.1016/j.ces.2016.12.032
• Holland, J. H. Adaptation in natural and artificial systems. An introductory analysis with application to biology, control, and artificial intelligence. Ann Arbor, MI: University of Michigan Press, (1975).
• Hulhoven, X., Vande Wouwer, A.,Bogaerts, P. State observer scheme for joint kinetic parameter and state estimation. Chemical Engineering Science , 63 (19), 4810-4819, (2008). https://doi.org/10.1016/j.ces.2007.11.042
  » https://doi.org/10.1016/j.ces.2007.11.042
• Hulhoven, X., Wouwer, A. V.,Bogaerts, P. Hybrid extended Luenberger-asymptotic observer for bioprocess state estimation. Chemical Engineering Science , 61 (21), 7151-7160, (2006). https://doi.org/10.1016/j.ces.2006.06.018
  » https://doi.org/10.1016/j.ces.2006.06.018
• Hunag, W.-H., Shieh, G. S.,Wang, F.-S. Optimization of fed-batch fermentation using mixture of sugars to produce ethanol. Journal of the Taiwan Institute of Chemical Engineers, 43 (1), 1-8, (2012). https://doi.org/10.1016/j.jtice.2011.06.007
  » https://doi.org/10.1016/j.jtice.2011.06.007
• Imtiaz, U., Assadzadeh, A., Jamuar, S. S.,Sahu, J. N. Bioreactor temperature profile controller using inverse neural network (INN) for production of ethanol. Journal of Process Control , 23 (5), 731-742, (2013). https://doi.org/10.1016/j.jprocont.2013.03.005
  » https://doi.org/10.1016/j.jprocont.2013.03.005
• Imtiaz, U., Jamuar, S. S., Sahu, J.,Ganesan, P. Bioreactor profile control by a nonlinear auto regressive moving average neuro and two degree of freedom PID controllers. Journal of Process Control , 24 (11), 1761-1777, (2014). https://doi.org/10.1016/j.jprocont.2014.09.012
  » https://doi.org/10.1016/j.jprocont.2014.09.012
• Ismail, M., Moghavvemi, M.,Mahlia, T. Genetic algorithm based optimization on modeling and design of hybrid renewable energy systems. Energy Conversion and Management, 85 120-130, (2014). https://doi.org/10.1016/j.enconman.2014.05.064
  » https://doi.org/10.1016/j.enconman.2014.05.064
• Jin, H., Chen, X., Yang, J., Wu, L.,Wang, L. Hybrid intelligent control of substrate feeding for industrial fed-batch chlortetracycline fermentation process. ISA transactions , 53 (6), 1822-1837, (2014). https://doi.org/10.1016/j.isatra.2014.08.015
  » https://doi.org/10.1016/j.isatra.2014.08.015
• Karakuzu, C., Türker, M.,Öztürk, S. Modelling, on-line state estimation and fuzzy control of production scale fed-batch baker’s yeast fermentation. Control Engineering Practice , 14 (8), 959-974, (2006). https://doi.org/10.1016/j.conengprac.2005.05.007
  » https://doi.org/10.1016/j.conengprac.2005.05.007
• Kookos, I. K. Optimization of Batch and Fed-Batch Bioreactors using Simulated Annealing. Biotechnology progress, 20 (4), 1285-1288, (2004). https://doi.org/10.1021/bp034233g
  » https://doi.org/10.1021/bp034233g
• Kravaris, C., Sotiropoulos, V., Georgiou, C., Kazantzis, N., Xiao, M.,Krener, A. J. Nonlinear observer design for state and disturbance estimation. Systems & Control Letters , 56 (11-12), 730-735, (2007). https://doi.org/10.1016/j.sysconle.2007.05.001
  » https://doi.org/10.1016/j.sysconle.2007.05.001
• Lara-Cisneros, G., Aguilar-López, R., Dochain, D.,Femat, R. On-line estimation of VFA concentration in anaerobic digestion via methane outflow rate measurements. Computers & Chemical Engineering , 94 250-256, (2016). https://doi.org/10.1016/j.compchemeng.2016.07.005
  » https://doi.org/10.1016/j.compchemeng.2016.07.005
• Lehouche, H., Guéguen, H.,Mendil, B. Supervisory control based on closed‐loop adaptive control approach of nonlinear systems: application to CSTR process. Asian Journal of Control, 14 (1), 258-270, (2012). https://doi.org/10.1002/asjc.317
  » https://doi.org/10.1002/asjc.317
• Liu, C., Gong, Z., Shen, B.,Feng, E. Modelling and optimal control for a fed-batch fermentation process. Applied Mathematical Modelling , 37 (3), 695-706, (2013). https://doi.org/10.1016/j.apm.2012.02.044
  » https://doi.org/10.1016/j.apm.2012.02.044
• Logist, F., Houska, B., Diehl, M.,Van Impe, J. F. A Toolkit for Multi-Objective Optimal Control in Bioprocess Engineering. IFAC Proceedings Volumes, 43 (6), 281-286, (2010). https://doi.org/10.3182/20100707-3-BE-2012.0063
  » https://doi.org/10.3182/20100707-3-BE-2012.0063
• Mangesh, M. G.,Jana, A. K. A comparison of three sets of DSP algorithms for monitoring the production of ethanol in a fed-batch baker’s yeast fermenter. Measurement, 41 (9), 970-985, (2008). https://doi.org/10.1016/j.measurement.2008.01.008
  » https://doi.org/10.1016/j.measurement.2008.01.008
• Mohammadi, Y., Pakdel, A. S., Saeb, M. R.,Boodhoo, K. Monte Carlo simulation of free radical polymerization of styrene in a spinning disc reactor. Chemical Engineering Journal, 247 231-240, (2014). https://doi.org/10.1016/j.cej.2014.03.013
  » https://doi.org/10.1016/j.cej.2014.03.013
• Mohammadzaheri, M.,Chen, L. Intelligent control of a nonlinear tank reactor. Asian Journal of Control , 13 (3), 439-444, (2011). https://doi.org/10.1002/asjc.347
  » https://doi.org/10.1002/asjc.347
• Motwani, R., Raghavan, P. . Randomized algorithms (cambridge international series on parallel computation). Cambridge University Press, (1995). https://doi.org/10.1017/CBO9780511814075
  » https://doi.org/10.1017/CBO9780511814075
• Ochoa, S. A new approach for finding smooth optimal feeding profiles in fed-batch fermentations. Biochemical Engineering Journal, 105 177-188, (2016). https://doi.org/10.1016/j.bej.2015.09.004
  » https://doi.org/10.1016/j.bej.2015.09.004
• Oliveira-Esquerre, K., Mori, M.,Bruns, R. Simulation of an industrial wastewater treatment plant using artificial neural networks and principal components analysis. Brazilian Journal of Chemical Engineering , 19 (4), 365-370, (2002). https://doi.org/10.1590/S0104-66322002000400002
  » https://doi.org/10.1590/S0104-66322002000400002
• Omar, M., Soliman, M., Ghany, A. A.,Bendary, F. Optimal tuning of PID controllers for hydrothermal load frequency control using ant colony optimization. International journal on electrical engineering and informatics, 5 (3), 348, (2013). https://doi.org/10.15676/ijeei.2013.5.3.8
  » https://doi.org/10.15676/ijeei.2013.5.3.8
• Öztürk, S., Çalık, P.,Özdamar, T. H. Fed-Batch Biomolecule Production by Bacillus subtilis: A State of the Art Review. Trends in Biotechnology, (2016). https://doi.org/10.1016/j.tibtech.2015.12.008
  » https://doi.org/10.1016/j.tibtech.2015.12.008
• Pantano, M. N., Fernández, M. C., Serrano, M. E., Ortíz, O. A.,Scaglia, G. J. E. Trajectory Tracking Controller for a Nonlinear Fed-batch Bioprocess. Revista Ingeniería Electrónica, Automática y Comunicaciones, 38 (1), 78, (2017a).
• Pantano, M. N., Serrano, M. E., Fernández, M. C., Rossomando, F. G., Ortiz, O. A.,Scaglia, G. J. Multivariable Control for Tracking Optimal Profiles in a Nonlinear Fed-Batch Bioprocess Integrated with State Estimation. Industrial & Engineering Chemistry Research, 56 (20), 6043-6056, (2017b). https://doi.org/10.1021/acs.iecr.7b00831
  » https://doi.org/10.1021/acs.iecr.7b00831
• Preuβ, K., Le Lann, M.-V., Cabassud, M.,Anne-Archard, G. Implementation procedure of an advanced supervisory and control strategy in the pharmaceutical industry. Control engineering practice, 11 (12), 1449-1458, (2003). https://doi.org/10.1016/S0967-0661(03)00098-4
  » https://doi.org/10.1016/S0967-0661(03)00098-4
• Roeva, O.,Slavov, T. PID Controller Tuning based on metaheuristic algorithms for bioprocess control. Biotechnology & Biotechnological Equipment, 26 (5), 3267-3277, (2012). https://doi.org/10.5504/BBEQ.2012.0065
  » https://doi.org/10.5504/BBEQ.2012.0065
• Sadatsakkak, S. A., Ahmadi, M. H., Bayat, R., Pourkiaei, S. M.,Feidt, M. Optimization density power and thermal efficiency of an endoreversible Braysson cycle by using non-dominated sorting genetic algorithm. Energy Conversion and Management , 93 31-39, (2015). https://doi.org/10.1016/j.enconman.2014.12.088
  » https://doi.org/10.1016/j.enconman.2014.12.088
• Salau, N. P., Trierweiler, J. O.,Secchi, A. R.. State estimators for better bioprocesses operation. Computer Aided Chemical Engineering . Elsevier. (2012). https://doi.org/10.1016/B978-0-444-59520-1.50112-3
  » https://doi.org/10.1016/B978-0-444-59520-1.50112-3
• Scaglia, G., Rosales, A., Quintero, L., Mut, V., Agarwal, R. A linear-interpolation-based controller design for trajectory tracking of mobile robots. Control Engineering Practice , 18 (3), 318-329, (2010). https://doi.org/10.1016/j.conengprac.2009.11.011
  » https://doi.org/10.1016/j.conengprac.2009.11.011
• Strang, G. Linear algebra and its applications, Thomson, Brooks/Cole, Belmont, CA. ISBN 0-030-10567-6. (2006)
• Tempo, R.,Ishii, H. Monte Carlo and Las Vegas Randomized Algorithms for Systems and Control*: An Introduction. European journal of control, 13 (2), 189-203, (2007). https://doi.org/10.3166/ejc.13.189-203
  » https://doi.org/10.3166/ejc.13.189-203
• Wechselberger, P., Seifert, A.,Herwig, C. PAT method to gather bioprocess parameters in real-time using simple input variables and first principle relationships. Chemical Engineering Science , 65 (21), 5734-5746, (2010). https://doi.org/10.1016/j.ces.2010.05.002
  » https://doi.org/10.1016/j.ces.2010.05.002
• Wood, B. J. Microbiology of fermented foods, Springer Science & Business Media. (2012)
• Xu, L., Ding, F., Gu, Y., Alsaedi, A.,Hayat, T. A multi-innovation state and parameter estimation algorithm for a state space system with d-step state-delay. Signal Processing, (2017). https://doi.org/10.1016/j.sigpro.2017.05.006
  » https://doi.org/10.1016/j.sigpro.2017.05.006
• Ye, J., Xu, H., Feng, E.,Xiu, Z. Optimization of a fed-batch bioreactor for 1, 3-propanediol production using hybrid nonlinear optimal control. Journal of Process Control , 24 (10), 1556-1569, (2014). https://doi.org/10.1016/j.jprocont.2014.08.002
  » https://doi.org/10.1016/j.jprocont.2014.08.002
• Yu, W., Li, B., Jia, H., Zhang, M.,Wang, D. Application of multi-objective genetic algorithm to optimize energy efficiency and thermal comfort in building design. Energy and Buildings , 88 135-143, (2015). https://doi.org/10.1016/j.enbuild.2014.11.063
  » https://doi.org/10.1016/j.enbuild.2014.11.063
• Zhou, K., Doyle, J.,Glover, K. Robust and Optimal Control, Prentice Hall. (1996)

APPENDIX A

This appendix aims to demonstrate the error convergence to zero. As was explained in the article, the values of the k_σ parameters must be between zero and one (0 < k_σ < 1), because that makes the tracking error tends to zero when n tends to infinity. The steps to demonstrate this fact are:

Once the sacrificed variable is determined, the matrix form of the equations system could be generically expressed as:

$$\underset{A}{\underbrace{\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\\ a_5\end{array}\right]}}{U}_n=\underset{b}{\underbrace{\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ b_4\\ b_5\end{array}\right]}}$$ (A.1)

Then, solving Eq. (A.1) with least squares:

$$U_n={\left(A^{T}A\right)}^{-1}{A}^{T}b=\frac{a_1{b}_1+a_2{b}_2+a_3{b}_3+a_4{b}_4+a_5{b}_5}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}$$ (A.2)

From Eq. (A.1):

$$\begin{array}{l}\frac{a_1}{a_2}=\frac{b_1}{b_2}\to b_2=\frac{a_2}{a_1}{b}_1\\ \frac{a_1}{a_3}=\frac{b_1}{b_3}\to b_3=\frac{a_3}{a_1}{b}_1\\ \frac{a_1}{a_4}=\frac{b_1}{b_4}\to b_4=\frac{a_4}{a_1}{b}_1\\ \frac{a_1}{a_5}=\frac{b_1}{b_5}\to b_5=\frac{a_5}{a_1}{b}_1\end{array}$$ (A.3)

Substituting Eq. (A.3) in Eq. (A.2):

$$\begin{array}{l}{U}_n=\frac{a_1{b}_1+(a_2^2{b}_1)/a_1+(a_3^2{b}_1)/a_1+(a_4^2{b}_1)/a_1+(a_5^2{b}_1)/a_1}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}\\ =\frac{(b_1/a_1)(a_1^2+a_2^2+a_3^2+a_4^2+a_5^2)}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}=\frac{b_1}{a_1}\\ =\frac{V_n(((X_{refn+1}-k_X(X_{refn}-X_n)-X_n)/T_S)-({\mu }_1(S_{1ezn},P_{1n})+{\mu }_2(S_{2n},P_{1n}))X_n)}{-X_n}\end{array}$$ (A.4)

Replacing Eq. (A.4) in Eq. (7):

$$X_{n+1}=X_{refn+1}-k_X(X_{refn}-X_n)+T_S[{\mu }_1(S_{1ezn},P_{1n})-{\mu }_1(S_{1n},P_{1n})]X_n$$ (A.5)

Then, as it was mentioned, the tracking error definition is:

$$e_{Xn+1}=X_{refn+1}-X_{n+1}$$ (A.6)

Introducing Eq. (A.5) in Eq. (A.6):

$$e_{Xn+1}=k_X(X_{refn}-X_n)-T_S[{\mu }_1(S_{1ezn},P_{1n})-{\mu }_1(S_{1n},P_{1n})]X_n$$ (A.7)

The Taylor approximation of µ ₁ (S _1n , P _1n ) in the desired value µ ₁ (S _{1ez n} , P _1n ) is:

$$\begin{array}{l}{\mu }_1(S_{1n},P_{1n})={\mu }_1(S_{1ezn},P_{1n})+\left.\frac{d{\mu }_1(S_1,P_{1n})}{d{S}_1}\right|\underset{{}_{S_1=S_{ezn}+\theta (S_{1n}-S_{1ezn})=S_{\theta }}}{(S_{1n}-S_{1ezn})}\\ where\to 0<\theta <1\end{array}$$ (A.8)

Replacing Eq. (A.8) in Eq. (A.7):

$$\begin{array}{l}{e}_{Xn+1}=k_X(X_{refn}-X_n)-T_S\left[{\mu }_1(S_{1ezn},P_{1n})+{\left.\frac{d{\mu }_1(S_1,P_{1n})}{d{S}_1}\right|}_{S_{\theta }}\underset{-e_{S_{1n}}}{\underbrace{(S_{1n}-S_{1ezn})}}-{\mu }_1(S_{1ezn},P_{1n})\right]X_n\\ e_{Xn+1}=k_X(X_{refn}-X_n)+T_S{\left.\frac{d{\mu }_1(S_1,P_{1n})}{d{S}_1}\right|}_{S_{\theta }}{e}_{S_{1n}}{X}_n\end{array}$$ (A.9)

In the same way, it is demonstrated for S ₁ :

From Eq. (A.1):

$$\begin{array}{l}\frac{a_1}{a_2}=\frac{b_1}{b_2}\to b_1=\frac{a_1}{a_2}{b}_2\\ \frac{a_3}{a_2}=\frac{b_3}{b_2}\to b_3=\frac{a_3}{a_2}{b}_2\\ \frac{a_4}{a_2}=\frac{b_4}{b_2}\to b_4=\frac{a_4}{a_2}{b}_2\\ \frac{a_5}{a_2}=\frac{b_5}{b_2}\to b_5=\frac{a_5}{a_2}{b}_2\end{array}$$ (A.10)

Substituting Eq. (A.10) in Eq. (A.2):

$$\begin{array}{l}{U}_n=\frac{(a_1^2{b}_2)/a_2+a_2{b}_2+(a_3^2{b}_2)/a_2+(a_4^2{b}_2)/a_2+(a_5^2{b}_2)/a_2}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}\\ =\frac{(b_2/a_2)(a_1^2+a_2^2+a_3^2+a_4^2+a_5^2)}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}=\frac{b_2}{a_2}\\ =\frac{V_n(((S_{1ezn+1}-k_{S_1}(S_{1ezn}-S_{1n})-S_{1n})/T_S)+\left(q_{S_1/P_1}(S_{1n},P_{1n})/Y_{P_1/S_1}+q_{S_1/P_2}(S_{1n},P_{2n})/Y_{P_2/S_1}\right)X_n)}{\lambda S_f-S_{1n}}\end{array}$$ (A.11)

Replacing Eq. (A.11) in Eq. (A.7):

$$S_{1n+1}=S_{1ezn+1}-k_{S_1}(S_{1ezn}-S_{1n})$$ (A.12)

Then, the tracking error definition is:

$$e_{S_{1n+1}}=S_{1refn+1}-S_{1n+1}$$ (A.13)

Introducing Eq. (A.12) in Eq. (A.13):

$$\begin{array}{l}{e}_{S_{1}n+1}=S_{1ezn+1}-(S_{1ezn+1}+k_{S_1}(S_{1ezn}-S_{1n}))\\ e_{S_{1}n+1}=k_{S_1}(S_{1ezn}-S_{1n})=k_{S_1}{e}_{S_{1}n}\end{array}$$ (A.14)

The same procedure is followed for S ₂ and the result is:

$$e_{S_{2}n+1}=k_{S_2}(S_{2refn}-S_{2n})=k_{S_2}{e}_{S_{2}n}$$ (A15)

In the same way, it is demonstrated for P ₁ :

From Eq. (A.1):

$$\begin{array}{l}\frac{a_1}{a_4}=\frac{b_1}{b_4}\to b_1=\frac{a_1}{a_4}{b}_4\\ \frac{a_2}{a_4}=\frac{b_2}{b_4}\to b_2=\frac{a_2}{a_4}{b}_4\\ \frac{a_3}{a_4}=\frac{b_3}{b_4}\to b_3=\frac{a_3}{a_4}{b}_4\\ \frac{a_5}{a_4}=\frac{b_5}{b_4}\to b_5=\frac{a_5}{a_4}{b}_4\end{array}$$ (A.16)

Substituting Eq. (A.16) in Eq. (A.2):

$$\begin{array}{l}{U}_n=\frac{(a_1^2{b}_4)/a_4+(a_2^2{b}_4)/a_4+(a_3^2{b}_4)/a_4+a_4{b}_4+(a_5^2{b}_4)/a_4}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}\\ =\frac{(b_4/a_4)(a_1^2+a_2^2+a_3^2+a_4^2+a_5^2)}{a_1^2+a_2^2+a_3^2+a_4^2+a_5^2}=\frac{b_4}{a_4}\\ =\frac{V_n(((P_{1ref,n+1}-k_{P_1}(P_{1ref,n}-P_{\mathrm{1,}n})-P_{\mathrm{1,}n})/T_S)-(q_{S_1/P_1,n}(S_{1ez,n},P_{\mathrm{1,}n})+q_{S_2/P_1}(S_{\mathrm{2,}n},P_{\mathrm{1,}n}))X_n)}{-X_n}\end{array}$$ (A.17)

Substituting Eq. (A.17) in Eq. (7):

$$P_{\mathrm{1,}n+1}=P_{1ref,n+1}-k_{P_1}(P_{1ref,n}-P_{\mathrm{1,}n})+T_S[q_{S_1/P_1,n}(S_{1ez,n},P_{\mathrm{1,}n})-q_{S_1/P_1,n}(S_{\mathrm{1,}n},P_{\mathrm{1,}n})]X_n$$ (A.18)

Then, the tracking error definition is:

$$e_{P_1,n+1}=P_{1ref,n+1}-P_{\mathrm{1,}n+1}$$ (A.19)

Substituting Eq. (A.18) in Eq. (A.19):

$$e_{P_1,n+1}=k_{P_1}(P_{1ref,n}-P_{\mathrm{1,}n})-T_S[q_{S_1/P_1,n}(S_{1ez,n},P_{\mathrm{1,}n})-q_{S_1/P_1,n}(S_{\mathrm{1,}n},P_{\mathrm{1,}n})]X_n$$ (A.20)

The Taylor approximation of q _S1/P1 (S _1n , P _1n ) in the desired value q _S1/P1 (S _1ezn , P _1n ) is:

$$\begin{array}{l}{q}_{S_1/P_1,n}(S_{\mathrm{1,}n},P_{\mathrm{1,}n})=q_{S_1/P_1,n}(S_{1ez,n},P_{\mathrm{1,}n})+\left.\frac{d{q}_{S_1/P_1,n}(S_1,P_{\mathrm{1,}n})}{d{S}_1}\right|\underset{{}_{S_1=S_{1ez,n}+\psi (S_{\mathrm{1,}n}-S_{1ez,n})=S\psi }}{\stackrel{-e_{S_{1}n}}{\overbrace{(S_{1n}-S_{1ezn})}}}\\ where\to 0<\psi <1\end{array}$$ (A.21)

Substituting Eq. (A.21) in Eq. (A.20):

$$\begin{array}{l}{e}_{P_1,n+1}=k_{P_1}(P_{1ref,n}-P_{\mathrm{1,}n})-T_S\left[q_{S_1/P_1,n}(S_{1ez,n},P_{\mathrm{1,}n})-q_{S_1/P_1}(S_{1ez,n},P_{\mathrm{1,}n})+{\left.\frac{d{q}_{S_1/P_1,n}(S_1,P_{\mathrm{1,}n})}{d{S}_1}\right|}_{{}_{S\psi }}{e}_{S_1,n}\right]X_n\\ e_{P_1,n+1}=k_{P_1}(P_{1ref,n}-P_{\mathrm{1,}n})-T_S{\left.\frac{d{q}_{S_1/P_1,n}(S_1,P_{\mathrm{1,}n})}{d{S}_1}\right|}_{{}_{S\psi }}{e}_{S_1,n}{X}_n\end{array}$$ (A.22)

The demonstration is similar for P ₂ and the result is:

$$\begin{array}{l}{e}_{P_2,n+1}=k_{P_2}(P_{2ref,n}-P_{\mathrm{2,}n})-T_S{\left.\frac{d{q}_{S_1/P_2,n}(S_1,P_{\mathrm{2,}n})}{d{S}_1}\right|}_{{}_{S\alpha }}{e}_{S_1,n}{X}_n\\ where\to 0<\alpha <1\end{array}$$ (A.23)

Finally, joining Eq. (A.9), (A.14), (A.15), (A.22) and (A.23):

(24) \[\left[ {\begin{array}{*{20}{c}} {{e_{Xn + 1}}}\\ {{e_{{S_1},n + 1}}}\\ {{e_{{S_2},n + 1}}}\\ {{e_{{P_1},n + 1}}}\\ {{e_{{P_2},n + 1}}} \end{array}} \right] = \underbrace {\left[ {\begin{array}{*{20}{c}} {{k_X}}&0&0&0&0\\ 0&{{k_{{S_1}}}}&0&0&0\\ 0&0&{{k_{{S_2}}}}&0&0\\ 0&0&0&{{k_{{P_1}}}}&0\\ 0&0&0&0&{{k_{{P_2}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{e_{Xn}}}\\ {{e_{{S_1},n}}}\\ {{e_{{S_2},n}}}\\ {{e_{{P_1},n}}}\\ {{e_{{P_2},n}}} \end{array}} \right]}_L + \underbrace {{T_S}{X_n}\left[ {\begin{array}{*{20}{c}} {{{\left. {\frac{{d{\mu _{1,n}}({S_1},{P_{1,n}})}}{{d{S_1}}}} \right|}_{_{{S_\theta }}}}}\\ 0\\ 0\\ {{{\left. { - \frac{{d{q_{{S_1}/{P_1},n}}({S_1},{P_{1,n}})}}{{d{S_1}}}} \right|}_{_{S\varphi }}}}\\ { - {{\left. {\frac{{d{q_{{S_1}/{P_2},n}}({S_1},{P_{2,n}})}}{{d{S_1}}}} \right|}_{_{S\alpha }}}} \end{array}} \right]}_{NL}{e_{{S_1},n}}\]

In the last expression, L is a linear system that tends to zero when k _σ values are between zero and one. NL is a bounded nonlinearity that tends to zero, because e _Sn has a tendency to zero; moreover, it is multiplied by a limited term (T _S, X _n ). Eq. (A.24) evidences that the tracking error tends to zero (Scaglia et al., 2010).

APPENDIX B

Remark 2: The optimal approximation error ε should be constant. ||ε|| ≤ $\overline{\epsilon }$

Proof: Defining a candidate function of Laypunov, with L _m a positive definite function given by:

$$L_m\left(k\right)={\tilde{x}}_n^T\Lambda {\tilde{x}}_n+{\gamma }^{-1}tr\left({\tilde{W}}_{n-1}^T{\tilde{W}}_{n-1}\right)$$ (B.1)

where tr(.) is a trace operator. Now, taking the difference in discrete time ΔL _m as follows,

$$\Delta L_m=\left({\tilde{x}}_{n+1}^T\Lambda {\tilde{x}}_{n+1}-{\tilde{x}}_n^T\Lambda {\tilde{x}}_n\right)+{\gamma }^{-1}tr\left({\tilde{W}}_n^T{\tilde{W}}_n-{\tilde{W}}_{n-1}^T{\tilde{W}}_{n-1}\right)$$ (B.2)

Defining ΔW as:

$$\Delta W={\gamma }_1^{-1}tr\left({\tilde{W}}_n^T{\tilde{W}}_n-{\tilde{W}}_{(n-1)}^T{\tilde{W}}_{\left(n-1\right)}\right)$$ (B.3)

Then, organizing Eq. (B.2) and considering $\Delta {\tilde{x}}_n={\tilde{x}}_{n+1}-{\tilde{x}}_n$, it is obtained,

$$\begin{array}{c}\Delta L_m=\left({\tilde{x}}_{n+1}^T\Lambda {\tilde{x}}_{n+1}-{\tilde{x}}_n^T\Lambda {\tilde{x}}_n\right)+\Delta W=\mathrm{...}\\ =\left({\left({\tilde{x}}_n+\Delta {\tilde{x}}_n\right)}^T\Lambda \left({\tilde{x}}_n+\Delta {\tilde{x}}_n\right)-{\tilde{x}}_n^T\Lambda {\tilde{x}}_n\right)+\Delta W=\mathrm{...}\\ =2{\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+\Delta {\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+\Delta W\end{array}$$ (B.4)

Next, substituting Eq. (20) into Eq. (B.4), ΔL _m is written as:

$$\begin{array}{c}\Delta L_m=2{\tilde{x}}_n^T\Lambda \left({\tilde{x}}_{n+1}-{\tilde{x}}_n\right)+\Delta {\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+\Delta W=\mathrm{...}\\ =2{\tilde{x}}_n^T\Lambda \left({\tilde{W}}_n^T\xi \left({\upsilon }_n\right)+{\epsilon }_n-{\tilde{x}}_n\right)+\Delta {\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+\Delta W\end{array}$$ (B.5)

From (B.3), and re-writing ΔW, it yields:

$$\begin{array}{c}\Delta W={\gamma }^{-1}tr\left({\tilde{W}}_n^T{\tilde{W}}_n\left.-{\left[{\tilde{W}}_n-\Delta {\tilde{W}}_n\right]}^T\left[{\tilde{W}}_n-\Delta {\tilde{W}}_n\right]\right)=\mathrm{...}\right.\\ =2{\gamma }^{-1}tr\left({\tilde{W}}_n^T\Delta {\tilde{W}}_n\right)-{\gamma }^{-1}tr\left(\Delta {\tilde{W}}_n^T\Delta {\tilde{W}}_n\right)\end{array}$$ (B.6)

Then, Eq. (B.6), is substituted into Eq. (B.5) where re-organizing terms, it is obtained,

$$\begin{array}{c}\Delta L_m=-2{\tilde{x}}_n^T\Lambda {\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\tilde{W}}_n^T\xi \left({\upsilon }_n\right)+\Delta {\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\epsilon }_n+\Delta W=\mathrm{...}\\ =-2{\tilde{x}}_n^T\Lambda {\tilde{x}}_n+\Delta {\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\epsilon }_n+\mathrm{...}\\ +tr\left(2{\tilde{W}}_n^T\left[\xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda +{\gamma }^{-1}\Delta {\tilde{W}}_n\right]\right)-{\gamma }^{-1}tr\left(\Delta {\tilde{W}}_n^T\Delta {\tilde{W}}_n\right)\end{array}$$ (B.7)

Replacing the adjustment laws (22), $\Delta \tilde{W}$, in , ΔL _m is represented by,

$$\Delta L_m=-2{\tilde{x}}_n^T\Lambda {\tilde{x}}_n+\Delta {\tilde{x}}_n^T\Lambda \Delta {\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\epsilon }_n-{\gamma }^{-1}tr\left({\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)}^T\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)\right)$$ (B.8)

In Eq. (B.8) the increment of the model error of $\Delta {\tilde{x}}_n$ is unknown, and it can be approximated by the following equation:

$$\Delta {\tilde{x}}_{ni}={\left(\frac{\partial {\tilde{x}}_{(n+1)i}}{\partial {\widehat{W}}_{ni}}\right)}^T\Delta {\widehat{W}}_{ni}$$ (B.9)

where the sub-index i denotes the i-th row of $\Delta \tilde{W}$. The partial derivative of the model error depends only on the weights of the neural network and can be rewritten as:

$$\Delta {\tilde{x}}_{ni}=-{\left(\frac{\partial {\widehat{x}}_{(n+1)i}}{\partial {\widehat{W}}_{ni}}\right)}^T\Delta {\widehat{W}}_{ni}$$ (B.10)

Changing the values of the weights according to (B.6), and considering that W* θ^* _f θ^* _f is constant ideal weight vectors required only for analytical purpose. Now re-organizing $\Delta {\tilde{W}}_{ni}$, it can be rewritten as:

$$\begin{array}{l}\Delta {\tilde{W}}_{ni}={\tilde{W}}_{ni}-{\tilde{W}}_{\left(n-1\right)i}=\left(W_{\mathrm{i}}^{*}-{\widehat{W}}_{ni}\right)\\ -\left(W_{\mathrm{i}}^{*}-{\widehat{W}}_{\left(n-1\right)i}\right)=-{\widehat{W}}_{ni}+{\widehat{W}}_{\left(n-1\right)i}=-\Delta {\widehat{W}}_{ni}\end{array}$$ (B.11)

Equation (B.10) can be written as:

$$\Delta {\tilde{x}}_{ni}=-{\left(\frac{\partial {\widehat{x}}_{(n+1)i}}{\partial {\widehat{W}}_{ni}}\right)}^T\left(\gamma {\tilde{x}}_{ni}\xi \left({\upsilon }_n\right){\lambda }_i\right)$$ (B.12)

Considering the value of the partial derivatives of the neural network (22) and replacing in (B.10), yields:

$$\Delta {\tilde{x}}_{ni}=-\gamma {\tilde{x}}_i{\xi }^T\left({\upsilon }_n\right)\xi \left({\upsilon }_n\right){\lambda }_i$$ (B.13)

The increase in the model error is defined as:

$$\Delta {\tilde{x}}_{ni}\le \gamma \left|{\tilde{x}}_{ni}\right|$$ (B.14)

Where

$$\begin{array}{c}\gamma =\gamma {‖\xi \left({\upsilon }_n\right)‖}^2\\ max\left(\xi \left({\upsilon }_n\right)\right)\le 1\end{array}$$ (B.15)

where ||ξ(υ_n))|| is a bounded function and λ_i is a constant. The value γ is a learning factor of the neural network (0<γ<1) and can be arbitrarily defined.

Expressing (B.14) in vector form.

$$‖\Delta {\tilde{x}}_n‖\le \gamma ‖{\tilde{x}}_n‖$$ (B.16)

Substituting the increment value of the model error in Eq. (B.16), the Lyapunov discrete difference is defined as:

$$\begin{array}{c}\Delta L_m={\gamma }^2{\tilde{x}}_n^T\Lambda {\tilde{x}}_n-2{\tilde{x}}_n^T\Lambda {\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\epsilon }_n-{\gamma }^{-1}tr\left({\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)}^T\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)\right)=\\ \Delta L_m=-{\tilde{x}}_n^T\left(\Lambda \left(1-{\gamma }^2\right)\right){\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\epsilon }_n-{\gamma }^{-1}tr\left({\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)}^T\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)\right)\end{array}$$ (B.17)

Now replacing the learning rule in Eq. (B.17) and applying norm, yields:

$$\begin{array}{c}\Delta L_m=-{\tilde{x}}_n^T\Lambda {\tilde{x}}_n+{\gamma }^2{\tilde{x}}_n^T\Lambda {\tilde{x}}_n+2{\tilde{x}}_n^T\Lambda {\epsilon }_n-{\gamma }^{-1}tr\left({\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)}^T\left(\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda \right)\right)\\ \Delta L_m\le -‖\Lambda ‖{‖{\tilde{x}}_n‖}^2+{\gamma }^2‖\Lambda ‖{‖{\tilde{x}}_n‖}^2+2{\lambda }_{\max }\left(\Lambda \right)‖{\tilde{x}}_n‖\overline{\epsilon }-{\gamma }^2\Lambda {‖\xi \left({\upsilon }_n\right)‖}^2{‖{\tilde{x}}_n‖}^2\\ \Delta L_m\le ‖\Lambda ‖{‖{\tilde{x}}_n‖}^2+2‖\Lambda ‖‖{\tilde{x}}_n‖\overline{\epsilon }\end{array}$$ (B.18)

where ||ξ (υ _n )|| ≤ 1, then Eq. (B.18) can be expressed as:

$$\begin{array}{c}\Delta L_m\le -‖\Lambda ‖{‖{\tilde{x}}_n‖}^2+2‖\Lambda ‖‖{\tilde{x}}_n‖\overline{\epsilon }=\mathrm{...}\\ =-{\lambda }_{\min }\left(\Lambda \right)\left(‖{\tilde{x}}_n‖-2\overline{\epsilon }\right)‖{\tilde{x}}_n‖\end{array}$$ (B.19)

Eq. (B.19) allows one to state that the control error is bounded in terms of the NN-RBF approximation error, this estimator is Uniformly Ultimately Bounded (UUB). Only in the case when ε=0, then x _n tends to zero when n tends to infinity.

If the error norm is $||{\tilde{x}}_n||<2\overline{\epsilon }$, it can occur that L>0 and the neural weights could tend to infinity. To prevent the above situation, the next update rule is used.

$$\Delta W=\left\{\begin{array}{c}\gamma \xi \left({\upsilon }_n\right){\tilde{x}}_n^T\Lambda if‖{\tilde{x}}_n‖\ge {\tilde{x}}_0\mathrm{and}‖W‖\le M\\ 0\mathrm{somewhere}else\end{array}\right.$$ (B.20)

where ${\tilde{x}}_o$ and M are design parameters.

Publication Dates

• Publication in this collection
  15 July 2019
• Date of issue
  Jan-Mar 2019

History

• Received
  17 July 2017
• Reviewed
  08 Aug 2018
• Accepted
  13 Aug 2018