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

Abstract

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 flow 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. Specifically, a novel and efficient 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.

1. 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 distinct features such as different material composition, efficiency, and operating temperatures. Some of the most common types that are currently in use and also heavily researched are the proton-exchange membrane fuel cell (PEMFC) [9], the solid oxide fuel cell (SOFC) [11], and the phosphoric acid fuel cell (PAFC) [12].

It is worth noting that modeling and simulation play a critical role in fuel cell design. A high-fidelity model gives the designer the capability to modify parameters accordingly and use them in the final implementation, therefore leading to an optimal desired design. One such PEMFC model was originally developed in [9]. This model has been investigated in a variety of domains over the years. For example, in [13], the authors consider dynamic and controlled operation of a PEMFC integrated with a natural gas fuel processor system (FPS) and a catalytic burner (CB). Optimization is performed to generate the air and fuel flow intake setpoints to the fuel processor system for different load levels, and then linear quadratic regulator techniques are used to develop a controller to mitigate hydrogen starvation in the fuel cell. In [9], analysis and design of air flow controllers that prevent oxygen deprivation of the fuel cell stack as the current changes was presented. Improvements in transient oxygen regulation when the fuel processing system voltage was included in the model were demonstrated. A combined PEMFC-FPS model was studied in [14], where the authors used order-reduction techniques to show that similar controller performance can be obtained if the model size is reduced.

Unlike the work of [9] and other aforementioned references, this research differs in that an augmented singularly perturbed model of the PEMFC-FPS is obtained. While the latter has been defined in [14], it has not been studied in a singularly perturbed fashion. The technological context as well as the motivation behind this work lies in the fact that fuel cell systems can contain inherent parasitics such as small impedances or load perturbations. A standard dynamic model does not include such parasitics and that can lead to inaccuracies in simulations which, in turn, would affect the efficiency of the overall system design. Singular perturbation methods are typically used for such models and are very important from a practical aspect since they capture the model’s slow and fast dynamics, hence providing a more accurate model. Besides the singularly perturbed modeling of the fuel cell system, case studies of optimal controller design are presented. Namely, linear quadratic (LQ) controllers are designed for the individual time-scales to meet operational requirements. It is important to note that in this paper, the term “fuel cell system” strictly refers to a setup consisting of a PEMFC and an FPS and does not include additional components such as filters, converters, etc.

The contributions of this paper are threefold. First, the theory that enables one to convert an arbitrary linear system with two or more time-scales into a singularly perturbed model is developed. The algorithm leverages an ordered Schur decomposition that arranges the system’s eigenvalues along the main diagonal. Furthermore, the algorithm is computationally efficient, which is essential from a practical perspective (e.g., parametric studies). Next, the singularly perturbed model of the overall PEMFC-FPS is derived and formalized. Lastly, LQ controllers for each of the lower order models are designed so that the response meets operational requirements.

The paper is organized as follows. In Section 2, the preliminaries of PEMFCs and singularly perturbed systems are introduced. In Section 3, theoretical details on time-scale decoupling via the ordered Schur decomposition are presented. Section 4 covers the singularly perturbed modeling of the PEMFC-FPS, followed by Section 5 where LQ controller design is illustrated. The remaining sections, namely, Section 6 and Section 7, correspond to the Results and Discussion, and Conclusions, respectively.

2. Preliminaries

2.1. PEM Fuel Cells

The purpose of this subsection is to familiarize the reader with the PEMFC. For more details on the topic, the reader is referred to [9,10,13] and the references therein. A simplified fuel cell diagram is 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.

2.2. 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

$$\begin{array}{cc}\hfill {\dot{x}}_1\left(t\right)& =A_1{x}_1\left(t\right)+A_2{x}_2\left(t\right)+B_{1}u\left(t\right)\hfill \\ \hfill \epsilon {\dot{x}}_2\left(t\right)& =A_3{x}_1\left(t\right)+A_4{x}_2\left(t\right)+B_{2}u\left(t\right)\hfill \\ \hfill y\left(t\right)& =C_1{x}_1\left(t\right)+C_2{x}_2\left(t\right),\hfill \end{array}$$

