Optimal Control Systems Using Evolutionary Algorithm-Control Input Range Estimation

: The closed-loop optimal control systems using the receding horizon control (RHC) structure make predictions based on a process model (PM) to calculate the current control output. In many applications, the optimal prediction over the current prediction horizon is calculated using a metaheuristic algorithm, such as an evolutionary algorithm (EA). The EAs, as other population-based metaheuristics, have large computational complexity. When integrated into the controller, the EA is carried out at each sampling moment and subjected to a time constraint: the execution time should be smaller than the sampling period. This paper proposes a software module integrated into the controller, called at each sampling moment. The module estimates using the PM integration the future process states, over a short time horizon, for different control input values covering the given technological interval. Only a narrower interval is selected for a ‘good’ evolution of the process, based on the so-called ‘state quality criterion’. The controller will consider only a shrunk control output range for the current sampling period. EA will search for its best prediction inside a smaller domain that does not cause the convergence to be affected. Simulations prove that the computational complexity of the controller will decrease.


Introduction
Dynamic system optimization subjected to a performance index is a common task in process engineering. Theoretical control laws can be applied when the process has sufficient mathematical properties and generates controllers with moderate computational complexity. Sometimes the process has profound nonlinearities or imprecise, uncertain, incomplete knowledge. In these cases, a metaheuristic algorithm, integrated within an adequate control structure, can be a realistic solution.
Metaheuristic algorithms (see [1][2][3]) have been employed intensively in control engineering (see [4][5][6][7][8]) due to their robustness and capacity to deal with high complexity problems. However, in these cases, the computational effort is very large and can be crucial for the controller to meet time constraints.
The RHC is a closed-loop control structure used in many works. The key factor is a process model, allowing the prediction of the process evolution. The controller organizes the moving of the prediction horizon in a specific way (see, for example, [9][10][11]). The model predictive control is a well-known particular case of the RHC, which makes at each sampling period a specific action: the prediction error's minimization.
Many papers have considered metaheuristics (genetic algorithm, simulated annealing, particle swarm optimization etc.) in conjunction with the RHC to implement closed-loop structures, which have been used successfully in real-time control. Reference [6] has a section that surveys this kind of work. Reference [12] introduced evolutionary predictive control, a technique for designing predictive controllers. This technique generates and evaluates a family of optimum predictive controllers using evolutionary algorithms. They have design parameters that are adapted at every sampling period. Finally, the best performer is selected using this adaptive process.
A systematic design procedure using the RHC (theoretic approach in [13]) to generate a closed-loop structure has been proposed in [14]. The optimal control structure integrates a metaheuristic algorithm slightly modified to play the role of controller. Obviously, an EA may be used with that end in view.
This work addresses the Park-Ramirez problem (PRP) (see for example [15,16]) that deals with the optimal production of secreted protein in a fed-batch reactor and is a kind of benchmark optimal control problem (OCP).
In a previous work [17], the authors have proposed a closed-loop solution for this problem using the RHC (see [18,19]). This structure implements the closed-loop control that makes optimal predictions based on a PM at each sampling period. The controller presented in [17] used an EA to make predictions and calculate the quasi-optimal control output. However, that paper's objective was to develop an analysis method for the optimal character degradation when the process is different from the PM.
The EAs, as other population-based metaheuristics, have high computational complexity. When integrated into the controller, the EA is carried out at each sampling moment and subjected to a time constraint: the execution time should be smaller than the sampling period. That is why computational complexity decrease of the EA in control applications is an important topic. For this purpose, the paper [20] proposes a method to encode the predicted sequences, which reduces the computational complexity.
The present work has the same context: the same problem, control structure and a similar EA within the controller. Still, it is devoted to another desideratum: to propose another method that improves the computational complexity.
Our objective has not been to implement a more effective control algorithm than others reported in previous papers, which would give better solutions for this particular problem.
The contribution of this paper is a software module (called 'ESTIMATOR') integrated into the controller, which could reduce the control output range of the controller at each sampling period. The process control input is, at the same time, the control output of the controller (it is the direct connection between the controller and the process). The ESTIMATOR will set the shrunk range by numerical integration of the PM for a short while (a sampling period at most). In this way, the EA will search for its best prediction inside a smaller domain that does not cause the convergence to be affected. Therefore, the computational complexity of the controller will decrease. From the design perspective, the time constraints have more chances to be met. Section 2 describes the scientific framework of this paper and introduces the notations to the lecturer. The objectives of this section are:

