# Towards Reduced-Order Models of Solid Oxide Fuel Cell

Authors: Maciej Ławryńczuk
Journal: Complexity
 (2018)

---

## Abstract

The objective of this work is to find precise reduced-order discrete-time models of a solid oxide fuel cell, which is a multiple-input multiple-output dynamic process. At first, the full-order discrete-time model is found from the continuous-time first-principle description. Next, the discrete-time submodels of hydrogen, oxygen, and water pressures (intermediate variables) are reduced. Two model reduction methods based on observability and controllability Grammians are compared: the state truncation method and reduction by residualisation. In all comparisons, the second method gives better results in terms of dynamic and steady-state errors as well as Nyquist plots. Next, the influence of the order of the pressure models on the errors of the process outputs (the voltage and the pressure difference) is studied. It is found that the number of pressure model parameters may be reduced from 25 to 19 without any deterioration of model accuracy. Two suboptimal reduced models are also discussed with only 14 and 11 pressure parameters, which give dynamic trajectories and steady-state characteristics that are very similar to those obtained from the full-order structure.


---

## Body

## 1. Introduction

There are three important reasons why renewable energy sources are becoming more and more popular. Firstly, burning of fossil fuels leads to air pollution and significant climate changes. Secondly, both mentioned phenomena badly affect public health, which has a very negative effect on the economy. Thirdly, fossil fuels are located in some countries only, whereas sources of renewable energy actually exist practically in all countries (of course not all of them are possible in all locations). Access to energy sources naturally increases national security. The countries that do not have fossil fuels may switch to renewable energy and become independent of other countries. As a result, the role of renewable energy is important and it is expected to grow fast in the future [1, 2]. Usually, renewable energy is obtained from wind turbines [3], geothermal systems [4], solar collectors [5], marine systems [6], and biofuels [7]. Additionally, solid oxide fuel cells (SOFCs) are very promising sources of energy. SOFCs are electrochemical devices that are able to directly convert the chemical energy stored in hydrocarbon fuels into electrical energy [8, 9]. They have many advantages, namely, high electrical efficiency, fuel flexibility, low emission, quiet operation, and relatively low cost. That is why SOFCs are expected to become sound alternatives to conventional power generation schemes not only for domestic but also for commercial and industrial sectors.Economically efficient and technologically safe operation of SOFCs requires well-designed control algorithms. Control of SOFCs is a challenging task, since they are nonlinear dynamic systems and it is essential to precisely satisfy some technological constraints that must be imposed on process variables [10]. Hence, for controlling SOFCs, advanced Model Predictive Control (MPC) algorithms are preferred rather than classical Proportional-Integral-Derivative (PID) controllers. An important feature of MPC is the fact that a mathematical model of the controlled SOFC is used online to successively calculate the best possible sequence of manipulated variable(s). In the simplest case, for prediction in MPC, linear models of the process may be used [11, 12] and the resulting control quality is better than that of the classical PID. However, in order to obtain very good control quality, a nonlinear dynamic model of the SOFC must be used in MPC. Different variants of nonlinear MPC algorithms for the SOFC are discussed in [13–16] (different model structures and optimisation algorithms are possible). Important applications of the mathematical model of the SOFC also include process optimisation [17], fault tolerant control [15], and estimation [16, 18].In MPC, optimisation, fault tolerant control, and estimation, different model structures may be used. Firstly, the first-principle model based on technological laws may be used [13, 16]. Secondly, empirical (black-box) models are possible, for example, neural networks [14, 19] or fuzzy systems [20, 21]. Although it may be easier to use empirical models than rigorous first-principle ones, it is necessary to point out three important disadvantages of black-box structures. Empirical structures make it possible to predict the sequence of the output variable(s) for a given sequence of the input and disturbance variable(s), but some intermediate process variables are usually not modelled. Moreover, frequently used black-bock models typically have numerous parameters, much more than the fundamental ones. Finally, accuracy of black-box models may be good in typical operating conditions, but for other ones they are likely to generate output value(s) far from those calculated by the fundamental models (and the real process).This work is concerned with the fundamental model of the SOFC, which is a multiple-input multiple-output nonlinear dynamic process. The objective is to find precise reduced-order discrete-time models. To achieve this goal, the discrete-time submodels of hydrogen, oxygen, and water pressures (intermediate variables) are reduced by means of two methods. In the first approach, the balanced Grammian of the state-space realisation is found and the state variables corresponding to small entries of the Grammian are removed. In the second reduction method, the model parameters are additionally adjusted in such a way that the steady-state gain of the reduced model is equal to that of the full-order one. The reduced models are compared in terms of dynamic and steady-state errors as well as Nyquist plots. Next, the influence of the order of the pressure models on the errors of the process outputs (the voltage and the pressure difference) is studied. The ideal reduced model and suboptimal ones, which give good compromise between accuracy and complexity, are discussed and compared with the full-order structure.This work is structured as follows. The SOFC is shortly described in Section2 and its full-order continuous-time model is detailed in Section 3. Section 4 derives the full-order discrete-time model. Section 5 discusses two model reduction methods. The main part of the paper, presented in Section 6, at first details reduction of the hydrogen, oxygen, and water pressure models and next studies the influence of the reduced pressure models on modelling accuracy of two process outputs. Finally, Section 7 concludes the paper.


## 2. SOFC System Description

The literature concerned with first-principle modelling of SOFCs is rich, for example, [9, 22–24]. The fundamental model of the SOFC introduced in [25] and next discussed in [12, 26] is considered here. The following assumptions are made:(1)
The gases are ideal.(2)
The stack is fed with hydrogen and air.(3)
The channels that transport gases along the electrodes have a fixed volume, but their lengths are small. Hence, it is only necessary to define one single pressure value in their interior.(4)
The exhaust of each channel is via a single orifice. The ratio of pressures between the interior and exterior of the channel is large enough to consider that the orifice is choked.(5)
The temperature is stable at all times.(6)
Because the working conditions are not close to the upper and lower extremes of current, the only source of losses is ohmic losses.(7)
The Nernst equation can be applied.The considered SOFC has two manipulated variables (the inputs of the process): the input gas flow rateqfin (mols-1) and the input oxygen flow rate qO2in (mols-1); one disturbance (the uncontrolled input) Iin which is the external current load (A) and four controlled variables (the outputs of the process): the stack output voltage Vr (V), fuel utilisation Uf (−), the fuel cell pressure difference pdif (atm) between the hydrogen and oxygen passing through the anode and cathode gas compartments, and the ratio RHO (−) between hydrogen and oxygen flow rates. The partial pressures of hydrogen, oxygen, and water are denoted by pH2, pO2, and pH2O, respectively (atm). The input hydrogen flow and the hydrogen reacted flow are denoted by qH2in and qH2r, respectively (mols-1).


## 3. Continuous-Time Model

Figure1 depicts the structure of the fundamental continuous-time model of the SOFC system. The pressure of hydrogen is(1)pH2=1/KH2τH2s+1qH2in-2KrIr,where the input hydrogen flow is(2)qH2in=1τfs+1qfin.Hence, the pressure of hydrogen is(3)pH2=1/KH2τH2s+11τfs+1qfin-2Kr1τes+1Iin,where KH2, τH2, τf, and τe denote the valve molar constant for hydrogen, the response time of hydrogen flow, the fuel processor time constant, and the electrical time constant, respectively. The pressure of oxygen is(4)pO2=1/KO2τO2s+1qO2in-KrIr=1/KO2τO2s+1qO2in-Kr1τes+1Iin,where KO2 and τO2 denote the valve molar constant for oxygen and the response time of oxygen flow, respectively. The pressure of water is(5)pH2O=1/KH2OτH2Os+1qH2r, where the hydrogen flow that reacts is(6)qH2r=2KrIr=2Krτes+1Iin.Hence, the pressure of water is(7)pH2O=2Kr/KH2Oτes+1τH2Os+1Iin,where KH2O and τH2O denote the valve molar constant for water and the response time of water flow, respectively. Finally, outputs of the process are defined. Applying Nernst’s equation and taking into account ohmic losses, the stack output voltage is(8)Vr=N0E0+R0T02F0lnpH2pO2pH2O-Irr,where N0, E0, R0, T0, and F0 denote the number of cells in series in the stack, the ideal standard potential, the universal gas constant, the absolute temperature, and Faraday’s constant, respectively. Fuel utilisation is defined as the ratio between the hydrogen flow that reacts and the input hydrogen flow.(9)Uf=qH2rqH2in.The ratio between hydrogen and oxygen flow rates is(10)RHO=qH2inqO2in.The pressure difference between the hydrogen and oxygen passing through the anode and cathode gas compartments is(11)pdif=pH2-pO2.All things considered, the continuous-time fundamental model consists of (2), (3), (4), (6), (7), (8), (9), (10), and (11). The values of the model parameters are given in Table 1. Table 2 gives values of process variables for the initial operating point. The values of process inputs are constrained: 0mols-1≤qfin≤1.7023mols-1 and 0mols-1≤qO2in≤1.6134mols-1; the value of the disturbance is also limited: 200A≤Iin≤300A.Table 1
Parameters of the fundamental continuous-time model of the SOFC system.
Parameter Value Unit Description E0 1.18 V Ideal standard potential F0 96.485 C  m o l - 1 Faraday’s constant KH2 8.43 × 1 0 - 1 m o l  s - 1  a t m - 1 Valve molar constant for hydrogen KH2O 2.81 × 1 0 - 1 m o l  s - 1  a t m - 1 Valve molar constant for water KO2 2.52 m o l  s - 1  a t m - 1 Valve molar constant for oxygen Kr 0.996 × 1 0 - 3 m o l  s - 1  A - 1 Constant,Kr=N0/4F0 N0 384 − Number of cells in series in the stack R0 8.314 J  m o l - 1  K - 1 Universal gas constant T0 1273 K Absolute temperature r 0.126 Ω Ohmic loss τe 0.8 s Electrical time constant τf 5 s The fuel processor time constant τH2 26.1 s Response time of hydrogen flow τH2O 78.3 s Response time of water flow τO2 2.91 s Response time of oxygen flowTable 2
The values of process variables for the initial operating point.
Variable Value Unit qfin 0.7023 m o l  s - 1 qO2in 0.6134 m o l  s - 1 Iin 300 A qH2in 0.7023 m o l  s - 1 qH2r 5.976000 × 1 0 - 1 m o l  s - 1 pH2 1.241993 × 1 0 - 1 a t m pO2 1.248413 × 1 0 - 1 a t m pH2O 2.126690 a t m Vr 3.335865 × 1 0 2 V Uf 8.509184 × 1 0 - 1 - RHO 1.144930 - pdif - 6.419816 × 1 0 - 4 a t mFigure 1
The structure of the fundamental continuous-time model of the SOFC system.


