A Flexible Tool for Modeling and Optimal Dispatch of Resources in Agri-Energy Hubs

: The dispatch of energy and resources in agricultural systems often involves the deﬁnition and resolution of optimization problems. This paper presents a novel tool composed of a set of MATLAB R (cid:13) and Simulink R (cid:13) ﬁles that has been developed to ease such tasks. In contrast to other alternatives, it allows the consideration of multiple kinds of resources in the problem and the relationships between the inputs and outputs of the system; its parametrization can be deﬁned graphically in Simulink R (cid:13) without requiring third party software, and the entire package is freely available on Github. The package can generate the constraints in MATLAB R (cid:13) code and can get the optimal dispatch schedule for the deterministic mixed-integer linear problem that represents the deﬁned system. Its main functions and blocks as well as a case study based on a traditional Mediterranean greenhouse and a photovoltaic parking lot located in Almeria (Spain) are included to demonstrate its use and clarify how the problem is formulated. The simulation performed validates the tool as being useful for decision-making (schedule irrigation and CO 2 enrichment, as well as managing storage systems) in these and similar environments. Future implementations are intended to incorporate the interconnection of agents with opposed interests and robust optimization strategies for uncertain scenarios.


Introduction
Today's realization that the energy sector must move towards a more efficient, affordable, reliable, safe, and sustainable system is encouraging research on the generation, transport, consumption, and storage processes [1]. Many authors postulate that decentralized and distributed approaches may exploit the local renewable energy sources [2] and combine different energy carriers, which results in a synergy that contributes to the sector's development [3]. In this regard, the energy hub (EH) concept [4] arose in 2006 for designating production systems with the latter characteristic (although most works include RES as well, but this is not something inherent to the EH's definition). Geidl and Andersson defined it as "the interface between producers, consumers, and transport infrastructure, which, from the point of view of a system, provides the relationship between the inputs, outputs, conversion, and storage of different energy carriers" [5]. Over the last years, the number of papers related to EHs has been constantly growing, as supported by some of the reviews that give an overview on this topic [6,7], and the term has settled down as well as others, such as microgrid (MG), virtual power plant (VPP), or multi-energy system (MES), which tend to appear in similar energy dispatch problems [8].
A quick look at some of the latest works related to distributed generation shows that most of them rely on solving optimization problems either for design or management purposes. For instance, particle swarm optimization (PSO) was used to obtain the optimal scheduling of a hybrid (thermal and electric) system [9] and to size ten grid-connected hybrid blocks aimed at reducing the supply deficit in Sierra Leone [10]. Model predictive control (MPC) is also associated with control problems in this kind of system: It has been applied to the optimal operation of hybrid power systems in the sugar cane industry [11], and more recently to networked microgrids [12,13] and desalination processes [14]. This is germane to the field of precision agriculture and agro-industrial production districts, where the growing concern about overexploitation of exhaustible resources and the need to maintain the modern economy and quality of life have raised new challenges for scientists. In terms of energy waste, the economic impact of the underlying ground during the production period has been analyzed in facilities such as broiler houses [15]. Another concern, particularly in greenhouses, is that energy consumption seems to have settled down as a significant factor that hinders their development [16], and many proposals have started to tackle dispatch problems from a wider perspective, that is, integrating all the resources required for crop development (water, energy, CO 2 , etc.); just to cite a few examples: optimally using combined heating and power facilities [17], determining the operation scheduling of a greenhouse with CO 2 enrichment and artificial lighting [18], or managing water production and its use [19,20].
The ODEHubs (Optimal Dispatch for Energy Hubs) tool developed during the CHROMAE project [21] by the authors of this paper-from the Automatic Control, Robotics, and Mechatronics Research Group (TEP-197) [22]-aims to contribute to the optimal management of these resources in a way that ensures equitable access, efficiency, and sustainability in agro-industrial districts. To put ODEHubs's features into context, the following paragraphs outline some related software packages as well as their strengths and weaknesses.
At present, many of the multifarious tools of interest that cope with these issues have been encompassed in reviews on district-scale [23] and urban [24] energy systems, where the spotlight is mainly on consolidated and specific engineering software for design and analysis. Other reviews focused on broader contexts [25,26] comprise less mature applications that rely on general-purpose programming environments or languages-such as CPLEX [27], GAMS [28], AMPL [29], AIMMS [30], MATLAB R [31], R [32], or Python [33]-and even, sometimes, on the self-learning capacity of the users. To the authors' knowledge, MATLAB R is extensively employed in this field of research, both as a standalone version and combined with other optimization-assisted modeling tools. It happens to be one the more common programming languages as well, together with Python and GAMS [25], but despite its potential, to date, very few specific tools for dealing with the most common problems in EHs have been developed and made freely available.
1. Ehub Modeling Tool [34]. This consists of a set of MATLAB R scripts for creating input case study data, which are sent and then executed in the optimization package AIMMS. It also incorporates R code for visualizing the results of the simulations. A subsequent version was released in Python, which includes a graph-based modeling tool. Software available at: https://github.com/hues-platform/ehub-modeling-tool. 2. EHCM Toolbox [35]. An object-oriented programming tool that was conceived to deal with the control necessities of a building and entirely in MATLAB R . Consequently, it might lack flexibility in certain aspects, which is counterbalanced by a sophisticated model for batteries and the possibility to integrate detailed building dynamics from another toolbox. Software available at: https://control.ee.ethz.ch/software/BRCM-Toolbox1.html. 3. MATPOWER Optimal Scheduling Tool [36] (MOST). This is an extension of the MATPOWER package, a free and open-source set of MATLAB R files for solving steady-state electric power scheduling problems. MOST can be used to solve from simple deterministic unconstrained problems to stochastic, security-constrained, combined unit-commitment, and multi-period optimal power flow problems. The current implementation is limited to DC power flow modeling of the network, so it is a bit limited when it comes to multi-energy systems (MESs) and EHs. Software available at: https://github.com/MATPOWER/most. 4. IMAKUS [37]. This is the name given to a deterministic model of the German power system aimed at making investment decisions on storage and conventional generation assets. Although it was implemented in MATLAB R code using a generic formulation, its authors do not clarify its possible use on other regions or with other resources than electricity. No support website. 5. LUSYM [38]. A tool quite similar to MOST in the sense that it is focused on the electricity sector and allows the introduction of nearly the same constraints in the problem. The main difference resides in the use of GAMS together with MATLAB R , which are interfaced so that GAMS is used for solving the optimization problem and MATLAB R for processing the input and output data. No support website. 6. EUPowerDispatch [39]. This is another tool for the optimal dispatch of the power system, but particularized at the European level, and it uses GAMS and MATLAB R in the same way that LUSYM does. More information available at: https://ses.jrc.ec.europa.eu/eupowerdispatch-model. 7. Software Library for Multi-Energy Systems (an official name has not been given yet) [40]. This is probably the most flexible and complete of the alternatives, is based on object-oriented programming, and is derived from the generalized modeling framework for MESs proposed by Long. It is intended to be used for MPC applications rather than scheduling problems. The same version was built both in Python and MATLAB R , although none of them is freely available. No support website.
All of these tools have contributed in one way or another to establishment of the basis from which developers can keep providing new algorithms and procedures to the field of energy and resource management. Nonetheless, there is still room for improvement in certain aspects that the researchers might demand: Less experienced users might struggle with coding, and there exists no graph-based modeling tool for the MATLAB R environment that allows one to formulate and solve resource dispatch problems; four out of the seven packages are focused on purely electric systems (3,4,5,6), only three are freely available online and properly documented (1,2,3), and some of them require third-party software that is subject to a license (1,5,6).
Thus, given the current alternatives, ODEHubs constitutes a tool with significant improvements. First, multiple carriers can be considered in the dispatch problem thanks to the general model for EHs and MESs [41] that has been incorporated, and, what is more important, the system can be defined via Simulink R blocks, by analogy with its input-output graph representation, or with MATLAB R code. Second, it allows the users to choose between two optimization modes and different solvers: on the one hand, MPC-based (receding horizon strategy) or scheduling (fixed horizon); on the other, the internal MATLAB R function intlinprog, so no third-party software installation is required, or any compatible external solver through YALMIP [42] (for instance, CPLEX [27]). Finally, the set of MATLAB R and Simulink R files that compose this tool are freely available on Github (https://github.com/ual-arm/odehubs), together with the 'ready-to-use' data-set of the case study presented in Section 4. Other contributions of this work are the extension of the prior general model [41] in order to add device-dependent variable loads and providing an example of the use of the tool in a case study with real-world data of a traditional Mediterranean greenhouse.
The remainder of this paper is organized as follows: Section 2 presents the model coded to solve the optimization problem; Section 3 contains a description of the main modules and their configuration parameters; Section 4 shows an example case where the ODEHubs toolbox is employed to obtain the operation schedule for the greenhouse and its facilities; Section 5 analyzes the current drawbacks of ODEHubs and future enhancements.

Conversion and Storage Model
The preliminary version of the ODEHubs model was formulated and published in a previous study [41], to which readers are referred for clarification, and has been kept without meaningful changes. In broad terms, the models of most energy hubs represent energy and mass balances between a set of input resources that can be transformed into another set of output resources. The basic relationship proposed by [5] is based on a coupling matrix C that represents the conversion process between multiple energy carriers in steady state (L − CP = 0), where P are the input powers or sources and L the output powers or loads. The coefficients of this matrix are given by the product of the efficiencies of all the conversion units or devices that intervene in the conversion and the so-called dispatch factors that need to be introduced whenever one energy carrier splits up to several flows inside the EH. Additional constraints may consider power limits, the presence of storage systems, and necessary energy balances among these dispatch factors. This has led to formulations that have progressively introduced elements to represent conversion and storage processes more accurately, but share a similar basis, as [41] discussed in the reference study [41]. In contrast to other proposals in the literature, the so-called "path vector" was defined to consider a decision variable for each possible path between inputs and outputs and to avoid the non-linear relationships of [5]'s model when two dispatch factors are multiplied in the coupling matrix [5]. The input and internal flows are obtained from the elements of the path vector and suitably defined coupling matrices. Some outputs or loads can be considered dependent on the fact that a certain device is turned on or off, which was taken into account by arranging [5]'s equation L − CP = 0. All the equations from [41] that encompass these and other circumstances are presented below, although, for the sake of simplicity, some of them will appear here in matrix notation. In addition, the term "power" will be superseded by "flow" or just omitted, since not only energy carriers are considered, but also material resources.

Prior Formulation
Consider a general energy hub with N i input flows, N o output flows, N p paths between these inputs and outputs, N d conversion devices with a total amount of N di input flows and N do output flows, and discrete time models using a uniform sample time T = t(k + 1) − t(k), where k constitutes any discrete time instant.
The conversion and storage processes need to satisfy balance conditions at each time step and between consecutive time steps, respectively, as stated in Equations (1) and (2): where P is the N p × 1 vector of flows between inputs and outputs, M is the N o × 1 vector of market sales flows, O is the N o × 1 vector of output flows, Q ch is the N o × 1 vector of charge flows, Q dis is the Additionally, the relationship between the internal P vector flows and both the input flows to the EH itself must be established, as defined in Equation (3), as well as the input-output flows in each device, as given in Equations (4) and (5): where I is the N i × 1 vector of input flows, D i is the N di × 1 vector of devices' flows (referred to their inputs), D o is the N do × 1 vector of devices' flows (referred to their outputs), C i is the N i × N p input coupling matrix, C di is the N di × N p device coupling matrix (referred to their inputs), and C do is the N do × N p device coupling matrix (referred to their outputs). Note that in Ref. [41], only one equation was defined, instead of two, to limit the flows through the devices of the EH. Except vector O, which a priori does not contain decision variables (see next subsection for the opposite case), and vector P, which is indirectly limited by the production capacity of the conversion devices from the following equations, the remaining vectors need to be constrained in order to ensure that the resource conversion and storage devices do not exceed their defined capacities. That is the raison d'être of Equations (6)- (12): where the superscripts max and min refer to the lower or upper limit of each vector (for instance, S max determines the maximum amount that each storage system may store), and the rest of the δ symbols define diagonal binary variables that indicate if the flow is allowed or not. In particular, δ M is the Additional constraints are required in order to take into account certain processes that cannot occur simultaneously. Three situations can be distinguished: preventing the charge and discharge of the same storage systems, by means of Equation (13): where 1 is the identity matrix (note that I already denotes a variable of the model); avoiding the use of the same infrastructure for acquiring and selling resources, as expressed in Equation (14): where i and o respectively denote the elements of δ M and δ I with common infrastructure; and impeding that some devices operate simultaneously, as stated in Equation (15): where d (1) ...d (n) refer to the elements of δ D i or δ D o corresponding to non-simultaneous converters.
Owing to the definition of the path vector, the variables corresponding to the outputs of a device that produces different resources at the same time (e.g., a combined heat and power system or a boiler with carbon capture) must equal the input amount of that device. Equation (16) satisfies that balance: where p (1) ...p (n) ...p (m) correspond to simultaneous converter outputs. These situations will be exemplified in Section 4, although readers are encouraged to consult Ref. [41] for a better understanding. Finally, the cost function of the optimization problem can be written as in Equation (17): where c(k) is the 1 × N i vector containing the price of each input, s(k) is the 1 × N o vector containing the price of sold resources (in terms of energy, mass, or volume), and H is the length (in samples) of the control horizon. Note that, depending on the units of the sample time T and the vectors, additional terms might need to be introduced in Equations (2) and (17) to convert between time units.

Device-Dependent Fixed and Variable Loads
Note that, in the arrangement proposed in Ref. [41], to consider loads coming from the use of a particular device consisted in adding any necessary element to vector O and using δ O in Equation (1). However, this solution is only suitable in the case of fixed demands, since the model does not contemplate the possibility of adding variable demands proportionally to the flow of the device. One straightforward way to iron out this situation consists in redefining the elements of O so that the demand equals the proportional flow of the device or, in what turns out to be equivalent, redefining matrices C and δ O . In the following paragraphs, both situations are exemplified by using a hypothetical EH in order to elucidate how the previous formulation needs to be adapted, which is relevant for programming.
Let Figure 1 represent a simplified example of an EH where an electric pump is in charge of keeping the water pressure of an irrigation facility. In order to ease the understanding of this example, free water discharge and no electricity storage or market sales flows are considered. The path vector only has two elements, that is, the number of possible paths in the EH: I 1 → D 1 → O 1 (P 1 ) and I 2 → O 2 (P 2 ), so, in this concrete case, C i becomes the identity matrix and P = I. Thus, the balance equations of the conversion processes and their limits, i.e., Equations (1), (3)-(5), (8) and (9), adopt the form of Equations (18)-(23): where the non-zero elements of the coupling matrix C and C do would be equal to one because no transformation of resources exits in the EH, and δ O contains the binary variable related to the activation of the pump (δ D,1 ). Now, if the electric pump is assumed to have a constant consumption, O 2 would be equal to that value instead of constituting a decision variable in the optimization problem, which is the case explained in Ref. [41]. However, if a consumption proportional to the pumped water flow is considered instead, with a constant of proportionality κ, Equation (18) needs to be reformulated as follows. First, it must be defined whether the consumption is proportional to the inlet or the outlet flow, which would lead to two different expressions: Only the second expression will be used from now on, as it will allow generalizations to the first one. In any case, as both P 1 and δ D,1 are decision variables, their product would turn the optimization problem into a non-linear one when substituting O 2 into Equation (18). To avoid this, δ D,1 can be removed from δ O (making that element equal to one), since P 1 is already constrained by the on/off state of the pump through Equations (22) and (23), which leads to transformation of Equation (18) into Equation (24): or, equivalently, arranging matrix C into Equation (25): Just to recap all of the above and as a rule of thumb, in case of device-dependent variable demands, the following modifications must be performed in the model presented previously [41]-removing from δ O , by making equal to one, all the elements of its diagonal related to those demands; setting to zero the analogous rows in O; and including in the columns of C related to that device the suitable conversion/transformation coefficients multiplied by the negative constant of proportionality between the demand and the flow through the device. Part of the ODEHubs code was generated to perform these modifications in the model.

ODEHubs Toolbox
ODEHubs consists of a set of functions and files that allow one to obtain the optimal dispatch or scheduling of any system that can be modeled according to the above-presented equations. Users can manually define the characteristics of this one (number of elements, flow relationship, etc.) and configure certain parameters depending on their needs. In order to serve as a short guide for anyone interested in the tool, the main functions and settings are described in the following lines, considering the general flowchart of the software presented in Figure 2), although specific documentation will be made available on the website of the tool (https://github.com/ual-arm/odehubs). The first step required to use ODEHubs, and from which the model is generated and adapted to be introduced in the optimization solver, consists in defining the problem, which involves manipulating three or four files, depending on whether the problem is defined graphically or by programming (see their names next to yellow trapezia in Figure 2). Except for Simulation_main.m, which is placed in the root directory, the rest of them can be found in the definition folder. Note also that although the manual procedures, which are described below, have been included sequentially in Figure 2, they are actually interchangeable.

•
Simulation_main.m is the base of the flowchart and is where most functions are called. At the beginning of the file, a set of parameters, whose use is clarified below, are declared, and users might need to modify these lines prior to simulation. A struct variable called EH is created during its execution and used to store variables and communicate them to other functions.

•
DataEH.mat is a file in the binary data container format that the MATLAB R program uses, with a single timetable variable called dataEH stored in it. This type of variable allows the grouping of column-oriented data in a table where each row represents a date and time. Thus, the columns consist of all the time-dependent variables required to simulate the system, such as the demanded outputs or weather conditions. Users are responsible for including the data that they need in this format, with the freedom to choose the name of each column/variable. Columns may be empty if no data are required to be loaded, but the timetable needs to be generated according to the simulation start time and horizon (see these parameters' definition below).

•
The folder user_models contains any function that users might need to define time-variable parameters in Equations (1)- (17). They need to be arranged following the syntax parameter = function (data, date, samples, tm), so that parameter is a vector of values obtained from data (which is an automatically filtered version of dataEH) starting at the time declared in date in each tm period of time during the horizon samples.
• EH_definition_code.m is the function employed to set the properties of the energy hubs, and it includes the parameters related to Equations (1)- (17) via MATLAB R code. The file provided by the authors of this work for the case study presented in Section 4 can be used as a template for other systems just by adapting its content.

•
ODEHubs_components.slx is the Simulink R library, which, with the support of the user_models folder, incorporates the predefined models of several storage and conversion devices developed during the CHROMAE project [21] (Figure 3). Its blocks, which are described in the next subsection, are employed to build Simulink R models that represent any EH either by copying and pasting from the library or by dragging from the library browser if ODEHubs_components.slx is added. As in any other library, the file is locked when first opened, but users, and even the developers in future releases of ODEHubs  Figure 3. ODEHubs's library components defined in Simulink R , which can be used to define customized systems either by including their elements in Simulink's browser or by copying them from the library file. The name that appears under each block can be changed when the block is included in a Simulink R model, as well as the labels within Input and Output blocks.
The remainder functions are not intended to interact with the users directly, but they are part of the automatic process of calculating the dynamic of the EH and storing and displaying the final results. They are distributed in different folders according to their purpose, but all of them appear in Simulation_main.m and are responsible for calling other specific functions. EH_definition_Simulink.m, features.m, and preset.m arrange the input data from the users to be processed by other functions; load_data.m loads the content of DataEH.mat and filters it by time according to values specified by the users; MPC_predictions.m and MPC.m provide an updated dispatch of the EH's resources at each sample time specified for the control actions; EH_system.m simulates the response of the EH at each sample time specified for the system's dynamic; and process_results.m aggregates the flows during the simulation horizon and displays the scheduled dispatch. At this point, note that including MPC_predictions.m and EH_system.m is a mere formal aspect because they currently make use of the same models and variables; thus, the control layer's dynamic will always match the one of the systems (except if different sample times are chosen for them). However, they will be useful in future releases of ODEHubs, which are expected to incorporate capabilities to simulate non-deterministic scenarios as well as different control strategies that deal with uncertainty.

ODEHubs's Block Library
The model presented in Section 2 can be configured both graphically and via code, and, in either case, the parameters that need to be introduced in ODEHubs are related to the elements of the above matrices, that is, the technical features of the devices that compose the EH. ODEHubs's block library will be used in this section to elaborate on that, since it might be more intuitive than programming code. Figure 3 shows the types of blocks that need to be used in the definition diagram, which were designed based on the elements found in the model presented in Section 2.2 and the appearance of the descriptive diagrams in many works related to EHss (see Figure 1, for example). They can be classified as follows.  • convergence Node and Divergence Node blocks ( Figure 5) do not currently have any utility and, in practice, could be omitted, but, because of the way Simulink R treats the signals between blocks, they will be needed in future releases of ODEHubs to allow users to merge and split flows numerically in simulations via Simulink R . For this reason, the diagrams must incorporate these blocks at any time that several flows converge into a single flow (so the sum of the original flows equals the resulting flow) or whenever a flow diverges into different sub-flows (so the sum of the sub-flows equals the original flow). Their only configurable parameter is the number of converging or diverging flows, respectively.
(a) (b) Figure 5. Parameters of the Convergence Node (a) and Divergence Node (b) blocks in Simulink R .

•
N o or less Storage System blocks (Figure 6a) need to be used to characterize the storage devices of the EH. The output of each of these blocks must be connected any Output block where energy, mass, or volume is stored. Their parameters, which are arranged in tabs, allow the definition of the charge and discharge flows and the capacity of each storage system, so they are elements of the vectors Q min ch , Q max ch , Q min dis , Q max dis , S min , and S max ; they also allow the definition of the efficiency of each storage process, taking part of the matrices C ch , C dis , and C s . Additionally, the initial state of the storage system must be defined S(0) in the Storage tab.
• N d blocks (Figure 6b) corresponding to devices need to be included in the Simulink R model employed to define the EH graphically. As depicted in Figure 3, ODEHubs's library consists of some predefined blocks with the same parameters: limits for the input(s) and output(s) of each device related to vectors D min di , D max di , D min do , and D max do , and the conversion factor(s) between them, which are required to build matrices C, C di , and C do . From the above parameters, the ones related to the model can be defined in one of the following ways: as a numeric data type in the case of static coefficients, as a string containing the name of a variable of the dataEH timetable corresponding to the value of such parameter at each sample time, as a string containing the name of a function in the folder user_models, which provides the value of this parameter at each sample time, or as a cell array composed of one or several of the latter ones, which is mandatory for those devices with more than one output, such as the Reversible Heat Pump, the Biomass Boiler, and the Combined Heat and Power blocks.

Main Configuration Parameters
Apart from setting the model's parameters, ODEHubs's users might want to customize the output of the program in order to suit their needs. Among the configuration parameters, one can select if third-party solvers are invoked, define the solver's options, select the start date for simulation (dataEH needs to be accordingly defined) and which results are presented (system response or control actions), and if they are exported as images or not. However, in the authors' opinion, the parameters with which new users could more easily struggle are those related to sample time and horizon, which are described below based on the interpretation of Figure 7.
• EH.simparam.nm defines the simulation horizon and, thus, depending on EH.simparam.tm, how many times the loop is executed to obtain the evolution of the system according to its dynamic, that is, updating the state of the storage systems through EH_system.m (which is still in its beta phase

Case Study
The EH employed in this how-to section on using ODEHubs consists of two of the plants from he CHROMAE project [21] (Figure 8 [43]. On the one hand, the parral-type greenhouse [44] has a total surface area of 877 m 2 (37.8 m × 23.2 m), a polyethylene cover, automated ventilation with lateral windows in the northern and southern walls, flap roof windows in each span, and fine-mesh anti-trip nets. The greenhouse orientation is east-west, the crop rows are aligned north-south, and the cropping conditions are very similar to the ones in commercial greenhouses, continuously monitoring several climatic parameters within the greenhouse. Outside, a weather station measures air temperature, relative humidity, solar radiation, liquid precipitations, and wind speed. In addition, the greenhouse has a heating system consisting of a GP 95 propane heater and a Missouri 150,000 multi-fuel boiler, which is part of the CO 2 storage and enrichment system as well [45].
On the other, the photovoltaic parking lot consists of an array of photovoltaic solar panels for production purposes, with 4830 Conergy PA 240P modules connected to 10 Fronius Agilo 100 inverters, and three other arrays for research, connected to their respective Fronius IG Plus 55v3 inverters and distributed as follows: 24

Parameters and Modeling
Both facilities are physically separate, but they will be considered a single ensemble that can be modeled by adapting Equations (1)- (16) to the structure represented in Figure 9. In fact, this is automatically done by ODEHubs, which generates the constraints of the optimization problem from the definition given by the users, but they have been included in Appendix A to clarify how the general model is particularized for a specific EHs. In addition, the input parameters required for the toolbox's blocks are described below before presenting the simulation results from ODEHubs.
Note that in Figure 9, the electric consumption of its 4.5 kW irrigation pump is considered separately from other devices' consumption in order to match it to the on/off state of the pump. Electricity, as well as the crop needs (water, heat, and CO 2 ), can be obtained from processing the data acquired by the sensor system of the greenhouse [44], whereas either technical data provided by the manufacturer of each device or one's own experimental data are employed to define most conversion and storage constraints of the EH. The production model of the photovoltaic parking was obtained by the authors in a former work, considering an isotropic sky model, the equivalent circuit for PV cells, and weather conditions [46]. This allows the determination of the instant availability of solar radiation and, consequently, the setting of the maximum input flow to the field (device 1), which has been scaled to its tenth part in order to analyze a realistic generation-demand situation and the conversion factor (performance) of the same one (Tables 1 and 2). Only the water, thermal, and electric storage systems in Figure 9 are hypothetical elements, whose parameters were estimated based on a former study [47]. Once it has finished entering the configuration parameters in the Simulink R model (Figure 10), ODEHubs will prepare the problem's setup and get the optimization results, as shown in Figure 2, which includes redefining Equations (1)- (16), as clarified in Appendix A. The parameters that need to be substituted in these equations come from various sources, and they can be classified into time-invariant or variable parameters. Tables 1-3 collect all of them, including their constant values for the first ones and the references of the models employed to periodically get the second ones. Note that the infinite symbol has been employed for unconstrained variables; thereby it is assumed that there is an inexhaustible supply of input resources. The same reasoning applies to the conversion devices, whose upper limits have been defined only for their inputs or for their outputs. There are some elements deliberately missing in Figure 9 (M 1 , M 2 , M 4 , M 5 , Q ch,5 , Q dis,5 , S 5 ), which implies that those flows do not take place in the EH, and the corresponding upper limits are set to zero by ODEHubs (see the fourth column in Table 1 and the fifth row in Table 3).

Eo
Ei (1) Eo (1) Rpk (2) Pi (3) Bi (4) Wi (5) Hgh (2) Go (3) Wo (4) Esgh ( Figure 10. Case study EH defined through ODEHubs's blocks in Simulink R . The labels within some blocks are customizable for representation purposes, although by default, ODEHubs names each input and output to present the results according to the number between parentheses, matching the notation in Figure 9. ODEHubs also renames the Storage system blocks according to the label of the Output block to which they are connected. In this diagram, the labels within Input and Output blocks and the name of each device have been chosen to identify the kind of resources and device, but they are related to Figure 9 as follows. Devices: D 1 → PVpk (photovoltaic parking lot), D 2 → PBgh (propane heater), D 3 → BBgh (biomass boiler), D 4 → IPgh (irrigation pump). Inputs: I 1 → Ei, I 2 → Rpk, I 3 → Pi, Note that M 3 is missing because the block Go has been configured to allow sale of the resources attached to that output (see Figure 4b).
Regarding the remaining parameters of the model, propane (c 3 = 1.694 EUR/kg), biomass (c 4 = 0.255 EUR/kg), water (c 5 = 0.547 EUR/m 3 ), and electricity prices are obtained from local supply companies. In the case of the electricity (c 1 ), three time periods are distinguished to define the cost: 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. In contrast to the old regulatory framework [41], solar radiation is not taxed (c 2 = 0). Moreover, as the only resource "sold" is the CO 2 freely released to the atmosphere (M 3 ), the price of all the sold resources is null: s(k) = 0 0 0 0 0 .  11.54 η s,2 0.94 η ch,2 0.9 η dis,2 0.9 η D, 3 4.25 η s,3 1 η ch,3 1 η dis,3 1 η D, 4 1.76 η s,4 1 η ch,4 1 η dis,4 1 η D, 5 1 η s,5 1 η ch,5 1 η dis,5 1 Table 3. Storage systems' upper and lower limits. Finally, the weather and greenhouse data correspond to the data collected at the said experimental station on 17 December 2018. The choice of that day was made based on the availability of solar radiation (sunny intervals and temperatures between 11 and 21 • C, as shown in Figure 11) and the fact that CO 2 enrichment and heating support were active during part of the day, which results in a scenario where the outputs of the EH have a non-null demand profile (electricity and water are demanded daily, so they are not critical to select the dataset).  Figure 11. Weather conditions on 17 December 2018. Apart from global horizontal irradiance (GHI), direct normal irradiance (DNI) and diffuse horizontal irradiance (DHI) determine the amount of power produced by the photovoltaic field, and temperature (a) allows the estimation of the performance of the cells. DNI requires a different scale (a), with respect to GHI and DHI (b), owing to the direction of the normal beam of radiation, which results in higher values of the irradiance around midday. As data are recorded every minute, ODEHubs calculates the hourly average to simulate the case study.

Simulation Results
The optimal dispatch for the greenhouse was obtained by running ODEHubs with the solver intlinprog [48] and scheduling mode (fixed horizon) on an Intel R Core TM i7-6700K 4GHz CPU, taking 0.3219 s to solve the optimization problem. That value means that the tool would probably work as well with shorter sampling times or using MPC mode (receding horizon), which would increase the computational burden. ODEHubs outputs a structure variable that contains all the data defined by the user and post-processed variables, for instance: general parameters of the EH, vectors and matrices of the model, constraints for the solver, and the values of the decision variables.
Although anyone can consult the aforementioned structure, by default, the results are presented in as many graphics as the number of EH's outputs (N o ), arranged in the same way (see Figure 12). Scaled to the left axis and expressed in terms of mass flow, volume flow, or power, the elements of O (demand profiles) are depicted by a thick line together with the elements of M, which are stacked over these latter ones and depicted by a thin line. This avoids overlapping between them, as they constitute the left term of Equation (1), which needs to equal the right term, represented with colored stacked bars. These bars actually do not represent any of the variables of the model, but they are obtained from the product of C and P with the purpose of indicating the contribution of each input flow to meeting the demand of resources (O). Scaled to the right axis and expressed in terms of mass, volume, or energy, the evolution of the storage systems (elements of S) is represented by a dashed line. In this regard, bear in mind that neither the charge (Q c ) nor discharge flows (Q d ) explicitly appear in the charts, but they can be deducted through either the difference between the bars and the stacked solid lines or the slope of the dashed lines. By default, ODEHubs sets the upper limit set equal to the maximum capacity of each storage device, but, for the sake of clarity, in this case, they have been customized to show each variable with an appropriate scale. Note that the storage systems are charged whenever the bars are over the stacked solid lines (demand plus sales profiles), or, equivalently, whenever the slope of the dashed line is positive and, by similar reasoning, they are discharged whenever the bars are under the stacked solid lines or the slope is negative. Nonetheless, when analyzing a negative slope, the presence of degradation must be also taken into account, since it has a similar effect. The color choice for each kind of resource is made by the users, and the appearance of the horizontal axis is determined by the sample time, as explained in Section 3.
From Figure 12, which shows the charts of the case study previously presented and the cumulative variables included in ODEHubs's output structure, the results can be summarized as follows: The 7.02 kWh of electric energy demanded (greenhouse and pump) is covered by purchasing 3.81 kWh from the power grid (0.39 EUR) and the photovoltaic field, which produces energy whenever there is enough solar radiation (see Figure 11); the 18.47 kWh of thermal energy is entirely generated in the biomass boiler because of its lower cost (1.46 EUR) in comparison with propane and its simultaneous use for producing the required 4.69 kg of CO 2 ; the 0.87 m 3 of water (0.48 EUR) is purchased in a single period, between 9:00 and 10:00, in which the availability of solar energy makes up for the slightly higher cost of the electricity at that time (0.1127 versus 0.0892 EUR/kWh at off-peak hours).

Conclusions
In view of the results, it can be concluded that ODEHubs is a solid alternative to consider for the formulation and analysis of the dispatch problem in agri-energy hubs, such as that of CHROMAE. Although, for the sake of conciseness, only a specific case study has been evaluated in this paper, it is expected that the tool will be useful in studying resource management strategies with different sets of real data in agro-industrial districts and, because of its conception as a flexible tool, even in industrial, commercial, or home environments. There is still a tradeoff between an accurate representation of some nonlinear processes and the actual formulation as a mixed-integer linear programming problem, which ensures the convergence to a global optimum while keeping the computational burden low. The current solution, however, allows the introduction of variable conversion coefficients, but the constraints could be adapted to consider nonlinear efficiencies due to partial loading [49].
To the best of the authors' knowledge, many of the control systems proposed in the literature come from the field of chemical engineering, which has traditionally used a hierarchy of levels or layers to distinguish different tasks performed in industry [50]. These are defined attending to the frequency at which decisions are made (i.e., how often the control actions are calculated), the kind and number of equipment involved, or the level of abstraction. Many contributions employ hierarchical control strategies in which an MPC can act in either the regulatory or the supervisory layers as well [51], and, in this regard, there is a dichotomy between Economic MPC and Real-Time Optimization with MPC [52]. Long [40] discusses both approaches in his thesis [40], as well as how the existence of uncertainty affecting the system behavior is sometimes handled by using certainty equivalent model predictive control [53]. From his ideas, one can conclude that, when it comes to certain energy management applications, it is acceptable to assume (a) that the lower-level regulatory controllers have a dynamic fast enough with respect to the sampling interval of the upper layers, and they do not need to be considered for optimization purposes (the corresponding steady states are instead considered), (b) that due to their presence, stability is not a primary concern, and (c) that for demonstration purposes, a general uncertain optimization problem can be converted into a deterministic one. This is the basis of many of the articles found in the literature about energy dispatch that deal with the energetic management of a campus [54] or the energy dispatch of an MG [55], to name a few, and is also the basis of ODEHubs in its current form, whose output corresponds to the computations performed in upper control layers.
Future implementations might incorporate: the interconnection of EHss, which is under development [56], so agents with opposing interests can be considered in the dispatch problem; uncertainty in some of the variables (especially the demanded resources) and arrangement of the consequent probabilistic formulation of the problem; integration of the dispatch problem into the design of facilities, for example, using a bi-level approach in which a genetic algorithm generates possible configurations of the EH [57]; addition of capabilities to incorporate elements from Simscape Simulink R 's library, which includes components to model different energy networks, such as power grids [58]; or migration of the package to other languages, such as Python or Modelica R . Python programs are executed from a command line, similarly to MATLAB R scripts from the command prompt, and there are already some tools built in that language [25], and Modelica R is an object-oriented modeling language that has already been used to implement models of some of the components of the CHROMAE project [59]. Funding: This work has been funded by the National R+D+i Plan Project DPI2017-85007-R of the Spanish Ministry of Science, Innovation, and Universities and ERDF funds.

Acknowledgments:
The authors wish to thank the Experimental Station of the Cajamar Foundation and the CIESOL Research Center for providing real data of their facilities. They thank as well the reviewers for their valuable comments, which have helped to improve the quality of this paper.

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

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Greenhouse Model in ODEHubs
The energy hub in Figure 9 counts with five inputs (N i = 5), five outputs (N o = 5), eight paths (N p = 8) between inputs and outputs (see Table A1), four conversion devices (virtually divided into five because of the two outputs of the biomass boiler) with a total number of four inputs (virtually divided into five because of the two outputs of the biomass boiler) N di = 4, and five outputs N do = 5. Thus, Equations (1)- (17) are defined according to Equations (A1)-(A18): s(k) = s 1 (k) s 2 (k) s 3 (k) s 4 (k) s 5 (k) .
In addition, note that Equation (13) remains unaltered, Equations (14) and (15) are removed from the constraints since there are no resources that can be acquired and sold nor devices that cannot operate at the same time, and Equation (16) turns into Equation (A19): where P 6 and P 7 correspond to the input flow of the biomass boiler. Note that the capacities' constraints are not explicitly included in this appendix, as they can be easily deduced from previous equations and Tables 1-3.