•
To recall the OCP's elements that are exemplified with the help of PRP.

•
To recall the RHC structure and optimal predictions. • To introduce the EA as the core of the predictor.
The paper's contribution is presented in Section 3: the proposed method to reduce the control input range by an ESTIMATOR, the new structure of the controller and its algorithm and the repercussions of this method to the EA. Sections 4 and 5 describe the results of the two closed-loop simulations: without and with shrunk control output ranges. Section 6 compares the results and shows that the computational complexity decreases by reducing the control input range.

Optimal Control Problem
In this subsection, we recall the PRP, a kind of benchmark problem describing an OCP. It is about a nonlinear process concerning the secreted protein production in a fed-batch reactor. Many papers reported quasi-optimal solutions obtained with different numerical integration techniques, among those we recall here [11]. We present hereafter the three parts defining the PRP as an OCP.

Process Model
The protein production process is modeled by algebraic and ordinary differential equations like in many papers (see [4,5]). PRP has the following dynamic model: The five state variables have the following meanings: The process has a single control input variable: • u(t): the nutrient feed rate.

Constraints
The ordinary differential equations are accompanied by a set of constraints to precise the process evolution. In the case of the PRP, there is a minimum set of constraints that have to be met: where t f denotes the final moment of the control horizon initialstate : To solve the PRP using an AE, we can discretize the control input and try to find the optimal 'command profile' (see [4]) that will be a chromosome. Hence, we consider a supplementary constraint: the control input is constant during each sampling period. The control horizon and the variable u are discretized as follows: We consider a new constraint (5) that replaces (4) Finally, the control profile with n = 15 constant values is u = (u 1 , u 2 , . . . , u n )

Performance Index
The solution of PRP is the control profile that maximizes the objective function (6) Automation 2022, 3

98
The performance index determined by this solution is Taking into account that constraint (4) can be replaced by (5), the PRP's optimal solution is the control profile u = (u 1 , u 2 , . . . , u n ) that meets all the constraints and maximizes the objective function (6).

Controller Using Predictions and An Evolutionary Algorithm
Generally speaking, PRP is one of the problems that can be solved through the RHC structure.
NB: To keep the presentation self-contained, we recall in Appendix A (Supplementary Materials) only the elements that define the well-known RHC structure. The controller corresponding to this control structure has the pseudocode description given below through Algorithm 1. The latter is rather a generic way to describe the receding horizon controller in real-time applications.
Algorithm 1 calls a generic function (called "Predictor_EA") that calculates the optimal prediction for the current sampling period, in line #2. The prediction can be calculated using any sound method, including metaheuristics, adapted to the PM and its constraints. In this work, we used a metaheuristic algorithm, namely the EA. The "Predictor_EA" function has two parts (detailed in Section 3.2):

•
The first part sets the bounds of the control inputs according to the control input range estimation, which is the contribution of this work.

•
The second part is an EA, whose details are given in Section 2.3, calculating the optimal prediction.
We denote by [0, H] the control horizon from our problem (t f = H · T). A prediction is an optimal command profile u(k) over the prediction horizon [k, H], k = 0, 1, · · · , H − 1]. A candidate prediction has the structure u(k) =< u(k|k), · · · , u(k + i|k), · · · , u(H − 1|k) >, (8) where, u(k + i|k), i = 0, . . . , H − k − 1 is the predicted value for the control input u(k + i) based on knowledge up to the moment k. Using Equations (1) and (8), the PM returns the predicted state sequence x(k) as where x(k + i|k ), i = 0, . . . , H − k is the predicted value of the state x(k + i) based on knowledge up to moment k. The EA computes the objective function value over the prediction horizon for each chosen sequence (individual of the population). After convergence, the EA returns the optimal control sequence (denoted u * ) Note that this optimal prediction is made at the beginning of the sampling period [k, k + 1]. The controller will send the value of the first element to the process as a real control input value, namely u * (k) = u * (k|k).
Algorithm 1 outlines the controller's algorithm carried out for every sampling period. The prediction is achieved using the function "Predictor_EA" in Line #2. Its input parameters are the current sampling period and the process state vector, and it returns the value u * (k) according to Equation (10). Line #3 sets the optimal control output u * (k) at the value u * (k|k) .

