Size Optimization of a Solar-wind Hybrid Energy System Using Two Simulation Based Optimization Techniques

1. Introduction

Many of the parts in this chapter have been reproduced from Ekren and Ekren, 2008; Ekren and Ekren, 2009; Ekren et al., 2009 with permission.

Energy is an important fact for a country for both its socio-economic development and economic growth. Main energy source on earth is the fossil fuels. However, the usage of fossil fuels causes global warming whose negative effects have recently been felt by all over the world. Also, because it is limited on earth, increased energy demand and high energy prices increase concerns on fossil fuels. For instance, petroleum reservoirs are present only in few countries therefore the other countries mostly purchase petroleum from these countries more than their production amount. Hence, decrease of fossil fuel reservoirs may create a new ‘energy crisis’ and/or “energy wars” as in 1970s in near future.

For a sustainable world the usage of fossil fuels must be decreased, in fact ended. Instead, the usage of renewable energy sources must be increased. As it is known, the interest in renewable energy sources has increased because it does not cause greenhouse effect in contrary to the fossil fuels. These energy sources are indigenous, environmental friendly and, they help to reduce the usage of fossil fuels. Solar, wind, wave, biomass and geothermal energies are renewable energy sources. Sun is the source of all energies. Solar energy is usually used in two aims: for thermal applications and for electricity production. Wind is the indirect form of solar energy and is always replenished. The past and the predictedamount of global renewable energy source statusby 2040 are presented in Table 1. As seen in this table the renewable energy consumption is predicted to increase in the future.

There are several studies on renewable energy sources and hybrid combination of these sources for electricity production. Electricity has high cost mostly due to centralized energy systems which operate mostly on fossil fuels and require large investments for establishing transmission and distribution of grids that can penetrate remote regions (Deepak, 2009). Unlike the centralized energy systems, decentralized energy systems are mostly based on renewable energy sources. They operate at lower scales (a few kWh scale) both in the presence and absence of grid, and easily accessible to remote locations due to generation of

	2001	2010	2020	2030	2040
Total energy consumption (million tons oil equivalent)	10,038	10,549	11,425	12,352	13,310
Biomass	1080	1313	1791	2483	3271
Large hydro	22.7	266	309	341	358
Geothermal	43.2	86	186	333	493
Small hydro	9.5	19	49	106	189
Wind	4.7	44	266	542	688
Solar thermal	4.1	15	66	244	480
Photovoltaic	0.1	2	24	221	784
Solar thermal electricity	0.1	0.4	3	16	68
Marine(tidal/wave/ocean)	0.05	0.1	0.4	3	20
Total renewable energy consumption	1,365.5	1,745.5	2,964.4	4,289	6,351
Renewable energy source consumption (%)	13.6	16.6	23.6	34.7	47.7

Table 1.

World global renewable energy sources scenario by 2040 (adapted from Kralova andSjöblom, 2010)

power in the propinquity of demand site. They are also called “stand-alone energy systems” and produce power independently from the utility grid. Because these systems are not connected to the utility grid, they usually need batteries for storage of electricity produced during off-peak demand periods.

Solar and wind energies are usually available for most of the remote areas as renewable sources. However, it is prudent that neither a standalone solar energy nor a wind energy system can provide a continuous supply of energy due to seasonal and periodical variations. Simultaneous utilization of multiple energy resources greatly enhances the certainty of meeting demands. These systems are called hybrid energy systems (see Figure 1). Because they are good complementary energy sources of each other, solar and wind energies have been widely used as hybrid combination for electricity supply in isolated locations far from the distribution network. However, they suffer from the fluctuating characteristics of available solar and wind energy sources. Therefore, properly sized wind turbine, photovoltaic panel and storage unit provides high reliability and low initial investment cost. In this chapter, we aim to show two simulation based size optimization procedures for a solar-wind hybrid energy system providing minimum cost. The case study is completed to meet the electricity demand of a GSM base station located near a seaside region in western of Turkey. As in Figure 1, the studied hybrid system’s electricity is also produced via photovoltaic array and wind turbine which are regulated by voltage regulator components and, the excess electricity produced is stored by the battery banks to be used for later lacking loads. Here, the amount of the electricity produced via the solar energy and the wind depends on the total solar radiation on horizontal surface and the wind speed respectively.

2. Background and motivation

Several researchers have studied hybrid renewable energy sources. Panwar et al. (2011) reviewed the renewable energy sources to define the role of the renewable energy sources for environmental protection. In their study, it is emphasized that renewable technologies are clean energy sources and optimal use of these resources minimize negative environmental

Figure 1.

A solar-wind hybrid energy system

impacts, produce minimum secondary wastes and provide a sustainable world. In another study, Angelis-Dimakis et al. (2011) evaluated the availabilityof renewable energy sources such as solar, wind, wave, biomass and geothermal energy. In their research, a detailed survey including existing methods and tools to determine the potential energy in renewable resources is presented. Also, tendency of using the renewable energy by the most developed countries in order to reduce the concentration of carbon dioxide in the atmosphere is emphasized. Their study also mentions the usability of hybrid energy system by mixing different renewable sources.

