Singularly Perturbed Modeling and LQR Controller Design for a Fuel Cell System

: Modeling and control of proton-exchange membrane fuel cells (PEMFC) has become a very popular research topic lately due to the increasing use of renewable energy. Despite this fact, most of the work in the current literature only studies standard dynamical models without taking into consideration possible parasitics such as small gas ﬂow perturbations that could be available in the system. This paper addresses this issue by elaborating on time-scale modeling of an augmented eighteenth-order PEMFC-reformer system via singular perturbation theory. The latter captures time scales that arise in the model due to the presence of small perturbations. Speciﬁcally, a novel and efﬁcient algorithm that helps identify the presence of different time-scales is developed. In addition, the method converts an implicit singularly perturbed model into an explicit equivalent where the time-scales are evident. Using this algorithm, a complete singularly perturbed dynamic model of the augmented eighteenth-order PEMFC-reformer system is obtained. Modeling of the PEMFC-reformer system is followed by linear quadratic regulator (LQR) design for the individual time-scales present in the system.


Introduction
An increase in worldwide energy demand and the quest for sustainability have led to major investments in renewable resources. The latter are eco-friendly, generally reliable, and have lower operating costs. While there are several disadvantages associated with renewable sources such as vulnerability or the inability to generate power in large quantities, the consensus among researchers is that energy in the future will primarily be generated by renewable sources. One of the latter is the fuel cell. This source has become very popular due to a higher efficiency compared to other renewable sources, low maintenance costs, and ubiquity of hydrogen. Fuel cells have been widely used in industrial settings (e.g., factories) [1], residential units [2], and the automotive industry (e.g., electric vehicles) [3,4], and are also commonplace in smart grid ecosystems where they are used as supplementary renewable sources in microgrids [5,6]. A major enabler that has aided in hydrogen mobility, and subsequently in the widespread use of fuel cells, has been the investment in infrastructure by several countries in North America, Europe, and Asia [7,8]. Government statistics in these countries [8] show that demand for a hydrogen economy is increasing, hence possibly putting the fuel cell at the forefront of the clean energy revolution.
The device works by converting the chemical energy of hydrogen and oxygen into electricity through a pair of redox reactions [9]. In addition to electricity, the chemical reaction produces water and heat [9,10]. While the method of operation is standard for all varieties of fuel cells, they have shown in Figure 1. The PEMFC contains an electrolyte (center unit in Figure 1), an anode, and a cathode. Hydrogen and oxygen gas are the primary fuels used in the redox chemical reaction. The products of the latter are energy, excess fuel (if any), water, and heat. Due to several advantages, such as no harmful gas emissions, light weight, small size, and low operating temperature, the PEMFC has become a very popular energy alternative.
The chemical reaction starts at the anode, where hydrogen fuel is processed and electrons are separated from protons on the surface of a catalyst. The electrons travel via a circuit to generate electricity while the protons pass through the membrane to the cathode of the cell. At the cathode, existing protons and electrons coming from the circuit are combined with oxygen to produce water, which together with heat are the only waste products. It is important to note that instead of oxygen gas in purified form, air can also be used as a fuel. An additional process at the electrode extracts oxygen from air which is then used in the chemical reaction. The hydrogen gas going into the fuel cell on the other hand has to be in pure form. It can either be supplied from a dedicated tank or obtained from a fuel processing system (FPS) also known as a gas reformer. The FPS typically uses natural gas as raw material and produces hydrogen via a series of chemical reactions.
This research investigates the modeling and controller design of a PEMFC and FPS originally studied in [9]. Unlike in [9], this work explores the singularly perturbed modeling of the PEMFC-FPS augmented system. On that end, an overview of singular perturbation methods in control is provided in the following subsection to familiarize the reader with the theory.

