Modeling, state estimation and nonlinear model predictive control of cathode exhaust gas mass flow for PEM fuel cells

Control Engineering Practice 49 (2016) 76–86

Contents lists available at ScienceDirect

Control Engineering Practice journal homepage: www.elsevier.com/locate/conengprac

Modeling, state estimation and nonlinear model predictive control of cathode exhaust gas mass flow for PEM fuel cells M. Schultze n, J. Horn Institute of Control Engineering, Department of Electrical Engineering, Helmut-Schmidt-University/University of the Federal Armed Forces Hamburg, Holstenhofweg 85, 22043 Hamburg, Germany

art ic l e i nf o

a b s t r a c t

Article history: Received 19 June 2015 Received in revised form 21 January 2016 Accepted 22 January 2016 Available online 3 February 2016

Polymer electrolyte membrane fuel cells are efficient energy converters and provide electrical energy, water and oxygen depleted air with a low oxygen content as exhaust gas if fed with air. Due to their low emission of greenhouse gases and noise they are investigated as replacement for auxiliary power units currently used for electrical power supply on aircraft. Oxygen depleted air, called ODA-gas, with an oxygen concentration of 10–11% and a low humidity can be used for tank-inerting on aircraft. A challenging task is controlling the fuel cell system for generation of dehumidified ODA-gas mass flow while simultaneously keeping bounds and gradients on control inputs. This task is attacked by a nonlinear model predictive control. Not all system states can be measured and some states measured exhibit a significant time delay. A nonlinear state estimation strategy builds the entire system state and compensates for the delay. The nonlinear model predictive control and the state estimation are derived from the system model, which is presented. Simulation and experimental results are shown. & 2016 Elsevier Ltd. All rights reserved.

Keywords: State estimation Nonlinear model predictive control PEM fuel cell

1. Introduction The use of electrically powered systems has been growing in civil aircraft during the last decades (McLoughlin, 2009). On aircraft electrical power supply during ground operation is done by an auxiliary power unit (APU) that significantly emits CO2, NOx and noise. Polymer electrolyte membrane (PEM) fuel cells are very efficient energy converters, suitable for dynamic applications and offer the potential of drastically reducing these emissions. Therefore, PEM fuel cells are studied as replacement of an aircraft's APU. In addition to supply of electrical energy, PEM fuel cells deliver water and oxygen depleted exhaust gas if provided with air. They are suitable for supply of electrical power, of water and of oxygen depleted air (ODA-gas) for tank-inerting on aircraft during ground operation as well as during flight. This multifunctional use on aircraft offers a great benefit as size of current systems for power and water supply as well as tank-inerting can be drastically reduced, which offers great savings on total aircraft weight. The multifunctional fuel cell system (MFFCS) and its functions on aircraft are shown in Fig. 1. For tank-inerting the MFFCS must provide ODA-gas with an oxygen content of 10–11% (Aviation Rulemaking Advisory Committee, 1998; Kallo et al., 2010) to prevent ignition of n

Corresponding author. E-mail addresses: [email protected] (M. Schultze), [email protected] (J. Horn). URL: http://www.hsu-hh.de (J. Horn). http://dx.doi.org/10.1016/j.conengprac.2016.01.006 0967-0661/& 2016 Elsevier Ltd. All rights reserved.

fuel vapors and a low humidity to prevent icing and contamination inside the tanks (Tomlinson, Barker, Venn, Hickson & Lam, 2011). PEM fuel cell systems have been studied for electrical power supply of autonomous robots (Niemeyer, 2009) or for automotive applications (Karnik, Sun, Stefanopoulou & Buckland, 2009; Pukrushpan, Stefanopoulou & Peng, 2004). Operation of PEM fuel cells for inerting is a novel research topic and has not yet been studied in detail.

Fig. 1. Multifunctional fuel cell system (MFFCS) for electrical power supply, water supply and tank-inerting on a civil aircraft.

Proper fuel cell system operation such as keeping the membrane well hydrated and to proper supply fuel and air as oxygen carrier is a central aspect (Borup et al., 2007; Pukrushpan et al., 2004). In the

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

multifunctional context, controlling the fuel cell system for ODA-gas mass flow further requires satisfaction of limitations on stack current, its gradient and exhaust gas dehumidifier cooling temperature. Due to its inherent capability of incorporating constraints model predictive control is an advantageous approach as compared to classical control techniques. Model predictive control (MPC) has already been successfully applied to fuel cell systems for a timeoptimal warm-up (Müller, Stefanopoulou & Guzzella, 2007), an optimal load-sharing between fuel cell and short-term energy storage (Vahidi, Stefanopoulou & Peng, 2004, 2006) or for electrical power supply (Niemeyer, 2009). Nonlinear Model Predictive Control (NMPC) of a MFFCS for exhaust gas is a novel topic. NMPC requires the current system state, which is gained by a nonlinear state estimation strategy based on the Sigma-Point-Kalman-Filter. NMPC strategies designed for the entire system (Schultze & Horn, 2013a) or assuming an underlying state estimation (Schultze, Hähnel & Horn, 2014) are part of previous work. Previous works have not presented the entire process of developing a NMPC for a MFFCS. This paper presents all steps in synthesis of the NMPC for exhaust gas mass flow, state estimation strategy and modeling of a MFFCS. The entire system model is presented in Section 2. This entire model is used for the design of the state estimation strategy presented in Section 3 and was used in designing the underlying fuel cell stack cooling temperature controller. The state estimation strategy supplies the NMPC with the entire state vector. Exploiting the underlying temperature control, the NMPC presented in Section 4 is based on a subset of the entire system model to save computational time. Simulation and experimental results are shown. The NMPC is compared to basic controls to show the advantages of this novel control strategy.

77

The system model described in the following is derived from the nonlinear model (Schultze & Horn, 2013b) and is part of previous research (Schultze & Horn, 2013c). The model consists of ordinary differential equations, has 11 states and covers the MFFCS shown in Fig. 2. Air is modeled consisting of 21% (vol.) oxygen (O2) and 79% nitrogen (N2). Electrochemical as well as electric processes are very fast as compared to gas-transport and thermal processes and are modeled as static. The underlying stack cooling temperature control as well as the system's high thermal mass cause the thermal process to have only little influence on the gasdynamic process. Gas dynamics are therefore considered being decoupled from thermal dynamics. Fig. 3 shows the entire system model structure. Due to a long pipe exhaust gas transport from the stack to the mass flow sensor underlies a significant transport delay. In the following, the delay is not stated explicitly. For the model, delayed measurements are derived by delaying the corresponding values by Td . 2.1. Gas dynamics The gas dynamics model describes the behavior of ODA-gas mass flow and its oxygen content. The vector of inputs ug = [Ist, λO2 ]T to the gas dynamics model comprises stack current Ist and stoichiometry λO2 . Feed air Wmfc is provided by a controlled valve (MFC), which is modeled as a

