Multi-Timescale-Based Partial Optimal Control of a Proton-Exchange Membrane Fuel Cell

: This paper presents a Proton-Exchange Membrane Fuel Cell (PEMFC) transient model in stack current cycling conditions and its partial optimal control. The derived model is used for a speciﬁc application of the recently published multistage control technique developed by the authors. The presented control-oriented transient PEMFC model is an extension of the steady-state control-oriented model previously established by the authors. The new model is experimentally validated for transient operating conditions on the Greenlight Innovation G60 testing station where the comparison of the experimental and simulation results is presented. The derived ﬁve-state nonlinear control-oriented model is linearized, and three clusters of eigenvalues can be clearly identiﬁed. This speciﬁc feature of the linearized model is known as the three timescale system. A novel multistage optimal control technique is particularly suitable for this class of systems. It is shown that this control technique enables the designer to construct a local LQR, pole-placement or any other linear controller type at the subsystem level completely independently, which further optimizes the performance of the whole non-decoupled system.


Introduction
Proton-Exchange Membrane Fuel Cells (PEMFC) are devices that convert the energy of electrochemical reactions to electrical energy. PEMFC is seen as a promising power source of the future, as an alternative to internal combustion engines and batteries, mainly because the only byproduct of the electrochemical reaction is water. Wide commercialization of such devices could reduce the greenhouse gas emissions, particularly associated with transportation. Those devices have a promising future in the automotive applications [1,2]. Even though the PEMFC have been extensively studied, they are still not even close to concur the market of power sources. The major bottlenecks for wider commercialization of PEMFC have cost and durability nature. Such systems operate sufficiently well on a fraction of their lifetime and then the performance decreases. Major reason for this is limited understanding of complex phenomena occurring during PEMFC operation and not controlling the system sufficiently well. There is an emerging need for taking a constructive control-oriented approach to PEMFC modeling and developing new control strategies.
PEMFC modeling and control is a multidisciplinary research area, because successful and efficient PEMFC operation implies deep understanding of various electrochemical, transport and material science phenomena. Over the years PEMFC have been studied from various aspects. Quite broad ongoing PEMFC research may be divided into two groups: (i) Examination of the fundamental electrochemical processes happening during PEMFC operation.
PEMFC modeling is the research area which represents the coupling between (i) and (ii). Physical PEMFC modeling uses the fundamental results of the group (i) research and lights the path which leads to a better design of cells, stacks, as outcome of the group (ii) PEMFC research. Numerous experimental, as well as theoretical studies have been conducted with a goal to explore the system's performance and enhance efficiency. Review of the models present in the literature can be found in [17,18]. Some researchers focused on the system design and modeling the water content present in the PEMFC stack [19]. Study [20] focused on the effect of membrane water content during cold start. Studies such as [21] carried out a 3D non-isothermal multi-phase modeling and analysis. Common approaches in modeling are mechanistic approaches and empirically-based analysis [7,[22][23][24][25]. In general, mechanistic approaches require specific parameters identification and empirically-based analysis relies heavily on experimental data which may be valid in a very narrow operating range. Also, those systems seem to be overparametrized [26], complicated and very challenging to control. In the current literature, to the authors' best knowledge, there is not present a transient control-oriented model, suitable for further study and control systems design. Such a model is derived and experimentally validated in this paper.
PEMFC systems age with exploitation and performance of the stack deteriorates [27][28][29]. One of the main reasons that causes the aging effect is the oxygen starvation phenomenon. Study [30] summarizes about 15 possible reasons for gas starvation and suggests that more research effort should be put in the diagnostics of starvation and adequate control in the operating regime. Furthermore, PEMFC material degradation, and local heterogeneity, caused by nonuniform fueling, were discussed in [31], where some possible mitigation strategies were suggested. For example, study [31] states that hydrogen starvation problems could be reduced by shortening the hydrogen presence in front of anode bipolar plate channels, which will help in uniform distribution along the stack channels. Stack behavior under fuel starvation conditions is explored and experimentally documented in [32]. PEMFC assembly aging was experimentally shown [33], where analysis has shown that frequent gas starvation operating conditions present especially during transients, have significant impact on performance. For successful PEMFC operation there is a need for an adequate control strategy and the controller design. Comprehensive studies [17,34] review most recent research efforts and attempts to define and use of various control techniques for PEMFC systems. Despite all the efforts, the community needs more control-oriented models and classical control schemes such as [35,36]. On the other hand, physical variables (states in the control-oriented state space models) coupling and dependence among them is not explored well enough. Dynamic phenomena present in the PEMFC have been studied on the improved model presented in [37], where it is shown that they have an important influence on the PEMFC fuel cell control system design. Coupling analysis showed that some sort of the decoupling control strategy would optimize the system performance. The recent approach presented in [38][39][40], tackle the coupling interference between physical variables represented as states in the model, by investigating PEMFC subsystems as systems with multiple time scales. It was suggested and demonstrated in [41] that PEMFC has at least three distinct characteristic time scales, processes with response time less than 1s, approximately 100s and more than 1000s. This unique feature of PEMFC systems makes them particularly suitable for the application of the multistage multi-timescale controller technique.
This paper provides a control-oriented framework to study PEMFC and it is consisted of two parts:

