New Approach for Stochastic Modelling of Microgrid Containing CHP

The focus in this article is on long term load forecasting and expansion generation planning of the microgrids containing multiple load type consumption. The novel approach of planning cogeneration systems by introducing the function of simultaneous heat and electricity duration is presented. The function is extended form of load duration curve in such a way as to capture the simultaneity of different type of energy load. The paper contains seven chapters. In the second chapter the load forecasting methods are summarized. In continuation as an introduction to the main goal of this article, definition, application and restriction of load duration curve is presented. In the fourth chapter, some outlines about cogeneration system planning are given. The new approach for planning of microgrid contacting CHP based on probabilistic approach is presented in the fifth chapter. In the sixth chapter are demonstrated some application functions of simultaneous heat and electricity duration as a tool for planning multi energy systems.


INTRODUCTION
The technical system lifetime cycle can be divided into several phases such as Planning, Design, Constructing, Exploitation and Decommission phase.Each of these phases can be further subdivided into more sub phases.The lifetime system quality mostly depends on planning phase quality.Planning can be defined as a complex intellectual process of determining system structure and operation strategy to accomplish the objectives in the future.
Planning requires knowledge about performance of the system under observation and prediction of the environment outside the system.There are numbers of scientific and engineering papers devoted to planning mostly in electric power systems both large scale and microgrids [1,2].A microgrid can be defined as a low or medium voltage distribution network comprising various DGs, storage devices, and responsive loads that can be operated in both grid-connected and islanded modes [3].The planning phase starts with problem definition (defining customer demands, the goals, inputs, outputs, system constrains, key performance indicators, investment evaluation criteria, etc.).In most cases, technical system under planning is going to be modelled (real world is simplified and presented by various mathematical formsmodels) [4].A special planning and modelling technics are developed for such system under uncertainty [6][7][8][9][10][11].The level of uncertainty of inputs and outputs significantly shapes models and determines the quality of planning.Unlike the large-scale power electric system, small scale system such as microgrid has dominantly uncertain outputs.Also, unlike the non-renewable power systems the renewable power systems are examples of the systems with uncertain inputs [5].The planning of such power systems involves various stochastic modelling and simulation methods [4,6,7].The survey of stochastic modelling of microgrid is presented in [14] and [15].A powerful mathematical tool for stochastic system analysing is the Monte Carlo Methods [16,17].
An important phase in the modelling and planning is gathering the historical data and its usage for demand forecasting [18].In the case of the planning of power supply systems (mechanical, electricity or/and heat) the balancing between power production (supply) and power consumption (demand) should be in generally satisfied during all lifecycle time but it is very important to optimize flows of various energy vectors [19].In reality despite all, power flow unbalance appears.One of the planning goals is to find commercially acceptable system architecture (usually depending on commercial tariffs [20]), component size, system operation rules and control strategies considering uncertainty [21].The reason for power flow unbalance is inherently time unbalance between capacity and demand (that is the crucial problem in power systems based on renewables and microgrids in general).The source of unbalance is of complex nature and generally it is hard to model (prediction of the sources capacity and availability as well as load demand).Prediction of both supply and demand side of the power supply system is a key element in planning phase, [13] and [20].Power load forecasting or demand forecasting is initial and vital process of planning.The forecasting assumes that historical analysis up to some extent can reveal the future behaviour.The process of load forecasting is always being improved since power consumption generally is increasing in magnitude and variety -so one can say that certain load forecasting of several years ago is not valid any more, [18].
Planning of multi-energy system is a more demanding task than planning of single energy systems.A subject of special interest is cogeneration facilities where two useful forms of energy are produced in the same technological process.While production of both forms of energy is coupled, their consumption generally is not!That introduces additional dimension in a planning process: simultaneity of energy products.In [23] a comprehensive two-stage optimal planning for CHP microgrid is presented.The first stage is for determination of optimal type and capacity while the second one is for operation optimization [24].

