FPGA-Based Implementation of an Optimization Algorithm to Maximize the Productivity of a Microbial Electrolysis Cell

: In this work, the design of the hardware architecture to implement an algorithm for optimizing the Hydrogen Productivity Rate (HPR) in a Microbial Electrolysis Cell (MEC) is presented. The HPR in the MEC is maximized by the golden section search algorithm in conjunction with a super-twisting controller. The development of the digital architecture in the implementation step of the optimization algorithm was developed in the Very High Description Language (VHDL) and synthesized in a Field Programmable Gate Array (FPGA). Numerical simulations demonstrated the feasibility of the proposed optimization strategy embedded in an FPGA Cyclone II. Results showed that only 21% of the total logic elements, 5.19% of dedicated logic registers, and 64% of the total eight-bits multipliers of the FPGA were used. On the other hand, the estimated power consumption required by the FPGA-embedded optimization algorithm was only 146 mW.


Introduction
Nowadays, biotechnological systems represent a very attractive option for hydrogen production. The degradation of organic matter through the use of bacteria has gained great interest in the scientific community because hydrogen can be produced in a clean way [1,2]. In contrast to current industrial methods, in which unfortunately 90% of the hydrogen produced requires the use of fossil fuels generating a large amount of CO 2 (10 tonnes of CO 2 per ton of H 2 ) [3], Microbial Electrolysis Cells (MEC) represent a great alternative to produce hydrogen because they require less energy compared to the classic techniques to produce hydrogen, such as the electrolysis of water [4,5].
A MEC is an electrochemical device which uses electroactive microorganisms as catalysts to convert the organic matter to hydrogen and provides a novel approach for producing economically viable hydrogen from a wide range of renewable biomass sources [6,7]. Furthermore, a waste biorefinery based on MECs to produce clean and renewable electrofuel and valuable chemical compounds holds the flexible potentials for pollutants removal and CO 2 capture [8]. Broadly speaking, unlike a Microbial Fuel Cell, a MEC requires the induction of a constant voltage supply generating a potential difference between the electrodes to produce a flow of hydrogen as a result of the degradation of the organic matter that is fed to the MEC.
Other widely biological approaches used for the production of hydrogen in a clean way include Dark Fermentation (DF) in which bioreactors are fed by wastewater with a high concentration of organic matter from domestic and industrial origin. However, its efficiency to produce hydrogen compared to a MEC is relatively low (40% or less) [9]. best of the authors knowledge, there is not FPGA-based control implementations applied to bioprocesses.
In the present work the optimization problem of maximizing the Hydrogen Production Rate (HPR) in a MEC is addressed. The productivity function is approximated from the MEC model in steady state, for which, a point of maximum performance in a well-defined operating region is ensured. Using the golden section search optimization algorithm coupled to a robust super-twisting controller, the MEC is online brought to its maximum hydrogen production performance. The proposed optimization strategy is embedded in an FPGA throughout different digital architectures that are executed in parallel without hardware sharing. The resulting digital architecture has mainly two advantages, first, the portability to be synthesized in an FPGA card from any manufacturer, and second, the low power consumption compared to a personal computer. The implementation of the optimization algorithm in an FPGA has the great advantage of being described in hardware. This allows an easy adaptation in the use of communication protocols with external devices.
The rest of the paper is organized as follows: in Section 2 the mathematical model of the MEC is described, and the objective function as the HPR is presented. A description of the optimization problem is described in detail in Section 3. In Section 4 the optimization problem is addressed by using the Golden Section Search algorithm coupled to the discrete time super-twisting controller. In addition, the maximum HPR numerically computed is verified analytically. The FPGA-based implementation of the optimization algorithm is presented in Section 5 including numerical algorithms for the implementation of arithmetic operations like division, multiplication and square root. The results are presented in Section 6, where numerical simulations are carried out in an FPGA to verify the performance, including both the truncation error and the synthesis report of the digital architecture. Finally some conclusions are pointed out in Section 7.

