Optimal Water Management in Agro-Industrial Districts: An Energy Hub's Case Study in the Southeast of Spain

In this work, the optimal management of the water grid belonging to a pilot agro-industrial district, based on greenhouse cultivation, is analyzed. Different water supply plants are considered in the district, some of them using renewable energies as power sources, i.e., a solar thermal desalination plant and a nanofiltration facility powered up by a photovoltaic field. Moreover, the trade with the water public utility network is also taken into account. As demanding agents, a greenhouse and an office building are contemplated. Due to the different water necessities, demand profiles, and the heterogeneous nature of the different plants considered as supplier agents, the management of the whole plant is not trivial. In this way, an algorithm based on the energy hubs approach, which takes into account economic terms and the optimal use of the available resources in its formulation, is proposed for the pilot district with a cropping area of 616 m2. Simulation results are provided in order to evidence the benefits of the proposed technique in two cases: Case 1 considers the flexible operation of the desalination plant, whereas in Case 2 the working conditions are forced to equal the plant’s maximum capacity (Case 2). A flexible operation results in a weekly improvement of 4.68% in profit, an optimized use of the desalination plant, and a reduction of the consumption of water from the public grid by 58.1%.


Introduction
Climate change has a great impact on the water supply sector of many regions worldwide, such as in southern Europe, where the hydrological stress is expected to shortly increase due to this phenomenon and some areas are already facing serious water problems indeed [1]. One of these is the province of Almería (southeast of Spain), which is identified as one of the driest regions in the continent, but has one of the largest agricultural production systems. What is more, the main driving force of the economy in this province is agriculture, with around 30,000 ha of greenhouse crop production [2]. The development of this industry has been tied, for a long time, to the depletion of freshwater reservoirs (despite the efficient management of this resource that has been performed in this sector [3]). This fact made the case for the inclusion of alternative water sources such as desalination plants, thus enhancing the availability of fresh water for the sustainability of the Almerian agricultural system [4]. As a consequence, the current panorama in Almería can be visualized as a distributed water network composed of (i) producers, based on traditional (wells and water public utility network) and non-traditional sources (desalination facilities and other non-conventional systems), and (ii) consumers, such as industries related to agriculture and greenhouses.
As in any other kind of distribution network or multi-agent system, this requires integral and optimal management [5][6][7] and, until now, different approaches dealing with this issue have been introduced and formulated in the literature. For example, in the studies of Ocampo Martinez et al. [8] and Pascual et al. [9], Model Predictive Control (MPC) formulations were proposed for the efficient management of the urban water cycle of Barcelona (Spain) in terms of operating costs. A similar approach was employed by Lopez Farias et al. [10] but, in this case, to improve the forecasting of the control method. In the work presented in Ref. [11], a distributed MPC approach, which was aimed at decreasing the required resolution time and computational cost, was put forward for the same problem. Another interesting strategy was presented in Ref. [12], where an MPC controller was in charge of optimizing the energy-water nexus in urban water networks. Furthermore, a scheduling method was formulated by Zhang et al. [13], in which the distributed water network was modeled and divided into three levels: water supply source, water station, and water user. Then, the developed model was used to implement effective scheduling strategies for the guidance of regional water distribution systems.
Although all the aforementioned studies presented satisfactory results, the proposed techniques are focused on optimizing the transport water network by minimizing the operational costs. However, in most cases, not only the water needs to be optimally managed but also energy trading. The presence of desalination plants and other non conventional water sources introduces new agents of intermittent nature in the problem, as they are normally powered by renewable energy to improve their efficiency [14,15] and to reduce water costs [16]. In addition, adequate storage systems [17] and the use of the public utility network as a backup are required for their continuous operation. These elements, together with the related industries and greenhouses' necessities, make the whole system to constitute an agro-industrial district [18] in which the generation, storage, and distribution of heterogeneous resources must be optimally performed for its efficient proper exploitation.
Over the last few years, some studies dealing with the optimal management of multiagents systems that include non-conventional water sources, such as desalination plants, have been published. Most of these analyses are based on the typical approach for an energy hub (EH), which relies on a simplified modeling methodology that represents the interactions given inside manifold systems, attending to its input-output configuration [19]. Gharffarpur et al. [20] proposed a scheduling method to manage an isolated district, including a Reverse Osmosis (RO) plant, renewable energy generation, and electricity, heat, and water demands. However, their formulation did not consider the connection to the public utility network and the objective of the scheduling method was not aimed at improving operating costs, but it was devoted to maintaining some levels of resilience. Authors in Ref. [21] introduced an alternative modeling methodology to integrate energy and water systems that include multiple energy sources. This was then used to pose a Mixed Integer Nonlinear optimization Problem (MINLP) and to perform the optimal management of the different systems comprising a shale-gas production plant, but no connection to the public network was taken into account either. In the work presented in Ref. [22], a similar procedure was also proposed to manage an isolated micro-grid located in an island, considering the optimal dispatch of water and energy. In this case, the planning method was tasked with minimizing the environmental pollutants as well as the operation and investment costs. Once more, the approach is only valid for off-grid regions (without public utility network connections) such as islands. A more complete approach was recently published by the authors in Ref. [23], in which instance the multiagent system also included issues related to energy and water, but both conventional (i.e., a well) and non-conventional water sources (i.e., a desalination plant) were considered. The management method was also based on other paradigms of energy hubs' and the planning problem included a multi-objective optimization algorithm that minimized, at the same time, the energy and the fresh water extracted from underground reservoirs. Although in this case a connection to the public utility network was included, the sale of surplus resources was not considered in the problem, which might prove to be a determining factor for the correct economical exploitation of agro-industrial districts. Accordingly, based on the above review, the main gaps identified in the literature that are the basis of this manuscript are listed as follows: • At present, the management strategies presented for distributed water networks are mainly focused on the optimal performance of the transportation system, without contemplating the sources of water, and without considering any other resource apart from water in the problem. • Regarding the strategies based on multi-agent systems, they are mainly focused on isolated plants. In this way, trading with the public utility network is not usually considered, neither are the sale of surplus resources nor backup systems to carry out continuous operations, which can be fundamental in agro-industrial districts to maintain the desired level of quality and productivity.
To cover the above issues, this paper extends the case study addressed in our previous work Ref. [18] by including all the plants related to water issues (as well as the exchange of resources among them) contemplated in the CHROMAE research project ("Control and Optimal Management of Heterogeneous Resources in Agroindustrial Production Districts Integrating Renewable Energies", www2.ual.es/chromae (accessed on 25 January 2021)). The resulting system is in fact an agro-industrial district in which several water sources and resources of different nature must be managed. In addition, the connection to the water public utility network and the sale of surpluses is considered. Unlike the approach presented in Ref. [18], which was focused on the development of an enabling platform for the management of the distributed facilities, and presented a basic MPC strategy for solving the problem, in this work we use the approach presented in Ref. [24] but adding new capabilities to the energy hub model. For example, variable dependent loads are considered by introducing a term in the model that reflects the variable production of the facilities. This is especially important to the operation of desalination facilities in these kind of environments as it allows the management strategy to adapt their production to the operating conditions at each sampling time (water needs, availability of resources, etc.), and therefore to perform an optimal dispatch. By using real historical meteorological data and a tool developed for such problems (ODEHubs), a simulation is presented to exemplify this fact and then compared to a manual operation, evidencing the benefits achieved by the proposed strategy.
The rest of the paper is arranged as follows. Section 2 presents the case study, the energy hub modeling approach, and the particularization of this model for the analyzed system. Section 3 presents the simulation results and the comparison with the manual strategy. Finally, Section 4 discusses about the main outcomes of the work.