(1)

where $x_1\left(t\right)\in {\mathbb{R}}^n$ and $x_2\left(t\right)\in {\mathbb{R}}^m$ are the slow and fast state variables respectively, $u\left(t\right)\in {\mathbb{R}}^p$ are the system control inputs, $y\left(t\right)\in {\mathbb{R}}^q$ are the system measurements, and $\epsilon$ is the small singular perturbation parameter $0<\epsilon \ll 1$. In a real physical system, $\epsilon$ 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 $\mathcal{O}\left(1\right)$ eigenvalues close to the imaginary axis and fast $\mathcal{O}\left(\frac{1}{\epsilon }\right)$ eigenvalues far from it. The following standard assumption is imposed [19].

Assumption 1.

Matrix $A_4$ is nonsingular.

When $\epsilon =0$ (a common strategy for order reduction), the following reduced-order system corresponding to the slow dynamics is obtained:

$$\begin{array}{cc}\hfill {\dot{\overline{x}}}_1\left(t\right)& =A_0{\overline{x}}_1\left(t\right)+B_{0}u\left(t\right)\hfill \\ \hfill \overline{y}\left(t\right)& =C_0{\overline{x}}_1\left(t\right)+D_{0}u\left(t\right),\hfill \end{array}$$

(2)

where

$$\begin{array}{c}{A}_0=A_1-A_2{A}_4^{-1}{A}_3{B}_0=B_1-A_2{A}_4^{-1}{B}_2\\ C_0=C_1-C_2{A}_4^{-1}{A}_3{D}_0=-C_2{A}_4^{-1}{B}_2.\end{array}$$

(3)

The approximated fast subsystem is [19]

$$\begin{array}{cc}\hfill \epsilon {\dot{\overline{x}}}_2\left(\tau \right)& =A_4{\overline{x}}_2\left(\tau \right)+B_{2}u\left(\tau \right)\hfill \\ \hfill {\overline{y}}_f\left(\tau \right)& =C_2{\overline{x}}_2\left(\tau \right),\hfill \end{array}$$

(4)

where $\tau =t/\epsilon$. According to the theory of singular perturbations [15,19], the approximation obtained in (2)–(4) satisfies

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill x_1\left(t\right)& ={\overline{x}}_1\left(t\right)+\mathcal{O}\left(\epsilon \right),\forall t\ge t_0\hfill \\ \hfill x_2\left(t\right)& ={\overline{x}}_2\left(t\right)-A_4^{-1}(A_3{\overline{x}}_1\left(t\right)+B_{2}u\left(t\right))+O\left(\epsilon \right),\hfill \\ & \forall t\ge t_0+O\left(\epsilon \right).\hfill \end{array}\end{array}$$

(5)

Hence, the smaller $\epsilon$, the better the approximation. Furthermore, $\lambda \left(A_0\right)$ and $\lambda (A_4/\epsilon )$ are $\mathcal{O}\left(\epsilon \right)$ 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.

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right]& =\left[\begin{array}{cc}I-\epsilon HL& \epsilon H\\ L& I\end{array}\right]\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]=T\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]& =\left[\begin{array}{cc}I& \epsilon H\\ -L& I-\epsilon LH\end{array}\right]\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right]=T^{-1}\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right].\hfill \end{array}$$

(6b)

After (6) is applied to (1), the decoupled system is then

$$\begin{array}{cc}\hfill {\dot{z}}_1\left(t\right)& =A_s{z}_1\left(t\right)+B_{s}u\left(t\right)\hfill \\ \hfill \epsilon {\dot{z}}_2\left(t\right)& =A_f{z}_2\left(t\right)+B_{f}u\left(t\right)\hfill \\ \hfill y\left(t\right)& =C_s{z}_1\left(t\right)+C_f{z}_2\left(t\right),\hfill \end{array}$$

(7)

where