Algorithm 1
Receding Horizon Controller's algorithm using predictions based on EA. 1 Get the current value of the state vector, x(k); /* Initialize k and Send u * (k) towards the dynamic system 5 Update the prediction horizon and wait for the next sampling period Line #5 makes preparations for the next prediction horizon. Owing to performance index (7), which is a terminal penalty, the right limit of the prediction horizon must be H even if the integral term is missing. The prediction horizon 'recedes' at each sampling period but keeps the final extremity. Its length decreases by one unit at each sampling period.
The prediction is essentially made by the EA included in the "Predictor_EA" function. The prediction horizon is [k, H] for the current sampling period, and its length is H − k. Therefore, the EA' computational complexity is variable and depends on the k value. We have chosen the discretization step for encoding the sequence u(k) to be equal to the sampling period T. Therefore, the length of the sequence u(k) is variable along the control horizon and decreases at each sampling period.
The execution time must be smaller than the sampling period; this is the more restrictive time constraint for the algorithm presented in Algorithm 1. The length of the prediction horizon could be very large, especially in the first sampling period; this is why we consider the decreasing of the computational complexity of EA to sometimes be crucial. The control system designer must verify through simulation this constraint to be met. The receding horizon controller, which uses an EA, is suitable for relatively large time constants processes.

Description of the Evolutionary Algorithm
Our objective has not been to implement a more effective control algorithm than others reported in previous papers, which would give better solutions for this particular problem.
The EA ( [3,15,16]) uses direct encoding, as mentioned previously. A chromosome having H − k genes represents the predicted control sequence over the interval [k, H], which comprises the control input's values that are supposed constant during each sampling period.
NB: The implemented EA has part of characteristics similar to those the authors proposed in the paper [14], except the crossover operator and the stop criteria. These characteristics are listed below: • Each population has µ individuals (predictions); • The number of children generated in each generation: λ (λ < µ); • The selection strategy: stochastic universal sampling; the individuals are ranked using the selection pressure (s) (see [3]); • A single point crossover operator; • The mutation step size is adapted; its global variance is modified according to the "1/5 success rule"; see [1] (pp. 245-274).

•
The replacement strategy causes the λ worst generation individuals to be replaced by the children.

•
The first stop criterion is that the population would evolve during NGen generations.

•
The second stop criterion is that the best individual has its objective function's value greater than or equal to J 0 . The value of J 0 is preestablished (see Section 6).
Other details concerning the implementation of the EA used in this work are given in Appendix C.

Potential Reduction of the Control Input Range
To keep the presentation simple, we consider the control input has only one dimension (like in our example), but the following considerations are also valid for the multidimensional case. Therefore, the set of control input values U is a real interval (like in (4)).
For each k, we denote by U k the control input domain, namely the set of all values that can be taken by u(k). The set U k generally corresponds to the technological limits. Hence, the control input value meets the constraint In real applications, the control input value is subjected to the constraint The minimum and maximum bounds u min and u max are mainly technological limits. The controller output can take values between these limits. The process control input u(k) is, at the same time, the control output given by the controller because it is the direct connection between the controller and the process.
In our example, the minimum and maximum limits of the controller output imposed by constraint (4) are The EA, which searches the ocs, generates an initial population consisting of control input values sequences. For fault of other information, these values fulfil the constraint (11) for all k. Hence, we have Despite the stochastic character of the EA, this choice will cause an important computational effort to find out the best control sequence. Therefore, the algorithm's computational complexity making the predictions will be very big.
The main contribution of our work is to propose an approach to obtain additional information that allows the control input range at moment k (denoted by R k ) to be reduced. We expect to have Remark 1: The following points briefly show how to reduce the control input range.