Mathematical Model
One of the most used Microbial Electrolysis Cell (MEC) configurations currently consists mainly of two chambers that are separated by a cathode membrane (see Figure 1). In the anode chamber, the anode is covered by a biofilm where the existence of anodophilic and methanogenic bacteria is considered. The degradation of VFAs in the MEC takes place in the anode chamber, where hydrogen protons and electrons are produced. Protons pass through a ionic membrane to the cathodic chamber where the production of hydrogen occurs. A relatively small voltage is supplied to the system generating a potential difference between the two electrodes, which allows the electrons released in the anode by the anodophilic bacteria to circulate and pass to the cathode to combine with the hydrogen protons. In the degradation process there is a competition between two types of microorganisms, anodophilic and methanogenic, to decide who will consume the substrate. This behavior is modeled by the following system of Ordinary Differential Equations (ODEs) [23]: where s is the acetate concentration (mg/L −1 ), while x a and x m are the concentration of the anodophilic and acetoclastic methanogenic microorganisms, respectively (mg/L −1 ); D in is the dilution rate, where F in is the input flow rate (Ld −1 ) and V reac is the reactor volume (L); α a and α m are the dimensionless biofilm retention constants. µ a and µ m are the growth rates (d −1 ) for anodophilic and acetoclastic methanogenic microorganisms, respectively, which are defined as follows: where µ max,a and µ max,m are the maximum grown rates (d −1 ), k s,a and k s,m are the half-rate Monod constants (mg (s) L −1 ), F is the Faraday constant (C mol −1 e −1 ), R is the ideal gas constant (J mol −1 K −1 ), T is the temperature (K), η = E anode − E Ka is the local potential, where E anode is the anode potential (V) and E Ka is the half-maximum-rate anodic Electron Aceptor (EA) potential (V) i.e., the potential that occurs when S = k S,a and the rate is half of the maximum rate [24].

Anode Cathode
Power supply

MEC Productivity
The hydrogen flow rate in the MEC is modeled by Equation (6), where it can be seen that the hydrogen produced is closely related to the current generated from the flow of electrons between the electrodes.
where Y H 2 is the dimensionless cathode efficiency, A a is the anode area (m 2 ), m is the electrons per mol specie (mol e − mol −1 M −1 ) and P is the pressure inside the cathodic chamber (atm). In the Equation (6) the methanogenic microorganisms consumption is neglected and it is considered that only anodophilic microorganisms are responsible for acetate degradation. The current in the MEC is modeled as: where γ s and γ x (mFM −1 W −1 s ) are the yield coefficients related to the number of coulombs that it is possible to obtain from W s (g mol −1 ) and W x (g mol −1 ), i.e., the substrate, and the biomass respectively; f 0 s is the dimensionless fraction of electrons used for cell synthesis, b is the endogenous decay coefficient (d −1 ) and L f is the biofilm thickness (m).
The hydrogen production rate (HPR) Q H 2,p is defined as the hydrogen flow rate produced per volume of reactor (L[H 2 ] L −1 d −1 ): where Q H 2 is the hydrogen flow rate defined by Equation (6).

Problem Statement
The HPR is function of both, the dilution rate D in and the inlet acetate concentration s in . D in is the optimization variable, while s in is considered as a disturbance. As it can be seen in Figure 2, the HPR presents a maximum hydrogen productivity point related to an optimal dilution rate (Q H 2,p max , D in,opt ) within a range of concentrations for the inlet acetate s in [2000, 6000] mL −1 . Therefore, the optimization problem consists in calculating the value of the optimal dilution rate D in,opt that ensures the maximum performance Q H 2,p max in the MEC. Assumption 1. The function Q H 2,p is twice continuously differentiable in Γ with respect to D in such that: Assumption 2. The function Q H 2,p is convex, unimodal and any D in,opt is a global maximizer for each s in in the operating region.
The optimization problem to maximize the hydrogen production rate in the MEC is proposed as: such that: where x = [s, x a , x m ] T is the state vector, f (x, D in , s in ) is defined by Equations (1)-(5) and the measured output Q H 2,p (x) is the hydrogen production rate defined by Equations (6)- (8).
As it is shown in Figure 2, only a maximum Q H 2,p max can be observed for each maximizer D in,opt in the operating region.
The optimization problem is online solved by the GSS algorithm coupled to a supertwisting controller. The GSS algorithm calculates the value Q H 2,p max using a hydrogen productivity function in relation to both, the dilution rate and the inlet acetate concentration of the MEC. The super-twisting controller uses Q H 2,p max as a reference to track the MEC productivity to the maximum value. The optimization scheme described before is depicted in Figure 3.  In order to optimize the hardware resources and to reduce the power consumption, the optimization strategy to maximize the HPR of the MEC is embedded in an FPGA. This way, the energy cost required to bring the MEC to its maximum HPR can be considerably reduced.

