4264
Ind. Eng. Chem. Res. 1997, 36, 4264-4272
Parameter and State Estimation of a Proton-Exchange Membrane Fuel Cell Using Sequential Quadratic Programming G. E. Suares and K. A. Kosanovich* Department of Chemical Engineering, University of South Carolina, Columbia, South Carolina 29208
All mathematical models contain parameters that must be determined for the model to represent accurately the behavior of the system. The parameter estimation problem is usually solved as an unconstrained optimization problem independent of the model equations. However, by integrating the parameter estimation problem with the generation of the model’s state profiles, constraints can be embedded directly into the optimizer, and an infeasible path solution approach can be used. Nonlinear programming is the ideal framework for formulating constrained optimization problems. The model is introduced into this framework as constraints using orthogonal collocation on finite elements. The resulting nonlinear programming problem is then solved using sequential quadratic programming. We demonstrate this approach on a mathematical model of a proton-exchange-membrane fuel cell in which four parameters are estimated and nine states’ profiles are determined. 1. Introduction There is an ever increasing demand for cost effective, alternative power sources that are both energy efficient and environmentally benign. There has been considerable interest in fuel cell systems, particularly the proton-exchange membrane (PEM) fuel cell. It has been reported that these fuel cells can be used for transportation applications and local power generation because of their low-temperature operation and ease of construction. However, the current designs are limited, for example, by their inability to be operated efficiently at a high power density. In order to optimize and improve fuel cell design, repeated experimentation must be carried out, which may be costly and time consuming. One viable alternative is to use a mathematical model to uncover the present design limitations and to determine where improvements can be made. Modeling of the PEM fuel cell can be done either mechanistically or empirically; regardless, it is a complex and difficult task (Amphlett et al., 1995a,b). The open literature contains several different models investigating different aspects of the PEM fuel cell. For example, Bernardi and Verbrugge focused on studying the transport mechanisms across the membrane while Nguyen and White and Fuller and Newman studied the importance of water and thermal management of the fuel cell (Bernardi and Verbrugge, 1992; Nguyen and White, 1993; Fuller and Newman, 1993). All of these models contain parameters, some of which can be measured readily while others must be estimated from experiments. It is attractive from the viewpoint of model validation to have a minimal number of parameters that may be adjusted to improve the fit of the model to the available experimental data. However, a parameter value obtained under ideal conditions may not be correct under normal operating conditions. This is often the case for parameters such as the diffusion coefficient, which is often replaced by an effective diffusion coefficient (Bernardi and Verbrugge, 1991). There is, therefore, a need for an efficient methodology * To whom correspondence should be addressed. Phone: (803) 777-0143. Fax: (803) 777-8265. E-mail: kosanoka@ sun.che.sc.edu. S0888-5885(97)00251-0 CCC: $14.00
to determine these parameters such that the model’s predictions accurately represent the behavior of a fuel cell. This paper is organized as follows. First, the particular fuel cell model used in this work is presented and compared to a number of alternative models. Next, the parameter estimation methodology and nonlinear programming formulation are presented. Last, simulation results are discussed when up to four adjustable parameters are estimated simultaneously with the state solution trajectories of the PEM fuel cell model. 2. PEM Model The model used in this work was developed by Nguyen and White. This “pseudo” two-dimensional (one-dimensional with coupled equations between the channels), steady-state model accounts for heat and water effects in a PEM fuel cell. It includes water transport across the membrane by electro-osmosis and diffusion, heat transfer from the solid to gas phase, and latent heat associated with water evaporation and condensation in the flow channels (see Figure 1). The model consists of material balances for four components that change due to variations in the normal flux across the membrane caused by the change in the current density along the length of the flow channels. The model (see Table 1) contains nine governing equations for the nine state (dependent) variables representing molar flow rates, temperatures, and the local current l v ˙ w,a , M ˙ w,a , M ˙ O2, M ˙ lw,c, M ˙ vw,c, Ta, Tc, and density: M ˙ H 2, M I(x) (Nguyen and White, 1993). The definition of the terms used can be found in the symbol section of the paper. The feed stream components, hydrogen fuel and air/ oxygen, are usually humidified to maintain proper membrane conditions. Because of differences between the partial and vapor pressures of water, it is necessary to account for water in both the liquid and vapor phases. The condensation and evaporation of water is dependent upon kc, the homogeneous rate constant. This is one parameter that can be adjusted so that the model matches the experimental data. It was assumed in this model that gas diffusion through the electrode porous layer can be neglected; hence the flux of hydrogen and oxygen can be character© 1997 American Chemical Society
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4265
Figure 1. Schematic of modeled regions of a PEM fuel cell. Table 1. Governing Equations for the Nguyen and White Model dM ˙ H2 dx
)-
wI(x) 2F
(
)(
dM ˙ lw,a kcwh M ˙ vw,a P - Psat ) w,a v dx R(Ta + 273) M ˙ w,a + M ˙ H2
)
dM ˙ vw,a dM ˙ lw,a RI(x)w )dx dx F
∑ (M˙ C i
p,i)
dTa dx
i
dM ˙ O2 dx
)-
) (Hvw,a - Hlw,a)
dM ˙ lw,a dx
+ Ua(Ts - Ta)
wI(x) 4F
(
)(
l v dM ˙ w,c kcwh M ˙ w,c sat P - Pw,c ) v dx R(Tc + 273) M ˙ w,c + M ˙ O2 + M ˙ N2
)
v dM ˙ w,c dM ˙ lw,a (1 + 2R)I(x)w )+ dx dx 2F
∑
(M ˙ iCp,i)
dTc
i
Vcell ) Voc -
dx
v l - Hw,c ) ) (Hw,c
l dM ˙ w,c
dx
+ Ua(Ts - Tc)
( )
R(273 + Ts) I(x)tm I(x) ln o 0.5F σm(x) I PO2(x)
ized completely by the corresponding local current density. However, the flux of water vapor depends on both migration and diffusion. In this model, the diffusion coefficient of water in the membrane is dependent on both temperature and the water content in the membrane as well as a diffusion parameter. Since this formulation was never validated, the diffusion parameter, D°, can be considered an adjustable parameter. Temperature effects, which are often overlooked, have been modeled for the anode and cathode flow streams along the length of the channel. The temperature of
the solid components of the fuel cell is assumed to be constant due to high thermal conductivities of the materials. The energy balances account for condensation or evaporation and conduction between the solid phase and the corresponding flow channel of each electrode. The value of the overall heat-transfer coefficient, U, used by Nguyen and White is based on an iron/water vapor system. This value is an approximation and may be used as an adjustable parameter to fit the model to experimental data. The cell potential is calculated as a function of electrode polarization and membrane resistance. The cell overpotential is dependent upon the oxygen exchange current density and the partial pressure of oxygen in the cathode stream. The exchange current density, I°, is temperature dependent. This implies that for different operating conditions I° will vary. Thus, even if I° is known for one operating condition, it must be estimated for others. 2.1. Comparison to Other Models. The model proposed by Fuller and Newman is also two-dimensional, with a focus on water and thermal management as well as fuel utilization (Fuller and Newman, 1993). Concentrated solution theory is used to represent the transport phenomena in the polymer electrolyte. Three of the four adjustable parameters discussed above appear in a similar fashion in this model. The fourth parameter, D°, is lumped into a constant that may be treated as a fitting parameter. Other fuel cell models have been developed to emphasize different operational and design aspects. In the model by Bernardi and Verbrugge, the focus is on the identification of cell performance limitations and model validation of the transport mechanisms (Bernardi and Verbrugge, 1991, 1992). In contrast to the Nguyen and White model, this model is isothermal and one-dimensional across the membrane. Only one adjustable parameter, the exchange current density, is modified to fit experimental data at one operating condition. Using this fixed value, the model is used to predict cell performance at other operating conditions. Comparisons of the model predictions to experimental data, using the parameter fixed at one value, indicate that model accuracy may be improved by adjustment of this or other parameters. This is corroborated by the work of Springer et al., who advocate model validation and parameter estimation that satisfy a family of operating conditions simultaneously (Springer et al., 1993). In the Springer et al. model, the focus is on the modeling of the water content and transport phenomena across the PEM fuel cell (Springer et al., 1991). They provide an isothermal, one-dimensional, steady-state model. This model and the one by Nguyen and White both use empirical equations to describe the electroosmotic drag and diffusion coefficient of water as a function of water vapor activity. D° can be considered an adjustable parameter in the Springer et al. model because it would appear as a parameter in the expression to determine the diffusion coefficient. The adjustable parameter I° also appears in their overpotential expression. Amphlett et al. present a mechanistic model for which there are a number of parameters that are not readily available (Amphlett et al., 1995a,b). To circumvent this, they advocate empirical modeling where appropriate because empiricism usually fits the available experimental data well. By combining both mechanistic and empirical model development, they show that extrapola-
4266 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997
Figure 2. Black box approach.
tion beyond the experimental data range is possible. Alternatively, there are methods that can provide parameter estimates without generating the empirical relationships. The parameter I° is included in their model and appears as a lumped parameter in the overpotential expression. 3. Methodology The conventional approach to the parameter estimation problem is to formulate it as an unconstrained optimization problem with a least-squares objective function, and to apply a black-box method that employs an ordinary differential equation (ODE) solver and an optimizer (see Figure 2) (Tjoa and Biegler, 1991; Dennis and Schnabel, 1983). The choice of the ODE solver is independent of the optimizer and is based solely on the behavior of the ODEs (e.g., stiffness). The updates to the parameters are found on the basis of the minimization of a suitably defined objective function. The usual choice of the update rule is a gradient-based method that uses numerical derivatives. In most cases, this technique is effective; however, because there is minimal communication between the solver and the optimizer, there is a high probability that some implicit constraint may be violated, leading to erroneous solutions. Another disadvantage is that process constraints cannot be included in the formulation (Tjoa and Biegler, 1991). An alternative to the black-box method is the all-atonce approach, which eliminates these problems by directly integrating the ODEs into the optimizer. The method of orthogonal collocation, along with polynomial finite element basis functions (orthogonal collocation on finite elements), is used to reduce the governing differential equations to a set of nonlinear algebraic residual equations (NAEs) containing state variables. The discretized model can then be embedded directly into the optimizer as constraints. For a given set of parameters, the NAEs may be solved for the unknown collocation coefficients, but this would only amount to substituting one black-box method for another. The black-box method follows what is known as a feasible path. For each iteration the model equations are satisfied for a corresponding set of parameters. Consider a very nonlinear constraint manifold and suppose the current point p0 is far from the solution p*. The feasible path solution requires that updates of p remain on the constraint manifold. However, it is easy to see, in the simple illustration shown in Figure 3, that by allowing infeasible movements off of the manifold, p* may be reached more efficiently. The all-at-once technique incorporates an infeasible path approach by treating the collocation coefficients as independent parameters. In this formulation, the optimizer not only determines the next iterate of parameters but also the next iterate of collocation coefficients. Thus, the model, written as equality constraints, is required to be valid only at the final solution. This formulation is known as a nonlinear programming (NLP) problem (Tjoa and Biegler, 1991; Edgar and Himmelblau, 1988; Luenberger, 1989) and can be solved by sequential quadratic programming (SQP) or reduced gradient approaches.
Figure 3. Hypothetical nonlinear constraint manifold.
In this work, the parameter estimation and state trajectory solution is formulated as an NLP and the SQP method is used to solve it. 4. Solution Procedure The parameters to be estimated are kc, the evaporation and condensation homogeneous rate constant; D°, a parameter in the diffusion coefficient expression for water; U, the overall heat-transfer coefficient; and I°, the oxygen exchange current density. There are eight states that are represented by ODEs and an algebraic expression for cell potential, for a total of nine state trajectories to be determined. The optimization problem is the minimization of an objective function defined as the sum squared errors between the data and the model’s predictions, r
min Φ(θ) ) θ
eTµ Weµ ∑ µ)1
(1)
subject to satisfying the model equations. θ represents the vector of parameters, e represents the vector of errors, W is a weighting matrix among the states, and Φ is the objective function. The minimization is achieved by selecting the optimal set of parameters. 4.1. Collocation. To use the NLP formulation, which allows the constraints to be embedded into the optimization, the ODEs are converted into AEs by orthogonal collocation on finite elements (Tjoa and Biegler, 1991). The state profiles are approximated by Lagrange polynomials of the form K
ziK(x) ) K
Ψ[ij](x) )
z[ij]Ψ[ij](x) ∑ j)0
(x - x[ik])
∏ k)0,j (x
K
) [ij]
- x[ik])
(τ - τk)
∏ k)0,j (τ
j
(2)
- τk)
x[ij] ) νi + τj(νi+1 - νi) ) νi + τj∆νi where z are the polynomial coefficients that are equal to the states at the collocation points, νi is the location of the ith finite element, τj is a collocation point found from the shifted roots of the orthogonal Legendre polynomial of order K, and x is the independent coordinate that represents the length along the electrode. The ODEs can be discretized and represented as residual equations, K
R(x[ij]) )
∑ z[ik]Ψ˙ k(τj) - ∆νiF(z[ij],θ,x[ij]) ) 0
k)0
i ) 1, 2, ..., NE
j ) 1, 2, ..., K
(3)
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4267
where NE is the number of finite elements and K is the number of collocation points in each element. This defines a set of optimization constraints that must be satisfied at the collocation points. Continuity of the piecewise polynomial approximation at the interior knots (between adjacent elements) is achieved by imposing the following constraints K
z[l0] )
zl-1,jΨj(τ)1) ∑ j)0
l ) 2, 3, ..., NE
(4)
3hTd ) -h
where H is the the Hessian (the symmetric matrix of second derivatives) of the Lagrange function and d is the step direction vector. However, evaluation of the Hessian is often computationally intensive. Alternatively, eq 7 can be solved in an iterative fashion by using a local quadratic approximation of the Lagrange subject to linearized constraints. The quadratic is minimized at each iteration to find the appropriate search direction for the next step.
4.2. NLP Problem Formulation. Define
y)
z
1 T T min 3Φ (yi)d + d Bid 2 d
θ
h(yi) + 3hT(yi)d ) 0
[]
which is an augmented vector of polynomial coefficients and parameters. The NLP problem is r
min Φ(y) ) y
eTµ Weµ ∑ µ)1
yL e y e yU
Bi+1 ) Bi +
r
∑ eTµ eµ µ)1
h(y) ) 0
(6)
4.3. SQP Strategy. One way to solve the NLP problem is to use the efficient optimization strategy SQP, which consists of three major steps. The first step involves the approximation of the Hessian of the Lagrangian function, the second solves a quadratic programming (QP) subproblem, and the third finds the step length in the search direction obtained in step two. The first-order necessary conditions for finding the optimum of a constrained NLP problem are known as the Kuhn-Tucker conditions for optimality. These require, at the solution, that both the equality constraints hold and the Lagrangian function
L ) Φ(y) + λTh(y)
(7)
is at an extremum. That is,
3L ) 3Φ(y) + λT3h(y) ) 0 h(y) ) 0
(8)
where λ is the vector of Lagrange multipliers. Application of Newton’s method, which minimizes a secondorder Taylor series expansion, to solve eq 8 yields
Hd + 3hλ ) -3Φ
qiqTi qTi qi
BipipTi Bi -
(11)
pTi Bipi
(5)
where eµ t zdata - zµ(y) represents the vector of errors µ between the data, zdata, and the model predictions, z(y) at the µth experimental measurement, h(y) is the set of equality constraints that includes the discretized model, continuity constraints, and boundary conditions at the end points, and (yU) and (yL) are upper and lower constraints on the parameters. In this work, W is selected to be the identity matrix and no bounds are placed on the parameters. This is expressed by the reduced NLP problem
min Φ(y) )
(10)
where B is an approximation of the Hessian. The Hessian is approximated with a quasi-Newton update formula, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method
h(y) ) 0
y
(9)
where
pi ) yi+1 - yi and m
qi ) 3Φ(yi+1) +
λj3hj(yi+1) - 3Φ(yi) + ∑ j)1 m
λj3hj(yi) ∑ j)1 because these methods guarantee superlinear convergence as they accumulate second-order information. The solution obtained from the QP subproblem is the search direction for the next iteration,
yi+1 ) yi + sidi
(12)
where s is the step length parameter. However, the length of this step may be beyond the valid range of the quadratic approximation leading to an insufficient improvement in the objective function or the constraints. Therefore, it is required that the step used at the next iteration reduces the value of a merit function. Let this merit function be the augmented Lagrange function
1 Φ(y) + λTh(y) + c||h(y)||2 2
(13)
that contains the objective function, the constraints, and a term that drives the Lagrange multipliers to their true values. c is a penalty parameter and ||‚|| is the Euclidean norm. The merit function is minimized for s in the direction obtained from the QP subproblem by a line search procedure. Common line search procedures include polynomial interpolations or curve fitting. 5. Results The PEM model by Nguyen and White is first normalized and then discretized by orthogonal collocation (Nguyen and White, 1993).
4268 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997
Figure 4. Location of collocation knots based on the profile of the molar flow rate of water vapor in the anode.
5.1. Orthogonal Collocation on Finite Elements. The collocation procedure is dependent upon the values of the collocation parameters, namely the number of elements, the location of the collocation points, and the other of the polynomials. These parameters affect the total number of variables, n ) m + p, that must be optimized, subject to m ) (NS)(NE)(K + 1) equality constraints and p parameters. To make the problem tractable, it is prudent to limit the polynomial order to four or less with three or four elements. For example, in the PEM model with four unknown parameters, n ) 184 when the polynomial order and the number of elements are each four. The locations of the knots, boundary points between finite elements, are made on the basis of the experimental state profiles. Shown in Figure 4 is a typical profile of v the molar flow rate of water vapor at the anode, M ˙ w,a . It is advantageous to use more elements in the steep gradient region than where the trajectory is flat. For this model, the knot points are placed at 10 and 30% of the channel length for three elements and an additional knot is placed at 5% for four elements. The placement of the knot points was not optimized. The accuracy of the orthogonal collocation method for different combinations of elements and polynomial orders is determined at the selected knot locations by solving the system of equations with known parameters and comparing the solution profiles with those obtained using the commercial differential and algebraic solver DDASAC (Double precision Differential-Algebraic Sensitivity Analysis Code) (Caracotsios and Stewart, 1985). The collocation solution is obtained using the fsolve routine of MATLAB. This routine finds roots of nonlinear equations by formulating the problem in a leastsquares structure and solving with the Gauss-Newton method. The DDASAC solution is taken to be the reference and is assumed to be sufficiently close to the actual solution. Comparisons are made by evaluating the solutions at 201 equally spaced points. The error, or difference, between the solution profiles obtained using collocation and DDASAC at each of these points is squared and then summed to obtain the sum squared error. The results of this comparison are shown in Figure 5 where second-, third-, and fourth-order polynomials are
Figure 5. Sum squared error comparison of DDASAC and collocation solutions.
Figure 6. Comparison of the effect of collocation parameters on the % error for the molar flow rate of water vapor in the cathode from 0 to 10% of the channel length.
each evaluated for three and four elements, respectively. From this figure, we see that the accuracy of the collocation solution improves at the expense of increased computational time as more elements and/or higher order polynomials are used. The sum squared errors are reduced by greater than 90% as the order of the polynomial increases and by greater than 65% as the number of elements increases. This result is also representative of the individual solution profiles. Most notable is the molar flow rate of water vapor at the cathode, which constitutes the largest error on an order of magnitude equal to the total error. Since the error of the collocation solution with respect to the reference solution is dominated by the error of the solution for the molar flow rate of water vapor at the cathode, this is the state that has been selected for closer inspection. Figures 6 and 7 provide a point by point perspective into the effects of our three collocation parameters: the number of elements (NE), the order of the polynomials (K), and the knot placements. The curves designate the percentage error of the collocation solutions obtained for combinations of three and four elements with third- and fourth-order polynomials with respect to the reference solution (DDASAC). Figure 6 shows the percentage error for the first element when NE ) 3 and the first and second elements
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4269 Table 3. Estimation of Single Parameters with Fixed Polynomial Coefficients elements
order
10-2I°, A cm-2
10-3U, J s-1 cm-2 °C-1
10-7D°, cm2 s-1
kc, s-1
3 4 3 4 3 4
2 2 3 3 4 4
1.0010 1.0013 1.0027 1.0014 1.0010 1.0022 1.0015
2.5000 2.5033 2.5035 2.5008 2.5005 2.5000 2.5003
5.5000 6.1167 5.9143 5.5045 5.5340 5.6083 5.5623
1.0000 0.9995 0.9993 0.9997 0.9998 1.0002 1.0000
Table 4. Parameter Estimation Results for Ideal Data 10-2I°, A cm-2
10-3U, J s-1 cm-2 °C-1
10-7D°, cm2 s-1
kc, s-1
1.0010 1.0010
2.5000
5.5000
1.0000
2.4938 5.5220 Figure 7. Comparison of the effect of collocation parameters on the % error for the molar flow rate of water vapor in the cathode from 10 to 100% of the channel length.
0.9979 1.0018 error 1.0027
max. error max. % error ∑error2
max. error max. % error ∑error2
M ˙ lw,a
M ˙ vw,a
Ta
10-4 0.11 10-6
10-14 0 10-28
10-3 1.2 10-5
10-12 0 10-27
M ˙ O2
l M ˙ w,c
v M ˙ w,c
Tc
10-4
10-5
10-3
0.05 10-6
78 10-7
0.98 10-4
10-3 0.14 10-5
when NE ) 4 with a knot placed at 5%. Notice that an additional element in the region of greatest error reduces the error more than increasing the polynomial order. Recall that this region of greatest error, near the entrance of the flow channel, corresponds to the steep region of the state profiles. However, the reduction in error in the remaining two elements, as seen in Figure 7, is significant enough to result in the smaller sum squared error of 3 × 10-5 (NE ) 3, K ) 4) as compared to 1 × 10-4 (NE ) 4, K ) 3). One may also note that the collocation points correspond to the stationary points of the percentage error curves and that there is almost no error at the elemental boundaries. If the aim of this work is to obtain a high degree of accuracy from the collocation solution, one would do well to use four or more elements and fourth-order polynomials; however, this is not the objective. The complete parameter estimation problem requires accuracy to the extent that the error between the collocation and actual solution is inconsequential in comparison to the error between the model and the data. Thus, for the same model parameters, the sum squared error for the data with respect to the actual and discretized model should be equivalent to significant digits. For different combinations of polynomial order and number of finite elements, it is found that an acceptable performance is obtained for at least three finite elements and third-order polynomials. A greater number of collocation points marginally improves the collocation solution, but the complete parameter estimation solution would not be affected appreciably. Table 2 summarizes the maximum errors, maximum percentage errors, and sum squared errors for eight states at the 201 equally spaced points using these collocation parameters (NE ) 3, K ) 3). The maximum percent error
error 0.9981 2.4951 2.4954
Table 2. Order of Error Obtained for Three Elements and Third-Order Polynomials M ˙ H2
2.4927
0.9966 1.0004 error 0.9967 0.44%
2.4943 2.4952 2.4970 2.4967 0.29%
5.5271 5.5741 5.5340 error 5.5281 5.5340 1.30%
0.9988 1.0001 0.9986 error 0.9994 0.9989 0.21%
of 78 for the molar flow rate of liquid water at the cathode is somewhat misleading. This error occurs at the entrance to the flow channel where the molar flow rate is extremely small, thus even a small discrepancy translates into a seemingly large percentage error. Additionally, an analysis where the collocation coefficients are fixed and at least one parameter is unknown is carried out to evaluate the choice of collocation parameters on the parameter estimation problem. Table 3 shows the results obtained for a single unknown parameter. The top line represents the actual parameter values, that is, the value of each parameter used in the model to obtain the simulated data. Observe that the parameters are effectively estimated, but because the coefficients are fixed, there is no apparent benefit to increasing either the number of finite elements or the polynomial order. 5.2. Parameter Estimation. The following results on parameter estimation are obtained using the collocation parameters (NE ) 3, K ) 3). The simultaneous solution and parameter estimation problem require an initial guess to the unknown parameters, which include the polynomial coefficients. Since these coefficients are the values of the states at the collocation points, they may be estimated from experimental data. This is achieved by interpolation from a fourth-order fit of the data. Likewise, a good initial guess of the parameters is determined by solving an unconstrained least-squares problem similar to what was done to generate Table 3. Experimental data on the flow rates are difficult to obtain. In practice, the data are of the inlet-outlet type or, more often, voltage and current. Although the parameter estimation technique may be applied to either of these situations, the ideal scenario of flow rate and temperature data along the channels is more instructive to demonstrate the solution framework. As such, the data used here are of this type, simulated using DDASAC, at two points in each element with one additional end point.
4270 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997
Figure 8. Ideal state profiles along the channel length, x (cm), for the flow rate of hydrogen, oxygen, liquid water, water vapor, and temperature.
Figure 9. Nonideal state profiles along the channel length, x (cm), for the flow rate of hydrogen, oxygen, liquid water, water vapor, and temperature.
The objective of the parameter estimation task is to determine the parameters of a model so that the model’s results best match the experimental data. Here, the minimum least-squares criterion is used as the merit for this objective. For, the case where simulated data are used, the actual parameter values are known; thus, one has the further benefit of being able to compare the estimates to the actual parameter values. The SQP strategy is used to solve the simultaneous solution and parameter estimation problem under ideal conditions. This is done when one or more of the four adjustable parameters are assumed to be unknown. The parameter estimation results are shown in Table 4, where the top row contains the actual values and the
bottom row shows the maximum percentage error for each parameter. Table 4 contains results for every possible combination of the parameters from one unknown to four unknowns. Observe that the largest maximum error obtained is a mere 1.3%. Figure 8 v shows the solution profiles of M ˙ H2, M ˙ w,a , M ˙ O2, M ˙ lw,c, v M ˙ w,c, and Tc when all four adjustable parameters are unknown. For comparison, the solution profiles using the true parameter values are also shown. Having demonstrated the capabilities of the solution procedure for idealized simulated data, it is natural to investigate the more realistic case of nonideal experimental data. In order to simulate an expected level of
Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 4271 Table 5. Parameter Estimation Results for Nonideal Data
actual value ideal value av value av % abs error max % abs error
10-2I°, A cm-2
10-3U, J s-1 cm-2 °C-1
10-7D°, cm2 s-1
kc, s-1
1.0010 0.9967 1.0299 11.39 29.71
2.5000 2.4967 2.4793 6.50 15.68
5.5000 5.5340 5.3462 6.93 22.67
1.0000 0.9989 0.9988 3.64 9.93
experimental error, Gaussian noise with a mean of zero and a variance of 0.01 is added to the data. Although this may not represent the possibility of process-model mismatch, the simulated nonideal data do present the challenge of convergence to a nonzero solution. Figure 9 shows the solution profiles for the same states as in Figure 8 when all four adjustable parameters are treated as unknowns. Also shown are the ideal profiles and the noisy data (×). Because the signal to noise ratio (0.143) is 1 order of magnitude smaller for Tc as compared to M ˙ H2 (1.43), one would expect larger discrepancies between the data points and the actual profiles of Tc. The solution profiles obtained using the SQP strategy are a closer fit of the noisy data than the ideal profile because the objective function minimizes the error between the model’s predictions and the noisy data. Thus, for the given objective, the solution parameters are a better solution than the parameters used to obtain the data. In the limit of infinite data subject to only Gaussian noise, the nonideal solution will approach the ideal one. Table 5 shows the parameter estimation results obtained when all four adjustable parameters are treated as unknowns. The results are for 25 sets of nonideal data all having a variance of 0.01. The first row represents the actual parameter values for the ideal data; for comparison, the second row contains the parameters estimated using ideal data, and the third provides the average results obtained for nonideal data. The fourth row represents the average absolute error, which is the average of the absolute value of the percentage error between the actual and the predicted parameter value for each data set. The fifth row shows the worst case scenario of the predicted parameters for the 25 sets of data. In the worst case, the error is no larger than 30%. However, this solution is still the best attainable as defined by the objective function. What can be concluded is the quality of the model’s prediction and the accuracy of the parameters estimated are only as good as the data supplied. 6. Summary and Conclusions This work has demonstrated parameter estimation by means of a nonlinear programming (NLP) formulation. The solution was obtained using sequential quadratic programming (SQP). The approach was demonstrated on a proton-exchange-membrane fuel cell model described by differential and algebraic equations. Nine states and four parameters were estimated, subject to satisfying an objective function that minimized the difference between known data and the model’s predictions. Both simulated ideal data and data corrupted with Gaussian random noise were used to validate the solution technique. Under ideal conditions, the worst estimate had no more than 1.3% error; in the case of nonideal conditions, this error may be as large as 30%. Preliminary results obtained when experimental voltage-current data are used corroborate these findings
(Suares and Kosanovich, 1996). On the basis of this and similar results on simpler systems found in the open literature (Tjoa and Biegler, 1991), it is concluded that this framework holds a great deal of promise for the simultaneous solution of both the states and parameters for systems of these types. Acknowledgment We acknowledge financial support from the U.S. Department of Energy, Cooperative Agreement No. DEFG02-91ER75666. List of Symbols a ) heat transfer area per unit length, cm B ) approximation of the Hessian c ) penalty parameter Cp ) heat capacity of gas, (J/mol)/°C d ) step direction vector D° ) a parameter used in the expression for the diffusion coefficient of water, cm2/s e ) vector of errors F ) Faraday constant, 96 487 C/equiv F ) algebraic right-hand side of the ODE h ) channel height, cm h(y) ) set of equality constraints H ) enthalpy of water, J/mol H ) Hessian of the Lagrange function I(x) ) local current density, A/cm2 I° ) exchange current density for the oxygen reaction, A/cm2 K ) order of the polynomial kc ) evaporation and condensation rate constant, s-1 L ) Lagrangian function M ˙ i ) molar flow rate of species i, mol/s m ) number of equality constraints n ) number of variables to be optimized NE ) number of finite elements NS ) number of states p ) number of parameters P ) cell total pressure, atm Psat ) vapor pressure, atm R ) gas constant, 82.06 cm3 atm mol-1 K-1 R ) residual equation s ) step length t ) thickness, cm T ) temperature, °C U ) overall heat-transfer coefficient, J s-1 cm-2 °C-1 V ) voltage, V w ) channel width, cm W ) weighting matrix among the state variables x ) direction along the channel length, cm y ) augmented vector of states and parameters z ) polynomial coefficient z(x) ) state approximation Greek Symbols R ) net water flux per proton flux λ ) vector of Lagrange multipliers νi ) location of the ith finite element Φ ) objective function Ψ ) Lagrange polynomial Ψ ˙ ) derivative of the Lagrange polynomial with respect to τ σ ) conductivity, 1/(Ω/cm) τ ) collocation point θ ) vector of parameters Subscripts and Superscripts a ) anode c ) cathode
4272 Ind. Eng. Chem. Res., Vol. 36, No. 10, 1997 H2 ) hydrogen l ) liquid L ) lower bound m ) membrane N2 ) nitrogen O2 ) oxygen s ) solid phase T ) transpose U ) upper bound v ) vapor w ) water
Literature Cited Amphlett, J. C.; Baumert, R. M.; Mann, R. F.; Peppley, B. A.; Roberge, P. R.; Harris, T. J. Performance modeling of the ballard mark iv solid polymer electrolyte fuel cell, i. mechanistic model development. J. Electrochem. Soc. 1995a, 142 (1), 1-8. Amphlett, J. C.; Baumert, R. M.; Mann, R. F.; Peppley, B. A.; Roberge, P. R.; Harris, T. J. Performance modeling of the ballard mark iv solid polymer electrolyte fuel cell, ii. empirical model development. J. Electrochem. Soc. 1995b, 142 (1), 9-15. Bernardi, D. M.; Verbrugge, M. W. Mathematical model of a gas diffusion electrode bonded to polymer electrolyte. AIChE J. 1991, 37 (8), 1151-1163. Bernardi, D. M.; Verbrugge, M. W. A mathematical model of the solid-polymer-electrolyte fuel cell. J. Electrochem. Soc. 1992, 139 (9), 2477-2491. Caracotsios, M.; Stewart, W. E. Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Comput. Chem. Eng. 1985, 9, 359-365. Dennis, J. E.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; Prentice-Hall, Inc.: Englewoods Cliffs, NJ, 1983.
Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw-Hill: New York, 1988. Fuller, T. F.; Newman, J. Water and thermal management in solidpolymer-electrolyte fuel cells. J. Electrochem. Soc. 1993, 140 (5), 1218-1225. Luenberger, D. G. Linear and Nonlinear Programming, 2nd ed.; Addison-Wesley: Reading, MA, 1989. Nguyen, T. V.; White, R. E. A water and heat management model for proton-exchange-membrane fuel cells. J. Electrochem. Soc. 1993, 140 (8), 2178-2186. Springer, T. E.; Zawodzinski, T. A.; Gottesfeld, S. Polymer electrolyte fuel cell model. J. Electrochem. Soc. 1991, 138, 2334-2342. Springer, T. E.; Wilson, M. S.; Gottesfeld, S. Modeling and experimental diagnostics in polymer electrolyte fuel cells. J. Electrochem. Soc. 1993, 140 (12), 3513-3526. Suares, G. E.; Kosanovich, K. A. Parameter estimation of a protonexchange membrane fuel cell using an NLP approach. Submitted for publication in J. Electrochem. Soc., 1996. Tjoa, I.-B.; Biegler, L. T. Simultaneous solution and optimization strategies for parameter estimation of differential-algebraic equation systems. Ind. Eng. Chem. Res. 1991, 30 (2), 376-385.
Received for review March 24, 1997 Revised manuscript received July 14, 1997 Accepted July 15, 1997X IE970251P
X Abstract published in Advance ACS Abstracts, September 15, 1997.