1.
A software module called ESTIMATOR will evaluate the process statesx(k + ∆t) in the following conditions: the initial state x(k) is the current state of the process that is known the evaluation is made by integration over the interval I = [k, k + ∆t], where ∆t is lesser than or equal to the sampling period 2.
The integrations consider a representative list of control inputs u(k) covering the set U k . For example, in this work, we have considered the list 3.
The set R k is defined as the interval covering the values belonging to L that transfers the process to a final state fulfilling a certain quality criterion.

4.
The quality criterion expresses that the final statex(k + ∆t) has a property, which is mandatory for the desired system behavior.

5.
The ESTIMATOR returns the two values defining the interval R k = [u m , u M ].
These aspects are illustrated in the sequel for Process (1). The approach described before can be applied if we have identified a quality criterion for the final state x(k + ∆t). Analyzing the PRP's good solutions, we have noted that the culture cell density x 3 (t) is monotonically increasing as a function of t. Therefore, we have considered the following quality criterion After process integration for all the values belonging to L, the ESTIMATOR must select only the control input values u(k) that produce final states x(k + ∆t) meeting the quality criterion (16). Finally, the ESTIMATOR returns the narrower interval that covers all these values u(k). (16) is insufficient to assure the monotony, but the value u(k) at hand can be eliminated if not met.

Remark 2: Constraint
The value x 3 (k) is the initial state for the current sampling period, and it is known at the moment k. Because it is the initial state for all integrations, we denote this value x 30 : The value x 3 (k + ∆t) could be considered as a function of control input u(k). Figure 1 shows the dependence of the estimated state quality Automation 2022, 3, FOR PEER REVIEW 8 Remark 3: A small value for alfa  will cause the efficiency of the new control range to be very small. The range k R is too large, and the method's advantage disappears. On the other hand, a big value for alfa  could involve the EA not to converge to the optimal solution. The range k R is too narrow, and the EA does not find the optimal control input values.
Simulations of the controller functioning for the PRP have proven an appropriate value Figure 1 is obtained using the ESTIMATOR and other programs, which the lecturer can use. Process (1) has an optimal evolution according to performance index (7), and the current state is ( ) x k with 5 k  . The ESTIMATOR evaluates the future state (6) x ( t  equals the sampling period) for the control input values belonging to L. Because the state quality (17)  The blue segment corresponds to the shrunk control input range, that is The new range length is 30% of the initial one and entails a smaller computational complexity that will be analyzed in Section 6.
If a quality criterion can be stated, then the current state of the Process entails, for certain physical reasons, a narrower interval can be used to diminish the computational complexity. The control input values will be looked for within a smaller interval. Hence, the constraint (11) could be replaced by In our analysis, the minimum bound is always 0.05, which is not useful for having a shrunk control input range. That is why, for u m -the minimum value of R k -we have adopted the rule Remark 3: A small value for alfa α will cause the efficiency of the new control range to be very small. The range R k is too large, and the method's advantage disappears. On the other hand, a big value for alfa α could involve the EA not to converge to the optimal solution. The range R k is too narrow, and the EA does not find the optimal control input values.
Simulations of the controller functioning for the PRP have proven an appropriate value α = 0.2. Hence the ESTIMATOR returns Figure 1 is obtained using the ESTIMATOR and other programs, which the lecturer can use. Process (1) has an optimal evolution according to performance index (7), and the current state is x(k) with k = 5. The ESTIMATOR evaluates the future state x(6) (∆t equals the sampling period) for the control input values belonging to L. Because the state quality (17) has to be positive, the maximum control input value is u M = 0.75.
The blue segment corresponds to the shrunk control input range, that is The new range length is 30% of the initial one and entails a smaller computational complexity that will be analyzed in Section 6.
If a quality criterion can be stated, then the current state of the Process entails, for certain physical reasons, a narrower interval R k = [u m , u M ] such that u(k) ∈ R k . This fact can be used to diminish the computational complexity. The control input values will be looked for within a smaller interval. Hence, the constraint (11) could be replaced by Remark 4: The approach described before can be applied if we have identified a quality criterion for the final statex(k + ∆t). We can identify specific quality criteria considering certain physical and chemical reasons or good solutions' properties. Moreover, the quality criterion has to reduce the control input range significantly. We emphasize that the new control input range depends on the current process state.

