Energy Packet Networks With Energy Harvesting

We investigate the cooperation among energy prosumers (unified energy provider and consumer) through the energy packet network (EPN) paradigm, which represents both the flow of work that requires energy, and the flow of energy itself, in terms of discrete units. This paper details a stochastic model of EPNs, which is inspired from a branch of queuing theory called G-networks. The model allows us to compute the equilibrium state of a system that includes energy storage units, energy transmission networks, and energy consumers, together with the intermittent energy sources. The model is then used to show how the flow of work and energy in the system can be optimized for certain utility functions that consider both the needs of the consumers, and the desire to maintain some reserve energy for potential future needs.


I. INTRODUCTION
In the future, billions of computer and communication devices may work through the Internet of Things (IoT) to manage and solve locally the major challenges of cities and human communities [1].Applications will include energy provisioning and the smart grid [2], physical safety and security [3], health [4] and the environment [5].In particular, the smart grid will itself require innovative communication solutions [6]- [10], but distributed data sources require synchronisation [11] while communication protocols create additional overhead [12] and hence energy consumption.
However energy and workload monitoring and management systems will also consume energy and add to societal issues such as CO 2 impact [13].Thus in order to move forward towards ubiquitous, self-sustainable and optimised systems such as smart cities, sustainable agriculture, energy microgrids, and local energy supply systems, with desirable properties such as reduced emissions and energy transmission losses, it will be necessary to power the IoT as much as possible through sources of renewable and clean energy and meet the energy needs of the communication networks, servers and data centres that are needed to support them [14].
In the meanwhile, the worldwide increase in electrical energy consumption for ICT [15] worldwide in the order of 5-7% per year, has also motivated research in energy harvesting [16]- [19] for computer and communications systems.Furthermore power is being increasingly generated in a distributed manner with consumers becoming ''prosumers'', who may also be producers and/or ''storers'' of energy (e.g. in home-based battery systems, or in individual vehicles).Thus there is an increasing motivation to dynamically manage electrical energy consumption, in conjunction with its generation and distribution.Such dynamic management requires the integration of consumers with communication networks and data centres so that energy management software may dynamically optimise the power grid.
These evolutions encourage us to study integrated models for renewable energy systems in the presence of dynamic and flexible energy consumption [20].In this respect, a variety of emerging technologies such as Intelligent Power Switches, PowerComm Interfaces [21] and wireless energy transfer [22] can assist in dynamically managing, storing and conveying electrical energy.Furthermore, energy storage in conjunction with harvesting, energy sharing between prosumers, and workload migration between certain energy consumers such as data centres, can be used to improve the economic cost and CO 2 impact of the data centres and computer networks that are required to support the energy systems of the future.Since energy harvesting will in general be intermittent, it will useful to store energy dispatch it dynamically to optimise overall system behaviour [23], [24].

A. CONTRIBUTIONS OF THIS PAPER
Thus this paper addresses the use of the ''Energy Packet Network'' (EPN) paradigm [25]- [27] to model and optimise the dynamic generation, storage and consumption, as well VOLUME 4, 2016 This work is licensed under a Creative Commons Attribution 3.0 License.For more information, see http://creativecommons.org/licenses/by/3.0/as the dynamic workload, of highly dynamic and distributed energy systems.Specifically, it addresses the use of the EPN, to model and optimise the generation, storage and consumption, as well as the dynamic workload, of highly dynamic and distributed energy systems, and presents several new results: • A novel case of practical interest is first developed, when we take a specific ''high level'' decision regarding the distribution of energy.In this decision, all energy stores should, on average, contain the same amount of energy, so that if any of them fails the others can best support the system, and the backlog of work at each of the workstations should on average be the same.
• Then we introduce specific utility functions for the general case, that take into account both the needs of the consumers, and the desire to maintain some reserve energy for future needs.
• The control variables considered are the rates at which we can move energy packets from energy stores to workstations, and also move energy between different energy stores.The model also allows us to represent the movement of jobs between workstations.The utility functions contain contradictory terms, so that as we vary the control variables, some of these terms increase while others decrease, leading to clear optima.
• In order to optimise these utility functions, we develop an algebraic formulation of the partial derivatives of the quantities of interest with respect to the control variables.This algebraic formulation is used to derive the computational algorithms that obtain the optimum operating points of the utility functions.
• A significant numerical example is also presented.It concerns a remote sensing/processing centre operating a (radio) telescope or radar system, together with a sensing and life support system, which constiute three ''workstations''.Electrical energy is provided by a photovoltaic source and by wind.Within the context of this system, the previously defined utility functions are considered and detailed optimisation results are illustrated numerically in several tables and figures, regarding the optimum flow of energy to each of the three workstations.The paper ends with conclusions including suggestion for further work.