Case Study Description
The case study adopted in this work comprises all the plants related to water issues included in the CHROMAE research project. This project aims to contribute, from a control point of view, to the optimal dispatch of heterogeneous resources in a way that ensures equitable access, sustainability, and efficiency. Specifically, in this work, we focus on a pilot agro-industrial cluster composed of an office building, a greenhouse, a solar thermal desalination plant, a nanofiltration plant, and a photovoltaic facility. This case study encompasses all the problems that can be found in real agro-industrial districts in terms of integral resources management, as we include all the resources shared among the different agents, renewable energy generation, and the connection to the public utility network. On the one hand, we consider the water and electricity demands of the greenhouse and the office building, and the thermal energy demand and generation (through a boiler) of the greenhouse. On the other, we take into account the generation of electricity, through the photovoltaic plant, and water, through the nanofiltration and the solar desalination plants, which in turn requires electricity, as well as thermal energy supplied by a solar thermal field for the latter one. The reader is referred to Figure 1 which visually reflects the exchange of resources among the different agents. In addition, all the plants are based on real facilities located in Almería (Southeast of Spain): • Solar Membrane Distillation (SMD) plant: The SMD plant is based on the pilot facility located at Plataforma Solar de Almería (PSA, www.psa.es (accessed on 25 January 2021)). This plant was totally described in the work in Ref. [25] and it is composed of: a solar thermal field, with a nominal capacity of 7 kW th at 90 • C, a slope angle of 7 • and facing south (azimuthal angle of 0 • ); a storage tank, with a volume of 1.5 m 3 ; and a Membrane Distillation (MD) module. Although there are several MD modules available at PSA, in this work we use the Aquastill module which has a total membrane surface area of 24 m 2 providing a maximum distillate production of around 30 L/h under the best operating conditions [26]. • Nanofiltration plant: This facility is also based on a pilot plant located at PSA. It consists of three membranes which can work in parallel or series; that is, the main stream of wastewater can be divided, in order to be filtered in each membrane, or entirely pass from one to the next. In either case, this results in a permeate stream (with low content of micropollutants that are able to break through the membranes) and a concentrated stream (rich in them), part of which is fed back to the tank containing the wastewater (recirculation stream). The system is configured to control the flow and electrical conductivity of each stream (permeate, concentrate, and recirculated), as well as the transmembrane pressure, thanks to the pressure, temperature, and flow meters and to the valves installed at certain points in the circuit. As the operation range of the plant depends on the amount of pollutants in the water, indirectly determined by the electrical conductivity of the stream, readers are referred to the most recent publication on this plant for further information, where the static behavior was tested using a feed stream (H2O + NaCl) with a conductivity of 2300 µS/cm [27]. • Greenhouse: The greenhouse environment used as reference in this work is also located at the Experimental Station of the Cajamar Foundation (40 km away the PSA). This environment consists of a parral-type greenhouse with a total surface area of 877 m 2 (37.8 m × 23.2 m), 616 m 2 of these being useful for cultivation. The heating system of this facility is a GP 95 propane heater and a Missouri 150,000 multi-fuel boiler fed with bio-mass. Further information about this system can be found elsewhere [28]. • Office building: The office building included in the agro-industrial district is based on the Centro de Investigación de la Energía Solar (CIESOL) building (www.ciesol.es (accessed on 25 January 2021)). This building is placed at the Almería's University campus, 20 km far from PSA, and it has a total surface area of 1071.91 m 2 divided into two floors. Although it can provide both electricity and thermal energy through its photovoltaic and solar thermal fields (see Ref. [29] for further information), in this work, we assume that it behaves as an agent demanding electricity and water, trying to imitate the role of a real office building in an agro-industrial district. • The photovoltaic parking of the University of Almeria (36.83 • N, 2.40 • W) consists of ten arrays of 483 Conergy PA 240P modules (twenty-one in series and twenty-three in parallel) connected to a Fronius Agilo 100 inverter (4830 modules and ten inverters in total) for production purposes, as well as another three arrays, for experimental testing, connected to their respective Fronius IG Plus 55v3 inverters, and distributed as follows: twenty-four Conergy PA 240P modules (twelve in series and two in parallel) in the first group, twenty-four Conergy Power Plus 240M modules (twelve in series and two in parallel) in the second and seventy-two First Solar FS-380 modules (eight in series and nine in parallel) in the third. All the modules have a slope angle of 7 • and are facing south (azimuthal angle of 21 • east) [30].