Controller Structure with ESTIMATOR
A preliminary analysis-based on Remark 4-can conclude whether the method to reduce the control input range can be applied to improve the EA's computational complexity. After the quality criterion is checked through simulation, the controller could include the ESTIMATOR module. Figure 2 presents the receding horizon controller with the structure described, for example, in [9,13,14]. The controller gathers mainly the prediction, optimization, and PM modules. There are certainly auxiliary modules that assure the process state x(k) to be acquired, the optimal control output to be sent to the process, and other actions to be achieved. structure described, for example, in [9,13,14]. The controller gathers mainly the prediction, optimization, and PM modules. There are certainly auxiliary modules that assure the process state ( ) x k to be acquired, the optimal control output to be sent to the process, and other actions to be achieved.
In our work, the optimization is accomplished by the EA that is included in the prediction module.  . The EA considers this new shrunk range and determines the optimal prediction for the current moment. Algorithm 2 is a version of Algorithm 1 that integrates the ESTIMATOR and can also be used in applications, which are not in real-time (for example, simulations). It calculates the optimal control output *( ) u k for the current sampling period. We recall that the predictions are made by the EA integrated into the "Predictor_EA" function (the EA is not a distinct function).

x(k)
Rk= [um, uM] x(k+Δt) Prediction and optimization Evolutionary Algorithm In our work, the optimization is accomplished by the EA that is included in the prediction module.

Receding Horizon Controller
The ESTIMATOR interacts with the PM and prediction modules and calculates the control input range R k = [u m , u M ]. The EA considers this new shrunk range and determines the optimal prediction for the current moment.
Algorithm 2 is a version of Algorithm 1 that integrates the ESTIMATOR and can also be used in applications, which are not in real-time (for example, simulations). It calculates the optimal control output u * (k) for the current sampling period. We recall that the predictions are made by the EA integrated into the "Predictor_EA" function (the EA is not a distinct function). If If the calling program wants to use the shrunk control input range, line #6 calls the ESTIMATOR function, only one time for each prediction horizon [k, H] to calculate R k = [u m , u M ]. This estimation is based on the state vector x(k) that is known because it is an input parameter.
Line #7 of Algorithm 2 calls the prediction function slightly modified because of the two new parameters (u m , u M ). Therefore, for the current predicted sequence, it holds

Predictions Using the EA
This subsection gives some details concerning the "Predictor_EA" function, whose main part is an EA. The latter one returns the optimal control sequence that maximizes the performance index over the prediction horizon ([k, H]) given the initial state x(k) u * (k) = arg max u(k) J(k, x(k), u(k)). The predictor sets the bounds (20) for each component of the predicted sequence and then uses the EA (we recall that the EA is a part of the "Predictor_EA" function). The EA generates the initial chromosome population, i.e., the control sequences ( ( ) u k ). Each one meets the constraint (20). After the random selection of the predicted sequence components inside the admissible domains, it holds The EA can now calculate all predicted state sequences (9) and evaluate their objective function values to determine the optimum.
Remark 5 states why the method that reduces the control input range can be applied only to the first component of the control input sequence. Hence, the improvement of the computational complexity is limited but quite important, as the next sections will prove.
The general structure of the Predictor_EA function is given in Algorithm 3. The main task of this function is to pave the way for the EA that will treat the bounds of the control input ranges differently. The range corresponding to the first sampling period ([k,k + 1]) will be defined by the two parameters m x and M x received from the ESTIMATOR. The generation of the initial population is achieved by lines #6-11.
The remaining part of this function is a general description of a usual EA, which implementation is not a key factor of our work. Nevertheless, the lecturer will find in Appendix C some details concerning a simple implementation, including the function "eval_PR.m" that calculates the objective function J (Equation (6)).  The predictor sets the bounds (20) for each component of the predicted sequence and then uses the EA (we recall that the EA is a part of the "Predictor_EA" function). The EA generates the initial chromosome population, i.e., the control sequences (u(k)). Each one meets the constraint (20). After the random selection of the predicted sequence components inside the admissible domains, it holds The EA can now calculate all predicted state sequences (9) and evaluate their objective function values to determine the optimum.
Remark 5 states why the method that reduces the control input range can be applied only to the first component of the control input sequence. Hence, the improvement of the computational complexity is limited but quite important, as the next sections will prove.
The general structure of the Predictor_EA function is given in Algorithm 3. The main task of this function is to pave the way for the EA that will treat the bounds of the control input ranges differently. The range corresponding to the first sampling period ([k, k + 1]) will be defined by the two parameters x m and x M received from the ESTIMATOR. The generation of the initial population is achieved by lines #6-11.
The remaining part of this function is a general description of a usual EA, which implementation is not a key factor of our work. Nevertheless, the lecturer will find in Appendix C some details concerning a simple implementation, including the function "eval_PR.m" that calculates the objective function J (Equation (6)). The implementation of the Predictor_EA function and the simulations were carried out using the MATLAB language and system. We have employed special functions devoted to integrating the differential equations to determine the process evolution. The value of H is a global constant, defined in Section 2.2 (see also Figure 3).
In our implementation, the EA also counts the objective function's calls number (Ncalls) to measure the calculation effort in generating the current optimal prediction. This variable is updated inside the function "eval_PR.m".