2. Multifunctional fuel cell system model The flow of electrical energy from the MFFCS, additional short time energy storages and the main aircraft generators to the electrical consumers would be managed by an electrical power management strategy (PM). The PM is not scope of this work. In the following it is assumed that the PM applies a fuel cell stack current as requested. Moreover, the PM is assumed to react ideally and immediately. As a simplification, the electric load in Fig. 2 represents the PM and the aircraft's electric system. Fig. 2 shows the experimental setup and a schematic of the multifunctional fuel cell system. The MFFCS fuel cell stack has an anode recirculation loop for efficient hydrogen use. The polymer membrane inside the stack separates the cathode and anode volumes. The membrane is electrically not conductible, however, allows for water and proton transport. Its thinness facilitates high water diffusivity and hence promotes a stack-internal membrane humidification. Stack cooling temperature is limited by manufacturer specifications, so that the stack does not dry out or experiences flooding. Supply of air and fuel, especially during transients, must be kept to minimize the risk of oxygen and fuel starvation. Limitation of stack current rate decreases these risks as well as the risk of dynamic voltage losses and therefore protects against stack damage. The stack cooling system consists of two cooling loops that are coupled by an intercooler. Cooler and cooling valve are in the external loop. Exhaust gas dehumidification is performed by water separators and a condenser that is connected to a temperature-controlled cooling unit. The fuel cell stack is electrically connected to a controllable ohmic electric load. System inputs are stack current Ist , stoichiometry λO2 , condenser cooling reference temperature Tcond,ref and cooling valve position uv . ODA-gas mass flow and O2 content are measured by a hot-wire anemometer and a lambda probe. Relative humidity, temperature and pressure measurements determine ODA-gas humidity.

Fig. 2. Experimental setup and schematic of the multifunctional fuel cell system (MFFCS) with fuel cell stack, controllable electric load, stack cooling system and exhaust gas dehumidification.

Fig. 3. Main structure of the multifunctional fuel cell system (MFFCS) model: gas dynamics and thermal dynamics.

78

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

PT1-element with time constant τmfc . Dynamics of the ODA-gas mass flow Woda through the stack cathode volume and the downstream pipes is modeled as a PT1-element with time constant τfc . Due to its low humidity after dehumidification ODA-gas is assumed a dry gas. The disturbance mass flow z at the cathode inlet covers tolerances in feed air supply and is modeled as a constant additive mass flow. Measurement of the ODA-gas oxygen content is modeled as a PT1-element with time constant τ O2:

τmfc Ẇ mfc = Wmfc,ref − Wmfc

τfc Ẇ oda = Wmfc + z − Ist

(1)

ncells MO2 − Woda 4F

Fig. 4. Stationary ODA-gas oxygen content vs. stoichiometry.

(2)

τ O2 cȮ 2 = cO2,ca − cO2

(3)

ż = 0

(4)

MO2 and MN2 are the molar masses of O2 and N2, F is the Faraday constant and ncells is the number of cells in the fuel cell stack. Feed air is assumed to consist of oxygen and nitrogen with the mass fractions χO2 = 0.23 of oxygen and χN2 = 0.77 of nitrogen. The chemical reaction taking place in the fuel cell stack consumes an oxygen mass flow that is proportional to the stack current. An immediate chemical reaction being coupled to the stack current is assumed (Pukrushpan et al., 2004). Hence, the cathode exhaust gas O2 content cO2,ca is gained by the stationary molar outlet flows of oxygen NO2,ca and nitrogen NN2,ca . Generation of ODA-gas is assumed to take place during a continuous chemical reaction inside the fuel cell:

Wmfc,ref = Ist λ O2

ncells MO2 + 4F

(

0.79 M 0.21 N2

)

cO2,ca = NO2,ca/( NO2,ca + NN2,ca )

NO2,ca = χO2 /MO2

(

)( W

+ z) −

(

)( W

+ z)

NN2,ca = χN2 /MN2

mfc

mfc

(5)

(6)

ncells Ist 4F

(7)

(8)

The gas dynamics model with the state vector x g and the output vector y g (without delay) is as follows:

xġ = f g ( xg , ug )

T xg = ⎡⎣ Wmfc Woda cO2 z⎤⎦

(9)

(10)

Fig. 5. Simulation results of the MFFCS gas dynamics model (ODA-gas mass flow delayed by Td ) compared with experimental results for a change of stack current from a low to a high (A) and a high to a low value (B) at a current rate ΔImax = 5, 10, 20 A/s.

simulation and experimental results for stack current rates of 5, 10 and 20 A/s. The overlaid plots show that the model describes the real dynamics very well. The experimental results show the significant time delay Td in ODA-gas mass flow measurement with Td = 4 s. 2.2. Thermal dynamics

T yg = ⎡⎣ Woda cO2 ⎤⎦

(11)

Analysis of the stationary model 0 = f g ( x g , ug ) assuming z¼ 0 reveals that ODA-gas oxygen content depends on stoichiometry λO2 as shown in Fig. 4. The fuel cell system used in this study can be operated at λO2 as shown in Fig. 4. When setting a certain λO2 the air mass flow provided to the stack only depends on stack current as follows from Wmfc,ref in Eq. (5). Parameter τfc has been identified by a least square error minimization using transient experimental data. τmfc and τ O2 stem from component documentation. Fig. 5 shows a comparison of

The model of the thermal dynamics describes the system pressure, fuel cell stack voltage and thermal behavior, the stack as well as condenser cooling system and water loading of dehumidified ODA-gas. The input vector uth = [Tcond,ref uv ]T comprises the condenser cooling system reference temperature Tcond,ref and the stack cooling system cooling valve position uv . 2.2.1. Pressure model Pressure establishes very fast as compared to the remaining system dynamics and therefore is described by a static model (Niemeyer, 2009; Pukrushpan et al., 2004). The values pca , p2, p3, p4

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

and p5 describe the pressures in the cathode volume and the downstream pipes 2–5 as depicted in Fig. 2. The pressures establish according to the dry ODA-gas mass flow Woda gained by the gas dynamics model. Pressure plays an important role in determining fuel cell heat flows, stack voltage and ODA-gas water loading. Water loading X of a humid gas is defined in the following as the fraction X = mv /mg = (mass of vapor)/ (mass of dry gas). The nonlinear mass flow equations for the pressure are given in the following. An iterative algorithm solves the set of equations. To save computational time the pressure equations are solved starting at pipe 5. The separators remove liquid water from the gas flow. Due to the high separation rate their efficiency is assumed 100%. Exhaust gas flow in pipe 5 is turbulent and has a low water loading due to dehumidification. Therefore, it is modeled as a dry gas with flow constant k5u and ambient pressure pamb as downstream pressure. Inversion of Eq. (12) leads to pressure p5:

Woda = k5u p5 − pamb

(12)

To solve for pressures p3 and p2 Eqs. (16) and (17) are rearranged leading to Eq. (20): 2

p2 = p3 + ( k34/k23 ) (p3 − p4 )