## 4. Model Discretisation

Discrete-time versions of the continuous-time transfer functions (3), (4), and (7) must be found. For this purpose, the continuous-time polynomials of the proper order of the variable s must be determined. From (3), one has(12)pH2=1/KH2τH2τfs2+τH2+τfs+1qfin-2Kr/KH2τH2τes2+τH2+τes+1IinFrom (4),(13)pO2=1/KO2τO2s+1qO2in-Kr/KO2τO2τes2+τO2+τes+1IinFrom (7),(14)pH2O=2Kr/KH2OτeτH2Os2+τe+τH2Os+1Iin.Next, model equations are discretised using the first-order holder with the sampling period equal to 1 second. The discrete-time version of (12) is(15)pH2k=β11q-1+β21q-21+α11q-1+α21q-2qfink+β12q-1+β22q-21+α21q-1+α22q-2Iink,where k denotes the discrete sampling instant and q-1 is a unit delay operator. The discrete-time version of (13) is(16)pO2k=β13q-11+α13q-1qO2ink+β14q-1+β24q-21+α14q-1+α24q-2Iink.The discrete-time version of (14) is(17)pH2Ok=β15q-1+β25q-21+α15q-1+α25q-2Iink.Additionally, the discrete-time version of (2)(18)qH2ink=β16q-11+α16q-1qfinkand the discrete-time version of (6)(19)qH2rk=β17q-11+α17q-1Iinkmust be taken into account. The final form of the discrete-time fundamental model is as follows. From the discrete transfer functions (15), (16), (17), (18), and (19), direct dependence of current values of process variables (i.e., for the current sampling instant k) as functions of corresponding variables in the previous instants is found. From (15), one has(20)pH2k=∑i=1npH2bi1,1qfink-i+∑i=1npH2bi1,2Iink-i-∑i=1npH2ai1pH2k-i,where b11,1=β11, b21,1=β21+α12β11, b31,1=α12β21+α22β11, b41,1=α22β21, b11,2=β12, b21,2=β22+α11β12, b31,2=α11β22+α21β12, b41,2=α21β22, a11=α11+α12, a21=α21+α22+α11α12, a31=α11α22+α21α12, and a41=α21α22. From (16), one has(21)pO2k=∑i=1npO2bi2,1qO2ink-i+∑i=1npO2bi2,2Iink-i-∑i=1npO2ai2pO2k-i,where b12,1=β13, b22,1=α14β23, b32,1=α24β13, b12,2=β14, b22,2=β24+α13β14, b32,2=α13β24, a12=α13+α14, a22=α24+α13α14, and a32=α13α24. From (17), one has(22)pH2Ok=∑i=1npH2Obi3Iink-i-∑i=1npH2Oai3pH2Ok-i,where b13=β15, b23=β25, a13=α15, and a23=α25. From (18), one has(23)qH2ink=b14qfink-i-a14qH2ink-i,where b14=β16 and a14=α16. From (19), one has(24)qH2rk=b15Iink-i-a15qH2rk-i,where b15=β17 and a15=α17. The discrete-time version of the Nernst equation (8) is(25)Vrk=N0E0+R0T02F0lnpH2kpO2kpH2Ok-Irkr.In the discrete-time domain fuel utilisation (9) is(26)Ufk=qH2rkqH2inkand the ratio between hydrogen and oxygen flow rates (10) is(27)RHOk=qH2inkqO2inkand the pressure difference (11) is(28)pdifk=pH2k-pO2k.All things considered, the full-order discrete-time model is defined by (20), (21), (22), (23), (24), (25), (26), (27), and (28). The order of model dynamics is defined by the integer numbers: npH2=4, npO2=3, and npH2O=2. Parameters of the full-order fundamental discrete-time model are given in Table 3.Table 3
Parameters of the full-order discrete-time fundamental model.
Equation Parameters (20) a 1 1 = - 3.030057 b 1 1,1 = 4.201819 × 1 0 - 3 b11,2=-3.832050×10-5 a 2 1 = 3.288185 b 2 1,1 = - 1.366717 × 1 0 - 3 b21,2=4.319957×10-5 a 3 1 = - 1.475212 b 3 1,1 = - 3.688449 b31,2=1.443104×10-5 a 4 1 = 2.172670 × 1 0 - 1 b 4 1,1 = 1.070128 × 1 0 - 3 b41,2=-1.974194×10-5  (21) a 1 2 = - 1.704879 b 1 2,1 = 1.154038 × 1 0 - 1 b12,2=-5.160259×10-5 a 2 2 = 9.093077 × 1 0 - 1 b 2 2,1 = - 1.149061 × 1 0 - 1 b22,2=6.187518×10-6 a 3 2 = - 1.440946 × 1 0 - 1 b 3 2,1 = 2.344823 × 1 0 - 2 b32,2=2.156491×10-5  (22) a 1 3 = - 1.273815 b 1 3 = 3.867789 × 1 0 - 5 a 2 3 = 2.828690 × 1 0 - 1 b 2 3 = 2.550836 × 1 0 - 5  (23) a 1 4 = - 8.187307 × 1 0 - 1 b 1 4 = 1.812692 × 1 0 - 1  (24) a 1 5 = - 2.865048 × 1 0 - 1 b 1 5 = 1.421282 × 1 0 - 3To demonstrate effectiveness of the considered model discretisation method, two full-order dynamic models are simulated: the continuous-time one and the discrete-time one. The continuous-time model is simulated in Simulink; the differential equations are solved by the Runge-Kutta method of order 45 with a variable step size. The discrete-time model is implemented in MATLAB; the differential equations are solved with the constant step (equal to the sampling time). Both models are excited by the same series of step changes of the manipulated and disturbance variables depicted in Figure2. Obtained trajectories of both dynamic systems (i.e., the sequences of three pressures and four process outputs) are compared in Figure 3. The discrete-time model gives practically the same responses as the continuous-time one; no important differences are present.Figure 2
The sequences of the manipulated variables and the disturbance used for dynamic simulations.Figure 3
The comparison of the pressures and the outputs calculated from continuous-time and discrete-time models of full order.


## 5. Model Reduction Methods