A great deal of research has been carried out on hybrid energy systems with respect to performance and optimization of these systems (Kellogg et al., 1996; Protegeropoulos et al., 1997; Seeling-Hochmuth, 1997; Markvart, 1996; Bagul et al., 1996; Borowy and Salameh, 1996; Morgan et al., 1997; Celik, 2002; Yang et al., 2003; Ashok, 2007; Bernal-Agustin and Dufo-Lo´ pez, 2009; Yang et al., 2009; Bilal et al., 2010; Zhou et al., 2010).Hybrid energy system studies in the past mostly based upon a particular design scenario with a certain set of design values yielding the near optimum design solution only (Kellogg et al., 1996; Protegeropoulos et al., 1997; Seeling-Hochmuth, 1997; Bagul et al., 1996; Morgan et al., 1997; Celik, 2002; Yang et al., 2003; Ashok, 2007). Such an approach, although providing the near optimum solution, unfortunately lacks the ability to provide a general understanding about how the total system cost changes with the size of the design parameters. A graphical optimization technique to optimize the size of the solar-wind hybrid energy system is studied by Markvart (1996). He considered monthly average solar and wind energy data in optimization. On the other hand, unlike the methods based on hourly, daily, and monthly average basis, a statistical approach for optimizing the size of PV arrays and the number of batteries for a standalone solar-wind hybrid system is presented by Bagul et al. (1996). They proposed a three-event probabilistic approach to overcome the limitations of the conventional two-event approach in matching the actual distribution of the energy generated by hybrid systems. Borowy and Salameh (1996) developed an algorithm to optimize a photovoltaic-array with battery bank for a standalone solar-wind hybrid energy system. Their model is based on a long-term hourly solar radiation and peak load demand data from the site. In this study, direct cost of the solar-wind hybrid energy system is considered. Different from these studies, Morgan et al. (1997) studied performance of battery units in a standalone hybrid energy system at various temperatures by taking into account the state of voltage (SOV) instead of the state of charge (SOC). Their algorithm is able to predict the performance of a hybrid energy system at various battery temperatures. This study is important for efficiency of a hybrid energy system because temperature affects the performance of a PV array and battery unit. Celik (2002) carried out techno-economic analysis and optimization of a solar-wind hybrid energy system. Yang et al. (2003) proposed an optimization technique by considering the loss of power supply probability (LPSP) model of a solar-wind hybrid system. They demonstrated the utility of the model for a hybrid energy system for a telecommunication system. Ashok (2007) presented a model based on different components of a hybrid energy system and developed a general model to define the optimal combination of renewable energy source components for a typical rural community. Recently, Bernal-Agustin and Dufo-Lo´ pez (2009) have analyzed usability of renewable energy systems. They use simulation to optimize the generation of electricity by hybrid energy systems. According to the authors stand-alone hybrid renewable energy systems are usually more suitable than only photovoltaic (PV) or wind systems in terms of lower cost and higher reliability. On the other hand, the design, control, and optimization of the hybrid systems are usually very complex tasks because of the high number of variables and the non-linearity. They come up with that the most reliable system is composed of PV–Wind–Battery. Although PV–Diesel–Battery is also reliable, this system uses fossil fuel. Optimal design and techno-economic analysis of a hybrid solar–wind power generation system is studied by Yang et al. (2009). They show sizing, optimizing and selecting the most suitable renewable source coupleswhich are crucial for a hybrid system. They also come up with that solar and wind energies are the most suitable renewable energy resources. The authors also propose an optimal design model for solar–wind hybrid energy system. Obtaining the optimum configurations ensures decreased annual cost and, the loss of power supply probability (LPSP) is satisfied. Bilal et al. (2010) studied size optimizationof a solar-wind-battery hybrid system for Potou which is an isolated site located in the northern coast of Senegal. This area is far away from the electricity supply. The methodology used in their study consists of sizing and optimization of a hybrid energy system by multi-objective genetic algorithm and the influence of the load profiles on the optimal configuration. Optimal configurations are examined for three profiles. Profile 1 is for the load for operation of refrigerators, domestic mill, welding machines, and other equipment in the village. Profile 2 is for the load for operation of a desalination and water pumping system, commercial refrigerators and domestic equipment, etc. The third load illustrates low consumption during the day (population working in the fields in the morning), and high power demand at night. Zhou et al. (2010) presented current status of researches on optimum sizing of stand-alone solar–wind hybrid power generation systems. The authors use two renewable energy sources because of the fact that their availability and topological advantages for local power generations.Also, these combinations of hybrid energy systems allow improving the system efficiency and power reliability and reduce the energy storage requirements for stand-alone applications. It is concluded that continued research and development in this area are still needed to improve these systems’ performance.

Most of the existing studies use historical data and/or intervals for the input variables – solar-wind energies and electricity demand - of the system. Different from the existing studies, we utilize probabilistic distributions in order to carry out random input simulation. A detailed studied carried out to fit the input variables to probabilistic distributions. We fit the probabilistic distributions based on hours for each month for the solar radiation and the wind speed values. The electricity consumption of the GSM base station is also fit hourly basis. Different from the previous studies, we use two different simulation based optimization techniques – RSM and OptQuest - to compare their results. We implement a case study to model a stand-alone solar-wind hybrid energy system at a remote location from the grid system located in western of Turkey (Ekren and Ekren, 2008; Ekren and Ekren, 2009; Ekren et al., 2009).

