Particle Swarm Optimization-Based Power and Temperature Control Scheme for Grid-Connected DFIG-Based Dish-Stirling Solar-Thermal System

: Variable-speed operation of a dish-Stirling (DS) concentrated solar-thermal power generating system can achieve higher energy conversion efﬁciency compared to the conventional ﬁxed-speed operation system. However, tuning of the controllers for the existing control schemes is cumbersome due to the presence of a large number of control parameters. This paper proposes a new control system design approach for the doubly-fed induction generator (DFIG)-based DS system to achieve maximum power point tracking and constant receiver temperature regulation. Based on a developed thermo-electro-pneumatic model, a coordinated torque and mean pressure control scheme is proposed. Through steady-state analysis, the optimal torque is calculated using the measured insolation and it serves as the tracking reference for direct torque control of the DFIG. To minimize the tracking error due to temperature variation and the compressor loss of the hydrogen supply system, four optimal control parameters are determined using particle swarm optimization (PSO). Model-order reduction and the process of the pre-examination of system stability are incorporated into the PSO algorithm, and it effectively reduces the search effort for the best solution to achieve maximum power point tracking and maintain the temperature around the set-point. The results from computational simulations are presented to show the efﬁcacy of the proposed scheme in supplying the grid system with smoothened maximum power generation as the solar irradiance varies. search efﬁcacy of the proposed in the grid with performance as of detailed of as due to can play impact on efﬁciency,


Introduction
Modern electricity supply systems have experienced significant changes due to the rapid increase of renewable energy generation in the past decades. Amongst various renewable energy technologies, dish-Stirling (DS) concentrated solar-thermal power generating system uses a parabolic dish-like reflector to collect and concentrate sunlight to a small heat receiver located at the focal point of the reflector. The high temperature achieved at the focal point, in which the receiver is located, is used as a heat source for a closed-cycle external heat engine called Stirling engine [1]. The DS system is capable of operating at a high energy conversion efficiency and releases no emissions, and thus it is considered to be a promising sustainable technology for future power grid structures such as microgrid and distributed generation [2].
For these reasons, in recent years, grid-connected DS system has received the attention of researchers aiming toward improving its modelling, design, control, and simulation, in order to increase the energy conversion efficiency and lower the overall cost of the system. A mathematical model intended for the optimal opto-geometric design of the concentrator and heat receiver is proposed in [3] and it is experimentally verified in [4]. Thermal efficiency of the heat receiver is studied and validated in [5]. A thermal model of a cavity receiver is developed for optimization in [6] where radiative exchange is evaluated using the radiosity method. The authors in [7] show that the global energy efficiency ranges between 19% and 26%, which can be achieved through optimal design based on a modified Iwamoto model of a Stirling engine and a solar prediction model. Further increase of utilization and efficiency can be achieved by incorporating novel components such as a free-piston Stirling engine [8], high-temperature superconducting linear synchronous generator [9], or integrated thermal energy storage [10]. These works mainly focus on the optimization and design for a specific component. System-level modelling, control, and simulation have been studied in References [11,12], where an ideal adiabatic model of a double-acting kinematic Stirling engine was used. In order to achieve acceptable accuracy, this Stirling engine model requires a very small step size (about 2 µs) for time-domain simulation as it has to be solved numerically as a cyclic boundary problem [13]. Due to the high complexity of the model, the control system can only be designed using empirical methods. A simple controller with proportional gain has to be used for the receiver temperature control, which causes lower thermal efficiency during the low insolation period. In [14], an averaged-value model of a solar Stirling engine is developed for the purpose of grid-integration studies and model-based control system design. The required step size for simulation is greatly increased to be compatible with that of electromechanical transient simulation. The control system has been improved with reduced impact on the interconnected grid system during insolation variation and/or grid fault. However, the scale of the conventional grid-connected DS system using squirrel-cage induction generators (SCIGs) is restricted due to the lack of control capability to meet the grid code requirements such as to provide frequency and voltage support. Although frequency support, such as a fixed-speed generation system, could be achieved by changing the working temperature [15], the significant reactive power injected to the grid system can still cause a severe voltage stability problem especially in a weak grid, as has been investigated in works [16,17].
In contrast, a variable-speed DS system is equipped with a doubly-fed induction generator (DFIG) [18,19] or permanent magnet synchronous generator (PMSG) [20,21]. With the back-to-back ac/dc/ac converters and advanced control techniques, independent regulation of the active power and reactive power becomes feasible. Maximum power point tracking (MPPT) and ancillary services, such as the frequency and voltage regulation, can be achieved for additional technological and economic benefits. Compared to PMSG, DFIG is more attractive as it is cheaper than PMSG and its power converter devices have a lower power rating. In [18], the authors propose two effective MPPT schemes by using direct power control, and perturbation and observation method, but the behaviors and limitations of thermodynamics and heat transfer of the prime mover have not been considered. In [19], it is found that the detailed modelling of the prime mover is important for a DS system as the control strategy of the prime mover can affect system stability when it is connected to the utility grid. However, the analysis and control system design presented is based on a simplified electrical model of DFIG, power converter, and grid system. In order to reduce the negative impact of temperature variation, a Takagi-Sugeno fuzzy supervisor is designed for several local state-feedback temperature controller to handle the nonlinearity of the system. As the control system is complex, the tuning of control parameters is cumbersome and no detailed discussion has been given.
To this end, in order to achieve the variable-speed operation of DS system with reduced effort in the control system design, this paper develops a new coordinated torque and mean pressure control scheme as well as a control parameter tuning method for MPPT and constant temperature regulation. The contribution of the work includes three aspects. First, a multi-physics dynamic model for a DFIG-based DS system has been established. For the first time, the detailed pneumatic model of the hydrogen supply system has been incorporated with the capability to consider the loss due to a gas compressor and limitation of external tanks. Second, a control scheme for a DFIG-based DS system Energies 2019, 12, 1300 3 of 23 is proposed to achieve simultaneous maximum power point tracking and temperature regulation with less complexity. Through steady-state analysis of the developed model, the optimal torque to for maximum power harnessing is calculated from the measured insolation, and then implemented using direct torque control, while a simple PI controller is used for temperature regulation. Third, to minimize the tracking error due to the disturbance of the temperature variation and the consumed power by the compressor, the optimal control parameters are determined using particle swarm optimization (PSO). With the incorporation of model order reduction and pre-examination of the system stability, the search effort can be greatly reduced.
The rest of the paper is organized as follows. Section 2 presents the overall system configuration and the modeling of the grid-connected DFIG-based DS system. Section 3 develops the control scheme with a simple control structure using direct torque control for active power regulation and using mean pressure control for temperature regulation. Section 4 proposes a method to use PSO to obtain the optimal control parameters with the consideration of system stability. Numerical examples are included in Section 5 to show the validity of the developed model and the effectiveness of the proposed control system. Section 6 concludes the main findings of the work.