II. ENERGY PACKET NETWORKS
The EPN model for networks of energy prosumers combines: • The intermittent flow of harvested of energy, represented by a random flow of arriving ''energy packets (EPs)''.An energy packet is a discrete unit of energy, e.g.X Joules, which is represented as arriving in ''one chunk''.It is simply a simplified version of the actual flow of energy which is continuous.
• The storage of energy (say in batteries), again in discrete units of EPs.An energy storage unit (ES) such as a battery is modelled as a queue of EPs that are waiting to be used.The ES is replenished by a flow of EPs from some external source including an energy harvesting unit, and it can be depleted both when energy is forwarded to a consumer, and through losses that represent leakage or line losses.
• The sources of power, the ESs and the consumers are interconnected by Power Switches (PSs) which dynamically connects the sources of power to the ESs, and the sources of power and ESs to the consumers.
• The consumers will request for EPs from either a PS or ESs, and these requests will typically be intermittent, since they are a function of the work that these consumers accomplish with the energy.A typical example of such consumers are ICT systems which intermittently receive computational work to accomplish (for instance in terms of programs that need to be executed, or in terms of data packets that need to be transmitted), and which in turn require energy to accomplish this work.Thus a consumer is also represented by: • A queue of work that it has to accomplish which may originate from some outside source or from some other consumer of energy which accomplished a prior work step, and • One or more servers that have job execution times which have a random duration which may depend on the nature of the job and the speed of the work unit, and also on the flow of power to the server.An EPN will contain multiple work units with external arrivals of jobs or data packets, the jobs or data packets may move from one consumer queue to the other, and there will be multiple queues of EPs, with external arrivals of renewable or steadily flowing energy, and energy itself may flow between ES, or it may leak from the energy stores, or move to consumers where it is used to execute jobs or transmit packets.Such EPNs can be analysed using a branch of queueing network theory called G-Networks [28], [29].Typically in queueing theory, jobs move around a set of servers where they queue up to receive service until they complete their needs and leave the system.In G-Networks on the other hand, some of the flows of ''customers'', in this case the EPs, can queue up (as in an ES), but they can move to act as servers or activators for other customers, in this case the jobs that need to be executed or the data packets that need to be transmitted, that queue up at the servers where work is accomplished.Thus the G-Network can be used to represent job flows that consume energy and energy flows that allow jobs to be executed.
Electronic systems that capture a similar physical behaviour have been recently built and experimentally tested under the name of ''power packet'' systems.They are designed for the smart dispatching of electrical energy [30], [31], and in some approaches the forwarding of an EP takes the form of a pulse of current with a fixed voltage and duration, and each EP is equipped with an encoded header compromising the destination information for the EP [32] that allows PSs to determine where the energy needs to be sent.
In previous work [33] we have used the EPN model to study the architectures that interconnect energy prosumer systems so that energy losses and the response time to service request are minimised.
In this paper we develop optimisation algorithms that optimise utility functions that include a linear combination of the throughput or useful, and the the probability that the system does not run out of energy.Thus the type of utility that we seek to maximise includes two sets of terms: the sum of the job processing rates of the system, and the sum of the probabilities that there is reserve energy in the system.Clearly if we try to just maximise the work done we may run out of energy, while if we try to just save energy then the system may delay the job execution times.