Modeling and Formulation of the Problem
Managing the system presented in Figure 1 involves deciding how and when each device should operate, which, for certain small-scale or pilot plants, can be manually done, based on operators' expertise, but it might better rely on making these decision automatically, by means of optimization techniques. As is the case with this study, the ODEHubs toolbox allows these kinds of problems to be analyzed in MATLAB ® and Simulink ® . The general model implemented in it, as well as the library itself, were presented by some of the authors of this paper formerly [24]. Before going deeper into the configuration parameters and the particular equations that would represent the above-mentioned district, it is necessary to identify a series of elements that give rise to what the authors call conceptual model of the EH ( Figure 2 and Table 1): the resources demanded by the plant that need to be produced within the EH (denoted with O), the ones that can be traded as purchases (denoted with I) and sales (denoted with M), and the devices available for their conversion (denoted with D) and storage (denoted with S).  Table 1. Note that although I 2 represents a flow of radiant power, it is actually calculated as the sum of flows corresponding to the solar facilities, hence different values of irradiance, solar angles, and geometries can be considered for each field. Water sold through the public utility network m 3 /h Notice that, in Figure 2, the arrows indicate the flows of resources within the EH which can have as start and end points either devices or convergence and diverge flows. This layout has been broadly used in other occasions since it eases the transition from the so-called conceptual model to the corresponding diagram in Simulink ® by representing as a remainder two important features, as exemplified in Ref. [24]: on the one hand, ODEHubs's users should care about the nodes where flows are merged and splitted, which need to be included in Simulink ® as specific blocks (CN, which stands for "convergence node" and DN, which stands for "divergence node", respectively, in Figure 3); on the other hand, the missing arrows indicate that no conversion exists for these flows, but they just represent the dependence of loads on the devices that constitute their end point (in ODEHubs they actually bypass these devices). The following subsections provide more information on the mathematical model auto-generated by ODEHubs for this district and the configuration parameters employed in the simulations.