The state-space representation of a discrete-time dynamic system is considered:(29)xk+1=Axk+Buk,yk=Cxk.It is assumed that matrix A is asymptotically stable. The controllability, Wc, and observability, Wo, Grammians are [27](30)Wc=∑i=0∞AiBBTATi,Wo=∑i=0∞ATiCTCAi. To reduce the order of the dynamic system (29), a balanced state-space representation must be used. For such a system, both Grammians are equal and diagonal:(31)Wc=Wo=Σ=σ10⋯00σ2⋯0⋮⋮⋱⋮00⋯σn. The balanced representation of system (29) is found by setting matrix A to T-1AT, setting matrix B to T-1B, and setting matrix C to CT, where the transformation matrix T may be found by the procedure described in [28, 29].In this work, two model reduction algorithms [30] are considered: model reduction by state truncation (algorithm 1) and model reduction by residualisation (algorithm 2). In algorithm 1, the states for which the corresponding coefficients σi are relatively “small” are simply removed. Model (29) is reformulated:(32)x1k+1x2k+1=A11A12A21A22x1kx2k+B1B2uk,yk=C1C2x1kx2k,where the states x1 and x2 are associated with “big” and “small” values of σi, respectively. In the reduced model, the first state variables are only considered:(33)x1k+1=A11x1k+B1uk,yk=C1x1k. When some of state variables are simply removed from the model, the steady-state properties of the resulting model may significantly differ from those of the original one. To solve the problem, in algorithm 2, the states of model (32) for which the corresponding coefficients σi are “small” are also removed, but afterwards the model matrices are modified to guarantee that the steady-state gains of the full-order and reduced models are the same. Since state variables x2 are assumed to be significantly faster than x1 ones, it follows that x2(k+1)=x2(k) and(34)x2k=I-A22-1A21x1k+B2uk. Hence, the reduced model becomes(35)x1k+1=A11+A12I-A22-1A21x1k+B1+A12I-A22-1B2uk,yk=C1+C2I-A22-1A21x1k+C2I-A22-1B2uk.In the second algorithm in the output equation, there is a possible direct influence of the process input(s) on the output(s), which is not necessary in the first algorithm. Therefore, the second algorithm leads to the difference equation for hydrogen pressure.(36)pH2k=∑i=0npH2bi1,1qfink-i+∑i=0npH2bi1,2Iink-i-∑i=1npH2ai1pH2k-i.In contrast to the original full-order model and the reduced one obtained in the first algorithm (see (20)), the coefficients b01,1 and b01,2 may be necessary. Similarly, the second algorithm leads to the difference equation for oxygen pressure.(37)pO2k=∑i=0npO2bi2,1qO2ink-i+∑i=0npO2bi2,2Iink-i-∑i=1npO2ai2pO2k-i.In contrast to the original full-order model and the reduced one obtained in the first algorithm (see (21)), the coefficients b02,1 and b02,2 may be necessary. Finally, the second algorithm leads to the difference equation for water pressure.(38)pH2Ok=∑i=0npH2Obi3Iink-i-∑i=1npH2Oai3pH2Ok-iin which the coefficient b03 may be necessary in contrast to the original full-order model and the reduced one generated by the first algorithm (see (22)).


## 6. Reduction of the Solid Oxide Fuel Cell Model

At first, separate reduction of hydrogen, oxygen, and water pressure models is considered. Next, the problem of finding the most appropriate orders of dynamics of these models is discussed taking into account two process outputs (the voltage and the pressure difference) whose modelling accuracy depends on the pressure models. Models of the remaining two process outputs, that is, fuel utilisation and the ratio between hydrogen and oxygen flow rates, are not reduced, since they are calculated from very simple equations, (26) and (27).


### 6.1. Reduction of the Hydrogen Pressure Model

The diagonal of the balanced Grammian of the state-space realisation corresponding to the discrete-time difference equation (20) describing hydrogen pressure is(39)g=6.7211×10-17.8946×10-23.7688×10-51.3656×10-17. The last two diagonal entries of the Grammian are small in comparison with the first two ones, which suggests that the last two state variables may be reduced. However, for a thorough analysis, the reduced models of the first, the second, and the third orders are considered; in each case, two reduction algorithms are used. The parameters of the reduced-order hydrogen pressure model of order npH2=1,2,3 calculated by the first algorithm are given in Table 4, whereas the results of the second algorithm are given in Table 5. The general form of the reduced models calculated by the first algorithm is the same as that of the full-order model (see (20)), whereas the reduced model determined by the second method is characterised by (36). In the latter case, the coefficients b01,1 and b01,2 may be necessary, which manifest a direct influence of the model input qfin and the disturbance Iin on the model output pH2. When the reduced model has the third order, the additional coefficients are not necessary to guarantee that the steady-state gain of the reduced model is the same as that of the full-order one. When the reduced model has the second order, only one of them is used, but in the case of the first-order model two of them are necessary.Table 4
Parameters of the reduced-order hydrogen pressure model of ordernpH2 calculated by the first algorithm.
n p H 2 Parameters 1 a 1 1 = - 9.768754 × 1 0 - 1 b 1 1,1 = 3.054446 × 1 0 - 2 b 1 1,2 = - 6.812056 × 1 0 - 5  2 a 1 1 = - 1.781141 b 1 1,1 = 4.201816 × 1 0 - 3 b 1 1,2 = - 6.557908 × 1 0 - 5 a 2 1 = 7.879551 × 1 0 - 1 b 2 1,1 = 3.881002 × 1 0 - 3 b 2 1,2 = 4.914773 × 1 0 - 5  3 a 1 1 = - 2.067646 b 1 1,1 = 4.201819 × 1 0 - 3 b 1 1,2 = - 3.832050 × 1 0 - 5 a 2 1 = 1.298261 b 2 1,1 = 2.677158 × 1 0 - 3 b 2 1,2 = 6.319519 × 1 0 - 6 a 3 1 = - 2.257529 × 1 0 - 1 b 3 1,1 = - 1.111925 × 1 0 - 3 b 3 1,2 = 2.051302 × 1 0 - 5Table 5
Parameters of the reduced-order hydrogen pressure model of ordernpH2 calculated by the second algorithm.
n p H 2 Parameters 1  b 0 1,1 = - 1.345838 × 1 0 - 1 b 0 1,2 = - 3.381314 × 1 0 - 5 a 1 1 = - 9.651884 × 1 0 - 1 b 1 1,1 = 1.758786 × 1 0 - 1 b 1 1,2 = - 4.844625 × 1 0 - 5  2   b 0 1,2 = 3.820428 × 1 0 - 5 a 1 1 = - 1.781141 b 1 1,1 = 4.201816 × 1 0 - 3 b 1 1,2 = - 1.519383 × 1 0 - 4 a 2 1 = 7.879551 × 1 0 - 1 b 2 1,1 = 3.881002 × 1 0 - 3 b 2 1,2 = 9.763310 × 1 0 - 5  3 a 1 1 = - 2.067646 b 1 1,1 = 4.201819 × 1 0 - 3 b 1 1,2 = - 3.832050 × 1 0 - 5 a 2 1 = 1.298261 b 2 1,1 = 2.677158 × 1 0 - 3 b 2 1,2 = 6.319519 × 1 0 - 6 a 3 1 = - 2.257529 × 1 0 - 1 b 3 1,1 = - 1.111925 × 1 0 - 3 b 3 1,2 = 2.051302 × 1 0 - 5To compare full-order and reduced models of the hydrogen pressure, it is convenient to consider their Nyquist plots. It is because Nyquist plots show both steady-state and dynamic behaviours of models. Figure4 depicts the Nyquist plots of reduced-order and full-order models. As expected, the first algorithm does not guarantee that the stead-state gain of the reduced model is exactly the same as that of the original one (the gain may be read for ω=0). Furthermore, for other frequencies, the Nyquist plot of the reduced models is different from that of the full-order one. This observation is true for the simplest first-order model and, for the disturbance Iin to output channel for the second-order structure, increasing the model order makes it possible to obtain precise steady-state gain. The second algorithm is much more efficient, since all reduced models, even the first-order one, have the same gain as the original one. Furthermore, the Nyquist plots of the reduced models are very similar to those of the full-order ones.Figure 4
The Nyquist plots for the hydrogen pressure models of reduced order (npH2=1,2,3) and full order (npH2=4): reduction algorithm 1 (a) and reduction algorithm 2 (b).
(a) (b)The models are also compared in a quantitative way. For this purpose, the dynamic model error is defined:(40)EpH2dyn=1ndyn∑k=1ndynpH2fullk-pH2reducedk2.The outputs of the full-order and reduced models for the sampling instant k are denoted by pH2full(k) and pH2reduced(k), respectively. Both models are excited by the same series of step changes of the manipulated variables and steps of the disturbance depicted in Figure 2; that is, ndyn=2000. Additionally, the steady-state model error is defined:(41)EpH2ss=1nss∑p=1nsspH2fullp-pH2reducedp2.In this case, pH2full(p) and pH2reduced(p) denote outputs of the full-order and reduced steady-state models for the data point p, respectively. They are derived easily from the dynamic models. For the second algorithm, using (36), the dynamic model for hydrogen pressure is (42)pH2=∑i=0npH2bi1,1qfin+∑i=0npH2bi1,2Iin1+∑i=1npH2ai1-1.In the case of the reduced model obtained by the first algorithm (see (20)), b01,1=b01,2=0. For the steady-state models, nss=4116. Dynamic and steady-state errors of the reduced models of hydrogen pressure are given in Table 6. As observed from the Nyquist plots, the second algorithm gives much better results (lower model errors). In particular, always EpH2ss=0 for the second algorithm, since it always guarantees that the steady-state properties of the original and reduced models are the same. When the errors are lower than the machine accuracy (2.2204×10-16), they are treated as 0. Both reduction algorithms give perfect models of order three.Table 6
The dynamic errorEpH2dyn and the steady-state error EpH2ss of the reduced hydrogen pressure models (npH2=1,2,3).
n p H 2 E p H 2 d y n E p H 2 s s Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 2.5687 × 1 0 - 3 1.1652 × 1 0 - 4 1.7068 × 1 0 - 2 0.0 2 2.0740 × 1 0 - 4 3.5651 × 1 0 - 8 4.1544 × 1 0 - 4 0.0 3 0.0 0.0 0.0 0.0Finally, Figure5 compares the hydrogen pressure trajectories (model outputs) calculated by reduced-order (npH2=1,2,3) and full-order (npH2=4) models determined by both algorithms; all models are excited by the same series of step changes of the manipulated variables and steps of the disturbance depicted in Figure 2. As expected, taking into account the Nyquist plots and model errors, the second algorithm gives much better results for the same model order. For the sampling instants k=230 and k=231, the simplest first-order model obtained by the second algorithm gives negative pressure (of order 10-2). Both reduction algorithms give perfect models of order three; that is, there is no difference between their outputs and the output of the full-order model. Reduced models of order one and two calculated by algorithm 2 give very precise trajectories, whereas the corresponding models found by algorithm 1 have bigger errors but they are still acceptable.Figure 5
The hydrogen pressure trajectories calculated by reduced order (npH2=1,2,3) and full order (npH2=4) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)