III. THE MATHEMATICAL MODEL
We study an EPN consisting of a finite set of nodes.It contains m energy storage (ES) nodes denoted E l , 1 ≤ l ≤ m store and dispatch EPs, while n workstation (WS) nodes W i , 1 ≤ i ≤ n are ''places of work'' where jobs are executed.The system diagram is described schematically in Figure 1.Energy packets enter the energy store E l at rate γ l from external renewable sources.Also, the energy stored in E l may be lost at a rate π l due to leakage and standby losses, and E i sends energy packets out of the store at rate r l either to other energy stores or to workstations.
The external arrival rate of jobs to the work station W i is λ i , while the workstation may also receive jobs from another workstation W j with probability P ji after workstation W j finishes processing that job.Moreover, uncompleted jobs stored at the work station W i might also be lost at a time-out rate β i .
However processing at W i can only occur if an energy packet arrives at the workstation from some ES E l .This can happen in two ways: • With probability s(l, i), E l sends an energy packet to W i , and it may send EPs to E k with probability S(l, k), with n l=1 s(l, i) + m k=1 S(l, k) = 1 for any l.We will write w(l, i) = r l s(l, i) and W (l, k) = r l S(l, k) so that r l = n l=1 w(l, i) + m k=1 W (l, k).• With probability P ij the workstation W i sends the job that has just completed to node W j for more processing, or with probability d i = 1 − n j=1 P ij the workstation W i will simply remove the job that terminates at W i from the system without forwarding it to another workstation.In fact a more complex model in which a single job can consume more than one energy packet at a time can also be used, and will be considered in future work.
As a consequence of the above assumptions, using G-Network theory [28], [34] the probabilities Q i that W i has at least one job in its queue, and q l that E l has at least one EP in its ES, will satisfy: where the rate of service at W i is determined by the rate at which it receives energy, and q l is the probability that energy storage station E l has at least one EP in store: The above equations also allow us to state the following result [34]: Theorem 1: Assume that the external arrivals of jobs, and the arrivals of EPs from renewable energy sources, are independent Poisson processes with rates λ i , and γ l .Suppose that the r l , π l and β i are the parameters of mutually independent exponentially distributed random variables.Then if (K 1 , . . ., K n ) ≥ (0, . . ., 0) and (k 1 , . . ., k m ) ≥ (0, . . ., 0) represent the backlogs of jobs to be executed at the workstations, and the number of energy packets stored at the ESs, respectively.Then Corollary 1: Because in steady-state the joint probability distribution of the number of jobs waiting, and energy packets stored, at each of the WS and ES, is the product of the marginal distributions, we have: A. VECTOR REPRESENTATION Denote by Q and q, respectively, the row vectors whose elements are (1) and ( 2).We can then write: where • w q is the n × n diagonal matrix whose diagonal elements are m l=1 q l w(l, j), • and the , are, respectively, the n and m dimensional vectors of elements The expression (3) can then be readily written as:

IV. AN INTERESTING SPECIAL CASE
A case of practical interest arises when we take a specific ''high level'' optimisation decision, namely: • All energy stores should, on average, contain the same amount of energy, so that if any of them fails the others can best support the system, so that: • The backlog of work at each of the workstations should on average be the same, i.e.
Then we have: and we know that w + (i) − n j=1 w + (j)P ji ≥ 0, and we recall that Q * , q * must satisfy 0 < Q * < 1 and 0 < q * < 1.If q * and Q * are fixed, we obtain: so that we must have: Then writing the m-vector c with elements: we have: which has a unique solution because the Markov chain P is transient.Since all the q k = q * are identical, we set so that all the w(k, i) are now determined.As a consequence, for each 1 ≤ k ≤ m we have: Turning now to the set of weights W (k, l) we have: where: and for 1 ≤ k ≤ m we have:

V. UTILITY FUNCTIONS AND OPTIMISATION
Simpler EPN optimisation problems, than the ones studied here were considered earlier in [35].The optimisation would be conducted on the assumption that: • The r l are variable, but they have an upper bound that represents the maximum amount of power that can be delivered by the energy stores.
• The P ij are fixed and represent the sequence of job steps that the jobs have to go through.
• The S(l, j) and s(l, k) are the ''control variables'': they are modified or selected so as to maximise U .
• Cases can be considered where the S(l, k) are also fixed.Obviously, we wish to limit the backlog of work at the workstations, while we also want to have some reserve of energy since there may be unpredictable needs.Thus a sensible cost or utility function would have the form: which is the weighted probability that the backlogs of work exceed the values K i at each workstation i, and the weighted probabilities that there are not at least k l energy packets at ES l, where the a i and b l are non-negative real numbers (the weights).A simpler utility that may be minimised is the weighted probabilities that there is some backlog at the workstations, plus that there is no energy in reserve, which is U 1 when K i = 1 and k l = 1 for all i, l: A cost function other than U 1 defined above that we may wish to minimise is: which differs from U 1 only in the first sum which is simply the weighted sum of the average response time of jobs waiting to be served at each of the workstations.On the other hand, the following utility function U 3 needs to be maximised since it establishes a balance between the throughput of the system (the first term) and the probability that some energy is kept in reserve:

VI. COMPUTING PARTIAL DERIVATIVES OF THE UTILITIES
The optimisation of the such utilities, where the utilities are continuous and differentialble functions of the control variables, will require the computation of the derivatives of the utility functions with respect to the control variables.In this case, the control variables are the routing probabilities s(l, k) and S(l, i) for the EPs.In fact because it is easier to work with non-negative real numbers of arbitrary size rather than with probabilities, we will consider that the control variables are the quantities w(l, k) and W (l, i).Thus, for some real valued function V of the w(l, i) and W (k, l), l, k ∈ {1, . . ., m} and i ∈ {1, . . ., n}, let us use the notation:

A. THE CASE WHERE THE r l ARE CONSTANT
In a certain number of circumstances, we can consider that the r l are constant, for instance when each of the ESs has a fixed and constant rate at which each of them can be emptied of its EPs.Here we will first focus on this case.= 0 for any l, a, b, it is easy to see that: and As a direct application of Remark 1, we construct the following lemmas that will be used later in the optimising algorithms.
Lemma 1: The derivative of q l with respect to any w(a, b) for l, a ∈ {1, . . ., m}, b ∈ {1, . . ., n} is given by: or in vector form: where • q w(a,b) is the m-vector each whose elements are the derivatives of the q l , • W is the m × m matrix of elements W (l, k)/(r k + π k ) and The derivative of q l with respect to any W (a, b) for l, a, b ∈ {1, . . ., m} is given by: where δ b = [0, . . ., 0, 1 r b +π b , 0, . . ., 0] is the m-vector with 0 elements everywhere except in the b-th position which has the value (r b + π b ) −1 .
Lemma 3: The derivative of Q i with respect to any w(a, b) for a, b ∈ {1, . . ., m}, i, b ∈ {1, . . ., n} is given by: Lemma 4: The derivative of the Q i with respect to any W (a, b) for a, b ∈ {1, . . ., m} is given by:

B. THE CASE WHERE THE r l ARE NOT CONSTANT
In other circumstances the r l will not be constant, so that when we modify any of the w(l, i) or W (l, k) we also change the r l without affecting the other parameters.
This may occur for instance when the ESs are organised as a stack of batteries in parallel so that their output power flow can be varied as a function of the number of storage units that are switched to the power output bus.
However in practical circumstances there will be a maximuml value of r M l ≥ r l which cannot be exceeded because of the physical limitations of the devices that are used to store energy and of the power circuits that carry the power flows.= 0 for any l, a, b, it is easy to see that: and As a consequence, when the r l < r M l , l = 1, . . ., m we have the following lemmas.
Lemma 5: The derivative of q l with respect to any w(a, b) for l, a ∈ {1, . . ., m}, b ∈ {1, . . ., n} is given by: or in vector form: where as before W is the m × m matrix of elements W (l, k)/(r k + π k ), and 1 l is the m × m matrix which is zero everywhere except in the diagonal terms which have the value 1 r l +π l .Lemma 6: The derivative of q l with respect to any W (a, b) for l, a, b ∈ {1, . . ., m} is given by: where δ b = [0, . . ., 0, 1 r b +π b , 0, . . ., 0] is the m-vector with 0 elements everywhere except in the b-th position which has the value (r b + π b ) −1 .
Lemma 7: The derivative of Q i with respect to any w(a, b) for a, b ∈ {1, . . ., m}, i, b ∈ {1, . . ., n} is given by: The derivative of the Q i with respect to any W (a, b) for a, b ∈ {1, . . ., m} is given by:

VII. OPTIMAL SOLUTION WITH GRADIENT DESCENT
Since the optimisation of the EPN optimisation can be expressed as the minimization or maximisation of utility functions such as those defined in Section V, it can be achieved by selecting the appropriate EP flow rates w(a, b) and W (a, b) when other system parameters are fixed.Since we are delaing with continuous and differentiable utility functions, the gradient descent algorithm is a useful tool in this case.
At a given operating point of the EPN X = [λ, γ , r, P, π, β], the gradient descent algorithm at its t th computational step is: where |η| is the rate of the gradient descent, we set η < 0 fo rthe minimisation of utility functions U 1 , U 0 1 , U 2 , while η > 0 forthe maximisation of utility function U 3 .
In practice, we will be interested in a gradual optimisation of the system, where we modify parameters progressively, in a system that should operate normally and adapt the ongoing work work and energy flows.
Thus we compute the partial derivative of the utility functions as: similarly, another cost function which we would like to minimise is: where Moreover, the partial derivative of the utility function which needs to be maximised is computed: The partial derivative with respect to the W (a, b) can be computed similarly.The w(l, i) w (a,b)  b) , q W (a,b) are detailed, for both constant and variable r l , in the previous section.Thus, the steps of the gradient algorithm are: • First initialize all the values w 0 (a, b), W 0 (a, b) and choose η, • Solve the system of non-linear equations given in ( 1) and ( 2) to obtain the steady state probabilities, • Calculate the partial derivatives as given in Section VI using the steady state probabilities computed in the previous step.

