Control-Oriented Modeling of the Cooling Process of a PEMFC-Based $\mu$ -CHP System

Micro-combined heat and power systems (<inline-formula> <tex-math notation="LaTeX">$\mu $ </tex-math></inline-formula>-CHP) based on proton exchange membrane fuel cell stacks (PEMFC) are capable of supplying electricity and heat for the residential housing sector with a high energy efficiency and a low level of <inline-formula> <tex-math notation="LaTeX">$CO_ {2}$ </tex-math></inline-formula> emissions. For this reason, they are regarded as a promising technology for coping with the current environmental challenges. In these systems, the temperature control of the stack is crucial, since it has a direct impact on its durability and electrical efficiency. In order to design a good temperature control, however, a dynamic model of the <inline-formula> <tex-math notation="LaTeX">$\mu $ </tex-math></inline-formula>-CHP cooling system is required. In this paper, we present a model of the cooling system of a PEMFC-based <inline-formula> <tex-math notation="LaTeX">$\mu $ </tex-math></inline-formula>-CHP system, which is oriented to the design of the temperature control of the stack. The model has been developed from a <inline-formula> <tex-math notation="LaTeX">$\mu $ </tex-math></inline-formula>-CHP system located in the laboratory of our research team, the predictive control and heuristic optimization group (CPOH). It is based on first principles, dynamic, non-linear, and has been validated against the experimental data. The model is implemented in Matlab/Simulink and the adjustment of its parameters was carried out using evolutionary optimization techniques. The methodology followed to obtain it is also described in detail. Both the model and the test data used for its adjustment and validation are accessible to anyone who wants to consult them. The results show that the model is able to faithfully represent the dynamics of the <inline-formula> <tex-math notation="LaTeX">$\mu $ </tex-math></inline-formula>-CHP cooling system, so it is appropriate for the design of the stack temperature control.


c w
Specific heat of water, C, J /(kg·K ). T ts Logarithmic mean temperature difference in the heat exchanger, V, K .

T ts in
Temperature difference between heat exchanger inlets (tube-shell), V, K . T ts out Temperature difference between heat exchanger outlets (tube-shell),V, K . F a Air flow rate, V, m 3 /s.

F a max
Maximum value of F a , C, l/min. F a min Minimum value of F a , C, l/min. F w 1 Water flow rate (primary circuit), V, m 3 /s. F w 1max Maximum value of F w 1 , C, l/min. F w 1min Minimum value of F w 1 , C, l/min. F w 2 Water flow rate (secondary circuit), V, m 3 /s. F w 2max Maximum value of F w 2 , C, l/min. F w 2min Minimum value of F w 2 , C, l/min. h a Air side heat transfer parameter in the stack, V, W /K . h a max Value of h a at F a max , P, W /K . h a min Value of h a at F a min , P, W /K . h aw Water-air interface heat transfer parameter in the stack, P, W /K . h e Conduction heat transfer parameter in the heat exchanger (wall), C, W /K . h fc loss Stack-ambient heat transfer parameter, P, W /K . h fc 1 Overall water-air heat transfer parameter in the stack, V, W /K . h fc 2 Water-housing heat transfer parameter in the stack, V, W /K . h fc 2max Value of h fc 2 at F w 1max , P, W /K . h fc 2min Value of h fc 2  Product ρ a c a , P, J /(m 3 ·K ). k w Product ρ w c w , C, J /(m 3 ·K ). n c Number of cells in the stack, C, ud. n ch Number of channels in the heat exchanger, C, ud. Q a Heat difference (stack air), V, W . Q fc loss Heat flow rate (stack-ambient), V, W .

Q heat
Heat flow rate generated by the reaction in the stack, V, W .

Q s
Heat difference (heat exchanger shell side), V, W . Q t Heat difference (heat exchanger tube side), V, W .

Q ts
Heat flow rate (tube-shell), V, W .

Q wa
Heat flow rate (water-air, in the stack), V, W .

Q wfc
Heat flow rate (water-housing, in the stack), V, W .

T a
Average air temperature in the stack, V, K .

T a in
Air temperature (stack air inlet), V, K . T amb Ambient temperature, V, K . T amb r Ambient temperature for the radiator, P, K .

T a out
Air temperature (stack air outlet), V, K . T fc Stack housing temperature, V, K .

T p1
Average water temperature in pipe 1, V, K .

T p4
Average water temperature in pipe 4, V, K . T p4 in Water temperature (pipe 4 inlet), V, K . T p4 out Water temperature (pipe 4 outlet), V, K . T r Average water temperature in radiator, V, K .

T r in
Water temperature (radiator inlet), V, K . T r out Water temperature (radiator outlet), V, K .

T s in
Water temperature (heat exchanger shell inlet), V, K .

T s out
Water temperature (heat exchanger shell outlet), V, K .

T t in
Water temperature (heat exchanger tube inlet), V, K .

T t out
Water temperature (heat exchanger tube outlet), V, K .

T t2 in
Water temperature (tank 2 inlet), V, K . T w Average water temperature in the stack, V, K .

T w in
Water temperature (stack water inlet), V, K . T w out Water temperature (stack water outlet), V, K . v Stack voltage, V, V . V a Volume of air in the stack, P, m 3 .
V p3 Volume of pipe 3, C, m 3 . V p4 Volume of pipe 4, P, m 3 . V r Volume of water in the radiator, P, m 3 .

V s
Volume of the heat exchanger shell, C, m 3 .

V t
Volume of the heat exchanger tube, C, m 3 . V t1 Volume of tank 1, P, m 3 . V t2 Volume of tank 2, P, m 3 .

V w
Volume of water in the stack, P, m 3 .
NOTE: In this paper, the symbol h (W /K ) is named ''heat transfer parameter'' and denotes the product of the overall heat transfer coefficient U in W /(K ·m 2 ) and the heat transfer area A in m 2 , i.e. h = UA.