### 6.2. Reduction of the Oxygen Pressure Model

The diagonal of the balanced Grammian of the state-space realisation corresponding to the discrete-time difference equation (21) describing oxygen pressure is(43)g=2.3217×10-16.1788×10-51.2064×10-20T. Although the last two diagonal entries of the Grammian are small in comparison with the first ones, which suggests that the last two state variables may be reduced, the reduced models of the first and the second orders are considered. The parameters of the reduced-order oxygen pressure model of order npH2=1,2 calculated by the first algorithm are given in Table 7, whereas the results of the second algorithm are given in Table 8. The first algorithm finds the models of a general structure the same as the full-order model (see (21)), whereas the reduced model of the first order determined by the second method is characterised by (37). In this case, the coefficient b02,2 may be necessary, which manifests a direct influence of the disturbance Iin on the model output pO2. For the second-order reduced model, the additional coefficient is not used, but it is necessary for the first-order model to modify its steady-state gain.Table 7
Parameters of the reduced-order oxygen pressure model of ordernpO2 calculated by the first algorithm.
n p O 2 Parameters 1 a 1 2 = - 7.091825 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 9.182485 × 1 0 - 5 2 a 1 2 = - 9.956871 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 5.160259 × 1 0 - 5 a 2 2 = 2.031841 × 1 0 - 1 b 2 2,1 = - 3.306375 × 1 0 - 2 b 2 2,2 = - 3.040813 × 1 0 - 5Table 8
Parameters of the reduced-order oxygen pressure model of ordernpO2 calculated by the second algorithm.
n p O 2 Parameters 1   b 0 2,2 = 5.637353 × 1 0 - 5 a 1 2 = - 7.091825 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 1.713157 × 1 0 - 4 2 a 1 2 = - 9.956871 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 5.160259 × 1 0 - 5 a 2 2 = 2.031841 × 1 0 - 1 b 2 2,1 = - 3.306375 × 1 0 - 2 b 2 2,2 = - 3.040813 × 1 0 - 5Figure6 depicts the Nyquist plots of reduced-order and full-order models of oxygen pressure. Both reduction methods give very precise models for the input qO2in to output channel; for all examined model orders, Nyquist plots are practically the same. Very important differences are present for the disturbance Iin to output channel. For the first-order model obtained by the first algorithm, the steady-state gain is different from that of the full-order model, but this discrepancy is not present for the second-order structure. For the second algorithm, modelling of the steady-state gain is excellent, but the first-order model gives wrong results as frequency increases. Similar to the first algorithm, this discrepancy is not present for the second-order structure.Figure 6
The Nyquist plots for the oxygen pressure models of reduced order (npO2=1,2) and full order (npO2=3): reduction algorithm 1 (a) and reduction algorithm 2 (b).
(a) (b)For the oxygen pressure approximator, the dynamic model error,EpO2dyn, and steady-state model error, EpO2ss, are defined similar to what is done in the case of hydrogen pressure (see (40) and (41), resp.). Errors of the reduced models of hydrogen pressure are given in Table 9. As suggested by the Nyquist plots, the second algorithm gives much better results (lower dynamic and steady-state errors). For the second algorithm, the steady-state error is 0. Both reduction algorithms give perfect models of order two.Table 9
The dynamic errorEpO2dyn and the steady-state error EpO2ss of the reduced oxygen pressure models (npO2=1,2).
npO2 E p O 2 d y n EpO2ss Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 5.5587×10-4 9.5827×10-8 1.1166×10-3 0.0 2 0.0 0.0 0.0 0.0Finally, Figure7 compares the hydrogen pressure trajectories calculated by reduced-order (npO2=1,2) and full-order (npO2=3) models determined by both algorithms. As expected, taking into account the Nyquist plots and model errors, the second algorithm gives much better results for the first-order model; that is, the model obtained by the first algorithm does not have the correct steady-state properties, but this error is not significant. All models of the second order are perfect.Figure 7
The oxygen pressure trajectories calculated by reduced order (npO2=1,2) and full order (npO2=3) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)


### 6.3. Reduction of the Water Pressure Model

The diagonal of the balanced Grammian of the state-space realisation corresponding to the discrete-time difference equation (22) describing water pressure is(44)g=3.6014×10-35.4358×10-5T.Although there is some difference in the order of the diagonal elements of the balanced Grammian, it is not so huge as in the case of the two previously considered pressure models. For the water pressure model, the only option is to reduce the second state variable. The parameters of the reduced-order hydrogen pressure model of order npH2O=1 calculated by the first algorithm are given in Table 10, whereas the results of the second algorithm are given in Table 11. The general form of the reduced models calculated by the first algorithm is the same as that of the model of the full order (see (22)), whereas the reduced model determined by the second method is characterised by (38). In the latter case, the coefficient b03 is necessary, which manifests a direct influence of the disturbance Iin on the model output pH2O.Table 10
Parameters of the reduced-order water pressure model of ordernpH2O=1 calculated by the first algorithm.
n p H 2 O Parameters 1 a 1 3 = - 9.877933 × 1 0 - 1 b 1 3 = 8.736752 × 1 0 - 5Table 11
Parameters of the reduced-order water pressure model of ordernpH2O calculated by the second algorithm.
n p H 2 O Parameters 1 a 1 3 = - 9.877933 × 1 0 - 1 b 1 3 = 1.580941 × 1 0 - 4 b 0 3 = - 6.819480 × 1 0 - 5Figure8 depicts the Nyquist plots of reduced-order and full-order models. As it is observed in the case of hydrogen and oxygen models, the first algorithm gives slightly worse results for all frequencies. In particular, it does not guarantee that the steady-state gain of the reduced model is exactly the same as that of the original one. The second algorithm leads to a very good model.Figure 8
The Nyquist plots for the water pressure model of reduced order (npH2O=1) and full order (npH2O=2): reduction algorithm 1 (a) and reduction algorithm 2 (b).
(a) (b)For the water pressure approximator, the dynamic model error,EpH2Odyn, and the steady-state model error, EpH2Oss, are defined similar to what is done in the case of hydrogen pressure (see (40) and (41), resp.). Errors of the reduced model of water pressure are given in Table 12. As observed from the Nyquist plots, the model obtained by the second algorithm is characterised by much lower errors. For the second algorithm, the steady-state error is 0, but due to a relatively low dynamic order of the original model (only two), the dynamic error is greater than 0, although it is very small.Table 12
The dynamic errorEpH2Odyn and the steady-state error EpH2Oss of the reduced water pressure model (npH2O=1).
npH2O E p H 2 O d y n EpH2Oss Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 4.0891×10-4 7.4056×10-8 8.2667×10-4 0.0Finally, Figure9 compares the water pressure trajectories calculated by reduced-order (npH2O=1) and full-order (npH2O=2) models determined by both algorithms. As expected, taking into account the Nyquist plots and model errors, the second algorithm gives excellent results; the model generated by the first one does not have correct steady-state properties, but the error is relatively small.Figure 9
The water pressure trajectories calculated by reduced order (npH2O=1) and full order (npH2O=2) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)


### 6.4. Reduction of the Voltage and the Pressure Difference Models