System Configuration
The energy conversion of a typical DS system consists of three processes: First, the parabolic dish collects, reflects, and concentrates the direct normal irradiance of the sunlight into a small hollow chamber called the receiver. Second, a Stirling engine converts the high temperature heat accumulated in the receiver into mechanical energy. Finally, an electrical generator, driven by the Stirling engine, converts mechanical energy into electricity. Figure 1 shows the schematic of the DS system as well as the power and mass flows between/with the sub-systems.
Energies 2019, 12, x 3 of 23 system is proposed to achieve simultaneous maximum power point tracking and temperature regulation with less complexity. Through steady-state analysis of the developed model, the optimal torque to for maximum power harnessing is calculated from the measured insolation, and then implemented using direct torque control, while a simple PI controller is used for temperature regulation. Third, to minimize the tracking error due to the disturbance of the temperature variation and the consumed power by the compressor, the optimal control parameters are determined using particle swarm optimization (PSO). With the incorporation of model order reduction and preexamination of the system stability, the search effort can be greatly reduced. The rest of the paper is organized as follows. Section 2 presents the overall system configuration and the modeling of the grid-connected DFIG-based DS system. Section 3 develops the control scheme with a simple control structure using direct torque control for active power regulation and using mean pressure control for temperature regulation. Section 4 proposes a method to use PSO to obtain the optimal control parameters with the consideration of system stability. Numerical examples are included in Section 5 to show the validity of the developed model and the effectiveness of the proposed control system. Section 6 concludes the main findings of the work.

System Configuration
The energy conversion of a typical DS system consists of three processes: First, the parabolic dish collects, reflects, and concentrates the direct normal irradiance of the sunlight into a small hollow chamber called the receiver. Second, a Stirling engine converts the high temperature heat accumulated in the receiver into mechanical energy. Finally, an electrical generator, driven by the Stirling engine, converts mechanical energy into electricity. Figure 1 shows the schematic of the DS system as well as the power and mass flows between/with the sub-systems. In contrast to the use of a fixed-speed SCIG in a conventional DS system [1], the variable-speed DFIG-based DS system studied in References [18,19] incorporates a wound-rotor induction generator and a back-to-back ac/dc/ac converter to interconnect the rotor windings of the generator to the grid. Although with similar structure to commercially-available DFIG-based wind turbine generator (WTG), the power rating of the prime mover of DS system, i.e., the Stirling engine, is much lower than that of a wind turbine. Many such generating units would have to be aggregated in some suitable manner to realize a utility-scale power station. Figure 2 shows the schematic of a possible In contrast to the use of a fixed-speed SCIG in a conventional DS system [1], the variable-speed DFIG-based DS system studied in References [18,19] incorporates a wound-rotor induction generator and a back-to-back ac/dc/ac converter to interconnect the rotor windings of the generator to the grid. Although with similar structure to commercially-available DFIG-based wind turbine generator (WTG), the power rating of the prime mover of DS system, i.e., the Stirling engine, is much lower than that of a wind turbine. Many such generating units would have to be aggregated in some suitable Energies 2019, 12, 1300 4 of 23 manner to realize a utility-scale power station. Figure 2 shows the schematic of a possible configuration of such a power generating station. It is based on the use of a multi-string-converter, consisting of a common dc-link and a common grid-side converter (GSC) similar to that described in [22] for a photovoltaic plant.
Energies 2019, 12, x 4 of 23 configuration of such a power generating station. It is based on the use of a multi-string-converter, consisting of a common dc-link and a common grid-side converter (GSC) similar to that described in [22] for a photovoltaic plant.
. . .  Figure 2. Proposed configuration of a DFIG-based DS solar-thermal power plant using a common dc bus and a common GSC unit.

==
In Figure 2, the DFIG-based DS power plant is shown to consist of a number of clusters in which each cluster is consists of N power conversion units (PCU) and one common GSC unit. Each PCU contains a prime mover (the DS), a DFIG, a rotor-side converter (RSC), and a corresponding local control system. The dc terminals of RSCs from the N PCUs are connected to a common dc bus. The common GSC unit consists of a three-phase voltage source inverter, dc-link capacitor, the associated line filter, and a GSC controller. Similar to the observations made in [22], there are several distinct advantages in the proposed configuration. First, power flows can still be controlled through the individual PCU. However, as a GSC with a larger power-rating tends to have a lower specific cost ($/kW) than that of a lower-rating GSC in this power level, the multi-string-converter configuration is therefore economically attractive. Moreover, a smoothing effect due to the expected diversity of the insolation level (Ik) among the PCU in a cluster would yield smaller power variations on the common dc bus. This would facilitate the control of the common GSC. Enhanced functionality of GSC, such as with regard to power quality and reactive power control, can also be achieved more economically using the larger-capacity converter. Thus, with the configuration shown in Figure 2, each RSC is to achieve active power control of its PCU and to provide reactive power to the external grid. As long as the dc-link voltage is effectively maintained by the common GSC, the RSC can be individually controlled. The main drawback of the common GSC is that should it malfunction, the dc-link voltage may increase to some unacceptable level as the cluster is unable to export power to the grid. If this occurs, a braking resistor in the common brake chopper circuit in Figure 2 is inserted to dissipate the surplus energy and maintain the dc voltage. Other possible solution includes the use of an AC crowbar to short-circuit the rotor windings, and the DFIG then operates as a SCIG. However, such a contingency condition is not discussed in this article as this is considered a rare occurrence.
This work mainly studies the modeling and the local control of the PCU. The modeling of the power converters, dc bus, and grid system can be found in various literature such as [23].