Objectives and Simulation's Hypothesis
In the context of using RHC to achieve closed-loop optimal control of a given process, the optimal predictions can be implemented through a metaheuristic algorithm. The solution is realistic but has a major inconvenience: the computational complexity of the controller.
In the previous sections, we proposed the ESTIMATOR subsystem integrated into the controller as a tool that can decrease the latter's computational complexity. Therefore, the major objective of our simulation study is to validate this hypothesis.
We have chosen for a case study the PRP, which is a well-known control problem reported in many papers to be solved in open loop and closed loop as well. As we are interested in control implementation aspects-not only computational intelligence and numerical aspects-we have treated the closed-loop problem, which is more complex. The fact that the PRP is already solved and analyzed presents a great advantage. Now, we can compare its known solution with the new one obtained using the RHC control structure, EA (for the predictions), and the proposed ESTIMATOR.
The simulation study presented in the sequel had a few objectives listed below:

1.
To implement and simulate the predictor function specific to the PRP, which uses the EA and considers the shrunk control input range.

2.
To implement and simulate a closed-loop based on RHC to solve the PRP problem.

3.
To compare the quasi-optimal closed-loop solutions produced by the controller without and with the ESTIMATOR module.

4.
To evaluate the computational complexity of the closed loop over the control horizon without and with the ESTIMATOR module.
The simulation study will be made by an application that emulates the closedloop functioning.

•
The real process is also simulated using the same PM. The red lines in Figure 2 show the simulated connections that create the closed loop.

•
The sequence of sampling moments is simulated.
The objective function's calls number (Ncalls) for each sampling period is cumulated with the help of variable Ncalls_C (line #7 and line #12). The Ncalls_C's value will be necessary to measure the computational complexity for the entire control horizon (all the sampling periods).
Function "RealProcessStep" calculates by numerical integration the process state at the next sampling moment (k + 1), using its arguments: the optimal control input and the initial state at the moment k. The next state becomes the initial one for the moment (k + 1). Despite its prefix, "RealProcess", this function will integrate the PM because actually, we are not interested in simulating the closed loop with a process different from its model. Still, we want to analyze the behavior of the closed loop both with and without ESTIMATOR.