In the previous sections, hydrogen, oxygen, and water pressure models have been reduced separately. The pressures are intermediate model variables; the actual outputs variables are the voltageVr and the pressure difference pdif. Hence, it is necessary to select order of the pressure models in such a way that not only pressures but also outputs are approximated precisely. For this purpose, the dynamic and steady-state errors of both outputs are defined similar to what is done in the case of hydrogen pressure (see (40) and (41), resp.). For the output Vr, the dynamic and steady-state errors are denoted by EVrdyn and EVrss, respectively. Next, for the output pdif, the errors are denoted by Epdifdyn and Epdifss. Table 13 gives values of the dynamic and steady-state errors of the voltage for all possible combinations of orders of hydrogen, oxygen, and water pressure models (npH2, npO2, and npH2O). As suggested by the partial results obtained for reduced models of the pressures, also when the errors of the output variable Vr are taken into account, the second algorithm gives much precise models. In practice, all models obtained by the second algorithm have the steady-state error equal to 0. Considering the dynamic error, three models are selected for further analysis. The perfect model (all errors are 0) that has the lowest number of parameters is characterised by the orders npH2=3, npO2=2, and npH2O=2. Furthermore, two suboptimal models (their errors are greater than 0) are chosen: in the first case, npH2=2, npO2=1, and npH2O=1; in the second case, npH2=npO2=npH2O=1. The first suboptimal model offers a reasonable compromise between accuracy complexities; the second one has the minimal order. Figure 10 compares dynamic trajectories of the voltage Vr calculated by three reduced-order models and the full-order one in case of two compared algorithms. In general, all models work correctly. The perfect reduced model gives no errors, whereas the suboptimal structures lead to small errors. From the comparison of the voltage trajectories, similar to the observations made so far for pressure models, it is clear that the second algorithm gives much more precise models.Table 13
The influence of the order of hydrogen, oxygen, and water pressure models (npH2, npO2, and npH2O) on the dynamic error EVrdyn and the steady-state error EVrss of the output variable Vr.
npH2 n p O 2 n p H 2 O E V r d y n E V r s s Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 1 1 1.7031 × 1 0 1 1.3307 1.2006 × 1 0 2 6.8080 × 1 0 - 10 1 1 2 1.7026 × 1 0 1 1.3307 1.198117 × 1 0 2 6.8080 × 1 0 - 10 1 2 1 1.6021 × 1 0 1 1.3247 1.084092 × 1 0 2 0.0 1 2 2 1.5278 × 1 0 2 1.3248 1.073773 × 1 0 2 0.0 1 3 1 1.6021 × 1 0 1 1.3247 1.084092 × 1 0 2 0.0 1 3 2 1.5278 × 1 0 2 1.3248 1.073773 × 1 0 2 0.0 2 1 1 4.4605 4.5855 × 1 0 - 3 6.2967 × 1 0 1 6.8380 × 1 0 - 10 2 1 2 4.8173 4.3832 × 1 0 - 3 6.3211 × 1 0 1 6.8380 × 1 0 - 10 2 2 1 1.9113 5.7115 × 1 0 - 4 4.5223 × 1 0 1 0.0 2 2 2 1.5297 5.7115 × 1 0 - 4 4.4634 × 1 0 1 0.0 2 3 1 1.5297 6.4325 × 1 0 - 4 4.5223 × 1 0 1 0.0 2 3 2 1.5297 5.7115 × 1 0 - 4 4.4634 × 1 0 1 0.0 3 1 1 5.6368 2.8609 × 1 0 - 3 2.3427 × 1 0 1 6.8080 × 1 0 - 10 3 1 2 6.3333 2.7232 × 1 0 - 3 2.4225 × 1 0 1 6.8080 × 1 0 - 10 3 2 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0 3 2 2 0.0 0.0 0.0 0.0 3 3 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0 3 3 2 0.0 0.0 0.0 0.0 4 1 1 5.6368 2.8609 × 1 0 - 3 2.3427 × 1 0 1 6.8080 × 1 0 - 10 4 1 2 6.3333 2.7232 × 1 0 - 3 2.4225 × 1 0 1 6.8080 × 1 0 - 10 4 2 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0 4 2 2 0.0 0.0 0.0 0.0 4 3 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0Figure 10
The voltageVr trajectories calculated by three reduced-order and full-order (npH2=4, npO2=3, and npH2O=2) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)Next, the second output variable is taken into account. Table14 gives values of the dynamic and steady-state errors of the pressure difference for all possible combinations of orders of hydrogen and pressure models (npH2 and npO2); the order of the water pressure model does not influence the pressure difference. All the errors are very small, since modelling the pressure difference is much easier than modelling of the voltage. It is because the Nernst equation (25) is nonlinear, whereas the pressure difference depends linearly only on hydrogen and oxygen pressures (see (28)). The second observation is that the second algorithm in comparison with the first one gives smaller errors. In particular, all steady-state errors of the model generated by the second method are equal to 0. The same three reduced-order models are selected taking into account the errors of the pressure difference: the perfect reduced structure (with all errors equal to 0) and two suboptimal ones. Figure 11 compares dynamic trajectories of the pressure difference pdif calculated by three reduced-order models and the full-order one in case of two compared algorithms. As in the case of the voltage model, also all pressure difference models work correctly. The perfect reduced model gives no errors, whereas the suboptimal reduced structures lead to small errors. As always, in this study, the second algorithm gives more precise models.Table 14
The influence of the order of hydrogen and oxygen pressure models (npH2 and npO2) on the dynamic error Epdifdyn and the steady-state error Epdifss of the output variable pdif.
n p H 2 n p O 2 E p d i f d y n E p d i f s s Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 1 4.4275 × 1 0 - 3 1.1629 × 1 0 - 4 2.3936 × 1 0 - 2 0.0 1 2 2.5494 × 1 0 - 3 1.1652 × 1 0 - 4 1.7068 × 1 0 - 2 0.0 1 3 2.5494 × 1 0 - 3 1.1652 × 1 0 - 4 1.7068 × 1 0 - 2 0.0 2 1 1.4349 × 1 0 - 3 2.1408 × 1 0 - 8 2.8942 × 1 0 - 3 0.0 2 2 2.0677 × 1 0 - 4 3.5651 × 1 0 - 8 4.1543 × 1 0 - 4 0.0 2 3 2.0677 × 1 0 - 4 3.5651 × 1 0 - 8 4.1543 × 1 0 - 4 0.0 3 1 5.5417 × 1 0 - 4 9.5827 × 1 0 - 8 1.1166 × 1 0 - 3 0.0 3 2 0.0 0.0 0.0 0.0 3 3 0.0 0.0 0.0 0.0 4 1 5.5417 × 1 0 - 4 9.5827 × 1 0 - 8 1.1166 × 1 0 - 3 0.0 4 2 0.0 0.0 0.0 0.0Figure 11
The pressure differencepdif trajectories calculated by three reduced-order and full-order (npH2=4 and npO2=3) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)To show good steady-state properties of the simplest reduced model (npH2=npO2=npH2O=1), its steady-state characteristics obtained by two considered algorithms are compared with the ideal characteristics of the full-order model in Figure 12. When the voltage is considered, there are very small differences between the ideal characteristic and that of the reduced model generated by algorithm 1. The characteristic of the model produced by the second algorithm is very precise, as demonstrated numerically in Table 13. When the pressure difference is considered, the differences between the model obtained by the first algorithm and the full-order one are practically not visible; the second algorithm gives perfect characteristic, as demonstrated numerically in Table 14.Figure 12
The steady-state voltageVr and partial difference pdif steady-state characteristics calculated by the full-order (npH2=4, npO2=3, and npH2O=2) model (a) and the simplest reduced model (npH2=npO2=npH2O=1) (b, c).
(a) (b) (c)Finally, Table15 gives a comparison of dynamic and steady-state model errors for the ideal full-order model and the chosen three reduced ones. The errors for two output variables are given and only the second algorithm is considered, since it always gives better results than the first one. The number of pressure model parameters is given. Although the perfect reduced-order model gives all errors equal to 0, it only reduces the overall number of parameters from 25 to 19. That is why in practical applications it may be more attractive to use suboptimal structures. The first suboptimal reduced model makes it possible to reduce the number of parameters to 14 but at the same time its dynamic and steady-state behaviours are practically the same as that of the full-order model. It is interesting to point out that even the simplest first-order reduced model, which makes it possible to reduce the number of parameters to 11, gives quite good results. Of course, selection of the models depends on specifications of the considered application.Table 15
Comparison of the ideal full-order model (npH2= 4, npO2=3, and npH2O=2) and the chosen reduced ones calculated by the second algorithms in terms of the number of pressure parameters as well as the dynamic and steady-state errors for two process output variables (algorithm 2).
n p H 2 n p O 2 n p H 2 O Parameters E V r d y n E V r s s E p d i f d y n E p d i f s s 4 3 2 25 0.0 0.0 0.0 0.0 3 2 2 19 0.0 0.0 0.0 0.0 2 1 1 14 4.3832 × 1 0 - 3 6.8380 × 1 0 - 10 2.1408 × 1 0 - 8 0.0 1 1 1 11 1.3307 6.8080 × 1 0 - 10 1.1629 × 1 0 - 4 0.0