Prime Mover Model
The prime mover of the PCU consists of the concentrator (i.e., the dish), the heat receiver/absorber, and the Stirling engine. It is shown as "DS" in Figures 1 and 2. The dish and the receiver are modelled using a normalized heat transfer equation [24]: In Figure 2, the DFIG-based DS power plant is shown to consist of a number of clusters in which each cluster is consists of N power conversion units (PCU) and one common GSC unit. Each PCU contains a prime mover (the DS), a DFIG, a rotor-side converter (RSC), and a corresponding local control system. The dc terminals of RSCs from the N PCUs are connected to a common dc bus. The common GSC unit consists of a three-phase voltage source inverter, dc-link capacitor, the associated line filter, and a GSC controller. Similar to the observations made in [22], there are several distinct advantages in the proposed configuration. First, power flows can still be controlled through the individual PCU. However, as a GSC with a larger power-rating tends to have a lower specific cost ($/kW) than that of a lower-rating GSC in this power level, the multi-string-converter configuration is therefore economically attractive. Moreover, a smoothing effect due to the expected diversity of the insolation level (I k ) among the PCU in a cluster would yield smaller power variations on the common dc bus. This would facilitate the control of the common GSC. Enhanced functionality of GSC, such as with regard to power quality and reactive power control, can also be achieved more economically using the larger-capacity converter. Thus, with the configuration shown in Figure 2, each RSC is to achieve active power control of its PCU and to provide reactive power to the external grid. As long as the dc-link voltage is effectively maintained by the common GSC, the RSC can be individually controlled. The main drawback of the common GSC is that should it malfunction, the dc-link voltage may increase to some unacceptable level as the cluster is unable to export power to the grid. If this occurs, a braking resistor in the common brake chopper circuit in Figure 2 is inserted to dissipate the surplus energy and maintain the dc voltage. Other possible solution includes the use of an AC crowbar to short-circuit the rotor windings, and the DFIG then operates as a SCIG. However, such a contingency condition is not discussed in this article as this is considered a rare occurrence.
This work mainly studies the modeling and the local control of the PCU. The modeling of the power converters, dc bus, and grid system can be found in various literature such as [23].

Prime Mover Model
The prime mover of the PCU consists of the concentrator (i.e., the dish), the heat receiver/absorber, and the Stirling engine. It is shown as "DS" in Figures 1 and 2. The dish and the receiver are modelled using a normalized heat transfer equation [24]: where I is the insolation level, T a is the ambient temperature, Q h is the rate of the heat flow absorbed by the Stirling engine, and the receiver temperature is assumed to be the same as the absorber/heater temperature T h of the Stirling engine. Constants T rec , K rec , and K con are as defined in [24], and their parametric values are determined according to the thermal characteristics of the receiver/absorber materials. For the study of the variable-speed operation of the DS system, the behaviors of a double-acting kinematic Stirling engine with four series-connected cylinders can be described by an average-value adiabatic model with single-cylinder equivalence and the correction of losses: where P m , τ m , and ω m are mechanical power, mechanical torque, and rotating angular speed of the Stirling engine, respectively. p mean and m se are the average mean pressure and the total mass of the hydrogen in the engine cylinder, respectively. p mean is proportional to m se , and the gain K p is a constant. gA se is the mass flow rate of the hydrogen from engine cylinder to hydrogen tanks. The constants a 00 , a 01 , a 10 , a 11 , b 00 , b 01 , b 10 , b 11 , b 02 , and b 12 are multivariable coefficients to represent the steady-state performance of the Stirling engine. Constants A and C are also defined in [24], and their parametric values are governed by the dimension and performance characteristics of the Stirling engine. These parameters shall also be affected by the working temperature T h . However, as T h is assumed to be well maintained at a maximum value T h,max , the effect of temperature variation on the these parameters is ignored. Note that all the parameters and variables given in Equations (1)-(4) are expressed in per unit, and the base values are selected and given in Table 1. The per unit parameters of the prime mover model used in this work are given in Table 2.

Hydrogen Supply System Model
The hydrogen supply system considered in this work consists of a high-pressure tank, a low-pressure tank, a gas compressor, and two solenoid valves, denoted as HV and LV, respectively. The thermodynamic behavior of the hydrogen stored in each of the pressure tanks is governed by Beattie-Bridgeman equation [25]. In Appendix A, the Beattie-Bridgeman equation is normalized to: where the gas pressure p X in the tank is expressed as a fourth degree polynomial function f X (·) of the hydrogen mass m X . The subscript X can be H or L, which denote the high-or low-pressure tank respectively. V X denotes the tank volume, and temperature T X is assumed to be the same as the cold sink temperature T k of the Stirling engine. Fast-responding solenoid valves HV and LV are used to regulate the mean pressure p mean of the hydrogen in the engine cylinder by controlling the mass exchange rate gA HV and gA LV respectively [11]. The control law, described by Equations (6) and (7), is such that if a positive valve command c is given, HV will open to allow hydrogen to flow into the cylinder while LV will close. This will increase both the mass and mean pressure of the hydrogen in the cylinder. Conversely, the reverse actions are taken for a negative command c. Thus, In Equations (6) and (7), K v and T v are the gain and the time constant of the solenoid valves, respectively. As p H > p mean > p L , hydrogen can only flow from the high-to the low-pressure tanks via the cylinder. As this flow continues, p H will be decreasing while p L will correspondingly increase. This is an unsustainable situation. A compressor is therefore used to maintain the gas pressures in the tanks at suitable levels. Considering the thermodynamic process to be a two-state polytropic compression, the power P comp consumed by the compressor can be expressed as [25]: where η comp is the overall efficiency of the compressor, n c is the polytropic constant, and the compressor temperature T comp is assumed to be the same as the cold sink temperature T k . gA comp is the mass flow rate of the gas in the compressor. With the compressor in operation, the mass balance equations are: where M tot is the total mass of the hydrogen in the DS system. The block diagram of the pneumatic model, containing Equations (5)- (12), of the hydrogen supply system is shown in Figure 3, and the parameters used in this work are given in Table 2. Note that the operation of the valves is limited by the maximum mass flow rate gA max , as indicated in Figure 3. that the operation of the valves is limited by the maximum mass flow rate gAmax, as indicated in Figure 3.

DFIG Model
For the time-domain simulation and control system design, a fifth-order induction machine model has been widely used [18,26]. Expressed in the synchronous reference frame, the voltage equations, flux linkage equations, and electromagnetic equation of this per-unit DFIG model are given as Equations (13)-(15), respectively: where vds and vqs are d-and q-axis stator voltages, respectively; vdr and vqr are d-and q-axis rotor voltages, respectively; ids and iqs are d-and q-axis stator currents, respectively; idr and iqr are d-and q-