Optimization of the MEC Productivity
with The objective function Q H 2,p (D in , s in ), defining the input-output map in steady state, is therefore expressed as: In this work, the acetate concentration in the inlet s in is assumed as known.

The Golden Section Search Algorithm
Golden ratio (ϕ) has been of a great interest to mathematicians, physicists, philosophers and artists. In antiquity, civilizations like Egyptians used the ϕ number as the main criterion for the construction of the Great Pyramids. The Parthenon in Greece was also built based on ϕ [26].
In relation to nature, the golden ratio is considered a natural constant that sets the standard in reproduction, growth patterns of living beings such as plants and animals. Their geometric relationship is described in Figure 4, where a line A-C is divided into two segments l and r by a point B where l is greater than r such that the ratio l/r is equal to the ratio (l + r)/l. The Golden Section Search (GSS) algorithm is an iterative process suggested to optimize one-dimensional, unimodal and well behaved functions [27], taking into account that the optimum value must be into a search region defined by a lower bound (A) and an upper bound (C), as described in Figure 4.
l r A C B The optimization of the HPR in the MEC begins defining the search region of Equation (15). In this case, the search region is defined by D in,A = 1 d −1 and D in,C = 3 d −1 . Then, two inner evaluation points D in,1 and D in,2 are selected as function of ϕ.
The error used by the GSS algorithm to stop its operation is defined as: The complete GSS algorithm to calculate the optimum point (D in,opt , Q H 2,p max ) is presented in Algorithm 1.

GSS Validation
First, let us analyze the stability of the nonlinear system (1)-(3) by calculating the eigenvalues (λ i ) of its linear approximation. The indirect Lyapunov method establish conditions that allow us to obtain conclusions about the stability of the nonlinear system in an operating point. Consider the nonlinear system (1)-(3) with the operating point x * = [s * , x * a , x * m ] that has the following linear approximatioṅ where x = x − x * , A, B u and B w are the Jacobian matrices of the system, u = D in − D * in and w = s in − s * in respectively. The indirect Lyapunov method states that the nonlinear system (1)-(3) is asymptotically stable if and only if Re(λ i ) < 0 of the matrix A, ∀λ i , i = 1, 2, 3, defined as: As it can be seen in Figure 5 the eigenvalues of the matrix A are Hurwitz in the operating region of the MEC. It must be pointed out that the closer the dilution rate is to the value 3 d −1 , the more the eigenvalues λ 1 and λ 2 approach the origin.

Algorithm 1: GSS algorithm description
The optimum value D in,opt is then obtained by differentiating the objective function (15) with respect to D in and equating the result to zero (first-order optimally condition), which leads to where   Figure 6 shows the Q H 2,p max value calculated both, by the GSS Algorithm 1 and by substituting D in,opt , calculated by setting the Equation (23) equal to zero (see Figure 7), in Equation (15). As it can be seen, the results of the GSS algorithm match the results obtained analytically.

Super-Twisting Controller
The MEC model (1)-(3) can be rewritten as follows: where γ(x) and g(x) are vector functions defined as: The relative degree σ of System (27) and (28) is computed by differentiating the output with respect to time as [28]: . Hence, the relative degree of the system (27) and (28) is σ = 1.
In this work the super-twisting controller, Equations (32) and (33), is therefore considered to track the maximum hydrogen flow rate computed by the GSS algorithm with the sliding variable defined as the tracking error [29].
In the super-twisting controller (32) and (33), the tracking error is defined as: ρ 1 and ρ 2 are the controller gains that ensure the finite-time stability of the tracking error, while D in,c is the control input necessary to bring the MEC to the maximum value Q H 2,p max . For implementation purposes in an FPGA, the discrete time super-twisting controller (DT-STC) is considered. The representative numerical solution showed in the Equations (35) and (36) is obtained from Equations (32) and (33) using the Euler's method. The controller uses the value Q H 2,p max as a reference to carry the real productivity to its maximum value in finite time.