Simulation of the Closed Loop without Control Input Range Estimation
This section presents the simulation results of the RHC structure solving the PRP. Here, we present the input data for the problem we have solved.
The EA has the parameters presented in Appendix C and tuned for this application.
Because the prediction is based on a stochastic algorithm, the "RHC_EA_new.m" program-that implements the RHC_EA algorithm-was carried out 30 times to analyze the results statistically. Practical details about this operation are given in Appendix D.1.
The objective function's calls number over the control horizon is a realistic measure of the computational complexity. The calls number is relevant, especially when the objective function involves numerical integrations, which have significant complexity and are timeconsuming. This measure is also adequate for comparison among versions of applications. The RHC_EA algorithm totalizes the calls' number for each sampling period, and the script "STATISTIC_DRAW30_0" calculates the average value for a sampling period, denoted Ncalls, in Table 1. The main results of this simulation series without range estimation are given in Table 1. For each execution, four values are displayed: the optimum criterion's value (J), the final value of x 1 , the final value of x 5 , and the average number of calls (Ncalls). The latter equals the number of calls divided by the sampling periods' number (H = 15). The sampling periods have very different calls' numbers: the greater the value of k, the smaller the number of calls. However, it is easier to compare the values of Ncalls. Thus, the total number of calls over the control horizon is 15*Ncalls.
A statistic of this simulation series; the minimum, average, and maximum values; and standard deviation of the objective function (J), is given in Table 2. We consider a simulation as typical execution if it produces the closest value to the average optimum criterion. Using the optimum criterion's value, we have determined the typical run among the 30 simulations.
In our simulation series, a specific execution yields Jtypical = 31.900. The other results associated with this run (line 30 in Table 1)  The controller calculates the best control output value as the first element of its predicted sequence. These 15 values of the quasi-optimal control output are recorded by the RHC_EA algorithm within the uRHC vector, tracing the controller's activity over the control horizon. Hence, this vector is not the result of a prediction made at the moment k = 0; in other words, it is not the sequence u * (0). This vector allows a final simulation of the closed loop.
The vector uRHC for the typical execution is depicted in Figure 4.
We consider a simulation as typical execution if it produces the closest value to t average optimum criterion. Using the optimum criterion's value, we have determined t typical run among the 30 simulations.
In our simulation series, a specific execution yields Jtypical = 31.900. The other resu associated with this run (line 30 in Table 1)  The controller calculates the best control output value as the first element of its pr dicted sequence. These 15 values of the quasi-optimal control output are recorded by t RHC_EA algorithm within the uRHC vector, tracing the controller's activity over the co trol horizon. Hence, this vector is not the result of a prediction made at the moment k = in other words, it is not the sequence *(0) u . This vector allows a final simulation of t closed loop.
The vector uRHC for the typical execution is depicted in Figure 4. The control outputs of uRHC sent to the Process involve the typical state evoluti depicted in Figure 5. The control outputs of uRHC sent to the Process involve the typical state evolution depicted in Figure 5.
The computational complexity of the first sampling period is the largest because the prediction horizon has H sampling periods. Despite the sequence u * (0) having an H length, the controller calculates the control variable in a few tens of seconds (depending on processor speed). This duration is quite satisfactory for a sampling period of 1 h.  The computational complexity of the first sampling period is the largest because the prediction horizon has H sampling periods. Despite the sequence *(0) u having an H length, the controller calculates the control variable in a few tens of seconds (depending on processor speed). This duration is quite satisfactory for a sampling period of 1 h.

Simulation of the Closed Loop with Control Input Range Estimation
In this implementation, ESTIMATOR is called within RHC_EA to reduce the control input range. The "RHC_EA_new.m" program was carried out 30 times for the same instance of the PRP as in the previous section. Practical details about how the simulation results are generated can be found in Appendix D2.
The results are presented in Table 3, having the same columns as Table 1. The new statistic over the 30 runs of the simulation series is given in Table 4.

Simulation of the Closed Loop with Control Input Range Estimation
In this implementation, ESTIMATOR is called within RHC_EA to reduce the control input range. The "RHC_EA_new.m" program was carried out 30 times for the same instance of the PRP as in the previous section. Practical details about how the simulation results are generated can be found in Appendix D.2.
The results are presented in Table 3, having the same columns as Table 1. The new statistic over the 30 runs of the simulation series is given in Table 4. The vector uRHC for the typical execution is depicted in Figure 6. These control output values generate the state evolution illustrated in Figure 7.  The vector uRHC for the typical execution is depicted in Figure 6. These control output values generate the state evolution illustrated in Figure 7.

Discussion
We first recall that this work has the same context as that presented in [12]. The basic version of algorithm RHC_EA was created similarly to the previous algorithm (with small differences concerning the EA). We proposed the control input range estimation and compared the closed-loop function in the two algorithm versions, without and with control input range estimation.
The simulations have proved that the RHC_EA program converges for all executions.