VIII. A SET OF NUMERICAL EXAMPLES
To illustrate the proposed approach, we have constructed a numerical example for a facility, such as a remote sensing and monitoring facility such as an astronomical observatory on a mountain peak, or a border control station in a remote part of the world, that is powered by energy harvesting devices such as solar panels and wind turbines.The arrival rate to the workstations is in units, each of which can be executed or carried out with exactly one EP.Thus the relative workload for the different workstations is expressed by the arrival rate of work to each workstations.In the example we will consider there are three distinct workstations carrying out distinct work so that they can transfer work from one workstation to another and the W (a, b) = 0 for all a = b.The work that uses energy is carried out in the three workstations: • W 1 is meant to model a radio (or optical) telescope or radar station that is operating continuously and requires substantial processing, represented by a flow of work arriving at rate λ 1 .We assign to it a priority a 1 ≥ 0.
• W 2 represents the infrastructure energy needs such as lighting, air-conditioning etc., that receives work requests at rate λ 2 , with a 2 > 0. Since this is a lifesupport system we would expect it to have high priority.
• W 3 is mean to model the perimeter security and surveillance, and the monitoring of weather conditions, with workload arriving at rate λ 3 .Since this represents the life-support with a priority a 3 > 0, which may be comparable to that of a 2 , while the priority a 1 may be lower.Radars or telescopes provide a substantial amount of imaging data per unit time to be processed and transmitted, thus λ 1 can be considerably larger than λ 2 and λ 3 .Moreover, security and atmospheric monitoring devices typically forward limited motion sensing, video camera, temperature and wind data, thus λ 1 > λ 2 > λ 3 .Likewise, the priorities of these different work stations are chosento be a 3 ≥ a 2 > a 1 .To reflect different forms of energy harvesting, such as photovoltaic and wind, we have two ESs, E 1 and E 2 which store energy from two different renewable sources at a rate γ 1 and γ 2 EPs per unit time.The system model for this example is summarised in Figure 2. To obtain significant results and to remain within the stability region for the model, we must operate with parameters that represent a balance between the energy flow and the flow of workload.Also, the impatience of jobs, represented by the parameters β i , can help avoid instability that occurs when the job queues increase in size and the waiting time of jobs at the execution queues become very large.Similarly, the energy leakage will also reduce the amount of energy that is stored when there is no work to perform.Note, however, that we will typically have π i r M i to represent the fact that leakage rates are only a small fraction of the maximum power (energy/time) that the ESs can offer.A plausible scenario is constructed with a set of numerical parameters shown in Table 1 in order to represent a plausible scenario.We first consider the case where the r l values are not constant, and the case where they are constant.
A. WHEN r 1 , r 2 ARE NOT CONSTANT Consider the case where the r l values are not constant, but a maximum practical value r M l is chosen as in Table 1.We set P(a, b) = 0 for a, b ∈ {1, 2, 3} so that jobs are not moved FIGURE 2. We use the EPN to model the energy and work flows in a remotely located sensing facility, such as an astronomical observatory or radar station, with two energy storage units, two renewable energy sources (e.g. one photovoltaic and one wind) and three main energy consuming workstations: the major instrument (e.g. a radar), the life support systems (e.g.air conditioning, lighting), and the atmospheric sensing and perimeter security sensors.between the different workstations.Also in the numerical examples we set s(l, i) ≥ 0.05 , i = {1, . . ., 3}, l = {1, 2} and 3 i=1 s(l, i) ≤ 1, l = {1, . . ., 2} to make sure that each WS receives some minimum amount of energy to avoid ''starvation'' that would lead to an infinite backlog of work at some workstations, and n i=1 w(l, i) ≤ r M l , l = {1, . . ., m}.The unconstrained optimisation problem is solved numerically for the four utilities U 1 , U 0 1 , U 2 , U 3 of Section V, and the resulting values for the optimum power flows are given in Tables 2, 3 Table 5 shows that in order to maximise U 3 and minimise U 1 , U 0 1 , U 2 , the ESs should be providing energy at their maximum total rate, i.e.This is illustrated in Tables 6, 7, 8, 9 where we see that depending on the requirements of optimisation, a fraction of the energy sent out by the ESs can be forwarded for storage to the other ES.Note however than in this example, the energy loss during transfer, i.e. the ohmic loss, is not taken into account.Clearly, if energy transfer losses become significant then an optimised system will shy away from making additional energy transfers which do not lead to direct consumption.