FPGA-Embedded Optimization Algorithm
The FPGA-based implementation of the optimization algorithm is depicted in Figures 8 and 9. Following the scheme presented in Figure 3, the implementation block diagram is integrated by the GSS algorithm digital architecture coupled to the DTSTC digital architecture. A finite state machine (FSM) and a down programmable counter are used to ensure the proper operation of the optimization algorithm embedded in the FPGA.   x, x, x, x, 1 x, x, x, x, x x, x, 1, x, x x, x, x, x, 0 x, x, x, x, x x, x, x, 1, x x, x, x, x, x x, x, 0, x, x x, x, x, 0, x The digital architecture of the optimization algorithm uses a fixed point format (16,24) to represent all the input-output signals and inner operations. The hardware description used to develop the digital architecture was VHDL and the target board used was the Cyclone II EP2C35F672C6 integrated in the ALTERA DE2 educational board with a clock frequency f CLK = 50 MHz.
The modules GSS_MEC and ST_CONTROLLER were designed for an easy interaction with the FSM_CONT_MEC module and any other external device through the STG, EOG, STCS and EOCS signals. When the input signals STG and STCS are assigned to the logical value '1' by the FSM_VCONT_MEC module, they will produce a busy mode of their respective modules due to the latency time in the calculation of their final results. The busy mode is indicated by the output signals EOG ='0' and EOCS = '0'. On the other hand, when EOG = '1' and EOCS = '1', it means that the modules GSS_MEC and ST_CONTROLLER have finished and the results are ready to be read.

Operation of the FPGA-Embedded Optimization Algorithm
The FSM depicted in the Figure 9 is a great help for understanding the operation of the digital architecture. The FPGA execution can be divided in two steps, the initialization step, which is controlled by the states S0 to S2, and the normal operation, which is controlled by the remanding states of the FSM_CONT_MEC module. The initialization is executed when the FPGA is energized and the INI signal has a binary value '1' . Otherwise, the FPGA remains in standby mode until an external source changes the value of that signal. In such case, the initialization is started by a push-button (see the state S0). When INI = '1' the FSM changes to the state S1 where STG = '1' and SEL = '0' in the GSS_MEC module and the two-one multiplexer. This will start the calculation of Q H 2,p ,max with the initial value s in,0 . In the next clock cycle, the EOMEC signal in the GSS_MEC module will change from logic '1' to logic '0' indicating that this module is in the process of calculating Q H 2,p ,max . At the same time, without any condition, a transition is made to the state S2 where the FSM is waiting by the logic value '1' in the EOMEC signal indicating that the result is ready. When Q H 2,p ,max is ready to be used by the ST_CONTROLLER module, the FSM make a transition to the state S3 where the initialization step is done, and the system now is in the normal operation where SEL = '1' and it is waiting for an external device to set the value STOMEC = '1'. During the initialization step, the down counter is loaded with an initial value decreased by one every sampling period until reaching the optimization period.
In the normal operation, the ST_CONTROLLER module and the down counter are executed every sampling period with the aim of controlling the HPR in the MEC, and decreasing the initial value of the counter. When the down counter reaches the value zero, this means that the optimization period has expired and the GSS_MEC module is executed to generate a new Q H 2,p ,max , after that, the down counter is reloaded with the initial value.
The normal operation starts in the state S3 and the digital architecture reads s in by SEL = '1' in the multiplexer. When the signal STOMEC = '1', the FPGA-based optimization algorithm generates the control input D in,c of the MEC after a latency time, otherwise, the system is in standby. The execution of the ST_CONTROLLER and the down counter are managed by the states S5 to S7 in the FSM every sampling period, while the states S8 and S9 manage the GSS_MEC MODULE and the reinitialization of the down counter when the optimization period has been reached. In order to know when the GSS_MEC module should be executed, the FSM reads the signal Z from the down counter in the state S4. When Z = '0' this means that the optimization period has not yet elapsed and the FSM is currently executing the ST_CONTROLLER module, otherwise, when Z = '1' the FSM executes one more time the GSS_MODULE and generates a new Q H 2,p ,max in function of the current value s in . The down counter is reinitialized as well.
The most used arithmetic operations in the optimization algorithm are product, addition, division and square root. The hardware description was developed using standard VHDL and therefore the designs presented in this work do not belong to any manufacturer.