I. INTRODUCTION
Currently, there are serious environmental problems: global warming, air pollution and depletion of fossil fuel reserves. Governments are aware of it and they are taking measures to address them. New technologies in the energy sector can help to achieve these objectives. One of these technologies is µ-CHP systems [1]- [4]. A µ-CHP system produces electricity and heat (cogeneration) for the energy supply in the residential sector. The main advantage of these systems over traditional power generation systems is that they take advantage of the heat which is produced during the process of electrical energy generation, which results in a total energy efficiency that can exceed 90%, much higher than that of conventional systems, which is around 40%. Their high efficiency means lower fuel consumption, which leads to a reduction in the operating cost and a reduction in the amount of greenhouse gases emitted into the atmosphere [5]- [7].
There are several types of µ-CHP systems, depending on the technology on which they are based: internal combustion engine, Stirling engine, turbines or fuel cells [8], [9]. Among them, the one based on proton exchange membrane fuel cell (PEMFC) is a promising option, because this first mover has some advantages over the rest: 1) high electrical efficiency, 2) low heat-to-power ratio, 3) low level of noise and 4) low level of CO 2 emissions [5]- [7], [10].
There are already some commercial µ-CHP systems based on PEMFC and the trend is upwards, but mass production has not been reached yet, because their prices are high [6], [7]. In order to overcome this frontier, it is necessary to advance in several technical aspects, with the aim of improving the performance of these systems and reducing their costs [10], [11]. One of the most important work areas and where there is a greater opportunity for improvement is the temperature control of the PEMFC. In effect, the correct design of the stack cooling system, which includes its control, plays a crucial role in the durability, cost, reliability and energy efficiency of the stack [11]- [18]. However, to design a good temperature control of the stack, it is necessary to have a good model of the cooling system.
In the literature, there are several models of the cooling system of a PEMFC-based µ-CHP system [19]- [26], [31]. There are also thermal models of PEMFCs which are interesting in the context of this work, although they do not belong to µ-CHP applications [13], [27]- [30], [32], i.e. they do not contemplate a heat recovery subsystem, which is an essential unit of a µ-CHP system [7]. However, most of the models of the first group are static, not dynamic, and to develop a temperature control a dynamic model is needed. Moreover, few have been validated against test data from a real system. Therefore, there is a need for dynamic models that have been experimentally validated [5], [33]. For example, Ahn and Choe [13] developed a control-oriented model of a PEMFC stack cooling system. This model is dynamic and based on first principles, but it was not empirically validated. Cozzolino et al. [19] obtained a model of the cooling system of a PEMFC-based µ-CHP system and they validated it experimentally against a real system. However, their model is static, not dynamic. Barelli et al. [20] presented a dynamic model of a complete PEMFC-based µ-CHP system, including the cooling system. They used it to analyze the behavior of their µ-CHP system in response to variable electrical and thermal demands, but it is a linear model that represents the system only around an operating point and it was not validated experimentally. Asensio et al. [21] proposed a control-oriented model of the cooling system of a 600W PEMFC stack. The stack is cooled by water and is integrated into a laboratory equipment that emulates a real µ-CHP system. This model is able to provide the stack temperature as a function of the water flow rate, the stack inlet water temperature and the electric current demanded, but this model is static, not dynamic, and it was experimentally validated only for a value of the water flow rate (1 l/min), i.e. at an operating point. For an overview of the thermal models found in the literature and their main characteristics, see Appendix B.
In this paper, we present a model of the cooling system of a PEMFC-based µ-CHP system, which includes both the stack and the heat recovery unit. The purpose of this model is to be useful for the design of the temperature control of the PEMFC stack. The model is based on a real system which is located in our laboratory. It has been adjusted by means of an evolutionary optimization algorithm, in particular a genetic algorithm (GA). It is a first principles model, non-linear, dynamic and has been experimentally validated. A model of these characteristics (dynamic and experimentally validated) is helpful for the design of the stack temperature control and, according to the bibliographic evidence supplied in this introduction and in Appendix B, it represents a relevant contribution to the existing body of knowledge in the subject area of µ-CHP Modeling.
The rest of the article is structured as follows: In Section II, we present the complete µ-CHP system that we have in our laboratory. In Section III, we describe the part of the µ-CHP system that is going to be modeled (the cooling system), the modeling methodology, the tests performed for the adjustment and validation of the model, the genetic optimization FIGURE 1. The µ-CHP system. The PEMFC stack is fed with hydrogen and air, and generates electricity and heat. The electronic load emulates electrical energy demands and the fan-cooled radiator emulates thermal energy demands. The hot water tank serves as a buffer where heat can be stored. The whole system is controlled by a control unit, which is based on a CompactRIO from National Instruments. algorithm used for its adjustment, and the equations of the model. In Section IV, we present the results and discuss them. Finally, in Section V, we conclude.

II. DESCRIPTION OF THE MICRO-CHP SYSTEM
In this section, we describe the µ-CHP system on which the tests were carried out and whose cooling system we model. Fig. 1 shows the complete system. The core of the µ-CHP system is the PEMFC stack, which is from the company Nedstack, model 2.0HP. The stack consumes hydrogen and oxygen, and produces power and heat. It is able to produce up to 2 kW of electrical energy and 3.3 kW of thermal energy. The electrical energy is consumed by an electronic load and the thermal energy is removed from the stack by means of a water cooling system and then transferred, through the heat exchanger, to a hot water tank that serves as a buffer of heat. This is the part of the system that we model, the cooling system (i.e. the stack and the heat recovery unit, see Fig. 3), and will be described in detail in Section III-A.
The supply of oxygen comes from an external air compressor and the supply of hydrogen from a pressurized bottle. A mass flow meter and a proportional valve are installed in each supply line. By means of these devices the supply flow rates can be measured and varied depending on the electrical energy demanded by the electronic load. Both air and hydrogen must enter the stack under certain conditions of temperature and humidity. In order to condition them, the equipment includes a heating and humidification system. This system takes part of the humidity and heat present in the air at the stack air outlet and transfers them to the supply gases.
The electronic load demands electrical energy from the stack, by requiring an electrical current (i). In order to generate this current the stack must be fed with hydrogen and oxygen (which is taken from the air supply). The amount of these gases to be supplied to the stack depends on the electrical current demanded. The hydrogen supply flow rate is F H 2 , which enters the stack at a temperature of T H 2 in ; the air supply flow rate is F a , which enters the stack at a temperature of T a in . Both flows enter the stack at a relative humidity of 40% (calculated at stack temperature). Periodically, the circuit of hydrogen is purged, in order to remove impurities inside the stack and avoid its flooding. This is done by a valve which is activated by the signal U H 2 purge . The exhaust air from the stack leaves the stack at a temperature of T a out and at a relative humidity of 100%.
It is possible to program profiles of electrical demand in the electronic load and, in this way, to emulate the electrical consumption of a house (lighting and appliances). The equipment also has a fan-cooled radiator through which it is possible to emulate the thermal energy consumption (heating and hot water). This radiator, however, works only in two states, on and off. When on, it extracts heat from the secondary circuit of the cooling system.
Several magnitudes have to be controlled to ensure the correct operation of the system: inlet and outlet water temperature of the stack, air supply flow rate, hydrogen pressure at the stack anode inlet, and the relative humidity and the temperature of air and hydrogen supply flows. For this, the equipment includes several sensors and actuators, which are connected to a control unit based on a CompactRIO from National Instruments. The controllers of these magnitudes are implemented within this control unit. There are a total of six control loops.
The CompactRIO communicates with a PC via Ethernet. A Supervisory Control and Data Acquisition (SCADA) system runs on the PC, by which it is possible to monitor the process, to log its signals, to change the set point of any control loop and to parameterize the controllers. Fig. 2 shows its main screen. The SCADA and the control system that runs on the CompactRIO have been implemented in LabView.
The facility that has been described (µ-CHP system, SCADA and control unit) is an open system and, consequently, very flexible when it comes to making modifications in it or performing tests. The system has been conceived so that research and development activities aimed at optimizing the performance of µ-CHP systems can be conducted, with an emphasis on energy efficiency improvement, and this in two areas: control and energy management. The model that we present in this article is the first result of this line of work.

III. MODELING
In this section we describe, first, the part of the µ-CHP that we model, i.e. the stack cooling system (Section III-A), the modeling methodology (Section III-B), the tests conducted to obtain the model and to validate it (Section III-C), and the optimization algorithm used to adjust the parameters of the model (Section III-D). Then, we introduce the four submodels that together constitute the complete model, as well as their equations and the assumptions assumed in their conception (Sections III-E, III-F, III-G and III-H). Finally, we present the complete model (Section III-I).