IX. CONCLUSIONS
In this paper, we have discussed the main properties of the Energy Packet Network model as a framework for a system where intermittent distributed workloads and intermittent sources of energy operate jointly with multiple energy consuming workstations, and with interconnected energy storage units, where both work and energy circulate.We obtain the equilibrium probability distribution for both the backlog of work and the backlog of energy throughout the system.An interesting special case is also considered where the work and energy flows have been designed so that all energy stores and work backlogs are the same.Other stochastic models that are more adapted to the study of wireless sensors with small scale energy harvesting [6], [16], [17], [36] are considered elsewhere [37].
Here we have studied several utility functions that describe the overall system performance, and an algorithmic method based on gradient descent (or ascent) of these utility functions is developed to seek energy and work distribution mechanisms that optimise the utility functions.Also, the gradient descent approach we have described requires matrix inversions to compute Q w (a,b) and q w(a,b) which are of computational complexity O(n 3 ) and O(m 3 ), respectively.Simpler numerical schemes, such as power expansions of the matrices, can result in useful approximations, especially when the energy and work flow rates vary frequently with time, or are imprecisely known and will be studied in future work.
Another direction that will be considered includes more complex models in which a single job can consume more than one energy packet at a time, or where an energy packet may service several jobs at a time.The approach studied in this paper can also be extended to multiple classes of EPs, with power colouring which can be used to represent the origin of the EPs, for instance different types of harvested energy from wind or photovoltaic, as well as from non-renewable or fossil sources, or with different economic cost.
Multiple classes of workloads depending on their energy needs, or their priorities, or their importance with regard to the income they produce, can also be considered.We will also examine ''good'' but suboptimal solutions which can be obtained at a lower computational cost.Finally, we will also study the dynamic behaviour of such systems when the flow of work and of energy is regulated in real-time, using techniques such as reinforcement learning.

FIGURE 1 .
FIGURE 1. Schematic representation of an EPN with two types of nodes and ''queues'': the ''positive'' nodes contain the work to be done while the ''negative'' nodes are the ESs which store energy and are replenished by harvesting.ESs provide power in the form of EPs to the ''positive'' nodes.Work can move from one positive node to another where they are processed when energy becomes available, and finally the work leaves the network after completing a certain number of work steps.EPs can leak from the ESs (negative nodes), or be transferred to another ES or they can be transferred to a positive node to accomplish work.
i w(l, i) * = r M l , when b 1 , b 2

TABLE 1 .
Parameters used in the numerical examples.

TABLE 2 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 100, the r 1 , r 2 are not constant and the parameters are given in

Table 1 .TABLE 3 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 10, the r 1 , r 2 are not constant and the other numerical parameters are in Table1.

TABLE 4 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 1, the r 1 , r 2 are not constant and the other numerical parameters are given in Table1.

TABLE 5 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 0.1, the r 1 , r 2 are not constant and the other numerical parameters are given in Table1.
are considerably smaller than a 1 , a 2 , a 3 .However when the b 1 , b 2 are larger, we observe that the optima occur when the sum i w(l, i) * lies between 0.15 * r M l and r M l depending

TABLE 6 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 100, r 1 , r 2 are constant and the other numerical parameters are given in Table1.

TABLE 7 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 10, r 1 , r 2 are constant and the other numerical parameters are given in Table1.

TABLE 8 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 1, r 1 , r 2 are constant and the other numerical parameters are given in

Table 1 .
on the values of b 1 , b 2 .Note that the second terms in the utility functions concern the reserve energy contained in the ESs.
B. WHEN r 1 , r 2 ARE CONSTANTWhen the r 1 and r 2 are constant, energy leaves each of the ESs at a constant rate.However not to waste energy, we will allow some of it to be sent to another ES rather than to the WS, if the optimisation requires it, i,e, the P(a, b), a, b ∈ {1, 2} are now unconstrained and can be positive.

TABLE 9 .
Optimised EP flow rates for the different utility functions when b 1 , b 2 = 0.1, r l are constant and the other numerical parameters are given in