(

p5 = ( Woda/k5u ) + pamb .

(13)

Pressure p4 is gained in analogy to Eq. (13) with flow constant k 45 and downstream pressure p5. Water loading after the condenser is low such that it is neglected for the flow equation for pressure p4 as follows: 2

p4 = ( Woda/k 45 ) + p5

(14)

In pipe 3 exhaust gas has a temperature of T3 = T2 and is fully saturated with vapor. Compared to pipes 3–5 that are assumed being perfectly insulated, pipe 2 shows a significant heat loss that is modeled by a constant temperature drop ΔT2 leading to T2 = Tst − ΔT2. Exhaust gas water loading X3 = X2 in pipe 3 is modeled as water loading X2 in pipe 2 as no significant evaporation or condensation is assumed to take place in the separator connecting pipes 2 and 3. Water loading X2 is gained as follows with the gas constant of ODA-gas Roda and of vapor RV . Saturated vapor pressure depends on temperature T in °C and is given as

(

pVsat (T ) = exp 17.2799 −

4102.99 T + 237.431

)

0.611657·103 (Baehr & Kabelac,

2012):

X2 =

( ) Roda p2 − pVsat ( T2 ) RV pVsat T2

(15)

The exhaust gas mass flows for pipes 2 and 3 are modeled as turbulent humid flows with water loading X2. Flow through the condenser connecting pipes 3 and 4 is gained by Eq. (16) with downstream pressure p4 and flow constant k34 . Flow through the separator connection pipes 2 and 3 is gained by Eq. (17) with downstream pressure p3 and flow constant k23:

(1 + X2 ) Woda = k34 p3 − p4

(16)

(1 + X2 ) Woda = k23 p2 − p3

(17)

At the cathode volume the gas flow is fully saturated. The humid gas mass flow in Eq. (19) depends on cathode water loading Xca , flow constant kca and downstream pressure p2:

X ca =

( ) Roda ( ) RV

pVsat Tst

pca − pVsat Tst

(1 + X ca ) Woda = k ca ( pca − p2 )

(20)

Taking the normal form of Eq. (16) and inserting Eqs. (15) and (20) leads to an equation with the unknown p3. This is done in analogy for Eq. (19) with Eq. (18) being inserted to lead to pca . The normal forms are solved by a Regula-Falsi-Algorithm of the Illinois-method (Ford, 1995). The algorithm starts with an initial interval such that the function values at the interval bounds have different signs. The interval is reduced iteratively. The Illinoismethod allows for a fast convergence. Due to the fuel cell stack's operating pressure and temperature liquid water resides in the cathode. Therefore, a fully saturated cathode gas is assumed. Based on the ideal gas law, cathode pressure pca is the sum of saturated vapor pressure pVsat (Tst ) and dry gas pressure. O2 partial pressure pO2,ca is gained as the fraction cO2,ca of dry gas pressure as follows:

( ))

pO2,ca = cO2,ca pca − pVsat Tst 2

79

(21)

