Scenario-based stochastic optimal operation of wind/PV/FC/CHP/boiler/tidal/energy storage system considering DR programs and uncertainties

Abstract

Background

Micro-grid (MG) can be described as a group of controllable loads and distributed energy resources that can be connected and disconnected from the main grid and utilized in grid-connected or islanded modes considering certain electrical constraints.

Methods

The objective of this article are as follows: (1) predict the uncertainties through the hybrid method of WT-ANN-ICA and (2) determine the optimal generation strategy of a MG containing wind farms (WFs), photovoltaic (PV), fuel cell (FC), combined heat and power (CHP) units, tidal steam turbine (TST), and also boiler and energy storage devices (ESDs). The uncertainties include wind speed, tidal steam speed, photovoltaic power generation (PVPG), market price, power, and thermal load demand. For modeling uncertainties, an effort has been made to predict uncertainties through the hybrid method of wavelet transform (WT) in order to reduce fluctuations in the historical input data. An improved artificial neural network (ANN) based on the nonlinear structure is applied for better training and learning. Furthermore, the imperialist competitive algorithm (ICA) is applied to find the best weights and biases for minimizing the mean square error of predictions.

Result

The scenario-based stochastic optimization problem is proposed to determine the optimal points for the energy resources generation and to maximize the expected profit considering demand response (DR) programs and uncertainties.

Conclusions

In this study, three cases are assessed to confirm the performance of the proposed method. In the first case study programming, MG is isolated from grid. In the second case study, which is grid-connected mode, the WT-ANN-ICA and WT-ANN uncertainty prediction methods are compared. In the third case, which is grid-connected mode, the effect of DR programs on the expected profit of energy resources is assessed.

Background

The micro-grid (MG) concept has recently attracted significant public attention. Integration of renewable sources, combined heat and power (CHP) systems, and energy storage technologies in the MGs will result in environmental friendly, low cost, and reliable energy. Recently, using CHP systems in MGs has attracted more attention. The primary motivation for incorporating CHP units is providing electrical and thermal energy, simultaneously. During electricity generation process of CHP systems, waste heat is employed to provide thermal energy. This process will result in the improvement of overall system efficiency as well as a significant reduction in the cost of thermal energy generation. It should be mentioned that, in a CHP unit, the power generation boundaries depend upon the heat generation of unit and the heat generation boundaries depend on the power generation of the unit [1].

Owners of renewable resources need to predict the uncertainties for optimal planning such as photovoltaic voltage/power generation [2], market price [3], and load forecasting [4], wind farm power generation/wind speed (WS) [5,6,7,8,9]. In [7], firstly, historical data of WF is decomposed using wavelet transform (WT) and then WF power generation is predicted by artificial neural network (ANN). This method is tested in two regions of china. Afterwards, comparing WT-ANN, ANN, and auto regressive moving average (ARMA) methods revealed that WT-ANN can significantly reduce the error in spite of ANN and ARMA methods. In [8], the optimal weights and biases of ANN are determined by genetic algorithm (GA), imperialist competitive algorithm (ICA), and ICA-GA methods; then, they are tested on six specified databases. In the end, the obtained results confirmed that ICA has higher capabilities. Similarly, ANN is employed to predict WF power generation, and then, ICA, GA, and particle swarm optimization (PSO) are chosen to determine the optimal weights and biases [9]. The prediction results were more satisfactory when ICA algorithm was utilized.

The second solution for uncertainty reduction in renewable units including renewable resources is to coordinate other energy resources which are quite expensive, but available and more reliable, such as pump-storage unit, hydro unit, gas turbines, combined cycle power plants, and energy storage batteries [10,11,12,13,14,15,16,17,18,19,20,21]. However, the share of these energy sources should diminish for many reasons [10]. In [11], the coordinated planning of WF, pump-storage unit, and thermal units is presented by the multi-stage stochastic planning and solved by scenario decreasing algorithm of PSO. In [12], the required reserve level is estimated in the presence of high-level WF penetration. In [13], the optimal strategy of WF is determined in the real-time market. The wind speed and market price are predicted by ARMA. Moreover, the expected profit is limited by FR and the required reserve is determined due to the error prediction in WF power generation. In [14], the coordinated planning problem of WF and thermal power plants are solved by artificial immune optimization method. This optimization method is implemented on a system including ten thermal power plants and two wind farms (WFs). A mixed integer programming algorithm is adopted for period planning of operation startup/shutdown and generating/pumping mode of pump-storage unit to maximize the profit in coordinated operation of WF and pump-storage unit [15]. A scenario-based and chance constrained optimization method is hired to consider the WF power generation prediction error. A rolling optimization method for WF coordination with the energy-storage systems in the day-ahead market is presented to increase the profit of these power plants.

The optimal scenario-based operation management of MG including WF, photovoltaic, micro-turbine/fuel cell, and energy storage devices are studied in [16]. In this paper, the considered uncertainties are load, WF power generation, photovoltaic power generation, and market price. In [17], the optimal biding strategy model in an electricity distributed company is considered in order to make the maximum profit in the day-ahead market. In [18], the modified particle swarm optimization algorithm is used to optimize energy in MG. Moreover, in this study, uncertainty of data is checked using Hong method. In [19], like [18], Hong method is applied for covering uncertainties; however, the modified firefly algorithm is utilized for optimization. In [20], studies on utilization of micro-network are made in the presence of generating resources of thermal and electrical energy and also Proton Exchange Membrane Fuel cell power plant along with the hydrogen storage. The modified algorithm of self-adaptive charge search algorithm is applied for optimization. In [21], the objective function is considered to maximize the profit of wind farm, fuel cell, boiler, CHP units, electrical power generation unit, and energy storage devices (ESDs) connecting to a MG regarding uncertainties. The uncertainties are predicted by time series methods.

In this paper, the presented issue can be shortly explained as follows:

• 1. Prediction of uncertainties via hybrid method (HM) of WT-ANN-ICA. According to the studies in [7,8,9], prediction of uncertainties using the proposed method can lessen errors of prediction of WS in comparison to ARMA, ANN, WT-ANN, WT-ANN-PSO, and WT-ANN-GA methods. Therefore, this approach may generate scenarios closer to reality and lead to the optimal programming.

• 2. Generating the scenarios of WS, tidal steam speed (TSS), photovoltaic power generation (PVPG), market price, and power/thermal load demand, decreasing the scenarios with the scenario-reduction backward method, and modeling them through the tree scenario method.

• 3. The programming of MG, including WFs, photovoltaic (PV), tidal steam turbine (TST), fuel cell (FC), CHP units, boiler, and electrical and thermal ESDs, considering constraints and the uncertainties of WS, TSS, PVG, market price, and power/thermal load demand.

• 4. Studying the expected profit of energy resources with and without DR program.

Methods

An algorithm is proposed for programming generation and unit commitment of an MG including three WFs, PV, TST, FC, two CHP units, boiler, and ESDs with and without considering DR program shown in Fig. 1.

Fig. 1

The flowchart of the proposed method

Scenario-based stochastic modeling

As a result of extending renewable resources and uncertainty in the nature of such resources, the modern complicated power systems should be analyzed in uncertain environment so that operating point and reliability of energy supply occur in approximation with the optimal point in reality. Therefore, having access to powerful tools is necessary for transition from uncertain environments with random variables, including their probability contributions, to the certain problems with certain variables. In the modern deregulated power supply markets, the most important random variables are load demands, wind speed, PVPG, and market price. The origin of the abovementioned uncertainties are found in issues like weather conditions, temperature variations, and government decision.

This proposed method for prediction of uncertainties is proposed in Fig. 1 as a flowchart. As observe here, it is assumed that the prediction for dth day can be made and the historical data extracted for every single hour of 24 h beginning 100 days ago.

Stage 1: data homogenization

The historical data are recalled and normalized to improve data homogenization.

Stage 2: data processing using wavelet theory

The components and features of data can be extracted through mathematical equations. More specifically, the components and features of time and frequency domain of data signal can be extracted through wavelet technique. The basic equations of WT are Eqs. (1 and 2) [22].

$$ \mathrm{WT}\left(a,b\right)=\frac{1}{\sqrt{a}}\underset{-\infty }{\overset{+\infty }{\int }}f(t)\psi \left(\frac{t-b}{a}\right)\mathrm{dt}\kern1em a={2}^{-j},b=k{2}^{-j}\in R,a\ne 0 $$

(1)

$$ f(t)=\frac{1}{\underset{2}{c\uppsi}}\int \underset{\int }{}{\uppsi}_{\psi}\frac{1}{a^2}\psi \left(\frac{t-b}{a}\right)\mathrm{dadb} $$

(2)

where ψ (a, b) is the wavelet function and f (t) is the input signal on which wavelet function is implemented until the WT (a, b) signal is yield. Furthermore, a and b are the parameters related to the WT which depend on the type of wavelet function. The approximated values are decomposed once more after some iterations; therefore, the signal is decomposed into smaller parts [7, 23]. WT is essential here to suppress the disturbances in historical data and to alleviate the fluctuation of input data. The input data is decomposed into three approximated components of (Dh1, Dh2, and Dh3) with lower accuracy together with a more precise component (Ah) which plays the most important role in prediction process [23].

Stage 3: ANN

McCulloch and Pitts tried to simulate the ANN by a logical model for the first time which now is widely applied in many fields. Here, the chosen ANN consist of three perception layers: the output layer with one neuron, the input layer with five neurons, and the hidden layer with three neurons. This ANN can predict the information of hours d (t + 1,…,t + 24) for the output signals of WT as the initial data.

Stage 4: ICA

ICA is a new optimization strategy based on political and social evolution of human. Basically, to determine the best solution, GA and PSO are inspired with biological evolutions, chromosomes, and particles. However, the source of inspiration in ICA is the social-political evolution, and it applies colonies (countries) as the variable for finding the optimal solution [8, 9]. The steps of ICA are briefed as follows:

• 1. Developing initial colonies: the ANN consist of input signals (Ah, Dh1, Dh2, and Dh3), the five neurons in the input layer (IL), the three neurons in hidden layer (HL), and the one neuron in output layer (OL). The matrixes of wrights (W) and biases (B) consist of ILW = [5 × 4], ILB = [5 × 1], HLW = [5 × 3], HLB = [3 × 1], OLW = [3 × 1], and OLB = [1 × 1]. Hence, each colony constitutes 47 variables. Initial colonies are selected through specific range based on initial training of ANN on a random basis. Regarding the cost function based on decreasing the prediction error, the optimization of weights and biases are performed within the neural network for better training. The cost function here is the mean square error which is proposed as Eq. (3).

$$ \min\;\mathrm{FunctionCost}\kern1em \mathrm{MSE}=\frac{1}{\mathrm{IN}}\sum \limits_{\mathrm{in}=1}^{\mathrm{IN}}{\left|{\widehat{y}}_{\mathrm{in}}-{y}_{\mathrm{in}}\right|}^2 $$

(3)

where ${\widehat{y}}_{\mathrm{in}}$ and y_in are predicted and real value for inth input.

• 2. Selecting the imperialist: in this stage, the colonies with minimum cost are selected as the imperialists.

• 3. Allocating the other countries as the colony to the imperialists: in this step, some colonies are allocated to each of imperialists and empires. This allocation is done according to imperialists fitness (fewer cost) by stochastic universal sampling method. The stages of 1–3 are the initialization stages of ICA.

• 4. Performing the act of assimilation or absorption policy: in this stage, each of the colonies is moved towards the imperialist in each empire. This stage proceeds to improve the exploitation of algorithm.

• 5. Performing the act of revolution: in this stage, the random changes are applied on each of the colonies. This action can improve the exploration of algorithm and prevent from involving the optimization in the local optimal points.

• 6. Computing the cost of colonies and imperialists.

• 7. Comparing the cost of colonies with imperialist in each empire: if a colony holds a lower cost than the imperialist, it will take its place.

• 8. Evaluating the empires: the cost for each empire is computed according to Eq. (4).

$$ {\mathrm{Cost}}_{\mathrm{empire}}={\mathrm{Cost}}_{\mathrm{imperialist}}+\frac{0.1}{N_{\mathrm{COL}}}\sum \limits_{n=1}^{N_{\mathrm{COL}}}\left({\mathrm{Cost}}_n\right) $$

(4)

where N_CLO is the number of colonies.

• 9. Decreasing the colonies: in this stage, a colony is omitted from the weakest empire and transmitted to another empire by roulette wheel method. According to this method, the empire with the lower cost has more chance to seize the colony.

• 10. Omitting the empire: if the weakest empire has no colony, the related imperialist will be transmitted to another empire as a colony.

Stage 5

Studying the termination condition, the stop condition is set based on the number of iterations obtained by trial and error method. If the stop condition of program is satisfactory, the results are moved to the scenario generation stage; otherwise, the algorithm returns to (4) to generate new colonies. ICA flowchart is illustrated in Fig. 2. The uncertainty prediction curves are shown in Fig. 3a–f. The comparison of mean square error of proposed method with ARMA, WT-ANN, and WT-ANN-PSO methods are listed in Table 1.

Fig. 2

Flowchart of ICA

Fig. 3

The uncertainty prediction curves: a Electrical demand. b Energy price. c PV power generation. d Thermal demand. e Wind speed. f Tidal steam speed

Table 1 Mean square error of four prediction method of WS

Generating scenarios and backward method scenarios reduction

According to the stated issues, the determination of optimal strategy for resources connected to the MG is analyzed randomly. To reach this goal, at first, a probability density function is defined for each variable. In this study, the applied probability density function is adopted for power/thermal load demand, TSS, PVPG, and market price with normal distributed functions. In the case of WS, the statistical model is not coordinated with the normal distribution but more harmonized with Weibull distribution function.

To accomplish this, first, a probability density function is defined for each variable. In this study, the applied probability density function is adopted for power/thermal load demand, TSS, PVPG, and market price with normal distributed functions. In case of WS, the statistical model is not coordinated with normal distribution but more harmonized with Weibull distribution function. First, this distribution function is portioned into N parts with the mean of zero from the center with the width of α. Next, each portion is allocated to the occurrence probability and specific error percentage at each level, Fig. 4, [16]. Then, the probability of each occurrence is normalized in a sense that their accumulated distribution function is equal to 1. Finally, a number is selected for each uncertainty variable and each time interval by roulette wheel method in a random matter; hence, the intended scenario is yield.

Fig. 4

The probability density function of uncertainties

The rate of each scenario is obtained by the sum of the error and predicted amount of variable [16]. Eq. (5) shows the amount of scenario for the WS. Consequently, 500 scenarios are generated for each uncertainty.

For modeling, all uncertainty parameters, including WS, TSS, PVG, market price, and power/thermal load demand, and many scenarios are generated. However, the huge number of scenarios makes it burdensome to solve the stochastic problem. In order to solve this problem, the number of scenarios should be declined by the backward method. The basis of this method is to merge the scenarios with close probability into one. This process could continue until reaching the favorable numbers [16, 21]. In this research, the number of scenarios abates down to 10 for each state.

$$ {P}_G^W\left(w,s,t\right)={P}_{G\mathrm{forecasted}}^W+\Delta {P}_G^W\left(w,s,t\right)\kern1em t=1,.\dots, 24,\kern1em s=1,\dots \dots, {s}_w,\kern1em w=1,\dots \dots, {W}_N $$

(5)

Objective function of WFs, PV, TST, FC, CHP units, boiler, and ESDs

In this study, the optimal scheduling of MG including WFs, PV, TST, FC, CHP units, boiler, and ESDs is examined with the 24-hour time horizon as well as considering uncertainties and DR programs in order to maximize the expected profit. The multi-stage stochastic programming is applied to deal with uncertainties. Since the generation power of units should be determined before applying stochastic processes, they are the first stages or here-and-now decisions and are not dependent to the scenarios. Other variables such as buy or sell power from the market and charge or discharge of storage devices are at the second stage or wait-and-see decisions. This mixed integer nonlinear optimization problem is solved through GAMS/COUENNE software. GAMS/COUENNE is a GAMS solver that allows users to combine the high level modeling capabilities of GAMS with the power of COUENNE optimizers which are designed to solve large and difficult problems quickly and with minimal user intervention. This solver is a general one which can be used to solve all scheduling problems. In these conditions, COUENNE solver can be used to solve the proposed optimization problem. COUENNE uses a branch and cut algorithm which solves a series of linear programming and sub-problems [24].

Problem modeling

In this section, an optimal bidding strategy is modeled and analyzed. The objective function of this optimization problem utilized for the first time is as Eq. (6). where, ρ _s is the probability of scenario s. P_sale(s, t) and P_buy(s, t) are the amount of power sold and bought to/from the market, respectively. According to Eq. (7), the probability of sth scenario is obtained by multiplying the probabilities of WS, TSS, PVPG, market price, power, and thermal load demand in each other. The function C _T is the total operation cost of units defined in Eq. (9). The objective is to maximize the expected profit of units considering constraints related to unit usages; the difference between selling excess energy to the market in grid-connected mode and costs. The costs include buying energy from the market in grid-connected mode, the expenditure of operation, startup and shutdown cost, the cost of charge of electrical ESD, and batteries of PV resource proposed as Eqs. (9–12). Eq. (15) is the power balancing constraint of MG.

$$ \operatorname{MAX}\kern0.5em {\mathrm{ER}}_T=\sum \limits_{t=1}^T\left[\sum \limits_{s\in S}{\rho}_s\left({P}_{\mathrm{sale}}\left(s,t\right).{E}_P\Big({s}_p,t\left)-{P}_{\mathrm{buy}}\right(s,t\left).{E}_P\right({s}_P,t\left)-{C}_T\right(s,W, PV, TST, FC, CHP,K,B,t\Big)\right)\right] $$

(6)

$$ {\rho}_s={p}_P\times {\rho}_W\times {\rho}_{TSS}\times {\rho}_{PV}\times {\rho}_{PL}\times {\rho}_{HL} $$

(7)

$$ s={s}_P,{s}_W,{s}_{PV},{s}_{PL},{s}_{HL},{s}_{TSS} $$

(8)

$$ {\displaystyle \begin{array}{l}{C}_T\left(s,W, PV, TST, FC, CHP,K,B,t\right)=\\ {}\left(\sum \limits_{W=1}^{W_N}{A}_W.M\left(W,t\right)+\sum \limits_{CHP=1}^{CHP_N}\begin{array}{l}\left({C}_{CHP}\left(P,H\right).M\left( CHP,t\right)\right)+\left({A}_{PV}.+{BATT}_{\mathrm{COST}}\left( PV,t\right)\right).M\left( PV,t\right)+\left({A}_{FC}+{B}_{FC}.{P}_G^{FC}\left(s,t\right)\right).M\left( FC,t\right)\\ {}+{BATT}_{\mathrm{COST}}\left(K,t\right).M\left(K,t\right)+{A}_B.{P}_G^B\left(s,t\right).M\left(B,t\right)+{A}_{TST}.M\left( TST,t\right)+\sum \limits_{i\in CHP, FC,B}\left\{{U}_{\mathrm{COST}}\left(i,t\right). SU\Big(i,t\left)+{D}_{\mathrm{COST}}\right(i,t\left). SD\right(i,t\Big)\right\}\end{array}\right)\end{array}} $$

(9)

$$ {C}_{CHP}\left(P,H\right)={A}_{CHP}+{B}_{CHP}.{P}_{G, CHP}(t)+{C}_{CHP}.{P}_{G, CHP}^2(t)+{D}_{CHP}.{H}_{G, CHP}^2(t)+{E}_{CHP}.{H}_{G, CHP}(t)+{F}_{CHP}.H(t).P(t) $$

(10)

$$ {BATT}_{\mathrm{COST}}\left(K,t\right)=\left[{a}^{CH}(K){Z}_{BATT}^{CH}\left(K,t\right)+{b}^{CH}(K){P}_{BATT}^{CH}\left(s,K,t\right)\right]+\left[{a}^{DCH}(K){Z}_{BATT}^{DCH}+{b}^{DCH}{P}_{BATT}^{DCH}\left(s,K,t\right)\right]+ CC(K) $$

(11)

$$ \begin{aligned}{BATT}_{\mathrm{COST}}\left(PV,t\right)= &\left[{a}^{CH}(PV){Z}_{BATT}^{CH}\left(PV,t\right)+{b}^{CH}(PV){P}_{BATT}^{CH}\left(s_{PV},PV,t\right)\right]\\ &\quad +\left[{a}^{DCH}(PV){Z}_{BATT}^{DCH}+{b}^{DCH}{P}_{BATT}^{DCH}\left(s_{PV},PV,t\right)\right]\\ &\quad + CC(PV)\end{aligned} $$

(12)

$$ SU\left(i,t\right)=M\left(i,t\right)\times \left(1-M\left(i,t-1\right)\right)\kern1em i\in CHP, FC,B $$

(13)

$$ SD\left(i,t\right)=\left(1-M\left(i,t\right)\right)\times M\left(i,t-1\right)\kern1em i\in CHP, FC,B $$

(14)

$$ {\displaystyle \begin{array}{l}{P}_{\mathrm{buy}}\left(s,t\right)+\sum \limits_{CH P=1}^{{CH P}_N}{P}_{G, CHP}(t)+{P}_G^{FC}(t)+\sum \limits_{W=1}^{W_N}{P}_G^W\left({s}_W,W,t\right)+{P}_{BATT}^{DCH}\left({s}_{PV}, PV,t\right)+{P}_G^{TST}\left({s}_{TSS}, TST,t\right)+{P}_{BATT}^{DCH}\left(s,K,t\right)=\\ {}{P}_{\mathrm{sale}}\left(s,t\right)+{P}_{BATT}^{CH}\left({s}_{PV}, PV,t\right)+{P}_{BATT}^{CH}\left(s,K,t\right)+\left\{\left(1- DR\left(s,t\right).{L}_0\right(s,t\left)+{L}_{\mathrm{shift}}\right(s,t\Big)\right\}\end{array}} $$

(15)

where i is the index of each energy resources; W, CHP, PV, FC, TST, K, and B are the index of wind farm, combined heat and power, photovoltaic, fuel cell, tidal steam turbine, electrical energy storage device, and boiler; A _W , A _TST , A _PV , A _CHP  − F _CHP , A _FC , and B _FC are the cost coefficients of wind farm, tidal steam turbine, photovoltaic, combined heat and power, and fuel cell; T/t is the total number/index of time intervals; S _P , S _TSS , S _W , S _PV , S _PL , and S _HL are the index of scenarios for market price, tidal steam speed, wind speed, photovoltaic power generation, power, and thermal load demand, respectively; Y is the sufficient large number; U_COST(i, t) and D_COST(i, t) are the startup/shutdown cost of ith generation unit at hour t; M(i, t) is the commitment state of ith generation unit at hour t; ρ _s is the probability of the s _W th wind speed; s _TSS th is the tidal steam speed; s _PV th is the photovoltaic generation; s _p th is the scenario of market price; s _PL th is the scenario of power load demand; s _HL th is the scenario of thermal load demand; C _T (i, t) is the value of total generation cost of ith generation unit at hour t; E _P (s _p , t) is the price of the market ($/MW) for energy for s _p th scenario of price at hour t, respectively; P_sale(s, t) and P_buy(s, t) are the amount of power sold and bought to/from the market at hour t in MW; ${P}_G^W\left({s}_W,t\right)$, ${P}_G^{TST}\left({s}_{TSS},t\right)$, P_{G, CHP}(t), and ${P}_G^{FC}(t)$ are the power generation of wind farm, tidal steam turbine, heat and power, and fuel cell at hour t in MW, respectively; SU(i, t)/SD(i, t) is the startup/shutdown status of ith unit at hour t.

Demand response program constrains

The aim of demand response programs is shifting the load of MG from high consumption hours (where the energy prices are high) to the low consumption hours. It should be noted that planning for load shifting is just able to change a part or percentage of load from an hour to another [21].

The final load after applying DR program:

$$ L\left(s,t\right)=\left(1- DR\left(s,t\right)\right)\times {L}_0\left(s,t\right)+{L}_{\mathrm{shift}}\left(s,t\right) $$

(16)

The maximum amount of movable load:

$$ DR\left(s,t\right)\le {DR}_{\mathrm{max}} $$

(17)

Maximum limit of load in each of the intervals:

$$ 0\le {L}_{\mathrm{increased}}\left(s,t\right)\le {\varepsilon}_{\mathrm{increased}}\left(s,t\right)\times {L}_0\left(s,t\right) $$

(18)

Load in hour t after applying DR program:

$$ {L}_{\mathrm{increased}}\left(s,t\right)={L}_{\mathrm{shift}}\left(s,t\right)-\left( DR\left(s,t\right)\times {L}_0\left(s,t\right)\right) $$

(19)

Coefficient limit of increasing load:

$$ {\varepsilon}_{\mathrm{increased}}\left(s,t\right)\le {\varepsilon}_{\mathrm{max}} $$

(20)

Since in day-ahead markets, in general, clearing is performed for 24 h ahead, it is assumed that the daily energy expenditure of MG is fixed based on Eq. (21).

$$ \sum \limits_{t=1}^T{L}_{\mathrm{increased}}\left(s,t\right)=\sum \limits_{t=1}^T\left( DR\left(s,t\right)\times {L}_0\Big(s,t\Big)\right) $$

(21)

where L₀(s, t)/L(s, t), DR(s, t), DR_max, L_shift(s, t), ε_increased(s, t), ε_max are load before/after applying demand response program, percentage of load shifting from hour t, maximum load which can be shifted, shifted load from other hours to hour t for sth scenario, amount of increased load at hour t, maximum amount of load which can be increased at hour t, respectively.

CHP units constraints

As observed in Fig. 5, the electrical power generations of CHP units are not independent of their thermal power and these two powers cannot be controlled in a separate manner [25]. The electrical-thermal characteristics of CHP units are shown Fig. 5. The operation constraints of CHP units can be extracted from Fig. 5. The area under curve is formulated through Eq. (22). Eqs. (23 and 24) represent the models for areas above and curves, respectively. Both electrical and thermal powers are equal zero in the case of non-participating CHP units in energy generation according to Eqs. (25 and 26), respectively [26].

$$ {P}_{G, CHP}(t)-{P}_{G, CHP}\left(\alpha \right)-\frac{P_{G, CHP}\left(\alpha \right)-{P}_{G, CHP}\left(\beta \right)}{H_{G, CHP}\left(\alpha \right)-{H}_{G, CHP}\left(\beta \right)}\left({H}_{G, CHP}(t)-{H}_{G, CHP}\left(\alpha \right)\right)\le 0 $$

(22)

Fig. 5

The electrical-thermal characteristic of CHP units

$$ {P}_{G, CHP}(t)-{P}_{G, CHP}\left(\beta \right)-\frac{P_{G, CHP}\left(\beta \right)-{P}_{G, CHP}\left(\theta \right)}{H_{G, CHP}\left(\beta \right)-{H}_{G, CHP}\left(\theta \right)}\left({H}_{G, CHP}(t)-{H}_{G, CHP}\left(\beta \right)\right)\ge -\left(1-M\left( CHP,t\right)\right)\times Y $$

(23)

$$ {P}_{G, CHP}(t)-{P}_{G, CHP}\left(\theta \right)-\frac{P_{G, CHP}\left(\theta \right)-{P}_{G, CHP}\left(\lambda \right)}{H_{G, CHP}\left(\theta \right)-{H}_{G, CHP}\left(\lambda \right)}\left({H}_{G, CHP}(t)-{H}_{G, CHP}\left(\theta \right)\right)\ge -\left(1-M\left( CHP,t\right)\right)\times Y $$

(24)

$$ 0\le {H}_{G, CHP}(t)\le {H}_{G, CHP}\left(\beta \right)\times M\left( CHP,t\right) $$

(25)

$$ 0\le {P}_{G, CHP}(t)\le {P}_{G, CHP}\left(\beta \right)\times M\left( CHP,t\right) $$

(26)

where α, β, θ, and λ are the four marginal points of the electrical-thermal characteristic of combined heat and power.

Heat energy storage device constraints

n general, the heat buffer tank is added to CHP units and boiler and serves as thermal storage. The total amount of generated thermal energy is obtained from Eq. (27). The delivered thermal energy to the buffer tank influenced is by losses (η_loss) and generated excess heat (η_gain), respectively. The shutdown and startup modes of CHP units and boiler at hour t can be extracted through Eq. (28) [21]. Therefore, the thermal power available for buffer tank at hour t is computed through Eq. (29). In Eq. (30), the heat storage capacity of buffer tank is calculated. The gradient rate of increasing and decreasing thermal energy are calculated through Eqs. (31 and 32), respectively. The heat capacity limits of boiler are expressed as Eq. (33).

$$ {H}_T^{\prime }(t)=\sum \limits_{CHP=1}^{CHP_N}{H}_{G, CHP}(t)+{H}_{G,B}(t) $$

(27)

$$ {H}_T(t)={H}_T^{\prime }(t)-{\eta}_{\mathrm{loss}} SU\left(i,t\right)+{\eta}_{\mathrm{gain}} SD\left(i,t\right)\kern1em i\in CHP,B $$

(28)

$$ AH(t)=\left(1-\sigma \right) AH\left(t-1\right)+{H}_T(t)-{H}_{HL}(t) $$

(29)

$$ {AH}_{\mathrm{min}}\le AH(t)\le {AH}_{\mathrm{max}} $$

(30)

$$ AH(t)- AH\left(t-1\right)\le {AH}_{\mathrm{charge}}^{\mathrm{max}} $$

(31)

$$ AH\left(t-1\right)- AH(t)\le {AH}_{\mathrm{dischage}}^{\mathrm{max}} $$

(32)

$$ {H}_{G\min}^B(B).M\left(B,t\right)\le {H}_G^B\left(B,t\right)\le {H}_{G\max}^B(B).M\left(B,t\right) $$

(33)

where H^′(t), H(t), AH(t), and σ are the total produced heat in combined heat and power and boiler at hour t, real heat which the buffer tank could be supplied at hour t, available heat in the buffer tank, and heat loss rate for heat buffer tank, respectively.

Fuel cell constraints

The power capacity limits of boiler:

$$ {P}_{G\min}^{FC}(FC).M\left( FC,t\right)\le {P}_G^{FC}\left( FC,t\right)\le {P}_{G\max}^{FC}(FC).M\left( FC,t\right) $$

(34)

Minimum up time constraint:

$$ \left[{t}_{\mathrm{on}}\Big( FC,\left(t-1\right)\Big)-{U}_{\mathrm{min}}(FC)\right]\times \left[U\Big( FC,\left(t-1\right)-U\Big( FC,t\Big)\right]\ge 0 $$

(35)

Minimum down time constraint:

$$ \left[{t}_{\mathrm{OFF}}\Big( FC,\left(t-1\right)\Big)-{D}_{\mathrm{min}}(FC)\right]\times \left[U\left( FC,t\right)-U\Big( FC,\left(t-1\right)\Big)\right]\ge 0 $$

(36)

Generation rate constraints: [14].

$$ {R}^{\mathrm{up}}\left( FC,t\right)=\min \left\{{R}_{\mathrm{max}}^{\mathrm{up}}\Big( FC,{P}_{G\max}^{FC}(FC)-{P}_G^{FC}\Big( FC,t\Big)\right\} $$

(37)

$$ {R}^{\mathrm{down}}\left( FC,t\right)=\min \left\{{R}_{\mathrm{max}}^{\mathrm{down}}(FC),{P}_G^{FC}\Big( FC,t\Big)-{P}_{G\min}^{FC}(FC)\right\} $$

(38)

where ${t}_{\mathrm{on}}\left( FC,t\right)/{t}_{\mathrm{off}}\left( FC,t\right),{U}_{\mathrm{min}}/{D}_{\mathrm{min}},{R}^{\mathrm{up}}\left( FC,t\right)/{R}^{\mathrm{down}}\left( FC,t\right),\mathrm{and}\;{R}_{\mathrm{max}}^{\mathrm{up}}/{R}_{\mathrm{max}}^{\mathrm{down}}$ are the duration for which fuel cell had been continuously up/down till period t, minimum up/down time of fuel cell, ramp up/down capacity of fuel cell at hour t, and maximum ramp up/down rate of fuel cell, respectively.

WFs constraints

Characteristics of power generation for wind farms are nonlinear according to the wind speed which varies under the influence of type, dimension, and design of turbine. Generally speaking, the generation power of wind unit can be obtained by Eq. (39).

Moreover, A and B are determined by Eqs. (40 and 41) [23].

$$ {\displaystyle \begin{array}{l}{P}_G^W\left({s}_W,W,t\right)=\\ {}\left\{\begin{array}{cc}0& 0\le WS\left({s}_W,W,t\right)<{WS}_{ci}(W)\\ {}{P}_{WN}(W).\left(A(w).{WS}^3\left({s}_W,W,t\right)-B(W)\right)& {WS}_{ci}(W)\le WS\left({s}_W,W,t\right)<{WS}_n(W)\\ {}{P}_{WN}(W)& {WS}_n(W)\le WS\left({s}_W,W,t\right)<{WS}_{co}(W)\\ {}0& {WS}_{co}(W)\le WS\left({s}_W,W,t\right)\end{array}\right\}\end{array}} $$

(39)

$$ A(W)=\frac{1}{WS_n^3(W)-{WS}_{ci}^3(W)} $$

(40)

$$ B(W)=\frac{WS_{ci}^3(W)}{WS_n^3(W)-{WS}_{ci}^3(W)} $$

(41)

where WS(s _W , W, t), P _WN (W), WS _ci (W), WS _n (W), and WS _co (W) are the amounts of wind speed of Wth wind farm for s _W th wind speed scenario at time t, rated power of Wth wind farm, the minimum wind speed required to start power efficiency in wind farm (cut-in speed), rated speed, and the cut-out wind speed (the wind speed by which turbine puts the blades parallel to wind to prevent damages), respectively.

PV constraints

Limits on the ESD of the PV unit while getting charged and discharged:

$$ 0\le {P}_{BATT}^{CH}\le \left( PV,{s}_{PV},t\right)\le {P}_{BATT}^{\max \kern0.24em CH}(PV){Z}_{CH}\left( PV,{s}_{PV},t\right) $$

(42)

$$ 0\le {P}_{BATT}^{DCH}\le \left( PV,{s}_{PV},t\right)\le {P}_{BATT}^{\max \kern0.24em DCH}(PV){Z}_{DCH}\left( PV,{s}_{PV},t\right) $$

(43)

Charge/discharge switching constraint:

$$ 0\le {Z}_{CH}\left( PV,t\right)+{Z}_{DCH}\left( PV,t\right)\le 1 $$

(44)

Initial/terminal energy of the battery:

$$ ENR\left( PV,t=1\right)={ENR}_{\mathrm{ini}}(PV), ENR\left( PV,t=24\right)\ge {ENR}_{\mathrm{end}} $$

(45)

Amount of saved energy in the battery:

$$ ENR\left({s}_{PV}, PV,t\right)= ENR\left({s}_{PV}, PV,t-1\right)+{P}_{BATT}^{CH}\left({s}_{PV}, PV,t-1\right)-{P}_{BATT}^{DCH}\left({s}_{PV}, PV,t-1\right) $$

(46)

The power generation of the PV unit:

$$ {P}_{\mathrm{sale}}^{PV}\left({s}_{PV}, PV,t\right)={P}_G^{PV}\left({s}_{PV}, PV,t\right)-{P}_{BATT}^{CH}\left({s}_{PV}, PV,t\right)+\delta {P}_{BATT}^{DCH}\left({s}_{PV}, PV,t\right) $$

(47)

where Z _CH (PV, t)/Z _DCH (PV, t), BATT _COST , ${P}_{BATT}^{CH}\left( PV,t\right)$, ${P}_{BATT}^{\max CH}(PV)$, ${P}_{BATT}^{DCH}\left( PV,t\right)$, ${P}_{BATT}^{\max DCH}(PV)$, ENR(PV, t), and δ are the charge/discharge state of energy saving device of photovoltaic unit at hour t, the cost of buying energy for battery charging, the amount of charging power of PVth photovoltaic unit at hour t and its maximum limit, the amount of power delivered while discharging energy-saving device of photovoltaic unit at hour t and its maximum limit, the amount of saved energy in PVth energy-saving device of photovoltaic unit at hour t, efficiency factor of electrical energy-saving device, respectively.

Electrical energy storage device constraints

The constraints of electrical energy storage devices correspond to Eq. (42–47). This difference is in charging and discharging of these devices and other constraints depend on scenarios pertaining to WS, PVPG, power/thermal load demand, while the PV constraints are just affected by PVPG.

Tidal turbine

In order to extract the tidal energy and generate power electricity, the following two systems are applied:

• 1. Tidal stream system where kinetic energy of the free-flowing water is consumed and

• 2. Tidal barrage system that consumes potential energy of the ocean at ebb and flow. Usually, this method is not adopted due to the environmental conditions [27].

The generation power of tidal stream turbine is obtained through Eq. (48).

$$ {\displaystyle \begin{array}{l}{P}_G^{TST}\left({s}_{TSS}, TST,t\right)=\\ {}\left\{\begin{array}{cc}0& 0\le TSS\left({s}_{TSS}, TST,t\right)<{TSS}_{ci}(TST)\\ {}{P}_{TST}^N(TST).{\left(\frac{TSS\left({s}_{TSS}, TST,t\right)-{TSS}_{ci}(TST)}{TSS_n(TST)-{TSS}_{ci}(TST)}\right)}^3& {TSS}_{ci}(TST)\le TSS\left({s}_{TSS}, TST,t\right)<{TSS}_n(TST)\\ {}{P}_{TST}^N(TST)& {TSS}_n(TST)\le TSS\left({s}_{TSS}, TST,t\right)<{TSS}_{co}(TST)\\ {}0& {TSS}_{co}(TST)\le TSS\left({s}_{TSS}, TST,t\right)\end{array}\right\}\end{array}} $$

(48)

where

$$ {P}_{TST}^N(TST)=0.5.{\rho}_{\mathrm{water}}.A.{C}_P.{TSS}_n^3(TST) $$

TSS(s _TST , TST, t), ${P}_{TST}^N(TST)$, TSS _ci (TST), TSS _n (TST), TSS _co (TST), ρ_water, A, and C _P are the amounts of steam speed of TSTth tidal steam turbine for s _TST th steam speed scenario at time t, the rated power of TSTth tidal steam turbine, cut-in steam speed, rated steam speed, cut-out speed of tidal steam turbine, fluid density ($\raisebox{1ex}{$\mathrm{kg}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.$), the cross-sectional area of the tidal steam turbine (m²), and the power coefficient, respectively.

Results and discussions

In this part, firstly, the structure of MG and numerical data concerned with energy resources are studied, and then, simulation results of optimal operation for the stochastic problem are analyzed.

Configuration of MG

In this article, three case studies will be assessed:

• 1. Planning isolated MG by predicting uncertainties by hybrid method of WT-ANN-ICA

• 2. Planning and determining the optimal strategy of MG energy resources connected to grid and comparing hybrid prediction methods of WT-ANN and WT-ANN-ICA, in order to examine the influence of predicting uncertainties upon the profit amount of MG

• 3. Programming and determining the optimal strategy of MG connected to the grid and applying hybrid prediction method of WT-ANN-ICA for predicting uncertainties and exploring the effect of DR problem on the profit of MG

The case studies are run on three WFs, two CHP units, one TST, one PV, one boiler, one low-temperature fuel cell (PAFC), one electrical energy storage device, one heat buffer tank together with the fixed and responsive electrical, and the related fixed thermal loads. The startup and shutdown costs of units are tabulated in Table 1. The heat buffer tank data and cost coefficients of CHP units are tabulated in Table 2. The heat buffer tank data and cost coefficients of CHP units are shown in Table 3. Both DR_max and ε_max are assumed 30%. The electrical-thermal characteristics of CHP units are displayed in Fig. 5. The parameters of WFs include WS_co(i) = 25^m/s, WS _n (i) = 11^m/s, WS_ci(i) = 2.5^m/s, and the rated output power are equal to P_WN1 = 1.5^MW, P_{WN2, 3} = 2.4^MW. Historical data pertaining to the WS, electrical demand and market price, electrical energy storage devices data, and photovoltaic power generation are proposed in [21, 28, 29] and [16], respectively. The PV nominal power generation is ${P}_{\mathrm{max}}^{MW}=4.68,{P}_{\mathrm{min}}^{MW}=0$ and δ = 0.75. Table 4 lists the parameters used for the tidal steam turbine [30].

Table 2 The startup and shutdown cost of units

Table 3 The heat buffer tank data and cost coefficients of CHP units

Table 4 The tidal steam turbine data

Case studies

Case study 1: planning of MG in the grid-isolated mode

In this case, P_buy = P_sale = 0 and the objective function is just concluded in cost terms of energy resources. The uncertainties are predicted by hybrid method of WT-ANN-ICA. The simulation results obtained from the first case study are listed in Table 5. According to Table 5, the generation cost is equal to the amount of objective function and $2032.76. In this state, the WF, TST, and PV are not working with their maximum capacity, while the cost of generating them is zero and this occurs due to the thermal load of MG. The CHP units produce heat to provide thermal load. The electrical-thermal characteristics of CHP units are the generation factor for both heat and electrical power (Fig. 6).

Table 5 Case studies results

Fig. 6

Typical MG under study

2. Case study 2

The effect of exchanging electrical energy with grid in connected mode and also the effect of more accurate prediction of random parameters on MG planning are studied by comparing the hybrid methods of WT-ANN-ICA and WT-ANN. The MG planning problem in the presence of all economic and technical constraints and the DR problems will be solved. The results are tabulated in Table 5, where, by WT-ANN-ICA prediction method, the generation cost increases by 13.77% in comparison with the first case. The profit of MG resulted from taking part in the market is $792.64 and $909.93 for WT-ANN and WT-ANN-ICA, respectively. This profit is due to the sale of power to the main grid.

The generated power of resources in the planning horizon is shown in Fig. 7. The expected amount of buying and selling powers, WFs, TST, and PV power generation are illustrated in Fig. 7. According to Fig. 7, WFs, TST, and PV power generations are readjusted at their maximum capacity.

Fig. 7

The generated power of resources in the planning horizon

Therefore, the excess WFs, TST, and PV power generations can be saved in devices of electrical energy storage and shifted to the hours with more demand (higher price) or can be on offer to the market. In Fig. 8c, d, the generated thermal power of CHP units and boiler are shown, respectively. According to Fig. 8c, d), considering costly expenses of generation, the boiler is less involved in supplying thermal demand, as compared to CHP units.

Fig. 8

a–e Energy generation of resources with and without DR

Due to the stochastic nature of WFs, TST and PV power generations, market price, and power/thermal load demand, more accurate prediction of these random parameters led to generating scenarios proximate to reality and with greater possibility. Consequently, a more detailed planning can be achievable. In Table 5, the obtained results of MG planning are compared by prediction hybrid method of WT-ANN and WT-ANN-ICA. The results indicate 14.79% increased profits of MG in WT-ANN-ICA method in comparison with WT-ANN method.

Case study 3: to explore the effect of DR program on determination of generation strategy of MG

In the third case, the MG planning is studied without DR program in order to investigate its effects on the expected profit of MG by hybrid prediction method of WT-ANN-ICA. Table 5 shows the results of case study 3 versus 2. The expected profit decreased while the cost of generation dropped slightly. According to Table 5, the expected profit of MG in the third case study is $712.05 which approximately decreased 27.7% as compared to the second case study (applied DR program). This reduction in profit of MG indicates the efficacy of DR program on the optimal planning of these units.

Regarding Table 5, the objective function values of both cases 2 and 3 are negative. That is because of high energy demand inside the MG and lack of possibility to offer excess power to the main grid.

The difference between energy generation of resources with and without DR is illustrated in Fig. 8. The received and delivered power to grid with and without DR is depicted in Fig. 8a. Regarding this figure, the movable loads can be shifted from peak time to other hours when the energy price is lower and makes profit.

Conclusions

In this paper, an algorithm was suggested to ascertain optimal strategy of a MG including WFs, PV, TST, fuel cell, CHP units, boiler, and ESDs, by considering economic and technical constraints and DR program. This research aimed to present an optimization program to maximize the profit of MG in grid–connected mode, and to minimize the cost of energy resources in grid- isolated mode. The uncertainties are WS, TSS, PVPG market price and power/ thermal load demand which are predicted by hybrid prediction methods of WT-ANN and WT- ANN-ICA, and related scenarios are generated by probability density functions appropriate to each uncertainty and scenario reduction method. The simulation results represent that applying more accurate prediction method, some scenarios proximate to reality and with greater possibility are generated. Hence, a more detailed and precise planning is achieved and the expected profit of MG might be increased. If the method of WT-ANN-ICA is used rather than WT-ANN, the expected profit of MG will be increased by 14.79%. Furthermore, the expected profit can be raised by applying DR program. According to the studied cases, although the DR program increases the generation cost by 7%, the expected profit rises more than 27.7% and it goes to $909.93, while without considering DR program, this profit would be $712.05.

In addition, the outcome of case studies illustrate that by implementing this proposed framework, the MG can obtain a meaningful profit in the grid-connected mode in comparison with the isolated mode, as well as supplying total electrical and heat demand.