$$\begin{array}{cc}\hfill A_s& =A_1-A_{2}L{B}_s=B_1-H{B}_2-\epsilon HL{B}_1\hfill \\ \hfill A_f& =A_4+\epsilon L{A}_2{B}_f=B_2+\epsilon L{B}_1\hfill \\ \hfill C_s& =C_1-C_{2}L{C}_f=C_2-\epsilon C_{2}LH+\epsilon C_{1}H.\hfill \end{array}$$

(8)

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\left(t\right)$ while $\overline{y}\left(t\right)$ from (2) includes a $D_{0}u\left(t\right)$ term.

Matrices L and H are obtained by solving the following equations.

$$\begin{array}{c}{A}_{4}L-A_3-\epsilon L(A_1-A_{2}L)=\mathbf{0}\\ H{A}_4-A_2+\epsilon (HL{A}_2-A_{1}H+A_{2}LH)=\mathbf{0}.\end{array}$$

(9)

The reader can refer to [23,24,25] for methods on finding the solution of L and H.

3. 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.

3.1. 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 ${\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_N$ are either known or can be easily determined. ${\epsilon }_1$ is associated with the slowest state variable and ${\epsilon }_N$ is associated with the fastest, namely, $0<{\epsilon }_N\ll {\epsilon }_{N-1}\ll \dots \ll {\epsilon }_1=1$.

Consider a general implicit multiple time-scale system without inputs, as shown in (10).

$$\dot{x}\left(t\right)=A\left(\epsilon \right)x\left(t\right),$$

(10)

where $x\left(t\right)={[x_1\left(t\right)x_2\left(t\right)\dots x_N\left(t\right)]}^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\in {\mathbb{R}}^{n\times n}$, there exits a unitary matrix $T\in {\mathbb{R}}^{n\times n}$ such that $T^{T}AT=\tilde{A}$ is upper quasi-triangular.

$$\tilde{A}=\left[\begin{array}{cccc}{R}_{11}& R_{12}& \cdots & R_{1N}\\ \mathbf{0}& R_{22}& \cdots & R_{2N}\\ ⋮& \ddots & \ddots & ⋮\\ \mathbf{0}& \cdots & \mathbf{0}& R_{NN}\end{array}\right].$$

(11)

Matrix blocks $R_{ij},i=j$ can be $1\times 1$ or $2\times 2$. $1\times 1$ blocks correspond to real eigenvalues while $2\times 2$ blocks correspond to complex eigenvalues.

The eigenvalues appearing along the diagonal of $\tilde{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\left(t\right)=Tx\left(t\right)$ [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

$$\dot{z}\left(t\right)=\tilde{A}z\left(t\right),$$

(12)

where

$$\begin{array}{c}\hfill \tilde{A}=\left[\begin{array}{cccc}{A}_{11}& {\mathsf{\Psi }}_{12}& \dots & {\mathsf{\Psi }}_{1N}\\ \mathbf{0}& A_{22}& \dots & {\mathsf{\Psi }}_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \dots & A_{NN}\end{array}\right].\end{array}$$

(13)

Diagonal block matrices $A_{ij}$, $i=j$, contain the system’s eigenvalues in descending order. $z{\left(t\right)}_1,z{\left(t\right)}_2,\dots ,z{\left(t\right)}_N$ are vectors each representing each time scale and ${\mathsf{\Psi }}_{ij}$ are matrix blocks of appropriate dimensions. Note that blocks $A_{ij}$, $i=j$, represent individual time-scales rather than individual eigenvalues.

3.2. 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), $\epsilon$ is commonly evaluated as $\epsilon =\max \mathrm{Re}\left\{{\lambda }_s\right\}/\min \mathrm{Re}\left\{{\lambda }_f\right\}$ [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

$${\epsilon }_i=\frac{\max |{\lambda }_{slowest}|}{\min |{\lambda }_{i_{fast}}|},i=2,3,\cdots ,N.$$

(14)

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:

$$\left[\begin{array}{c}{\epsilon }_1{\dot{z}}_1\left(t\right)\\ {\epsilon }_2{\dot{z}}_2\left(t\right)\\ ⋮\\ {\epsilon }_N{\dot{z}}_N\left(t\right)\end{array}\right]=\left[\begin{array}{cccc}{\tilde{A}}_{11}& {\tilde{\mathsf{\Psi }}}_{12}& \dots & {\tilde{\mathsf{\Psi }}}_{1N}\\ \mathbf{0}& {\tilde{A}}_{22}& \dots & {\tilde{\mathsf{\Psi }}}_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \dots & {\tilde{A}}_{NN}\end{array}\right]\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\\ ⋮\\ z_N\left(t\right)\end{array}\right],$$

(15)

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 (6).

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.

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\dot{\mathrm{Z}}}_1\left(t\right)\\ {\epsilon }_N{\dot{\mathrm{Z}}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{\tilde{A}}_1& {\tilde{A}}_2\\ {\tilde{A}}_3& {\tilde{A}}_4\end{array}\right]\left[\begin{array}{c}{\mathrm{Z}}_1\left(t\right)\\ {\mathrm{Z}}_2\left(t\right)\end{array}\right].\end{array}$$

(16)

In (16), matrices ${\tilde{A}}_3$ and ${\tilde{A}}_4$ represent the last row of the system’s matrix in (13) with ${\epsilon }_N$ extracted. ${\tilde{A}}_4$ represents the fastest time-scale, while ${\tilde{A}}_3$ contains the rest of the matrix blocks which happen to be all zero in this case. ${\tilde{A}}_1$ and ${\tilde{A}}_2$ are matrices of appropriate dimensions containing the rest of the system matrix in (12).

Utilizing (6), 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.

$$\begin{array}{cc}\hfill \dot{\xi }\left(t\right)& =({\tilde{A}}_1-{\tilde{A}}_2{L}_N)\xi \left(t\right)\hfill \end{array}$$

(17a)

$$\begin{array}{cc}\hfill {\epsilon }_N{\dot{\eta }}_N\left(t\right)& =({\tilde{A}}_4+{\epsilon }_N{L}_N{\tilde{A}}_2){\eta }_N\left(t\right).\hfill \end{array}$$

(17b)

As a reminder to the reader, matrices L and H satisfy (9). A modified Newton’s method developed in [25] follows.

$$\begin{array}{c}{D}_1^{\left(i\right)}{L}_N^{(i+1)}+L_N^{(i+1)}{D}_2^{\left(i\right)}=Q^{\left(i\right)}\\ D_1^{\left(i\right)}={\tilde{A}}_4+\epsilon L_N^{\left(i\right)}{\tilde{A}}_2{D}_2^{\left(i\right)}=-\epsilon ({\tilde{A}}_1-{\tilde{A}}_2{L}_N^{\left(i\right)})\\ Q^{\left(i\right)}={\tilde{A}}_3+\epsilon L_N^{\left(i\right)}{\tilde{A}}_2{L}_N^{\left(i\right)}{L}_N^{\left(0\right)}={\tilde{A}}_4^{-1}{\tilde{A}}_3.\end{array}$$

(18)

The new slow subsystem (17a) is partitioned again, as in (19), where ${\epsilon }_{N-1}$ is now extracted from the second fastest time-scale.

$$\left[\begin{array}{c}{\dot{Y}}_1\left(t\right)\\ {\epsilon }_{N-1}{\dot{Y}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{\tilde{G}}_1& {\tilde{G}}_2\\ {\tilde{G}}_3& {\tilde{G}}_4\end{array}\right]\left[\begin{array}{c}{Y}_1\left(t\right)\\ Y_2\left(t\right)\end{array}\right].$$

(19)

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

$$\left[\begin{array}{c}{\eta }_1\left(t\right)\\ {\eta }_2\left(t\right)\\ ⋮\\ {\eta }_N\left(t\right)\end{array}\right]=T\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\\ ⋮\\ z_N\left(t\right)\end{array}\right],$$

(20)

where T is defined as $T={\overline{T}}_2{\overline{T}}_3\cdots {\overline{T}}_{N-1}{T}_N$ and ${\overline{T}}_i$ is

$${\overline{T}}_i\triangleq \left[\begin{array}{cc}{T}_i& \mathbf{0}\\ \mathbf{0}& I\end{array}\right].$$

(21)

Matrix $T_i$ represents the linear transformation for each time-scale i. On the other hand, matrix ${\overline{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), ${\tilde{A}}_3$ in (16) is $\mathbf{0}$. Then, equations for the solution of matrices L and H in (9) simplify to the following.

$$\begin{array}{cc}\hfill & {\tilde{A}}_4{L}_N-\epsilon L_N({\tilde{A}}_1-{\tilde{A}}_2{L}_N)=\mathbf{0}\hfill \end{array}$$

(22a)

$$\begin{array}{cc}\hfill & \epsilon (H_N{L}_N{\tilde{A}}_2-{\tilde{A}}_1{H}_N+{\tilde{A}}_2{L}_N{H}_N)-{\tilde{A}}_2+H{\tilde{A}}_4=\mathbf{0}.\hfill \end{array}$$

(22b)

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}_N^{\left(1\right)}=\mathbf{0}.$$

(23)

Matrix M is defined as $(I_n\otimes {\tilde{A}}_4-\epsilon {\tilde{A}}_1^T\otimes I_n)$ and is full rank. Therefore, $ker\left(M\right)=0$, which implies $L_N^{\left(1\right)}=0$. Since ${L_N}^{\left(0\right)}=0$ and $L_N^{\left(1\right)}=0$, then by induction, for all other iterations $L_N^{\left(i\right)}=0$. This applies to matrix $L_i,i=N-1,N-2,\dots ,1$, for the remaining time-scales. □

In the absence of matrix L, (22b) then becomes just a Sylvester equation.

$$H_N{\tilde{A}}_4-{\tilde{A}}_2-\epsilon {\tilde{A}}_1{H}_N=\mathbf{0}.$$

(24)

Likewise, transformation (6a) simplifies to (25). Note that the state variables in (25) are arbitrary.

$$\left[\begin{array}{c}\xi \left(t\right)\\ \eta \left(t\right)\end{array}\right]=\left[\begin{array}{cc}I& -\epsilon H\\ \mathbf{0}& I\end{array}\right]\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\end{array}\right].$$