GSS Implementation
The digital architecture of the GSS optimization strategy, described in Algorithm 1, is depicted in Figure 10. The digital architecture of such algorithm is made up of registers, full adders, 8-bit embedded multipliers, multiplexers and full comparators using the previously mentioned fixed point format. Notice that the objective function shown in Equation (15) was programmed in the block Q H 2,p . Its implementation needed a simplified representation with the objective to calculate the hydrogen productivity with few hardware resources and small latency time. By precalculating constant parameters and making a separation by variables the following objective function is obtained: where the values of constants β 1 , β 2 and β 3 are defined in Table 1.  The complete comparator that determines if f 1 > f 2 , in Algorithm 1, was designed taking into account that the operation involves real numbers and therefore the classical definition of a complete comparator of binary numbers is not sufficient for this implementation.

DTSTC Implementation
The digital architecture of the DTSTC (see Figure 11) is simpler than that one of the GSS algorithm. Although only combinational elements are required, its response speed is quite fast to generate the control action compared to the speed of change to generate the reference computed by the GSS algorithm. The controller correction term σ 1 | | requires a digital circuit capable of computing the square root of the tracking error. Particularly in this work, the Pencil and Paper algorithm [30] proved to be very useful as a basis for the design of the SQRT arithmetic circuit.
The arithmetic circuit of the multiplier in the DTSTC architecture is based on the Coordinate Digital Computer Algorithm (CORDIC) with its rotating linear version (see Figure 12) [31], i.e.,: The results obtained after a sequence of fixed micro-rotations are given in the following way: The resulting operation y f in Equation (40) has the necessary shape to implement the DTSTC. As it can be seen in Equation (35), D in,c can be calculated from the final result y f by these two arithmetic operations; i.e., the product and the addition. The CORDIC-based Multiplier Digital Circuit presented in the Figure 12 has the shape necessary to implement DTSTC without the need of using embedded multipliers in the FPGA and it has a short latency time.

Results
The feasibility of the FPGA-embedded optimization algorithm was demonstrated through numerical simulations. The MEC model (1)-(6) was simulated in Matlab, the ODEs were solved by the stiff solver ode15s. The parameters used in the numerical simulations are listed Table 2. In order to demonstrate the robustness of the optimization strategy proposed, modified parameters between ±30% from their nominal value were considered. The hardware required for the verification test is depicted in the Figure 13. As it can be seen, a serial communication was used to communicate the FPGA with Matlab, which was executed in a personal computer with Windows 10, Intel Core i7 and memory RAM DDR3 of 32 GB. In these conditions, six hours were needed to perform the verification test of the optimization algorithm in a Cyclone II FPGA running at f clock = 50 MHz and a reception-transmission data rate of 70 Mbps. The operation time of the MEC simulated in the computer was of 200 d with a sampling period τ = 0.004 d. The hardware resources in the target board are summarized in Table 3.   The inlet acetate concentration s in used to feed the MEC in the numerical simulations is depicted in Figure 14. The digital architecture verification test of the MEC optimization algorithm consists mainly in comparing the results obtained from the FPGA working with the fixed point format (16,24) with the results of the same algorithm executed in Matlab in a floating point representation format. The resulting HPR obtained by executing the optimization algorithm both in the FPGA and in Matlab is shown in Figure 15. The green dashed-line represents the HPR by the MEC model, the red line represents the maximum HPR computed in Matlab, while the blue dashed-line represents the maximum HPR computed by the FPGA. On the other hand, the dilution rate computed by the optimization algorithm both in the FPGA and in Matlab is shown in Figure 16. The green dashed-line represents the optimum dilution rate computed by the GSS algorithm, the red line represents the dilution rate computed by the DTSTC in Matlab, while the blue dashed-line represents the dilution rate computed by the DTSTC in the FPGA. As it can be seen, the numerical representation format used to design the optimizer's digital architecture reduces properly the truncation error due to the finite number of bits. Initially, the optimization algorithm requires eighteen days to reach the optimal point, as shown in Figure 15. The super-twisting controller requires this period for the control error to converge to zero using the gains specified in Table 4. In this transitory period, the GSS algorithm is initialized with 105 mL[H 2 ] mL −1 d −1 and this value was taken as the initial reference for the DTSTC. Once the tracking error has converged to zero, the GSS algorithm reads the inlet acetate concentration value s in , every optimization period equivalent to D −1 in,max = 0.33 d to update the maximum productivity value Q H 2,p max used as reference by the DTSTC.
The acetate concentration in the MEC is showed in the Figure 17. It is easy to see in Figures 18 and 19 that the most of the acetate used to feed the MEC is consumed by the anodophilic bacteria x a because there is a inhibition process in the methanogenic bacteria growing x m . As expected, the current between the MEC electrodes is closely related to the HPR (see Figure 20).