Section 3 explains the simulation modeling of the hybrid energy system. In this section, the measured values of the system’s inputs - solar radiation, wind speed and GSM base station’s electricity consumption – are also provided. Besides, hourly fitted distributions for each month of the solar radiation and the wind speed and, hourly fitted distributions of the GSM base station’s electricity consumption are also presented. In Section 4, the first optimization methodology - RSM – is introduced and used to optimize the hybrid system. In this section, we also provide the conducted experiments. In Section 5, the second optimization methodology – OptQuest – is explained. Last, we conclude the study.

3. Simulation modeling of the solar-wind hybrid energy system

The hybrid system under study relies on solar and wind energies as the primary power resources, and it is backed up by the batteries (see Figure 1). Batteries are used because of the stochastic characteristics of the system inputs. Namely, it is used to meet the electricity demand while the solar and wind energies are not adequate. The basic input variables of the hybrid model are: solar radiation, wind speed, and the electricity consumption of the GSM base station. Because the characteristics of these variables are non-deterministic, we fit to probability distributions to carry out a Monte Carlo (MC) simulation (see Tables 3-5). The probability distributions are specified in the input analyzer tool of the ARENA simulation software. Random data for solar radiation, wind speed, and the electricity consumption are generated using these distributions in ARENA (Kelton et al., 2004). Because in the simulation model hourly data are used, one of the system’s assumptions is that the input variables do not change throughout an hour. This means that the solar radiation and the wind speed input values are constant e.g. from 12:00 pm to 1:00 pm. in any month in the model. The length of each simulation run is considered as twenty years of the economical life which consists of 365 days/year, 24 hours/day, in total 175,200 hours. For each run, 5 independent replications are completed. In the simulation model, since it is a popular and useful variance reduction technique to compare two or more alternative configurations, the common random numbers (CRN) variance reduction technique is used. And since a steady state analysis is needed to analyze a long time period non-terminating system, the warm-up period is decided as 12,000 hours (Law, 2007). Hourly mean solar radiation and wind speed data for the period of 2001-2003 (26,280 data = 24hours*365days*3years) are recorded at a meteorological station where the suggested hybrid energy system is to be established. Technical specifications of the meteorological station are given in Table 2 (Ekren, 2003).

Instrument	Specification/Description
Pyranometer (CM11)	Wieving Angle : 2π; Irradiance : 0 - 1400 W / m2; Sensitive : 5.11*10-6 Volts per W/ m2 (+/- 0.5 % at 20 oC and 500 W/m2 ); Expected Signal Output : 0 – 10 mV; Response time for 95 % response : < 15 sec.
Data Logger	Module capacity:192896 bytes; 12 signal inputs
Anemometer (for speed); Wind Vane (for direction); Thermometer; Hygrometer; Barometer	Measurement RangeRecording ResolutionAccuracy; 0.3 to 50 m/s ≤0.1 m/s ± 0.3 m/s; 00 - 3600 ≤1º ± 20; (-30)–(+70) 0C≤0.1ºC ± 0.2 K; 0-100 % RH 1 % RH ±2 % RH; 800 to 1600 kPa ≤ 1kPa -

Table 2.

Main Characteristics of the Meteorological Station

3.1. Solar radiation

Figure 2presents average measured hourly total solar radiation on horizontal surface, H, based on months in a year. Hourly total solar radiation on tilted surface, IT, is calculated using H and optimum tilted angle of the PV panel, β is taken as 380(Eke et al., 2005).

Figure 2.

Average hourly total solar radiation on horizontal surface (Measured)

Table 3 presents the fitted hourly total solar radiation distributions, for example, for three months, in June, July, and August. Each solar radiation distribution is different at each hour of the month. This means that solar radiation may vary in accordance with its distribution at each hour of the months which is not known in advance. Here, the existences of the zero values are due to the sunset.

Hours	Months
June	July	August
00:00-01:00	0	0	0
01:00-02:00	0	0	0
02:00-03:00	0	0	0
03:00-04:00	0	0	0
04:00-05:00	0	0	0
05:00-06:00	0	0	0
06:00-07:00	-0.001 + 37 * BETA(0.0269, 0.0857)	0	0
07:00-08:00	48 + 123 * BETA(2.71, 0.787)	NORM(113, 12.9)	18 + 80 * BETA(2.22, 1.22)
08:00-09:00	NORM(313, 16.6)	TRIA(240, 278, 305)	80 + 189 * BETA(3.49, 1.22)
09:00-10:00	NORM(482, 19.9)	270 + 209 * BETA(3.26, 0.724)	183 + 246 * BETA(1.57, 0.542)
10:00-11:00	481 + 173 * BETA(2.39, 0.602)	NORM(602, 17.2)	NORM(545, 55.8)
11:00-12:00	558 + 206 * BETA(2.51, 0.601)	NORM(714, 28.1)	265 + 457 * BETA(2.76, 0.465)
12:00-13:00	538 + 297 * BETA(1.34, 0.431)	517 + 337 * BETA(3.95, 1.11)	221 + 617 * BETA(3.05, 0.62)
13:00-14:00	472 + 399 * BETA(1.5, 0.42)	672 + 187 * BETA(3.27, 0.858)	145 + 704 * BETA(2.09, 0.385)
14:00-15:00	229 + 617 * BETA(1.13, 0.313)	400 + 436 * BETA(1.32, 0.415)	144 + 655 * BETA(0.978, 0.366)
15:00-16:00	359 + 411 * BETA(1.02, 0.402)	453 + 311 * BETA(0.931, 0.405)	334 + 393 * BETA(1.45, 0.507)
16:00-17:00	198 + 450 * BETA(1.09, 0.303)	174 + 481 * BETA(1.34, 0.407)	193 + 464 * BETA(2.78, 1.04)
17:00-18:00	178 + 326 * BETA(1.28, 0.418)	46 + 452 * BETA(1.16, 0.347)	76 + 371 * BETA(1.47, 0.483)
19:00-20:00	103 + 233 * BETA(1.55, 0.387)	69 + 258 * BETA(2.46, 0.527)	TRIA(78, 266, 278)
20:00-21:00	63 + 91 * BETA(1.54, 0.604)	-0.001 + 148 * BETA(3.43, 1.02)	-0.001 + 101 * BETA(0.974, 1.13)
22:00-23:00	0	0	0
23:00-00:00	0	0	0

Table 3.

Hourly solar radiation on horizontal surface distributions for months June, July, August (W/m2)

ARENA simulation software uses nine different theoretical distributions to fit data to a theoretical distribution. These are: Exponential, Gamma, Lognormal, Normal, Triangular, Uniform, Weibull, Erlang, and Beta distributions. Each of the distribution has its own probabilistic characteristics in creating random variables in a stochastic model.

3.2. Wind speed

In order to measure the wind speed and the prevailing wind direction, a three-cup anemometer and a wind vane are used. Hourly average wind speed at 10-meter-height for all months of the year, can be seen in Figure 3. These average hourly measured solar radiation and wind speed data figures are given for a general idea of the energy potential of the area. Otherwise, in the simulation model we use long-term dynamic hourly data. The height used to measure the wind speed is the universally standard meteorological measurement height (AWS Scientific, 1997). Table 4 shows hourly wind speed distributions of three months as an example.

Figure 3.

Average hourly average wind speeds (measured)

Hours	Months
January	February	March
00:00-01:00	TRIA(2, 3.08, 18)	2 + 13 * BETA(1.07, 1.57)	1 + 12 * BETA(1.46, 1.79)
01:00-02:00	1 + LOGN(7.02, 5.28)	WEIB(7.75, 1.88)	1 + 11 * BETA(1.54, 1.66)
02:00-03:00	2 + 14 * BETA(0.881, 1.23)	GAMM(2.24, 3.31)	TRIA(2, 3.95, 13)
03:00-04:00	NORM(8.22, 3.96)	1 + LOGN(7.1, 5.85)	TRIA(1, 4.52, 14)
04:00-05:00	2 + WEIB(6.7, 1.52)	1 + LOGN(6.61, 5.58)	NORM(6.93, 2.66)
05:00-06:00	TRIA(2, 4.21, 18)	1 + LOGN(6.3, 5.69)	2 + WEIB(5.3, 1.78)
06:00-07:00	1 + 15 * BETA(1.37, 1.6)	2 + WEIB(5.5, 1.2)	TRIA(1, 7.62, 12)
07:00-08:00	1 + LOGN(7.23, 5.21)	2 + WEIB(5.87, 1.48)	TRIA(1, 7, 13)
08:00-09:00	2 + WEIB(6.84, 1.66)	2 + 17 * BETA(1.14, 2.28)	TRIA(1, 7.69, 10.9)
09:00-10:00	1 + WEIB(8.11, 2.12)	2 + WEIB(5.98, 1.41)	1 + GAMM(1.13, 4.69)
10:00-11:00	2 + WEIB(6.71, 1.8)	1 + WEIB(7.1, 1.61)	2 + WEIB(5.15, 2.25)
11:00-12:00	2 + WEIB(6.78, 1.87)	1 + ERLA(2.11, 3)	NORM(6.77, 2.54)
12:00-13:00	3 + ERLA(2.63, 2)	2 + 22 * BETA(0.955, 2.76)	NORM(6.68, 2.81)
13:00-14:00	2 + GAMM(2.45, 2.54)	1 + WEIB(7.24, 1.56)	1 + ERLA(1.98, 3)
14:00-15:00	2 + WEIB(6.97, 1.58)	1 + WEIB(6.92, 1.57)	1 + WEIB(6.97, 1.76)
15:00-16:00	2 + WEIB(6.54, 1.6)	1 + WEIB(7.06, 1.61)	NORM(7.17, 3.63)
16:00-17:00	2 + WEIB(6.41, 1.63)	1 + WEIB(6.86, 1.6)	2 + 17 * BETA(1.37, 3.1)
17:00-18:00	NORM(7.35, 3.12)	1 + ERLA(2.8, 2)	NORM(6.98, 2.9)
19:00-20:00	TRIA(2, 4.7, 16)	NORM(6.84, 3.82)	1 + LOGN(6.07, 3.79)
20:00-21:00	2 + WEIB(5.76, 1.57)	1 + WEIB(6.46, 1.57)	2 + 12 * BETA(1.06, 1.56)
22:00-23:00	1 + WEIB(7.01, 1.73)	1 + WEIB(6.4, 1.62)	1 + GAMM(1.84, 3.1)
23:00-00:00	NORM(7.41, 4.06)	1 + WEIB(6.4, 1.62)	TRIA(2, 4.39, 13)

Table 4.

Hourly average wind speed distributions for months, January, February, March (m/sec)

3.3. Electricity consumption

The third stochastic data is the electricity consumption of the GSM base station. The data are collected from the base station for every hour of the day. In this study, the existence of a seasonal effect on the GSM base station’s electricity consumption is ignored. The statistical data are collected in 15 random days in each season, fall, winter, spring, and summer. Hence, totally 15*4 = 60 electricity consumption data are collected for an hour (e.g. for 1.00 pm). Then these data are fit to theoretical distributions without considering seasonal effects. Fig. 4 illustrates the hourly mean electricity consumption values of the GSM base station. And the fitted distributions are given in Table 5.

The output of the wind generator and PV panels are DC power and, inverter converts it to the AC power.

Figure 4.

Hourly average demand of the GSM base station

3.4. Formulations

If the hybrid energy systems are well designed, they provide a reliable service for an extended period of time. In the optimization procedure, the sizes of system components are decision variables, and their costs are objective function. The cost of the system is considered as the total cost of PV, wind turbine rotor, battery, battery charger, installation, maintenance, and engineering. A solar and wind hybrid energy system with the sizes of as and aw, respectively, can be defined by (1)–(2):

$$a_s=\eta .A_s$$ E1

where, η is the PV module efficiency, Asis the PV array area and:

$$a_w=C_p.\left(\pi .r^2\right)$$ E2

where, Cp is the power coefficient, and r is the rotor radius. Here, π.r2 represents Aw, rotor swept area. η value is taken as a variable value depending on PV module type and module temperature. In the study, mono-crystal silicon PV module type, rated output of 75 W at

Hours	Demand
00:00-01:00	NORM(1.21, 0.602)
01:00-02:00	1.72 * BETA(1.13, 1.55)
02:00-03:00	WEIB(0.645, 1.28)
03:00-04:00	GAMM(0.565, 1.11)
04:00-05:00	EXPO(0.854)
05:00-06:00	3.31 * BETA(1.28, 3.34)
06:00-07:00	2.86 * BETA(1.15, 1.96)
07:00-08:00	WEIB(1.7, 2.36)
08:00-09:00	TRIA(0, 0.827, 3.84)
09:00-10:00	3.84 * BETA(2.46, 3.25)
10:00-11:00	0.15 + 3.34 * BETA(2.54, 2.87)
11:00-12:00	0.28 + LOGN(1.91, 0.951)
12:00-13:00	1 + ERLA(0.828, 2)
13:00-14:00	1 + LOGN(1.82, 1.43)
14:00-15:00	NORM(2.27, 0.626)
15:00-16:00	1 + GAMM(0.669, 2.4)
16:00-17:00	1 + GAMM(0.417, 3.45)
17:00-18:00	NORM(2.47, 0.685)
19:00-20:00	1 + ERLA(0.377, 4)
20:00-21:00	1.11 + ERLA(0.436, 4)
22:00-23:00	1.29 + LOGN(2.2, 1.59)
23:00-00:00	1 + GAMM(0.435, 3.48)

Table 5.

Hourly electricity demands (kW)

1000 W/m2 is used. Here, η value changes between 7% and 17% based on module surface temperatures which are between 100C and 700C (Kemmoku, 2004; Muselli, 1999). In the simulation model, for December, January and February the temperature and the η values are assumed to be 100C and 17%, respectively. For March, April and May the temperature and the η values are assumed to be 500C and 10%, respectively. For June, July, August the temperature and the η values are assumed to be 700C and 7%, respectively. And, for September, October, November the temperature and the η values are assumed to be 300C and 13%, respectively. These values are obtained from a manufacturer firm. Cp value is also taken from a manufacturer firm as a graphic value which changes according to the wind speed value. Wind energy density, W, is calculated by (3):

$$W=\text{1}/\text{2}.\rho .V^3.D$$ E3

D is length of period. Because of hourly operating state, it is taken as 1 hour. ρ is air density which is considered as 1.225, and V is hourly average wind velocity whose distribution is shown in Table 4.

Solar radiation on tilted plate is calculated for isotropic sky assumption by (4a)-(4b):

$$\begin{array}{c}{I}_T =I_b R_B +I_d (\text{1}+\text{cos}\beta )/\text{2} +\left(I_b+I_d \right)r (\text{1}-\text{cos}\beta )/\text{2}\\ R_b=\text{cos}\theta /\text{cos}\theta z\end{array}$$ E4

where, IT is total solar radiation on tilted surface, Ib is horizontal beam radiation, Id is horizontal diffuse radiation, Rb is ratio of beam radiation on tilt factor, θ is incidence angle, θz is zenith angle, ρis surface reflectivity, βis tilted angle of the plate.

The supply of hourly solar and wind energies must meet the hourly demand, d. This expression can be formulated as:

$$S.a_s+W.a_w\ge d$$ E5

where S is the solar energy density on tilted surface, HT (kWh/m2) and W is the wind energy density (kWh/m2). If (5) is not realized, the stored energy in the battery will be used. The battery’s efficiency is assumed as 85%. The inverter’s efficiency is considered 90% here. If the total of solar, wind, and battery energies still cannot meet the demand, the energy shortage will be supplied by an auxiliary energy source, whose unit cost, herein, is considered $0.5 per kWh electricity. In this study, the location of the hybrid system is assumed in such a place where the unit cost of the auxiliary energy is more expensive than the electricity produced by the hybrid system. Therefore, the cost of extra energy here is decided such that it is three times as much of the average unit cost of the electricity produced by the hybrid system. Otherwise, if the unit cost of the auxiliary energy was less than that of the electricity produced by the hybrid system, the shortage could be supplied from the auxiliary energy source all the time without the need for such a hybrid system (Celik, 2002). As seen in (6), the auxiliary energy cost is also a part of the total hybrid system cost. Cs, Cw, CB, Csh and CT are the unit cost of photovoltaic, wind energy generator, battery, shortage electricity, and total hybrid energy system, respectively. BC denotes the battery capacity. Ei is the total amount of the electricity energy shortage because of not meeting the demand during an hour i. Cs are $5.8/Wp (mono-crystal silicon, rated output of 75 W at 1000 W/m2), $5.5/Wp (multi-crystal silicon, rated output of 75 W at 1000 W/m2) whereas Cw is US $3/W (rated output of 5000 W at 10 m/s), and CB is $180/kWh (200 Ah 12V lead acid battery) (Eke et al., 2005). Csh is $0.5 per kWh as explained above, and n is the simulation time period, 175,200 hours. In this study, the inflation rate and time value of money are not considered.

$$C_T= C_s.a_s+C_w.a_w+C_B.B_C +{\displaystyle \sum _{i=1}^n{C}_{sh}.E_i}$$ E6

In addition, a total of US $500 battery charger cost, and 5% installation, maintenance and engineering cost of the initial hardware is also added into the total system cost for an assumed 20-year-lifetime.

4. Response surface methodolgy

RSM consists of a group of mathematical and statistical techniques useful for developing, improving and optimizing processes(Myers et al., 2009). It is also useful in the design, development and formulation of new systems as well as improvement of existing systems. In RSM, usually there are several input variables that affect some performance measure which is also called response. Most real world applications generally involve more than one response. The input variables are sometimes called independent variables which can be controlled by the controller for purposes of a test or an experiment (Myers et al., 2009).

Figure 5 shows an example for the relationship between the response variable (bean yield) and the two input variables (PhosAcid and Nitrogen) in a chemical process graphically. Since there is a response lying above the two input variables’ plane the term “response surface” comes from this reason.

Figure 5.

Response surface of a chemical process

RSM requires developing an approximating model also called metamodel for the true response surface. The metamodel is based on observed data from the system and is an empirical model. Metamodels are developed to obtain a better understanding of the nature of the true relationship between the input variables and the responses of the system. This approximate formula could be used as a proxy for the full-blown simulation itself in order to get at least a rough idea of what would happen for a large number of input-parameter combinations. Multiple regression is a technique useful for building the types of empirical models required in RSM (Kleijnen, 1987; Friedman, 1996; Kleijnen and Sargent, 2000).

A first-order response surface metamodel is given by (7) as an example.

$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\epsilon$$ E7

where y represents the response, x1 and x2 represent the input variables. is a fixed variable β0 and β1 and β2 are the coefficients of x1 and x2,respectively. Hence, this is a multiple linear- regression with two independent (input) variables.

In RSM, first the significant input variables are identified. While identifying the significant input variables, it is also important to determine the current level of these variables. Here, the objective is to select the levels that provide the response close to its optimum. An RSM begins when the process is near optimum. At this phase it is aimed to predict the model that estimates the true response with a small region around the optimum. Because there is usually a curvature near the optimum a second-order model is mostly used. After determining the metamodel of the system, RSM uses steepest ascent method to optimize it(Myers et al., 2009).

The sequential nature of RSM provides the experimenter to have an idea about these five steps: a) the optimum point’s location region b) the type of metamodel required c) the proper choice of experimental designs d) how much replication is necessary e) whether or not transformations on responses or any of process variables are required (Myers et al., 2009). The RSM process typically involves taking observations in a starting region, usually according to an experimental design such as a factorial (2k) or fractional (2k-p). Here, k is the number of factors (input variables) and 1/2p is called the degree of fractionation, because it represents the fraction of observations from a 2k design that is required. For example, if a problem having 5 factors is studied, the necessary number of runs in the experiment would be 25=32. Because each run may require time-consuming and costly setting and resetting of machinery, it is often not feasible to require many different production runs for the experiment. In these conditions, fractional factorials (2k-p) are used that “sacrifice” interaction effects so that main effects may still be computed correctly (Law, 2007).

After obtaining the metamodel, it is necessary to examine whether the fitted model provides an adequate estimate on the true system or not. Also, it is necessary to verify whether the pre-defined regression assumptions are satisfied or not. There are several techniques to check the model adequacy (Myers et al., 2009). In this study, since it is easy and popular we use residual analysis to check the model adequacy.

4.1. Residual analysis

Normality assumption should be checked in metamodel fitting. This assumption is satisfied if the normality plot of the residuals is as in Figure 6a and it is not if it is as in Figure 6b. Namely, if the residuals plot approximately along a straight line then the normality assumption is satisfied (Myers et al., 2009).

Figure 6.

Probability plot of normal data b. Probability plot of nonnormal data

4.2. Lack of fit test

In RSM, it is important to check whether the regression function is a good fit or not. We use the lack of fit test to decide this. This test uses F test to reject or accept the hypothesis “the fitted model adequately describes the data”. If the F value is high so, the p value is small (< 0.05) then we reject the hypothesis. Under this condition, the model does not reflect the data adequately (Myers et al., 2009).

4.3. Design of experiment and results

Design of Experiment (DOE) is a design tool that makes changes to the factors (input variables) to determine their effect on the responses. It not only identifies the significant factors that affect the response, but also how these factors affect the response (Montgomery, 2001).

In this study, a three-level, three-factorial Box–Behnken experimental design is used to evaluate the effects of selected independent variables on the responses to characterize the hybrid energy system and to optimize the procedure (Myers et al., 2009). This design is suitable for exploration of second-order (quadratic) response surfaces and for construction of second-order polynomial models, thus helping to optimize them by a small number of experimental runs (Box and Draper, 1987; Myers et al., 2009). Box–Behnken experimental design is an orthogonal design. Therefore, the factor levels are evenly spaced and coded for low, medium, and high settings, as −1, 0 and +1 (Montgomery, 2001; Myers et al., 2009). For the three-level, three-factorial Box–Behnken experimental design, a total of 15 experimental runs, shown in Table 6, are needed. Here, x1, x2and x3 are the factors that could affect the cost function as given in (6). -1, 0 and +1 show the coded variables of these factor levels. -1 is the low level, 0 is the middle level and +1 is the high level of the factors. In the middle level, three independent replications, namely 13th, 14th, and 15th experiments in Table 6 are needed.

Experiment	Factor and factor level
x1	x2	x3
1	−1	−1	0
2	−1	+1	0
3	+1	−1	0
4	+1	+1	0
5	−1	0	−1
6	−1	0	+1
7	+1	0	−1
8	+1	0	+1
9	0	−1	−1
10	0	−1	+1
11	0	+1	−1
12	0	+1	+1
13	0	0	0
14	0	0	0
15	0	0	0

Table 6.

Experiments and coded levels for Box–Behnken design

In the design model, three factors are chosen as PV size, As, wind turbine rotor swept area, Aw, and the battery capacity, BC. A list of factors and their levels are provided in Table 7. The levels are 0 and 10, for -1 and +1 levels of PV size; 17 and 37 for -1 and +1 levels of wind turbine rotor swept area; and 10 and 50 for -1 and +1 levels of battery capacity. Since the amount of wind energy in this area is greater than the solar energy, the interval between the low and the high level of the wind turbine rotor swept area is greater than the interval of the low and the high level of the PV size. This also means that at the optimum point wind turbine rotor swept area tends to be greater than the PV size.

Factor	Level
−1	0	+1
x1 : PV (m^2)	0	5	10
x2 : wind (m^2)	17	27	37
x3 : BC (kWh)	10	30	50

Table 7.

Factors and factor levels used in Box–Behnken experimental design

4.4. Second-order polynomial (quadratic) model fitting and RSM results

The general representation of second order regression model of the design is shown by (8).

$$Y=b_0+b_1{x}_1+b_2{x}_2+b_3{x}_3+b_4{x}_1{x}_2+b_5{x}_2{x}_3+b_6{x}_1{x}_3+b_7{x}_1^2+b_8{x}_2^2+b_9{x}_3^2$$ E8

where Y is the selected response; b0–b9 are the regression coefficients and x1-x3 are the factors. This model is used to estimate the relationship between the cost function, Y, and the three independent factors, PV size, x1, wind turbine rotor swept area, x2, and battery capacity, x3. Here, b1, b2, b3 coefficients denote the main effect of factors x1, x2and x3, respectively. Besides, b4 denotes the interaction between factors x1, x2; b5 denotes the interaction between factors x2, x3, and b6 denotes the interaction between factors x1, x3. Finally, b7, b8, b9 denote the quadratic effect of factors x1, x2and x3, respectively.

In order to fit the metamodel given by (8), 15 experiments with three independent replications in the middle are utilized as illustrated in Table 6. The validity of the fitted model is tested by computing a lack-of-fit, F-test and residual analysis. For fitting the metamodel and the optimization procedure, a software called Design Expert 7.1 a popular statistic software package is used (Montgomery, 2001; Myers et al., 2009).

As a result, the fitted standardized metamodel for the hybrid system is given by (8). The model is found to be significant at 95% confidence level by the F-test. In addition, the model does not exhibit lack-of-fit (p>0.05). The lack-of-fit test measures the failure of the model to represent data in the experimental domain at points that are not included in the regression. If a model is significant, meaning that the model contains one or more important terms, and the model does not suffer from lack-of-fit, does not necessarily mean that the model is a good one. If the experimental environment is quite noisy or some important variables are left out of the experiment, then it is possible that the portion of the variability in the data not explained by the model, also called the residual, could be large. Thus, a measure of the model’s overall performance referred to as the coefficient of determination and denoted by R2 must be considered. The value R2 quantifies goodness of fit. It is a fraction between 0 and 1, and has no units. Higher values indicate that the model fits the data better. At the same time, adjusted R2 allowing for the degrees of freedom associated with the sums of the squares is also considered in the lack-of-fit test, which should be an approximate value of R2.

$$\widehat{Y}=38812.37-5008.24x_2+2816x_1{x}_2+10638.18x_2^2+3393.62x_3^2$$ E9

R2 and adjusted R2 are calculated as 0.984 and 0.956, respectively. If adjusted R2 is significantly lower than R2, it normally means that one or more explanatory variables are missing. Here, the two R2 values are not significantly different, and the normal probability plots of residuals do not show evidence of strong departures from normality as depicted in Figure 7. Therefore, the overall second-order metamodel, as expressed by (9), for the response measure is significant and adequate.

In general, estimated standardized metamodel coefficients provide two types of information. The magnitude of a coefficient indicates how important that particular effect is and its sign indicates whether the factor has a positive or a negative effect on the response.

The estimated cost function, $$\widehat{Y}$$ , of the hybrid energy system obtained by RSM in (9) indicates that there is no main effect of PV size and battery capacity on the cost, while wind turbine rotor swept area, x2, significantly affects the response variable (p<0.05). Instead of the main effect of PV size, there is a significant interaction term of PV size and wind turbine

Figure 7.

Normal probability plot of residuals

rotor swept area, x1x2. Besides, the quadratic effects of x2 and x3 are statistically significant (p<0.05). It should be noticed that only the second factor affects the cost negatively while the others affect it positively.

Design Expert 7.1 RSM result indicates that the near optimum value of the fitted response surface, X0 = (x1, x2, x3), is (-0.21, 0.24, 0.096) as coded variable and (3.95, 29.4, 37.55) as the real value yielding a predicted mean response of $$Y_0=$$ $37,033.9 (CT), a minimum in the experimental region. When we run the simulation model at the optimum point the total energy met by the auxiliary energy source in 175,200 hours is obtained as 2,043.85 kWh, which means a shortage cost of $1,021.925. This shortage cost can be interpreted as a small portion in the whole hybrid system cost. Figure 8, illustrates the 3D surface graph of cost response at the optimum battery capacity, 31.92 kWh, as a function of two factors, PV size, x1, and wind turbine rotor swept area, x2.

Figure 8.

Response surface of cost at fixed (optimum) battery capacity as a function of PV size, x1, and wind turbine rotor swept area, x2.

From this figure, while the battery capacity is fixed at 37.55 kWh, the trade-offs between PV size and wind turbine rotor swept area can be seen easily. And also, by this figure, the coded optimum values can be seen around the minimum point of the convex shape.

5. OptQuest optimization

The OptQuest tool is provided with the student version of the ARENA software with free of charge (Kelton et al., 2004). Like all practical simulation optimization methods, OptQuest is also an iterative heuristic. It treats the simulation model as a black box by observing only the input and the output of the simulation model (Kleijnen and Wan, 2007). It is an optimization tool that combines the meta-heuristics of tabu search, neural networks, and scatter search into a single search heuristic.

OptQuest requires the specification of lower, suggested, and upper values for the variables that are to be optimized. The suggested values determine the starting point of the decision variables. This choice affects the efficiency and effectiveness of the search. In practice, the simulation analyst may run the optimization tool twice so that he/she can use the current solution as a suggested value in the second run. We also run the OptQuest twice in order to determine the region of the initial points. Then, the obtained optimum values by the first run are used as suggested values for the second run of the optimization procedure. We thus narrow the search space based on the suggested values obtained from the first run.

Moreover, OptQuest requires the selection of the (random) simulation output that is the goal or objective variable to be minimized. In our problem, we select the minimization of the hybrid system cost as objective. And, OptQuest allows the user to explicitly define integer and linear constraints on the deterministic simulation inputs. Here, we do not specify any constraints, on the decision variables.

OptQuest allows different precision criteria for both the objective and the constrained simulation outputs. For example, it selects the number of replicates such that the halfwidth of the 95% confidence interval for the average output is within a user-selected percentage of the true mean. Because the pre-defined half-width could be obtained by, we use a fixed number of replications, m = 5 (Law, 2007). And, we select 95% confidence interval. Although, OptQuest allows alternative stopping criteria, for example, stop the search either after 300 min (5 h) or after 500 ‘nonimproving solutions, we use a tolerance value of 0.1 to determine when two solutions are equal and to stop the optimization automatically. Figure 9 shows an option example for OptQuest tool.

We optimize the same system using OptQuest. The optimum total cost of the system is obtained as $32,962.5 (CT) which corresponds to (3.04, 32.5, 33.67) for PV size, wind turbine rotor swept area, and battery capacity, respectively.

It should be noted that we obtained better result by OptQuest than RSM. This is probably because OptQuest integrates several heuristics into a single heuristic. Also, running the simulation twice might narrow the search space by getting closer to the optimum point. The only disadvantage using OptQuest is its time consuming property compared to RSM. Also, because the problem is a continuous type optimization in another word, the decision variables - PV size, wind turbine rotor swept area, and battery capacity- can get any numerical values between their levels, the OptQuest tries enormous number of alternative designs while running to find out the near optimum. Therefore, RSM is much faster than the OptQuest.

Figure 9.

OptQuest tool options

6. Conclusion

In this study, size of a PV-wind hybrid energy conversion system with battery storage is optimized using two simulation based optimization methods - RSM and OptQuest. The hybrid system is based on an hourly operating cost to meet the electricity demand of a GSM base station located in western Turkey. First, we fit the probabilistic distributions for solar radiation, wind speed and electricity consumption of the GSM base station using the input analyzer in ARENA – a commercial simulation software. Second, we complete a MC simulation and model a detailed hybrid energy system in ARENA. Third, the system’s cost is optimized by considering the size of PV, wind turbine rotor swept area and battery capacity as decision variables, using RSM and the OptQuest tool in ARENA. Finally, we noticed that we could obtain better result by OptQuest than RSM corresponding to $32,962.5 cost with (3.04, 32.5, 33.67) PV size(m2), wind turbine rotor swept area(m2), and battery capacity (kWh), respectively.