## 6.1. Reduction of the Hydrogen Pressure Model

The diagonal of the balanced Grammian of the state-space realisation corresponding to the discrete-time difference equation (20) describing hydrogen pressure is(39)g=6.7211×10-17.8946×10-23.7688×10-51.3656×10-17. The last two diagonal entries of the Grammian are small in comparison with the first two ones, which suggests that the last two state variables may be reduced. However, for a thorough analysis, the reduced models of the first, the second, and the third orders are considered; in each case, two reduction algorithms are used. The parameters of the reduced-order hydrogen pressure model of order npH2=1,2,3 calculated by the first algorithm are given in Table 4, whereas the results of the second algorithm are given in Table 5. The general form of the reduced models calculated by the first algorithm is the same as that of the full-order model (see (20)), whereas the reduced model determined by the second method is characterised by (36). In the latter case, the coefficients b01,1 and b01,2 may be necessary, which manifest a direct influence of the model input qfin and the disturbance Iin on the model output pH2. When the reduced model has the third order, the additional coefficients are not necessary to guarantee that the steady-state gain of the reduced model is the same as that of the full-order one. When the reduced model has the second order, only one of them is used, but in the case of the first-order model two of them are necessary.Table 4
Parameters of the reduced-order hydrogen pressure model of ordernpH2 calculated by the first algorithm.
n p H 2 Parameters 1 a 1 1 = - 9.768754 × 1 0 - 1 b 1 1,1 = 3.054446 × 1 0 - 2 b 1 1,2 = - 6.812056 × 1 0 - 5  2 a 1 1 = - 1.781141 b 1 1,1 = 4.201816 × 1 0 - 3 b 1 1,2 = - 6.557908 × 1 0 - 5 a 2 1 = 7.879551 × 1 0 - 1 b 2 1,1 = 3.881002 × 1 0 - 3 b 2 1,2 = 4.914773 × 1 0 - 5  3 a 1 1 = - 2.067646 b 1 1,1 = 4.201819 × 1 0 - 3 b 1 1,2 = - 3.832050 × 1 0 - 5 a 2 1 = 1.298261 b 2 1,1 = 2.677158 × 1 0 - 3 b 2 1,2 = 6.319519 × 1 0 - 6 a 3 1 = - 2.257529 × 1 0 - 1 b 3 1,1 = - 1.111925 × 1 0 - 3 b 3 1,2 = 2.051302 × 1 0 - 5Table 5
Parameters of the reduced-order hydrogen pressure model of ordernpH2 calculated by the second algorithm.
n p H 2 Parameters 1  b 0 1,1 = - 1.345838 × 1 0 - 1 b 0 1,2 = - 3.381314 × 1 0 - 5 a 1 1 = - 9.651884 × 1 0 - 1 b 1 1,1 = 1.758786 × 1 0 - 1 b 1 1,2 = - 4.844625 × 1 0 - 5  2   b 0 1,2 = 3.820428 × 1 0 - 5 a 1 1 = - 1.781141 b 1 1,1 = 4.201816 × 1 0 - 3 b 1 1,2 = - 1.519383 × 1 0 - 4 a 2 1 = 7.879551 × 1 0 - 1 b 2 1,1 = 3.881002 × 1 0 - 3 b 2 1,2 = 9.763310 × 1 0 - 5  3 a 1 1 = - 2.067646 b 1 1,1 = 4.201819 × 1 0 - 3 b 1 1,2 = - 3.832050 × 1 0 - 5 a 2 1 = 1.298261 b 2 1,1 = 2.677158 × 1 0 - 3 b 2 1,2 = 6.319519 × 1 0 - 6 a 3 1 = - 2.257529 × 1 0 - 1 b 3 1,1 = - 1.111925 × 1 0 - 3 b 3 1,2 = 2.051302 × 1 0 - 5To compare full-order and reduced models of the hydrogen pressure, it is convenient to consider their Nyquist plots. It is because Nyquist plots show both steady-state and dynamic behaviours of models. Figure4 depicts the Nyquist plots of reduced-order and full-order models. As expected, the first algorithm does not guarantee that the stead-state gain of the reduced model is exactly the same as that of the original one (the gain may be read for ω=0). Furthermore, for other frequencies, the Nyquist plot of the reduced models is different from that of the full-order one. This observation is true for the simplest first-order model and, for the disturbance Iin to output channel for the second-order structure, increasing the model order makes it possible to obtain precise steady-state gain. The second algorithm is much more efficient, since all reduced models, even the first-order one, have the same gain as the original one. Furthermore, the Nyquist plots of the reduced models are very similar to those of the full-order ones.Figure 4
The Nyquist plots for the hydrogen pressure models of reduced order (npH2=1,2,3) and full order (npH2=4): reduction algorithm 1 (a) and reduction algorithm 2 (b).
(a) (b)The models are also compared in a quantitative way. For this purpose, the dynamic model error is defined:(40)EpH2dyn=1ndyn∑k=1ndynpH2fullk-pH2reducedk2.The outputs of the full-order and reduced models for the sampling instant k are denoted by pH2full(k) and pH2reduced(k), respectively. Both models are excited by the same series of step changes of the manipulated variables and steps of the disturbance depicted in Figure 2; that is, ndyn=2000. Additionally, the steady-state model error is defined:(41)EpH2ss=1nss∑p=1nsspH2fullp-pH2reducedp2.In this case, pH2full(p) and pH2reduced(p) denote outputs of the full-order and reduced steady-state models for the data point p, respectively. They are derived easily from the dynamic models. For the second algorithm, using (36), the dynamic model for hydrogen pressure is (42)pH2=∑i=0npH2bi1,1qfin+∑i=0npH2bi1,2Iin1+∑i=1npH2ai1-1.In the case of the reduced model obtained by the first algorithm (see (20)), b01,1=b01,2=0. For the steady-state models, nss=4116. Dynamic and steady-state errors of the reduced models of hydrogen pressure are given in Table 6. As observed from the Nyquist plots, the second algorithm gives much better results (lower model errors). In particular, always EpH2ss=0 for the second algorithm, since it always guarantees that the steady-state properties of the original and reduced models are the same. When the errors are lower than the machine accuracy (2.2204×10-16), they are treated as 0. Both reduction algorithms give perfect models of order three.Table 6
The dynamic errorEpH2dyn and the steady-state error EpH2ss of the reduced hydrogen pressure models (npH2=1,2,3).
n p H 2 E p H 2 d y n E p H 2 s s Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 2.5687 × 1 0 - 3 1.1652 × 1 0 - 4 1.7068 × 1 0 - 2 0.0 2 2.0740 × 1 0 - 4 3.5651 × 1 0 - 8 4.1544 × 1 0 - 4 0.0 3 0.0 0.0 0.0 0.0Finally, Figure5 compares the hydrogen pressure trajectories (model outputs) calculated by reduced-order (npH2=1,2,3) and full-order (npH2=4) models determined by both algorithms; all models are excited by the same series of step changes of the manipulated variables and steps of the disturbance depicted in Figure 2. As expected, taking into account the Nyquist plots and model errors, the second algorithm gives much better results for the same model order. For the sampling instants k=230 and k=231, the simplest first-order model obtained by the second algorithm gives negative pressure (of order 10-2). Both reduction algorithms give perfect models of order three; that is, there is no difference between their outputs and the output of the full-order model. Reduced models of order one and two calculated by algorithm 2 give very precise trajectories, whereas the corresponding models found by algorithm 1 have bigger errors but they are still acceptable.Figure 5
The hydrogen pressure trajectories calculated by reduced order (npH2=1,2,3) and full order (npH2=4) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)


## 6.2. Reduction of the Oxygen Pressure Model