Discussion
We first recall that this work has the same context as that presented in [12]. The basic version of algorithm RHC_EA was created similarly to the previous algorithm (with small differences concerning the EA). We proposed the control input range estimation and compared the closed-loop function in the two algorithm versions, without and with control input range estimation.
The simulations have proved that the RHC_EA program converges for all executions. Aiming to obtain good solutions for PRP and a satisfactory computation time, the predictor stops the iterative process when for the best solution encountered it holds The stop criterion is either to fulfil constraint (21) or to evolve during Ngen generations. The reported optimum value of J is 32.6, but for algorithms that solve PRP standing alone, not in a closed-loop implementation. As EA now works in a closed loop, its stop-criterion must assure good solutions and adequate execution time for the controller. The main findings of our simulation series are that Figures 5 and 7 ascertain the process has the same quasi-optimal behavior. In the two cases (without and with control input range estimation), the average values of the performance index J are 31.9 and 31.88, respectively. Hence, the ESTIMATOR does not worsen the closed-loop behavior:

•
The optimal behavior is not degraded; the two J values are practically identical.

•
The process evolutions are very alike.

•
The convergence is not impeded.
On the contrary, the RHC_EA algorithm with control input range estimation has a beneficial influence on the computational complexity. Table 5 shows this improvement. Hence, the computational complexity was diminished by 18.46% using the control input range estimation. This result was expected, but the reduction is large enough considering that range reduction applies only to the first sampling period of the predicted sequences.
The controller meets the time constraint all the more so since the version without ESTIMATOR already fulfils it. Hence, the control structure could be a solution even for real-time control.
In our simulations, we decided that the process is identical to the PM, which is generally not realistic. However, we were interested in seeing the net contribution of the range estimation to the decrease of the computational complexity. If we consider simultaneously perturbative factors (e.g., process = PM), the computational complexity will increase to ensure the controller's optimal character. On the other hand, the papers [9] and [12] treated to a large extent this situation, when process = PM, and how simulations can estimate the optimal character degradation.
The ESTIMATOR worked in these simulations having the current state given by the PM. Therefore, it could be useful, at least in simulation studies. In real-time applications, the current state would be given by the real process. Data Availability Statement: The study did not report any data.

Conflicts of Interest:
The authors declare no conflict of interest.

Receding Horizon Control
Because the RHC is already well-known, we recall below only the elements that define the RHC strategy:

•
The next controller's output is calculated by looking ahead for several steps in terms of a given cost criterion, but it is only implemented by one step; • The controller makes predictions for the number of steps taken into account (prediction horizon), using a dynamic PM; • The controller's output is sent to the process, and a new decision is made by taking updated information into account and looking ahead for a new prediction horizon.

•
The prediction horizon 'recedes' at the next sampling period but keeps its final extremity. This is the most difficult situation when the objective function has a terminal penalty term [17]. Hence, the prediction horizon length decreases by one unit at each sampling period.
The RHC structure is depicted in Figure A1, where the notations have the usual meaning. The objective function is denoted in a simplified form by J(u(k), k), where u(k) is the sequence of control input values over the prediction horizon that begins at the current moment k. The value u*(k) is the controller's output sent to the process at the moment k.
The receding horizon controller can use a metaheuristic algorithm to optimize the objective function at hand. A slightly modified EA is integrated into the closed loop in our case. Model predictive control is a special case of RHC when the prediction error is minimized to determine some controller parameters (see [18]). Only the RHC mechanism generates a quasi-optimal solution over the control horizon in this paper. Useful applications concerning the RHC structure are presented in papers [19,20].
meaning. The objective function is denoted in a simplified form by J( ( ) u k , k ), where is the sequence of control input values over the prediction horizon that begins at the rent moment k. The value u*(k) is the controller's output sent to the process at the mo k. The receding horizon controller can use a metaheuristic algorithm to optimiz objective function at hand. A slightly modified EA is integrated into the closed loop i case. Model predictive control is a special case of RHC when the prediction error is mized to determine some controller parameters (see [18]). Only the RHC mechanism erates a quasi-optimal solution over the control horizon in this paper. Useful applica concerning the RHC structure are presented in papers [19,20].