Particularization of the General Conversion and Storage Model
As the full generic model that defines of the optimization problem employed by ODEHubs to get the dispatch of resources was presented in Ref. [24], only the brief definition of the vectors and matrices, which have been adapted to the EH presented in in Figure 2, will be provided here. The relationships between them can be consulted in Equations (1)-(17) of the cited paper, and the dispatch of resources is given by the following optimization problem: where k constitutes any discrete time instant, c is the vector containing the price of each input, s is the vector containing the price of sold resources, δ I is the binary diagonal matrix of input flow activation, δ M is the N o × N o binary diagonal matrix of market sales flow activation, and H is the length (in samples) of the control horizon. Several constraints, which define mass and energy balances as well as production and transportation limits, are imposed on the vectors whose elements constitute decision variables: P, I, M, Q ch , Q dis S, Let T and diag be the transpose operator and the operator that returns a square diagonal matrix with the elements of a vector, respectively, and consider any variable written in bold font as the distinctive form of matrices and vectors. Notice also that the energy hub in Figure 2 counts with six inputs, seven outputs, 15 paths between inputs and outputs (see Table 2), and five conversion devices with a single input and output.  or stored by a certain device and the purchase or sales of them. The parameters for each block that define their minimum and maximum limits, as well as their efficiencies, will be presented in the following section.
In addition, the coupling matrices that define the efficiency of the charge and discharge of resources and the degradation of the stored resources are defined as diagonal matrices, where η represents each of the devices' efficiencies in these processes: C ch = diag([η ch,1 · · · η ch,7 ]), C dis = diag([η dis,1 · · · η dis,7 ]), C s = diag([η s,1 · · · η s,7 ]). In contrast, when it comes to the rest of the coupling matrices, which depend on the internal structure of the EH, their definition is not straightforward and for the case of Figure 2 this results in Equations (2)-(5): In order to model all the constraints of the optimization problem, the remaining diagonal matrices are required to define the state of operation of devices and distribution networks (on/off), where δ represents each of the binary variables: . To define δ O , one must however consider the existence of dependent loads, as further explained in Ref. [24]. Figure 2 represents O 2 , O 4 , O 5 , and O 6 as fixed dependent loads, since they appear multiplied by the binary variables of the device to which these loads are attached (δ D,2 , δ D,3 , δ D,3 , and δ D,4 , respectively). Actually, the amount of heat required by the distillation units is proportional to the flow of water, hence O 5 constitute a variable dependent load defined as follows: O 5 = κη D,3 P 13 . As P 13 is already dependent on η D,3 through the constraints on D i and D o , that binary variable does not need to be included in δ O , which yields: δ O = diag([1 δ D,2 1 δ D,3 1 δ D,4 1]). Note that C has already been arranged to consider this issue by introducing the term −κη D,3 in Equation (2).
Lastly, c = [c 1 · · · c 6 ] is the vector containing the price of each input and s = [s 1 · · · s 7 ] is the vector containing the price of sold resources (in terms of energy, mass, or volume) and owing to the use of the public utility networks for both importing and exporting electricity and water, Equations (6) and (7) are considered as constraints in ODEHubs: δ I,6 + δ M,7 ≤ 1. Figure 3 shows the Simulink ® model of the case study built from ODEHubs's blocks, which will be used to simulate the dispatch problem and to obtain the optimization results according to the procedure presented in [24] and the elements found in Figure 2. Tables 3-5 summarize the parameters that define the case study and that must be substituted in the equations presented in the above section. In those tables, the infinite symbol is used to denote unconstrained variables, i.e., those variables related to inexhaustible resources. Note also that regarding the elements not included in Figure 2 (M 2 , M 3 , M 4 , M 5 , M 6 , Q ch,2 , Q dis,2 , S 4 ,Q ch,4 , Q dis,4 , S 4 , Q ch,6 , Q dis,6 , S 6 ), the corresponding limits are fixed to zero by the ODEHubs library as shown in Tables 3 and 5.