The diagonal of the balanced Grammian of the state-space realisation corresponding to the discrete-time difference equation (21) describing oxygen pressure is(43)g=2.3217×10-16.1788×10-51.2064×10-20T. Although the last two diagonal entries of the Grammian are small in comparison with the first ones, which suggests that the last two state variables may be reduced, the reduced models of the first and the second orders are considered. The parameters of the reduced-order oxygen pressure model of order npH2=1,2 calculated by the first algorithm are given in Table 7, whereas the results of the second algorithm are given in Table 8. The first algorithm finds the models of a general structure the same as the full-order model (see (21)), whereas the reduced model of the first order determined by the second method is characterised by (37). In this case, the coefficient b02,2 may be necessary, which manifests a direct influence of the disturbance Iin on the model output pO2. For the second-order reduced model, the additional coefficient is not used, but it is necessary for the first-order model to modify its steady-state gain.Table 7
Parameters of the reduced-order oxygen pressure model of ordernpO2 calculated by the first algorithm.
n p O 2 Parameters 1 a 1 2 = - 7.091825 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 9.182485 × 1 0 - 5 2 a 1 2 = - 9.956871 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 5.160259 × 1 0 - 5 a 2 2 = 2.031841 × 1 0 - 1 b 2 2,1 = - 3.306375 × 1 0 - 2 b 2 2,2 = - 3.040813 × 1 0 - 5Table 8
Parameters of the reduced-order oxygen pressure model of ordernpO2 calculated by the second algorithm.
n p O 2 Parameters 1   b 0 2,2 = 5.637353 × 1 0 - 5 a 1 2 = - 7.091825 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 1.713157 × 1 0 - 4 2 a 1 2 = - 9.956871 × 1 0 - 1 b 1 2,1 = 1.154038 × 1 0 - 1 b 1 2,2 = - 5.160259 × 1 0 - 5 a 2 2 = 2.031841 × 1 0 - 1 b 2 2,1 = - 3.306375 × 1 0 - 2 b 2 2,2 = - 3.040813 × 1 0 - 5Figure6 depicts the Nyquist plots of reduced-order and full-order models of oxygen pressure. Both reduction methods give very precise models for the input qO2in to output channel; for all examined model orders, Nyquist plots are practically the same. Very important differences are present for the disturbance Iin to output channel. For the first-order model obtained by the first algorithm, the steady-state gain is different from that of the full-order model, but this discrepancy is not present for the second-order structure. For the second algorithm, modelling of the steady-state gain is excellent, but the first-order model gives wrong results as frequency increases. Similar to the first algorithm, this discrepancy is not present for the second-order structure.Figure 6
The Nyquist plots for the oxygen pressure models of reduced order (npO2=1,2) and full order (npO2=3): reduction algorithm 1 (a) and reduction algorithm 2 (b).
(a) (b)For the oxygen pressure approximator, the dynamic model error,EpO2dyn, and steady-state model error, EpO2ss, are defined similar to what is done in the case of hydrogen pressure (see (40) and (41), resp.). Errors of the reduced models of hydrogen pressure are given in Table 9. As suggested by the Nyquist plots, the second algorithm gives much better results (lower dynamic and steady-state errors). For the second algorithm, the steady-state error is 0. Both reduction algorithms give perfect models of order two.Table 9
The dynamic errorEpO2dyn and the steady-state error EpO2ss of the reduced oxygen pressure models (npO2=1,2).
npO2 E p O 2 d y n EpO2ss Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 5.5587×10-4 9.5827×10-8 1.1166×10-3 0.0 2 0.0 0.0 0.0 0.0Finally, Figure7 compares the hydrogen pressure trajectories calculated by reduced-order (npO2=1,2) and full-order (npO2=3) models determined by both algorithms. As expected, taking into account the Nyquist plots and model errors, the second algorithm gives much better results for the first-order model; that is, the model obtained by the first algorithm does not have the correct steady-state properties, but this error is not significant. All models of the second order are perfect.Figure 7
The oxygen pressure trajectories calculated by reduced order (npO2=1,2) and full order (npO2=3) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)


## 6.3. Reduction of the Water Pressure Model

The diagonal of the balanced Grammian of the state-space realisation corresponding to the discrete-time difference equation (22) describing water pressure is(44)g=3.6014×10-35.4358×10-5T.Although there is some difference in the order of the diagonal elements of the balanced Grammian, it is not so huge as in the case of the two previously considered pressure models. For the water pressure model, the only option is to reduce the second state variable. The parameters of the reduced-order hydrogen pressure model of order npH2O=1 calculated by the first algorithm are given in Table 10, whereas the results of the second algorithm are given in Table 11. The general form of the reduced models calculated by the first algorithm is the same as that of the model of the full order (see (22)), whereas the reduced model determined by the second method is characterised by (38). In the latter case, the coefficient b03 is necessary, which manifests a direct influence of the disturbance Iin on the model output pH2O.Table 10
Parameters of the reduced-order water pressure model of ordernpH2O=1 calculated by the first algorithm.
n p H 2 O Parameters 1 a 1 3 = - 9.877933 × 1 0 - 1 b 1 3 = 8.736752 × 1 0 - 5Table 11
Parameters of the reduced-order water pressure model of ordernpH2O calculated by the second algorithm.
n p H 2 O Parameters 1 a 1 3 = - 9.877933 × 1 0 - 1 b 1 3 = 1.580941 × 1 0 - 4 b 0 3 = - 6.819480 × 1 0 - 5Figure8 depicts the Nyquist plots of reduced-order and full-order models. As it is observed in the case of hydrogen and oxygen models, the first algorithm gives slightly worse results for all frequencies. In particular, it does not guarantee that the steady-state gain of the reduced model is exactly the same as that of the original one. The second algorithm leads to a very good model.Figure 8
The Nyquist plots for the water pressure model of reduced order (npH2O=1) and full order (npH2O=2): reduction algorithm 1 (a) and reduction algorithm 2 (b).
(a) (b)For the water pressure approximator, the dynamic model error,EpH2Odyn, and the steady-state model error, EpH2Oss, are defined similar to what is done in the case of hydrogen pressure (see (40) and (41), resp.). Errors of the reduced model of water pressure are given in Table 12. As observed from the Nyquist plots, the model obtained by the second algorithm is characterised by much lower errors. For the second algorithm, the steady-state error is 0, but due to a relatively low dynamic order of the original model (only two), the dynamic error is greater than 0, although it is very small.Table 12
The dynamic errorEpH2Odyn and the steady-state error EpH2Oss of the reduced water pressure model (npH2O=1).
npH2O E p H 2 O d y n EpH2Oss Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 4.0891×10-4 7.4056×10-8 8.2667×10-4 0.0Finally, Figure9 compares the water pressure trajectories calculated by reduced-order (npH2O=1) and full-order (npH2O=2) models determined by both algorithms. As expected, taking into account the Nyquist plots and model errors, the second algorithm gives excellent results; the model generated by the first one does not have correct steady-state properties, but the error is relatively small.Figure 9
The water pressure trajectories calculated by reduced order (npH2O=1) and full order (npH2O=2) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)


## 6.4. Reduction of the Voltage and the Pressure Difference Models