Error Analysis
The truncation error in the digital architecture of the optimization algorithm is mainly due to the bits fixed quantity in the representation format established in this work. If the resolution in the intermediate operations required to run the optimization algorithm on the FPGA is not sufficient, the truncation error will propagate in such a way that the results obtained are greatly affected. Figures 21 and 22 show the behavior of the truncation error throughout the simulation process. It can be seen that the error is small enough to determine that the (16,24) format is sufficient to implement the optimization algorithm architecture in the FPGA.

Hardware Report
The FPGA hardware resources needed for embedding the digital architecture of the optimization algorithm on Figure 8 are summarized in Table 5. Only a 21% of the total logic elements (L.E.), 5.19% of dedicated logic registers (D.L.R.) and 64% of total eight-bits multipliers (8b-Mult.) in the chip Cyclone II were used. The input to output delay in the implementation was of 150 µs. The estimated power consumption required by the EP2C35F672C6 device using the aforementioned hardware resources is 146 mW. This estimate was calculated by the PowerPlay Early Power Estimator spreadsheet for Cyclone II family v8.0 SP1. The hardware resources used by the most important functional blocks in the optimization algorithm are summarized in Tables 6-8. As it can be seen in Tables 6 and 7, 64% of the total 8-b multipliers in the FPGA are used in the GSS algorithm, where 48.57% is destined to the Q H 2,p functional block where the objective function defined by Equation (15) is processed. It should be pointed out that the Q H 2,p block is part of the GSS algorithm functional block (see Figure 10). The GSS algorithm needs at least 4 cycles in the worse of the cases to reach the tolerance error (err = 0.001) defined by Equation (20). Therefore, embedded multipliers most be used in the GSS algorithm digital architecture to have a short latency time. On the other hand, the hardware resources used in the DTSTC and its inner functional block, the CORDIC Multiplier, are summarized in Tables 8 and 9. Most DTSTC inner operations are implemented using a CORDIC-based multiplier that has a latency time of 1.48 µs in the worse of the cases, before the tracking error converges to zero. After that, the multiplier is executed faster than 1.48 µs. It should be noted that the CORDIC-based multiplier internally uses an 8-bit expansion in the fractional part to substantially improve the truncation error generated by the fixed-point format considered (see Figure 22). Finally, the arithmetic operation | | in the DTSTC is processed by the SQRT functional block, which is based on the Pencil and Paper algorithm. Its digital architecture is primarily based on bit additions and shifts. Table 10 shows the hardware resources needed.

Conclusions
In this work an FPGA-embedded optimization algorithm to maximize the hydrogen production rate (HPR) of a microbial electrolysis cell (MEC) using the golden section search (GSS) algorithm coupled to a discrete-time super-twisting controller (DTSTC) was presented. The correct performance of the GSS algorithm was analyzed analytically. Furthermore, it was proven that the relative degree of the MEC model is one, a necessary condition to use the DTSTC to bring the HPR to its maximum performance point in finite time.
To reduce the power consumption required to bring the MEC to its maximum performance, a digital architecture of the optimization algorithm was designed and embedded in an FPGA. Although the FPGA used in this work was the Cyclone II of ALTERA, the digital architectures presented in this work were designed to be implemented in any FPGA, regardless of the manufacturer.
The results of the FPGA-embedded optimization algorithm showed a correct performance with low hardware resources and low power consumption compared with a personal computer. Besides, the truncation error generated by the fixed point format used in this work was practically negligible.
Such results allow us to conclude that the implementation of control and optimization algorithms in FPGAs represents an excellent alternative to replace personal computers. Particularly, as demonstrated in the previous section, the FPGA-embedded optimization algorithm proposed to maximize the HPR in the MEC, represents a lower cost alternative in terms of consumed power and resources.