Experimental Setup
An experimental setup present in a Villanova University lab is consisted of: PEMFC testing station Greenlight Innovation G60 ( Figure 1) and a Greenlight Innovation TP50 fuel cell stack ( Figure 2). Model derivation in this paper assumes that the stack temperature during transients is constant (stack endplate heaters and external fan ensure operating temperature), so the thermal control loops are purposely omitted on the block diagram. Block diagram in Figure 3, shows only flow loops present in the system: fuel (Hydrogen) flow loop, oxidant (Air) flow loop, purging gas (nitrogen) flow loop, deionized (DI) water flow loop. As it can be seen input gases are driven through external humidifiers prior to entering the PEMFC stack.   In the model derivation the following assumptions were made:

•
Inlet pressure is maintained constant during transient periods.

•
All gases present are assumed to be single phase gases.

•
Flow velocity is the same on the cathode (CA) and anode (AN) sides at all points where the gas flow is considered laminar, despite the fact that on the anode side during low stoic and pure H 2 condition this assumption might be violated. Since this assumption is used only for pressure calculation, non-uniformity of the flow on the anode side can be neglected in the developed model.

•
Internal stack relative humidity never exceeds 100(%). In this model, liquid water content (known as flooding) is not considered.

•
CA and AN volumes are calculated as single lumped volume.

•
Stack temperature is constant.

Control-Oriented Transient PEMFC Model: Derivation
The control-oriented PEMFC model is structured as In the rest of the paper, for convenience, explicit input, output, disturbance, and state time dependency will be omitted. The state space vector has five states, three in cathode (CA) and two in anode (AN), as outlined in Figure 4. CA and AN volumes are lumped volumes of serpentine distributed flow channels in CA and AN. The state vector is defined as  Particular state equations take the following form CA: CA: AN: AN: CA: where partial species present in CA or AN mass flow are denoted as W with subscripts in, out, rct, gen, mem stand for inlet and outlet mass flows, reacted (depleted in electrochemical reaction) mass flow, generated species mass flow and mass flow through membrane respectively. Schematic of the mass flows present in the equations are provided in Figure 5. The derivations are procedure is split in four parts cathode, anode, outlet pressures and stack voltage.

Cathode
This part focuses on Equations (3), (4) and (7). Particular mass flows considered are given in Figure 6. Cathode inlet components, mass flows of oxygen W O 2 ,in , nitrogen W N 2 ,in , and water vapor can be represented as mass fractions of the total CA inlet mass flow W a,ca,in .