In the previous sections, hydrogen, oxygen, and water pressure models have been reduced separately. The pressures are intermediate model variables; the actual outputs variables are the voltageVr and the pressure difference pdif. Hence, it is necessary to select order of the pressure models in such a way that not only pressures but also outputs are approximated precisely. For this purpose, the dynamic and steady-state errors of both outputs are defined similar to what is done in the case of hydrogen pressure (see (40) and (41), resp.). For the output Vr, the dynamic and steady-state errors are denoted by EVrdyn and EVrss, respectively. Next, for the output pdif, the errors are denoted by Epdifdyn and Epdifss. Table 13 gives values of the dynamic and steady-state errors of the voltage for all possible combinations of orders of hydrogen, oxygen, and water pressure models (npH2, npO2, and npH2O). As suggested by the partial results obtained for reduced models of the pressures, also when the errors of the output variable Vr are taken into account, the second algorithm gives much precise models. In practice, all models obtained by the second algorithm have the steady-state error equal to 0. Considering the dynamic error, three models are selected for further analysis. The perfect model (all errors are 0) that has the lowest number of parameters is characterised by the orders npH2=3, npO2=2, and npH2O=2. Furthermore, two suboptimal models (their errors are greater than 0) are chosen: in the first case, npH2=2, npO2=1, and npH2O=1; in the second case, npH2=npO2=npH2O=1. The first suboptimal model offers a reasonable compromise between accuracy complexities; the second one has the minimal order. Figure 10 compares dynamic trajectories of the voltage Vr calculated by three reduced-order models and the full-order one in case of two compared algorithms. In general, all models work correctly. The perfect reduced model gives no errors, whereas the suboptimal structures lead to small errors. From the comparison of the voltage trajectories, similar to the observations made so far for pressure models, it is clear that the second algorithm gives much more precise models.Table 13
The influence of the order of hydrogen, oxygen, and water pressure models (npH2, npO2, and npH2O) on the dynamic error EVrdyn and the steady-state error EVrss of the output variable Vr.
npH2 n p O 2 n p H 2 O E V r d y n E V r s s Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 1 1 1.7031 × 1 0 1 1.3307 1.2006 × 1 0 2 6.8080 × 1 0 - 10 1 1 2 1.7026 × 1 0 1 1.3307 1.198117 × 1 0 2 6.8080 × 1 0 - 10 1 2 1 1.6021 × 1 0 1 1.3247 1.084092 × 1 0 2 0.0 1 2 2 1.5278 × 1 0 2 1.3248 1.073773 × 1 0 2 0.0 1 3 1 1.6021 × 1 0 1 1.3247 1.084092 × 1 0 2 0.0 1 3 2 1.5278 × 1 0 2 1.3248 1.073773 × 1 0 2 0.0 2 1 1 4.4605 4.5855 × 1 0 - 3 6.2967 × 1 0 1 6.8380 × 1 0 - 10 2 1 2 4.8173 4.3832 × 1 0 - 3 6.3211 × 1 0 1 6.8380 × 1 0 - 10 2 2 1 1.9113 5.7115 × 1 0 - 4 4.5223 × 1 0 1 0.0 2 2 2 1.5297 5.7115 × 1 0 - 4 4.4634 × 1 0 1 0.0 2 3 1 1.5297 6.4325 × 1 0 - 4 4.5223 × 1 0 1 0.0 2 3 2 1.5297 5.7115 × 1 0 - 4 4.4634 × 1 0 1 0.0 3 1 1 5.6368 2.8609 × 1 0 - 3 2.3427 × 1 0 1 6.8080 × 1 0 - 10 3 1 2 6.3333 2.7232 × 1 0 - 3 2.4225 × 1 0 1 6.8080 × 1 0 - 10 3 2 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0 3 2 2 0.0 0.0 0.0 0.0 3 3 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0 3 3 2 0.0 0.0 0.0 0.0 4 1 1 5.6368 2.8609 × 1 0 - 3 2.3427 × 1 0 1 6.8080 × 1 0 - 10 4 1 2 6.3333 2.7232 × 1 0 - 3 2.4225 × 1 0 1 6.8080 × 1 0 - 10 4 2 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0 4 2 2 0.0 0.0 0.0 0.0 4 3 1 4.1917 × 1 0 - 2 7.5057 × 1 0 - 6 4.0895 × 1 0 - 2 0.0Figure 10
The voltageVr trajectories calculated by three reduced-order and full-order (npH2=4, npO2=3, and npH2O=2) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)Next, the second output variable is taken into account. Table14 gives values of the dynamic and steady-state errors of the pressure difference for all possible combinations of orders of hydrogen and pressure models (npH2 and npO2); the order of the water pressure model does not influence the pressure difference. All the errors are very small, since modelling the pressure difference is much easier than modelling of the voltage. It is because the Nernst equation (25) is nonlinear, whereas the pressure difference depends linearly only on hydrogen and oxygen pressures (see (28)). The second observation is that the second algorithm in comparison with the first one gives smaller errors. In particular, all steady-state errors of the model generated by the second method are equal to 0. The same three reduced-order models are selected taking into account the errors of the pressure difference: the perfect reduced structure (with all errors equal to 0) and two suboptimal ones. Figure 11 compares dynamic trajectories of the pressure difference pdif calculated by three reduced-order models and the full-order one in case of two compared algorithms. As in the case of the voltage model, also all pressure difference models work correctly. The perfect reduced model gives no errors, whereas the suboptimal reduced structures lead to small errors. As always, in this study, the second algorithm gives more precise models.Table 14
The influence of the order of hydrogen and oxygen pressure models (npH2 and npO2) on the dynamic error Epdifdyn and the steady-state error Epdifss of the output variable pdif.
n p H 2 n p O 2 E p d i f d y n E p d i f s s Algorithm 1 Algorithm 2 Algorithm 1 Algorithm 2 1 1 4.4275 × 1 0 - 3 1.1629 × 1 0 - 4 2.3936 × 1 0 - 2 0.0 1 2 2.5494 × 1 0 - 3 1.1652 × 1 0 - 4 1.7068 × 1 0 - 2 0.0 1 3 2.5494 × 1 0 - 3 1.1652 × 1 0 - 4 1.7068 × 1 0 - 2 0.0 2 1 1.4349 × 1 0 - 3 2.1408 × 1 0 - 8 2.8942 × 1 0 - 3 0.0 2 2 2.0677 × 1 0 - 4 3.5651 × 1 0 - 8 4.1543 × 1 0 - 4 0.0 2 3 2.0677 × 1 0 - 4 3.5651 × 1 0 - 8 4.1543 × 1 0 - 4 0.0 3 1 5.5417 × 1 0 - 4 9.5827 × 1 0 - 8 1.1166 × 1 0 - 3 0.0 3 2 0.0 0.0 0.0 0.0 3 3 0.0 0.0 0.0 0.0 4 1 5.5417 × 1 0 - 4 9.5827 × 1 0 - 8 1.1166 × 1 0 - 3 0.0 4 2 0.0 0.0 0.0 0.0Figure 11
The pressure differencepdif trajectories calculated by three reduced-order and full-order (npH2=4 and npO2=3) models determined by algorithm 1 (a) and algorithm 2 (b).
(a) (b)To show good steady-state properties of the simplest reduced model (npH2=npO2=npH2O=1), its steady-state characteristics obtained by two considered algorithms are compared with the ideal characteristics of the full-order model in Figure 12. When the voltage is considered, there are very small differences between the ideal characteristic and that of the reduced model generated by algorithm 1. The characteristic of the model produced by the second algorithm is very precise, as demonstrated numerically in Table 13. When the pressure difference is considered, the differences between the model obtained by the first algorithm and the full-order one are practically not visible; the second algorithm gives perfect characteristic, as demonstrated numerically in Table 14.Figure 12
The steady-state voltageVr and partial difference pdif steady-state characteristics calculated by the full-order (npH2=4, npO2=3, and npH2O=2) model (a) and the simplest reduced model (npH2=npO2=npH2O=1) (b, c).
(a) (b) (c)Finally, Table15 gives a comparison of dynamic and steady-state model errors for the ideal full-order model and the chosen three reduced ones. The errors for two output variables are given and only the second algorithm is considered, since it always gives better results than the first one. The number of pressure model parameters is given. Although the perfect reduced-order model gives all errors equal to 0, it only reduces the overall number of parameters from 25 to 19. That is why in practical applications it may be more attractive to use suboptimal structures. The first suboptimal reduced model makes it possible to reduce the number of parameters to 14 but at the same time its dynamic and steady-state behaviours are practically the same as that of the full-order model. It is interesting to point out that even the simplest first-order reduced model, which makes it possible to reduce the number of parameters to 11, gives quite good results. Of course, selection of the models depends on specifications of the considered application.Table 15
Comparison of the ideal full-order model (npH2= 4, npO2=3, and npH2O=2) and the chosen reduced ones calculated by the second algorithms in terms of the number of pressure parameters as well as the dynamic and steady-state errors for two process output variables (algorithm 2).
n p H 2 n p O 2 n p H 2 O Parameters E V r d y n E V r s s E p d i f d y n E p d i f s s 4 3 2 25 0.0 0.0 0.0 0.0 3 2 2 19 0.0 0.0 0.0 0.0 2 1 1 14 4.3832 × 1 0 - 3 6.8380 × 1 0 - 10 2.1408 × 1 0 - 8 0.0 1 1 1 11 1.3307 6.8080 × 1 0 - 10 1.1629 × 1 0 - 4 0.0


## 7. Conclusions

This work discusses reduction of the dynamic model of the solid oxide fuel cell. Firstly, the discrete-time representation of the continuous-time model is found. It is shown that the discrete-time model gives process trajectories practically the same as the original continuous-time one. Next, the discrete-time submodels of hydrogen, oxygen, and water pressures (intermediate variables) are reduced by means of state truncation and residualisation methods. The reduced pressure models are assessed taking into account both dynamic and steady-state errors as well as Nyquist plots. Next, the influence of the order of the pressure models on the dynamic and steady-state errors of the process outputs (the voltage and the pressure difference) is studied. Three reduced models are chosen: with 19, 14, and 11 pressure parameters (the original one has 25 of them). In general, in all cases, the second reduction method is more precise.It is important to summarise the advantages of the described reduced models:(1)
The reduced models calculate not only the output variables (as the black-box models) but also the intermediate variables (the pressures).(2)
The reduced models are selected taking into account dynamic and static errors as well as Nyquist plots.(3)
The first reduced model (with 19 pressure parameters) does not have any steady-state and dynamic errors; it is able to calculate pressures and process outputs exactly in the same way it is done by the full-order model.(4)
Although the second and the third reduced models (with 14 and 11 pressure parameters, resp.) have some errors, they are quite small and the reduced models give dynamic trajectories and steady-state characteristics that are very similar to those obtained by the full-order structure.The reduced-order model of the solid oxide fuel cell may be further used in all model-based methods. Firstly, it may be used for design and implementation of MPC algorithms [31], in particular computationally efficient MPC schemes [32]. In addition to classical MPC algorithms, robust versions with guaranteed stability may be considered [33]. Secondly, it may be also used in online process optimisation [31] cooperating with MPC. Thirdly, the model may be used in fault diagnosis of the process and fault-tolerant control [34]. The presented model reduction approach may be also used for different variants of dynamic models of fuel cells [9, 22–24, 35, 36].



---

*Source: 6021249-2018-07-04.xml*