(25)

After the process is repeated for all the N time-scales available in the system, the individual subsystems are then given as follows.

$$\begin{array}{cc}\hfill {\epsilon }_1{\dot{\eta }}_1\left(t\right)& ={\tilde{A}}_1{\eta }_1\left(t\right)\hfill \\ \hfill {\epsilon }_2{\dot{\eta }}_2\left(t\right)& ={\tilde{A}}_2{\eta }_2\left(t\right)\hfill \\ \hfill & ⋮\hfill \\ \hfill {\epsilon }_N{\dot{\eta }}_N\left(t\right)& ={\tilde{A}}_N{\eta }_N\left(t\right).\hfill \end{array}$$

(26)

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 $B_s=B_1-H{B}_2-\epsilon HL{B}_1$ and $B_f=B_2+\epsilon L{B}_1$, respectively.

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    3:       apply swapping algorithm [27] or [29]    4: Evaluate ${\epsilon }_i$ for each time-scale and form explicit system    5: Use Chang transformation [22] to get individual time-scales    6: Output: Completely decoupled system

4. 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 ${}^{\circ }$C to 100 ${}^{\circ }$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_1^o=40$ kW and oxygen excess ratio is $z_2^o=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.

4.1. 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 $\delta x\in {\mathbb{R}}^8,\delta x=\left[\begin{array}{cccccccc}{m}_{O_2}& m_{H_2}& m_{N_2}& {\omega }_{cp}& p_{sm}& m_{sm}& m_{w,an}& p_{rm}\end{array}\right]$. Corresponding descriptions of state variables is provided in