2.2.2. Stack voltage Identical cells in the fuel cell stack are assumed. Stack voltage Ust is the sum of all ncells cell voltages and is gained by reversible cell voltage Urev (O'Hayre, Cha, Colella & Prinz, 2009), activation voltage losses Uact (Amphlett, Baumert, Mann & Peppley, 1995) and ohmic voltage losses Uohm as Ust = ncells (Urev − Uact − Uohm ). The anode pressure control leads to the assumption of a constant anode hydrogen pressure pH2. Due to its thinness and the limited stack current gradient the polymer membrane is assumed to be perfectly hydrated. Therefore, the membrane's ohmic resistance, which is influenced by membrane humidity, is modeled as a constant resistance Rohm . By parameters ζ1, ζ2, … , ζ4 and Rohm the voltage model is adapted to the fuel cell stack:

Urev = 1.229 − 0.85· 10−3 ( Tst − 298.15) ⎛ pH pO ,ca ⎞ 1 2 + 4.3·10−5 ⎜ ln + ln 2 ⎟ 2 p0 p0 ⎠ ⎝

(22)

⎛ pO2,ca 498/ T ⎞ st ⎟ + ζ T ln ( I ) Uact = ζ1 + ζ2 Tst + ζ3 Tst ln ⎜ e 4 st st ⎝ 5.08· 106 ⎠

(23)

Uohm = R ohm Ist

(24)

2.2.3. Stack heat flow and stack cooling mass flow The heat flow Q̇ st emit by the fuel cell stack is determined by a stationary energy balance across the stack involving the chemical reaction (thermoneutral voltage UHHV = ncells 1.48 V ) and fluid flows that are fed to the stack and leave the stack. Temperature of feed air Tair is constant. Water is generated during the chemical reaction and is transported outside the stack as vapor WV and liquid WL ( MH2O is the molar mass of water). Evaporation is assumed to occur instantaneously. Vapor mass flow is determined as the lower value of the maximum vapor mass flow Xca Woda that can be transported by the cathode exhaust gas and mass flow of reaction water. The constants coda , cair , cL , cV and h0 are the specific heat capacities of ODA-gas, air, liquid water, vapor and enthalpy of evaporation of water:

WV = min (Ist ncells MH2O/(2F ), X ca Woda )

(25)

WL = Ist ncells MH2O/(2F ) − WV

(26)

(18)

(19)

80

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

Q̇ st = (UHHV − Ust ) Ist + (Wmfc + z ) cair Tair − Woda coda Tst − WV (h0 + cV Tst ) − WL cL Tst

̇ = uv We ( Te2 − Tcw ) − mcw Tcw (27)

The cooling temperature difference across the stack ΔTst is controlled for the reference value ΔTst,ref by an underlying controller actuating the cooling mass flow Wi , which has faster dynamics than the system's thermal dynamics. Therefore, Wi is determined by the stationary energy balance in Eq. (28) and is limited to its minimum and maximum values Wmin and Wmax :

(

(

)

Wi = min max Q̇ st/(ci ΔTst,ref ), Wmin , Wmax

)

(28)

2.2.4. Temperature model and ODA-gas water loading All pipes of the cooling system are modeled as continuously stirred tanks with perfect insulation. Conductance UAcw and temperature difference between cooler temperature Tcw and ambient temperature Tu determine the heat flow dissipated by the cooler. Specific heat capacities of the coolants in the internal and external circuit are ci and ce . Coolant masses mi1, mi2 , me1, me2 and mcw in the cooling system are constant. The cooling mass flow in the external circuit We is set constant. The H2-pump for hydrogen recirculation is being cooled at temperature Te1 with a cooling mass flow of 6% We . The H2-pump's heat flow is small as compared to the fuel cell stack's one. The intercooler is a counterflow heat exchanger and exhibits a fast heat transfer due to high coolant mass flows as well as a low thermal mass as compared to the remaining cooling system. Hence, the intercooler's behavior is described by a static model using the Epsilon-NTU-method (Shah & Sekulic, 2003). NTU (Number of Transfer Units) is a coefficient for heat exchanger design. Outlet temperatures are determined explicitly and are gained by a stationary energy balance. The heat flow Q̇ ic exchanged inside the intercooler is gained as the product of the cooling flows inlet temperature difference with the effectiveness ϵNTU and of the lower heat capacity flow Cmin = min (Ci, Ce ). The two heat capacity flows are Ci = Wi ci and Ce = 0.94We ce . The effectiveness ϵNTU is determined by NTU ¼ UA ic/Cmin, by the heat capacity flows Cmax = max (Ci, Ce ) and C⁎ = Cmin/Cmax as well as the conductance UAic :

Q̇ ic = ϵNTU Cmin ( Ti2 − Te1)

(29)

UA cw ( Tcw − Tu ) ce

(35)

Cst Tsṫ = Q̇ st − Wi ci (Tst − Ti1)

(36)

̇ τc Tcond = Tcond,ref − Tcond

(37)

ODA-gas water loading Xoda is determined for a fully saturated flow and the mean pressure in the condenser:

X oda =

(

pVsat Tcond 1 (p 2 3

)

(

+ p4 ) − pVsat Tcond

)

R oda RV

(38)

The thermal dynamics model with the state vector x th and the output vector y th is as follows. The state vector includes stack temperature Tst , stack cooling system temperatures Ti1, Ti2, Te1, Te2, Tcw and condenser cooling temperature Tcond. The output vector comprises all values measured, whereas Ti2 and Te2 are not available through measurement.

̇ = f th ( xg , xth , ug , uth ) xth

(39)

T xth = ⎡⎣ Tst Ti1 Ti2 Te1 Te2 Tcw Tcond ⎤⎦

(40)

T yth = ⎡⎣ Tst Ti1 Te1 Tcw Tcond X oda ⎤⎦

(41)

Model parameters kca , k23, k34 , k 45 and k5u of the pressure model and the parameters ζ1, ζ2, … , ζ4 and Rohm of the voltage model have been identified by a least square error minimization using stationary experimental data with stack current, stoichiometry and stack temperature varied. Model parameters UAic and UAcw of the temperature model have been identified by a least square error minimization using stationary experimental data with coolant temperature and coolant mass flows varied. Parameters Cst and τc have been identified by a least square error minimization using transient experimental data. Experiments have been conducted with the MFFCS test setup shown in Fig. 2. Experimental data for identification is different to verification data. 2.3. Simulation of the MFFCS model

ϵNTU

⎧ NTU C* = 1 , ⎪ ⎪ 1 + NTU =⎨ ⎪ 1 − e−NTU ( 1 − C *) , otherwise ⎪ ⎩ 1 − C *e−NTU ( 1 − C⁎)

(30)

The heat flow from the fuel cell stack into the cooling system is Wi ci (Tst − Ti1) and τc is the time constant of the condenser cooling system that is modeled as a PT1-element. Stack and condenser cooling system temperatures are gained as follows:

⎛⎛ ⎞ Q̇ ⎞ mi1Ti1̇ = Wi ⎜ ⎜ Ti2 − ic ⎟ − Ti1⎟ Ci ⎠ ⎝⎝ ⎠

(31)

mi2 Ti2̇ = Wi ( Tst − Ti2 )

(32)

̇ = We ( uv Tcw + (1 − uv ) Te2 − Te1) me1Te1

(33)

⎛ ⎛ ⎞ ̇ ⎞ ̇ = We ⎜ 94 ⎜ Te1 + Q ic ⎟ + 6 Te1 − Te2 ⎟ me2 Te2 Ce ⎠ 100 ⎝ 100 ⎝ ⎠

(34)

For operation of the MFFCS the fuel cell stack cooling temperature Ti1 is controlled for a reference temperature of 58 °C by an underlying temperature controller operating the cooling valve uv in the external cooling loop (Schultze, Kirsten, Helmker & Horn, 2012). A MFFCS simulation was run at a stoichiometry of λO2 = 1.70 to gain an ODA-gas oxygen content of 10% and a stack current profile as shown in Fig. 6 for stack currents of 25%, 50%, 75% and 100%. Condenser cooling temperature was controlled for 10 °C. Fig. 6 further compares the simulation results with the experimental results. As shown in Fig. 6 the simulation model agrees well with the experimental results: the simulated stack voltage is close to the experimental stack voltage and the simulated temperatures for stack in- and outlet are close to the experimental ones. Temperatures such as the intercooler inlet temperature deviate marginally, however, show the same trend as in the experimental case. The simulated ODA-gas mass flow is a dry gas flow as compared to the experimental one, which is a wet gas flow with low water loading. However, for a water loading of 20 g/kg a deviation between assuming a dry and a humid mass flow is less than 2%. Another reason for deviations in the ODA-gas mass flow are the tolerances in the air mass flow supply by the MFC. An air mass flow slightly higher than requested leads to a higher ODA-

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

81

systems (Haseltine & Rawlings, 2005; Haykin, 2001). An Extended Kalman Filter runs at a low computational cost, however, approximates the model only up to first order. An EKF may not converge to the true system state especially if initial conditions are poorly chosen (Haseltine & Rawlings, 2005). A MHE has the advantage of considering the nonlinear model and to incorporate state constraints, however, has a fairly high computational cost (Allesandri, Baglietto, Battistelli & Zavala, 2010). In contrast to MHE the SPKF works at nearly the same computational cost as the EKF and leads to an approximation of second order. For control of electrical power of a fuel cell system state estimation based on a SPKF has already been applied by Niemeyer (2009). State estimation of a multifunctional fuel cell system using a SPKF algorithm is new. Also, using the SPKF algorithm for a nonlinear model predictive control of the multifunctional fuel cell system is a new topic. To compensate for the significant time delay in the MFFCS Birk (2000) proposes operating the estimator in the past and running a prediction to compensate for the delay. Applying this strategy to the whole system would neglect the available measurement of the thermal system and would only consider the past evolution through the artificially delayed signals y th (t − Td ). The state estimation strategy presented in the following takes advantage of the system model structure. State estimation is performed using y g (t − Td ) and y th (t ) in three consecutive steps (Schultze & Horn, 2013c). 3.1. Sigma-Point-Kalman-Filter The Sigma-Point-Kalman-Filter for nonlinear systems was proposed by Julier, Uhlmann and Durrant-Whyte (1995) and is based on the Kalman-Filter (Kalman, 1960) for linear systems. Weighted sigma points are mapped through the nonlinear and time-discrete system equations F . By this mapping the covariance matrices are approximated. A detailed description of the SPKF is further given by Haykin (2001). The system model with Eq. (42) is assumed to exhibit additive process and measurement noise vk and w k , respectively. Order of the system is n. The noise is assumed to be normally distributed, with zero mean, white and uncorrelated and is described by the covariances Q k and Rk . System equations F result from the continuous system equations f by a numerical time-integration with sampling time Ts using a Runge–Kutta-scheme of fourth order:

x k + 1 = F ( x k , uk ) + v k y k = H ( x k , uk ) + w k

Fig. 6. Simulation results of the multifunctional fuel cell system model compared with experimental results.

gas mass flow and a higher ODA-gas oxygen content as shown. ODA-gas water loading of simulation and experiment are very close. The simulation model describes the dynamics and stationary behavior of the MFFCS well and therefore is suitable for design of a state estimation and a nonlinear model predictive control of dehumidified ODA-gas mass flow.

3. State estimation strategy Due to the MFFCS model's nonlinearity only a few algorithms for state estimation come into consideration. Extended and SigmaPoint-Kalman-Filters (EKF, SPKF) as well as Moving Horizon Estimators (MHE) are candidates for state estimation of nonlinear

(42)

In the original algorithm of the SPKF (Julier et al., 1995) the covariance matrices Pk − 1 and Pk− are decomposed into their square roots Sk − 1 and Sk−. Numerical inaccuracies due to the final representation of numbers can lead to instability of the filter. The Square-Root-Implementation developed by Van Der Merwe and Wan (2001) guarantees positive definiteness and propagates the square roots of the covariance matrices. The Square-Root-Implementation uses three powerful linear algebra techniques. ‘qr’

decomposes a matrix A into an orthogonal matrix Q and an upper triangular matrix R with AT = QR . ‘cholup’ for P ± uuT the rank-1-update of the Cholesky factor S (with P = SST ) is cholup (S , u, ± 1). ‘ /’ is an efficient technique for matrix inversion. S vv and S ww are the square roots of the noise covariances Q k = (S vv )(S vv )T and Rk = (S ww )(S ww )T . Starting from the current estimate x^k − 1 (2n + 1) sigma points are calculated using the square root Sk − 1 of the error covariance matrix. S |i is the i-th column of this matrix:

82

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

⎡⎣ X 0, k − 1… X 2n, k − 1⎤⎦ = ⎡⎣ x^k − 1 x^k − 1 + γS|i x^k − 1 − γS|i ⎤⎦

(43)

Afterwards, the sigma points are mapped through the nonlinear system equations (42). Based on the weights W i(m) and W i(c ) the a − priori-estimate x^ as well as the square root S − of the a priori

by 0 < α ≤ 1 through λ = α 2 (n + κ ) − n and γ = n + λ . Parameters used are κ = 0, α = 0.5 and β = 2 (Haykin, 2001; Niemeyer, 2009):

λ , n+λ

W0(c ) =

error covariance matrix is determined. The case i¼0 is being handled separately by a rank-1 Cholesky factor update due to the negative weight W0(c ):

Wi(m) = Wi(c ) =

1 2 (n + λ )

X i−, k = F ( X i, k − 1, uk − 1)

(44)

3.2. State estimation strategy

(45)

The state estimation strategy is based on the SPKF and utilizes the system model structure as shown in Fig. 3: the gas dynamics f g is independent of the thermal dynamics f th , which on the other

k

k

2n

− x^k =

∑ Wi(m) X i−, k i=0

^− Sk = qr

Sk−

{ ⎡⎣

⎤ − ^− Wi(c ) ( X 1:2 n, k − x k ) S vv ⎦

= cholup

{

^− Sk ,

−W0(c ) ( X 0,− k

}

− − x^k ), − 1

(46)

}

⎡ x^ − ⎣ k

− x^k

+

− γSk−|i x^k

− γSk−|i ⎤⎦

(47)

(48)

−

The estimate of the outputs y^k is gained by mapping the sigma points through the output-equations H into the values Y i−, k and summing these up:

Y i−, k = H ( X i−, k , uk )

− y^k =

(49)

Wi(m) Y i−, k

(50)

i=0

The correction step requires the cross-covariance matrix P −^^ as xy − of the output-error-covariance well as the Cholesky factor S yy matrix. The case i¼ 0 is handled separately by a rank-1-update:

^− S yy ^ ^ = qr

{ ⎡⎣

S −yy ^ ^ = cholup

⎤ − ^− Wi(c ) ( Y 1:2 n, k − y k ) S ww ⎦

{

^− ^ ^, S yy

}

− −W0(c ) ( Y 0,− k − y^k ), − 1

2n

− = P xy ^^

(

}

(52)

− T

−

∑ Wi(c) ( X i−, k − x^k )( Y i−, k − y^k ) i=0

− /S − T P xy ^^ yy ^^

)

/S −yy ^^

(53)

{

x^k =

(58)

hand depends on the gas dynamics. State estimation involves three steps in order to compensate for the time delay. (1) State estimation of gas dynamics (τ = t − Td ):

ficially delayed input ug (τ ) = ug (t − Td ), which is available through an internal storage:

ẋg = f g ( xg , ug (τ )) yg = hg ( xg , ug (τ ))

(59)

The nonlinear and time-discrete system equations F g for the SPKF are gained through numerical time-integration of Eq. (59) for one sampling time step Ts based on a Runge–Kutta-scheme of fourth order with constant step size. State estimation is performed for dimension n ¼4 of x g and delivers x^ g (t − Td ).

(

−

+ Lk yk − y^k

)

}

^ An estimate x^ g (τ = t ) of the gas states at time t is gained through prediction over the time interval t − Td < τ ≤ t of the delay using the input ug . Initial condition for the numerical timeintegration is the state x^ (τ ) estimated in step one for τ = t − T : g

0

0

d

^̇ ^ x^ g (τ ) = f g ( x^ g (τ ), ug (τ ))

(60)

^ x^ g (τ 0 ) = x^g (τ 0 )

(61)

Piecewise constant inputs [ ug, Nk , … ug,1] are stored with ug, k = ug (t − kTs ≤ τ < t − (k − 1) Ts ) over the time interval. This vector consists of Nk = (Td/Ts ) elements. Ts is chosen such that Nk is an integer number. System state at time τ = t is gained by Nk piecewise numerical time-integrations of equations (60) by a ^ Runge–Kutta-scheme of fourth order. The predicted state is x^ g (t ). (3) State estimation of thermal dynamics (τ = t ):

(55)

The model for thermal dynamics state estimation is based on Eq. (62) with input vector u = [ ug , uth ] and output vector y th . Due to the dependency of the thermal dynamics on the gas dynamics, the latter ones are time-integrated along with the thermal dynamics by incorporating the state vector ζ g . As gas dynamics states are not estimated in this step, they are labeled as ζ g :

(56)

⎡ ζ ̇ ⎤ ⎡ f g ( ζ g , u) ⎤ ⎥ ⎢ g⎥ = ⎢ ⎢⎣ xth ̇ ⎥⎦ ⎢⎣ f th ( ζ g, xth , u)⎥⎦

(54)

Sk = cholup Sk−, Lk S −yy ^ ^, − 1

− x^k

(i = 1, 2, …, 2n)

(51)

The state estimate x^k is gained using the Kalman-gain Lk as well as the deviation between current and estimated measurement. The Cholesky factor Sk is updated by a rank-1-update:

Lk =

(57)

(2) Prediction of gas dynamics (t − Td < τ ≤ t ):

2n

∑

λ + (1 − α 2 + β ) n+λ

Gas dynamics state estimation is done at time τ being shifted into the past by the delay Td . The system model f g has the arti-

For a higher accuracy in case of additive process noise the sigma points are redrawn:

⎡⎣ X 0,− k … X 2−n, k ⎤⎦ =

W0(m) =

The distance of the sigma points to the current estimate is defined

yth = hth ( ζ g, xth , u)

(62)

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

83

Simulink and was successfully applied experimentally. The experiments were conducted with a constant stoichiometry of λO2 = 1.70. Condenser cooling temperature was controlled for 10 °C to gain a low water loading. A stack current profile was applied such that the current changes with 5 A/s. Fig. 7 shows the input signals to the fuel cell system as well as the simulation and experimental results. As shown, the SPKF estimates the ODA-gas mass flow and oxygen concentration very well. The mass flow disturbance z is different from zero showing that the model and the plant deviate slightly from each other. This deviation is due to tolerances in air mass flow supply, which can be seen by the ODA-gas oxygen content that is slightly higher than the expected value of circa 10%. An approximation of the current ODA-gas mass flow and therefore a compensation of the time delay is obtained by the prediction in step two of the state estimation strategy. However, model/plant mismatches are not updated during the prediction, which can lead to a slight inaccuracy between measurement and state estimation as shown during transients (e.g. for mass flow z). As further shown the temperatures are estimated very well. ODAgas water loading is estimated very well, too. Its estimation deviates only slightly from the measurement. As shown the state estimation strategy delivers an estimate that keeps track of the system evolution. The state vector estimate is handed over to the nonlinear model predictive control that requires the current system state.

4. NMPC of ODA-gas mass flow

Fig. 7. State estimation results showing the stack current as input and the experimental measurement as well as the estimation of ODA-gas mass flow, oxygen content as well as ODA-gas water loading and temperatures and the estimation of disturbance mass flow.

The nonlinear time-discrete system model for the SPKF is gained through numerical time-integration of Eq. (62) for one sampling time step Ts based on a Runge–Kutta-scheme of fourth order with constant step size. The time evolution of the gas dynamics is being considered starting from the previous step. State estimation of the thermal dynamics is performed for dimension n ¼ 7 of x th and ^ delivers x^ th (t ). ^T ^T The system state estimate is [ x^ g x^ th ]T . 3.3. SPKF simulation and experimental results The state estimation strategy has been implemented in Matlab/

In model predictive control (MPC) a finite horizon open-loop optimal control problem is solved online for every sampling time step and yields an optimal control sequence from which the first element is applied. The current system state, which is obtained by measurement or estimation, serves as initial state. MPC's major advantage is handling constraints on states and inputs (Mayne, Rawlings, Rao & Scokaert, 2000). MPC is based on linear constraints and a linear model describing the system dynamics. Nonlinear model predictive control (NMPC) refers to MPC schemes being applied to nonlinear constraints and nonlinear system dynamics (Allgöwer, Findeisen & Nagy, 2004). Generally, NMPC problems cannot be solved analytically, calling for numerical algorithms. Optimal control problems can be solved by dynamic programming, direct and indirect methods. In dynamic programming the computational cost grows with the problem's dimension (Kirk, 2004). Its application is limited to problems with low order. Indirect methods are based on the first order optimality conditions of calculus of variations and involve solving a two-point boundary value problem, which may require significant effort. NMPC problems are widely solved by direct methods that convert the optimal control problem into a nonlinear programming problem, which is solved iteratively by sequential quadratic programming (SQP). Inequality constraints can easily be included. A comprehensive overview over NMPC is given in Grüne and Pannek (2011). The NMPC of ODA-gas mass flow (Schultze et al., 2014) actuates stack current Ist and condenser cooling reference temperature Tcond,ref . For satisfying operational limits and to prevent stack damage by cathode flooding or oxygen starvation Ist and its gradient are limited. The value Tcond,ref shall not fall below the user defined value Tcond,l and shall not be below 5 °C to prevent icing. The control objective and the limitations are translated into the cost function and the constraints, respectively. The NMPC is run with a sampling time of Ts = 0.5 s.

84

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

4.1. NMPC system model

2

2

J ( uD ) = q11 ( x2, N − x2,ref ) + q12 ( yN )

The NMPC system model has a low dimension to reduce the computational cost. Due to the underlying stack cooling temperature control actuating the cooling valve uv the thermal model can be reduced such that the stack cooling system is neglected. However, fuel cell stack temperature dynamics are important and remain in the NMPC system model. The NMPC system model describes the gas dynamics as given by Eqs. (1)–(6) as well as the thermal dynamics of the fuel cell stack and the condenser as given by Eqs. (12)–(28) and (36)–(38). Ti1 is handled as a constant input over the NMPC horizon. The model input vector is unmpc (t ) = [u1 u2 u3 ]T = [Ist λO2 Tcond,ref ]T and the output is ODA-gas water loading Xoda . The NMPC system equations are f nmpc with states x nmpc . The entire NMPC system model is given by (64)–(65) as follows: T xnmpc = ⎡⎣ Wmfc Woda cO2 z Tst Tcond ⎤⎦

(63)

ẋnmpc = f nmpc ( xnmpc , unmpc )

(64)

ynmpc = hnmpc ( xnmpc , unmpc ) = X oda

(65)

4.2. Prediction of current system state The NMPC receives x nmpc from the SPKF state estimation at time t, whereas the NMPC result is applied at time t + Ts . This would cause a delay of one sampling time period. An approximation of the system state x^ (t + T ) at time t + T is generated nmpc

s

s

by prediction of the system model through numerical time-integration with a Runge–Kutta-scheme of fourth order and constant step size over the time interval [t t + Ts ] based on the current system input. The predicted system state is x = x^ (t + T ) . 0

nmpc

N−1

+ Ts

∑ k=0 N−1

+ Ts

∑ k=0

2 2⎤ ⎡q x − x 2,ref ) + q22 ( yk ) ⎦ ⎣ 21 ( 2, k 2⎤ 2 ⎡q u − u 1,ss ) + q24 ( u3 ) ⎦ ⎣ 23 ( 1, k

(66)

J ( uD ) combines tracking control of ODA-gas mass flow for reference value x2,ref = Woda,ref by simultaneously requiring a low ODA-gas water loading. Input variables are penalized to prevent their heavy use. Weighing factors q11, q12, q21, … , q24 influence the balance between control performance and accuracy in steady state. They are adjusted such that q11 and q21 are chosen larger than q12 and q22 to achieve close reference tracking of the ODA-gas mass flow. Weighing factors q23 and q24 influence reference tracking as well. They are chosen to achieve smooth control variable-signals and to gain close reference tracking of the ODA-gas mass flow. The weighing factors have been adjusted through an iterative procedure that includes simulation studies. This procedure started with a default set of weighing factors. In consecutive iterations the weighing factors were adjusted aiming for a smooth stack current curve and close reference tracking. Due to model/ plant-deviations and disturbances steady state accuracy cannot be achieved (Pannocchia & Rawlings, 2003). A concept that allows for steady state accuracy is presented by Morari and Maeder (2012). A stationary analysis of the system model under influence of the disturbances leads to stationary inputs that are incorporated into the cost function such that J ( uD ) rates the deviation between this stationary value and the input. In the following, disturbance z is not considered to prevent nonlinear effects due to noise. A stationary analysis of Eqs. (1)–(2) using the ODA-gas mass flow reference value x2,ref leads to u1,ss that is incorporated in the cost function to improve steady state accuracy for reaching x2,ref :

u1,ss =

x2,ref 4F ⎛ ⎞ 0.79 ncells ⎜ (λ O2 − 1) MO2 + λ O2 MN ⎟ ⎝ 0.21 2 ⎠

(67)

s

4.3. Nonlinear optimization problem An introduction to the following NMPC algorithm is given in Grüne and Pannek (2011). In NMPC the continuous system model is discretized assuming that the inputs uk are constant over one sampling time period Ts . The piecewise time-integration of the continuous model over the time intervals [kTs (k + 1) Ts ] delivers the states x k + 1 for k = 0, … , N − 1. The prediction horizon has a length of NTs . Numerical time-integration is performed by a Runge–Kutta-scheme of fourth order and leads to the time-discrete model equations x k + 1 = F D ( x k , uk ) and the output yk = hnmpc ( x k , uk ). Starting with the initial state x 0 and the initial

The nonlinear optimization problem of minimizing cost function J ( uD ) subject to system dynamics (69) and constraints (70)–(74) is as follows. The linear inequality constraints (72)–(74) are derived by approximation as finite differences. The input value u˜1 is the input applied in the previous sampling time step. The optimal sequence uD⁎ leads to the minimal cost J⁎ ( uD⁎ ) and is the solution of ⁎ and u3⁎ of the solution the optimization problem. The elements u1,0 ⁎ ⁎ vector and λO2 are applied as input signal unmpc = [u1,0 λO2 u3⁎ ]T . The optimizationproblem (68)–(74) is solved by sequential quadratic programming using the solver solnp (Ye, 1989):

min J ( uD ) uD

(68)

x k + 1 = F D ( x k , uk )

(69)

Ist,min ≤ u1, k ≤ Ist,max

(70)

max ( 5, Tcond,l ) ≤ u3 ≤ Tcond,max

(71)

u˜1 − ΔImax Ts ≤ u1,0 ≤ u˜1 + ΔImax Ts

(72)

−ΔImax Ts ≤ u1,1 − u1,0 ≤ + ΔImax Ts

(73)

]T