DFIG Model
For the time-domain simulation and control system design, a fifth-order induction machine model has been widely used [18,26]. Expressed in the synchronous reference frame, the voltage equations, flux linkage equations, and electromagnetic equation of this per-unit DFIG model are given as Equations (13)-(15), respectively: T e = L m (i qr i ds − i dr i qs ) (15) where v ds and v qs are dand q-axis stator voltages, respectively; v dr and v qr are dand q-axis rotor voltages, respectively; i ds and i qs are dand q-axis stator currents, respectively; i dr and i qr are dand q-axis rotor currents, respectively; ψ ds and ψ qs are dand q-axis stator flux linkages, respectively; ψ dr and ψ qr are dand q-axis rotor flux linkages, respectively; R s and R r are the stator and rotor resistances; L s , L r , and L m are the stator self-inductance, rotor self-inductance, and mutual inductance, respectively; ω r is the shaft angular speed; and T e is the electromagnetic torque T e of the DFIG, respectively. Note that the generator convention is used and all the quantities are expressed in per-unit except the synchronous speed ω e = 2πf, where f is the nominal system frequency. Moreover, the motion equation of the DFIG is: where H is the inertia constant of a one-mass lumped-parameter shaft model of the DFIG-based DS system and D is the damping constant. T m is the mechanical torque of the DFIG. Note that the selected base values for angular speed and torque quantities of the DFIG per-unit system are different from those given in Table 1, and the following equations are used for unit conversion: where S N is the nominal apparent power rating of the DFIG, p n is the number of pole pairs of the DFIG, and r ω represents the gearbox ratio. Finally, taking the power P comp consumed by the compressor into consideration, the total active power generated by the PCU is thus: where P s and P r are the stator and rotor powers of the DFIG, respectively. The model parameters of DFIG used in this work are given in Table 3.

Local Control System of the PCU
The design of the local control system for each PCU is presented in this section. There are three main control objectives under consideration, including (1) active power control, (2) receiver/absorber temperature control, and (3) tank pressure control. The study of the interaction between PCU and the grid system, including the reactive power and voltage control, is not the focus of the present investigation. Hence, we assume the DS system is connected to a single-machine infinite bus system and the terminal voltage at the point of connection is ideally maintained at the nominal voltage.

Direct Torque Control
For active power regulation, direct torque control (DTC) can simplify the procedure of controller design and reduce the parameter dependence of the system compared to the conventional voltage-modulated control schemes [27]. This approach is employed by taking advantage of the main features of this control concept, i.e., lower machine parameters dependency, simpler structure, no coordinate transformation, high dynamic performances, and minimal response time, which achieves better torque control conditions than field-oriented control.
The schematic of the DTC is shown in Figure 4. First, the measured three-phase rotor and stator currents are transformed into the components i αr , i βr , i αs , and i βs in the two-phase stationary reference frame using a Clark transformation. The electromagnetic torque T e , the magnitude Ψ r and angle θ Ψr of the rotor flux linkage are then estimated using Equations (20)- (22), respectively: where ψ αs and ψ βs are calculated using: The estimated T e is compared with the torque reference T e * obtained from the optimal torque control (OTC), which will be discussed in Section 3.2. The difference ∆T e is used as the input of a three-level hysteresis comparator to generate the torque index K T . Similarly, the estimated Ψ r is compared with the flux reference Ψ r * to meet the requirement of reactive power control (which is not specifically discussed in this article). The difference ∆Ψ r is sent to a two-level hysteresis comparator to generate the flux linkage index K Ψ . Furthermore, θ Ψr is used to determine the sector number K s . The six gate signals for the RSC will be selected using a switching table based on the calculated combination of indices K T , K Ψ , and K s . The detail of the implementation of DTC can be found in [27] and it is not elaborated on in this article.
where ψαs and ψβs are calculated using: The estimated Te is compared with the torque reference Te * obtained from the optimal torque control (OTC), which will be discussed in Section 3.2. The difference ΔTe is used as the input of a three-level hysteresis comparator to generate the torque index KT. Similarly, the estimated Ψr is compared with the flux reference Ψr * to meet the requirement of reactive power control (which is not specifically discussed in this article). The difference ΔΨr is sent to a two-level hysteresis comparator to generate the flux linkage index KΨ. Furthermore, θΨr is used to determine the sector number Ks. The six gate signals for the RSC will be selected using a switching table based on the calculated combination of indices KT, KΨ, and Ks. The detail of the implementation of DTC can be found in [27] and it is not elaborated on in this article.  Figure 4. Schematic of the DFIG control.

Optimal Torque Control
The DTC introduced in the previous subsection decouples the regulation of the generated active power and the reactive power. By setting the torque reference Te * , the active power of the DFIG-based DS system can be regulated. In order to achieve MPPT, the steady-state analysis is conducted in this subsection to obtain the relationship between the optimal operating speed and the optimal torque.
First, using Equation (3), the steady-state mechanical torque τm of the engine can be expressed as:

Optimal Torque Control
The DTC introduced in the previous subsection decouples the regulation of the generated active power and the reactive power. By setting the torque reference T e * , the active power of the DFIG-based DS system can be regulated. In order to achieve MPPT, the steady-state analysis is conducted in this subsection to obtain the relationship between the optimal operating speed and the optimal torque. First, using Equation (3), the steady-state mechanical torque τ m of the engine can be expressed as: Substituting Equation (2) into Equation (1) to eliminate Q h , and setting dT h /dt = 0 and T h = 1 p.u. in the resulting equation, the steady-state mean pressure p mean of the working gas can be expressed as a function of the insolation level I and the engine speed ω m : Next, substituting Equation (25) into Equation (24) to eliminate p mean yields an expression of τ m as a function of insolation level I and engine speed ω m : Using Equation (26), a family of τ m versus ω m curves at different insolation levels can be obtained and plotted, as shown in Figure 5. It can be seen that at the same insolation level, there is a decreasing trend in τ m as the engine speed ω m increases. In order to achieve the maximum energy harness, one can evaluate the engine power P m = τ m ω m at every insolation level, and the theoretical MPPT or OTC curve on the τ m versus ω m plane can then be obtained and shown as H-K in Figure 5.
Using Equation (26), a family of τm versus ωm curves at different insolation levels can be obtained and plotted, as shown in Figure 5. It can be seen that at the same insolation level, there is a decreasing trend in τm as the engine speed ωm increases. In order to achieve the maximum energy harness, one can evaluate the engine power Pm = τmωm at every insolation level, and the theoretical MPPT or OTC curve on the τm versus ωm plane can then be obtained and shown as H-K in Figure 5. However, the practical operation of the DS system is restricted by several physical limits. In Figure 5, these limits are indicated by: a) the maximum mean pressure curve A-H (pmean = pmax = 1.0 p.u.), b) the minimum mean pressure curve D-G (pmean = pmin), c) the maximum insolation curve A-B (I = 1.0 p.u.), d) the minimum speed limit F-G (ωm = ωm,min), and e) the maximum speed limit B-C (ωm = ωm,max). Consequently, the feasible OTC curve becomes A-H-K-G, in which A-H and K-G are restricted by the maximum mean pressure of the working gas and the minimum engine speed, respectively. As such, it is expected that the DS system can track A-H-K-G in order to extract as much power when insolation varies.
When implementing A-H-K-G in the control system, the torque reference shall be determined using the measured speed quantities. However, it can be seen that the torque is not uniquely determined by the speed when the system is operating on the minimum speed line G-K. Therefore, respectively. As such, it is expected that the DS system can track A-H-K-G in order to extract as much power when insolation varies.
When implementing A-H-K-G in the control system, the torque reference shall be determined using the measured speed quantities. However, it can be seen that the torque is not uniquely determined by the speed when the system is operating on the minimum speed line G-K. Therefore, the tracking curve A-H-K-G is modified by replacing the constant speed line G-K with a ramp line G-T as shown in Figure 5. The revised feasible MPPT curve is thus A-H-T-G instead. The feasible optimal torque is thus expressed as a function f opt of engine speed: Using Equations (16)- (18) and Equation (27), one can obtain the steady-state T e versus ω r relationship that will be used in the DFIG control system, as shown in Figure 4:

Receiver/Absorber Temperature Control
Temperature control of a DS system is essential because it ensures the safe and efficient operation of the PCU. The receiver/absorber temperature T h is anticipated to be as high as possible to increase the thermal efficiency (which is defined as the energy conversion efficiency from the absorbed heat transfer rate Q h to the mechanical power P m of the Stirling engine) of the Stirling engine, while the material of the receiver/absorber will be damaged if the maximum tolerable temperature is exceeded [11].
As explained in [19], to overcome the problems caused by speed variation and system nonlinearity in a variable-speed DS system, a fuzzy supervisory control can be designed by taking advantage of the relatively slow variation of the insolation. Several local temperature controllers are designed using the pole-placement method by linearizing the system model at several selected operating points. The supervisor determines the control decision c of the solenoid valves. Although effective, a large set of control parameters have to be tuned. The selection of suitable positions of poles are empirical due to the complex control structure. In this subsection, we propose a simple temperature control scheme and only three control parameters need to be determined. The schematic of the proposed temperature control system is shown in Figure 6. increase the thermal efficiency (which is defined as the energy conversion efficiency from the absorbed heat transfer rate Qh to the mechanical power Pm of the Stirling engine) of the Stirling engine, while the material of the receiver/absorber will be damaged if the maximum tolerable temperature is exceeded [11]. As explained in [19], to overcome the problems caused by speed variation and system nonlinearity in a variable-speed DS system, a fuzzy supervisory control can be designed by taking advantage of the relatively slow variation of the insolation. Several local temperature controllers are designed using the pole-placement method by linearizing the system model at several selected operating points. The supervisor determines the control decision c of the solenoid valves. Although effective, a large set of control parameters have to be tuned. The selection of suitable positions of poles are empirical due to the complex control structure. In this subsection, we propose a simple temperature control scheme and only three control parameters need to be determined. The schematic of the proposed temperature control system is shown in Figure 6. As shown in Figure 6, instead of using the pole placement method for the design of a state feedback controller, a mean pressure control (MPC) scheme for the conventional fixed-speed DS system is used here, which is expressed as [14]: where * mean p is the set-point of the mean pressure of the working gas, and KMPC is the control gain of MPC. In order to reduce the negative impact of the disturbances of the varying insolation and speed, in the proposed temperature control scheme, the pressure set-point * mean p is expressed as: fb, fb, T ω p p Figure 6. Schematic of the proposed prime mover control.
As shown in Figure 6, instead of using the pole placement method for the design of a state feedback controller, a mean pressure control (MPC) scheme for the conventional fixed-speed DS system is used here, which is expressed as [14]: where p * mean is the set-point of the mean pressure of the working gas, and K MPC is the control gain of MPC. In order to reduce the negative impact of the disturbances of the varying insolation and speed, in the proposed temperature control scheme, the pressure set-point p * mean is expressed as: Here p * mean consists of three components: p fb,T is the output of the temperature PI controller, p OP is the term for feedforward compensation of the disturbance due to the insolation variation, and p fb,ω is the term to compensate for the disturbance due to the speed variation.
The optimal mean pressure (OP) function p OP (I) and the optimal speed (OS) function ω OS (I) depend on the steady-state p mean versus I and ω m versus I relationships, respectively, which can be predetermined analytically as follows. First, substituting Equation (27) into Equation (26) yields the ω m versus I relationship. Next, using this ω m versus I relationship and Equation (25) gives the p mean versus I relationship. These relationships are shown in Figure 7 and they can be implemented using two look-up tables.

Tank Pressure Control
Moreover, as shown in Figure 3, the compressor shall be used to maintain the tank pressure. The control law is expressed as: Hence, it can be seen that only four control parameters K1, K2, K3, and K4 are required to be determined in Equations (30) and (31). The selection of the optimal control parameters will be determined using the PSO method in the next section.

Constraints on Temperature Control Parameters
To determine the optimal control parameters, model order reduction is done first to simplify the stability analysis and to reduce the simulation time. First, to study the stability of the temperature control system, the effects of the compressor and pressure tanks are ignored by setting gAcomp = 0 and Pcomp = 0. In doing so, the hydrogen supply system model described using Equations (5)- (12) collapses to that shown in [19]: With Equations (32), (33), and (29), the transfer function GMPC(s) of the MPC is:

Tank Pressure Control
Moreover, as shown in Figure 3, the compressor shall be used to maintain the tank pressure. The control law is expressed as: Hence, it can be seen that only four control parameters K 1 , K 2 , K 3 , and K 4 are required to be determined in Equations (30) and (31). The selection of the optimal control parameters will be determined using the PSO method in the next section.

Constraints on Temperature Control Parameters
To determine the optimal control parameters, model order reduction is done first to simplify the stability analysis and to reduce the simulation time. First, to study the stability of the temperature control system, the effects of the compressor and pressure tanks are ignored by setting gA comp = 0 and P comp = 0. In doing so, the hydrogen supply system model described using Equations (5)- (12) collapses to that shown in [19]: where T p = 1/(K MPC K v K p ) = 4ξ 2 T v and ξ is the damping ratio. For ξ = 0.707, T p = 2T v . As T v << 1, the effect of the dynamics of the MPC system is ignored in a temperature control system design (i.e., p mean = p * mean ). Second, as the DTC of the DFIG provides very fast response to active power regulation, the torque is assumed ideally controlled by τ e = τ e * .
The resulting reduced order model made up of Equations (1) Figure 8, where the gains are defined as:

∂I
As can be seen in Figure 8, the plant for the temperature control design (including the prime mover, DFIG, MPC, and OTC) is a third-order single-input-single-output (SISO) system with one disturbance (insolation I).
Energies 2019, 12, x 13 of 23 As can be seen in Figure 8, the plant for the temperature control design (including the prime mover, DFIG, MPC, and OTC) is a third-order single-input-single-output (SISO) system with one disturbance (insolation I).  The characteristic equation of the system can be readily obtained as: a s a s a s   The characteristic equation of the system can be readily obtained as: where: Applying the Routh-Hurwitz stability criterion to Equation (35) gives the following Constraints (36)-(38): Constraints (36)-(38) will be used as constraints in the next subsection for particle swarm optimization to determine the optimal control parameters.

Particle Swarm Optimization
In the present investigation, the control objective is to maximize the energy generation, while maintain the operating temperature at a suitable level. Moreover, the hydrogen pressure in the tanks shall be maintained. For a given insolation profile I(t), we formulate the selection of control parameters as a finite horizon constrained optimization problem: subject to system model using Equations (1)- (19) and Constraints (36)-(38). In Constraint (39), the first term represents the total generated energy, ϕ 1 is used to implement a state constraint to limit the temperature within a narrow band|T h − T h,max | ≤ ∆T h,max , and ϕ 2 is used to introduce an additional cost to make sure that the final value of the tank pressure p H restores to its set-point. They are expressed as: where α and β are the penalty coefficients with large values. To obtain the cost J for each combination of control parameters (i.e., each particle), time-domain simulation has to be carried out. In order to reduce the searching time, the system model of Equations (1)- (19) will be simplified by assuming that the electromagnetic torque is ideally controlled using a DTC and no converter dynamics are considered, which yields τ e = τ * e . Next, representing Equation (16) in the per unit system of the prime mover and considering Equation (27), thus gives: where H and D are converted from H and D as given in Table 3. The PCU power Equation (19) is therefore modified to give: where η is a coefficient considering the losses in the DFIG windings, and it is expressed as a function of insolation I. Hence, the reduced order model for optimization algorithm design shall consists of Equations (1) In Equation (39), t max is the time horizon of the insolation profile. Parameters α and β are penalties for the evaluation of the performance of active power regulation and temperature control, respectively.
PSO is a simple and effective metaheuristic optimization technique that is initialized with some population of random candidates and a search algorithm to reach the global optimum for the defined objective function [28]. In an n-dimensional optimization problem, the position vector x ∈ R n and velocity vector v ∈ R n of the i-th particle in the iteration number j are expressed as: Each position vector represents a potential solution to the problem, while the velocity vector indicates the extent of variation in the updating procedure. The equations with which the velocity and the positions of the particles are updated are: where p is the recorded best solution for specific particle and g is the global optimal solution among all p. r and s are two random numbers between 0 and 1. w, c 1 , and c 2 are referred as the inertia constant, cognitive scaling parameter, and social scaling parameter, respectively. The recommended values c 1 = c 2 = 2 in [29] are used in this article and w is selected to be 0.8. These coefficients determine how much of each component from the previous iteration affects the next movement of the algorithm. For the problem in hand, the position vector is a four-dimensional vector: For the control scheme proposed in [19], there are 20 control parameters that need to be determined: 18 for the temperature controller and 2 for the speed controller. Hence, the number of control parameters have been significantly reduced.
In this article, the values of K 1 , K 2 , and K 3 shall satisfy Constraints (36)-(38). Hence, during the step to randomly select the particle, an additional step will be included to examine the stability requirement. Unaccepted particles will be abandoned before conducting the time-domain simulation.
However, as most of the parametric values of the gains in the linearized model depends on the insolation level I, the feasible range of K 1 , K 2 , and K 3 varies. To illustrate this point, in Figure 9 we plot an example of the feasible search subspace (K 1 , K 2 ) ∈ R 2 when K 3 = 0. It can be observed that as the insolation increases, the maximum K 1 also increases, which can be derived from Constraint (36). The other boundary also moves which is governed by Constraints (37) and (38). To guarantee system stability, the search subspace shall be shared area under all insolation levels. However, considering that all the gains change smoothly, only the cases for the maximum and the minimum insolation are investigated for simplicity. The feasible search subspace is therefore obtained and shown as the shaded area in Figure 9. insolation increases, the maximum K1 also increases, which can be derived from Constraint (36). The other boundary also moves which is governed by Constraints (37) and (38). To guarantee system stability, the search subspace shall be shared area under all insolation levels. However, considering that all the gains change smoothly, only the cases for the maximum and the minimum insolation are investigated for simplicity. The feasible search subspace is therefore obtained and shown as the shaded area in Figure 9.

Model Validation and Comparison
In order to validate the accuracy of the proposed model, in this subsection, the simulation results from the developed model were compared with that obtained using the existing models. All the models were implemented and solved using MATLAB R2016a/Simulink 8.7 (The MathWorks, Inc., Natick, MA, USA). The parameters of the system model used in the simulation are given in Tables 1-3. The optimal control parameters obtained in the latter subsection are used.
First, for the Stirling engine, a multi-cylinder adiabatic model of the Stirling engine was used as the benchmark to evaluate the model accuracy of proposed model of Equations (1)- (19). This Figure 9. Feasible search subspace for PSO to determine K 1 and K 2 for K 3 = 0.

Model Validation and Comparison
In order to validate the accuracy of the proposed model, in this subsection, the simulation results from the developed model were compared with that obtained using the existing models. All the models were implemented and solved using MATLAB R2016a/Simulink 8.7 (The MathWorks, Inc., Natick, MA, USA). The parameters of the system model used in the simulation are given in Tables 1-3. The optimal control parameters obtained in the latter subsection are used.
First, for the Stirling engine, a multi-cylinder adiabatic model of the Stirling engine was used as the benchmark to evaluate the model accuracy of proposed model of Equations (1)- (19). This adiabatic model describes the thermodynamic cycles of four series-connected cylinders, and has been studied intensively in existing literature [11,12,21]. In this case, detailed representation of the converters and the grid system were established using the SymPowerSystems Toolbox. Figure 10 shows the simulation results when a step increase of insolation from 0.6 to 0.65 p.u. occurred at 0.5 s. The purpose was to compare the dynamic response of the models that could affect the stability of the system. It can be seen that using the average-value Stirling engine model, the dynamic behavior of the benchmark model could be successfully captured, while the high-frequency oscillating components in Figure 10c had an ignorable effect on other variables. This is in line with the conclusion of Reference [24]. Furthermore, the simulation results using an assumed prime mover model used in [18] are also presented. This model assumes the prime mover of the DS system is a first-order system and it ignores the detailed description of the thermodynamics of the Stirling engine and temperature control system. The mechanical torque was calculated according to the measured rotational speed of DFIG and insolation, according to the relationship presented in Figure 5. Corresponding temperature and absorbed heat transfer rate are not shown in Figure 10a,c as these variables are not considered in this assumed model. Clearly, it cannot reflect the dynamic characteristics of the DS system for the design of the temperature control system.
The effect of the tank volume V X was investigated next and the proposed hydrogen supply system model was compared with the simplified model used in References [12,19] (i.e., Equations (32) and (33)). The simplified model ignores the detailed representation of the tank and compressor, which leads to V X = +∞ and P comp = 0. In this case, the initial insolation level was 0.3 p.u., and it started to increase at 5 s with a ramp rate of 30 W/(m 2 ·s) or 0.03 p.u./s. This value has been considered to be very fast for solar power generation studies [14]. Figure 11 shows the simulations results. A low-pass filter has been applied to the PCU power to remove the high-frequency components. From Figure 11b,d, it can be seen that for V X = 5.5 p.u., the gas pressure in the low pressure tank reached zero during the transient, and a large amount of power was consumed by the compressor from 5 s to 20 s. By increasing the tank volume, the P comp profile became flatter. For V X = 500 p.u., the simulation results were close to that of the simplified model. Hence, when designing the optimal control scheme, the effect of the pressure tank and gas compressor should be carefully considered to improve the model accuracy.
the benchmark model could be successfully captured, while the high-frequency oscillating components in Figure 10c had an ignorable effect on other variables. This is in line with the conclusion of Reference [24]. Furthermore, the simulation results using an assumed prime mover model used in [18] are also presented. This model assumes the prime mover of the DS system is a first-order system and it ignores the detailed description of the thermodynamics of the Stirling engine and temperature control system. The mechanical torque was calculated according to the measured rotational speed of DFIG and insolation, according to the relationship presented in Figure 5. Corresponding temperature and absorbed heat transfer rate are not shown in Figure 10a,c as these variables are not considered in this assumed model. Clearly, it cannot reflect the dynamic characteristics of the DS system for the design of the temperature control system. The effect of the tank volume VX was investigated next and the proposed hydrogen supply system model was compared with the simplified model used in References [12,19] (i.e., Equations (32) and (33)). The simplified model ignores the detailed representation of the tank and compressor, which leads to VX = +∞ and Pcomp = 0. In this case, the initial insolation level was 0.3 p.u., and it started to increase at 5 s with a ramp rate of 30 W/(m 2 ·s) or 0.03 p.u./s. This value has been considered to be very fast for solar power generation studies [14]. Figure 11 shows the simulations results. A low-pass filter has been applied to the PCU power to remove the high-frequency components. From Figure  11b,d, it can be seen that for VX = 5.5 p.u., the gas pressure in the low pressure tank reached zero during the transient, and a large amount of power was consumed by the compressor from 5 s to 20 s. By increasing the tank volume, the Pcomp profile became flatter. For VX = 500 p.u., the simulation results were close to that of the simplified model. Hence, when designing the optimal control scheme, the effect of the pressure tank and gas compressor should be carefully considered to improve the model accuracy. Furthermore, the proposed control scheme described in Section 3 was compared with an existing scheme proposed in [19]. In [19], MPPT was achieved using a classic double-loop control structure, with the outer feedback speed control loop and the inner current control loop. Space-vector pulse width modulation (SVPWM) was used to generate the gate signal for the converter. Detail of such a scheme can be found in most textbooks and thus it is not elaborated on here. Figure 12 shows the simulation results. Again, low-pass filters have been applied to the PCU power and electromagnetic torque. Furthermore, the proposed control scheme described in Section 3 was compared with an existing scheme proposed in [19]. In [19], MPPT was achieved using a classic double-loop control structure, with the outer feedback speed control loop and the inner current control loop. Space-vector pulse width modulation (SVPWM) was used to generate the gate signal for the converter. Detail of such a scheme can be found in most textbooks and thus it is not elaborated on here. Figure 12 shows the simulation results. Again, low-pass filters have been applied to the PCU power and electromagnetic torque. different tank volumes: (a) power generated by the PCU, (b) power consumed by the gas compressor, (c) gas pressure in the high pressure tank, and (d) gas pressure in the low pressure tank. Furthermore, the proposed control scheme described in Section 3 was compared with an existing scheme proposed in [19]. In [19], MPPT was achieved using a classic double-loop control structure, with the outer feedback speed control loop and the inner current control loop. Space-vector pulse width modulation (SVPWM) was used to generate the gate signal for the converter. Detail of such a scheme can be found in most textbooks and thus it is not elaborated on here. Figure 12 shows the simulation results. Again, low-pass filters have been applied to the PCU power and electromagnetic torque.  It can be seen from Figure 12b that for both schemes, the receiver temperature can be effectively regulated around the set-point 1033 K. However, Figure 12a shows that the proposed control scheme provides a smoothened power profile as insolation changes. The reason is explained briefly as follows. For the feedback speed control, in order to track the maximum power points during the insolation increase, an accelerating torque has to be provided. In Figure 12c, a higher ramp rate of speed is observed from 20 s to 30 s, compared to that before 20 s. This was due to the steady-state characteristics of the MPPT curve given in Figure 5. Hence, during this period, the control system had to reduce the electromagnetic torque of the DFIG so that more accelerating torque could be provided, and this could be observed from Figure 12d. Meanwhile, for the proposed DTC scheme, the torque is directly regulated according to the measured insolation, such that an unexpected decrease in electromagnetic torque can be effectively avoided, resulting in a flatter speed profile, as shown in Figure 12c.