Table 1. PEMFC nomenclature.

Symbol	Variable Definition	Units
$m_{O_2}$	Mass of oxygen	kg
$m_{H_2}$	Mass of hydrogen	kg
$m_{N_2}$	Mass of nitrogen	kg
${\omega }_{cp}$	Compressor speed	rad/s
$p_{sm}$	Pressure of gas in supply manifold	kPa
$m_{sm}$	Mass of gas in supply manifold	kg
$m_{w,an}$	Mass of water in anode channel	kg
$p_{rm}$	Pressure in return manifold	kPa

.

Likewise, for the FPS, masses and pressures of gases are essential for modeling. The state vector in this case is $\delta x\in {\mathbb{R}}^{10},\delta x=\left[\begin{array}{cccccccccc}{T}_{cpox}& p_{H_2}^{an}& p^{an}& p^{hex}& {\omega }_{blower}& p^{hds}& p_{{CH}_4}^{mix}& p_{air}^{mix}& p_{H_2}^{wrox}& p_{wrox}\end{array}\right]$.

Table 2. Fuel processing system nomenclature.

Symbol	Variable Definition	Units
$T_{cpox}$	Catalyst temperature	K
$p_{H_2}^{an}$	Pressure of hydrogen in anode	kPa
$p^{an}$	Anode pressure	kPa
$p^{hex}$	Heat-exchanger pressure	kPa
${\omega }_{blower}$	Speed of the blower (rad/sec)	kPa
$p^{hds}$	Pressure of hydrodesulfurizer	kPa
$p_{{CH}_4}^{mix}$	Pressure of $C{H}_4$ in mixer	kPa
$p_{air}^{mix}$	Pressure of air in mixer	kPa
$p_{H_2}^{wrox}$	Hydrogen pressure in water gas shift converter (WROX)	kPa
$p_{wrox}$	Total pressure in WROX	kPa