output y0 there are N additional state vectors x k = [x1, k x2, k…x6, k and outputs yk . Stoichiometry is set constant such that the independent system inputs are u1 and u3. These are the optimization variables of the nonlinear optimization problem. Furthermore, the slow dynamics of condenser cooling are exploited to reduce the number of optimization variables and hence computational cost. Input variable u3 is set constant over the entire prediction horizon, so there are N identical input variables u3. In total, there are N values of u1, k and a scalar u3 as optimization variables. Instead of 2N optimization variables, their number has been reduced to N + 1. These variables are grouped in uD = [u1,0… u1, N − 1 u3 ]T . To reach the control objective the cost function J ( uD ) rates inputs and states:

… − ΔImax Ts ≤ u1, N − 1 − u1, N − 2 ≤ + ΔImax Ts.

(74)

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

85

Fig. 9. Simulation and experimental results for the NMPC of the MFFCS for stack current gradients of ΔImax = 5, 10, 20 A/s.

Fig. 8. Simulation and experimental results of the NMPC.

4.4. NMPC simulation and experimental results The NMPC has been run in combination with the SPKF at the entire system model in Matlab/Simulink as well as at the MFFCS experimental setup. Stoichiometry was set λO2 = 1.70 to gain an ODA-gas oxygen content close to 10%. ODA-gas mass flow reference value x2,ref = Woda,ref was varied between 0.10 and 0.167 g/s/cell and stack current gradient was limited to ΔImax = 5 A/s. For condenser cooling a value of Tcond,l = 5 °C was set. The simulation was run with z¼0 and a prediction horizon of N¼ 9 leading to a time horizon of NTs = 4.5 s. This prediction horizon was chosen as a good compromise between computational burden and accuracy of the controller. The state estimation algorithm has a computation time of less than 5 ms and the NMPC routine takes less than 50 ms to return its result. The entire algorithm of state estimation and NMPC finishes in less than a sampling time of 500 ms and hence shows the real time capability of the entire algorithm. Fig. 8 shows simulation as well as experimental results of the NMPC for ODA-gas mass flow of the MFFCS. ODA-gas mass flow is controlled for with high accuracy. Furthermore, limitation of the stack current gradient is satisfied.

ODA-gas water loading is low and matches between simulation and experiment. Deviations occur between simulation and experiment for ODA-gas oxygen content and stack current. This is due to the disturbance mass flow occurring in the experiments. An increased supply of feed air mass flow due to tolerances in the air supply leads to an increase in ODA-gas oxygen content and results in a lower stack current to reach the reference value. Fig. 9 shows simulation and experimental results for the NMPC run at stack current gradients of ΔImax = 5, 10, 20 A/s. Obviously, the mass flow reference value is reached quicker for an increased gradient ΔImax . Due to incorporation as constraint the NMPC adapts to the value of ΔImax and leads to the quick response as shown in Fig. 5 with the advantage of closely reaching the reference value. Control of the MFFCS ODA-gas mass flow could be realized by classical linear control approaches such as proportional-integral control or state space control (Schultze & Horn, 2013b) as well. For controller design a linear model would be derived from the nonlinear model. However, the linear model looses its validity over the entire operating region. Further, to satisfy the constraints on the stack current and its gradient a limiter on the controller output would have to be added. Parametrization of the linear controller for variable maximum stack current gradients would be difficult. Parametrization for a certain stack current gradient may lead to the case that the dynamics of the controlled system saturate at a parameter-specific maximum stack current gradient and higher gradients could not be reached. NMPC eases the incorporation of constraints significantly and has several advantages as compared to classical linear control techniques. Due to the use of the nonlinear model, NMPC is valid over the entire operating range of the MFFCS. Further, the NMPC of the MFFCS shows very good results even if the maximum stack current gradient is varied as it incorporates constraints and solves the optimization problem with updated constraints every sampling time period. In combination with the state estimation the NMPC is a successful algorithm for controlling a MFFCS for ODA-gas mass flow with a low water loading by simultaneously satisfying constraints on the inputs. 5. Conclusions For a multifunctional fuel cell system a nonlinear model predictive control of the dehumidified exhaust gas mass flow with

86

M. Schultze, J. Horn / Control Engineering Practice 49 (2016) 76–86

low oxygen content was presented. The paper showed all stages involved in designing the controller. These stages cover system modeling, state estimation and controller synthesis. Nonlinear model predictive control requires the current system state that is gained by a nonlinear state estimation based on the Sigma-PointKalman-Filter. The algorithm presented further compensates the significant delay in the measurement. The fuel cell system model presented serves for design of the control algorithm and the state estimation strategy.

Acknowledgments This study as part of the project “Cabin technology and multifunctional fuel cell systems” has been supported by Airbus and the German Federal Ministry of Education and Research (support code: 03CL03A).

References Allesandri, A., Baglietto, M., Battistelli, G. & Zavala, V. (2010). Advances in moving horizon estimation for nonlinear systems. In 2010 49th IEEE conference on decision and control, Atlanta, Georgia (pp. 5681–5688). Allgöwer, F., Findeisen, R., & Nagy, Z. K. (2004). Nonlinear model predictive control: From theory to application. Journal of Chinese Institute of Chemical Engineers, 35 (3), 299–315. Amphlett, J. C., Baumert, R. M., Mann, R. F., Peppley, B. A., et al. (1995). Performance modeling of the Ballard Mark IV solid polymer electrolyte fuel cell I. Mechanistic model development. Journal of the Electrochemical Society, 142(1), 1–8. Aviation Rulemaking Advisory Committee (1998). Fuel tank inerting – Task group 3. Fuel Tank Harmonization Working Group, Final report. Baehr, H. D., & Kabelac, S. (2012). Thermodynamik: Grundlagen und technische Anwendungen [English translation – Thermodynamics: Fundamentals and technical applications] (15th ed.). Berlin Heidelberg: Springer-Verlag (in German). Birk, J. (2000). Zustandsschätzung für nichtlineare Prozesse mit zeitdiskreten und totzeitbehafteten Messgrössen [English title: State estimation for nonlinear processes with discrete-time and time-delayed measurement]. at-Automatisierungstechnik, 48 (5), 235–239 (in German). Borup, R., Meyers, J., Pivovar, B., Kim, Y. S., Mukundan, R., Garland, N., et al. (2007). Scientific aspects of polymer electrolyte fuel cell durability and degradation. Chemical Reviews, 107 (10), 3904–3951. Ford, J. A. (1995). Improved algorithms of Illinois-type for the numerical solution of nonlinear equations. Wivenhoe Park, Colchester, Essex, UK: University of Essex. Grüne, L., & Pannek, J. (2011). Nonlinear model predictive control: Theory and algorithms. London: Springer-Verlag. Haseltine, E. L., & Rawlings, J. B. (2005). Critical evaluation of extended Kalman filtering and moving-horizon estimation. Industrial & Engineering Chemistry Research, 44(8), 2451–2460. Haykin, S. (2001). Kalman filtering and neural networks. New York: John Wiley & Sons. Julier, S. J., Uhlmann, J. K. & Durrant-Whyte, H. F. (1995). A new approach for filtering nonlinear systems. In Proceedings of the 1995 American control conference, Seattle, Washington (Vol. 3, pp. 1628–1632). IEEE. Kallo, J., Renouard-Vallet, G., Saballus, M., Schmithals, G., Schirmer, J., Friedrich, K. A., et al. (2010). Fuel cell system development and testing for aircraft applications. In Proceedings of the 18th world hydrogen energy conference 2010 – WHEC

2010, Essen, Germany. Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35–45. Karnik, A. Y., Sun, J., Stefanopoulou, A. G. & Buckland, J. H. (2009). Humidity and pressure regulation in a PEM fuel cell using a gain-scheduled static feedback controller. IEEE Transactions on Control Systems Technology, 17 (6), 283–297. Kirk, D. E. (2004). Optimal control theory: An introduction. Mineola, New York: Dover Publications. Mayne, D. Q., Rawlings, J. B., Rao, C. V., & Scokaert, P. O. M. (2000). Constrained model predictive control: Stability and optimality. Automatica, 36(6), 789–814. McLoughlin, A. (2009). More electric – ready for take off? In EPE'09, 13th European conference on power electronics and applications 2009, Barcelona, Spain (pp. 1–7). IEEE. Morari, M., & Maeder, U. (2012). Nonlinear offset-free model predictive control. Automatica, 48(9), 2059–2067. Müller, E. A., Stefanopoulou, A. G. & Guzzella, L. (2007). Optimal power control of hybrid fuel cell systems for an accelerated system warm-up. Transactions on Control Systems Technology, IEEE, 15 (2), 290–305. Niemeyer, J. (2009). Modellprädiktive Regelung eines PEM-Brennstoffzellensystems [English translation: Model predictive control of a PEM fuel cell system]. Dissertation, University of Karlsruhe (in German). O'Hayre, R., Cha, S.-W., Colella, W., & Prinz, F. B. (2009). Fuel cell fundamentals. Hoboken, NJ: John Wiley & Sons. Pannocchia, G., & Rawlings, J. B. (2003). Disturbance models for offset-free modelpredictive control. AIChE Journal, 49(2), 426–437. Pukrushpan, J. T., Stefanopoulou, A. G., & Peng, H. (2004). Control of fuel cell power systems: Principles, modeling, analysis and feedback design. London: SpringerVerlag. Schultze, M. & Horn, J. (2013a). Nonlinear model predictive control of PEM fuel cell systems for generation of exhaust gas with low oxygen content. In 19th IFAC symposium on automatic control in aerospace, Würzburg, Germany (pp. 470– 475). IFAC. Schultze, M. & Horn, J. (2013b). State estimation with time delay and state feedback control of cathode exhaust gas mass flow for PEM fuel cell systems. In European control conference (ECC) 2013, Zurich, Switzerland (pp. 3560–3565). IEEE. Schultze, M. & Horn, J. (2013c). State estimation for PEM fuel cell systems with time delay by an unscented Kalman filter and predictor strategy. In 21st Mediterranean conference on control & automation (MED) 2013, Barcelona, Spain (pp. 104– 112). IEEE. Schultze, M., Kirsten, M., Helmker, S. & Horn, J. (2012). Modeling and simulation of a coupled double-loop-cooling system for PEM-fuel cell stack cooling. In 2012 UKACC international conference on control, Cardiff, UK (pp. 857–863). Schultze, M., Hähnel, C. & Horn, J. (2014). Nonlinear model predictive control of a PEM fuel cell system for cathode exhaust gas generation. In 19th IFAC world congress 2014, Cape Town, South Africa (pp. 9432–9437). Shah, R. K., & Sekulic, D. P. (2003). Fundamentals of heat exchanger design. Hoboken, New Jersey: Jon Wiley & Sons. Tomlinson, S., Barker, M., Venn, D., Hickson, L. & Lam, J.K.-W. (2011). Mathematical model of water contamination in aircraft fuel tanks. SAE technical paper, 201101-2540. Vahidi, A., Stefanopoulou, A. G. & Peng, H. (2004). Model predictive control for starvation prevention in a hybrid fuel cell system. In Proceedings of the 2004 American control conference, Boston, Massachusetts (Vol. 1, pp. 834–839). IEEE. Vahidi, A., Stefanopoulou, A. G. & Peng, H. (2006). Current management in a hybrid fuel cell power system: A model-predictive control approach. Transactions on Control Systems Technology, IEEE, 14 (6), 1047–1057. Van Der Merwe, R. & Wan, E. A. (2001). The square-root unscented Kalman filter for state and parameter-estimation. In Proceedings of the IEEE international conference on acoustics, speech and signal processing, Salt Lake City, Utah, 2001 (Vol. 6, pp. 3461–3464). IEEE. Ye, Y. (1989). SOLNP user's guide – A nonlinear optimization program in MATLAB. Department of Management Sciences, University of Iowa.