Case Study
This section will present a case study where the simulation results were obtained from MATLAB/Simulink software. In this case study, the following 100-second insolation profile I(t) was used for PSO: the initial insolation was 0.25 p.u. and the insolation was ramped up to the maximum insolation level 1.0 p.u. and then ramped down back to 0.25 p.u. The global optimal solution and the corresponding cost for each iteration are shown in Figure 13. Here, the particle size was 20. From Figure 13, it can be seen that the result converged and became stable at the 18th iteration, and the optimal control parameters were found to be K 1 = 11.1, K 2 = 11.8, K 3 = 0.18, and K 4 = 25.5. Figure 14 shows the positions of the randomly generated particles using PSO at iteration 1, 2, 11, 12, 20, and 21. Note that only two dimensions (K 1 , K 2 ) are shown. It can be seen that the generated particles were approaching the optimal position as iterations increased, and all the particles fell within the feasible search subspace shown in Figure 9.
MATLAB/Simulink software. In this case study, the following 100-second insolation profile I(t) was used for PSO: the initial insolation was 0.25 p.u. and the insolation was ramped up to the maximum insolation level 1.0 p.u. and then ramped down back to 0.25 p.u. The global optimal solution and the corresponding cost for each iteration are shown in Figure 13. Here, the particle size was 20. From Figure 13, it can be seen that the result converged and became stable at the 18th iteration, and the optimal control parameters were found to be K1 = 11.1, K2 = 11.8, K3 = 0.18, and K4 = 25.5.  Figure 14 shows the positions of the randomly generated particles using PSO at iteration 1, 2, 11, 12, 20, and 21. Note that only two dimensions (K1, K2) are shown. It can be seen that the generated particles were approaching the optimal position as iterations increased, and all the particles fell within the feasible search subspace shown in Figure 9. In order to show the efficacy of the obtained optimal control parameters, in Figure 15 we compare the simulation waveforms of the temperature Th, PCU power PPCU, compressor power Pcomp, and the tank pressures pH and pL using four groups of different control parameters, including the optimal one and three other groups of parameters. As each group of the parameters can be considered as a particle in the PSO search space, we denoted the control parameters as Particle 1, Particle 2, and Particle 3, respectively. The parametric values of the three particles were selected to be:  In order to show the efficacy of the obtained optimal control parameters, in Figure 15 we compare the simulation waveforms of the temperature T h , PCU power P PCU , compressor power P comp , and the tank pressures p H and p L using four groups of different control parameters, including the optimal one and three other groups of parameters. As each group of the parameters can be considered as a particle in the PSO search space, we denoted the control parameters as Particle 1, Particle 2, and Particle 3, respectively. The parametric values of the three particles were selected to be: • Particle 1: K 1 = 20, K 2 = 11.8, K 3 = 0.18, and K 4 = 25.5. In order to show the efficacy of the obtained optimal control parameters, in Figure 15 we compare the simulation waveforms of the temperature Th, PCU power PPCU, compressor power Pcomp, and the tank pressures pH and pL using four groups of different control parameters, including the optimal one and three other groups of parameters. As each group of the parameters can be considered as a particle in the PSO search space, we denoted the control parameters as Particle 1, Particle 2, and Particle 3, respectively. The parametric values of the three particles were selected to be:  From Figure 15a, it can be seen that the temperature profiles from all four selected particles were limited within a narrow band between 1030 K (0.997 p.u) to 1036 K (1.003 p.u.). For the first particle, the system was unstable and the temperature was oscillating after the system reached the steady state after about 77 s. In the steady state, the compressor did not work and no power was be consumed by the compressor. However, it can be seen from Figure 15c that for the first particle, there is a persistent compressor power. The pressure in the high-pressure tank decreased and the pressure in the low-pressure tank increased, which can be observed from Figure 15d,e, respectively. Figure 15c shows the second particle consumed the least amount of compressor power, which helped to increase the total generated power from the PCU. However, from Figure 10e, it can be seen that at 100 s, the gas pressure in the high-pressure tank had not reached the set-point. For the third particle, although it provided a better control on the high pressure tank, it consumed more compressor power during the period of investigation.

Conclusions
This paper proposes an optimal control system design approach for a DFIG-based dish-Stirling concentrated solar-thermal generating system, aiming to achieve maximum power point tracking and constant temperature regulation. The overall system model with multi-physics nature has been developed for time-domain simulation studies in a grid-connected condition. The model incorporates an average-valued thermodynamic model of Stirling engine, a heat transfer model of the dish/receiver, and a pneumatic model of the hydrogen supply system. A coordinated torque and mean pressure control scheme were proposed to provide smoothened active power tracking and temperature regulation. This scheme requires tuning for only four control parameters, which is a significant reduction from 20 in the existing literature. The optimal control parameters were determined using the PSO method, and the process of the pre-examination of system stability was incorporated into the algorithm to effectively reduce the search effort. The results from the time-domain simulation were presented to show the efficacy of the proposed scheme in supplying the grid system with improved performance of active power and temperature regulation as solar irradiance varied. It also showed the importance of the detailed modeling of the hydrogen supply system as the parasitic loss due to the gas compressor and the limited tank volume can play negative impact on the energy conversion efficiency, which should be taken into consideration during the design of the DS system.

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

D, H
damping and inertia constants of DS system (p.u. in electrical system) D , H damping and inertia constants of DS system (p.u. in mechanical system) I(t) solar insolation (p.u.) K 1 , K 2 proportional and integral gains of the temperature controller K 3 gain of the speed feedforward controller K 4 gain of the hydrogen compressor controller K v , T v gain and time constant of solenoid valve K con gain of the concentrator K rec , T rec gain and time constant of receiver K p engine pressure-mass gain M tot total mass in the engine cylinders P e (t), P m (t) electromagnetic power, mechanical power (p.u.) P comp (t) power consumed by the hydrogen compressor (p.u.) P PCU (t) power generated from PCU (p.u.) Q h (t) heater and cooler heat transfer rate (p.u.). T a ambient temperature (p.u.) T e (t), T m (t) electromagnetic and mechanical torques (p.u. in electrical system) T h (t) heater/receiver temperature (p.u.) a 00 , a 01 , a 10 , a 11 b 00 , b 01 , b 10 , b 11 , b 02 , b 12 steady-state efficiency correction coefficients of Stirling engine c(t) control command of solenoid valve (p.u.) gA se (t) total mass flow rate from tanks to engine cylinder (p.u.) m se (t) mass of the hydrogen in the engine cylinder m H (t), m L (t) masses of the hydrogen in the high-and low-pressure tanks (p.u.) p mean (t) mean pressure of the working gas in the Stirling engine cylinder (p.u.) p max , p min maximum and minimum pressure in the Stirling engine cylinder (p.u.) p H (t), p L (t) pressures of the hydrogen in the high-and low-pressure tanks (p.u.) τ e (t), τ m (t) electromagnetic and mechanical torques (p.u. in mechanical system) ω m (t) shaft speed of engine/generator (p.u. in mechanical system) ω r (t) shaft speed of engine/generator (p.u. in electrical system)

Appendix A
The original Beattie-Bridgeman equation is given in [25] is: where n = m/M is the amount of the substance, and m is the mass and M is the molar mass of hydrogen. The units of the states of the equation are: p (atm), n (mol), m (kg), V (L), T (K). By using the base values given in Table 1, Equation (A1) can be normalized to: where the hydrogen constants are: Re-arranging Equation (A2) yields Equation (5), where the coefficients β k are: The hydrogen constants used in Equation (A1)