A. THE COOLING SYSTEM
The cooling system of the µ-CHP plant which was presented in Section II is shown in Fig. 3 and is described in detail in this section. It consists of two water circuits, the primary circuit and the secondary circuit, coupled by a heat exchanger. The water of the primary circuit, which passes through the stack, is demineralized water, with a conductivity lower than 10 µS/cm.
The stack generates the electrical energy required by the electronic load and, as a by-product, produces heat within itself. The stack voltage (v) depends on the electrical current demanded (i) and on the temperature of the stack, among other variables, and is determined by the characteristic curve of the stack. The heat generated by the stack is removed by the water flow of the primary circuit (F w 1 ) and transferred to the secondary circuit through the heat exchanger. This heat finally ends up in the hot water tank (tank 2), which serves as a heat storage.
Pump 1 propels the water in the primary circuit and operates in steady state (U P1 = 100%). The water flow that passes through the stack (F w 1 ) is regulated by the motorized valve installed after pump 1 (U v ). There is an internal control loop of F w 1 , so a user (or a controller) can directly assign set points to F w 1 without having to worry about the motorized valve. This water flow rate is measured with a flow meter. By varying F w 1 it is possible to vary the amount of heat extracted from the stack. If F w 1 increases, more heat is extracted and the stack cools down. If it decreases, less heat is removed and the stack heats up. The outlet water temperature of the stack (T w out ) is used as a measure of the temperature of the stack.
The radiator emulates the demand of thermal energy. This radiator is activated by the signal R and has two states, on and off. When the radiator is on, it extracts heat from the secondary circuit and evacuates it to the environment, as a result the water temperature in tank 2 (T t2 ) decreases. When the radiator is off, the heat losses in the system are residual, although they exist. These residual losses are due to the transfer of heat from all the elements of the system to the environment. The ambient temperature is, therefore, a relevant variable in the process. A temperature sensor is installed for its measurement (T amb ).
The water flow rate of the secondary circuit (F w 2 ) is varied by means of the control signal of pump 2 (U P2 ), which is driven by a variable-speed drive. Like F w 1 , F w 2 is also internally controlled. The flow rate F w 2 is estimated from a measurement of the pressure at the inlet of the heat exchanger shell (P s in ). If this flow rate decreases, the amount of heat transferred through the heat exchanger also decreases, as a result less heat passes from the primary circuit to the secondary circuit, and this causes the water temperature at the stack inlet (T w in ) to increase; and vice versa, if the flow rate F w 2 increases, the amount of heat transferred from the primary to the secondary circuit increases and the temperature of the water entering the stack decreases.
For the stack to work properly, there must be a temperature control system that maintains T w out and T w in at its set points, minimizing the effect of the disturbances to which these signals will be subjected as a result of changes in the demands of electrical and thermal energy. The set points of T w out and T w in are, on the recommendation of the stack manufacturer, 65 • C and 60 • C, respectively. This control system will manipulate the flows F w 1 and F w 2 , which, as we have seen, have an influence on those temperatures and, consequently, serve as manipulable variables. This temperature control system, as noted in the Introduction, is of great importance, since the electrical efficiency and lifetime of the stack depend on it. Currently, a relatively simple temperature control system is implemented. This control consists of two proportional-integral (PI) controllers, one for each temperature, which were adjusted using a linear model of the process. The non-linear model that we present in this work would enable a control engineer to develop a more sophisticated temperature control and, consequently, to improve the efficiency and lifetime of the stack.

B. MODELING METHODOLOGY
The methodology used for the development of the model consists of the following stages: A) To delimit the part of the µ-CHP system that we want to model (the cooling system), specifying its inputs and outputs. B) Definition of the range of validity of the model. C) Conceptual division of the cooling system into subsystems, specifying the inputs and outputs of each subsystem. D) Modeling and validation of each subsystem separately. E) Interconnection of the submodels obtained in the previous stage to build the complete model. F) Fine adjustment of the complete model. G) Validation of the complete model. The modeling and validation of each subsystem (stage D) is divided, in turn, into the following steps: 1) Identification of the main elements of the subsystem.
2) Identification of the natural phenomena involved in the operation of these elements as well as the physical principles that govern them. 3) Introduction of simplifying hypotheses. 4) Determination of the structure of the submodel (set of equations) and its parameters, which is based on the physical principles and simplifying hypotheses. 5) Adjustment of the parameters of the submodel. 6) Evaluation of the results of the adjustment. If OK, validation (step 7); if not OK, revision of the simplifying hypotheses (step 3). 7) Validation of the submodel. If OK, the submodel is released; if not OK, revision of the simplifying hypotheses (step 3). The part of the µ-CHP plant that we want to model is the cooling system. This system was described in Section III-A and is represented in Fig. 3. Note that here we consider this system as a thermal process. This is important, because the selection of input and output signals that we detail below can only be understood in light of that fact. The signals that we consider as inputs to the system and as outputs are indicated in Fig. 3, the inputs in red and the outputs in blue. The electric current (i) and the stack voltage (v) provided by the stack are, considering the system as a thermal process, inputs, because the voltage and the current of the stack determine the amount of thermal energy it produces. The signals F H 2 and T H 2 in are also, from this point of view, inputs to the system (U H 2 purge only affects F H 2 ). However, we decided to ignore them, because the influence of the hydrogen flow rate and its temperature on the system turned out to be negligible, in comparison with the effect of the air flow rate (F a ) and its temperature (T a in ), which we do take into account. The signals F w 1 and F w 2 are inputs to the system. In fact, these are the variables which any temperature control will have to manipulate to control the temperature of the stack. The ambient temperature T amb is also considered as an input, because it affects the thermal behavior of the system, just like R, the activation signal of the radiator. The signals chosen as outputs are: T w out , T w in , T t2 , T a out , T s in and T s out .
In general, the range of validity of a model, i.e. the range of operation within which the model is intended to represent the process, depends on the purpose of the model. In our case, the purpose of the model is to serve for the design of the stack temperature control. By taking into account this basic consideration, we established that the model had to be valid for current demands (i) between 140 A and 200 A, and for water temperatures in tank 2 (T t2 ) between 53 • C and 56 • C. Outside these ranges, no controller will be able to keep T w out and T w in at their set points, because its control actions (F w 1 and F w 2 ) saturate if these ranges are exceeded. This is simply a physical restriction of the system. Since our model is oriented towards control design, further increasing its range of validity would increase the complexity of the model without adding any additional value to its purpose.
The adjustments of the parameters of the models, both those of the submodels (step 5) and the fine adjustment of the complete model (stage F), were made using the data set from the test for modeling (setM). The validations (stage G and step 7) were performed using the data set from the test for validation (setV). These two tests are described in detail in Section III-C.
Note that the process of modeling is iterative: if the results of the adjustment of a submodel (step 6) or its validation (step 7) are unsatisfactory, one goes back to step 3 to revise the simplifying hypotheses. This is because it is not always easy to determine a priori which physical phenomena are relevant and which are not. In the same way, it is not self-evident to what extent the equations representing these physical phenomena can be simplified without this entailing a deterioration in the capacity of the submodel for representing reality. What we were looking for was a model in first principles, suitable for the design of a temperature control, but not unnecessarily complicated. Therefore, in practice, for each submodel we started with a simple structure and evaluated the goodness of its adjustment. If it was insufficient, we revised the simplifying assumptions and incorporated more complexity, until the adjustment and the validation satisfied our expectations.
For example, when modeling the radiator, we initially assumed that the heat transfer parameters to the ambient for the states ON and OFF (h r ON and h r OFF ) were constant. Under this hypothesis, the adjustment achieved by the optimization algorithm displayed an intolerable discrepancy with respect to the test data. This fact led us to revise that hypothesis, assuming now that those heat transfer parameters depended on the water flow rate of the secondary circuit (F w 2 ), which is the water flow that passes through the radiator. The new adjustment improved substantially and, in addition, this improvement could be quantified.
It is interesting to note that in some cases the opposite happened, that is, the results of the adjustment indicated that a certain physical phenomenon that had been considered as relevant was not. For example, it was initially assumed that there was a loss of heat between the water in tank 2 and the environment, and that the corresponding heat transfer parameter could be approximated by a constant. When adjusting submodel SM4, to which tank 2 belongs, the optimization algorithm returned a value close to 0 W /K for that coefficient. This informed us that that phenomenon, which we had assumed as relevant, was in fact negligible.
Additionally, during the analysis of the adjustment of each submodel, we paid attention to the values of the parameters provided by the optimization algorithm and we checked whether these values were plausible or not. We could do this because the submodels are expressed in first principles, so their parameters have a direct physical sense. If the value of any of the parameters was unlikely, this was a clue that possibly there was some other relevant phenomenon that we were not considering, whose effect was being absorbed, incorrectly, by that parameter.
The fact that the model is formulated in first principles represented still another advantage during the modeling process: since the parameters of the model have a clear physical sense, it is simple and intuitive to estimate an initial search range for each of them. Once again, the analysis of the adjustment guided us in this task: when the search of a parameter tried to exit its search range, this suggested an increase in the limits of the search range of that parameter.
As can be seen, the analysis of the adjustment provided by the optimization algorithm gives many clues about which physical phenomena are relevant and are not being considered, and vice versa, which physical phenomena that were incorporated into the model are not relevant. For this reason, it serves as a starting point for the revision of the simplifying hypotheses. This is a useful help for the modeler, because it allows them to find, through this iterative process, the desired trade-off between model complexity and goodness of fit.
The adjustment of the parameters of each submodel (step 5) consists in finding the set of values of these parameters that provides the ''best'' fit between the submodel response and the experimental data from the test for modeling (setM), which means to solve an optimization problem. To quantify the goodness of the fit, i.e. what ''best'' means, we used an index. The formulation of the optimization problem and the index will be described in detail in Section III-D. As in any optimization problem, here there arises the difficulty that the cost function to be optimized (i.e. the index) presents local minimums and that the optimization algorithm, instead of finding the global optimum, is trapped in one of these local minimums. In our case, since we are dealing with the adjustment of a non-linear model with a lot of parameters, the probability that there are local minimums is high. To cope with this problem we did two things. First, we employed a genetic algorithm as optimization algorithm, because genetic algorithms have proven to be effective as global optimizers. Moreover, our research team has experience in their use. Second, we launched the algorithm three times (in parallel) for each adjustment. In this way, since the genetic algorithm begins with a random population, if the three executions of the algorithm returned the same solution, this was a reasonable proof that this solution was precisely the global optimum.
It should be added that the search space was discretized. That is, for each parameter, apart from its range of values, which defines the search space for that parameter, a minimum step was defined. For example, for a given heat transfer parameter the range could be from 0 to 10 W/K and the step 0.1 W/K; this would result in a set of 101 possible values for that parameter. The discretization of the search space makes optimization more efficient. The determination of the minimum steps is intuitive, thanks to the fact that the model is expressed in first principles. In addition, the minimum steps can be revised if the accuracy of the adjustment is unsatisfactory.
In order to evaluate the adjustment of each submodel (step 6), we used, on the one hand, the aforementioned cost function, which gives an objective measure of the goodness of the adjustment and, on the other hand, a visual comparison of the submodel response with the test data. This comparison allowed an assessment of how well the submodel was able to capture the dynamics of the plant, regardless of its accuracy, i.e. regardless of the value of the cost function.
The process of modeling, adjustment and validation of each submodel ends with a set of values, namely, the values of the parameters of the submodel that minimize the cost function. Once the four submodels were adjusted and validated independently, they were merged to form the complete model (see Fig. 11). A fine adjustment of the parameters of the complete model was then carried out (step F), by using again the optimization algorithm. On this occasion, the set of values of the parameters of the submodels already adjusted was included, as a starting solution, in the initial population of the algorithm. In addition, search ranges were narrowed and minimum steps were also reduced. This fine adjustment was launched three times too. In this way, we obtained the final set of values for all the parameters of the complete model. Finally, once the complete model was adjusted, it was validated against the data set from the validation test.

C. MODELING AND VALIDATION TESTS
The tests were conducted on the µ-CHP system by using the SCADA implemented in LabView (Fig. 2). Two tests were carried out: 1) one test for the modeling and for the adjustment of the parameters and 2) one test for the validation of the model. They were conducted on different days, the test for modeling on 10/13/2017 and the test for validation on 10/17/2017, and both have an approximate duration of 2.5 hours. In both tests, the process was excited by VOLUME 7, 2019 manipulating the variables i, F w 1 , F w 2 and R, that is, the electrical energy demand, the water flow rates of the primary and secondary circuits, and the thermal energy demand. These signals are represented in Fig. 4 (test for modeling) and in Fig. 5 (test for validation). In addition, the following signals were measured and recorded: on the one hand, T w out , T w in , T t2 , T a out , T s in and T s out (which are the output signals of the process and will be presented in Section IV, where the adjustment and validation results are displayed) and, on the other hand, T amb , v, T a in and F a . These last signals, remember, are also inputs to the process, but not manipulable. For reasons of simplicity and clarity, these last signals have not been included in the figures.
The tests are step tests. This is relevant for two reasons. First, because in a real µ-CHP system the electrical and thermal energy demands (in our system, i and R) have usually the form of steps. Second, because a temperature controller may apply also high frequency control actions (in our system, F w 1 and F w 2 ) in response to disturbances. Therefore, the fast dynamics of the process must be captured by the model. Both tests are analogous in nature (step tests), although quantitatively different (they differ in the values of the steps and in their sequence). The tests were carried out in open loop. This was done to prevent the PI controllers from attenuating the mutual effects produced by the programmed steps in the manipulated variables, which would have partially concealed the dynamics of the system. For example, a change in F w 1 affects T w out , but also T w in , so the controller of T w in will act on F w 2 to counteract the disturbance produced by F w 1 ; as a result the effect of F w 1 on T w in will be attenuated and it will be difficult to identify. Again, it is important that the model (since it is a model oriented to control) captures the dynamics of the process well, because if the control designed is based on a model that does not satisfy that requirement, such a control could make the process unstable, especially if it is an aggressive control.

D. THE OPTIMIZATION ALGORITHM
The parameters of each of the four submodels are adjusted using the data set from the modeling test (step 5). This adjustment task is formulated as an optimization problem and this problem is solved by an optimization algorithm.
Given the mathematical structure of a submodel, expressed in its general form as a set of non-linear differential equations: where u(t) are the inputs, y(t) the outputs, x(t) the state variables and θ = [θ 1 , ..., θ n ] the vector of parameters to be adjusted, the optimization problem consists in finding the set of values of θ that minimize a certain cost function (J ), which measures the goodness of fit of the model against the real data. This is to say: subject to the following restrictions: Basically, these restrictions define the search space, in other words, θ i and θ i are the lower and upper limits of the search range of the parameter θ i . Solving that optimization problem is precisely what the optimization algorithm does. There are several ways in which the cost function can be defined. We defined it as: where: where n is the number of outputs,ŷ k (t) is the output k of the process, y k (t) the output k of the submodel and T is the duration of the test. That is to say, what is minimized (J ) is the sum of the average errors in the outputs of each submodel. Note that, due to the way in which the index has been defined, all the errors in the submodel outputs have the same relative importance.
In short, the optimization algorithm takes as inputs: 1) the mathematical structure of the submodel, 2) the modeling test (setM), 3) the cost function (J ) and 4) the search space; and it returns the values of the parameters of the submodel that minimize the cost function. The optimization algorithm used is a genetic algorithm. In particular, we use the genetic algorithm described in [40], available in Matlab Central, at https://es.mathworks.com/matlabcentral/fileexchange/ 39021-basic-genetic-algorithm?s_tid=prof_contriblnk.

E. SUBMODEL 1 (SM1)
The subsystem (SS1) corresponding to this submodel (SM1) consists of a single element, the PEMFC stack (see Fig. 3). The inputs of this submodel are F a , T amb , T a in , v, i, T w in and F w 1 ; and its outputs are: T w out and T a out (see Fig. 11). The parameters to be adjusted in this submodel are: V w , V a , k a , h fc 2max , h fc 2min , h fc loss , h a max , h a min , h w max , h w min , h aw , cal T w out and cal T a out . These signals and parameters will be explained below, as the equations of this submodel are introduced. But first, a clarifying note: In this section and in the three that follow, we describe the equations of each element of the submodels. The binding equations, i.e. the equations that model the connections between elements within a submodel or between submodels within the complete model, VOLUME 7, 2019 for example T w in = T p1 out , we do not specify them, since they are trivial just by looking at Fig. 11, where they are represented by arrows. Moreover, since a lot of the equations are based on the same modeling principles, we have included the Appendix A', where the five most common modeling principles are numbered and summarized in a table (Table 5). In order to simplify the introduction of the equations and to avoid unnecessary repetitions, we will refer to these modeling principles (MP) throughout the following four sections.
The PEMFC stack working principle diagram that we have assumed is displayed in Fig. 6. It shows the physical phenomena that have been considered as relevant, as well as the main variables, constants and parameters involved. Two electrochemical reactions take place at each cell of the stack. At the anode: At the cathode: The maximum electrical energy that a cell can produce is given by the Gibbs free energy: And the ideal potential of the cell is: where n is the number of electrons released in the reaction (two) and F is Faraday's constant (96485 As/mol). At 25 • C and atmospheric pressure E c = 1.23 V (ideal standard potential). Therefore, the maximum electrical energy that the stack can generate, assuming that there are no losses of any kind (ideal situation), is: where n c is the number of cells of the stack (16 cells) and i is the electric current generated. However, the electric power actually available between the terminals of the stack is lower and it is given by: where v is the voltage between the terminals of the stack. The difference between the ideal electrical energy and the electrical energy actually available corresponds to the energy that the stack produces in the form of heat, which is Both signals v and i are inputs to this submodel and they were measured during the tests. Part of this heat is extracted by the water flow of the primary circuit, another part passes to the air that runs through the stack (Q wa ) and the rest is transferred to the stack shell (Q wfc ). Thus, the evolution of the water temperature at the stack outlet can be expressed as: The constant k w is the product c w ρ w ; c w is the specific heat of the water and ρ w the density of the water, both of them for a water temperature of 60 • C, which we assume as a representative value of the water temperature in both the primary and the secondary circuits. The values of these two parameters were taken from [39] and their product gives k w = 4112120.21 J /m 3 K . The volume V w is the volume of water inside the stack and it is a parameter to be adjusted by the optimization algorithm.
The heat difference Q w is calculated according to MP2 (F = F w 1 , T in = T w in , T out = T w out , k = k w ) and corresponds to the amount of heat that is removed from the stack by the water flow of the primary circuit.
The amount of heat transferred from the water to the stack shell is: which depends, on the one hand, on the temperature gradient between the temperature of the water inside the stack (T w ) and the temperature of the stack shell (T fc ), and on the other hand, on the heat transfer parameter h fc 2 . As a representative value of the temperature of the water inside the stack, we chose the average temperature between the temperature of the water at the stack inlet and the temperature of the water at the stack outlet (MP4: T in = T w in , T out = T w out ). The heat transfer parameter h fc2 is assumed to depend on the water flow, according to MP5 (h min = h fc 2min , h max = h fc 2max , F min = F w 1min , F max = F w 1max , F = F w 1 ). That is to say, it has been assumed (as an approximation) that there is a linear relationship between the water flow rate F w 1 and the heat transfer parameter h fc 2 , i.e. the higher the water flow rate, the more heat will be transferred from the water to the stack shell, for a given temperature gradient. Note that this relation between h fc2 and F w 1 is the result of a trial and error process, as we explained in Section III-B. F w 1min and F w 1max are the minimum and the maximum values of F w 1 , 2.9005·10 −5 m 3 /s (1.74 liters/minute) and 1.1657·10 −4 m 3 /s (6.99 liters/minute), respectively. The constants h fc 2min and h fc 2max are the minimum and the maximum values of h fc 2 , and they are parameters to be adjusted by the optimization algorithm.
The heat flow rate transferred from water to air is given by: which depends, on the one hand, on the temperature gradient between the temperature of the water inside the stack (T w ) and the temperature of the air inside the stack (T a ), and on the other hand, on the heat transfer parameter h fc 1 . Here also, a representative value of the temperature of the air inside the stack (T a ) is calculated after MP4 (T in = T a in , T out = T a out ). The overall heat transfer parameter between the water and the air (h fc 1 ) is calculated from three independent heat transfer parameters (h w , h a and h aw ), according to the following equation: It is also assumed, as an approximation, that the heat transfer parameter h a depends linearly on the air flow rate (F a ), according to MP5 (h min = h a min , h max = h a max , F min = F a min , F max = F a max , F = F a ). F a min and F a max are the minimum and the maximum values of F a , 0.0020 m 3 /s (120 liters/min) and 0.0027 m 3 /s (162 liters/minute), respectively. The constants h a min and h a max are the minimum and maximum values of h a , and they are parameters to be adjusted by the optimization algorithm. Similarly, we assume that the heat transfer parameter h w is linearly dependent on the water flow rate (MP5: h min = h w min , h max = h w max , F min = F w 1min , F max = F w 1max , F = F w 1 ). The constants h w min and h w max are the minimum and maximum values of h w . Both are parameters to be adjusted by the optimization algorithm. Finally, the heat transfer parameter h aw is assumed to be constant, and it is also a parameter to be adjusted.
The heat flow rate transferred from the water to the air (Q wa ) will increase the temperature of the air volume within the stack. Thus, the air temperature at the stack outlet (T a out ) can be calculated by the following expression: The constant k a is the product c a ρ a ; c a is the specific heat of the air and ρ a the density of the air inside the stack. Since the air density within the stack depends on the pressure, it is not easy to estimate k a theoretically, as we did for k w . For this reason, it will be adjusted by the optimization algorithm. The volume V a is the volume of air inside the stack and is also a parameter to be adjusted. The heat difference Q a is calculated according to MP2 The heat transferred from the water to the stack shell (Q wfc ) will cause the temperature of the stack shell (T fc ) to increase, according to the following expression: C fc is the stack heat capacity. This constant is calculated theoretically as the product of the mass of the stack, 14.3 kg, and the specific heat of aluminum, 905 J /kg · K (value taken from [39]). That product gives C fc = 12941.5 J /K . The heat flow rate Q fc loss is the amount heat that is transferred from the stack shell to the environment, that is, the heat lost to the environment due to natural convection. It is given by: That is to say, it depends, on the one hand, on the temperature gradient between the temperature of the stack shell (T fc ) and the ambient temperature (T amb ), and on the other hand, on the heat transfer parameter h fc loss , which is assumed to be constant and is a parameter to be adjusted. During the adjustment of this submodel, when comparing its output T w out with the output of the real system, we observed that their dynamics were very similar, but there was a constant difference between the two signals. From this fact, we inferred that this discrepancy could be due to a calibration error in the temperature sensor that measures the water temperature at the stack outlet. Although that sensor was calibrated before it was assembled, the assembly itself could have caused a slight decalibration. We also thought that the same problem could be taking place in the temperature sensor that measures the air temperature at the stack outlet, although with a less noticeable effect. So, in order to model these hypothetical calibration errors, we included two additional parameters in the submodel, i.e. two correction coefficients of the calibration errors. These coefficients are cal T w out and cal T a out , they are constant and are added to the outputs of the submodel. These two parameters were adjusted by the optimization algorithm together with the rest of parameters. Their search space was set between −2 • C and 2 • C. Note that if the hypothesis were incorrect, the optimization algorithm would have returned null values for these coefficients. In a subsequent test, we confirmed that there was indeed a calibration error in the measurement of the water temperature at the stack outlet and the value of this error corresponded approximately to the value provided by the optimization algorithm for the parameter cal T w out .
Note that our goal is to develop a thermal model of the cooling system. Consequently, the modeling of the electrochemical phenomena involved in the functioning of the PEMFC stack does not belong to the scope of this work. Moreover, this kind of modeling have already been thoroughly studied in the existing literature. In particular, note that the submodel of the PEMFC stack (SM1) does not include the typical polarization curve (which provides the stack voltage (v) as a function of the electric current (i)), because this polarization curve would be part of an electrical model, not a thermal one. For this reason, the stack voltage (v) is considered as an input to the model and it is measured. This measured signal (v) is used in (12) and (13) to calculate P e and Q heat , respectively.

F. SUBMODEL 2 (SM2)
The subsystem (SS2) that corresponds to this submodel (SM2), consists of a single element, pipe 4 (see Fig. 3). The inputs of this submodel are T p4 in , F w 2 and T amb ; and its output is: T p4 out (see Fig. 11). This submodel is quite simple and has only three parameters to be adjusted: V p4 , h p4 loss and cal T p4 out . A diagram with the thermal phenomena assumed in pipe 4 is shown in Fig. 7. The heat difference in pipe 4 ( Q p4 ) is calculated according to MP2 (F = F w 2 , T in = T p4 in , T out = T p4 out , k = k w ), where T p4 in is the water temperature at pipe 4 inlet, T p4 out is the water temperature at pipe 4 outlet, the constant k w was explained in the previous section and F w 2 is the water flow rate of the secondary circuit. It is assumed that there is a loss of heat from the water to the environment, which is given by: h p4 loss is the heat transfer parameter from the water within the pipe to the environment, which is assumed to be constant and is a parameter to be adjusted. The temperatureT p4 is the average value between the inlet water temperature and the outlet water temperature (MP4: T in = T p4 in , T out = T p4 out ), value that is assumed to be representative of the water temperature in the pipe. By considering the incoming and outgoing heat flow rates in pipe 4, the variation of the water temperature at its outlet is given by: where V p4 is the volume of the pipe and a parameter to be adjusted. Also in this submodel, as we did in submodel 1, in order to correct a possible calibration error in the measurement of the outlet water temperature, a correction coefficient was included (cal T p4 out ). This coefficient is constant and a parameter to be adjusted. Its search range is [−2 • C, 2 • C].

G. SUBMODEL 3 (SM3)
The subsystem (SS3) that corresponds to this submodel (SM3) consists of three elements: heat exchanger, tank 1 and pipe 1 (see Fig. 3). The inputs of this submodel are T t in , T s in , F w 1 , F w 2 and T amb ; and its outputs are: T s out and T p1 out (see Fig. 11). The parameters to be adjusted are: h s min , h s max , V p1 and h p1 loss . A working principle diagram of the heat exchanger is shown in Fig. 8 and of tank 1 in Fig. 9. The diagram of pipe 1 is identical to pipe 4's, which was introduced in the previous section. Let us begin with the heat exchanger. The transfer of heat in the heat exchanger takes place from the tube side to the shell side. The water flow of the primary circuit (F w 1 ) circulates through the tube and the water flow of the secondary circuit (F w 2 ) circulates through the shell. Thus, the heat is transferred from the primary circuit to the secondary circuit. The heat exchanger is configured in a parallel-flow arrangement.
The primary water flow enters the tube at a temperature T t in and leaves it at a temperature T t out . Similarly, the secondary water flow enters the shell at a temperature T s in and leaves it at a temperature T s out . The heat difference in the tube side ( Q t ) and the heat difference in the shell side ( Q s ) are calculated after MP2 for F = F w 1 , T in = T t in , T out = T t out , k = k w , and The amount of heat transferred from the tube side to the shell side is given by: where h ts is the overall heat transfer parameter of the heat exchanger and T ts the temperature gradient between tube and shell. The coefficient h ts depends on three heat transfer phenomena: 1) forced convection in the tube side, 2) conduction through the wall between tube and shell and 3) forced convection in the shell side. It can be calculated as: where n ch is the number of channels in the heat exchanger (9 channels), h t is the heat transfer parameter due to forced convection in the tube side, h e is the conduction heat transfer parameter, which is constant, and h s is the heat transfer parameter due to forced convection in the shell side. h t is assumed to depend linearly on F w 1 according to where h t min and h t max are the minimum and maximum values of h t . Both are parameters to be adjusted by the optimization algorithm. Similarly, the coefficient h s is assumed to depend linearly on F w 2 according to MP5 (h min = h s min , h max = h s max , F min = F w 2min , F max = F w 2max , F = F w 2 ), where F w 2min and F w 2max are the minimum and maximum values of F w 2 , 4.2137·10 −5 m 3 /s (2.53 liters/minute) and 1.5493·10 −4 m 3 /s (9.29 liters/minute), respectively. h s min and h s max are the minimum and maximum values of h s . Both are parameters to be adjusted. h e was calculated theoretically from the technical data of the heat exchanger (area, thickness and thermal conductivity of the plates), its value is h e = 216.72 W /K .
On the other hand, the temperature gradient T ts in (23) is the so-called logarithmic mean temperature difference (LMTD) [39], which is calculated according to: where, for a parallel-flow configuration: and By considering the incoming and outgoing heat flow rates in tube and shell sides, it follows that their outlet water temperatures have to satisfy: and where V t is the volume of water in the tube side and V s is the volume of water in the shell side. These volumes were calculated theoretically from the dimensions of the heat exchanger, both resulting in 1.8 · 10 −4 m 3 (0.18 liters). Now let us look at the equations of tank 1. The heat difference in tank 1 ( Q t1 ) is calculated after MP2 (F = F w 1 , T in = T t1 in , T out = T t1 out , k = k w ), where T t1 in is the temperature of the water entering the tank and T t1 is the temperature of the water leaving the tank. It was not necessary to assume any lost of heat to the environment, so the water temperature at tank 1 outlet is given by: where V t1 is the volume of water in the tank, which is constant and is also a parameter to be adjusted. Finally, the equations that model pipe 1 are exactly the same as those of pipe 4 (see Section III-F). The heat difference in pipe 1 ( Q p1 ) is based on MP2 (F = F w 1 , T in = T p1 in , T out = T p1 out , k = k w ), where T p1 in is the water temperature at pipe 1 inlet and T p1 out is the water temperature at pipe 1 outlet. The loss of heat to the ambient is: where h p1 loss is the heat transfer parameter to the environment andT p1 is calculated after MP4 (T in = T p1 in , T out = T p1 out ). So, the balance of heat flows in the pipe 1 results: where V p1 is the volume of pipe 1, which is unknown and, therefore, a parameter to be adjusted.

H. SUBMODEL SM4
The subsystem (SS4) corresponding to this submodel (SM4) consists of four elements: pipe 2, radiator, pipe 3 and tank 2 (see Fig. 3). The inputs of submodel are T p2 in , F w 2 and R; and its output is: T t2 (see Fig. 11). The parameters to be adjusted are: h r OFFmin , h r OFFmax , h r ONmin , h r ONmax , V t2 , V r and T amb r . The diagrams of pipes 2 and 3 are similar to that of pipe 4, which has already been introduced in Section III-F, with the only difference that for pipes 2 and 3 it was not necessary to model the losses of heat to the environment, because they are negligible compared to the amount of heat that the radiator extracts from the secondary circuit, in both states, on and off. The diagram of tank 2 is identical to that of tank 1, presented in Section III-G. The diagram of the radiator is shown in Fig. 10. Let us begin with pipes 2 and 3. The heat difference in pipe 2 ( Q p2 ) is calculated after MP2 (F = F w 2 , T in = T p2 in , T out = T p2 out , k = k w ), where T p2 in is the water temperature at the pipe inlet and T p2 out is the water temperature at the pipe outlet. Since there are no heat losses to the environment, the temperature of the water leaving the pipe is given by: where V p2 is the volume of pipe 2. This constant is calculated theoretically from the length and the cross-section area of the pipe, its value is 7.0686 · 10 −4 m 3 (0.7069 liters). The heat difference in pipe 3 ( Q p3 ) is calculated after MP2 (F = F w 2 , T in = T p3 in , T out = T p3 out , k = k w ), where T p3 in is the water temperature at the pipe inlet and T p3 out is the water temperature at the outlet of the pipe. Thus, the water temperature at the outlet of the pipe is given by: where V p3 is the volume of pipe 3. Also theoretically calculated, its value is 3.5343 · 10 −4 m 3 (0.3534 liters). Regarding tank 2, the heat difference Q t2 is calculated where T t2 in is the water temperature at the tank inlet and T t2 is the water temperature of the tank. As heat losses to the environment are not considered, T t2 is given by: where V t2 is the volume of water in tank 2 and a parameter to be adjusted. Finally, we address the radiator, which is the most important element of this submodel. The water flow of the secondary circuit (F w 2 ) circulates through the radiator. It enters the radiator at a temperature T r in and leaves it at T r out . The heat difference in the radiator ( Q r ) is calculated after MP2 (F = F w 2 , T in = T r in , T out = T r out , k = k w ). The amount of heat transferred from the hot water within the radiator to the environment depends on whether the fan is on (R = ON ) or off (R = OFF), since in the first case the transfer of heat occurs by natural convection, whereas in the second by forced convection. This results in two different heat transfer parameters, h r ON and h r OFF , one for each state (ON, OFF). Therefore, the heat transfer parameter of the radiator can be expressed as follows: Each of these coefficients depends, in turn, on the water flow rate F w 2 . As we did before, we also assume here that these dependencies are linear. h r ON is calculated according to MP5 , where h r OFFmin and h r OFFmax are the minimum and the maximum values of h r OFF and also parameters to be adjusted. The amount of heat transferred to the environment depends, as we said, on the heat transfer parameter h r , but it also depends on the temperature gradient between the temperature of the water within the radiator and the temperature of the air that surrounds it: T r is the average value of the water temperatures at the inlet and outlet of the radiator, and it is the value that has been assumed as representative of the water temperature inside the radiator, calculated after MP4 (T in = T r in , T out = T r out ). T amb r is the temperature of the air surrounding the radiator. This temperature is different from the temperature assumed as ambient temperature for the rest of submodels (T amb ), for whose measurement a sensor was installed. The reason for this difference is that the radiator is mounted next to a window that is partially open (it can not be closed because through it the purge of hydrogen is evacuated). Thus, the temperature of the air surrounding the radiator is not the one inside the laboratory (T amb ), but, presumably, a weighted average of two variables, the temperature inside the laboratory and the temperature outside the building (considerably lower than the former, since the tests were conducted in October). As we had not installed any sensor to measure the air temperature surrounding the radiator (T amb r ), we decided to estimate it. We assumed for this temperature a constant and unknown value, a value within the temperature range between the inside temperature (T amb ) and the outside temperature, and then we adjusted its value as one more parameter of the model, with the data set from the modeling test (that is why T amb r is listed as a parameter of submodel SM4 at the beginning of this section). Moreover, since the validation test was performed only four days after the modeling test and during the same time span, we initially considered that it would be plausible to assume that the ambient temperature of the radiator (T amb r ) on that day was approximately equal to the one on the day of the modeling test. We will return to this hypothesis in Section IV-B, when we comment the results of the model validation.
Lastly, the balance of heat flows in the radiator determines the water temperature variation at the radiator outlet:   the six outputs of the model are compared with the six outputs of the real system, having been both, model and system, excited with the same input signals (Fig. 4). The values of the parameters that achieve that adjustment are shown in Table 1.
As can be seen from these results, the adjustment is satisfactory. First, the accuracy of the model is good: the overall average error (the average of the average errors of the six outputs) is 0.17 • C (see Table 2). Second, the adjustment is balanced, that is, all the outputs of the model match their corresponding process outputs approximately equally well (T s in is the best adjusted and T t2 the worst, and there is a difference of 0.089 • C between their average errors). Third, the model captures the dynamics of the system well. Fourth, the values of the parameters (Table 1) are realistic, which is a good evidence that the model remains close to reality, from a conceptual point of view.

B. RESULTS OF THE MODEL VALIDATION
The results of the validation of the complete model are shown in Fig. 13. In this figure, the six outputs of the model are compared, again, with the six outputs of the real system, having been both, model and system, excited with the same input signals, in this case the ones corresponding to the validation test (see Fig. 5).
The first thing that can be observed is that the response of the model against the validation test data, in comparison with the results of the model adjustment (see Fig. 12), deteriorates. The accuracy (overall average error) achieved in this case is 0.438 • C (see Table 3), compared to 0.17 • C that was achieved in the adjustment with the modeling test data. A certain degree of deterioration in accuracy when validating, with respect to the adjustment in the modeling phase, is usual and normal. However, in our case, we think that this deterioration has an additional cause, namely, that the hypothesis assumed, according to which the ambient temperature for the radiator (T amb r ) on the modeling test day was equal to the one on the validation test day, is false. We believe this for the following reason. If T amb r is readjusted (maintaing the same values for the remaining 29 parameters) with the validation test data, the optimization algorithm returns a new value of T amb r , 27.86 • C (0.75 • C higher than the one which resulted from the adjustment), and when this new value of T amb r is used, the fit between the model outputs and the real signals of the validation test substantially improves (see Fig. 14): the overall average error shifts from 0.438 • C to 0.296 • C (see TABLE 3. Average error of each output of the model with respect to the outputs of the process when it is tested against the validation test data. Table 4), which means an improvement of 0.142 • C. Note also that the temperature T amb r represents one of the two reference points of the thermal process (the other one is T amb ), i.e. it is something like the ground in an electrical circuit (using the analogy between thermal and electrical systems). Consequently, it only determines the operating point of the system in steady state. This is why the responses of the model before and after the adjustment of T amb r are practically identical in shape, differing only in their average value (compare Fig. 13  and 14). Certainly, the ideal would have been to dispose, from the beginning, an additional temperature sensor next to the radiator, in order to have a real measurement of T amb r on both days.
A second observation from Fig. 13 is that the model, despite the aforementioned deterioration in its accuracy, continues to faithfully capture the dynamics of the process, exactly the same as it did for the modeling test data (Fig. 12). This means that, in terms of the ability of the model to represent the dynamics of the process, there has been no deterioration.
The ability of a model to faithfully represent the dynamics of the process is precisely the most important characteristic that a control-oriented model must fulfill, because any temperature control of a PEMFC stack will always incorporate an integral action to compensate for errors in steady state and, therefore, it will absorb any constant difference that may exist between the model and the process. Thus, after weighing up both aspects, accuracy and ability to faithfully represent the dynamics of the process, and given that the former is not critical and the latter is crucial, we decided to release the model, because it fulfills its purpose, which is to be useful for a subsequent design of a temperature control.

V. CONCLUSION
In this article, we have presented a model of the cooling process of a µ-CHP system based on a PEMFC stack. The model has been developed to be used for the design of the temperature control of the stack. It is based on first principles, dynamic, non-linear and has been validated empirically. The results show that the model is able to faithfully represent the process. Especially, it is remarkable the accuracy with which the model captures the dynamics of the real process, which is the most important feature that a control-oriented model TABLE 4. Average error of each output of the model with respect to the outputs of the process when it is tested against the validation test data, after the readjustment of T amb r . must fulfill. Therefore, the model is useful for a subsequent temperature control design.
As far as we know, there is no model in the literature that satisfies these characteristics (see Appendix B) and its need has been pointed out in several reviews (see [5], [33]). Therefore, the main contribution of this work consists of developing and publishing a model of the cooling process of a µ-CHP system based on a PEMFC stack (i.e. stack and heat recovery unit) which is dynamic and has been experimentally validated. The advantage of such a model is that it enables the design of a more sophisticated temperature control of the stack and, consequently, with a better performance than the one that a control based on a simple model (such as, for example, a First Order Plus Dead Time (FOPDT) model) can achieve.
The model has been implemented in Matlab/Simulink and it is accessible, together with the test data that were used for its development, at http://hdl.handle.net/10251/118336. Moreover, the methodology employed to develop the model has been described in detail. This, added to the fact that the model is expressed in first principles (so it has a close relation with the structure of the real system and its parameters have physical meaning), allows other modelers to take our model as a starting point and, by making the necessary modifications, to build models of their µ-CHP systems.
In a future work we plan to design, by using the model we have presented here, a temperature control for the stack of our µ-CHP system and to validate it experimentally, in order to improve the performance of the system with respect to the performance that the current control (based on a linear and simpler model) offers. Table 5 summarizes the main common modeling principles on which the four submodels (Sections III-E, III-F, III-G and III-H) are based. Table 6 summarizes the main characteristics of the thermal models found in the literature review that we carried out VOLUME 7, 2019   Note that this table refers only to thermal models. In some of the works listed, the authors develop other models as well (for example, an electrochemical model of the PEMFC), apart from a thermal model. In these cases, the attributes ''dynamic or static'' and ''experimentally validated'' refer to the thermal models exclusively.

APPENDIX B THERMAL MODELS IN THE LITERATURE
As Table 6 shows, there are only three models that satisfy being dynamic and having been experimentally validated (Musio et al. [27], Real et al. [28] and San Martín et al. [29]). However, in these works 1) the stack is refrigerated by air (by means of a fan) and 2) the heat removed from the stack is evacuated to the atmosphere, i.e. the systems that these authors model do not include a heat recovery subsystem, which is the feature that makes an energy generator system a CHP system (cogeneration), (see [7]). In other words, the systems modelled in these works are substantially different from the system that we model in our work. Consequently, a direct comparation is not possible. Nonetheless, it is worth saying some words about these models.
In Musio et al. [27], the validation results show a discrepancy between real and simulated stack temperatures of almost 4 • C (maximum error) and the model dynamics do not represent well the dynamics of the real system. Moreover, the stack air outlet temperature is not validated. In Real et al. [28], the validation results are remarkable. The stack temperature provided by this model follows the real stack temperature quite well within a wide range of temperatures. However, the model is validated only for changes in the electric current demanded to the stack and not for changes in the manipulable variable (fan voltage), which stays constant during the validation test. In order to design a dynamic control of the stack temperature, the model should represent not only the dynamic relation between the output (stack temperature) and the perturbation (electric current) but also the dynamic relation between the output and the manipulable variable, because this manipulable variable will be, once the control is designed, the control action. The model presented in Martín et al. [29] has exactly the same limitation. This model achieves a satisfactory matching between model and process for changes in the electric current demanded to the stack (perturbation), but there is no validation for changes in the manipulable variable (fan voltage). The reason why these two last works do not validate their models for changes in the manipulable variable (fan voltage) may be the following: in both cases the system modeled is a closed system (Nexa 1200) and, therefore, it probably hinders from accessing to the fan voltage directly, so that it can be freely manipulated, and independently from the electric current demanded.