MEA
where M a and M v are air and water vapor molar masses. Mass fraction coefficients x O 2 and x N 2 are defined as where y O 2 = 0.21 is the oxygen mole fraction of the incoming air and M a , M N 2 are molar masses. In (10) externally controlled CA inlet pressure p ca,in is assumed to be constant, while the partial water vapor pressure depends on the incoming gas relative humidity φ ca,in and saturation pressure p T ca,in sat calculated at cathode inlet temperature T ca,in [42] as The saturation pressure defines an upper limit of water vapor partial pressure in the water-gas mixture. If this pressure is exceeded, condensation will occur. In transient conditions, delay between the stack current change and the net mass flow of oxygen W O 2 ,net can be represented by the first order differential equation where the net flow is the difference between inlet and outlet oxygen mass flows, Coefficient τ ca is a fitting parameter, with usual value of a couple of seconds. Next, CA outlet components are defined as where m ca = m O 2 + m N 2 + m v,ca is the total mass present in the CA. Total outlet CA mass flow can be determined from W ca,out = k ca (p ca − p ca,out ), where k ca is the empirically determined orifice constant. The total pressure in the cathode is sum of all components partial pressures p ca = p O 2 + p N 2 + p v,ca , and p ca,out is the CA outlet pressure defined in subsection outlet pressures. Partial pressures at stack temperature T st present in the cathode can be calculated by using gas constants R O 2 , R N 2 and CA volume V ca as: where φ ca is the relative humidity in CA. Relative humidity inside CA depends on the state m v,ca as where M v is the water molar mass and R is the universal gas constant. Due to chemical reaction some portion of the oxygen is depleted on the CA side (represented by −W O 2 ,rct in (3)) and some water vapor is generated. For stack current I st and n cells present in the stack, the use of Faraday's Law [43], leads to where F is the Faraday's constant. Two phenomena dominate the membrane transport of water vapor: electro-osmotic drag and water back diffusion. The electro-osmotic drag occurs when H + ions dragged through membrane carry some water molecules as well. Also, due to new water vapor molecules excessive generation on the CA side a pressure gradient is present between CA and AN, so some molecules are being transferred back from CA to AN side. This process is called water back diffusion. On higher current densities, in a typical PEMFC working regime, water back diffusion is dominated by electro-osmotic drag, so the resulting net mass transfer of water vapor molecules is from AN to CA side. Mass transfer through membrane with active area A f c for stack with n cells and current density i has the following form [44] W where water concentrations in mol . Water content coefficients λ ca , λ an and λ m are fitted from experimental data taken from [45]. They are dependent on water activities a ca and a an , which are calculated as ratio of the partial water vapor pressure p v,ca or p v,an and the saturation pressure at stack temperature as a ca = p v,ca where: k = ca, an, m (average water content) In Equation (21) electro-osmotic drag coefficient n d is defined as in [46]. It is a polynomial function of average water content λ m Water diffusion coefficient D w in Equation (21) is given in [44] D w = D λ e where

Anode
This modeling part focuses on Equations (5) and (6). Particular mass flows considered are given in Figure 7. Anode inlet components, hydrogen mass flow W H 2 ,in and water vapor mass flow W v,an,in are where M H 2 is hydrogen molar mass. In (27), AN inlet pressure p an,in is controlled externally, by the station and assumed to be constant during transients. Partial water vapor pressure depends on the incoming gas relative humidity φ an,in and saturation pressure p In transient conditions, delay between the stack current change and the net mass flow of hydrogen W H 2 ,net can be represented by the first order differential equation where the net flow is the difference between inlet and outlet hydrogen mass flows, W H 2 ,net = W H 2 ,in − W H 2 ,out . Coefficient τ an is a fitting parameter, with usual value of a couple of seconds. The outlet hydrogen mass flow W H 2 ,out and the outlet water vapor mass flow W v,an,out are mass fractions of the total mass in AN, m an = m H 2 + m v,an , given as W v,an,out = m v,an m an W an,out where the total outlet mass flow W an,out is proportional to pressure difference as W an,out = k an (p an − p an,out ), between total AN pressure p an and AN outlet pressure p an,out defined in subsection outlet pressures. AN orifice constant k an is determined empirically. The total AN pressure is sum of partial hydrogen pressure p H 2 and partial water vapor pressure p v,an , p an = p H 2 + p v,an , where partial pressures are given as p H 2 = R H 2 T st V an m H 2 and p v,an = φ an p Tst sat , with the hydrogen gas constant R H 2 and AN lumped volume V an . Relative humidity inside AN, φ an is a function of the state variable m v,an as In chemical reaction on the AN side some hydrogen is depleted, which is represented by term −W H 2 ,rct in (5) and defined by Faraday's Law [43] as The water vapor mass flow through membrane W v,mem in Equation (6) is explained previously, but here has a negative sign −W v,mem , because the direction of the flow is from AN side to CA side.