describes the state variables of the tenth order FPS model.

Using information from the PEMFC and FPS models, an augmented linear state-space model is created as in (27) [14].

$$\begin{array}{cc}\hfill \left[\begin{array}{c}\delta {\dot{x}}_{PEMFC}\left(t\right)\\ \delta {\dot{x}}_{FPS}\left(t\right)\end{array}\right]& =\left[\begin{array}{cc}{A}_{PEMFC}& \mathbf{0}\\ \mathbf{0}& A_{FPS}\end{array}\right]\left[\begin{array}{c}\delta x_{PEMFC}\left(t\right)\\ \delta x_{FPS}\left(t\right)\end{array}\right]+\left[\begin{array}{cc}{B}_{PEMFC}& \mathbf{0}\\ B_{{FPS}_1}& B_{{FPS}_2}\end{array}\right]\left[\begin{array}{c}\delta u_{blow}\\ \delta u_{val}\end{array}\right]\hfill \\ \hfill \left[\begin{array}{c}{y}_{PEMFC}\left(t\right)\\ y_{FPS}\left(t\right)\end{array}\right]& =\left[\begin{array}{cc}{C}_{PEMFC}& \mathbf{0}\\ \mathbf{0}& C_{FPS}\end{array}\right]\left[\begin{array}{c}\delta x_{PEMFC}\left(t\right)\\ \delta x_{FPS}\left(t\right)\end{array}\right]+\left[\begin{array}{cc}{D}_{PEMFC}^{blow}& \mathbf{0}\\ D_{FPS}^{blow}& D_{FPS}^{val}\end{array}\right]\left[\begin{array}{c}\delta u_{blow}\\ \delta u_{val}\end{array}\right].\hfill \end{array}$$

(27)

State-space system in (27) can be written as follows:

$$\begin{array}{cc}\hfill \delta \dot{x}\left(t\right)& =A_{aug}\delta x\left(t\right)+B_{aug}\delta u\left(t\right)\hfill \\ \hfill y\left(t\right)& =C_{aug}\delta x\left(t\right)+D_{aug}\delta u\left(t\right),\hfill \end{array}$$

(28)

where $A_{aug}\in {\mathbb{R}}^{18\times 18}$, $B_{aug}\in {\mathbb{R}}^{18\times 2}$, $C_{aug}\in {\mathbb{R}}^{5\times 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}=\mathbf{0}$.

In the following subsection, a singularly perturbed linear model of (27) is derived using methods from Section 3.

4.2. 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 ${\epsilon }_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 $|{\lambda }_9|=3.33$. The magnitude of the smallest eigenvalues of the second and third time-scales are $|{\lambda }_{10}|=12.17$ and $|{\lambda }_{16}|=157.90$. Hence, ${\epsilon }_1$, ${\epsilon }_2$, and ${\epsilon }_3$ are as follows.

$${\epsilon }_1=1,{\epsilon }_2=\frac{3.33}{12.17}\approx 0.27,{\epsilon }_3=\frac{3.33}{157.90}\approx 0.02.$$

(29)

The singularly perturbed model then becomes

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{z}}_1\left(t\right)\\ {\epsilon }_2{\dot{z}}_2\left(t\right)\\ {\epsilon }_3{\dot{z}}_3\left(t\right)\end{array}\right]& =\left[\begin{array}{ccc}{\tilde{A}}_{11}& {\tilde{\mathsf{\Psi }}}_{12}& {\tilde{\mathsf{\Psi }}}_{13}\\ \mathbf{0}& {\tilde{A}}_{22}& {\tilde{\mathsf{\Psi }}}_{23}\\ \mathbf{0}& \mathbf{0}& {\tilde{A}}_{33}\end{array}\right]\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\\ z_3\left(t\right)\end{array}\right]+\left[\begin{array}{c}{\tilde{B}}_1\\ {\tilde{B}}_2\\ {\tilde{B}}_3\end{array}\right]\left[\begin{array}{c}\delta u_{blow}\\ \delta u_{val}\end{array}\right]\hfill \\ \hfill y\left(t\right)& =\left[\begin{array}{ccc}{\tilde{C}}_1& {\tilde{C}}_2& {\tilde{C}}_3\end{array}\right]\left[\begin{array}{c}{z}_1\left(t\right)\\ z_2\left(t\right)\\ z_3\left(t\right)\end{array}\right].\hfill \end{array}$$

(30)

State variables belonging to each time-scale are presented in

Table 3. State variable distribution among time-scales.

	State Variables
Time-scale 1	$p_{{CH}_4}^{mix}$   $m_{H_2}$   ${\omega }_{blower}$   $m_{O_2}$   $p^{hex}$   $p_{air}^{mix}$   $p_{H_2}^{wrox}$   $p^{an}$   $p_{wrox}$
Time-scale 2	$p^{hds}$   $p_{rm}$   $m_{w,an}$   $p_{sm}$   $p_{H_2}^{an}$   ${\omega }_{cp}$
Time-scale 3	$T_{cpox}$   $m_{N_2}$   $m_{sm}$

.

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}=T{B}_{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 $\epsilon$ of the corresponding time-scale.

Lastly, the Chang transformation [22] is applied sequentially to obtain a completely decoupled system (refer to Algorithm 1).

5. 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.

5.1. LQ Control Overview

In LQ control, our objective is to drive the state vector $x\left(t\right)$ 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\left(t\right)$ 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\left(t\right)$ will be slow and a small control action will be required. This trade-off between the state-vector $x\left(t\right)$ and control action $u\left(t\right)$ can be represented by the following cost function.

$$J=\frac{1}{2}{\int }_0^{\infty }\left[x^T\left(t\right)Qx\left(t\right)+u^T\left(t\right)Ru\left(t\right)\right]dt.$$

(31)

The LQ problem is then to find a control $u\left(t\right)$ such that the cost function (31) is minimized. An important assumption here is that $R>0$ and $Q\ge 0$. The solution to this infinite-horizon optimal control problem is a control input given by [32]

$$u\left(t\right)=-Fx\left(t\right),$$

(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].

$$A^{T}P+PA+Q-PB{R}^{-1}{B}^{T}P=\mathbf{0}.$$

(33)

The state vectors $x\left(t\right)$ for each time-scale are given in Table 3.

5.2. 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 (6) to get a decoupled three time-scale model as follows.

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\dot{\eta }}_1\left(t\right)\\ {\epsilon }_2{\dot{\eta }}_2\left(t\right)\\ {\epsilon }_3{\dot{\eta }}_3\left(t\right)\end{array}\right]& =\left[\begin{array}{ccc}{\widehat{A}}_{11}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\widehat{A}}_{22}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\widehat{A}}_{33}\end{array}\right]\left[\begin{array}{c}{\eta }_1\left(t\right)\\ {\eta }_2\left(t\right)\\ {\eta }_3\left(t\right)\end{array}\right]+\left[\begin{array}{c}{\widehat{B}}_1\\ {\widehat{B}}_2\\ {\widehat{B}}_3\end{array}\right]\left[\begin{array}{c}\delta u_{blow}\\ \delta u_{val}\end{array}\right].\hfill \end{array}$$

(34)

In (34), ${\eta }_1\left(t\right)\in {\mathbb{R}}^9$, ${\eta }_2\left(t\right)\in {\mathbb{R}}^6$, and ${\eta }_3\left(t\right)\in {\mathbb{R}}^3$. The overall similarity transformation to arrive to (34) from original model (27) is $T=T_{ordSchur}{T}_{Chang}$, where $T_{Chang}$ is in itself two different linear transformations, namely, $T_{Chang}={\overline{T}}_{Chang2}{T}_{Chang3}$. Note that ${\overline{T}}_{Chang2}$ has the form given in (21).

The control input vector for system (34) now becomes as follows.

$$u\left(t\right)=\left[\begin{array}{c}\delta u_{blow}\\ \delta u_{val}\end{array}\right]=-F\left[\begin{array}{c}{\eta }_1\left(t\right)\\ {\eta }_2\left(t\right)\\ {\eta }_3\left(t\right)\end{array}\right].$$

(35)

Note that it is important to convert the control input vector into the original coordinates of the system [33]. To do so, the transformation ${\eta }_i\left(t\right)=T_{Chang}{z}_i\left(t\right)$ 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\left(t\right)=T_{ordSchur}{x}_i\left(t\right)$. The control input them becomes

$$u\left(x\left(t\right)\right)=-F{T}_{Chang}{T}_{ordSchur}x\left(t\right).$$

(36)

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.

5.3. 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. Fuel cell system operating time requirements.

Process	Operating Time
Electrochemistry	$\mathcal{O}({10}^{-19}$ s)
Hydrogen & air manifolds	$\mathcal{O}({10}^{-1}$ s)
Flow control/supercharging devices	$\mathcal{O}({10}^0$ s)

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 $\mathcal{O}\left({10}^0\right)$ 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).

$$Q=diag\left\{\begin{array}{ccccccccc}3613& 7417& 7059& 7009& 62& 3743& 9015& 3183& 5971\end{array}\right\},R=298.$$

(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 $\mathcal{O}\left({10}^{-1}\right)$ 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.

The parameters for this desired design are given as follows.

$$Q=diag\left\{\begin{array}{cccccc}\mathrm{8,223,940}& \mathrm{251,505}& \mathrm{4,144,289}& \mathrm{7,314,075}& \mathrm{7,813,740}& \mathrm{3,672,859}\end{array}\right\},R=1.$$

(38)

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.

Matrices Q and R for this design are given in (39).

$$Q=diag\left\{\begin{array}{ccc}131,830& 123,500& 190,903\end{array}\right\},R=1.$$

(39)

6. Discussion of Results

Results of this research can be divided into two main parts:

• Introduction of an algorithm that converts implicit singularly perturbed systems into explicit ones.
• Singularly perturbed modeling of the PEMFC-FPS augmented model followed by LQ controller design.

6.1. 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.

6.2. 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.

6.3. 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.

7. 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.