Characterization of Devices and Demands via ODEHubs
In the case of conversion devices (Tables 3 and 4), their upper limits are imposed only for their inputs or outputs, obtaining a proper representation of each device according to their physical constraints or limitations on production. The biomass boiler's parameters correspond to the ones found in its data-sheet, whereas the behavior of the other devices has been analyzed in former studies. For example, the isotropic sky model allows ODEHubs to calculate the incident radiation on sloped surfaces for the respective parameters (latitude, longitude, slope angle, area, etc.) of the parking and the solar collectors [30], which has been encapsulated in internal functions to get the maximum input flow of both devices (D max i,1 and D max i,2 ). Similarly, the equivalent circuit for photovoltaic cells provides the field's performance (η D,1 ) via another user-defined function [30]. The desalination and nanofiltration plants have been an object of analysis to determine their features as well [26,27] and only the storage systems are hypothetical elements, as justified in Ref. [31].     In respect to the resources considered as outputs (demands) of the energy hub, the electricity consumption of the pumping systems of the solar collectors (0.1943 kW), the MD plant (0.1720 kW), and the nanofiltration plant (0.2502 kW) were taken into account in the problem separately, as device-dependent outputs (O 2 , O 4 , and O 6 , respectively), allowing us to relate them to the on/off state of each pump. In addition, the value of κ in O 5 = κη D,3 P 13 is fixed at κ = 121.25 kW/m 3 , considering the experimental data and model provided in Ref. [26]. The rest of the demands (O 1 , O 3 , and O 7 ) were calculated by processing the real historical data measured by the sensor systems of each of the plants considered in the agro-industrial district (see Figure 4). At this point, note that both CIESOL's consumption profile and the parking's size have been ten times reduced for a more realistic test.
6.5  Regarding the parameters that define the cost function of the optimization problem solved by ODEHubs, water (c 6 = 0.9024 EUR/m 3 ), biomass (c 3 = 0.255 EUR/kg), and electricity current prices are obtained from local suppliers, with their corresponding valueadded tax included. In the case of the electricity (c 1 ), the cost during winter time depends on the period of the day as follows: 0.0892 EUR/kWh from 0:00 h to 8:00 h, 0.2044 EUR/kWh from 18:00 h to 22:00 h, and 0.1127 EUR/kWh for the rest of the day. The rest of resources are assumed to be freely available, hence neither solar radiation nor untreated water are taxed (c 2 = c 4 = c 5 = 0).
On the other hand, all the elements of vector s except s 1 and s 7 are equal to zero because there is no possibility to sell the resources of outputs 2 to 6 (as represented in Figure 2). Water is considered to be sold at a much lower price (s 7 = 0.2 EUR/m 3 ) because many Almerian farmers, who usually make use of private wells, have been claiming to the local councils that the price of desalinated water should be subsidized to achieve competitiveness by keeping it around 0.2-0.3 EUR/m 3 . In addition, the actual hourly electricity prices from intraday markets, together with the generation fee of 0.5 EUR/MWh and the 7 % tax on the value of electricity generation, are considered to obtain s 1 (see Figure 5), as further detailed in Ref. [19]. They show an overall tendency to fall during midnight or early morning and correspond to period of time selected for this study, between the 16-22 March 2014, which was chosen because of the combination of several given conditions: relatively high consumption of water, simultaneous heat demand at the greenhouse, and presence of both overcast and clear sky. Finally, the meteorological variables attached to the solar production models are presented in Figure 6. Over the week of data, one can observe that on 16-19 and on 21 March, it was sunny with some passing clouds on 18 and 21, whereas the two other days (on 20 and 22) were nearly overcast. This phenomena is noticeable in the rising and falling peaks of both the direct and global irradiances, and also in the diffuse radiation, which tend to increase due to cloudiness. The temperatures demonstrate the typical warm climate of the region, with a minimum of 8.