Outlet Pressures
Greenlight TP50 PEMFC present in the experimental setup has a serpentine flow pattern. This pattern causes significant pressure loss as the gas is traveling through the channels. Pressure difference between inlet and outlet pressures p ca,in and p ca,out on the CA side and p an,in and p an,out on the AN side is caused by three dominant losses: • Friction losses. As flow is considered to be laminar at all times, this loss ∆P can be modeled by Darcy's Law as where Darcy's friction factor is given as f d,j = 64 R e j , with Reynolds number R e j = ρ j L j µ j v j and flow velocity v j = ρ j A c W j,in , which depends on the inlet mass flow W j,in . In the previous expressions ρ j is fluid density, µ j is fluid viscosity, L j is channel chord length, D H is hydraulic diameter of the rectangular flow channel, and A c stands for the channel cross sectional area.
• Channel bending losses. This pressure loss ∆P bend,j is a consequence of bends present in the serpentine channel pattern [47], dependable on flow velocity v j where k j = k loc + k f r,j are friction coefficients defined as k loc = 0.21 • Chemical reaction pressure losses ∆P rct dependable on flow velocity v j : Coefficient values can be found in Appendix A, Tables A2 and A3.
Subtracting all three defined losses from inlet pressures p j,in leads to formula for outlet pressures p j,out on the CA side (j = ca) and the AN side (j = an)

Stack Voltage
The actual stack output voltage consisted of n cells is the open circuit voltage E cell reduced by three quantities, in this work considered, voltage losses: [48], or ohmic voltage loss v ohm , activation voltage loss v act , and concentration voltage loss v conc By following the principles of thermodynamic potentials [42,49], it is known that a negative change in Gibbs free energy directly corresponds to the maximum electrical work, seen as an open circuit voltage, at standard state conditions. Under nonstandard state conditions, the open circuit voltage E cell is a function of variations of pressures of the reactant species, in PEMFC case, oxygen pressure p O 2 and hydrogen pressure p H 2 and it is represented by the Nernst equation as where ∆ŝ is the heat capacity change of the substance. At stack temperature T = T st ,with standard temperature T 0 = 298. 15[K] and with ∆ŝ taken from thermodynamic tables in [49], the final form of the open circuit voltage equation is For a membrane thickness t m and membrane conductivity σ m , ohmic voltage loss v ohm is Membrane conductivity takes the following form [45] where λ m is water content provided in Equation (22). Portion of the energy is consumed in the chemical reaction, where it is used to support chemical bonds breakage and electron transport. This energy consumption is the main reason for activation voltage losses v act . This loss can be described by Tafel where v 0 is the voltage drop at zero current density, v a values are taken from [44] and c 1 = 20.8605 is fitted to our experimental data. The voltage drop at zero current density is in the form The concentration voltage drops dominate the voltage response of the stack at higher current densities. At stack current density i the voltage concentration drop takes the form where c 2 is taken from [44] while parameters c 3 = 4.4469 and i max = 0.96 are fitted to match experimental data. Membrane Electrolyte Assembly (MEA) separates metal plates in space and allows positive H + ions to pass while electron flow is blocked. On the AN side there is a negative charge build up and on the CA side there is a positive charge build up. In this way, an electrochemical double layer is formed, and it is acting as a regular electrical capacitor, which can store some electrical energy. This phenomenon is incorporated in the model by introducing the equivalent electrical circuit [50], where equivalent electrical capacitance C needs to be determined. The stack voltage becomes where governing equation for capacitor voltage V C for circuit provided in Figure 8 is where R act and R conc are equivalent resistances of the activation voltage loss v act given in (43) and the concentration voltage loss v conc given in (45).

Multistage Control Technique
The multistage control technique applied to a linear system in state space form can be split in two parts [40]: part 1 of the design transforms original system into the upper triangular form, and part 2 implements the multistage procedure where in each stage reduced order subsystem is decoupled from full-order system and a local controller is designed. A linear time invariant (LTI) continuous time dynamic system in a three-system partitioned form looks like   ẋ where state vector variables are x(t) ∈ R n×1 , x 1 (t) ∈ R n 1 ×1 , x 2 (t) ∈ R n 2 ×1 and x 3 (t) ∈ R n 3 ×1 , n = n 1 + n 2 + n 3 , the control input vector u(t) ∈ R 1×1 . State space matrices A and B are constant and submatrices A ij and B ij have appropriate dimensions, i, j = 1, 2, 3. The diagonal blocks of matrix A, submatrices A ii define subsystems x 1 (t), x 2 (t), x 3 (t), with coupling between them determined by the non-diagonal elements A ij , i = j. In the following presentation explicit time dependency of the state and input variables is omitted to simplify notation. Part 1 derivation steps explained in [40] lead to a similarity transformation L, which maps the original system (49) in x coordinates to the upper triangular form in η coordinates as outlined bellow where transformation matrices T 2 , T 1 are defined as L and P transformations are similarity transformations, so when the system is transformed back to original coordinates x, the system eigenvalues are preserved. Carrying out derivations can be cumbersome and might blur the outline of the multistage design. For a quick reference, those equations are summarized in a table [p. 194] [40]. Solving nonlinear algebraic equations for L 1 , L 2 is a very challenging task, if possible at all. To overcome this issue, this approach is specialized for three timescale singularly perturbed systems in the form with small positive singular perturbation parameters ε 1 ε 2 ≥ 0. MatricesÂ ij andB i1 are singularly perturbed matrix elements corresponding to original ones A ij and B i1 , i, j = 1..3. State vector variables x ∈ R n×1 , n = n 1 + n 2 + n 3 are split as slow state variables: x 1 ∈ R n 1 ×1 , fast state variables: x 2 ∈ R n 2 ×1 , and very fast state variables: x 3 ∈ R n 3 ×1 . For singularly perturbed system (53), resulting governing equations can be solved by iteration algorithm [40]. L similarity transformation (50) coefficients L 1 , L 2 and L 3 can be determined from the following iterative equations where L f 1 and L f 2 are final iteration values determined in the first two equations. P similarity transformation (52) is determined by solving for P 1 , P 2 and P 3 coefficients, from the following iterative equations [40] where L f 3 , P f 2 are final iteration values obtained from previous equations. TermsÂ 1 ,Â 2 ,Â 3 ,B 2 ,B 3 , α 12 ,α 13 andβ 2 are spelled out in the Appendix A.2. Comments on the required computational power needed to solve those iterative equations and achieved precision and convergence rate are provided in [39].

Algorithm 1: Three stage technique implementation
input : Singularly perturbed matrices A sp ,B sp , eqn. (53) output : Equivalent gain for feedback G eq solve for L 1 , L 2 matrices (54) ; solve for L 3 matrix (54) ; set G 3 control gain law with very fast desired eigenvalues λ d very f ast ; solve for P 3 matrix (55) ; set G 2 control gain law with fast desired eigenvalues 1 ε 1 λ d f ast ; solve for P 2 matrix (55) ; solve for P 1 matrix (55) ; set G 1 control gain law with slow desired eigenvalues 1 ε 2 λ d slow ; set equivalent gain G eq = G 1 G 2 G 3 T f 2 T f 1 L f ; → feedback control: u = −G eq x ; The presented approach assumes that the following matrices are non-singular (invertible): • MatrixÂ 33 . This is a standard assumption in singular perturbation analysis [51,52].
If this cannot be satisfied then L 3 equation in (54) can be solved iteratively as a sequence of algebraic Sylvester like equationŝ Those matrices are asymptotically stable, so they are always invertible [54]. Feedback controller gains G 3 and G 2 are chosen such that this is satisfied.

Control-Oriented Transient PEMFC Model: Experimental Validation
Nonlinear model simulation is run for stack current step from I st = 25A to I st = 30A and compared with experimental data-log for the same stack current step. Model transient response provides a good match overall with the experimental data, as observed in Figures 9 and 10.

Multistage Control Technique Application: Simulation Results and Analysis
The developed five-state nonlinear model is linearized around nominal operating point with the following nominal inputs and disturbance values  The corresponding linear model takes the following standard state space form The transient model is structured to have two inputs, but the multistage control technique is applicable to single input systems. A dual input system can be reduced for simulation purposes to a single input system as δẋ = Aδx + B u 1 δu 1 δy = C y δx + D u 1 δu 1 (59) System input u 2 is calculated from relation u 2 = c f it · u 1 , where c f it is the fitting constant. Simulation results in Figures 11-13 are obtained from the following model Small perturbation parameters are defined as Permutation matrix V puts the original system matrix A (given in the Appendix A) in the explicit singularly perturbed form provided in (53), where A sp = VAV T . One randomly chosen permutation matrix V solution has the following form Solving (55)

Control-Oriented Transient PEMFC Model
Control-oriented state space PEMFC model development enables us to better understand the PEMFC modeled dynamics, get insight what are the operating challenges and to do in depth analysis. Model is set up in a standard state space form so that wider engineering and scientific community can use it as a framework for design and implement new control techniques that may lead to more efficient and more durable PEMFC systems. The control-oriented model features to be highlighted are: • A complex physical system is represented by simple, but a realistic model, without overwhelming details. Following this modeling approach, scientist equipped with such a model can have full focus on designing a control strategies and controllers for the most dominant effects and study them theoretically and experimentally with less efforts. • In this model setup stack current is treated as a measured disturbance. The stack current is the fastest physical variable and it can be used in various control designs, with accent on feedforward control approaches, which one possible future research direction.

•
Since the inlet relative humidity is included in the model, it can be manipulated (controlled) to improve the water management in the stack, which is crucial for the correct operation of the fuel cell.

•
The presented logic model can be easily expanded. Flexibility is important, because for example, if a researcher wants to focus on the analysis in the higher current densities operational range, one might want to consider adding extra states which represent liquid phase generated water.

•
The system is in state space form, which enables us to perform controllability and observability analysis.
Future work: Future research directions would be derivation of aging dynamics and incorporating it in the presented model with the intention of developing a useful diagnostic tool for PEMFC stack health monitoring. Possible future work includes integration of the developed PEMFC stack model with the hydrogen gas processing system.

Multistage Control Technique Application
With the proposed multistage control technique developed and elaborated in detail in [40], different types of state feedback controllers may be designed to control different subsystems. These local feedback controllers can stabilize certain subsystems of interest only, while the rest of the subsystems are seamlessly not affected at all. The notable features of the proposed technique are:

•
Similarity transformations L and P preserve the controllability/observability property of the full system. If the original full-order system is controllable/observable then local systems will also be controllable/observable.

•
The original full-order system is decoupled in several subsystems, each with the unique timescale. A composite feedback controller applied to the full-order system is constructed with different state feedback controllers: eigenvalue assignment, optimal, robust, etc.

•
Computational costs are significantly reduced because all calculations are done with reduced order matrices where numerical ill-conditioning of higher order matrices is completely avoided. Instead, very high accuracy can be achieved by using iteration algorithms to solve nonlinear algebraic equations. • Partial full-state feedback and partial output feedback (for subsystem for which only the output is available for feedback) can be implemented through this design. If certain subsystem does not need to be controlled, the gain for that subsystem is set to zero. The gains of the other subsystems are set accordingly and combined in the composite feedback matrix.
The multistage control technique is a state feedback technique, so the states need to be available for feedback. The PEMFC derived system model has states which are not measurable, but observable and the state space form enables us to do the observability analysis. Due to multiple time scales nature it is a great challenge to design an adequate observer to estimate unavailable states, if the system is observable at all. Research in that direction, where linear observers for systems with slow and fast modes are considered, could be found in [55]. On the other hand, if multistage design technique is used, then observers could be constructed at subsystems level, which simplifies the design. In the review paper [30], some possible mitigation strategies are suggested: system intake control, nitrogen governing, purging strategy, etc. This approach in controlling large scale systems with multiple time scales is applied to nuclear reactors [56]. Future work: Our group specialized this design for 2 timescale systems [57] and expanded the work to 3 timescale systems [39,58]. The multistage design framework is quite universal and future research directions could include the possible expansion to 4 timescale systems, which is already outlined in [59]. Also, it will be interesting to study how observers affect optimality when optimal controllers are used at subsystem levels [60].
Author Contributions: The authors contributed equally to this work. All authors have read and agree to the published version of the manuscript.