Overview of Singularly Perturbed Methods in Control
While singular perturbation methods in control have been used extensively since the sixties [15], they are still an important research topic [16][17][18]. For a detailed treatment of singular perturbation methods, the reader is referred to [15,[19][20][21] and the references therein.
A linear time-invariant (LTI) strictly proper singularly perturbed system in mathematical form is represented asẋ where x 1 (t) ∈ R n and x 2 (t) ∈ R m are the slow and fast state variables respectively, u(t) ∈ R p are the system control inputs, y(t) ∈ R q are the system measurements, and ε is the small singular perturbation parameter 0 < ε 1. In a real physical system, ε can represent parasitic capacitance in a circuit, moment of inertia in a mechanical system, etc. All the matrices in (1) are constant and of appropriate dimensions. Two time-scale LTI singularly perturbed systems have eigenvalues located in two disjoint groups: slow O(1) eigenvalues close to the imaginary axis and fast O( 1 ε ) eigenvalues far from it. The following standard assumption is imposed [19].
When ε = 0 (a common strategy for order reduction), the following reduced-order system corresponding to the slow dynamics is obtained:x where The approximated fast subsystem is [19] where τ = t/ε. According to the theory of singular perturbations [15,19], the approximation obtained in (2)-(4) satisfies Hence, the smaller ε, the better the approximation. Furthermore, λ(A 0 ) and λ(A 4 /ε) are O(ε) perturbation from the slow and fast eigenvalues of the original model, respectively.
Another effective method to obtain exact dynamic decoupling is by utilizing the transformation developed in [22]. Referred to as the Chang transformation, it diagonalizes the system by exposing the slow and fast dynamics. The transformation is given as follows.
After (6a) is applied to (1), the decoupled system is theṅ where A noticeable difference between the reduced-order slow model defined in (2)-(3) and the decoupled system (7) is that the measurement in (7) lacks an input u(t) whileȳ(t) from (2) includes a D 0 u(t) term.
Matrices L and H are obtained by solving the following equations.
The reader can refer to [23][24][25] for methods on finding the solution of L and H.

Time-Scale Decoupling via the Ordered Schur Decomposition
This section serves as the basis for the singularly perturbed modeling of the PEMFC-FPS augmented system. Namely, the algorithm that achieves complete decoupling of an arbitrary singularly perturbed model is developed. This method is then applied in the following section. Theoretical techniques presented here were briefly introduced in [26]. This section contains a thorough explanation of the method.

Ordered Schur Transformation
As noted earlier, many real physical systems contain small parasitics when they are modeled. This forces parts of the system to operate in different time-scales. Quite often, it may be difficult to distinguish between the time-scales. Methods such as permutation matrices or other similarity transformations [19] have been proven successful to obtain a standard singularly perturbed form for two time-scale systems but it is challenging when additional time-scale dynamics are present.
In this paper, the aforementioned issue is addressed by developing a method that brings an implicit singularly perturbed system into its explicit form, where the perturbation parameters ε 1 , ε 2 , · · · , ε N are either known or can be easily determined. ε 1 is associated with the slowest state variable and ε N is associated with the fastest, namely, 0 < ε N ε N−1 . . . ε 1 = 1. Consider a general implicit multiple time-scale system without inputs, as shown in (10). where T is the state vector. To simplify the problem, an ordered Schur decomposition is employed to transform the model into a well-conditioned form. This is followed by the extraction of the perturbation parameter and sequential decoupling to obtain the individual time scales. The Schur decomposition is an efficient method used to find the system's eigenvalues by utilizing the QR algorithm [27]. For a matrix A ∈ R n×n , there exits a unitary matrix T ∈ R n×n such that T T AT = A is upper quasi-triangular.
The eigenvalues appearing along the diagonal of A can be arbitrarily ordered. An additional transformation has to be employed to achieve desired reordering (descending) of the system matrix [27][28][29]. A transformation z(t) = Tx(t) [30] can be found such that the unitary matrix T decomposes the system into the Schur form. Upon applying the ordered Schur algorithm, the dynamic Equation (10) takes the following formż (t) = Az(t), where Diagonal block matrices A ij , i = j, contain the system's eigenvalues in descending order. z(t) 1 , z(t) 2 , . . . , z(t) N are vectors each representing each time scale and Ψ ij are matrix blocks of appropriate dimensions. Note that blocks A ij , i = j, represent individual time-scales rather than individual eigenvalues.

Parameter Extraction and Time-Scale Decoupling
Prior to decoupling the transformed system, it is essential to convert it to an explicit singularly perturbed form by extracting the perturbation parameters from the system matrix. This is achieved by defining the perturbation parameters. For two time-scale stable systems with clearly separated eigenvalues (real or complex with small imaginary parts), ε is commonly evaluated as ε = max Re{λ s }/min Re{λ f } [31]. However, referring to examples such as the PEMFC-FPS under consideration, it can be inferred that the imaginary parts of the eigenvalues are indispensable for calculating the perturbation parameter. Hence, the latter is evaluated as the ratio of the magnitudes of the largest eigenvalue of the slowest time-scale with the smallest eigenvalue of the next fastest time-scale Since the system is in ordered Schur form, the perturbation parameters can be easily evaluated using (14) and extracted to put the system into the explicit singularly perturbed form. The explicit multi-time-scale system now looks as follows: where the elements of the system matrix have been scaled in accordance with parameter extraction.
The explicit system can now be decoupled into multiple distinct time-scale systems by successively applying the Chang transformation in (6a).
To initiate the decoupling, the perturbation parameter is extracted from the fastest time-scale and (12) is rewritten as a standard two time-scale singularly perturbed form.
In (16), matrices A 3 and A 4 represent the last row of the system's matrix in (13) with ε N extracted. A 4 represents the fastest time-scale, while A 3 contains the rest of the matrix blocks which happen to be all zero in this case. A 1 and A 2 are matrices of appropriate dimensions containing the rest of the system matrix in (12).
Utilizing (6a), the system in (16) is initially decoupled into two subsystems, where the fast subsystem represents the fastest time-scale available and the slow subsystem contains the rest of the time-scales.ξ As a reminder to the reader, matrices L and H satisfy (9). A modified Newton's method developed in [25] follows. D The new slow subsystem (17a) is partitioned again, as in (19), where ε N−1 is now extracted from the second fastest time-scale.
The algorithm is applied sequentially till all the perturbation parameters have been extracted and the system is in explicit singularly perturbed form. The relation between the system in original coordinated z(·) and the new one is given as where T is defined as T =T 2T3 · · ·T N−1 T N andT i is Matrix T i represents the linear transformation for each time-scale i. On the other hand, matrix T i is an augmentation of T i with an identity matrix I of appropriate dimensions so that its size is the same as T. Unlike in [31], the system matrix in ordered Schur form simplifies the computations. For a quasi-triangular system such as (12), A 3 in (16) is 0. Then, equations for the solution of matrices L and H in (9) simplify to the following.
An additional simplification comes due to the new system matrix structure.
Theorem 1. Due to the structure of (15), matrix L N of the Chang transformation evaluates to a zero matrix.
Proof. Upon applying the recursive algorithm (18), it is easy to show that matrix L N in (22a) evaluates to a zero matrix by solving the Sylvester equation for the first iteration

Mvec L
(1) and is full rank. Therefore, ker(M) = 0, which implies L In the absence of matrix L, (22b) then becomes just a Sylvester equation.
After the process is repeated for all the N time-scales available in the system, the individual subsystems are then given as follows.
System (26) is now completely decoupled into a standard explicit singularly perturbed form. Note that if control is considered, the input matrix for each time-scale would be obtained sequentially in a similar fashion. The slow and fast input matrices for a two time-scale system are evaluated as The overall process discussed in this section can be summarized as follows.
Algorithm 1 Time-Scale Decoupling of Implicit Singularly Perturbed Systems 1: Input: Implicit singularly perturbed system 2: Apply Schur decomposition 3: If Eigenvalues are not ordered then apply swapping algorithm [27] or [29] 4: Evaluate ε i for each time-scale and form explicit system 5: Use Chang transformation [22] to get individual time-scales 6: Output: Completely decoupled system

Singularly Perturbed Modeling of the PEMFC-FPS
The general nonlinear models of the PEMFC and FPS systems studied in this paper were originally derived and used for controller design in [9]. The fuel cell considered in [9] has been used for research in the automotive industry and has a stack size of approximately 40 kW. Its operation temperature varies from 50 • C to 100 • C. The associated FPS is typically used for PEMFCs with stack size of 100 kW and is made up of a hydrodesulfurizer, a catalytic partial oxidation reactor, a water gas shift reactor, and a preferential oxidation reactor [9]. Original PEMFC and FPS models have been linearized around a nominal operation point, where the system net power is z o 1 = 40 kW and oxygen excess ratio is z o 2 = 2 [9]. Both linear and nonlinear model responses are compared in [9] and it is validated that the errors are insignificant.
This paper does not elaborate any further on the nonlinear systems but instead uses the linearized system for modeling and controller design. The reader is referred to [9,10] for additional details.

PEMFC and FPS Linear Model
Variables that are used in determining the model of the PEMFC fall into two main categories: masses of the gases that are used in the chemical reaction and the pressure of these gases. The state vector for the linear PEMFC model is defined as δx ∈ R 8 , δx = m O 2 m H 2 m N 2 ω cp p sm m sm m w,an p rm . Corresponding descriptions of state variables is provided in Table 1.   Using information from the PEMFC and FPS models, an augmented linear state-space model is created as in (27) [14].
State-space system in (27) can be written as follows: where A aug ∈ R 18×18 , B aug ∈ R 18×2 , C aug ∈ R 5×18 , and D aug represent matrices in (27). Numerical values of A PEMFC , A FPS , B PEMFC , B FPS , C PEMFC , and C FPS are available in [9] and also included in Appendix A for reference. Note that the feedthrough matrix D is zero in both the PEMFC and FPS, therefore D aug = 0. In the following subsection, a singularly perturbed linear model of (27) is derived using methods from Section 3.

Singularly Perturbed Model of the PEMFC-FPS
The eigenvalues of the augmented system are key in determining the time-scale composition of the model. Figure 2 depicts graphically the sorted eigenvalues of the PEMFC-FPS. Recall from Section 3.2 that the magnitude of the eigenvalues is used to evaluate the singular perturbation parameter ε i . Based on the distribution shown in Figure 2, a three time-scale model is chosen. The first nine modes represent the slow subsystem, the following six represent the second time-scale, and the last three constitute the fastest time-scale. Note that this decision is not necessarily strict. At this point, the singular perturbation parameters of the three time-scale model can be calculated using (14). The magnitude of the largest eigenvalue of the slow subsystem is |λ 9 | = 3.33. The magnitude of the smallest eigenvalues of the second and third time-scales are |λ 10 | = 12.17 and |λ 16 | = 157.90. Hence, ε 1 , ε 2 , and ε 3 are as follows.
The singularly perturbed model then becomes   ż State variables belonging to each time-scale are presented in Table 3. After applying the algorithms developed in Section 3 to the original augmented PEMFC-FPS model in (27), the matrices that help form the singularly perturbed model (30) are obtained. The state matrix A schur = T −1 A aug T, control matrix B schur = TB aug , and output matrix C schur = C aug T −1 are available in Appendix B. Linear mapping T represents the transformation from the original augmented model (27) to the ordered Schur model. It is also important to note that matrices A schur and B schur , specified in Appendix B, do not represent the state and control matrices in (30) because the singular perturbation parameters have not been extracted yet. To obtain model (27), A schur and B schur are simply multiplied by the ε of the corresponding time-scale.
Lastly, the Chang transformation [22] is applied sequentially to obtain a completely decoupled system (refer to Algorithm 1).

LQ Control System Design for the PEMFC-Reformer System
In this section, LQ controller design is demonstrated for the decoupled subsystems of the PEMFC-FPS model obtained in the previous section. Specifically, all three time-scales are considered independent models and LQ control is applied to each of them. An overview of LQ control is first presented followed by simulation results.

LQ Control Overview
In LQ control, our objective is to drive the state vector x(t) to the origin from any nonzero initial state [32]. If only state-feedback control is used and the closed-loop eigenvalues are assigned far inside the left half of the s-plane, the state x(t) will die out quickly but elements of the feedback gain can be high in magnitude, therefore requiring a high control cost. Alternatively, if the closed-loop eigenvalues are assigned closer to the imaginary axis, the rate of decay of x(t) will be slow and a small control action will be required. This trade-off between the state-vector x(t) and control action u(t) can be represented by the following cost function.
The LQ problem is then to find a control u(t) such that the cost function (31) is minimized. An important assumption here is that R > 0 and Q ≥ 0. The solution to this infinite-horizon optimal control problem is a control input given by [32] u(t) = −Fx(t), (32) where the gain is given as F = R −1 B T P and matrix P is the solution of the following algebraic Riccati equation [32].
The state vectors x(t) for each time-scale are given in Table 3.

Feedback Controller Design
Controller design starts with the decoupled PEMFC-FPS model presented in (27). Initially, the ordered Schur transformation T schur is applied to obtain model (30) followed by sequential Chang transformation (6a) to get a decoupled three time-scale model as follows.
Note that it is important to convert the control input vector into the original coordinates of the system [33]. To do so, the transformation η i (t) = T Chang z i (t) is used to put the system in z(·) coordinates, where z(·) represents the ordered Schur model. Another similarity transformation to get to the original coordinates x(·) is necessary, namely, z i (t) = T ordSchur x i (t). The control input them becomes u(x(t)) = −FT Chang T ordSchur x(t).
In the following subsection, the aforementioned principles are applied to the decoupled PEMFC-FPS model. It is worth noting that pairs (A, B) and (A, C) for all the state-space models considered are stabilizable and detectable, respectively.

LQ Control Simulation Results
LQ controller design in this section is based on requirements of an actual physical system as noted in [10]. Table 4 shows operating time specifications for different processes that occur in a PEFMC-reformer system. In this section, LQ controllers for one state in each time-scale are designed. Scaled reference step inputs are considered in each of the three cases. The selection of matrices Q and R in this section is based on parametric studies. This method ensures that (31) is minimized and requirements in Table 4 are met. An example of LQ control for various values of Q and R is shown in Figure 3 for m O 2 . It is evident that different values of Q (shown as the trace of Q in this case) and R can have a significant impact on the controller's response.  Recall that the first subsystem corresponding to the slowest time-scale is a ninth-order model. For this case, only the mass of oxygen, namely m O 2 is selected for control. From Table 4 (see [10]), it is observed that this process has a time constant of O(10 0 ) seconds. Ideally, for a satisfactory performance, it is preferred that the response settles in less than ten seconds. Likewise, a relatively small overshoot would be ideal so that the mass of oxygen in the chamber does damage the hardware. Figure 4 shows the response for a reference input set to 1 kg. The simulation is run for fifteen seconds. It can be observed that the settling time is around eight seconds (within the requirements). There is a noticeable overshoot but given the time duration, it does not pose any risks to the hardware. The values of the matrices Q and R used for this design are presented in (37).   Figure 5 represents the second fastest subsystem. Specifically, the state variable corresponding to the pressure in the supply manifold has been used to demonstrate controller design. A reference input of 7 kPa has been used in the simulation. According to the requirements, the operational time for processes occurring in the manifold are of O(10 −1 ) seconds. Therefore, the design entails a very short settling time and a small overshoot.
Simulation results of the design are presented in Figure 5. It is clear that this design provides a short rise time, a very minimal settling time which falls within the requirements, and no overshoot.
Lastly, response for the fastest subsystem is shown in Figure 6. In the latter, results for the state corresponding to the mass of nitrogen (m N 2 ) are shown, for a reference input of 0.5 kg. The simulation is run for 0.05 seconds as the reader can observe. It is clear from Figure 6 that the mass of m N 2 goes to the desired setpoint in a very short period of time. In addition, the design ensures a short rise time and no overshoot.

Discussion of Results
Results of this research can be divided into two main parts: 1.
Introduction of an algorithm that converts implicit singularly perturbed systems into explicit ones.

2.
Singularly perturbed modeling of the PEMFC-FPS augmented model followed by LQ controller design.

From Implicit to Explicit Singularly Perturbed Systems
The theoretical basis behind this research lies in the development of an algorithm that converts a linear system with known two or more time-scales (this could be determined from the eigenvalue distribution for example) into a singularly perturbed model. The algorithm relies on an ordered Schur decomposition. Singular perturbation parameters are determined via the ratio of the eigenvalues, and their extraction leads to a final singularly perturbed form. Then, exact decoupling of the obtained singularly perturbed system is sequentially performed.
The developed algorithm is an efficient method that help obtain an explicit form where the time-scales are clearly defined. While there are different methods to put an implicit singularly perturbed system into its explicit equivalent (for example, using permutation matrices), our method offers more flexibility. For instance, if the model is very large scale, using permutation matrices may not be feasible and it could be computationally expensive.

Singularly Perturbed Modeling of the PEMFC-FPS and Controller Design
The main novelty of this research lies in the formulation of the singularly perturbed model of the PEMFC-FPS. Starting with an augmented PEMFC-FPS model, the theoretical basis to obtain a model that is completely decoupled into different time-scales was presented. The enablers behind the modeling are the eigenvalues of the system. Large magnitude differences among the eigenvalues show that the augmented system consists of multiple time-scales, and since singular perturbation theory is commonly used to study such systems, it was applied to (27).
Using this method, the explicit singularly perturbed model (30) was obtained. The latter model was decoupled and controllers were designed for each sub-subsystem. Results showed that a performance within operational requirements can be attained.

Real-World Implications
Modeling and simulation are key tools used for a successful fuel cell system design. In this paper, a modeling methodology was initially developed by using singular perturbation techniques. The method helps obtain an explicit singularly perturbed model and also offer a way to decouple the model into individual time-scales which can be used separately to design controllers. Having a complete model and efficient algorithms is very important not only because it enables the designer to rigorously analyze the system but also because it helps alleviate computational burdens if, for example, parametric studies or Monte Carlo runs are considered.

Conclusions
A lack of PEMFC-FPS time-scale modeling served as motivation for this research. Initially, the theoretical basis that enables the designer convert a linear model with multiple time-scales into an explicit singular perturbation model was introduced. These techniques were applied to a popular augmented PEMFC-FPS model and it was shown that it can be represented as a three time-scale singularly perturbed model. The latter was then decoupled into individual submodels, each corresponding to a different time-scale, and LQ control was applied in each case. The methods developed in this paper capture all the dynamics present in the system including small parasitics. From a practical standpoint, the work presented here is indispensable to the designer since, in addition to a complete model, the efficient algorithm helps decouple the full-order model into individual time-scales for separate analysis such as controller design, as it was shown via case studies.
Author Contributions: Conceptualization, methodology, investigation, writing-original draft preparation, and software, K.K.; conceptualization, validation, formal analysis, investigation, and writing-review and editing, N.Z. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript ( The corresponding control and output matrices are The corresponding control and output matrices are

Appendix B. Singularly Perturbed Model Matrices
Matrices A ordschur , B ordschur , and C ordschur after the ordered Schur transformation.