Results and Discussion
The analysis presented in the following subsections consists of a comparison between the operation of the EH, over a week in March 2014, under two different conditions placed on the desalination plant: Case 1 considers the operation point, in terms of the desalinized flow of water (D o,3 ), as a decision variable; whereas in Case 2 that variable is constrained to force that the plant work at its maximum capacity whenever it is activated (see Table 3), which is representative of what a human operator would do. ODEHubs was used to solve the optimization problem presented in Ref. [30] and above adapted for the district described in Section 2.1. The toolbox was configured to use a scheduling mode strategy (fixed horizon, with deterministic and measurable disturbances) and the solver intlinprog [32] and it was executed on an Intel ® Core TM i7-6700K 4GHz CPU, taking around 96 s to simulate each of the two cases considered. Figures 7 and 8 correspond to Case 1 and Figures 9 and 10 correspond to Case 2, in which the results are arranged in the same way, that is, using a think solid line to represent each of the hourly demand profiles (outputs of the energy hub that form O) together with a colored stacked bar graph showing the sources (elements of I) that meet those demands. The solid thin lines depicting the evolution of the sold resources represent in fact the accumulated values of M 1 + O 1 and M 7 + O 7 , respectively, hence the actual values of M 1 and M 7 can be obtained by subtraction. All these elements are expressed in terms of power or flow according to the scale of the left vertical axis, whereas the right one is only employed for the dashed lines that symbolize the state of charge of each of the storage systems (elements of S), in terms of energy, volume, or mass. Note that both the charge (Q c ) and discharge flows (Q d ) are deliberately missing in the said figures, since they can be deducted from either the slope of the dashed lines or the difference between the colored stacked bar graph and the stacked solid lines. The results have been split into non-dependent outputs (Figures 7 and 9), i.e., those ones related to loads or demands that cannot be controlled and would need to be met ad hoc; and dependent outputs (Figures 8  and 10), which hinge on the on/off state of some of the devices.
Regarding the similarities in both cases, the broad strategy to schedule the dispatch of resources over the week consists in making use of the solar energy to yield either energy or water from dawn to dusk at a lower cost. This is close to zero, since radiation was considered to be freely available (c 2 = 0), but not exactly zero because sometimes using electricity or biomass compensates the cost of acquiring water from the public network, which is the case for the nanofiltration plant on 16-20 March (Figures 8 and 10). In Case 1 (but not in Case 2), this also happens to the desalination plant on 18 March night ( Figure 9) and on 17's dawn and dusk, when the modules were allowed to operate at partial load. As a result that the electricity price tends to be higher around midday, except for 20 March (see Figure 5), and the electricity demanded by the facilities (CIESOL and the greenhouse) remains quite stable at about 6 kW, most of the electricity produced by the photovoltaic field is either directly sold or stored to be sold straightaway after noon, when the price is still profitable (see Figures 7 and 9). That is also the reason why no electricity is kept stored at night, given the abundance of solar radiation the following days and the lower price at that time. O 3 and O 7 are managed analogously, but they show more variability in the demand. Water (O 7 ) is required especially during the working time and when the evapotranspiration rises (as outdoor temperature and irradiance do), which happens to be higher on 17-20 March (hence, storage is required during the previous night); whereas the thermal needs (O 3 ), which are nil on 16 March, exhibit a quite irregular pattern.
On the contrary, the differences found between Case 1 and Case 2 are mainly justified from Figures 8 and 10. Since the desalination plant is no longer able to work at a partial load in Case 2, the amount of excess water that would be produced in comparison to Case 1 does not compensate the additional energy required for it. In other words, as the levels of irradiance and stored energy close to the dawn and the dusk are not enough to cover the needs of the plant operating at its maximum capacity, owing to physical constraints, it is preferable to use water from either the public network or the nanofiltration plant. This fact is noticeable in the increase of water coming from these sources in Figure 8 and in the amount of electric energy supplied to O 6 on 20-22 March (Figure 10). For the same reason, Figure 10 shows a minor usage of S 5 , in comparison to the same variable in Figure 9. Finally, why a first glance might look like a baffling performance on 18 March, when S 3 is no longer used to store heat, is explained because in Case 1 the biomass boiler is used to feed the desalination plant (see Figure 9), whereas in Case 2 it is used to feed the facilities, but the amount of consumed biomass is altogether the same, as summarized in Table 6.    In addition to the above discussion and figures, the accumulated amounts of resources involved in each case are summarized in Table 6 for vectors I, O, and M, as well as the total cost of purchasing or selling resources according to vectors c and s. Note that Case 1 and Case 2 differ in O 2 , O 4 , O 5 , and O 6 , which depend on the decision variables that minimize the operation cost of the plant. As shown in Figures 7 and 8, the constraint placed on the desalination plant in Case 2 produces a shift of supply sources and more water is acquired via nanofiltration or the public network. Thus, the flexible operation of Case 1 results in a higher profit of 4.68%, a decrease in the water consumption from the public grid of 1.8 m 3 (58.1%), and an increase in the desalination plant's production of 20.3%, in comparison to Case 2. Note that, in addition to the economic benefits, if this is extrapolated to the annual operation of the entire Almerian region, it would have a significant impact on the environment owing to the amount of water that would not be extracted from the sweetwater reservoirs.