LOAD FORECASTING OVERVIEW
There are several approaches for load forecasting presented in literature, [1,18,[25][26][27].Load forecasting techniques are divided into two main groups: Statistically based and artificial intelligence based methods.The first group of methods is further divided in: Regression-based methods, Time-series methods, Similar-day approach, and Simulation methods.
Regression based methods (also called econometric) are used for short and mid-term load forecast.The method is based on assumption that the load depends on several external factors.The factors can be weather depending (temperature, moisture, rainfall, wind speed and direction, solar irradiation etc.) and socio-economic factors (population, GDP) [25].
Time series group of methods is dominant and the oldest method used for load forecasting.Those methods are based on presumption that load as data has internal structure which can be decomposed into one or more components (depends on forecasting horizon) such as: trend component, cyclic component, seasonal component and noise (stochastical variation).Various types of filtering are used for getting separate components of load time series [1,25].The simple time series method is ARMA (auto regression moving average) method, which is, based on presumption that fluctuation pattern of load can be assessed from previous measurements.ARMA model is used for short term load forecasting and in the case where time series of load data can be treated as stationary process.The method is further expended on non-stationary process introducing ARIMA methods.
Besides standard, on statistical date based methods, by the development of computational technology and efficient software algorithms, an Artificial Intelligence -based group of forecasting method is appearing.The following methods belong to this group: Artificial neural networks (ANN), Expert systems, Fuzzy logics systems, Support Vector Machines, Ant Colony, and Particle Swarm Optimization.All these methods are generally based on learning process using historical data for estimation of current one [1,[25][26][27][28].

LOAD DURATION CURVE, DEFINITION, APPLICATION, RESTRICTION 3.1 Definition and Presentation of Load Duration Curve
Classical load forecasting approach (above-explained regression or time series) is used for short or middle term rather than for long term planning where chronological sequence of load data is less important than information about volatility across time span.The regression approach is unpractical on long scale, because of lack of available and reliable information about the external factors (temperature, wind…) used in regression formula.The weather depending factors on long scale became a stochastic value and regressions on such external factors have dominant stochastic nature.On the other side, time series approach requires an estimation of model parameters (ARMA, ARIMA, exponential smoothing).By increasing the time horizon for modelling of the load by the time series, the error of estimator (e.g.least square method) increases.Generally, the single set of model parameters should fit the load curve over a quite different period (season).Instead of load date being described by the deterministic expressions (load curves, LC) for long term planning, (especially for generation expansion or reliability evaluation) a probability description of the load via load duration curve LDC is more practical for usage [32,33].By using LDC additional information about how many hours (time) of a given period loads will have at a given value [32] i.e. how the power is used in time (frequency distribution) is available [33].In this approach load is treated as random variable [34] instead of deterministic one (obtained in regression or time series method).This approach is often used in stochastic modelling where besides the randomness of load, unpredictable behaviour of the power source as well as probabilistic nature of discreet events (such as a failure) are properly modelled [16].A practical way of describing the cumulative distribution function (CDF) of the load as random variables (Pt) over long period is to use inverse function of load duration curve [4,[33][34][35].
Load duration curve is defined as a fraction of time, the demand is greater or equal to a certain level [32,33].The easiest way to evaluate LDC is to plot the load demand on a system in decreasing order of magnitude over the period of interest [35,36].The fraction of time for a certain event (load is greater than a certain value) has meaning the probability for such event.Thus, LDC present cumulative frequency distribution of system loads.Sometimes instead of fraction of time, total time is used.In that case, LDC is time function, where time presents cumulative time when system load is greater than given load level.An example of LDC is shown in Fig. 1, [35].Actually, that is inverse function of (1) were time is replaced by value of fraction of total time (time horizon under consideration).The function directly corresponds to cumulative distribution function of the load random variable (2).Relation between alternative presentation of LDC and CDF of load as a stochastic value is presented in Fig. 2. From above elaboration, it can be written: 1 ( ) ( ) 1 Mathematical expression in (1) is more symbolic than practical guide for getting anexact analytic form of the LDC.In several papers, a various estimation (approximation) method for analytic form of LDC are explored.Unlike the load curve, the load duration curve is much easier for approximation.In [21] an example of forecasting of LDC is described.As a first step load time series Xt are forecasted.Instead of modelling X t directly logarithmic transformation of Y t = log(X t ) is modelled as Gaussian structural time series [21].Using such approach a smooth analytic LDC is forecasted by a weighted sum of transformed normal distribution functions.
Similar approach is in [34] where MONA (mixtures of normal approximation) technique is used for approximation of load duration curve.The approximation is based on the segmentation of load range into load classes (k = 1, 2,…, K).For each class parameters of corresponding normal distribution is evaluated.Approximation of LDC by MONA technique is as follows: Where ξ(p) is alternative presentation of LDC, Φ(p, μ k , σ k ) is normal cumulative distribution function, corresponds to fraction of time, μ k and σ k are mean and variance respectively for k class of load.
In [35] the LDC is approximated by Hills function which commonly used in biochemistry and ecology.It was originally put forward by Archibald Hill in 1910 to cover the dissociation curves of hemoglobin.Hill functions is generally representing by (4).The model parameters a, b, c and m in (4) are estimated by Levenberg-Marquardt method.[34] .
There are many other attempts for approximation of LDC by appropriate analytic curve, such like fifth order polynomial function across several times blocks in the year (mostly four) which is also used in WASP program package or combination of exponential function to capture LDC [31].All above presented technics results in mathematically well defend function representing a smooth curve as continues approximation for LDC.
By the rising of computational capability, the number of parameters describing the LDC function became less relevant factor.Normally, LDC can be easily derived out from the series of recoded value of load across the time horizon in a way that load date are reordered by the decreasing value.Since the load is usually recorded on hourly basis, thus both LC and LDC vector contains 8760 date on one year horizon.That huge number of data involve more detail than is merited in planning models.So, LDC are generally approximated over the year by a finite number (much less than 8760) of average loads [32].In such a discrete approximation, a time span is divided into N time segments.Approximation assumes that load value across each time segment is constant and equal to total energy consumed per time segment divided by duration of time segment.In many application of LDC discrete approximation is good enough.A further improvement of discrete approximation of LDC is piecewise linear function where LDC in each time segment is replaced by linear function.The time duration of certain segments should be chosen such that load share of all types of generation (peak, intermediate, base) is well presented [33].Determination of the number of time series and duration of each segments is special challenging task in such approximation.

Usage of Load Duration Curve
Load duration curve directly brings crucial information about how power is used in system, such as total energy consumed, both maximum and minimum load value across a time span.LDC is vital part production simulation, generation expansion (long term planning) reliability evaluation and electricity market analysis [1,36,37].

Figure 3 Simulation by inverse transform method
Since LDC is form of cumulative distribution of load, it can be used for load simulation in the Monte Carlo method in the case where reliability or financial indices cannot be simple derive out using an analytic platforms.The simulation is based on the facts that inverse cumulative distribution function of uniformly distributed random variable is distributed on the same way as random variable which cumulative distribution is subject of interest.The method is called inverse transform sampling (ITS) [38][39][40][41], (5), Fig. 3.
Where u is uniformly distributed random number on interval [0, 1], X t random variable of interest and F Xt distribution function of that random variable.Operator ← denotes "distributed as".
The Method has certain restriction related to inverse transform of distribution function.The issue can be involving LDC (or CDF) approximation (linear, polynomial, MONA or numbers of other approximation technique).
Most common usage of LDC is determination of optimal power technology as an initial stage of expansion planning known as a Screening curve method [12,21,29], Fig. 4. In that phase, the best candidate or plant alternatives are selected as an entry for the second phase of planning where candidates are being compered regarding certain criteria.The method is based on the operating cost function of various power production technology.Each technology is presented via function of total cost regarding capacity factor (or total time of usage) [30].The total cost functions are liner functions which intersections divide the total load area into sections.Each section corresponds to certain power production technology (PT) with the lowest total cost.Besides LDC provides important measures of reliability such as the loss of load probability and expected unserved demand [18].Also, LDC is very often used in power system simulation based on probability where availability of production resources is combined with the load forming equivalent load [1,4,13,16,22].Analysing load duration curve of equivalent load several power system indices can be derived out such as: loss of load probability (LOLP) and expected energy not supplied (EENS).

PLANNING OF COGENERATION SYSTEMS
The focus of this article is on long term load forecasting for expansion generation planning of the microgrids containing multiple load type consumption.Cogeneration is simultaneous generation of two different useful energy forms from single primary energy source.In this article, cogeneration of electricity and heat is considered.In that case cogeneration is also called combined heat and power production (CHP).Unlike utility based power supply, electricity and heat supply from the microgrid cogeneration system is strictly site related process.The exemption are large cogeneration plants for district heating.
Planning of cogeneration system is similar to oneenergy form generation system, but with remarkable difference, certain rate of simultaneity of both production and consumption of both energy components should be taken into consideration.Generally, cogeneration facility always involves generation mix of different ways for electricity or heat supply (or production), Fig. 5.

Figure 5 Diagram of microgrid electricity and heat supply system
Having in mind complexity of microgrids with CHP, integrated resources planning (IRP) should be applied.The aims of IRP are identifying an optimum energy supply mix between a simultaneous electric/thermal production system and utility energy supply [41].Planning of microgrid is more challenging issue than planning other utility connected supply system.Considering both stand-alone (SA) and grid connected (GC) mode of operation of microgrid with CHP and different type of electrical and heat loads that may be employed, it is impossible to maintain simultaneous balance in both form of the energy consumed by using CHP only.Certain level of unmatched load will appear and that may have an impact on operational cost due loss of production or due the cost of alternative supply system (from the utility or the storage).
There are several factors, which influence the microgrid planning.The most important is heat to power ration (HPR).Traditional approach for planning of standalone cogeneration facility is based on HPR as selection criteria (for certain time slot HPR of CHP should match with the HPR of microgrid loads).Typical HPR of certain cogeneration technology is presented in Table 1, [42].The load characteristic can also be quantified via HPR factor.HPR vary between different types of the microgrid application.Typical HPR for industrial microgrid is presented in Tab. 2.
HPR may vary during different times of the day and seasons of the year.To extend the range of application, several engineering methods is used to modify HPR of given cogeneration facility such as: additional firing (increase HPR) or, combining of different technology (gas turbine + steam turbine or ORC to decrease of HPR).Quality of heat (temperature and pressure of transferring media) is another important factor which strongly influence to microgrid planning.Unlike the electricity, which can be easily and efficiently transforms across the wide range of use (AC, DC, different voltage and frequency) available heat (e.g.waste heat of engine) cannot so easily being used for given heat load (free heat transfer goes from higher to lower temperature only).
The load matching principle is also important influencing factor.Depending on the type of the load and application, the planning may be based on electric load or heat load matching.In case of bottoming cycle, CHP is used to meet local need for heat at first.Deviation in electric power supply is settled via utility or other back up power supply or storage.In case of topping cycle the loading, priority is inversed and deviation in heat supply is settled via additional boiler or heat storage.

PLANNING OF MICROGRID CONTACTING CHP BASED ON PROBABILISTIC APPROACH
Planning based on HPR criteria is simple and fast for rough analysis of application certain cogeneration technology in given microgrid, which contain both electric and thermal load, especially in SA mode of operation.The method assumes that HPR of the load does not change or have small variation in each time slot or correlation between heat and electricity is strong.In reality, the microgrid contain several types of electrical or heat devices or devices which consume both energy form in same time where the heat load does not follow electricity load so strictly (HPR vary trough the time in wide range) or does not have a known pattern.The variety of the load, as well as the availability of the production unit on macro scale is further emphasized on micro scale -microgrid.So the planning of the microgrid containing both heat and electricity load cannot be based on load duration curves for heat or/and electricity only.Figs. 6 and 7 present two examples (I and II) of time series of heat and electricity respectively.For the sake of illustration, grid connected mode of microgrid is assumed where the deviation in electricity supply is settled via utility and for heat via additional heat source.In addition, it is assumed the loading follow heat priority principle.Fig. 7 presents the same load profile for heat and the different one for electricity load presented in Fig. 7.In both cases, the load duration curves for both heat and electricity demand are the same.
As can be concluded from the Fig. 6 and Fig. 7 the same load duration curves and HPR index on given time slot, do not produce the same performance of certain CHP.To capture this occurrence, the load demand and unit's capacity should be compered in the same space and manner where the time interdependence of load and CHP capacity is properly treated.As it was emphasized in chapter 3.1.,for long term analysis of microgrid, stochastic simulation and analysis is more appropriate then deterministic approach as it was proposed in [13-17, 22, 29, 32, 38, 44-51].Once developed stochastic models for loads, power supply system various reliability analysis, sizing, optimization of operation and control of microgrid can be performed using a powerful tool such as Monte Carlo Method [14,40,41].In following subchapters, a new approach for stochastic modelling of both load and power production as well as its usage in stochastic analysis is described.

Stochastic Model of Two-Dimensional Load
From the reason explained in chapter 3.1.,stochastic model of load demand vector of the certain microgrid LE,H = (e, h) is more appropriate than time series format.In that case L E,H is two-dimensional random vector, presented by its load probability density function LPDF E,H .Many authors, like [45][46][47][48][49][50][51], used this approach.The meaning of LPDF E,H = f E,H (e, h) of random energy vector L E,H (heat, electricity) is defined by (6).

E ,H f e,h e h P E e,e e H h,h h
Based on given LPDF E,H and outage rate of individual energy sources, the equivalent load probability density function ELPDF E,H is calculated [45,[47][48][49][50].The ELPDF's present alternative load probability functions where random outage of the source unit has been considered.Once LPDF's and/or ELPDF's are known in analytic form, a various availability indicators such as Loss of Load Probability (LOLP), Expected unserved energy (EUE) and others [1,6,45,49,50] can be calculated by analytic integration.
The analytic form of LPDFE,H plays crucial role for stochastic analysis CHP in general.Special challenge is appropriate approximation of LPDF from raw date to get appropriate analytic and ease to integrate form.
The LPDF gotten from recorded data is rather "spikey" than smooth curve so large errors in approximation can be expected.Furthermore, integration of such function to get certain expected figure of merits became quite difficult.
Instead of LPDFE,H cumulative distribution function of two dimensional energy vector CDF E,H is proposed as basis of stochastic analysis.Moreover, this approach is based on two graph-analytic analogies with power only system) LDC to CDF analogy explained in the chapter 3. 2 and b) the analogy of the value of the surface below the CDF and value of expected value energy used.
Since multi-dimensional nature of the energy vector under consideration, it is impractical to use a term"curve" so the term "function of simultaneous heat and electricity duration" (FSDE,H) is proposed.
FSD E,H (e, h) is two dimensional function of electricity and heat with the meaning of total number of hours in corresponding time range where electricity load demand is greater than given value e and heat load demand is greater than given value h.Note that one-dimensional load duration curve (LDC) mostly used in form offunction of time p = p D (τ), while FSD E,H is alternative presentation as function of energy vector (7).

T E ,H t T e,h t,e,h t e,h LR
Where The function of FSD E,H presents stochastic description of given two-dimensional load demand on certain location.Obviously, such characteristic can be presented via surface in E-H space.FSD E,H can be derived out from the recorded data and appropriate approximation.

Determination and Approximation of FSD E,H
The procedure goes in several steps as follows: First step is collecting a data regarding heat and electricity on given location during a certain time span, e.g. a year.For the sake of long-term analysis, at least hourly average value for electricity and heat is recommended.Very often, the load data for electricity and heat was count separately.Thus, a certain form of data synchronization is needed.The raw data should be prepared and written in a form of load matrix L3×K = {(t k ,e k ,h k ); k = 1, 2,…, K}, where t k , are corresponding time interval of length ΔT(h), e k and h k are average load of electricity and heat in that time interval.All data can be also presented as a dots in set of recorded data S r = {(e k , h k ), k = 1, 2,…, K} in E-H space, Fig. 8.
Second step.The load range LR is divided into load subranges (LR p,q ): ) .
Where p and q are integers in range (1, N), N is the fineness of approximation.Note that N<< .K Further dE and dH are load increments: dE = (E max − E min )/N; dH = (H max − H min )/N.
Third step.For each p and q to determine subset S p,q which contain a records (e k , h k ), which fulfil the following conditions: 1, 2,..., : Where E p = E min + pdE and H q = H min + qdH presents starting coordinates of LR p,q subrange, Fig. 8.
Forth step.Determination of the number of elements C p,q in each subset S p,q ; C p,q = card(S p,q ) and formationof a . The values in p row and q column of the matrix present a total time of simultaneous duration of electricity and heat load above E p and H q value respectively.By appropriate selection of number N according to total recorded data, every load vector (E,H) can beapproximated with (E p , H q ), so it can be concluded that T p,q is approximation of FSD E,H in E-H space.For the sake of demonstration of the new approach, the load data for electricity and heat consumption in one industrial microgrid in time range July 2015 up to May 2016 are collected and used for evaluation function of simultaneous heat and electricity duration, Fig. 9. Technical Gazette 24, 5(2017), 201-212

Model of CHP device in E-H space
In CHP device, the rate of produced heat and electricity are strongly correlated.The correlation function can be presented by following form of HPR expression.

amb ( , , , , ). h e HPR f T y t Constraints = ⋅ (10)
Where h and e are value of produced heat and produced electricity respectively.Generally, HPR depends on several parameters.Depending on the type of CHP units, the HPR is independent from a rate of partial load f in most of range of its value.Partial load f corresponds to level of fuel flow.At low partial load values the HPR function swings toward to certain point (0, H r ) Fig. 10a.That is because of very low electricity efficiency at small partial load level.In most cases, operation of CHP unit below certain partial level is not recommended.
Ambient temperature T amb has a non-neglecting influence on CHP performance.There are common rules that electric efficiency decrease and thermal efficiency increase by the increase of air temperature of CHP unit.Consequently, HPR increase with ambient temperature.For the sake of demonstration, in this article the temperature dependency will be neglected.
Besides fuel control, the values of produced heat can be controlled by flow control of heating inter-media (like in extraction steam turbines based CHP) or by control of the flow of media for heat removal (like in gas engine or gas turbine based CHP).In some cases, electricity is also affected by this way.Such control is represented by nondimensional value y.If there are more than one control options, y became a vector, Fig. 10b.
The performance of CHP also depends on the state of thermal engine: cold or hot state.That is represented by the time in operation t from cold condition.Running in cold state is operation under poor efficiency and high operational cost unless of running in hot state.For the sake of simplicity and demonstration of the CHP unit modelling this feature can be neglected in long term planning.
There are several constrains to be taken into consideration when CHP unit is under planning for certain microgrid.First of all, technical constrains are related to prevention severe damage caused by overheating either thermal units (turbines, engine, pipes and other thermal installation) or generators (overcurrent).Technical constrains are represented by Emax and H max .Besides there are economical constrains which are related to minimal partial load at which the operational costs are too high (poor efficiency).At the end, there are legal constrains which are related to the limits defined by legal authority such as limitation of emission or limitation for liable producer (e.g.primary energy saving rate > 10%).
Depending on variables and constraints explained above, the correlation of the produced heat and produced electricity can be presented as an area of points (h, e) in E-H space, where each point presents allowable working condition, which can be reached by certain control.This area is called working area (WA) or area of feasibility operation [45,47,48,[52][53][54][55], Fig. 11.
The shape of WA can vary from straight line for simple CHP unit up to complex non-convex area.The working area is used to established limitation for e and h at various optimizations algorithm (11), (12), [49].max min ( ) ( ).
The function min ( ), chp H e in (11) and ( 12) are determined by the boundary of working area, Fig. 12a.In such approach, working area is treated as a set of points (h,e) in E-H space bordered by several boundary lines.This approach is preferable at short-term optimizations and economic dispatching.The alternative approach, more suitable for long term CHP plaining in combination with FSD E,H , is to treat the WA as a family of line h=e•HPR(y), Fig. 11b.

USAGE OF FUNCTION OF SIMULTANEOUS HEAT AND ELECTRICITY DURATION
In Chapter 5.1 the presentation of the multiple load for given period through the duration function giving the information about amplitude and frequency of occurrence of certain load level was described.Such presentation, FSDE,H function, can be used for sampling the stochastic value of electricity and heat load.Following the analogy with the power-only load duration function and sampling in chapter 3.2, the inverse transform sampling (ITS) method can be applied on FSD E,H too.Generally, such approach assumes analytic form of two-dimensional distribution and its inverse function what can be demanded task even for power-only system.
In case of two-dimensional distribution presented by FSD E,H , the value of the distribution are only known in the nodes of the mesh created by discrete value of E p and H q , where p,q = 1, 2,…, N.
The creation of the analytic form of the FSD E,H means determination continues two dimensional function which best describes (interpolates) the value of FSD E,H in whole load range based on the discreet value in mesh of nodes S p,q in (9).For similar purpose, in case of determination of the distribution function, the interpolation by some wellknown distribution is very often used [17,18,31,34].For the sake of robust and fast stochastic analysis of the microgrid the piecewise linearization is suggested.
The value of the FSD E,H in arbitrary point T(e, h) can determined by known values in subsequent nodes of S p,q as follows:

E e T T e,h e E T E
In case of e p < e < e p+1 ; h In case of e = e 1 ; In case of E p < e < E p+1 : H q < h < H q+1 ; 1 < p < N; 1 < q < N; where ; p,q p,q p,q p ,q p ,q T T D T T In case of e = E N or h = H N .Note that (13) presents marginal distribution for electricity with alternative notation T E (e) = T(e, h = h 1 ) while (14) presents marginal distribution for heat T H (h) = T(e = e 1 , h).

Heat and Electricity Sampling
Since the function of simultaneous heat and electricity duration FSD E,H has the meaning of heat and electricity distribution where the interdependency between these two system variables are properly described, it can be used for heat and electricity sampling.Such obtained pseudo random system variables (load) is input set for stochastic simulation.The procedure of heat and electricity sampling using FSD E,H and its numerical approximation is shown in Fig. 12.
3. Sampling of heat load according to marginal distribution of heat distribution: Where q satisfy inequality T q < t 1 < T q+1 .

Determination of conditional electric load distribution: sampl ( )
E|h T e , using (15) with substitution h = h sampl ; 5. Determination of inverse function 1 sampl ( ) Where ˆˆˆp,q q p,q q a T H h T h H 6. Scaling of uniformly distributed random variables for conditional electric distribution.
Alternative procedure is similar to above explained with the difference that electric load is sampled at first and then heat load.

Calculation of Total Heat and Electricity Energy Produced in Cogeneration Process
Besides of the application for heat and electricity load sampling, the approximation of FSDE,H is used for determination of total heat and electricity energy produced in certain cogeneration facility.The procedure is based on combining the working area of certain cogeneration facility and FSD E,H in the same E-H space, Fig. 13.Let assume that working area is presented as a set of family of line h = e•HPR(y), Fig. 11, where y is heat control vector (e.g.: rate of steam extraction in case of cogeneration based on extraction steam turbine).There is certain trade-off between operational cost and rate of heat control.In most cases the least cost per produced energy is at y = 1.Other heat control, if are required often, leads to higher cost and lower efficiency).
If y is assumed as fixed value, heat and electricity production are generally linked by linear HPR correlation: are limitation of CHP unit.In Fig. 14 the FSD E,H of certain load and Working area/line are presented.When the working area is projected on the FSD E,H surface, certain curve K C is obtained.That curve represents the rate of heat and electricity load to be satisfied by single cogeneration unit.Index C means that the rate of simultaneous heat and electricity load feed by single cogeneration unit depends on the way of heat control.
The practical meaning of the applied procedure is when curve Ψ C is mapped on the heat and electrical marginal distribution plane T E (e) and T H (h) respectively.The curves Ψ C|H and Ψ C|E presents distribution function of for heat and electricity produced in certain cogeneration unit (system) and used in certain consumption area (microgrid).The surface below the curves represents the value of both heat and electricity energy supplied by the unit.
The curve Ψ C can be determined as orthogonal mapping of the heat to electricity curve (HPR) of cogeneration unit (22) on the surface of FSD E,H (13)-( 16 Where T(e, h) is ( 15) rewritten in canonical form with following coefficients: Also, matrix D p,q are according to (15).Since T(e, h) is defined on the whole area for each combination of indexes p and q, the ψ C depends on p and q which are now functions of the value of electricity power or heat.To avoid confusion, they will be written as p' and q': here INT is function which rounds an argument down to the nearest integer.By ( 29) and (30) and for a certain CHP unit a particular Unit -Load subrange (ULS) as a subset of Load Range LR is created (yellow rectangle in Fig. 15).
Finally, the amount of electricity produced by certain cogeneration unit and used in consumption area on yearly basis can be determined by integration of Where, e f and e b are the value of the electricity on the edge of Unit Load Subrange ULR p',q' .In similar way, the amount of heat produced can be determined (by the integration of

Calculation of CHP Potential in Grid Connected and Stand-Alone Mode of Operation of Microgrid
Without going in deeper elaboration, a method presented in 6.2 can be used for evaluation of efficiency of usage of certain CHP.Previously described method, brings the capacity of certain CHP unit in simultaneously production of electricity and heat in stand-alone mode without any storage capacity on the site.
In case of grid connected mode of Microgrid operation or in stand-alone with some energy storage, a certain amount of heat or electricity can be taken over or delivered aside of local consumption.
The procedure starts with mutually mapping of marginal distributions TE(e) and T H (h) on the opposite planes, Fig. 15: ( ) ( ) Where h(e) is HPR correlation in general form.Four characteristic area are obtained: Area A is the amount of energy simultaneously produced and used in local area (stand-alone, no storage capacity of CHP, (39)).
Area B is the amount of energy, which cannot be produced in CHP unit because of power of load, exceeds maximum capacity of CHP unit or it is below the technical minimum of CHP unit.That amount should be purchased in any other way or it presents amount of the loss.
Area C corresponds to potential additional amount of heat or electricity produced above the amount of A in case the CHP production follows electricity or heat load respectively.That amount presents the surpluses which need to be delivered to utility or storage.
Area D presents the amount of heat or electricity missing for direct load of local consumption in island mode of operation.

CONCLUSION
Planning is initial phase of each life cycle of certain technical system.It determines main figures such as system structure and operation strategy.Small power system containing multi-energy load and as well distributive generation operable in both grid connected and stand-alone mode -microgrid, is special concern of planning.The quality of various planning method in most cases depends on the quality input data.Load forecasting is vital part of each planning analysis, especially where stochastically nature of load is emphasized (what is the case in microgrids).Besides of load curves, which is often used as a form of load representation, a load duration curve is also used when stochastic nature of the load need to be captured.By the using of LDC instead of LC some additional information, such as frequency distribution of the load, are available at planning.Planning of the microgrid containing cogeneration facility is demanded task since both mode of operation and emphasized unpredictability of the load should be taken into consideration.Various stochastic model has been developed for such purpose.The article derived out a new approach in stochastic modelling of microgrid, which is equipped with cogeneration facility.The approach is based on the function of simultaneous heat and electricity duration given by bilinear interpolation of load date and appropriate presentation of the two-dimensional characteristic of the cogeneration units in the same space.Such approach can be used for numerical evaluation of the capacity of certain cogeneration unit in given load conditions on long-term basis.Besides, heat and electricity sampling can be easily derived out using such twodimensional distribution function.
In further research, the described approach can be used for various efficiency evaluation and optimization methods of CHP unit used in microgrid.

Figure 1
Figure 1 Load curve and Load duration curve [35]If p L (t) is load curve function, the load duration curve p D (t) can be mathematically expressed as follows:

Figure 2 a
Figure 2 a) and b) Alternative presentation of LDC, c) CDF of P t

Figure 4
Figure 4 Screening Curve Method

Figure 6 Figure 7
Figure 6 An example I of time series of heat and electricity

Figure 8
Figure 8 Recorded data in load range, E-H space

Figure 9
Figure 9 An approximation FSD E,H

Figure 10 Figure 11
Figure 10 a) Types of Working area a) single curve in case without of flow control heating media; b) in case of flow control of heating media

Figure 12
Figure 12 Heat and Electricity sampling by FSD E,HSampling procedure: 1.Two independent, uniformly distributed, random variables u 1 and u 2 was sampled; 2. Scaling of uniformly distributed random variable for heat distribution by using u 1

Figure 13
Figure 13 Heat and Electricity sampling by FSD E,H

Figure 14
Figure 14 Overlap of CHP Unit Working area and Load Range

Figure 15
Figure 15 CHP potential in Microgird in Stand-alone mode Finally, B + D is gross value of directly not-supplied energy from the CHP unit to the local consumption and B + D − C is net shortage or surplus in grid connected mode or stand-alone mode with hypothetical infinitive storage.

Table 1
HPR of different cogeneration technology