Digital twin of a Fresnel solar collector for solar cooling

Applied Energy 339 (2023) 120944 Available online 17 March 2023 0306-2619/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/bync/4.0/). Contents lists available at ScienceDirect Applied Energy journal homepage: www.elsevier.com/locate/apenergy Digital twin of a Fresnel solar collector for solar cooling Diogo Ortiz Machado a,b,c, William D. Chicaizac, Juan M. Escañoc, Antonio J. Gallego c, Gustavo A. de Andrade b, Julio E. Normey-Ricob, Carlos Bordons c,d,∗, Eduardo F. Camachoc,d aIFRS - Instituto Federal de Educação, Ciencia e Tecnolgia do Rio Grande do Sul, Rua Alfredo Huch, 475, Rio Grande, 96201 460, Rio Grande do Sul, Brazil bUFSC - Universidade Federal de Santa Catarina. Departamento de Automação e Sistemas, R. Eng. Agronômico Andrei Cristian Ferreira, Florianópolis, 88040 900, Santa Catarina, Brazil cUS - Universidad de Sevilla. Departamento de Ingeniería de Sistemas y Automática, Camino de los Descubrimientos, Sevilla, 41092, Andalucía, Spain dENGREEN - Laboratory of Engineering for Energy and Environmental Sustainability, Universidad de Sevilla, Spain ARTICLE INFO Dataset link: https://data.mendeley.com/datas ets/rzggrvczf6 Keywords: Adaptive neuro fuzzy inference system Distributed parameter dynamic model Dynamic modeling Validation Absorption plant ABSTRACT This work develops digital entities of a commercial Fresnel Solar Collector (FSC) installed in an absorption cooling plant. The objective is to create and validate models that describe the FSC dynamics across its whole operation range during the day and the night. Thus, the temperatures range between operation temperature of 180 ◦C and almost ambient temperature due to overnight heat losses. In the same sense, the flow range between zero to 13m3∕h. The idea is that the digital twin will aid start-up and shut-down optimization and control design reliability. The paper employs two modeling approaches, then evaluates their twinning/adaptation time and performance validation. One model uses phenomenological modeling through Partial Differential Equations (PDE) and parameters identification, and another uses a data-driven technique with Adaptive Neuro-Fuzzy Inference Systems (ANFIS). The available measurement data sets comprise 25 days of operation with a sampling time of 20 s which, after outlier removal, filtering and treatment, resulted in 108416 samples. The validation considers six separate operating days. Results show that both models can twinning/adapt considering measured data. The models present pretty good results and are suitable for control and optimization. Besides, this is the first paper considering the FSC mirror defocus action on dynamic modeling and validation. 1. Introduction The solar power that daily strikes planet earth is the driving energy that sustains all life that evolved to an intricate and delicate equilibrium. Therefore, knowing how to harness this power source in a practical and versatile manner will pave the way for a sustainable future. This work contributes to developing two adaptive models of a Fresnel Solar Collector (FSC) in the framework of digital twins. One uses data-driven neuro-fuzzy (NF) networks, and the other uses phenomenological Partial Differential Equations (PDE) with parameter identification. Both digital entities will plan, integrate, and control the absorption plant installed at the Escuela Tecnica Superior de Ingenieria - ETSI, Seville, Spain, to increase renewable energy use in future studies. Policies have been increasing investments and knowledge development to solve the energy problem, considering fossil fuel burning continues to worsen the climate crisis. The reason is that CO2emissions from the energy and industry sectors have increased by 60% despite the United Nations Framework Convention on Climate Change in 1992 [1]. The energy sector is the primary source of global emissions, accounting ∗Corresponding author at: US - Universidad de Sevilla. Departamento de Ingeniería de Sistemas y Automática, Camino de los Descubrimientos, Sevilla, 41092, Andalucía, Spain. E-mail address: [email protected] (C. Bordons). for approximately 60% of global greenhouse gas emissions [2], and the building sector uses 40% of the world’s energy production, which 5% is for cooling. Furthermore, the cooling demand tends to grow due to a hotter climate [3]. In this line, investments of 7.4 trillion euros are estimated, in the next 25 years, for the deployment of technologies that eliminate net CO2emissions [3]. Therefore, solar cooling technologies development is a clever way to reduce the CO2emissions while enhancing cooling plants’ technological maturity and economic viability [4]. Thus, solar cooling technology has received much attention [5]. Solar absorption plants produce cold from a solar-heated source through an absorption thermodynamic cycle [6]. It has the feature of having solar energy availability following the cooling demand. The problem is that using solar energy as the primary source adds the complexity of having a process that is intermittent during the day, night, and year seasons. Besides, solar irradiance that strikes the absorber is subject to clouds that strongly affect the plant’s dynamic behavior. Thus, generating a stable heat source for the absorption chiller is critical. https://doi.org/10.1016/j.apenergy.2023.120944 Received 4 August 2022; Received in revised form 9 December 2022; Accepted 5 March 2023 Applied Energy 339 (2023) 120944 2 D.O. Machado et al. Several concentrating solar collector types harvest thermal solar energy. Parabolic Trough Collectors (PTC) and Fresnel Solar Collectors (FSC) are among the line focus systems. Despite PTC’s technological maturity and higher efficiency, the FSC has advantages compared to PTC [7], mainly cheaper production resulting in equivalent Levelized Cost of Electricity. Besides, FSC has room for further technological development [5,7]. In this way, improving the technological/economic performance of these systems [8] is highly desirable [9]. This work seeks to improve FSC operation through its Digital Twin formulation for future control and optimization studies. From the whole plant control perspective, a hierarchical structure typically solves the scheduling (integration), optimizing (plan), and controlling problems, separating them into layers considering the computation effort and operation horizon of each one [10]. The upper optimization layers minimize a given cost function and send the optimal schedule or operation points solution to the control layer. The control layers seek stable operation, reference tracking, and disturbance rejection. The optimization layers typically use the Real-Time Optimization (RTO) technique [11], while the multivariable control layer uses Model-based Predictive Control (MPC) techniques [12]. Both run models continuously in a feedback loop. Therefore, the model’s accuracy directly impacts the control, optimization, thus, the process operation. Specifically for concentrating solar collectors that operate intermittently, with a wide operational range, and present slow and fast dynamics, a typical RTO with non-linear steady-state models and a linear MPC with dynamic models presents limitations [13–15]. With this, the dynamic models resulting from this work are accurate for both steady-state and dynamic optimizations with low computational burden, enabling assets for control and optimization [16,17]. These models facilitate more advanced techniques realizations, feasible in shorter sampling times, such as DRTO [18], Economic MPC [19], or hybrid approaches [20]. A Digital Twin (DT) has many definitions. In a review study, Rasheed et al. [21] state that a Digital Twin is a ‘‘virtual representation of a physical asset enabled through data and simulators for realtime prediction, optimization, monitoring, controlling, and improved decision making’’. Although the concept, its definition, and the related studies continue to evolve in many areas [22]. Digital Twin is one of the most promising enabling technologies for Industry 4.0 viability, already having applications in the industry through publications, patents, and best practices of leading companies. The applications cover product lifecycle, product design, reliable/flexible production, prognostics, and health management [23]. The application of DT in the energy sector is relatively scarce. Rasheed et al. [21] compiles DT contributions to the energy sector, having studies on asset maintenance, energy saving, improving efficiency in smart factories aiming to reduce both production costs and greenhouse gas emissions, design, construction, and performance of residential buildings. Although, only one study published a full-scale DT of a district heating and cooling network. As can be seen, a core DT enabling technology is the modeling, simulation, verification, and validation. This paper develops both physical and data-driven models of the Fresnel Solar Collector (FSC) of Escuela Técnica Superior de Ingeniería de Seville (ETSI), Spain, to compose its DT. A validation states their performances considering the computational execution time and statistical indexes. The plant location is on the roof of the ETSI building. It aims to supplement the air conditioning system with chilled water to reduce electric consumption, CO2emissions, and operating costs [24]. The plant is a multi-energy system once it transforms solar irradiance into internal thermal energy, hot water into chilled water (thermal to thermal), and gas chemical energy into internal thermal energy. This system is also complex, highly non-linear, and dynamic, switching between electric, gas, and solar resources according to meteorology and demand profiles. The objective is to integrate and enhance the energy system with control [25], and optimization techniques [26]. It is worth noting that the ETSI plant is a physical entity. It has computers, servers, and industrial communication structures to connect virtual and physical spaces. Thus, the only asset not available for the ETSI absorption plant DT is its adaptive virtual entity or its adaptive dynamic model. This work seeks to expand the ETSI FSC model with the following specifications: 1. To consider the transitory regime, the continuous focus/defocus of the mirrors, and the cleanliness factor. 2. To have defined accuracy and precision indexes through validation with a representative amount of data. 3. To run fast enough to be used in Model Predictive Control (MPC) techniques. 4. To describe transients, and part-load, during the day and night for control, optimization, and what-if decisions. 5. To adapt and cope with the aging of the FSC for updated plant operation, long-term, and lifecycle assessment. These specifications have the following justification. Describing the solar collector considering day and night is especially important because solar plants are intrinsically intermittent and have appreciable heat losses during the night. Thus, this model feature can aid start-up/shut-down decision-making, contributing to optimal thermal storage strategy formulation [24]. The focus/defocus action is critical for the safety [27] and control [28] of FSC. However, to the best of the authors’ knowledge, there is no FSC validated model with the defocus feature. The cleanliness factor is related to thermal efficiency and varies over time. Thus, a model that can adapt such parameters is necessary to predict when to stop the FSC for maintenance and cleaning [27]. Stating the model’s accuracy and precision with mean error and standard deviation, among other statistical indexes, gives information about the models’ confidence and the fidelity level of DT’s simulations and what-if analysis [21,29–31]. The updating/adapting procedure of the virtual entities/models is called twinning. The twining is synchronizing the virtual and physical entities [29], measuring the physical entity state, and updating the virtual entity to reduce the difference between them. This process can occur either from physicalto-virtual or virtual-to-physical spaces, in a closed loop, with a twinning rate. To the authors’ best knowledge, no published FSC dynamic model is available considering such features. A specific literature review about the available modeling of FSC on the ETSI absorption plant leads to the following works: Robledo et al. [32] develop a phenomenological lumped parameter dynamic model considering the optical and thermal models and an ordinary differential equation. The authors used least-squares methods to identify the model’s parameters. The validation considered four operation days data, with a 7.5 h duration each. The validation outlet temperature range is from 85 to 165 ◦C, and the results are qualitatively presented, contrasting the actual data and model output plots. Spoladore et al. [33] develop a phenomenological, distributed parameters, dynamic model, considering the optical and thermal models and Partial Differential Equations (PDE). The authors used leastsquared methods to identify the global heat loss and the metal-fluid heat transmission coefficients considering second-degree polynomial functions. The validation considered two operation days data from 11:00 to 18:00, on May 27, 2010, and November 17, 2009. The outlet temperature ranges from 100 to 180 ◦C, with a maximum error of 8 ◦C. The authors present qualitative comparison plots between model prediction and actual data. Pino et al. [27] develop phenomenological, steady-state, thermal, and optical models using algebraic equations that account for the focus of the mirrors. The validation considers the mirror’s rows inclinations, absorbed heat, and outlet temperature during operation between 13:00 and 15:30 on May 27, 2009. The validated outlet temperature range is from 147 to 168.5 ◦C, with relative errors of less than 1%. Chicaiza et al. [34] develop data-driven dynamic models based on Adaptive Neuro-Fuzzy Inference System (ANFIS) to compose a DT. The Applied Energy 339 (2023) 120944 3 D.O. Machado et al. ANFIS training considers two days of actual and artificial data from PDE models. The validated outlet temperature range is from 40 to 180 ◦Cin three days of actual data. The authors present qualitative plots between measured and model outlet temperatures and calculate the Root Mean Squared Error (RMSE) and the Mean Absolute Percentage Error (MAPE) for each day. The worst-case RMSE =12.96 ◦Cand MAPE =8.54%, while the best-case RMSE =3.67 ◦Cand MAPE =2.27%. Considering the previous works on the ETSI FSC modeling, this paper has the following contributions. 1. It develops dynamic models using both ANFIS [35] and PDE [36] to generate a transparent model with explicit equations and expand FSC models capabilities. This work uses 25 days of operation data with a sampling time of 𝑡𝑠=20 s, with continuous data between day and night. 2. The ANFIS twinning (training and checking) and validation consider 19 days of measured data. On the other hand, identifying PDE model parameters uses four days of measurement data resampled every 5 s to avoid numerical integration instability. The massive amount of data in twinning results in generalized dynamic models of both ANFIS and PDE. 3. The validation procedure considers six days of actual operation data, three in June and three in October. The validation results indicate that the models can adequately represent the behavior of the FSC outlet in a wide temperature range, from 40 to 180 ◦C. 4. The models describe operations continuously during day and night, part-load operations, or overnight, with heat losses. 5. To the best authors’ knowledge, this work presents the first FSC-validated dynamic models that describe the mirror’s focus/defocus action. 6. Lastly, this work contributes to testing the computational speed of resulting dynamic models during twinning and simulation, setting the execution time limits of the models regarding MPC and dynamic optimization techniques. The organization of the rest of the paper is as follows. Section 2 defines the FSC process and presents the actual data preparation and its correlation analysis, closing with the FSC operation description. Section 3presents the FSC phenomenological distributed parameters modeling and its identification. Section 4states the ANFIS architecture and defines its training procedure. Section 5defines the simulation planning with the model structure definition and data sets preparation. Section 6states the results divided into two parts, the first is the twinning time of both models (training and identification), and the second is devoted to validation results and indexes considering actual data. Lastly, Section 7shows the findings of this work and closes the paper. 2. Fresnel solar collector process The Fresnel solar collector (FSC) installed at ETSI was constructed by PSE AG, which activities have passed to Mirrox Gmbh, and today is Industrial Solar [37]. Fig. 1 depicts the referred Fresnel concentrating solar collector that generates heat for an absorption chiller that, in its turn, supplements the building air conditioning system with a renewable primary energy source. The solar field has an 18◦orientation in the east–west direction, a total area of 352 m2, where 11 mirror rows and 16 modules compose the optical system, summing 64 m of length, which focus the solar irradiance to the receiver as depicted in Fig. 2. The receiver, in its turn, is composed of a secondary reflector and an absorber tube, where glass thermal insulation equips both. Steel DNI 1.4541 (AISI321) composes the SCHOTT PTR70 [38] absorber tube, installed 4 meters above the mirror’s plane. Water is the heat transfer fluid and flows inside the absorber tube, with 13 bar and 180 ◦Cof nominal pressure and temperature, respectively. Table 1 presents the main characteristics of the Fresnel solar collector. Table 1 ETSI Fresnel solar collector characteristics [24,38]. Parameter (symbol) Value Unit Total area 512 m2 Total mirror area (𝐴𝑡) 352 m2 Absorber length 64 m Absorber height 4 m Mirror unitary length 4 m Mirror unitary width 0.5 m Total Mirror aperture (𝐺) 5.5 m Mirrors rows 11 – Number of mirrors 176 – Absorber tube specific mass (𝜌𝑚) 8027 kg∕m3 Absorber tube external diameter 0.07 m Absorber internal diameter 0.066 m Absorber tube specific heat (𝑐𝑚) 500 J∕(kg ◦C) Receiver cavity aperture 0.3 m Heat transfer fluid Water – Nominal temperature 180 ◦C Nominal pressure 13 bar Nominal thermal power 120 kW The FSC objective is to generate a hot outlet temperature flow. The incident irradiance and the inlet temperature are critical disturbances of the solar collector process because they fastly and widely affect the outlet temperature. The plant start-up occurs when the solar irradiance reaches a minimum, starting the pump, mirror’s sun tracking, and focusing the solar beam into the absorber tube. After this beginning, the controller manipulates the flow and the mirror’s focus to regulate the outlet temperature of the collector with the objectives of reference tracking and disturbance rejection. The plant’s shut-down occurs when irradiance reaches a given minimum value. 2.1. Data preparation Neuro-fuzzy (NF) and Partial Differential Equations (PDE) models are identified and validated using actual data for predicting the FSC dynamic behavior. Therefore, an identified model is only as good as the data used to fit/train it [39]. The impact of feeding the identification system and validating with poor-quality data is low accuracy and veracity, leading to pointless work and unuseful results that, if used, would result in economic and safety losses [40]. Thus, data preparation means the model is built right [41]. Data quality assessment uses accuracy, completeness, consistency, and timeliness classifications [42]. This work executes the following preparation seeking these classifications: filtering instrumentation noise, homogenizing sampling times, excluding data garbage, and carefully choosing variables that affect the process. The preparation strongly affects the model’s accuracy in the results section. This work employs raw data from the solar absorption plant’s SCADA system that stores measurements in Excel data sheets daily. Each file has 18 variables (columns) imported to Matlab as timetables with the respective headers or variable names. An algorithm tests each Excel sheet header to maintain the same variables’ order names to avoid mixing different variables in one array. Then each available day of operation was concatenated by time, considering consecutive days. Data inspection indicates incomplete data, missing measurement values, and irregularly sampled data. The data preparation consists of re-sampling the initial timetables using linear interpolation, defining the new sample as the statistic mode of raw data sampling times, which is 𝑡𝑠=20 s. It is worth noting that the sampling time of 20 s repeats in more than 96% of the raw data. Therefore, the plant’s automation system has this sampling time which the authors maintained for homogenization. After the data import and re-sample procedure, it passes through a concatenation, resulting in eight regular timetables: from June 9th to 12th, 23rd to 26th, 29th to 30th, August 11st to 12th, 14th to 17th, September 2nd to 3rd, 22nd to 25th, and October 14th to 16th, totaling 25 days of operational data. Applied Energy 339 (2023) 120944 4 D.O. Machado et al. Fig. 1. Fresnel solar collector installed at ETSI absorption solar plant [34]. Fig. 2. Fresnel solar collector model installed at ETSI absorption solar plant [37]. Next, a moving mean filter with a 6 min window was applied to remove outliers from data sets. The window was chosen by trial and error for proper outliers detection and substitution considering a series of visual inspections of the raw data and the outliers’ clean data. In addition, a one-by-one evaluation substituted negative pressures, flows, and inconsistent variables by the minimum or maximum absorption solar collector operation range considering the manufacturer’s datasheet [37]. The variable-by-variable inspection also led to the discovery of empty columns that were further deleted. Furthermore, pressure and flow variables were smoothed using Gaussian interpolation with a 15 sampling window. The sampling window was chosen by trial and error, considering a good compromise between lag and dynamic information loss (long window) and noise (short window). As a result, 14 variables compose the final data set with 108 416 samples. It is worth saying that the previous literature on the modeling of the FSC on the ETSI does not discuss the experimental data quality, the preprocessing, or the data preparation procedure for simulation reproducibility nor considers the intermittent operation of the plant Applied Energy 339 (2023) 120944 5 D.O. Machado et al. with the same amount of data points. These differences highlight an important improvement in FSC modeling in this work. The following section presents the data selection based on correlation maps. 2.2. Correlation analysis of FSC field I/O variables The correlation coefficient measures the association between variables, the most widely used is the linear correlation coefficient. Pearson’s correlation coefficient (𝜌)for pairs of variables (𝐱,𝐲)with 𝑛 samples 𝐱= [𝑥1,1,…, 𝑥𝑛,1]and 𝐲= [𝑦1,1,…, 𝑦𝑛,1]is given by 𝜌(𝑥, 𝑦) = 1 𝑛− 1 𝑛 ∑ 𝑖=1 (𝑥𝑖,1−𝜇𝑥 𝜎𝑥)(𝑦𝑖,2−𝜇𝑦 𝜎𝑦),(1a) 𝜇𝑥=1 𝑛 𝑛 ∑ 𝑖=1 𝑥1,𝑖,(1b) 𝜎𝑥=√ √ √ √1 𝑛− 1 𝑛 ∑ 𝑖=1 (𝑥1,𝑖 −𝜇𝑥)2,(1c) where 𝜇𝑥and 𝜎𝑥are the mean and standard deviation of the samples of 𝐱, respectively, and 𝜇𝑦and 𝜎𝑦are the mean and standard deviation of 𝐲. The values that 𝜌can take are between [−1,1], where 𝜌= −1 represents a complete negative correlation, 𝜌= 1 represents a complete positive correlation, and a value of 𝜌= 0 indicates that the variables (𝐱,𝐲)are uncorrelated. The correlation coefficient matrix (𝐑)of the random variables 𝑀is 𝐑∈ℜ𝑀×𝑀for each combination of pairwise variables. Only for 𝑥and 𝑦, result in the following matrix 𝐑∈ℜ2×2: 𝐑=(𝜌(𝑥, 𝑥)𝜌(𝑥, 𝑦) 𝜌(𝑦, 𝑥)𝜌(𝑦, 𝑦))=(1𝜌(𝑥, 𝑦) 𝜌(𝑦, 𝑥) 1 ).(2) The next step is sorting the input variables in a degree of relationship. The sort evaluates the degree of correlation for each input regarding the desired output to analyze which variables impact the model output concerning a threshold. In the case of FSC, the desired output is the outlet temperature 𝑇 𝑓2.Fig. 3 shows the correlation coefficient matrix of the actual Fresnel solar field data. Where 𝑓is the mirrors’ focus, varying between 0 for full defocus to 1 for complete focus; 𝑞is the volumetric flow, 𝑃is the hydraulic pressure, 𝑇 𝑓1and 𝑇 𝑓2are the inlet and outlet temperatures, respectively; 𝑤𝑠 is the wind speed, 𝐼 is the solar irradiance, 𝐻is the humidity, 𝑇 𝑎𝑚𝑏 is the environment temperature. 𝐸𝑓𝑓𝑂𝑝𝑡 is the optical efficiency accordingly to Brandão et al. [43], and 𝑇 𝑓1𝑑𝑒𝑙𝑎𝑦 is the delayed inlet temperature calculated accordingly to Normey-Rico et al. [44]. 2.3. FSC operation Fig. 4 depicts the main operation variables from the prepared data considering two operation days: from June 24, 06:00, to June 26, 00:00, 2009. It is worth saying that the validation section uses the same exact actual data. By inspecting Fig. 4.a, one can see two days of irradiance profile. The irradiance profile (continuous yellow line) of June, 24, is smooth, with the sunrise at about 06:00, a peak at 15:00, and the sunset at 22:00. The operation stages follow the irradiance profile, where the start-up occurs when there is enough irradiance, 𝐼≥200 (W/m2), at almost 09:00 turning on the pump, see Fig. 4.b. Note that the flow goes from 0 to 13 m3/h together with the mirror’s focus from 0 to 100%. Next, the heating phase occurs with the FSC increasing its temperature from 09:00 to 12:00. Then, the absorption chiller consumes the FSC thermal power. It generates a sudden temperature drop followed by a temperature increase return between 12:00 and 14:00. The temperature drop results from the absorption chiller’s lower internal temperature injection in the FSC hydraulic loop. After a recovery period, the temperature increases again due to the gas boiler operation. After the absorption chiller start-up, aided by the backup gas boiler from 12:00 to 14:00, the boiler is shut-down. The plant operates roughly at the nominal point with solar irradiance only, producing chilled water for the ESTI air conditioning system. Note that the temperatures are highly oscillatory from 14:00 to 20:00 due to the absorption chiller on–off controller operation. The two-position valve effect is evident in the hydraulic loop flow (blue dotted line) in Fig. 4.b, where 𝐹varies between two levels during operation. This bi-stable event happens because when the controller feeds the chiller High Temperature Generator, a heat exchanger with long tubes generates a pressure drop reflecting the flow changes in Fig. 4.b. For further information on the whole plant control logic, refer to [28]. When the irradiance is too low, that is I ≤200 (W/m2), the flow and mirror’s focus go to zero. That is why the plant shuts down at 20:00. Then, from June 24, 22:00, to 25, 09:00, the FSC temperatures decrease overnight, see Fig. 4.a. After these heat losses, the operation is re-started on June 25, at 09:00, where the sequence of the same events as the previous day happens. It is worth noting that the irradiance profile of June 25, unlike the previous day, is very oscillatory. Fig. 4 depicts the operation data of the ETSI absorption plant. The process is complex, highly non-linear, and intermittent, with a broad operational range of the variables, such as temperatures, flows, and solar irradiance. Thus, developing the FSC digital twin is not trivial, and the respective models must cope with and describe such a wide range of operations and phenomena to reflect the dynamic behavior of the FSC. After the FSC digital twin creation, it allows for optimizing the start-up, operation, and shut-down, considering changing meteorological and plant conditions while offering accurate predictions for model-based predictive control techniques. The following sections present the two model structures used in this work. 3. Phenomenological - PDE modeling Partial differential Eqs. (3) and (4) describe mathematically the Fresnel’s temperature distributed in time and space [36]. 𝜌𝑚𝑐𝑚𝐴𝑚 𝜕𝑇𝑚 𝜕𝑡 (𝑡, 𝑥) =  𝑄𝑠𝑢𝑛(𝑡) −  𝑄𝑎(𝑡, 𝑥) −  𝑄𝑓(𝑡, 𝑥),(3) 𝜌𝑓𝑐𝑓𝐴𝑓 𝜕𝑇𝑓 𝜕𝑡 (𝑡, 𝑥) + 𝜌𝑓𝑐𝑓𝑞(𝑡)𝜕𝑇𝑓 𝜕𝑥 (𝑡, 𝑥) =  𝑄𝑓(𝑡, 𝑥),(4) where the sub-indexes 𝑚,𝑓,𝑠𝑢𝑛, and 𝑎refer to the metal absorber tube, fluid, the sun, and the environment, respectively. The variable 𝜌is the specific mass (kg/m3), 𝑐is the specific heat (J/(kg ◦C)) and 𝑞is the volumetric flow (𝑚3∕ℎ). Variable 𝑇, is the temperature (◦C), 𝑡is the time (s), 𝑥∈ [0, 𝐿], with 𝐿 > 0, is the space (m),  𝑄𝑠𝑢𝑛 is the sun heat rate (W) that flows from the sun to the solar collector,  𝑄𝑎represents the thermal losses (W) to the ambient, and  𝑄𝑓is the heat rate (W) that relates to the mass flow. Eq. (5) gives the boundary condition of Eq. (3), 𝑇(𝑡, 0) = 𝑇 𝑓1(𝑡),(5) where 𝑇 𝑓1is the Fresnel inlet temperature. Finally, the initial conditions of the system are given by Eq. (6) 𝑇𝑚(0, 𝑥) = 𝑇0 𝑚(𝑥), 𝑇𝑓(0, 𝑥) = 𝑇0 𝑓(𝑥),(6) where 𝑇0 𝑚and 𝑇0 𝑓are functions that satisfy the steady state condition of Eqs. (3) and (4). Papers in literature use the system of Eqs. (3)–(6) for both parabolic through and Fresnel collectors. The difference among the different solar collector’s equations is mainly in calculating the optical efficiency 𝜂𝑜𝑝𝑡 once each type of solar collector has a specific primary and secondary mirrors scheme [45]. The optical efficiency is in the inlet solar heat term,  𝑄𝑠𝑢𝑛, which Eq. (7) describes,  𝑄𝑠𝑢𝑛 =𝜂𝑇𝜂𝑜𝑝𝑡𝐴𝑡𝐼(𝑡),(7) Applied Energy 339 (2023) 120944 6 D.O. Machado et al. Fig. 3. Correlation coefficients matrix of the prepared ETSI FSC field data. Fig. 4. FSC consecutive days operation actual data. a. The FSC outlet temperature 𝑇 𝑓2is the controlled variable, while the inlet temperature 𝑇 𝑓1and Irradiance 𝐼are the prominent disturbances. b. The FSC installed at ETSI has two manipulated variables, the typical flow 𝑞, and the mirror’s focus 𝑓. where 𝜂𝑇and 𝜂𝑜𝑝𝑡 are the thermal and optical efficiency, respectively, 𝐴𝑡(m2) is the total area of collector’s mirrors, and 𝐼is the solar irradiance per mirror length in W/m2. A Fresnel efficiency has both variable and constant parameters in time. For example, having a given collector’s orientation and solar time enables one to calculate the solar beam’s incidence angles and the respective reflection cosine losses. However, the reflexivity and other characteristics of the mirrors vary with plant aging, dirt accumulation, and water condensations in the mirrors. Thus, this work considers that the variable 𝜂𝑜𝑝𝑡 contains the deterministic and constant parameters, which are calculated accordingly to Brandão et al. [43], considering the geometrical relation between the mirrors and absorber. With 𝜂𝑜𝑝𝑡 is possible to calculate the equivalent irradiance that arrives in the absorber tube 𝐼𝑒𝑞 =𝐼(𝑡)𝜂𝑜𝑝𝑡. The time-varying and, a priori, unknown efficiency-related parameters are embedded in one unique parameter called thermal efficiency 𝜂𝑇, which this work estimates using an identification technique. This work applies Euler’s finite differences discrete approximation, according to Fig. 5 schematic, to solve Eqs. (3) and (4), that are continuous in time and space. The spatial derivative in a given time instant is given by 𝜕𝑇𝑓 𝜕𝑥 (𝑡, 𝑥) ≈ 𝑇𝑓(𝑡, 𝑛) − 𝑇𝑓(𝑡, 𝑛 − 1) 𝛥𝑥 , Applied Energy 339 (2023) 120944 7 D.O. Machado et al. Fig. 5. Absorber tube discretization considering S =16. where 𝛥𝑥 =𝐿∕𝑆is the length of the spatial discretization, 𝑆is the number of points, and 𝑛∈ {1,…, 𝑆}is the respective given volume. Note that 𝑇𝑡,0=𝑇 𝑓1(𝑡), as Eq. (5) presents. Constant time intervals, 𝑀, composes the whole time horizon 𝑡∈ [0, 𝑡𝑓], where 𝑡𝑓is the final time, such that 𝜏𝑘∈ [0, 𝑡𝑓],𝑘= {1,…, 𝑀}, are points in the discretization mesh concerning the time, with 0 = 𝜏1<⋯< 𝜏𝑀−1 < 𝜏𝑀=𝑡𝑓. This discretization is considered equidistant for simplicity. Thus, 𝛥𝑡 =𝑡𝑓 𝑀− 1 , 𝜏𝑘= (𝑘− 1)𝛥𝑡, 𝑘 ∈ {1,…, 𝑀}. and the time derivatives approximations are 𝜕𝑇𝑚 𝜕𝑡 (𝑡, 𝑥) ≈ 𝑇𝑚(𝑘+ 1, 𝑥) − 𝑇𝑚(𝑘, 𝑥) 𝛥𝑡 , 𝜕𝑇𝑓 𝜕𝑡 (𝑡, 𝑥) ≈ 𝑇𝑓(𝑘+ 1, 𝑥) − 𝑇𝑓(𝑘, 𝑥) 𝛥𝑡 . The derivative approximations in time and space above transform Eqs. (3)–(4) in algebraic discrete Eqs. (8) and (9) 𝑇𝑚(𝑘+ 1, 𝑛) = 𝑇𝑚(𝑘, 𝑛) +𝛥𝑡 ( 𝑄𝑠𝑜𝑙(𝑘) 𝜌𝑚𝑐𝑚𝐴𝑚 −  𝑄𝑎(𝑘, 𝑛) 𝜌𝑚𝑐𝑚𝐴𝑚 −  𝑄𝑓(𝑘, 𝑛) 𝜌𝑚𝑐𝑚𝐴𝑚),(8) 𝑇𝑓(𝑘+ 1, 𝑛) = 𝑇𝑓(𝑘, 𝑛) +𝛥𝑡 (𝑞(𝑘) 𝐴𝑓 𝑇𝑓(𝑘, 𝑛) − 𝑇𝑓(𝑘, 𝑛 − 1) 𝛥𝑥 +  𝑄𝑓(𝑘, 𝑛) 𝜌𝑓𝑐𝑓𝐴𝑓),(9) where  𝑄𝑎=𝐷𝑚𝜋[𝑎(𝑇𝑚(𝑘, 𝑛) − 𝑇𝑎(𝑘))3+𝑏(𝑇𝑚(𝑘, 𝑛) − 𝑇𝑎(𝑘))],(10)  𝑄𝑓=𝐷𝑓𝜋𝐻𝑡(𝑇𝑚(𝑘, 𝑛) − 𝑇𝑓(𝑘, 𝑛)),(11) 𝐷𝑚(𝑚)is the equivalent diameter of the tube walls, 𝐷𝑓(𝑚)the internal tube diameter, 𝑎and 𝑏are the coefficients of heat losses of the absorber polynomial function, and 𝐻𝑡(W∕(m ◦C)) the coefficient of heat transfer between the tube walls and the fluid. It is worth noting that this paper proposes a third-order heat loss polynomial function to describe the heat losses coefficient, Eq. (10), instead of fourth and second-order polynomials that are typically used [32,33]. The third-order polynomial is employed because it becomes possible to change the sign of  𝑄𝑎of Eq. (10), and consequently, the heat losses term of Eq. (8). That is, by using a third-order polynomial, the tube model can also represent the case where the heat enters the tube instead of exiting. The case where the tube receives heat occurs if the plant does not operate for days. Therefore, the tube’s temperature tends to be ambient temperature. In this case, it is possible that 𝑇𝑚< 𝑇𝑎. Thus, the model can reproduce the plant’s dynamic behavior in realtime, becoming a tool for what-if analysis of the start-up and shut-down of the plant independently of the state of the process after day and night of operation. To integrate Eqs. (8) and (9), it is necessary to iterate 𝑇𝑚(𝑘+ 1, 𝑛) and 𝑇𝑓(𝑘+1, 𝑛)from the initial and boundary conditions given in 𝑇(0, 𝑛) and 𝑇(𝑘, 0), respectively, from n =1 to n =S, and from k =1 until k = M. This work employs the ode45 package of MATLAB [46] to integrate the phenomenological model. 3.1. Identification The previous section has presented the FSC’s model considering phenomenological concepts and, therefore, the main dynamics that occur in the plant. The model has thermodynamical characteristics such as specific mass, specific heat, and heat transfer coefficients between the heat transfer fluid, the metallic tube, and the environment. Despite knowing these parameters for pure substances and materials references, the values can vary sensibly due to the plant aging, corrosion, mirror soiling, and internal pipe walls crusting. Thus, these parameters are time-varying; thus, it is necessary to identify their values to plug in the model for representing the process in a given time regarding the data used for identification. One advantage of Eqs. (8) and (9) is the possibility to estimate the parameters’ values by comparing the model output with experimental data. The estimation employs a quadratic non-linear minimization algorithm through Eq. (12) to adjust the proposed model parameters. The cost function is the sum of the normalized quadratic error between the Fresnel model predicted outlet temperature, 𝑇𝑓(𝑘, 64), and the measured outlet temperature, 𝑇∗ 𝑓(𝑘). The model’s parameters are the decision variables of the optimization problem, with the initial points and maximum bounds set considering experimental and manufacturer data. For 𝑛∈ {1,…, 𝑆}, Eq. (12) defines the optimization problem: min 𝜂𝑇,𝜌𝑓,𝑐𝑓,𝑎,𝑏,𝐻𝑡 𝑀 ∑ 𝑘=1 (𝑇𝑓(𝑘, 64) − 𝑇∗ 𝑓(𝑘))2 𝑇∗ 𝑓(𝑘)2 subject to, Eq. (8), Eq. (9), 𝑇𝑚(0, 𝑗) = 𝑇0 𝑚(𝑗), 𝑇𝑓(0, 𝑗) = 𝑇0 𝑓(𝑗), 𝑇𝑓(𝑘, 0) = 𝑇 𝑓1(𝑘), 0≤𝜂𝑇≤1, 800 ≤𝜌𝑓≤1000, 4200 ≤𝑐𝑓≤4500, 0≤𝑎≤∞, 0≤𝑏≤∞, 353 ≤𝐻𝑡 ≤2500, (12) where 𝑀= 207 605 is the number of measurements, the lower and upper bounds were chosen based on the materials’ physical parameters and properties tables. Each Eq. (12) iteration integrates Eqs. (8) and (9) along the tube, from 1to 𝑆, and the time, from 1to 𝑀. The code execution continues until the stopping criteria condition. Such conditions are the maximum number of iterations and the objective function derivative convergence to a minimum constant value. The fmincon [47] Matlab’s algorithm solves the model parameter identification problem. The decision variables are the thermal efficiency 𝜂𝑇, the specific mass of the fluid 𝜌𝑓, the polynomial coefficients 𝑎 and 𝑏of thermal losses regarding Eq. (10). Lastly, 𝐻𝑡 is the heat transfer coefficient between the metal and water. The identification procedure Applied Energy 339 (2023) 120944 8 D.O. Machado et al. Fig. 6. Adaptive Neuro-fuzzy Inference System (ANFIS) architecture [35]. considers that all parameters are equal along the tube length. Besides, the simulation considers steps of 5 s to avoid numerical instability; therefore, it is necessary to interpolate the input data accordingly. Section 6.1 discusses the identification results of the FSC parameters. 4. Neuro-fuzzy modeling The modeling of a system to faithfully represent its behavior has a certain level of complexity due to the random dynamics of the unknown nature of the process, making it a challenge to describe its behavior utilizing mathematical equations. Often an accurate model will be represented by several mathematical equations. However, a model supported by mathematical tools (e.g., differential equations) is only sometimes adequate to deal with uncertain systems. In addition, developing control and optimization strategies using nonlinear models to obtain the system prediction implies a high computational burden. It is the principal hindrance if the problem must be solved in a given constrained time, considering the sampling time, so a fast model is indispensable. Fuzzy models have proven to be an effective technique for modeling and controlling nonlinear systems, successfully expressing the original nonlinear model as a set of local linear models interpolated by a membership function containing the nonlinearities of the original model. In addition, fuzzy inference systems (FIS) describe the behavior of a process based on rules with linguistic labels from human language. Thus, FIS represents the qualitative aspects of human knowledge and reasoning processes, avoiding precise quantitative analyses. The ANFIS architecture [35], also called the Neuro-Fuzzy (NF) system, is an artificial intelligence (AI) technique. One ANFIS has five layers, as depicted in Fig. 6, where the nodes of the first (I) and fourth layers (IV) have adaptive (square blocks) parameters, and the remaining layers have fixed parameters (circle blocks). The first layer contains the fuzzification interface, which transforms a crisp input into linguistic labels with a certain degree of membership, forming the fuzzy sets [48,49] characterized by the membership functions (MF). The parameters are adaptive and are called antecedent parameters. The second layer outputs the product of the incoming signals from each fixed node and represents the firing strength of each node rule. The third layer normalizes the output of each node, calculated as the ratio of the firing strength of the node’s rule to the sum of all the firing strengths of each node rule. In the fourth layer is the defuzzification interface, and each node function provides the weighting of a first-order polynomial crisp function, whose parameters are called consequent parameters. Lastly, the fifth layer contains a single node aggregating all rule outputs. For a detailed description of the ANFIS architecture network, refer to [35,50]. 4.1. Training of ANFIS The resulting FIS, after training, contains a set of rules of Takagi– Sugeno [51] type as the following: 𝙸𝙵 𝑥1𝚒𝚜 𝐹1𝑗𝚊𝚗𝚍 𝑥2𝚒𝚜 𝐹2𝑗𝚊𝚗𝚍 𝑥𝑖𝚒𝚜 𝐹𝑖𝑗 , 𝚃𝙷𝙴𝙽 ∶𝑓𝑗(𝑥) = 𝑔0𝑗+𝑔1𝑗𝑥1+⋯+𝑔𝑖𝑗 𝑥𝑖. Each rule has an antecedent and consequent parameter. The fuzzy sets 𝐹𝑖𝑗 of each crisp entry 𝑥𝑖consist of Gaussian membership functions of the type: 𝜇𝐹𝑖𝑗 (𝑥𝑖) = 1 1 + [(𝑥𝑖−𝑐𝑖𝑗 𝑎𝑖𝑗 )2]𝑏𝑖𝑗 (13) where {𝑎𝑖𝑗 , 𝑏𝑖𝑗 , 𝑐𝑖𝑗 }are the antecedents parameters that define the mean, height, and width of the Gaussian used to vary the MFs,1and the terms 𝑔𝑖𝑗 ∈ℜof each first-order polynomial function are the consequent parameters. The learning process of the ANFIS neuronal architecture adapts both parameters. The output of each rule 𝑓𝑗is a linear combination of input variables added to a constant term. The final output of the fuzzy inference system is the weighted average of each rule’s output. The learning process of the ANFIS network architecture uses the training and checking sets to capture better the system’s dynamics, which allows for an acceptable model that predicts its behavior. The ANFIS uses as inputs the variables {𝐼, 𝐸𝑓𝑓𝑂𝑝𝑡, 𝑓, 𝑇𝑎𝑚𝑏, 𝑞, 𝑇 𝑓1, 𝑤𝑠, 𝐻, 𝑃 , 𝑇 𝑓1𝑑𝑒𝑙𝑎𝑦}and as output 𝑇 𝑓2, thus constructing a mapping of the input– output variables that represent the behavior of the solar field. Initially, the subtractive clustering (SC) method [52] is used to estimate the number and initialization centers of the Gaussian MF of the fuzzy rules. In addition, learning employs a hybrid learning method. The method runs a back-propagation algorithm [53] to obtain the parameters defining the MF of each fuzzy set (antecedents parameters). Next, the learning also executes a least-squares to estimate the terms of the first-order polynomial function (consequent parameters) of the output of each rule at each epoch. An epoch, or sweep, is one forward and backward parameter update. The checking procedure evaluates the error between the ANFIS output and the actual output of the checking data set (a new data set not used in training). The checking runs after each epoch during the training and aims to evaluate if the ANFIS training results in generalized learning. If the ANFIS output has low errors with unknown inputs, then it is said that the ANFIS model had general learning. Typically, the checking considers the root mean squared error (RMSE) given by Eq. (14) 𝑅𝑀𝑆𝐸 =√∑𝑁 𝑖=1(𝑥𝑖,𝑗 −𝑥𝑖,𝑗 )2 𝑁,(14) 1The value that the function 𝜇𝐹𝑖𝑗 takes for a given 𝑥𝑖is known as the degree of membership of 𝑥𝑖for the fuzzy set 𝐹𝑖𝑗 . Applied Energy 339 (2023) 120944 9 D.O. Machado et al. Fig. 7. a - SFC outlet heat-transfer fluid temperature PDE model. b - SFC outlet heat-transfer fluid temperature ANFIS model. Table 2 ANFIS parameters. Description ANFIS MF type: Gaussian Number MFs: 4 Number rules: 4 Influence range 0.8 Epoch number: 1500 Table 3 Simulation planning. Schematic Fig. 7.a Fig. 7.b Fig. 7.b Structure PDE(dt =5) ANFIS(dt =5) ANFIS(dt =20) noinputs variables 6 10 10 nooutputs variables 𝑇𝑓and 𝑇𝑚𝑇 𝑓2𝑇 𝑓2 Twinning Par. identification Learning Learning Method LSaLS and GDbLS and GD Data-set 𝐺𝑇 𝑤𝑖𝑛𝑛𝑖𝑛𝑔 𝑃 𝐷𝐸 𝐺𝑇 𝑤𝑖𝑛𝑛𝑖𝑛𝑔 𝑁𝐹 𝐺𝑇 𝑤𝑖𝑛𝑛𝑖𝑛𝑔 𝑁𝐹 Data-set size 4 days 19 days 19 days Data-set samples 14 188 311 666 78 032 Sampling time (s) 5 5 20 Stopping criteria Eq. (12) nRMSEcnRMSE Simulation ODE45 1◦polynomial 1◦polynomial Data-set 𝐺𝑉 𝑎𝑙1,𝐺𝑉 𝑎𝑙2𝐺𝑉 𝑎𝑙1,𝐺𝑉 𝑎𝑙2𝐺𝑉 𝑎𝑙1,𝐺𝑉 𝑎𝑙2 Data-set size 6 days of 2009 6 days of 2009 6 days of 2009 Data-set samples 95 282 95 282 23 822 Integration step (dt) 5 5 20 Validation Table 5 Table 5 Table 5 aLS - Least-squares. bGD - Gradient Descent. cnRMSE - Normalized Root Mean Squared Error. where 𝑥𝑖,𝑗 is a given actual variable j with 𝑁samples, and 𝑥𝑖,𝑗 is the output of the predicted variable. This work considers normalized outputs (𝑧𝑖,𝑗 )for training and checking. Therefore, the normalized RMSE, given by Eq. (15), is used. 𝑛𝑅𝑀𝑆𝐸 =√∑𝑁 𝑖=1(𝑧𝑖,𝑗 −𝑧𝑖,𝑗 ) 𝑁.(15) Table 2 presents the parameters of the ANFIS that capture the fluid outlet temperature behavior. 5. Simulations planning This section describes the FSC PDE, and NF models structure with inputs/outputs, their twinning/updating characteristics, and validation specifications. The input selection considers the prepared data sorted accordingly to the outlet temperature 𝑇 𝑓2correlation coefficients, described in the Correlation Matrix in Fig. 3. Thus, the PDE and NF models have the structure, inputs, and characteristics summarized in 7.a and .b, respectively. Table 3 compares Fig. 7.a PDE model and Fig. 7.b model with two different integrations steps of 5 s and 20 s. The PDE model has six inputs, while the NF model has ten inputs. This difference happens because, on the one hand, the PDE model does not describe meteorological variables such as wind speed 𝑤𝑠, humidity 𝐻, and hydraulic circuit pressure 𝑃. On the other hand, it is easy to add these input variables to the NF model. The wind speed and air humidity affect the heating process of the Heat Transfer Fluid, changing the convective heat coefficient between the absorber tube and the environment, refer to Fig. 3. In addition, the PDE model intrinsically describes the FSC inlet temperature dead time; thus, the delayed inlet temperature 𝑇 𝑓1𝑑𝑒𝑙𝑎𝑦 is employed only in the NF model. Besides, the model outputs are different. While the PDE model has two outputs vectors 𝑇𝑓(𝑘)and 𝑇𝑚(𝑘) that describe the temperature gradient along the absorber tube length, the NF model has a scalar outlet temperature output 𝑇 𝑓2. It is worth noting that the PDE outputs are arrays with 16 values, where the last value represents the outlet temperature 𝑇 𝑓2 = 𝑇𝑓(16, 𝑘); After defining the models’ structure, the next step is twinning them to adjust their parameters to follow the plant’s actual data. The PDE model is updated using least squares (LS), while the NF model uses LS and Gradient Descent (GD). The PDE model employs 𝐺𝑇 𝑤𝑖𝑛𝑛𝑖𝑛𝑔 𝑃 𝐷𝐸 as twinning data set with four days of measurement, totaling 14 188 samples. While the NF model twinning data set is 𝐺𝑇 𝑤𝑖𝑛𝑛𝑖𝑛𝑔 𝑁𝐹 , which uses 19 days of data, totaling 311666 and 78032 samples, for sampling times of 5 s and 20 s, respectively. Refer to the middle section of Table 3 for details. The twinning rate of a given digital twin is the updating rate of a model. Thus, knowing the twinning time of a given model is essential to defining the twinning rate of the virtual entity once the twinning rate must be lower than the models’ twinning time. The twinning time is the processing time the model’s twinning takes to converge the model error to a given minimum stopping criteria tolerance. The twining data reduction for the PDE model is because the computational burden of the PDE parameter identification is high. In preliminary tests, the identification leads to more than one week of processing using the whole available twinning data set, which results in impractical twinning rates of the FSC. In addition, the sampling time was reduced to five seconds to avoid numerical instability of the PDE integration, increasing the total number of samples—see Table 3. The re-sampling of 𝐺𝑇 𝑤𝑖𝑛𝑛𝑖𝑛𝑔 𝑃 𝐷𝐸 , further PDE model updating and simulation would lead to unfair execution times compared to the NF model. Note that the latter model was initially updated and simulated, considering the prepared data with a sampling time of 20 s. Therefore, to compare the PDE and NF models updating performance, the NF model twinning is run with sampling times of 5 s, which is the same re-sampled data of the PDE model. Each column of Table 3 describes each situation, with the first column representing the PDE model with integration steps of dt =5 s in four days of actual data, while the second and third columns have the NF model with dt =5 s and the NF model with dt =20 s, specifications, respectively. Both NF models use the same 19 days of actual data operation. The performance indexes for twinning are the total twinning time 𝑡𝑡𝑤, which is the total time that it takes to update the model, the twining time per sample, 𝑡𝑡𝑤∕𝑠𝑎𝑚𝑝𝑙𝑒,