Conclusions
This paper compares two cases of water management where the operation of a desalination plant is constrained to work similarly to a manual operation mode, in contrast to its flexible use, in which the amount of distilled flow is adapted to the consumption needs. The simulations performed on a realistic test-bed plant, which can be defined as an energy hub, show that those constraints actually make the system to economically underperform. It also proves ODEHubs to be suitable tool for resource scheduling problems, since the results for both cases are coherent with respect to the intuitively expected behavior of the controller regarding the economic objective and they constitute a paradigm of device-dependent variable loads, whose theoretical framework was presented in Ref. [24] but not exemplified by a real-world study as the one above-presented. Although the scheduling strategy was used in this paper to demonstrate the differences in the flexible operation-deterministic scenario where all the disturbances are beforehand known, without a feedback loop-ODEHubs allows the users to simulate MPC-based strategies with receding horizon, which would take into account the presence of uncertainty in more realistic scenarios.
Some related studies have shown a similar performance when operating distillation facilities in terms of economic costs, as in Ref. [33], in which the cost per unit of demanded water was 0.44 EUR/m 3 . If compared with the results obtained, considering the same electricity and water prices as in Ref. [33] together with the material and energetic needs in Table 1, cost results in 0.12 EUR/m 3 for Case 1 and 0.14 EUR/m 3 for Case 2. This difference is due to the fact that in this study the distillation unit is enforced to yield water at the optimal operation point. Note that, in both studies, the thermal cost is neglected because it is assumed to be met by the field of solar collectors.
On the other hand, these analyses and methods are applicable to the wider Almerian ecosystem, in order to manage sets of greenhouse and water producers distributed over certain area. The modifications to be performed on the ODEHubs's model and on the code itself have been already posed and are under implementation within a framework based on game theory in which each of the prosumers constitute different players with, at times, opposing interests.

Acknowledgments:
The authors would like to thank Plataforma Solar de Almería, Experimental Station of Cajamar foundation, and CIESOL Center for facilitating real data of their facilities.

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

Abbreviations
The following abbreviations are used in this manuscript: