Enhancing the efficiency and lifetime of a proton exchange membrane fuel cell using nonlinear model-predictive control with nonlinear observation

Enhancing the Efficiency and Lifetime of a Proton Exchange Membrane Fuel Cell using Nonlinear Model Predictive Control with Nonlinear Observation Julio Luna, Student Member,IEEE, Elio Usai, Member,IEEE, Attila Husar, and Maria Serra Abstract—The aim of this research is to develop and test in a simulation environment an advanced model-based control solution for a Proton Exchange Membrane Fuel Cell (PEMFC) system. A Nonlinear Model Predictive Control (NMPC) strategy is proposed to maximise the active catalytic surface area at the Cathode Catalyst Layer (CCL) to increase the available reaction area of the stack and to avoid starvation at the catalyst sites. The PEMFC stack model includes a spatial discretisation that permits the control strategy to take into account the internal conditions of the system. These internal states are estimated and fed to the NMPC via a Nonlinear Distributed Parameters Observer (NDPO). The air-fed cathode of the PEMFC simulation model includes a two-phase water model for better representation of the stack voltage. The stack temperature is regulated through the use of an active cooling system. The control strategy is evaluated in an automotive application using a driving cycle based on the New European Driving Cycle (NEDC) profile as the case study. Index Terms—Electrochemically active surface area, nonlinear model predictive control, nonlinear observation, proton exchange membrane fuel cells, degradation, starvation. NOMENCLATURE Throughout this paper, spatially distributed systems are treated, denoting the spatial variables as x,yand z. Subscripts iand jare associated to the reactant and discretisation volume respectively. For instance, ci,j refers to the concentration value of the i-th gas at the jdiscretisation volume. Column vectors are denoted by bold style, e.g. x. Matrices are denoted by bold upper case, e.g., A. Scalars are denoted by non-bold style, This work has been partially supported by the Spanish national project MICAPEM (ref. DPI2015-69286-C3-2-R, MINECO/FEDER), the Regione Autonoma della Sardegna project CRP-7733 and the Fuel Cells and Hydrogen 2 Joint Undertaking under grant agreement No 735969. This Joint Undertaking receives support from the European Union’s Horizon 2020 research and innovation programme and Hydrogen Europe and N.ERGHY. J. Luna, A. Husar and M. Serra are with the Institut de Rob` otica i Inform` atica Industrial (CSIC-UPC), Llorens i Artig` as 4-6, 08028 Barcelona, Spain. E-mail: {jluna, ahusar, maserra}@iri.upc.edu E. Usai is with the Department of Electrical and Electronic Engineering (DIEE), University of Cagliari, Cagliari 09123, Italy. E-mail: [email protected] e.g., B. The set of real numbers is denoted by R. Observed variables are denoted by the caret symbol, e.g., ˆx. I. INTRODUCTION AS energy consumption increases, society, industry and governments have become aware of the necessities to invest in sustainable energies that can decrease the problems associated with the use of fossil fuels and nuclear energy. Recent studies [1] show that the use of hydrogen as an energy vector can aid to satisfy the present and future energy demands without additional carbon emissions. In this context, Proton Exchange Membrane Fuel Cells (PEMFCs), which use hydrogen as fuel and provide high power densities while operating at low temperatures, are one of the most promising technologies for both stationary and mobile applications. To guarantee the optimal operation of PEMFC-based systems when designing novel control strategies, Balance of Plant (BoP) auxiliary components such as compressors, pumps, heat exchanger, etc., have to be taken into account. To compete with other power generation systems, such as internal combustion engines for automotive applications, PEMFCs have to achieve a similar efficiency and cost. Cost reduction by means of materials improvement has already been achieved during the last decade [2]. Nevertheless, there is still room for improvement regarding the efficiency and durability of PEMFCs. Efficiency and durability are associated with the operating conditions of the system which are subject to changes due to cycling and current demand. Furthermore, the internal conditions of the PEMFC also affect the performance and durability of the system. Lifetime of PEMFCs is related mainly with catalyst degradation, caused either by Platinum (Pt) dissolution or carbonsupport corrosion [3, 4]. During the normal operation of a PEMFC, degradation can mainly occur due to three mechanisms: baseline degradation, cycling degradation and incidentinduced degradation. In hybrid systems, the lifetime of the fuel cell is also affected by the power distribution between the battery and the fuel cell. However, hybridisation and its effects on the PEMFC durability are out of the scope of the present research. Designing proper control strategies can reduce the degradation rate of the PEMFC through the use of the available manipulable inputs of the system to avoid starvation at the Catalyst Layers (CLs). Moreover, it is possible to mitigate the effect of unavoidable degradation mechanisms and improve the PEMFC efficiency by means of proper water management in the CLs to maximise the available active surface. Quantifying degradation is a challenging task. An approach proposed in the literature [5] is to model the effective area where the reaction can occur. This area is known as Electrochemically active Surface Area (AECSA). The AECSA is a measure of the total active Pt available in the carbon-support layer at the Cathode CL (CCL) and it depends on the Pt loading of the CCL, the pore distribution, the CCL hydration state and the degradation condition of the stack. To maximise the AECSA, the only available manipulable variable is the hydration state of the system, which can be actively controlled mainly by modifying the temperature of the stack and the inlet cathode Relative Humidity (RH). The in-situ characterisation of AECSA has improved in the last years [6, 7]. However, determining its value while the system is being operated is not yet possible with the current technology. In this sense, modelling and estimating AECSA is an important step forward to actively control this parameter. Internal conditions greatly affect the performance and degradation of PEMFCs. However, most of the published control solutions have been based on models that do not consider the spatial dynamics of the PEMFC [8], providing an overly simplified lumped description of the system when advanced control strategies need to be designed. Existing sensor technology is not capable of measuring internal variables due to the enclosed nature of the system. In the literature, modelbased observation with nonlinear distributed models has been proposed to tackle this drawback [9]. Regarding control strategies, the range of control techniques used in PEMFC-based systems is wide: unfalsified controllers, predictive controllers and variable structure controllers are some of the most used control strategies as analysed in recent review works [10, 11]. As any real system, PEMFCs have plenty of fast dynamic behaviours and variables bounded by physical limits that should be considered when designing a control law. Moreover, the definition of several operational constraints, in the same way as variable bounds, should be taken into account when formulating a closed-loop control scheme. In this sense, Nonlinear Model Predictive Control (NMPC) [12] is a promising option since it has the ability to handle state and input constraints. Moreover, NMPC is able to deal with the nonlinearities that are present in fuel cells [13]. An additional advantage of NMPC is its intrinsic capability of considering multiple manipulable variables and control objectives as a single multi-objective control problem. However, NMPC requires a reliable prediction model and output-feedback of the current system state, which includes an additional computational burden to the controller. The main contribution of this paper relies on the combination of a NMPC strategy that considers the estimation of the AECSA along the CCL in the direction of the gas channels. To achieve this, a nonlinear distributed parameters model of a PEMFC [14] and its BoP auxiliaries is implemented along with a Nonlinear Distributed Parameters Observer (NDPO) to estimate the unmeasured states that are injected into the NMPC. The control strategy uses a prediction model to optimise the operating point of the fuel cell, computing a set of optimal control actions for a given cost function at each time instant. The cost function is selected to guarantee the fuel cell performance in terms of AECSA maximisation (and therefore, fuel cell efficiency) and the controller restrictions will be in charge of the lifetime enhancement of the PEMFC, avoiding starvation at the catalyst sites. The simulation model used as plant incorporates a water transport model that considers the macroscopic two-phase flow of water with mesoscopic pore filling effects in the cathode diffusion and catalyst layers to represent the voltage drop at each single cell [5]. Regarding the BoP, the cathode is air-fed with a compressor and the stack temperature is controlled with an active cooling system. This paper is organised as follows. In Section II the system description and simulation model based on distributed parameters are presented. The NDPO used to observe the internal states of the PEMFC is portrayed in Section III. In Section IV, the NMPC strategy is stated and developed based on the model presented in Section II. Simulation results for a given case study are discussed and analysed in Section V. Finally, Section VI recaps the conclusions of this paper and proposes some research lines for future work. II. SIMULATION MODEL A. System Description Figure 1 shows the plant scheme. It contains four subsystems: the PEMFC stack and load, the hydrogen delivery and recirculation auxiliaries, the air delivery and humidification auxiliaries and the cooling system. In this paper, it is considered that all the power is delivered by the fuel cell stack: no additional power sources or batteries are connected to the system. The hydrogen is stored in a high-pressure container and delivered to the anode through a pressure regulating valve. The cathode is air-fed with a compressor. Hydrogen is passively humidified through recirculation and the air is actively humidified at the humidity exchanger shown in Figure 1.    ʹ   ZĞĐŝƌĐƵůĂƚŝŽŶ ,LJĚƌŽŐĞŶ ƐƵƉƉůLJ     ŝƌƐƵƉƉůLJ ,ĞĂƚ ĞdžĐŚĂŶŐĞƌ ŵďŝĞŶƚ ƉƌĞƐƐƵƌĞ          ʹ ൉     ൉ WƵƌŐĞ ǀĂůǀĞ Fig. 1: PEMFC stack and BoP The PEMFC stack is an assembly of nfc single-channel cells identical as the one displayed in Figure 2a. Each cell has a channel length of 0.4 m, a channel width of 1 mm and a channel depth of 0.7 mm. The main fuel cell parameters are included in Table I. a) Cathode discretisation detail (y-direction) Membrane GDL Channel CCL    Vapour transport Liquid transport Condensation Evaporation Diffusion Anode Cathode Water production Electro-osmotic drag ȟ Ɂ  Ɂ  Ɂ  Ɂ  Ɂ  Ɂ  Ɂ  Anode gas channel Cathode gas channel +  +  2 1   2  +  2 PEMFC discretisation detail (z-direction) … … ȟ Membrane b) Fig. 2: (a) Single-channel PEMFC scheme and (b) detail of the discretisation volumes and water transport mechanisms B. Fuel cell stack model Each cell in the stack is modelled with a 1+1D or quasitwo dimensional parameters model. The gas flow transports along the z-direction are described with partial derivatives and the transports along the y-direction are considered lumped parameters [14]. The main model assumptions are the following: The gas species behave as ideal gases in all the simulation domain. The anode over-voltage is considered negligible compared to the cathode [15] and thus, activation and concentration losses are only considered at the cathode side. Correspondingly, liquid water formation is only considered in the cathode. Water mechanisms in the membrane include electro-osmotic drag and back-diffusion [15]. The CCL is a porous structure consisting of a number of pores (np) with fixed radius τp. At each pore there is a number of platinum particles (nP t) with a fixed radius τP t. The temperature of each single-cell (Tcell) is described according to the thermal model included in Section II-D. The load current is an input of the electrochemical model. 1) Electrochemical model: The fuel cell operating voltage Vfc is calculated as the total sum of the nfc single-cell TABLE I: Fuel cell parameters Parameter Description Units ciConcentration of i-th gas mol m−3 DiDiffusion coefficient of i-th gas m2s−1 ErIdeal potential voltage V FFaraday constant C mol−1 iCurrent density A m−2 KPressure drop coefficient m2s−1Pa−1 MMolar mass kg mol−1 ˙niy-direction flux of i-th gas mol m−2s−1 nfc Number of fuel cells in the stack — n(y,z)(y, z)-direction discretisation volumes — pPressure Pa RGas constant J mol−1K−1 SWater source term mol m−3s−1 TTemperature K VElectrical potential V vFlow velocity m s−1 δy-axis thickness m ∆G∗Gibbs activation energy J mol−1 ρVolumetric mass density kg m−3 τRadius m voltages Vfc = nfc X h=1 Vfc,cell,h,(1) being hthe index for each individual single-cell. Each singlecell voltage is modelled with the Butler-Volmer equation Vfc,cell,h =Er− RTcell,h α2Flog i i0−log pO2 pO2,ref −iRohm, (2) where Eris the ideal potential voltage of the fuel cell, αis the cathode charge transfer coefficient, Rohm is the internal resistance of the membrane, that depends on its water content value and pO2and pO2,ref are the oxygen pressure and oxygen pressure reference at the CCL respectively. The exchange current density at the cathode i0is a function of the fuel cell temperature Tcell,h, oxygen pressure at the CCL and the Electrochemically active Surface Area at the CCL (AECSA) [15, 16] i0=i0,ref AECSA Ageo pO2 pO2,ref 0.5 eh−∆G∗ RTcell,h 1−Tcell,h Tref i, (3) where i0,ref is the intrinsic catalytic Pt activity at normal conditions (Tref and pO2,ref ), Ageo is the total surface area of the electrode and ∆G∗is the Gibbs activation energy for the oxygen reduction reaction at the CCL. 2) Two-phase water model: As shown in Equation (2), the fuel cell voltage depends on the exchange current density i0which is greatly determined by the fuel cell temperature and the AECSA. In this paper, the AECSA in the CCL is modelled following a mesoscopic pore structure that considers only primary pores with a fixed pore size [5] AECSA =4πτ2 P tnpnP t,if 2τP t < τp1−3 √1−sCCL, 2πτP tτpnpnP t 1−3 √1−sCCL,else, (4) where npis the number of pores in the CCL volume with a fixed radius of τp= 10 nm and nP t= 1 is the number of Pt particles per pore, each one with a fixed radius of τP t= 2 nm. The equivalent number of pores npis computed considering the total CCL volume, the volume of a single pore and the specific porosity of the CCL. In the case where a pore-size or a particle-size distribution exists, Equation (4) is not valid and a new expression for AECSA would be required as proposed in [5]. The ratio of liquid volume to the total volume of void space in the porous structure of the CCL is defined by sCCL. Figure 3 shows an example of a pore in the CCL structure and the active area of a single Pt particle with different levels of water. At the right side of Figure 3, the layer of water covers the Pt particle and therefore, the active area is the total surface of the sphere, denoted by 4πτ2 P t. When the layer of water does not cover the Pt particle (2τP t > τp−τc), only the area in contact with the water is considered active. Following the geometrical relation between τcand sCCL proposed and experimentally validated in [5], the expression for the active surface in the second case of Equation (4) is obtained. ߬ ௉௧ ߬ ௉ ߬ ௖ ߬ ௉௧ ߬ ௉ ߬ ௖ ĐƚŝǀĞWƚĂƌĞĂ Fig. 3: Variation of the estimated Pt active area (depicted by the red line) in Equation (4) A two-phase (liquid and vapour) water model is implemented at the CCL and the GDLs. The partial differential equation that defines the ratio of liquid volume sin the CCL and cathode Gas Diffusion Layer (GDL) is expressed as follows [5]: ∂s ∂t = Sl H2O+Ds∂2s ∂y2+∂2s ∂z2 Ksorpερl H2O .(5) where Dsdenotes the liquid water diffusivity throughout a layer with specific porosity ε,ρl H2Ois the liquid water density and Ksorp is the time constant for the sorption of water into the porous layers. To compute the liquid water source term Sl H2Oin Equation (5), the water generation (Sgen H2O), water evaporation (Sevap H2O) and water transport through the membrane (SM H2O) terms are needed (see first row of Table II). The expressions for these terms can be found in [5]. Consequently, to obtain sCCL, the source terms for the CCL are the ones introduced in Equation (5). Then, sCCL is used to obtain the value of AECSA in Equation (4). TABLE II: Water source terms Anode Cathode GDL CL GDL CL Sl H2O0 0 −Sevap H2O−Sevap H2O+Sgen H2O+SM H2O Sg H2O0 0 Sevap H2OSevap H2O 3) Gas flow model: The gas species flow dynamics are described by mass balance equations along the PEMFC gas channels (see Figure 2b) [14]: ∂ci ∂t =∂ ∂z (vci)−˙ni δ+Sg i,(6a) v=−κ∂ ∂z ,(6b) p=RTcell,h Pici,(6c) being cithe concentration of the i-th gas, where subscript istands for the gaseous species, namely i=H2for the hygrogen, i=O2for the oxygen, i=N2for the nitrogen and i=H2Ofor the vapour water. The reaction and water molar transports from the Membrane Electrode Assembly (MEA) are modelled in ˙niand they are defined as lumped parameters perpendicular to the gas channels in the y-direction [14]. The y-direction thickness of the anode and cathode gas channels is represented by δ.Sg iis a function of the evaporation rate of the liquid water in the fuel cell when i=H2O(Sg i= 0 for any other gaseous species). In this paper the effect of liquid water is only considered in the cathode side of the PEMFC as depicted by the source terms in Table II. 4) Diffusion model: Hydrogen, oxygen and water diffuse from the gas channels through the GDLs and CLs because of a concentration gradient. The effect of nitrogen diffusion is not considered since it does not react. In this paper, diffusion follows Fick’s first law [15] and the concentrations at the CLs after diffusing are equal to: ∂ci,CL ∂t =Di ∂ci ∂y2,(7) being Dithe diffusion coefficient of the i-th gas species. C. Air supply system model In this paper, the cathode side of the fuel cell stack is fed with compressed air. A nonlinear compressor model is included in the simulation model. The nonlinear dynamics for the compressor, that include the dynamics of the oxygen (pO2), nitrogen (pN2) input pressures and the compressor angular velocity (ωcmp), are described in [17]. The input of the compressor model is the compressor motor current Icmp. The fuel cell current Ifc affects the oxygen partial pressure: the oxygen consumption in the PEMFC alters the cathode input manifold pressure [15]. The main parasitic loss in the system configuration described in Figure 1 is the power consumption of the compressor which can be expressed as Pcmp =τcmpωcmp,(8) being τcmp the torque of the compressor motor given a Icmp current. D. Thermal model 1) Stack temperature model: In Figure 1 the active cooling system model implemented is represented. It consists of a pump that circulates the coolant fluid (water) through a heat exchanger. The temperature of the coolant is reduced in the heat exchanger by means of forced air convection using an electric fan. The energy balance of the PEMFC enables the computation of the stack temperature Tfc as follows [18]: MCfc dTfc dt =˙ Etot −˙ Egross −Qcool −Qconv,(9) being MCfc the thermal mass of the fuel cell stack, Etot the total energy available, Egross the gross electrical energy supplied by the fuel cell, Qcool is the thermal energy dissipated by the coolant in the heat exchanger circuit and Qconv is the convective heat transfer to the environment. To complete Equation (9), the following set of equations is employed: ˙ Etot =nIfc 2F∆H, (10a) ˙ Egross =VfcIfc,(10b) Qcool = ˙mcoolCp,cool∆Tcool,(10c) Qconv =kconvAstack (Tfc −Tamb),(10d) where ˙mcool and Cp,cool are the mass flow and heat capacity of the coolant fluid respectively. The enthalpy variation during the reaction is ∆H,kconv is the thermal conductivity coefficient for the stack assembly and Astack is its surface area. The coolant temperature drop ∆Tcool is a design parameter of the model that is assumed to be constant. 2) Single-cell temperature model: The thermal model presented in the previous section represents the temperature value of the stack (Tfc). Internally, fuel cells have temperature gradients as a result of the internal chemical reactions and heat generation at the CLs. In [19] different thermal gradients are shown for different simulation cases. After achieving steady state, the thermal variation between the end plates and the middle point of the stack appears to be approximately 10 K. In this paper a probabilistic approach is applied to Tfc in order to represent the temperature gradient at the stack without increasing the complexity of the model. Considering the normal probability distribution f(x|µ, σ2), where µis the mean of the distribution and σis the standard deviation, the fuel cell temperature Tfc from Equation (9) is then used to compute each cell temperature as follows: Tcell,h =Tfc +kdistf(x|µ, σ2),(11) being kdist the constant to model the total temperature gradient inside the stack. For the simulations in Section V, kdist = 50 (which gives a ∆T= 10 Kbetween the middle and border cells, as proposed by [19]) is selected. E. Finite-difference discretisation The yand xspatial derivatives in the simulation model are discretised to numerically solve them and to take advantage of the boundary conditions of the problem (e.g. input molar fluxes and the external ambient pressure at the end of the gas channels). The details of the discretisation procedure were studied in [9]. III. OBSERVER MODEL In this Section, an improved version of the High-Order Sliding-Mode (HOSM) observer presented in [9] is developed to estimate the PEMFC full gas concentrations profiles. The main improvement is the extension of the estimation procedure to the GDLs and CLs, which allows the dynamic and robust estimation of relevant variables for the control strategy: AECSA, water content (Λ), Rohm and ci,CCL. This observer model only includes the gas species flow dynamics in Equation (6) and the diffusion model in Equation (7). A. HOSM state observer 1) Block controllable structure: The first step to implement the HOSM observer for the gas concentrations profiles (ˆci,j) is to express Equation (6) in block-controllable form [20] for nzfinite discretisation volumes along the gas channel. The complete mathematical procedure was studied in [9]. 2) HOSM back-stepping algorithm: To estimate the full concentrations profiles, a back-stepping observation algorithm is employed. The observation procedure guarantees that the estimation is achieved even in the presence of model uncertainties. It follows an r-step algorithm where r∈[1,2, . . . , nz], as presented in [9]. 3) Gas diffusion estimation: Once the concentrations estimation along the gas channels is achieved, the diffusion through the GDLs and CLs in the y-direction is carried out using the diffusion law in Equation (7). The estimated gas concentrations at the GDLs and CLs can be expressed as follows: ˙ ˆci,(j,n)=     Di ∆y2ˆci,j −2ˆci,(j,n)+ ˆci,(j,n+1),if n= 1, Di ∆y2ˆci,(j,n−1) −ˆci,(j,n),if n=ny, Di ∆y2ˆci,(j,n−1) −2ˆci,(j,n)+ ˆci,(j,n+1),else, (12) where j,[1, . . . nz]is the discretisation volume subscript for the z-direction and n,[1, . . . ny]is the discretisation volume subscript for the y-direction. Note that the term ˆci,j when n= 1 in Equation (12) is the observed variable at the gas channel. B. Estimation of the AECSA Once the full gas concentrations profile ˆci,j is obtained, AECSA is computed using the estimates and the PEMFC model presented in Section II. The discretised variables and their associated dependencies needed to obtain AECSA are the following: ˆ Λj=f(ˆcH2O,j),(13) ˆ Rohm,j =f(ˆ Λj),(14) ˆci,(j,CCL)=f(Equation (12),ˆci,j).(15) The amount of active surface at each discretised volume of the CCL, AECSA,j is isolated from Equation (3) AECSA,j = i0,j Ageo (i0,ref )−1pO2,j pO2,ref −0.5 e−∆G∗1−Tcell,h Tref R−1Tcell,h−1.(16) And assuming that the temperature of the individual cells can be measured, the only unknown in Equation (16) is the current exchange at each discretisation volume i0,j. IV. NMPC STRATEGY The fuel cell stack in Figure 1 delivers the demanded power. The control strategy has to guarantee that this is done under the proper operating conditions to prevent the accelerated degradation of the system. In this section the proposed control strategy to achieve this objective is described. A. Control objectives As mentioned beforehand, degradation of the PEMFC derives into a reduction of the AECSA [5] and therefore, less available area for the chemical reaction to take place (lower PEMFC efficiency). In this paper, one of the objectives of the control strategy is to use the available control inputs of the system to maximise the AECSA of the stack. Moreover, it is critical to maintain suitable amounts of fuel and oxidant at the catalyst sites to avoid starvation, which causes permanent damage in PEMFCs [4]. The control strategy has to maintain a safe amount of hydrogen and oxygen concentrations along the catalyst sites during the operation to guarantee that the degradation rate of the fuel cell is not accelerated. Finally, the amount of water in the system has to remain between certain boundaries guarantee that the membrane is humidified without flooding the fuel cell, in order to avoid the acceleration of degradation mechanisms [4]. B. Prediction model The NMPC uses Equations (2), (6), (7) and (16) as prediction model over a discrete-time variable k∈R. Moreover, the prediction model does not include the full complexity of the simulation model. Actuating over the fuel cell temperature Tfc and the oxygen pressure pO2, it is possible to find an optimal AECSA over the prediction horizon and guarantee the first control objective defined in Section IV-A. Note that Equation (16) refers to each jdiscretisation volume along the y-direction. The optimisation will be done for the total sum of the active surface. C. Input and state constraints The input constraints are fixed by the physical characteristics of the equipment employed. They are set as follows: 1.2≤uIcmp ≤3.5A,(17a) 0≤uRHC ref ≤1,(17b) 0≤u˙mcool ≤7kg s−1,(17c) 0≤uH2,in ≤250 mol m−2s−1.(17d) Regarding the state constraints, the lower bounds of the hydrogen and oxygen concentrations are set to be higher than zero at the catalyst sites, guaranteeing that the system does not operate under starvation conditions. The higher bounds are fixed taking into account nominal values from models reported in the literature [14]. Moreover the water concentrations constraints are set in order to guarantee the correct humidification of the membrane without flooding the system. The state constraints are summarised by the following set of equations: 0< cH2,j,CCL ≤70 mol m3,∀j, (18a) 0< cO2,j,CCL ≤25 mol m3,∀j, (18b) 0< cH2O,j,CCL ≤10 mol m3,∀j. (18c) D. Cost function According to the control objectives (Section IV-A), the performance indicator for the AECSA is defined for all the nzdiscretisation volumes at instant kas `e k,−||AECSA||2.(19) Notation || · ||2indicates the 2-norm (Euclidean norm) [21] of AECSA = [AECSA,1, . . . , AECSA,j]. Because of the minimisation problem described in Section IV-E, `e khas a negative sign to maximise the AECSA during the optimisation procedure. Regarding the control objective to operate under smooth control actions, it can be expressed at each time instant kas follows: `∆u k=||∆u||Wu=p∆uTWu∆u(20) being ∆u(k),u(k)−u(k−1) the slew-rate of the control signals in Equation (17). The slew-rate terms are weighted with a diagonal weight matrix Wuof suitable dimensions: Wu=diag(Wu1;Wu2;Wu3;Wu4). Given Equations (19) and (20), the resultant control function that has to be minimised is the following one: Jk=λ1`e k+λ2`∆u k,(21) being λ1and λ2weights to prioritise between control objectives `e kand `∆u k. E. Optimisation problem The cost function stated in Section IV-D is minimised using the receding horizon principle for NMPC [12], solving an optimisation problem at each step of the prediction horizon. The result is an optimal input sequence that aims to minimise Equation (21) at each sampling time. Let u(k),(u(0|k), . . . , u(Hp−1|k)) (22) be the sequence of control inputs over a fixed-time prediction horizon Hp(Hp≥2), depending also on the initial condition x(0|k),x0. The NMPC Finite-Time Open-loop Optimization Problem (FTOOP) is formulated as Problem 1 (NMPC FTOOP): min u(k)∈Rm×HpJ(x0,u(k)),(23) subject to •predicted states from the HOSM observer at time k, •system model in (2), (6), (7) and (16) over Hp, •input constraints in (17) over Hp, •state constraints in (18) over Hp, where J(·) : Um×Hp×RHp7→ Rin (21) is the cost function, with m= 4 denoting the number of control inputs. Assuming that the FTOOP (23) is feasible, there will be an optimal solution for the sequence of control inputs u∗(k),(u∗(0|k), u∗(1|k), . . . , u∗(Hp−1|k)) (24) and then, according to the receding horizon philosophy, u∗ i(0|k)is applied to the system, while the process is repeated for the next time instant k∈Z. V. SIMULATION RESULTS The control strategy has been tested by simulation using the driving cycle described in Section V-A. The mesh for the simulation model consists of 5 elements equally distributed along the z-direction, 5 elements for the GDL and 5 elements for the CCL along the y-direction. The mesh of the observation model for the full concentrations profile consists of 5 elements along the z-direction and the initial state for the i-th gas concentrations values are ˆci(t= 0) = 0∈R5×nz. Simulations have been carried out using MATLAB R2011a (32 bits), running in a PC Intel Core i7-3770 at 3.40 GHz with 8 GB of RAM. A. Case study A synthetic driving cycle is used to test the control strategy in a simulation framework. It is based on the New European Driving Cycle (NEDC). The NEDC is a speed profile that represents urban and highway scenarios to evaluate pollutant emissions and energy management strategies for different types of engines (i.e. gasoline, electric, etc.). In [22] the NEDC speed profile was converted into a demanded current density profile (see Figure 5) for a fuel cell powered car. This is going to be the case study profile used to analyse the performance of the control strategy. Figure 4 shows the driving cycle used to test the proposed controller. The fuel cell has to provide the total power demanded by the driving cycle, the compressor and the secondary auxiliaries. This is denoted by the following: Pfc,elec =PDC +Pcmp +Paux,(25) 0 500 1000 1500 2000 2500 3000 0 0.25 0.5 0.75 1 1.25 1.5 Current density [A cm−2] Time [s] 0 500 1000 1500 2000 2500 3000 0 50 100 150 200 250 300 350 Power [W] Demanded NEDC current density Auxiliaries power demand Fig. 4: Synthetic NEDC profile and experimental determination of the auxiliaries power demand where PDC is the driving cycle power demand. Moreover, the power losses of the secondary auxiliary BoP subsystems Paux have been experimentally characterised in a test bench [23] as presented in Figure 4. For comparison purposes, in the following sections the NMPC controller is compared with a classic constant stoichiometry controller. The stoichiometries for CCS are 1.3 at the anode and 2.0 at the cathodic side. Table III shows CNMP C weights, sampling time and simulation times for the case study. Moreover, the stack is composed of nfc = 6 identical singlechannel PEMFCs with a total surface area Ageo = 25 cm2. TABLE III: CNMP C configuration Parameter Description Value λiObjective prioritisation weights [1,1] Wuiuislew-rate weight [1,1,1,1] HpPrediction horizon 2 ∆tSampling time 1 s B. Results and discussion The behaviour of the modelled AECSA and the observed variable \ AECSA for the controllers CNMP C and CCS is presented in Figure 5. At the beginning of the simulation, the low current demand produces an accumulation of condensed liquid water due to a lower temperature and therefore, an increase of AECSA in both cases. After the current demand is increased, the liquid water is being evaporated from the CCL, which is represented by the decreasing values of AECSA. The controller CNMP C , through the use of the cooling circuit reduces Tfc until the AECSA starts increasing, maximising its value even in the presence of sudden current demand variations. The active area using CCS is clearly lower than in the case of CNMP C during all the simulation time. As expected, the dynamic response of AECSA is slow due to its dependency to slow dynamic effects such as the temperature and evaporation of water. 0 500 1000 1500 2000 2500 3000 0 5 10 15 20 25 AECSA [cm2] Time [s] 0 500 1000 1500 2000 2500 3000 0.3 0.6 0.9 1.2 1.5 Current density [A cm−2] AECSA for CNMP C b AECSA for CNMP C AECSA for CCS b AECSA for CCS Current density Fig. 5: Evolution of modelled and observed AECSA for CNMP C and CCS and NEDC stack current density demand 0 500 1000 1500 2000 2500 3000 47.5 50 52.5 55 57.5 60 Efficiency [%] Time [s] 0 500 1000 1500 2000 2500 3000 0 0.25 0.5 0.75 1 1.25 Current density [A cm−2] CNMP C : max(ECSA) CCS :λa= 1.3; λc= 2.0 Demanded NEDC current density Fig. 6: Fuel cell efficiency for CNMP C and CCS Regarding the observation of AECSA, it is shown in Figure 5 that the observer tracks the real value properly throughout the NEDC cycle. This is an important contribution of this research, since as mentioned in Section I, AECSA can not be physically measured while the PEMFC is operating. Regarding the controllers efficiency, the fuel cell efficiency ηfc is higher when using controller CNMP C as presented in Figure 6. This is because higher values of AECSA produce an increase of Vfc,cell in Equation (2), and therefore, better fuel cell efficiency [15]. Moreover, as shown in Table IV, the proposed strategy uses a lower quantity of injected hydrogen for the same driving cycle, contributing to the increase in the global efficiency of the system. However, CNMP C makes use of higher RH values and therefore, more quantity of injected water as presented in Table IV. The dynamic behaviour of the manipulable inputs applied to the system using CNMP C are shown in Figure 7. These optimal control inputs maximise the AECSA while maintaining suitable operating conditions for the PEMFC, such as avoiding starvation situations, which would accelerate the degradation of the fuel cell. 0 500 1000 1500 2000 2500 3000 1 1.2 1.4 1.6 1.8 Compressor current [A] a) 0 500 1000 1500 2000 2500 3000 0 2 4 6 8 mCool [kg s−1] b) 0 500 1000 1500 2000 2500 3000 50 51 52 53 54 55 Cathode RH [%] c) 0 500 1000 1500 2000 2500 3000 0 50 100 150 200 H2 Input [mol m−2 s−1] Time [s] d) Fig. 7: Control actions supplied by the NMPC The HOSM observer feeds the NMPC controller with the estimated state vector for the gas concentrations. While in Figure 8 only the observation of the concentrations in the middle point of the gas channels is presented, the observation is performed in all of the discretisation volumes. This is done to facilitate the reading of the results. Figures 8a and b, refer to the concentration estimation at the anode side. On the other hand, Figure 8c and d, refer to the cathode gas channel concentrations estimation. In the case of Figures 8a and c, denoting the hydrogen and oxygen concentrations, the controller maintains these values between certain boundaries with the objective of avoiding local starvation. For the water concentrations in the anode and cathode sides of the PEMFC, depicted in Figures 8b and d, the controller guarantees that the humidification of the fuel cell is adequate without flooding the system. As pointed out previously, the implementation of the NMPC using a nonlinear distributed parameters prediction model 0 500 1000 1500 2000 2500 3000 30 40 50 60 70 80 Concentrations [mol m−3] cH2 ˆcH2,NDP O a) 0 500 1000 1500 2000 2500 3000 0 1 2 3 4 Concentrations [mol m−3] cH2O ˆcH2O,NDP O b) 0 500 1000 1500 2000 2500 3000 7 7.2 7.4 7.6 7.8 8 Concentrations [mol m−3] cO2 ˆcO2,NDP O c) 0 500 1000 1500 2000 2500 3000 0 0.5 1 1.5 2 Time [s] Concentrations [mol m−3] cH2O ˆcH2O,NDP O d) Fig. 8: Behaviour of the NDPO versus the plant states in the middle discretisation volume of the anode (a and b) and the cathode (c and d) introduces a high computational effort. Nevertheless, the total accumulated computation time remains below the total simulation time as shown in Figure 9, making it feasible for implementation in future revisions of this work. TABLE IV: Results for CNMP C and CCS ηfc [%] ¯ AECSA [cm2]. H2,inj [gr] H2Oinj [gr] CNMP C 57.56 11.57 777.66 383.82 CCS 55.27 5.30 814.50 280.39 VI. CONCLUSION In this paper, an NMPC strategy has been proposed to guarantee the maximum active area in the CCL while avoiding fuel and oxidant starvation at both sides of the fuel cell. The 0 500 1000 1500 2000 2500 3000 0 0.5 1 1.5 2 2.5 Time [s] Iteration computation time [s] 0 500 1000 1500 2000 2500 3000 0 100 200 300 400 500 Accumulated computation time [s] Fig. 9: NMPC computation time performance of the NMPC has been evaluated, obtaining satisfactory results considering an automotive driving cycle in a simulation scenario. A NDPO estimates the unmeasured states and injects this information into the controller. The simulation model in Section II includes complex dynamics. However, the observation and prediction models are a reduced form of the simulation model, simplifying future real implementation of the controller. The dynamic behaviour of the AECSA is orders of magnitude slower than the chemical reactions that take place in the PEMFC. As shown in Figure 5, the control strategy guarantees that the AECSA improves during the simulation time when compared to other control strategies that do not consider the maximisation of the active surface. Meanwhile, the control optimiser avoids starvation scenarios that could harm the fuel cell and reduce its lifespan. Therefore, the combination of the control objectives with the system constraints provides an enhancement of the lifetime of the PEMFC, mitigating the degradation mechanisms that naturally occur in these systems when they operate. In this work, the analytical development of a novel NMPC strategy with nonlinear observation in a PEMFC-based system has been studied. A forthcoming study regarding the experimental validation of the solution in a real PEMFC-based system is in progress. REFERENCES [1] O. Z. Sharaf and M. F. Orhan, “An overview of fuel cell technology: Fundamentals and applications,” Renewable and Sustainable Energy Reviews, vol. 32, pp. 810–853, 2014. [2] A. de Frank Bruijn and G. J. Janssen, “PEM fuel cell materials: Costs, performance and durability,” in Fuel Cells, pp. 249–303. Springer, 2013. [3] T. Jahnke et al., “Performance and degradation of proton exchange membrane fuel cells: State of the art in modeling from atomistic to system scale,” Journal of Power Sources, vol. 304, pp. 207–233, 2016.